{"cells":[{"metadata":{},"cell_type":"markdown","source":"There are two data sets with molecular properties provided for the train molecules but not the test molecules: dipole moments and potential energy. It is possible that these could be useful for predicting the scalar coupling constant. So, here I use some molecular properties (derived here and in a previous kernel) to fit the magnitude of the dipole moment and the potential energy and store the predicted values on the train and test sets for future work."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"from numpy.random import permutation\nfrom sklearn import metrics\nimport lightgbm\nfrom matplotlib import pyplot as plt\nimport math","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/champs-scalar-coupling/train.csv')\ntest = pd.read_csv('../input/champs-scalar-coupling/test.csv')\nsub = pd.read_csv('../input/champs-scalar-coupling/sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# previously constructed features\ntrain_dist = pd.read_csv('../input/distance-features/train_dist.csv')\ntest_dist = pd.read_csv('../input/distance-features/test_dist.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get xyz data for each atom\nstructures = pd.read_csv('../input/champs-scalar-coupling/structures.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# dipole moments of molecules\ndipole_moments = pd.read_csv('../input/champs-scalar-coupling/dipole_moments.csv')\ntrain = pd.merge(train, dipole_moments, how='left',\n              left_on='molecule_name',\n              right_on='molecule_name')\n# calculate the magnitude (attempts at fitting the components using the below properties failed utterly) note that this is in QM9 (mu)\ntrain['dipole_moment'] = train.apply(lambda row: math.sqrt(row['X']**2 + row['Y']**2 + row['Z']**2), axis=1)\ntrain.drop(['X','Y','Z'], axis=1, inplace=True)\ndel dipole_moments","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# potential energy of molecules\npotential_energy = pd.read_csv('../input/champs-scalar-coupling/potential_energy.csv')\ntrain = pd.merge(train, potential_energy, how='left',\n              left_on='molecule_name',\n              right_on='molecule_name')\ndel potential_energy","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# dipole_moment and potential_energy are molecular features,\n# so let's get some molecular features and use them to estimate\n# count the total no. atoms in the molecule and the number of each type\natom_cnt = structures['molecule_name'].value_counts().reset_index(level=0)\natom_cnt.rename({'index': 'molecule_name', 'molecule_name': 'atom_count'}, axis=1, inplace=True)\ntrain = pd.merge(train, atom_cnt, how='left', on='molecule_name')\ntest = pd.merge(test, atom_cnt, how='left', on='molecule_name')\ndel atom_cnt\n\n# (the following should be put in a loop)\nH_cnt = structures['molecule_name'][structures['atom']=='H'].value_counts().reset_index(level=0)\nH_cnt.rename({'index': 'molecule_name', 'molecule_name': 'H_count'}, axis=1, inplace=True)\ntrain = pd.merge(train, H_cnt, how='left', on='molecule_name')\ntest = pd.merge(test, H_cnt, how='left', on='molecule_name')\ndel H_cnt\n\nC_cnt = structures['molecule_name'][structures['atom']=='C'].value_counts().reset_index(level=0)\nC_cnt.rename({'index': 'molecule_name', 'molecule_name': 'C_count'}, axis=1, inplace=True)\ntrain = pd.merge(train, C_cnt, how='left', on='molecule_name')\ntest = pd.merge(test, C_cnt, how='left', on='molecule_name')\ndel C_cnt\n\nO_cnt = structures['molecule_name'][structures['atom']=='O'].value_counts().reset_index(level=0)\nO_cnt.rename({'index': 'molecule_name', 'molecule_name': 'O_count'}, axis=1, inplace=True)\ntrain = pd.merge(train, O_cnt, how='left', on='molecule_name')\ntest = pd.merge(test, O_cnt, how='left', on='molecule_name')\ndel O_cnt\n\nN_cnt = structures['molecule_name'][structures['atom']=='N'].value_counts().reset_index(level=0)\nN_cnt.rename({'index': 'molecule_name', 'molecule_name': 'N_count'}, axis=1, inplace=True)\ntrain = pd.merge(train, N_cnt, how='left', on='molecule_name')\ntest = pd.merge(test, N_cnt, how='left', on='molecule_name')\ndel N_cnt\n\nF_cnt = structures['molecule_name'][structures['atom']=='F'].value_counts().reset_index(level=0)\nF_cnt.rename({'index': 'molecule_name', 'molecule_name': 'F_count'}, axis=1, inplace=True)\ntrain = pd.merge(train, F_cnt, how='left', on='molecule_name')\ntest = pd.merge(test, F_cnt, how='left', on='molecule_name')\ndel F_cnt\n\ntrain = train.fillna(0)\ntest = test.fillna(0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get molecular distance props from previous kernel\n# https://www.kaggle.com/robertburbidge/distance-features\nmol_props = ['molecule_dist_mean', 'molecule_dist_std', 'molecule_dist_skew',\n             'molecule_dist_kurt', 'meanx', 'meany', 'meanz', 'meanxH', 'meanyH', 'meanzH',\n             'meanxC', 'meanyC', 'meanzC','meanxN', 'meanyN', 'meanzN','meanxO', 'meanyO', 'meanzO',\n             'meanxF', 'meanyF', 'meanzF']\n\ntrain[mol_props] = train_dist[mol_props]\ntest[mol_props] = test_dist[mol_props]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# electronegativity (this could be better informed by physical chemistry)\natoms = ['H', 'C', 'N', 'O', 'F']\nelectronegativity = {'H': 2.2, 'C': 2.55, 'N': 3.04, 'O': 3.44, 'F': 3.98}\nfor atom in atoms:\n    meanx = (structures[structures['atom']==atom].groupby('molecule_name')['x'].apply(np.mean) *\n             electronegativity[atom]).reset_index()\n    meanx.rename({'x': 'meanxe' + atom}, axis=1, inplace=True)\n    train = pd.merge(train, meanx, how='left', on='molecule_name')\n    test = pd.merge(test, meanx, how='left', on='molecule_name')\n\n    meany = (structures[structures['atom']==atom].groupby('molecule_name')['y'].apply(np.mean) *\n             electronegativity[atom]).reset_index()\n    meany.rename({'y': 'meanye' + atom}, axis=1, inplace=True)\n    train = pd.merge(train, meany, how='left', on='molecule_name')\n    test = pd.merge(test, meany, how='left', on='molecule_name')\n\n    meanz = (structures[structures['atom']==atom].groupby('molecule_name')['z'].apply(np.mean) *\n             electronegativity[atom]).reset_index()\n    meanz.rename({'z': 'meanze' + atom}, axis=1, inplace=True)\n    train = pd.merge(train, meanz, how='left', on='molecule_name')\n    test = pd.merge(test, meanz, how='left', on='molecule_name')\n    \ndel meanx, meany, meanz, structures","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# train, predict at molecular level (could have done this earlier but you may want to pull in extra features \n# and aggregate them before this stage, so I left this step to just prior to modelling)\ntrain_mol = train.drop_duplicates(subset=['molecule_name'])\ntest_mol = test.drop_duplicates(subset=['molecule_name'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# molecular features for predicting mu\npred_vars = [v for v in train.columns if v not in ['id', 'molecule_name', 'atom_index_0', 'atom_index_1',\n                                                   'type', 'dipole_moment', 'potential_energy',\n                                                   'scalar_coupling_constant']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# train-val split by molecule_name (since train and test data have disjoint molecules)\nmolecule_names = pd.DataFrame(permutation(train_mol['molecule_name'].unique()),columns=['molecule_name'])\nnm = molecule_names.shape[0]\nntrn = int(0.9*nm)\nnval = int(0.1*nm)\n\ntmp_train = pd.merge(train_mol, molecule_names[0:ntrn], how='right', on='molecule_name')\ntmp_val = pd.merge(train_mol, molecule_names[ntrn:nm], how='right', on='molecule_name')\n\nX_train = tmp_train[pred_vars]\nX_val = tmp_val[pred_vars]\ny_train = tmp_train['dipole_moment']\ny_val = tmp_val['dipole_moment']\ndel tmp_train, tmp_val","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# data for LightGBM\ntrain_data = lightgbm.Dataset(X_train, label=y_train)\nval_data = lightgbm.Dataset(X_val, label=y_val)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# heuristic parameters for LightGBM\nparams = { 'objective': 'regression_l1',\n           'learning_rate': 0.1,\n           'num_leaves': 255,\n           'num_threads': -1,\n           'bagging_fraction': 0.9,\n           'bagging_freq': 10,\n           'feature_fraction': 0.9,\n           'lambda_l1': 10.0,\n           'max_bin': 255,\n           'min_child_samples': 50,\n           }","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# training\nmodel = lightgbm.train(params,\n                       train_data,\n                       valid_sets=[train_data, val_data], verbose_eval=500,\n                       num_boost_round=4000,\n                       early_stopping_rounds=100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# validation\npred_val = model.predict(X_val)\npred_median = np.full(y_val.shape, np.median(y_val))\nprint(metrics.mean_absolute_error(y_val, pred_val) / metrics.mean_absolute_error(y_val, pred_median))\nplt.scatter(y_val, pred_val)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# train & test predictions\npred_train = model.predict(train_mol[pred_vars])\npred_test = model.predict(test_mol[pred_vars])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# save these for inputs for further modelling\ntrain_dipole_moment = pd.DataFrame(train_mol['molecule_name'])\ntrain_dipole_moment['dipole_moment_pred'] = pred_train\ntrain_dipole_moment.to_csv('train_dipole_moment.csv', index=False)\ntest_dipole_moment = pd.DataFrame(test_mol['molecule_name'])\ntest_dipole_moment['dipole_moment_pred'] = pred_test\ntest_dipole_moment.to_csv('test_dipole_moment.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# molecular features for predicting potential energy\npred_vars = [v for v in train.columns if v not in ['id', 'molecule_name', 'atom_index_0', 'atom_index_1',\n                                                   'type', 'dipole_moment', 'potential_energy',\n                                                   'scalar_coupling_constant']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# train-val split by molecule_name\nmolecule_names = pd.DataFrame(permutation(train_mol['molecule_name'].unique()),columns=['molecule_name'])\nnm = molecule_names.shape[0]\nntrn = int(0.9*nm)\nnval = int(0.1*nm)\n\ntmp_train = pd.merge(train_mol, molecule_names[0:ntrn], how='right', on='molecule_name')\ntmp_val = pd.merge(train_mol, molecule_names[ntrn:nm], how='right', on='molecule_name')\n\nX_train = tmp_train[pred_vars]\nX_val = tmp_val[pred_vars]\ny_train = tmp_train['potential_energy']\ny_val = tmp_val['potential_energy']\ndel tmp_train, tmp_val","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# data for LightGBM\ntrain_data = lightgbm.Dataset(X_train, label=y_train)\nval_data = lightgbm.Dataset(X_val, label=y_val)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# heuristic parameters for LightGBM\nparams = { 'objective': 'regression_l1',\n           'learning_rate': 0.1,\n           'num_leaves': 255,\n           'num_threads': -1,\n           'bagging_fraction': 0.9,\n           'bagging_freq': 10,\n           'feature_fraction': 0.9,\n           'lambda_l1': 10.0,\n           'max_bin': 255,\n           'min_child_samples': 50,\n           }","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# training\nmodel = lightgbm.train(params,\n                       train_data,\n                       valid_sets=[train_data, val_data], verbose_eval=100,\n                       num_boost_round=1000,\n                       early_stopping_rounds=100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# validation\npred_val = model.predict(X_val)\npred_median = np.full(y_val.shape, np.median(y_val))\nprint(metrics.mean_absolute_error(y_val, pred_val) / metrics.mean_absolute_error(y_val, pred_median))\nplt.scatter(y_val, pred_val)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# train & test predictions\npred_train = model.predict(train_mol[pred_vars])\npred_test = model.predict(test_mol[pred_vars])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# save these for inputs for further modelling\ntrain_potential_energy = pd.DataFrame(train_mol['molecule_name'])\ntrain_potential_energy['potential_energy_pred'] = pred_train\ntrain_potential_energy.to_csv('train_potential_energy.csv', index=False)\ntest_potential_energy = pd.DataFrame(test_mol['molecule_name'])\ntest_potential_energy['potential_energy_pred'] = pred_test\ntest_potential_energy.to_csv('test_potential_energy.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# example of using these features to improve (slightly) on previous performance with dist features\ntrain.drop(['potential_energy', 'dipole_moment'], axis=1, inplace=True)\ntrain = pd.merge(train, train_dipole_moment, how='left', on='molecule_name')\ntrain = pd.merge(train, train_potential_energy, how='left', on='molecule_name')\ntest = pd.merge(test, test_dipole_moment, how='left', on='molecule_name')\ntest = pd.merge(test, test_potential_energy, how='left', on='molecule_name')\ndel train_dipole_moment, train_potential_energy, test_dipole_moment, test_potential_energy","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get remaining features from previous kernel\naddfeats = [v for v in train_dist.columns if v not in mol_props]\ntrain[addfeats] = train_dist[addfeats]\ntest[addfeats] = test_dist[addfeats]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# https://www.kaggle.com/artgor/artgor-utils\ndef reduce_mem_usage(df, verbose=True):\n    numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']\n    start_mem = df.memory_usage().sum() / 1024 ** 2\n    for col in df.columns:\n        col_type = df[col].dtypes\n        if col_type in numerics:\n            c_min = df[col].min()\n            c_max = df[col].max()\n            if str(col_type)[:3] == 'int':\n                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n                    df[col] = df[col].astype(np.int8)\n                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n                    df[col] = df[col].astype(np.int16)\n                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n                    df[col] = df[col].astype(np.int32)\n                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n                    df[col] = df[col].astype(np.int64)\n            else:\n                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:\n                    df[col] = df[col].astype(np.float16)\n                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:\n                    df[col] = df[col].astype(np.float32)\n                else:\n                    df[col] = df[col].astype(np.float64)\n    end_mem = df.memory_usage().sum() / 1024 ** 2\n    if verbose: print('Mem. usage decreased to {:5.2f} Mb ({:.1f}% reduction)'.format(end_mem, 100 * (\n                start_mem - end_mem) / start_mem))\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = reduce_mem_usage(train)\ntest = reduce_mem_usage(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# features for prediction\npred_vars = [v for v in train.columns if v not in ['id', 'molecule_name', 'scalar_coupling_constant']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# train-val split by molecule_name\nmolecule_names = pd.DataFrame(permutation(train['molecule_name'].unique()),columns=['molecule_name'])\nnm = molecule_names.shape[0]\nntrn = int(0.9*nm)\nnval = int(0.1*nm)\n\ntmp_train = pd.merge(train, molecule_names[0:ntrn], how='right', on='molecule_name')\ntmp_val = pd.merge(train, molecule_names[ntrn:nm], how='right', on='molecule_name')\n\nX_train = tmp_train[pred_vars]\nX_val = tmp_val[pred_vars]\ny_train = tmp_train['scalar_coupling_constant']\ny_val = tmp_val['scalar_coupling_constant']\ndel tmp_train, tmp_val","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# heuristic parameters for LightGBM\nparams = { 'objective': 'regression_l1',\n           'learning_rate': 0.1,\n           'num_leaves': 1023,\n           'num_threads': -1,\n           'bagging_fraction': 0.5,\n           'bagging_freq': 1,\n           'feature_fraction': 0.9,\n           'lambda_l1': 10.0,\n           'max_bin': 255,\n           'min_child_samples': 15,\n           }","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# categorical features (note that these are already integer-coded)\ncat_feats = ['type', 'type_0', 'type_1', 'atom_0l', 'atom_0r', 'atom_1l', 'atom_1r']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# data for lightgbm\ntrain_data = lightgbm.Dataset(X_train, label=y_train, categorical_feature=cat_feats)\nval_data = lightgbm.Dataset(X_val, label=y_val, categorical_feature=cat_feats)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# training\nmodel = lightgbm.train(params,\n                       train_data,\n                       valid_sets=[train_data, val_data], verbose_eval=500,\n                       num_boost_round=4000,\n                       early_stopping_rounds=100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# evaluation metric for validation\n# https://www.kaggle.com/abhishek/competition-metric\ndef metric(df, preds):\n    df[\"prediction\"] = preds\n    maes = []\n    for t in df.type.unique():\n        y_true = df[df.type==t].scalar_coupling_constant.values # column 1 is the target whatever the y_var\n        y_pred = df[df.type==t].prediction.values\n        mae = np.log(metrics.mean_absolute_error(y_true, y_pred))\n        maes.append(mae)\n    return np.mean(maes)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# validation metric\npreds = model.predict(X_val)\nmetric(pd.concat([X_val, y_val], axis=1), preds)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# submission\npreds_sub = model.predict(test[pred_vars])\nsub['scalar_coupling_constant'] = preds_sub\nsub.to_csv('submission_feats_dist_mol01.csv', index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}