Merge pull request #38 from RMeli/develop

Increase robustness and number of tests

.appveyor.yml

@@ -53,7 +53,7 @@ build: false
 
 test_script:
     # Run test and get coverage report
-  - pytest -v --cov=spyrmsd tests --runslow
+  - pytest -v --cov=spyrmsd tests --benchmark --large
 
     # Run checks
   - flake8

.gitignore

@@ -70,6 +70,10 @@ ENV/
 # mypy
 .mypy_cache/
 
+# tests
+tests/docking.zip
+tests/docking
+
 # duecredit
 .duecredit.p
 

.travis.yml

@@ -120,7 +120,7 @@ install:
 
 script:
     # Run test and get coverage report
-  - pytest -v --cov=spyrmsd tests/ --runslow
+  - pytest -v --cov=spyrmsd tests/  --benchmark --large
 
     # Run checks
   - flake8

README.md

@@ -11,7 +11,7 @@
 [![Documentation Status](https://readthedocs.org/projects/spyrmsd/badge/?version=develop)](https://spyrmsd.readthedocs.io/en/develop/?badge=develop)
 
 [![License](https://img.shields.io/github/license/RMeli/pyrmsd?color=%2333BBFF)](https://opensource.org/licenses/MIT)
-[![PyPI](https://img.shields.io/badge/PyPI-v0.3.5%20-ff69b4)](https://pypi.org/project/spyrmsd/)
+[![PyPI](https://img.shields.io/badge/PyPI-v0.4.0%20-ff69b4)](https://pypi.org/project/spyrmsd/)
 [![Conda Version](https://img.shields.io/conda/vn/conda-forge/spyrmsd.svg)](https://anaconda.org/conda-forge/spyrmsd)
 
 [![J. Cheminform.](https://img.shields.io/badge/J.%20Cheminform.-10.1186%2Fs13321--020--00455--2-blue)](https://doi.org/10.1186/s13321-020-00455-2)

devtools/README.md

@@ -1,6 +1,6 @@
 # Development, testing, and deployment tools
 
-This directory contains a collection of tools for running Continuous Integration (CI) tests, 
+This directory contains a collection of tools for running Continuous Integration (CI) tests,
 conda installation, and other development tools not directly related to the coding process.
 
 ## Continuous Integration

docs/source/api/spyrmsd.graphs.rst

@@ -10,6 +10,7 @@ Submodules
 ----------
 
 .. toctree::
+   :maxdepth: 4
 
    spyrmsd.graphs.gt
    spyrmsd.graphs.nx

docs/source/api/spyrmsd.optional.rst

@@ -10,6 +10,7 @@ Submodules
 ----------
 
 .. toctree::
+   :maxdepth: 4
 
    spyrmsd.optional.obabel
    spyrmsd.optional.rdkit

docs/source/api/spyrmsd.rst

@@ -10,6 +10,7 @@ Subpackages
 -----------
 
 .. toctree::
+   :maxdepth: 4
 
    spyrmsd.graphs
    spyrmsd.optional
@@ -18,6 +19,7 @@ Submodules
 ----------
 
 .. toctree::
+   :maxdepth: 4
 
    spyrmsd.graph
    spyrmsd.hungarian

spyrmsd/io.py

@@ -47,12 +47,8 @@
         "bonds",
     ]
 
-from typing import List
 
-from spyrmsd import molecule
-
-
-def loadmol(fname: str, adjacency: bool = True) -> molecule.Molecule:
+def loadmol(fname: str, adjacency: bool = True):
     """
     Load molecule from file.
 
@@ -72,7 +68,7 @@ def loadmol(fname: str, adjacency: bool = True) -> molecule.Molecule:
     return to_molecule(mol, adjacency=adjacency)
 
 
-def loadallmols(fname: str, adjacency: bool = True) -> List[molecule.Molecule]:
+def loadallmols(fname: str, adjacency: bool = True):
     """
     Load molecules from file.
 

spyrmsd/optional/obabel.py

@@ -94,7 +94,7 @@ def adjacency_matrix(mol) -> np.ndarray:
     return A
 
 
-def to_molecule(mol, adjacency: bool = True) -> molecule.Molecule:
+def to_molecule(mol, adjacency: bool = True):
     """
     Transform molecule to `pyrmsd` molecule.
 

spyrmsd/optional/rdkit.py

@@ -81,7 +81,7 @@ def adjacency_matrix(mol) -> np.ndarray:
     return Chem.rdmolops.GetAdjacencyMatrix(mol)
 
 
-def to_molecule(mol, adjacency: bool = True) -> molecule.Molecule:
+def to_molecule(mol, adjacency: bool = True):
     """
     Transform molecule to `pyrmsd` molecule.
 
@@ -98,6 +98,10 @@ def to_molecule(mol, adjacency: bool = True) -> molecule.Molecule:
         `spyrmsd` molecule
     """
 
+    if mol is None:
+        # Propagate RDKit parsing failure
+        return None
+
     atoms = mol.GetAtoms()
 
     n = len(atoms)

spyrmsd/qcp.py

@@ -196,12 +196,33 @@ def dP(x):
         """
         return 4 * x ** 3 + 2 * c2 * x + c1
 
-    x0 = (Ga + Gb) / 2.0
+    x0 = (Ga + Gb) * 0.5
 
-    return optimize.newton(P, x0, fprime=dP)
+    lmax = optimize.newton(P, x0, fprime=dP)
 
+    return lmax
 
-def qcp_rmsd(A: np.ndarray, B: np.ndarray) -> float:
+
+def _lambda_max_eig(K: np.ndarray) -> float:
+    """
+    Find largest eigenvalue of :math:`K`.
+
+    Parameters
+    ----------
+    K: np.ndarray
+        Symmetric key matrix
+
+    Returns
+    -------
+    float
+        Largest eigenvalue of :math:`K`, :math:`\\lambda_\\text{max}`
+    """
+    e, _ = np.linalg.eig(K)
+
+    return max(e)
+
+
+def qcp_rmsd(A: np.ndarray, B: np.ndarray, atol: float = 1e-9) -> float:
     """
     Compute RMSD using the quaternion polynomial method.
 
@@ -211,6 +232,8 @@ def qcp_rmsd(A: np.ndarray, B: np.ndarray) -> float:
         Coordinates of structure A
     B: numpy.ndarray
         Coordinates of structure B
+    atol: float
+        Absolute tolerance parameter (see notes)
 
     Returns
     -------
@@ -232,7 +255,7 @@ def qcp_rmsd(A: np.ndarray, B: np.ndarray) -> float:
     This means that :math:`s = G_a + G_bb - 2 * \\lambda_\\text{max}` can become
     negative because of numerical errors and therefore :math:`\\sqrt{s}` fails.
     In order to avoid this problem, the final RMSD is set to :math:`0`
-    if :math:`|s| < 10^{-12}`.
+    if :math:`|s| < atol`.
     """
 
     assert A.shape == B.shape
@@ -247,11 +270,18 @@ def qcp_rmsd(A: np.ndarray, B: np.ndarray) -> float:
 
     c2, c1, c0 = coefficients(M, K)
 
-    l_max = lambda_max(Ga, Gb, c2, c1, c0)
+    try:
+        # Fast calculation of the largest eigenvalue of K as root of the characteristic
+        # polynomial.
+        l_max = lambda_max(Ga, Gb, c2, c1, c0)
+    except RuntimeError:  # Newton method fails to converge; see GitHub Issue #35
+        # Fallback to (slower) explicit calculation of the largest eigenvalue of K
+        l_max = _lambda_max_eig(K)
 
     s = Ga + Gb - 2 * l_max
-    if abs(s) < 1e-12:  # Avoid numerical errors when Ga + Gb = 2 * l_max
-        rmsd = 0
+
+    if abs(s) < atol:  # Avoid numerical errors when Ga + Gb = 2 * l_max
+        rmsd = 0.0
     else:
         rmsd = np.sqrt(s / N)
 

spyrmsd/rmsd.py

@@ -12,6 +12,7 @@ def rmsd(
     atomicn2: np.ndarray,
     center: bool = False,
     minimize: bool = False,
+    atol: float = 1e-9,
 ) -> float:
     """
     Compute RMSD
@@ -30,6 +31,8 @@ def rmsd(
         Center molecules at origin
     minimize: bool
         Compute minimum RMSD (with QCP method)
+    atol: float
+        Absolute tolerance parameter for QCP method (see :func:`qcp_rmsd`)
 
     Returns
     -------
@@ -54,7 +57,7 @@ def rmsd(
     c2 = utils.center(coords2) if center or minimize else coords2
 
     if minimize:
-        rmsd = qcp.qcp_rmsd(c1, c2)
+        rmsd = qcp.qcp_rmsd(c1, c2, atol)
     else:
         n = coords1.shape[0]
 
@@ -121,6 +124,7 @@ def _rmsd_isomorphic_core(
     center: bool = False,
     minimize: bool = False,
     isomorphisms: Optional[List[Tuple[List[int], List[int]]]] = None,
+    atol: float = 1e-9,
 ) -> Tuple[float, List[Tuple[List[int], List[int]]]]:
     """
     Compute RMSD using graph isomorphism.
@@ -145,6 +149,8 @@ def _rmsd_isomorphic_core(
         Compute minized RMSD
     isomorphisms: Optional[List[Dict[int,int]]]
         Previously computed graph isomorphism
+    atol: float
+        Absolute tolerance parameter for QCP (see :func:`qcp_rmsd`)
 
     Returns
     -------
@@ -186,7 +192,7 @@ def _rmsd_isomorphic_core(
             result = np.sum((c1i - c2i) ** 2)
         else:
             # Compute minimized RMSD using QCP
-            result = qcp.qcp_rmsd(c1i, c2i)
+            result = qcp.qcp_rmsd(c1i, c2i, atol)
 
         min_result = result if result < min_result else min_result
 
@@ -208,6 +214,7 @@ def symmrmsd(
     center: bool = False,
     minimize: bool = False,
     cache: bool = True,
+    atol: float = 1e-9,
 ) -> Any:
     """
     Compute RMSD using graph isomorphism for multiple coordinates.
@@ -230,6 +237,10 @@ def symmrmsd(
         Centering flag
     minimize: bool
         Minimum RMSD
+    cache: bool
+        Cache graph isomorphisms
+    atol: float
+        Absolute tolerance parameter for QCP (see :func:`qcp_rmsd`)
 
     Returns
     -------
@@ -267,6 +278,7 @@ def symmrmsd(
                 center=center,
                 minimize=minimize,
                 isomorphisms=isomorphism,
+                atol=atol,
             )
 
             RMSD.append(srmsd)
@@ -283,6 +295,7 @@ def symmrmsd(
             center=center,
             minimize=minimize,
             isomorphisms=None,
+            atol=atol,
         )
 
     return RMSD

tests/conftest.py

@@ -9,21 +9,33 @@
 
 def pytest_addoption(parser):
     parser.addoption(
-        "--runslow", action="store_true", default=False, help="run slow tests"
+        "--benchmark", action="store_true", default=False, help="run benchmark"
+    )
+
+    parser.addoption(
+        "--large",
+        action="store_true",
+        default=False,
+        help="run large number of randomly selected tests",
     )
 
 
 def pytest_configure(config):
-    config.addinivalue_line("markers", "slow: mark test as slow to run")
+    config.addinivalue_line("markers", "benchmark: mark test as benchmark")
+    config.addinivalue_line("markers", "large: mark test as large")
 
 
 def pytest_collection_modifyitems(config, items):
-    if config.getoption("--runslow"):
-        # --runslow given in cli: do not skip slow tests
-        return
+    if not config.getoption("--benchmark"):
+        skip_benchmark = pytest.mark.skip(reason="need --benchmark option to run")
+
+        for item in items:
+            if "benchmark" in item.keywords:
+                item.add_marker(skip_benchmark)
 
-    skip_slow = pytest.mark.skip(reason="need --runslow option to run")
+    if not config.getoption("--large"):
+        skip_large = pytest.mark.skip(reason="need --large option to run")
 
-    for item in items:
-        if "slow" in item.keywords:
-            item.add_marker(skip_slow)
+        for item in items:
+            if "large" in item.keywords:
+                item.add_marker(skip_large)