Comprehensive_Molecular_Featurizer_V2.0   .ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CMF  \n",
    "## 分子综合化特征器  \n",
    "## Comprehensive Molecular Featurizer   \n",
    "Writed by Yiheng Dai and Yuming Su\n",
    "### 目标：  \n",
    "综合多种特征化方式，搭建一体式的化学分子特征化工具  \n",
    "将容纳多种由后续探索加入的描述符  \n",
    "### 方法简介：  \n",
    "内含ECFP(改)方法，版本为6.0\n",
    "内含RDKit描述符方法\n",
    "内含自创的共轭结构描述符方法(92个描述符，含5套较复杂的原子级描述符)\n",
    "内含溶剂与波长合并模块 \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:09:04.610382Z",
     "start_time": "2022-06-27T13:09:04.602904Z"
    }
   },
   "outputs": [],
   "source": [
    "# 特征器整体参数：\n",
    "VERSION = '2.02sym'\n",
    "import time\n",
    "# 输入文件：应当包含两列，第一列标题为'smiles'，第二列标题任意，为标签值\n",
    "INPUT_NAME = 'TPA_856_0307.csv'\n",
    "LOG_NAME = 'Log_CMF_V'+VERSION+'_'+time.strftime(\"%Y-%m-%d_%H-%M-%S\", time.localtime())+'.txt'\n",
    "SMILES_CHECK = True                       # 进行SMILES检查的开关，开启后若存在错误的smiles，将会终止运行\n",
    "VALUES_DIV_MOLWT = True                   # 将标签值除以摩尔质量的开关\n",
    "VALUES_LN = True                          # 将标签值取对数值的开关\n",
    "PANDAS_DATASET_GENERATE = True            # 生成Pandas可用的数据集的开关"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:09:04.705030Z",
     "start_time": "2022-06-27T13:09:04.613395Z"
    }
   },
   "outputs": [],
   "source": [
    "# =====================ECFP part=====================\n",
    "\n",
    "# ECFP特征器参数：\n",
    "ECFP_SWITCH = True                        # ECFP特征器的总开关\n",
    "# 输出文件：一共6个：X，y，特征列标题Title，无标题的SMILES，两个整合的数据集\n",
    "RADIUS = 4                                # 整数，设置ECFP提取特征的半径\n",
    "\n",
    "# V1.5加入'特征列填充率保留模块'：\n",
    "FEATURE_RESERVE = 0.015                  # 浮点数，特征矩阵内非零数据的比例大于等于这个一般保留阈值，这个特征就会被保留\n",
    "EXACT_PIECE_RESERVE = False                # 布尔值，表示是否开启使用特定基团保留片段的功能\n",
    "SMILES_PATT = []                      # 字符串列表，每个元素都是一段SMILES值(其实用的是SMART，这样可以匹配类似'ccc'的芳环残片)，用于在保留时作为匹配模式(仅在EXACT_PIECE_RESERVE为True时生效)\n",
    "EXACT_PIECE_RESERVE_THRESHOLD = 0.015    # 浮点数，针对特定基团的阈值，和上面的类似，若想要全部保留含有某特定基团的特征，此项请用0.0(仅在EXACT_PIECE_RESERVE为True时生效)\n",
    "# 假如只想保留含有特定基团的片段，可以把FEATURE_RESERVE设为1.0，而EXACT_PIECE_RESERVE设为0.0，因为特定基团筛选的优先级高于一般保留阈值的筛选\n",
    "\n",
    "# V2.0为针对'CO'的处理而加入，V5.0修改为更具有普遍性的'布尔值转化模块'：\n",
    "# 作用：将一个max值为n的特征列转化为n或(n+1)列布尔值特征列，分别表示拥有(0),1,2,3,...,n个特征片段\n",
    "# 注意：若进行全部数据的转化以求获得一个只由{0, 1}构成的特征矩阵，不应和加入RDKit描述符的功能一起使用\n",
    "# V5.1：添加了自动删除空特征列的功能\n",
    "# V6.0：整合了两种转化为Bool值的模式，即a版(n列)与b版(n+1列)\n",
    "TURN_TO_BOOL = False                      # 布尔值，作为布尔值转化模块的总开关使用\n",
    "ALL_TURN_TO_BOOL = True                   # 布尔值，若为True，则将所有非{0, 1}整数数据全部转化为多个布尔值特征列，最后的特征值矩阵将完全由{0, 1}构成\n",
    "BOOL_PATT_LIST = []                       # 字符串列表，仅在ALL_TURN_TO_BOOL为False时生效，表示将列表中指定的基团转变为布尔值特征内(注意：必须特征列的标题和这里面的字符串一样才进行转换操作)\n",
    "BOOL_TURN_MODE = 'a'                      # 字符，'a'或'b'，代表两种布尔值转化的模式，a为n列，b为n+1列\n",
    "\n",
    "# V3.0加入'RDKit描述符合并模块'，现独立出来作为另一种特征化方法处理\n",
    "# INCLUDE_DESC = False                      # 布尔值，表示是否需要计算并加入RDKit的描述符\n",
    "\n",
    "# V4.0加入'片段原子数量调控模块'：\n",
    "ATOM_COUNT_CONTROL = True                 # 布尔值，作为“原子数量调控模块”的总开关\n",
    "SMARTS_MIN_LENGTH = 2                     # 整数值，表示特征SMARTS内非氢原子的最小数目，包含非氢原子数目大于等于该数目的基团片段才会被保留，设为0就会跳过这一步骤\n",
    "SMARTS_MAX_LENGTH = 50                    # 整数值，表示特征SMARTS内非氢原子的最小数目，包含非氢原子数目小于等于该数目的基团片段才会被保留，若不想在最大值上设限，可设为5*RADIUS或更大的数值\n",
    "ATOM_COUNT_CONTROL_OMIT_PATTERN = []      # 字符串列表，表示若特征片段内含有此列表中的SMARTS片段，此列将不受“原子数量调控模块”的调控\n",
    "\n",
    "# V4.1加入：新输出一个带Title和Value的完整的数据矩阵\n",
    "\n",
    "SIMILARITY_SWITCH = False                  # 布尔值，作为“相似度特征化模块”的总开关，使用时请确保TURN_TO_BOOL为False\n",
    "SIMILARITY_METHOD = 'Tanimoto'            # 字符串，'Tanimoto'或'MACCS'，作为相似度计算的方法\n",
    "SIMILARITY_MODE = 'mean'                  # 字符串，'max'或'mean'，作为标签值为0的特征的转化方法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:09:04.852698Z",
     "start_time": "2022-06-27T13:09:04.708429Z"
    }
   },
   "outputs": [],
   "source": [
    "# =====================RDKit Descriptor part=====================\n",
    "\n",
    "# RDKit描述符特征器参数：\n",
    "# 使用前请先使用Smiles检查程序对数据包中的smiles进行检查，确保\n",
    "RDKIT_DESC_SWITCH = True           # RDKit描述符特征器的总开关\n",
    "allowedDescriptors = ['NOCount', 'VSA_EState4', 'NumHDonors', 'SlogP_VSA12', 'NumRadicalElectrons', \n",
    "                      'SlogP_VSA4', 'Kappa2', 'Chi2n', 'PEOE_VSA3', 'PEOE_VSA7', 'PEOE_VSA4', 'Chi1', \n",
    "                      'MolWt', 'SMR_VSA1', 'SlogP_VSA9', 'VSA_EState9', 'MaxAbsPartialCharge', 'NumSaturatedHeterocycles', \n",
    "                      'MaxPartialCharge', 'VSA_EState1', 'PEOE_VSA6', 'EState_VSA11', 'SMR_VSA4', 'EState_VSA7', \n",
    "                      'VSA_EState2', 'NHOHCount', 'SlogP_VSA10', 'SMR_VSA7', 'PEOE_VSA9', 'NumAliphaticRings', 'EState_VSA8', \n",
    "                      'PEOE_VSA5', 'BertzCT', 'SlogP_VSA1', 'SlogP_VSA6', 'PEOE_VSA1', 'VSA_EState7', 'MinAbsPartialCharge', \n",
    "                      'LabuteASA', 'SlogP_VSA2', 'EState_VSA4', 'MolMR', 'Kappa1', 'NumHAcceptors', 'EState_VSA9', 'MolLogP', \n",
    "                      'NumAromaticHeterocycles', 'BalabanJ', 'FractionCSP3', 'SMR_VSA3', 'RingCount', 'NumSaturatedRings', \n",
    "                      'PEOE_VSA2', 'MaxAbsEStateIndex', 'Kappa3', 'Chi3n', 'NumRotatableBonds', 'Chi4n', 'VSA_EState3', \n",
    "                      'SMR_VSA8', 'MinPartialCharge', 'EState_VSA6', 'SMR_VSA9', 'PEOE_VSA13', 'NumValenceElectrons', \n",
    "                      'MaxEStateIndex', 'SMR_VSA6', 'VSA_EState8', 'EState_VSA2', 'NumAromaticCarbocycles', 'SMR_VSA10', \n",
    "                      'SlogP_VSA3', 'HallKierAlpha', 'PEOE_VSA14', 'HeavyAtomCount', 'VSA_EState10', 'SlogP_VSA11', \n",
    "                      'ExactMolWt', 'MinAbsEStateIndex', 'TPSA', 'PEOE_VSA10', 'SMR_VSA2', 'Chi1v', 'Chi4v', 'PEOE_VSA8', \n",
    "                      'EState_VSA5', 'Chi1n', 'VSA_EState5', 'SlogP_VSA7', 'HeavyAtomMolWt', 'MinEStateIndex', \n",
    "                      'NumAliphaticHeterocycles', 'VSA_EState6', 'Chi0v', 'SlogP_VSA5', 'SMR_VSA5', 'Chi0', 'Chi2v', \n",
    "                      'NumSaturatedCarbocycles', 'NumAromaticRings', 'Chi0n', 'PEOE_VSA12', 'Chi3v', 'NumAliphaticCarbocycles', \n",
    "                      'EState_VSA10', 'EState_VSA3', 'EState_VSA1', 'NumHeteroatoms', 'SlogP_VSA8', 'PEOE_VSA11']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:09:05.011577Z",
     "start_time": "2022-06-27T13:09:04.856095Z"
    }
   },
   "outputs": [],
   "source": [
    "# =====================Conjugation Descriptor part=====================\n",
    "\n",
    "# 共轭特征器模块参数：\n",
    "CONJU_DESC_SWITCH = True                  # 共轭描述符特征器的总开关\n",
    "KEEP_TYPE = 'max'           # 保留方式，设为'max'、'mean'或'acc-mean'\n",
    "patt_list_d = ['C=C', 'C#C', 'C#N', 'C=O', 'C=S', 'C=N', 'N=N', '[N+]([O-])=O']\n",
    "patt_list_m = ['N', 'O', 'S', 'F', 'Cl', 'Br', 'I', 'P']\n",
    "one_list = ['C']\n",
    "two_list = ['N', 'O', 'S', 'P', 'F', 'Cl', 'Br', 'I']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:09:05.111954Z",
     "start_time": "2022-06-27T13:09:05.031924Z"
    }
   },
   "outputs": [],
   "source": [
    "# 附加溶剂模块参数：\n",
    "SOLVENT_SWITCH = True                  # 附加溶剂\n",
    "SOLVENT_IN = 'Solvent_856.csv'\n",
    "SOLVENT_TITLE = ['ET(30) (Solvent)', 'Dielectic Constant (Solvent)', 'Dipole Moment (Solvent)']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:09:05.220836Z",
     "start_time": "2022-06-27T13:09:05.131262Z"
    }
   },
   "outputs": [],
   "source": [
    "# 附加波长模块参数：\n",
    "WAVE_SWITCH = True                  # 附加波长\n",
    "WAVE_IN = 'Wave_856.csv'\n",
    "WAVE_TITLE = ['Wavelength (Exp nm)']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:09:10.019488Z",
     "start_time": "2022-06-27T13:09:05.225596Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/jyb/.conda/envs/dyhpy/lib/python3.7/site-packages/requests/__init__.py:91: RequestsDependencyWarning: urllib3 (1.26.4) or chardet (4.0.0) doesn't match a supported version!\n",
      "  RequestsDependencyWarning)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import deepchem as dc\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem import AllChem\n",
    "from rdkit.Chem import Descriptors\n",
    "from rdkit.Chem import rdMolDescriptors\n",
    "from rdkit import DataStructs\n",
    "from rdkit.Chem import MACCSkeys\n",
    "from rdkit.Chem.EState.EState import EStateIndices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:09:10.030346Z",
     "start_time": "2022-06-27T13:09:10.024552Z"
    }
   },
   "outputs": [],
   "source": [
    "# 创建一个目录来保存生成的数据包\n",
    "import os\n",
    "from pathlib import Path\n",
    "DIR = 'CMF_'+VERSION+'_'+time.strftime(\"%Y-%m-%d_%H-%M-%S\", time.localtime())\n",
    "os.mkdir(DIR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:09:10.155910Z",
     "start_time": "2022-06-27T13:09:10.033381Z"
    }
   },
   "outputs": [],
   "source": [
    "# 输出数据的创建\n",
    "X_out = []\n",
    "title_out = []\n",
    "smiles_out = []\n",
    "values_out = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:09:10.325072Z",
     "start_time": "2022-06-27T13:09:10.158555Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "deepchem.utils.save has been deprecated.\n",
      "The utilities in save.py are moved to deepchem.utils.data_utils or deepchem.utils.genomics_utils.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Columns of dataset: ['smiles' 'values']\n",
      "Number of examples in dataset: 856\n"
     ]
    }
   ],
   "source": [
    "# 数据集的读取，供DeepChem相关的特征器使用\n",
    "from deepchem.utils.save import load_from_disk\n",
    "dataset_file= INPUT_NAME\n",
    "dataset = load_from_disk(dataset_file)\n",
    "print(\"Columns of dataset: %s\" % str(dataset.columns.values))\n",
    "print(\"Number of examples in dataset: %s\" % str(dataset.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:09:16.996851Z",
     "start_time": "2022-06-27T13:09:10.327545Z"
    }
   },
   "outputs": [],
   "source": [
    "# 数据集的读取，供RDKit相关的特征器使用\n",
    "newdata = []\n",
    "with open(INPUT_NAME) as file:\n",
    "    data = file.readlines()\n",
    "del data[0]\n",
    "for i in range(len(data)):\n",
    "    newdata.append(data[i].split(','))\n",
    "    newdata[i][1] = eval(newdata[i][1])\n",
    "newdata = np.array(newdata)\n",
    "smiles_rd = newdata[:, 0].flatten().tolist()\n",
    "values_rd = newdata[:, 1].astype(float).flatten().tolist()\n",
    "# 检查Smiles\n",
    "if SMILES_CHECK:\n",
    "    error_list = []\n",
    "    for smi in smiles_rd:\n",
    "        try:\n",
    "            mol = Chem.MolFromSmiles(smi)\n",
    "        except:\n",
    "            print('error:', smi)\n",
    "# 将标签值除以分子质量\n",
    "if VALUES_DIV_MOLWT:\n",
    "    values_rd_dwt = []\n",
    "    for descriptor, function in Descriptors.descList:\n",
    "        if descriptor=='MolWt':\n",
    "            for i in range(len(smiles_rd)):\n",
    "                smile = smiles_rd[i]\n",
    "                mol = Chem.MolFromSmiles(str(smile))\n",
    "                v = values_rd[i]/function(mol)\n",
    "                if v==0:\n",
    "                    print(values_rd[i], function(mol))\n",
    "                values_rd_dwt.append(v)\n",
    "    values_dwt_out = np.array(values_rd_dwt).reshape(len(values_rd_dwt), 1)\n",
    "smiles_out = np.array(smiles_rd).reshape(len(smiles_rd), 1)\n",
    "values_out = np.array(values_rd).reshape(len(values_rd), 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:10:06.043223Z",
     "start_time": "2022-06-27T13:09:17.000492Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "smiles_field is deprecated and will be removed in a future version of DeepChem.Use feature_field instead.\n",
      "/home/jyb/.conda/envs/dyhpy/lib/python3.7/site-packages/deepchem/data/data_loader.py:162: FutureWarning: featurize() is deprecated and has been renamed to create_dataset().featurize() will be removed in DeepChem 3.0\n",
      "  \"featurize() will be removed in DeepChem 3.0\", FutureWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sparse:\n",
      "((856,), (856, 1), (856, 1), (856,))\n",
      "(856,)\n",
      "unique hash: 65272\n",
      "unique pieces 7845\n",
      "7844\n",
      "Before delete sparse coloumns: length of title: 7844\n",
      "Before delete sparse coloumns: shape of data: (856, 7844)\n"
     ]
    }
   ],
   "source": [
    "if ECFP_SWITCH:\n",
    "    # 定义‘特征器’(Featurizer)来将SMILES转化为某些特征\n",
    "    # featurizer_fix = dc.feat.CircularFingerprint(size=1024)\n",
    "    featurizer_sparse = dc.feat.CircularFingerprint(size=1024, sparse=True, smiles=True, radius=RADIUS)\n",
    "    # if INCLUDE_DESC:\n",
    "        # featurizer_desc = dc.feat.RDKitDescriptors()\n",
    "    # 对SMILES进行特征化，输出为哈希值以及SMILES块\n",
    "    # 因为hash值与smiles块一一对应，我们只需挑出unique的SMILES/hash作为特征即可\n",
    "    loader_sparse = dc.data.CSVLoader(\n",
    "          tasks=[dataset.columns.values[1]], smiles_field=dataset.columns.values[0],\n",
    "          featurizer=featurizer_sparse)\n",
    "    dataset_sparse = loader_sparse.featurize(dataset_file)\n",
    "    print('Sparse:')\n",
    "    print(dataset_sparse.get_shape())\n",
    "    # 此时的特征“矩阵”应当为一个大的字典\n",
    "    print(dataset_sparse.get_shard(0)[0].shape)\n",
    "    # 取出所有的分子，这一步得到一个原数据集长度[0]的列表，每个元素为1个字典，字典中为hash:{'smiles':'str', 'count':, int}\n",
    "    fps_sparse = list(dataset_sparse.get_shard(0)[0])\n",
    "    smiles = []\n",
    "    # count = []\n",
    "    for i in range(dataset_sparse.get_shard(0)[0].shape[0]):\n",
    "        # 这一步得到的依然是1个字典\n",
    "        fps_dict = fps_sparse[i]\n",
    "        for hash_value, sub_dict in fps_dict.items():\n",
    "            # 这一步把所有的smiles和count挂进列表，count会影响某种特征的权重，也有必要记录下来\n",
    "            smiles.append(sub_dict['smiles'])\n",
    "            # count.append(sub_dict['count'])\n",
    "    print('unique hash:', len(smiles))\n",
    "    # 这一步得到所有特异的smiles，发现1个smiles可能对应多个hash\n",
    "    smiles_unique = np.unique(smiles)\n",
    "    print('unique pieces', len(smiles_unique))\n",
    "    # 目标是把特征矩阵缩到smiles_unique的长度，这样每一列代表一种特征，特征值即为count的权重\n",
    "    # 如果这样，可以选择跳过hash，直接统计特征\n",
    "    title = []\n",
    "    data = np.zeros((dataset_sparse.get_shard(0)[0].shape[0], len(smiles_unique)-1))\n",
    "    idx = 0\n",
    "    index = 0\n",
    "    for i in range(dataset_sparse.get_shard(0)[0].shape[0]):\n",
    "        # 这一步得到的依然是1个字典\n",
    "        fps_dict = fps_sparse[i]\n",
    "        # 数据提炼去重：\n",
    "        \n",
    "        smi = dataset_sparse.get_shard(0)[3][index]\n",
    "        m = Chem.MolFromSmiles(smi)\n",
    "            \n",
    "        for hash_value, sub_dict in fps_dict.items():\n",
    "            SMILES = sub_dict['smiles']\n",
    "            patt = Chem.MolFromSmarts(SMILES)\n",
    "            atomids = m.GetSubstructMatches(patt)\n",
    "            COUNT = len(atomids)\n",
    "            # COUNT = sub_dict['count']\n",
    "            # 跳过'这个空SMILES\n",
    "            if SMILES=='':\n",
    "                continue\n",
    "            if SMILES not in title:\n",
    "                title.append(SMILES)\n",
    "                data[i][idx] += COUNT\n",
    "                idx +=1\n",
    "            else:\n",
    "                if data[i][title.index(SMILES)]==0:\n",
    "                    data[i][title.index(SMILES)] += COUNT\n",
    "        index += 1\n",
    "    print(idx)\n",
    "    print('Before delete sparse coloumns: length of title:', len(title))\n",
    "    print('Before delete sparse coloumns: shape of data:', data.shape)\n",
    "    data = data.astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:10:30.242252Z",
     "start_time": "2022-06-27T13:10:06.046229Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "After delete sparse coloumns: length of title: 564\n",
      "After delete sparse coloumns: shape of data: (856, 564)\n"
     ]
    }
   ],
   "source": [
    "if ECFP_SWITCH:\n",
    "    # 对于过于稀疏的特征列，设置FEATURE_RESERVE特征保留阈值进行筛选，列填充比例大于此的才会被保留\n",
    "    # 同时还要保留含有特定基团的片段，此操作的优先级高于按保留阈值筛选\n",
    "    delete = []\n",
    "    for i in range(data.shape[1]):\n",
    "        filled_rate = sum(data[:, i]!=0)/data.shape[0]\n",
    "        if filled_rate<FEATURE_RESERVE:\n",
    "            if (EXACT_PIECE_RESERVE):                               # 这一块是在挑出含有指定集团的片段\n",
    "                m = Chem.MolFromSmarts(title[i])\n",
    "                for patt_temp in SMILES_PATT:\n",
    "                    patt = Chem.MolFromSmarts(patt_temp)\n",
    "                    flag = m.HasSubstructMatch(patt)\n",
    "                    if flag:\n",
    "                        break\n",
    "                if flag and (filled_rate>=EXACT_PIECE_RESERVE_THRESHOLD):\n",
    "                    continue\n",
    "            delete.append(i)\n",
    "    data = np.delete(data, delete, axis=1)\n",
    "    title = np.array(title).reshape(1, len(title))\n",
    "    title = np.delete(title, delete, axis=1)\n",
    "    title = title.flatten().tolist()\n",
    "    print('After delete sparse coloumns: length of title:', len(title))\n",
    "    print('After delete sparse coloumns: shape of data:', data.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:10:30.289953Z",
     "start_time": "2022-06-27T13:10:30.245493Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "After atom-count control: length of title: 564\n",
      "After atom-count control: shape of data: (856, 564)\n"
     ]
    }
   ],
   "source": [
    "# 原子数量调控模块：\n",
    "if ECFP_SWITCH and ATOM_COUNT_CONTROL:\n",
    "    delete = []\n",
    "    for i in range(data.shape[1]):\n",
    "        flag = False\n",
    "        for patt_temp in ATOM_COUNT_CONTROL_OMIT_PATTERN:\n",
    "            patt = Chem.MolFromSmarts(patt_temp)\n",
    "            flag = m.HasSubstructMatch(patt)\n",
    "            if flag:\n",
    "                break\n",
    "        if flag:\n",
    "            continue\n",
    "        m = Chem.MolFromSmarts(title[i])\n",
    "        atom_count = len(m.GetAtoms())\n",
    "        if atom_count<SMARTS_MIN_LENGTH or atom_count>SMARTS_MAX_LENGTH:\n",
    "            delete.append(i)\n",
    "    data = np.delete(data, delete, axis=1)\n",
    "    title = np.array(title).reshape(1, len(title))\n",
    "    title = np.delete(title, delete, axis=1)\n",
    "    title = title.flatten().tolist()\n",
    "    print('After atom-count control: length of title:', len(title))\n",
    "    print('After atom-count control: shape of data:', data.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:10:30.385080Z",
     "start_time": "2022-06-27T13:10:30.292777Z"
    }
   },
   "outputs": [],
   "source": [
    "# 这一块里，会删除某一指定的特征列，将其转化为表示含有0,1,...,n个基团的布尔表示列，其中n为所有样本中含基团的最大值，最后会自动把空列删除\n",
    "if ECFP_SWITCH and TURN_TO_BOOL:\n",
    "    if ALL_TURN_TO_BOOL:\n",
    "        BOOL_PATT_LIST = title\n",
    "    for bool_patt in BOOL_PATT_LIST:\n",
    "        print('Processing feature: ', bool_patt)\n",
    "        print('title length(before turn to bool):', len(title))\n",
    "        print('data shape(before turn to bool):', data.shape)\n",
    "        index_t = title.index(bool_patt)\n",
    "        print('Target conlumn index:', index_t)\n",
    "        feature_max = int(max(data[:, index_t]))\n",
    "        print('max number of target patt:', feature_max)\n",
    "        # 制作新的bool数据数组\n",
    "        if BOOL_TURN_MODE=='a':\n",
    "            matrix_bool = np.zeros((data.shape[0], feature_max))\n",
    "        else:\n",
    "            matrix_bool = np.zeros((data.shape[0], feature_max+1))\n",
    "        for i in range(data.shape[0]):\n",
    "            if BOOL_TURN_MODE=='a':\n",
    "                if data[i, index_t]>=1:\n",
    "                    matrix_bool[i, int(data[i, index_t]-1)] = 1\n",
    "            else:\n",
    "                matrix_bool[i, int(data[i, index_t])] = 1\n",
    "        # 删除对应的数据列和标题列\n",
    "        data = np.delete(data, [index_t], axis=1)\n",
    "        title = np.array(title).reshape(1, len(title))\n",
    "        title = np.delete(title, [index_t], axis=1)\n",
    "        title = title.flatten().tolist()\n",
    "        # 制作额外的标题\n",
    "        new_title = []\n",
    "        if BOOL_TURN_MODE=='a':\n",
    "            for i in range(feature_max):\n",
    "                new_title.append('With_'+str(i+1)+'_'+bool_patt)\n",
    "        else:\n",
    "            for i in range(feature_max+1):\n",
    "                new_title.append('With_'+str(i)+'_'+bool_patt)\n",
    "        data = np.hstack((data, matrix_bool))\n",
    "        title = title+new_title\n",
    "        print('title length(after turn to bool):', len(title))\n",
    "        print('data shape(after turn to bool):', data.shape, '\\n')\n",
    "    # 将空列自动删除，V5.1加入\n",
    "    delete = []\n",
    "    for i in range(data.shape[1]):\n",
    "        if sum(data[:, i])==0:\n",
    "            delete.append(i)\n",
    "    print(data.shape)\n",
    "    print(len(title))\n",
    "    data = np.delete(data, delete, axis=1)\n",
    "    title = np.array(title).reshape(1, len(title))\n",
    "    title = np.delete(title, delete, axis=1)\n",
    "    title = title.flatten().tolist()\n",
    "    print(data.shape)\n",
    "    print(len(title))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:10:30.517412Z",
     "start_time": "2022-06-27T13:10:30.387927Z"
    }
   },
   "outputs": [],
   "source": [
    "# ECFP配合相似性方法使用\n",
    "if ECFP_SWITCH and SIMILARITY_SWITCH:\n",
    "    data = data.astype(float)\n",
    "    for i in range(data.shape[0]):\n",
    "        emp = []\n",
    "        emp_list = []\n",
    "        fill = []\n",
    "        fill_list = []\n",
    "        for j in range(len(title)):\n",
    "            if data[i, j]>=1:\n",
    "                fill.append(title[j])\n",
    "                fill_list.append(j)\n",
    "            elif data[i, j]==0:\n",
    "                emp.append(title[j])\n",
    "                emp_list.append(j)\n",
    "        for smi1_idx in emp_list:\n",
    "            similarity = []\n",
    "            smi1 = title[smi1_idx]\n",
    "            fps1 = Chem.RDKFingerprint(Chem.MolFromSmarts(smi1))\n",
    "            for smi2_idx in fill_list:\n",
    "                smi2 = title[smi2_idx]\n",
    "                fps2 = Chem.RDKFingerprint(Chem.MolFromSmarts(smi2))\n",
    "                if SIMILARITY_METHOD=='Tanimoto':\n",
    "                    sm = DataStructs.FingerprintSimilarity(fps1,fps2)\n",
    "                elif SIMILARITY_METHOD=='MACCS':\n",
    "                    sm = DataStructs.FingerprintSimilarity(fps1,fps2,metric=DataStructs.DiceSimilarity)\n",
    "                sm *= data[i, smi2_idx]\n",
    "                similarity.append(sm)\n",
    "            if SIMILARITY_MODE=='max':\n",
    "                try:\n",
    "                    sim = max(similarity)\n",
    "                except:\n",
    "                    sim = 0.0\n",
    "            elif SIMILARITY_MODE=='mean':\n",
    "                try:\n",
    "                    sim = np.mean(similarity)\n",
    "                except:\n",
    "                    sim = 0.0\n",
    "            data[i, smi1_idx] += sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:10:30.668316Z",
     "start_time": "2022-06-27T13:10:30.520026Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(856, 564)\n",
      "564\n"
     ]
    }
   ],
   "source": [
    "# ECFP数据集的整理\n",
    "if ECFP_SWITCH:\n",
    "    print(data.shape)\n",
    "    print(len(title))\n",
    "    X = data\n",
    "    y = dataset_sparse.get_shard(0)[1]\n",
    "    title = np.array(title).reshape(X.shape[1], 1)\n",
    "    X_out.append(X)\n",
    "    title_out.append(title)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### =======MFF=======\n",
    "### =======分割=======\n",
    "### =======RDKit======="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:11:44.134087Z",
     "start_time": "2022-06-27T13:10:30.670935Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(856, 110)\n"
     ]
    }
   ],
   "source": [
    "# RDKit描述符特征器模块\n",
    "if RDKIT_DESC_SWITCH:\n",
    "    descriptors = []\n",
    "    descList = []\n",
    "    for descriptor, function in Descriptors.descList:\n",
    "        if descriptor in allowedDescriptors:\n",
    "            descriptors.append(descriptor)\n",
    "            descList.append((descriptor, function))\n",
    "    title_rd = np.array(descriptors).reshape(len(descriptors), 1)\n",
    "    def _featurize(mol):\n",
    "        rval = []\n",
    "        for desc_name, function in descList:\n",
    "            rval.append(function(mol))\n",
    "        return rval\n",
    "    Mol = []\n",
    "    for i in range(len(smiles_rd)):\n",
    "        mol = Chem.MolFromSmiles(smiles_rd[i])\n",
    "        Mol.append(_featurize(mol))\n",
    "    Mol = np.array(Mol)\n",
    "    print(Mol.shape)\n",
    "    X_out.append(Mol)\n",
    "    title_out.append(title_rd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ======RDKit=======\n",
    "### ======Conju======="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:11:44.158312Z",
     "start_time": "2022-06-27T13:11:44.137193Z"
    }
   },
   "outputs": [],
   "source": [
    "# 初始化原子半径、键长矩阵，用于计算VSA：\n",
    "if CONJU_DESC_SWITCH:\n",
    "    dij_m = np.zeros((54, 54))\n",
    "    ri_m = np.zeros((54, 1))\n",
    "    ri_m[6, 0] = 1.950  # C\n",
    "    ri_m[7, 0] = 1.950  # N\n",
    "    ri_m[8, 0] = 1.779  # O\n",
    "    ri_m[9, 0] = 1.496  # F\n",
    "    ri_m[15, 0] = 2.287  # P\n",
    "    ri_m[16, 0] = 2.185  # S\n",
    "    ri_m[17, 0] = 2.044  # Cl\n",
    "    ri_m[35, 0] = 2.166  # Br\n",
    "    ri_m[53, 0] = 2.358  # I\n",
    "    dij_m[6, 35] = 1.970  # C-Br\n",
    "    dij_m[35, 6] = 1.970\n",
    "    dij_m[7, 35] = 1.840  # N-Br\n",
    "    dij_m[35, 7] = 1.840\n",
    "    dij_m[6, 6] = 1.540  # C-C\n",
    "    dij_m[7, 7] = 1.450  # N-N\n",
    "    dij_m[8, 8] = 1.470  # N-N\n",
    "    dij_m[6, 17] = 1.800  # C-Cl\n",
    "    dij_m[17, 6] = 1.800\n",
    "    dij_m[6, 9] = 1.350  # C-F\n",
    "    dij_m[9, 6] = 1.350\n",
    "    dij_m[6, 53] = 2.120  # C-I\n",
    "    dij_m[53, 6] = 2.120\n",
    "    dij_m[6, 7] = 1.470  # C-N\n",
    "    dij_m[7, 6] = 1.470\n",
    "    dij_m[6, 8] = 1.430  # C-N\n",
    "    dij_m[8, 6] = 1.430\n",
    "    dij_m[6, 15] = 1.850  # C-P\n",
    "    dij_m[15, 6] = 1.850\n",
    "    dij_m[6, 16] = 1.810  # C-S\n",
    "    dij_m[16, 6] = 1.810\n",
    "    dij_m[7, 8] = 1.460  # N-O\n",
    "    dij_m[8, 7] = 1.460\n",
    "    dij_m[7, 15] = 1.600  # N-P\n",
    "    dij_m[15, 7] = 1.600\n",
    "    dij_m[7, 16] = 1.760  # N-S\n",
    "    dij_m[16, 7] = 1.760\n",
    "    dij_m[8, 15] = 1.570  # O-P\n",
    "    dij_m[15, 8] = 1.570\n",
    "    dij_m[8, 16] = 1.570  # O-S\n",
    "    dij_m[16, 8] = 1.570"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:11:44.264565Z",
     "start_time": "2022-06-27T13:11:44.160879Z"
    }
   },
   "outputs": [],
   "source": [
    "def find_conju(mol, a_m, d_m):\n",
    "    global patt_list_d, patt_list_m\n",
    "    ring_list = []\n",
    "    f_a = []\n",
    "    patt = Chem.MolFromSmarts('c')\n",
    "    atomids = mol.GetSubstructMatches(patt)\n",
    "    atoms = mol.GetAtoms()\n",
    "    temp_list = []\n",
    "    \n",
    "    def find_ring(atom_id, found_atoms):\n",
    "        nonlocal a_m, ring_list, f_a, atoms, temp_list\n",
    "        flag = False\n",
    "        c_list = np.argwhere(a_m[atom_id] == 1).flatten().tolist()\n",
    "        for atom in c_list:\n",
    "            if atom not in f_a:\n",
    "                a = atoms[atom]\n",
    "                if a.IsInRing() and str(a.GetHybridization())!='SP3':\n",
    "                    found_atoms.append(atom)\n",
    "                    f_a.append(atom)\n",
    "                    find_ring(atom, found_atoms)\n",
    "                    flag = True\n",
    "        if not flag:\n",
    "            temp_list.append(found_atoms)\n",
    "    \n",
    "    for atom in atomids:\n",
    "        a = atom[0]\n",
    "        if a not in f_a:\n",
    "            find_ring(a, [])\n",
    "        if len(temp_list)>0:\n",
    "            max_ring = temp_list[0]\n",
    "            for l in temp_list:\n",
    "                if len(l)>len(max_ring):\n",
    "                    max_ring = l\n",
    "            ring_list.append(max_ring)\n",
    "        temp_list = []\n",
    "    for patt in patt_list_d:\n",
    "        f = Chem.MolFromSmarts(patt)\n",
    "        atomids = mol.GetSubstructMatches(f)\n",
    "        if len(atomids)>0:\n",
    "            for pair in atomids:\n",
    "                n_l = []\n",
    "                flag_f = False\n",
    "                for a in pair:\n",
    "                    if a in f_a:\n",
    "                        flag_f = True\n",
    "                        break\n",
    "                    neighbors = atoms[a].GetNeighbors()\n",
    "                    for na in neighbors:\n",
    "                        n_l.append(na.GetIdx())\n",
    "                if flag_f:\n",
    "                    continue\n",
    "                temp = []\n",
    "                temp_r_id = []\n",
    "                # 到这里，找到了双、叁键的邻接原子\n",
    "                for n in n_l:\n",
    "                    if atoms[n].GetAtomicNum() in [6, 7, 8]:\n",
    "                        for i in range(len(ring_list)):\n",
    "                            ring = ring_list[i]\n",
    "                            if n in ring:\n",
    "                                temp.append(ring)\n",
    "                                temp_r_id.append(i)\n",
    "                if len(temp)==1:\n",
    "                    ring_list[temp_r_id[0]].append(pair[0])\n",
    "                    ring_list[temp_r_id[0]].append(pair[1])\n",
    "                    f_a.append(pair[0])\n",
    "                    f_a.append(pair[1])\n",
    "                else:\n",
    "                    # 合并多个列表\n",
    "                    t_r = []\n",
    "                    for r in temp:\n",
    "                        t_r += r\n",
    "                    # 加上双、叁键两端\n",
    "                    t_r.append(pair[0])\n",
    "                    t_r.append(pair[1])\n",
    "                    # 删掉原来的环\n",
    "                    temp_r_id.sort()\n",
    "                    temp_r_id = np.unique(temp_r_id)\n",
    "                    for i in reversed(temp_r_id):\n",
    "                        del ring_list[i]\n",
    "                    # 加上新环\n",
    "                    ring_list.append(t_r)\n",
    "                    f_a.append(pair[0])\n",
    "                    f_a.append(pair[1])\n",
    "    for patt in patt_list_m:\n",
    "        f = Chem.MolFromSmarts(patt)\n",
    "        atomids = mol.GetSubstructMatches(f)\n",
    "        if len(atomids)>0:\n",
    "            for atom in atomids:\n",
    "                a = atom[0]\n",
    "                if a not in f_a:\n",
    "                    neighbors = atoms[a].GetNeighbors()\n",
    "                    n_l = []\n",
    "                    for na in neighbors:\n",
    "                        n_l.append(na.GetIdx())\n",
    "                    temp = []\n",
    "                    temp_r_id = []\n",
    "                    # 到这里，找到了杂原子的邻接原子\n",
    "                    for n in n_l:\n",
    "                        for i in range(len(ring_list)):\n",
    "                            ring = ring_list[i]\n",
    "                            if (n in ring) and (i not in temp_r_id):\n",
    "                                temp.append(ring)\n",
    "                                temp_r_id.append(i)\n",
    "                    if len(temp)==1:\n",
    "                        ring_list[temp_r_id[0]].append(a)\n",
    "                        f_a.append(a)\n",
    "                    else:\n",
    "                        # 合并多个列表\n",
    "                        t_r = []\n",
    "                        for r in temp:\n",
    "                            t_r += r\n",
    "                        # 加上杂原子\n",
    "                        t_r.append(a)\n",
    "                        # 删掉原来的环\n",
    "                        temp_r_id.sort()\n",
    "                        for i in reversed(temp_r_id):\n",
    "                            del ring_list[i]\n",
    "                        # 加上新环\n",
    "                        if len(t_r)>1:\n",
    "                            ring_list.append(t_r)\n",
    "                            f_a.append(a)\n",
    "    for i in range(len(atoms)):\n",
    "        if i not in f_a:\n",
    "            aa = atoms[i]\n",
    "            if aa.GetSymbol()!='C' or str(aa.GetHybridization())!='SP2':\n",
    "                continue\n",
    "            aa_n = aa.GetNeighbors()\n",
    "            flag = False\n",
    "            for aaa in aa_n:\n",
    "                if aaa.GetIdx() in f_a:\n",
    "                    flag = True\n",
    "                    break\n",
    "            if flag:\n",
    "                a = i\n",
    "                neighbors = atoms[a].GetNeighbors()\n",
    "                n_l = []\n",
    "                for na in neighbors:\n",
    "                    n_l.append(na.GetIdx())\n",
    "                temp = []\n",
    "                temp_r_id = []\n",
    "                # 到这里，找到了杂原子的邻接原子\n",
    "                for n in n_l:\n",
    "                    for i in range(len(ring_list)):\n",
    "                        ring = ring_list[i]\n",
    "                        if (n in ring) and (i not in temp_r_id):\n",
    "                            temp.append(ring)\n",
    "                            temp_r_id.append(i)\n",
    "                if len(temp)==1:\n",
    "                    ring_list[temp_r_id[0]].append(a)\n",
    "                    f_a.append(a)\n",
    "                else:\n",
    "                    # 合并多个列表\n",
    "                    t_r = []\n",
    "                    for r in temp:\n",
    "                        t_r += r\n",
    "                    # 加上杂原子\n",
    "                    t_r.append(a)\n",
    "                    # 删掉原来的环\n",
    "                    temp_r_id.sort()\n",
    "                    for i in reversed(temp_r_id):\n",
    "                        del ring_list[i]\n",
    "                    # 加上新环\n",
    "                    if len(t_r)>1:\n",
    "                        ring_list.append(t_r)\n",
    "                        f_a.append(a)\n",
    "    # 最后核对共轭结构是否相连\n",
    "    if len(ring_list)>1:\n",
    "        temp_count = 0\n",
    "        flag = True\n",
    "        while flag:\n",
    "            t_temp = int(len(ring_list)*(len(ring_list)-1)/2)\n",
    "            temp = 0\n",
    "            break_flag = False\n",
    "            for i in range(len(ring_list)-1):\n",
    "                for j in range(len(ring_list)-i-1):\n",
    "                    r_1 = ring_list[i]\n",
    "                    r_2 = ring_list[i+j+1]\n",
    "                    if np.sum(a_m[r_1, :][:, r_2]) == 0:\n",
    "                        temp += 1\n",
    "                    else:\n",
    "                        # 需要进行合并\n",
    "                        for k in r_2:\n",
    "                            ring_list[i].append(k)\n",
    "                        ring_list[i] = np.unique(ring_list[i])\n",
    "                        del ring_list[i+j+1]\n",
    "                        break_flag = True\n",
    "                        break\n",
    "                if break_flag:\n",
    "                    break\n",
    "            if temp == t_temp:\n",
    "                flag = False\n",
    "    for i in range(len(ring_list)):\n",
    "        ring_list[i] = np.unique(ring_list[i]).flatten().tolist()\n",
    "    return (ring_list, f_a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:11:44.363009Z",
     "start_time": "2022-06-27T13:11:44.267152Z"
    }
   },
   "outputs": [],
   "source": [
    "def find_elec_num(kind, hyb):\n",
    "    global one_list, two_list\n",
    "    if kind in one_list:\n",
    "        return 1\n",
    "    elif kind in two_list:\n",
    "        if hyb=='SP' and kind in ['N', 'P', 'O']:\n",
    "            return 1\n",
    "        elif hyb=='SP2' and kind in ['N', 'O', 'S']:\n",
    "            return 1\n",
    "        else:\n",
    "            return 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:11:44.441851Z",
     "start_time": "2022-06-27T13:11:44.365754Z"
    }
   },
   "outputs": [],
   "source": [
    "def find_acc_mean(l_in, conju_size_list):\n",
    "    temp = 0\n",
    "    for i in range(len(l_in)):\n",
    "        temp += l_in[i] * conju_size_list[i]\n",
    "    temp /= sum(conju_size_list)\n",
    "    return temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:50:11.185431Z",
     "start_time": "2022-06-27T13:45:49.886989Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# 共轭描述符特征器模块\n",
    "if CONJU_DESC_SWITCH:\n",
    "    CONJU_TITLE = []\n",
    "    CONJU_PRE = ['Apperant-Elec-Count', 'PEOE-Charge', 'EState-Indice', 'Atomic-LogP', 'Atomic-MR']\n",
    "    CONJU_NUM = 5 + 7 + 16 * len(CONJU_PRE)\n",
    "    data_conju = np.zeros((len(smiles_rd), CONJU_NUM))\n",
    "    descList = []\n",
    "    allowedDescriptors = ['MolWt']\n",
    "    for descriptor, function in Descriptors.descList:\n",
    "        if descriptor in allowedDescriptors:\n",
    "            descList.append((descriptor, function))\n",
    "    mff_title = title_out[0].flatten().tolist()\n",
    "    for _ in range(len(smiles_rd)):\n",
    "        smi = smiles_rd[_]\n",
    "        mol = Chem.MolFromSmiles(Chem.MolToSmiles(Chem.MolFromSmiles(smi)))\n",
    "        atoms = mol.GetAtoms()\n",
    "        a_m = Chem.rdmolops.GetAdjacencyMatrix(mol)\n",
    "        d_m = Chem.rdmolops.GetDistanceMatrix(mol)\n",
    "        res = find_conju(mol, a_m, d_m)\n",
    "        ring_list = res[0]\n",
    "        lstmax = sorted(ring_list, key = len)[-1]\n",
    "        f_a = res[1]\n",
    "        conju_size_list = [len(r) for r in ring_list]\n",
    "        ## 计算原子级描述符列表：\n",
    "        # 计算表观共轭电荷\n",
    "        app_elec = []\n",
    "        for a in range(len(atoms)):\n",
    "            if a in f_a:\n",
    "                a_kind = atoms[a].GetSymbol()\n",
    "                hyb = str(atoms[a].GetHybridization())\n",
    "                app_elec.append(find_elec_num(a_kind, hyb))\n",
    "            else:\n",
    "                app_elec.append(0)\n",
    "        # 计算PEOE电荷\n",
    "        AllChem.ComputeGasteigerCharges(mol, nIter=25)\n",
    "        peoe_charge = [mol.GetAtomWithIdx(i).GetDoubleProp('_GasteigerCharge') for i in range(mol.GetNumAtoms())]\n",
    "        # 计算电子拓扑系数EState\n",
    "        estate_index = EStateIndices(mol)\n",
    "        # 计算LogP与MR的原子贡献\n",
    "        contribs = rdMolDescriptors._CalcCrippenContribs(mol)\n",
    "        logp = [contribs[i][0] for i in range(len(contribs))]\n",
    "        mr = [contribs[i][1] for i in range(len(contribs))]\n",
    "        # 结合各原子级描述符\n",
    "        atom_props = [app_elec, peoe_charge, estate_index, logp, mr]\n",
    "        \n",
    "        ## 共轭结构的总和特征：\n",
    "        if _ == 0:\n",
    "            CONJU_TITLE.append('Num-Conju-Stru')  # 1.\n",
    "            CONJU_TITLE.append('Conju-Num-Atom-All')  # 2.\n",
    "            CONJU_TITLE.append('Conju-Num-Atom-Ratio')  # 3.\n",
    "            CONJU_TITLE.append('Conju-Wt-Ave')  # 4.\n",
    "            CONJU_TITLE.append('Full-Mol Wiener Index')  # 5.\n",
    "            CONJU_TITLE.append('Conju-Num-Atom-Individual')  # 6.\n",
    "            CONJU_TITLE.append('Conju-Wt-Part')  # 7.\n",
    "            CONJU_TITLE.append('Conju-Wt-Ave')  # 8.\n",
    "            CONJU_TITLE.append('Conju-Max-Distance')  # 9.\n",
    "            CONJU_TITLE.append('Conju-Branch-Ratio')  # 10.\n",
    "            CONJU_TITLE.append('Conju-Stru-Wiener-Index')  # 11.\n",
    "            CONJU_TITLE.append('Conju-Stru-VSA')  # 12.\n",
    "        # 1.共轭结构数量\n",
    "        data_conju[_, 0] = len(ring_list)\n",
    "        # 2.共轭原子总数\n",
    "        data_conju[_, 1] = len(f_a)\n",
    "        # 3.共轭结构数量占比\n",
    "        data_conju[_, 2] = len(f_a)/len(atoms)\n",
    "        # 4.共轭结构质量占比\n",
    "        rval = []\n",
    "        for desc_name, function in descList:\n",
    "            rval.append(function(mol))\n",
    "        wt_list = []\n",
    "        mwt_list = []\n",
    "        for r in ring_list:\n",
    "            tt = 0\n",
    "            for a in r:\n",
    "                tt += atoms[a].GetMass()\n",
    "            wt_list.append(tt)\n",
    "            mwt_list.append(tt/len(r))\n",
    "        data_conju[_, 3] = sum(wt_list)/rval[0]\n",
    "        # 5.全分子维纳指数\n",
    "        if '.' in smi:\n",
    "            mol22 = Chem.MolFromSmiles(smi.split('.')[0])\n",
    "            dm22 = Chem.rdmolops.GetDistanceMatrix(mol22)\n",
    "            data_conju[_, 4] = np.sum(dm22)/(2*dm22.shape[0]*(dm22.shape[0]-1))\n",
    "        else:\n",
    "            data_conju[_, 4] = np.sum(d_m)/(2*d_m.shape[0]*(d_m.shape[0]-1))\n",
    "        \n",
    "        ## 共轭结构独立特征\n",
    "        conju_props = []\n",
    "        # 6.共轭结构独立原子数\n",
    "        size_l = []\n",
    "        for r in ring_list:\n",
    "            size_l.append(len(r))\n",
    "        conju_props.append(size_l)\n",
    "        \n",
    "        # 7.共轭结构独立质量\n",
    "        # 8.共轭结构独立原子平均质量\n",
    "        conju_props.append(wt_list)\n",
    "        conju_props.append(mwt_list)\n",
    "        # 9. 共轭结构长度\n",
    "        conju_max_dis = []\n",
    "        for r in ring_list:\n",
    "            conju_max_dis.append(np.max(d_m[r, :][:, r]))\n",
    "        conju_props.append(conju_max_dis)\n",
    "        \n",
    "        # 10.共轭结构分支系数\n",
    "#         branch_l = []\n",
    "#         for i in range(len(wt_list)):\n",
    "#             branch_l.append(np.sum(a_m[ring_list[i], :][:, ring_list[i]])/(2*size_l[i]))\n",
    "    \n",
    "        for i in range(len(ring_list)):\n",
    "            d_mring=d_m[ring_list[i], :][:, ring_list[i]] \n",
    "        dlist=np.where(d_mring==np.max(d_mring))\n",
    "        l_dlist=len(dlist[0])/2\n",
    "        ltemp=dlist[0][0:round(l_dlist)]\n",
    "        distemp=np.max(ltemp)-np.min(ltemp)\n",
    "        conjuratio=[distemp/max(conju_max_dis)]        \n",
    "        conju_props.append(conjuratio)\n",
    "        # 11.共轭结构维纳指数\n",
    "        wi_l = []\n",
    "        for r in ring_list:\n",
    "            d_m_temp = d_m[r, :][:, r]\n",
    "            wi_l.append(np.sum(d_m_temp)/(2*d_m_temp.shape[0]*(d_m_temp.shape[0]-1)))\n",
    "        conju_props.append(wi_l)\n",
    "        \n",
    "        # 12.共轭结构VSA\n",
    "        conju_vsa_l = []\n",
    "        for r in ring_list:\n",
    "            vsa_t = 0\n",
    "            for i in range(len(r)):\n",
    "                vsa_tt = 0\n",
    "                atom = atoms[r[i]]\n",
    "                n_l = atom.GetNeighbors()\n",
    "                aid_1 = atom.GetAtomicNum()\n",
    "                ar_1 = ri_m[aid_1, 0]\n",
    "                for j in range(len(n_l)):\n",
    "                    aid_2 = n_l[j].GetAtomicNum()\n",
    "                    ar_2 = ri_m[aid_2, 0]\n",
    "                    dij_i = dij_m[aid_1, aid_2]\n",
    "                    dij = min(max(abs(ar_1-ar_2), dij_i), ar_1+ar_2)\n",
    "                    vsa_tt += (ar_2**2-(ar_1-dij)**2)/dij\n",
    "                vsa_t += 4*np.pi*ar_1**2 - np.pi*ar_1*vsa_tt\n",
    "            conju_vsa_l.append(vsa_t)\n",
    "        conju_props.append(conju_vsa_l)\n",
    "        \n",
    "        # 接下来计算原子描述符的性质\n",
    "        for __ in range(len(atom_props)):\n",
    "            PRE = 'Conju-'+CONJU_PRE[__]+'-'\n",
    "            PRE = CONJU_PRE[__]+'-'\n",
    "            END_P = ' (MFF-Conju)'\n",
    "            atom_props_list = atom_props[__]\n",
    "            if _ == 0:\n",
    "                CONJU_TITLE.append(PRE+'Sum')  # 13.1.\n",
    "                CONJU_TITLE.append(PRE+'Ave')  # 13.2.\n",
    "                CONJU_TITLE.append(PRE+'Max')  # 13.3.\n",
    "                CONJU_TITLE.append(PRE+'Min')  # 13.4.\n",
    "                CONJU_TITLE.append(PRE+'Delta')  # 13.5. \n",
    "                CONJU_TITLE.append(PRE+'Weighted')  # 13.6.\n",
    "                CONJU_TITLE.append(PRE+'Weighted-Ave')  # 13.7.\n",
    "                CONJU_TITLE.append(PRE+'PositiveDisCoef')  # 13.8.\n",
    "                CONJU_TITLE.append(PRE+'PositiveDisCoef-PairMean')  # 13.9.\n",
    "                CONJU_TITLE.append(PRE+'NegativeDisCoef')  # 13.10.\n",
    "                CONJU_TITLE.append(PRE+'NegativeDisCoef-PairMean')  # 13.11.\n",
    "                CONJU_TITLE.append(PRE+'GradSum')  # 13.12.\n",
    "                CONJU_TITLE.append(PRE+'GradSum-PairMean')  # 13.13. \n",
    "                CONJU_TITLE.append(PRE+'MaxMinDisRatio')  # 13.14.\n",
    "                # CONJU_TITLE.append(PRE+'AppearantU')  # 13.15.\n",
    "                CONJU_TITLE.append(PRE+'LaplaceSum')  # 13.16.\n",
    "                CONJU_TITLE.append(PRE+'Laplace-PairMean')  # 13.17.\n",
    "            # 13.1. 不含碎片的求和\n",
    "            x_count_l = []\n",
    "            # 13.2. 不含碎片的按原子平均的平均值\n",
    "            x_atom_mean_l = []\n",
    "            # 13.3. 含碎片的最大值\n",
    "            x_max_l = []\n",
    "            # 13.4. 含碎片的最小值\n",
    "            x_min_l = []\n",
    "            # 13.5. 含碎片的差值\n",
    "            x_delta_l = []\n",
    "            # 13.6. 影响力\n",
    "            x_infl_l = []\n",
    "            # 13.7. 按原子平均的影响力\n",
    "            x_atom_infl_l = []\n",
    "            # 13.8. 一阶距离影响力系数\n",
    "            x_pos_dis_coef_l = []\n",
    "            # 13.9. 按二元原子对数目平均的一阶距离影响力系数\n",
    "            x_pair_mean_pos_dis_coef_l = []\n",
    "            # 13.10. 负一阶距离影响力系数\n",
    "            x_neg_dis_coef_l = []\n",
    "            # 13.11. 按二元原子对数目平均的负一阶距离影响力系数\n",
    "            x_pair_mean_neg_dis_coef_l = []\n",
    "            # 13.12. 总一阶梯度\n",
    "            x_sum_grad_l = []\n",
    "            # 13.13. 按二元原子对数目平均的一阶梯度\n",
    "            x_pair_mean_grad_l = []\n",
    "            # 13.14. 含碎片的极值-极值距离占比\n",
    "            x_dis_ratio_l = []\n",
    "            # 13.15. 按电荷性质与图论距离计算的势能（与负一阶一样，删除掉）\n",
    "            # x_u_l = []\n",
    "            # 13.16. 总二阶梯度\n",
    "            x_sum_laplace_l = []\n",
    "            # 13.17. 按二元原子对数目平均的二阶梯度\n",
    "            x_pair_mean_laplace_l = []\n",
    "            # 开始计算性质\n",
    "            for r in [lstmax]:\n",
    "                # 整合1.2.\n",
    "                a_p_l = [atom_props_list[a] for a in r]\n",
    "                x_count_l.append(sum(a_p_l))  # 1.\n",
    "                x_atom_mean_l.append(sum(a_p_l)/len(r))  # 2.\n",
    "                # 生成碎片性质列表\n",
    "                frag_x = []\n",
    "                frag_atom_id = []\n",
    "                for i in range(len(mff_title)):\n",
    "                    patt = mff_title[i]\n",
    "                    f = Chem.MolFromSmarts(patt)\n",
    "                    atomids = mol.GetSubstructMatches(f)\n",
    "                    if len(atomids) > 0:\n",
    "                        for j in range(len(atomids)):\n",
    "                            peoe_flag = True\n",
    "                            for k in atomids[j]:\n",
    "                                if k not in r:\n",
    "                                    peoe_flag = False\n",
    "                                    break\n",
    "                            if peoe_flag:\n",
    "                                frag_atom_id.append(atomids[j])\n",
    "                                x_temp = 0\n",
    "                                for k in atomids[j]:\n",
    "                                    x_temp += atom_props_list[k]\n",
    "                                frag_x.append(x_temp)\n",
    "#                     for j in range(len(r)):\n",
    "#                         atom_id = r[j]\n",
    "#                         frag_x.append(atom_props_list[atom_id])\n",
    "#                         frag_atom_id.append([atom_id])\n",
    "                if len(frag_x)==0:\n",
    "                    frag_x=[0]\n",
    "                    frag_atom_id=[()]\n",
    "                x_max_l.append(max(frag_x))  # 3.\n",
    "                x_min_l.append(min(frag_x))  # 4.\n",
    "                x_delta_l.append(max(frag_x) - min(frag_x))  # 5.\n",
    "                f_1 = [i for i in frag_atom_id[frag_x.index(min(frag_x))]]\n",
    "                f_2 = [i for i in frag_atom_id[frag_x.index(max(frag_x))]]\n",
    "                s = np.max(d_m[f_1, :][:, f_2])/np.max(d_m[r, :][:, r])\n",
    "                x_dis_ratio_l.append(s)  # 14.\n",
    "                # 计算影响力\n",
    "                temp_infl = 0\n",
    "                for a in r:\n",
    "                    temp_infl += atom_props_list[a] * atoms[a].GetDegree()\n",
    "                x_infl_l.append(temp_infl)  # 6.\n",
    "                x_atom_infl_l.append(temp_infl/len(r))  # 7.\n",
    "                # 计算距离系数与一阶、二阶梯度\n",
    "                ttp = 0\n",
    "                tti = 0\n",
    "                grad_temp = 0\n",
    "                laplace_temp = 0\n",
    "                count = 0\n",
    "                for i in range(len(r)):\n",
    "                    for j in range(len(r)-i-1):\n",
    "                        a = r[i]\n",
    "                        b = r[i+j+1]\n",
    "                        count += 1\n",
    "                        a_e = atom_props_list[a]\n",
    "                        b_e = atom_props_list[b]\n",
    "                        ttp += a_e*b_e*d_m[a, b]\n",
    "                        tti += a_e*b_e/d_m[a, b]\n",
    "                        grad_temp += abs(a_e-b_e)/d_m[a, b]\n",
    "                        laplace_temp += abs(a_e-b_e)/(d_m[a, b]**2)\n",
    "                x_pos_dis_coef_l.append(ttp)  # 8.\n",
    "                x_pair_mean_pos_dis_coef_l.append(ttp/count)  # 9.\n",
    "                x_neg_dis_coef_l.append(tti)  # 10.\n",
    "                x_pair_mean_neg_dis_coef_l.append(tti/count)  # 11.\n",
    "                x_sum_grad_l.append(grad_temp)  # 12.\n",
    "                x_pair_mean_grad_l.append(grad_temp/count)  # 13.\n",
    "                x_sum_laplace_l.append(laplace_temp)  # 16.\n",
    "                x_pair_mean_laplace_l.append(laplace_temp/count)  # 17.\n",
    "            conju_props.append(x_count_l)\n",
    "            conju_props.append(x_atom_mean_l)\n",
    "            conju_props.append(x_max_l)\n",
    "            conju_props.append(x_min_l)\n",
    "            conju_props.append(x_delta_l)\n",
    "            conju_props.append(x_infl_l)\n",
    "            conju_props.append(x_atom_infl_l)\n",
    "            conju_props.append(x_pos_dis_coef_l)\n",
    "            conju_props.append(x_pair_mean_pos_dis_coef_l)\n",
    "            conju_props.append(x_neg_dis_coef_l)\n",
    "            conju_props.append(x_pair_mean_neg_dis_coef_l)\n",
    "            conju_props.append(x_sum_grad_l)\n",
    "            conju_props.append(x_pair_mean_grad_l)\n",
    "            conju_props.append(x_dis_ratio_l)\n",
    "            # conju_props.append(x_u_l)\n",
    "            conju_props.append(x_sum_laplace_l)\n",
    "            conju_props.append(x_pair_mean_laplace_l)\n",
    "        # 开始填入特征\n",
    "        for i in range(len(conju_props)):\n",
    "            index = 5 + i\n",
    "            if KEEP_TYPE == 'max':\n",
    "                data_conju[_, index] = max(conju_props[i])\n",
    "            elif KEEP_TYPE == 'mean':\n",
    "                data_conju[_, index] = mean(conju_props[i])\n",
    "            elif KEEP_TYPE == 'acc-mean':\n",
    "                temp = 0\n",
    "                for j in range(len(conju_size_list)):\n",
    "                    temp += conju_size_list[j] * conju_props[i][j]\n",
    "                data_conju[_, index] = temp / sum(conju_size_list)\n",
    "    X_out.append(data_conju)\n",
    "    title_out.append(np.array(CONJU_TITLE).reshape(len(CONJU_TITLE), 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:50:18.667322Z",
     "start_time": "2022-06-27T13:50:18.645831Z"
    }
   },
   "outputs": [],
   "source": [
    "if SOLVENT_SWITCH:\n",
    "    sol_d = np.loadtxt(SOLVENT_IN, dtype=float, delimiter=',')\n",
    "    X_out.append(sol_d)\n",
    "    title_out.append(np.array(SOLVENT_TITLE).reshape(3, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:50:18.922195Z",
     "start_time": "2022-06-27T13:50:18.853510Z"
    }
   },
   "outputs": [],
   "source": [
    "if WAVE_SWITCH:\n",
    "    wave_d = np.loadtxt(WAVE_IN, dtype=float, delimiter=',')\n",
    "    data_w = np.zeros((len(smiles_rd), 1))\n",
    "    for i in range(len(smiles_rd)):\n",
    "        data_w[i, 0] = wave_d[i, ]\n",
    "    X_out.append(data_w)\n",
    "    title_out.append(np.array(WAVE_TITLE).reshape(1, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:50:19.029472Z",
     "start_time": "2022-06-27T13:52:09.991Z"
    }
   },
   "outputs": [],
   "source": [
    "# 预留给更多的特征化模块的空间"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:50:31.233255Z",
     "start_time": "2022-06-27T13:50:19.145848Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(856, 770) (1, 770)\n",
      "(857, 775)\n"
     ]
    }
   ],
   "source": [
    "# 将所有对应的数据包叠加起来，得到一个大数据包\n",
    "X_init = X_out[0]\n",
    "if len(X_out)>=2:\n",
    "    for i in range(len(X_out)-1):\n",
    "        X_init = np.hstack((X_init, X_out[i+1]))\n",
    "title_init = title_out[0]\n",
    "if len(title_out)>=2:\n",
    "    for i in range(len(title_out)-1):\n",
    "        title_init = np.vstack((title_init, title_out[i+1]))\n",
    "title_init = np.transpose(title_init)\n",
    "print(X_init.shape, title_init.shape)\n",
    "title_extra = ['smiles', 'values']\n",
    "data_out = np.hstack((X_init, smiles_out))\n",
    "data_out = np.hstack((data_out, values_out))\n",
    "if VALUES_DIV_MOLWT:\n",
    "    title_extra.append('values_dwt')\n",
    "    data_out = np.hstack((data_out, values_dwt_out))\n",
    "if VALUES_LN:\n",
    "    title_extra.append('values_ln')\n",
    "    data_out = np.hstack((data_out, np.log(values_out)))\n",
    "    if VALUES_DIV_MOLWT:\n",
    "        title_extra.append('values_dwt_ln')\n",
    "        data_out = np.hstack((data_out, np.log(values_dwt_out)))\n",
    "title_extra = np.array(title_extra).reshape(1, len(title_extra))\n",
    "title_out = np.hstack((title_init, title_extra))\n",
    "full_data = np.vstack((title_out, data_out))\n",
    "print(full_data.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:52:45.940069Z",
     "start_time": "2022-06-27T13:52:27.590791Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[]\n",
      "['NumRadicalElectrons' 'SMR_VSA8' 'SlogP_VSA9']\n",
      "(857, 772)\n"
     ]
    }
   ],
   "source": [
    "# 删除含有nan或infi的行，以及只有0的列\n",
    "del_list = []\n",
    "for i in range(full_data.shape[0]-1):\n",
    "    j = i+1\n",
    "    if (np.isnan(full_data[j, :-5].astype(float)).sum()>0)or(np.inf in full_data[j, :-5].astype(float)):\n",
    "        del_list.append(j)\n",
    "        continue\n",
    "print(full_data[0, del_list])\n",
    "full_data = np.delete(full_data, del_list, axis=0)\n",
    "if X_init.shape[0]>1:\n",
    "    del_list = []\n",
    "    for i in range(X_init.shape[1]):\n",
    "        if max(full_data[1:, i].astype(float))==min(full_data[1:, i].astype(float)):\n",
    "            del_list.append(i)\n",
    "            continue\n",
    "    print(full_data[0, del_list])\n",
    "    full_data = np.delete(full_data, del_list, axis=1)\n",
    "print(full_data.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:52:59.484161Z",
     "start_time": "2022-06-27T13:52:52.985509Z"
    }
   },
   "outputs": [],
   "source": [
    "# 将大数据包拆分开来\n",
    "OUT_NAME_FULL = 'Full_'+str(full_data.shape[0])+'_'+str(full_data.shape[1])+'.csv'\n",
    "OUT_NAME_FULL = Path('.', DIR, OUT_NAME_FULL)\n",
    "np.savetxt(OUT_NAME_FULL, full_data, fmt='%s', delimiter=',')\n",
    "data_t = full_data[1:, :]\n",
    "full_t = full_data\n",
    "if VALUES_LN:\n",
    "    if VALUES_DIV_MOLWT:\n",
    "        values_out_wt_ln = data_t[:, -1]\n",
    "        OUT_NAME_VALUES_WT_LN = 'Values_True_w_ln_'+str(values_out_wt_ln.shape[0])+'.csv'\n",
    "        OUT_NAME_VALUES_WT_LN = Path('.', DIR, OUT_NAME_VALUES_WT_LN)\n",
    "        np.savetxt(OUT_NAME_VALUES_WT_LN, values_out_wt_ln, fmt='%s', delimiter=',')\n",
    "        data_t = data_t[:, :-1]\n",
    "        full_t = full_t[:, :-1]\n",
    "    values_out_ln = data_t[:, -1]\n",
    "    OUT_NAME_VALUES_LN = 'Values_True_ln_'+str(values_out_ln.shape[0])+'.csv'\n",
    "    OUT_NAME_VALUES_LN = Path('.', DIR, OUT_NAME_VALUES_LN)\n",
    "    np.savetxt(OUT_NAME_VALUES_LN, values_out_ln, fmt='%s', delimiter=',')\n",
    "    data_t = data_t[:, :-1]\n",
    "    full_t = full_t[:, :-1]\n",
    "if VALUES_DIV_MOLWT:\n",
    "    values_out_wt = data_t[:, -1]\n",
    "    OUT_NAME_VALUES_WT = 'Values_True_w_'+str(values_out_wt.shape[0])+'.csv'\n",
    "    OUT_NAME_VALUES_WT = Path('.', DIR, OUT_NAME_VALUES_WT)\n",
    "    np.savetxt(OUT_NAME_VALUES_WT, values_out_wt, fmt='%s', delimiter=',')\n",
    "    data_t = data_t[:, :-1]\n",
    "    full_t = full_t[:, :-1]\n",
    "values_out = data_t[:, -1]\n",
    "OUT_NAME_VALUES = 'Values_True_'+str(values_out.shape[0])+'.csv'\n",
    "OUT_NAME_VALUES = Path('.', DIR, OUT_NAME_VALUES)\n",
    "np.savetxt(OUT_NAME_VALUES, values_out, fmt='%s', delimiter=',')\n",
    "data_t = data_t[:, :-1]\n",
    "full_t = full_t[:, :-1]\n",
    "smiles_out = data_t[:, -1]\n",
    "OUT_NAME_SMILES = 'Smiles_'+str(values_out.shape[0])+'.csv'\n",
    "OUT_NAME_SMILES = Path('.', DIR, OUT_NAME_SMILES)\n",
    "np.savetxt(OUT_NAME_SMILES, smiles_out, fmt='%s', delimiter=',')\n",
    "X_out = data_t[:, :-1]\n",
    "full_t = full_t[:, :-1]\n",
    "OUT_NAME_X = 'Features_'+str(X_out.shape[0])+'_'+str(X_out.shape[1])+'.csv'\n",
    "OUT_NAME_X = Path('.', DIR, OUT_NAME_X)\n",
    "np.savetxt(OUT_NAME_X, X_out, fmt='%s', delimiter=',')\n",
    "OUT_NAME_TITLE = 'Title_'+str(full_t.shape[1])+'.csv'\n",
    "OUT_NAME_TITLE = Path('.', DIR, OUT_NAME_TITLE)\n",
    "np.savetxt(OUT_NAME_TITLE, np.transpose(full_t[0, :]), fmt='%s', delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:52:59.560067Z",
     "start_time": "2022-06-27T13:54:46.419Z"
    }
   },
   "outputs": [],
   "source": [
    "# 输出综合的数据集，供pandas使用\n",
    "if PANDAS_DATASET_GENERATE:\n",
    "    dataset_out = np.hstack((full_t, np.hstack((np.array(['values']), values_out)).reshape(full_t.shape[0], 1)))\n",
    "    OUT_NAME_D1 = 'Dataset_'+str(X_out.shape[0])+'_'+str(X_out.shape[1])+'.csv'\n",
    "    OUT_NAME_D1 = Path('.', DIR, OUT_NAME_D1)\n",
    "    np.savetxt(OUT_NAME_D1, dataset_out, fmt='%s', delimiter=',')\n",
    "    if VALUES_DIV_MOLWT:\n",
    "        dataset_out = np.hstack((full_t, np.hstack((np.array(['values_wt']), values_out_wt)).reshape(full_t.shape[0], 1)))\n",
    "        OUT_NAME_D2 = 'Dataset_w_'+str(X_out.shape[0])+'_'+str(X_out.shape[1])+'.csv'\n",
    "        OUT_NAME_D2 = Path('.', DIR, OUT_NAME_D2)\n",
    "        np.savetxt(OUT_NAME_D2, dataset_out, fmt='%s', delimiter=',')\n",
    "    if VALUES_LN:\n",
    "        dataset_out = np.hstack((full_t, np.hstack((np.array(['values_ln']), values_out_ln)).reshape(full_t.shape[0], 1)))\n",
    "        OUT_NAME_D3 = 'Dataset_ln_'+str(X_out.shape[0])+'_'+str(X_out.shape[1])+'.csv'\n",
    "        OUT_NAME_D3 = Path('.', DIR, OUT_NAME_D3)\n",
    "        np.savetxt(OUT_NAME_D3, dataset_out, fmt='%s', delimiter=',')\n",
    "        if VALUES_DIV_MOLWT:\n",
    "            dataset_out = np.hstack((full_t, np.hstack((np.array(['values_wt_ln']), values_out_wt_ln)).reshape(full_t.shape[0], 1)))\n",
    "            OUT_NAME_D4 = 'Dataset_w_ln_'+str(X_out.shape[0])+'_'+str(X_out.shape[1])+'.csv'\n",
    "            OUT_NAME_D4 = Path('.', DIR, OUT_NAME_D4)\n",
    "            np.savetxt(OUT_NAME_D4, dataset_out, fmt='%s', delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-27T13:52:59.577650Z",
     "start_time": "2022-06-27T13:54:47.129Z"
    }
   },
   "outputs": [],
   "source": [
    "# 输出Log文件\n",
    "LOG_NAME = Path('.', DIR, LOG_NAME)\n",
    "f1 = open(LOG_NAME, 'w+')\n",
    "f1.write('Log for CMF\\n')\n",
    "f1.write('Version: V'+VERSION+'\\n\\n')\n",
    "f1.write('Log generation time: '+time.strftime(\"%Y.%m.%d-%H:%M:%S\", time.localtime())+'\\n\\n')\n",
    "if SMILES_CHECK:\n",
    "    f1.write('Smiles check is on.\\n')\n",
    "if VALUES_DIV_MOLWT:\n",
    "    f1.write('Data with values divided by molecular weight is generated.\\n')\n",
    "if VALUES_LN:\n",
    "    f1.write('Data with ln(values) is generated.\\n')\n",
    "if PANDAS_DATASET_GENERATE:\n",
    "    f1.write('Dataset(s) which can be used by pandas is(are) generated.\\n')\n",
    "f1.write('\\n\\n')\n",
    "f1.write('Parameters:\\n\\n')\n",
    "if ECFP_SWITCH:\n",
    "    f1.write('ECFP(modified) featurizer is on.\\n')\n",
    "    f1.write('ECFP radius: '+str(RADIUS)+'\\n')\n",
    "    f1.write('Rate of feature reservation: '+str(FEATURE_RESERVE)+'\\n')\n",
    "    if EXACT_PIECE_RESERVE:\n",
    "        f1.write('   Exact smiles piece reservation is on.\\n')\n",
    "        f1.write('   Smiles piece(s): '+str(SMILES_PATT)+'\\n')\n",
    "        f1.write('   Rate of feature reservation for exact smiles piece(s): '+str(EXACT_PIECE_RESERVE_THRESHOLD)+'\\n')\n",
    "    if TURN_TO_BOOL:\n",
    "        f1.write('Turn-to-bool module is on.\\n')\n",
    "        f1.write('Turn-to-bool mode: '+BOOL_TURN_MODE+'\\n')\n",
    "        if ALL_TURN_TO_BOOL:\n",
    "            f1.write('All features will be turned to bool values.\\n')\n",
    "        else:\n",
    "            if len(BOOL_PATT_LIST)>0:\n",
    "                f1.write('Smiles pieces to be turned into bool: '+str(BOOL_PATT_LIST)+'\\n')\n",
    "    if ATOM_COUNT_CONTROL:\n",
    "        f1.write('Atom-count-control module is on.\\n')\n",
    "        f1.write('Minimum atom (except H) number of allowed smarts piece: '+str(SMARTS_MIN_LENGTH)+'\\n')\n",
    "        f1.write('Maximum atom (except H) number of allowed smarts piece: '+str(SMARTS_MAX_LENGTH)+'\\n')\n",
    "        if len(ATOM_COUNT_CONTROL_OMIT_PATTERN)>0:\n",
    "            f1.write('Excepption smarts pieces of Atom-count-control module: '+str(ATOM_COUNT_CONTROL_OMIT_PATTERN)+'\\n')\n",
    "    if SIMILARITY_SWITCH:\n",
    "        f1.write('Fingerprint-Similarity module is on.\\n')\n",
    "        f1.write('Fingerprint-Similarity method: '+SIMILARITY_METHOD+'\\n')\n",
    "        f1.write('Fingerprint-Similarity mode: '+SIMILARITY_MODE+'\\n')\n",
    "    f1.write('\\n')\n",
    "if RDKIT_DESC_SWITCH:\n",
    "    f1.write('RDKit Descriptor featurizer is on.\\n')\n",
    "    f1.write('Length of allowed descriptor set: '+str(len(allowedDescriptors))+'\\n')\n",
    "    f1.write('Allowed descriptors:\\n')\n",
    "    for i in range(len(allowedDescriptors)):\n",
    "        if (i)%5 == 0:\n",
    "            f1.write('    ')\n",
    "        f1.write(allowedDescriptors[i]+'   ')\n",
    "        if (i+1)%5 == 0:\n",
    "            f1.write('\\n')\n",
    "    f1.write('\\n')\n",
    "f1.write('Some samples may be deleted beacause of the presence of nan or infinite.\\n')\n",
    "f1.write('Some columns may be deleted beacause of only 0 is contained.\\n')\n",
    "f1.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "PythonForDYH",
   "language": "python",
   "name": "dyhpy"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "position": {
    "height": "363.797px",
    "left": "2500.03px",
    "right": "20px",
    "top": "120px",
    "width": "359.977px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}