add 1dqsar

1d-qsar/cuda/qssar1d.py

@@ -0,0 +1,152 @@
+#!/usr/bin/env python
+# -*- encoding: utf-8 -*-
+'''
+@file: RFE_cuml_permutation.py
+@Description: 使用GPU cuML实现特征选择，手动实现RFECV，并通过置换重要性计算来剔除最不重要的特征
+@Date: 2025/03/01
+@Author: lyzeng
+'''
+
+import numpy as np
+import pandas as pd
+
+# ---------------------------
+# 数据加载与特征构造
+# ---------------------------
+data = pd.read_csv("../data_smi.csv", index_col="Entry ID")
+target_data = data[["SMILES", "S.aureus ATCC25923"]]
+target_data.columns = ["SMILES", "TARGET"]
+
+from rdkit import Chem
+from rdkit.Chem import Descriptors
+
+def getMolDescriptors(smiles: str, missingVal=None):
+    mol = Chem.MolFromSmiles(smiles)
+    res = {}
+    for nm, fn in Descriptors._descList:
+        try:
+            val = fn(mol)
+        except Exception as e:
+            val = missingVal
+        res[nm] = val
+    return res
+
+# 构建特征矩阵与目标变量
+X_df = pd.DataFrame([getMolDescriptors(smi) for smi in target_data['SMILES']])
+y = target_data['TARGET'].values
+
+# 剔除无效特征（取值唯一的列）
+invalid_features = [col for col in X_df.columns if X_df[col].nunique() <= 1]
+X_df.drop(columns=invalid_features, inplace=True)
+X = X_df.values
+
+print(f"训练样本数: {len(y)}, 特征数量: {X.shape[1]}")
+
+# ---------------------------
+# 转换为GPU数据格式（cupy数组）
+# ---------------------------
+import cupy as cp
+X_gpu = cp.asarray(X)
+y_gpu = cp.asarray(y)
+
+# ---------------------------
+# 定义辅助函数：计算模型的MSE
+# ---------------------------
+def model_mse(model, X_data, y_data):
+    y_pred = model.predict(X_data)
+    mse = cp.mean((y_data - y_pred) ** 2)
+    return mse.get()
+
+# ---------------------------
+# 定义置换重要性计算函数（在GPU上）
+# ---------------------------
+def permutation_importance_gpu(model, X_data, y_data, n_repeats=3, random_state=0):
+    """
+    对于给定的模型，计算每个特征的置换重要性：
+      重要性 = (打乱特征后的MSE均值) - (原始MSE)
+    随机打乱每个特征的值，并观察模型性能的下降程度来衡量特征的重要性。 如果打乱某个特征导致模型性能显著下降，说明该特征对模型很重要；反之，如果打乱某个特征对模型性能没有明显影响，说明该特征不太重要。
+    """
+    # 先将数据复制到CPU进行打乱操作（数据量较小时不会有太大开销）
+    X_cpu = cp.asnumpy(X_data)
+    baseline = model_mse(model, X_data, y_data)
+    importances = np.zeros(X_cpu.shape[1])
+    rng = np.random.RandomState(random_state)
+    for j in range(X_cpu.shape[1]):
+        scores = []
+        for _ in range(n_repeats):
+            X_permuted = X_cpu.copy()
+            rng.shuffle(X_permuted[:, j])
+            X_permuted_gpu = cp.asarray(X_permuted)
+            score = model_mse(model, X_permuted_gpu, y_data)
+            scores.append(score)
+        importances[j] = np.mean(scores) - baseline
+    return importances
+
+# ---------------------------
+# 使用 cuML 随机森林和置换重要性手动实现 RFECV
+# ---------------------------
+from cuml.ensemble import RandomForestRegressor as cuRF
+from sklearn.model_selection import KFold
+
+def cross_val_score_gpu(X_data, y_data, cv=5):
+    """
+    使用 KFold 分割数据，利用 cuML 随机森林评估当前特征子集的CV均方误差
+    """
+    kf = KFold(n_splits=cv, shuffle=True, random_state=0)
+    mse_scores = []
+    # 使用 CPU 索引进行分割
+    X_np = cp.asnumpy(X_data)
+    for train_index, test_index in kf.split(X_np):
+        X_train = X_data[train_index, :]
+        X_test = X_data[test_index, :]
+        y_train = y_data[train_index]
+        y_test = y_data[test_index]
+        model = cuRF(n_estimators=100, random_state=0)
+        model.fit(X_train, y_train)
+        mse = model_mse(model, X_test, y_test)
+        mse_scores.append(mse)
+    return np.mean(mse_scores)
+
+def manual_rfecv(X_data, y_data, cv=5):
+    """
+    手动实现 RFECV:
+      - 在当前特征子集上计算CV MSE
+      - 使用置换重要性评估各特征的重要性
+      - 移除置换重要性最低的特征
+      - 记录使CV MSE最小的特征组合
+    找到一个最佳的特征组合，用于训练最终的模型。 这个子集的大小在每次迭代中都在变化，最终目标是找到一个使交叉验证 MSE 最小的子集
+     cuml.ensemble.RandomForestRegressor 的 max_features 默认值通常是 sqrt (特征总数的平方根)。 这与 scikit-learn 的 RandomForestRegressor 的默认值相同。
+    """
+    n_features = X_data.shape[1]
+    features = list(range(n_features))
+    best_score = float('inf')
+    best_features = features.copy()
+    scores_history = []
+    
+    while len(features) > 0:
+        current_score = cross_val_score_gpu(X_data[:, features], y_data, cv=cv)
+        scores_history.append((features.copy(), current_score))
+        print(f"当前特征数: {len(features)}, CV MSE: {current_score:.4f}")
+        if current_score < best_score:
+            best_score = current_score
+            best_features = features.copy()
+        if len(features) == 1:
+            break
+        # 训练模型并计算置换重要性
+        model = cuRF(n_estimators=100, random_state=0) #  cuml.ensemble.RandomForestRegressor 的 max_features 默认值通常是 sqrt (特征总数的平方根)。 这与 scikit-learn 的 RandomForestRegressor 的默认值相同。
+        model.fit(X_data[:, features], y_data)
+        imp = permutation_importance_gpu(model, X_data[:, features], y_data, n_repeats=3, random_state=0)
+        idx_to_remove = int(np.argmin(imp))
+        removed_feature = features[idx_to_remove]
+        print(f"剔除特征索引: {removed_feature}，置换重要性: {imp[idx_to_remove]:.4f}")
+        del features[idx_to_remove]
+    
+    return best_features, best_score, scores_history
+
+# 执行手动 RFECV，使用5折交叉验证
+cv_folds = 5
+best_features_rfecv, best_mse_rfecv, history = manual_rfecv(X_gpu, y_gpu, cv=cv_folds)
+
+print(f"\n手动RFECV选择了 {len(best_features_rfecv)} 个特征，CV MSE: {best_mse_rfecv:.4f}")
+selected_feature_names = [X_df.columns[i] for i in best_features_rfecv]
+print("选定特征名称：", selected_feature_names)