{"cells":[{"metadata":{},"cell_type":"markdown","source":"There is a kernel [here](https://www.kaggle.com/borisdee/predicting-mulliken-charges-with-acsf-descriptors) that describes predicting the Mulliken charges on the test set using atom-centered symmetry functions (ACSFs) and also gives an introduction and references to ACSFs. Estimated Mulliken charges may well be useful for predicting the scalar coupling constant as applied in this [kernel](https://www.kaggle.com/robertburbidge/using-estimated-mulliken-charges).\n\nHere, I use ACSFs to directly predict the scalar coupling constant. I use the packages Dscribe and ASE to handle the chemistry. These are worth investigating for this competition.\n\nhttps://github.com/SINGROUP/dscribe\n\nhttps://wiki.fysik.dtu.dk/ase/"},{"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":{"trusted":true},"cell_type":"code","source":"!pip install dscribe","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"from sklearn import metrics\nimport lightgbm\nfrom sklearn.model_selection import KFold\nfrom sklearn.preprocessing import LabelEncoder\nfrom dscribe.descriptors import ACSF\nfrom ase.io import read as ase_read","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# settings for cross-validation and LightGBM training\nnfolds = 5\nniters = 4000","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/train.csv')\ntest = pd.read_csv('../input/test.csv')\nsub = pd.read_csv('../input/sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Setting up the ACSF descriptor\n# see the kernel referenced above as well as the Dscribe docs:\n# https://github.com/SINGROUP/dscribe/blob/master/docs/tutorials/acsf.html\n# chemists can add more functions and tune the parameters here\nacsf = ACSF(\n    species=['H', 'C', 'N', 'O', 'F'],\n    rcut=3.0,\n    g2_params=[(0.4, 0.2),(0.4, 0.5),(0.4, 1.0),(0.5, 2.0),(0.5, 3.0),(0.5, 4.0)],\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# calculate ACSFs\ndef calc_acsf(df):\n    acsf_arr = np.zeros([df.shape[0], 70])\n    for i in range(df.shape[0]):\n        molecule_name = df['molecule_name'].iloc[i]\n        atoms = ase_read('../input/structures/' + molecule_name + '.xyz')\n        acsfi = acsf.create(atoms, positions = [df['atom_index_0'].iloc[i], df['atom_index_1'].iloc[i]], n_jobs=4)\n        acsf_arr[i, :] = np.reshape(acsfi, [70])\n    df = pd.concat([df.reset_index(drop=True), pd.DataFrame(acsf_arr)], axis=1)\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = calc_acsf(train)\ntest = calc_acsf(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# predictive vars for LightGBM\npred_vars = [v for v in train_sample.columns if v not in ['id', 'molecule_name', 'scalar_coupling_constant',\n                                                   'atom_index_0', 'atom_index_1']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# encode type\ncat_feats = ['type']\nfor f in cat_feats:\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': 255,\n           'min_data_in_leaf': 100,\n           'max_depth': 10,\n           'num_threads': -1,\n           'bagging_fraction': 0.5,\n           'bagging_freq': 1,\n           'feature_fraction': 0.9,\n           'lambda_l1': 10.0,\n           'lambda_l2': 10.0,\n           'max_bin': 255,\n           'verbosity': -1\n           }","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":"# KFold cross-validation by molecule, separate models for each type\nkf = KFold(n_splits=nfolds)\nmolecule_names = train_sample['molecule_name'].unique()\n\npreds = np.empty([train_sample.shape[0]])\npreds_sub = np.zeros([sub.shape[0]])\ntypes = train_sample['type'].unique()\n\nfor train_mol_idx, val_mol_idx in kf.split(molecule_names):\n    #\n    train_idx = pd.merge(train_sample[['molecule_name']].reset_index(),\n                         pd.DataFrame(molecule_names[train_mol_idx]), how='right',\n                         left_on='molecule_name', right_on=0)['index']\n    val_idx = pd.merge(train_sample[['molecule_name']].reset_index(),\n                         pd.DataFrame(molecule_names[val_mol_idx]), how='right',\n                       left_on='molecule_name', right_on=0)['index']\n    #\n    for type in types:\n        train_data = lightgbm.Dataset(train_sample.iloc[train_idx,:].loc[train_sample.iloc[train_idx,:]['type']==type][pred_vars],\n                                      label=train_sample.iloc[train_idx,:].loc[train_sample.iloc[train_idx,:]['type']==type]['scalar_coupling_constant'],\n                                      categorical_feature=cat_feats)\n        val_data = lightgbm.Dataset(train_sample.iloc[val_idx,:].loc[train_sample.iloc[val_idx,:]['type']==type][pred_vars],\n                                      label=train_sample.iloc[val_idx,:].loc[train_sample.iloc[val_idx,:]['type']==type]['scalar_coupling_constant'],\n                                      categorical_feature=cat_feats)\n        #\n        # training\n        model = lightgbm.train(params,\n                               train_data,\n                               valid_sets=[train_data, val_data], verbose_eval=int(niters/8),\n                               num_boost_round=niters,\n                               early_stopping_rounds=int(niters/40))\n        #\n        tmp_idx = val_idx.values[train_sample.iloc[val_idx, :]['type'] == type]\n        preds[tmp_idx] =\\\n            model.predict(train_sample.iloc[val_idx,:].loc[train_sample.iloc[val_idx,:]['type']==type][pred_vars])\n        #\n        preds_sub[test['type']==type] = preds_sub[test['type']==type] + \\\n            model.predict(test.loc[test['type']==type,:][pred_vars])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# validation performance\nprint(metric(pd.concat([train_sample[pred_vars], train_sample['scalar_coupling_constant']], axis=1), preds))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# submission\nsub['scalar_coupling_constant'] = preds_sub\nsub.to_csv('submission_acsf01.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}