Add polarizability

skfp/fingerprints/_new_mordred/calculator.py

@@ -14,6 +14,7 @@
     cpsa,
     extended_topochemical_atom,
     morse,
+    polarizability,
     rdkit_descriptors,
     ring_count,
     rotatable_bond,
@@ -120,6 +121,7 @@ def compute(mol: Mol, use_3D: bool) -> np.ndarray:
         aromatic.calc(mol_regular),
         topological_charge.calc(adjacency_matrix_regular, distance_matrix_regular),
         cpsa_2d,
+        polarizability.calc(mol_hydrogens),
     ]
 
     for values, feature_names in descriptors_2d:

skfp/fingerprints/_new_mordred/descriptors/polarizability.py

@@ -0,0 +1,32 @@
+import numpy as np
+from rdkit.Chem import Mol
+
+from skfp.fingerprints._new_mordred.utils.atomic_properties import get_polarizability
+
+"""
+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 = ["apol", "bpol"]
+
+
+def calc(mol_hydrogens: Mol) -> tuple[np.ndarray, list[str]]:
+    atom_polarizability = sum(
+        get_polarizability(atom) for atom in mol_hydrogens.GetAtoms()
+    )
+    bond_polarizability = sum(
+        abs(
+            get_polarizability(bond.GetBeginAtom())
+            - get_polarizability(bond.GetEndAtom())
+        )
+        for bond in mol_hydrogens.GetBonds()
+    )
+
+    values = np.array(
+        [atom_polarizability, bond_polarizability],
+        dtype=np.float32,
+    )
+    return values, FEATURE_NAMES

tests/fingerprints/_new_mordred/polarizability.py

@@ -0,0 +1,42 @@
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+
+from skfp.fingerprints._new_mordred.descriptors.polarizability import (
+    FEATURE_NAMES,
+    calc,
+)
+from skfp.fingerprints._new_mordred.utils.mol_preprocess import preprocess_mol
+
+"""
+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.
+"""
+
+
+@pytest.mark.parametrize(
+    "name, expected",
+    [
+        ("Hexane", [19.355102, 14.044898]),
+        ("Benzene", [14.020758, 6.019242]),
+        ("Caffeine", [26.03193, 18.03807]),
+        ("Lycopene", [104.140408, 56.179592]),
+        ("Epicatechin", [39.197102, 15.780898]),
+        ("Limonene", [27.368688, 16.051312]),
+        ("Allicin", [23.28993, 14.59007]),
+        ("Glutathione", [39.047481, 23.556519]),
+        ("Digoxin", [122.372752, 77.225248]),
+        ("Capsaicin", [51.569411, 30.260589]),
+        ("EllagicAcid", [33.796758, 11.227242]),
+        ("Astaxanthin", [104.681236, 53.902764]),
+    ],
+)
+def test_polarizability_values(name, expected, mordred_test_mols):
+    mol_hydrogens = preprocess_mol(mordred_test_mols[name], explicit_hydrogens=True)
+
+    values, feature_names = calc(mol_hydrogens)
+
+    assert feature_names == FEATURE_NAMES
+    assert_allclose(values, np.asarray(expected, dtype=np.float32), rtol=1e-5)

tests/fingerprints/new_mordred_fp.py

@@ -14,7 +14,14 @@ def test_new_mordred_fingerprint(smallest_mols_list):
     # temporary mask - will be eventually removed
     mask = ~(np.isnan(X_new) | np.isnan(X_old)).all(axis=0)
 
-    assert_allclose(X_new[:, mask], X_old[:, mask], equal_nan=True, atol=1e-3)
+    feature_names = new_mordred_fp.get_feature_names_out()
+    assert_allclose(
+        X_new[:, mask],
+        X_old[:, mask],
+        equal_nan=True,
+        atol=1e-3,
+        err_msg=_mismatched_features_msg(X_new, X_old, mask, feature_names, atol=1e-3),
+    )
     assert_equal(X_new.shape, (len(smallest_mols_list), 1613))
     assert X_new.dtype == np.float32
 
@@ -26,11 +33,19 @@ def test_new_mordred_sparse_fingerprint(smallest_mols_list):
     mordred_fp = MordredFingerprint(sparse=True, n_jobs=-1)
     X_old = mordred_fp.transform(smallest_mols_list)
 
+    X_new = X_new.toarray()
+    X_old = X_old.toarray()
+
     # temporary mask - will be eventually removed
-    mask = ~(np.isnan(X_new.toarray()) | np.isnan(X_old.toarray())).all(axis=0)
+    mask = ~(np.isnan(X_new) | np.isnan(X_old)).all(axis=0)
 
+    feature_names = new_mordred_fp.get_feature_names_out()
     assert_allclose(
-        X_new[:, mask].toarray(), X_old[:, mask].toarray(), equal_nan=True, atol=1e-3
+        X_new[:, mask],
+        X_old[:, mask],
+        equal_nan=True,
+        atol=1e-3,
+        err_msg=_mismatched_features_msg(X_new, X_old, mask, feature_names, atol=1e-3),
     )
     assert_equal(X_new.shape, (len(smallest_mols_list), 1613))
     assert X_new.dtype == np.float32
@@ -46,7 +61,14 @@ def test_new_mordred_3D_fingerprint(mols_conformers_list, smallest_mols_list):
     # temporary mask - will be eventually removed
     mask = ~(np.isnan(X_new) | np.isnan(X_old)).all(axis=0)
 
-    assert_allclose(X_new[:, mask], X_old[:, mask], equal_nan=True, atol=1e-2)
+    feature_names = new_mordred_fp.get_feature_names_out()
+    assert_allclose(
+        X_new[:, mask],
+        X_old[:, mask],
+        equal_nan=True,
+        atol=1e-2,
+        err_msg=_mismatched_features_msg(X_new, X_old, mask, feature_names, atol=1e-2),
+    )
     assert_equal(X_new.shape, (len(mols_conformers_list), 1826))
     assert X_new.dtype == np.float32
 
@@ -58,11 +80,19 @@ def test_new_mordred_3D_sparse_fingerprint(mols_conformers_list, smallest_mols_l
     mordred_fp = MordredFingerprint(use_3D=True, sparse=True, n_jobs=-1)
     X_old = mordred_fp.transform(smallest_mols_list)
 
+    X_new = X_new.toarray()
+    X_old = X_old.toarray()
+
     # temporary mask - will be eventually removed
-    mask = ~(np.isnan(X_new.toarray()) | np.isnan(X_old.toarray())).all(axis=0)
+    mask = ~(np.isnan(X_new) | np.isnan(X_old)).all(axis=0)
 
+    feature_names = new_mordred_fp.get_feature_names_out()
     assert_allclose(
-        X_new[:, mask].toarray(), X_old[:, mask].toarray(), equal_nan=True, atol=1e-2
+        X_new[:, mask],
+        X_old[:, mask],
+        equal_nan=True,
+        atol=1e-2,
+        err_msg=_mismatched_features_msg(X_new, X_old, mask, feature_names, atol=1e-2),
     )
     assert_equal(X_new.shape, (len(mols_conformers_list), 1826))
     assert X_new.dtype == np.float32
@@ -118,3 +148,20 @@ def test_new_mordred_3D_feature_names():
     feature_names_old = feature_names_old[~changed_name]
 
     assert_equal(feature_names_new, feature_names_old)
+
+
+def _mismatched_features_msg(
+    X_new: np.ndarray,
+    X_old: np.ndarray,
+    mask: np.ndarray,
+    feature_names: np.ndarray,
+    atol: float,
+) -> str:
+    # build an error message listing the feature names whose values
+    # differ beyond the tolerance
+    close = np.isclose(X_new[:, mask], X_old[:, mask], equal_nan=True, atol=atol)
+    mismatched_cols = ~close.all(axis=0)
+    mismatched_names = feature_names[mask][mismatched_cols]
+    if len(mismatched_names) == 0:
+        return ""
+    return "Mismatched feature names:\n" + "\n".join(mismatched_names)