{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# XGBoost Pyramid [CV 0.7968]\n\n## Pipeline:\n1. Feature engineering: https://www.kaggle.com/code/roberthatch/amex-feature-engg-gpu-or-cpu-process-in-chunks\n2. Training: (this notebook). Can also do test predictions if DO_SUBMIT is set to True.\n3. Test Predictions: https://www.kaggle.com/roberthatch/xgboost-pyramid-test-predictions\n\nNote: main reason for pipeline is to streamline GPU use and simplify experimentation.\n\n### Acknowledgements:\n1. https://www.kaggle.com/code/cdeotte/xgboost-starter-0-793\n2. https://www.kaggle.com/code/jiweiliu/rapids-cudf-feature-engineering-xgb \n\n\n### Training: Personal touches:\n1. The pyramid! \n2. Early stopping against auc, not logloss nor amex score.\n3. Lower sampling and learning rate. \n\n\n## About the Pyramid:\nMain goal is to reduce over-specialization without the performance impact of DART on XGBoost GPU.\n\nStart with a large forest with fewer number of rounds and higher learning rate to promote diversity that all the boosting will rely on.\nBetween early layers, rescale the predictions to allow the next layer more room to impact. Again promoting more diversity. Rescale by less and less each layer.\nSlightly lower true learning rate per layer. True learning rate = learning_rate / num_parallel_tree.\n\nPyramid terminology and theory: Purely my own homebrew theory-crafting after reading about DART. If there's existing scholorly articles or prior work in this area, or if blending/stacking/whatever is a more correct term for what I'm doing, please let me know in the comments!","metadata":{}},{"cell_type":"markdown","source":"# Load Libraries","metadata":{}},{"cell_type":"code","source":"# LOAD LIBRARIES\nimport pandas as pd, numpy as np # CPU libraries\nimport matplotlib.pyplot as plt, gc, os\n\nGPU = True\ntry:\n    import cupy, cudf\nexcept ImportError:\n    GPU = False\n\nif GPU:\n    print('RAPIDS version',cudf.__version__)\nelse:\n    print(\"Disabling cudf, using pandas instead\")\n    cudf = pd","metadata":{"execution":{"iopub.status.busy":"2022-07-18T23:24:48.059269Z","iopub.execute_input":"2022-07-18T23:24:48.059704Z","iopub.status.idle":"2022-07-18T23:24:52.820784Z","shell.execute_reply.started":"2022-07-18T23:24:48.059623Z","shell.execute_reply":"2022-07-18T23:24:52.819948Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# VERSION NAME FOR SAVED MODEL FILES\nVER = 1\nFEATURE_VER = 111\n\n# RANDOM SEED\nSEED = 108+5*VER+100*FEATURE_VER\n\n# FOLDS PER MODEL\nFOLDS = 5\n\n# NOTEBOOK PATH\nFEATURE_PATH = '../input/amex-feature-engg-gpu-or-cpu-process-in-chunks/'\n\nDO_SUBMIT = False\n\nprint(\"VER:\", VER)\nprint(\"fVER:\", FEATURE_VER)","metadata":{"execution":{"iopub.status.busy":"2022-07-18T23:24:52.822653Z","iopub.execute_input":"2022-07-18T23:24:52.823256Z","iopub.status.idle":"2022-07-18T23:24:52.829728Z","shell.execute_reply.started":"2022-07-18T23:24:52.823218Z","shell.execute_reply":"2022-07-18T23:24:52.828933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load Data\nFeature engineering is all done in: https://www.kaggle.com/code/roberthatch/amex-feature-engg-gpu-or-cpu-process-in-chunks","metadata":{}},{"cell_type":"code","source":"print('Reading train data...')\nTRAIN_PATH = f'{FEATURE_PATH}train_fe_v{FEATURE_VER}.parquet'\ntrain = pd.read_parquet(TRAIN_PATH)\nprint(train.shape)\n\ntrain = train.sample(frac=1, random_state=SEED)\ntrain = train.reset_index(drop=True)\ntrain.head()","metadata":{"execution":{"iopub.status.busy":"2022-07-18T23:24:52.831279Z","iopub.execute_input":"2022-07-18T23:24:52.831632Z","iopub.status.idle":"2022-07-18T23:25:25.812659Z","shell.execute_reply.started":"2022-07-18T23:24:52.831598Z","shell.execute_reply":"2022-07-18T23:25:25.811741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train XGB\nWe will train using `DeviceQuantileDMatrix`. This has a very small GPU memory footprint.","metadata":{}},{"cell_type":"code","source":"# LOAD XGB LIBRARY\nfrom sklearn.model_selection import KFold\nimport xgboost as xgb\nprint('XGB Version',xgb.__version__)\n\n\n# XGB MODEL PARAMETERS\nBASE_LEARNING_RATE = 0.01\nxgb_params = { \n    'max_depth': 7,\n    'subsample':0.75,\n    'colsample_bytree': 0.35,\n    'gamma':1.5,\n    'lambda':70,\n    'min_child_weight':8,\n\n    'objective':'binary:logistic',\n    'eval_metric':['logloss', 'auc'],  ## Early stopping is based on the last metric listed.\n    'tree_method':'gpu_hist',\n    'predictor':'gpu_predictor',\n    'random_state':SEED,\n\n    'num_parallel_tree':1\n}","metadata":{"execution":{"iopub.status.busy":"2022-07-18T23:25:25.815027Z","iopub.execute_input":"2022-07-18T23:25:25.815457Z","iopub.status.idle":"2022-07-18T23:25:25.930380Z","shell.execute_reply.started":"2022-07-18T23:25:25.815421Z","shell.execute_reply":"2022-07-18T23:25:25.929586Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# NEEDED WITH DeviceQuantileDMatrix BELOW\nclass IterLoadForDMatrix(xgb.core.DataIter):\n    def __init__(self, df=None, features=None, target=None, batch_size=256*1024):\n        self.features = features\n        self.target = target\n        self.df = df\n        self.it = 0 # set iterator to 0\n        self.batch_size = batch_size\n        self.batches = int( np.ceil( len(df) / self.batch_size ) )\n        super().__init__()\n\n    def reset(self):\n        '''Reset the iterator'''\n        self.it = 0\n\n    def next(self, input_data):\n        '''Yield next batch of data.'''\n        if self.it == self.batches:\n            return 0 # Return 0 when there's no more batch.\n        \n        a = self.it * self.batch_size\n        b = min( (self.it + 1) * self.batch_size, len(self.df) )\n        dt = cudf.DataFrame(self.df.iloc[a:b])\n        input_data(data=dt[self.features], label=dt[self.target]) #, weight=dt['weight'])\n        self.it += 1\n        return 1","metadata":{"execution":{"iopub.status.busy":"2022-07-18T23:25:25.931549Z","iopub.execute_input":"2022-07-18T23:25:25.932380Z","iopub.status.idle":"2022-07-18T23:25:25.941821Z","shell.execute_reply.started":"2022-07-18T23:25:25.932344Z","shell.execute_reply":"2022-07-18T23:25:25.941003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## TODO, replace with newer version.\n\n# https://www.kaggle.com/kyakovlev\n# https://www.kaggle.com/competitions/amex-default-prediction/discussion/327534\ndef amex_metric_mod(y_true, y_pred):\n\n    labels     = np.transpose(np.array([y_true, y_pred]))\n    labels     = labels[labels[:, 1].argsort()[::-1]]\n    weights    = np.where(labels[:,0]==0, 20, 1)\n    cut_vals   = labels[np.cumsum(weights) <= int(0.04 * np.sum(weights))]\n    top_four   = np.sum(cut_vals[:,0]) / np.sum(labels[:,0])\n\n    gini = [0,0]\n    for i in [1,0]:\n        labels         = np.transpose(np.array([y_true, y_pred]))\n        labels         = labels[labels[:, i].argsort()[::-1]]\n        weight         = np.where(labels[:,0]==0, 20, 1)\n        weight_random  = np.cumsum(weight / np.sum(weight))\n        total_pos      = np.sum(labels[:, 0] *  weight)\n        cum_pos_found  = np.cumsum(labels[:, 0] * weight)\n        lorentz        = cum_pos_found / total_pos\n        gini[i]        = np.sum((lorentz - weight_random) * weight)\n\n    print(\"  4%  :\", top_four)\n    print(\"  Gini:\", gini[1]/gini[0])\n    print(\"Kaggle:\", 0.5 * (gini[1]/gini[0] + top_four))\n    return 0.5 * (gini[1]/gini[0] + top_four)","metadata":{"execution":{"iopub.status.busy":"2022-07-18T23:25:25.944892Z","iopub.execute_input":"2022-07-18T23:25:25.945165Z","iopub.status.idle":"2022-07-18T23:25:25.959580Z","shell.execute_reply.started":"2022-07-18T23:25:25.945141Z","shell.execute_reply":"2022-07-18T23:25:25.958728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"importances = []\nPYRAMID_W = [0.5, 2/3, 0.75, 0.875, 1, 0]\n\ndef run_training(train, features):\n    oof = []\n\n    skf = KFold(n_splits=FOLDS)\n    for fold,(train_idx, valid_idx) in enumerate(skf.split(\n                train, train.target )):\n        print('#'*25)\n        print('### Fold',fold+1)\n    \n        # TRAIN, VALID, TEST FOR FOLD K\n        Xy_train = IterLoadForDMatrix(train.loc[train_idx], features, 'target')\n        X_valid = train.loc[valid_idx, features]\n        y_valid = train.loc[valid_idx, 'target']\n\n        print('### Train size',len(train_idx),'Valid size',len(valid_idx),'Valid positives',y_valid.sum())\n        print(f'### Training with all of fold data...')\n        print('#'*25)\n\n        dtrain = xgb.DeviceQuantileDMatrix(Xy_train, max_bin=256)\n        dvalid = xgb.DMatrix(data=X_valid, label=y_valid)\n\n        # TRAIN MODEL FOLD K\n        # PYRAMID: Smoothly go from diverse forest of early trees into focused boosted trees correcting residuals.\n        #   final layer must have w==0\n        #   columns:    forest|boost|adj_eta|w\n        pyramid_layers = [(100,  10,  1.56,  0.5),\n                          ( 20,  50,  1.3,   2/3),\n                          (  1,1000,  1.25,  0.75),\n                          (  1,1000,  1.125, 0.875),\n                          (  1,3000,  1.0,   1),\n                          (  1,9000,  0.5,   0)]\n        assert(PYRAMID_W == [layer[-1] for layer in pyramid_layers])\n        for (layer, (n_trees, n_rounds, adj_learning, w)) in enumerate(pyramid_layers):\n            ## Load the manual parameters from the pyramid layer\n            xgb_params['num_parallel_tree'] = n_trees\n            xgb_params['learning_rate'] = n_trees*adj_learning*BASE_LEARNING_RATE\n            xgb_params['random_state'] += 1\n            \n            ## No early stopping except on final round. This is important since the weighting causes the model to go backwards for a time at the start of the next layer.\n            early_stop = None\n            if w == 0:\n                early_stop = 300\n\n            print(\"Learning Rate:\", xgb_params['learning_rate'])\n            model = xgb.train(xgb_params, \n                        dtrain=dtrain,\n                        evals=[(dtrain,'train'),(dvalid,'valid')],\n                        num_boost_round=n_rounds,\n                        early_stopping_rounds=early_stop,\n                        verbose_eval=100//n_trees)\n            ## save model layer here\n            model.save_model(f'XGB_v{VER}_fold{fold}_layer{layer}.xgb')\n\n            ## predict to load the predictions on the next model layer\n            ## Don't set base margin on final layer. w = 0 is used as an encoded way to skip this step.\n            if (w != 0):\n                ptrain = model.predict(dtrain, output_margin=True)\n                pvalid = model.predict(dvalid, output_margin=True)\n\n                ## reduce the impact of all model layers so far by w. This should be another way to reduce over-specialization, without the computational cost of DART\n                if (w < 1.0):\n                    ptrain = ptrain * w\n                    pvalid = pvalid * w\n\n                ## This set_base_margin on the DMatrix data is what informs the next layer of the prior training.\n                ## See code example from official demos: https://github.com/dmlc/xgboost/blob/master/demo/guide-python/boost_from_prediction.py\n                dtrain.set_base_margin(ptrain)\n                dvalid.set_base_margin(pvalid)\n\n                plt.hist(pvalid, bins=100)\n                plt.title(f'Layer {layer} OOF Predictions')\n                plt.show()\n\n                del model, ptrain, pvalid\n                gc.collect()\n\n        # GET FEATURE IMPORTANCE FOR FOLD K\n        dd = model.get_score(importance_type='weight')\n        df = pd.DataFrame({'feature':dd.keys(),f'importance_{fold}':dd.values()})\n        importances.append(df)\n\n        # INFER OOF FOLD K\n        # Note: Not necessary with current implementation having final pyramid layer with num_parallel_tree == 1, but more robust to divide best_ntree_limit\n        #   by num_parallel_tree. Oddly, iteration range is based only on num_boost_rounds, but best_ntree_limit is stored as num_boost_rounds * num_parallel_trees\n        print(\"Best_ntree_limit:\", model.best_ntree_limit//xgb_params['num_parallel_tree'])\n        oof_preds = model.predict(dvalid, iteration_range=(0,model.best_ntree_limit//xgb_params['num_parallel_tree']))\n        print('For this fold:')\n        ## TODO: update metric. Fork this notebook to confirm the latest version of the numpy implementation from author is even faster and equally accurate.\n        ## https://www.kaggle.com/code/rohanrao/amex-competition-metric-implementations\n        amex_metric_mod(y_valid.values, oof_preds)\n    \n        # SAVE OOF\n        df = train.loc[valid_idx, ['customer_ID','target'] ].copy()\n        df['oof_pred'] = oof_preds\n        oof.append( df )\n\n        del dtrain, Xy_train, dd, df\n        del X_valid, y_valid, dvalid, model\n        gc.collect()\n\n    print('#'*25)\n    print('OVERALL CV:')\n    oof = pd.concat(oof,axis=0,ignore_index=True).set_index('customer_ID')\n    amex_metric_mod(oof.target.values, oof.oof_pred.values)\n    return oof\n","metadata":{"scrolled":true,"execution":{"iopub.status.busy":"2022-07-18T23:25:25.962636Z","iopub.execute_input":"2022-07-18T23:25:25.962940Z","iopub.status.idle":"2022-07-18T23:25:25.988773Z","shell.execute_reply.started":"2022-07-18T23:25:25.962917Z","shell.execute_reply":"2022-07-18T23:25:25.987961Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"features = train.columns[1:-1]\nprint(f'There are {len(features)} features!')\nprint(train.shape)\n\noof = run_training(train, features)","metadata":{"execution":{"iopub.status.busy":"2022-07-18T23:25:25.993730Z","iopub.execute_input":"2022-07-18T23:25:25.997434Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# CLEAN RAM\ndel train\n_ = gc.collect()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Save OOF Preds","metadata":{}},{"cell_type":"code","source":"oof_xgb = pd.read_parquet(TRAIN_PATH, columns=['customer_ID']).drop_duplicates()\noof_xgb = oof_xgb.set_index('customer_ID')\noof_xgb = oof_xgb.merge(oof, left_index=True, right_index=True)\noof_xgb = oof_xgb.sort_index().reset_index(drop=True)\noof_xgb.to_csv(f'oof_xgb_v{VER}.csv',index=False)\noof_xgb.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# PLOT OOF PREDICTIONS\nplt.hist(oof_xgb.oof_pred.values, bins=100)\nplt.title('OOF Predictions')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# CLEAR VRAM, RAM FOR INFERENCE BELOW\ndel oof_xgb, oof\n_ = gc.collect()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Feature Importance","metadata":{}},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n\ndf = importances[0].copy()\nfor k in range(1,FOLDS): df = df.merge(importances[k], on='feature', how='left')\ndf['importance'] = df.iloc[:,1:].mean(axis=1)\ndf = df.sort_values('importance',ascending=False)\ndf.to_csv(f'xgb_feature_importance_v{VER}.csv',index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"NUM_FEATURES = 30\nplt.figure(figsize=(10,5*NUM_FEATURES//10))\nplt.barh(np.arange(NUM_FEATURES,0,-1), df.importance.values[:NUM_FEATURES])\nplt.yticks(np.arange(NUM_FEATURES,0,-1), df.feature.values[:NUM_FEATURES])\nplt.title(f'XGB Feature Importance - Top {NUM_FEATURES}')\nplt.show()\n\ndel df, importances","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Infer Test","metadata":{"_kg_hide-output":true,"_kg_hide-input":true}},{"cell_type":"code","source":"if DO_SUBMIT:\n    gc.collect()\n\n    # INFER TEST DATA IN PARTS\n\n    TEST_SECTIONS = 2\n    TEST_SUB_SECTIONS = 2\n\n    test_preds = []\n    customers = False\n    for k in range(TEST_SECTIONS):\n        for i in range(TEST_SUB_SECTIONS):    \n            # READ PART OF TEST DATA\n            print(f'\\nReading test data...')\n            test = cudf.read_parquet(f'{FEATURE_PATH}test{k}_fe_v{FEATURE_VER}.parquet')\n            if i == 0:\n                print(f'=> Test part {k+1} has shape', test.shape )\n                if k == 0:\n                    customers = test.index.copy()\n                else:\n                    customers = customers.append(test.index)\n\n            # TEST DATA FOR XGB\n            X_test = test[features]\n            n_rows = len(test.index)//TEST_SUB_SECTIONS\n            print(\".\")\n            if i+1 < TEST_SUB_SECTIONS:\n                X_test = X_test.iloc[i*n_rows:(i+1)*n_rows, :].copy()\n            elif TEST_SUB_SECTIONS > 1:\n                X_test = X_test.iloc[i*n_rows:, :].copy()\n            print(f'=> Test piece {k+1}, {i+1} has shape', X_test.shape )\n            del test\n            gc.collect()\n            dtest = xgb.DMatrix(data=X_test)\n            del X_test\n            gc.collect()\n            ## Need to reset to level 0 between folds.\n            reset_margin = dtest.get_base_margin()\n\n            # INFER XGB MODELS ON TEST DATA\n            print(\".\")\n            for f in range(FOLDS):\n                if (f > 0):\n                    dtest.set_base_margin(reset_margin)\n                for (layer, w) in enumerate(PYRAMID_W[:-1]):\n                    model = xgb.Booster()\n                    model.load_model(f'XGB_v{VER}_fold{f}_layer{layer}.xgb')\n                    print(f'Loaded fold{f}, layer{layer}')\n                    ptest = model.predict(dtest, output_margin=True)\n\n                    ## reduce the impact of all model layers so far by w. This should be another way to reduce over-specialization, without the computational cost of DART\n                    if (w < 1.0):\n                        ptest = ptest * w\n\n                    ## This set_base_margin is what informs the next layer of the prior training.\n                    ## See code example from official demos: https://github.com/dmlc/xgboost/blob/master/demo/guide-python/boost_from_prediction.py\n                    dtest.set_base_margin(ptest)\n\n                layer = len(PYRAMID_W) - 1\n                model = xgb.Booster()\n                model.load_model(f'XGB_v{VER}_fold{f}_layer{layer}.xgb')\n                print(\"Best_ntree_limit\", model.best_ntree_limit//xgb_params['num_parallel_tree'])\n                if f == 0:\n                    preds = model.predict(dtest, output_margin=True, iteration_range=(0,model.best_ntree_limit//xgb_params['num_parallel_tree']))\n                else:\n                    preds += model.predict(dtest, output_margin=True, iteration_range=(0,model.best_ntree_limit//xgb_params['num_parallel_tree']))\n            preds /= FOLDS\n            test_preds.append(preds)\n\n            # CLEAN MEMORY\n            del dtest, model, reset_margin\n            _ = gc.collect()","metadata":{"_kg_hide-output":false,"_kg_hide-input":false,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create Submission CSV","metadata":{"_kg_hide-output":false,"_kg_hide-input":false}},{"cell_type":"code","source":"if DO_SUBMIT:\n    # WRITE SUBMISSION FILE\n    test_preds = np.concatenate(test_preds)\n    test = cudf.DataFrame(index=customers,data={'prediction':test_preds})\n    sub = cudf.read_csv('../input/amex-default-prediction/sample_submission.csv')[['customer_ID']]\n    sub['customer_ID_hash'] = sub['customer_ID'].str[-16:].str.hex_to_int().astype('int64')\n    sub = sub.set_index('customer_ID_hash')\n    sub = sub.merge(test[['prediction']], left_index=True, right_index=True, how='left')\n    sub = sub.reset_index(drop=True)\n\n\n    # DISPLAY PREDICTIONS\n    sub.to_csv(f'submission.csv',index=False)\n    print('Submission file shape is', sub.shape )\n    sub.head()\n\n    # PLOT PREDICTIONS\n    plt.hist(sub.to_pandas().prediction, bins=100)\n    plt.title('Test Predictions')\n    plt.show()","metadata":{"_kg_hide-output":false,"_kg_hide-input":false,"trusted":true},"execution_count":null,"outputs":[]}]}