{"cells":[{"metadata":{"trusted":true,"collapsed":true,"_uuid":"6707fffe3c7a8c42bdb3f906977b6f342259f4c7"},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport lightgbm as lgb\nimport os\n\nfrom sklearn.ensemble import ExtraTreesRegressor, ExtraTreesClassifier, RandomForestClassifier\nfrom sklearn.model_selection import train_test_split\n\nfrom hep_ml.gradientboosting import UGradientBoostingClassifier\nfrom hep_ml.losses import BinFlatnessLossFunction","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"6b66dc605e27d4c6dfdc712578560d4d59581911"},"cell_type":"code","source":"import numpy\nfrom sklearn.metrics import roc_curve, auc\n\n\ndef __rolling_window(data, window_size):\n    \"\"\"\n    Rolling window: take window with definite size through the array\n\n    :param data: array-like\n    :param window_size: size\n    :return: the sequence of windows\n\n    Example: data = array(1, 2, 3, 4, 5, 6), window_size = 4\n        Then this function return array(array(1, 2, 3, 4), array(2, 3, 4, 5), array(3, 4, 5, 6))\n    \"\"\"\n    shape = data.shape[:-1] + (data.shape[-1] - window_size + 1, window_size)\n    strides = data.strides + (data.strides[-1],)\n    return numpy.lib.stride_tricks.as_strided(data, shape=shape, strides=strides)\n\n\ndef __cvm(subindices, total_events):\n    \"\"\"\n    Compute Cramer-von Mises metric.\n    Compared two distributions, where first is subset of second one.\n    Assuming that second is ordered by ascending\n\n    :param subindices: indices of events which will be associated with the first distribution\n    :param total_events: count of events in the second distribution\n    :return: cvm metric\n    \"\"\"\n    target_distribution = numpy.arange(1, total_events + 1, dtype='float') / total_events\n    subarray_distribution = numpy.cumsum(numpy.bincount(subindices, minlength=total_events), dtype='float')\n    subarray_distribution /= 1.0 * subarray_distribution[-1]\n    return numpy.mean((target_distribution - subarray_distribution) ** 2)\n\n\ndef compute_cvm(predictions, masses, n_neighbours=200, step=50):\n    \"\"\"\n    Computing Cramer-von Mises (cvm) metric on background events: take average of cvms calculated for each mass bin.\n    In each mass bin global prediction's cdf is compared to prediction's cdf in mass bin.\n\n    :param predictions: array-like, predictions\n    :param masses: array-like, in case of Kaggle tau23mu this is reconstructed mass\n    :param n_neighbours: count of neighbours for event to define mass bin\n    :param step: step through sorted mass-array to define next center of bin\n    :return: average cvm value\n    \"\"\"\n    predictions = numpy.array(predictions)\n    masses = numpy.array(masses)\n    assert len(predictions) == len(masses)\n\n    # First, reorder by masses\n    predictions = predictions[numpy.argsort(masses)]\n\n    # Second, replace probabilities with order of probability among other events\n    predictions = numpy.argsort(numpy.argsort(predictions, kind='mergesort'), kind='mergesort')\n\n    # Now, each window forms a group, and we can compute contribution of each group to CvM\n    cvms = []\n    for window in __rolling_window(predictions, window_size=n_neighbours)[::step]:\n        cvms.append(__cvm(subindices=window, total_events=len(predictions)))\n    return numpy.mean(cvms)\n\n\ndef __roc_curve_splitted(data_zero, data_one, sample_weights_zero, sample_weights_one):\n    \"\"\"\n    Compute roc curve\n\n    :param data_zero: 0-labeled data\n    :param data_one:  1-labeled data\n    :param sample_weights_zero: weights for 0-labeled data\n    :param sample_weights_one:  weights for 1-labeled data\n    :return: roc curve\n    \"\"\"\n    labels = [0] * len(data_zero) + [1] * len(data_one)\n    weights = numpy.concatenate([sample_weights_zero, sample_weights_one])\n    data_all = numpy.concatenate([data_zero, data_one])\n    fpr, tpr, _ = roc_curve(labels, data_all, sample_weight=weights)\n    return fpr, tpr\n\n\ndef compute_ks(data_prediction, mc_prediction, weights_data, weights_mc):\n    \"\"\"\n    Compute Kolmogorov-Smirnov (ks) distance between real data predictions cdf and Monte Carlo one.\n\n    :param data_prediction: array-like, real data predictions\n    :param mc_prediction: array-like, Monte Carlo data predictions\n    :param weights_data: array-like, real data weights\n    :param weights_mc: array-like, Monte Carlo weights\n    :return: ks value\n    \"\"\"\n    assert len(data_prediction) == len(weights_data), 'Data length and weight one must be the same'\n    assert len(mc_prediction) == len(weights_mc), 'Data length and weight one must be the same'\n\n    data_prediction, mc_prediction = numpy.array(data_prediction), numpy.array(mc_prediction)\n    weights_data, weights_mc = numpy.array(weights_data), numpy.array(weights_mc)\n\n    assert numpy.all(data_prediction >= 0.) and numpy.all(data_prediction <= 1.), 'Data predictions are out of range [0, 1]'\n    assert numpy.all(mc_prediction >= 0.) and numpy.all(mc_prediction <= 1.), 'MC predictions are out of range [0, 1]'\n\n    weights_data /= numpy.sum(weights_data)\n    weights_mc /= numpy.sum(weights_mc)\n\n    fpr, tpr = __roc_curve_splitted(data_prediction, mc_prediction, weights_data, weights_mc)\n\n    Dnm = numpy.max(numpy.abs(fpr - tpr))\n    return Dnm\n\n\ndef roc_auc_truncated(labels, predictions, tpr_thresholds=(0.2, 0.4, 0.6, 0.8),\n                      roc_weights=(4, 3, 2, 1, 0)):\n    \"\"\"\n    Compute weighted area under ROC curve.\n\n    :param labels: array-like, true labels\n    :param predictions: array-like, predictions\n    :param tpr_thresholds: array-like, true positive rate thresholds delimiting the ROC segments\n    :param roc_weights: array-like, weights for true positive rate segments\n    :return: weighted AUC\n    \"\"\"\n    assert numpy.all(predictions >= 0.) and numpy.all(predictions <= 1.), 'Data predictions are out of range [0, 1]'\n    assert len(tpr_thresholds) + 1 == len(roc_weights), 'Incompatible lengths of thresholds and weights'\n    fpr, tpr, _ = roc_curve(labels, predictions)\n    area = 0.\n    tpr_thresholds = [0.] + list(tpr_thresholds) + [1.]\n    for index in range(1, len(tpr_thresholds)):\n        tpr_cut = numpy.minimum(tpr, tpr_thresholds[index])\n        tpr_previous = numpy.minimum(tpr, tpr_thresholds[index - 1])\n        area += roc_weights[index - 1] * (auc(fpr, tpr_cut, reorder=True) - auc(fpr, tpr_previous, reorder=True))\n    tpr_thresholds = numpy.array(tpr_thresholds)\n    # roc auc normalization to be 1 for an ideal classifier\n    area /= numpy.sum((tpr_thresholds[1:] - tpr_thresholds[:-1]) * numpy.array(roc_weights))\n    return area\n\ndef feature_importance(forest, X_train):\n    ranked_list = []\n    \n    importances = forest.feature_importances_\n\n    std = np.std([tree.feature_importances_ for tree in forest.estimators_], axis=0)\n    indices = np.argsort(importances)[::-1]\n\n    # Print the feature ranking\n    print(\"Feature ranking:\")\n\n    for f in range(X_train.shape[1]):\n        print(\"%d. feature %d (%f)\" % (f + 1, indices[f], importances[indices[f]]) + \" - \" + X_train.columns[indices[f]])\n        ranked_list.append(X_train.columns[indices[f]])\n    \n    return ranked_list","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"f2f6f087244bc0f707f9bb91df08210cb8f5cd11"},"cell_type":"code","source":"data_path = \"../input\"\ntrain = pd.read_csv(os.path.join(data_path, 'training.csv'), index_col='id')\ntest = pd.read_csv(os.path.join(data_path, 'test.csv'), index_col='id')\ncheck_agreement = pd.read_csv(os.path.join(data_path, 'check_agreement.csv'), index_col='id')\n\ntrainids = train.index.values\ntestids = test.index.values\ncaids = check_agreement.index.values\ntrainsignals = train.signal.ravel()\nsignal = train.signal","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"e5910d9e9a35fa07c8ebf91cde36ff94d004fd72"},"cell_type":"code","source":"def add_features(df):\n    # features used by the others on Kaggle\n    df['NEW_FD_SUMP']=df['FlightDistance']/(df['p0_p']+df['p1_p']+df['p2_p'])\n    df['NEW5_lt']=df['LifeTime']*(df['p0_IP']+df['p1_IP']+df['p2_IP'])/3\n    df['p_track_Chi2Dof_MAX'] = df.loc[:, ['p0_track_Chi2Dof', 'p1_track_Chi2Dof', 'p2_track_Chi2Dof']].max(axis=1)\n    #df['flight_dist_sig'] = df['FlightDistance']/df['FlightDistanceError'] # modified to:\n    df['flight_dist_sig2'] = (df['FlightDistance']/df['FlightDistanceError'])**2\n    # features from phunter\n    df['flight_dist_sig'] = df['FlightDistance']/df['FlightDistanceError']\n    df['NEW_IP_dira'] = df['IP']*df['dira']\n    df['p0p2_ip_ratio']=df['IP']/df['IP_p0p2']\n    df['p1p2_ip_ratio']=df['IP']/df['IP_p1p2']\n    df['DCA_MAX'] = df.loc[:, ['DOCAone', 'DOCAtwo', 'DOCAthree']].max(axis=1)\n    df['iso_bdt_min'] = df.loc[:, ['p0_IsoBDT', 'p1_IsoBDT', 'p2_IsoBDT']].min(axis=1)\n    df['iso_min'] = df.loc[:, ['isolationa', 'isolationb', 'isolationc','isolationd', 'isolatione', 'isolationf']].min(axis=1)\n    # My:\n    # new combined features just to minimize their number;\n    # their physical sense doesn't matter\n    df['NEW_iso_abc'] = df['isolationa']*df['isolationb']*df['isolationc']\n    df['NEW_iso_def'] = df['isolationd']*df['isolatione']*df['isolationf']\n    df['NEW_pN_IP'] = df['p0_IP']+df['p1_IP']+df['p2_IP']\n    df['NEW_pN_p']  = df['p0_p']+df['p1_p']+df['p2_p']\n    df['NEW_IP_pNpN'] = df['IP_p0p2']*df['IP_p1p2']\n    df['NEW_pN_IPSig'] = df['p0_IPSig']+df['p1_IPSig']+df['p2_IPSig']\n    #My:\n    # \"super\" feature changing the result from 0.988641 to 0.991099\n    df['NEW_FD_LT']=df['FlightDistance']/df['LifeTime']\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"a35c0b238f4e152958ffa4491f2e5778b9f9e5e6"},"cell_type":"code","source":"train = add_features(train)\ntest = add_features(test)\ncheck_agreement = add_features(check_agreement)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a59bff8e34bae91fbfd827eed2f9de9e0af19d72"},"cell_type":"markdown","source":"## Five Lines Model"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"01bd7494d33a52fb0e5b2e7879ab478c439942b0"},"cell_type":"code","source":"p1 = 11.05855369567871094\np2 = 0.318310\np3 = 1.570796\n\ndef Output(p):\n    return 1/(1.+np.exp(-p))\n\ndef GP(data):\n    return Output(  1.0*np.tanh(((((((((data[\"IPSig\"]) + (data[\"ISO_SumBDT\"]))) - (np.minimum(((-2.0)), ((data[\"ISO_SumBDT\"])))))) / (data[\"ISO_SumBDT\"]))) / (np.minimum((((-1.0*((data[\"ISO_SumBDT\"]))))), ((data[\"IPSig\"])))))) +\n                    1.0*np.tanh((-1.0*((((data[\"iso\"]) + (((((((((((((data[\"VertexChi2\"]) + ((3.0)))) / (data[\"ISO_SumBDT\"]))) * (data[\"IP\"]))) * 2.0)) / (data[\"ISO_SumBDT\"]))) * (((((((((((data[\"VertexChi2\"]) + ((3.0)))) / (data[\"ISO_SumBDT\"]))) * (data[\"IP\"]))) * 2.0)) / (data[\"ISO_SumBDT\"])))))))))) +\n                    1.0*np.tanh((-1.0*(((((((((data[\"IPSig\"]) * ((((data[\"iso\"]) + (((data[\"IP\"]) * 2.0)))/2.0)))) + (np.tanh((data[\"p0_IsoBDT\"]))))/2.0)) * ((((data[\"p0_IsoBDT\"]) + (data[\"IPSig\"]))/2.0))))))) +\n                    1.0*np.tanh(((np.minimum(((np.cos((((np.cos((((data[\"p0_track_Chi2Dof\"]) * (np.cos((data[\"p0_track_Chi2Dof\"]))))))) * (np.log((data[\"IP_p0p2\"])))))))), ((np.cos((data[\"p0_track_Chi2Dof\"])))))) * (data[\"p0_track_Chi2Dof\"]))) +\n                    1.0*np.tanh((((((((((p1)) / (((((p1)) + (((((data[\"SPDhits\"]) / 2.0)) / 2.0)))/2.0)))) - (data[\"IP\"]))) - (((data[\"SPDhits\"]) / (data[\"p1_pt\"]))))) * 2.0)) +\n                    1.0*np.tanh((((((((((((((data[\"CDF3\"]) / (data[\"dira\"]))) > (data[\"CDF3\"]))*1.)) > (data[\"CDF3\"]))*1.)) / 2.0)) + ((-1.0*((((((data[\"CDF3\"]) * (data[\"p2_track_Chi2Dof\"]))) * (((data[\"CDF3\"]) * (data[\"p2_track_Chi2Dof\"])))))))))/2.0)) +\n                    1.0*np.tanh((((-1.0*((((data[\"DOCAthree\"]) / (data[\"CDF2\"])))))) + (np.minimum(((((data[\"p2_pt\"]) / (data[\"p0_p\"])))), ((np.minimum(((data[\"CDF2\"])), ((((np.sin((p3))) / 2.0)))))))))) +\n                    1.0*np.tanh(np.minimum((((-1.0*(((((((data[\"FlightDistance\"]) < (data[\"IPSig\"]))*1.)) / 2.0)))))), ((((np.minimum(((np.cos((np.log((data[\"p0_pt\"])))))), ((np.cos((data[\"p1_track_Chi2Dof\"])))))) / (p2)))))) +\n                    1.0*np.tanh(((np.sin((np.where(data[\"iso\"]>0, ((((data[\"iso\"]) - ((-1.0*((((data[\"IPSig\"]) / 2.0))))))) / 2.0), ((((3.0) * (data[\"IP\"]))) * 2.0) )))) / 2.0)) +\n                    1.0*np.tanh(((((np.cos(((((data[\"ISO_SumBDT\"]) + (p2))/2.0)))) - (np.sin((np.log((data[\"p1_eta\"]))))))) - ((((((data[\"ISO_SumBDT\"]) + (np.cos((data[\"p2_IsoBDT\"]))))/2.0)) * ((((data[\"ISO_SumBDT\"]) + (np.cos((data[\"p2_IsoBDT\"]))))/2.0)))))))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"fd4221664ef29a8841606bffceab9fc7ab884a0d"},"cell_type":"code","source":"tr_preds_1 = GP(train).values\ntest_preds_1 = GP(test).values\nca_preds_1 = GP(check_agreement).values\n\ntest_predictions = pd.DataFrame({'id':testids,'predictions_1':test_preds_1})\ntrain_predictions_all = pd.DataFrame({'id':trainids,'predictions_1':tr_preds_1})\nca_predictions = pd.DataFrame({'id':caids,'predictions_1':ca_preds_1})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"1e03455e724c7786865e9b3b24046ca76e700d3b"},"cell_type":"code","source":"# since the target is not used for this model we can add the feature to our data without any leakage\ntrain['lines'] = tr_preds_1\ncheck_agreement['lines'] = ca_preds_1\ntest['lines'] = test_preds_1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"627684ba2e115d589db593a8076addce874da233","collapsed":true},"cell_type":"code","source":"agreement_probs = ca_predictions.predictions_1\n\nks = compute_ks(\n    agreement_probs[check_agreement['signal'].values == 0],\n    agreement_probs[check_agreement['signal'].values == 1],\n    check_agreement[check_agreement['signal'] == 0]['weight'].values,\n    check_agreement[check_agreement['signal'] == 1]['weight'].values)\n\nprint('KS metric', ks, ks < 0.09)\n# print(roc_auc_truncated(y_cv, cv_predictions.predictions_1))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2e35733acf8fe23490ee0f3b1badddbf6ccfa3b3"},"cell_type":"markdown","source":"### Split Data"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"4ebbcee9ca7bf25da3d4ad581673316e70fbfa27"},"cell_type":"code","source":"# split data into train and cv\nX_tr, X_cv, y_tr, y_cv, X_tr_id, X_cv_id, train_predictions, cv_predictions = train_test_split(train, signal, trainids, train_predictions_all, random_state=100, test_size=0.25, shuffle=True)\n\n# copy our predictions so they are not slices and we won't get errors\n# train_predictions = train_predictions.copy()\ncv_predictions = cv_predictions.copy()\nX_cv = X_cv.copy()\n\n# train on whole data set now\ntrain_predictions = train_predictions_all.copy()\nX_tr = train.copy()\ny_tr = signal\nX_tr = X_tr.copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"18d0c0fa8a8752061d6b854e46185eb057a0fa57"},"cell_type":"code","source":"def add_lines(data):\n    data['line1'] = 1.0*np.tanh(((((((((data[\"IPSig\"]) + (data[\"ISO_SumBDT\"]))) - (np.minimum(((-2.0)), ((data[\"ISO_SumBDT\"])))))) / (data[\"ISO_SumBDT\"]))) / (np.minimum((((-1.0*((data[\"ISO_SumBDT\"]))))), ((data[\"IPSig\"]))))))\n    data['line2'] = 1.0*np.tanh((-1.0*((((data[\"iso\"]) + (((((((((((((data[\"VertexChi2\"]) + ((3.0)))) / (data[\"ISO_SumBDT\"]))) * (data[\"IP\"]))) * 2.0)) / (data[\"ISO_SumBDT\"]))) * (((((((((((data[\"VertexChi2\"]) + ((3.0)))) / (data[\"ISO_SumBDT\"]))) * (data[\"IP\"]))) * 2.0)) / (data[\"ISO_SumBDT\"]))))))))))\n    data['line3'] = 1.0*np.tanh((-1.0*(((((((((data[\"IPSig\"]) * ((((data[\"iso\"]) + (((data[\"IP\"]) * 2.0)))/2.0)))) + (np.tanh((data[\"p0_IsoBDT\"]))))/2.0)) * ((((data[\"p0_IsoBDT\"]) + (data[\"IPSig\"]))/2.0)))))))\n    data['line4'] = 1.0*np.tanh(((np.minimum(((np.cos((((np.cos((((data[\"p0_track_Chi2Dof\"]) * (np.cos((data[\"p0_track_Chi2Dof\"]))))))) * (np.log((data[\"IP_p0p2\"])))))))), ((np.cos((data[\"p0_track_Chi2Dof\"])))))) * (data[\"p0_track_Chi2Dof\"])))\n    data['line5'] = 1.0*np.tanh((((((((((p1)) / (((((p1)) + (((((data[\"SPDhits\"]) / 2.0)) / 2.0)))/2.0)))) - (data[\"IP\"]))) - (((data[\"SPDhits\"]) / (data[\"p1_pt\"]))))) * 2.0))\n    data['line6'] = 1.0*np.tanh((((((((((((((data[\"CDF3\"]) / (data[\"dira\"]))) > (data[\"CDF3\"]))*1.)) > (data[\"CDF3\"]))*1.)) / 2.0)) + ((-1.0*((((((data[\"CDF3\"]) * (data[\"p2_track_Chi2Dof\"]))) * (((data[\"CDF3\"]) * (data[\"p2_track_Chi2Dof\"])))))))))/2.0))\n    data['line7'] = 1.0*np.tanh((((-1.0*((((data[\"DOCAthree\"]) / (data[\"CDF2\"])))))) + (np.minimum(((((data[\"p2_pt\"]) / (data[\"p0_p\"])))), ((np.minimum(((data[\"CDF2\"])), ((((np.sin((p3))) / 2.0))))))))))\n    data['line8'] = 1.0*np.tanh(np.minimum((((-1.0*(((((((data[\"FlightDistance\"]) < (data[\"IPSig\"]))*1.)) / 2.0)))))), ((((np.minimum(((np.cos((np.log((data[\"p0_pt\"])))))), ((np.cos((data[\"p1_track_Chi2Dof\"])))))) / (p2))))))\n    data['line9'] = 1.0*np.tanh(((np.sin((np.where(data[\"iso\"]>0, ((((data[\"iso\"]) - ((-1.0*((((data[\"IPSig\"]) / 2.0))))))) / 2.0), ((((3.0) * (data[\"IP\"]))) * 2.0) )))) / 2.0))\n    data['line10'] = 1.0*np.tanh(((((np.cos(((((data[\"ISO_SumBDT\"]) + (p2))/2.0)))) - (np.sin((np.log((data[\"p1_eta\"]))))))) - ((((((data[\"ISO_SumBDT\"]) + (np.cos((data[\"p2_IsoBDT\"]))))/2.0)) * ((((data[\"ISO_SumBDT\"]) + (np.cos((data[\"p2_IsoBDT\"]))))/2.0))))))\n    \n    return data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"dabcec1aa901f20ada118fc9b76116f1eb40cd10"},"cell_type":"code","source":"X_tr = add_lines(X_tr)\nX_cv = add_lines(X_cv)\ntest = add_lines(test)\ntrain = add_lines(train)\ncheck_agreement = add_lines(check_agreement)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d1613b2a69d7123a644af9403b63d25f16d3b4c3"},"cell_type":"markdown","source":"### LightGBM"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"cc1a5d6dbef010741bb1c25785536acd8048b512"},"cell_type":"code","source":"feature_names = ['LifeTime',\n 'dira',\n 'FlightDistance',\n 'FlightDistanceError',\n 'IP',\n 'IPSig',\n 'VertexChi2',\n 'pt',\n 'iso',\n 'ISO_SumBDT',\n 'NEW_FD_SUMP',\n 'NEW5_lt',\n 'p_track_Chi2Dof_MAX',\n 'flight_dist_sig2',\n 'flight_dist_sig',\n 'NEW_IP_dira',\n 'p0p2_ip_ratio',\n 'p1p2_ip_ratio',\n 'DCA_MAX',\n 'iso_bdt_min',\n 'iso_min',\n 'NEW_iso_abc',\n 'NEW_iso_def',\n 'NEW_pN_IP',\n 'NEW_pN_p',\n 'NEW_IP_pNpN',\n 'NEW_pN_IPSig',\n 'NEW_FD_LT',\n                \n                'line1', 'line2', 'line3', 'line4',\n                 'line6',\n                 'line7',\n                 'line9',\n                 'line10',\n                 'line8',\n                ]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0c37e744e572e1d52d9494839b128fa112155d95","collapsed":true},"cell_type":"code","source":"# use the full training set for our cv, and then train only on training set so we can validate\n# train_all = lgb.Dataset(train[feature_names],signal)\ntrain_all = train_set = lgb.Dataset(X_tr[feature_names],y_tr)\n\nparams = {\n    'task': 'train',\n    'boosting_type': 'gbdt',\n    'objective': 'binary',\n    'num_leaves': 2**8,\n    'metric': {'auc'},\n    'min_data_in_leaf': 31,\n    'max_depth': 12,\n    'learning_rate': 0.05,\n    'bagging_fraction': 0.5,\n    'lambda': 0.1,\n    'feature_fraction': 0.5,\n}\n\ncv_output = lgb.cv(\n    params,\n    train_all,\n    num_boost_round=450,\n    nfold=5,\n)\n\nbest_niter = np.argmax(cv_output['auc-mean'])\nbest_score = cv_output['auc-mean'][best_niter]\nprint('Best number of iterations: {}'.format(best_niter))\nprint('Best CV score: {}'.format(best_score))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"818ebd0cec79490fcab9ca3348487ae975d5b7a4"},"cell_type":"code","source":"model = lgb.train(params, train_set, num_boost_round=best_niter)\n\ntrain_predictions['predictions_2'] = model.predict(X_tr[feature_names])\ncv_predictions['predictions_2'] = model.predict(X_cv[feature_names])\ntest_predictions['predictions_2'] = model.predict(test[feature_names])\nca_predictions['predictions_2'] = model.predict(check_agreement[feature_names])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"96356ba552bd96474817af64fdac64de7d5b82e0","collapsed":true},"cell_type":"code","source":"agreement_probs = ca_predictions.predictions_2.values\n\nks = compute_ks(\n    agreement_probs[check_agreement['signal'].values == 0],\n    agreement_probs[check_agreement['signal'].values == 1],\n    check_agreement[check_agreement['signal'] == 0]['weight'].values,\n    check_agreement[check_agreement['signal'] == 1]['weight'].values)\n\nprint('KS metric', ks, ks < 0.09)\nprint(roc_auc_truncated(y_cv, cv_predictions['predictions_2']))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"4711751ec0fb32276c580bc84ebe757d4fe807de","collapsed":true},"cell_type":"code","source":"# depth 12\nprint('KS metric', ks, ks < 0.09)\nprint(roc_auc_truncated(y_cv, cv_predictions['predictions_2']))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c29d76e2d45b3ef65b4175480718c781bb0f185c"},"cell_type":"markdown","source":"## ExtraTrees"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"f3b25cf5cb90edf9b5f1a65a86de3f7f678dbd3d"},"cell_type":"code","source":"et_features = ['LifeTime', 'dira', 'FlightDistance', 'FlightDistanceError', 'IP',\n       'VertexChi2', 'pt', 'DOCAone', 'DOCAtwo', 'DOCAthree',\n       'IP_p0p2', 'IP_p1p2', 'isolationa', 'isolationb', 'isolationc',\n       'isolationd', 'isolatione', 'isolationf', 'iso',\n       'ISO_SumBDT', 'p0_IsoBDT', 'p1_IsoBDT', 'p2_IsoBDT', 'p0_track_Chi2Dof',\n       'p1_track_Chi2Dof', 'p2_track_Chi2Dof', 'p0_IP', 'p1_IP', 'p2_IP',\n       'p0_IPSig', 'p1_IPSig', 'p2_IPSig', 'p0_pt', 'p1_pt', 'p2_pt', 'p0_p',\n       'p1_p', 'p2_p', 'p0_eta', 'p1_eta', 'p2_eta', 'NEW_FD_SUMP',\n       'NEW5_lt', 'p_track_Chi2Dof_MAX', 'flight_dist_sig2', 'flight_dist_sig',\n       'NEW_IP_dira', 'p0p2_ip_ratio', 'p1p2_ip_ratio', 'DCA_MAX',\n       'iso_bdt_min', 'iso_min', 'NEW_iso_abc', 'NEW_iso_def', 'NEW_pN_IP',\n       'NEW_pN_p', 'NEW_IP_pNpN', 'NEW_pN_IPSig', 'NEW_FD_LT', \n       'lines', 'line1', 'line2', 'line3', 'line4',  'line6', 'line7', # 'line5',\n       'line8', 'line9', 'line10']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"41ef040f497e537baa8ca51464fed25f9f593b42","collapsed":true},"cell_type":"code","source":"et = ExtraTreesClassifier(n_estimators=100, random_state=0, max_depth=22, min_impurity_decrease=1e-8, min_samples_leaf=15, n_jobs=-1, verbose=1)\net.fit(X_tr[et_features], y_tr)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"20bb2904f1521d090cf5ab5cbfdee5022663c2b8","collapsed":true},"cell_type":"code","source":"tr_predictions_3 = et.predict_proba(X_tr[et_features])[:,1]\ncv_predictions_3 = et.predict_proba(X_cv[et_features])[:,1]\ntest_predictions_3 = et.predict_proba(test[et_features])[:,1]\nagreement_predictions_3 = et.predict_proba(check_agreement[et_features])[:,1]\n\nprint(\"Train Max:\", np.max(tr_predictions_3))\nprint(\"Test Max:\", np.max(test_predictions_3))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6b7441c25e7f97354c23b7a27df1a680eab4168d","collapsed":true},"cell_type":"code","source":"agreement_probs = agreement_predictions_3\n\nks = compute_ks(\n    agreement_probs[check_agreement['signal'].values == 0],\n    agreement_probs[check_agreement['signal'].values == 1],\n    check_agreement[check_agreement['signal'] == 0]['weight'].values,\n    check_agreement[check_agreement['signal'] == 1]['weight'].values)\n\nnoise = np.random.normal(0,0.01,len(cv_predictions_3))\n\nprint('KS metric', ks, ks < 0.09)\nprint(roc_auc_truncated(y_cv, cv_predictions_3))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ab61c461a91e395fa8a10d504dbe989ae78aeec4","collapsed":true},"cell_type":"code","source":"print('KS metric', ks, ks < 0.09)\nprint(roc_auc_truncated(y_cv, cv_predictions_3))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"eb7ad33ad68c43dfa14d6c3a3ee8b77f45b5b0f0"},"cell_type":"code","source":"train_predictions['predictions_3'] = tr_predictions_3\ncv_predictions['predictions_3'] = cv_predictions_3\ntest_predictions['predictions_3'] = test_predictions_3\nca_predictions['predictions_3'] = agreement_predictions_3","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c5f0e37b11fb82926ad46c9fde61e30d7858e498","collapsed":true},"cell_type":"code","source":"ranked = feature_importance(et, test[et_features])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a2404631fe10bc9213bdf53b4d3e10093278264c"},"cell_type":"markdown","source":"### RandomForest"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"b323b0f067cfed33eeba159a1f50587c7bff00db"},"cell_type":"code","source":"rf_features = ['lines',\n 'line3',\n 'IPSig',\n 'p0p2_ip_ratio',\n 'line2',\n 'IP',\n 'line8',\n 'p_track_Chi2Dof_MAX',\n 'dira',\n 'p1p2_ip_ratio',\n 'VertexChi2',\n 'line9',\n 'p0_track_Chi2Dof',\n 'iso_bdt_min',\n 'ISO_SumBDT',\n 'NEW_FD_SUMP',\n 'DCA_MAX',\n 'p0_IP',\n 'p0_IPSig',\n 'flight_dist_sig2',\n 'NEW_IP_dira',\n 'LifeTime',\n 'flight_dist_sig',\n 'NEW_pN_IPSig',\n 'p0_IsoBDT',\n 'NEW_pN_IP',\n 'line1',\n 'line4',\n 'NEW_pN_p',\n 'p2_IPSig',\n 'p1_track_Chi2Dof',\n 'p1_IsoBDT',\n 'NEW_FD_LT',\n 'IP_p1p2',\n 'pt',\n 'NEW5_lt',\n 'line7',\n 'iso',\n 'p2_track_Chi2Dof',\n 'IP_p0p2',\n 'p1_p',\n 'p1_IPSig',\n 'p2_IsoBDT',\n 'p0_p',\n 'p1_eta',\n 'line10',\n 'p2_IP',\n 'NEW_IP_pNpN',\n 'DOCAone',\n 'p0_pt',\n 'FlightDistance',\n 'DOCAthree',\n 'p1_pt',\n 'p0_eta',\n 'p1_IP',\n 'FlightDistanceError',\n 'DOCAtwo',\n 'p2_pt',\n 'p2_eta',\n 'p2_p',\n 'CDF1',\n 'CDF2']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"43583c174b50c1f8d25228416ac46fdd810737c1","collapsed":true},"cell_type":"code","source":"rf = RandomForestClassifier(n_estimators=200, random_state=0, max_depth=15, min_impurity_decrease=1e-6, min_samples_leaf=20, n_jobs=-1, verbose=1)\nrf.fit(X_tr[rf_features], y_tr)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"aed6a12ba2753097cdc06d31ec00e7e21d53ca65","collapsed":true},"cell_type":"code","source":"tr_predictions_4 = rf.predict_proba(X_tr[rf_features])[:,1]\ncv_predictions_4 = rf.predict_proba(X_cv[rf_features])[:,1]\ntest_predictions_4 = rf.predict_proba(test[rf_features])[:,1]\nagreement_predictions_4 = rf.predict_proba(check_agreement[rf_features])[:,1]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"916e2f21c1905f7e1677aa34c69e3ad313870d9d","collapsed":true},"cell_type":"code","source":"agreement_probs = agreement_predictions_4\n\nks = compute_ks(\n    agreement_probs[check_agreement['signal'].values == 0],\n    agreement_probs[check_agreement['signal'].values == 1],\n    check_agreement[check_agreement['signal'] == 0]['weight'].values,\n    check_agreement[check_agreement['signal'] == 1]['weight'].values)\n\nprint('KS metric', ks, ks < 0.09)\nprint(roc_auc_truncated(y_cv, cv_predictions_4))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"765e9e26401d7e508cbbc8a441676d300b6a73cc","collapsed":true},"cell_type":"code","source":"print('KS metric', ks, ks < 0.09)\nprint(roc_auc_truncated(y_cv, cv_predictions_4))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"e0e23420214afa8497cf65c8e6a3049ed2d4dc2c"},"cell_type":"code","source":"train_predictions['predictions_4'] = tr_predictions_4\ncv_predictions['predictions_4'] = cv_predictions_4\ntest_predictions['predictions_4'] = test_predictions_4\nca_predictions['predictions_4'] = agreement_predictions_4","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ff2303a2c497b8f6d21b6a48c89d7f7dd041e4e1","collapsed":true},"cell_type":"code","source":"ranked = feature_importance(rf, test[rf_features])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"dba1fbb964787abd4cd5653b0ce87d24c3ebbaa3"},"cell_type":"markdown","source":"### UGBC"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"46f3687376cbd0b323a8735a75137b4d354e9ba3"},"cell_type":"code","source":"ugbc_features = ['LifeTime',\n 'dira',\n 'FlightDistance',\n 'FlightDistanceError',\n 'IP',\n 'IPSig',\n 'VertexChi2',\n 'pt',\n 'iso',\n 'ISO_SumBDT',\n 'NEW_FD_SUMP',\n 'NEW5_lt',\n 'p_track_Chi2Dof_MAX',\n 'flight_dist_sig2',\n 'flight_dist_sig',\n 'NEW_IP_dira',\n 'p0p2_ip_ratio',\n 'p1p2_ip_ratio',\n 'DCA_MAX',\n 'iso_bdt_min',\n 'iso_min',\n 'NEW_iso_abc',\n 'NEW_iso_def',\n 'NEW_pN_IP',\n 'NEW_pN_p',\n 'NEW_IP_pNpN',\n 'NEW_pN_IPSig',\n 'NEW_FD_LT', \n \n 'lines',\n 'line3',\n]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b432b67ea88d972a5d07611d00ed2494199abb52","collapsed":true},"cell_type":"code","source":"loss = BinFlatnessLossFunction(['mass'], n_bins=15, uniform_label=0 , fl_coefficient=15, power=2)\nugbc = UGradientBoostingClassifier(loss=loss, n_estimators=250,\n                                 max_depth=8,\n                                 learning_rate=0.15,\n                                 train_features=ugbc_features,\n                                 subsample=0.7,\n                                 random_state=123)\n\nugbc.fit(X_tr[ugbc_features + ['mass']], y_tr)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"cc94d1cea00ff7b9777b221481f0f31f8e00b6dd"},"cell_type":"code","source":"tr_predictions_5 = ugbc.predict_proba(X_tr[ugbc_features])[:,1]\ncv_predictions_5 = ugbc.predict_proba(X_cv[ugbc_features])[:,1]\ntest_predictions_5 = ugbc.predict_proba(test[ugbc_features])[:,1]\nagreement_predictions_5 = ugbc.predict_proba(check_agreement[ugbc_features])[:,1]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"efc92f0cb8fb41d94e5025eeaf4f4d7289419778","collapsed":true},"cell_type":"code","source":"agreement_probs = agreement_predictions_5\n\nks = compute_ks(\n    agreement_probs[check_agreement['signal'].values == 0],\n    agreement_probs[check_agreement['signal'].values == 1],\n    check_agreement[check_agreement['signal'] == 0]['weight'].values,\n    check_agreement[check_agreement['signal'] == 1]['weight'].values)\n\nprint('KS metric', ks, ks < 0.09)\nprint(roc_auc_truncated(y_cv, cv_predictions_5))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b85a4988329500fd9659f7c081c759d0b8d1bab5","collapsed":true},"cell_type":"code","source":"# without lines\nprint('KS metric', ks, ks < 0.09)\nprint(roc_auc_truncated(y_cv, cv_predictions_5))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"4b51c6ddce93015fc2a90be37aef1400867a7af3"},"cell_type":"code","source":"train_predictions['predictions_5'] = tr_predictions_5\ncv_predictions['predictions_5'] = cv_predictions_5\ntest_predictions['predictions_5'] = test_predictions_5\nca_predictions['predictions_5'] = agreement_predictions_5","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"8821a33b1949fef426e0543b7bec7ba4f156dd54"},"cell_type":"code","source":"test_predictions[['id', 'predictions_5']].to_csv(\"20180720_ugbc_1.csv\", index=False, header=[\"id\", \"prediction\"])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"86393ade4056b0a5a6033616874d72674b16eee1"},"cell_type":"markdown","source":"### Average Our Predictions"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"11677ba98dd83c263c09717ea981e3dc316df681"},"cell_type":"code","source":"avg_drop_cols = [\"id\"]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"01efad3eb3b858cd9cb7a44d39c1a39ddaec1a33"},"cell_type":"code","source":"test_predictions['avg_preds'] = test_predictions.drop(avg_drop_cols, axis=1).mean(axis=1)\nca_predictions['avg_preds'] = ca_predictions.drop(avg_drop_cols, axis=1).mean(axis=1)\ncv_predictions['avg_preds'] = cv_predictions.drop(avg_drop_cols, axis=1).mean(axis=1)\ntrain_predictions['avg_preds'] = train_predictions.drop(avg_drop_cols, axis=1).mean(axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"aa531895bbe63b2703f5b69c5a62e0a84ecd4ca9","collapsed":true},"cell_type":"code","source":"agreement_probs = ca_predictions['avg_preds']\n\nks = compute_ks(\n    agreement_probs[check_agreement['signal'].values == 0],\n    agreement_probs[check_agreement['signal'].values == 1],\n    check_agreement[check_agreement['signal'] == 0]['weight'].values,\n    check_agreement[check_agreement['signal'] == 1]['weight'].values)\n\nprint('KS metric', ks, ks < 0.09)\nprint(roc_auc_truncated(y_cv, cv_predictions['avg_preds']))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9063681465e291477e7939dc2e2d429b540b7b57","collapsed":true},"cell_type":"code","source":"print('KS metric', ks, ks < 0.09)\nprint(roc_auc_truncated(y_cv, cv_predictions['avg_preds']))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"collapsed":true,"_uuid":"09f76e0cf7e7de6545d467ff3e19c97e7e01084d"},"cell_type":"code","source":"test_predictions[['id', 'avg_preds']].to_csv(\"20180720_averaged_1.csv\", index=False, header=[\"id\", \"prediction\"])","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"collapsed":true,"_uuid":"ca020fdc6b8834c0cba983fc76b617eb885d53c9"},"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}