Constitutional descriptor added

skfp/fingerprints/_new_mordred/calculator.py

@@ -11,6 +11,7 @@
     barysz_matrix,
     bond_count,
     carbon_types,
+    constitutional,
     extended_topochemical_atom,
     morse,
     rdkit_descriptors,
@@ -97,6 +98,7 @@ def compute(mol: Mol, use_3D: bool) -> np.ndarray:
         atom_count.calc(mol_regular),
         bond_count.calc(mol_hydrogens, mol_kekulized_hydrogens),
         carbon_types.calc(mol_kekulized),
+        constitutional.calc(mol_hydrogens),
         rotatable_bond.calc(mol_regular),
         vertex_adjacency_info.calc(mol_regular),
         ring_count.calc(mol_regular),

skfp/fingerprints/_new_mordred/descriptors/constitutional.py

@@ -0,0 +1,72 @@
+"""Constitutional descriptors."""
+
+import numpy as np
+from rdkit.Chem import Mol
+from rdkit.Chem.rdchem import Atom
+
+from skfp.fingerprints._new_mordred.utils.atomic_properties import (
+    PROPERTY_FUNCS,
+)
+
+"""
+This code has been adapted from the BSD-licensed mordred-community library.
+https://github.com/JacksonBurns/mordred-community
+
+See skfp/fingerprints/data/mordred-community_bsd_license.txt for the license text.
+"""
+
+FEATURE_NAMES = [
+    "SZ",
+    "Sm",
+    "Sv",
+    "Sse",
+    "Spe",
+    "Sare",
+    "Sp",
+    "Si",
+    "MZ",
+    "Mm",
+    "Mv",
+    "Mse",
+    "Mpe",
+    "Mare",
+    "Mp",
+    "Mi",
+]
+
+_NUM_ELEMENTS = 119
+
+
+def _get_normalized_property_table() -> np.ndarray:
+    """Precompute each atomic property divided by its value for carbon."""
+    atoms = tuple(Atom(atomic_num) for atomic_num in range(_NUM_ELEMENTS))
+    property_values = np.asarray(
+        [
+            [property_func(atom) for atom in atoms]
+            for property_func in PROPERTY_FUNCS.values()
+        ],
+        dtype=np.float64,
+    )
+    carbon_values = property_values[:, [6]]
+    return property_values / carbon_values
+
+
+_CARBON_NORMALIZED_PROPERTIES = _get_normalized_property_table()
+
+
+def calc(mol_hydrogens: Mol) -> tuple[np.ndarray, list[str]]:
+    """
+    Compute the Mordred constitutional descriptors.
+
+    For each atomic property, the property values of all atoms, including explicit
+    hydrogens, are normalized by the corresponding value for carbon. The `S*`
+    descriptors are their sums and the `M*` descriptors are their means.
+    """
+    atomic_numbers = np.fromiter(
+        (atom.GetAtomicNum() for atom in mol_hydrogens.GetAtoms()),
+        dtype=np.intp,
+        count=mol_hydrogens.GetNumAtoms(),
+    )
+    sums = _CARBON_NORMALIZED_PROPERTIES[:, atomic_numbers].sum(axis=1)
+    means = sums / atomic_numbers.size if atomic_numbers.size else np.full(8, np.nan)
+    return np.concatenate((sums, means)).astype(np.float32), FEATURE_NAMES

tests/fingerprints/_new_mordred/constitutional.py

@@ -0,0 +1,37 @@
+import numpy as np
+import pytest
+from mordred import Calculator, Constitutional
+from numpy.testing import assert_allclose
+from rdkit import Chem
+
+from skfp.fingerprints._new_mordred import calculator
+from skfp.fingerprints._new_mordred.descriptors import constitutional
+from skfp.fingerprints._new_mordred.utils.mol_preprocess import preprocess_mol
+
+
+@pytest.mark.parametrize(
+    "smiles",
+    ["CCO", "c1ccccc1", "C(F)(Cl)(Br)I", "[Na+]", "[H][H]", ""],
+)
+def test_constitutional_matches_mordred(smiles):
+    mol = Chem.MolFromSmiles(smiles)
+    mol_hydrogens = preprocess_mol(mol, explicit_hydrogens=True)
+
+    values, feature_names = constitutional.calc(mol_hydrogens)
+    mordred_calc = Calculator(Constitutional)
+    expected = np.asarray(list(mordred_calc(mol)), dtype=np.float32)
+
+    assert feature_names == [str(desc) for desc in mordred_calc.descriptors]
+    assert_allclose(values, expected, rtol=1e-6, equal_nan=True)
+
+
+def test_constitutional_is_connected_to_calculator():
+    mol = Chem.MolFromSmiles("CCO")
+    result = calculator.compute(mol, use_3D=False)
+    names = calculator.get_feature_names(use_3d=False)
+    expected, feature_names = constitutional.calc(
+        preprocess_mol(mol, explicit_hydrogens=True)
+    )
+
+    indices = [np.flatnonzero(names == name)[0] for name in feature_names]
+    assert_allclose(result[indices], expected)