{"cells":[{"metadata":{},"cell_type":"markdown","source":"The magnetic shielding tensor is provided for the train set but not the test set. Perhaps if we can impute the tensors on the test set we can obtain better predictions of the scalar coupling constant. The magnetic shielding tensor is described here:\nhttps://www.ncbi.nlm.nih.gov/pmc/articles/PMC3058154/#S3title. "},{"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 sklearn.preprocessing import LabelEncoder","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":"# 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":"# previously constructed distance features\n# https://www.kaggle.com/robertburbidge/distance-features\ntrain_dist = pd.read_csv('train_dist.csv')\ntest_dist = pd.read_csv('test_dist.csv')\ntrain = pd.merge(train.drop(['atom_index_0', 'atom_index_1', 'type'], axis=1), train_dist, how='left', on='id')\ntest = pd.merge(test.drop(['atom_index_0', 'atom_index_1', 'type'], axis=1), test_dist, how='left', on='id')\ndist_columns = list(train_dist.columns)\ndel train_dist, test_dist","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# atomic properties\n# https://www.lenntech.com/periodic-chart-elements/\natomic_radius = {'H': 0.38, 'C': 0.77, 'N': 0.75, 'O': 0.73, 'F': 0.71, np.nan: 0}\natomic_number = {'H': 1, 'C': 6, 'N': 7, 'O': 8, 'F': 9, np.nan: 0}\natomic_mass = {'H': 1.0079, 'C': 12.0107, 'N': 14.0067, 'O': 15.9994, 'F': 18.9984, np.nan: 0}\nvanderwaalsradius = {'H': 120, 'C': 185, 'N': 154, 'O': 140, 'F': 135, np.nan: 0}\ncovalenzradius = {'H': 30, 'C': 77, 'N': 70, 'O': 66, 'F': 58, np.nan: 0}\nelectronegativity = {'H': 2.2, 'C': 2.55, 'N': 3.04, 'O': 3.44, 'F': 3.98, np.nan: 0}\nionization_energy = {'H': 13.5984, 'C': 11.2603, 'N': 14.5341, 'O': 13.6181, 'F': 17.4228, np.nan: np.inf}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def atom_props(df, suffix):\n    df['atomic_radius' + suffix] = df['atom_' + suffix].apply(lambda x: atomic_radius[x])\n    df['atomic_protons' + suffix] = df['atom_' + suffix].apply(lambda x: atomic_number[x])\n    df['atomic_mass' + suffix] = df['atom_' + suffix].apply(lambda x: atomic_mass[x])\n    df['vanderwaalsradius' + suffix] = df['atom_' + suffix].apply(lambda x: vanderwaalsradius[x])\n    df['covalenzradius' + suffix] = df['atom_' + suffix].apply(lambda x: covalenzradius[x])\n    df['electronegativity' + suffix] = df['atom_' + suffix].apply(lambda x: electronegativity[x])\n    df['ionization_energy' + suffix] = df['atom_' + suffix].apply(lambda x: ionization_energy[x])\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# atom_0, atom_1\ntrain = pd.merge(train, structures[['molecule_name', 'atom_index', 'atom']], how='left',\n                 left_on=['molecule_name', 'atom_index_0'],\n                 right_on=['molecule_name', 'atom_index']).rename({'atom': 'atom_0'}, axis=1)\ntrain = pd.merge(train, structures[['molecule_name', 'atom_index', 'atom']], how='left',\n                 left_on=['molecule_name', 'atom_index_1'],\n                 right_on=['molecule_name', 'atom_index']).rename({'atom': 'atom_1'}, axis=1)\ntest = pd.merge(test, structures[['molecule_name', 'atom_index', 'atom']], how='left',\n                 left_on=['molecule_name', 'atom_index_0'],\n                 right_on=['molecule_name', 'atom_index']).rename({'atom': 'atom_0'}, axis=1)\ntest = pd.merge(test, structures[['molecule_name', 'atom_index', 'atom']], how='left',\n                 left_on=['molecule_name', 'atom_index_1'],\n                 right_on=['molecule_name', 'atom_index']).rename({'atom': 'atom_1'}, axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# neighbours by index (replace the integer codes from {train|test}_dist with the original atom chars)\ndef lr(df):\n    df.drop(['atom_0l', 'atom_0r', 'atom_1l', 'atom_1r'], axis=1, inplace=True)\n    df['atom_index_0l'] = df['atom_index_0'].apply(lambda i: max(i - 1, 0))\n    tmp = df[['atom_index_0', 'atom_count']]\n    df['atom_index_0r'] = tmp.apply(lambda row: min(row['atom_index_0'] + 1, row['atom_count']), axis=1)\n    df = pd.merge(df, structures[['molecule_name', 'atom_index', 'atom']], how='left',\n                     left_on=['molecule_name', 'atom_index_0l'],\n                     right_on=['molecule_name', 'atom_index']).rename({'atom': 'atom_0l'}, axis=1)\n    df = pd.merge(df, structures[['molecule_name', 'atom_index', 'atom']], how='left',\n                     left_on=['molecule_name', 'atom_index_0r'],\n                     right_on=['molecule_name', 'atom_index']).rename({'atom': 'atom_0r'}, axis=1)\n    df['atom_index_1l'] = df['atom_index_1'].apply(lambda i: max(i - 1, 0))\n    tmp = df[['atom_index_1', 'atom_count']]\n    df['atom_index_1r'] = tmp.apply(lambda row: min(row['atom_index_1'] + 1, row['atom_count']), axis=1)\n    df = pd.merge(df, structures[['molecule_name', 'atom_index', 'atom']], how='left',\n                     left_on=['molecule_name', 'atom_index_1l'],\n                     right_on=['molecule_name', 'atom_index']).rename({'atom': 'atom_1l'}, axis=1)\n    df = pd.merge(df, structures[['molecule_name', 'atom_index', 'atom']], how='left',\n                     left_on=['molecule_name', 'atom_index_1r'],\n                     right_on=['molecule_name', 'atom_index']).rename({'atom': 'atom_1r'}, axis=1)\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = lr(train)\ntest = lr(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get atomic properties of both atoms and their neighbours\ntrain = atom_props(train, '0')\ntrain = atom_props(train, '0l')\ntrain = atom_props(train, '0r')\ntrain = atom_props(train, '1')\ntrain = atom_props(train, '1l')\ntrain = atom_props(train, '1r')\ntest = atom_props(test, '0')\ntest = atom_props(test, '0l')\ntest = atom_props(test, '0r')\ntest = atom_props(test, '1')\ntest = atom_props(test, '1l')\ntest = atom_props(test, '1r')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"atom_feats = [v for v in train.columns if v not in (dist_columns +\n              ['molecule_name', 'scalar_coupling_constant', 'atom_index_x', 'atom_index_y'])]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# save these for future use\ntrain[['id'] + atom_feats].to_csv('train_atom_feats.csv', index=False)\ntest[['id'] + atom_feats].to_csv('test_atom_feats.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# drop duplicate columns\ntrain.drop(['atom_index_x', 'atom_index_y'], axis=1, inplace=True)\ntest.drop(['atom_index_x', 'atom_index_y'], axis=1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# https://www.kaggle.com/c/champs-scalar-coupling/discussion/96655#latest-558745\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                c_prec = df[col].apply(lambda x: np.finfo(x).precision).max()\n                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max and c_prec == np.finfo(np.float16).precision:\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 and c_prec == np.finfo(np.float32).precision:\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 * (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":"# magnetic_shielding_tensors\nmagnetic_shielding_tensors = pd.read_csv('../input/magnetic_shielding_tensors.csv')\ntrain = pd.merge(train, magnetic_shielding_tensors, how='left',\n              left_on=['molecule_name', 'atom_index_0'],\n              right_on=['molecule_name', 'atom_index'])\n\ntrain.rename(columns = { \"XX\": \"XX_0\", \"YX\": \"YX_0\", \"ZX\": \"ZX_0\",\n                         \"XY\": \"XY_0\", \"YY\": \"YY_0\", \"ZY\": \"ZY_0\",\n                         \"XZ\": \"XZ_0\", \"YZ\": \"YZ_0\", \"ZZ\": \"ZZ_0\" }, inplace=True)\n\ntrain = pd.merge(train, magnetic_shielding_tensors, how='left',\n              left_on=['molecule_name', 'atom_index_1'],\n              right_on=['molecule_name', 'atom_index'])\n\ntrain.rename(columns = { \"XX\": \"XX_1\", \"YX\": \"YX_1\", \"ZX\": \"ZX_1\",\n                         \"XY\": \"XY_1\", \"YY\": \"YY_1\", \"ZY\": \"ZY_1\",\n                         \"XZ\": \"XZ_1\", \"YZ\": \"YZ_1\", \"ZZ\": \"ZZ_1\" }, inplace=True)\ndel magnetic_shielding_tensors","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y_vars = ['XX_0', 'YX_0', 'ZX_0', 'XY_0', 'YY_0', 'ZY_0',\n             'XZ_0', 'YZ_0', 'ZZ_0', 'XX_1', 'YX_1', 'ZX_1', 'XY_1',\n             'YY_1', 'ZY_1', 'XZ_1', 'YZ_1', 'ZZ_1']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# features for prediction of magnetic shielding tensor\npred_vars = [v for v in train.columns if v not in ['id', 'molecule_name', 'atom_0', 'atom_1',\n                                                   'scalar_coupling_constant', 'atom_index_x', 'atom_index_y'] + y_vars]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# encode categorical features as integers for LightGBM\ncat_feats = ['type', 'type_0', 'type_1', 'atom_0l', 'atom_0r', 'atom_1l', 'atom_1r']\ncat_feats_to_encode = ['atom_0l', 'atom_0r', 'atom_1l', 'atom_1r']\nfor f in cat_feats_to_encode:\n    lbl = LabelEncoder()\n    lbl.fit(list(train[f].values) + list(test[f].values))\n    train[f] = lbl.transform(list(train[f].values))\n    test[f] = lbl.transform(list(test[f].values))","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":"# loop through the magnetic shielding tensor components, model, predict\n# this takes a long time and the quality of the predictions is low\nfor i in np.arange(0, len(y_vars)):\n    #\n    # train-val split by molecule_name\n    molecule_names = pd.DataFrame(permutation(train['molecule_name'].unique()),columns=['molecule_name'])\n    nm = molecule_names.shape[0]\n    ntrn = int(0.9*nm)\n    nval = int(0.1*nm)\n    #\n    tmp_train = pd.merge(train, molecule_names[0:ntrn], how='right', on='molecule_name')\n    tmp_val = pd.merge(train, molecule_names[ntrn:nm], how='right', on='molecule_name')\n    #\n    X_train = tmp_train[pred_vars]\n    X_val = tmp_val[pred_vars]\n    y_train = tmp_train[y_vars[i]]\n    y_val = tmp_val[y_vars[i]]\n    del tmp_train, tmp_val\n    #\n    train_data = lightgbm.Dataset(X_train, label=y_train, categorical_feature=cat_feats)\n    val_data = lightgbm.Dataset(X_val, label=y_val, categorical_feature=cat_feats)\n    #\n    # training\n    model = 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=20)\n    #\n    # validation\n    pred_val = model.predict(X_val)\n    pred_median = np.full(y_val.shape, np.median(y_val))\n    print(metrics.mean_absolute_error(y_val, pred_val) / metrics.mean_absolute_error(y_val, pred_median))\n    #\n    train[y_vars[i]] = model.predict(train[pred_vars])\n    test = pd.concat([test, pd.DataFrame(model.predict(test[pred_vars]))], axis=1)\n    test.rename(columns={0: y_vars[i]}, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# save features for future use\ntrain[['id'] + y_vars].to_csv('train_mst.csv', index=False)\ntest[['id'] + y_vars].to_csv('test_mst.csv', index=False)","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":"# example use of the imputed tensors to attempt to improve on my previous kernel:\n# https://www.kaggle.com/robertburbidge/using-estimated-mulliken-charges\n# hence, we drop the atom_feats for now as they weren't used previously and we want to see if the tensors\n# have any predictive power\ntrain.drop(atom_feats, axis=1, inplace=True)\ntest.drop(atom_feats, axis=1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# previously imputed molecular features\n# https://www.kaggle.com/robertburbidge/imputing-molecular-features\ntrain_dipole_moment = pd.read_csv('../input/imputing-molecular-features/train_dipole_moment.csv')\ntest_dipole_moment = pd.read_csv('../input/imputing-molecular-features/test_dipole_moment.csv')\ntrain = pd.merge(train, train_dipole_moment, how='left', on='molecule_name')\ntest = pd.merge(test, test_dipole_moment, how='left', on='molecule_name')\ntrain_potential_energy = pd.read_csv('../input/imputing-molecular-features/train_potential_energy.csv')\ntest_potential_energy = pd.read_csv('../input/imputing-molecular-features/test_potential_energy.csv')\ntrain = pd.merge(train, train_potential_energy, how='left', on='molecule_name')\ntest = pd.merge(test, test_potential_energy, how='left', on='molecule_name')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# mulliken charges with Open Babel\n# https://www.kaggle.com/asauve/v7-estimation-of-mulliken-charges-with-open-babel\ntrain_ob_charges = pd.read_csv('../input/v7-estimation-of-mulliken-charges-with-open-babel/train_ob_charges.csv')\ntest_ob_charges = pd.read_csv('../input/v7-estimation-of-mulliken-charges-with-open-babel/test_ob_charges.csv')\ntrain = pd.merge(train, train_ob_charges[['molecule_name', 'atom_index', 'eem']], how='left',\n         left_on=['molecule_name', 'atom_index_0'], right_on=['molecule_name', 'atom_index']).\\\n    rename({'eem': 'eem0'}, axis=1)\ntrain = pd.merge(train, train_ob_charges[['molecule_name', 'atom_index', 'eem']], how='left',\n         left_on=['molecule_name', 'atom_index_1'], right_on=['molecule_name', 'atom_index']).\\\n    rename({'eem': 'eem1'}, axis=1)\ntest = pd.merge(test, test_ob_charges[['molecule_name', 'atom_index', 'eem']], how='left',\n         left_on=['molecule_name', 'atom_index_0'], right_on=['molecule_name', 'atom_index']).\\\n    rename({'eem': 'eem0'}, axis=1)\ntest = pd.merge(test, test_ob_charges[['molecule_name', 'atom_index', 'eem']], how='left',\n         left_on=['molecule_name', 'atom_index_1'], right_on=['molecule_name', 'atom_index']).\\\n    rename({'eem': 'eem1'}, axis=1)","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 for lightgbm (already integer-encoded)\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\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\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_mol_mc_mst01.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}