{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom sklearn.preprocessing import OneHotEncoder\nfrom sklearn.model_selection import StratifiedKFold\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport lightgbm as lgb\nimport warnings\nfrom itertools import product\nimport re\n\n# import torch\n# import torch.nn.functional as F\n# from torch.autograd import grad\nimport tensorflow as tf\ntf.enable_eager_execution()\ntfe = tf.contrib.eager\n\nimport pdb\nimport gc\n\nsns.set_style('whitegrid')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ad6704059e30aae6115cf479777247ac3a9bf393"},"cell_type":"markdown","source":"Let us see if we can improve the basic model in several public kernels by replacing the training objective with the actual objective!\nUpdates: \n* removed some nonsense features (like mjd_max)\n* changed mjd_diff to det_mjd_diff (detected observations window length)\n* used tensorflow eager mode as an alternative to PyTorch. Hessian is also working now. Note that Hessian is only used as weights in LightGBM/XGBoost, so its scale does not matter. However, the hyperparamter 'min_child_weight'/'min_sum_hessian_in_leaf' can mess with the Hessian scale. It is easier to simply set it to 0. Using Hessian does not seem to impact prediction quality much.\n* learning rate scale is in line with what you can expect from vanilla LightGBM now.\n* There does not seem to be any difference between using clipping or simple log softmax in the objective."},{"metadata":{"trusted":true,"_uuid":"d2af2160ac84987f607f38b2c4f58679d9aca736"},"cell_type":"code","source":"# # this is a reimplementation of the above loss function using pytorch expressions.\n# # Alternatively this can be done in pure numpy (not important here)\n# # note that this function takes raw output instead of probabilities from the booster\n# # Also be aware of the index order in LightDBM when reshaping (see LightGBM docs 'fobj')\n# def wloss_metric(preds, train_data):\n#     y_t = torch.tensor(train_data.get_label(), requires_grad=False).type(torch.LongTensor)\n#     y_h = torch.zeros(\n#         y_t.shape[0], len(classes), requires_grad=False).scatter(1, y_t.reshape(-1, 1), 1)\n#     y_h /= y_h.sum(dim=0, keepdim=True)\n#     y_p = torch.tensor(preds, requires_grad=False).type(torch.FloatTensor)\n#     if len(y_p.shape) == 1:\n#         y_p = y_p.reshape(len(classes), -1).transpose(0, 1)\n#     ln_p = torch.log_softmax(y_p, dim=1)\n#     wll = torch.sum(y_h * ln_p, dim=0)\n#     loss = -torch.dot(weight_tensor, wll) / torch.sum(weight_tensor)\n#     return 'wloss', loss.numpy() * 1., False\n\n\n# # with autograd or pytorch you can pretty much come up with any loss function you want\n# # without worrying about implementing the gradients yourself\n# def wloss_objective(preds, train_data):\n#     y_t = torch.tensor(train_data.get_label(), requires_grad=False).type(torch.LongTensor)\n#     y_h = torch.zeros(\n#         y_t.shape[0], len(classes), requires_grad=False).scatter(1, y_t.reshape(-1, 1), 1)\n#     ys = y_h.sum(dim=0, keepdim=True)\n#     y_h /= ys\n#     y_p = torch.tensor(preds, requires_grad=True).type(torch.FloatTensor)\n#     y_r = y_p.reshape(len(classes), -1).transpose(0, 1)\n#     ln_p = torch.log_softmax(y_r, dim=1)\n#     wll = torch.sum(y_h * ln_p, dim=0)\n#     loss = -torch.dot(weight_tensor, wll)\n#     grads = grad(loss, y_p, create_graph=True)[0]\n#     grads *= float(len(classes)) / torch.sum(1 / ys)  # scale up grads\n#     hess = torch.ones(y_p.shape)  # haven't bothered with properly doing hessian yet\n#     return grads.detach().numpy(), \\\n#         hess.detach().numpy()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"classes = [6, 15, 16, 42, 52, 53, 62, 64, 65, 67, 88, 90, 92, 95]\nclass_weight = {6: 1, 15: 2, 16: 1, 42: 1, 52: 1, 53: 1, 62: 1,\n                64: 2, 65: 1, 67: 1, 88: 1, 90: 1, 92: 1, 95: 1}\nweight_tensor = tf.convert_to_tensor(list(class_weight.values()), dtype=tf.float32)\nclass_dict = {c: i for i, c in enumerate(classes)}\n\ndef label_to_code(labels):\n    return np.array([class_dict[c] for c in labels])\n\n# this is the simplified original loss function by Olivier. It works excellently as an\n# evaluation function, but we won't be able to use it in training\ndef multi_weighted_logloss(y_true, y_preds):\n    if len(np.unique(y_true)) > 14:\n        classes.append(99)\n        class_weight[99] = 2\n    enc = OneHotEncoder(sparse=False, categories='auto')\n    if len(y_true.shape) == 1:\n        y_true = np.expand_dims(y_true, 1)\n    y_ohe = enc.fit_transform(y_true)\n    y_p = np.clip(a=y_preds, a_min=1e-15, a_max=1 - 1e-15)\n    if y_p.shape[0] > y_true.shape[0]:\n        y_p = y_p.reshape(y_true.shape[0], len(classes), order='F')\n        if y_p.shape[0] != y_true.shape[0]:\n            raise ValueError(\n                'Dimension Mismatch for y_p {0} and y_true {1}!'.format(\n                    y_p.shape, y_true.shape))\n    y_p_log = np.log(y_p)\n    y_log_ones = np.sum(np.multiply(y_ohe, y_p_log), axis=0)\n    nb_pos = np.sum(y_ohe, axis=0).astype(float)\n    class_arr = np.array([class_weight[k] for k in sorted(class_weight.keys())])\n    y_w = y_log_ones * class_arr / nb_pos\n    loss = - np.sum(y_w) / np.sum(class_arr)\n    return loss\n\n# we may implement the above loss function using Tensorflow then use automatic \n# differentiation for the grad and hess in the objective function\ndef wloss_metric(preds, train_data):\n    y_t = tf.convert_to_tensor(train_data.get_label())\n    y_h = tf.one_hot(y_t, depth=14, dtype=tf.float32)\n    y_h /= tf.reduce_sum(y_h, axis=0, keepdims=True)\n    y_p = tf.convert_to_tensor(preds, dtype=tf.float32)\n    if len(y_p.shape) == 1:\n        y_p = tf.transpose(tf.reshape(y_p, (len(classes), -1)), perm=(1, 0))\n#     ln_p = tf.nn.log_softmax(y_p, axis=1)\n    ln_p = tf.log(tf.clip_by_value(tf.nn.softmax(y_p, axis=1), 1e-15, 1-1e-15))\n    wll = tf.reduce_sum(y_h * ln_p, axis=0)\n    loss = -tf.reduce_sum(weight_tensor * wll) / tf.reduce_sum(weight_tensor)\n    return 'wloss', loss.numpy(), False\n\ndef grad(f):\n    return lambda x: tfe.gradients_function(f)(x)[0]\n\ndef wloss_objective(preds, train_data):\n    y_t = tf.convert_to_tensor(train_data.get_label())\n    y_h = tf.one_hot(y_t, depth=14, dtype=tf.float32)\n    ys = tf.reduce_sum(y_h, axis=0, keepdims=True)\n    y_h /= ys\n    y_p = tf.convert_to_tensor(preds, dtype=tf.float32)\n    def loss(y_p):\n        if len(y_p.shape) == 1:\n            y_p = tf.transpose(tf.reshape(y_p, (len(classes), -1)), perm=(1, 0))\n        ln_p = tf.nn.log_softmax(y_p, axis=1)\n#         ln_p = tf.log(tf.clip_by_value(tf.nn.softmax(y_p, axis=1), 1e-15, 1-1e-15))\n        wll = tf.reduce_sum(y_h * ln_p, axis=0)\n        return -tf.reduce_sum(weight_tensor * wll) * len(train_data.get_label())\n    grads = grad(loss)(y_p)\n#     hess = grad(grad(loss))(y_p)\n#     hess /= tf.reduce_mean(hess)\n    hess = tf.ones(y_p.shape)\n    return grads.numpy(), hess.numpy()\n\ndef softmax(x, axis=1):\n    z = np.exp(x)\n    return z / np.sum(z, axis=axis, keepdims=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"78522fd6a1aea6671b1f98b89a1bfc0d1aff5fd1"},"cell_type":"code","source":"# we use some synthetic data to verify that the loss functions are correct:\nmock_y_true = np.array(classes + [6] * (100 - 14))\nmock_pred_score = np.zeros((100, 14))\nmock_pred_score[:, 0] = 10\nmock_pred_score[:, 1:] = 5\nmock_preds = softmax(mock_pred_score)\nmulti_weighted_logloss(mock_y_true, mock_preds)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"583e4832bd387d8dda4bd5f0a3f7eac9bea9e90f"},"cell_type":"code","source":"wloss_metric(np.reshape(mock_pred_score, (-1), order='F'),\n             lgb.Dataset(None, label_to_code(mock_y_true)))[1]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"591b53b8ca95fbb35e33a76daec4dd77de12231e"},"cell_type":"code","source":"train_meta = pd.read_csv('../input/training_set_metadata.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"835ab7afdff4773be5e9bd715ea49595fe6c3c92"},"cell_type":"code","source":"gc.enable()\n\ntrain = pd.read_csv('../input/training_set.csv')\ntrain['flux_ratio_sq'] = np.power(train['flux'] / train['flux_err'], 2.0)\ntrain['flux_by_flux_ratio_sq'] = train['flux'] * train['flux_ratio_sq']\n\naggs = {\n#     'mjd': ['min', 'max', 'size'],\n#     'passband': ['min', 'max', 'mean', 'median', 'std','skew'],\n    'flux': ['min', 'max', 'mean', 'median', 'std','skew'],\n    'flux_err': ['min', 'max', 'mean', 'median', 'std','skew'],\n    'detected': ['mean','std'],\n    'flux_ratio_sq':['sum','skew'],\n    'flux_by_flux_ratio_sq':['sum','skew'],\n}\n\ntrain_feats = train.groupby('object_id').agg(aggs)\nnew_columns = [\n    k + '_' + agg for k in aggs.keys() for agg in aggs[k]\n]\ntrain_feats.columns = new_columns\ndetected_groups = train[train['detected'] == 1].groupby('object_id')\ntrain_feats['det_mjd_diff'] = detected_groups['mjd'].transform('max') \\\n    - detected_groups['mjd'].transform('min')\ntrain_feats['flux_diff'] = train_feats['flux_max'] - train_feats['flux_min']\ntrain_feats['flux_dif2'] = (train_feats['flux_max'] - \n                            train_feats['flux_min']) / train_feats['flux_mean']\ntrain_feats['flux_w_mean'] = train_feats['flux_by_flux_ratio_sq_sum'] \\\n    / train_feats['flux_ratio_sq_sum']\ntrain_feats['flux_dif3'] = (train_feats['flux_max'] - \n                            train_feats['flux_min']) / train_feats['flux_w_mean']\n\n# del train_feats['mjd_max'], train_feats['mjd_min']\n\ndel train\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9ba5257bd616772e9bb7c61fa49a6c18334da637"},"cell_type":"code","source":"full_features = train_feats.reset_index().merge(\n    right=train_meta,\n    how='outer',\n    on='object_id'\n)\n\nif 'target' in full_features:\n    target = full_features['target']\n    del full_features['target']\nclasses = sorted(target.unique())\n\n# Taken from Giba's topic : https://www.kaggle.com/titericz\n# https://www.kaggle.com/c/PLAsTiCC-2018/discussion/67194\n# with Kyle Boone's post https://www.kaggle.com/kyleboone\nclass_weight = {\n    c: 1 for c in classes\n}\nfor c in [64, 15]:\n    class_weight[c] = 2\n\nprint('Unique classes : ', classes)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9a9ebe0ae8eab4158d7528da367e3203b69d0de9"},"cell_type":"code","source":"if 'object_id' in full_features:\n    oof_df = full_features[['object_id']]\n    del full_features['object_id'], full_features['distmod'], full_features['hostgal_specz']\n    \n    \ntrain_mean = full_features.mean(axis=0)\nfull_features.fillna(train_mean, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1023d9d4a8a4872bedd9af33bf7048ab5ad22f82"},"cell_type":"code","source":"full_features.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6923f01a83c4bc3769b4d0996656a59ad91edd9f"},"cell_type":"code","source":"folds = StratifiedKFold(n_splits=5, shuffle=True, random_state=1111)\nclf_params = {\n    'boosting_type': 'gbdt',\n    'objective': 'multiclass',\n    'num_class': 14,\n    'metric': 'None',\n    'learning_rate': 0.01,\n    'subsample': .9,\n    'colsample_bytree': .75,\n    'reg_alpha': 1.0e-2,\n    'reg_lambda': 1.0e-2,\n    'min_split_gain': 0.01,\n#     'min_child_weight': 10,\n    'min_child_samples': 20,\n#     'n_estimators': 2000,\n#     'silent': -1,\n#     'verbose': -1,\n    'max_depth': 3,\n    'importance_type': 'gain',\n    'n_jobs': -1\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"befea89e6bf5f58c7933a31130aebdf8edfc7e13","scrolled":false},"cell_type":"code","source":"boosters = []\nimportances = pd.DataFrame()\noof_preds = np.zeros((full_features.shape[0], target.unique().shape[0]))\n\nwarnings.simplefilter('ignore', FutureWarning)\nfor fold_id, (train_idx, validation_idx) in enumerate(folds.split(full_features, target)):\n    print('processing fold {0}'.format(fold_id))\n    X_train, y_train = full_features.iloc[train_idx], target.iloc[train_idx]\n    X_valid, y_valid = full_features.iloc[validation_idx], target.iloc[validation_idx]\n\n    train_dataset = lgb.Dataset(X_train, label_to_code(y_train))\n    valid_dataset = lgb.Dataset(X_valid, label_to_code(y_valid))\n    \n    booster = lgb.train(clf_params.copy(), train_dataset, \n                        num_boost_round=2000,\n                        fobj=wloss_objective, \n                        feval=wloss_metric,\n                        valid_sets=[train_dataset, valid_dataset],\n                        verbose_eval=100,\n                        early_stopping_rounds=100\n                       )\n    oof_preds[validation_idx, :] = booster.predict(X_valid)\n\n    imp_df = pd.DataFrame()\n    imp_df['feature'] = full_features.columns\n    imp_df['gain'] = booster.feature_importance('gain')\n    imp_df['fold'] = fold_id\n    importances = pd.concat([importances, imp_df], axis=0, sort=False)\n\n    boosters.append(booster)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f98b80a226b178b12e8b0dd13e3a0fc69654b152"},"cell_type":"code","source":"loss = multi_weighted_logloss(y_true=target, y_preds=softmax(oof_preds))\n_, loss2, _ = wloss_metric(oof_preds, lgb.Dataset(full_features, label_to_code(target)))\nprint(f'OG wloss : {loss:.5f}, Re-implemented wloss: {loss2:.5f} ')  # 0.93526","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0f23ca923962fbb2913b015f4826194290446ebc"},"cell_type":"code","source":"mean_gain = importances[['gain', 'feature']].groupby('feature').mean()\nimportances['mean_gain'] = importances['feature'].map(mean_gain['gain'])\nplt.figure(figsize=(14, 14), facecolor='w')\nsns.barplot(x='gain', y='feature', \n            data=importances.sort_values('mean_gain', ascending=False).iloc[:5 * 40])\nplt.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5d0423683e3f947816aa7f278ca8d44d23169c87"},"cell_type":"code","source":"meta_test = pd.read_csv('../input/test_set_metadata.csv')\n\nimport time\n\nstart = time.time()\nchunks = 5000000\nfor i_c, df in enumerate(pd.read_csv('../input/test_set.csv', chunksize=chunks, iterator=True)):\n    df['flux_ratio_sq'] = np.power(df['flux'] / df['flux_err'], 2.0)\n    df['flux_by_flux_ratio_sq'] = df['flux'] * df['flux_ratio_sq']\n    # Group by object id\n    agg_test = df.groupby('object_id').agg(aggs)\n    agg_test.columns = new_columns\n#     agg_test['mjd_diff'] = agg_test['mjd_max'] - agg_test['mjd_min']\n    detected_groups = df[df['detected'] == 1].groupby('object_id')\n    agg_test['det_mjd_diff'] = detected_groups['mjd'].transform('max') \\\n        - detected_groups['mjd'].transform('min')\n    agg_test['flux_diff'] = agg_test['flux_max'] - agg_test['flux_min']\n    agg_test['flux_dif2'] = (agg_test['flux_max'] - agg_test['flux_min']) / agg_test['flux_mean']\n    agg_test['flux_w_mean'] = agg_test['flux_by_flux_ratio_sq_sum'] / agg_test['flux_ratio_sq_sum']\n    agg_test['flux_dif3'] = (agg_test['flux_max'] - agg_test['flux_min']) / agg_test['flux_w_mean']\n\n#     del agg_test['mjd_max'], agg_test['mjd_min']\n    \n    full_test = agg_test.reset_index().merge(\n        right=meta_test,\n        how='left',\n        on='object_id'\n    )\n    full_test = full_test.fillna(train_mean)\n    \n    # Make predictions\n    preds = None\n    for booster in boosters:\n        if preds is None:\n            preds = softmax(booster.predict(full_test[full_features.columns])) / folds.n_splits\n        else:\n            preds += softmax(booster.predict(full_test[full_features.columns])) / folds.n_splits\n    \n   # Compute preds_99 as the proba of class not being any of the others\n    # preds_99 = 0.1 gives 1.769\n    preds_99 = np.ones(preds.shape[0])\n    for i in range(preds.shape[1]):\n        preds_99 *= (1 - preds[:, i])\n    \n    # Store predictions\n    preds_df = pd.DataFrame(preds, columns=['class_' + str(s) for s in classes])\n    preds_df['object_id'] = full_test['object_id']\n    preds_df['class_99'] = 0.14 * preds_99 / np.mean(preds_99) \n    \n    if i_c == 0:\n        preds_df.to_csv('predictions.csv',  header=True, mode='a', index=False, float_format='%.6f')\n    else: \n        preds_df.to_csv('predictions.csv',  header=False, mode='a', index=False, float_format='%.6f')\n        \n    del agg_test, full_test, preds_df, preds\n    gc.collect()\n    \n    if (i_c + 1) % 10 == 0:\n        print('%15d done in %5.1f' % (chunks * (i_c + 1), (time.time() - start) / 60))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"614f3a57dd1177580e906ca9424edc839dba1f1c"},"cell_type":"code","source":"z = pd.read_csv('predictions.csv')\n\nprint(z.groupby('object_id').size().max())\nprint((z.groupby('object_id').size() > 1).sum())\n\nz = z.groupby('object_id').mean()\n\nz.to_csv('single_predictions.csv', index=True)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}