Clean up code for tdims

README.md

@@ -52,6 +52,7 @@ Each folder represents a different chemical representation and architecture type
 | 3DGrid-VQGAN [\[HuggingFace\]](https://huggingface.co/ibm-research/materials.3dgrid_vqgan) | 3DGrid-VQGAN is an encoder-decoder chemical foundation model for representing 3D electron density grids that efficiently encodes high-dimensional data into compact latent representations. This approach could significantly reduce reliance on computationally intensive quantum chemical simulations, offering simulation-grade data derived directly from learned representations.|
 | SMILESDFT-CLIP [\[HuggingFace\]](https://huggingface.co/ibm-research/materials.smilesdft-clip) | SMILESDFT-CLIP is a multimodal foundation model that jointly train a continuous 3D-field encoder, based on a vector-quantized generative adversarial network (VQGAN), and a SMILES-based transformer encoder on a dataset of 855,000 molecules, each represented by a DFT-computed electron density grid and a corresponding canonical SMILES string.|
 | STR-Bamba [\[HuggingFace\]](https://huggingface.co/ibm-research/materials.str-bamba) | STR-Bamba is an encoder-decoder chemical foundation model based on the IBM Bamba architecture, a hybrid of Transformers and Mamba-2 layers, designed to support multi-representational molecular string inputs.|
+| TDiMS [\[Code\]](models/tdims) | TDiMS (Topological Distance of intraMolecular Substructures) is a molecular descriptor that captures non-local intramolecular interactions by summarizing enumerated pairwise topological distances between molecular substructures. Unlike the pre-trained models above, TDiMS is computed directly from SMILES and requires no pre-training, making it lightweight and interpretable for small-data property prediction.|
 
 
 

models/tdims/examples/example_notebook.ipynb

@@ -879,16 +879,19 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from experiment_sample import main\n",
+    "exp_path = (Path.cwd().resolve().parent / \"experiments\").resolve()\n",
+    "if str(exp_path) not in sys.path:\n",
+    "    sys.path.insert(0, str(exp_path))\n",
+    "from run_nested_cv_experiment import main\n",
     "\n",
     "main(\n",
-    "    database=\"./data/cmpCl3_200.csv\",\n",
+    "    database=\"../data/cmpCl3_200.csv\",\n",
     "    prop=\"def_EmAbs\",\n",
     "    desc_name=\"TDiMS\",\n",
     "    outer_random_state=0,\n",
     "    outer_n_repeats=1,\n",
     "    n_jobs=-1,\n",
-    "    out_dir_sub=\"experimentcode\",\n",
+    "    out_dir=\"result\",\n",
     "    etc=\"notebook_example\",\n",
     "    save_cv_results=True,\n",
     "    save_joblib=False,\n",
@@ -900,9 +903,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python (TDiMS_Nature_v6)",
+   "display_name": "Python 3",
    "language": "python",
-   "name": "tdims_nature_v6"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {

models/tdims/requirements-notebook.txt

@@ -0,0 +1,27 @@
+# Notebook dependencies for examples/example_notebook.ipynb
+#
+# This file EXTENDS requirements.txt (installed via the -r line below),
+# so you do not need to install requirements.txt separately.
+#
+# Recommended setup (see README.md):
+#   conda create -n NCS python=3.10 -y
+#   conda activate NCS
+#   conda install -c conda-forge numba llvmlite shap
+#   pip install -r requirements-notebook.txt
+#
+# Installing shap from conda-forge first avoids numba / llvmlite build
+# failures on some platforms. If shap is already present, pip keeps the
+# existing (conda) build instead of rebuilding it.
+
+-r requirements.txt
+
+# Jupyter environment
+jupyterlab>=4.0,<5.0
+notebook>=7.0,<8.0
+ipykernel>=6.20,<7.0
+ipywidgets>=8.0,<9.0
+
+# Feature importance analysis (SHAP section of the notebook)
+shap>=0.44,<0.46
+numba>=0.57,<0.61
+llvmlite>=0.40,<0.44

models/tdims/src/tdims/load.py

@@ -12,7 +12,6 @@
 
 import pandas as pd
 from collections import defaultdict
-import types
 
 import warnings
 warnings.simplefilter('ignore')
@@ -32,7 +31,10 @@ def __init__(self, sm_list, radius=1, func_dis=-2, func_merge=sum, fragment_set=
             radius (int, optional): radius of finger print. Default to 1.
             func_dis (int, float, function, optional) : Calculation method for computing the feature value from bonds distance. Default to -2 (inverse square)
             func_merge (function, optional) :  Calculation method for merging feature values of distance in the same set of substructures. Default to sum.
-            fragment_set (bool, optional) : True if you want to include this substructure type to extract the distance. Default to True.
+            fragment_set (bool or list, optional):
+                If True, use the default CEP fragment list.
+                If False, do not use fragment-based features.
+                If a list/tuple/set of SMILES is given, use those fragments instead.
             atom_set (bool, optional) : Types True if you want to include this substructure type to extract the distance. Default to True
             fingerprint_set (bool, optional) : True if you want to include this substructure type to extract the distance. Default to True
         """
@@ -45,6 +47,7 @@ def __init__(self, sm_list, radius=1, func_dis=-2, func_merge=sum, fragment_set=
         self.atom_set = atom_set
         self.fingerprint_set = fingerprint_set
         self.nBit = nBit
+        self.failed_smiles = []
 
     def mfp_subset(self, mol):
 
@@ -71,52 +74,69 @@ def mfp_subset(self, mol):
         return mfp_set
 
     def extract_mol_features(self, sm):
-        
+
         topological_distance = dict()
 
         mol = Chem.MolFromSmiles(sm)
         if mol is None:
-            print(f"RDKit failed to read SMILES: {sm}")
+            logger.debug(f"RDKit failed to read SMILES: {sm}")
+            self.failed_smiles.append(sm)
             return topological_distance
 
         all_dic = dict()
         mol_distance_matrix = Chem.rdmolops.GetDistanceMatrix(mol)
         eps = 1.0e-10
 
-         # atom index of Fragment
-        if self.fragment_set:
-            sm_list = ['C1C=CC=C1',
-                        'S1N=C2C=CN=CC2=N1',
-                        'S1N=C2C=CC=CC2=N1',
-                        'O1C=CC2=CSC=C12',
-                        'C1=CC=NC=C1',
-                        'C1=NC=NC=N1',
-                        'O1C=CC=C1',
-                        '[SiH2]1C=CC=C1',
-                        '[SiH2]1C=C2C=CC=CC2=C1',
-                        'N1C=CC2=CSC=C12',
-                        'C1=CC2=CC=CC=C2C=C1',
-                        'S1C=CC=C1',
-                        'S1C=C2SC=CC2=C1',
-                        'O1C=C2C=CC=CC2=C1',
-                        '[Se]1C=CC=C1',
-                        'N1C=C2C=CC=CC2=C1',
-                        'N1C=CC=C1',
-                        'S1C=C2C(=C1)C1=CC=CC=C1C1=CC=CC=C21',
-                        'C1=CC=CC=C1',
-                        'S1C=C2N=CC=NC2=C1',
-                        'S1C=CN=C1',
-                        '[SiH2]1C=CC2=CSC=C12',
-                        'S1C=C2C=CC=CC2=C1',
-                        'S1C=C2[Se]C=CC2=C1',
-                        'C1C=C2C=CC=CC2=C1',
-                        'C1C=CC2=CSC=C12']
-
+        # atom index of Fragment
+        cep_fragment_smiles = [
+            'C1C=CC=C1',
+            'S1N=C2C=CN=CC2=N1',
+            'S1N=C2C=CC=CC2=N1',
+            'O1C=CC2=CSC=C12',
+            'C1=CC=NC=C1',
+            'C1=NC=NC=N1',
+            'O1C=CC=C1',
+            '[SiH2]1C=CC=C1',
+            '[SiH2]1C=C2C=CC=CC2=C1',
+            'N1C=CC2=CSC=C12',
+            'C1=CC2=CC=CC=C2C=C1',
+            'S1C=CC=C1',
+            'S1C=C2SC=CC2=C1',
+            'O1C=C2C=CC=CC2=C1',
+            '[Se]1C=CC=C1',
+            'N1C=C2C=CC=CC2=C1',
+            'N1C=CC=C1',
+            'S1C=C2C(=C1)C1=CC=CC=C1C1=CC=CC=C21',
+            'C1=CC=CC=C1',
+            'S1C=C2N=CC=NC2=C1',
+            'S1C=CN=C1',
+            '[SiH2]1C=CC2=CSC=C12',
+            'S1C=C2C=CC=CC2=C1',
+            'S1C=C2[Se]C=CC2=C1',
+            'C1C=C2C=CC=CC2=C1',
+            'C1C=CC2=CSC=C12'
+        ]
+
+        fragment_smiles_list = None
+        if self.fragment_set is True:
+            fragment_smiles_list = cep_fragment_smiles
+        elif isinstance(self.fragment_set, (list, tuple, set)):
+            fragment_smiles_list = list(self.fragment_set)
+
+        if fragment_smiles_list:
             fragment_atomidx_dic = defaultdict(list)
-            for smiles in sm_list:
-                atm_set=mol.GetSubstructMatches(Chem.MolFromSmiles(smiles))
+            for frag_smiles in fragment_smiles_list:
+                frag_mol = Chem.MolFromSmiles(frag_smiles)
+                if frag_mol is None:
+                    logger.warning(f"Invalid fragment SMILES was skipped: {frag_smiles}")
+                    continue
+
+                atm_set = mol.GetSubstructMatches(frag_mol)
                 if atm_set != ():
-                    fragment_atomidx_dic[f'{smiles}_CEPfrag']=[list(atm_idx) for atm_idx in atm_set]
+                    suffix = "CEPfrag" if self.fragment_set is True else "frag"
+                    fragment_atomidx_dic[f"{frag_smiles}_{suffix}"] = [
+                        list(atm_idx) for atm_idx in atm_set
+                    ]
             all_dic.update(fragment_atomidx_dic)
 
         # atom index of HeavyAtom
@@ -170,10 +190,10 @@ def extract_mol_features(self, sm):
 
             all_dic.update(sub_atomidx_dic_slc)
 
-
         dis_dic = defaultdict(list)
 
         for (sm1, sm2) in itertools.combinations_with_replacement(sorted(all_dic.keys(), key=len, reverse=True), 2):
+
             sub_pair = f'{sm1} & {sm2}'
 
             # collect each pair of substracture distance
@@ -188,9 +208,8 @@ def extract_mol_features(self, sm):
                                 distance_tmp.append(mol_distance_matrix[x][y])
                         av_dis = sum(distance_tmp)/len(distance_tmp)
 
-                        if isinstance(self.func_dis, types.FunctionType):
+                        if callable(self.func_dis):
                             calc_dis = self.func_dis(av_dis)
-
                         else:
                             calc_dis = av_dis**self.func_dis
 
@@ -205,24 +224,20 @@ def extract_mol_features(self, sm):
                                 distance_tmp.append(mol_distance_matrix[x][y])
                         av_dis = sum(distance_tmp) / len(distance_tmp)
 
-                        if isinstance(self.func_dis,types.FunctionType):
+                        if callable(self.func_dis):
                             calc_dis = self.func_dis(av_dis)
-
                         else:
                             calc_dis = av_dis ** self.func_dis
 
                         dis_dic[sub_pair].append(calc_dis)
-
 
             if dis_dic[sub_pair] == []:
                 distance_final = 0
             else:
                 distance_final = self.func_merge(dis_dic[sub_pair])
 
-
             if distance_final-eps > 0:
                 topological_distance[sub_pair] = distance_final
-
             elif distance_final < 0:
                 logger.error(f'Minus feature for mol:{Chem.MolToSmiles(mol)} feature:{sub_pair}')
 
@@ -232,11 +247,18 @@ def feature_extraction(self):
 
         key_all={}
         dic_all=[]
+        self.failed_smiles = []
         for sm in self.sm_list:
             tdims_dic = self.extract_mol_features(sm)
             dic_all.append(tdims_dic)
             key_all = {**key_all, **tdims_dic}
-
+
+        if self.failed_smiles:
+            logger.warning(
+                f"{len(self.failed_smiles)} of {len(dic_all)} SMILES could not be parsed by RDKit "
+                f"and produced all-zero feature rows. See TDiMS.failed_smiles for the list."
+            )
+
         X=[]
         for dic in dic_all:
             x = [dic.get(key, 0) for key in key_all]

models/tdims/src/tdims/sparse_transformers.py

@@ -12,16 +12,6 @@
 
 
 class SparseFeatureFilterVCSR(BaseEstimator, TransformerMixin):
-    """
-    非ゼロサンプル数が少ない特徴量を除去（dense/sparse 両対応）
-
-    Parameters
-    ----------
-    min_samples : int
-        各特徴量が非ゼロとなるサンプル数が min_samples 未満なら除去
-    verbose : bool
-        fit時に特徴量数の変化をprintするか
-    """
     def __init__(self, min_samples: int = 30, verbose: bool = True):
         self.min_samples = min_samples
         self.verbose = verbose
@@ -52,9 +42,6 @@ def get_support(self, indices: bool = False):
 
 
 class ToCSR(BaseEstimator, TransformerMixin):
-    """
-    入力を CSR sparse matrix に変換（すでに sparse なら .tocsr()）
-    """
     def fit(self, X, y=None):
         return self
 
@@ -64,7 +51,6 @@ def transform(self, X):
         return csr_matrix(X)
 
 class ClipGreaterThanOneToZero(BaseEstimator, TransformerMixin):
-    """CSR（疎行列）を想定し、値が threshold より大きい要素を 0 にする。"""
     def __init__(self, threshold: float = 1.0):
         self.threshold = threshold
 
@@ -74,16 +60,14 @@ def fit(self, X, y=None):
     def transform(self, X):
         # sparse
         if issparse(X):
-            X = X.tocsr(copy=True)  # 元を壊さない
+            X = X.tocsr(copy=True) 
             if X.nnz == 0:
                 return X
             mask = X.data > self.threshold
             if np.any(mask):
                 X.data[mask] = 0.0
                 X.eliminate_zeros()
             return X
-
-        # denseが来た場合の保険（通常ここには来ない想定）
         X = np.array(X, copy=True)
         X[X > self.threshold] = 0.0
         return csr_matrix(X)

models/tdims/src/tdims/tdims_ext.py

@@ -51,6 +51,9 @@ def get_representation_with_fs_selection(sm_list, y, model="tdims", radius=1, fu
 
         return x_slc, list(key_slc), key_all
 
+    else:
+        raise ValueError("Invalid model input")
+
 
 
 def run_tdims_regression_cv(