Commit c6c87c5f authored by nd-02110114's avatar nd-02110114
Browse files

✨ add cgcnn featurizer

parent 1b7d83bd
Loading
Loading
Loading
Loading
+1 −1
Original line number Diff line number Diff line
@@ -22,4 +22,4 @@ from deepchem.feat.atomic_coordinates import AtomicCoordinates
from deepchem.feat.atomic_coordinates import NeighborListComplexAtomicCoordinates
from deepchem.feat.adjacency_fingerprints import AdjacencyFingerprint
from deepchem.feat.smiles_featurizers import SmilesToSeq, SmilesToImage
from deepchem.feat.materials_featurizers import ElementPropertyFingerprint, SineCoulombMatrix, StructureGraphFeaturizer
from deepchem.feat.materials_featurizers import ElementPropertyFingerprint, SineCoulombMatrix, CGCNNFeaturizer
+39 −38
Original line number Diff line number Diff line
@@ -2,8 +2,8 @@ from typing import Optional, Sequence
import numpy as np


class MoleculeGraphData:
  """MoleculeGraphData class
class GraphData:
  """GraphData class

  This data class is almost same as `torch_geometric.data.Data 
  <https://pytorch-geometric.readthedocs.io/en/latest/modules/data.html#torch_geometric.data.Data>`_.
@@ -14,8 +14,6 @@ class MoleculeGraphData:
    Node feature matrix with shape [num_nodes, num_node_features]
  edge_index: np.ndarray, dtype int
    Graph connectivity in COO format with shape [2, num_edges]
  targets : np.ndarray
    Graph or node targets with arbitrary shape
  edge_features: np.ndarray, optional (default None)
    Edge feature matrix with shape [num_edges, num_edge_features]
  graph_features: np.ndarray, optional (default None)
@@ -26,7 +24,7 @@ class MoleculeGraphData:
    The number of features per node in the graph
  num_edges: int
    The number of edges in the graph
  num_edges_features : int, , optional (default None)
  num_edges_features: int, optional (default None)
    The number of features per edge in the graph
  """

@@ -34,7 +32,6 @@ class MoleculeGraphData:
      self,
      node_features: np.ndarray,
      edge_index: np.ndarray,
      targets: np.ndarray,
      edge_features: Optional[np.ndarray] = None,
      graph_features: Optional[np.ndarray] = None,
  ):
@@ -45,8 +42,6 @@ class MoleculeGraphData:
      Node feature matrix with shape [num_nodes, num_node_features]
    edge_index: np.ndarray, dtype int
      Graph connectivity in COO format with shape [2, num_edges]
    targets : np.ndarray
      Graph or node targets with arbitrary shape
    edge_features: np.ndarray, optional (default None)
      Edge feature matrix with shape [num_edges, num_edge_features]
    graph_features: np.ndarray, optional (default None)
@@ -55,20 +50,21 @@ class MoleculeGraphData:
    # validate params
    if isinstance(node_features, np.ndarray) is False:
      raise ValueError('node_features must be np.ndarray.')

    if isinstance(edge_index, np.ndarray) is False:
      raise ValueError('edge_index must be np.ndarray.')
    elif edge_index.dtype != np.int:
      raise ValueError('edge_index.dtype must be np.int')
    elif edge_index.shape[0] != 2:
      raise ValueError('The shape of edge_index is [2, num_edges].')
    if isinstance(targets, np.ndarray) is False:
      raise ValueError('y must be np.ndarray.')

    if edge_features is not None:
      if isinstance(edge_features, np.ndarray) is False:
        raise ValueError('edge_features must be np.ndarray or None.')
      elif edge_index.shape[1] != edge_features.shape[0]:
        raise ValueError('The first dimension of edge_features must be the \
                    same as the second dimension of edge_index.')

    if graph_features is not None and isinstance(graph_features,
                                                 np.ndarray) is False:
      raise ValueError('graph_features must be np.ndarray or None.')
@@ -77,37 +73,42 @@ class MoleculeGraphData:
    self.edge_index = edge_index
    self.edge_features = edge_features
    self.graph_features = graph_features
    self.targets = targets
    self.num_nodes, self.num_node_features = self.node_features.shape
    self.num_edges = edge_index.shape[1]
    if self.node_features is not None:
      self.num_edge_features = self.edge_features.shape[1]

  def to_pyg_data(self):
  def to_pyg_data(self, target):
    """Convert to PyTorch Geometric Data instance

    Parameters
    ----------
    target: np.ndarray
      Graph or node targets with arbitrary shape

    Returns
    -------
    torch_geometric.data.Data
      Molecule graph data for PyTorch Geometric
      Graph data for PyTorch Geometric
    """
    try:
      import torch
      from torch_geometric.data import Data
    except ModuleNotFoundError:
      raise ValueError("This class requires PyTorch Geometric to be installed.")
      raise ValueError(
          "This function requires PyTorch Geometric to be installed.")

    return Data(
      x=torch.from_numpy(self.node_features),
      edge_index=torch.from_numpy(self.edge_index),
      edge_attr=None if self.edge_features is None \
        else torch.from_numpy(self.edge_features),
      y=torch.from_numpy(self.targets),
      y=torch.from_numpy(target),
    )


class BatchMoleculeGraphData(MoleculeGraphData):
  """Batch MoleculeGraphData class
class BatchGraphData(GraphData):
  """Batch GraphData class

  Attributes
  ----------
@@ -115,37 +116,34 @@ class BatchMoleculeGraphData(MoleculeGraphData):
    This vector indicates which graph the node belongs with shape [num_nodes,]
  """

  def __init__(self, molecule_graphs: Sequence[MoleculeGraphData]):
  def __init__(self, graphs: Sequence[GraphData]):
    """
    Parameters
    ----------
    molecule_graphs : Sequence[MoleculeGraphData]
      List of MoleculeGraphData
    graphs: Sequence[GraphData]
      List of GraphData
    """
    # stack features and targets
    batch_node_features = np.vstack(
        [graph.node_features for graph in molecule_graphs])
    batch_targets = np.vstack([graph.targets for graph in molecule_graphs])
    # stack features
    batch_node_features = np.vstack([graph.node_features for graph in graphs])

    # before stacking edge_features or graph_features,
    # we should check whether these are None or not
    if molecule_graphs[0].edge_features is not None:
      batch_edge_features = np.vstack(
          [graph.edge_features for graph in molecule_graphs])
    if graphs[0].edge_features is not None:
      batch_edge_features = np.vstack([graph.edge_features for graph in graphs])
    else:
      batch_edge_features = None

    if molecule_graphs[0].graph_features is not None:
    if graphs[0].graph_features is not None:
      batch_graph_features = np.vstack(
          [graph.graph_features for graph in molecule_graphs])
          [graph.graph_features for graph in graphs])
    else:
      batch_graph_features = None

    # create new edge index
    num_nodes_list = [graph.num_nodes for graph in molecule_graphs]
    num_nodes_list = [graph.num_nodes for graph in graphs]
    batch_edge_index = np.hstack(
      [graph.edge_index + prev_num_node for prev_num_node, graph \
        in zip([0] + num_nodes_list[:-1], molecule_graphs)]
        in zip([0] + num_nodes_list[:-1], graphs)]
    ).astype(int)

    # graph_index indicates which nodes belong to which graph
@@ -157,30 +155,33 @@ class BatchMoleculeGraphData(MoleculeGraphData):
    super().__init__(
        node_features=batch_node_features,
        edge_index=batch_edge_index,
        targets=batch_targets,
        edge_features=batch_edge_features,
        graph_features=batch_graph_features,
    )

    @staticmethod  # type: ignore
    def to_pyg_data(molecule_graphs: Sequence[MoleculeGraphData]):
    def to_pyg_data(graphs: Sequence[GraphData], targets: Sequence[np.ndarray]):
      """Convert to PyTorch Geometric Batch instance

      Parameters
      ----------
      molecule_graphs : Sequence[MoleculeGraphData]
        List of MoleculeGraphData
      graphs: Sequence[GraphData]
        List of GraphData
      targets: Sequence[np.ndarray]
        List of graph or node targets with arbitrary shape

      Returns
      -------
      torch_geometric.data.Batch
        Batch data of molecule graph for PyTorch Geometric
        Batch data of graphs for PyTorch Geometric
      """
      try:
        from torch_geometric.data import Batch
      except ModuleNotFoundError:
        raise ValueError(
            "This class requires PyTorch Geometric to be installed.")
            "This function requires PyTorch Geometric to be installed.")

      data_list = [mol_graph.to_pyg_data() for mol_graph in molecule_graphs]
      data_list = [
          graph.to_pyg_data(target) for graph, target in zip(graphs, targets)
      ]
      return Batch.from_data_list(data_list=data_list)
+0 −298
Original line number Diff line number Diff line
"""
Featurizers for inorganic crystals.
"""

import numpy as np

from deepchem.feat import MaterialStructureFeaturizer, MaterialCompositionFeaturizer
from deepchem.utils import pad_array


class ElementPropertyFingerprint(MaterialCompositionFeaturizer):
  """
  Fingerprint of elemental properties from composition.

  Based on the data source chosen, returns properties and statistics
  (min, max, range, mean, standard deviation, mode) for a compound
  based on elemental stoichiometry. E.g., the average electronegativity
  of atoms in a crystal structure. The chemical fingerprint is a 
  vector of these statistics. For a full list of properties and statistics,
  see ``matminer.featurizers.composition.ElementProperty(data_source).feature_labels()``.

  This featurizer requires the optional dependencies pymatgen and
  matminer. It may be useful when only crystal compositions are available
  (and not 3D coordinates).

  See references [1]_ [2]_ [3]_ [4]_ for more details.

  References
  ----------
  .. [1] MagPie data: Ward, L. et al. npj Comput Mater 2, 16028 (2016).
	 https://doi.org/10.1038/npjcompumats.2016.28

  .. [2] Deml data: Deml, A. et al. Physical Review B 93, 085142 (2016).
	 10.1103/PhysRevB.93.085142

  .. [3] Matminer: Ward, L. et al. Comput. Mater. Sci. 152, 60-69 (2018).

  .. [4] Pymatgen: Ong, S.P. et al. Comput. Mater. Sci. 68, 314-319 (2013). 

  """

  def __init__(self, data_source='matminer'):
    """
    Parameters
    ----------
    data_source : {"matminer", "magpie", "deml"}
      Source for element property data.

    """

    self.data_source = data_source

  def _featurize(self, composition):
    """
    Calculate chemical fingerprint from crystal composition.

    Parameters
    ----------
    composition: pymatgen.Composition object
      Composition object.

    Returns
    -------
    feats: np.ndarray
      Vector of properties and statistics derived from chemical
      stoichiometry. Some values may be NaN.

    """
    try:
      from matminer.featurizers.composition import ElementProperty
    except ModuleNotFoundError:
      raise ValueError("This class requires matminer to be installed.")

    ep = ElementProperty.from_preset(self.data_source)

    try:
      feats = ep.featurize(composition)
    except:
      feats = []

    return np.array(feats)


class SineCoulombMatrix(MaterialStructureFeaturizer):
  """
  Calculate sine Coulomb matrix for crystals.

  A variant of Coulomb matrix for periodic crystals.

  The sine Coulomb matrix is identical to the Coulomb matrix, except
  that the inverse distance function is replaced by the inverse of
  sin**2 of the vector between sites which are periodic in the 
  dimensions of the crystal lattice.

  Features are flattened into a vector of matrix eigenvalues by default
  for ML-readiness. To ensure that all feature vectors are equal
  length, the maximum number of atoms (eigenvalues) in the input
  dataset must be specified.

  This featurizer requires the optional dependencies pymatgen and
  matminer. It may be useful when crystal structures with 3D coordinates 
  are available.

  See [1]_ for more details.

  References
  ----------
  .. [1] Faber et al. Inter. J. Quantum Chem. 115, 16, 2015.

  """

  def __init__(self, max_atoms, flatten=True):
    """
    Parameters
    ----------
    max_atoms : int
      Maximum number of atoms for any crystal in the dataset. Used to
      pad the Coulomb matrix.
    flatten : bool (default True)
      Return flattened vector of matrix eigenvalues.

    """

    self.max_atoms = int(max_atoms)
    self.flatten = flatten

  def _featurize(self, struct):
    """
    Calculate sine Coulomb matrix from pymatgen structure.

    Parameters
    ----------
    struct : pymatgen.Structure
      A periodic crystal composed of a lattice and a sequence of atomic
      sites with 3D coordinates and elements.
      
    Returns
    -------
    features: np.ndarray
      2D sine Coulomb matrix with shape (max_atoms, max_atoms),
      or 1D matrix eigenvalues with shape (max_atoms,). 

    """

    try:
      from matminer.featurizers.structure import SineCoulombMatrix as SCM
    except ModuleNotFoundError:
      raise ValueError("This class requires matminer to be installed.")

    # Get full N x N SCM
    scm = SCM(flatten=False)
    sine_mat = scm.featurize(struct)

    if self.flatten:
      eigs, _ = np.linalg.eig(sine_mat)
      zeros = np.zeros((1, self.max_atoms))
      zeros[:len(eigs)] = eigs
      features = zeros
    else:
      features = pad_array(sine_mat, self.max_atoms)

    features = np.asarray(features)

    return features


class StructureGraphFeaturizer(MaterialStructureFeaturizer):
  """
  Calculate structure graph features for crystals.

  Based on the implementation in Crystal Graph Convolutional
  Neural Networks (CGCNN). The method constructs a crystal graph
  representation including atom features (atomic numbers) and bond
  features (neighbor distances). Neighbors are determined by searching
  in a sphere around atoms in the unit cell. A Gaussian filter is
  applied to neighbor distances. All units are in angstrom.  

  This featurizer requires the optional dependency pymatgen. It may
  be useful when 3D coordinates are available and when using graph 
  network models and crystal graph convolutional networks.

  See [1]_ for more details.

  References
  ----------
  .. [1] T. Xie and J. C. Grossman, Phys. Rev. Lett. 120, 2018.

  """

  def __init__(self, radius=8.0, max_neighbors=12, step=0.2):
    """
    Parameters
    ----------
    radius : float (default 8.0)
      Radius of sphere for finding neighbors of atoms in unit cell.
    max_neighbors : int (default 12)
      Maximum number of neighbors to consider when constructing graph.
    step : float (default 0.2)
      Step size for Gaussian filter.

    """

    self.radius = radius
    self.max_neighbors = int(max_neighbors)
    self.step = step

  def _featurize(self, struct):
    """
    Calculate crystal graph features from pymatgen structure.

    Parameters
    ----------
    struct : pymatgen.Structure
      A periodic crystal composed of a lattice and a sequence of atomic
      sites with 3D coordinates and elements.

    Returns
    -------
    feats: np.array
      Atomic and bond features. Atomic features are atomic numbers 
      and bond features are Gaussian filtered interatomic distances.

    """

    features = self._get_structure_graph_features(struct)
    features = np.array(features)

    return features

  def _get_structure_graph_features(self, struct):
    """
    Calculate structure graph features from pymatgen structure.

    Parameters
    ----------
    struct : pymatgen.Structure
      A periodic crystal composed of a lattice and a sequence of atomic
      sites with 3D coordinates and elements.

    Returns
    -------
    feats: tuple[np.array]
      atomic numbers, filtered interatomic distance tensor, and neighbor ids
    
    """

    atom_features = np.array([site.specie.Z for site in struct], dtype='int32')

    neighbors = struct.get_all_neighbors(self.radius, include_index=True)
    neighbors = [sorted(n, key=lambda x: x[1]) for n in neighbors]

    # Get list of lists of neighbor distances
    neighbor_features, neighbor_idx = [], []
    for neighbor in neighbors:
      if len(neighbor) < self.max_neighbors:
        neighbor_idx.append(
            list(map(lambda x: x[2], neighbor)) +
            [0] * (self.max_neighbors - len(neighbor)))
        neighbor_features.append(
            list(map(lambda x: x[1], neighbor)) +
            [self.radius + 1.] * (self.max_neighbors - len(neighbor)))
      else:
        neighbor_idx.append(
            list(map(lambda x: x[2], neighbor[:self.max_neighbors])))
        neighbor_features.append(
            list(map(lambda x: x[1], neighbor[:self.max_neighbors])))

    neighbor_features = np.array(neighbor_features)
    neighbor_idx = np.array(neighbor_idx)
    neighbor_features = self._gaussian_filter(neighbor_features)
    neighbor_features = np.vstack(neighbor_features)

    return (atom_features, neighbor_features, neighbor_idx)

  def _gaussian_filter(self, distances):
    """
    Apply Gaussian filter to an array of interatomic distances.

    Parameters
    ----------
    distances : np.array
      Matrix of distances of dimension (num atoms) x (max neighbors). 

    Returns
    -------
    expanded_distances: np.array 
      Expanded distance tensor after Gaussian filtering. Dimensionality
      is (num atoms) x (max neighbors) x (len(filt))
    
    """

    filt = np.arange(0, self.radius + self.step, self.step)

    # Increase dimension of distance tensor and apply filter
    expanded_distances = np.exp(
        -(distances[..., np.newaxis] - filt)**2 / self.step**2)

    return expanded_distances
+6 −0
Original line number Diff line number Diff line
"""
Featurizers for inorganic crystals.
"""
from deepchem.feat.materials_featurizers.element_property_fingerprint import ElementPropertyFingerprint
from deepchem.feat.materials_featurizers.sine_coulomb_matrix import SineCoulombMatrix
from deepchem.feat.materials_featurizers.cgcnn_featurizer import CGCNNFeaturizer
+179 −0
Original line number Diff line number Diff line
import os
import json
import numpy as np
from typing import Tuple

from deepchem.utils import download_url, get_data_dir
from deepchem.utils.typing import PymatgenStructure
from deepchem.feat import MaterialStructureFeaturizer
from deepchem.feat.graph_data import GraphData

# FIXME: it is better to add this json to DeepChem AWS
ATOM_JSON_URL = 'https://raw.githubusercontent.com/txie-93/cgcnn/master/data/sample-regression/atom_init.json'


class CGCNNFeaturizer(MaterialStructureFeaturizer):
  """
  Calculate structure graph features for crystals.

  Based on the implementation in Crystal Graph Convolutional
  Neural Networks (CGCNN). The method constructs a crystal graph
  representation including atom features and bond features (neighbor
  distances). Neighbors are determined by searching in a sphere around
  atoms in the unit cell. A Gaussian filter is applied to neighbor distances.
  All units are in angstrom.

  This featurizer requires the optional dependency pymatgen. It may
  be useful when 3D coordinates are available and when using graph
  network models and crystal graph convolutional networks.

  See [1]_ for more details.

  References
  ----------
  .. [1] T. Xie and J. C. Grossman, Phys. Rev. Lett. 120, 2018.

  Note
  ----
  This class requires Pymatgen to be installed.
  """

  def __init__(self,
               radius: float = 8.0,
               max_neighbors: float = 8,
               step: float = 0.2):
    """
    Parameters
    ----------
    radius: float (default 8.0)
      Radius of sphere for finding neighbors of atoms in unit cell.
    max_neighbors: int (default 8)
      Maximum number of neighbors to consider when constructing graph.
    step: float (default 0.2)
      Step size for Gaussian filter. This value is used when building edge features.
    """

    self.radius = radius
    self.max_neighbors = int(max_neighbors)
    self.step = step

    # load atom_init.json
    data_dir = get_data_dir()
    download_url(ATOM_JSON_URL, data_dir)
    atom_init_json_path = os.path.join(data_dir, 'atom_init.json')
    with open(atom_init_json_path, 'r') as f:
      atom_init_json = json.load(f)

    self.atom_features = {
        int(key): np.array(value, dtype=np.float32)
        for key, value in atom_init_json.items()
    }
    self.valid_atom_number = set(self.atom_features.keys())

  def _featurize(self, struct: PymatgenStructure) -> GraphData:
    """
    Calculate crystal graph features from pymatgen structure.

    Parameters
    ----------
    struct: pymatgen.Structure
      A periodic crystal composed of a lattice and a sequence of atomic
      sites with 3D coordinates and elements.

    Returns
    -------
    graph: GraphData
      A crystal graph with CGCNN style features.
    """

    node_features = self._get_node_features(struct)
    edge_index, edge_features = self._get_edge_features_and_index(struct)
    graph = GraphData(node_features, edge_index, edge_features)
    return graph

  def _get_node_features(self, struct: PymatgenStructure) -> np.ndarray:
    """
    Get the node feature from `atom_init.json`. The `atom_init.json` was collected
    from `data/sample-regression/atom_init.json` in the CGCNN repository.

    Parameters
    ----------
    struct: pymatgen.Structure
      A periodic crystal composed of a lattice and a sequence of atomic
      sites with 3D coordinates and elements.

    Returns
    -------
    node_features: np.ndarray
      A numpy array of shape `(num_nodes, 92)`.
    """
    node_features = []
    for site in struct:
      # check whether the atom feature exists or not
      assert site.specie.number in self.valid_atom_number
      node_features.append(self.atom_features[site.specie.number])
    node_features = np.vstack(node_features).astype(np.float)
    return node_features

  def _get_edge_features_and_index(
      self, struct: PymatgenStructure) -> Tuple[np.ndarray, np.ndarray]:
    """
    Calculate the edge feature and edge index from pymatgen structure.

    Parameters
    ----------
    struct: pymatgen.Structure
      A periodic crystal composed of a lattice and a sequence of atomic
      sites with 3D coordinates and elements.

    Returns
    -------
    edge_idx np.ndarray, dtype int
      A numpy array of shape with `(2, num_edges)`.
    edge_features: np.ndarray
      A numpy array of shape with `(num_edges, filter_length)`. The `filter_length` is
      (self.radius / self.step) + 1. The edge features were built by applying gaussian
      filter to the distance between nodes.
    """

    neighbors = struct.get_all_neighbors(self.radius, include_index=True)
    neighbors = [sorted(n, key=lambda x: x[1]) for n in neighbors]

    # construct bi-directed graph
    src_idx, dest_idx = [], []
    edge_distances = []
    for node_idx, neighbor in enumerate(neighbors):
      neighbor = neighbor[:self.max_neighbors]
      src_idx.extend([node_idx] * len(neighbor))
      dest_idx.extend([site[2] for site in neighbor])
      edge_distances.extend([site[1] for site in neighbor])

    edge_idx = np.array([src_idx, dest_idx], dtype=np.int)
    edge_distances = np.asarray(edge_distances)
    edge_features = self._gaussian_filter(edge_distances)
    return edge_idx, edge_features

  def _gaussian_filter(self, distances: np.ndarray) -> np.ndarray:
    """
    Apply Gaussian filter to an array of interatomic distances.

    Parameters
    ----------
    distances : np.ndarray
      A numpy array of the shape `(num_edges, )`.

    Returns
    -------
    expanded_distances: np.ndarray
      Expanded distance tensor after Gaussian filtering.
      The shape is `(num_edges, filter_length)`. The `filter_length` is
      (self.radius / self.step) + 1.
    """

    filt = np.arange(0, self.radius + self.step, self.step)

    # Increase dimension of distance tensor and apply filter
    expanded_distances = np.exp(
        -(distances[..., np.newaxis] - filt)**2 / self.step**2)

    return expanded_distances
Loading