{"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":"# Feature Engineer","metadata":{}},{"cell_type":"code","source":"# LOAD LIBRARIES\nimport pandas as pd, numpy as np # CPU libraries\nimport 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-28T12:09:40.485485Z","iopub.execute_input":"2022-07-28T12:09:40.486230Z","iopub.status.idle":"2022-07-28T12:09:44.202019Z","shell.execute_reply.started":"2022-07-28T12:09:40.486133Z","shell.execute_reply":"2022-07-28T12:09:44.200231Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"PROCESS_TEST_DATA = True\n\n# VERSION NAME FOR SAVED PARQUET FILES\nVER = 111\n\n# FILL NAN VALUE\nNAN_VALUE = -127 # will fit in int8\n\nif GPU:\n    TRAIN_NUM_PARTS = 2\n    TEST_SECTIONS = 2\n    TEST_NUM_PARTS = 2\nelse:\n    TRAIN_NUM_PARTS = 6\n    TEST_SECTIONS = 2\n    TEST_NUM_PARTS = 6\n\nprint(\"VER:\", VER)\nif not PROCESS_TEST_DATA:\n    print(\"NOT processing test data!\")\n    \n    \nLABEL_DATA_PATH = \"../input/boolart-user-default-prediction/train_labels.csv\"\nTRAIN_PATH = '../input/boolart-user-default-prediction/train.parquet'\nTEST_PATH = '../input/boolart-user-default-prediction/test.parquet'\nOUTPUT_PATH = \"./feat/\"\nfrom pathlib import Path\nPath(OUTPUT_PATH).mkdir(parents=True, exist_ok=True)","metadata":{"execution":{"iopub.status.busy":"2022-07-28T12:09:44.204111Z","iopub.execute_input":"2022-07-28T12:09:44.204444Z","iopub.status.idle":"2022-07-28T12:09:44.212277Z","shell.execute_reply.started":"2022-07-28T12:09:44.204408Z","shell.execute_reply":"2022-07-28T12:09:44.211177Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Basic reading and formatting from initial raddar parquet file\n\ndef process_customer_columns(df):\n    # REDUCE DTYPE FOR CUSTOMER AND DATE\n    if GPU:\n        df['customer_ID'] = df['customer_ID'].str[-16:].str.hex_to_int().astype('int64')\n    else:\n        df['customer_ID'] = df['customer_ID'].str[-16:].apply(int, base=16).astype('int64')\n    year = cudf.to_numeric(df['S_2'].str[:4])\n    month = cudf.to_numeric(df['S_2'].str[5:7])\n    df['S_2'] = year.mul(12).add(month).sub(24207).astype('int8')\n    return df\n\ndef read_file(path='', usecols=None):\n    # LOAD DATAFRAME\n    if usecols is not None:\n        df = cudf.read_parquet(path, columns=usecols)\n        df = process_customer_columns(df)\n    else:\n        df = cudf.read_parquet(path)\n\n    print('Shape of data:', df.shape)\n    \n    return df","metadata":{"execution":{"iopub.status.busy":"2022-07-28T12:09:44.214029Z","iopub.execute_input":"2022-07-28T12:09:44.214465Z","iopub.status.idle":"2022-07-28T12:09:44.229342Z","shell.execute_reply.started":"2022-07-28T12:09:44.214423Z","shell.execute_reply":"2022-07-28T12:09:44.228449Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## CALCULATE SIZE OF EACH SEPARATE PART\ndef get_rows(customers, df, NUM_PARTS=4, verbose=''):\n    chunk = len(customers)//NUM_PARTS\n    if verbose != '':\n        print(f'We will process {verbose} data as {NUM_PARTS} separate parts.')\n        print(f'There will be {chunk} customers in each part (except the last part).')\n        print('Below are number of rows in each part:')\n    rows = []\n\n    for k in range(NUM_PARTS):\n        if k==NUM_PARTS-1: cc = customers[k*chunk:]\n        else: cc = customers[k*chunk:(k+1)*chunk]\n        s = df.loc[df.customer_ID.isin(cc)].shape[0]\n        rows.append(s)\n    if verbose != '': print( rows )\n    return rows,chunk\n\ndef getAndProcessDataInChunks(filename, is_train=False, NUM_PARTS=4, NUM_SECTIONS=1, split_k=0, verbose=''):\n    gc.collect()\n\n    print(f'Reading customer_IDs from {verbose} data...')\n    df = read_file(path = filename, usecols = ['customer_ID','S_2'])\n    customers = df[['customer_ID']].drop_duplicates().sort_index().values.flatten()\n    rows,num_cust = get_rows(customers, df[['customer_ID']], NUM_PARTS=NUM_PARTS*NUM_SECTIONS, verbose=verbose)\n\n    # INFER DATA IN PARTS\n    skip_rows = 0\n    skip_cust = 0\n    allData = []\n\n    del df\n    gc.collect()\n\n    print(f'\\nReading {verbose} data...')\n    df_file = read_file(path = filename)\n\n    if is_train:\n        assert(NUM_SECTIONS == 1) ## Splitting not implemented for target labels\n        targets = cudf.read_csv(LABEL_DATA_PATH)\n        if GPU:\n            targets['customer_ID'] = targets['customer_ID'].str[-16:].str.hex_to_int().astype('int64')\n        else:\n            targets['customer_ID'] = targets['customer_ID'].str[-16:].apply(int, base=16).astype('int64')\n        targets = targets.set_index('customer_ID')\n        targets.target = targets.target.astype('int8')\n\n    if NUM_SECTIONS > 1:\n        startRow = 0\n        for i in range(NUM_SECTIONS):\n            if i == split_k:\n                startRow = skip_rows\n            for k in range(NUM_PARTS):\n                skip_rows += rows[i*NUM_PARTS + k]\n            if i == split_k:\n                df_file = df_file.iloc[startRow:skip_rows].reset_index(drop=True)\n                rows = rows[i*NUM_PARTS:(i+1)*NUM_PARTS]\n                gc.collect()\n                skip_rows = 0\n                break\n    for k in range(NUM_PARTS):\n        # READ PART OF DATA\n        df = df_file.iloc[skip_rows:skip_rows+rows[k]].reset_index(drop=True)\n        skip_rows += rows[k]\n        print(f'=> {verbose} part {k+1} has shape', df.shape )\n\n        # PROCESS AND FEATURE ENGINEER PART OF DATA\n        df = process_and_feature_engineer(df)\n\n        if is_train:\n            ## Relies on assumption that initial train data has customer IDs in same sorted order as train_labels.csv\n            if k==NUM_PARTS-1: targetSlice = targets.iloc[skip_cust:]\n            else: targetSlice = targets.iloc[skip_cust:skip_cust+num_cust]\n            skip_cust += num_cust\n\n            print(\"|...\")\n            df = cudf.concat([df, targetSlice], axis=1)\n            print(\" ...|\")\n\n        if GPU:\n            print(\"|...\")\n            df = df.to_pandas()\n            print(\" ...|\")\n\n        allData.append(df)\n        gc.collect()\n\n    print(\".\", end='')\n    del df_file\n    gc.collect()\n    allData = pd.concat(allData, axis=0)\n    del df\n    gc.collect()\n    if is_train:\n        print(\".\", end='')\n        allData = allData.sort_index()\n        gc.collect()\n        print(\".\", end='')\n        allData = allData.reset_index()\n    print(\"|\")\n    return allData","metadata":{"execution":{"iopub.status.busy":"2022-07-28T12:09:44.231341Z","iopub.execute_input":"2022-07-28T12:09:44.231658Z","iopub.status.idle":"2022-07-28T12:09:44.260604Z","shell.execute_reply.started":"2022-07-28T12:09:44.231626Z","shell.execute_reply":"2022-07-28T12:09:44.259651Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"##\n## Easy to replace or modify this function with your own custom feature engineering,\n## and keep the other boilerplate untouched to allow switching from GPU to CPU based on GPU quota.\n##\ndef process_and_feature_engineer(df):\n    print(\".\", end = '')\n\n    ## Save space on customer ID, and encode S_2 based on month and year as 0-12 for train set, 13-25 or 19-31 for test set.\n    df = process_customer_columns(df)\n\n    print(\".\", end = '')\n\n    df.drop(\"B_29\",inplace=True,axis=1)\n\n    # compute \"after pay\" features\n    for bcol in [f'B_{i}' for i in [11,14,17]]+['D_39','D_131']+[f'S_{i}' for i in [16,23]]:\n        for pcol in ['P_2','P_3']:\n            if bcol in df.columns:\n                result = df[bcol] - df[pcol]\n                df[f'{bcol}-{pcol}'] = result.fillna(0)\n\n    # FEATURE ENGINEERING heavily modified, started from: \n    all_cols = [c for c in list(df.columns) if c not in ['customer_ID']]\n    cat_features = [\"B_30\",\"B_38\",\"D_114\",\"D_116\",\"D_117\",\"D_120\",\"D_126\",\"D_63\",\"D_64\",\"D_66\",\"D_68\"]\n    num_features = [col for col in all_cols if col not in (cat_features + [\"S_2\"])]\n\n    print(\".\", end = '')\n\n    ## For each customer, count all NaN in any row, and count all NaN in the last row. Later we will add it as two columns\n    df_nan = (df.mul(0) + 1).fillna(0)\n    df_nan['customer_ID'] = df['customer_ID']\n    nan_sum = df_nan.groupby(\"customer_ID\").sum().sum(axis=1)\n    nan_last = df_nan.groupby(\"customer_ID\").last().sum(axis=1)\n    del df_nan\n    print(\".\", end = '')\n\n    groups = df.groupby(\"customer_ID\")\n    test_num_agg = groups[num_features].agg(['mean', 'std', 'last'])\n    test_num_agg.columns = ['_'.join(x) for x in test_num_agg.columns]\n\n    print(\"+\", end = '')\n\n    ## TODO: One-hot encode or convert to non-numeric? I believe raddar's clean dataset (and original data?) stores as numeric,\n    ## and XGBoost doesn't(?) have a way to explicitly call out categorical columns\n    test_cat_agg = groups[cat_features].last()\n    test_cat_agg.columns = [x + \"_last\" for x in test_cat_agg.columns]\n\n    print(\".\", end = '')\n\n    ## S_2 test data has different range of values for S_2, normalize 'min' by subtracting max-12. Aka min = min + 12 - max.\n    ## S_2 max will always be the same as last, and (after normalization), the same for every customer.\n    ## Min tells us if they are a shorter term customer or not. Count under 13 tells us they're EITHER a shorter term customer OR they have some gap months.\n    ## TODO: Min is more relevant (much higher default rate for short term customers than for gap customers), but might be worth encoding both into single column.\n    ##   Basically, count + 143 minus 12*min (normal customer gets 156, long term gap customer gets 145-155, short term customer gets 0-143)\n    ##   If not combining into one, could arguably use count+min instead of count, to directly highlight gap customers.\n    test_s2_agg = groups[[\"S_2\"]].agg(['min', 'count', 'max'])\n    test_s2_agg.columns = ['_'.join(x) for x in test_s2_agg.columns]\n    test_s2_agg['S_2_min'] = test_s2_agg['S_2_min'] + 12 - test_s2_agg['S_2_max']\n    test_s2_agg.drop(['S_2_max'],inplace=True,axis=1)\n\n    ## Quick sanity sorting check: confirm last value of each group is from the last (max) statement by checking S_2:\n    assert_s2_check = groups[\"S_2\"].last() - groups[\"S_2\"].max()\n    assert((assert_s2_check == 0).all())\n\n    ## Drop delta for now. I haven't had success getting 1500+ features without running out of memory.\n    ##   Out of memory: Not only while feature engineering (probably solvable), but again while running the model from pre-processed feature dataset (harder to solve)\n    ## TODO: try swapping out other feature to add this one in?\n    ##   Feature selection on subsets, then combine best features?\n    ##   Reduce float64, int64, etc, to lower precision and fit more that way?\n#     ### Add delta\n#     test_num_agg2 = groups[num_features].nth(-1) - groups[num_features].nth(-2)\n#     test_num_agg2 = test_num_agg2.fillna(0)\n#     test_num_agg2.columns = [x + '_delta' for x in test_num_agg2.columns]\n\n    ## Note for optimization: several of these are probably super inefficient to re-calculate it from the group data. Should be able to re-use test_num_agg in some way.\n\n    ### Add current level: range from 0.0-1.0. For example: 0.0 means last = min; 0.5 means last = (max+min)/2; 1.0 means last = max.\n    test_num_agg2 = (groups[num_features].last() - groups[num_features].min()) / (groups[num_features].max() - groups[num_features].min())\n    test_num_agg2 = test_num_agg2.fillna(0)\n    test_num_agg2.columns = [x + '_curLevel' for x in test_num_agg2.columns]\n\n    ### Add magnitude: max - min\n    test_num_agg3 = groups[num_features].max() - groups[num_features].min()\n    test_num_agg3 = test_num_agg3.fillna(0)\n    test_num_agg3.columns = [x + '_magnitude' for x in test_num_agg3.columns]\n\n    ### Add last-mean\n    test_num_agg4 = groups[num_features].last() - groups[num_features].mean()\n    test_num_agg4 = test_num_agg4.fillna(0)\n    test_num_agg4.columns = [x + '_last-mean' for x in test_num_agg4.columns]\n\n    ### Add match for categorical: 1 if last and next to last are the same, 0 if not. If both are nan, or next to last is nan, treat as the same via fillna(1)\n    ## TODO: in case it matters, move below forwards/backwards fill, just below.\n    test_cat_agg2 = 1 + groups[cat_features].nth(-1) - groups[cat_features].nth(-2)\n    test_cat_agg2 = test_cat_agg2.fillna(1)\n    test_cat_agg2[test_cat_agg2 != 1] = 0\n    test_cat_agg2.columns = [x + '_match' for x in test_cat_agg2.columns]\n\n    print(\".\", end = '')\n\n    ## Forward fill, and then backward fill to remove all nans before using \"nth\" to calc hma\n    groups = groups.ffill()\n    groups[\"customer_ID\"] = df[\"customer_ID\"]\n    groups = groups.groupby(\"customer_ID\").bfill()\n    groups[\"customer_ID\"] = df[\"customer_ID\"]\n    groups = groups.groupby(\"customer_ID\")\n\n    print(\".\", end = '')\n\n    ### Add HMA. Hull moving average is a smoothed moving average sometimes used on time series data. (e.g. stock price).\n    ## I happen to like it. As usual, I can't actually say whether it helps a lot, a little, or not at all.\n    ##   Especially with the inherent CV variance, I haven't done nearly enough (any) experiments to see whether it helps a lot, a little, or not at all.\n    ##   It didn't obviously help compared with other XGB public notebooks until I lowered column subsampling and learning rate.\n    ##   With those hyper parameter changes helping a ton, I haven't checked which personal feature engg touches are actually contributing (if any).\n    ## For calculation simplicity, this version of HMA is a bit simpler, while keeping the key idea that mean(range(10)) = 5, but hma(range(10)) = 10.\n    ## TODO: Maybe consider calculating reverse HMA, for example last - reverse HMA, or HMA-reverseHMA to more directly look at slope of the data.\n    ##\n    ## Calculate the biggest down to the smallest, then take the biggest that isn't nan, to handle varying group sizes\n    h13 = hma13(groups, num_features)\n    h11 = hma11(groups, num_features)\n    h9 = hma9(groups, num_features)\n    h7 = hma7(groups, num_features)\n    h5 = hma5(groups, num_features)\n    h3 = hma3(groups, num_features)\n    h1 = hma1(groups, num_features)\n    hma_df = cudf.concat([h1, h3, h5, h7, h9, h11, h13], axis=0)\n    del h1, h3, h5, h7, h9, h11, h13, assert_s2_check\n    gc.collect()\n    print(\".\", end = '')\n    hma_df = hma_df.sort_index()\n    hma_df = hma_df.reset_index()\n    hma_df = hma_df.groupby(\"customer_ID\").last()\n    hma_df.columns = [x + \"_hma\" for x in hma_df.columns]\n\n    print(\".\", end = '')\n\n    df = cudf.concat([test_s2_agg, test_num_agg, test_cat_agg, test_num_agg2, test_cat_agg2, hma_df, test_num_agg3, test_num_agg4], axis=1)\n\n    print(\".\")\n\n    ## Finally add NaN counts from earlier\n    df[\"total_data_count\"] = nan_sum\n    df[\"total_data_last\"] = nan_last\n\n    ## Per discussion I forgot to save the link, maybe on XGBoost Starter notebook, there's two columns with actual numbers often going below '-127', the default fillna.\n    ## TODO is to handle some categories of nan differently anyways (per raddar's work), and agg often ignores them (which I think is a fine way) anyways.\n    ##   The remaining nans: it's probably better in XGBoost to just leave them there, and let XGBoost decide!\n    nan_col = ['D_50_mean', 'D_50_std', 'D_50_last', 'S_23_mean', 'S_23_std', 'S_23_last']\n    df[nan_col] = df[nan_col].fillna(-32783)\n    df = df.fillna(NAN_VALUE)\n\n    print('shape after engineering', df.shape )\n    for col in df.columns:\n        if len(df[col].unique()) == 1:\n            print(\"Consider dropping column!?\", col)\n\n    del test_s2_agg, test_num_agg, test_cat_agg\n    del test_num_agg2, test_cat_agg2, test_num_agg3, test_num_agg4\n\n    return df\n\n\ndef hma13(groups, columns):\n    return (5/14)*groups[columns].nth(-1) + (27/91)*groups[columns].nth(-2) + (43/182)*groups[columns].nth(-3) + (16/91)*groups[columns].nth(-4) + (3/26)*groups[columns].nth(-5) + (5/91)*groups[columns].nth(-6) - (1/182)*groups[columns].nth(-7) - (6/91)*groups[columns].nth(-8) - (5/91)*groups[columns].nth(-9) - (4/91)*groups[columns].nth(-10) - (3/91)*groups[columns].nth(-11) - (2/91)*groups[columns].nth(-12) - (1/91)*groups[columns].nth(-13)\ndef hma11(groups, columns):\n    return (17/42)*groups[columns].nth(-1) + (25/77)*groups[columns].nth(-2) + (113/462)*groups[columns].nth(-3) + (38/231)*groups[columns].nth(-4) + (13/154)*groups[columns].nth(-5) + (1/231)*groups[columns].nth(-6) - (5/66)*groups[columns].nth(-7) - (2/33)*groups[columns].nth(-8) - (1/22)*groups[columns].nth(-9) - (1/33)*groups[columns].nth(-10) - (1/66)*groups[columns].nth(-11)\ndef hma9(groups, columns):\n    return (7/15)*groups[columns].nth(-1) + (16/45)*groups[columns].nth(-2) + (11/45)*groups[columns].nth(-3) + (2/15)*groups[columns].nth(-4) + (1/45)*groups[columns].nth(-5) - (4/45)*groups[columns].nth(-6) - (1/15)*groups[columns].nth(-7) - (2/45)*groups[columns].nth(-8) - (1/45)*groups[columns].nth(-9)\ndef hma7(groups, columns):\n    return (11/20)*groups[columns].nth(-1) + (27/70)*groups[columns].nth(-2) + (31/140)*groups[columns].nth(-3) + (2/35)*groups[columns].nth(-4) - (3/28)*groups[columns].nth(-5) - (1/14)*groups[columns].nth(-6) - (1/28)*groups[columns].nth(-7)\ndef hma5(groups, columns):\n    return (2/3)*groups[columns].nth(-1) + (2/5)*groups[columns].nth(-2) + (2/15)*groups[columns].nth(-3) - (2/15)*groups[columns].nth(-4) - (1/15)*groups[columns].nth(-5)\ndef hma3(groups, columns):\n    return (5/6)*groups[columns].nth(-1) + (1/3)*groups[columns].nth(-2) - (1/6)*groups[columns].nth(-3)\ndef hma1(groups, columns):\n    return groups[columns].nth(-1)","metadata":{"execution":{"iopub.status.busy":"2022-07-28T12:09:44.262148Z","iopub.execute_input":"2022-07-28T12:09:44.262428Z","iopub.status.idle":"2022-07-28T12:09:44.312925Z","shell.execute_reply.started":"2022-07-28T12:09:44.262405Z","shell.execute_reply":"2022-07-28T12:09:44.312038Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = getAndProcessDataInChunks(TRAIN_PATH, is_train=True, NUM_PARTS=TRAIN_NUM_PARTS, verbose='train')\n\nprint(train.shape)\nprint(train.head())\ntrain.to_parquet(OUTPUT_PATH+f'train_fe_v{VER}.parquet')\nprint(\"done\")\n\n## ~2 minutes GPU execution time with 1 part, versus\n## ~2.5-3 minutes GPU execution time with 4 parts\n\n## ~9 minutes CPU execution time with 4 parts","metadata":{"execution":{"iopub.status.busy":"2022-07-28T12:09:44.316358Z","iopub.execute_input":"2022-07-28T12:09:44.316882Z","iopub.status.idle":"2022-07-28T12:10:26.608645Z","shell.execute_reply.started":"2022-07-28T12:09:44.316854Z","shell.execute_reply":"2022-07-28T12:10:26.607626Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if PROCESS_TEST_DATA:\n    del train\n    gc.collect()\n\n    \n    for k in range(TEST_SECTIONS):\n        test = getAndProcessDataInChunks(TEST_PATH, NUM_PARTS=TEST_NUM_PARTS, NUM_SECTIONS=TEST_SECTIONS, split_k=k, verbose='test')\n\n        print(test.shape)\n        print(test.head())\n        test.to_parquet(OUTPUT_PATH+f'test{k}_fe_v{VER}.parquet')\n        print(\"done\")\n\n        del test\n        gc.collect()\n\n## ~4 minutes GPU: 2*1 parts\n## ~5-6 minutes GPU: 2*4 parts\n\n## ~20 minutes CPU: 2*4","metadata":{"execution":{"iopub.status.busy":"2022-07-28T12:10:26.610219Z","iopub.execute_input":"2022-07-28T12:10:26.610582Z","iopub.status.idle":"2022-07-28T12:11:28.485397Z","shell.execute_reply.started":"2022-07-28T12:10:26.610546Z","shell.execute_reply":"2022-07-28T12:11:28.484263Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Xgboost pyramid layers","metadata":{}},{"cell_type":"code","source":"# LOAD LIBRARIES\nimport os\nos.environ['CUDA_VISIBLE_DEVICES'] = \"1\"\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-28T12:11:28.488581Z","iopub.execute_input":"2022-07-28T12:11:28.488896Z","iopub.status.idle":"2022-07-28T12:11:28.495167Z","shell.execute_reply.started":"2022-07-28T12:11:28.488871Z","shell.execute_reply":"2022-07-28T12:11:28.494118Z"},"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 = './feat/'\n\nDO_SUBMIT = True\n\nprint(\"VER:\", VER)\nprint(\"fVER:\", FEATURE_VER)","metadata":{"execution":{"iopub.status.busy":"2022-07-28T12:11:28.496560Z","iopub.execute_input":"2022-07-28T12:11:28.497144Z","iopub.status.idle":"2022-07-28T12:11:28.508788Z","shell.execute_reply.started":"2022-07-28T12:11:28.497108Z","shell.execute_reply":"2022-07-28T12:11:28.507868Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-28T12:11:28.510520Z","iopub.execute_input":"2022-07-28T12:11:28.510880Z","iopub.status.idle":"2022-07-28T12:11:29.861597Z","shell.execute_reply.started":"2022-07-28T12:11:28.510845Z","shell.execute_reply":"2022-07-28T12:11:29.860625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-28T12:11:29.862850Z","iopub.execute_input":"2022-07-28T12:11:29.863846Z","iopub.status.idle":"2022-07-28T12:11:29.964130Z","shell.execute_reply.started":"2022-07-28T12:11:29.863805Z","shell.execute_reply":"2022-07-28T12:11:29.963249Z"},"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-28T11:55:39.024009Z","iopub.execute_input":"2022-07-28T11:55:39.024987Z","iopub.status.idle":"2022-07-28T11:55:39.034737Z","shell.execute_reply.started":"2022-07-28T11:55:39.024939Z","shell.execute_reply":"2022-07-28T11:55:39.033576Z"},"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_train = train.loc[train_idx, features]\n        y_train = train.loc[train_idx, 'target']\n        X_valid = train.loc[valid_idx, features]\n        y_valid = train.loc[valid_idx, 'target']\n        print(X_train.shape, y_train.isnull().sum())\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.DMatrix(data=X_train, label=y_train)\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    \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, X_train, y_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#     auc(oof.target.values, oof.oof_pred.values)\n    return oof","metadata":{"execution":{"iopub.status.busy":"2022-07-28T12:15:26.271972Z","iopub.execute_input":"2022-07-28T12:15:26.272320Z","iopub.status.idle":"2022-07-28T12:15:26.294542Z","shell.execute_reply.started":"2022-07-28T12:15:26.272290Z","shell.execute_reply":"2022-07-28T12:15:26.293604Z"},"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-28T12:15:28.049151Z","iopub.execute_input":"2022-07-28T12:15:28.049810Z","iopub.status.idle":"2022-07-28T12:30:22.912328Z","shell.execute_reply.started":"2022-07-28T12:15:28.049776Z","shell.execute_reply":"2022-07-28T12:30:22.911421Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# CLEAN RAM\ndel train\n_ = gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-07-28T12:30:22.914169Z","iopub.execute_input":"2022-07-28T12:30:22.914831Z","iopub.status.idle":"2022-07-28T12:30:23.029396Z","shell.execute_reply.started":"2022-07-28T12:30:22.914794Z","shell.execute_reply":"2022-07-28T12:30:23.028301Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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":{"execution":{"iopub.status.busy":"2022-07-28T12:30:23.031203Z","iopub.execute_input":"2022-07-28T12:30:23.031994Z","iopub.status.idle":"2022-07-28T12:30:23.187093Z","shell.execute_reply.started":"2022-07-28T12:30:23.031955Z","shell.execute_reply":"2022-07-28T12:30:23.185968Z"},"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":{"execution":{"iopub.status.busy":"2022-07-28T12:30:23.190777Z","iopub.execute_input":"2022-07-28T12:30:23.191167Z","iopub.status.idle":"2022-07-28T12:30:23.496388Z","shell.execute_reply.started":"2022-07-28T12:30:23.191131Z","shell.execute_reply":"2022-07-28T12:30:23.495526Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# CLEAR VRAM, RAM FOR INFERENCE BELOW\ndel oof_xgb, oof\n_ = gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-07-28T12:30:23.497882Z","iopub.execute_input":"2022-07-28T12:30:23.498212Z","iopub.status.idle":"2022-07-28T12:30:23.627742Z","shell.execute_reply.started":"2022-07-28T12:30:23.498179Z","shell.execute_reply":"2022-07-28T12:30:23.626768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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":{"execution":{"iopub.status.busy":"2022-07-28T12:30:23.629056Z","iopub.execute_input":"2022-07-28T12:30:23.629968Z","iopub.status.idle":"2022-07-28T12:30:23.663791Z","shell.execute_reply.started":"2022-07-28T12:30:23.629932Z","shell.execute_reply":"2022-07-28T12:30:23.662965Z"},"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":{"execution":{"iopub.status.busy":"2022-07-28T12:30:23.666403Z","iopub.execute_input":"2022-07-28T12:30:23.667196Z","iopub.status.idle":"2022-07-28T12:30:23.978836Z","shell.execute_reply.started":"2022-07-28T12:30:23.667160Z","shell.execute_reply":"2022-07-28T12:30:23.977842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ngc.collect()\n\n# INFER TEST DATA IN PARTS\n\nTEST_SECTIONS = 2\nTEST_SUB_SECTIONS = 2\n\ntest_preds = []\ncustomers = False\nfor 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":{"execution":{"iopub.status.busy":"2022-07-28T12:30:23.980308Z","iopub.execute_input":"2022-07-28T12:30:23.980896Z","iopub.status.idle":"2022-07-28T12:30:40.343182Z","shell.execute_reply.started":"2022-07-28T12:30:23.980861Z","shell.execute_reply":"2022-07-28T12:30:40.341982Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n# WRITE SUBMISSION FILE\ntest_preds = np.concatenate(test_preds)\ntest = cudf.DataFrame(index=customers,data={'prediction':test_preds})\nsub = cudf.read_csv('../input/boolart-user-default-prediction/sample_submission.csv')[['customer_ID']]\nsub['customer_ID_hash'] = sub['customer_ID'].str[-16:].str.hex_to_int().astype('int64')\nsub = sub.set_index('customer_ID_hash')\nsub = sub.merge(test[['prediction']], left_index=True, right_index=True, how='left')\nsub = sub.reset_index(drop=True)\n\n\n# DISPLAY PREDICTIONS\nsub.to_csv(f'submission.csv',index=False)\nprint('Submission file shape is', sub.shape )\nsub.head()\n\n# PLOT PREDICTIONS\nplt.hist(sub.to_pandas().prediction, bins=100)\nplt.title('Test Predictions')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-28T12:30:40.345012Z","iopub.execute_input":"2022-07-28T12:30:40.345360Z","iopub.status.idle":"2022-07-28T12:30:40.832174Z","shell.execute_reply.started":"2022-07-28T12:30:40.345326Z","shell.execute_reply":"2022-07-28T12:30:40.831270Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}