{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom sklearn import metrics\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import train_test_split\n\nimport xgboost as xgb","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def metric(y_true, y_pred, types, floor=1e-9):\n    maes = (y_true-y_pred).abs().groupby(types).mean()\n    return np.log(maes.map(lambda x: max(x, floor))).mean()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 读取数据\ndf = dict()\n\ndf['train'] = pd.read_csv('../input/train.csv', index_col='id')\ndf['infer'] = pd.read_csv('../input/test.csv', index_col='id')\ndf['structure'] = pd.read_csv('../input/structures.csv')\ndf['charges'] = pd.read_csv('../input/mulliken_charges.csv')\n\n# merge both datasets\ndf['full'] = pd.concat([df['train'], df['infer']], axis=0, sort=False)\n\n# free up some memory\ndel df['train'], df['infer']","metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"key=\"full\"\nfor atom in [0,1]:\n\n    # adding atom coordinates\n    df[key] = pd.merge(df[key], df['structure'], how = 'left',\n                        left_on  = ['molecule_name', f'atom_index_{atom}'],\n                        right_on = ['molecule_name',  'atom_index'], \n                    )\n\n    df[key] = df[key].drop('atom_index', axis=1)\n    df[key] = df[key].rename(columns={\n                            'atom': f'A{atom}',\n                            'x': f'x{atom}',\n                            'y': f'y{atom}',\n                            'z': f'z{atom}'})\n\nm_data = df[key] \n# 分子中的原子数\natom_num = pd.DataFrame(df['structure'].groupby('molecule_name')['atom_index'].max())\nm_data = pd.merge(m_data, atom_num, how='left',left_on='molecule_name', right_on='molecule_name').rename(columns={'atom_index' : 'nAtoms'})\n# 原子键的类型\nm_data.loc[:, 'bond'] = [value[:2] for value in m_data['type'].values]\n# 原子对\nm_data.loc[:, 'pair'] = [value[2:] for value in m_data['type'].values]\n# 构成一个原子对两个原子的L1 L2距离\nm_data['L2'] = np.linalg.norm(m_data[['x0', 'y0', 'z0']].values - m_data[['x1', 'y1', 'z1']].values, ord=2, axis=1)\nm_data['L1'] = np.linalg.norm(m_data[['x0', 'y0', 'z0']].values - m_data[['x1', 'y1', 'z1']].values, ord=1, axis=1)\n# 构成一个原子对的两个原子的距离在三个维度上的分量\nm_data.loc[:, 'dx'] = np.abs(m_data.loc[:, 'x0'] - m_data.loc[:, 'x1'])\nm_data.loc[:, 'dy'] = np.abs(m_data.loc[:, 'y0'] - m_data.loc[:, 'y1'])\nm_data.loc[:, 'dz'] = np.abs(m_data.loc[:, 'z0'] - m_data.loc[:, 'z1'])\n# 各原子对在一个分子中的相对L1,L2距离(该原子对的L1/L2距离 比上 该分子中最大的L1/L2距离)\nm_data.loc[:, 'L1_RM'] = m_data.loc[:, 'L1'] / m_data.groupby('molecule_name')['L1'].transform('max')\nm_data.loc[:, 'L2_RM'] = m_data.loc[:, 'L2'] / m_data.groupby('molecule_name')['L2'].transform('max')\n# 各原子对在所有由atom1组成的原子对中的相对L1 L2距离\nm_data.loc[:, 'L1_RA'] = m_data.loc[:, 'L1'] / m_data.groupby('A1')['L1'].transform('max')\nm_data.loc[:, 'L2_RA'] = m_data.loc[:, 'L2'] / m_data.groupby('A1')['L2'].transform('max')\n# 各原子对在所有同bond原子对中的相对L1、L2距离\nm_data.loc[:, 'L1_RB'] = m_data.loc[:, 'L1'] / m_data.groupby('bond')['L1'].transform('max')\nm_data.loc[:, 'L2_RB'] = m_data.loc[:, 'L2'] / m_data.groupby('bond')['L2'].transform('max')\n# 各原子对在所有同type原子对中的相对L1，l2距离\nm_data.loc[:, 'L1_RT'] = m_data.loc[:, 'L1'] / m_data.groupby('type')['L1'].transform('max')\nm_data.loc[:, 'L2_RT'] = m_data.loc[:, 'L2'] / m_data.groupby('type')['L2'].transform('max')\n\n# 分子的体积 \nvolumes = np.prod(df['structure'].groupby('molecule_name')[['x','y','z']].max() - df['structure'].groupby('molecule_name')[['x','y','z']].min(), axis=1)\nvolumes = pd.DataFrame(volumes) \n\nm_data = pd.merge(m_data, volumes, how='left', left_on='molecule_name', right_index=True)\nm_data = m_data.rename(columns={0 : 'volume'})\n\ndel volumes\n\n# 分子中原子长宽高的最远距离\nsize = df['structure'].groupby('molecule_name')[['x','y','z']].max() - \\\n       df['structure'].groupby('molecule_name')[['x','y','z']].min()\nsize = pd.DataFrame(size) \nm_data = pd.merge(m_data, size, how='left', left_on='molecule_name', right_index=True)\nm_data = m_data.rename(columns={'x':'Lx', 'y':'Ly', 'z':'Lz'})\n\ndel size\n\n# 增加相对原子质量\natomic_weights = {\n    'H' : 1.008,\n    'C' : 12.011,\n    'N' : 14.007,\n    'O' : 15.999,\n    'F' : 18.998\n}\n\ndef weight(element):\n    return atomic_weights[element]\n\ndf['structure'].loc[:, 'atomic_weight'] = list(map(weight, df['structure']['atom']))\n\nweights = df['structure'].groupby('molecule_name')['atomic_weight'].sum()\nweights = pd.DataFrame(weights) \n\nm_data = pd.merge(m_data, weights, how='left', left_on='molecule_name', right_index=True)\nm_data = m_data.rename(columns={'atomic_weight' : 'weight'})\n\ndel weights\n\n# 特征交叉\nm_data.loc[:, 'weightAtom'] = list(map(weight, m_data['A1']))\nm_data.loc[:, 'wA/w'] = m_data.loc[:, 'weightAtom'] / m_data.loc[:, 'weight'] \nm_data.loc[:, 'L1/Lx'] = m_data.loc[:, 'L1'] / m_data.loc[:, 'Lx']\nm_data.loc[:, 'L2/Lx'] = m_data.loc[:, 'L2'] / m_data.loc[:, 'Lx']\nm_data.loc[:, 'L1/Ly'] = m_data.loc[:, 'L1'] / m_data.loc[:, 'Ly']\nm_data.loc[:, 'L2/Ly'] = m_data.loc[:, 'L2'] / m_data.loc[:, 'Ly']\nm_data.loc[:, 'L1/Lz'] = m_data.loc[:, 'L1'] / m_data.loc[:, 'Lz']\nm_data.loc[:, 'L2/Lz'] = m_data.loc[:, 'L2'] / m_data.loc[:, 'Lz']\nm_data.loc[:, 'dx/Lx'] = m_data.loc[:, 'dx'] / m_data.loc[:, 'Lx']\nm_data.loc[:, 'dy/Ly'] = m_data.loc[:, 'dy'] / m_data.loc[:, 'Ly']\nm_data.loc[:, 'dz/Lz'] = m_data.loc[:, 'dz'] / m_data.loc[:, 'Lz']\nm_data.loc[:, 'dx/Lx_RB'] = m_data.loc[:, 'dx'] / m_data.groupby('bond')['Lx'].transform('max')\nm_data.loc[:, 'dy/Ly_RB'] = m_data.loc[:, 'dy'] / m_data.groupby('bond')['Ly'].transform('max')\nm_data.loc[:, 'dz/Lz_RB'] = m_data.loc[:, 'dz'] / m_data.groupby('bond')['Lz'].transform('max')\n\nm_data.head(10)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 将字符串特征 type bond pair 热编码\nencoder = LabelEncoder()\nfor feature in ['type', 'bond', 'pair']:\n    m_data[feature] = encoder.fit_transform(m_data[feature])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 训练集测试集\nX, y = dict(), dict()\n\nto_drop = [\n    'A0', 'A1', 'x0', 'y0', 'z0', 'x1', 'y1', 'z1', 'weightAtom', 'scalar_coupling_constant'\n]\n\nX['train'] = m_data[m_data['scalar_coupling_constant'].isna() == False].drop(to_drop, axis=1)\nX['infer'] = m_data[m_data['scalar_coupling_constant'].isna() == True].drop(to_drop, axis=1)\n\ny['train'] = m_data[m_data['scalar_coupling_constant'].isna() == False]['scalar_coupling_constant']\n\nX['train'], X['valid'], y['train'], y['valid'] = train_test_split(X['train'], y['train'], \n                                                                  random_state=0, test_size=0.2)\n\ndel m_data","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for key in ['train', 'valid']:\n    # atom0的电荷M0\n    X[key] = pd.merge(X[key], df['charges'], how='left',\n             left_on=['molecule_name', 'atom_index_0'], right_on=['molecule_name', 'atom_index'])\n    \n    X[key] = X[key].rename(columns={'mulliken_charge':'M0'})\n    \n    # atom1的电荷M1\n    X[key] = pd.merge(X[key], df['charges'], how='left',\n             left_on=['molecule_name', 'atom_index_1'], right_on=['molecule_name', 'atom_index'])\n\n    X[key] = X[key].rename(columns={'mulliken_charge':'M1'})\n\n    X[key].drop(['atom_index_x', 'atom_index_y'], axis=1, inplace=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('{}, {}, {}'.format(X['train'].shape, X['valid'].shape, X['infer'].shape))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for key in ['train', 'valid', 'infer']:\n    X[key].drop(['molecule_name', 'atom_index_0', 'atom_index_1'], axis=1, inplace=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# building model for M0\n# M0训练\nmodel0 = xgb.XGBRegressor(\n    n_estimators=1000, \n    learning_rate=0.2,\n    max_depth = 9,\n    n_jobs=-1,\n    random_state=0,\n    subsample=0.8,\n    tree_method='gpu_hist'\n)\n\nto_drop = ['M0', 'M1']\n\nmodel0.fit(\n    X['train'].drop(to_drop,axis=1), X['train']['M0'], \n    eval_set=[(X['valid'].drop(to_drop,axis=1), X['valid']['M0'])], \n    eval_metric='mae',early_stopping_rounds=20)","metadata":{"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# building model for M1\n# 训练M1预测器\nmodel1 = xgb.XGBRegressor(\n    n_estimators=1000, \n    learning_rate=0.2,\n    max_depth = 9,\n    n_jobs=-1,\n    random_state=0,\n    subsample=0.8,\n    tree_method='gpu_hist'\n)\n\nto_drop = ['M0', 'M1']\n\nmodel1.fit(\n    X['train'].drop(to_drop,axis=1), X['train']['M1'], \n    eval_set=[(X['valid'].drop(to_drop,axis=1), X['valid']['M1'])], \n    eval_metric='mae',early_stopping_rounds=20)","metadata":{"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('{}, {}, {}'.format(X['train'].shape, X['valid'].shape, X['infer'].shape))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 预测模型\nmodel = xgb.XGBRegressor(\n    n_estimators=10000, \n    learning_rate=0.2,\n    max_depth = 9,\n    n_jobs=-1,\n    random_state=0,\n    subsample=0.8,\n    tree_method='gpu_hist'\n)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 训练\nmodel.fit(\n    X['train'], y['train'], \n    eval_set=[(X['valid'], y['valid'])], \n    eval_metric='mae',early_stopping_rounds=20)","metadata":{"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"f, ax = plt.subplots(1,1, figsize=(17,8))\nxgb.plot_importance(model, ax)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions = dict()\n\npredictions['train'] = model.predict(X['train'])\npredictions['valid'] = model.predict(X['valid'])\n# predictions['infer'] = model.predict(X['infer'])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"metric(y['train'], predictions['train'], X['train']['type'])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"metric(y['valid'], predictions['valid'], X['valid']['type'])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = {\n    'id' : pd.read_csv('../input/test.csv', index_col='id').index,\n    'scalar_coupling_constant' : predictions['infer']\n}","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = pd.DataFrame(submission)\nsubmission.to_csv('submision_12.csv', index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}