Explicit hydrogen removal parity #67

Author: LouieSlocombe

README.md

@@ -171,7 +171,7 @@ supplied in full. Options may appear before or after the input and use
 | `--accept-palindromes=<0\|1>` | `0` | Identify a string fragment with its reversal. |
 | `--parallel=<auto\|on\|off>` | `off` | Select parallel search automatically, require it, or disable it. |
 | `--threads=<auto\|N>` | `auto` | Set the OpenMP thread count per process; `N` must be positive. |
-| `--remove-hydrogens=<0\|1>` | `1` | Remove explicit hydrogens from molfiles. |
+| `--remove-hydrogens=<0\|1>` | `1` | Remove explicit hydrogens from MOL/SDF and native graph inputs. |
 | `--verbose=<0\|1>` | `0` | Print the parsed input graph. |
 | `--compensate-disjoint=<0\|1>` | `0` | Subtract one per processed component after the first. |
 | `--memory-report=<0\|1>` | `0` | Write Linux peak virtual memory to `memUsage`. |

unitTests/graphioTester.cpp

@@ -120,6 +120,7 @@ namespace
 
 int main()
 {
+    removeHydrogens = true;
     verbose = false;
 
     molGraph graph = parse(validGraph);
@@ -135,6 +136,39 @@ int main()
     assert(graph.atoms[1].bonds[0].bondType == 1);
     assert(graph.atoms[1].bonds[1].bondType == 2);
 
+    const string explicitHydrogenGraph =
+        "explicit hydrogens\n"
+        "4\n"
+        "1 3 2 3 3 4\n"
+        "H H C C\n"
+        "1 1 1\n";
+    const molGraph hydrogensRemoved = parse(explicitHydrogenGraph);
+    assert(hydrogensRemoved.atoms.size() == 2);
+    assert(hydrogensRemoved.totalBonds == 1);
+    assert(hydrogensRemoved.atoms[0].atomType == "C");
+    assert(hydrogensRemoved.atoms[1].atomType == "C");
+    assert(hydrogensRemoved.atoms[0].bonds[0].neighbourAtomIndex == 1);
+
+    removeHydrogens = false;
+    const molGraph hydrogensRetained = parse(explicitHydrogenGraph);
+    assert(hydrogensRetained.atoms.size() == 4);
+    assert(hydrogensRetained.totalBonds == 3);
+    assert(hydrogensRetained.atoms[0].atomType == "H");
+    assert(hydrogensRetained.atoms[1].atomType == "H");
+
+    removeHydrogens = true;
+    const molGraph unrestrictedLabels = parse(
+        "unrestricted labels\n"
+        "2\n"
+        "1 2\n"
+        "COLLAPSE He\n"
+        "1\n"
+    );
+    assert(unrestrictedLabels.atoms.size() == 2);
+    assert(unrestrictedLabels.totalBonds == 1);
+    assert(unrestrictedLabels.atoms[0].atomType == "COLLAPSE");
+    assert(unrestrictedLabels.atoms[1].atomType == "He");
+
     molGraph replacementTarget = sentinelGraph();
     istringstream validInput(validGraph);
     graphio(validInput, replacementTarget);

unitTests/libraryBatchTester.cpp

@@ -91,6 +91,30 @@ int main(int argc, char **argv)
         )
     ) return 1;
 
+    const std::string explicitHydrogenGraph =
+        "explicit hydrogens\n"
+        "6\n"
+        "1 3 2 3 3 4 4 5 4 6\n"
+        "H H C C H H\n"
+        "1 1 1 1 1\n";
+    std::istringstream filteredGraphStream(explicitHydrogenGraph);
+    const assemblycpp::CalculationResult filteredGraph =
+        assemblycpp::calculateGraph(filteredGraphStream);
+    assemblycpp::CalculationOptions retainedHydrogenOptions;
+    retainedHydrogenOptions.removeHydrogens = false;
+    std::istringstream retainedGraphStream(explicitHydrogenGraph);
+    const assemblycpp::CalculationResult retainedGraph =
+        assemblycpp::calculateGraph(
+            retainedGraphStream,
+            retainedHydrogenOptions
+        );
+    if (
+        !require(filteredGraph.succeeded, "filtered graph calculation failed") ||
+        !require(filteredGraph.assemblyIndex == 0, "filtered graph index mismatch") ||
+        !require(retainedGraph.succeeded, "retained graph calculation failed") ||
+        !require(retainedGraph.assemblyIndex == 3, "retained graph index mismatch")
+    ) return 1;
+
     std::istringstream invalidGraph(
         "invalid graph\n2\n1 3\nC C\n1\n"
     );

unitTests/unitTester.py

@@ -1275,6 +1275,7 @@ def run_cli_checks(executable: Path) -> int:
                 5,
             ),
         )
+        molecular_hydrogen_graphs = {}
         for (
             name,
             hydrogen_options,
@@ -1309,20 +1310,43 @@ def run_cli_checks(executable: Path) -> int:
                 "atoms or bonds",
                 completed,
             )
+            molecular_hydrogen_graphs[tuple(hydrogen_options)] = graph
             scenarios += 1
 
-        native_hydrogen_graph = "native-hydrogens\n4\n1 3 2 3 3 4\nH H C C\n1 1 1\n"
-        expected_native_graph = None
-        for name, hydrogen_option in (
-            ("native-hydrogens-on", "--remove-hydrogens=1"),
-            ("native-hydrogens-off", "--remove-hydrogens=0"),
-        ):
+        native_hydrogen_graph = (
+            "native-hydrogens\n"
+            "6\n"
+            "1 3 2 3 3 4 4 5 4 6\n"
+            "H H C C H H\n"
+            "1 1 1 1 1\n"
+        )
+        native_hydrogen_cases = (
+            ("native-hydrogens-default", [], ["C", "C"], 1),
+            (
+                "native-hydrogens-on",
+                ["--remove-hydrogens=1"],
+                ["C", "C"],
+                1,
+            ),
+            (
+                "native-hydrogens-off",
+                ["--remove-hydrogens=0"],
+                ["H", "H", "C", "C", "H", "H"],
+                5,
+            ),
+        )
+        for (
+            name,
+            hydrogen_options,
+            expected_colours,
+            expected_edge_count,
+        ) in native_hydrogen_cases:
             case_directory = working_directory / name
             case_directory.mkdir()
             (case_directory / "input").write_text(native_hydrogen_graph)
             completed = run_cli_command(
                 executable,
-                ["input", "--pathway=1", hydrogen_option],
+                ["input", "--pathway=1", *hydrogen_options],
                 case_directory,
             )
             require_cli(
@@ -1338,20 +1362,18 @@ def run_cli_checks(executable: Path) -> int:
             )
             graph = parse_pathway_document(pathway_path)["file_graph"][0]
             require_cli(
-                graph["VertexColours"] == ["H", "H", "C", "C"]
-                and len(graph["Edges"]) == 3,
-                "--remove-hydrogens should not transform native graph inputs",
+                graph["VertexColours"] == expected_colours
+                and len(graph["Edges"]) == expected_edge_count,
+                f"native-graph hydrogen scenario {name!r} transformed the "
+                "wrong atoms or bonds",
+                completed,
+            )
+            require_cli(
+                graph == molecular_hydrogen_graphs[tuple(hydrogen_options)],
+                f"native-graph hydrogen scenario {name!r} should match its "
+                "MOL equivalent",
                 completed,
             )
-            if expected_native_graph is None:
-                expected_native_graph = graph
-            else:
-                require_cli(
-                    graph == expected_native_graph,
-                    "native graph output should be identical with hydrogen "
-                    "removal on or off",
-                    completed,
-                )
             scenarios += 1
 
         all_hydrogen_mol = (

v5/assemblycpp.h

@@ -24,6 +24,7 @@ struct CalculationOptions
 {
     std::uint64_t runtimeTicks = std::numeric_limits<std::uint64_t>::max();
     int enumerationLimit = 50000000;
+    /** Remove explicit H vertices from MOL/SDF and native graph inputs. */
     bool removeHydrogens = true;
     bool compensateDisjoint = false;
     bool verbose = false;

v5/graphio.h

@@ -236,6 +236,7 @@ inline void graphio(std::istream &inputStream, molGraph &molecule)
         );
     }
 
+    if (removeHydrogens) parsed.removeExplicitHydrogens();
     if (verbose)
     {
         const std::vector<std::string> nameFields = graphioDetail::tokens(nameLine);

v5/ioflag.h

@@ -109,7 +109,7 @@ const vector<InputFlagDefinition>& inputFlagDefinitions()
             "remove-hydrogens",
             "0|1",
             "1",
-            "Remove explicit hydrogens from molfile inputs.",
+            "Remove explicit hydrogens from MOL/SDF and native graph inputs.",
             {"removeHydrogens"}
         },
         {

v5/molGraph.h

@@ -143,17 +143,18 @@ struct molGraph
 
 private:
     /**
-     * @brief Rebuild the graph while dropping removed atoms and zero-order bonds.
+     * @brief Rebuild the graph while dropping selected atoms and zero-order bonds.
      */
-    void rebuild(bool removeMarkedAtoms)
+    template <typename AtomRemovalPredicate>
+    void rebuild(AtomRemovalPredicate shouldRemoveAtom)
     {
         const size_t removed = atoms.size();
         vector<size_t> reverseMap(atoms.size(), removed);
         molGraph output;
 
         for (size_t i = 0; i < atoms.size(); i++)
         {
-            if (removeMarkedAtoms && atoms[i].atomType == "COLLAPSE") continue;
+            if (shouldRemoveAtom(atoms[i])) continue;
             reverseMap[i] = output.atoms.size();
             output.addAtom(atoms[i].atomType);
         }
@@ -190,26 +191,36 @@ struct molGraph
      */
     void collapse()
     {
-        rebuild(false);
+        rebuild([](const atom &) { return false; });
     }
 
-    /**
-     * @brief For explicit hydrogen removal
-     *
-     */
+    /** Mark one atom for removal by removeAndCollapse(). */
     void removeAtom(size_t i)
     {
         if (i >= atoms.size()) return;
         atoms[i].atomType = "COLLAPSE";
     }
 
-    /**
-     * @brief For explicit hydrogen removal
-     *
-     */
+    /** Remove atoms marked by removeAtom(), then rebuild the graph. */
     void removeAndCollapse()
     {
-        rebuild(true);
+        rebuild(
+            [](const atom &candidate)
+            {
+                return candidate.atomType == "COLLAPSE";
+            }
+        );
+    }
+
+    /** Remove every explicitly represented hydrogen atom and its bonds. */
+    void removeExplicitHydrogens()
+    {
+        rebuild(
+            [](const atom &candidate)
+            {
+                return candidate.atomType == "H";
+            }
+        );
     }
 
     /**

v5/molfileParser.h

@@ -166,14 +166,7 @@ void molfileParser(std::istream &molfile, molGraph &molecule)
             );
         parsed.addBond(atomA - 1, atomB - 1, bondOrder);
     }
-    if (removeHydrogens)
-    {
-        for (size_t atomIndex = 0; atomIndex < parsed.atoms.size(); atomIndex++)
-        {
-            if (parsed.atomType(atomIndex) == "H") parsed.removeAtom(atomIndex);
-        }
-        parsed.removeAndCollapse();
-    }
+    if (removeHydrogens) parsed.removeExplicitHydrogens();
     if (verbose) parsed.printToCout();
     molecule = std::move(parsed);
 }