Merge pull request #66 from RMeli/develop

Version 0.5.2
RMeli · Feb 23, 2022 · eb121af · eb121af
2 parents 9e9cd9a + f5183a9
commit eb121af
Show file tree

Hide file tree

Showing 14 changed files with 250 additions and 98 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -3,23 +3,46 @@
 
 ------------------------------------------------------------------------------
 
+## Version 0.5.2
+
+Date:            23/02/2022
+Contributors:    @RMeli
+
+### Fixed
+
+* Inconsistent number of nodes with `graph-tool` for disconnected graphs [PR #61 | @RMeli]
+
+### Improved
+
+* Support for more types of node properties (including strings) with `graph-tool` [PR #64 | @RMeli]
+
+### Changed
+
+* `ValueError` exception into `NonIsomorphicGraphs(ValueError)` exception [PR #65 | @RMeli]
+
+### Added
+
+* Warning for disconnected graphs [PR #61| @RMeli]
+
+------------------------------------------------------------------------------
+
 ## Version 0.5.1
 
 Date:            21/09/2021
 Contributors:    @RMeli
 
 ### Fixed
 
-* Fixed wrong covalent radius in `graph.adjacency_matrix_from_atomic_coordinates()` [PR #58 | @RMeli]
+* Wrong covalent radius in `graph.adjacency_matrix_from_atomic_coordinates()` [PR #58 | @RMeli]
 
 ### Added
 
-* Added [pre-commit](https://pre-commit.com/) configuration file [PR #57 | @RMeli]
-* Added support for gzip-compressed files (`.gz`) [PR #56 | @RMeli]
+* [pre-commit](https://pre-commit.com/) configuration file [PR #57 | @RMeli]
+* Support for gzip-compressed files (`.gz`) [PR #56 | @RMeli]
 
 ### Removed
 
-* Removed dependency `QCElemental` [PR #58 | @RMeli]
+* Dependency `QCElemental` [PR #58 | @RMeli]
 
 ------------------------------------------------------------------------------
 
@@ -31,15 +54,15 @@ Contributors:    @RMeli
 
 ### Added
 
-* Added `molecule.Molecule` constructor from RDKit molecule [PR #50 | @RMeli]
-* Added `molecule.Molecule` constructor from Open Babel molecule [PR #50 | @RMeli]
-* Added `--n-tests` option for `pytest` [PR #44 | @RMeli]
+* `molecule.Molecule` constructor from RDKit molecule [PR #50 | @RMeli]
+* `molecule.Molecule` constructor from Open Babel molecule [PR #50 | @RMeli]
+* `--n-tests` option for `pytest` [PR #44 | @RMeli]
 
 ### Improved
 
-* Improved `spyrmsd.rmsdwrapper` to deal with single molecule [PR #51 | @RMeli]
-* Improved issue template [PR #46 | @RMeli]
-* Improved speed of computation of squared pairwise distances [PR #45 | @RMeli]
+* `spyrmsd.rmsdwrapper` to deal with single molecule [PR #51 | @RMeli]
+* Issue template [PR #46 | @RMeli]
+* Speed of computation of squared pairwise distances [PR #45 | @RMeli]
 
 ### Changed
 
@@ -50,8 +73,8 @@ Contributors:    @RMeli
 
 ### Removed
 
-* Removed `spyrmsd` module [PR #52 | @RMeli]
-* Removed Travis CI and AppVeyor bindings [PR #44 | @RMeli]
-* Removed `--long` option for `pytest` [PR #44 | @RMeli]
+* `spyrmsd` module [PR #52 | @RMeli]
+* Travis CI and AppVeyor bindings [PR #44 | @RMeli]
+* `--long` option for `pytest` [PR #44 | @RMeli]
 
 ------------------------------------------------------------------------------
diff --git a/README.md b/README.md
@@ -122,14 +122,16 @@ The function  `rmsd.rmsd` computes RMSD without symmetry correction. The atoms a
 def rmsd(
     coords1: np.ndarray,    # Coordinates of molecule 1
     coords2: np.ndarray,    # Coordinates of molecule 2
-    atomicn1: np.ndarray,   # Atomic number of molecule 1
-    atomicn2: np.ndarray,   # Atomic number of molecule 2
+    aprops1: np.ndarray,    # Atomic properties of molecule 1
+    aprops2: np.ndarray,    # Atomic properties of molecule 2
     center: bool = False,   # Flag to center molecules at origin
     minimize: bool = False, # Flag to compute minimum RMSD
     atol: float = 1e-9,     # Numerical tolerance for QCP method
 )
 ```
 
+Note: atomic properties (`aprops`) can be any Python object when using [NetworkX](https://networkx.github.io/), or integers, floats, or strings when using [graph-tool](https://graph-tool.skewed.de/).
+
 #### Symmetry-Corrected RMSD
 
 The function `rmsd.symmrmsd` computes symmetry-corrected RMSD using molecular graph isomorphisms. Symmetry correction requires molecular adjacency matrices describing the connectivity but needs not the atoms to be in the same order.
@@ -140,8 +142,8 @@ Atom matching is performed according to the molecular graph. This function shoul
 def symmrmsd(
     coordsref: np.ndarray,                       # Reference coordinated
     coords: Union[np.ndarray, List[np.ndarray]], # Coordinates (one set or multiple sets)
-    atomicnumsref: np.ndarray,                   # Reference atomic numbers
-    atomicnums: np.ndarray,                      # Atomic numbers
+    apropsref: np.ndarray,                       # Reference atomic properties
+    aprops: np.ndarray,                          # Atomic properties
     amref: np.ndarray,                           # Reference adjacency matrix
     am: np.ndarray,                              # Adjacency matrix
     center: bool = False,                        # Flag to center molecules at origin
@@ -151,6 +153,8 @@ def symmrmsd(
 )
 ```
 
+Note: atomic properties (`aprops`) can be any Python object when using [NetworkX](https://networkx.github.io/), or integers, floats, or strings when using [graph-tool](https://graph-tool.skewed.de/).
+
 ## Development
 
 To ensure code quality and consistency the following tools are used during development:

diff --git a/docs/source/api/spyrmsd.exceptions.rst b/docs/source/api/spyrmsd.exceptions.rst
@@ -0,0 +1,7 @@
+spyrmsd.exceptions module
+=========================
+
+.. automodule:: spyrmsd.exceptions
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs/source/api/spyrmsd.rst b/docs/source/api/spyrmsd.rst
@@ -23,6 +23,7 @@ Submodules
 
    spyrmsd.constants
    spyrmsd.due
+   spyrmsd.exceptions
    spyrmsd.graph
    spyrmsd.hungarian
    spyrmsd.io

diff --git a/spyrmsd/exceptions.py b/spyrmsd/exceptions.py
@@ -0,0 +1,6 @@
+class NonIsomorphicGraphs(ValueError):
+    """
+    Raised when graphs are not isomorphic
+    """
+
+    pass
diff --git a/spyrmsd/graph.py b/spyrmsd/graph.py
@@ -40,15 +40,15 @@
 
 
 def adjacency_matrix_from_atomic_coordinates(
-    atomicnums: np.ndarray, coordinates: np.ndarray
+    aprops: np.ndarray, coordinates: np.ndarray
 ) -> np.ndarray:
     """
     Compute adjacency matrix from atomic coordinates.
 
     Parameters
     ----------
-    atomicnums: numpy.ndarray
-        Atomic numbers
+    aprops: numpy.ndarray
+        Atomic properties
     coordinates: numpy.ndarray
         Atomic coordinates
 
@@ -74,17 +74,17 @@ def adjacency_matrix_from_atomic_coordinates(
        (1991).
     """
 
-    n = len(atomicnums)
+    n = len(aprops)
 
     assert coordinates.shape == (n, 3)
 
     A = np.zeros((n, n))
 
     for i in range(n):
-        r_i = constants.anum_to_covalentradius[atomicnums[i]]
+        r_i = constants.anum_to_covalentradius[aprops[i]]
 
         for j in range(i + 1, n):
-            r_j = constants.anum_to_covalentradius[atomicnums[j]]
+            r_j = constants.anum_to_covalentradius[aprops[j]]
 
             distance = np.sqrt(np.sum((coordinates[i] - coordinates[j]) ** 2))
 

diff --git a/spyrmsd/graphs/_common.py b/spyrmsd/graphs/_common.py
@@ -0,0 +1,8 @@
+warn_disconnected_graph: str = "Disconnected graph detected. Is this expected?"
+warn_no_atomic_properties: str = (
+    "No atomic property information stored on nodes. Node matching is not performed..."
+)
+
+error_non_isomorphic_graphs: str = (
+    "Graphs are not isomorphic.\nMake sure graphs have the same connectivity."
+)
diff --git a/spyrmsd/graphs/gt.py b/spyrmsd/graphs/gt.py
@@ -5,10 +5,53 @@
 import numpy as np
 from graph_tool import generation, topology
 
+from spyrmsd.exceptions import NonIsomorphicGraphs
+from spyrmsd.graphs._common import (
+    error_non_isomorphic_graphs,
+    warn_disconnected_graph,
+    warn_no_atomic_properties,
+)
+
+
+# TODO: Implement all graph-tool supported types
+def _c_type(numpy_dtype):
+    """
+    Get C type compatible with graph-tool from numpy dtype
+
+    Parameters
+    ----------
+    numpy_dtype: np.dtype
+        Numpy dtype
+
+    Returns
+    -------
+    str
+        C type
+
+    Raises
+    ------
+    ValueError
+        If the data type is not supported
+
+    Notes
+    -----
+    https://graph-tool.skewed.de/static/doc/quickstart.html#sec-property-maps
+    """
+    name: str = numpy_dtype.name
+
+    if "int" in name:
+        return "int"
+    elif "float" in name:
+        return "double"
+    elif "str" in name:
+        return "string"
+    else:
+        raise ValueError(f"Unsupported property type: {name}")
+
 
 def graph_from_adjacency_matrix(
     adjacency_matrix: Union[np.ndarray, List[List[int]]],
-    atomicnums: Optional[Union[np.ndarray, List[int]]] = None,
+    aprops: Optional[Union[np.ndarray, List[Any]]] = None,
 ):
     """
     Graph from adjacency matrix.
@@ -17,8 +60,8 @@ def graph_from_adjacency_matrix(
     ----------
     adjacency_matrix: Union[np.ndarray, List[List[int]]]
         Adjacency matrix
-    atomicnums: Union[np.ndarray, List[int]], optional
-        Atomic numbers
+    aprops: Union[np.ndarray, List[Any]], optional
+        Atomic properties
 
     Returns
     -------
@@ -33,13 +76,27 @@ def graph_from_adjacency_matrix(
     # Get upper triangular adjacency matrix
     adj = np.triu(adjacency_matrix)
 
+    assert adj.shape[0] == adj.shape[1]
+    num_vertices = adj.shape[0]
+
     G = gt.Graph(directed=False)
+    G.add_vertex(n=num_vertices)
     G.add_edge_list(np.transpose(adj.nonzero()))
 
-    if atomicnums is not None:
-        vprop = G.new_vertex_property("short")  # Create property map (of C type short)
-        vprop.a = atomicnums  # Assign atomic numbers to property map array
-        G.vertex_properties["atomicnum"] = vprop  # Set property map
+    # Check if graph is connected, for warning
+    cc, _ = topology.label_components(G)
+    if set(cc.a) != {0}:
+        warnings.warn(warn_disconnected_graph)
+
+    if aprops is not None:
+        if not isinstance(aprops, np.ndarray):
+            aprops = np.array(aprops)
+
+        assert aprops.shape[0] == num_vertices
+
+        ptype: str = _c_type(aprops.dtype)  # Get C type
+        vprop = G.new_vertex_property(ptype, vals=aprops)  # Create property map
+        G.vertex_properties["aprops"] = vprop  # Set property map
 
     return G
 
@@ -62,7 +119,7 @@ def match_graphs(G1, G2) -> List[Tuple[List[int], List[int]]]:
 
     Raises
     ------
-    ValueError
+    NonIsomorphicGraphs
         If the graphs `G1` and `G2` are not isomorphic
     """
 
@@ -71,26 +128,19 @@ def match_graphs(G1, G2) -> List[Tuple[List[int], List[int]]]:
             G1,
             G2,
             vertex_label=(
-                G1.vertex_properties["atomicnum"],
-                G2.vertex_properties["atomicnum"],
+                G1.vertex_properties["aprops"],
+                G2.vertex_properties["aprops"],
             ),
             subgraph=False,
         )
-    except KeyError:  # No "atomicnum" vertex property
-        warnings.warn(
-            "No atomic number information stored on nodes. "
-            + "Node matching is not performed..."
-        )
+    except KeyError:
+        warnings.warn(warn_no_atomic_properties)
 
         maps = topology.subgraph_isomorphism(G1, G2, subgraph=False)
 
     # Check if graphs are actually isomorphic
     if len(maps) == 0:
-        # TODO: Create a new exception
-        raise ValueError(
-            "Graphs are not isomorphic."
-            "\nMake sure graphs have the same connectivity."
-        )
+        raise NonIsomorphicGraphs(error_non_isomorphic_graphs)
 
     n = num_vertices(G1)