Merge pull request #437 from seoklab/fmt/pdb/no-guess

feat(fmt/pdb): don't guess in reader

include/nuri/algo/guess.h

@@ -6,6 +6,8 @@
 #ifndef NURI_ALGO_GUESS_H_
 #define NURI_ALGO_GUESS_H_
 
+#include <absl/base/attributes.h>
+
 #include "nuri/core/molecule.h"
 
 namespace nuri {
@@ -26,8 +28,9 @@ constexpr double kDefaultThreshold = 0.5;
  * If connectivity information is already present and is correct, consider using
  * guess_all_types().
  */
-extern bool guess_everything(MoleculeMutator &mut, int conf = 0,
-                             double threshold = kDefaultThreshold);
+ABSL_MUST_USE_RESULT extern bool
+guess_everything(MoleculeMutator &mut, int conf = 0,
+                 double threshold = kDefaultThreshold);
 
 /**
  * @brief Guess connectivity information of a molecule.
@@ -53,7 +56,7 @@ extern void guess_connectivity(MoleculeMutator &mut, int conf = 0,
  * and all atom/bond types and implicit hydrogen counts are incorrect. The
  * information present in the molecule could be overwritten by this function.
  */
-extern bool guess_all_types(Molecule &mol, int conf = 0);
+ABSL_MUST_USE_RESULT extern bool guess_all_types(Molecule &mol, int conf = 0);
 
 /**
  * @brief Guess formal charges of a molecule.
@@ -84,6 +87,10 @@ extern void guess_hydrogens_2d(Molecule &mol);
  * guess_hydrogens() in sequence.
  */
 extern void guess_fcharge_hydrogens_2d(Molecule &mol);
+
+namespace internal {
+  ABSL_MUST_USE_RESULT extern bool guess_update_subs(Molecule &mol);
+}  // namespace internal
 }  // namespace nuri
 
 #endif /* NURI_ALGO_GUESS_H_ */

include/nuri/fmt/base.h

@@ -129,7 +129,7 @@ class MoleculeReader {
    */
   virtual Molecule parse(const std::vector<std::string> &block) const = 0;
 
-  virtual bool sanitized() const = 0;
+  virtual bool bond_valid() const = 0;
 
   MoleculeStream<MoleculeReader> stream() { return { *this }; }
 };

include/nuri/fmt/mol2.h

@@ -31,7 +31,7 @@ class Mol2Reader: public DefaultReaderImpl<read_mol2> {
 
   bool getnext(std::vector<std::string> &block) override;
 
-  bool sanitized() const override { return false; }
+  bool bond_valid() const override { return true; }
 
 private:
   bool read_mol_header_ = false;

include/nuri/fmt/pdb.h

@@ -31,7 +31,7 @@ class PDBReader: public DefaultReaderImpl<read_pdb> {
 
   bool getnext(std::vector<std::string> &block) override;
 
-  bool sanitized() const override { return true; }
+  bool bond_valid() const override { return false; }
 
 private:
   std::vector<std::string> header_;

include/nuri/fmt/sdf.h

@@ -31,7 +31,7 @@ class SDFReader: public DefaultReaderImpl<read_sdf> {
 
   bool getnext(std::vector<std::string> &block) override;
 
-  bool sanitized() const override { return false; }
+  bool bond_valid() const override { return true; }
 };
 
 class SDFReaderFactory: public DefaultReaderFactoryImpl<SDFReader> {

include/nuri/fmt/smiles.h

@@ -33,7 +33,7 @@ class SmilesReader: public DefaultReaderImpl<read_smiles> {
 
   bool getnext(std::vector<std::string> &block) override;
 
-  bool sanitized() const override { return false; }
+  bool bond_valid() const override { return true; }
 };
 
 class SmilesReaderFactory: public DefaultReaderFactoryImpl<SmilesReader> {

python/docs/conf.py.in

@@ -8,9 +8,9 @@
 # For the full list of built-in configuration values, see the documentation:
 # https://www.sphinx-doc.org/en/master/usage/configuration.html
 
+import doctest
 import re
 
-import doctest
 from sphinx.ext import autodoc
 
 # -- Project information -----------------------------------------------------
@@ -61,10 +61,10 @@ html_static_path = ["@CMAKE_CURRENT_LIST_DIR@/_static"]
 if "@DOXYGEN_OUTPUT_DIR@":
     html_extra_path = ["@DOXYGEN_OUTPUT_DIR@/html"]
 html_css_files = [
-    "css/rtd-property.css" # workaround https://github.com/readthedocs/sphinx_rtd_theme/issues/1301
+    "css/rtd-property.css"  # workaround https://github.com/readthedocs/sphinx_rtd_theme/issues/1301
 ]
 html_theme_options = {
-    'navigation_depth': -1,
+    "navigation_depth": -1,
 }
 
 intersphinx_mapping = {
@@ -83,7 +83,8 @@ doctest_global_setup = "import nuri"
 class UnneededBaseStripDocumenter(autodoc.ClassDocumenter):
     _strip_re = re.compile(
         r"\s*Bases:\s*:py:class:"
-        r"`([^`]*?\bpybind11_[A-Za-z0-9_\.]+|object)`\s*")
+        r"`([^`]*?\bpybind11_[A-Za-z0-9_\.]+|object)`\s*"
+    )
 
     def add_line(self, line: str, source: str, *lineno: int):
         if self._strip_re.fullmatch(line):

python/src/nuri/fmt/fmt_module.cpp

@@ -22,6 +22,7 @@
 #include <pybind11/pytypes.h>
 #include <pybind11/stl/filesystem.h>
 
+#include "nuri/algo/guess.h"
 #include "nuri/core/molecule.h"
 #include "nuri/fmt/base.h"
 #include "nuri/fmt/mol2.h"
@@ -38,7 +39,8 @@ class PyMoleculeReader {
 public:
   PyMoleculeReader(std::unique_ptr<std::istream> is, std::string_view fmt,
                    bool sanitize, bool skip_on_error)
-      : stream_(std::move(is)), skip_on_error_(skip_on_error) {
+      : stream_(std::move(is)), sanitize_(sanitize),
+        skip_on_error_(skip_on_error) {
     if (!*stream_)
       throw py::value_error(absl::StrCat("Invalid stream object"));
 
@@ -51,7 +53,7 @@ class PyMoleculeReader {
     if (!reader_)
       throw py::value_error(absl::StrCat("Failed to create reader for ", fmt));
 
-    sanitize_ = sanitize && !reader_->sanitized();
+    guess_ = sanitize && !reader_->bond_valid();
   }
 
   auto next() {
@@ -78,7 +80,16 @@ class PyMoleculeReader {
         continue;
       }
 
-      if (sanitize_ && !MoleculeSanitizer(mol).sanitize_all()) {
+      if (guess_ && mol.is_3d()) {
+        if (!internal::guess_update_subs(mol)) {
+          log_or_throw("Failed to guess molecule atom/bond types");
+          continue;
+        }
+      } else if (sanitize_ && !MoleculeSanitizer(mol).sanitize_all()) {
+        ABSL_LOG_IF(WARNING, guess_ && !mol.is_3d())
+            << "Reader might produce molecules with invalid bonds, but the "
+               "molecule is missing 3D coordinates; guessing is disabled.";
+
         log_or_throw("Failed to sanitize molecule");
         continue;
       }
@@ -102,6 +113,7 @@ class PyMoleculeReader {
   std::vector<std::string> block_;
   bool sanitize_;
   bool skip_on_error_;
+  bool guess_;
 };
 
 template <class F, class... Args>
@@ -146,9 +158,10 @@ Read a molecule from a file.
 
 :param fmt: The format of the file.
 :param path: The path to the file.
-:param sanitize: Whether to sanitize the produced molecule. Note that if the
-  underlying reader produces a sanitized molecule, this option is ignored and
-  the molecule is always sanitized.
+:param sanitize: Whether to sanitize the produced molecule. For formats that is
+  known to produce molecules with insufficient bond information (e.g. PDB), this
+  option will trigger guessing based on the 3D coordinates
+  (:func:`nuri.algo.guess_everything()`).
 :param skip_on_error: Whether to skip a molecule if an error occurs, instead of
   raising an exception.
 :raises OSError: If any file-related error occurs.
@@ -171,9 +184,10 @@ Read a molecule from string.
 
 :param fmt: The format of the file.
 :param data: The string to read.
-:param sanitize: Whether to sanitize the produced molecule. Note that if the
-  underlying reader produces a sanitized molecule, this option is ignored and
-  the molecule is always sanitized.
+:param sanitize: Whether to sanitize the produced molecule. For formats that is
+  known to produce molecules with insufficient bond information (e.g. PDB), this
+  option will trigger guessing based on the 3D coordinates
+  (:func:`nuri.algo.guess_everything()`).
 :param skip_on_error: Whether to skip a molecule if an error occurs, instead of
   raising an exception.
 :raises ValueError: If the format is unknown or sanitization fails, unless

src/algo/guess/3d.cpp

@@ -2170,4 +2170,19 @@ bool guess_all_types(Molecule &mol, int conf) {
   reset_bonds(mol);
   return guess_types_common(mol, mol.confs()[conf]);
 }
+
+namespace internal {
+  bool guess_update_subs(Molecule &mol) {
+    bool ret;
+    {
+      auto mut = mol.mutator();
+      ret = guess_everything(mut);
+    }
+
+    for (auto &sub: mol.substructures())
+      sub.refresh_bonds();
+
+    return ret;
+  }
+}  // namespace internal
 }  // namespace nuri