{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Basic feature generation with molmod library\n## Inroduction / prerequisites\n\nThe [molmod](http://molmod.github.io/molmod/index.html) library is \n> > a Python library with many compoments that are useful to write molecular modeling programs.\n\nThe library provides fairly easy way to extract data from `*.xyz` files, starting from basic features like bonds and distances between atoms up to advanced features like dihedral angles and advanced pattern matching using graphs.\n\nHere I present an example how to extract some basic features from [CHAMPS competition](https://www.kaggle.com/c/champs-scalar-coupling).\n\nIn order to use run code from this kernel You need molmod library installed.\n\n## Preparations\n\nLet's load libraries"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport re\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.model_selection import GroupKFold\nfrom sklearn.ensemble import RandomForestRegressor\nimport matplotlib.pyplot as plt\nimport lightgbm as lgb\nimport gc\nimport seaborn as sns\nimport molmod\nimport warnings\nimport multiprocessing\nwarnings.simplefilter(action='ignore', category=FutureWarning)\nimport os\nprint(os.listdir('../input'))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Check wehere our competition data is stored."},{"metadata":{"trusted":true},"cell_type":"code","source":"data_dir = '../input/champs-scalar-coupling/' if 'champs-scalar-coupling' in os.listdir('../input') else '../input/'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Evaluation function for predictions."},{"metadata":{"trusted":true},"cell_type":"code","source":"def eval_fn(y_pred, y_test, c_type):\n    diff = np.full((y_pred.shape[0], 2), 1e-9)\n    diff[:,1] = np.abs(y_pred - y_test)\n    step_1 = pd.DataFrame(\n        {'diff': np.amax(diff, 1), 'type': c_type}\n    ).groupby(\n        'type'\n    ).mean()\n    return np.sum(np.log(step_1['diff'])) / step_1.shape[0]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Load and preprocess data"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ntrain = pd.read_csv(data_dir + 'train.csv')\ntest = pd.read_csv(data_dir + 'test.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Add basic features - atomic number for atom_1 (atom_0 is always hydrogen) and number of bonds between target atoms."},{"metadata":{"trusted":true},"cell_type":"code","source":"atoms = {'H': 1, 'C': 6, 'N': 7, 'O': 8, 'F': 9}\ndef map_atom(atom):\n    return atoms[atom]\ntrain['atom_1'] = train['type'].astype(str).str[3].apply(map_atom)\ntest['atom_1'] = test['type'].astype(str).str[3].apply(map_atom)\ntrain['n_bonds'] = train['type'].astype(str).str[0].astype(np.int)\ntest['n_bonds'] = test['type'].astype(str).str[0].astype(np.int)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Some helper functions"},{"metadata":{"trusted":true},"cell_type":"code","source":"def pad_1d_array(a, length, value = np.nan):\n    '''This function will right pad a numpy array `a` to `length` with `value`'''\n    if a.shape[0] >= length:\n        return a\n    return np.pad(a, (0, length - a.shape[0]), 'constant', constant_values = value)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_neighbor_indices(molecule, index, limit = 4):\n    '''Get neighbor indices from molmod.molecules.Molecules object for `index`.\n    As sometimes the max bond limit is exceeded add limiting value'''\n    neighbor_indices = list(molecule.graph.neighbors[index])\n    while len(neighbor_indices) > limit:\n        distances_tmp = molecule.distance_matrix[index, neighbor_indices]\n        neighbor_indices.remove(neighbor_indices[np.argmax(distances_tmp)])\n    return neighbor_indices","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Feature generation\n\nMain feature generating function"},{"metadata":{"trusted":true},"cell_type":"code","source":"def gen_features(row):\n    '''Generate features for train/test entry, using molmod.molecules.Molecule class''' \n    # This is some hashing in order not to repeat loading molecule data from file\n    # We expect the following indices in `row` 0: molecule name; 1: atom index 0; 2: atom index 1\n    if row[0] != gen_features.molecule_name:\n        # Load molecule in hashed var\n        gen_features.molecule = molmod.molecules.Molecule.from_file(f'{data_dir}/structures/{row[0]}.xyz')\n        # Remember loaded molecule name\n        gen_features.molecule_name = row[0]\n        # Generate graph\n        gen_features.molecule.set_default_graph()\n    # Distance between target atoms\n    distance = gen_features.molecule.distance_matrix[row[1], row[2]]\n    # Neighbor indices for both target atoms (up to 4 - default limit)\n    neighbors_idxs_0 = get_neighbor_indices(gen_features.molecule, row[1], 1)\n    neighbors_idxs_1 = get_neighbor_indices(gen_features.molecule, row[2], 4)\n    # Get atomic numbers for neighbor atoms\n    # Atom 0 is always hydrogen and has no more than 1 neighbor\n    neighbor_atoms_0 = pad_1d_array(gen_features.molecule.numbers[neighbors_idxs_0].astype(np.float), 1)\n    # Atom 1 may be anyone and may have up to 4 neighbors\n    neighbor_atoms_1 = pad_1d_array(gen_features.molecule.numbers[neighbors_idxs_1].astype(np.float), 4)\n    # Get distances to neighboring atoms\n    neighbor_distances_0 = pad_1d_array(\n        gen_features.molecule.distance_matrix[row[1], neighbors_idxs_0]\n        , 1\n    )\n    neighbor_distances_1 = pad_1d_array(\n        gen_features.molecule.distance_matrix[row[2], neighbors_idxs_1]\n        , 4\n    )\n    # Get normalized graph to first order neighborhood, may be used directly as factor or to match some patterns\n    neighborhood = gen_features.molecule.graph.get_subgraph(\n        list(set([*neighbors_idxs_0, *neighbors_idxs_1])) + [row[1], row[2]]\n        , True\n    ).blob\n    # Put all together in a single array\n    return np.hstack((\n        distance, neighbor_atoms_0, neighbor_atoms_1, neighbor_distances_0, neighbor_distances_1, neighborhood\n    ))\n\n# Variables for hashing\ngen_features.molecule = gen_features.molecule_name = None\n\n# Names for newly generated features\n# new_columns = ['distance'] + [\n#     f'{feature}_{i}_{j}' for feature in ['n_atom', 'n_distance'] for i in range(2) for j in range(4)\n# ] + ['neighborhood']\nnew_columns = [\n    'distance', 'n_atom_0_0', 'n_atom_1_0', 'n_atom_1_1', 'n_atom_1_2', 'n_atom_1_3'\n    , 'n_distance_0_0', 'n_distance_1_0', 'n_distance_1_1', 'n_distance_1_2', 'n_distance_1_3', 'neighborhood']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Add new features to train and test sets. We will not use `pandas.DataFrame.apply()` for our non-trivial `gen_features()` as people on the internet [say](https://ys-l.github.io/posts/2015/08/28/how-not-to-use-pandas-apply/) it is not memory and performance friendly. We split our data frame in chunks and generate features for each chunk separately. To speed things up we use some parallel processing."},{"metadata":{"trusted":true},"cell_type":"code","source":"def index_marks(nrows, chunk_size):\n    '''Get indices where to split df'''\n    return range(1 * chunk_size, (nrows // chunk_size + 1) * chunk_size, chunk_size)\n\ndef split_df(df, chunk_size):\n    '''Split df into chunks not larger than `chunk_size`'''\n    indices = index_marks(df.shape[0], chunk_size)\n    return np.split(df, indices)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"At first run on small part of data to check that everything is OK"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_head = train.head(13).copy()\ntrain_parts = split_df(train_head, 2)\nfor part in train_parts:\n    with multiprocessing.Pool(4) as pool:\n        part[new_columns] = pd.DataFrame( pool.map(gen_features, np.array(part[['molecule_name', 'atom_index_0', 'atom_index_1']])), index = part.index)\n    gc.collect()\ntrain_head = pd.concat(train_parts)\ndel train_parts\ngc.collect()\ntrain_head","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ntrain_parts = split_df(train, 100000)\nfor part in train_parts:\n    with multiprocessing.Pool(4) as pool:\n        part[new_columns] = pd.DataFrame( pool.map(gen_features, np.array(part[['molecule_name', 'atom_index_0', 'atom_index_1']])), index = part.index)\n    gc.collect()\ntrain = pd.concat(train_parts)\ndel train_parts\ngc.collect()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head(20)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ntest_parts = split_df(test, 100000)\nnew_vals = []\nfor part in test_parts:\n    with multiprocessing.Pool(4) as pool:\n        part[new_columns] = pd.DataFrame( pool.map(gen_features, np.array(part[['molecule_name', 'atom_index_0', 'atom_index_1']])), index = part.index)\n    gc.collect()\ntest = pd.concat(test_parts)\ndel test_parts\ngc.collect()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for col in [col for col in new_columns if not re.search('(neighborhood|type)', col)]:\n    train[col] = train[col].astype(np.float)\n    test[col] = test[col].astype(np.float)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Missing value imputation"},{"metadata":{"trusted":true},"cell_type":"code","source":"atom_cols = [c for c in train.columns if re.search('_atom', c)]\ndistance_cols = [c for c in train.columns if re.search('distance', c)]\n\ndef summary(series):\n    return {\n        'mean': np.mean(series)\n        , 'sd': np.std(series)\n        , 'min': np.amin(series)\n        , 'max': np.amax(series)\n        , 'NAs': np.sum(np.isnan(series)) / series.shape[0]\n    }\n\nfor col in atom_cols:\n    print(f'{col}: {summary(train[col])}')\n\nfor col in distance_cols:\n    print(f'{col}: {summary(train[col])}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Seems that `-10` should be fine for missing values.\n\nSet all atomic numbers to integers."},{"metadata":{"trusted":true},"cell_type":"code","source":"for col in atom_cols + distance_cols:\n    train[col] = train[col].fillna(-10)\n\nfor col in atom_cols:\n    train[col] = train[col].astype(np.int)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for f in [c for c in train.columns if re.search('(neighborhood|type)', c)]:\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))\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Select columns for model."},{"metadata":{"trusted":true},"cell_type":"code","source":"model_cols = [c for c in train.columns if not re.search('(^id$|molecule_name|_index|coupling_constant)', c)]\nmodel_cols","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Split training set for modelling."},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train, X_test, y_train, y_test = train_test_split(\n    train[model_cols]\n    , train['scalar_coupling_constant']\n    , test_size = 0.15, random_state = 0\n)\nprint(X_train.shape)\nprint(X_test.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Do some default LightGBM modeling."},{"metadata":{"trusted":true},"cell_type":"code","source":"params = {'boosting': 'gbdt', 'colsample_bytree': 1, \n          'learning_rate': 0.1, 'max_depth': 200, 'metric': 'mae',\n          'min_child_samples': 50, 'num_leaves': 500, \n          'objective': 'regression', 'reg_alpha': 0.5, \n          'reg_lambda': 0.8, 'subsample': 0.5,\n          'n_jobs': 4\n     }\n\nlgtrain = lgb.Dataset(X_train, label = y_train)\nlgval = lgb.Dataset(X_test, label = y_test)\n\nmodel_lgb = lgb.train(params, lgtrain, 1000, valid_sets = [lgtrain, lgval], early_stopping_rounds = 250, verbose_eval = 500)\n\ndef eval_fn(y_pred, y_test, c_type):\n    diff = np.full((y_pred.shape[0], 2), 1e-9)\n    diff[:,1] = np.abs(y_pred - y_test)\n    step_1 = pd.DataFrame(\n        {'diff': np.amax(diff, 1), 'type': c_type}\n    ).groupby(\n        'type'\n    ).mean()\n    return np.sum(np.log(step_1['diff'])) / step_1.shape[0]\n\ny_pred = model_lgb.predict(X_test)\nscore = eval_fn(y_pred, y_test, X_test['atom_1'])\nprint(f'Evaluation score: {score}')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Feature importance."},{"metadata":{"trusted":true},"cell_type":"code","source":"def plotImp(model, X , num = 20):\n    feature_imp = pd.DataFrame(sorted(zip(model.feature_importance(),X.columns)), columns = ['Value','Feature'])\n    plt.figure(figsize = (40, 20))\n    sns.set(font_scale = 5)\n    sns.barplot(x = \"Value\", y = \"Feature\", data=feature_imp.sort_values(by = \"Value\", ascending = False)[0:num])\n    plt.title('LightGBM Features (avg over folds)')\n    plt.tight_layout()\n    plt.xscale('log')\n    plt.show()\n\nplotImp(model_lgb, X_train, 20)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Predict and save."},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions = model_lgb.predict(test[model_cols])\nresult = pd.DataFrame({'id': test['id'], 'scalar_coupling_constant': predictions})\nprint(result.shape)\nresult.to_csv('prediction.csv', index = False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Conclusions\n\nHere we extracted just some very basic features with molmod library and achieved pretty good result. The library is certainly not limited to this, it provides vast opurtunities to generate geometric and molecular pattern features.\n\nThe feature extraction is quite slow as it is done for each row indivually. However it could be speeded up by parallel processing and/or by extracting only the required atom indices with molmod and then calculating distances and angles in vectorized way.\n\nWhile using it for a particular `type` group I am able to get close to **-2** scores and I'm not quite done yet."}],"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}