{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# Import Libraries\nimport pandas as pd\nimport numpy as np\nfrom numba import jit\nimport tsfresh\nfrom tsfresh.feature_extraction import extract_features\nimport sys\nimport gc; gc.enable()\nimport time\n\nfrom datetime import datetime\nfrom functools import partial\nfrom lightgbm import LGBMClassifier\nfrom sklearn.model_selection import StratifiedKFold\n\nimport os\nprint(os.listdir(\"../input\"))","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"@jit\ndef haversine_plus(lon1, lat1, lon2, lat2):\n    \"\"\"\n    Calculate the great circle distance between two points on the earth (specified in decimal degrees) from\n    #https://stackoverflow.com/questions/4913349/haversine-formula-in-python-bearing-and-distance-between-two-gps-points\n    \"\"\"\n    # Convert decimal degrees to Radians:\n    lon1 = np.radians(lon1)\n    lat1 = np.radians(lat1)\n    lon2 = np.radians(lon2)\n    lat2 = np.radians(lat2)\n\n    # Haversine Formula:\n    dlon = np.subtract(lon2, lon1)\n    dlat = np.subtract(lat2, lat1)\n\n    a = np.add(np.power(np.sin(np.divide(dlat, 2)), 2),\n               np.multiply(np.cos(lat1),\n                           np.multiply(np.cos(lat2),\n                                       np.power(np.sin(np.divide(dlon, 2)), 2))))\n\n    haversine = np.multiply(2, np.arcsin(np.sqrt(a)))\n    return {\n        'haversine': haversine,\n        'latlon1': np.subtract(np.multiply(lon1, lat1), np.multiply(lon2, lat2)),\n    }\n\n\n@jit\ndef process_flux(df):\n    flux_ratio_sq = np.power(df['flux'].values / df['flux_err'].values, 2.0)\n\n    df_flux = pd.DataFrame({\n        'flux_ratio_sq': flux_ratio_sq,\n        'flux_by_flux_ratio_sq': df['flux'].values * flux_ratio_sq\n        },\n        index=df.index\n    )\n\n    return pd.concat([df, df_flux], axis=1)\n\n\n@jit\ndef process_flux_agg(df):\n    flux_w_mean = df['flux_by_flux_ratio_sq_sum'].values / df['flux_ratio_sq_sum'].values\n    flux_diff = df['flux_max'].values - df['flux_min'].values\n\n    df_flux_agg = pd.DataFrame({\n        'flux_w_mean': flux_w_mean,\n        'flux_diff1': flux_diff,\n        'flux_diff2': flux_diff / df['flux_mean'].values,\n        'flux_diff3': flux_diff / flux_w_mean,\n        },\n        index=df.index\n    )\n\n    return pd.concat([df, df_flux_agg], axis=1)\n\n\nAGGS = {\n    'flux': ['min', 'max', 'mean', 'median', 'std', 'skew'],\n    'flux_err': ['min', 'max', 'mean', 'median', 'std', 'skew'],\n    'detected': ['min', 'max', 'mean', 'median', 'std', 'skew'],\n    'flux_ratio_sq': ['sum', 'skew'],\n    'flux_by_flux_ratio_sq': ['sum', 'skew'],\n}\n\n# tsfresh features\nFCP = {\n    'flux': {\n        'longest_strike_above_mean': None,\n        'longest_strike_below_mean': None,\n        'mean_change': None,\n        'mean_abs_change': None,\n        'length': None,\n    },\n    'flux_by_flux_ratio_sq': {\n        'longest_strike_above_mean': None,\n        'longest_strike_below_mean': None,\n    },\n    'flux_passband': {\n        'fft_coefficient': [\n            {'coeff': 0, 'attr': 'abs'},\n            {'coeff': 1, 'attr': 'abs'}\n        ],\n        'kurtosis': None,\n        'skewness': None,\n    },\n    'mjd': {\n        'maximum': None,\n        'minimum': None,\n        'mean_change': None,\n        'mean_abs_change': None,\n    },\n}\n\n\ndef featurize(df, df_meta, n_jobs=4):\n    \"\"\"\n    Extracting Features from train set\n    Features from olivier's kernel https://www.kaggle.com/ogrellier/plasticc-in-a-kernel-meta-and-data\n    Very smart and powerful feature that is generously given here https://www.kaggle.com/c/PLAsTiCC-2018/discussion/69696#410538\n    Per passband features with tsfresh library.\n    fft features added to capture periodicity https://www.kaggle.com/c/PLAsTiCC-2018/discussion/70346#415506\n    \"\"\"\n\n    df = process_flux(df)\n\n    # agg features\n    aggs = AGGS\n    # tsfresh features\n    fcp = FCP\n\n    agg_df = df.groupby('object_id').agg(aggs)\n    agg_df.columns = ['{}_{}'.format(k, agg) for k in aggs.keys() for agg in aggs[k]]\n    agg_df = process_flux_agg(agg_df)  # new feature to play with tsfresh\n\n    # Add more features with\n    agg_df_ts_flux_passband = extract_features(df,\n                                               column_id='object_id',\n                                               column_sort='mjd',\n                                               column_kind='passband',\n                                               column_value='flux',\n                                               default_fc_parameters=fcp['flux_passband'],\n                                               n_jobs=n_jobs\n                                               )\n\n    agg_df_ts_flux = extract_features(df,\n                                      column_id='object_id',\n                                      column_value='flux',\n                                      default_fc_parameters=fcp['flux'], n_jobs=n_jobs)\n\n    agg_df_ts_flux_by_flux_ratio_sq = extract_features(df,\n                                                       column_id='object_id',\n                                                       column_value='flux_by_flux_ratio_sq',\n                                                       default_fc_parameters=fcp['flux_by_flux_ratio_sq'],\n                                                       n_jobs=n_jobs)\n\n    # Add smart feature that is suggested here https://www.kaggle.com/c/PLAsTiCC-2018/discussion/69696#410538\n    # dt[detected==1, mjd_diff:=max(mjd)-min(mjd), by=object_id]\n    df_det = df[df['detected'] == 1].copy()\n    agg_df_mjd = extract_features(df_det,\n                                  column_id='object_id',\n                                  column_value='mjd',\n                                  default_fc_parameters=fcp['mjd'], n_jobs=n_jobs)\n    agg_df_mjd['mjd_diff_det'] = agg_df_mjd['mjd__maximum'].values - agg_df_mjd['mjd__minimum'].values\n    del agg_df_mjd['mjd__maximum'], agg_df_mjd['mjd__minimum']\n\n    agg_df_ts_flux_passband.index.rename('object_id', inplace=True)\n    agg_df_ts_flux.index.rename('object_id', inplace=True)\n    agg_df_ts_flux_by_flux_ratio_sq.index.rename('object_id', inplace=True)\n    agg_df_mjd.index.rename('object_id', inplace=True)\n    agg_df_ts = pd.concat([agg_df,\n                           agg_df_ts_flux_passband,\n                           agg_df_ts_flux,\n                           agg_df_ts_flux_by_flux_ratio_sq,\n                           agg_df_mjd], axis=1).reset_index()\n\n    result = agg_df_ts.merge(right=df_meta, how='left', on='object_id')\n    return result\n\n\ndef process_meta(filename):\n    meta_df = pd.read_csv(filename)\n\n    meta_dict = dict()\n    # distance\n    meta_dict.update(haversine_plus(meta_df['ra'].values, meta_df['decl'].values, meta_df['gal_l'].values, meta_df['gal_b'].values))\n\n    meta_dict['hostgal_photoz_certain'] = np.multiply(\n        meta_df['hostgal_photoz'].values,\n        np.exp(meta_df['hostgal_photoz_err'].values))\n\n    meta_df = pd.concat([meta_df, pd.DataFrame(meta_dict, index=meta_df.index)], axis=1)\n    return meta_df\n\ndef multi_weighted_logloss(y_true, y_preds, classes, class_weights):\n    \"\"\"\n    refactor from\n    @author olivier https://www.kaggle.com/ogrellier\n    multi logloss for PLAsTiCC challenge\n    \"\"\"\n    y_p = y_preds.reshape(y_true.shape[0], len(classes), order='F')\n    # Trasform y_true in dummies\n    y_ohe = pd.get_dummies(y_true)\n    # Normalize rows and limit y_preds to 1e-15, 1-1e-15\n    y_p = np.clip(a=y_p, a_min=1e-15, a_max=1 - 1e-15)\n    # Transform to log\n    y_p_log = np.log(y_p)\n    # Get the log for ones, .values is used to drop the index of DataFrames\n    # Exclude class 99 for now, since there is no class99 in the training set\n    # we gave a special process for that class\n    y_log_ones = np.sum(y_ohe.values * y_p_log, axis=0)\n    # Get the number of positives for each class\n    nb_pos = y_ohe.sum(axis=0).values.astype(float)\n    # Weight average and divide by the number of positives\n    class_arr = np.array([class_weights[k] for k in sorted(class_weights.keys())])\n    y_w = y_log_ones * class_arr / nb_pos\n\n    loss = - np.sum(y_w) / np.sum(class_arr)\n    return loss\n\n\ndef lgbm_multi_weighted_logloss(y_true, y_preds):\n    \"\"\"\n    refactor from\n    @author olivier https://www.kaggle.com/ogrellier\n    multi logloss for PLAsTiCC challenge\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\n    galactic_classes = [6, 16, 53, 65, 92]\n    extragalactic_classes = [15, 42, 52, 62, 64, 67, 88, 90, 95]\n    galactic_classes_weights = {c: 1 for c in galactic_classes}\n    extragalactic_classes_weights = {c: 1 for c in extragalactic_classes}\n    extragalactic_classes_weights.update({c: 2 for c in [64, 15]})\n    if len(y_preds) == 5 * y_true.shape[0]:\n        classes = galactic_classes\n        class_weights = galactic_classes_weights\n    else:\n        classes = extragalactic_classes\n        class_weights = extragalactic_classes_weights\n\n    loss = multi_weighted_logloss(y_true, y_preds, classes, class_weights)\n    return 'wloss', loss, False\n\n\ndef xgb_multi_weighted_logloss(y_predicted, y_true, classes, class_weights):\n    loss = multi_weighted_logloss(y_true.get_label(), y_predicted, classes, class_weights)\n    return 'wloss', loss","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"462dec2438e13323fe65da1ae7dc77c43d0febc0"},"cell_type":"code","source":"galactic_classes = [6, 16, 53, 65, 92]\nextragalactic_classes = [15, 42, 52, 62, 64, 67, 88, 90, 95]\n\ndef predict_chunk(df_, clfs_gal_, clfs_ext_, meta_, features):\n    # processing all features\n    full_test = featurize(df_, meta_)\n    full_test.fillna(0, inplace=True)\n\n    galactic_cut = full_test['hostgal_photoz'] == 0\n    gal_test = full_test[galactic_cut]\n    ext_test = full_test[~galactic_cut]\n\n    # predictions\n    preds_gal = None\n    if not gal_test.empty:\n        for clf in clfs_gal_:\n            if preds_gal is None:\n                preds_gal = clf.predict_proba(gal_test[features])\n            else:\n                preds_gal += clf.predict_proba(gal_test[features])\n\n        preds_gal = preds_gal / len(clfs_gal_)\n\n        preds_99_gal = np.ones(preds_gal.shape[0])\n        for i in range(preds_gal.shape[1]):\n            preds_99_gal *= (1 - preds_gal[:, i])\n\n        # Create DataFrame from predictions\n        preds_gal = pd.DataFrame(preds_gal,\n                                 columns=['class_{}'.format(s) for s in clfs_gal_[0].classes_])\n        assert preds_gal.shape[0] == gal_test.shape[0], 'len of preds={}, test={}'.format(preds_gal.shape[0], gal_test.shape[0])\n        preds_gal['object_id'] = gal_test['object_id'].values\n        for c in ['class_{}'.format(s) for s in extragalactic_classes]:\n            preds_gal.insert(0, c, 0.0)\n        preds_gal['class_99'] = 0.017 * preds_99_gal / np.mean(preds_99_gal)\n\n    preds_ext = None\n    if not ext_test.empty:\n        for clf in clfs_ext_:\n            if preds_ext is None:\n                preds_ext = clf.predict_proba(ext_test[features])\n            else:\n                preds_ext += clf.predict_proba(ext_test[features])\n\n        preds_ext = preds_ext / len(clfs_ext_)\n\n        preds_99_ext = np.ones(preds_ext.shape[0])\n        for i in range(preds_ext.shape[1]):\n            preds_99_ext *= (1 - preds_ext[:, i])\n\n        # Create DataFrame from predictions\n        preds_ext = pd.DataFrame(preds_ext,\n                                 columns=['class_{}'.format(s) for s in clfs_ext_[0].classes_])\n        assert preds_ext.shape[0] == ext_test.shape[0], 'len of preds={}, test={}'.format(preds_ext.shape[0], ext_test.shape[0])\n        preds_ext['object_id'] = ext_test['object_id'].values\n        for c in ['class_{}'.format(s) for s in galactic_classes]:\n            preds_ext.insert(0, c, 0.0)\n        preds_ext['class_99'] = 0.17 * preds_99_ext / np.mean(preds_99_ext)\n\n    preds_df_ = pd.concat([preds_gal, preds_ext], ignore_index=True, sort=False)\n\n    return preds_df_\n\n\ndef process_test(clfs_gal, clfs_ext,\n                 features,\n                 filename='predictions.csv',\n                 chunks=5000000):\n    start = time.time()\n\n    meta_test = process_meta('../input/test_set_metadata.csv')\n\n    remain_df = None\n    for i_c, df in enumerate(pd.read_csv('../input/test_set.csv', chunksize=chunks, iterator=True)):\n        # Check object_ids\n        # I believe np.unique keeps the order of group_ids as they appear in the file\n        unique_ids = np.unique(df['object_id'])\n\n        new_remain_df = df.loc[df['object_id'] == unique_ids[-1]].copy()\n        if remain_df is None:\n            df = df.loc[df['object_id'].isin(unique_ids[:-1])]\n        else:\n            df = pd.concat([remain_df, df.loc[df['object_id'].isin(unique_ids[:-1])]], axis=0)\n        # Create remaining samples df\n        remain_df = new_remain_df\n\n        preds_df = predict_chunk(df_=df,\n                                 clfs_gal_=clfs_gal,\n                                 clfs_ext_=clfs_ext,\n                                 meta_=meta_test,\n                                 features=features)\n\n        if i_c == 0:\n            preds_df.to_csv(filename, header=True, mode='a', index=False)\n        else:\n            preds_df.to_csv(filename, header=False, mode='a', index=False)\n\n        del preds_df\n        gc.collect()\n        print('{:15d} done in {:5.1f} minutes'.format(\n            chunks * (i_c + 1), (time.time() - start) / 60), flush=True)\n\n    # Compute last object in remain_df\n    preds_df = predict_chunk(df_=remain_df,\n                             clfs_gal_=clfs_gal,\n                             clfs_ext_=clfs_ext,\n                             meta_=meta_test,\n                             features=features)\n\n    preds_df.to_csv(filename, header=False, mode='a', index=False)\n    return\n\n\n\ndef lgbm_modeling_cross_validation(params,\n                                   full_train,\n                                   y,\n                                   classes,\n                                   class_weights,\n                                   id,\n                                   part,\n                                   nr_fold=5,\n                                   random_state=99,\n                                   ):\n    # Compute weights\n    w = y.value_counts()\n    weights = {i: np.sum(w) / w[i] for i in w.index}\n\n    clfs = []\n    importances = pd.DataFrame()\n    folds = StratifiedKFold(n_splits=nr_fold,\n                            shuffle=True,\n                            random_state=random_state)\n\n    oof_preds = np.zeros((len(full_train), np.unique(y).shape[0]))\n\n    for fold_, (trn_, val_) in enumerate(folds.split(y, y)):\n        trn_x, trn_y = full_train.iloc[trn_], y.iloc[trn_]\n        val_x, val_y = full_train.iloc[val_], y.iloc[val_]\n\n        clf = LGBMClassifier(**params)\n        clf.fit(\n            trn_x, trn_y,\n            eval_set=[(trn_x, trn_y), (val_x, val_y)],\n            eval_metric=lgbm_multi_weighted_logloss,\n            verbose=100,\n            early_stopping_rounds=50,\n            sample_weight=trn_y.map(weights)\n        )\n        clfs.append(clf)\n\n        oof_preds[val_, :] = clf.predict_proba(val_x, num_iteration=clf.best_iteration_)\n        print('no {}-fold loss: {}'.format(fold_ + 1,\n                                           multi_weighted_logloss(val_y, oof_preds[val_, :],\n                                                                  classes, class_weights)))\n\n        imp_df = pd.DataFrame({\n            'feature': full_train.columns,\n            'gain': clf.feature_importances_,\n            'fold': [fold_ + 1] * len(full_train.columns),\n        })\n        importances = pd.concat([importances, imp_df], axis=0, sort=False)\n\n    score = multi_weighted_logloss(y_true=y, y_preds=oof_preds,\n                                   classes=classes, class_weights=class_weights)\n    print('MULTI WEIGHTED LOG LOSS: {:.5f}'.format(score))\n\n    oof_preds_pd = pd.DataFrame(data=oof_preds, columns=['class_{}'.format(s) for s in classes])\n    pd.concat([id, oof_preds_pd], axis=1).to_csv('lgbm_train_oof_preds_{}.csv'.format(part), index=False)\n    df_importances = save_importances(importances_=importances)\n    df_importances.to_csv('lgbm_importances_{}.csv'.format(part), index=False)\n\n    return clfs, score\n\ndef save_importances(importances_):\n    mean_gain = importances_[['gain', 'feature']].groupby('feature').mean()\n    importances_['mean_gain'] = importances_['feature'].map(mean_gain['gain'])\n    return importances_","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"dbcd9a90342bc48a82062c63d94d25041ca7ae2d"},"cell_type":"code","source":"xgb_params = {\n    'objective': 'multi:softprob',\n    'eval_metric': 'mlogloss',\n    'silent': True,\n    'num_class': 14,\n\n    'booster': 'gbtree',\n    'n_jobs': 4,\n    'n_estimators': 1000,\n    'tree_method': 'hist',\n    'grow_policy': 'lossguide',\n    'base_score': 0.25,\n    'max_depth': 7,\n    'max_delta_step': 2,  # default=0\n    'learning_rate': 0.03,\n    'max_leaves': 11,\n    'min_child_weight': 64,\n    'gamma': 0.1,  # default=\n    'subsample': 0.7,\n    'colsample_bytree': 0.68,\n    'reg_alpha': 0.01,  # default=0\n    'reg_lambda': 10.,  # default=1\n    'seed': 1234\n\n}\n\nlgbm_params = {\n    'device': 'cpu',\n    'objective': 'multiclass',\n    'num_class': 14,\n    'boosting_type': 'gbdt',\n    'n_jobs': -1,\n    'max_depth': 7,\n    'n_estimators': 500,\n    'subsample_freq': 2,\n    'subsample_for_bin': 5000,\n    'min_data_per_group': 100,\n    'max_cat_to_onehot': 4,\n    'cat_l2': 1.0,\n    'cat_smooth': 59.5,\n    'max_cat_threshold': 32,\n    'metric_freq': 10,\n    'verbosity': -1,\n    'metric': 'multi_logloss',\n    'xgboost_dart_mode': False,\n    'uniform_drop': False,\n    'colsample_bytree': 0.5,\n    'drop_rate': 0.173,\n    'learning_rate': 0.0267,\n    'max_drop': 5,\n    'min_child_samples': 10,\n    'min_child_weight': 100.0,\n    'min_split_gain': 0.1,\n    'num_leaves': 7,\n    'reg_alpha': 0.1,\n    'reg_lambda': 0.00023,\n    'skip_drop': 0.44,\n    'subsample': 0.75\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1a97586f152916750f7037ce9b4d102d7ff3e0c1"},"cell_type":"code","source":"def train_model(full_train, classes, class_weights, part):\n    if 'target' in full_train:\n        y = full_train['target']\n        del full_train['target']\n\n    if 'object_id' in full_train:\n        object_id = full_train['object_id']\n        del full_train['object_id']\n        del full_train['hostgal_specz']\n        del full_train['ra'], full_train['decl'], full_train['gal_l'], full_train['gal_b']\n        del full_train['ddf']\n\n    eval_func = partial(lgbm_modeling_cross_validation,\n                        full_train=full_train,\n                        y=y,\n                        classes=classes,\n                        class_weights=class_weights,\n                        id=object_id,\n                        part=part,\n                        nr_fold=5,\n                        random_state=99,\n                        )\n\n    lgbm_params.update({'n_estimators': 2000})\n\n    # modeling from CV\n    clfs, score = eval_func(lgbm_params)\n\n    return clfs, score","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d9025684d4356c7fa60a9a20d3f8918532c8ccc3"},"cell_type":"code","source":"def main(argc, argv):\n    meta_train = process_meta('../input/training_set_metadata.csv')\n\n    train = pd.read_csv('../input/training_set.csv')\n    full_train = featurize(train, meta_train)\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\n    galactic_classes_weights = {c: 1 for c in galactic_classes}\n    extragalactic_classes_weights = {c: 1 for c in extragalactic_classes}\n    extragalactic_classes_weights.update({c: 2 for c in [64, 15]})\n\n    full_train.fillna(0, inplace=True)\n    galactic_cut = full_train['hostgal_photoz'] == 0\n\n    clfs_gal, score_gal = train_model(full_train[galactic_cut],\n                            galactic_classes, galactic_classes_weights, 'gal')\n    clfs_ext, score_ext = train_model(full_train[~galactic_cut],\n                            extragalactic_classes, extragalactic_classes_weights, 'ext')\n\n    filename = 'subm_{:.6f}_{:.6f}_{}.csv'.format(\n            score_gal, score_ext,\n            datetime.now().strftime('%Y-%m-%d-%H-%M'))\n    print('save to {}'.format(filename))\n\n    # TEST\n    if 'target' in full_train:\n        y = full_train['target']\n        del full_train['target']\n\n    if 'object_id' in full_train:\n        object_id = full_train['object_id']\n        del full_train['object_id']\n        del full_train['hostgal_specz']\n        del full_train['ra'], full_train['decl'], full_train['gal_l'], full_train['gal_b']\n        del full_train['ddf']\n    process_test(clfs_gal, clfs_ext,\n                 features=full_train.columns,\n                 filename=filename,\n                 chunks=5000000)\n\n    z = pd.read_csv(filename)\n    print(\"Shape BEFORE grouping: {}\".format(z.shape))\n    z = z.groupby('object_id').mean()\n    print(\"Shape AFTER grouping: {}\".format(z.shape))\n    z.to_csv('single_{}'.format(filename), index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"4b7c981075773cce8b59d0f5baedd9b634f86916"},"cell_type":"code","source":"if __name__ == '__main__':\n    main(len(sys.argv), sys.argv)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"4359471259138987e55feef75286f11e9f9f804b"},"cell_type":"code","source":"","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}