{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Obtain a decent score with just 11 distances and 10 types\n\nInspired by the really cool kernel \"Distance - is all you need\" (https://www.kaggle.com/criskiev/distance-is-all-you-need-lb-1-481), I tried to devise my own system of simple, geomtric features. I started out with the simplest set of distance-features, I could think of.\n\n## Steps\n\n* Calculate the center-point on the axis between the two j-coupled atoms\n  (let's refer to them as atom_0 and and atom_1)\n* Get the 10 closet atoms and rank them by their distance to the center point (from closest to farthest away)\n\n\n## Features\nNow, or each of the 10-nearest neighbor atoms, we have 3 features:\n\n* atom-type\n* distance from the center-point\n   \nAdding the distance between atom_0 and atom_1 gives at total of 21 features.\n\n    \n## Intuition\n\nThese features provide three sources of information:\n\n* Which atom: H, C, N, O, or F\n* How far away from the j-coupled atoms\n\n\n## Model\n\nThe same LGBM used in the above mentioned kernel.  \nOnly more estimators are used for test-set prediction\n\n## Limitations\n\nThe decent score given the simple features suggests that geometry and atom type is all you need (it's also all you got ... :-D).\nHowever, this simple approach suffers from 3 main drawbacks:\n\n#### Lack of geometric information:\n\nDistance get's you surprisingly far, but has it's limitations. This features contain no information about whether a neighbor atom is on the side of atom_0, on the side of atom_1 or orthogonal to the atom_0-atom_1 axis. This is obviously a severe limation and it's surprising, the approch even reaches a relativels decent score.\n\n#### Permutation\n  A small difference in distances can lead to a different order of neighbor-atoms. This way very similar geometries can end up a different feature ordering which impedes learning.\n... that's why Message Passing Neural Networks or Graph Neural networks are used: the employ permutation invariant aggregation functions, such as $sum$ or $mean$.\n  \n#### No feature structure\n\nThe 21 features above are not simply 21 scalar features. Obviously there is some structure:\nFor each atom, the features type and distance belong togheter. However, the model can't account for this structure when it's simply provided a list of features.\n\n## Runtime\n\nIt takes forever.  \nThe computation was not optimized due to the inherent limitations of the approach explained above :-)"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"%%capture\n\nimport os\nfrom os.path import join\nimport sys\nfrom pprint import pprint\n\nimport numpy as np\nfrom numpy.linalg import norm\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom scipy.spatial import distance_matrix\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import mean_absolute_error\n\nfrom lightgbm import LGBMRegressor\n\nimport warnings\nwarnings.filterwarnings('ignore')\n\nDATA_DIR = '../input/champs-scalar-coupling'\nATOMIC_NUMBERS = {\n    'H': 1,\n    'C': 6,\n    'N': 7,\n    'O': 8,\n    'F': 9\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load Data"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"structures_dtypes = {\n    'molecule_name': 'category',\n    'atom_index': 'int8',\n    'atom': 'category',\n    'x': 'float32',\n    'y': 'float32',\n    'z': 'float32'\n}\nstructures_df = pd.read_csv(join(DATA_DIR, 'structures.csv'), dtype=structures_dtypes)\n\nstructures_df['molecule_index'] = structures_df.molecule_name.str.replace('dsgdb9nsd_', '').astype('int32')\nstructures_df['atom'] = structures_df['atom'].replace(ATOMIC_NUMBERS).astype('int8')\nstructures_df.index = structures_df.molecule_index\nstructures_df = structures_df[['atom_index', 'atom', 'x', 'y', 'z']]\n\nstructures_df[['x', 'y', 'z']] = structures_df[['x', 'y', 'z']] / 10  # puts all distances approx. in range [0, 1]\nprint(structures_df.shape)\ndisplay(structures_df.head())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\nThese two functions are from this cool kernel:\nhttps://www.kaggle.com/criskiev/distance-is-all-you-need-lb-1-481\n\"\"\"\n\ndef add_coordinates(base, structures, index):\n    df = pd.merge(base, structures, how='inner',\n                  left_on=['molecule_index', f'atom_index_{index}'],\n                  right_on=['molecule_index', 'atom_index']).drop(['atom_index'], axis=1)\n    df = df.rename(columns={\n        'atom': f'atom_{index}',\n        'x': f'x_{index}',\n        'y': f'y_{index}',\n        'z': f'z_{index}'\n    })\n    return df\n\n\ndef add_center(df):\n    df['x_c'] = ((df['x_1'] + df['x_0']) * np.float32(0.5))\n    df['y_c'] = ((df['y_1'] + df['y_0']) * np.float32(0.5))\n    df['z_c'] = ((df['z_1'] + df['z_0']) * np.float32(0.5))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_j_coupling_csv(file_path: str, train=True, verbose=False):\n    train_dtypes = {\n        'molecule_name': 'category',\n        'atom_index_0': 'int8',\n        'atom_index_1': 'int8',\n        'type': 'category',\n        'scalar_coupling_constant': 'float32'\n    }\n    df = pd.read_csv(file_path, dtype=train_dtypes)\n    df['molecule_index'] = df.molecule_name.str.replace('dsgdb9nsd_', '').astype('int32')\n    \n    if train:\n        cols = ['id', 'molecule_index', 'atom_index_0', 'atom_index_1', 'type', 'scalar_coupling_constant']\n    else: \n        cols = ['id', 'molecule_index', 'atom_index_0', 'atom_index_1', 'type']\n    df = df[cols]\n\n    df = add_coordinates(df, structures_df, 0)\n    df = add_coordinates(df, structures_df, 1)\n    add_center(df)\n    \n    if verbose:\n        print(df.shape)\n        display(df.head())\n        \n    return df\n\n\ntrain_df = load_j_coupling_csv(join(DATA_DIR, 'train.csv'), verbose=True)\n\nmol2distance_matrix = structures_df.groupby('molecule_index').apply(\n    lambda df: distance_matrix(df[['x','y', 'z']].values, df[['x','y', 'z']].values))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Calculate Features"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_knn_features_center(j_coupling: pd.Series,\n                            structures=structures_df,\n                            mol2dist=mol2distance_matrix,\n                            k=10) -> np.array:\n    \n    center      = j_coupling[['x_c', 'y_c', 'z_c']].values.reshape(1, 3)\n    mol_df      = structures.loc[j_coupling.molecule_index]\n    coordinates = mol_df[['x','y', 'z']].values\n    \n    center_distances = distance_matrix(center, coordinates).ravel()\n    knn = np.argsort(center_distances)[:(k + 2)]  # atom-indices of KNN-atoms\n    \n    a_0, a_1 = j_coupling.atom_index_0, j_coupling.atom_index_1\n    knn = np.array([x for x in knn if x not in (a_0, a_1)])\n    \n    distances = center_distances[knn]\n    types = mol_df.iloc[knn].atom\n    \n    distances = np.pad(distances, (0, k - len(distances)), 'constant')\n    types =     np.pad(types,     (0, k - len(types)),     'constant')\n    \n    dist_matrix = mol2distance_matrix[j_coupling.molecule_index]\n    d_a0_a1 = dist_matrix[a_0, a_1]\n\n    knn_features = np.concatenate([[d_a0_a1], distances, types])\n    return knn_features\n\n\n# this may take a while...\nid2center_knn = {row.id : get_knn_features_center(row) for _, row in train_df.iterrows()}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Make Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_data(df: pd.DataFrame, id2features: dict, random_state=128, split=True):\n    tmp_df = df.copy()\n    tmp_df['features'] = tmp_df.id.map(id2features)\n    \n    X = np.stack(tmp_df.features)\n    y = tmp_df.scalar_coupling_constant.values\n\n    if split:\n        X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=random_state)\n        return (X_train, y_train), (X_val, y_val)\n    else:\n        return X, y\n\n    \n# hyper-parameters like in this kernel:\n# https://www.kaggle.com/criskiev/distance-is-all-you-need-lb-1-481\nLGB_PARAMS = {\n        'objective': 'regression',\n        'metric': 'mae',\n        'verbosity': -1,\n        'boosting_type': 'gbdt',\n        'learning_rate': 0.2,\n        'num_leaves': 128,\n        'min_child_samples': 79,\n        'max_depth': 9,\n        'subsample_freq': 1,\n        'subsample': 0.9,\n        'bagging_seed': 11,\n        'reg_alpha': 0.1,\n        'reg_lambda': 0.3,\n        'colsample_bytree': 1.0\n        }\n\n    \ndef train_model(train, validation):\n    \n    X_train, y_train = train\n    X_val,   y_val   = validation\n\n    model = LGBMRegressor(**LGB_PARAMS, n_estimators=1500, n_jobs = -1)  # 6000 estimators would be better but take much longer\n    model.fit(X_train, y_train,\n            eval_set=[(X_train, y_train), (X_val, y_val)], eval_metric='mae',\n            verbose=100, early_stopping_rounds=200)\n\n    y_pred    = model.predict(X_val)\n    score     = np.log(mean_absolute_error(y_val, y_pred))\n    residuals = y_val - y_pred\n    \n    print(f'competition-metric score: {score}')\n    return model, score, residuals\n\n\ndef plot_residuals(residuals: np.array):\n    plt.hist(residuals, bins=50)\n    plt.title('residual distribution')\n    plt.show();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# test pipeline with smallest type:\n\nsub_train_df = train_df.query('type == \"1JHN\"')\n\n(X_train, y_train), (X_val, y_val) = make_data(sub_train_df, id2center_knn)\nprint(X_train.shape, X_val.shape, y_train.shape, y_val.shape)\n\nmodel, score, residuals = train_model((X_train, y_train), (X_val, y_val))\nplot_residuals(residuals)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Cross-validate"},{"metadata":{"trusted":true},"cell_type":"code","source":"def cross_validate(df, id2features):\n    type2model = {}\n    scores = {}\n\n    for type_, type_df in df.groupby('type'):\n        print(f'\\n\\n### {type_}')\n        train, validation = make_data(type_df, id2features)\n        model, score, residuals = train_model(train, validation)\n        type2model[type_] = model\n        scores[type_] = score\n        plot_residuals(residuals)\n\n    assert len(scores) == len(type2model) == 8\n    return  scores\n\n\nscores = cross_validate(train_df, id2center_knn)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'competition-metric: {np.mean(list(scores.values())):.2f}')\nprint('scores per type:')\npprint(scores, width=1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Predict and submit"},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_test_data(df: pd.DataFrame, id2features: dict, random_state=128):\n    tmp_df = df.copy()\n    tmp_df['features'] = tmp_df.id.map(id2features)\n    X = np.stack(tmp_df.features)\n    return X\n\n\ntest_df = load_j_coupling_csv(join(DATA_DIR, 'test.csv'), train=False, verbose=True)\nid2center_knn_test = {row.id : get_knn_features_center(row) for _, row in test_df.iterrows()}\nprediction_df = pd.DataFrame()\n\n\nfor type_ in sorted(train_df.type.unique()):\n    print(f'\\n### {type_}')\n    \n    train_type_df = train_df.query('type == @type_')\n    X_train, y_train = make_data(train_type_df, id2center_knn, split=False)\n    model = LGBMRegressor(**LGB_PARAMS, n_estimators=2000, n_jobs = -1)  # more estimators for test-set\n    model.fit(X_train, y_train, eval_metric='mae')\n    \n    test_type_df = test_df.query('type == @type_')\n    X_test = make_test_data(test_type_df, id2center_knn_test)\n    y_hat  = model.predict(X_test)\n    \n    type_pred_df  = pd.DataFrame({'id': test_type_df.id, 'scalar_coupling_constant': y_hat})\n    prediction_df = pd.concat([prediction_df, type_pred_df], ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"prediction_df = prediction_df.sort_values('id')\nprediction_df.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"assert len(prediction_df) == len(test_df)\nprint(prediction_df.shape)\ndisplay(prediction_df.head())","execution_count":null,"outputs":[]}],"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":4,"nbformat_minor":1}