{"cells":[{"metadata":{},"cell_type":"markdown","source":"# How to get a good score with only 31 features\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.\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* For each of the 10-nearest neighbor atoms calculate the angle between this atom, the center-point and atom_1. \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\n* angle between atom, center-point and atom_1 (using the cosine, the feature is automatically scaled between 1 and -1)\n   \nThrow in the distance between atom_0 and atom_1 and we have our 31 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* Which side of the interactions: \n    - close to 1:  side of atom_1\n    - close to -1: side of atom_0\n    - close to 0:  equal distance to both (in-between)\n    \nEssentially, what we got is a 2C polar-coordinate system centered around the center-point between atom_0 and atom_1.  \nThe only thing missing for a 3D poloar-coordinate system representing the whole molecular geometry would be to add the angle of rotion around the atom_0-atom_1 axis.  \nSee next section, why this wasn't done :-)\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 two main drawbacks:\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\nNote that the 31 features above are not simply 31 scalar features. Obviously there is some structure:\nFor each atom, the features type, distance and cosine belong togheter (actually it's more like a single, vector valued feature rather than 3 features).  \nHowever, the model can't account for this structure when it's simply provided a list of features. For this reason, adding even more features to this approach would probably not lead to substantial improvements.\n\n## Runtime\n\nIt takes forever.  \nThe computation was not optimized due to the inherent limitations of the approach explained above :-)\n  \n"},{"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":{"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":"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    if verbose:\n        print(df.shape)\n        display(df.head())\n        \n    return df\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":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"def get_center(p0: np.array, p1: np.array) -> np.array:\n    return p0 + (p1 - p0)/2\n\n\ndef get_vector(p0: np.array, p1: np.array) -> np.array:\n    return p1 - p0\n\n\ndef calc_cosines(coords: np.array, a0: int, a1: int) -> np.array:\n    \"\"\"\n    @coordinates:\n    atom    x     y     z\n    1       1     2     3\n    2       4     5     6\n    ...\n    \"\"\"\n    atom_axes = coords - get_center(coords[a0, :], coords[a1, :])\n    main_axis = get_vector(coords[a0, :], coords[a1, :])\n    \n    dot_products = np.dot(atom_axes, main_axis)\n    atom_axes_norms = np.apply_along_axis(norm, 1, atom_axes)\n    main_axis_norm = norm(main_axis)\n    \n    return dot_products / (atom_axes_norms * main_axis_norm)\n\n\ndef get_knn_coordinates(j_coupling: pd.Series,\n                        structures=structures_df,\n                        mol2dist=mol2distance_matrix,\n                        k=10) -> np.array:\n    \n    a_0, a_1 = j_coupling.atom_index_0, j_coupling.atom_index_1\n    mol_df   = structures.loc[j_coupling.molecule_index]\n    \n    coordinates = mol_df[['x','y', 'z']].values\n    a0_coords   = coordinates[a_0, :]\n    a1_coords   = coordinates[a_1, :]\n    center      = get_center(a0_coords, a1_coords)\n    \n    cosines = calc_cosines(coordinates, a_0, a_1)\n    \n    center_distances = distance_matrix(center.reshape(1, 3), coordinates).ravel()\n    knn = np.argsort(center_distances)[:(k + 2)]  # atom-indices of KNN-atoms\n    \n    knn = np.array([x for x in knn if x not in (a_0, a_1)])\n    \n    distances = center_distances[knn]\n    cosines   = cosines[knn]\n    types = mol_df.iloc[knn].atom\n    \n    distances = np.pad(distances, (0, k - len(distances)), 'constant')\n    cosines   = np.pad(cosines,   (0, k - len(cosines)),   'constant')\n    types     = np.pad(types,     (0, k - len(types)),     'constant')\n    \n    d_a0_a1 = norm(a1_coords - a0_coords)\n\n    return np.concatenate([[d_a0_a1], distances, cosines, types])\n\n\n# this may take a while...\nid2features = {row.id : get_knn_coordinates(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, id2features)\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, id2features)","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)\nid2features_test = {row.id : get_knn_coordinates(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, id2features, 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, id2features_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}