cosmolkit-core 0.2.7

//! Subgraph isomorphism matching (VF2) for molecule pattern matching.
//!
//! ## RDKit provenance (protocol: dev/source_reproduction_protocol.md)
//!
//! This module reproduces RDKit's substructure matching from:
//! - `third_party/rdkit/Code/GraphMol/Substruct/vf2.hpp` (~682 lines C++)
//! - `third_party/rdkit/Code/GraphMol/Substruct/SubstructMatch.cpp` (~735 lines C++)
//!
//! The VF2 algorithm implementation is adapted from vflib-2.0 by P. Foggia,
//! extensively modified by Greg Landrum, ported to Rust with depth-based
//! term_1/term_2 tracking (BackTrack decrements counters instead of
//! recomputing from scratch).
//!
//! ## Marker convention
//!
//! Each copied C++ block below uses the two-axis status marker:
//! - RDKit✔️✔️: fully reproduced behavior and performance
//! - RDKit✔️❌: functionally correct, but with a known performance gap
//! - RDKit❗✔️: approximately implemented, perf-equivalent
//! - RDKit❌❌: not yet ported

use crate::chemistry::forcefield::crystalff::build_crystalff_query_molecule;
use crate::search::query::{
    QueryMatchContext, atom_predicate_matches_with_context, bond_predicate_matches_with_context,
    build_query_match_context,
};
use crate::{Atom, AtomQueryPredicate, Bond, BondOrder, BondQueryPredicate, Molecule};
use std::collections::BTreeMap;
use std::env;
use std::time::Instant;

fn row61_torsion349_substruct_trace_enabled(
    mol: &Molecule,
    query: &Molecule,
    params: &SubstructMatchParams,
) -> bool {
    env::var("RDKIT_ROW61_TRACE").ok().as_deref() == Some("1")
        && mol.num_atoms() == 26
        && query.num_atoms() == 5
        && query.num_bonds() == 4
        && !params.uniquify
}

// ---------------------------------------------------------------------------
// Result types
// ---------------------------------------------------------------------------

/// Result of a single substructure match.
#[derive(Debug, Clone, PartialEq, Eq)]
pub struct SubstructMatchResult {
    /// Mapping from query atom index to molecule atom index.
    pub atom_mapping: Vec<usize>,
    /// Mapping from query bond index to molecule bond index.
    pub bond_mapping: Vec<usize>,
}

/// Parameters controlling substructure matching behaviour.
#[derive(Debug, Clone)]
pub struct SubstructMatchParams {
    /// Maximum number of matches to return (default: 1000).
    pub max_matches: usize,
    /// Whether to uniquify results (default: true).
    pub uniquify: bool,
}

type RecursiveQueryMatchCache = BTreeMap<String, Vec<bool>>;

impl Default for SubstructMatchParams {
    fn default() -> Self {
        Self {
            max_matches: 1000,
            uniquify: true,
        }
    }
}

// ---------------------------------------------------------------------------
// Internal minimum-degree graph representation
// ---------------------------------------------------------------------------

/// Minimal adjacency info needed for VF2.
///
/// `nbrs[i]` is a slice into `edges`.
#[derive(Debug, Clone)]
struct Vf2Graph {
    n_atoms: usize,
    n_bonds: usize,
    /// For each atom index, the neighbor indices and bond ids.
    adjacency: Vec<Vec<(usize, usize)>>, // (neighbor_atom_index, bond_index)
}

/// Build a VF2-compatible adjacency view from a molecule.
///
/// The C++ code iterates `out_edges` via Boost graph. We pre-build adjacency
/// once and use raw index lookups.
fn build_vf2_graph(mol: &Molecule) -> Vf2Graph {
    // RDKit source (implicit in vf2.hpp usage of out_edges):
    //   The VF2 state stores Graph *g1, *g2 and calls:
    //     boost::out_edges(node, *g)
    //     boost::out_degree(node, *g)
    //     boost::adjacent_vertices(node, *g)
    //   These are all O(1) in Boost adjacency_list.
    //
    // RDKit✔️❌: We build a flat adjacency Vec<(usize, usize)> per atom.
    //   This adds a one-time O(V+E) allocation vs the Boost inline storage,
    //   but lookups are O(degree) which matches the original hot-path cost.
    let n_atoms = mol.num_atoms();
    let mut adjacency: Vec<Vec<(usize, usize)>> = vec![Vec::new(); n_atoms];
    for (bond_idx, bond) in mol.bonds().iter().enumerate() {
        let b = bond.begin().index();
        let e = bond.end().index();
        adjacency[b].push((e, bond_idx));
        adjacency[e].push((b, bond_idx));
    }
    Vf2Graph {
        n_atoms,
        n_bonds: mol.num_bonds(),
        adjacency,
    }
}

impl Vf2Graph {
    fn out_degree(&self, node: usize) -> usize {
        self.adjacency[node].len()
    }

    fn out_edges(&self, node: usize) -> &[(usize, usize)] {
        &self.adjacency[node]
    }
}

// ---------------------------------------------------------------------------
// Atom and bond matching functors
// ---------------------------------------------------------------------------

// RDKit source (SubstructMatch.cpp):
//   class AtomLabelFunctor {
//    public:
//     AtomLabelFunctor(const ROMol &query, const ROMol &mol,
//                      const SubstructMatchParameters &ps)
//         : d_query(query), d_mol(mol), d_params(ps) {};
//     bool operator()(unsigned int i, unsigned int j) const {
//       bool res = false;
//       if (d_params.useChirality) {
//         const Atom *qAt = d_query.getAtomWithIdx(i);
//         if (qAt->getChiralTag() == Atom::CHI_TETRAHEDRAL_CW ||
//             qAt->getChiralTag() == Atom::CHI_TETRAHEDRAL_CCW) {
//           const Atom *mAt = d_mol.getAtomWithIdx(j);
//           if (!d_params.specifiedStereoQueryMatchesUnspecified &&
//               mAt->getChiralTag() != Atom::CHI_TETRAHEDRAL_CW &&
//               mAt->getChiralTag() != Atom::CHI_TETRAHEDRAL_CCW) {
//             return false;
//           }
//         }
//       }
//       res = atomCompat(d_query[i], d_mol[j], d_params);
//       return res;
//     }
//    private:
//     const ROMol &d_query;
//     const ROMol &d_mol;
//     const SubstructMatchParameters &d_params;
//   };
//
// RDKit✔️❌: AtomLabelFunctor ported as a plain function. useChirality
//   handling is preserved but chirality is not yet wired into the matching.
//   The atomCompat logic is inlined for the basic case (atomic number,
//   aromaticity, isotope, formal charge). Full QueryNode eval is deferred.

/// RDKit✔️❌: Basic atom compatibility check matching RDKit's atomCompat.
///
/// RDKit `atomCompat` checks: atomic number, aromaticity, isotope,
/// formal charge, and query atom queries. We implement the non-query
/// subset here. Full SMARTS/QueryNode evaluation is blocked on
/// MatchSubqueries porting.
fn atom_matches(query_atom: &Atom, query_mol: &Molecule, mol_atom: &Atom, mol: &Molecule) -> bool {
    // RDKit source (atomCompat, approx):
    //   if (query->hasQuery()) {
    //     return query->Match(mol);
    //   }
    //   if (query->getAtomicNum() != mol->getAtomicNum() &&
    //       query->getAtomicNum() > 0) return false;
    //   if (query->getIsAromatic() != mol->getIsAromatic() &&
    //       query->getIsAromatic()) return false;
    //   if (query->getIsotope() && query->getIsotope() != mol->getIsotope())
    //     return false;
    //   if (query->getFormalCharge() != mol->getFormalCharge() &&
    //       query->getFormalCharge() != 0) return false;
    //   return true;

    // RDKit❗✔️: If the query atom has a query predicate tree, evaluate it
    // through the non-cached path. The recursive-SMARTS cached path is wired
    // separately through `atom_matches_with_recursive_cache()`.
    if let Some(query_node) = query_atom.query() {
        let query_ctx = build_query_match_context(mol);
        return match query_node {
            crate::QueryNode::Predicate(pred) => {
                atom_query_predicate_matches_for_substruct(mol_atom, pred, mol, None, &query_ctx)
            }
            crate::QueryNode::And(children) => children.iter().all(|child| match child {
                crate::QueryNode::Predicate(pred) => atom_query_predicate_matches_for_substruct(
                    mol_atom, pred, mol, None, &query_ctx,
                ),
                _ => evaluate_atom_query(child, mol_atom, mol, None, &query_ctx),
            }),
            crate::QueryNode::Or(children) => children.iter().any(|child| match child {
                crate::QueryNode::Predicate(pred) => atom_query_predicate_matches_for_substruct(
                    mol_atom, pred, mol, None, &query_ctx,
                ),
                _ => evaluate_atom_query(child, mol_atom, mol, None, &query_ctx),
            }),
            crate::QueryNode::Not(child) => {
                !evaluate_atom_query(child, mol_atom, mol, None, &query_ctx)
            }
        };
    }

    let q_an = query_atom.atomic_number();
    let m_an = mol_atom.atomic_number();

    // If query specifies a (non-dummy) atomic number, it must match.
    if q_an != 0 && q_an != m_an {
        return false;
    }

    // Aromaticity check: if query atom is aromatic, mol atom must also be aromatic.
    if query_atom.is_aromatic() && !mol_atom.is_aromatic() {
        return false;
    }

    // Isotope check: if query specifies an isotope, it must match.
    if let Some(q_iso) = query_atom.isotope() {
        if mol_atom.isotope() != Some(q_iso) {
            return false;
        }
    }

    // Formal charge check: if query has a non-zero charge, it must match.
    let q_charge = query_atom.formal_charge();
    if q_charge != 0 && q_charge != mol_atom.formal_charge() {
        return false;
    }

    let _ = query_mol;
    true
}

fn recursive_smarts_root_matches(
    atom: &Atom,
    recursive_smarts: &str,
    mol: &Molecule,
    recursive_cache: Option<&RecursiveQueryMatchCache>,
) -> bool {
    // BEGIN RDKIT CPP FUNCTION detail::MatchSubqueries / detail::RecursiveMatcher (SubstructMatch.cpp:613-712 via QueryOps.h RecursiveStructureQuery)
    // RDKit❗✔️: RecursiveStructureQuery(ROMol const *query, unsigned int serialNumber = 0)
    // RDKit❗✔️:   unsigned int res = RecursiveMatcher(mol, *queryMol, matchStarts,
    // RDKit❗✔️:                                       subqueryMap, params, locked);
    // RDKit❗✔️:   if (res) {
    // RDKit❗✔️:     for (int &matchStart : matchStarts) {
    // RDKit❗✔️:       rsq->insert(matchStart);
    // RDKit❗✔️:     }
    // RDKit❗✔️:   }
    // RDKit❗✔️:   if (!query.hasProp(common_properties::_queryRootAtom)) {
    // RDKit❗✔️:     matches.push_back(pairs.begin()->second);
    // RDKit❗✔️:   }
    if let Some(cache) = recursive_cache
        && let Some(match_starts) = cache.get(recursive_smarts)
    {
        return match_starts
            .get(atom.id().index())
            .copied()
            .unwrap_or(false);
    }

    let Some(inner) = recursive_smarts
        .strip_prefix("$(")
        .and_then(|value| value.strip_suffix(')'))
    else {
        return false;
    };
    let Ok(query) = build_crystalff_query_molecule(inner) else {
        return false;
    };
    substruct_match_impl(
        mol,
        &query,
        &SubstructMatchParams {
            max_matches: 1000,
            uniquify: false,
        },
    )
    .into_iter()
    .any(|matched| matched.atom_mapping.first().copied() == Some(atom.id().index()))
}

fn atom_query_predicate_matches_for_substruct(
    atom: &Atom,
    pred: &AtomQueryPredicate,
    mol: &Molecule,
    recursive_cache: Option<&RecursiveQueryMatchCache>,
    query_ctx: &QueryMatchContext,
) -> bool {
    match pred {
        AtomQueryPredicate::RecursiveSmarts(recursive_smarts) => {
            recursive_smarts_root_matches(atom, recursive_smarts, mol, recursive_cache)
        }
        _ => atom_predicate_matches_with_context(atom, pred, mol, query_ctx),
    }
}

/// RDKit❗✔️: Partial evaluation of an atom query node.
///
/// Only the basic predicates (AtomicNumber, IsAromatic, FormalCharge,
/// Isotope) are implemented. Recursive, RGroup, and SMARTS queries
/// are not yet supported.
fn evaluate_atom_query(
    query: &crate::QueryNode<AtomQueryPredicate>,
    atom: &Atom,
    mol: &Molecule,
    recursive_cache: Option<&RecursiveQueryMatchCache>,
    query_ctx: &QueryMatchContext,
) -> bool {
    match query {
        crate::QueryNode::Predicate(pred) => {
            atom_query_predicate_matches_for_substruct(atom, pred, mol, recursive_cache, query_ctx)
        }
        crate::QueryNode::And(children) => children
            .iter()
            .all(|c| evaluate_atom_query(c, atom, mol, recursive_cache, query_ctx)),
        crate::QueryNode::Or(children) => children
            .iter()
            .any(|c| evaluate_atom_query(c, atom, mol, recursive_cache, query_ctx)),
        crate::QueryNode::Not(child) => {
            !evaluate_atom_query(child, atom, mol, recursive_cache, query_ctx)
        }
    }
}

// RDKit source (SubstructMatch.cpp):
//   class BondLabelFunctor {
//    public:
//     BondLabelFunctor(const ROMol &query, const ROMol &mol,
//                      const SubstructMatchParameters &ps)
//         : d_query(query), d_mol(mol), d_params(ps) {};
//     bool operator()(MolGraph::edge_descriptor i,
//                     MolGraph::edge_descriptor j) const {
//       if (d_params.useChirality) {
//         const Bond *qBnd = d_query[i];
//         if (qBnd->getBondType() == Bond::DOUBLE &&
//             qBnd->getStereo() > Bond::STEREOANY) {
//           const Bond *mBnd = d_mol[j];
//           if (mBnd->getBondType() == Bond::DOUBLE &&
//               !d_params.specifiedStereoQueryMatchesUnspecified &&
//               mBnd->getStereo() <= Bond::STEREOANY) {
//             return false;
//           }
//         }
//       }
//       bool res = bondCompat(d_query[i], d_mol[j], d_params);
//       return res;
//     }
//    private:
//     const ROMol &d_query;
//     const ROMol &d_mol;
//     const SubstructMatchParameters &d_params;
//   };

/// RDKit✔️❌: Basic bond compatibility check matching RDKit's bondCompat.
///
/// Handles bond order matching with aromatic fallback rules.
/// Query bonds via BondQueryPredicate are deferred.
fn bond_matches(query_bond: &Bond, query_mol: &Molecule, mol_bond: &Bond, mol: &Molecule) -> bool {
    // RDKit source (bondCompat, approx):
    //   if (qBnd->hasQuery()) return qBnd->Match(mBnd);
    //   if (qBnd->getIsAromatic() && mBnd->getIsAromatic()) return true;
    //   return qBnd->getBondType() == mBnd->getBondType();

    // RDKit❗✔️: If query bond has a query predicate tree, try to evaluate.
    if let Some(query_node) = query_bond.query() {
        let query_ctx = build_query_match_context(mol);
        return evaluate_bond_query(query_node, mol_bond, mol, &query_ctx);
    }

    let q_aromatic = query_bond.is_aromatic();
    let m_aromatic = mol_bond.is_aromatic();

    // Both aromatic: match regardless of bond order.
    if q_aromatic && m_aromatic {
        return true;
    }

    // Either could be aromatic, use order comparison.
    // RDKit rule: aromatic query bond matches aromatic mol bond (handled above)
    // or single bond in mol (for Kekule forms).
    // Also: query single bond matches aromatic mol bond.
    let q_order = query_bond.order();
    let m_order = mol_bond.order();

    if q_order == m_order {
        return true;
    }

    // Aromatic/Single interchange.
    if q_aromatic && m_order == BondOrder::Single {
        return true;
    }
    if m_aromatic && q_order == BondOrder::Single {
        return true;
    }

    let _ = query_mol;
    false
}

/// RDKit❗✔️: Partial evaluation of a bond query node.
fn evaluate_bond_query(
    query: &crate::QueryNode<BondQueryPredicate>,
    bond: &Bond,
    mol: &Molecule,
    query_ctx: &QueryMatchContext,
) -> bool {
    match query {
        crate::QueryNode::Predicate(pred) => {
            bond_predicate_matches_with_context(bond, pred, mol, query_ctx)
        }
        crate::QueryNode::And(children) => children
            .iter()
            .all(|c| evaluate_bond_query(c, bond, mol, query_ctx)),
        crate::QueryNode::Or(children) => children
            .iter()
            .any(|c| evaluate_bond_query(c, bond, mol, query_ctx)),
        crate::QueryNode::Not(child) => !evaluate_bond_query(child, bond, mol, query_ctx),
    }
}

// ---------------------------------------------------------------------------
// VF2 State Machine
// ---------------------------------------------------------------------------
//
// ## RDKit source reproduction: vf2.hpp
//
// The following section reproduces the VF2SubState class from vf2.hpp.
// The C++ code is shown as verbatim comments with RDKit markers.
//
// ### Key design differences from RDKit:
//
// 1. `core_1`/`core_2`: Same role — mapping from query atom idx → mol atom idx
//    and vice versa. Uses `Option<usize>` instead of NULL_NODE sentinel.
//
// 2. `term_1`/`term_2`: Stores the core_len *depth* at which each atom was
//    added to the terminal set, exactly as in vf2.hpp. BackTrack decrements
//    counters keyed by depth, not recomputes from scratch.
//
// 3. No shared_ptr copy semantics: VF2SubState in RDKit uses COW with
//    `share_count`. Rust's Clone+Vf2State avoids raw pointer sharing.
//    This means each VF2 recursive branch owns its state, which is
//    semantically correct but allocates O(depth * n) instead of
//    O(n) shared storage. For typical molecule sizes (<1000 atoms) this
//    is negligible; for very large searches the COW approach could be
//    reinstated with Arc<Vec<NodeId>>.
//
// 4. No boost graph: hand-rolled Vf2Graph adjacency.

// RDKit source (vf2.hpp):
//   typedef std::uint32_t node_id;
//   const node_id NULL_NODE = 0xFFFFFFFF;

type NodeId = usize;
const NULL_NODE: NodeId = usize::MAX;

// RDKit source (vf2.hpp):
//   template <class Graph>
//   struct Pair {
//     node_id n1, n2;
//     bool hasiter{false};
//     RDK_ADJ_ITER nbrbeg, nbrend;
//     Pair() : n1(NULL_NODE), n2(NULL_NODE) {}
//   };

#[derive(Debug, Clone)]
struct Vf2Pair {
    n1: NodeId,
    n2: NodeId,
    hasiter: bool,
    /// VF2+ source atom in the mol graph (g2) whose adjacency drives the
    /// neighbor iterator.
    nbr_node: NodeId,
    /// VF2+ neighbor iterator over mol graph (g2) neighbors.
    nbr_cursor: usize,
    nbr_end: usize,
}

impl Vf2Pair {
    fn new() -> Self {
        Self {
            n1: NULL_NODE,
            n2: NULL_NODE,
            hasiter: false,
            nbr_node: NULL_NODE,
            nbr_cursor: 0,
            nbr_end: 0,
        }
    }
}

// RDKit source (vf2.hpp):
//   /**
//    * The ordering by in/out degree
//    */
//   static bool nodeInfoComp1(const NodeInfo &a, const NodeInfo &b) {
//     if (a.out < b.out) { return true; }
//     if (a.out > b.out) { return false; }
//     if (a.in < b.in) { return true; }
//     if (a.in > b.in) { return false; }
//     return false;
//   }

#[derive(Debug, Clone, Copy)]
struct NodeInfo {
    id: usize,
    in_deg: usize,
    out_deg: usize,
}

// RDKit✔️✔️: nodeInfoComp1 — sort ascending by out-degree then in-degree.
fn node_info_cmp1(a: &NodeInfo, b: &NodeInfo) -> std::cmp::Ordering {
    // RDKit✔️✔️: if (a.out < b.out) { return true; }
    // RDKit✔️✔️: if (a.out > b.out) { return false; }
    // RDKit✔️✔️: if (a.in < b.in) { return true; }
    // RDKit✔️✔️: if (a.in > b.in) { return false; }
    // RDKit✔️✔️: return false;
    a.out_deg
        .cmp(&b.out_deg)
        .then_with(|| a.in_deg.cmp(&b.in_deg))
}

// RDKit source (vf2.hpp):
//   static int nodeInfoComp2(const NodeInfo &a, const NodeInfo &b) {
//     if (!a.in && b.in) return 1;
//     if (a.in && !b.in) return -1;
//     if (a.out < b.out) return -1;
//     if (a.out > b.out) return 1;
//     if (a.in < b.in) return -1;
//     if (a.in > b.in) return 1;
//     return 0;
//   }

// RDKit✔️✔️: nodeInfoComp2 — sort ascending by frequency (out=run count) then
//   valence (in=valence sum), except zero-valence nodes sort after non-zero.
fn node_info_cmp2(a: &NodeInfo, b: &NodeInfo) -> std::cmp::Ordering {
    // RDKit✔️✔️: if (!a.in && b.in) return 1;
    // RDKit✔️✔️: if (a.in && !b.in) return -1;
    if a.in_deg == 0 && b.in_deg != 0 {
        return std::cmp::Ordering::Greater;
    }
    if a.in_deg != 0 && b.in_deg == 0 {
        return std::cmp::Ordering::Less;
    }
    // RDKit✔️✔️: if (a.out < b.out) return -1;
    // RDKit✔️✔️: if (a.out > b.out) return 1;
    // RDKit✔️✔️: if (a.in < b.in) return -1;
    // RDKit✔️✔️: if (a.in > b.in) return 1;
    // RDKit✔️✔️: return 0;
    a.out_deg
        .cmp(&b.out_deg)
        .then_with(|| a.in_deg.cmp(&b.in_deg))
}

// RDKit source (vf2.hpp), SortNodesByFrequency:
//   Sorts the nodes of a graphs, returning a heap-allocated vector
//   with the node ids in the proper orders.
//   The sorting criterion takes into account:
//     1 - The number of nodes with the same in/out degree.
//     2 - The valence of the nodes.
//   The nodes at the beginning of the vector are the most singular,
//   from which the matching should start.

// RDKit✔️✔️: SortNodesByFrequency — returns query node ordering for VF2.
fn sort_nodes_by_frequency(g: &Vf2Graph) -> Vec<NodeId> {
    // RDKit✔️✔️: template <class Graph>
    // RDKit✔️✔️: node_id *SortNodesByFrequency(const Graph *g) {
    // RDKit✔️✔️:   std::vector<NodeInfo> vect;
    // RDKit✔️✔️:   vect.reserve(boost::num_vertices(*g));
    // RDKit✔️✔️:   typename Graph::vertex_iterator bNode, eNode;
    // RDKit✔️✔️:   boost::tie(bNode, eNode) = boost::vertices(*g);
    // RDKit✔️✔️:   while (bNode != eNode) {
    // RDKit✔️✔️:     NodeInfo t;
    // RDKit✔️✔️:     t.id = vect.size();
    // RDKit✔️✔️:     t.in = boost::out_degree(*bNode, *g);
    // RDKit✔️✔️:     t.out = boost::out_degree(*bNode, *g);
    // RDKit✔️✔️:     vect.push_back(t);
    // RDKit✔️✔️:     ++bNode;
    // RDKit✔️✔️:   }
    // RDKit✔️✔️:   std::sort(vect.begin(), vect.end(), nodeInfoComp1);
    let mut vect: Vec<NodeInfo> = (0..g.n_atoms)
        .map(|i| {
            let deg = g.out_degree(i);
            NodeInfo {
                id: i,
                in_deg: deg,
                out_deg: deg,
            }
        })
        .collect();
    vect.sort_unstable_by(node_info_cmp1);

    // RDKit✔️✔️:   unsigned int run = 1;
    // RDKit✔️✔️:   for (unsigned int i = 0; i < vect.size(); i += run) {
    // RDKit✔️✔️:     for (run = 1; i+run < vect.size() && vect[i+run].in == vect[i].in
    // RDKit✔️✔️:             && vect[i+run].out == vect[i].out; ++run) { ; }
    // RDKit✔️✔️:     for (unsigned int j = 0; j < run; ++j) {
    // RDKit✔️✔️:       vect[i+j].in += vect[i+j].out;  // in=valence sum
    // RDKit✔️✔️:       vect[i+j].out = run;            // out=frequency
    // RDKit✔️✔️:     }
    // RDKit✔️✔️:   }
    let mut i = 0;
    while i < vect.len() {
        let mut run = 1;
        while i + run < vect.len()
            && vect[i + run].in_deg == vect[i].in_deg
            && vect[i + run].out_deg == vect[i].out_deg
        {
            run += 1;
        }
        for j in 0..run {
            vect[i + j].in_deg += vect[i + j].out_deg; // valence sum
            vect[i + j].out_deg = run; // frequency
        }
        i += run;
    }

    // RDKit✔️✔️:   std::sort(vect.begin(), vect.end(), nodeInfoComp2);
    vect.sort_unstable_by(node_info_cmp2);

    // RDKit✔️✔️:   node_id *nodes = new node_id[vect.size()];
    // RDKit✔️✔️:   for (unsigned int i = 0; i < vect.size(); ++i) {
    // RDKit✔️✔️:     nodes[i] = vect[i].id;
    // RDKit✔️✔️:   }
    // RDKit✔️✔️:   return nodes;
    vect.iter().map(|ni| ni.id).collect()
}

// RDKit source (vf2.hpp), VF2SubState class:
//   template <class Graph, class VertexCompatible, class EdgeCompatible,
//             class MatchChecking>
//   class VF2SubState {
//    private:
//     Graph *g1, *g2;
//     VertexCompatible &vc;
//     EdgeCompatible &ec;
//     MatchChecking &mc;
//     unsigned int n1, n2;
//     unsigned int core_len;
//     unsigned int t1_len;
//     unsigned int t2_len;  // Core nodes are also counted by these...
//     node_id *core_1;
//     node_id *core_2;
//     node_id *term_1;
//     node_id *term_2;
//     node_id *order;
//     long *share_count;
//     int *vs_compared;

/// RDKit❗✔️: VF2 subgraph isomorphism state.
///
/// g1 = query graph, g2 = molecule graph.
/// core_1[i] = mapping from query atom i -> mol atom j (or None).
/// core_2[j] = mapping from mol atom j -> query atom i (or None).
/// term_1[i] = depth (core_len) when atom i entered terminal set (0 = not terminal).
/// term_2[j] = same for mol atoms.
struct Vf2SubState<'a> {
    g1: &'a Vf2Graph,
    g2: &'a Vf2Graph,
    n1: usize,
    n2: usize,
    core_len: usize,
    t1_len: usize,
    t2_len: usize,
    core_1: Vec<NodeId>,
    core_2: Vec<NodeId>,
    term_1: Vec<usize>,
    term_2: Vec<usize>,
    order: Option<Vec<NodeId>>,
}

impl<'a> Vf2SubState<'a> {
    // RDKit source (vf2.hpp):
    //   VF2SubState(Graph *ag1, Graph *ag2, VertexCompatible &avc,
    //               EdgeCompatible &aec, MatchChecking &amc, bool sortNodes=false)
    //       : g1(ag1), g2(ag2), vc(avc), ec(aec), mc(amc),
    //         n1(num_vertices(*ag1)), n2(num_vertices(*ag2)) {
    //     if (sortNodes) { order = SortNodesByFrequency(ag1); }
    //     else { order = nullptr; }
    //     core_len = 0; t1_len = 0; t2_len = 0;
    //     core_1 = new node_id[n1]; core_2 = new node_id[n2];
    //     term_1 = new node_id[n1]; term_2 = new node_id[n2];
    //     for (unsigned int i = 0; i < n1; i++) {
    //       core_1[i] = NULL_NODE; term_1[i] = 0;
    //     }
    //     for (unsigned int i = 0; i < n2; i++) {
    //       core_2[i] = NULL_NODE; term_2[i] = 0;
    //     }
    //   }

    /// RDKit✔️❌: Construct VF2 state. `sort_nodes` enables frequency-based
    ///   query node ordering. Performance gap: Rust Vec allocation vs C++
    ///   raw new[] — negligible for typical molecule sizes.
    fn new(g1: &'a Vf2Graph, g2: &'a Vf2Graph, sort_nodes: bool) -> Self {
        let n1 = g1.n_atoms;
        let n2 = g2.n_atoms;
        let order = if sort_nodes {
            Some(sort_nodes_by_frequency(g1))
        } else {
            None
        };

        // RDKit✔️✔️: core_len = 0; t1_len = 0; t2_len = 0;
        // RDKit✔️✔️: core_1[i] = NULL_NODE; term_1[i] = 0;
        // RDKit✔️✔️: core_2[j] = NULL_NODE; term_2[j] = 0;
        Self {
            g1,
            g2,
            n1,
            n2,
            core_len: 0,
            t1_len: 0,
            t2_len: 0,
            core_1: vec![NULL_NODE; n1],
            core_2: vec![NULL_NODE; n2],
            term_1: vec![0usize; n1],
            term_2: vec![0usize; n2],
            order,
        }
    }

    fn debug_order(&self) -> Option<&[NodeId]> {
        self.order.as_deref()
    }

    // RDKit source (vf2.hpp):
    //   bool IsGoal() { return core_len == n1; }
    //   bool IsDead() { return n1 > n2 || t1_len > t2_len; }

    /// RDKit✔️✔️: IsGoal — all query atoms matched.
    fn is_goal(&self) -> bool {
        self.core_len == self.n1
    }

    /// RDKit✔️✔️: IsDead — more query atoms than mol atoms, or more
    ///   terminal query atoms than terminal mol atoms.
    fn is_dead(&self) -> bool {
        self.n1 > self.n2 || self.t1_len > self.t2_len
    }

    // RDKit source (vf2.hpp):
    //   bool NextPair(Pair<Graph> &pair) {
    //     if (pair.n1 == NULL_NODE) { pair.n1 = 0; }
    //     if (pair.n2 == NULL_NODE) { pair.n2 = 0; }
    //     else { pair.n2++; }
    //     ...
    //     if (t1_len > core_len && t2_len > core_len) {
    //       while (pair.n1 < n1 &&
    //              (core_1[pair.n1] != NULL_NODE || term_1[pair.n1] == 0)) {
    //         pair.n1++; pair.n2 = 0;
    //       }
    //       ...
    //     } else if (pair.n1 == 0 && order != nullptr) {
    //       // Optimisation: ...
    //       unsigned int i = 0;
    //       while (i < n1 && core_1[pair.n1 = order[i]] != NULL_NODE) { i++; }
    //       ...
    //     } else {
    //       while (pair.n1 < n1 && core_1[pair.n1] != NULL_NODE) {
    //         pair.n1++; pair.n2 = 0;
    //       }
    //     }
    //     // VF2 Plus iterator ...
    //     if (pair.hasiter) { ... }
    //     else if (t1_len > core_len && t2_len > core_len) {
    //       while (pair.n2 < n2 &&
    //              (core_2[pair.n2] != NULL_NODE || term_2[pair.n2] == 0)) {
    //         pair.n2++;
    //       }
    //     } else {
    //       while (pair.n2 < n2 && core_2[pair.n2] != NULL_NODE) { pair.n2++; }
    //     }
    //     return pair.n1 < n1 && pair.n2 < n2;
    //   }

    /// RDKit✔️❌: NextPair — find the next candidate pair (n1 from query,
    ///   n2 from mol) to try matching.
    ///
    /// Uses terminal-set-based iteration from vf2.hpp, including the VF2+
    /// neighbor iterator that restricts mol-side candidates to neighbors of
    /// the already-mapped terminal predecessor.
    fn next_pair(&self, pair: &mut Vf2Pair) -> bool {
        // RDKit✔️✔️: if (pair.n1 == NULL_NODE) pair.n1 = 0;
        // RDKit✔️✔️: if (pair.n2 == NULL_NODE) pair.n2 = 0;
        // RDKit✔️✔️: else pair.n2++;
        if pair.n1 == NULL_NODE {
            pair.n1 = 0;
        }
        if pair.n2 == NULL_NODE {
            pair.n2 = 0;
        } else {
            pair.n2 += 1;
        }

        // --- Select query node (n1) ---
        // RDKit✔️✔️: if (t1_len > core_len && t2_len > core_len) {
        if self.t1_len > self.core_len && self.t2_len > self.core_len {
            // RDKit✔️✔️: while (pair.n1 < n1 &&
            // RDKit✔️✔️:   (core_1[pair.n1] != NULL_NODE || term_1[pair.n1] == 0)) {
            // RDKit✔️✔️:   pair.n1++; pair.n2 = 0;
            // RDKit✔️✔️: }
            while pair.n1 < self.n1
                && (self.core_1[pair.n1] != NULL_NODE || self.term_1[pair.n1] == 0)
            {
                pair.n1 += 1;
                pair.n2 = 0;
            }
            // RDKit✔️✔️: /* Initialize VF2 Plus neighbor iterator.
            // RDKit✔️✔️:  * The next query node (pair.n1) has been selected from the terminal
            // RDKit✔️✔️:  * set and is therefore adjacent to an already mapped atom (in
            // RDKit✔️✔️:  * core_1). Rather than select pair.n2 from all atoms (0...n2) we can
            // RDKit✔️✔️:  * select it from the neighbors of this mapped atom (0...deg(nbor))
            // RDKit✔️✔️:  * since it must also be adajcent to this mapped atom!
            // RDKit✔️✔️:  */
            // RDKit✔️✔️: if (!pair.hasiter) {
            // RDKit✔️✔️:   boost::tie(n1iter_beg, n1iter_end) =
            // RDKit✔️✔️:       boost::adjacent_vertices(pair.n1, *g1);
            // RDKit✔️✔️:   while (n1iter_beg != n1iter_end && core_1[*n1iter_beg] == NULL_NODE) {
            // RDKit✔️✔️:     ++n1iter_beg;
            // RDKit✔️✔️:   }
            // RDKit✔️✔️:   assert(n1iter_beg != n1iter_end);
            // RDKit✔️✔️:   boost::tie(pair.nbrbeg, pair.nbrend) =
            // RDKit✔️✔️:       boost::adjacent_vertices(core_1[*n1iter_beg], *g2);
            // RDKit✔️✔️:   pair.hasiter = true;
            // RDKit✔️✔️: }
            if !pair.hasiter {
                let mut mapped_terminal_neighbor = NULL_NODE;
                for &(query_neighbor, _) in self.g1.out_edges(pair.n1) {
                    if self.core_1[query_neighbor] != NULL_NODE {
                        mapped_terminal_neighbor = self.core_1[query_neighbor];
                        break;
                    }
                }
                debug_assert_ne!(mapped_terminal_neighbor, NULL_NODE);
                if mapped_terminal_neighbor != NULL_NODE {
                    pair.nbr_node = mapped_terminal_neighbor;
                    pair.nbr_cursor = 0;
                    pair.nbr_end = self.g2.out_edges(mapped_terminal_neighbor).len();
                    pair.hasiter = true;
                }
            }
        } else if pair.n1 == 0 {
            // RDKit✔️✔️: } else if (pair.n1 == 0 && order != nullptr) {
            if let Some(order) = &self.order {
                // RDKit✔️✔️:   unsigned int i = 0;
                // RDKit✔️✔️:   while (i < n1 && core_1[pair.n1 = order[i]] != NULL_NODE) { i++; }
                // RDKit✔️✔️:   if (i == n1) pair.n1 = n1;
                let mut i = 0;
                while i < self.n1 {
                    let candidate = order[i];
                    if self.core_1[candidate] == NULL_NODE {
                        pair.n1 = candidate;
                        break;
                    }
                    i += 1;
                }
                if i == self.n1 {
                    pair.n1 = self.n1;
                }
            } else {
                // RDKit✔️✔️: } else {
                // RDKit✔️✔️:   while (pair.n1 < n1 && core_1[pair.n1] != NULL_NODE) {
                // RDKit✔️✔️:     pair.n1++; pair.n2 = 0;
                // RDKit✔️✔️:   }
                while pair.n1 < self.n1 && self.core_1[pair.n1] != NULL_NODE {
                    pair.n1 += 1;
                    pair.n2 = 0;
                }
            }
        } else {
            // RDKit✔️✔️: } else {
            // RDKit✔️✔️:   while (pair.n1 < n1 && core_1[pair.n1] != NULL_NODE) {
            // RDKit✔️✔️:     pair.n1++; pair.n2 = 0;
            // RDKit✔️✔️:   }
            while pair.n1 < self.n1 && self.core_1[pair.n1] != NULL_NODE {
                pair.n1 += 1;
                pair.n2 = 0;
            }
        }

        // --- Select mol node (n2) ---
        // RDKit✔️✔️: if (pair.hasiter) { ... }
        if pair.hasiter {
            // RDKit✔️✔️: while (pair.nbrbeg < pair.nbrend && core_2[*pair.nbrbeg] != NULL_NODE) {
            // RDKit✔️✔️:   ++pair.nbrbeg;
            // RDKit✔️✔️: }
            let neighbors = self.g2.out_edges(pair.nbr_node);
            while pair.nbr_cursor < pair.nbr_end
                && self.core_2[neighbors[pair.nbr_cursor].0] != NULL_NODE
            {
                pair.nbr_cursor += 1;
            }
            // RDKit✔️✔️: if (pair.nbrbeg < pair.nbrend) {
            // RDKit✔️✔️:   pair.n2 = *pair.nbrbeg;
            // RDKit✔️✔️:   ++pair.nbrbeg;
            // RDKit✔️✔️: } else {
            // RDKit✔️✔️:   pair.n2 = n2;
            // RDKit✔️✔️: }
            if pair.nbr_cursor < pair.nbr_end {
                pair.n2 = neighbors[pair.nbr_cursor].0;
                pair.nbr_cursor += 1;
            } else {
                pair.n2 = self.n2;
            }
        } else if self.t1_len > self.core_len && self.t2_len > self.core_len {
            // RDKit✔️✔️: } else if (t1_len > core_len && t2_len > core_len) {
            // RDKit✔️✔️:   while (pair.n2 < n2 &&
            // RDKit✔️✔️:     (core_2[pair.n2] != NULL_NODE || term_2[pair.n2] == 0)) {
            // RDKit✔️✔️:     pair.n2++;
            // RDKit✔️✔️:   }
            while pair.n2 < self.n2
                && (self.core_2[pair.n2] != NULL_NODE || self.term_2[pair.n2] == 0)
            {
                pair.n2 += 1;
            }
        } else {
            // RDKit✔️✔️: } else {
            // RDKit✔️✔️:   while (pair.n2 < n2 && core_2[pair.n2] != NULL_NODE) { pair.n2++; }
            // RDKit✔️✔️: }
            while pair.n2 < self.n2 && self.core_2[pair.n2] != NULL_NODE {
                pair.n2 += 1;
            }
        }

        // RDKit✔️✔️: return pair.n1 < n1 && pair.n2 < n2;
        let ok = pair.n1 < self.n1 && pair.n2 < self.n2;
        if ok
            && env::var("RDKIT_ROW61_TRACE").ok().as_deref() == Some("1")
            && self.n1 == 5
            && self.n2 == 26
        {
            println!(
                "row61_substruct_next_pair core_len={} t1_len={} t2_len={} hasiter={} n1={} n2={}",
                self.core_len, self.t1_len, self.t2_len, pair.hasiter, pair.n1, pair.n2
            );
        }
        ok
    }

    // RDKit source (vf2.hpp), IsFeasiblePair:
    //   bool IsFeasiblePair(node_id node1, node_id node2) {
    //     assert(node1 < n1); assert(node2 < n2);
    //     assert(core_1[node1] == NULL_NODE); assert(core_2[node2] == NULL_NODE);
    //
    //     // O(1) check for adjacency list
    //     if (boost::out_degree(node1, *g1) > boost::out_degree(node2, *g2)) {
    //       return false;
    //     }
    //     if (!vc(node1, node2)) { return false; }
    //
    //     unsigned int other1, other2;
    //     // Check the out edges of node1
    //     typename Graph::out_edge_iterator bNbrs, eNbrs;
    //     boost::tie(bNbrs, eNbrs) = boost::out_edges(node1, *g1);
    //     while (bNbrs != eNbrs) {
    //       other1 = getOtherIdx(*g1, *bNbrs, node1);
    //       if (core_1[other1] != NULL_NODE) {
    //         other2 = core_1[other1];
    //         typename Graph::edge_descriptor oEdge;
    //         bool found;
    //         boost::tie(oEdge, found) = boost::edge(node2, other2, *g2);
    //         if (!found || !ec(*bNbrs, oEdge)) { return false; }
    //       }
    //       ++bNbrs;
    //     }
    //     return true;
    //   }

    /// RDKit✔️❌: IsFeasiblePair — check if (node1, node2) can be added.
    ///
    /// Performs degree check, vertex compatibility, and edge compatibility
    /// for already-matched neighbors. RDK_VF2_PRUNING (terminal count
    /// pre-check) is not enabled — the C++ code also has it behind an
    /// ifdef that is not defined at the top of vf2.hpp.
    fn is_feasible_pair(
        &self,
        node1: NodeId,
        node2: NodeId,
        atom_fn: &impl Fn(usize, usize) -> bool,
        bond_fn: &impl Fn(usize, usize) -> bool,
    ) -> bool {
        // RDKit✔️✔️: assert(node1 < n1); assert(node2 < n2);
        // RDKit✔️✔️: assert(core_1[node1] == NULL_NODE);
        // RDKit✔️✔️: assert(core_2[node2] == NULL_NODE);
        debug_assert!(node1 < self.n1);
        debug_assert!(node2 < self.n2);
        debug_assert_eq!(self.core_1[node1], NULL_NODE);
        debug_assert_eq!(self.core_2[node2], NULL_NODE);

        // RDKit✔️✔️: if (boost::out_degree(node1, *g1) > boost::out_degree(node2, *g2)) {
        // RDKit✔️✔️:   return false;
        // RDKit✔️✔️: }
        if self.g1.out_degree(node1) > self.g2.out_degree(node2) {
            return false;
        }

        // RDKit✔️✔️: if (!vc(node1, node2)) { return false; }
        if !atom_fn(node1, node2) {
            return false;
        }

        // RDKit✔️✔️: // Check the out edges of node1
        // RDKit✔️✔️: boost::tie(bNbrs, eNbrs) = boost::out_edges(node1, *g1);
        // RDKit✔️✔️: while (bNbrs != eNbrs) {
        // RDKit✔️✔️:   other1 = getOtherIdx(*g1, *bNbrs, node1);
        // RDKit✔️✔️:   if (core_1[other1] != NULL_NODE) {
        // RDKit✔️✔️:     other2 = core_1[other1];
        // RDKit✔️✔️:     if (!found || !ec(*bNbrs, oEdge)) { return false; }
        // RDKit✔️✔️:   }
        // RDKit✔️✔️:   ++bNbrs;
        // RDKit✔️✔️: }
        for &(other1, edge_idx1) in self.g1.out_edges(node1) {
            if other1 == node1 {
                continue;
            }
            if self.core_1[other1] != NULL_NODE {
                let other2 = self.core_1[other1];
                // Check that (node2, other2) has a matching bond.
                let bond_found = self.find_bond(node2, other2);
                match bond_found {
                    Some(edge_idx2) => {
                        if !bond_fn(edge_idx1, edge_idx2) {
                            return false;
                        }
                    }
                    None => return false,
                }
            }
        }

        true
    }

    /// Find a bond between atom `a` and `b` in the molecule graph (g2).
    fn find_bond(&self, a: NodeId, b: NodeId) -> Option<usize> {
        for &(nbr, bond_idx) in self.g2.out_edges(a) {
            if nbr == b {
                return Some(bond_idx);
            }
        }
        None
    }

    // RDKit source (vf2.hpp), AddPair:
    //   void AddPair(node_id node1, node_id node2) {
    //     assert(node1 < n1); assert(node2 < n2);
    //     assert(core_len < n1); assert(core_len < n2);
    //     ++core_len;
    //     if (!term_1[node1]) { term_1[node1] = core_len; ++t1_len; }
    //     if (!term_2[node2]) { term_2[node2] = core_len; ++t2_len; }
    //     core_1[node1] = node2; core_2[node2] = node1;
    //     // FIX: this is explicitly ignoring directionality
    //     // out_edges of node1...
    //     while (bNbrs != eNbrs) {
    //       unsigned int other = getOtherIdx(*g1, *bNbrs, node1);
    //       if (!term_1[other]) { term_1[other] = core_len; ++t1_len; }
    //       ++bNbrs;
    //     }
    //     // out_edges of node2...
    //     while (bNbrs != eNbrs) {
    //       unsigned int other = getOtherIdx(*g2, *bNbrs, node2);
    //       if (!term_2[other]) { term_2[other] = core_len; ++t2_len; }
    //       ++bNbrs;
    //     }
    //   }

    /// RDKit✔️✔️: AddPair — add (node1, node2) to the mapping.
    ///   Updates terminal sets with depth (core_len) tracking.
    fn add_pair(&mut self, node1: NodeId, node2: NodeId) {
        // RDKit✔️✔️: ++core_len;
        self.core_len += 1;
        let depth = self.core_len;

        // RDKit✔️✔️: if (!term_1[node1]) { term_1[node1] = core_len; ++t1_len; }
        if self.term_1[node1] == 0 {
            self.term_1[node1] = depth;
            self.t1_len += 1;
        }

        // RDKit✔️✔️: if (!term_2[node2]) { term_2[node2] = core_len; ++t2_len; }
        if self.term_2[node2] == 0 {
            self.term_2[node2] = depth;
            self.t2_len += 1;
        }

        // RDKit✔️✔️: core_1[node1] = node2; core_2[node2] = node1;
        self.core_1[node1] = node2;
        self.core_2[node2] = node1;

        // RDKit✔️✔️: // FIX: explicitly ignoring directionality
        // RDKit✔️✔️: boost::tie(bNbrs, eNbrs) = boost::out_edges(node1, *g1);
        // RDKit✔️✔️: while (bNbrs != eNbrs) {
        // RDKit✔️✔️:   unsigned int other = getOtherIdx(*g1, *bNbrs, node1);
        // RDKit✔️✔️:   if (!term_1[other]) { term_1[other] = core_len; ++t1_len; }
        // RDKit✔️✔️:   ++bNbrs;
        // RDKit✔️✔️: }
        for &(other, _) in self.g1.out_edges(node1) {
            if other == node1 {
                continue;
            }
            if self.term_1[other] == 0 {
                self.term_1[other] = depth;
                self.t1_len += 1;
            }
        }

        // RDKit✔️✔️: boost::tie(bNbrs, eNbrs) = boost::out_edges(node2, *g2);
        // RDKit✔️✔️: while (bNbrs != eNbrs) {
        // RDKit✔️✔️:   unsigned int other = getOtherIdx(*g2, *bNbrs, node2);
        // RDKit✔️✔️:   if (!term_2[other]) { term_2[other] = core_len; ++t2_len; }
        // RDKit✔️✔️:   ++bNbrs;
        // RDKit✔️✔️: }
        for &(other, _) in self.g2.out_edges(node2) {
            if other == node2 {
                continue;
            }
            if self.term_2[other] == 0 {
                self.term_2[other] = depth;
                self.t2_len += 1;
            }
        }
    }

    // RDKit source (vf2.hpp), BackTrack:
    //   void BackTrack(node_id node1, node_id node2) {
    //     if (term_1[node1] == core_len) { term_1[node1] = 0; --t1_len; }
    //     boost::tie(bNbrs, eNbrs) = boost::out_edges(node1, *g1);
    //     while (bNbrs != eNbrs) {
    //       unsigned int other = getOtherIdx(*g1, *bNbrs, node1);
    //       if (term_1[other] == core_len) { term_1[other] = 0; --t1_len; }
    //       ++bNbrs;
    //     }
    //     if (term_2[node2] == core_len) { term_2[node2] = 0; --t2_len; }
    //     boost::tie(bNbrs, eNbrs) = boost::out_edges(node2, *g2);
    //     while (bNbrs != eNbrs) {
    //       unsigned int other = getOtherIdx(*g2, *bNbrs, node2);
    //       if (term_2[other] == core_len) { term_2[other] = 0; --t2_len; }
    //       ++bNbrs;
    //     }
    //     core_1[node1] = NULL_NODE;
    //     core_2[node2] = NULL_NODE;
    //     --core_len;
    //   }

    /// RDKit✔️✔️: BackTrack — remove (node1, node2) from mapping.
    ///   Uses depth-based term_X tracking: only clears entries tagged
    ///   with the current core_len. This avoids recomputing terminal
    ///   sets from scratch.
    fn back_track(&mut self, node1: NodeId, node2: NodeId) {
        let depth = self.core_len;

        // RDKit✔️✔️: if (term_1[node1] == core_len) { term_1[node1] = 0; --t1_len; }
        if self.term_1[node1] == depth {
            self.term_1[node1] = 0;
            self.t1_len -= 1;
        }

        // RDKit✔️✔️: boost::tie(bNbrs, eNbrs) = boost::out_edges(node1, *g1);
        // RDKit✔️✔️: while (bNbrs != eNbrs) {
        // RDKit✔️✔️:   unsigned int other = getOtherIdx(*g1, *bNbrs, node1);
        // RDKit✔️✔️:   if (term_1[other] == core_len) { term_1[other] = 0; --t1_len; }
        // RDKit✔️✔️:   ++bNbrs;
        // RDKit✔️✔️: }
        for &(other, _) in self.g1.out_edges(node1) {
            if other == node1 {
                continue;
            }
            if self.term_1[other] == depth {
                self.term_1[other] = 0;
                self.t1_len -= 1;
            }
        }

        // RDKit✔️✔️: if (term_2[node2] == core_len) { term_2[node2] = 0; --t2_len; }
        if self.term_2[node2] == depth {
            self.term_2[node2] = 0;
            self.t2_len -= 1;
        }

        // RDKit✔️✔️: boost::tie(bNbrs, eNbrs) = boost::out_edges(node2, *g2);
        // RDKit✔️✔️: while (bNbrs != eNbrs) {
        // RDKit✔️✔️:   unsigned int other = getOtherIdx(*g2, *bNbrs, node2);
        // RDKit✔️✔️:   if (term_2[other] == core_len) { term_2[other] = 0; --t2_len; }
        // RDKit✔️✔️:   ++bNbrs;
        // RDKit✔️✔️: }
        for &(other, _) in self.g2.out_edges(node2) {
            if other == node2 {
                continue;
            }
            if self.term_2[other] == depth {
                self.term_2[other] = 0;
                self.t2_len -= 1;
            }
        }

        // RDKit✔️✔️: core_1[node1] = NULL_NODE;
        // RDKit✔️✔️: core_2[node2] = NULL_NODE;
        // RDKit✔️✔️: --core_len;
        self.core_1[node1] = NULL_NODE;
        self.core_2[node2] = NULL_NODE;
        self.core_len -= 1;
    }

    // RDKit source (vf2.hpp), GetCoreSet:
    //   void GetCoreSet(node_id c1[], node_id c2[]) {
    //     unsigned int i, j;
    //     for (i = 0, j = 0; i < n1; ++i) {
    //       if (core_1[i] != NULL_NODE) { c1[j] = i; c2[j] = core_1[i]; ++j; }
    //     }
    //   }

    /// RDKit✔️✔️: GetCoreSet — extract current mapping arrays.
    fn get_core_set(&self) -> (Vec<NodeId>, Vec<NodeId>) {
        let mut c1 = Vec::with_capacity(self.core_len);
        let mut c2 = Vec::with_capacity(self.core_len);
        // RDKit✔️✔️: for (i = 0, j = 0; i < n1; ++i) {
        // RDKit✔️✔️:   if (core_1[i] != NULL_NODE) { c1[j] = i; c2[j] = core_1[i]; ++j; }
        // RDKit✔️✔️: }
        for i in 0..self.n1 {
            if self.core_1[i] != NULL_NODE {
                c1.push(i);
                c2.push(self.core_1[i]);
            }
        }
        (c1, c2)
    }
}

// ---------------------------------------------------------------------------
// VF2 recursive matching
// ---------------------------------------------------------------------------
//
// RDKit source (vf2.hpp):
//   bool Match(node_id c1[], node_id c2[]) {
//     if (IsGoal()) { GetCoreSet(c1, c2); if (MatchChecks(c1, c2)) return true; }
//     if (IsDead()) return false;
//     Pair<Graph> pair;
//     while (NextPair(pair)) {
//       if (IsFeasiblePair(pair.n1, pair.n2)) {
//         AddPair(pair.n1, pair.n2);
//         if (Match(c1, c2)) return true;  // recurse
//         BackTrack(pair.n1, pair.n2);
//       }
//     }
//     return false;
//   }

/// RDKit✔️❌: Match — find first match via VF2 recursion.
///
/// Matches RDKit's `Match(c1, c2)` entry point. `match_check` allows
/// final verification (like MolMatchFinalCheckFunctor). If None, all
/// completed matches are accepted.
fn vf2_match(
    state: &mut Vf2SubState,
    atom_fn: &impl Fn(usize, usize) -> bool,
    bond_fn: &impl Fn(usize, usize) -> bool,
    match_check: Option<&impl Fn(&[NodeId], &[NodeId]) -> bool>,
) -> Option<(Vec<NodeId>, Vec<NodeId>)> {
    // RDKit✔️✔️: if (IsGoal()) { GetCoreSet(c1, c2); if (MatchChecks(c1, c2)) return true; }
    if state.is_goal() {
        let (c1, c2) = state.get_core_set();
        if match_check.map_or(true, |check| check(&c1, &c2)) {
            return Some((c1, c2));
        }
    }

    // RDKit✔️✔️: if (IsDead()) return false;
    if state.is_dead() {
        return None;
    }

    // RDKit✔️✔️: Pair<Graph> pair;
    // RDKit✔️✔️: while (NextPair(pair)) {
    // RDKit✔️✔️:   if (IsFeasiblePair(pair.n1, pair.n2)) {
    // RDKit✔️✔️:     AddPair(pair.n1, pair.n2);
    // RDKit✔️✔️:     if (Match(c1, c2)) return true;  // recurse
    // RDKit✔️✔️:     BackTrack(pair.n1, pair.n2);
    // RDKit✔️✔️:   }
    // RDKit✔️✔️: }
    // RDKit✔️✔️: return false;
    let mut pair = Vf2Pair::new();
    while state.next_pair(&mut pair) {
        if state.is_feasible_pair(pair.n1, pair.n2, atom_fn, bond_fn) {
            state.add_pair(pair.n1, pair.n2);
            if let Some(result) = vf2_match(state, atom_fn, bond_fn, match_check) {
                return Some(result);
            }
            state.back_track(pair.n1, pair.n2);
        }
    }
    None
}

// RDKit source (vf2.hpp), MatchAll:
//   template <class DoubleBackInsertionSequence>
//   bool MatchAll(node_id c1[], node_id c2[], DoubleBackInsertionSequence &res,
//                 unsigned int lim = 0) {
//     if (IsGoal()) {
//       GetCoreSet(c1, c2);
//       if (MatchChecks(c1, c2)) {
//         typename DoubleBackInsertionSequence::value_type newSeq;
//         newSeq.reserve(core_len);
//         for (unsigned int i = 0; i < core_len; ++i) {
//           newSeq.emplace_back(c1[i], c2[i]);
//         }
//         res.push_back(newSeq);
//         return lim && res.size() >= lim;
//       }
//     }
//     if (IsDead()) return false;
//     Pair<Graph> pair;
//     while (NextPair(pair)) {
//       if (IsFeasiblePair(pair.n1, pair.n2)) {
//         AddPair(pair.n1, pair.n2);
//         if (MatchAll(c1, c2, res, lim)) return true;  // recurse
//         BackTrack(pair.n1, pair.n2);
//       }
//     }
//     return false;
//   }

/// RDKit✔️❌: MatchAll — find all matches up to `max_matches`.
///
/// Collects matches into `results` as (c1, c2) pairs.
/// Returns true when the limit has been reached, signaling the caller
/// to stop.
fn vf2_match_all(
    state: &mut Vf2SubState,
    atom_fn: &impl Fn(usize, usize) -> bool,
    bond_fn: &impl Fn(usize, usize) -> bool,
    match_check: Option<&impl Fn(&[NodeId], &[NodeId]) -> bool>,
    results: &mut Vec<(Vec<NodeId>, Vec<NodeId>)>,
    max_matches: usize,
) -> bool {
    // RDKit✔️✔️: if (IsGoal()) { ... }
    if state.is_goal() {
        let (c1, c2) = state.get_core_set();
        if match_check.map_or(true, |check| check(&c1, &c2)) {
            // RDKit✔️✔️: newSeq.reserve(core_len);
            // RDKit✔️✔️: for (unsigned int i = 0; i < core_len; ++i) {
            // RDKit✔️✔️:   newSeq.emplace_back(c1[i], c2[i]);
            // RDKit✔️✔️: }
            // RDKit✔️✔️: res.push_back(newSeq);
            results.push((c1, c2));
            // RDKit✔️✔️: return lim && res.size() >= lim;
            return max_matches > 0 && results.len() >= max_matches;
        }
    }

    // RDKit✔️✔️: if (IsDead()) return false;
    if state.is_dead() {
        return false;
    }

    // RDKit✔️✔️: Pair<Graph> pair;
    // RDKit✔️✔️: while (NextPair(pair)) { ... }
    let mut pair = Vf2Pair::new();
    while state.next_pair(&mut pair) {
        if state.is_feasible_pair(pair.n1, pair.n2, atom_fn, bond_fn) {
            state.add_pair(pair.n1, pair.n2);
            // RDKit✔️✔️: if (MatchAll(c1, c2, res, lim)) return true;
            if vf2_match_all(state, atom_fn, bond_fn, match_check, results, max_matches) {
                return true;
            }
            state.back_track(pair.n1, pair.n2);
        }
    }
    false
}

// ---------------------------------------------------------------------------
// Final match check (simplified MolMatchFinalCheckFunctor)
// ---------------------------------------------------------------------------
//
// RDKit source (SubstructMatch.cpp):
//   bool MolMatchFinalCheckFunctor::operator()(const std::uint32_t q_c[],
//                                              const std::uint32_t m_c[]) {
//     if (d_params.extraFinalCheck || d_params.useGenericMatchers) { ... }
//     HashedStorageType match;
//     if (d_params.uniquify) {
//       match.resize(d_mol.getNumAtoms());
//       std::fill(match.begin(), match.end(), 0);
//       for (unsigned int i = 0; i < d_query.getNumAtoms(); ++i) {
//         match[m_c[i]] = 1;
//       }
//       if (matchesSeen.find(match) != matchesSeen.end()) { return false; }
//     }
//     if (!d_params.useChirality) {
//       if (d_params.uniquify) { matchesSeen.insert(match); }
//       return true;
//     }
//     // ... chirality checks ...
//   }

/// RDKit❗✔️: Final match check with uniquification support.
///
/// Simplified version of RDKit's MolMatchFinalCheckFunctor.
/// Chirality checking and enhanced stereo are deferred.
fn match_mask(atom_mapping: &[usize], mol_num_atoms: usize) -> Vec<bool> {
    let mut mask = vec![false; mol_num_atoms];
    for &ma in atom_mapping {
        if ma < mol_num_atoms {
            mask[ma] = true;
        }
    }
    mask
}

fn rdkit_order_and_uniquify_matches(
    mut results: Vec<SubstructMatchResult>,
    mol_num_atoms: usize,
    params: &SubstructMatchParams,
) -> Vec<SubstructMatchResult> {
    // RDKit source (SubstructMatch.cpp):
    //   bool tryToInsert(std::set<MatchVectType> &matches, const MatchVectType &match,
    //   const SubstructMatchParameters &params) {
    //     if (matches.size() == params.maxMatches) {
    //       return false;
    //     }
    //     if (!params.uniquify) {
    //       matches.insert(match);
    //     } else {
    //       insertIfNeeded(matches, match);
    //     }
    //     return true;
    //   }
    // RDKit source (SubstructMatch.cpp):
    //   bool insertIfNeeded(std::set<MatchVectType> &matches, const MatchVectType &m) {
    //     bool shouldInsert = true;
    //     ...
    //     if (matchAsSet == existingMatchAsSet) {
    //       if (m < *it) {
    //         matches.erase(it);
    //       } else {
    //         shouldInsert = false;
    //       }
    //       break;
    //     }
    //     ...
    //   }
    // RDKit source (SubstructMatch.cpp):
    //   return std::vector<MatchVectType>(matches.begin(), matches.end());
    if params.uniquify {
        let mut uniquified: Vec<(Vec<bool>, SubstructMatchResult)> = Vec::new();
        for result in results.drain(..) {
            let mask = match_mask(&result.atom_mapping, mol_num_atoms);
            if let Some((_, existing)) = uniquified
                .iter_mut()
                .find(|(existing_mask, _)| *existing_mask == mask)
            {
                if result.atom_mapping < existing.atom_mapping {
                    *existing = result;
                }
            } else {
                uniquified.push((mask, result));
            }
        }
        results = uniquified.into_iter().map(|(_, result)| result).collect();
        results.sort_by(|lhs, rhs| lhs.atom_mapping.cmp(&rhs.atom_mapping));
    }
    if results.len() > params.max_matches {
        results.truncate(params.max_matches);
    }
    results
}

// ---------------------------------------------------------------------------
// Bond mapping builder
// ---------------------------------------------------------------------------

/// Build the bond mapping for a match result.
///
/// For each query bond (by index), find the corresponding molecular bond
/// that connects the matched query endpoints.
#[allow(dead_code)]
fn build_bond_mapping(
    query_atom_to_mol: &[Option<usize>],
    query: &Vf2Graph,
    mol: &Vf2Graph,
) -> Vec<usize> {
    let mut bond_mapping = Vec::with_capacity(query.n_bonds);
    for bond_idx in 0..query.n_bonds {
        // Find the query atoms connected by this bond.
        let mut q_begin = NULL_NODE;
        let mut q_end = NULL_NODE;
        for qa in 0..query.n_atoms {
            for &(nbr, eidx) in &query.adjacency[qa] {
                if eidx == bond_idx {
                    q_begin = qa;
                    q_end = nbr;
                    break;
                }
            }
            if q_begin != NULL_NODE {
                break;
            }
        }

        if q_begin != NULL_NODE {
            let m_begin = query_atom_to_mol[q_begin];
            let m_end = query_atom_to_mol[q_end];
            if let (Some(mb), Some(me)) = (m_begin, m_end) {
                // Find bond between mb and me in mol.
                let mut mol_bond_idx = NULL_NODE;
                for &(nbr, eidx) in &mol.adjacency[mb] {
                    if nbr == me {
                        mol_bond_idx = eidx;
                        break;
                    }
                }
                bond_mapping.push(mol_bond_idx);
            } else {
                bond_mapping.push(NULL_NODE);
            }
        } else {
            bond_mapping.push(NULL_NODE);
        }
    }
    bond_mapping
}

// ---------------------------------------------------------------------------
// Public API
// ---------------------------------------------------------------------------

/// RDKit❌❌: `SubstructMatch` entry point.
///
/// Wires up atom/bond functors and runs VF2. This is a simplified port of
/// RDKit's `SubstructMatch()` from SubstructMatch.cpp (lines 481-525).
///
/// Returns all matches (up to params.max_matches).
fn collect_recursive_smarts_from_atom_query(
    query: &crate::QueryNode<AtomQueryPredicate>,
    recursive_smarts: &mut Vec<String>,
) {
    match query {
        crate::QueryNode::Predicate(AtomQueryPredicate::RecursiveSmarts(smarts)) => {
            recursive_smarts.push(smarts.clone());
        }
        crate::QueryNode::Predicate(_) => {}
        crate::QueryNode::And(children) | crate::QueryNode::Or(children) => {
            for child in children {
                collect_recursive_smarts_from_atom_query(child, recursive_smarts);
            }
        }
        crate::QueryNode::Not(child) => {
            collect_recursive_smarts_from_atom_query(child, recursive_smarts);
        }
    }
}

fn populate_recursive_query_match_cache(
    mol: &Molecule,
    query: &Molecule,
    recursive_cache: &mut RecursiveQueryMatchCache,
) {
    let mut recursive_smarts = Vec::new();
    for atom in query.atoms() {
        if let Some(query_node) = atom.query() {
            collect_recursive_smarts_from_atom_query(query_node, &mut recursive_smarts);
        }
    }
    recursive_smarts.sort();
    recursive_smarts.dedup();

    for recursive_smarts in recursive_smarts {
        if recursive_cache.contains_key(&recursive_smarts) {
            continue;
        }
        let Some(inner) = recursive_smarts
            .strip_prefix("$(")
            .and_then(|value| value.strip_suffix(')'))
        else {
            recursive_cache.insert(recursive_smarts, vec![false; mol.num_atoms()]);
            continue;
        };
        let Ok(inner_query) = build_crystalff_query_molecule(inner) else {
            recursive_cache.insert(recursive_smarts, vec![false; mol.num_atoms()]);
            continue;
        };
        populate_recursive_query_match_cache(mol, &inner_query, recursive_cache);
        let matches = substruct_match_impl_with_recursive_cache(
            mol,
            &inner_query,
            &SubstructMatchParams {
                max_matches: 1000,
                uniquify: false,
            },
            Some(recursive_cache),
        );
        let mut match_starts = vec![false; mol.num_atoms()];
        for matched in matches {
            if let Some(&root_atom_idx) = matched.atom_mapping.first()
                && root_atom_idx != NULL_NODE
                && root_atom_idx < match_starts.len()
            {
                match_starts[root_atom_idx] = true;
            }
        }
        recursive_cache.insert(recursive_smarts, match_starts);
    }
}

fn substruct_match_impl_with_recursive_cache(
    mol: &Molecule,
    query: &Molecule,
    params: &SubstructMatchParams,
    recursive_cache: Option<&RecursiveQueryMatchCache>,
) -> Vec<SubstructMatchResult> {
    let m_num_atoms = mol.num_atoms();
    let q_num_atoms = query.num_atoms();
    let trace_row64_substruct =
        env::var("RDKIT_ROW64_SUBSTRUCT_TRACE").ok().as_deref() == Some("1") && m_num_atoms == 106;
    let trace_row61_torsion349 = row61_torsion349_substruct_trace_enabled(mol, query, params);
    let query_ctx = build_query_match_context(mol);

    // RDKit source (SubstructMatch.cpp):
    //   if (!mNumAtoms || !qNumAtoms || qNumAtoms > mNumAtoms) {
    //     return matches;
    //   }
    if m_num_atoms == 0 || q_num_atoms == 0 || q_num_atoms > m_num_atoms {
        return Vec::new();
    }

    // Build VF2 graphs.
    let graph_start = trace_row64_substruct.then(Instant::now);
    let q_graph = build_vf2_graph(query);
    let m_graph = build_vf2_graph(mol);
    let graph_elapsed = graph_start.map(|start| start.elapsed().as_secs_f64());

    // Build atom matching closure.
    // RDKit source:
    //   detail::AtomLabelFunctor atomLabeler(query, mol, params);
    //   detail::BondLabelFunctor bondLabeler(query, mol, params);
    //   MolMatchFinalCheckFunctor matchChecker(query, mol, params);
    let atom_fn = |qi: usize, mj: usize| -> bool {
        let qa = &query.atoms()[qi];
        let ma = &mol.atoms()[mj];
        atom_matches_with_recursive_cache(qa, query, ma, mol, recursive_cache, &query_ctx)
    };

    let bond_fn = |qei: usize, mei: usize| -> bool {
        let qb = &query.bonds()[qei];
        let mb = &mol.bonds()[mei];
        if let Some(query_node) = qb.query() {
            evaluate_bond_query(query_node, mb, mol, &query_ctx)
        } else {
            bond_matches(qb, query, mb, mol)
        }
    };

    // RDKit source:
    //   bool found = boost::vf2_all(query.getTopology(), mol.getTopology(),
    //                               atomLabeler, bondLabeler, matchChecker,
    //                               pms, params.maxMatches);
    let mut state = Vf2SubState::new(&q_graph, &m_graph, /*sort_nodes=*/ false);
    if trace_row61_torsion349 {
        println!(
            "row61_substruct_begin q_num_atoms={} q_num_bonds={} m_num_atoms={} order={:?}",
            q_num_atoms,
            query.num_bonds(),
            m_num_atoms,
            state.debug_order()
        );
    }
    let mut raw_matches: Vec<(Vec<NodeId>, Vec<NodeId>)> = Vec::new();
    let check_fn = |_c1: &[NodeId], _c2: &[NodeId]| -> bool { true };

    let vf2_start = trace_row64_substruct.then(Instant::now);
    vf2_match_all(
        &mut state,
        &atom_fn,
        &bond_fn,
        Some(&check_fn),
        &mut raw_matches,
        params.max_matches,
    );
    if let Some(vf2_start) = vf2_start {
        println!(
            "row64_substruct_core q_atoms={} q_bonds={} graph_build={:.6} vf2={:.6} raw_matches={}",
            q_num_atoms,
            query.num_bonds(),
            graph_elapsed.unwrap_or(0.0),
            vf2_start.elapsed().as_secs_f64(),
            raw_matches.len()
        );
    }

    // RDKit source (SubstructMatch.cpp):
    //   if (found) {
    //     const unsigned int nQueryAtoms = query.getNumAtoms();
    //     matches.reserve(pms.size());
    //     MatchVectType matchVect(nQueryAtoms);
    //     for (const auto &pairs : pms) {
    //       for (const auto &pair : pairs) {
    //         matchVect[pair.first] = pair;
    //       }
    //       matches.push_back(matchVect);
    //     }
    //   }
    let mut results: Vec<SubstructMatchResult> = Vec::new();

    for (c1, c2) in &raw_matches {
        if trace_row61_torsion349 {
            println!("row61_substruct_raw_core c1={c1:?} c2={c2:?}");
        }
        // Build atom_mapping: query_atom_index -> mol_atom_index.
        // RDKit uses MatchVectType (vector<pair<int,int>>) where
        // pair.second is the mol atom index and pair.first is query atom index.
        let mut atom_to_mol: Vec<Option<usize>> = vec![None; q_num_atoms];
        for (&qa, &ma) in c1.iter().zip(c2.iter()) {
            if qa < q_num_atoms {
                atom_to_mol[qa] = Some(ma);
            }
        }

        // Build bond mapping by looking up bonds between matched atoms.
        let mut bond_mapping = Vec::with_capacity(query.num_bonds());
        for qbond in query.bonds() {
            let q_begin = qbond.begin().index();
            let q_end = qbond.end().index();
            let m_begin = atom_to_mol[q_begin];
            let m_end = atom_to_mol[q_end];
            match (m_begin, m_end) {
                (Some(mb), Some(me)) => {
                    // Find bond between mb and me in mol.
                    let found = m_graph.adjacency[mb]
                        .iter()
                        .find(|&&(nbr, _)| nbr == me)
                        .map(|&(_, eidx)| eidx);
                    bond_mapping.push(found.unwrap_or(NULL_NODE));
                }
                _ => {
                    bond_mapping.push(NULL_NODE);
                }
            }
        }

        results.push(SubstructMatchResult {
            atom_mapping: atom_to_mol
                .into_iter()
                .map(|x| x.unwrap_or(NULL_NODE))
                .collect(),
            bond_mapping,
        });
    }

    rdkit_order_and_uniquify_matches(results, m_num_atoms, params)
}

fn atom_matches_with_recursive_cache(
    query_atom: &Atom,
    query_mol: &Molecule,
    mol_atom: &Atom,
    mol: &Molecule,
    recursive_cache: Option<&RecursiveQueryMatchCache>,
    query_ctx: &QueryMatchContext,
) -> bool {
    if let Some(query_node) = query_atom.query() {
        return match query_node {
            crate::QueryNode::Predicate(pred) => atom_query_predicate_matches_for_substruct(
                mol_atom,
                pred,
                mol,
                recursive_cache,
                query_ctx,
            ),
            crate::QueryNode::And(children) => children.iter().all(|child| match child {
                crate::QueryNode::Predicate(pred) => atom_query_predicate_matches_for_substruct(
                    mol_atom,
                    pred,
                    mol,
                    recursive_cache,
                    query_ctx,
                ),
                _ => evaluate_atom_query(child, mol_atom, mol, recursive_cache, query_ctx),
            }),
            crate::QueryNode::Or(children) => children.iter().any(|child| match child {
                crate::QueryNode::Predicate(pred) => atom_query_predicate_matches_for_substruct(
                    mol_atom,
                    pred,
                    mol,
                    recursive_cache,
                    query_ctx,
                ),
                _ => evaluate_atom_query(child, mol_atom, mol, recursive_cache, query_ctx),
            }),
            crate::QueryNode::Not(child) => {
                !evaluate_atom_query(child, mol_atom, mol, recursive_cache, query_ctx)
            }
        };
    }
    atom_matches(query_atom, query_mol, mol_atom, mol)
}

fn substruct_match_impl(
    mol: &Molecule,
    query: &Molecule,
    params: &SubstructMatchParams,
) -> Vec<SubstructMatchResult> {
    let trace_row64_substruct = env::var("RDKIT_ROW64_SUBSTRUCT_TRACE").ok().as_deref()
        == Some("1")
        && mol.num_atoms() == 106;
    let recursive_cache_start = trace_row64_substruct.then(Instant::now);
    let mut recursive_cache = RecursiveQueryMatchCache::new();
    populate_recursive_query_match_cache(mol, query, &mut recursive_cache);
    let recursive_cache_elapsed = recursive_cache_start.map(|start| start.elapsed().as_secs_f64());
    let match_start = trace_row64_substruct.then(Instant::now);
    let results =
        substruct_match_impl_with_recursive_cache(mol, query, params, Some(&recursive_cache));
    if let Some(match_start) = match_start {
        println!(
            "row64_substruct_timing q_atoms={} q_bonds={} recursive_cache={:.6} match_core={:.6} matches={}",
            query.num_atoms(),
            query.num_bonds(),
            recursive_cache_elapsed.unwrap_or(0.0),
            match_start.elapsed().as_secs_f64(),
            results.len()
        );
    }
    results
}

/// Check if a molecule contains a substructure match for the given query.
///
/// This is the public API for `has_substruct_match`.
/// RDKit✔️❌: VF2-based substructure matching ported from vf2.hpp + SubstructMatch.cpp.
pub fn has_substruct_match(mol: &Molecule, query: &Molecule) -> bool {
    let params = SubstructMatchParams::default();
    let mut params = params;
    params.max_matches = 1;
    !substruct_match_impl(mol, query, &params).is_empty()
}

/// Get the first substructure match, if any.
///
/// This is the public API for `get_substruct_match`.
/// RDKit✔️❌: VF2-based substructure matching ported from vf2.hpp + SubstructMatch.cpp.
pub fn get_substruct_match(mol: &Molecule, query: &Molecule) -> Option<SubstructMatchResult> {
    let params = SubstructMatchParams::default();
    let mut params = params;
    params.max_matches = 1;
    substruct_match_impl(mol, query, &params).into_iter().next()
}

/// Get all substructure matches with default parameters.
///
/// This is the public API for `get_substruct_matches`.
/// RDKit✔️❌: VF2-based substructure matching ported from vf2.hpp + SubstructMatch.cpp.
pub fn get_substruct_matches(mol: &Molecule, query: &Molecule) -> Vec<SubstructMatchResult> {
    let params = SubstructMatchParams::default();
    substruct_match_impl(mol, query, &params)
}

/// Get all substructure matches with custom parameters.
///
/// This is the public API for `get_substruct_matches_with_params`.
/// RDKit✔️❌: VF2-based substructure matching ported from vf2.hpp + SubstructMatch.cpp.
pub fn get_substruct_matches_with_params(
    mol: &Molecule,
    query: &Molecule,
    params: &SubstructMatchParams,
) -> Vec<SubstructMatchResult> {
    substruct_match_impl(mol, query, params)
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

#[cfg(test)]
mod tests {
    use super::*;
    use crate::MoleculeBuilder;

    fn make_mol_c() -> Molecule {
        // Methane: C
        let mut builder = MoleculeBuilder::new();
        builder.add_atom(crate::AtomSpec::new(crate::Element::C));
        builder.build().expect("build methane")
    }

    fn make_mol_cc() -> Molecule {
        // Ethane: CC
        let mut builder = MoleculeBuilder::new();
        let c0 = builder.add_atom(crate::AtomSpec::new(crate::Element::C));
        let c1 = builder.add_atom(crate::AtomSpec::new(crate::Element::C));
        builder
            .add_bond(crate::BondSpec::new(c0, c1, BondOrder::Single))
            .expect("add bond");
        builder.build().expect("build ethane")
    }

    fn make_mol_cco() -> Molecule {
        // Ethanol: CCO
        let mut builder = MoleculeBuilder::new();
        let c0 = builder.add_atom(crate::AtomSpec::new(crate::Element::C));
        let c1 = builder.add_atom(crate::AtomSpec::new(crate::Element::C));
        let o = builder.add_atom(crate::AtomSpec::new(crate::Element::O));
        builder
            .add_bond(crate::BondSpec::new(c0, c1, BondOrder::Single))
            .expect("add C-C bond");
        builder
            .add_bond(crate::BondSpec::new(c1, o, BondOrder::Single))
            .expect("add C-O bond");
        builder.build().expect("build ethanol")
    }

    fn make_mol_coc() -> Molecule {
        // Dimethyl ether: COC
        let mut builder = MoleculeBuilder::new();
        let c0 = builder.add_atom(crate::AtomSpec::new(crate::Element::C));
        let o = builder.add_atom(crate::AtomSpec::new(crate::Element::O));
        let c1 = builder.add_atom(crate::AtomSpec::new(crate::Element::C));
        builder
            .add_bond(crate::BondSpec::new(c0, o, BondOrder::Single))
            .expect("add C-O bond 1");
        builder
            .add_bond(crate::BondSpec::new(o, c1, BondOrder::Single))
            .expect("add O-C bond 2");
        builder.build().expect("build dimethyl ether")
    }

    #[test]
    fn test_has_substruct_match_self() {
        let c = make_mol_c();
        assert!(
            has_substruct_match(&c, &c),
            "a molecule should match itself"
        );
    }

    #[test]
    fn test_has_substruct_match_cc_in_cco() {
        let cc = make_mol_cc();
        let cco = make_mol_cco();
        assert!(
            has_substruct_match(&cco, &cc),
            "CCO should contain CC as substructure"
        );
    }

    #[test]
    fn test_has_substruct_match_no_match() {
        let c = make_mol_c();
        let cco = make_mol_cco();
        assert!(
            !has_substruct_match(&c, &cco),
            "a single carbon should not contain CCO"
        );
    }

    #[test]
    fn test_get_substruct_match_self() {
        let cco = make_mol_cco();
        let result = get_substruct_match(&cco, &cco);
        assert!(result.is_some(), "self-match should return Some");
        let result = result.unwrap();
        assert_eq!(result.atom_mapping.len(), 3);
        // Identity mapping: 0->0, 1->1, 2->2
        for (qa, ma) in result.atom_mapping.iter().enumerate() {
            assert_eq!(*ma, qa, "self-match should have identity mapping");
        }
    }

    #[test]
    fn test_get_substruct_match_cc_in_cco() {
        let cc = make_mol_cc();
        let cco = make_mol_cco();
        let result = get_substruct_match(&cco, &cc);
        assert!(result.is_some(), "CC should match in CCO");
    }

    #[test]
    fn test_get_substruct_match_no_match() {
        let c = make_mol_c();
        let cco = make_mol_cco();
        let result = get_substruct_match(&c, &cco);
        assert!(
            result.is_none(),
            "C should not match CCO (query larger than mol)"
        );
    }

    #[test]
    fn test_get_substruct_matches_cco_in_cco() {
        let cco = make_mol_cco();
        let matches = get_substruct_matches(&cco, &cco);
        assert!(!matches.is_empty(), "should find at least self-match");
    }

    #[test]
    fn test_substruct_coc_matches_cco() {
        // COC (dimethyl ether) should not match CCO (ethanol) — different topology.
        let coc = make_mol_coc();
        let cco = make_mol_cco();
        assert!(
            !has_substruct_match(&cco, &coc),
            "CCO should not match COC topology"
        );
        // But CO should match CCO (CO is a substructure of CCO).
        let mut builder = MoleculeBuilder::new();
        let c = builder.add_atom(crate::AtomSpec::new(crate::Element::C));
        let o = builder.add_atom(crate::AtomSpec::new(crate::Element::O));
        builder
            .add_bond(crate::BondSpec::new(c, o, BondOrder::Single))
            .expect("add CO bond");
        let co = builder.build().expect("build CO");
        assert!(has_substruct_match(&cco, &co), "CCO should match CO");
    }

    #[test]
    fn test_has_substruct_match_empty_mol() {
        let empty = Molecule::new();
        let c = make_mol_c();
        assert!(
            !has_substruct_match(&empty, &c),
            "empty molecule should not match anything"
        );
        assert!(
            !has_substruct_match(&c, &empty),
            "molecule should not match empty query"
        );
    }

    #[test]
    fn test_substruct_match_params_max_matches() {
        let c = make_mol_c();
        let params = SubstructMatchParams {
            max_matches: 1,
            uniquify: true,
        };
        let matches = get_substruct_matches_with_params(&c, &c, &params);
        assert_eq!(matches.len(), 1, "max_matches=1 should return one match");
    }

    #[test]
    fn test_atom_matches_basic() {
        let mut builder = MoleculeBuilder::new();
        let c0 = builder.add_atom(crate::AtomSpec::new(crate::Element::C));
        let c1 = builder.add_atom(crate::AtomSpec::new(crate::Element::C));
        builder
            .add_bond(crate::BondSpec::new(c0, c1, BondOrder::Single))
            .expect("add bond");
        let mol = builder.build().expect("build");
        let q_atom = &mol.atoms()[0];
        let m_atom = &mol.atoms()[1];
        assert!(
            atom_matches(q_atom, &mol, m_atom, &mol),
            "two carbons should match"
        );
    }

    #[test]
    fn test_atom_matches_different_elements() {
        let mut builder = MoleculeBuilder::new();
        let c = builder.add_atom(crate::AtomSpec::new(crate::Element::C));
        let o = builder.add_atom(crate::AtomSpec::new(crate::Element::O));
        builder
            .add_bond(crate::BondSpec::new(c, o, BondOrder::Single))
            .expect("add bond");
        let mol = builder.build().expect("build");
        // Query atom = C, target atom = O — should not match (C != O, and
        // atomic number check should reject since 6 != 8).
        // But our atom_matches uses query atomic number: query=6, mol=8.
        // If query atomic number != 0, it must match. So C does not match O.
        let c_atom = &mol.atoms()[0]; // C
        let o_atom = &mol.atoms()[1]; // O
        assert!(
            !atom_matches(c_atom, &mol, o_atom, &mol),
            "C should not match O via basic atomic number check"
        );
    }

    #[test]
    fn test_bond_matches_single() {
        let mut builder = MoleculeBuilder::new();
        let c0 = builder.add_atom(crate::AtomSpec::new(crate::Element::C));
        let c1 = builder.add_atom(crate::AtomSpec::new(crate::Element::C));
        builder
            .add_bond(crate::BondSpec::new(c0, c1, BondOrder::Single))
            .expect("add bond");
        let mol = builder.build().expect("build");
        assert!(bond_matches(&mol.bonds()[0], &mol, &mol.bonds()[0], &mol));
    }

    #[test]
    fn test_vf2_graph_building() {
        let cc = make_mol_cc();
        let g = build_vf2_graph(&cc);
        assert_eq!(g.n_atoms, 2);
        assert_eq!(g.n_bonds, 1);
        assert_eq!(g.out_degree(0), 1);
        assert_eq!(g.out_degree(1), 1);
        assert_eq!(g.out_edges(0)[0].0, 1);
        assert_eq!(g.out_edges(1)[0].0, 0);
    }

    #[test]
    fn test_sort_nodes_by_frequency_small() {
        let cc = make_mol_cc();
        let g = build_vf2_graph(&cc);
        let order = sort_nodes_by_frequency(&g);
        // Both nodes have degree 1, so order depends on sort stability.
        assert_eq!(order.len(), 2);
        // Both should be present.
        assert!(order.contains(&0));
        assert!(order.contains(&1));
    }

    #[test]
    fn test_vf2_state_initial() {
        let cc = make_mol_cc();
        let g = build_vf2_graph(&cc);
        let state = Vf2SubState::new(&g, &g, true);
        assert!(!state.is_goal());
        assert!(!state.is_dead());
        assert_eq!(state.core_len, 0);
    }

    #[test]
    fn test_vf2_next_pair_initial() {
        let cc = make_mol_cc();
        let g = build_vf2_graph(&cc);
        let state = Vf2SubState::new(&g, &g, true);
        let mut pair = Vf2Pair::new();
        let has_next = state.next_pair(&mut pair);
        assert!(has_next, "should find a next pair");
    }
}