From 27ac684278a902c76b082a413c93c7d579f56ec4 Mon Sep 17 00:00:00 2001 From: Kevin Boyd Date: Fri, 24 Jul 2026 10:22:22 -0400 Subject: [PATCH 1/8] Add fMCS fast incremental matching --- src/mcs/fmcs_cuda/fmcs_match.cuh | 455 ++++++++++++++++++++++ tests/CMakeLists.txt | 6 + tests/test_fmcs_match.cu | 645 +++++++++++++++++++++++++++++++ 3 files changed, 1106 insertions(+) create mode 100644 src/mcs/fmcs_cuda/fmcs_match.cuh create mode 100644 tests/test_fmcs_match.cu diff --git a/src/mcs/fmcs_cuda/fmcs_match.cuh b/src/mcs/fmcs_cuda/fmcs_match.cuh new file mode 100644 index 00000000..afdc049b --- /dev/null +++ b/src/mcs/fmcs_cuda/fmcs_match.cuh @@ -0,0 +1,455 @@ +// SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +// SPDX-License-Identifier: Apache-2.0 +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#ifndef FMCS_CUDA_FMCS_MATCH_CUH +#define FMCS_CUDA_FMCS_MATCH_CUH + +#include + +#include "fmcs_cuda/fmcs_match_tables.cuh" +#include "fmcs_cuda/fmcs_seed.cuh" + +namespace mcs { +namespace fmcs { + +/// Bond endpoints are stored across the kernel as a single uint32 with +/// the u-endpoint atom index in the high 16 bits and the v-endpoint +/// atom index in the low 16 bits. Tiers cap maxAtoms at 128 so 16 bits +/// per index is plenty. +constexpr int kBondEndpointShift = 16; +constexpr std::uint32_t kBondEndpointMask = 0xFFFFu; + +/// Resolved target endpoints for a successful single-(query bond, target +/// bond, orientation) compatibility check. Populated by +/// @ref matchSingleBondWithinThread on success only; contents are +/// unspecified on failure. @c targetAtomU is the target atom that the +/// query bond's u endpoint was mapped to (likewise V). +struct SingleBondMatch { + uint8_t targetAtomU; + uint8_t targetAtomV; +}; + +template struct FmcsSubstructureScratch { + // Scratch for the RDKit checkIfMatchAndAppend fallback. This is deliberately + // caller-owned shared memory, not function-local state: tier-128 scratch is + // too large to risk compiler-created stack/local memory in the matcher. + std::uint8_t seedAtomList[maxAtoms]; + std::uint8_t seedAtoms[maxAtoms]; + std::uint8_t seedDegree[maxAtoms]; + std::uint8_t targetDegree[maxTargetAtoms]; + std::uint8_t orderedQueryAtom[maxAtoms]; + std::uint8_t queryOrderPos[maxAtoms]; + std::uint8_t targetAtomForQuery[maxAtoms]; + int currentCount; + int nextCount; + int found; + int overflowed; +}; + +template +__device__ __forceinline__ bool seedContainsBondWithinThread(const Seed& seed, + const int queryBondIdx) { + using SeedT = Seed; + using BondWord = typename SeedT::bond_word_type; + constexpr int kBondBitsPerWord = SeedT::kBondBitsPerWord; + if (queryBondIdx < 0 || queryBondIdx >= maxBonds) + return false; + const BondWord word = seed.bonds[queryBondIdx / kBondBitsPerWord]; + return ((word >> (queryBondIdx % kBondBitsPerWord)) & 1) != 0; +} + +template __host__ __device__ constexpr bool topologyHasAdjacencyBondIndices() { + if constexpr (requires { Topology::kHasAdjacencyBondIndices; }) { + return Topology::kHasAdjacencyBondIndices; + } else { + return false; + } +} + +/// Within-thread: per-lane single-(query bond, target bond, orientation) +/// compatibility check used by Phase 1 initial-seed enumeration. Writes +/// resolved target atom indices for the two endpoints of @p queryBondIdx +/// into @p outMatch and returns true on success; on false the caller +/// should not read @p outMatch. +/// +/// @p reversed selects the orientation: when false, the query bond's u +/// endpoint maps to the target bond's u endpoint; when true, to the +/// target bond's v endpoint. +/// +/// @p queryTopology and @p targetTopology must expose a +/// @c bondEndpoints array of packed (u<<16 | v) entries. +template +__device__ __forceinline__ bool matchSingleBondWithinThread(const int queryBondIdx, + const int targetBondIdx, + const bool reversed, + const QueryTopology& queryTopology, + const TargetTopology& targetTopology, + const PairMatchTablesDevice& tables, + SingleBondMatch& outMatch) { + // Cheap bond-table check first; if the bond labels are incompatible + // we never need to touch the atom table or decode endpoints. + if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) + return false; + + // Decode the packed (u<<16 | v) endpoint encoding for both bonds. + const std::uint32_t queryEndpoints = queryTopology.bondEndpoints[queryBondIdx]; + const int queryEndpointU = static_cast(queryEndpoints >> kBondEndpointShift); + const int queryEndpointV = static_cast(queryEndpoints & kBondEndpointMask); + + const std::uint32_t targetEndpoints = targetTopology.bondEndpoints[targetBondIdx]; + const int targetEndpointU = static_cast(targetEndpoints >> kBondEndpointShift); + const int targetEndpointV = static_cast(targetEndpoints & kBondEndpointMask); + + // Pick which target endpoint the query's u maps to; v gets the other. + // reversed=false: queryU -> targetU, queryV -> targetV. + // reversed=true: queryU -> targetV, queryV -> targetU. + const int targetForQueryU = reversed ? targetEndpointV : targetEndpointU; + const int targetForQueryV = reversed ? targetEndpointU : targetEndpointV; + + // Both atom-pairings must be label-compatible; either failure rejects + // this orientation. + if (!tables.atoms.testBit(queryEndpointU, targetForQueryU)) + return false; + if (!tables.atoms.testBit(queryEndpointV, targetForQueryV)) + return false; + + outMatch.targetAtomU = static_cast(targetForQueryU); + outMatch.targetAtomV = static_cast(targetForQueryV); + return true; +} + +/// Cooperative: extend @p match by every query bond in @p seed.bonds +/// whose @c match.targetBondIdx[q] is still @ref kUnmappedTargetIdx +/// (i.e., unmapped by the parent's recorded embedding). For each such +/// bond: +/// - Both endpoints already mapped -> ring-closing case. The lanes +/// of @p group scan target bonds in parallel for one whose endpoint +/// pair exactly matches the mapped (queryU, queryV) target atoms +/// and is unvisited, with the bond-match-table bit set. First +/// compatible target bond commits. +/// - Exactly one endpoint mapped -> atom-adding case. Lane-parallel +/// scan for a target bond incident to the mapped target atom whose +/// other end is unvisited, atom-table-compatible with the unmapped +/// query atom, and bond-table-compatible. First compatible +/// candidate commits both the new bond mapping and the new atom +/// mapping, and marks both visited. +/// - Both endpoints unmapped -> defensive fail (shouldn't occur on +/// well-formed seeds, where Phase 1 maps both initial atoms before +/// pushing). +/// Any bond that fails to extend causes the function to return false; +/// @p match is left in an unspecified state and the caller should +/// discard the seed. +template +__device__ __forceinline__ bool matchIncrementalFastCooperative(const GroupT& group, + const Seed& seed, + const QueryTopology& queryTopology, + const TargetTopology& targetTopology, + const PairMatchTablesDevice& tables, + MatchResult& match) { + using SeedT = Seed; + using MatchT = MatchResult; + using BondWord = typename SeedT::bond_word_type; + + constexpr int kBondBitsPerWord = SeedT::kBondBitsPerWord; + constexpr int kBondWords = SeedT::kBondWords; + constexpr int kTargetAtomBitsPerWord = MatchT::kTargetAtomBitsPerWord; + constexpr int kTargetBondBitsPerWord = MatchT::kTargetBondBitsPerWord; + using TargetAtomWord = typename MatchT::target_atom_word; + using TargetBondWord = typename MatchT::target_bond_word; + + const int laneRank = static_cast(group.thread_rank()); + const int laneCount = static_cast(group.num_threads()); + + // Outer loop walks set bits of seed.bonds via __ffs/__ffsll. Each + // iteration handles one query bond q; if q is already mapped (parent + // saw it) we skip, otherwise we extend the match by one bond + + // possibly one atom. + for (int wordIdx = 0; wordIdx < kBondWords; ++wordIdx) { + BondWord remainingBondBits = seed.bonds[wordIdx]; + while (remainingBondBits != 0) { + int bitPosInWord; + if constexpr (sizeof(BondWord) == 4) { + bitPosInWord = __ffs(static_cast(remainingBondBits)) - 1; + } else { + bitPosInWord = __ffsll(static_cast(remainingBondBits)) - 1; + } + const int queryBondIdx = wordIdx * kBondBitsPerWord + bitPosInWord; + remainingBondBits &= remainingBondBits - 1; // clear lowest set bit + + // Keep this control decision uniform across the group; diverging before + // the later ballot/shuffle would make the winning lane undefined. + int mappedTargetBond = kUnmappedTargetIdx; + if (laneRank == 0) + mappedTargetBond = match.targetBondIdx[queryBondIdx]; + mappedTargetBond = group.shfl(mappedTargetBond, 0); + if (mappedTargetBond != kUnmappedTargetIdx) + continue; + + // Decode this bond's query endpoints from the packed (u<<16 | v). + const std::uint32_t queryEndpoints = queryTopology.bondEndpoints[queryBondIdx]; + const int queryEndpointU = static_cast(queryEndpoints >> kBondEndpointShift); + const int queryEndpointV = static_cast(queryEndpoints & kBondEndpointMask); + + // Look up the parent's atom mapping for both endpoints; either + // may already be mapped (from an earlier bond) or still unmapped + // (this bond is the one bringing it in). + int targetForQueryU = kUnmappedTargetIdx; + int targetForQueryV = kUnmappedTargetIdx; + if (laneRank == 0) { + targetForQueryU = match.targetAtomIdx[queryEndpointU]; + targetForQueryV = match.targetAtomIdx[queryEndpointV]; + } + targetForQueryU = group.shfl(targetForQueryU, 0); + targetForQueryV = group.shfl(targetForQueryV, 0); + const bool queryUIsMapped = targetForQueryU != kUnmappedTargetIdx; + const bool queryVIsMapped = targetForQueryV != kUnmappedTargetIdx; + + // Both endpoints unmapped means the seed is missing an earlier + // bond that should have brought one of them in. Phase 1 always + // anchors both initial atoms before pushing, so this never hits + // for well-formed seeds; defensive return false. + if (!queryUIsMapped && !queryVIsMapped) + return false; + + // Each lane records its first compatible target bond, and (for + // the atom-adding case) the resulting new target atom. After + // the per-lane scan, the warp ballots and the lowest-rank winner + // is broadcast as the committed extension. + int chosenTargetBond = -1; + int chosenTargetAtomForUnmapped = -1; // -1 means ring-closing. + + if (queryUIsMapped && queryVIsMapped) { + // Ring-closing: the new bond connects two atoms that are both + // already in the seed's mapping. Find a target bond whose + // endpoints are exactly the pair { targetForQueryU, targetForQueryV }. + const int srcTargetAtom = targetForQueryU; + const int dstTargetAtom = targetForQueryV; + if constexpr (topologyHasAdjacencyBondIndices()) { + if (srcTargetAtom >= 0 && srcTargetAtom < targetTopology.numAtoms) { + const int begin = static_cast(targetTopology.rowOffsets[srcTargetAtom]); + const int end = static_cast(targetTopology.rowOffsets[srcTargetAtom + 1]); + for (int adjIdx = begin + laneRank; adjIdx < end && chosenTargetBond < 0; adjIdx += laneCount) { + const int otherTargetAtom = static_cast(targetTopology.colIndices[adjIdx]); + if (otherTargetAtom != dstTargetAtom) + continue; + const int targetBondIdx = static_cast(targetTopology.bondIndices[adjIdx]); + if (targetBondIdx < 0 || targetBondIdx >= targetTopology.numBonds || targetBondIdx >= maxTB) { + continue; + } + const TargetBondWord visitedBondsWord = match.visitedTargetBonds[targetBondIdx / kTargetBondBitsPerWord]; + if ((visitedBondsWord >> (targetBondIdx % kTargetBondBitsPerWord)) & 1) { + continue; + } + if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) + continue; + chosenTargetBond = targetBondIdx; + } + } + } else { + const bool scanAdjacency = targetTopology.rowOffsets != nullptr && targetTopology.colIndices != nullptr && + targetTopology.bondIndices != nullptr && srcTargetAtom >= 0 && + srcTargetAtom < targetTopology.numAtoms; + if (scanAdjacency) { + const int begin = static_cast(targetTopology.rowOffsets[srcTargetAtom]); + const int end = static_cast(targetTopology.rowOffsets[srcTargetAtom + 1]); + for (int adjIdx = begin + laneRank; adjIdx < end && chosenTargetBond < 0; adjIdx += laneCount) { + const int otherTargetAtom = static_cast(targetTopology.colIndices[adjIdx]); + if (otherTargetAtom != dstTargetAtom) + continue; + const int targetBondIdx = static_cast(targetTopology.bondIndices[adjIdx]); + if (targetBondIdx < 0 || targetBondIdx >= targetTopology.numBonds || targetBondIdx >= maxTB) { + continue; + } + const TargetBondWord visitedBondsWord = match.visitedTargetBonds[targetBondIdx / kTargetBondBitsPerWord]; + if ((visitedBondsWord >> (targetBondIdx % kTargetBondBitsPerWord)) & 1) { + continue; + } + if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) + continue; + chosenTargetBond = targetBondIdx; + } + } else { + for (int targetBondIdx = laneRank; targetBondIdx < targetTopology.numBonds && chosenTargetBond < 0; + targetBondIdx += laneCount) { + // Skip target bonds already used by the parent's match. + const TargetBondWord visitedBondsWord = match.visitedTargetBonds[targetBondIdx / kTargetBondBitsPerWord]; + if ((visitedBondsWord >> (targetBondIdx % kTargetBondBitsPerWord)) & 1) { + continue; + } + const std::uint32_t targetEndpoints = targetTopology.bondEndpoints[targetBondIdx]; + const int targetEndpointU = static_cast(targetEndpoints >> kBondEndpointShift); + const int targetEndpointV = static_cast(targetEndpoints & kBondEndpointMask); + // Match either orientation -- target bonds are undirected. + const bool endpointsMatch = (targetEndpointU == srcTargetAtom && targetEndpointV == dstTargetAtom) || + (targetEndpointU == dstTargetAtom && targetEndpointV == srcTargetAtom); + if (!endpointsMatch) + continue; + if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) + continue; + chosenTargetBond = targetBondIdx; + } + } + } + } else { + // Atom-adding: one query endpoint (`src`) is already mapped to + // a target atom; the other (`dst`) is what we're trying to + // place. Scan target bonds incident to srcTargetAtom for one + // whose other endpoint is unvisited, atom-table-compatible + // with the unmapped query atom, and bond-table-compatible. + const int unmappedQueryAtom = queryUIsMapped ? queryEndpointV : queryEndpointU; + const int srcTargetAtom = queryUIsMapped ? targetForQueryU : targetForQueryV; + if constexpr (topologyHasAdjacencyBondIndices()) { + if (srcTargetAtom >= 0 && srcTargetAtom < targetTopology.numAtoms) { + const int begin = static_cast(targetTopology.rowOffsets[srcTargetAtom]); + const int end = static_cast(targetTopology.rowOffsets[srcTargetAtom + 1]); + for (int adjIdx = begin + laneRank; adjIdx < end && chosenTargetBond < 0; adjIdx += laneCount) { + const int targetBondIdx = static_cast(targetTopology.bondIndices[adjIdx]); + if (targetBondIdx < 0 || targetBondIdx >= targetTopology.numBonds || targetBondIdx >= maxTB) { + continue; + } + const TargetBondWord visitedBondsWord = match.visitedTargetBonds[targetBondIdx / kTargetBondBitsPerWord]; + if ((visitedBondsWord >> (targetBondIdx % kTargetBondBitsPerWord)) & 1) { + continue; + } + const int candidateTargetAtom = static_cast(targetTopology.colIndices[adjIdx]); + if (candidateTargetAtom < 0 || candidateTargetAtom >= targetTopology.numAtoms || + candidateTargetAtom >= maxTA) { + continue; + } + const TargetAtomWord visitedAtomsWord = + match.visitedTargetAtoms[candidateTargetAtom / kTargetAtomBitsPerWord]; + if ((visitedAtomsWord >> (candidateTargetAtom % kTargetAtomBitsPerWord)) & 1) { + continue; + } + if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) + continue; + if (!tables.atoms.testBit(unmappedQueryAtom, candidateTargetAtom)) + continue; + chosenTargetBond = targetBondIdx; + chosenTargetAtomForUnmapped = candidateTargetAtom; + } + } + } else { + const bool scanAdjacency = targetTopology.rowOffsets != nullptr && targetTopology.colIndices != nullptr && + targetTopology.bondIndices != nullptr && srcTargetAtom >= 0 && + srcTargetAtom < targetTopology.numAtoms; + if (scanAdjacency) { + const int begin = static_cast(targetTopology.rowOffsets[srcTargetAtom]); + const int end = static_cast(targetTopology.rowOffsets[srcTargetAtom + 1]); + for (int adjIdx = begin + laneRank; adjIdx < end && chosenTargetBond < 0; adjIdx += laneCount) { + const int targetBondIdx = static_cast(targetTopology.bondIndices[adjIdx]); + if (targetBondIdx < 0 || targetBondIdx >= targetTopology.numBonds || targetBondIdx >= maxTB) { + continue; + } + const TargetBondWord visitedBondsWord = match.visitedTargetBonds[targetBondIdx / kTargetBondBitsPerWord]; + if ((visitedBondsWord >> (targetBondIdx % kTargetBondBitsPerWord)) & 1) { + continue; + } + const int candidateTargetAtom = static_cast(targetTopology.colIndices[adjIdx]); + if (candidateTargetAtom < 0 || candidateTargetAtom >= targetTopology.numAtoms || + candidateTargetAtom >= maxTA) { + continue; + } + const TargetAtomWord visitedAtomsWord = + match.visitedTargetAtoms[candidateTargetAtom / kTargetAtomBitsPerWord]; + if ((visitedAtomsWord >> (candidateTargetAtom % kTargetAtomBitsPerWord)) & 1) { + continue; + } + if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) + continue; + if (!tables.atoms.testBit(unmappedQueryAtom, candidateTargetAtom)) + continue; + chosenTargetBond = targetBondIdx; + chosenTargetAtomForUnmapped = candidateTargetAtom; + } + } else { + for (int targetBondIdx = laneRank; targetBondIdx < targetTopology.numBonds && chosenTargetBond < 0; + targetBondIdx += laneCount) { + const TargetBondWord visitedBondsWord = match.visitedTargetBonds[targetBondIdx / kTargetBondBitsPerWord]; + if ((visitedBondsWord >> (targetBondIdx % kTargetBondBitsPerWord)) & 1) { + continue; + } + const std::uint32_t targetEndpoints = targetTopology.bondEndpoints[targetBondIdx]; + const int targetEndpointU = static_cast(targetEndpoints >> kBondEndpointShift); + const int targetEndpointV = static_cast(targetEndpoints & kBondEndpointMask); + // Identify the candidate target atom on the far side of the + // bond from srcTargetAtom; skip bonds not incident to it. + int candidateTargetAtom; + if (targetEndpointU == srcTargetAtom) { + candidateTargetAtom = targetEndpointV; + } else if (targetEndpointV == srcTargetAtom) { + candidateTargetAtom = targetEndpointU; + } else { + continue; + } + // Candidate target atom must not already be in the embedding. + const TargetAtomWord visitedAtomsWord = + match.visitedTargetAtoms[candidateTargetAtom / kTargetAtomBitsPerWord]; + if ((visitedAtomsWord >> (candidateTargetAtom % kTargetAtomBitsPerWord)) & 1) { + continue; + } + if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) + continue; + if (!tables.atoms.testBit(unmappedQueryAtom, candidateTargetAtom)) + continue; + chosenTargetBond = targetBondIdx; + chosenTargetAtomForUnmapped = candidateTargetAtom; + } + } + } + } + + // Warp-wide pick: ballot for any lane with a candidate, broadcast + // the lowest-rank winner's choice to all lanes. No candidate + // anywhere -> this bond can't be extended -> abandon the seed. + const unsigned ballot = group.ballot(chosenTargetBond >= 0 ? 1u : 0u); + if (ballot == 0u) + return false; + const int firstWinningLane = __ffs(ballot) - 1; + const int committedTargetBond = group.shfl(chosenTargetBond, firstWinningLane); + const int committedTargetAtom = group.shfl(chosenTargetAtomForUnmapped, firstWinningLane); + if (committedTargetBond < 0 || committedTargetBond >= targetTopology.numBonds || committedTargetBond >= maxTB) + return false; + if (committedTargetAtom < -1 || committedTargetAtom >= targetTopology.numAtoms || committedTargetAtom >= maxTA) + return false; + + // Lane 0 commits the new mapping into the shared MatchResult. + // Subsequent bonds in this same call will read the updated + // visited bitsets after group.sync below. + if (laneRank == 0) { + match.targetBondIdx[queryBondIdx] = static_cast(committedTargetBond); + match.visitedTargetBonds[committedTargetBond / kTargetBondBitsPerWord] |= + static_cast(1) << (committedTargetBond % kTargetBondBitsPerWord); + match.matchedBondSize += 1; + match.empty = false; + if (committedTargetAtom >= 0) { + // Atom-adding case: also commit the new atom mapping. + const int unmappedQueryAtom = queryUIsMapped ? queryEndpointV : queryEndpointU; + match.targetAtomIdx[unmappedQueryAtom] = static_cast(committedTargetAtom); + match.visitedTargetAtoms[committedTargetAtom / kTargetAtomBitsPerWord] |= + static_cast(1) << (committedTargetAtom % kTargetAtomBitsPerWord); + match.matchedAtomSize += 1; + } + } + group.sync(); + } + } + return true; +} + +} // namespace fmcs +} // namespace mcs + +#endif // FMCS_CUDA_FMCS_MATCH_CUH diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f4480ba1..f4647a0c 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -81,6 +81,11 @@ target_compile_options(test_fmcs_grow add_executable(test_fmcs_foundations test_fmcs_foundations.cu) target_link_libraries(test_fmcs_foundations PRIVATE mcs_fmcs_foundations) +add_executable(test_fmcs_match test_fmcs_match.cu) +target_include_directories(test_fmcs_match PRIVATE ${CMAKE_SOURCE_DIR}/src/mcs) +target_link_libraries(test_fmcs_match PRIVATE mcs_fmcs_foundations + device_vector) + add_executable(test_openmp_helpers test_openmp_helpers.cpp) target_link_libraries(test_openmp_helpers PRIVATE openmp_helpers OpenMP::OpenMP_CXX) @@ -428,6 +433,7 @@ set(TEST_LIST test_fmcs_primitives test_fmcs_foundations test_fmcs_grow + test_fmcs_match test_work_splitting test_fire_minimizer test_fire_minimizer_permol) diff --git a/tests/test_fmcs_match.cu b/tests/test_fmcs_match.cu new file mode 100644 index 00000000..5e933e69 --- /dev/null +++ b/tests/test_fmcs_match.cu @@ -0,0 +1,645 @@ +// SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +// SPDX-License-Identifier: Apache-2.0 +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. +// +// Unit tests for the per-helper device functions in fmcs_cuda/. Each test +// launches a tiny __global__ driver that constructs inputs, calls one +// helper, and copies results back to host memory for assertion. This +// file is populated incrementally as Steps 1-5 land real implementations. + +#include +#include + +#include +#include +#include +#include + +#include "fmcs_cuda/fmcs_match.cuh" +#include "fmcs_cuda/fmcs_match_tables.cuh" +#include "fmcs_cuda/fmcs_seed.cuh" +#include "src/utils/device_vector.h" + +namespace { + +using nvMolKit::AsyncDevicePtr; +using nvMolKit::AsyncDeviceVector; + +using mcs::fmcs::MatchResult; +using mcs::fmcs::Seed; + +} // namespace + +// --------------------------------------------------------------------------- +// matchSingleBondWithinThread / matchIncrementalFastCooperative +// --------------------------------------------------------------------------- + +namespace { + +using mcs::fmcs::MatchTableDevice; +using mcs::fmcs::PairMatchTablesDevice; +using mcs::fmcs::SingleBondMatch; + +// Tiny CSR-view used by the match helpers via duck-typing; satisfies the +// QueryTopology / TargetTopology template requirement (bondEndpoints + +// numAtoms / numBonds), with optional CSR adjacency fields. +struct TestCsrView { + static constexpr bool kHasAdjacencyBondIndices = false; + + const std::uint32_t* bondEndpoints = nullptr; + int numAtoms = 0; + int numBonds = 0; + const std::uint32_t* rowOffsets = nullptr; + const std::uint32_t* colIndices = nullptr; + const std::uint32_t* bondIndices = nullptr; +}; + +// Match tables staged on the host and uploaded to device memory on +// demand. Both atom and bond tables are 32-bit row-packed bitmasks +// (one bit per (q, t) pair). Tests mutate the host staging vectors +// via the set*Bit helpers; device() uploads (if dirty) and returns the +// device view to pass to kernels. +struct ManagedMatchTables { + std::vector atomHost; + std::vector bondHost; + AsyncDeviceVector atomData; + AsyncDeviceVector bondData; + int qNumAtoms = 0; + int qNumBonds = 0; + + void allocate(int qAtoms, int tAtoms, int qBonds, int tBonds) { + qNumAtoms = qAtoms; + qNumBonds = qBonds; + const int atomWordsPerRow = (tAtoms + 31) / 32; + const int bondWordsPerRow = (tBonds + 31) / 32; + atomHost.assign(qAtoms * atomWordsPerRow, 0); + bondHost.assign(qBonds * bondWordsPerRow, 0); + dev_.atoms = MatchTableDevice{nullptr, qAtoms, tAtoms, atomWordsPerRow}; + dev_.bonds = MatchTableDevice{nullptr, qBonds, tBonds, bondWordsPerRow}; + dirty_ = true; + } + + void setAtomBit(int qAtom, int tAtom) { + atomHost[qAtom * dev_.atoms.wordsPerRow + tAtom / 32] |= (1u << (tAtom % 32)); + dirty_ = true; + } + void setBondBit(int qBond, int tBond) { + bondHost[qBond * dev_.bonds.wordsPerRow + tBond / 32] |= (1u << (tBond % 32)); + dirty_ = true; + } + void setAllAtomBits() { + for (int q = 0; q < dev_.atoms.nRows; ++q) + for (int t = 0; t < dev_.atoms.nCols; ++t) + setAtomBit(q, t); + } + void setAllBondBits() { + for (int q = 0; q < dev_.bonds.nRows; ++q) + for (int t = 0; t < dev_.bonds.nCols; ++t) + setBondBit(q, t); + } + + PairMatchTablesDevice device() { + if (dirty_) { + atomData.setFromVector(atomHost); + bondData.setFromVector(bondHost); + dev_.atoms.data = atomData.data(); + dev_.bonds.data = bondData.data(); + dirty_ = false; + } + return dev_; + } + + private: + PairMatchTablesDevice dev_{}; + bool dirty_ = true; +}; + +AsyncDeviceVector makeBondEndpointsDevice(const std::vector>& edges) { + std::vector host(edges.size()); + for (size_t i = 0; i < edges.size(); ++i) { + host[i] = (static_cast(edges[i].first) << 16) | static_cast(edges[i].second); + } + AsyncDeviceVector dev(edges.size()); + dev.copyFromHost(host); + return dev; +} + +// ---- matchSingleBondWithinThread ---- + +struct SingleBondTestOut { + bool ok; + SingleBondMatch match; +}; + +__global__ void matchSingleBondDriver(int qBondIdx, + int tBondIdx, + bool reversed, + const std::uint32_t* qBondEndpoints, + int qNumAtoms, + int qNumBonds, + const std::uint32_t* tBondEndpoints, + int tNumAtoms, + int tNumBonds, + PairMatchTablesDevice tables, + SingleBondTestOut* out) { + if (threadIdx.x != 0 || blockIdx.x != 0) + return; + TestCsrView qView{qBondEndpoints, qNumAtoms, qNumBonds}; + TestCsrView tView{tBondEndpoints, tNumAtoms, tNumBonds}; + SingleBondMatch sm{}; + out->ok = mcs::fmcs::matchSingleBondWithinThread(qBondIdx, tBondIdx, reversed, qView, tView, tables, sm); + out->match = sm; +} + +} // namespace + +TEST(FMCSUnit, MatchSingleBondForwardOrientation) { + // Query bond (0,1), target bond (0,1). All atoms / bonds compatible. + auto qBE = makeBondEndpointsDevice({ + {0, 1} + }); + auto tBE = makeBondEndpointsDevice({ + {0, 1} + }); + ManagedMatchTables tables; + tables.allocate(2, 2, 1, 1); + tables.setAllAtomBits(); + tables.setAllBondBits(); + + AsyncDevicePtr d_out; + matchSingleBondDriver<<<1, 1>>>(0, + 0, + /*reversed=*/false, + qBE.data(), + 2, + 1, + tBE.data(), + 2, + 1, + tables.device(), + d_out.data()); + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + SingleBondTestOut out{}; + d_out.get(out); + + EXPECT_TRUE(out.ok); + EXPECT_EQ(out.match.targetAtomU, 0u); // qU=0 -> tU=0 + EXPECT_EQ(out.match.targetAtomV, 1u); // qV=1 -> tV=1 +} + +TEST(FMCSUnit, MatchSingleBondReverseOrientation) { + auto qBE = makeBondEndpointsDevice({ + {0, 1} + }); + auto tBE = makeBondEndpointsDevice({ + {0, 1} + }); + ManagedMatchTables tables; + tables.allocate(2, 2, 1, 1); + tables.setAllAtomBits(); + tables.setAllBondBits(); + + AsyncDevicePtr d_out; + matchSingleBondDriver<<<1, 1>>>(0, + 0, + /*reversed=*/true, + qBE.data(), + 2, + 1, + tBE.data(), + 2, + 1, + tables.device(), + d_out.data()); + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + SingleBondTestOut out{}; + d_out.get(out); + + EXPECT_TRUE(out.ok); + EXPECT_EQ(out.match.targetAtomU, 1u); // qU=0 -> tV=1 (reversed) + EXPECT_EQ(out.match.targetAtomV, 0u); // qV=1 -> tU=0 +} + +TEST(FMCSUnit, MatchSingleBondBondTableRejection) { + auto qBE = makeBondEndpointsDevice({ + {0, 1} + }); + auto tBE = makeBondEndpointsDevice({ + {0, 1} + }); + ManagedMatchTables tables; + tables.allocate(2, 2, 1, 1); + tables.setAllAtomBits(); + // Deliberately leave bondData all zero -> bond-table reject. + + AsyncDevicePtr d_out; + matchSingleBondDriver<<<1, 1>>>(0, + 0, + /*reversed=*/false, + qBE.data(), + 2, + 1, + tBE.data(), + 2, + 1, + tables.device(), + d_out.data()); + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + SingleBondTestOut out{}; + d_out.get(out); + + EXPECT_FALSE(out.ok); +} + +TEST(FMCSUnit, MatchSingleBondAtomTableRejection) { + auto qBE = makeBondEndpointsDevice({ + {0, 1} + }); + auto tBE = makeBondEndpointsDevice({ + {0, 1} + }); + ManagedMatchTables tables; + tables.allocate(2, 2, 1, 1); + tables.setAllBondBits(); + // Atom 0 compatible with target atom 0, but atom 1 is incompatible + // with target atom 1 -> forward orientation rejects on second atom. + tables.setAtomBit(0, 0); + tables.setAtomBit(0, 1); + tables.setAtomBit(1, 0); + // Note: (1, 1) deliberately left unset. + + AsyncDevicePtr d_out; + matchSingleBondDriver<<<1, 1>>>(0, + 0, + /*reversed=*/false, + qBE.data(), + 2, + 1, + tBE.data(), + 2, + 1, + tables.device(), + d_out.data()); + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + SingleBondTestOut out{}; + d_out.get(out); + + EXPECT_FALSE(out.ok); +} + +// ---- matchIncrementalFastCooperative ---- + +namespace { + +// Helper: build a parent MatchResult that records the mapping +// {qAtom[i] -> tAtom[i]} and {qBond[i] -> tBond[i]} from caller-supplied +// parallel arrays. Used by the incremental tests to set up "what the +// parent already had matched" before adding new bonds. +template +__device__ __forceinline__ void buildParentMatch(mcs::fmcs::MatchResult& match, + const int* qAtoms, + const int* tAtoms, + int nAtomMaps, + const int* qBonds, + const int* tBonds, + int nBondMaps) { + using MatchT = mcs::fmcs::MatchResult; + mcs::fmcs::matchResultClearWithinThread(match); + for (int i = 0; i < nAtomMaps; ++i) { + match.targetAtomIdx[qAtoms[i]] = static_cast(tAtoms[i]); + const int t = tAtoms[i]; + match.visitedTargetAtoms[t / MatchT::kTargetAtomBitsPerWord] |= (typename MatchT::target_atom_word{1}) + << (t % MatchT::kTargetAtomBitsPerWord); + } + for (int i = 0; i < nBondMaps; ++i) { + match.targetBondIdx[qBonds[i]] = static_cast(tBonds[i]); + const int t = tBonds[i]; + match.visitedTargetBonds[t / MatchT::kTargetBondBitsPerWord] |= (typename MatchT::target_bond_word{1}) + << (t % MatchT::kTargetBondBitsPerWord); + } + match.matchedAtomSize = static_cast(nAtomMaps); + match.matchedBondSize = static_cast(nBondMaps); + match.empty = (nAtomMaps == 0 && nBondMaps == 0); +} + +} // namespace + +namespace mcs_fmcs_incremental_test { + +using QueuedT16 = mcs::fmcs::QueuedSeed<16, 16, 16, 16>; + +struct IncrementalTestOut { + bool ok; + QueuedT16 child; +}; + +// One-warp driver: builds parent match in shared mem, then constructs +// the child seed (parent + new bonds) and runs +// matchIncrementalFastCooperative on it. +__global__ void matchIncrementalAtomAddingDriver(const std::uint32_t* qBE, + int qNumAtoms, + int qNumBonds, + const std::uint32_t* tBE, + int tNumAtoms, + int tNumBonds, + PairMatchTablesDevice tables, + IncrementalTestOut* out) { + __shared__ QueuedT16 child; + if (threadIdx.x == 0) { + mcs::fmcs::seedClearWithinThread(child.seed); + // Parent: bond (0,1) mapped 0->0, 1->1, bond 0->target bond 0. + int qA[] = {0, 1}; + int tA[] = {0, 1}; + int qB[] = {0}; + int tB[] = {0}; + buildParentMatch(child.match, qA, tA, 2, qB, tB, 1); + // Child seed has bonds {0, 1} and atoms {0, 1, 2}. Bond 1 is the + // new atom-adding bond (qU=1 already mapped, qV=2 unmapped). + mcs::fmcs::seedAddBondWithinThread(child.seed, 0); + mcs::fmcs::seedAddBondWithinThread(child.seed, 1); + mcs::fmcs::seedAddAtomWithinThread(child.seed, 0); + mcs::fmcs::seedAddAtomWithinThread(child.seed, 1); + mcs::fmcs::seedBeginGrowStepWithinThread(child.seed); + mcs::fmcs::seedAddAtomWithinThread(child.seed, 2); + } + __syncthreads(); + + auto block = cooperative_groups::this_thread_block(); + auto warp = cooperative_groups::tiled_partition<32>(block); + TestCsrView qView{qBE, qNumAtoms, qNumBonds}; + TestCsrView tView{tBE, tNumAtoms, tNumBonds}; + bool ok = mcs::fmcs::matchIncrementalFastCooperative(warp, child.seed, qView, tView, tables, child.match); + __syncthreads(); + + if (threadIdx.x == 0) { + out->ok = ok; + out->child = child; + } +} + +__global__ void matchIncrementalRingClosingDriver(const std::uint32_t* qBE, + int qNumAtoms, + int qNumBonds, + const std::uint32_t* tBE, + int tNumAtoms, + int tNumBonds, + PairMatchTablesDevice tables, + IncrementalTestOut* out) { + __shared__ QueuedT16 child; + if (threadIdx.x == 0) { + mcs::fmcs::seedClearWithinThread(child.seed); + // 4-atom square query: bonds (0,1), (1,2), (2,3), (0,3). Parent + // has all 4 atoms mapped via 3 path bonds; child adds the + // ring-closing 4th bond (0,3). + int qA[] = {0, 1, 2, 3}; + int tA[] = {0, 1, 2, 3}; + int qB[] = {0, 1, 2}; + int tB[] = {0, 1, 2}; + buildParentMatch(child.match, qA, tA, 4, qB, tB, 3); + for (int b : {0, 1, 2, 3}) + mcs::fmcs::seedAddBondWithinThread(child.seed, b); + for (int a : {0, 1, 2, 3}) + mcs::fmcs::seedAddAtomWithinThread(child.seed, a); + } + __syncthreads(); + + auto block = cooperative_groups::this_thread_block(); + auto warp = cooperative_groups::tiled_partition<32>(block); + TestCsrView qView{qBE, qNumAtoms, qNumBonds}; + TestCsrView tView{tBE, tNumAtoms, tNumBonds}; + bool ok = mcs::fmcs::matchIncrementalFastCooperative(warp, child.seed, qView, tView, tables, child.match); + __syncthreads(); + + if (threadIdx.x == 0) { + out->ok = ok; + out->child = child; + } +} + +__global__ void matchIncrementalVisitedConflictDriver(const std::uint32_t* qBE, + int qNumAtoms, + int qNumBonds, + const std::uint32_t* tBE, + int tNumAtoms, + int tNumBonds, + PairMatchTablesDevice tables, + IncrementalTestOut* out) { + __shared__ QueuedT16 child; + if (threadIdx.x == 0) { + mcs::fmcs::seedClearWithinThread(child.seed); + // Query: atoms 0,1,2; bonds (0,1), (0,2). Target: atoms 0,1; one + // bond (0,1). Parent has bond (0,1) mapped (qU=0->t=0, qV=1->t=1). + // Now child wants atom-adding bond (0,2), but the only target bond + // out of t=0 is (0,1) and t=1 is already visited -> must fail. + int qA[] = {0, 1}; + int tA[] = {0, 1}; + int qB[] = {0}; + int tB[] = {0}; + buildParentMatch(child.match, qA, tA, 2, qB, tB, 1); + mcs::fmcs::seedAddBondWithinThread(child.seed, 0); + mcs::fmcs::seedAddBondWithinThread(child.seed, 1); + mcs::fmcs::seedAddAtomWithinThread(child.seed, 0); + mcs::fmcs::seedAddAtomWithinThread(child.seed, 1); + mcs::fmcs::seedBeginGrowStepWithinThread(child.seed); + mcs::fmcs::seedAddAtomWithinThread(child.seed, 2); + } + __syncthreads(); + + auto block = cooperative_groups::this_thread_block(); + auto warp = cooperative_groups::tiled_partition<32>(block); + TestCsrView qView{qBE, qNumAtoms, qNumBonds}; + TestCsrView tView{tBE, tNumAtoms, tNumBonds}; + bool ok = mcs::fmcs::matchIncrementalFastCooperative(warp, child.seed, qView, tView, tables, child.match); + __syncthreads(); + + if (threadIdx.x == 0) { + out->ok = ok; + out->child = child; + } +} + +__global__ void matchIncrementalTwoBondChainDriver(const std::uint32_t* qBE, + int qNumAtoms, + int qNumBonds, + const std::uint32_t* tBE, + int tNumAtoms, + int tNumBonds, + PairMatchTablesDevice tables, + IncrementalTestOut* out) { + __shared__ QueuedT16 child; + if (threadIdx.x == 0) { + mcs::fmcs::seedClearWithinThread(child.seed); + // 4-atom path 0-1-2-3. Parent has only bond (0,1) mapped. Child + // adds bonds (1,2) and (2,3) in one matchIncrementalFast call. + int qA[] = {0, 1}; + int tA[] = {0, 1}; + int qB[] = {0}; + int tB[] = {0}; + buildParentMatch(child.match, qA, tA, 2, qB, tB, 1); + mcs::fmcs::seedAddBondWithinThread(child.seed, 0); + mcs::fmcs::seedAddBondWithinThread(child.seed, 1); + mcs::fmcs::seedAddBondWithinThread(child.seed, 2); + mcs::fmcs::seedAddAtomWithinThread(child.seed, 0); + mcs::fmcs::seedAddAtomWithinThread(child.seed, 1); + mcs::fmcs::seedBeginGrowStepWithinThread(child.seed); + mcs::fmcs::seedAddAtomWithinThread(child.seed, 2); + mcs::fmcs::seedAddAtomWithinThread(child.seed, 3); + } + __syncthreads(); + + auto block = cooperative_groups::this_thread_block(); + auto warp = cooperative_groups::tiled_partition<32>(block); + TestCsrView qView{qBE, qNumAtoms, qNumBonds}; + TestCsrView tView{tBE, tNumAtoms, tNumBonds}; + bool ok = mcs::fmcs::matchIncrementalFastCooperative(warp, child.seed, qView, tView, tables, child.match); + __syncthreads(); + + if (threadIdx.x == 0) { + out->ok = ok; + out->child = child; + } +} + +} // namespace mcs_fmcs_incremental_test + +TEST(FMCSUnit, MatchIncrementalFastAtomAdding) { + using mcs_fmcs_incremental_test::IncrementalTestOut; + using mcs_fmcs_incremental_test::matchIncrementalAtomAddingDriver; + + // Query/target are both a 3-atom path 0-1-2 with bonds (0,1), (1,2). + auto qBE = makeBondEndpointsDevice({ + {0, 1}, + {1, 2} + }); + auto tBE = makeBondEndpointsDevice({ + {0, 1}, + {1, 2} + }); + ManagedMatchTables tables; + tables.allocate(3, 3, 2, 2); + tables.setAllAtomBits(); + tables.setAllBondBits(); + + AsyncDevicePtr d_out; + matchIncrementalAtomAddingDriver<<<1, 32>>>(qBE.data(), 3, 2, tBE.data(), 3, 2, tables.device(), d_out.data()); + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + IncrementalTestOut out{}; + d_out.get(out); + + EXPECT_TRUE(out.ok); + EXPECT_EQ(out.child.match.targetBondIdx[1], 1u); + EXPECT_EQ(out.child.match.targetAtomIdx[2], 2u); + EXPECT_EQ(out.child.match.matchedBondSize, 2); + EXPECT_EQ(out.child.match.matchedAtomSize, 3); +} + +TEST(FMCSUnit, MatchIncrementalFastRingClosing) { + using mcs_fmcs_incremental_test::IncrementalTestOut; + using mcs_fmcs_incremental_test::matchIncrementalRingClosingDriver; + + // 4-atom square with one diagonal-free closure. Bonds: (0,1) (1,2) + // (2,3) (0,3). + auto qBE = makeBondEndpointsDevice({ + {0, 1}, + {1, 2}, + {2, 3}, + {0, 3} + }); + auto tBE = makeBondEndpointsDevice({ + {0, 1}, + {1, 2}, + {2, 3}, + {0, 3} + }); + ManagedMatchTables tables; + tables.allocate(4, 4, 4, 4); + tables.setAllAtomBits(); + tables.setAllBondBits(); + + AsyncDevicePtr d_out; + matchIncrementalRingClosingDriver<<<1, 32>>>(qBE.data(), 4, 4, tBE.data(), 4, 4, tables.device(), d_out.data()); + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + IncrementalTestOut out{}; + d_out.get(out); + + EXPECT_TRUE(out.ok); + EXPECT_EQ(out.child.match.targetBondIdx[3], 3u); + EXPECT_EQ(out.child.match.matchedBondSize, 4); + // No new atoms added by the ring-closing bond. + EXPECT_EQ(out.child.match.matchedAtomSize, 4); +} + +TEST(FMCSUnit, MatchIncrementalFastVisitedConflictFails) { + using mcs_fmcs_incremental_test::IncrementalTestOut; + using mcs_fmcs_incremental_test::matchIncrementalVisitedConflictDriver; + + // Query: 3 atoms / 2 bonds. Target: 2 atoms / 1 bond. Parent has + // bond (0,1) mapped; trying to extend with bond (0,2) forces atom 2 + // onto target atom 1, which is already visited. + auto qBE = makeBondEndpointsDevice({ + {0, 1}, + {0, 2} + }); + auto tBE = makeBondEndpointsDevice({ + {0, 1} + }); + ManagedMatchTables tables; + tables.allocate(3, 2, 2, 1); + tables.setAllAtomBits(); + tables.setAllBondBits(); + + AsyncDevicePtr d_out; + matchIncrementalVisitedConflictDriver<<<1, 32>>>(qBE.data(), 3, 2, tBE.data(), 2, 1, tables.device(), d_out.data()); + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + IncrementalTestOut out{}; + d_out.get(out); + + EXPECT_FALSE(out.ok); +} + +TEST(FMCSUnit, MatchIncrementalFastTwoBondChain) { + using mcs_fmcs_incremental_test::IncrementalTestOut; + using mcs_fmcs_incremental_test::matchIncrementalTwoBondChainDriver; + + // Both sides are the 4-atom path 0-1-2-3 with bonds 0=(0,1), 1=(1,2), + // 2=(2,3). + auto qBE = makeBondEndpointsDevice({ + {0, 1}, + {1, 2}, + {2, 3} + }); + auto tBE = makeBondEndpointsDevice({ + {0, 1}, + {1, 2}, + {2, 3} + }); + ManagedMatchTables tables; + tables.allocate(4, 4, 3, 3); + tables.setAllAtomBits(); + tables.setAllBondBits(); + + AsyncDevicePtr d_out; + matchIncrementalTwoBondChainDriver<<<1, 32>>>(qBE.data(), 4, 3, tBE.data(), 4, 3, tables.device(), d_out.data()); + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + IncrementalTestOut out{}; + d_out.get(out); + + EXPECT_TRUE(out.ok); + EXPECT_EQ(out.child.match.targetBondIdx[1], 1u); + EXPECT_EQ(out.child.match.targetBondIdx[2], 2u); + EXPECT_EQ(out.child.match.targetAtomIdx[2], 2u); + EXPECT_EQ(out.child.match.targetAtomIdx[3], 3u); + EXPECT_EQ(out.child.match.matchedBondSize, 3); + EXPECT_EQ(out.child.match.matchedAtomSize, 4); +} From d982a2b26d4deaddcfc30af434e427885f722b19 Mon Sep 17 00:00:00 2001 From: Kevin Boyd Date: Fri, 24 Jul 2026 13:46:46 -0400 Subject: [PATCH 2/8] Use project-rooted includes in fMCS matching --- src/mcs/fmcs_cuda/fmcs_match.cuh | 4 ++-- tests/test_fmcs_match.cu | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/mcs/fmcs_cuda/fmcs_match.cuh b/src/mcs/fmcs_cuda/fmcs_match.cuh index afdc049b..956deb50 100644 --- a/src/mcs/fmcs_cuda/fmcs_match.cuh +++ b/src/mcs/fmcs_cuda/fmcs_match.cuh @@ -18,8 +18,8 @@ #include -#include "fmcs_cuda/fmcs_match_tables.cuh" -#include "fmcs_cuda/fmcs_seed.cuh" +#include "src/mcs/fmcs_cuda/fmcs_match_tables.cuh" +#include "src/mcs/fmcs_cuda/fmcs_seed.cuh" namespace mcs { namespace fmcs { diff --git a/tests/test_fmcs_match.cu b/tests/test_fmcs_match.cu index 5e933e69..9c2aa141 100644 --- a/tests/test_fmcs_match.cu +++ b/tests/test_fmcs_match.cu @@ -26,9 +26,9 @@ #include #include -#include "fmcs_cuda/fmcs_match.cuh" -#include "fmcs_cuda/fmcs_match_tables.cuh" -#include "fmcs_cuda/fmcs_seed.cuh" +#include "src/mcs/fmcs_cuda/fmcs_match.cuh" +#include "src/mcs/fmcs_cuda/fmcs_match_tables.cuh" +#include "src/mcs/fmcs_cuda/fmcs_seed.cuh" #include "src/utils/device_vector.h" namespace { From 2c147406da951834dc3278a53fa7a70eb528b4d6 Mon Sep 17 00:00:00 2001 From: Kevin Boyd Date: Tue, 28 Jul 2026 13:22:06 -0400 Subject: [PATCH 3/8] Drop the topology-generality concept from fMCS matching There is exactly one topology representation. The QueryTopology / TargetTopology template parameters and the kHasAdjacencyBondIndices trait made it look like there were two, but the only type that ever set the trait false was a test double -- production code has always passed a CSR view with all four arrays populated. The trait dates from when the CSR-adjacency scan was added on top of the original linear bond-endpoint scan and the old path was kept alive behind `if constexpr`. Both paths now land together, so the fallback was dead on arrival: each `else` branch held a verbatim copy of the adjacency loop under a runtime null-check plus a linear scan nothing reaches. Hoist DeviceCsrView into a new fmcs_topology.cuh (along with the bond endpoint pack/unpack constants, which fmcs_grow.cuh had been duplicating), take it concretely in matchSingleBondWithinThread, matchIncrementalFastCooperative, and fillNewBondsCooperative, and delete the trait with both fallback branches. fmcs_match.cuh drops from 453 to 321 lines with no behaviour change. The match tests relied on the fallback -- they never populated CSR at all -- so they now build real CSR arrays via a TestGraph helper and exercise the path that actually ships. --- src/mcs/fmcs_cuda/fmcs_grow.cuh | 12 +- src/mcs/fmcs_cuda/fmcs_match.cuh | 231 +++++---------------- src/mcs/fmcs_cuda/fmcs_topology.cuh | 52 +++++ tests/test_fmcs_grow.cu | 19 +- tests/test_fmcs_match.cu | 304 +++++++++++++++------------- 5 files changed, 272 insertions(+), 346 deletions(-) create mode 100644 src/mcs/fmcs_cuda/fmcs_topology.cuh diff --git a/src/mcs/fmcs_cuda/fmcs_grow.cuh b/src/mcs/fmcs_cuda/fmcs_grow.cuh index 14b6b418..d2e3998f 100644 --- a/src/mcs/fmcs_cuda/fmcs_grow.cuh +++ b/src/mcs/fmcs_cuda/fmcs_grow.cuh @@ -19,14 +19,12 @@ #include #include "src/mcs/fmcs_cuda/fmcs_seed.cuh" +#include "src/mcs/fmcs_cuda/fmcs_topology.cuh" #include "src/mcs/mcs_common/mcs_cooperative_copy.cuh" namespace mcs { namespace fmcs { -constexpr int kGrowBondEndpointShift = 16; -constexpr std::uint32_t kGrowBondEndpointMask = 0xFFFFu; - template __device__ __forceinline__ void seedAddNewBondWithinThread(Seed& seed, const NewBond& bond) { seedAddBondWithinThread(seed, bond.bondIdx); @@ -68,10 +66,10 @@ __device__ __forceinline__ void seedAddNewBondWithinThread(Seed +template __device__ __forceinline__ bool fillNewBondsCooperative(const GroupT& group, const Seed& seed, - const QueryTopology& queryTopology, + const DeviceCsrView& queryTopology, NewBond* outBonds, int* outCount, int maxNewBonds) { @@ -96,8 +94,8 @@ __device__ __forceinline__ bool fillNewBondsCooperative(const GroupT& // Decode this query bond's endpoints. const std::uint32_t queryEndpoints = queryTopology.bondEndpoints[q]; - const int queryEndpointU = static_cast(queryEndpoints >> kGrowBondEndpointShift); - const int queryEndpointV = static_cast(queryEndpoints & kGrowBondEndpointMask); + const int queryEndpointU = static_cast(queryEndpoints >> kBondEndpointShift); + const int queryEndpointV = static_cast(queryEndpoints & kBondEndpointMask); // Bond is a "new boundary" candidate iff at least one endpoint // was added in the most recent grow step. diff --git a/src/mcs/fmcs_cuda/fmcs_match.cuh b/src/mcs/fmcs_cuda/fmcs_match.cuh index 956deb50..4c033a9d 100644 --- a/src/mcs/fmcs_cuda/fmcs_match.cuh +++ b/src/mcs/fmcs_cuda/fmcs_match.cuh @@ -20,17 +20,11 @@ #include "src/mcs/fmcs_cuda/fmcs_match_tables.cuh" #include "src/mcs/fmcs_cuda/fmcs_seed.cuh" +#include "src/mcs/fmcs_cuda/fmcs_topology.cuh" namespace mcs { namespace fmcs { -/// Bond endpoints are stored across the kernel as a single uint32 with -/// the u-endpoint atom index in the high 16 bits and the v-endpoint -/// atom index in the low 16 bits. Tiers cap maxAtoms at 128 so 16 bits -/// per index is plenty. -constexpr int kBondEndpointShift = 16; -constexpr std::uint32_t kBondEndpointMask = 0xFFFFu; - /// Resolved target endpoints for a successful single-(query bond, target /// bond, orientation) compatibility check. Populated by /// @ref matchSingleBondWithinThread on success only; contents are @@ -42,9 +36,7 @@ struct SingleBondMatch { }; template struct FmcsSubstructureScratch { - // Scratch for the RDKit checkIfMatchAndAppend fallback. This is deliberately - // caller-owned shared memory, not function-local state: tier-128 scratch is - // too large to risk compiler-created stack/local memory in the matcher. + // Scratch for the RDKit checkIfMatchAndAppend fallback. std::uint8_t seedAtomList[maxAtoms]; std::uint8_t seedAtoms[maxAtoms]; std::uint8_t seedDegree[maxAtoms]; @@ -70,14 +62,6 @@ __device__ __forceinline__ bool seedContainsBondWithinThread(const Seed> (queryBondIdx % kBondBitsPerWord)) & 1) != 0; } -template __host__ __device__ constexpr bool topologyHasAdjacencyBondIndices() { - if constexpr (requires { Topology::kHasAdjacencyBondIndices; }) { - return Topology::kHasAdjacencyBondIndices; - } else { - return false; - } -} - /// Within-thread: per-lane single-(query bond, target bond, orientation) /// compatibility check used by Phase 1 initial-seed enumeration. Writes /// resolved target atom indices for the two endpoints of @p queryBondIdx @@ -87,15 +71,11 @@ template __host__ __device__ constexpr bool topologyHasAdjacenc /// @p reversed selects the orientation: when false, the query bond's u /// endpoint maps to the target bond's u endpoint; when true, to the /// target bond's v endpoint. -/// -/// @p queryTopology and @p targetTopology must expose a -/// @c bondEndpoints array of packed (u<<16 | v) entries. -template __device__ __forceinline__ bool matchSingleBondWithinThread(const int queryBondIdx, const int targetBondIdx, const bool reversed, - const QueryTopology& queryTopology, - const TargetTopology& targetTopology, + const DeviceCsrView& queryTopology, + const DeviceCsrView& targetTopology, const PairMatchTablesDevice& tables, SingleBondMatch& outMatch) { // Cheap bond-table check first; if the bond labels are incompatible @@ -151,11 +131,11 @@ __device__ __forceinline__ bool matchSingleBondWithinThread(const int /// Any bond that fails to extend causes the function to return false; /// @p match is left in an unspecified state and the caller should /// discard the seed. -template +template __device__ __forceinline__ bool matchIncrementalFastCooperative(const GroupT& group, const Seed& seed, - const QueryTopology& queryTopology, - const TargetTopology& targetTopology, + const DeviceCsrView& queryTopology, + const DeviceCsrView& targetTopology, const PairMatchTablesDevice& tables, MatchResult& match) { using SeedT = Seed; @@ -236,70 +216,24 @@ __device__ __forceinline__ bool matchIncrementalFastCooperative(const GroupT& // endpoints are exactly the pair { targetForQueryU, targetForQueryV }. const int srcTargetAtom = targetForQueryU; const int dstTargetAtom = targetForQueryV; - if constexpr (topologyHasAdjacencyBondIndices()) { - if (srcTargetAtom >= 0 && srcTargetAtom < targetTopology.numAtoms) { - const int begin = static_cast(targetTopology.rowOffsets[srcTargetAtom]); - const int end = static_cast(targetTopology.rowOffsets[srcTargetAtom + 1]); - for (int adjIdx = begin + laneRank; adjIdx < end && chosenTargetBond < 0; adjIdx += laneCount) { - const int otherTargetAtom = static_cast(targetTopology.colIndices[adjIdx]); - if (otherTargetAtom != dstTargetAtom) - continue; - const int targetBondIdx = static_cast(targetTopology.bondIndices[adjIdx]); - if (targetBondIdx < 0 || targetBondIdx >= targetTopology.numBonds || targetBondIdx >= maxTB) { - continue; - } - const TargetBondWord visitedBondsWord = match.visitedTargetBonds[targetBondIdx / kTargetBondBitsPerWord]; - if ((visitedBondsWord >> (targetBondIdx % kTargetBondBitsPerWord)) & 1) { - continue; - } - if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) - continue; - chosenTargetBond = targetBondIdx; - } - } - } else { - const bool scanAdjacency = targetTopology.rowOffsets != nullptr && targetTopology.colIndices != nullptr && - targetTopology.bondIndices != nullptr && srcTargetAtom >= 0 && - srcTargetAtom < targetTopology.numAtoms; - if (scanAdjacency) { - const int begin = static_cast(targetTopology.rowOffsets[srcTargetAtom]); - const int end = static_cast(targetTopology.rowOffsets[srcTargetAtom + 1]); - for (int adjIdx = begin + laneRank; adjIdx < end && chosenTargetBond < 0; adjIdx += laneCount) { - const int otherTargetAtom = static_cast(targetTopology.colIndices[adjIdx]); - if (otherTargetAtom != dstTargetAtom) - continue; - const int targetBondIdx = static_cast(targetTopology.bondIndices[adjIdx]); - if (targetBondIdx < 0 || targetBondIdx >= targetTopology.numBonds || targetBondIdx >= maxTB) { - continue; - } - const TargetBondWord visitedBondsWord = match.visitedTargetBonds[targetBondIdx / kTargetBondBitsPerWord]; - if ((visitedBondsWord >> (targetBondIdx % kTargetBondBitsPerWord)) & 1) { - continue; - } - if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) - continue; - chosenTargetBond = targetBondIdx; + if (srcTargetAtom >= 0 && srcTargetAtom < targetTopology.numAtoms) { + const int begin = static_cast(targetTopology.rowOffsets[srcTargetAtom]); + const int end = static_cast(targetTopology.rowOffsets[srcTargetAtom + 1]); + for (int adjIdx = begin + laneRank; adjIdx < end && chosenTargetBond < 0; adjIdx += laneCount) { + const int otherTargetAtom = static_cast(targetTopology.colIndices[adjIdx]); + if (otherTargetAtom != dstTargetAtom) + continue; + const int targetBondIdx = static_cast(targetTopology.bondIndices[adjIdx]); + if (targetBondIdx < 0 || targetBondIdx >= targetTopology.numBonds || targetBondIdx >= maxTB) { + continue; } - } else { - for (int targetBondIdx = laneRank; targetBondIdx < targetTopology.numBonds && chosenTargetBond < 0; - targetBondIdx += laneCount) { - // Skip target bonds already used by the parent's match. - const TargetBondWord visitedBondsWord = match.visitedTargetBonds[targetBondIdx / kTargetBondBitsPerWord]; - if ((visitedBondsWord >> (targetBondIdx % kTargetBondBitsPerWord)) & 1) { - continue; - } - const std::uint32_t targetEndpoints = targetTopology.bondEndpoints[targetBondIdx]; - const int targetEndpointU = static_cast(targetEndpoints >> kBondEndpointShift); - const int targetEndpointV = static_cast(targetEndpoints & kBondEndpointMask); - // Match either orientation -- target bonds are undirected. - const bool endpointsMatch = (targetEndpointU == srcTargetAtom && targetEndpointV == dstTargetAtom) || - (targetEndpointU == dstTargetAtom && targetEndpointV == srcTargetAtom); - if (!endpointsMatch) - continue; - if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) - continue; - chosenTargetBond = targetBondIdx; + const TargetBondWord visitedBondsWord = match.visitedTargetBonds[targetBondIdx / kTargetBondBitsPerWord]; + if ((visitedBondsWord >> (targetBondIdx % kTargetBondBitsPerWord)) & 1) { + continue; } + if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) + continue; + chosenTargetBond = targetBondIdx; } } } else { @@ -310,103 +244,34 @@ __device__ __forceinline__ bool matchIncrementalFastCooperative(const GroupT& // with the unmapped query atom, and bond-table-compatible. const int unmappedQueryAtom = queryUIsMapped ? queryEndpointV : queryEndpointU; const int srcTargetAtom = queryUIsMapped ? targetForQueryU : targetForQueryV; - if constexpr (topologyHasAdjacencyBondIndices()) { - if (srcTargetAtom >= 0 && srcTargetAtom < targetTopology.numAtoms) { - const int begin = static_cast(targetTopology.rowOffsets[srcTargetAtom]); - const int end = static_cast(targetTopology.rowOffsets[srcTargetAtom + 1]); - for (int adjIdx = begin + laneRank; adjIdx < end && chosenTargetBond < 0; adjIdx += laneCount) { - const int targetBondIdx = static_cast(targetTopology.bondIndices[adjIdx]); - if (targetBondIdx < 0 || targetBondIdx >= targetTopology.numBonds || targetBondIdx >= maxTB) { - continue; - } - const TargetBondWord visitedBondsWord = match.visitedTargetBonds[targetBondIdx / kTargetBondBitsPerWord]; - if ((visitedBondsWord >> (targetBondIdx % kTargetBondBitsPerWord)) & 1) { - continue; - } - const int candidateTargetAtom = static_cast(targetTopology.colIndices[adjIdx]); - if (candidateTargetAtom < 0 || candidateTargetAtom >= targetTopology.numAtoms || - candidateTargetAtom >= maxTA) { - continue; - } - const TargetAtomWord visitedAtomsWord = - match.visitedTargetAtoms[candidateTargetAtom / kTargetAtomBitsPerWord]; - if ((visitedAtomsWord >> (candidateTargetAtom % kTargetAtomBitsPerWord)) & 1) { - continue; - } - if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) - continue; - if (!tables.atoms.testBit(unmappedQueryAtom, candidateTargetAtom)) - continue; - chosenTargetBond = targetBondIdx; - chosenTargetAtomForUnmapped = candidateTargetAtom; + if (srcTargetAtom >= 0 && srcTargetAtom < targetTopology.numAtoms) { + const int begin = static_cast(targetTopology.rowOffsets[srcTargetAtom]); + const int end = static_cast(targetTopology.rowOffsets[srcTargetAtom + 1]); + for (int adjIdx = begin + laneRank; adjIdx < end && chosenTargetBond < 0; adjIdx += laneCount) { + const int targetBondIdx = static_cast(targetTopology.bondIndices[adjIdx]); + if (targetBondIdx < 0 || targetBondIdx >= targetTopology.numBonds || targetBondIdx >= maxTB) { + continue; } - } - } else { - const bool scanAdjacency = targetTopology.rowOffsets != nullptr && targetTopology.colIndices != nullptr && - targetTopology.bondIndices != nullptr && srcTargetAtom >= 0 && - srcTargetAtom < targetTopology.numAtoms; - if (scanAdjacency) { - const int begin = static_cast(targetTopology.rowOffsets[srcTargetAtom]); - const int end = static_cast(targetTopology.rowOffsets[srcTargetAtom + 1]); - for (int adjIdx = begin + laneRank; adjIdx < end && chosenTargetBond < 0; adjIdx += laneCount) { - const int targetBondIdx = static_cast(targetTopology.bondIndices[adjIdx]); - if (targetBondIdx < 0 || targetBondIdx >= targetTopology.numBonds || targetBondIdx >= maxTB) { - continue; - } - const TargetBondWord visitedBondsWord = match.visitedTargetBonds[targetBondIdx / kTargetBondBitsPerWord]; - if ((visitedBondsWord >> (targetBondIdx % kTargetBondBitsPerWord)) & 1) { - continue; - } - const int candidateTargetAtom = static_cast(targetTopology.colIndices[adjIdx]); - if (candidateTargetAtom < 0 || candidateTargetAtom >= targetTopology.numAtoms || - candidateTargetAtom >= maxTA) { - continue; - } - const TargetAtomWord visitedAtomsWord = - match.visitedTargetAtoms[candidateTargetAtom / kTargetAtomBitsPerWord]; - if ((visitedAtomsWord >> (candidateTargetAtom % kTargetAtomBitsPerWord)) & 1) { - continue; - } - if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) - continue; - if (!tables.atoms.testBit(unmappedQueryAtom, candidateTargetAtom)) - continue; - chosenTargetBond = targetBondIdx; - chosenTargetAtomForUnmapped = candidateTargetAtom; + const TargetBondWord visitedBondsWord = match.visitedTargetBonds[targetBondIdx / kTargetBondBitsPerWord]; + if ((visitedBondsWord >> (targetBondIdx % kTargetBondBitsPerWord)) & 1) { + continue; + } + const int candidateTargetAtom = static_cast(targetTopology.colIndices[adjIdx]); + if (candidateTargetAtom < 0 || candidateTargetAtom >= targetTopology.numAtoms || + candidateTargetAtom >= maxTA) { + continue; } - } else { - for (int targetBondIdx = laneRank; targetBondIdx < targetTopology.numBonds && chosenTargetBond < 0; - targetBondIdx += laneCount) { - const TargetBondWord visitedBondsWord = match.visitedTargetBonds[targetBondIdx / kTargetBondBitsPerWord]; - if ((visitedBondsWord >> (targetBondIdx % kTargetBondBitsPerWord)) & 1) { - continue; - } - const std::uint32_t targetEndpoints = targetTopology.bondEndpoints[targetBondIdx]; - const int targetEndpointU = static_cast(targetEndpoints >> kBondEndpointShift); - const int targetEndpointV = static_cast(targetEndpoints & kBondEndpointMask); - // Identify the candidate target atom on the far side of the - // bond from srcTargetAtom; skip bonds not incident to it. - int candidateTargetAtom; - if (targetEndpointU == srcTargetAtom) { - candidateTargetAtom = targetEndpointV; - } else if (targetEndpointV == srcTargetAtom) { - candidateTargetAtom = targetEndpointU; - } else { - continue; - } - // Candidate target atom must not already be in the embedding. - const TargetAtomWord visitedAtomsWord = - match.visitedTargetAtoms[candidateTargetAtom / kTargetAtomBitsPerWord]; - if ((visitedAtomsWord >> (candidateTargetAtom % kTargetAtomBitsPerWord)) & 1) { - continue; - } - if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) - continue; - if (!tables.atoms.testBit(unmappedQueryAtom, candidateTargetAtom)) - continue; - chosenTargetBond = targetBondIdx; - chosenTargetAtomForUnmapped = candidateTargetAtom; + const TargetAtomWord visitedAtomsWord = + match.visitedTargetAtoms[candidateTargetAtom / kTargetAtomBitsPerWord]; + if ((visitedAtomsWord >> (candidateTargetAtom % kTargetAtomBitsPerWord)) & 1) { + continue; } + if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) + continue; + if (!tables.atoms.testBit(unmappedQueryAtom, candidateTargetAtom)) + continue; + chosenTargetBond = targetBondIdx; + chosenTargetAtomForUnmapped = candidateTargetAtom; } } } diff --git a/src/mcs/fmcs_cuda/fmcs_topology.cuh b/src/mcs/fmcs_cuda/fmcs_topology.cuh new file mode 100644 index 00000000..dd26e14b --- /dev/null +++ b/src/mcs/fmcs_cuda/fmcs_topology.cuh @@ -0,0 +1,52 @@ +// SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +// SPDX-License-Identifier: Apache-2.0 +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#ifndef FMCS_CUDA_FMCS_TOPOLOGY_CUH +#define FMCS_CUDA_FMCS_TOPOLOGY_CUH + +#include + +namespace mcs { +namespace fmcs { + +/// Bond endpoints are stored across the kernel as a single uint32 with +/// the u-endpoint atom index in the high 16 bits and the v-endpoint +/// atom index in the low 16 bits. Tiers cap maxAtoms at 128 so 16 bits +/// per index is plenty. +constexpr int kBondEndpointShift = 16; +constexpr std::uint32_t kBondEndpointMask = 0xFFFFu; + +/// Non-owning view over one side's CSR + bond-endpoint arrays. Passed to +/// the matcher and grow helpers as the query- or target-side topology. +/// +/// All four arrays are required; helpers index them unconditionally. +/// @c rowOffsets has @c numAtoms + 1 entries; @c colIndices and +/// @c bondIndices are parallel and hold one entry per directed CSR edge +/// (the neighbour atom and the undirected bond id respectively). +/// @c bondEndpoints has one packed (u << 16 | v) entry per undirected +/// bond, ordered to match the bond dimension of the match tables. +struct DeviceCsrView { + const std::uint32_t* rowOffsets = nullptr; + const std::uint32_t* colIndices = nullptr; + const std::uint32_t* bondIndices = nullptr; + const std::uint32_t* bondEndpoints = nullptr; + int numAtoms = 0; + int numBonds = 0; +}; + +} // namespace fmcs +} // namespace mcs + +#endif // FMCS_CUDA_FMCS_TOPOLOGY_CUH diff --git a/tests/test_fmcs_grow.cu b/tests/test_fmcs_grow.cu index 6f23af04..0f9962fe 100644 --- a/tests/test_fmcs_grow.cu +++ b/tests/test_fmcs_grow.cu @@ -18,12 +18,6 @@ namespace { using nvMolKit::AsyncDevicePtr; using nvMolKit::AsyncDeviceVector; -struct TestCsrView { - const std::uint32_t* bondEndpoints = nullptr; - int numAtoms = 0; - int numBonds = 0; -}; - AsyncDeviceVector makeBondEndpointsDevice(const std::vector>& edges) { std::vector host(edges.size()); for (std::size_t i = 0; i < edges.size(); ++i) { @@ -74,10 +68,15 @@ __device__ __forceinline__ void fillNewBondsRun(SeedSetup&& setup, } __syncthreads(); - auto block = cooperative_groups::this_thread_block(); - auto warp = cooperative_groups::tiled_partition<32>(block); - TestCsrView qView{qBondEndpoints, qNumAtoms, qNumBonds}; - bool ok = mcs::fmcs::fillNewBondsCooperative(warp, seed, qView, bonds, &count, maxNewBonds); + auto block = cooperative_groups::this_thread_block(); + auto warp = cooperative_groups::tiled_partition<32>(block); + // fillNewBondsCooperative reads only the bond-endpoint side of the + // topology, so the CSR arrays stay null here. + mcs::fmcs::DeviceCsrView qView{}; + qView.bondEndpoints = qBondEndpoints; + qView.numAtoms = qNumAtoms; + qView.numBonds = qNumBonds; + bool ok = mcs::fmcs::fillNewBondsCooperative(warp, seed, qView, bonds, &count, maxNewBonds); __syncthreads(); if (threadIdx.x == 0) { diff --git a/tests/test_fmcs_match.cu b/tests/test_fmcs_match.cu index 9c2aa141..f79a9974 100644 --- a/tests/test_fmcs_match.cu +++ b/tests/test_fmcs_match.cu @@ -24,6 +24,7 @@ #include #include #include +#include #include #include "src/mcs/fmcs_cuda/fmcs_match.cuh" @@ -51,19 +52,7 @@ using mcs::fmcs::MatchTableDevice; using mcs::fmcs::PairMatchTablesDevice; using mcs::fmcs::SingleBondMatch; -// Tiny CSR-view used by the match helpers via duck-typing; satisfies the -// QueryTopology / TargetTopology template requirement (bondEndpoints + -// numAtoms / numBonds), with optional CSR adjacency fields. -struct TestCsrView { - static constexpr bool kHasAdjacencyBondIndices = false; - - const std::uint32_t* bondEndpoints = nullptr; - int numAtoms = 0; - int numBonds = 0; - const std::uint32_t* rowOffsets = nullptr; - const std::uint32_t* colIndices = nullptr; - const std::uint32_t* bondIndices = nullptr; -}; +using mcs::fmcs::DeviceCsrView; // Match tables staged on the host and uploaded to device memory on // demand. Both atom and bond tables are 32-bit row-packed bitmasks @@ -125,15 +114,68 @@ struct ManagedMatchTables { bool dirty_ = true; }; -AsyncDeviceVector makeBondEndpointsDevice(const std::vector>& edges) { - std::vector host(edges.size()); - for (size_t i = 0; i < edges.size(); ++i) { - host[i] = (static_cast(edges[i].first) << 16) | static_cast(edges[i].second); +// Owns the device buffers behind a DeviceCsrView. Built from an +// undirected edge list: bond i is edges[i], and the CSR is the symmetric +// expansion of that list (each edge contributes one entry to each +// endpoint's row), matching what the host-side graph builder uploads for +// real molecules. +class TestGraph { + public: + TestGraph(int numAtoms, const std::vector>& edges) { + const int numBonds = static_cast(edges.size()); + + std::vector bondEndpointsHost(numBonds); + for (int i = 0; i < numBonds; ++i) { + bondEndpointsHost[i] = + (static_cast(edges[i].first) << 16) | static_cast(edges[i].second); + } + + // Counting sort into CSR: degree pass, prefix sum, then scatter. + std::vector rowOffsetsHost(numAtoms + 1, 0); + for (const auto& edge : edges) { + ++rowOffsetsHost[edge.first + 1]; + ++rowOffsetsHost[edge.second + 1]; + } + for (int atom = 0; atom < numAtoms; ++atom) { + rowOffsetsHost[atom + 1] += rowOffsetsHost[atom]; + } + + std::vector colIndicesHost(2 * numBonds); + std::vector bondIndicesHost(2 * numBonds); + std::vector cursor(rowOffsetsHost.begin(), rowOffsetsHost.end() - 1); + for (int bond = 0; bond < numBonds; ++bond) { + const int u = edges[bond].first; + const int v = edges[bond].second; + colIndicesHost[cursor[u]] = static_cast(v); + bondIndicesHost[cursor[u]] = static_cast(bond); + ++cursor[u]; + colIndicesHost[cursor[v]] = static_cast(u); + bondIndicesHost[cursor[v]] = static_cast(bond); + ++cursor[v]; + } + + bondEndpoints_.setFromVector(bondEndpointsHost); + rowOffsets_.setFromVector(rowOffsetsHost); + colIndices_.setFromVector(colIndicesHost); + bondIndices_.setFromVector(bondIndicesHost); + + view_.bondEndpoints = bondEndpoints_.data(); + view_.rowOffsets = rowOffsets_.data(); + view_.colIndices = colIndices_.data(); + view_.bondIndices = bondIndices_.data(); + view_.numAtoms = numAtoms; + view_.numBonds = numBonds; } - AsyncDeviceVector dev(edges.size()); - dev.copyFromHost(host); - return dev; -} + + DeviceCsrView view() const { return view_; } + + private: + AsyncDeviceVector bondEndpoints_; + AsyncDeviceVector rowOffsets_; + AsyncDeviceVector colIndices_; + AsyncDeviceVector bondIndices_; + DeviceCsrView view_{}; +}; // ---- matchSingleBondWithinThread ---- @@ -145,18 +187,12 @@ struct SingleBondTestOut { __global__ void matchSingleBondDriver(int qBondIdx, int tBondIdx, bool reversed, - const std::uint32_t* qBondEndpoints, - int qNumAtoms, - int qNumBonds, - const std::uint32_t* tBondEndpoints, - int tNumAtoms, - int tNumBonds, + DeviceCsrView qView, + DeviceCsrView tView, PairMatchTablesDevice tables, SingleBondTestOut* out) { if (threadIdx.x != 0 || blockIdx.x != 0) return; - TestCsrView qView{qBondEndpoints, qNumAtoms, qNumBonds}; - TestCsrView tView{tBondEndpoints, tNumAtoms, tNumBonds}; SingleBondMatch sm{}; out->ok = mcs::fmcs::matchSingleBondWithinThread(qBondIdx, tBondIdx, reversed, qView, tView, tables, sm); out->match = sm; @@ -166,11 +202,13 @@ __global__ void matchSingleBondDriver(int qBondIdx, TEST(FMCSUnit, MatchSingleBondForwardOrientation) { // Query bond (0,1), target bond (0,1). All atoms / bonds compatible. - auto qBE = makeBondEndpointsDevice({ - {0, 1} + TestGraph query(2, + { + {0, 1} }); - auto tBE = makeBondEndpointsDevice({ - {0, 1} + TestGraph target(2, + { + {0, 1} }); ManagedMatchTables tables; tables.allocate(2, 2, 1, 1); @@ -181,12 +219,8 @@ TEST(FMCSUnit, MatchSingleBondForwardOrientation) { matchSingleBondDriver<<<1, 1>>>(0, 0, /*reversed=*/false, - qBE.data(), - 2, - 1, - tBE.data(), - 2, - 1, + query.view(), + target.view(), tables.device(), d_out.data()); ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); @@ -199,11 +233,13 @@ TEST(FMCSUnit, MatchSingleBondForwardOrientation) { } TEST(FMCSUnit, MatchSingleBondReverseOrientation) { - auto qBE = makeBondEndpointsDevice({ - {0, 1} + TestGraph query(2, + { + {0, 1} }); - auto tBE = makeBondEndpointsDevice({ - {0, 1} + TestGraph target(2, + { + {0, 1} }); ManagedMatchTables tables; tables.allocate(2, 2, 1, 1); @@ -214,12 +250,8 @@ TEST(FMCSUnit, MatchSingleBondReverseOrientation) { matchSingleBondDriver<<<1, 1>>>(0, 0, /*reversed=*/true, - qBE.data(), - 2, - 1, - tBE.data(), - 2, - 1, + query.view(), + target.view(), tables.device(), d_out.data()); ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); @@ -232,11 +264,13 @@ TEST(FMCSUnit, MatchSingleBondReverseOrientation) { } TEST(FMCSUnit, MatchSingleBondBondTableRejection) { - auto qBE = makeBondEndpointsDevice({ - {0, 1} + TestGraph query(2, + { + {0, 1} }); - auto tBE = makeBondEndpointsDevice({ - {0, 1} + TestGraph target(2, + { + {0, 1} }); ManagedMatchTables tables; tables.allocate(2, 2, 1, 1); @@ -247,12 +281,8 @@ TEST(FMCSUnit, MatchSingleBondBondTableRejection) { matchSingleBondDriver<<<1, 1>>>(0, 0, /*reversed=*/false, - qBE.data(), - 2, - 1, - tBE.data(), - 2, - 1, + query.view(), + target.view(), tables.device(), d_out.data()); ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); @@ -263,11 +293,13 @@ TEST(FMCSUnit, MatchSingleBondBondTableRejection) { } TEST(FMCSUnit, MatchSingleBondAtomTableRejection) { - auto qBE = makeBondEndpointsDevice({ - {0, 1} + TestGraph query(2, + { + {0, 1} }); - auto tBE = makeBondEndpointsDevice({ - {0, 1} + TestGraph target(2, + { + {0, 1} }); ManagedMatchTables tables; tables.allocate(2, 2, 1, 1); @@ -283,12 +315,8 @@ TEST(FMCSUnit, MatchSingleBondAtomTableRejection) { matchSingleBondDriver<<<1, 1>>>(0, 0, /*reversed=*/false, - qBE.data(), - 2, - 1, - tBE.data(), - 2, - 1, + query.view(), + target.view(), tables.device(), d_out.data()); ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); @@ -347,12 +375,8 @@ struct IncrementalTestOut { // One-warp driver: builds parent match in shared mem, then constructs // the child seed (parent + new bonds) and runs // matchIncrementalFastCooperative on it. -__global__ void matchIncrementalAtomAddingDriver(const std::uint32_t* qBE, - int qNumAtoms, - int qNumBonds, - const std::uint32_t* tBE, - int tNumAtoms, - int tNumBonds, +__global__ void matchIncrementalAtomAddingDriver(DeviceCsrView qView, + DeviceCsrView tView, PairMatchTablesDevice tables, IncrementalTestOut* out) { __shared__ QueuedT16 child; @@ -375,11 +399,9 @@ __global__ void matchIncrementalAtomAddingDriver(const std::uint32_t* qBE, } __syncthreads(); - auto block = cooperative_groups::this_thread_block(); - auto warp = cooperative_groups::tiled_partition<32>(block); - TestCsrView qView{qBE, qNumAtoms, qNumBonds}; - TestCsrView tView{tBE, tNumAtoms, tNumBonds}; - bool ok = mcs::fmcs::matchIncrementalFastCooperative(warp, child.seed, qView, tView, tables, child.match); + auto block = cooperative_groups::this_thread_block(); + auto warp = cooperative_groups::tiled_partition<32>(block); + bool ok = mcs::fmcs::matchIncrementalFastCooperative(warp, child.seed, qView, tView, tables, child.match); __syncthreads(); if (threadIdx.x == 0) { @@ -388,12 +410,8 @@ __global__ void matchIncrementalAtomAddingDriver(const std::uint32_t* qBE, } } -__global__ void matchIncrementalRingClosingDriver(const std::uint32_t* qBE, - int qNumAtoms, - int qNumBonds, - const std::uint32_t* tBE, - int tNumAtoms, - int tNumBonds, +__global__ void matchIncrementalRingClosingDriver(DeviceCsrView qView, + DeviceCsrView tView, PairMatchTablesDevice tables, IncrementalTestOut* out) { __shared__ QueuedT16 child; @@ -414,11 +432,9 @@ __global__ void matchIncrementalRingClosingDriver(const std::uint32_t* qBE, } __syncthreads(); - auto block = cooperative_groups::this_thread_block(); - auto warp = cooperative_groups::tiled_partition<32>(block); - TestCsrView qView{qBE, qNumAtoms, qNumBonds}; - TestCsrView tView{tBE, tNumAtoms, tNumBonds}; - bool ok = mcs::fmcs::matchIncrementalFastCooperative(warp, child.seed, qView, tView, tables, child.match); + auto block = cooperative_groups::this_thread_block(); + auto warp = cooperative_groups::tiled_partition<32>(block); + bool ok = mcs::fmcs::matchIncrementalFastCooperative(warp, child.seed, qView, tView, tables, child.match); __syncthreads(); if (threadIdx.x == 0) { @@ -427,12 +443,8 @@ __global__ void matchIncrementalRingClosingDriver(const std::uint32_t* qBE, } } -__global__ void matchIncrementalVisitedConflictDriver(const std::uint32_t* qBE, - int qNumAtoms, - int qNumBonds, - const std::uint32_t* tBE, - int tNumAtoms, - int tNumBonds, +__global__ void matchIncrementalVisitedConflictDriver(DeviceCsrView qView, + DeviceCsrView tView, PairMatchTablesDevice tables, IncrementalTestOut* out) { __shared__ QueuedT16 child; @@ -456,11 +468,9 @@ __global__ void matchIncrementalVisitedConflictDriver(const std::uint32_t* qBE, } __syncthreads(); - auto block = cooperative_groups::this_thread_block(); - auto warp = cooperative_groups::tiled_partition<32>(block); - TestCsrView qView{qBE, qNumAtoms, qNumBonds}; - TestCsrView tView{tBE, tNumAtoms, tNumBonds}; - bool ok = mcs::fmcs::matchIncrementalFastCooperative(warp, child.seed, qView, tView, tables, child.match); + auto block = cooperative_groups::this_thread_block(); + auto warp = cooperative_groups::tiled_partition<32>(block); + bool ok = mcs::fmcs::matchIncrementalFastCooperative(warp, child.seed, qView, tView, tables, child.match); __syncthreads(); if (threadIdx.x == 0) { @@ -469,12 +479,8 @@ __global__ void matchIncrementalVisitedConflictDriver(const std::uint32_t* qBE, } } -__global__ void matchIncrementalTwoBondChainDriver(const std::uint32_t* qBE, - int qNumAtoms, - int qNumBonds, - const std::uint32_t* tBE, - int tNumAtoms, - int tNumBonds, +__global__ void matchIncrementalTwoBondChainDriver(DeviceCsrView qView, + DeviceCsrView tView, PairMatchTablesDevice tables, IncrementalTestOut* out) { __shared__ QueuedT16 child; @@ -498,11 +504,9 @@ __global__ void matchIncrementalTwoBondChainDriver(const std::uint32_t* qBE, } __syncthreads(); - auto block = cooperative_groups::this_thread_block(); - auto warp = cooperative_groups::tiled_partition<32>(block); - TestCsrView qView{qBE, qNumAtoms, qNumBonds}; - TestCsrView tView{tBE, tNumAtoms, tNumBonds}; - bool ok = mcs::fmcs::matchIncrementalFastCooperative(warp, child.seed, qView, tView, tables, child.match); + auto block = cooperative_groups::this_thread_block(); + auto warp = cooperative_groups::tiled_partition<32>(block); + bool ok = mcs::fmcs::matchIncrementalFastCooperative(warp, child.seed, qView, tView, tables, child.match); __syncthreads(); if (threadIdx.x == 0) { @@ -518,13 +522,15 @@ TEST(FMCSUnit, MatchIncrementalFastAtomAdding) { using mcs_fmcs_incremental_test::matchIncrementalAtomAddingDriver; // Query/target are both a 3-atom path 0-1-2 with bonds (0,1), (1,2). - auto qBE = makeBondEndpointsDevice({ - {0, 1}, - {1, 2} + TestGraph query(3, + { + {0, 1}, + {1, 2} }); - auto tBE = makeBondEndpointsDevice({ - {0, 1}, - {1, 2} + TestGraph target(3, + { + {0, 1}, + {1, 2} }); ManagedMatchTables tables; tables.allocate(3, 3, 2, 2); @@ -532,7 +538,7 @@ TEST(FMCSUnit, MatchIncrementalFastAtomAdding) { tables.setAllBondBits(); AsyncDevicePtr d_out; - matchIncrementalAtomAddingDriver<<<1, 32>>>(qBE.data(), 3, 2, tBE.data(), 3, 2, tables.device(), d_out.data()); + matchIncrementalAtomAddingDriver<<<1, 32>>>(query.view(), target.view(), tables.device(), d_out.data()); ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); IncrementalTestOut out{}; d_out.get(out); @@ -550,17 +556,19 @@ TEST(FMCSUnit, MatchIncrementalFastRingClosing) { // 4-atom square with one diagonal-free closure. Bonds: (0,1) (1,2) // (2,3) (0,3). - auto qBE = makeBondEndpointsDevice({ - {0, 1}, - {1, 2}, - {2, 3}, - {0, 3} + TestGraph query(4, + { + {0, 1}, + {1, 2}, + {2, 3}, + {0, 3} }); - auto tBE = makeBondEndpointsDevice({ - {0, 1}, - {1, 2}, - {2, 3}, - {0, 3} + TestGraph target(4, + { + {0, 1}, + {1, 2}, + {2, 3}, + {0, 3} }); ManagedMatchTables tables; tables.allocate(4, 4, 4, 4); @@ -568,7 +576,7 @@ TEST(FMCSUnit, MatchIncrementalFastRingClosing) { tables.setAllBondBits(); AsyncDevicePtr d_out; - matchIncrementalRingClosingDriver<<<1, 32>>>(qBE.data(), 4, 4, tBE.data(), 4, 4, tables.device(), d_out.data()); + matchIncrementalRingClosingDriver<<<1, 32>>>(query.view(), target.view(), tables.device(), d_out.data()); ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); IncrementalTestOut out{}; d_out.get(out); @@ -587,12 +595,14 @@ TEST(FMCSUnit, MatchIncrementalFastVisitedConflictFails) { // Query: 3 atoms / 2 bonds. Target: 2 atoms / 1 bond. Parent has // bond (0,1) mapped; trying to extend with bond (0,2) forces atom 2 // onto target atom 1, which is already visited. - auto qBE = makeBondEndpointsDevice({ - {0, 1}, - {0, 2} + TestGraph query(3, + { + {0, 1}, + {0, 2} }); - auto tBE = makeBondEndpointsDevice({ - {0, 1} + TestGraph target(2, + { + {0, 1} }); ManagedMatchTables tables; tables.allocate(3, 2, 2, 1); @@ -600,7 +610,7 @@ TEST(FMCSUnit, MatchIncrementalFastVisitedConflictFails) { tables.setAllBondBits(); AsyncDevicePtr d_out; - matchIncrementalVisitedConflictDriver<<<1, 32>>>(qBE.data(), 3, 2, tBE.data(), 2, 1, tables.device(), d_out.data()); + matchIncrementalVisitedConflictDriver<<<1, 32>>>(query.view(), target.view(), tables.device(), d_out.data()); ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); IncrementalTestOut out{}; d_out.get(out); @@ -614,15 +624,17 @@ TEST(FMCSUnit, MatchIncrementalFastTwoBondChain) { // Both sides are the 4-atom path 0-1-2-3 with bonds 0=(0,1), 1=(1,2), // 2=(2,3). - auto qBE = makeBondEndpointsDevice({ - {0, 1}, - {1, 2}, - {2, 3} + TestGraph query(4, + { + {0, 1}, + {1, 2}, + {2, 3} }); - auto tBE = makeBondEndpointsDevice({ - {0, 1}, - {1, 2}, - {2, 3} + TestGraph target(4, + { + {0, 1}, + {1, 2}, + {2, 3} }); ManagedMatchTables tables; tables.allocate(4, 4, 3, 3); @@ -630,7 +642,7 @@ TEST(FMCSUnit, MatchIncrementalFastTwoBondChain) { tables.setAllBondBits(); AsyncDevicePtr d_out; - matchIncrementalTwoBondChainDriver<<<1, 32>>>(qBE.data(), 4, 3, tBE.data(), 4, 3, tables.device(), d_out.data()); + matchIncrementalTwoBondChainDriver<<<1, 32>>>(query.view(), target.view(), tables.device(), d_out.data()); ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); IncrementalTestOut out{}; d_out.get(out); From b12721c8a8096314a01a99f989c92ac33e7baacd Mon Sep 17 00:00:00 2001 From: Kevin Boyd Date: Tue, 28 Jul 2026 13:44:49 -0400 Subject: [PATCH 4/8] Clarify greedy fMCS match contract --- src/mcs/fmcs_cuda/fmcs_grow.cuh | 10 ++++++---- src/mcs/fmcs_cuda/fmcs_match.cuh | 34 ++++++++++++++++++-------------- tests/test_fmcs_match.cu | 22 ++++++++++----------- 3 files changed, 36 insertions(+), 30 deletions(-) diff --git a/src/mcs/fmcs_cuda/fmcs_grow.cuh b/src/mcs/fmcs_cuda/fmcs_grow.cuh index d2e3998f..e72d44d7 100644 --- a/src/mcs/fmcs_cuda/fmcs_grow.cuh +++ b/src/mcs/fmcs_cuda/fmcs_grow.cuh @@ -151,10 +151,12 @@ __device__ __forceinline__ bool fillNewBondsCooperative(const GroupT& /// across all lanes of @p group), and within each iteration the only /// truly group-parallel step is the @c warpCopy of @c parent into /// @p childWorkspace. Lane 0 then patches in the singleton bond / -/// atom; the per-bond match call is whatever @p matchFn does (when -/// the kernel passes @ref matchIncrementalFastCooperative the inner -/// target-bond scan IS lane-parallel); @p childSink is typically lane -/// 0 only (cache.insert + queue.push). +/// atom; the per-bond match call is whatever @p matchFn does. The +/// kernel's callback may try @ref tryMatchIncrementalGreedyCooperative, +/// but it must own the exact fallback too: a greedy miss alone must not +/// prune the bond. The greedy helper's inner target-bond scan is +/// lane-parallel; @p childSink is typically lane 0 only (cache.insert + +/// queue.push). /// /// In other words: this function is cooperative only for the /// parent->child copy and for delegating to its callbacks; it does diff --git a/src/mcs/fmcs_cuda/fmcs_match.cuh b/src/mcs/fmcs_cuda/fmcs_match.cuh index 4c033a9d..b800ca40 100644 --- a/src/mcs/fmcs_cuda/fmcs_match.cuh +++ b/src/mcs/fmcs_cuda/fmcs_match.cuh @@ -110,10 +110,10 @@ __device__ __forceinline__ bool matchSingleBondWithinThread(const int return true; } -/// Cooperative: extend @p match by every query bond in @p seed.bonds -/// whose @c match.targetBondIdx[q] is still @ref kUnmappedTargetIdx -/// (i.e., unmapped by the parent's recorded embedding). For each such -/// bond: +/// Cooperative greedy fast path: try to extend @p match by every query +/// bond in @p seed.bonds whose @c match.targetBondIdx[q] is still +/// @ref kUnmappedTargetIdx (i.e., unmapped by the parent's recorded +/// embedding). For each such bond: /// - Both endpoints already mapped -> ring-closing case. The lanes /// of @p group scan target bonds in parallel for one whose endpoint /// pair exactly matches the mapped (queryU, queryV) target atoms @@ -122,22 +122,26 @@ __device__ __forceinline__ bool matchSingleBondWithinThread(const int /// - Exactly one endpoint mapped -> atom-adding case. Lane-parallel /// scan for a target bond incident to the mapped target atom whose /// other end is unvisited, atom-table-compatible with the unmapped -/// query atom, and bond-table-compatible. First compatible +/// query atom, and bond-table-compatible. The first compatible /// candidate commits both the new bond mapping and the new atom -/// mapping, and marks both visited. +/// mapping, and marks both visited. This path does not backtrack +/// if that greedy choice prevents a later bond from matching. /// - Both endpoints unmapped -> defensive fail (shouldn't occur on /// well-formed seeds, where Phase 1 maps both initial atoms before /// pushing). -/// Any bond that fails to extend causes the function to return false; -/// @p match is left in an unspecified state and the caller should -/// discard the seed. +/// A true return proves that @p match was extended successfully. A false +/// return is inconclusive: the greedy choices may have missed another valid +/// embedding. The caller must run the exact substructure fallback before +/// rejecting the seed or marking a NewBond dead. On false, @p match is left +/// in an unspecified state and must not be reused as a valid embedding. template -__device__ __forceinline__ bool matchIncrementalFastCooperative(const GroupT& group, - const Seed& seed, - const DeviceCsrView& queryTopology, - const DeviceCsrView& targetTopology, - const PairMatchTablesDevice& tables, - MatchResult& match) { +__device__ __forceinline__ bool tryMatchIncrementalGreedyCooperative( + const GroupT& group, + const Seed& seed, + const DeviceCsrView& queryTopology, + const DeviceCsrView& targetTopology, + const PairMatchTablesDevice& tables, + MatchResult& match) { using SeedT = Seed; using MatchT = MatchResult; using BondWord = typename SeedT::bond_word_type; diff --git a/tests/test_fmcs_match.cu b/tests/test_fmcs_match.cu index f79a9974..004da3f0 100644 --- a/tests/test_fmcs_match.cu +++ b/tests/test_fmcs_match.cu @@ -43,7 +43,7 @@ using mcs::fmcs::Seed; } // namespace // --------------------------------------------------------------------------- -// matchSingleBondWithinThread / matchIncrementalFastCooperative +// matchSingleBondWithinThread / tryMatchIncrementalGreedyCooperative // --------------------------------------------------------------------------- namespace { @@ -326,7 +326,7 @@ TEST(FMCSUnit, MatchSingleBondAtomTableRejection) { EXPECT_FALSE(out.ok); } -// ---- matchIncrementalFastCooperative ---- +// ---- tryMatchIncrementalGreedyCooperative ---- namespace { @@ -374,7 +374,7 @@ struct IncrementalTestOut { // One-warp driver: builds parent match in shared mem, then constructs // the child seed (parent + new bonds) and runs -// matchIncrementalFastCooperative on it. +// tryMatchIncrementalGreedyCooperative on it. __global__ void matchIncrementalAtomAddingDriver(DeviceCsrView qView, DeviceCsrView tView, PairMatchTablesDevice tables, @@ -401,7 +401,7 @@ __global__ void matchIncrementalAtomAddingDriver(DeviceCsrView qView, auto block = cooperative_groups::this_thread_block(); auto warp = cooperative_groups::tiled_partition<32>(block); - bool ok = mcs::fmcs::matchIncrementalFastCooperative(warp, child.seed, qView, tView, tables, child.match); + bool ok = mcs::fmcs::tryMatchIncrementalGreedyCooperative(warp, child.seed, qView, tView, tables, child.match); __syncthreads(); if (threadIdx.x == 0) { @@ -434,7 +434,7 @@ __global__ void matchIncrementalRingClosingDriver(DeviceCsrView qView, auto block = cooperative_groups::this_thread_block(); auto warp = cooperative_groups::tiled_partition<32>(block); - bool ok = mcs::fmcs::matchIncrementalFastCooperative(warp, child.seed, qView, tView, tables, child.match); + bool ok = mcs::fmcs::tryMatchIncrementalGreedyCooperative(warp, child.seed, qView, tView, tables, child.match); __syncthreads(); if (threadIdx.x == 0) { @@ -470,7 +470,7 @@ __global__ void matchIncrementalVisitedConflictDriver(DeviceCsrView qVie auto block = cooperative_groups::this_thread_block(); auto warp = cooperative_groups::tiled_partition<32>(block); - bool ok = mcs::fmcs::matchIncrementalFastCooperative(warp, child.seed, qView, tView, tables, child.match); + bool ok = mcs::fmcs::tryMatchIncrementalGreedyCooperative(warp, child.seed, qView, tView, tables, child.match); __syncthreads(); if (threadIdx.x == 0) { @@ -506,7 +506,7 @@ __global__ void matchIncrementalTwoBondChainDriver(DeviceCsrView qView, auto block = cooperative_groups::this_thread_block(); auto warp = cooperative_groups::tiled_partition<32>(block); - bool ok = mcs::fmcs::matchIncrementalFastCooperative(warp, child.seed, qView, tView, tables, child.match); + bool ok = mcs::fmcs::tryMatchIncrementalGreedyCooperative(warp, child.seed, qView, tView, tables, child.match); __syncthreads(); if (threadIdx.x == 0) { @@ -517,7 +517,7 @@ __global__ void matchIncrementalTwoBondChainDriver(DeviceCsrView qView, } // namespace mcs_fmcs_incremental_test -TEST(FMCSUnit, MatchIncrementalFastAtomAdding) { +TEST(FMCSUnit, TryMatchIncrementalGreedyAtomAdding) { using mcs_fmcs_incremental_test::IncrementalTestOut; using mcs_fmcs_incremental_test::matchIncrementalAtomAddingDriver; @@ -550,7 +550,7 @@ TEST(FMCSUnit, MatchIncrementalFastAtomAdding) { EXPECT_EQ(out.child.match.matchedAtomSize, 3); } -TEST(FMCSUnit, MatchIncrementalFastRingClosing) { +TEST(FMCSUnit, TryMatchIncrementalGreedyRingClosing) { using mcs_fmcs_incremental_test::IncrementalTestOut; using mcs_fmcs_incremental_test::matchIncrementalRingClosingDriver; @@ -588,7 +588,7 @@ TEST(FMCSUnit, MatchIncrementalFastRingClosing) { EXPECT_EQ(out.child.match.matchedAtomSize, 4); } -TEST(FMCSUnit, MatchIncrementalFastVisitedConflictFails) { +TEST(FMCSUnit, TryMatchIncrementalGreedyVisitedConflictNeedsFallback) { using mcs_fmcs_incremental_test::IncrementalTestOut; using mcs_fmcs_incremental_test::matchIncrementalVisitedConflictDriver; @@ -618,7 +618,7 @@ TEST(FMCSUnit, MatchIncrementalFastVisitedConflictFails) { EXPECT_FALSE(out.ok); } -TEST(FMCSUnit, MatchIncrementalFastTwoBondChain) { +TEST(FMCSUnit, TryMatchIncrementalGreedyTwoBondChain) { using mcs_fmcs_incremental_test::IncrementalTestOut; using mcs_fmcs_incremental_test::matchIncrementalTwoBondChainDriver; From 8fb960b6b93a1557b09072b44dc27d4b8edb8b94 Mon Sep 17 00:00:00 2001 From: Kevin Boyd Date: Fri, 24 Jul 2026 10:26:44 -0400 Subject: [PATCH 5/8] Add fMCS full-substructure fallback matching --- src/mcs/fmcs_cuda/fmcs_match.cuh | 849 +++++++++++++++++++++++++++++++ tests/CMakeLists.txt | 7 + tests/test_fmcs_substructure.cu | 582 +++++++++++++++++++++ 3 files changed, 1438 insertions(+) create mode 100644 tests/test_fmcs_substructure.cu diff --git a/src/mcs/fmcs_cuda/fmcs_match.cuh b/src/mcs/fmcs_cuda/fmcs_match.cuh index b800ca40..a60c07d7 100644 --- a/src/mcs/fmcs_cuda/fmcs_match.cuh +++ b/src/mcs/fmcs_cuda/fmcs_match.cuh @@ -318,6 +318,855 @@ __device__ __forceinline__ bool tryMatchIncrementalGreedyCooperative( return true; } +template +__device__ __forceinline__ bool findTargetBondBetweenAtomsWithinThread(const int targetAtomA, + const int targetAtomB, + const int queryBondIdx, + const TargetTopology& targetTopology, + const PairMatchTablesDevice& tables, + int& outTargetBondIdx) { + if constexpr (topologyHasAdjacencyBondIndices()) { + if (targetAtomA < 0 || targetAtomA >= targetTopology.numAtoms) { + return false; + } + const int begin = static_cast(targetTopology.rowOffsets[targetAtomA]); + const int end = static_cast(targetTopology.rowOffsets[targetAtomA + 1]); + for (int adjIdx = begin; adjIdx < end; ++adjIdx) { + const int otherTargetAtom = static_cast(targetTopology.colIndices[adjIdx]); + if (otherTargetAtom != targetAtomB) + continue; + const int targetBondIdx = static_cast(targetTopology.bondIndices[adjIdx]); + if (targetBondIdx < 0 || targetBondIdx >= targetTopology.numBonds) { + continue; + } + if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) + continue; + outTargetBondIdx = targetBondIdx; + return true; + } + return false; + } else { + const bool scanAdjacency = targetTopology.rowOffsets != nullptr && targetTopology.colIndices != nullptr && + targetTopology.bondIndices != nullptr && targetAtomA >= 0 && + targetAtomA < targetTopology.numAtoms; + if (scanAdjacency) { + const int begin = static_cast(targetTopology.rowOffsets[targetAtomA]); + const int end = static_cast(targetTopology.rowOffsets[targetAtomA + 1]); + for (int adjIdx = begin; adjIdx < end; ++adjIdx) { + const int otherTargetAtom = static_cast(targetTopology.colIndices[adjIdx]); + if (otherTargetAtom != targetAtomB) + continue; + const int targetBondIdx = static_cast(targetTopology.bondIndices[adjIdx]); + if (targetBondIdx < 0 || targetBondIdx >= targetTopology.numBonds) { + continue; + } + if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) + continue; + outTargetBondIdx = targetBondIdx; + return true; + } + return false; + } + + for (int targetBondIdx = 0; targetBondIdx < targetTopology.numBonds; ++targetBondIdx) { + const std::uint32_t targetEndpoints = targetTopology.bondEndpoints[targetBondIdx]; + const int targetEndpointU = static_cast(targetEndpoints >> kBondEndpointShift); + const int targetEndpointV = static_cast(targetEndpoints & kBondEndpointMask); + const bool endpointsMatch = (targetEndpointU == targetAtomA && targetEndpointV == targetAtomB) || + (targetEndpointU == targetAtomB && targetEndpointV == targetAtomA); + if (!endpointsMatch) + continue; + if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) + continue; + outTargetBondIdx = targetBondIdx; + return true; + } + return false; + } +} + +template +__device__ __forceinline__ bool rebuildMatchFromSubstructureMappingWithinThread( + const Seed& seed, + const QueryTopology& queryTopology, + const TargetTopology& targetTopology, + const PairMatchTablesDevice& tables, + MatchResult& match, + FmcsSubstructureScratch& scratch) { + using SeedT = Seed; + using MatchT = MatchResult; + using BondWord = typename SeedT::bond_word_type; + using TargetBondWord = typename MatchT::target_bond_word; + using TargetAtomWord = typename MatchT::target_atom_word; + + constexpr int kBondBitsPerWord = SeedT::kBondBitsPerWord; + constexpr int kBondWords = SeedT::kBondWords; + constexpr int kTargetAtomBitsPerWord = MatchT::kTargetAtomBitsPerWord; + constexpr int kTargetBondBitsPerWord = MatchT::kTargetBondBitsPerWord; + + matchResultClearWithinThread(match); + for (int i = 0; i < seed.numAtoms; ++i) { + const int queryAtomIdx = scratch.seedAtomList[i]; + const int targetAtomIdx = scratch.targetAtomForQuery[queryAtomIdx]; + if (targetAtomIdx == kUnmappedTargetIdx) + return false; + match.targetAtomIdx[queryAtomIdx] = static_cast(targetAtomIdx); + match.visitedTargetAtoms[targetAtomIdx / kTargetAtomBitsPerWord] |= static_cast(1) + << (targetAtomIdx % kTargetAtomBitsPerWord); + } + match.matchedAtomSize = seed.numAtoms; + + int matchedBondCount = 0; + for (int wordIdx = 0; wordIdx < kBondWords; ++wordIdx) { + BondWord remaining = seed.bonds[wordIdx]; + while (remaining != 0) { + int bitPosInWord; + if constexpr (sizeof(BondWord) == 4) { + bitPosInWord = __ffs(static_cast(remaining)) - 1; + } else { + bitPosInWord = __ffsll(static_cast(remaining)) - 1; + } + const int queryBondIdx = wordIdx * kBondBitsPerWord + bitPosInWord; + remaining &= remaining - 1; + + const std::uint32_t queryEndpoints = queryTopology.bondEndpoints[queryBondIdx]; + const int queryEndpointU = static_cast(queryEndpoints >> kBondEndpointShift); + const int queryEndpointV = static_cast(queryEndpoints & kBondEndpointMask); + const int targetEndpointU = match.targetAtomIdx[queryEndpointU]; + const int targetEndpointV = match.targetAtomIdx[queryEndpointV]; + if (targetEndpointU == kUnmappedTargetIdx || targetEndpointV == kUnmappedTargetIdx) { + matchResultClearWithinThread(match); + return false; + } + + int targetBondIdx = -1; + if (!findTargetBondBetweenAtomsWithinThread(targetEndpointU, + targetEndpointV, + queryBondIdx, + targetTopology, + tables, + targetBondIdx)) { + matchResultClearWithinThread(match); + return false; + } + const TargetBondWord visitedWord = match.visitedTargetBonds[targetBondIdx / kTargetBondBitsPerWord]; + if ((visitedWord >> (targetBondIdx % kTargetBondBitsPerWord)) & 1) { + matchResultClearWithinThread(match); + return false; + } + match.targetBondIdx[queryBondIdx] = static_cast(targetBondIdx); + match.visitedTargetBonds[targetBondIdx / kTargetBondBitsPerWord] |= static_cast(1) + << (targetBondIdx % kTargetBondBitsPerWord); + ++matchedBondCount; + } + } + match.matchedBondSize = static_cast(matchedBondCount); + if (matchedBondCount != seed.numBonds) { + matchResultClearWithinThread(match); + return false; + } + match.empty = false; + return true; +} + +template +__device__ __forceinline__ void initializeSeedSubstructureScratchCooperative( + const GroupT& group, + const TargetTopology& targetTopology, + FmcsSubstructureScratch& scratch) { + const int laneRank = static_cast(group.thread_rank()); + const int laneCount = static_cast(group.num_threads()); + + for (int i = laneRank; i < maxAtoms; i += laneCount) { + scratch.seedDegree[i] = 0; + scratch.orderedQueryAtom[i] = 0; + scratch.queryOrderPos[i] = kUnmappedTargetIdx; + scratch.targetAtomForQuery[i] = kUnmappedTargetIdx; + } + + if constexpr (topologyHasAdjacencyBondIndices()) { + for (int targetAtomIdx = laneRank; targetAtomIdx < targetTopology.numAtoms; targetAtomIdx += laneCount) { + scratch.targetDegree[targetAtomIdx] = static_cast(targetTopology.rowOffsets[targetAtomIdx + 1] - + targetTopology.rowOffsets[targetAtomIdx]); + } + } else { + if (targetTopology.rowOffsets != nullptr) { + for (int targetAtomIdx = laneRank; targetAtomIdx < targetTopology.numAtoms; targetAtomIdx += laneCount) { + scratch.targetDegree[targetAtomIdx] = static_cast(targetTopology.rowOffsets[targetAtomIdx + 1] - + targetTopology.rowOffsets[targetAtomIdx]); + } + } else { + for (int i = laneRank; i < maxTA; i += laneCount) { + scratch.targetDegree[i] = 0; + } + group.sync(); + if (laneRank == 0) { + for (int targetBondIdx = 0; targetBondIdx < targetTopology.numBonds; ++targetBondIdx) { + const std::uint32_t targetEndpoints = targetTopology.bondEndpoints[targetBondIdx]; + const int targetEndpointU = static_cast(targetEndpoints >> kBondEndpointShift); + const int targetEndpointV = static_cast(targetEndpoints & kBondEndpointMask); + ++scratch.targetDegree[targetEndpointU]; + ++scratch.targetDegree[targetEndpointV]; + } + } + } + } + + if (laneRank == 0) { + scratch.currentCount = 0; + scratch.nextCount = 0; + scratch.found = 0; + scratch.overflowed = 0; + } + group.sync(); +} + +template +__device__ __forceinline__ bool prepareSeedSubstructureSearchWithinThread( + const Seed& seed, + const QueryTopology& queryTopology, + const TargetTopology& targetTopology, + const PairMatchTablesDevice& tables, + FmcsSubstructureScratch& scratch, + int& numSeedAtoms) { + using SeedT = Seed; + using AtomWord = typename SeedT::atom_word_type; + using BondWord = typename SeedT::bond_word_type; + + constexpr int kAtomBitsPerWord = SeedT::kAtomBitsPerWord; + constexpr int kAtomWords = SeedT::kAtomWords; + constexpr int kBondBitsPerWord = SeedT::kBondBitsPerWord; + constexpr int kBondWords = SeedT::kBondWords; + + if (seed.numAtoms > targetTopology.numAtoms || seed.numBonds > targetTopology.numBonds) { + return false; + } + + numSeedAtoms = 0; + for (int wordIdx = 0; wordIdx < kAtomWords; ++wordIdx) { + AtomWord remaining = seed.atoms[wordIdx]; + while (remaining != 0) { + int bitPosInWord; + if constexpr (sizeof(AtomWord) == 4) { + bitPosInWord = __ffs(static_cast(remaining)) - 1; + } else { + bitPosInWord = __ffsll(static_cast(remaining)) - 1; + } + const int queryAtomIdx = wordIdx * kAtomBitsPerWord + bitPosInWord; + remaining &= remaining - 1; + if (queryAtomIdx < queryTopology.numAtoms) { + scratch.seedAtomList[numSeedAtoms++] = static_cast(queryAtomIdx); + } + } + } + if (numSeedAtoms != seed.numAtoms) + return false; + + for (int wordIdx = 0; wordIdx < kBondWords; ++wordIdx) { + BondWord remaining = seed.bonds[wordIdx]; + while (remaining != 0) { + int bitPosInWord; + if constexpr (sizeof(BondWord) == 4) { + bitPosInWord = __ffs(static_cast(remaining)) - 1; + } else { + bitPosInWord = __ffsll(static_cast(remaining)) - 1; + } + const int queryBondIdx = wordIdx * kBondBitsPerWord + bitPosInWord; + remaining &= remaining - 1; + + const std::uint32_t queryEndpoints = queryTopology.bondEndpoints[queryBondIdx]; + const int queryEndpointU = static_cast(queryEndpoints >> kBondEndpointShift); + const int queryEndpointV = static_cast(queryEndpoints & kBondEndpointMask); + ++scratch.seedDegree[queryEndpointU]; + ++scratch.seedDegree[queryEndpointV]; + } + } + + for (int orderPos = 0; orderPos < numSeedAtoms; ++orderPos) { + int bestAtom = -1; + int bestMappedNeighborCount = -1; + int bestDegree = -1; + int bestCandidateCount = maxTA + 1; + + for (int atomListIdx = 0; atomListIdx < numSeedAtoms; ++atomListIdx) { + const int queryAtomIdx = scratch.seedAtomList[atomListIdx]; + if (scratch.orderedQueryAtom[queryAtomIdx]) + continue; + + int mappedNeighborCount = 0; + if (orderPos > 0) { + if constexpr (topologyHasAdjacencyBondIndices()) { + const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); + const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); + for (int adjIdx = begin; adjIdx < end; ++adjIdx) { + const int queryBondIdx = static_cast(queryTopology.bondIndices[adjIdx]); + if (queryBondIdx >= queryTopology.numBonds || + !seedContainsBondWithinThread(seed, queryBondIdx)) { + continue; + } + const int otherQueryAtom = static_cast(queryTopology.colIndices[adjIdx]); + if (otherQueryAtom >= 0 && otherQueryAtom < queryTopology.numAtoms && + scratch.orderedQueryAtom[otherQueryAtom]) { + ++mappedNeighborCount; + } + } + } else { + const bool scanQueryAdjacency = queryTopology.rowOffsets != nullptr && queryTopology.colIndices != nullptr && + queryTopology.bondIndices != nullptr && queryAtomIdx >= 0 && + queryAtomIdx < queryTopology.numAtoms; + if (scanQueryAdjacency) { + const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); + const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); + for (int adjIdx = begin; adjIdx < end; ++adjIdx) { + const int queryBondIdx = static_cast(queryTopology.bondIndices[adjIdx]); + if (queryBondIdx >= queryTopology.numBonds || + !seedContainsBondWithinThread(seed, queryBondIdx)) { + continue; + } + const int otherQueryAtom = static_cast(queryTopology.colIndices[adjIdx]); + if (otherQueryAtom >= 0 && otherQueryAtom < queryTopology.numAtoms && + scratch.orderedQueryAtom[otherQueryAtom]) { + ++mappedNeighborCount; + } + } + } else { + for (int wordIdx = 0; wordIdx < kBondWords; ++wordIdx) { + BondWord remaining = seed.bonds[wordIdx]; + while (remaining != 0) { + int bitPosInWord; + if constexpr (sizeof(BondWord) == 4) { + bitPosInWord = __ffs(static_cast(remaining)) - 1; + } else { + bitPosInWord = __ffsll(static_cast(remaining)) - 1; + } + const int queryBondIdx = wordIdx * kBondBitsPerWord + bitPosInWord; + remaining &= remaining - 1; + + const std::uint32_t queryEndpoints = queryTopology.bondEndpoints[queryBondIdx]; + const int queryEndpointU = static_cast(queryEndpoints >> kBondEndpointShift); + const int queryEndpointV = static_cast(queryEndpoints & kBondEndpointMask); + if (queryEndpointU == queryAtomIdx && scratch.orderedQueryAtom[queryEndpointV]) { + ++mappedNeighborCount; + } else if (queryEndpointV == queryAtomIdx && scratch.orderedQueryAtom[queryEndpointU]) { + ++mappedNeighborCount; + } + } + } + } + } + } + if (orderPos > 0 && mappedNeighborCount == 0) + continue; + + int candidateCount = 0; + for (int targetAtomIdx = 0; targetAtomIdx < targetTopology.numAtoms; ++targetAtomIdx) { + if (scratch.targetDegree[targetAtomIdx] < scratch.seedDegree[queryAtomIdx]) { + continue; + } + if (!tables.atoms.testBit(queryAtomIdx, targetAtomIdx)) + continue; + ++candidateCount; + } + if (candidateCount == 0) + return false; + + const int degree = scratch.seedDegree[queryAtomIdx]; + const bool better = bestAtom < 0 || mappedNeighborCount > bestMappedNeighborCount || + (mappedNeighborCount == bestMappedNeighborCount && degree > bestDegree) || + (mappedNeighborCount == bestMappedNeighborCount && degree == bestDegree && + candidateCount < bestCandidateCount) || + (mappedNeighborCount == bestMappedNeighborCount && degree == bestDegree && + candidateCount == bestCandidateCount && queryAtomIdx < bestAtom); + if (better) { + bestAtom = queryAtomIdx; + bestMappedNeighborCount = mappedNeighborCount; + bestDegree = degree; + bestCandidateCount = candidateCount; + } + } + + if (bestAtom < 0) { + for (int atomListIdx = 0; atomListIdx < numSeedAtoms; ++atomListIdx) { + const int queryAtomIdx = scratch.seedAtomList[atomListIdx]; + if (!scratch.orderedQueryAtom[queryAtomIdx]) { + bestAtom = queryAtomIdx; + break; + } + } + } + if (bestAtom < 0) + return false; + scratch.seedAtoms[orderPos] = static_cast(bestAtom); + scratch.orderedQueryAtom[bestAtom] = 1; + scratch.queryOrderPos[bestAtom] = static_cast(orderPos); + } + return true; +} + +__device__ __forceinline__ bool partialUsesTargetAtomWithinThread(const std::uint8_t* partial, + const int depth, + const int targetAtomIdx) { + for (int i = 0; i < depth; ++i) { + if (partial[i] == targetAtomIdx) + return true; + } + return false; +} + +template +__device__ __forceinline__ bool findMappedQueryNeighborWithinThread( + const Seed& seed, + const QueryTopology& queryTopology, + const FmcsSubstructureScratch& scratch, + const int depth, + const int queryAtomIdx, + int& outNeighborOrderPos) { + using SeedT = Seed; + using BondWord = typename SeedT::bond_word_type; + constexpr int kBondBitsPerWord = SeedT::kBondBitsPerWord; + constexpr int kBondWords = SeedT::kBondWords; + + outNeighborOrderPos = -1; + if constexpr (topologyHasAdjacencyBondIndices()) { + const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); + const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); + for (int adjIdx = begin; adjIdx < end; ++adjIdx) { + const int queryBondIdx = static_cast(queryTopology.bondIndices[adjIdx]); + if (queryBondIdx >= queryTopology.numBonds || + !seedContainsBondWithinThread(seed, queryBondIdx)) { + continue; + } + const int otherQueryAtom = static_cast(queryTopology.colIndices[adjIdx]); + if (otherQueryAtom < 0 || otherQueryAtom >= queryTopology.numAtoms) { + continue; + } + const int otherOrderPos = scratch.queryOrderPos[otherQueryAtom]; + if (otherOrderPos != kUnmappedTargetIdx && otherOrderPos < depth) { + outNeighborOrderPos = otherOrderPos; + return true; + } + } + return false; + } else { + const bool scanQueryAdjacency = queryTopology.rowOffsets != nullptr && queryTopology.colIndices != nullptr && + queryTopology.bondIndices != nullptr && queryAtomIdx >= 0 && + queryAtomIdx < queryTopology.numAtoms; + if (scanQueryAdjacency) { + const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); + const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); + for (int adjIdx = begin; adjIdx < end; ++adjIdx) { + const int queryBondIdx = static_cast(queryTopology.bondIndices[adjIdx]); + if (queryBondIdx >= queryTopology.numBonds || + !seedContainsBondWithinThread(seed, queryBondIdx)) { + continue; + } + const int otherQueryAtom = static_cast(queryTopology.colIndices[adjIdx]); + if (otherQueryAtom < 0 || otherQueryAtom >= queryTopology.numAtoms) { + continue; + } + const int otherOrderPos = scratch.queryOrderPos[otherQueryAtom]; + if (otherOrderPos != kUnmappedTargetIdx && otherOrderPos < depth) { + outNeighborOrderPos = otherOrderPos; + return true; + } + } + return false; + } + + for (int wordIdx = 0; wordIdx < kBondWords; ++wordIdx) { + BondWord remaining = seed.bonds[wordIdx]; + while (remaining != 0) { + int bitPosInWord; + if constexpr (sizeof(BondWord) == 4) { + bitPosInWord = __ffs(static_cast(remaining)) - 1; + } else { + bitPosInWord = __ffsll(static_cast(remaining)) - 1; + } + const int queryBondIdx = wordIdx * kBondBitsPerWord + bitPosInWord; + remaining &= remaining - 1; + + const std::uint32_t queryEndpoints = queryTopology.bondEndpoints[queryBondIdx]; + const int queryEndpointU = static_cast(queryEndpoints >> kBondEndpointShift); + const int queryEndpointV = static_cast(queryEndpoints & kBondEndpointMask); + int otherQueryAtom = -1; + if (queryEndpointU == queryAtomIdx) { + otherQueryAtom = queryEndpointV; + } else if (queryEndpointV == queryAtomIdx) { + otherQueryAtom = queryEndpointU; + } else { + continue; + } + + const int otherOrderPos = scratch.queryOrderPos[otherQueryAtom]; + if (otherOrderPos != kUnmappedTargetIdx && otherOrderPos < depth) { + outNeighborOrderPos = otherOrderPos; + return true; + } + } + } + return false; + } +} + +template +__device__ __forceinline__ bool substructurePartialEdgeConsistentWithinThread( + const Seed& seed, + const QueryTopology& queryTopology, + const TargetTopology& targetTopology, + const PairMatchTablesDevice& tables, + const FmcsSubstructureScratch& scratch, + const std::uint8_t* partial, + const int depth, + const int queryAtomIdx, + const int targetAtomIdx) { + using SeedT = Seed; + using BondWord = typename SeedT::bond_word_type; + constexpr int kBondBitsPerWord = SeedT::kBondBitsPerWord; + constexpr int kBondWords = SeedT::kBondWords; + + if constexpr (topologyHasAdjacencyBondIndices()) { + const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); + const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); + for (int adjIdx = begin; adjIdx < end; ++adjIdx) { + const int queryBondIdx = static_cast(queryTopology.bondIndices[adjIdx]); + if (queryBondIdx >= queryTopology.numBonds || + !seedContainsBondWithinThread(seed, queryBondIdx)) { + continue; + } + const int otherQueryAtom = static_cast(queryTopology.colIndices[adjIdx]); + if (otherQueryAtom < 0 || otherQueryAtom >= queryTopology.numAtoms) { + continue; + } + + const int otherOrderPos = scratch.queryOrderPos[otherQueryAtom]; + if (otherOrderPos == kUnmappedTargetIdx || otherOrderPos >= depth) { + continue; + } + const int otherTargetAtom = partial[otherOrderPos]; + + int targetBondIdx = -1; + if (!findTargetBondBetweenAtomsWithinThread(targetAtomIdx, + otherTargetAtom, + queryBondIdx, + targetTopology, + tables, + targetBondIdx)) { + return false; + } + } + return true; + } else { + const bool scanQueryAdjacency = queryTopology.rowOffsets != nullptr && queryTopology.colIndices != nullptr && + queryTopology.bondIndices != nullptr && queryAtomIdx >= 0 && + queryAtomIdx < queryTopology.numAtoms; + if (scanQueryAdjacency) { + const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); + const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); + for (int adjIdx = begin; adjIdx < end; ++adjIdx) { + const int queryBondIdx = static_cast(queryTopology.bondIndices[adjIdx]); + if (queryBondIdx >= queryTopology.numBonds || + !seedContainsBondWithinThread(seed, queryBondIdx)) { + continue; + } + const int otherQueryAtom = static_cast(queryTopology.colIndices[adjIdx]); + if (otherQueryAtom < 0 || otherQueryAtom >= queryTopology.numAtoms) { + continue; + } + + const int otherOrderPos = scratch.queryOrderPos[otherQueryAtom]; + if (otherOrderPos == kUnmappedTargetIdx || otherOrderPos >= depth) { + continue; + } + const int otherTargetAtom = partial[otherOrderPos]; + + int targetBondIdx = -1; + if (!findTargetBondBetweenAtomsWithinThread(targetAtomIdx, + otherTargetAtom, + queryBondIdx, + targetTopology, + tables, + targetBondIdx)) { + return false; + } + } + return true; + } + + for (int wordIdx = 0; wordIdx < kBondWords; ++wordIdx) { + BondWord remaining = seed.bonds[wordIdx]; + while (remaining != 0) { + int bitPosInWord; + if constexpr (sizeof(BondWord) == 4) { + bitPosInWord = __ffs(static_cast(remaining)) - 1; + } else { + bitPosInWord = __ffsll(static_cast(remaining)) - 1; + } + const int queryBondIdx = wordIdx * kBondBitsPerWord + bitPosInWord; + remaining &= remaining - 1; + + const std::uint32_t queryEndpoints = queryTopology.bondEndpoints[queryBondIdx]; + const int queryEndpointU = static_cast(queryEndpoints >> kBondEndpointShift); + const int queryEndpointV = static_cast(queryEndpoints & kBondEndpointMask); + int otherQueryAtom = -1; + if (queryEndpointU == queryAtomIdx) { + otherQueryAtom = queryEndpointV; + } else if (queryEndpointV == queryAtomIdx) { + otherQueryAtom = queryEndpointU; + } else { + continue; + } + + const int otherOrderPos = scratch.queryOrderPos[otherQueryAtom]; + if (otherOrderPos == kUnmappedTargetIdx || otherOrderPos >= depth) { + continue; + } + const int otherTargetAtom = partial[otherOrderPos]; + + int targetBondIdx = -1; + if (!findTargetBondBetweenAtomsWithinThread(targetAtomIdx, + otherTargetAtom, + queryBondIdx, + targetTopology, + tables, + targetBondIdx)) { + return false; + } + } + } + return true; + } +} + +template +__device__ __forceinline__ bool matchSeedSubstructureCooperative(const GroupT& group, + const Seed& seed, + const QueryTopology& queryTopology, + const TargetTopology& targetTopology, + const PairMatchTablesDevice& tables, + MatchResult& match, + FmcsSubstructureScratch& scratch, + std::uint8_t* partialStorage, + int partialCapacity, + bool* overflowedFlag) { + const int laneRank = static_cast(group.thread_rank()); + const int laneCount = static_cast(group.num_threads()); + + int numSeedAtoms = 0; + int prepared = 0; + if (laneRank == 0) { + matchResultClearWithinThread(match); + } + if (seed.numAtoms != 0) { + initializeSeedSubstructureScratchCooperative(group, targetTopology, scratch); + } + if (laneRank == 0) { + if (seed.numAtoms == 0) { + prepared = (seed.numBonds == 0) ? 2 : -1; + } else { + prepared = + prepareSeedSubstructureSearchWithinThread(seed, queryTopology, targetTopology, tables, scratch, numSeedAtoms) ? + 1 : + -1; + } + if ((partialStorage == nullptr || partialCapacity <= 0) && prepared == 1) { + scratch.overflowed = 1; + prepared = -1; + } + } + group.sync(); + prepared = group.shfl(prepared, 0); + numSeedAtoms = group.shfl(numSeedAtoms, 0); + if (prepared == 2) + return true; + if (prepared < 0) { + if (laneRank == 0 && scratch.overflowed && overflowedFlag != nullptr) { + *overflowedFlag = true; + } + return false; + } + + const int stride = numSeedAtoms; + const int halfBytes = partialCapacity * maxAtoms; + const int effectiveCapacity = halfBytes / stride; + std::uint8_t* currentPartials = partialStorage; + std::uint8_t* nextPartials = partialStorage + halfBytes; + + const int firstQueryAtom = scratch.seedAtoms[0]; + for (int targetAtomIdx = laneRank; targetAtomIdx < targetTopology.numAtoms; targetAtomIdx += laneCount) { + if (scratch.targetDegree[targetAtomIdx] < scratch.seedDegree[firstQueryAtom]) { + continue; + } + if (!tables.atoms.testBit(firstQueryAtom, targetAtomIdx)) + continue; + + const int slot = atomicAdd(&scratch.currentCount, 1); + if (slot < effectiveCapacity) { + currentPartials[slot * stride] = static_cast(targetAtomIdx); + } else { + atomicExch(&scratch.overflowed, 1); + } + } + group.sync(); + + if (numSeedAtoms == 1) { + if (scratch.currentCount > 0) { + if (laneRank == 0) { + scratch.targetAtomForQuery[firstQueryAtom] = currentPartials[0]; + prepared = + rebuildMatchFromSubstructureMappingWithinThread(seed, queryTopology, targetTopology, tables, match, scratch) ? + 1 : + -1; + } + group.sync(); + prepared = group.shfl(prepared, 0); + return prepared == 1; + } + return false; + } + + for (int depth = 1; depth < numSeedAtoms; ++depth) { + if (scratch.currentCount == 0 || scratch.found != 0) + break; + if (laneRank == 0) + scratch.nextCount = 0; + group.sync(); + + const int queryAtomIdx = scratch.seedAtoms[depth]; + const int numPartials = scratch.currentCount < effectiveCapacity ? scratch.currentCount : effectiveCapacity; + + for (int partialIdx = 0; partialIdx < numPartials && scratch.found == 0; ++partialIdx) { + const std::uint8_t* partial = currentPartials + partialIdx * stride; + int neighborOrderPos = -1; + const bool hasMappedNeighbor = + findMappedQueryNeighborWithinThread(seed, queryTopology, scratch, depth, queryAtomIdx, neighborOrderPos); + const bool scanAdjacency = + hasMappedNeighbor && targetTopology.rowOffsets != nullptr && targetTopology.colIndices != nullptr; + const int targetScanBegin = + scanAdjacency ? static_cast(targetTopology.rowOffsets[partial[neighborOrderPos]]) : 0; + const int targetScanEnd = scanAdjacency ? + static_cast(targetTopology.rowOffsets[partial[neighborOrderPos] + 1]) : + targetTopology.numAtoms; + + for (int targetScanIdx = targetScanBegin + laneRank; targetScanIdx < targetScanEnd && scratch.found == 0; + targetScanIdx += laneCount) { + const int targetAtomIdx = + scanAdjacency ? static_cast(targetTopology.colIndices[targetScanIdx]) : targetScanIdx; + if (scratch.targetDegree[targetAtomIdx] < scratch.seedDegree[queryAtomIdx]) { + continue; + } + if (!tables.atoms.testBit(queryAtomIdx, targetAtomIdx)) + continue; + if (partialUsesTargetAtomWithinThread(partial, depth, targetAtomIdx)) { + continue; + } + if (!substructurePartialEdgeConsistentWithinThread(seed, + queryTopology, + targetTopology, + tables, + scratch, + partial, + depth, + queryAtomIdx, + targetAtomIdx)) { + continue; + } + + if (depth == numSeedAtoms - 1) { + if (atomicCAS(&scratch.found, 0, 1) == 0) { + for (int orderPos = 0; orderPos < depth; ++orderPos) { + const int mappedQueryAtom = scratch.seedAtoms[orderPos]; + scratch.targetAtomForQuery[mappedQueryAtom] = partial[orderPos]; + } + scratch.targetAtomForQuery[queryAtomIdx] = static_cast(targetAtomIdx); + } + } else { + const int slot = atomicAdd(&scratch.nextCount, 1); + if (slot < effectiveCapacity) { + std::uint8_t* next = nextPartials + slot * stride; + for (int orderPos = 0; orderPos < depth; ++orderPos) { + next[orderPos] = partial[orderPos]; + } + next[depth] = static_cast(targetAtomIdx); + } else { + atomicExch(&scratch.overflowed, 1); + } + } + } + } + group.sync(); + + if (scratch.found != 0) + break; + if (laneRank == 0) + scratch.currentCount = scratch.nextCount; + std::uint8_t* tmp = currentPartials; + currentPartials = nextPartials; + nextPartials = tmp; + group.sync(); + } + + int found = scratch.found; + if (found != 0) { + if (laneRank == 0) { + found = + rebuildMatchFromSubstructureMappingWithinThread(seed, queryTopology, targetTopology, tables, match, scratch) ? + 1 : + 0; + } + group.sync(); + found = group.shfl(found, 0); + return found != 0; + } + + if (laneRank == 0) { + matchResultClearWithinThread(match); + if (scratch.overflowed && overflowedFlag != nullptr) { + *overflowedFlag = true; + } + } + group.sync(); + return false; +} + +template +__device__ __forceinline__ bool matchSeedWithSubstructureFallbackCooperative( + const GroupT& group, + const Seed& seed, + const QueryTopology& queryTopology, + const TargetTopology& targetTopology, + const PairMatchTablesDevice& tables, + MatchResult& match, + FmcsSubstructureScratch& scratch, + int* scratchLock, + std::uint8_t* partialStorage, + int partialCapacity, + bool* overflowedFlag) { + if (matchIncrementalFastCooperative(group, seed, queryTopology, targetTopology, tables, match)) { + return true; + } + if (group.thread_rank() == 0) { + while (atomicCAS(scratchLock, 0, 1) != 0) { + } + } + group.sync(); + const bool ok = matchSeedSubstructureCooperative(group, + seed, + queryTopology, + targetTopology, + tables, + match, + scratch, + partialStorage, + partialCapacity, + overflowedFlag); + group.sync(); + if (group.thread_rank() == 0) { + atomicExch(scratchLock, 0); + } + group.sync(); + return ok; +} + } // namespace fmcs } // namespace mcs diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f4647a0c..180bf3bd 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -86,6 +86,12 @@ target_include_directories(test_fmcs_match PRIVATE ${CMAKE_SOURCE_DIR}/src/mcs) target_link_libraries(test_fmcs_match PRIVATE mcs_fmcs_foundations device_vector) +add_executable(test_fmcs_substructure test_fmcs_substructure.cu) +target_include_directories(test_fmcs_substructure + PRIVATE ${CMAKE_SOURCE_DIR}/src/mcs) +target_link_libraries(test_fmcs_substructure + PRIVATE mcs_fmcs_foundations device_vector) + add_executable(test_openmp_helpers test_openmp_helpers.cpp) target_link_libraries(test_openmp_helpers PRIVATE openmp_helpers OpenMP::OpenMP_CXX) @@ -434,6 +440,7 @@ set(TEST_LIST test_fmcs_foundations test_fmcs_grow test_fmcs_match + test_fmcs_substructure test_work_splitting test_fire_minimizer test_fire_minimizer_permol) diff --git a/tests/test_fmcs_substructure.cu b/tests/test_fmcs_substructure.cu new file mode 100644 index 00000000..1e7925ae --- /dev/null +++ b/tests/test_fmcs_substructure.cu @@ -0,0 +1,582 @@ +// SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +// SPDX-License-Identifier: Apache-2.0 +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. +// +// Unit tests for the per-helper device functions in fmcs_cuda/. Each test +// launches a tiny __global__ driver that constructs inputs, calls one +// helper, and copies results back to host memory for assertion. This +// file is populated incrementally as Steps 1-5 land real implementations. + +#include +#include + +#include +#include +#include +#include + +#include "fmcs_cuda/fmcs_match.cuh" +#include "fmcs_cuda/fmcs_match_tables.cuh" +#include "fmcs_cuda/fmcs_seed.cuh" +#include "src/utils/device_vector.h" + +namespace { + +using nvMolKit::AsyncDevicePtr; +using nvMolKit::AsyncDeviceVector; + +using mcs::fmcs::MatchResult; +using mcs::fmcs::Seed; + +} // namespace + +// --------------------------------------------------------------------------- +// matchSingleBondWithinThread / matchIncrementalFastCooperative +// --------------------------------------------------------------------------- + +namespace { + +using mcs::fmcs::MatchTableDevice; +using mcs::fmcs::PairMatchTablesDevice; +using mcs::fmcs::SingleBondMatch; + +// Tiny CSR-view used by the match helpers via duck-typing; satisfies the +// QueryTopology / TargetTopology template requirement (bondEndpoints + +// numAtoms / numBonds), with optional CSR adjacency fields. +struct TestCsrView { + static constexpr bool kHasAdjacencyBondIndices = false; + + const std::uint32_t* bondEndpoints = nullptr; + int numAtoms = 0; + int numBonds = 0; + const std::uint32_t* rowOffsets = nullptr; + const std::uint32_t* colIndices = nullptr; + const std::uint32_t* bondIndices = nullptr; +}; + +// Match tables staged on the host and uploaded to device memory on +// demand. Both atom and bond tables are 32-bit row-packed bitmasks +// (one bit per (q, t) pair). Tests mutate the host staging vectors +// via the set*Bit helpers; device() uploads (if dirty) and returns the +// device view to pass to kernels. +struct ManagedMatchTables { + std::vector atomHost; + std::vector bondHost; + AsyncDeviceVector atomData; + AsyncDeviceVector bondData; + int qNumAtoms = 0; + int qNumBonds = 0; + + void allocate(int qAtoms, int tAtoms, int qBonds, int tBonds) { + qNumAtoms = qAtoms; + qNumBonds = qBonds; + const int atomWordsPerRow = (tAtoms + 31) / 32; + const int bondWordsPerRow = (tBonds + 31) / 32; + atomHost.assign(qAtoms * atomWordsPerRow, 0); + bondHost.assign(qBonds * bondWordsPerRow, 0); + dev_.atoms = MatchTableDevice{nullptr, qAtoms, tAtoms, atomWordsPerRow}; + dev_.bonds = MatchTableDevice{nullptr, qBonds, tBonds, bondWordsPerRow}; + dirty_ = true; + } + + void setAtomBit(int qAtom, int tAtom) { + atomHost[qAtom * dev_.atoms.wordsPerRow + tAtom / 32] |= (1u << (tAtom % 32)); + dirty_ = true; + } + void setBondBit(int qBond, int tBond) { + bondHost[qBond * dev_.bonds.wordsPerRow + tBond / 32] |= (1u << (tBond % 32)); + dirty_ = true; + } + void setAllAtomBits() { + for (int q = 0; q < dev_.atoms.nRows; ++q) + for (int t = 0; t < dev_.atoms.nCols; ++t) + setAtomBit(q, t); + } + void setAllBondBits() { + for (int q = 0; q < dev_.bonds.nRows; ++q) + for (int t = 0; t < dev_.bonds.nCols; ++t) + setBondBit(q, t); + } + + PairMatchTablesDevice device() { + if (dirty_) { + atomData.setFromVector(atomHost); + bondData.setFromVector(bondHost); + dev_.atoms.data = atomData.data(); + dev_.bonds.data = bondData.data(); + dirty_ = false; + } + return dev_; + } + + private: + PairMatchTablesDevice dev_{}; + bool dirty_ = true; +}; + +AsyncDeviceVector makeBondEndpointsDevice(const std::vector>& edges) { + std::vector host(edges.size()); + for (size_t i = 0; i < edges.size(); ++i) { + host[i] = (static_cast(edges[i].first) << 16) | static_cast(edges[i].second); + } + AsyncDeviceVector dev(edges.size()); + dev.copyFromHost(host); + return dev; +} + +} // namespace + +namespace { + +// Helper: build a parent MatchResult that records the mapping +// {qAtom[i] -> tAtom[i]} and {qBond[i] -> tBond[i]} from caller-supplied +// parallel arrays. Used by the incremental tests to set up "what the +// parent already had matched" before adding new bonds. +template +__device__ __forceinline__ void buildParentMatch(mcs::fmcs::MatchResult& match, + const int* qAtoms, + const int* tAtoms, + int nAtomMaps, + const int* qBonds, + const int* tBonds, + int nBondMaps) { + using MatchT = mcs::fmcs::MatchResult; + mcs::fmcs::matchResultClearWithinThread(match); + for (int i = 0; i < nAtomMaps; ++i) { + match.targetAtomIdx[qAtoms[i]] = static_cast(tAtoms[i]); + const int t = tAtoms[i]; + match.visitedTargetAtoms[t / MatchT::kTargetAtomBitsPerWord] |= (typename MatchT::target_atom_word{1}) + << (t % MatchT::kTargetAtomBitsPerWord); + } + for (int i = 0; i < nBondMaps; ++i) { + match.targetBondIdx[qBonds[i]] = static_cast(tBonds[i]); + const int t = tBonds[i]; + match.visitedTargetBonds[t / MatchT::kTargetBondBitsPerWord] |= (typename MatchT::target_bond_word{1}) + << (t % MatchT::kTargetBondBitsPerWord); + } + match.matchedAtomSize = static_cast(nAtomMaps); + match.matchedBondSize = static_cast(nBondMaps); + match.empty = (nAtomMaps == 0 && nBondMaps == 0); +} + +} // namespace +namespace mcs_fmcs_substructure_test { + +using QueuedT16 = mcs::fmcs::QueuedSeed<16, 16, 16, 16>; + +struct SubstructureTestOut { + bool ok; + bool overflowed; + QueuedT16 child; +}; + +constexpr int kTestSubstructurePartialCapacity = 64; + +__device__ __forceinline__ void addMaskSeed(QueuedT16& child, std::uint32_t atomMask, std::uint32_t bondMask) { + mcs::fmcs::seedClearWithinThread(child.seed); + mcs::fmcs::matchResultClearWithinThread(child.match); + for (int a = 0; a < 16; ++a) { + if ((atomMask >> a) & 1u) { + mcs::fmcs::seedAddAtomWithinThread(child.seed, a); + } + } + for (int b = 0; b < 16; ++b) { + if ((bondMask >> b) & 1u) { + mcs::fmcs::seedAddBondWithinThread(child.seed, b); + } + } +} + +__global__ void matchSubstructureMaskDriver(const std::uint32_t* qBE, + int qNumAtoms, + int qNumBonds, + const std::uint32_t* tBE, + int tNumAtoms, + int tNumBonds, + PairMatchTablesDevice tables, + std::uint32_t atomMask, + std::uint32_t bondMask, + std::uint8_t* partialStorage, + int partialCapacity, + SubstructureTestOut* out) { + __shared__ QueuedT16 child; + __shared__ mcs::fmcs::FmcsSubstructureScratch<16, 16> scratch; + if (threadIdx.x == 0) { + addMaskSeed(child, atomMask, bondMask); + } + __syncthreads(); + + auto block = cooperative_groups::this_thread_block(); + auto warp = cooperative_groups::tiled_partition<32>(block); + TestCsrView qView{qBE, qNumAtoms, qNumBonds}; + TestCsrView tView{tBE, tNumAtoms, tNumBonds}; + bool overflowed = false; + bool ok = mcs::fmcs::matchSeedSubstructureCooperative(warp, + child.seed, + qView, + tView, + tables, + child.match, + scratch, + partialStorage, + partialCapacity, + &overflowed); + __syncthreads(); + + if (threadIdx.x == 0) { + out->ok = ok; + out->overflowed = overflowed; + out->child = child; + } +} + +__global__ void matchFallbackBadParentDriver(const std::uint32_t* qBE, + int qNumAtoms, + int qNumBonds, + const std::uint32_t* tBE, + int tNumAtoms, + int tNumBonds, + PairMatchTablesDevice tables, + std::uint8_t* partialStorage, + int partialCapacity, + SubstructureTestOut* out) { + __shared__ QueuedT16 child; + __shared__ mcs::fmcs::FmcsSubstructureScratch<16, 16> scratch; + __shared__ int scratchLock; + if (threadIdx.x == 0) { + // Query seed is the 4-edge path inside the triangle-with-leaves + // repro: query bonds 1,2,3,4 and all five atoms. The stored parent + // match maps q bond 1 = (0,2) onto target bond 0 = (0,3), which is + // locally valid but blocks q bond 2 = (0,4) in the fast extender. + addMaskSeed(child, /*atomMask=*/0x1Fu, /*bondMask=*/0x1Eu); + mcs::fmcs::matchResultClearWithinThread(child.match); + using MatchT = decltype(child.match); + child.match.targetAtomIdx[0] = 0; + child.match.targetAtomIdx[2] = 3; + child.match.visitedTargetAtoms[0 / MatchT::kTargetAtomBitsPerWord] |= typename MatchT::target_atom_word{1} + << (0 % MatchT::kTargetAtomBitsPerWord); + child.match.visitedTargetAtoms[3 / MatchT::kTargetAtomBitsPerWord] |= typename MatchT::target_atom_word{1} + << (3 % MatchT::kTargetAtomBitsPerWord); + child.match.targetBondIdx[1] = 0; + child.match.visitedTargetBonds[0 / MatchT::kTargetBondBitsPerWord] |= typename MatchT::target_bond_word{1} + << (0 % MatchT::kTargetBondBitsPerWord); + child.match.matchedAtomSize = 2; + child.match.matchedBondSize = 1; + child.match.empty = false; + scratchLock = 0; + } + __syncthreads(); + + auto block = cooperative_groups::this_thread_block(); + auto warp = cooperative_groups::tiled_partition<32>(block); + TestCsrView qView{qBE, qNumAtoms, qNumBonds}; + TestCsrView tView{tBE, tNumAtoms, tNumBonds}; + bool overflowed = false; + bool ok = mcs::fmcs::matchSeedWithSubstructureFallbackCooperative(warp, + child.seed, + qView, + tView, + tables, + child.match, + scratch, + &scratchLock, + partialStorage, + partialCapacity, + &overflowed); + __syncthreads(); + + if (threadIdx.x == 0) { + out->ok = ok; + out->overflowed = overflowed; + out->child = child; + } +} + +} // namespace mcs_fmcs_substructure_test + +TEST(FMCSUnit, MatchSeedSubstructurePath) { + using namespace mcs_fmcs_substructure_test; + + auto qBE = makeBondEndpointsDevice({ + {0, 1}, + {1, 2}, + {2, 3} + }); + auto tBE = makeBondEndpointsDevice({ + {0, 1}, + {1, 2}, + {2, 3}, + {3, 4} + }); + ManagedMatchTables tables; + tables.allocate(4, 5, 3, 4); + tables.setAllAtomBits(); + tables.setAllBondBits(); + + AsyncDevicePtr d_out; + AsyncDeviceVector partials(2 * kTestSubstructurePartialCapacity * 16); + matchSubstructureMaskDriver<<<1, 32>>>(qBE.data(), + 4, + 3, + tBE.data(), + 5, + 4, + tables.device(), + /*atomMask=*/0xFu, + /*bondMask=*/0x7u, + partials.data(), + kTestSubstructurePartialCapacity, + d_out.data()); + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + SubstructureTestOut out{}; + d_out.get(out); + + EXPECT_TRUE(out.ok); + EXPECT_FALSE(out.overflowed); + EXPECT_EQ(out.child.match.matchedAtomSize, 4); + EXPECT_EQ(out.child.match.matchedBondSize, 3); + for (int q = 0; q < 4; ++q) { + EXPECT_NE(out.child.match.targetAtomIdx[q], mcs::fmcs::kUnmappedTargetIdx); + } + for (int q = 0; q < 3; ++q) { + EXPECT_NE(out.child.match.targetBondIdx[q], mcs::fmcs::kUnmappedTargetIdx); + } +} + +TEST(FMCSUnit, MatchSeedSubstructureRejectsNoMatch) { + using namespace mcs_fmcs_substructure_test; + + auto qBE = makeBondEndpointsDevice({ + {0, 1}, + {1, 2}, + {0, 2} + }); + auto tBE = makeBondEndpointsDevice({ + {0, 1}, + {1, 2} + }); + ManagedMatchTables tables; + tables.allocate(3, 3, 3, 2); + tables.setAllAtomBits(); + tables.setAllBondBits(); + + AsyncDevicePtr d_out; + AsyncDeviceVector partials(2 * kTestSubstructurePartialCapacity * 16); + matchSubstructureMaskDriver<<<1, 32>>>(qBE.data(), + 3, + 3, + tBE.data(), + 3, + 2, + tables.device(), + /*atomMask=*/0x7u, + /*bondMask=*/0x7u, + partials.data(), + kTestSubstructurePartialCapacity, + d_out.data()); + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + SubstructureTestOut out{}; + d_out.get(out); + + EXPECT_FALSE(out.ok); + EXPECT_FALSE(out.overflowed); + EXPECT_TRUE(out.child.match.empty); +} + +TEST(FMCSUnit, MatchSeedSubstructureRespectsAtomTable) { + using namespace mcs_fmcs_substructure_test; + + auto qBE = makeBondEndpointsDevice({ + {0, 1}, + {1, 2} + }); + auto tBE = makeBondEndpointsDevice({ + {0, 1}, + {1, 2} + }); + ManagedMatchTables tables; + tables.allocate(3, 3, 2, 2); + tables.setAtomBit(0, 0); + tables.setAtomBit(1, 1); + tables.setAtomBit(2, 2); + tables.setAllBondBits(); + + AsyncDevicePtr d_out; + AsyncDeviceVector partials(2 * kTestSubstructurePartialCapacity * 16); + matchSubstructureMaskDriver<<<1, 32>>>(qBE.data(), + 3, + 2, + tBE.data(), + 3, + 2, + tables.device(), + /*atomMask=*/0x7u, + /*bondMask=*/0x3u, + partials.data(), + kTestSubstructurePartialCapacity, + d_out.data()); + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + SubstructureTestOut out{}; + d_out.get(out); + + EXPECT_TRUE(out.ok); + EXPECT_FALSE(out.overflowed); + EXPECT_EQ(out.child.match.targetAtomIdx[0], 0u); + EXPECT_EQ(out.child.match.targetAtomIdx[1], 1u); + EXPECT_EQ(out.child.match.targetAtomIdx[2], 2u); + EXPECT_EQ(out.child.match.matchedAtomSize, 3); + EXPECT_EQ(out.child.match.matchedBondSize, 2); +} + +TEST(FMCSUnit, MatchSeedSubstructureRespectsBondTable) { + using namespace mcs_fmcs_substructure_test; + + auto qBE = makeBondEndpointsDevice({ + {0, 1}, + {1, 2} + }); + auto tBE = makeBondEndpointsDevice({ + {0, 1}, + {1, 2}, + {0, 2} + }); + ManagedMatchTables tables; + tables.allocate(3, 3, 2, 3); + tables.setAllAtomBits(); + tables.setBondBit(0, 0); + tables.setBondBit(1, 1); + + AsyncDevicePtr d_out; + AsyncDeviceVector partials(2 * kTestSubstructurePartialCapacity * 16); + matchSubstructureMaskDriver<<<1, 32>>>(qBE.data(), + 3, + 2, + tBE.data(), + 3, + 3, + tables.device(), + /*atomMask=*/0x7u, + /*bondMask=*/0x3u, + partials.data(), + kTestSubstructurePartialCapacity, + d_out.data()); + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + SubstructureTestOut out{}; + d_out.get(out); + + EXPECT_TRUE(out.ok); + EXPECT_FALSE(out.overflowed); + EXPECT_EQ(out.child.match.targetBondIdx[0], 0u); + EXPECT_EQ(out.child.match.targetBondIdx[1], 1u); + EXPECT_EQ(out.child.match.matchedAtomSize, 3); + EXPECT_EQ(out.child.match.matchedBondSize, 2); +} + +TEST(FMCSUnit, MatchSeedSubstructureFindsPathInsideTriangleWithLeaves) { + using namespace mcs_fmcs_substructure_test; + + auto qBE = makeBondEndpointsDevice({ + {0, 1}, + {0, 2}, + {0, 4}, + {1, 2}, + {1, 3} + }); + auto tBE = makeBondEndpointsDevice({ + {0, 3}, + {1, 2}, + {1, 4}, + {2, 3} + }); + ManagedMatchTables tables; + tables.allocate(5, 5, 5, 4); + tables.setAllAtomBits(); + tables.setAllBondBits(); + + AsyncDevicePtr d_out; + AsyncDeviceVector partials(2 * kTestSubstructurePartialCapacity * 16); + matchSubstructureMaskDriver<<<1, 32>>>(qBE.data(), + 5, + 5, + tBE.data(), + 5, + 4, + tables.device(), + /*atomMask=*/0x1Fu, + /*bondMask=*/0x1Eu, + partials.data(), + kTestSubstructurePartialCapacity, + d_out.data()); + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + SubstructureTestOut out{}; + d_out.get(out); + + EXPECT_TRUE(out.ok); + EXPECT_FALSE(out.overflowed); + EXPECT_EQ(out.child.match.matchedAtomSize, 5); + EXPECT_EQ(out.child.match.matchedBondSize, 4); + EXPECT_EQ(out.child.match.targetBondIdx[0], mcs::fmcs::kUnmappedTargetIdx); + for (int q : {1, 2, 3, 4}) { + EXPECT_NE(out.child.match.targetBondIdx[q], mcs::fmcs::kUnmappedTargetIdx); + } +} + +TEST(FMCSUnit, MatchSeedFallbackRebuildsAfterGreedyFailure) { + using namespace mcs_fmcs_substructure_test; + + auto qBE = makeBondEndpointsDevice({ + {0, 1}, + {0, 2}, + {0, 4}, + {1, 2}, + {1, 3} + }); + auto tBE = makeBondEndpointsDevice({ + {0, 3}, + {1, 2}, + {1, 4}, + {2, 3} + }); + ManagedMatchTables tables; + tables.allocate(5, 5, 5, 4); + tables.setAllAtomBits(); + tables.setAllBondBits(); + + AsyncDevicePtr d_out; + AsyncDeviceVector partials(2 * kTestSubstructurePartialCapacity * 16); + matchFallbackBadParentDriver<<<1, 32>>>(qBE.data(), + 5, + 5, + tBE.data(), + 5, + 4, + tables.device(), + partials.data(), + kTestSubstructurePartialCapacity, + d_out.data()); + ASSERT_EQ(cudaDeviceSynchronize(), cudaSuccess); + SubstructureTestOut out{}; + d_out.get(out); + + EXPECT_TRUE(out.ok); + EXPECT_FALSE(out.overflowed); + EXPECT_EQ(out.child.match.matchedAtomSize, 5); + EXPECT_EQ(out.child.match.matchedBondSize, 4); + for (int q = 0; q < 5; ++q) { + EXPECT_NE(out.child.match.targetAtomIdx[q], mcs::fmcs::kUnmappedTargetIdx); + } + for (int q : {1, 2, 3, 4}) { + EXPECT_NE(out.child.match.targetBondIdx[q], mcs::fmcs::kUnmappedTargetIdx); + } +} From 0a7916253c43856955e052ec458b0faeabc3de94 Mon Sep 17 00:00:00 2001 From: Kevin Boyd Date: Fri, 24 Jul 2026 13:47:48 -0400 Subject: [PATCH 6/8] Use project-rooted includes in fMCS substructure tests --- tests/test_fmcs_substructure.cu | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_fmcs_substructure.cu b/tests/test_fmcs_substructure.cu index 1e7925ae..cb612719 100644 --- a/tests/test_fmcs_substructure.cu +++ b/tests/test_fmcs_substructure.cu @@ -26,9 +26,9 @@ #include #include -#include "fmcs_cuda/fmcs_match.cuh" -#include "fmcs_cuda/fmcs_match_tables.cuh" -#include "fmcs_cuda/fmcs_seed.cuh" +#include "src/mcs/fmcs_cuda/fmcs_match.cuh" +#include "src/mcs/fmcs_cuda/fmcs_match_tables.cuh" +#include "src/mcs/fmcs_cuda/fmcs_seed.cuh" #include "src/utils/device_vector.h" namespace { From 5c31891c80351e1bb1015bd31bc482ceaef2dde0 Mon Sep 17 00:00:00 2001 From: Kevin Boyd Date: Tue, 28 Jul 2026 14:56:19 -0400 Subject: [PATCH 7/8] Drop the topology-generality concept from the substructure fallback Extends the cut in the parent commits to the exact substructure layer, which carried five more `if constexpr` sites on the same kHasAdjacencyBondIndices trait: findTargetBondBetweenAtomsWithinThread, the target-degree seeding, the mapped-neighbour count, the mapped-query- neighbour lookup, and the partial-mapping bond check. Every else branch was the same dead shape as before -- a copy of the CSR walk behind a runtime null-check, then a linear scan over all bonds that no production caller could reach. Only the test double set the trait false. Also renames the fast-path call to tryMatchIncrementalGreedyCooperative. The substructure fallback is exactly the caller the greedy contract describes: it runs the exact search when the greedy attempt returns false, rather than treating false as a rejection. test_fmcs_substructure.cu now builds real CSR via the same TestGraph helper as the match tests. fmcs_match.cuh drops another 431 lines. --- src/mcs/fmcs_cuda/fmcs_match.cuh | 431 +++++++------------------------ tests/test_fmcs_substructure.cu | 283 ++++++++++---------- 2 files changed, 241 insertions(+), 473 deletions(-) diff --git a/src/mcs/fmcs_cuda/fmcs_match.cuh b/src/mcs/fmcs_cuda/fmcs_match.cuh index a60c07d7..d58c4f1d 100644 --- a/src/mcs/fmcs_cuda/fmcs_match.cuh +++ b/src/mcs/fmcs_cuda/fmcs_match.cuh @@ -318,78 +318,38 @@ __device__ __forceinline__ bool tryMatchIncrementalGreedyCooperative( return true; } -template __device__ __forceinline__ bool findTargetBondBetweenAtomsWithinThread(const int targetAtomA, const int targetAtomB, const int queryBondIdx, - const TargetTopology& targetTopology, + const DeviceCsrView& targetTopology, const PairMatchTablesDevice& tables, int& outTargetBondIdx) { - if constexpr (topologyHasAdjacencyBondIndices()) { - if (targetAtomA < 0 || targetAtomA >= targetTopology.numAtoms) { - return false; - } - const int begin = static_cast(targetTopology.rowOffsets[targetAtomA]); - const int end = static_cast(targetTopology.rowOffsets[targetAtomA + 1]); - for (int adjIdx = begin; adjIdx < end; ++adjIdx) { - const int otherTargetAtom = static_cast(targetTopology.colIndices[adjIdx]); - if (otherTargetAtom != targetAtomB) - continue; - const int targetBondIdx = static_cast(targetTopology.bondIndices[adjIdx]); - if (targetBondIdx < 0 || targetBondIdx >= targetTopology.numBonds) { - continue; - } - if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) - continue; - outTargetBondIdx = targetBondIdx; - return true; - } + if (targetAtomA < 0 || targetAtomA >= targetTopology.numAtoms) { return false; - } else { - const bool scanAdjacency = targetTopology.rowOffsets != nullptr && targetTopology.colIndices != nullptr && - targetTopology.bondIndices != nullptr && targetAtomA >= 0 && - targetAtomA < targetTopology.numAtoms; - if (scanAdjacency) { - const int begin = static_cast(targetTopology.rowOffsets[targetAtomA]); - const int end = static_cast(targetTopology.rowOffsets[targetAtomA + 1]); - for (int adjIdx = begin; adjIdx < end; ++adjIdx) { - const int otherTargetAtom = static_cast(targetTopology.colIndices[adjIdx]); - if (otherTargetAtom != targetAtomB) - continue; - const int targetBondIdx = static_cast(targetTopology.bondIndices[adjIdx]); - if (targetBondIdx < 0 || targetBondIdx >= targetTopology.numBonds) { - continue; - } - if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) - continue; - outTargetBondIdx = targetBondIdx; - return true; - } - return false; - } - - for (int targetBondIdx = 0; targetBondIdx < targetTopology.numBonds; ++targetBondIdx) { - const std::uint32_t targetEndpoints = targetTopology.bondEndpoints[targetBondIdx]; - const int targetEndpointU = static_cast(targetEndpoints >> kBondEndpointShift); - const int targetEndpointV = static_cast(targetEndpoints & kBondEndpointMask); - const bool endpointsMatch = (targetEndpointU == targetAtomA && targetEndpointV == targetAtomB) || - (targetEndpointU == targetAtomB && targetEndpointV == targetAtomA); - if (!endpointsMatch) - continue; - if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) - continue; - outTargetBondIdx = targetBondIdx; - return true; + } + const int begin = static_cast(targetTopology.rowOffsets[targetAtomA]); + const int end = static_cast(targetTopology.rowOffsets[targetAtomA + 1]); + for (int adjIdx = begin; adjIdx < end; ++adjIdx) { + const int otherTargetAtom = static_cast(targetTopology.colIndices[adjIdx]); + if (otherTargetAtom != targetAtomB) + continue; + const int targetBondIdx = static_cast(targetTopology.bondIndices[adjIdx]); + if (targetBondIdx < 0 || targetBondIdx >= targetTopology.numBonds) { + continue; } - return false; + if (!tables.bonds.testBit(queryBondIdx, targetBondIdx)) + continue; + outTargetBondIdx = targetBondIdx; + return true; } + return false; } -template +template __device__ __forceinline__ bool rebuildMatchFromSubstructureMappingWithinThread( const Seed& seed, - const QueryTopology& queryTopology, - const TargetTopology& targetTopology, + const DeviceCsrView& queryTopology, + const DeviceCsrView& targetTopology, const PairMatchTablesDevice& tables, MatchResult& match, FmcsSubstructureScratch& scratch) { @@ -469,10 +429,10 @@ __device__ __forceinline__ bool rebuildMatchFromSubstructureMappingWithinThread( return true; } -template +template __device__ __forceinline__ void initializeSeedSubstructureScratchCooperative( const GroupT& group, - const TargetTopology& targetTopology, + const DeviceCsrView& targetTopology, FmcsSubstructureScratch& scratch) { const int laneRank = static_cast(group.thread_rank()); const int laneCount = static_cast(group.num_threads()); @@ -484,32 +444,9 @@ __device__ __forceinline__ void initializeSeedSubstructureScratchCooperative( scratch.targetAtomForQuery[i] = kUnmappedTargetIdx; } - if constexpr (topologyHasAdjacencyBondIndices()) { - for (int targetAtomIdx = laneRank; targetAtomIdx < targetTopology.numAtoms; targetAtomIdx += laneCount) { - scratch.targetDegree[targetAtomIdx] = static_cast(targetTopology.rowOffsets[targetAtomIdx + 1] - - targetTopology.rowOffsets[targetAtomIdx]); - } - } else { - if (targetTopology.rowOffsets != nullptr) { - for (int targetAtomIdx = laneRank; targetAtomIdx < targetTopology.numAtoms; targetAtomIdx += laneCount) { - scratch.targetDegree[targetAtomIdx] = static_cast(targetTopology.rowOffsets[targetAtomIdx + 1] - - targetTopology.rowOffsets[targetAtomIdx]); - } - } else { - for (int i = laneRank; i < maxTA; i += laneCount) { - scratch.targetDegree[i] = 0; - } - group.sync(); - if (laneRank == 0) { - for (int targetBondIdx = 0; targetBondIdx < targetTopology.numBonds; ++targetBondIdx) { - const std::uint32_t targetEndpoints = targetTopology.bondEndpoints[targetBondIdx]; - const int targetEndpointU = static_cast(targetEndpoints >> kBondEndpointShift); - const int targetEndpointV = static_cast(targetEndpoints & kBondEndpointMask); - ++scratch.targetDegree[targetEndpointU]; - ++scratch.targetDegree[targetEndpointV]; - } - } - } + for (int targetAtomIdx = laneRank; targetAtomIdx < targetTopology.numAtoms; targetAtomIdx += laneCount) { + scratch.targetDegree[targetAtomIdx] = static_cast(targetTopology.rowOffsets[targetAtomIdx + 1] - + targetTopology.rowOffsets[targetAtomIdx]); } if (laneRank == 0) { @@ -521,11 +458,11 @@ __device__ __forceinline__ void initializeSeedSubstructureScratchCooperative( group.sync(); } -template +template __device__ __forceinline__ bool prepareSeedSubstructureSearchWithinThread( const Seed& seed, - const QueryTopology& queryTopology, - const TargetTopology& targetTopology, + const DeviceCsrView& queryTopology, + const DeviceCsrView& targetTopology, const PairMatchTablesDevice& tables, FmcsSubstructureScratch& scratch, int& numSeedAtoms) { @@ -595,63 +532,18 @@ __device__ __forceinline__ bool prepareSeedSubstructureSearchWithinThread( int mappedNeighborCount = 0; if (orderPos > 0) { - if constexpr (topologyHasAdjacencyBondIndices()) { - const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); - const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); - for (int adjIdx = begin; adjIdx < end; ++adjIdx) { - const int queryBondIdx = static_cast(queryTopology.bondIndices[adjIdx]); - if (queryBondIdx >= queryTopology.numBonds || - !seedContainsBondWithinThread(seed, queryBondIdx)) { - continue; - } - const int otherQueryAtom = static_cast(queryTopology.colIndices[adjIdx]); - if (otherQueryAtom >= 0 && otherQueryAtom < queryTopology.numAtoms && - scratch.orderedQueryAtom[otherQueryAtom]) { - ++mappedNeighborCount; - } + const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); + const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); + for (int adjIdx = begin; adjIdx < end; ++adjIdx) { + const int queryBondIdx = static_cast(queryTopology.bondIndices[adjIdx]); + if (queryBondIdx >= queryTopology.numBonds || + !seedContainsBondWithinThread(seed, queryBondIdx)) { + continue; } - } else { - const bool scanQueryAdjacency = queryTopology.rowOffsets != nullptr && queryTopology.colIndices != nullptr && - queryTopology.bondIndices != nullptr && queryAtomIdx >= 0 && - queryAtomIdx < queryTopology.numAtoms; - if (scanQueryAdjacency) { - const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); - const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); - for (int adjIdx = begin; adjIdx < end; ++adjIdx) { - const int queryBondIdx = static_cast(queryTopology.bondIndices[adjIdx]); - if (queryBondIdx >= queryTopology.numBonds || - !seedContainsBondWithinThread(seed, queryBondIdx)) { - continue; - } - const int otherQueryAtom = static_cast(queryTopology.colIndices[adjIdx]); - if (otherQueryAtom >= 0 && otherQueryAtom < queryTopology.numAtoms && - scratch.orderedQueryAtom[otherQueryAtom]) { - ++mappedNeighborCount; - } - } - } else { - for (int wordIdx = 0; wordIdx < kBondWords; ++wordIdx) { - BondWord remaining = seed.bonds[wordIdx]; - while (remaining != 0) { - int bitPosInWord; - if constexpr (sizeof(BondWord) == 4) { - bitPosInWord = __ffs(static_cast(remaining)) - 1; - } else { - bitPosInWord = __ffsll(static_cast(remaining)) - 1; - } - const int queryBondIdx = wordIdx * kBondBitsPerWord + bitPosInWord; - remaining &= remaining - 1; - - const std::uint32_t queryEndpoints = queryTopology.bondEndpoints[queryBondIdx]; - const int queryEndpointU = static_cast(queryEndpoints >> kBondEndpointShift); - const int queryEndpointV = static_cast(queryEndpoints & kBondEndpointMask); - if (queryEndpointU == queryAtomIdx && scratch.orderedQueryAtom[queryEndpointV]) { - ++mappedNeighborCount; - } else if (queryEndpointV == queryAtomIdx && scratch.orderedQueryAtom[queryEndpointU]) { - ++mappedNeighborCount; - } - } - } + const int otherQueryAtom = static_cast(queryTopology.colIndices[adjIdx]); + if (otherQueryAtom >= 0 && otherQueryAtom < queryTopology.numAtoms && + scratch.orderedQueryAtom[otherQueryAtom]) { + ++mappedNeighborCount; } } } @@ -713,10 +605,10 @@ __device__ __forceinline__ bool partialUsesTargetAtomWithinThread(const std::uin return false; } -template +template __device__ __forceinline__ bool findMappedQueryNeighborWithinThread( const Seed& seed, - const QueryTopology& queryTopology, + const DeviceCsrView& queryTopology, const FmcsSubstructureScratch& scratch, const int depth, const int queryAtomIdx, @@ -727,92 +619,32 @@ __device__ __forceinline__ bool findMappedQueryNeighborWithinThread( constexpr int kBondWords = SeedT::kBondWords; outNeighborOrderPos = -1; - if constexpr (topologyHasAdjacencyBondIndices()) { - const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); - const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); - for (int adjIdx = begin; adjIdx < end; ++adjIdx) { - const int queryBondIdx = static_cast(queryTopology.bondIndices[adjIdx]); - if (queryBondIdx >= queryTopology.numBonds || - !seedContainsBondWithinThread(seed, queryBondIdx)) { - continue; - } - const int otherQueryAtom = static_cast(queryTopology.colIndices[adjIdx]); - if (otherQueryAtom < 0 || otherQueryAtom >= queryTopology.numAtoms) { - continue; - } - const int otherOrderPos = scratch.queryOrderPos[otherQueryAtom]; - if (otherOrderPos != kUnmappedTargetIdx && otherOrderPos < depth) { - outNeighborOrderPos = otherOrderPos; - return true; - } + const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); + const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); + for (int adjIdx = begin; adjIdx < end; ++adjIdx) { + const int queryBondIdx = static_cast(queryTopology.bondIndices[adjIdx]); + if (queryBondIdx >= queryTopology.numBonds || + !seedContainsBondWithinThread(seed, queryBondIdx)) { + continue; } - return false; - } else { - const bool scanQueryAdjacency = queryTopology.rowOffsets != nullptr && queryTopology.colIndices != nullptr && - queryTopology.bondIndices != nullptr && queryAtomIdx >= 0 && - queryAtomIdx < queryTopology.numAtoms; - if (scanQueryAdjacency) { - const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); - const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); - for (int adjIdx = begin; adjIdx < end; ++adjIdx) { - const int queryBondIdx = static_cast(queryTopology.bondIndices[adjIdx]); - if (queryBondIdx >= queryTopology.numBonds || - !seedContainsBondWithinThread(seed, queryBondIdx)) { - continue; - } - const int otherQueryAtom = static_cast(queryTopology.colIndices[adjIdx]); - if (otherQueryAtom < 0 || otherQueryAtom >= queryTopology.numAtoms) { - continue; - } - const int otherOrderPos = scratch.queryOrderPos[otherQueryAtom]; - if (otherOrderPos != kUnmappedTargetIdx && otherOrderPos < depth) { - outNeighborOrderPos = otherOrderPos; - return true; - } - } - return false; + const int otherQueryAtom = static_cast(queryTopology.colIndices[adjIdx]); + if (otherQueryAtom < 0 || otherQueryAtom >= queryTopology.numAtoms) { + continue; } - - for (int wordIdx = 0; wordIdx < kBondWords; ++wordIdx) { - BondWord remaining = seed.bonds[wordIdx]; - while (remaining != 0) { - int bitPosInWord; - if constexpr (sizeof(BondWord) == 4) { - bitPosInWord = __ffs(static_cast(remaining)) - 1; - } else { - bitPosInWord = __ffsll(static_cast(remaining)) - 1; - } - const int queryBondIdx = wordIdx * kBondBitsPerWord + bitPosInWord; - remaining &= remaining - 1; - - const std::uint32_t queryEndpoints = queryTopology.bondEndpoints[queryBondIdx]; - const int queryEndpointU = static_cast(queryEndpoints >> kBondEndpointShift); - const int queryEndpointV = static_cast(queryEndpoints & kBondEndpointMask); - int otherQueryAtom = -1; - if (queryEndpointU == queryAtomIdx) { - otherQueryAtom = queryEndpointV; - } else if (queryEndpointV == queryAtomIdx) { - otherQueryAtom = queryEndpointU; - } else { - continue; - } - - const int otherOrderPos = scratch.queryOrderPos[otherQueryAtom]; - if (otherOrderPos != kUnmappedTargetIdx && otherOrderPos < depth) { - outNeighborOrderPos = otherOrderPos; - return true; - } - } + const int otherOrderPos = scratch.queryOrderPos[otherQueryAtom]; + if (otherOrderPos != kUnmappedTargetIdx && otherOrderPos < depth) { + outNeighborOrderPos = otherOrderPos; + return true; } - return false; } + return false; } -template +template __device__ __forceinline__ bool substructurePartialEdgeConsistentWithinThread( const Seed& seed, - const QueryTopology& queryTopology, - const TargetTopology& targetTopology, + const DeviceCsrView& queryTopology, + const DeviceCsrView& targetTopology, const PairMatchTablesDevice& tables, const FmcsSubstructureScratch& scratch, const std::uint8_t* partial, @@ -824,124 +656,43 @@ __device__ __forceinline__ bool substructurePartialEdgeConsistentWithinThread( constexpr int kBondBitsPerWord = SeedT::kBondBitsPerWord; constexpr int kBondWords = SeedT::kBondWords; - if constexpr (topologyHasAdjacencyBondIndices()) { - const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); - const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); - for (int adjIdx = begin; adjIdx < end; ++adjIdx) { - const int queryBondIdx = static_cast(queryTopology.bondIndices[adjIdx]); - if (queryBondIdx >= queryTopology.numBonds || - !seedContainsBondWithinThread(seed, queryBondIdx)) { - continue; - } - const int otherQueryAtom = static_cast(queryTopology.colIndices[adjIdx]); - if (otherQueryAtom < 0 || otherQueryAtom >= queryTopology.numAtoms) { - continue; - } - - const int otherOrderPos = scratch.queryOrderPos[otherQueryAtom]; - if (otherOrderPos == kUnmappedTargetIdx || otherOrderPos >= depth) { - continue; - } - const int otherTargetAtom = partial[otherOrderPos]; - - int targetBondIdx = -1; - if (!findTargetBondBetweenAtomsWithinThread(targetAtomIdx, - otherTargetAtom, - queryBondIdx, - targetTopology, - tables, - targetBondIdx)) { - return false; - } + const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); + const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); + for (int adjIdx = begin; adjIdx < end; ++adjIdx) { + const int queryBondIdx = static_cast(queryTopology.bondIndices[adjIdx]); + if (queryBondIdx >= queryTopology.numBonds || + !seedContainsBondWithinThread(seed, queryBondIdx)) { + continue; } - return true; - } else { - const bool scanQueryAdjacency = queryTopology.rowOffsets != nullptr && queryTopology.colIndices != nullptr && - queryTopology.bondIndices != nullptr && queryAtomIdx >= 0 && - queryAtomIdx < queryTopology.numAtoms; - if (scanQueryAdjacency) { - const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); - const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); - for (int adjIdx = begin; adjIdx < end; ++adjIdx) { - const int queryBondIdx = static_cast(queryTopology.bondIndices[adjIdx]); - if (queryBondIdx >= queryTopology.numBonds || - !seedContainsBondWithinThread(seed, queryBondIdx)) { - continue; - } - const int otherQueryAtom = static_cast(queryTopology.colIndices[adjIdx]); - if (otherQueryAtom < 0 || otherQueryAtom >= queryTopology.numAtoms) { - continue; - } - - const int otherOrderPos = scratch.queryOrderPos[otherQueryAtom]; - if (otherOrderPos == kUnmappedTargetIdx || otherOrderPos >= depth) { - continue; - } - const int otherTargetAtom = partial[otherOrderPos]; - - int targetBondIdx = -1; - if (!findTargetBondBetweenAtomsWithinThread(targetAtomIdx, - otherTargetAtom, - queryBondIdx, - targetTopology, - tables, - targetBondIdx)) { - return false; - } - } - return true; + const int otherQueryAtom = static_cast(queryTopology.colIndices[adjIdx]); + if (otherQueryAtom < 0 || otherQueryAtom >= queryTopology.numAtoms) { + continue; } - for (int wordIdx = 0; wordIdx < kBondWords; ++wordIdx) { - BondWord remaining = seed.bonds[wordIdx]; - while (remaining != 0) { - int bitPosInWord; - if constexpr (sizeof(BondWord) == 4) { - bitPosInWord = __ffs(static_cast(remaining)) - 1; - } else { - bitPosInWord = __ffsll(static_cast(remaining)) - 1; - } - const int queryBondIdx = wordIdx * kBondBitsPerWord + bitPosInWord; - remaining &= remaining - 1; - - const std::uint32_t queryEndpoints = queryTopology.bondEndpoints[queryBondIdx]; - const int queryEndpointU = static_cast(queryEndpoints >> kBondEndpointShift); - const int queryEndpointV = static_cast(queryEndpoints & kBondEndpointMask); - int otherQueryAtom = -1; - if (queryEndpointU == queryAtomIdx) { - otherQueryAtom = queryEndpointV; - } else if (queryEndpointV == queryAtomIdx) { - otherQueryAtom = queryEndpointU; - } else { - continue; - } - - const int otherOrderPos = scratch.queryOrderPos[otherQueryAtom]; - if (otherOrderPos == kUnmappedTargetIdx || otherOrderPos >= depth) { - continue; - } - const int otherTargetAtom = partial[otherOrderPos]; - - int targetBondIdx = -1; - if (!findTargetBondBetweenAtomsWithinThread(targetAtomIdx, - otherTargetAtom, - queryBondIdx, - targetTopology, - tables, - targetBondIdx)) { - return false; - } - } + const int otherOrderPos = scratch.queryOrderPos[otherQueryAtom]; + if (otherOrderPos == kUnmappedTargetIdx || otherOrderPos >= depth) { + continue; + } + const int otherTargetAtom = partial[otherOrderPos]; + + int targetBondIdx = -1; + if (!findTargetBondBetweenAtomsWithinThread(targetAtomIdx, + otherTargetAtom, + queryBondIdx, + targetTopology, + tables, + targetBondIdx)) { + return false; } - return true; } + return true; } -template +template __device__ __forceinline__ bool matchSeedSubstructureCooperative(const GroupT& group, const Seed& seed, - const QueryTopology& queryTopology, - const TargetTopology& targetTopology, + const DeviceCsrView& queryTopology, + const DeviceCsrView& targetTopology, const PairMatchTablesDevice& tables, MatchResult& match, FmcsSubstructureScratch& scratch, @@ -1128,12 +879,12 @@ __device__ __forceinline__ bool matchSeedSubstructureCooperative(const GroupT& return false; } -template +template __device__ __forceinline__ bool matchSeedWithSubstructureFallbackCooperative( const GroupT& group, const Seed& seed, - const QueryTopology& queryTopology, - const TargetTopology& targetTopology, + const DeviceCsrView& queryTopology, + const DeviceCsrView& targetTopology, const PairMatchTablesDevice& tables, MatchResult& match, FmcsSubstructureScratch& scratch, @@ -1141,7 +892,7 @@ __device__ __forceinline__ bool matchSeedWithSubstructureFallbackCooperative( std::uint8_t* partialStorage, int partialCapacity, bool* overflowedFlag) { - if (matchIncrementalFastCooperative(group, seed, queryTopology, targetTopology, tables, match)) { + if (tryMatchIncrementalGreedyCooperative(group, seed, queryTopology, targetTopology, tables, match)) { return true; } if (group.thread_rank() == 0) { diff --git a/tests/test_fmcs_substructure.cu b/tests/test_fmcs_substructure.cu index cb612719..d50c768d 100644 --- a/tests/test_fmcs_substructure.cu +++ b/tests/test_fmcs_substructure.cu @@ -42,7 +42,7 @@ using mcs::fmcs::Seed; } // namespace // --------------------------------------------------------------------------- -// matchSingleBondWithinThread / matchIncrementalFastCooperative +// matchSingleBondWithinThread / tryMatchIncrementalGreedyCooperative // --------------------------------------------------------------------------- namespace { @@ -51,19 +51,7 @@ using mcs::fmcs::MatchTableDevice; using mcs::fmcs::PairMatchTablesDevice; using mcs::fmcs::SingleBondMatch; -// Tiny CSR-view used by the match helpers via duck-typing; satisfies the -// QueryTopology / TargetTopology template requirement (bondEndpoints + -// numAtoms / numBonds), with optional CSR adjacency fields. -struct TestCsrView { - static constexpr bool kHasAdjacencyBondIndices = false; - - const std::uint32_t* bondEndpoints = nullptr; - int numAtoms = 0; - int numBonds = 0; - const std::uint32_t* rowOffsets = nullptr; - const std::uint32_t* colIndices = nullptr; - const std::uint32_t* bondIndices = nullptr; -}; +using mcs::fmcs::DeviceCsrView; // Match tables staged on the host and uploaded to device memory on // demand. Both atom and bond tables are 32-bit row-packed bitmasks @@ -125,15 +113,68 @@ struct ManagedMatchTables { bool dirty_ = true; }; -AsyncDeviceVector makeBondEndpointsDevice(const std::vector>& edges) { - std::vector host(edges.size()); - for (size_t i = 0; i < edges.size(); ++i) { - host[i] = (static_cast(edges[i].first) << 16) | static_cast(edges[i].second); +// Owns the device buffers behind a DeviceCsrView. Built from an +// undirected edge list: bond i is edges[i], and the CSR is the symmetric +// expansion of that list (each edge contributes one entry to each +// endpoint's row), matching what the host-side graph builder uploads for +// real molecules. +class TestGraph { + public: + TestGraph(int numAtoms, const std::vector>& edges) { + const int numBonds = static_cast(edges.size()); + + std::vector bondEndpointsHost(numBonds); + for (int i = 0; i < numBonds; ++i) { + bondEndpointsHost[i] = + (static_cast(edges[i].first) << 16) | static_cast(edges[i].second); + } + + // Counting sort into CSR: degree pass, prefix sum, then scatter. + std::vector rowOffsetsHost(numAtoms + 1, 0); + for (const auto& edge : edges) { + ++rowOffsetsHost[edge.first + 1]; + ++rowOffsetsHost[edge.second + 1]; + } + for (int atom = 0; atom < numAtoms; ++atom) { + rowOffsetsHost[atom + 1] += rowOffsetsHost[atom]; + } + + std::vector colIndicesHost(2 * numBonds); + std::vector bondIndicesHost(2 * numBonds); + std::vector cursor(rowOffsetsHost.begin(), rowOffsetsHost.end() - 1); + for (int bond = 0; bond < numBonds; ++bond) { + const int u = edges[bond].first; + const int v = edges[bond].second; + colIndicesHost[cursor[u]] = static_cast(v); + bondIndicesHost[cursor[u]] = static_cast(bond); + ++cursor[u]; + colIndicesHost[cursor[v]] = static_cast(u); + bondIndicesHost[cursor[v]] = static_cast(bond); + ++cursor[v]; + } + + bondEndpoints_.setFromVector(bondEndpointsHost); + rowOffsets_.setFromVector(rowOffsetsHost); + colIndices_.setFromVector(colIndicesHost); + bondIndices_.setFromVector(bondIndicesHost); + + view_.bondEndpoints = bondEndpoints_.data(); + view_.rowOffsets = rowOffsets_.data(); + view_.colIndices = colIndices_.data(); + view_.bondIndices = bondIndices_.data(); + view_.numAtoms = numAtoms; + view_.numBonds = numBonds; } - AsyncDeviceVector dev(edges.size()); - dev.copyFromHost(host); - return dev; -} + + DeviceCsrView view() const { return view_; } + + private: + AsyncDeviceVector bondEndpoints_; + AsyncDeviceVector rowOffsets_; + AsyncDeviceVector colIndices_; + AsyncDeviceVector bondIndices_; + DeviceCsrView view_{}; +}; } // namespace @@ -198,12 +239,8 @@ __device__ __forceinline__ void addMaskSeed(QueuedT16& child, std::uint32_t atom } } -__global__ void matchSubstructureMaskDriver(const std::uint32_t* qBE, - int qNumAtoms, - int qNumBonds, - const std::uint32_t* tBE, - int tNumAtoms, - int tNumBonds, +__global__ void matchSubstructureMaskDriver(DeviceCsrView qView, + DeviceCsrView tView, PairMatchTablesDevice tables, std::uint32_t atomMask, std::uint32_t bondMask, @@ -217,12 +254,10 @@ __global__ void matchSubstructureMaskDriver(const std::uint32_t* qBE, } __syncthreads(); - auto block = cooperative_groups::this_thread_block(); - auto warp = cooperative_groups::tiled_partition<32>(block); - TestCsrView qView{qBE, qNumAtoms, qNumBonds}; - TestCsrView tView{tBE, tNumAtoms, tNumBonds}; - bool overflowed = false; - bool ok = mcs::fmcs::matchSeedSubstructureCooperative(warp, + auto block = cooperative_groups::this_thread_block(); + auto warp = cooperative_groups::tiled_partition<32>(block); + bool overflowed = false; + bool ok = mcs::fmcs::matchSeedSubstructureCooperative(warp, child.seed, qView, tView, @@ -241,12 +276,8 @@ __global__ void matchSubstructureMaskDriver(const std::uint32_t* qBE, } } -__global__ void matchFallbackBadParentDriver(const std::uint32_t* qBE, - int qNumAtoms, - int qNumBonds, - const std::uint32_t* tBE, - int tNumAtoms, - int tNumBonds, +__global__ void matchFallbackBadParentDriver(DeviceCsrView qView, + DeviceCsrView tView, PairMatchTablesDevice tables, std::uint8_t* partialStorage, int partialCapacity, @@ -278,12 +309,10 @@ __global__ void matchFallbackBadParentDriver(const std::uint32_t* qBE, } __syncthreads(); - auto block = cooperative_groups::this_thread_block(); - auto warp = cooperative_groups::tiled_partition<32>(block); - TestCsrView qView{qBE, qNumAtoms, qNumBonds}; - TestCsrView tView{tBE, tNumAtoms, tNumBonds}; - bool overflowed = false; - bool ok = mcs::fmcs::matchSeedWithSubstructureFallbackCooperative(warp, + auto block = cooperative_groups::this_thread_block(); + auto warp = cooperative_groups::tiled_partition<32>(block); + bool overflowed = false; + bool ok = mcs::fmcs::matchSeedWithSubstructureFallbackCooperative(warp, child.seed, qView, tView, @@ -308,16 +337,18 @@ __global__ void matchFallbackBadParentDriver(const std::uint32_t* qBE, TEST(FMCSUnit, MatchSeedSubstructurePath) { using namespace mcs_fmcs_substructure_test; - auto qBE = makeBondEndpointsDevice({ - {0, 1}, - {1, 2}, - {2, 3} + TestGraph query(4, + { + {0, 1}, + {1, 2}, + {2, 3} }); - auto tBE = makeBondEndpointsDevice({ - {0, 1}, - {1, 2}, - {2, 3}, - {3, 4} + TestGraph target(5, + { + {0, 1}, + {1, 2}, + {2, 3}, + {3, 4} }); ManagedMatchTables tables; tables.allocate(4, 5, 3, 4); @@ -326,12 +357,8 @@ TEST(FMCSUnit, MatchSeedSubstructurePath) { AsyncDevicePtr d_out; AsyncDeviceVector partials(2 * kTestSubstructurePartialCapacity * 16); - matchSubstructureMaskDriver<<<1, 32>>>(qBE.data(), - 4, - 3, - tBE.data(), - 5, - 4, + matchSubstructureMaskDriver<<<1, 32>>>(query.view(), + target.view(), tables.device(), /*atomMask=*/0xFu, /*bondMask=*/0x7u, @@ -357,14 +384,16 @@ TEST(FMCSUnit, MatchSeedSubstructurePath) { TEST(FMCSUnit, MatchSeedSubstructureRejectsNoMatch) { using namespace mcs_fmcs_substructure_test; - auto qBE = makeBondEndpointsDevice({ - {0, 1}, - {1, 2}, - {0, 2} + TestGraph query(3, + { + {0, 1}, + {1, 2}, + {0, 2} }); - auto tBE = makeBondEndpointsDevice({ - {0, 1}, - {1, 2} + TestGraph target(3, + { + {0, 1}, + {1, 2} }); ManagedMatchTables tables; tables.allocate(3, 3, 3, 2); @@ -373,12 +402,8 @@ TEST(FMCSUnit, MatchSeedSubstructureRejectsNoMatch) { AsyncDevicePtr d_out; AsyncDeviceVector partials(2 * kTestSubstructurePartialCapacity * 16); - matchSubstructureMaskDriver<<<1, 32>>>(qBE.data(), - 3, - 3, - tBE.data(), - 3, - 2, + matchSubstructureMaskDriver<<<1, 32>>>(query.view(), + target.view(), tables.device(), /*atomMask=*/0x7u, /*bondMask=*/0x7u, @@ -397,13 +422,15 @@ TEST(FMCSUnit, MatchSeedSubstructureRejectsNoMatch) { TEST(FMCSUnit, MatchSeedSubstructureRespectsAtomTable) { using namespace mcs_fmcs_substructure_test; - auto qBE = makeBondEndpointsDevice({ - {0, 1}, - {1, 2} + TestGraph query(3, + { + {0, 1}, + {1, 2} }); - auto tBE = makeBondEndpointsDevice({ - {0, 1}, - {1, 2} + TestGraph target(3, + { + {0, 1}, + {1, 2} }); ManagedMatchTables tables; tables.allocate(3, 3, 2, 2); @@ -414,12 +441,8 @@ TEST(FMCSUnit, MatchSeedSubstructureRespectsAtomTable) { AsyncDevicePtr d_out; AsyncDeviceVector partials(2 * kTestSubstructurePartialCapacity * 16); - matchSubstructureMaskDriver<<<1, 32>>>(qBE.data(), - 3, - 2, - tBE.data(), - 3, - 2, + matchSubstructureMaskDriver<<<1, 32>>>(query.view(), + target.view(), tables.device(), /*atomMask=*/0x7u, /*bondMask=*/0x3u, @@ -442,14 +465,16 @@ TEST(FMCSUnit, MatchSeedSubstructureRespectsAtomTable) { TEST(FMCSUnit, MatchSeedSubstructureRespectsBondTable) { using namespace mcs_fmcs_substructure_test; - auto qBE = makeBondEndpointsDevice({ - {0, 1}, - {1, 2} + TestGraph query(3, + { + {0, 1}, + {1, 2} }); - auto tBE = makeBondEndpointsDevice({ - {0, 1}, - {1, 2}, - {0, 2} + TestGraph target(3, + { + {0, 1}, + {1, 2}, + {0, 2} }); ManagedMatchTables tables; tables.allocate(3, 3, 2, 3); @@ -459,12 +484,8 @@ TEST(FMCSUnit, MatchSeedSubstructureRespectsBondTable) { AsyncDevicePtr d_out; AsyncDeviceVector partials(2 * kTestSubstructurePartialCapacity * 16); - matchSubstructureMaskDriver<<<1, 32>>>(qBE.data(), - 3, - 2, - tBE.data(), - 3, - 3, + matchSubstructureMaskDriver<<<1, 32>>>(query.view(), + target.view(), tables.device(), /*atomMask=*/0x7u, /*bondMask=*/0x3u, @@ -486,18 +507,20 @@ TEST(FMCSUnit, MatchSeedSubstructureRespectsBondTable) { TEST(FMCSUnit, MatchSeedSubstructureFindsPathInsideTriangleWithLeaves) { using namespace mcs_fmcs_substructure_test; - auto qBE = makeBondEndpointsDevice({ - {0, 1}, - {0, 2}, - {0, 4}, - {1, 2}, - {1, 3} + TestGraph query(5, + { + {0, 1}, + {0, 2}, + {0, 4}, + {1, 2}, + {1, 3} }); - auto tBE = makeBondEndpointsDevice({ - {0, 3}, - {1, 2}, - {1, 4}, - {2, 3} + TestGraph target(5, + { + {0, 3}, + {1, 2}, + {1, 4}, + {2, 3} }); ManagedMatchTables tables; tables.allocate(5, 5, 5, 4); @@ -506,12 +529,8 @@ TEST(FMCSUnit, MatchSeedSubstructureFindsPathInsideTriangleWithLeaves) { AsyncDevicePtr d_out; AsyncDeviceVector partials(2 * kTestSubstructurePartialCapacity * 16); - matchSubstructureMaskDriver<<<1, 32>>>(qBE.data(), - 5, - 5, - tBE.data(), - 5, - 4, + matchSubstructureMaskDriver<<<1, 32>>>(query.view(), + target.view(), tables.device(), /*atomMask=*/0x1Fu, /*bondMask=*/0x1Eu, @@ -535,18 +554,20 @@ TEST(FMCSUnit, MatchSeedSubstructureFindsPathInsideTriangleWithLeaves) { TEST(FMCSUnit, MatchSeedFallbackRebuildsAfterGreedyFailure) { using namespace mcs_fmcs_substructure_test; - auto qBE = makeBondEndpointsDevice({ - {0, 1}, - {0, 2}, - {0, 4}, - {1, 2}, - {1, 3} + TestGraph query(5, + { + {0, 1}, + {0, 2}, + {0, 4}, + {1, 2}, + {1, 3} }); - auto tBE = makeBondEndpointsDevice({ - {0, 3}, - {1, 2}, - {1, 4}, - {2, 3} + TestGraph target(5, + { + {0, 3}, + {1, 2}, + {1, 4}, + {2, 3} }); ManagedMatchTables tables; tables.allocate(5, 5, 5, 4); @@ -555,12 +576,8 @@ TEST(FMCSUnit, MatchSeedFallbackRebuildsAfterGreedyFailure) { AsyncDevicePtr d_out; AsyncDeviceVector partials(2 * kTestSubstructurePartialCapacity * 16); - matchFallbackBadParentDriver<<<1, 32>>>(qBE.data(), - 5, - 5, - tBE.data(), - 5, - 4, + matchFallbackBadParentDriver<<<1, 32>>>(query.view(), + target.view(), tables.device(), partials.data(), kTestSubstructurePartialCapacity, From a1f14966ebef16d5869da091d27a8ca37ffb3a6e Mon Sep 17 00:00:00 2001 From: Kevin Boyd Date: Fri, 31 Jul 2026 10:51:51 -0400 Subject: [PATCH 8/8] Formatting --- src/mcs/fmcs_cuda/fmcs_match.cuh | 13 ------------- tests/CMakeLists.txt | 4 ++-- 2 files changed, 2 insertions(+), 15 deletions(-) diff --git a/src/mcs/fmcs_cuda/fmcs_match.cuh b/src/mcs/fmcs_cuda/fmcs_match.cuh index 9ca6eca8..be6f679a 100644 --- a/src/mcs/fmcs_cuda/fmcs_match.cuh +++ b/src/mcs/fmcs_cuda/fmcs_match.cuh @@ -318,7 +318,6 @@ __device__ __forceinline__ bool tryMatchIncrementalGreedyCooperative( return true; } -<<<<<<< HEAD __device__ __forceinline__ bool findTargetBondBetweenAtomsWithinThread(const int targetAtomA, const int targetAtomB, const int queryBondIdx, @@ -614,11 +613,6 @@ __device__ __forceinline__ bool findMappedQueryNeighborWithinThread( const int depth, const int queryAtomIdx, int& outNeighborOrderPos) { - using SeedT = Seed; - using BondWord = typename SeedT::bond_word_type; - constexpr int kBondBitsPerWord = SeedT::kBondBitsPerWord; - constexpr int kBondWords = SeedT::kBondWords; - outNeighborOrderPos = -1; const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); @@ -652,11 +646,6 @@ __device__ __forceinline__ bool substructurePartialEdgeConsistentWithinThread( const int depth, const int queryAtomIdx, const int targetAtomIdx) { - using SeedT = Seed; - using BondWord = typename SeedT::bond_word_type; - constexpr int kBondBitsPerWord = SeedT::kBondBitsPerWord; - constexpr int kBondWords = SeedT::kBondWords; - const int begin = static_cast(queryTopology.rowOffsets[queryAtomIdx]); const int end = static_cast(queryTopology.rowOffsets[queryAtomIdx + 1]); for (int adjIdx = begin; adjIdx < end; ++adjIdx) { @@ -919,8 +908,6 @@ __device__ __forceinline__ bool matchSeedWithSubstructureFallbackCooperative( return ok; } -======= ->>>>>>> main } // namespace fmcs } // namespace mcs diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index bf7b291f..0445863a 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -92,8 +92,8 @@ target_link_libraries(test_fmcs_match PRIVATE mcs_fmcs_foundations add_executable(test_fmcs_substructure test_fmcs_substructure.cu) target_include_directories(test_fmcs_substructure PRIVATE ${CMAKE_SOURCE_DIR}/src/mcs) -target_link_libraries(test_fmcs_substructure - PRIVATE mcs_fmcs_foundations device_vector) +target_link_libraries(test_fmcs_substructure PRIVATE mcs_fmcs_foundations + device_vector) add_executable(test_openmp_helpers test_openmp_helpers.cpp) target_link_libraries(test_openmp_helpers PRIVATE openmp_helpers