update symbolic classifiers

Author: sfluegel05

README.md

@@ -10,8 +10,9 @@ Not all models can be installed automatically at the moment:
 - `chebai-graph` and its dependencies. To install them, follow
 the instructions in the [chebai-graph repository](https://github.com/ChEB-AI/python-chebai-graph).
 - `chemlog-extra` can be installed with `pip install git+https://github.com/ChEB-AI/chemlog-extra.git`
-- The automatically installed version of `c3p` may not work under Windows. If you want to run chebifier on Windows, we
-recommend using this forked version: `pip install git+https://github.com/sfluegel05/c3p.git`
+- `c3p` reads its generated programs assuming a UTF-8 locale and guards each of them with a
+SIGALRM-based timeout, neither of which holds on Windows. The `c3p` predictor works around both
+(see `_patch_c3p`), at the price of running the programs without a timeout there.
 
 
 You can get the package from PyPI:

chebifier/cli.py

@@ -234,7 +234,9 @@ def build(
     ensemble_param,
 ):
     """Build (calibrate) an ensemble on the ChEBI validation set."""
-    base_learners = build_base_learners(ensemble_config)
+    base_learners = build_base_learners(
+        ensemble_config, prediction_cache_dir=prediction_cache_dir, split="validation"
+    )
     ensemble_model = ENSEMBLES[ensemble_type](
         ensemble_dir, **parse_ensemble_params(ensemble_param)
     )

chebifier/prediction_models/c3p_predictor.py

@@ -1,10 +1,37 @@
+import functools
 from pathlib import Path
 from typing import List, Optional
 
 import tqdm
 
 from chebifier import modelwise_smiles_lru_cache
 from chebifier.prediction_models import BasePredictor
+from chebifier.utils import get_superclasses, to_smiles
+
+
+def _patch_c3p(c3p_classifier):
+    """Two things C3P 0.5.0 assumes that do not hold on Windows, both of which make it return no
+    classification at all rather than a wrong one:
+
+    - it reads its generated programs with `open(program, "r")`, i.e. in the platform default
+      encoding, while the programs are UTF-8. `open` is looked up in the module globals before the
+      builtins, so binding a UTF-8 `open` there fixes the reads without affecting any other module.
+    - it guards every program with timeout_decorator, which needs SIGALRM. Running the programs
+      without their 2s timeout is the only way to get predictions on a platform that has no
+      SIGALRM - timeout_decorator's signal-free mode forks a process per call, which is not
+      affordable for 300 programs per molecule.
+    """
+    import signal
+
+    if not hasattr(c3p_classifier, "open"):
+        c3p_classifier.open = functools.partial(open, encoding="utf-8")
+    if hasattr(signal, "SIGALRM"):
+        return
+    from c3p import learn
+
+    if hasattr(learn.eval_with_timeout, "__wrapped__"):
+        print("No SIGALRM on this platform, running C3P programs without a timeout.")
+        learn.eval_with_timeout = learn.eval_with_timeout.__wrapped__
 
 
 class C3PPredictor(BasePredictor):
@@ -28,6 +55,9 @@ def __init__(
     def predict_list(self, smiles_list: list[str]) -> list:
         from c3p import classifier as c3p_classifier
 
+        _patch_c3p(c3p_classifier)
+        # C3P only takes SMILES, while the evaluation datasets hand out RDKit molecules
+        smiles_list = [to_smiles(molecule) for molecule in smiles_list]
         result_list = []
         for batch_start in tqdm.tqdm(
             range(0, len(smiles_list), 32), desc="Classifying with C3P"
@@ -55,9 +85,7 @@ def predict_list(self, smiles_list: list[str]) -> list:
         for result in tqdm.tqdm(result_list, desc="Reformatting C3P results"):
             chebi_id = result.class_id.split(":")[1]
             if result.is_match and self.chebi_graph is not None:
-                parents = [
-                    str(parent) for parent in self.chebi_graph.predecessors(chebi_id)
-                ]
+                parents = get_superclasses(self.chebi_graph, chebi_id)
             else:
                 parents = []
             for idx in indices_by_smiles[result.input_smiles]:
@@ -74,6 +102,7 @@ def explain_smiles(self, smiles):
         """
         from c3p import classifier as c3p_classifier
 
+        _patch_c3p(c3p_classifier)
         highlights = []
         result_list = c3p_classifier.classify(
             [smiles], self.program_directory, self.chemical_classes, strict=False

chebifier/prediction_models/chebi_lookup.py

@@ -6,7 +6,7 @@
 
 from chebifier import modelwise_smiles_lru_cache
 from chebifier.prediction_models import BasePredictor
-from chebifier.utils import _smiles_to_mol, load_chebi_graph
+from chebifier.utils import _smiles_to_mol, get_superclasses, load_chebi_graph
 
 
 class ChEBILookupPredictor(BasePredictor):
@@ -61,7 +61,7 @@ def build_smiles_lookup(self):
                         smiles_lookup[canonical_smiles] = []
                     # if the canonical SMILES is already in the lookup, append "different interpretation of the SMILES"
                     smiles_lookup[canonical_smiles].append(
-                        (chebi_id, list(self.chebi_graph.predecessors(chebi_id)))
+                        (chebi_id, list(get_superclasses(self.chebi_graph, chebi_id)))
                     )
                 except Exception as e:
                     print(

chebifier/prediction_models/chemlog_predictor.py

@@ -4,6 +4,7 @@
 
 from chebifier import modelwise_smiles_lru_cache
 from chebifier.prediction_models.base_predictor import BasePredictor
+from chebifier.utils import CHEBI_VERSION, get_superclasses, to_mol
 
 AA_DICT = {
     "A": "L-alanine",
@@ -70,41 +71,41 @@ def predict_list(self, smiles_list: list[str]) -> list:
         return self._predict_smiles_list(smiles_list)
 
     def _predict_smiles_list(self, smiles_list: list[str]) -> list:
-        from chemlog.cli import _smiles_to_mol
-
-        mol_list = [_smiles_to_mol(smiles) for smiles in smiles_list]
+        mol_list = [to_mol(molecule) for molecule in smiles_list]
         res = self.classifier.classify(mol_list)
         if self.chebi_graph is not None:
             for sample in res:
                 sample_additions = dict()
                 for cls in sample:
                     if sample[cls] == 1:
-                        successors = list(self.chebi_graph.predecessors(cls))
-                        if successors:
-                            for succ in successors:
-                                sample_additions[str(succ)] = 1
+                        for superclass in get_superclasses(self.chebi_graph, cls):
+                            sample_additions[superclass] = 1
                 sample.update(sample_additions)
         return res
 
 
 class ChemlogXMolecularEntityPredictor(ChemlogExtraPredictor):
-    def __init__(self, model_name: str, **kwargs):
+    def __init__(self, model_name: str, chebi_version: int = CHEBI_VERSION, **kwargs):
         from chemlog_extra.alg_classification.by_element_classification import (
             XMolecularEntityClassifier,
         )
 
         super().__init__(model_name, **kwargs)
-        self.classifier = XMolecularEntityClassifier(chebi_graph=self.chebi_graph)
+        self.classifier = XMolecularEntityClassifier(
+            chebi_graph=self.chebi_graph, chebi_version=chebi_version
+        )
 
 
 class ChemlogOrganoXCompoundPredictor(ChemlogExtraPredictor):
-    def __init__(self, model_name: str, **kwargs):
+    def __init__(self, model_name: str, chebi_version: int = CHEBI_VERSION, **kwargs):
         from chemlog_extra.alg_classification.by_element_classification import (
             OrganoXCompoundClassifier,
         )
 
         super().__init__(model_name, **kwargs)
-        self.classifier = OrganoXCompoundClassifier(chebi_graph=self.chebi_graph)
+        self.classifier = OrganoXCompoundClassifier(
+            chebi_graph=self.chebi_graph, chebi_version=chebi_version
+        )
 
 
 class ChemlogLopsterPredictor(ChemlogExtraPredictor):
@@ -142,9 +143,9 @@ def __init__(self, model_name: str, **kwargs):
         print(f"Initialised ChemLog model {self.model_name}")
 
     def predict(self, smiles: str) -> Optional[dict]:
-        from chemlog.cli import _smiles_to_mol, strategy_call
+        from chemlog.cli import strategy_call
 
-        mol = _smiles_to_mol(smiles)
+        mol = to_mol(smiles)
         if mol is None:
             return None
         pos_labels = [
@@ -157,9 +158,9 @@ def predict(self, smiles: str) -> Optional[dict]:
         ]
         if self.chebi_graph:
             indirect_pos_labels = [
-                str(pr)
+                superclass
                 for label in pos_labels
-                for pr in self.chebi_graph.predecessors(label)
+                for superclass in get_superclasses(self.chebi_graph, label)
             ]
             pos_labels = list(set(pos_labels + indirect_pos_labels))
         return {
@@ -181,7 +182,7 @@ def _predict_smiles_list(self, smiles_list: list[str]) -> list:
 
         return results
 
-    def get_chemlog_result_info(self, smiles):
+    def get_chemlog_result_info(self, molecule):
         """Get classification for single molecule with additional information."""
         from chemlog.alg_classification.charge_classifier import get_charge_category
         from chemlog.alg_classification.peptide_size_classifier import (
@@ -194,10 +195,9 @@ def get_chemlog_result_info(self, smiles):
             is_diketopiperazine,
             is_emericellamide,
         )
-        from chemlog.cli import _smiles_to_mol
 
-        mol = _smiles_to_mol(smiles)
-        if mol is None or not smiles:
+        mol = to_mol(molecule) if molecule else None
+        if mol is None:
             return {"error": "Failed to parse SMILES"}
 
         charge_category = get_charge_category(mol)

chebifier/utils.py

@@ -3,11 +3,15 @@
 import os
 import pickle
 
+import networkx as nx
 import yaml
+from chebi_utils.obo_extractor import get_hierarchy_subgraph
 from rdkit import Chem
 
 from chebifier.hugging_face import download_model_files
 
+CHEBI_VERSION = 252
+
 
 def load_chebi_graph(filename=None):
     """Load ChEBI graph from Hugging Face (if filename is None) or local file"""
@@ -17,7 +21,7 @@ def load_chebi_graph(filename=None):
             {
                 "repo_id": "chebai/chebifier",
                 "repo_type": "dataset",
-                "files": {"f": "chebi_graph_v252.pkl"},
+                "files": {"f": f"chebi_graph_v{CHEBI_VERSION}.pkl"},
             }
         )["f"]
     else:
@@ -89,3 +93,40 @@ def _smiles_to_mol(smiles: str):
         except Chem.KekulizeException as e:
             print(f"Failed to Kekulize {smiles}: {e}")
     return mol
+
+
+def to_mol(molecule: str | Chem.Mol):
+    """Molecules reach a predictor either as SMILES or as RDKit molecules (the evaluation datasets
+    store the latter). Rule-based classifiers expect kekulised molecules, and Kekulize works in
+    place, so a molecule that is not ours to modify is copied first."""
+    if not isinstance(molecule, Chem.Mol):
+        return _smiles_to_mol(molecule)
+    molecule = Chem.Mol(molecule)
+    try:
+        Chem.Kekulize(molecule)
+    except Chem.KekulizeException as e:
+        print(f"Failed to Kekulize {Chem.MolToSmiles(molecule)}: {e}")
+    return molecule
+
+
+def to_smiles(molecule: str | Chem.Mol) -> str:
+    return Chem.MolToSmiles(molecule) if isinstance(molecule, Chem.Mol) else molecule
+
+
+@functools.lru_cache(maxsize=2)
+def _isa_graph(chebi_graph):
+    return get_hierarchy_subgraph(chebi_graph)
+
+
+@functools.lru_cache(maxsize=None)
+def get_superclasses(chebi_graph, chebi_id: str) -> tuple[str, ...]:
+    """All transitive superclasses of a ChEBI class.
+
+    is-a edges point from child to parent, and the graph also carries non-subsumption relations
+    (has role, conjugate acid/base, ...), so the superclasses of a node are neither its
+    predecessors nor all of its successors.
+    """
+    isa_graph = _isa_graph(chebi_graph)
+    if chebi_id not in isa_graph:
+        return ()
+    return tuple(str(cls) for cls in nx.descendants(isa_graph, chebi_id))