{"cells":[{"metadata":{"_uuid":"50dc3fd18d9644aaabfd59dfdf468b8c8b23a5f7"},"cell_type":"markdown","source":"Based on [UGBC GS](https://www.kaggle.com/sionek/ugbc-gs) and the Coursera course  [Addressing Large Hadron Collider Challenges by Machine Learning](https://www.coursera.org/learn/hadron-collider-machine-learning/home/welcome)  \n...and whatever else I borrow from other kernels"},{"metadata":{"trusted":true,"_uuid":"31c6faf233559970c671181c16a11226bcbd00f3","scrolled":true},"cell_type":"code","source":"%pylab inline","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6d9704060ba941a8c7e4baf4a7ba8c7e4f994c41"},"cell_type":"markdown","source":"**Imports**"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import warnings\nwarnings.filterwarnings(\"ignore\", category=DeprecationWarning)\nimport pandas as pd\nfrom sklearn.model_selection import train_test_split\nfrom hep_ml.gradientboosting import UGradientBoostingClassifier\nfrom hep_ml.losses import BinFlatnessLossFunction\nfrom sklearn.metrics import roc_curve, roc_auc_score, auc\nfrom hep_ml import metrics\nfrom sklearn.utils.validation import column_or_1d\nprint(\"Imports added...\")","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"## Load dataset and split into training / test"},{"metadata":{"trusted":true,"_uuid":"bc9b73d86f3dae5d9235b092287375f8345449db","scrolled":true},"cell_type":"code","source":"# -------------- loading data files -------------- #\nprint(\"Load the train/test/eval data using pandas\")\ntrain_ugbc = pd.read_csv(\"../input/training.csv\")\ntrain_ugbc = train_ugbc[train_ugbc['min_ANNmuon'] > 0.4]\ntest_ugbc  = pd.read_csv(\"../input/test.csv\")\ncheck_agreement = pd.read_csv('../input/check_agreement.csv', index_col='id')\n\ntrainids = train_ugbc.index.values\ntestids = test_ugbc.index.values\ncaids = check_agreement.index.values\ntrainsignals = train_ugbc.signal.ravel()\nsignal = train_ugbc.signal\nprint(\"Data loaded...\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"fe8cffe9f0e010af7a3e838aaa2dc047164a6bf9"},"cell_type":"code","source":"# control switches\nDO_5_LINES = True    # Do 5 lines model?\nDO_5_ENS   = False     # Ensemble 5 lines model\nDO_IMP = False       # Do feature permutation importances?\nDO_GRAMOLIN_IMP = False   # Basis for my current 'best'. Use this as ensemble for ensemble importances.\nDO_MASS_PLOT = True # Plot mass correlation?\nDO_MASS_CORR = True  # Check mass correlation?\nDO_NOISE = True     # Add noise to improve monte carlo vs. real. My rough understanding is this means you are training more to \n                     # predicting monte carlo rather than signal? \nDO_GRAMOLIN = True   # Include gramolin solution in final ensemble","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a0d67dc5e526bab0cb7c3ca4ebe75a431ca4e5ec"},"cell_type":"code","source":"# Constants\nMC = 0.002        # Competition allowed mass correlation","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"74c2d1cc9c635fb5b50ced88af007a230a144182"},"cell_type":"markdown","source":"## MonteCarlo vs Real difference"},{"metadata":{"_uuid":"35eb3a2cbe43930125d6de62b2aea725045c3a22"},"cell_type":"markdown","source":"**get_ks_metric**"},{"metadata":{"trusted":true,"_uuid":"b208ffe400d579fb0c4851ec286bb47915045bd5","_kg_hide-input":true},"cell_type":"code","source":"def get_ks_metric(df_agree, df_test):\n    sig_ind = df_agree[df_agree['signal'] == 1].index\n    bck_ind = df_agree[df_agree['signal'] == 0].index\n\n    mc_prob = numpy.array(df_test.loc[sig_ind]['prediction'])\n    mc_weight = numpy.array(df_agree.loc[sig_ind]['weight'])\n    data_prob = numpy.array(df_test.loc[bck_ind]['prediction'])\n    data_weight = numpy.array(df_agree.loc[bck_ind]['weight'])\n    val, agreement_metric = check_agreement_ks_sample_weighted(data_prob, mc_prob, data_weight, mc_weight)\n    return agreement_metric['ks']","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"00495f52651addba1a9b377b56a3001f4aede9af"},"cell_type":"markdown","source":"**check_agreement_ks_sample_weighted (code)**"},{"metadata":{"trusted":true,"_uuid":"a60feb55f3686e02a4984993e090dc50edf1c648","_kg_hide-input":true},"cell_type":"code","source":"def check_agreement_ks_sample_weighted (data_prediction, mc_prediction, weights_data, weights_mc):\n    data_prediction, weights_data = map(column_or_1d, [data_prediction, weights_data])\n    mc_prediction, weights_mc = map(column_or_1d, [mc_prediction, weights_mc])\n\n    assert numpy.all(data_prediction >= 0.) and numpy.all(data_prediction <= 1.), 'error in prediction'\n    assert numpy.all(mc_prediction >= 0.) and numpy.all(mc_prediction <= 1.), 'error in prediction'\n\n    weights_data = weights_data / numpy.sum(weights_data)\n    weights_mc = weights_mc / numpy.sum(weights_mc)\n\n    data_neg = data_prediction[weights_data < 0]\n    weights_neg = -weights_data[weights_data < 0]\n    mc_prediction = numpy.concatenate((mc_prediction, data_neg))\n    weights_mc = numpy.concatenate((weights_mc, weights_neg))\n    data_prediction = data_prediction[weights_data >= 0]\n    weights_data = weights_data[weights_data >= 0]\n\n    assert numpy.all(weights_data >= 0) and numpy.all(weights_mc >= 0)\n    assert numpy.allclose(weights_data.sum(), weights_mc.sum())\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    Dnm_part = numpy.max(numpy.abs(fpr - tpr)[fpr + tpr < 1])\n\n    result = {'ks': Dnm, 'ks_part': Dnm_part}\n    return Dnm_part < 0.03, result","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"869215aa528daa09d13a7bcaf650b6b58f0558c0"},"cell_type":"markdown","source":"**__roc_curve_splitted (code)**"},{"metadata":{"trusted":true,"_uuid":"f5fe5e5b6e40b50ba99d0d39124937f3ef9c1397","_kg_hide-input":false},"cell_type":"code","source":"def __roc_curve_splitted(data_zero, data_one, sample_weights_zero, sample_weights_one):\n    \"\"\"\n    Compute roc curve\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","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1f9b49aaa67b7bf431da603354640bc844c9aaa5"},"cell_type":"markdown","source":"**compute_ks (code)**"},{"metadata":{"trusted":true,"_uuid":"ec1e337713e6c8a2ff32ab63b87e37788d0437fe","_kg_hide-input":true},"cell_type":"code","source":"def 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    :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","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"05bc11d6ddfd40145580090b0572daa975db9501"},"cell_type":"markdown","source":"**roc_auc_truncated **(code)    From starter kit evaluation.py"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"021bb228e40bdc445be07089f1c381dc8b857da6"},"cell_type":"code","source":"def 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","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"503fba8a179e38d8dd355ad397264f1d3fdecbdf"},"cell_type":"markdown","source":"**check_correlation (code)**   For checking mass correlation with solution"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"894816c745f8e3978ff44149e99da82efe04ea97"},"cell_type":"code","source":"def check_correlation(probabilities, mass):\n    probabilities, mass = map(column_or_1d, [probabilities, mass])\n\n    y_pred = numpy.zeros(shape=(len(probabilities), 2))\n    y_pred[:, 1] = probabilities\n    y_pred[:, 0] = 1 - probabilities\n    y_true = [0] * len(probabilities)\n    df_mass = pd.DataFrame({'mass': mass})\n    cvm = metrics.BinBasedCvM(uniform_features=['mass'], uniform_label=0)\n    cvm.fit(df_mass, y_true)\n    return cvm(y_true, y_pred, sample_weight=None)\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7c20011107bb0579e0539c830a4d4bc1571eff5f"},"cell_type":"markdown","source":"**Get the evaluation data**"},{"metadata":{"trusted":true,"_uuid":"1739a848168b183675a0c4ca066b584f5ad203e7"},"cell_type":"code","source":"df_agreement = pd.read_csv('../input/check_agreement.csv')\ndf_corr_check = pd.read_csv(\"../input/check_correlation.csv\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b63c997eca50cf8d91cd6b825ea8d87d08ecc666"},"cell_type":"markdown","source":"**Gramolin 2nd place features from prior competition**  \n[Second-ranked solution to the Kaggle \"Flavours of Physics\" competition](https://github.com/gramolin/flavours-of-physics)"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"309bece879978e53ce04daac93497642de97d54d"},"cell_type":"code","source":"# Physical constants:\nc = 299.792458     # Speed of light\nm_mu = 105.6583715 # Muon mass (in MeV)\nm_tau = 1776.82    # Tau mass (in MeV)\n\n# List of the features for the first booster:\nlist1 = [\n# Original features:\n         'FlightDistance',\n         'FlightDistanceError',\n         'LifeTime',\n         'IP',\n         'IPSig',\n         'VertexChi2',\n         'dira',\n         'pt',\n         'DOCAone',\n         'DOCAtwo',\n         'DOCAthree',\n         'IP_p0p2',\n         'IP_p1p2',\n         'isolationa',\n         'isolationb',\n         'isolationc',\n         'isolationd',\n         'isolatione',\n         'isolationf',\n         'iso',\n         'CDF1',\n         'CDF2',\n         'CDF3',\n         'ISO_SumBDT',\n         'p0_IsoBDT',\n         'p1_IsoBDT',\n         'p2_IsoBDT',\n         'p0_track_Chi2Dof',\n         'p1_track_Chi2Dof',\n         'p2_track_Chi2Dof',\n         'p0_IP',\n         'p0_IPSig',\n         'p1_IP',\n         'p1_IPSig',\n         'p2_IP',\n         'p2_IPSig',\n# Extra features:\n         'E',\n         'FlightDistanceSig',\n         'DOCA_sum',\n         'isolation_sum',\n         'IsoBDT_sum',\n         'track_Chi2Dof',\n         'IP_sum',\n         'IPSig_sum',\n         'CDF_sum'\n        ]\n\n# List of the features for the second booster:\nlist2 = [\n# Original features:\n         'dira',\n         'pt',\n         'p0_pt',\n         'p0_p',\n         'p0_eta',\n         'p1_pt',\n         'p1_p',\n         'p1_eta',\n         'p2_pt',\n         'p2_p',\n         'p2_eta',\n# Extra features:\n         'E',\n         'pz',\n         'beta',\n         'gamma',\n         'beta_gamma',\n         'Delta_E',\n         'Delta_M',\n         'flag_M',\n         'E0',\n         'E1',\n         'E2',\n         'E0_ratio',\n         'E1_ratio',\n         'E2_ratio',\n         'p0_pt_ratio',\n         'p1_pt_ratio',\n         'p2_pt_ratio',\n         'eta_01',\n         'eta_02',\n         'eta_12',\n         't_coll'\n         ]\n\n# Function to add extra features:\ndef add_features_gramolin(df):\n  \n  # Number of events:\n  N = len(df)\n  \n  # Internal arrays:\n  p012_p = np.zeros(3)\n  p012_pt = np.zeros(3)\n  p012_z = np.zeros(3)\n  p012_eta = np.zeros(3)\n  p012_IsoBDT = np.zeros(3)\n  p012_track_Chi2Dof = np.zeros(3)\n  p012_IP = np.zeros(3)\n  p012_IPSig = np.zeros(3)\n  CDF123 = np.zeros(3)\n  isolation = np.zeros(6)\n  \n  # Kinematic features related to the mother particle:\n  E = np.zeros(N)\n  pz = np.zeros(N)\n  beta = np.zeros(N)\n  gamma = np.zeros(N)\n  beta_gamma = np.zeros(N)\n  M_lt = np.zeros(N)\n  M_inv = np.zeros(N)\n  Delta_E = np.zeros(N)\n  Delta_M = np.zeros(N)\n  flag_M = np.zeros(N)\n  \n  # Kinematic features related to the final-state particles p0, p1, and p2:\n  E012 = np.zeros((N,3))\n  E012_ratio = np.zeros((N,3))\n  p012_pt_ratio = np.zeros((N,3))\n  eta_01 = np.zeros(N)\n  eta_02 = np.zeros(N)\n  eta_12 = np.zeros(N)\n  t_coll = np.zeros(N)\n  \n  # Other extra features:\n  FlightDistanceSig = np.zeros(N)\n  DOCA_sum = np.zeros(N)\n  isolation_sum = np.zeros(N)\n  IsoBDT_sum = np.zeros(N)\n  track_Chi2Dof = np.zeros(N)\n  IP_sum = np.zeros(N)\n  IPSig_sum = np.zeros(N)\n  CDF_sum = np.zeros(N)\n  \n  for i in range(N):\n    # Read some of the original features:  \n    pt = df['pt'].values[i]\n    dira = df['dira'].values[i]\n    LifeTime = df['LifeTime'].values[i]\n    FlightDistance = df['FlightDistance'].values[i]\n    FlightDistanceError = df['FlightDistanceError'].values[i]\n    DOCAone = df['DOCAone'].values[i]\n    DOCAtwo = df['DOCAtwo'].values[i]\n    DOCAthree = df['DOCAthree'].values[i]\n    isolation[0] = df['isolationa'].values[i]\n    isolation[1] = df['isolationb'].values[i]\n    isolation[2] = df['isolationc'].values[i]\n    isolation[3] = df['isolationd'].values[i]\n    isolation[4] = df['isolatione'].values[i]\n    isolation[5] = df['isolationf'].values[i]\n    \n    for j in range(3):\n      p012_p[j] = df['p'+str(j)+'_p'].values[i]\n      p012_pt[j] = df['p'+str(j)+'_pt'].values[i]\n      p012_eta[j] = df['p'+str(j)+'_eta'].values[i]\n      p012_IsoBDT[j] = df['p'+str(j)+'_IsoBDT'].values[i]\n      p012_track_Chi2Dof[j] = df['p'+str(j)+'_track_Chi2Dof'].values[i]\n      p012_IP[j] = df['p'+str(j)+'_IP'].values[i]\n      p012_IPSig[j] = df['p'+str(j)+'_IPSig'].values[i]\n      CDF123[j] = df['CDF'+str(j+1)].values[i]\n    \n    # Differences between pseudorapidities of the final-state particles:\n    eta_01[i] = p012_eta[0] - p012_eta[1]\n    eta_02[i] = p012_eta[0] - p012_eta[2]\n    eta_12[i] = p012_eta[1] - p012_eta[2]\n    \n    # Transverse collinearity of the final-state particles (equals to 1 if they are collinear):\n    t_coll[i] = sum(p012_pt[:])/pt\n    \n    # Longitudinal momenta of the final-state particles:\n    p012_z[:] = p012_pt[:]*np.sinh(p012_eta[:])\n    \n    # Energies of the final-state particles:\n    E012[i,:] = np.sqrt(np.square(m_mu) + np.square(p012_p[:]))\n    \n    # Energy and momenta of the mother particle:\n    E[i] = sum(E012[i,:])\n    pz[i] = sum(p012_z[:])\n    p = np.sqrt(np.square(pt) + np.square(pz[i]))\n    \n    # Energies and momenta of the final-state particles relative to those of the mother particle:\n    E012_ratio[i,:] = E012[i,:]/E[i]\n    p012_pt_ratio[i,:] = p012_pt[:]/pt\n    \n    # Mass of the mother particle calculated from FlightDistance and LifeTime:\n    beta_gamma[i] = FlightDistance/(LifeTime*c)\n    M_lt[i] = p/beta_gamma[i]\n    \n    # If M_lt is around the tau mass then flag_M = 1 (otherwise 0):\n    if np.fabs(M_lt[i] - m_tau - 1.44) < 17: flag_M[i] = 1\n    \n    # Invariant mass of the mother particle calculated from its energy and momentum:        \n    M_inv[i] = np.sqrt(np.square(E[i]) - np.square(p))\n    if (np.isnan(M_inv[i])):      # mjh for about 11 records this is true\n        M_inv[i] = 0\n        gamma[i] = 0\n        beta[i] = 0\n    else:\n        # Relativistic gamma and beta of the mother particle:\n        gamma[i] = E[i]/M_inv[i]   \n        beta[i] = np.sqrt(np.square(gamma[i]) - 1.)/gamma[i]\n    \n    # Difference between M_lt and M_inv:\n    Delta_M[i] = M_lt[i] - M_inv[i]\n    \n    # Difference between energies of the mother particle calculated in two different ways:\n    Delta_E[i] = np.sqrt(np.square(M_lt[i]) + np.square(p)) - E[i]\n    \n    # Other extra features:\n    FlightDistanceSig[i] = FlightDistance/FlightDistanceError\n    DOCA_sum[i] = DOCAone + DOCAtwo + DOCAthree\n    isolation_sum[i] = sum(isolation[:])\n    IsoBDT_sum[i] = sum(p012_IsoBDT[:])\n    track_Chi2Dof[i] = np.sqrt(sum(np.square(p012_track_Chi2Dof[:] - 1.)))\n    IP_sum[i] = sum(p012_IP[:])\n    IPSig_sum[i] = sum(p012_IPSig[:])\n    CDF_sum[i] = sum(CDF123[:])\n  \n  # Kinematic features related to the mother particle:\n  df['E'] = E\n  df['pz'] = pz\n  df['beta'] = beta\n  df['gamma'] = gamma\n  df['beta_gamma'] = beta_gamma\n  df['M_lt'] = M_lt\n  df['M_inv'] = M_inv\n  df['Delta_E'] = Delta_E\n  df['Delta_M'] = Delta_M\n  df['flag_M'] = flag_M\n  \n  # Kinematic features related to the final-state particles:\n  df['E0'] = E012[:,0]\n  df['E1'] = E012[:,1]\n  df['E2'] = E012[:,2]\n  df['E0_ratio'] = E012_ratio[:,0]\n  df['E1_ratio'] = E012_ratio[:,1]\n  df['E2_ratio'] = E012_ratio[:,2]\n  df['p0_pt_ratio'] = p012_pt_ratio[:,0]\n  df['p1_pt_ratio'] = p012_pt_ratio[:,1]\n  df['p2_pt_ratio'] = p012_pt_ratio[:,2]\n  df['eta_01'] = eta_01\n  df['eta_02'] = eta_02\n  df['eta_12'] = eta_12\n  df['t_coll'] = t_coll\n  \n  # Other extra features:\n  df['FlightDistanceSig'] = FlightDistanceSig\n  df['DOCA_sum'] = DOCA_sum\n  df['isolation_sum'] = isolation_sum\n  df['IsoBDT_sum'] = IsoBDT_sum\n  df['track_Chi2Dof'] = track_Chi2Dof\n  df['IP_sum'] = IP_sum\n  df['IPSig_sum'] = IPSig_sum\n  df['CDF_sum'] = CDF_sum\n  \n  return df","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5799b9d6eafa305061bcb928c9248d1538784f30"},"cell_type":"markdown","source":"**Add Features (code)**"},{"metadata":{"trusted":true,"_uuid":"2e8c84654b62c34426c9c66753015d07fe1e201d","_kg_hide-input":true},"cell_type":"code","source":"#--------------- feature engineering -------------- #\ndef 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":{"_uuid":"f893c970ebd37ef7ea2182bd14b2405470df0ff3"},"cell_type":"markdown","source":"### Five Lines Model  \n\nOriginal kernel [Five Line Model](https://www.kaggle.com/scirpus/five-line-model)  \nThis is based on how it's used in this kernel [Ensemble with UGBC](https://www.kaggle.com/skooch/ensemble-with-ugbc)\n"},{"metadata":{"trusted":true,"_uuid":"7b5fb3a60b0bef75ff558ecf06f197d8395d845a","_kg_hide-input":true},"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,"_uuid":"31948bd87e773a4bb2841d2fe3a7581dbe7d93a9"},"cell_type":"code","source":"if DO_5_LINES or DO_5_ENS:\n    tr_preds_1 = GP(train_ugbc).values\n    test_preds_1 = GP(test_ugbc).values\n    ca_preds_1 = GP(check_agreement).values\n\n    test_predictions = pd.DataFrame({'preds_line5':test_preds_1})\n    train_predictions_all = pd.DataFrame({'id':trainids,'predictions_1':tr_preds_1})\n    ca_predictions = pd.DataFrame({'id':caids,'predictions_1':ca_preds_1})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8142c94e2ed2b3ab0bfe493989a0f605a7b2728e"},"cell_type":"code","source":"# since the target is not used for this model we can add the feature to our data without any leakage\nif DO_5_LINES:\n    train_ugbc['lines'] = tr_preds_1\n    check_agreement['lines'] = ca_preds_1\n    test_ugbc['lines'] = test_preds_1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"71078fe34e92443f9b69ea751054c8af2a15bda2"},"cell_type":"code","source":"if DO_5_LINES or DO_5_ENS:\n    agreement_probs = ca_predictions.predictions_1\n\n    ks = 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\n    print('5 line KS metric', ks, ks < 0.09)\n# print(roc_auc_truncated(y_cv, cv_predictions.predictions_1))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"37a85d95902e2c4daa04ea6140ba0e03805e0fef"},"cell_type":"markdown","source":"**add_lines (code)**"},{"metadata":{"trusted":true,"_uuid":"ca0a2fc4d07276ebab59a6c2dcc63bf84b40c3aa","_kg_hide-input":true},"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":{"_uuid":"afa348e1caeb0486195e934adf263a787be375b3"},"cell_type":"markdown","source":"## Train simple model using part of the training sample"},{"metadata":{"trusted":true,"_uuid":"3a9a5266e9f536fbaa400e0837e153b221b028e0"},"cell_type":"code","source":"if DO_5_LINES:\n    train, test, y_train, y_test, train_id, test_id, train_predictions, cv_predictions = train_test_split(train_ugbc, signal, trainids, train_predictions_all, random_state=100, test_size=0.25, shuffle=True)\n    cv_predictions = cv_predictions.copy()\n    train_predictions = train_predictions_all.copy()\nelse:\n    train, test, y_train, y_test, train_id, test_id = train_test_split(train_ugbc, signal, trainids, 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()\n\ntest = test.copy()\n\n# train on whole data set now\ntrain = train.copy()\ny_train = signal\n\n#train = X_tr.copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e67f9c1f0aa25ae08d98a9b3aea8641acf94bc22"},"cell_type":"code","source":"if DO_5_LINES:\n    train = add_lines(train)\n    test = add_lines(test)\n    test_ugbc = add_lines(test_ugbc)\n    train_ugbc = add_lines(train_ugbc)\n    check_agreement = add_lines(check_agreement)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c198b7538751445923bcfa32bb451637602dbc4b"},"cell_type":"markdown","source":"**Add features**"},{"metadata":{"trusted":true,"_uuid":"b8f5a5319fde7365680b966b916179e0230f641a"},"cell_type":"code","source":"print(\"Add features\")\ntrain_gramolin = add_features_gramolin(train.copy())\ntrain = add_features(train)\ntrain_ugbc = add_features(train_ugbc)\ntest_gramolin = add_features_gramolin(test)\ntest = add_features(test)\ntest_ugbc = add_features(test_ugbc)\ncheck_agreement = add_features(check_agreement)\nprint(\"features added...\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d22c201b27c150c2c6fd3455cc5d439294089db7"},"cell_type":"markdown","source":"**Eliminate features**"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"dbc6d80d3aeda31318b723d41c627d7b120c24c9"},"cell_type":"code","source":"print(\"Eliminate features\")\nfilter_out = ['id', 'min_ANNmuon', 'production', 'mass', 'signal',\n              'SPDhits','CDF1', 'CDF2', 'CDF3',\n              'isolationb', 'isolationc','p0_pt', 'p1_pt', 'p2_pt',\n              'p0_p', 'p1_p', 'p2_p', 'p0_eta', 'p1_eta', 'p2_eta',\n              'isolationa', 'isolationb', 'isolationc', 'isolationd', 'isolatione', 'isolationf',\n              'p0_IsoBDT', 'p1_IsoBDT', 'p2_IsoBDT',\n              'p0_IP', 'p1_IP', 'p2_IP',\n              'IP_p0p2', 'IP_p1p2',\n              'p0_track_Chi2Dof', 'p1_track_Chi2Dof', 'p2_track_Chi2Dof',\n              'p0_IPSig', 'p1_IPSig', 'p2_IPSig',\n              'DOCAone', 'DOCAtwo', 'DOCAthree',\n              'lines', 'line5', 'line6']\n              #'line10']","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d906bf2129e1bfa585d6c44993bca92e1d2ea954"},"cell_type":"markdown","source":"**Final features**"},{"metadata":{"trusted":true,"_uuid":"bea43ff3a55c559732d25d11824b7b597433f9d4"},"cell_type":"code","source":"features = list(f for f in train.columns if f not in filter_out)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b83daae3fbe5a8162835e81cde0018bbabf57342"},"cell_type":"markdown","source":"**UGBC model(s)**"},{"metadata":{"_uuid":"26936a27f0a4f6ebf990ad439361ac687fde8ca1"},"cell_type":"markdown","source":"**Gramolin models**"},{"metadata":{"trusted":true,"_uuid":"f391bec6e49e36e19feadef98669715cd6818446"},"cell_type":"code","source":"if DO_5_LINES:\n    train_gramolin = add_lines(train_gramolin)\n    test_gramolin = add_lines(test_gramolin)\n    # list1.extend(['line1','line2','line3','line4','line6','line7','line8']) # ,'line9','line10'])\n    list1.extend(['line1','line2','line3','line4','line6'])\n\n# mjh - looking to update my later importance functions to include the agreement and correlation metrics\n#       by feature. Maybe make this process of merging features, and passing metric constraint tests \n#       a little less trial and err. This version should be close to back to my best, except maybe that was just\n#       5 lines model lines1-4 only? It appeared that way when I did a quick diff. Version 59 was I think the best.\n# mjh try some of the original ugbc features and eliminate some that eliminates\n#train_gramolin = add_features(train_gramolin)\n#test_gramolin = add_features(test_gramolin)\n#list1 = list(f for f in list1 if f not in filter_out)\n#list1.extend(['NEW5_lt'])\n#list2 = list(f for f in list2 if f not in filter_out)\n#list2.extend(['NEW_FD_SUMP'])\n\nloss = BinFlatnessLossFunction(['mass'], n_bins=15, uniform_label=0, fl_coefficient=15, power=2)\ngramolin1 = UGradientBoostingClassifier(loss=loss, n_estimators=550,\n                                 max_depth=6,\n                                 learning_rate=0.15,\n                                 train_features=list1,\n                                 subsample=0.7,\n                                 random_state=123)\ngramolin1.fit(train_gramolin[list1 + ['mass']], train_gramolin['signal'])\ny_pred_gramolin1 = gramolin1.predict_proba(test_gramolin[list1])[:, 1]\nroc_auc_gramolin1 = roc_auc_score(test_gramolin['signal'], y_pred_gramolin1)\nprint(\"Gramolin 1 AUC:\",roc_auc_gramolin1)\ndf_agreement_gramolin = add_features_gramolin(df_agreement)\ndf_agreement_gramolin = add_features(df_agreement_gramolin)\nif DO_5_LINES:\n    df_agreement_gramolin = add_lines(df_agreement_gramolin)\nagreement_probs_gramolin1 = gramolin1.predict_proba(df_agreement_gramolin[list1])[:, 1]\nks = compute_ks(\n    agreement_probs_gramolin1[df_agreement_gramolin['signal'].values == 0],\n    agreement_probs_gramolin1[df_agreement_gramolin['signal'].values == 1],\n    df_agreement_gramolin[df_agreement_gramolin['signal'] == 0]['weight'].values,\n    df_agreement_gramolin[df_agreement_gramolin['signal'] == 1]['weight'].values)\ndf_corr_check_gramolin = add_features_gramolin(df_corr_check)\ndf_corr_check_gramolin = add_features(df_corr_check_gramolin)\nif DO_5_LINES:\n    df_corr_check_gramolin = add_lines(df_corr_check_gramolin)\ny_mc = gramolin1.predict_proba(df_corr_check_gramolin[list1])[:, 1]\nmc1 = check_correlation(y_mc, df_corr_check['mass'])\nprint ('Gramolin 1 KS metric:', ks, \"is OK:\", ks < 0.09,'MC metric:', mc1, \"is OK:\", mc1 < MC)\ntest_ugbc_gramolin = add_features_gramolin(test_ugbc)\ntest_ugbc_gramolin = add_features(test_ugbc_gramolin)\nif DO_5_LINES:\n    test_ugbc_gramolin = add_lines(test_ugbc_gramolin)\ntest_p_gramolin1 = gramolin1.predict_proba(test_ugbc_gramolin[list1])[:,1]\nresult1 = pd.DataFrame({'id': test_ugbc['id']})\nresult1['prediction'] = test_p_gramolin1 \nresult1.to_csv('gramolin1.csv', index=False, header=[\"id\", \"prediction\"], sep=',', mode='a')\ngramolin2 = UGradientBoostingClassifier(loss=loss, n_estimators=550,\n                                 max_depth=6,\n                                 learning_rate=0.15,\n                                 train_features=list2,\n                                 subsample=0.7,\n                                 random_state=123)\ngramolin2.fit(train_gramolin[list2 + ['mass']], train_gramolin['signal'])\ny_pred_gramolin2 = gramolin2.predict_proba(test_gramolin[list2])[:, 1]\nroc_auc_gramolin2 = roc_auc_score(test_gramolin['signal'], y_pred_gramolin2)\nprint(\"Gramolin 2 AUC:\",roc_auc_gramolin2)\nagreement_probs_gramolin2 = gramolin2.predict_proba(df_agreement_gramolin[list2])[:, 1]\n\nks = compute_ks(\n    agreement_probs_gramolin1[df_agreement_gramolin['signal'].values == 0],\n    agreement_probs_gramolin1[df_agreement_gramolin['signal'].values == 1],\n    df_agreement_gramolin[df_agreement_gramolin['signal'] == 0]['weight'].values,\n    df_agreement_gramolin[df_agreement_gramolin['signal'] == 1]['weight'].values)\ny_mc = gramolin2.predict_proba(df_corr_check_gramolin[list1])[:, 1]\nmc2 = check_correlation(y_mc, df_corr_check['mass'])\nprint ('Gramolin 2 KS metric:', ks, \"is OK:\", ks < 0.09,'MC metric:', mc2, \"is OK:\", mc2 < MC)\ntest_p_gramolin2 = gramolin2.predict_proba(test_ugbc_gramolin[list2])[:,1]\nresult2 = pd.DataFrame({'id': test_ugbc['id']})\nresult2['prediction'] = test_p_gramolin2 \nresult2.to_csv('gramolin2.csv', index=False, header=[\"id\", \"prediction\"], sep=',', mode='a')\n\np_weight = 0.94   \n# Weighted average of the predictions:\nresult = pd.DataFrame({'id': test_ugbc['id']})\nresult['prediction'] = 0.5*(p_weight*test_p_gramolin1 + (1 - p_weight)*test_p_gramolin2)\n# Write to the submission file:\nresult.to_csv('gramolin_sub.csv', index=False, header=[\"id\", \"prediction\"], sep=',', mode='a')\nif not 'test_predictions' in locals():\n    test_predictions = pd.DataFrame({'gramolin':result['prediction']})\nelse:\n    test_predictions['gramolin'] = result['prediction']\nprint(\"Gramolin Done...\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8bea62ebd08433d486b374e0c228af7e92b759ed"},"cell_type":"code","source":"#-------------------  UGBC model -------------------- #\nprint(\"Train a UGradientBoostingClassifier\")\nloss = BinFlatnessLossFunction(['mass'], n_bins=15, uniform_label=0 , fl_coefficient=15, power=2)\nugbc = UGradientBoostingClassifier(loss=loss, n_estimators=550,\n                                 max_depth=6,\n                                 learning_rate=0.15,\n                                 train_features=features,\n                                 subsample=0.7,\n                                 random_state=123)\nugbc.fit(train[features + ['mass']], train['signal'])\nprint(\"Done...\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"133930ea41518da524eaf9c3bde55e46f6ab5f05"},"cell_type":"markdown","source":"**Check model quality on the training sample**"},{"metadata":{"trusted":true,"_uuid":"79f6809f3c816c3e6d0c74dddf5dec136a32bd3a"},"cell_type":"code","source":"def plot_metrics(y_true, y_pred):\n    fpr, tpr, thresholds = roc_curve(y_true, y_pred)\n    roc_auc = roc_auc_score(y_true, y_pred)\n\n    plt.plot(fpr, tpr, label='ROC AUC=%f' % roc_auc)\n    plt.xlabel(\"FPR\")\n    plt.ylabel(\"TPR\")\n    plt.legend()\n    plt.title(\"ROC Curve\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"5197cfc506b17d615542e73bd5f3ddb83e5adb2a"},"cell_type":"code","source":"y_pred = ugbc.predict_proba(test[features])[:, 1]\n\nplot_metrics(test['signal'], y_pred)\ntrain.shape, y_pred.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6583802ce3397722395bbe7858419265b3df0064"},"cell_type":"markdown","source":"ROC AUC is just a part of the solution, you also have to make sure that\n* the classifier output is not correlated with the mass\n* classifier performs similarily on MC and real data of the normalization channel"},{"metadata":{"_uuid":"125bc97aa4d457d238f353de1872ff33f6295803"},"cell_type":"markdown","source":"## Mass correlation check\n"},{"metadata":{"trusted":true,"_uuid":"2da3cd75a3aeeab618cc22e05224e2b29887e2c5"},"cell_type":"code","source":"if DO_MASS_CORR:\n    df_corr_check = pd.read_csv(\"../input/check_correlation.csv\")\n    df_corr_check = add_features(df_corr_check)\n    df_corr_check['lines'] = GP(df_corr_check).values\n    df_corr_check = add_lines(df_corr_check)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"03cfae043968722ada711bd7072d5aed60cbcd7e"},"cell_type":"code","source":"if DO_MASS_PLOT:\n    df_corr_check.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"535cd01a4b5f1f850a0493d5b2de8343c2f3ad1a"},"cell_type":"code","source":"if DO_MASS_PLOT or DO_MASS_CORR:\n    y_pred = ugbc.predict(df_corr_check[features])","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_uuid":"1290e7bf5d360d8eefcff9b6435a0384224210e9"},"cell_type":"markdown","source":"**efficiencies (code)**"},{"metadata":{"trusted":true,"_uuid":"c3a9c48d2b91e4f4e14f9b5215e0b8645213f0a6","_kg_hide-input":true},"cell_type":"code","source":"def efficiencies(features, thresholds=None, mask=None, bins=30, labels_dict=None, ignored_sideband=0.0,\n                     errors=False, grid_columns=2):\n        \"\"\"\n        Efficiencies for spectators\n        :param features: using features (if None then use classifier's spectators)\n        :type features: None or list[str]\n        :param bins: bins for histogram\n        :type bins: int or array-like\n        :param mask: mask for data, which will be used\n        :type mask: None or numbers.Number or array-like or str or function(pandas.DataFrame)\n        :param list[float] thresholds: thresholds on prediction\n        :param bool errors: if True then use errorbar, else interpolate function\n        :param labels_dict: label -- name for class label\n            if None then {0: 'bck', '1': 'signal'}\n        :type labels_dict: None or OrderedDict(int: str)\n        :param int grid_columns: count of columns in grid\n        :param float ignored_sideband: (0, 1) percent of plotting data\n        :rtype: plotting.GridPlot\n        \"\"\"\n        mask, data, class_labels, weight = self._apply_mask(\n            mask, self._get_features(features), self.target, self.weight)\n        labels_dict = self._check_labels(labels_dict, class_labels)\n\n        plots = []\n        for feature in data.columns:\n            for name, prediction in self.prediction.items():\n                prediction = prediction[mask]\n                eff = OrderedDict()\n                for label, label_name in labels_dict.items():\n                    label_mask = class_labels == label\n                    eff[label_name] = utils.get_efficiencies(prediction[label_mask, label],\n                                                             data[feature][label_mask].values,\n                                                             bins_number=bins,\n                                                             sample_weight=weight[label_mask],\n                                                             thresholds=thresholds, errors=errors,\n                                                             ignored_sideband=ignored_sideband)\n\n                for label_name, eff_data in eff.items():\n                    if errors:\n                        plot_fig = plotting.ErrorPlot(eff_data)\n                    else:\n                        plot_fig = plotting.FunctionsPlot(eff_data)\n                    plot_fig.xlabel = feature\n                    plot_fig.ylabel = 'Efficiency for {}'.format(name)\n                    plot_fig.title = '{} flatness'.format(label_name)\n                    plot_fig.ylim = (0, 1)\n                    plots.append(plot_fig)\n\n        return plotting.GridPlot(grid_columns, *plots)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e20cc9ec39100579c4a5ed980f2f0e1037dcb6a3"},"cell_type":"markdown","source":"**check_arrays (code)**"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"79f8c767922d88cc973b821390a29b6a8b08f6aa"},"cell_type":"code","source":"def check_arrays(*arrays):\n    \"\"\"\n    Left for consistency, version of `sklearn.validation.check_arrays`\n    :param list[iterable] arrays: arrays with same length of first dimension.\n    \"\"\"\n    assert len(arrays) > 0, 'The number of array must be greater than zero'\n    checked_arrays = []\n    shapes = []\n    for arr in arrays:\n        if arr is not None:\n            checked_arrays.append(numpy.array(arr))\n            shapes.append(checked_arrays[-1].shape[0])\n        else:\n            checked_arrays.append(None)\n    assert numpy.sum(numpy.array(shapes) == shapes[0]) == len(shapes), 'Different shapes of the arrays {}'.format(\n        shapes)\n    return checked_arrays","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8450fbb21b1788118a78097eef6acd1e565677c9"},"cell_type":"markdown","source":"**get_efficiencies (code)**"},{"metadata":{"trusted":true,"_uuid":"686d1ef1994ed0ed4fd0e0454fabfb60cccb97a1","_kg_hide-input":true},"cell_type":"code","source":"def get_efficiencies(prediction, spectator, sample_weight=None, bins_number=20,\n                     thresholds=None, errors=False, ignored_sideband=0.0):\n    \"\"\"\n    Construct efficiency function dependent on spectator for each threshold\n    Different score functions available: Efficiency, Precision, Recall, F1Score,\n    and other things from sklearn.metrics\n    :param prediction: list of probabilities\n    :param spectator: list of spectator's values\n    :param bins_number: int, count of bins for plot\n    :param thresholds: list of prediction's threshold\n        (default=prediction's cuts for which efficiency will be [0.2, 0.4, 0.5, 0.6, 0.8])\n    :return:\n        if errors=False\n        OrderedDict threshold -> (x_values, y_values)\n        if errors=True\n        OrderedDict threshold -> (x_values, y_values, y_err, x_err)\n        All the parts: x_values, y_values, y_err, x_err are numpy.arrays of the same length.\n    \"\"\"\n    prediction, spectator, sample_weight = \\\n        check_arrays(prediction, spectator, sample_weight)\n\n    spectator_min, spectator_max = weighted_quantile(spectator, [ignored_sideband, (1. - ignored_sideband)])\n    mask = (spectator >= spectator_min) & (spectator <= spectator_max)\n    spectator = spectator[mask]\n    prediction = prediction[mask]\n    bins_number = min(bins_number, len(prediction))\n    sample_weight = sample_weight if sample_weight is None else numpy.array(sample_weight)[mask]\n\n    if thresholds is None:\n        thresholds = [weighted_quantile(prediction, quantiles=1 - eff, sample_weight=sample_weight)\n                      for eff in [0.2, 0.4, 0.5, 0.6, 0.8]]\n\n    binner = Binner(spectator, bins_number=bins_number)\n    if sample_weight is None:\n        sample_weight = numpy.ones(len(prediction))\n    bins_data = binner.split_into_bins(spectator, prediction, sample_weight)\n\n    bin_edges = numpy.array([spectator_min] + list(binner.limits) + [spectator_max])\n    xerr = numpy.diff(bin_edges) / 2.\n    result = OrderedDict()\n    for threshold in thresholds:\n        x_values = []\n        y_values = []\n        N_in_bin = []\n        for num, (masses, probabilities, weights) in enumerate(bins_data):\n            y_values.append(numpy.average(probabilities > threshold, weights=weights))\n            N_in_bin.append(numpy.sum(weights))\n            if errors:\n                x_values.append((bin_edges[num + 1] + bin_edges[num]) / 2.)\n            else:\n                x_values.append(numpy.mean(masses))\n\n        x_values, y_values, N_in_bin = check_arrays(x_values, y_values, N_in_bin)\n        if errors:\n            result[threshold] = (x_values, y_values, numpy.sqrt(y_values * (1 - y_values) / N_in_bin), xerr)\n        else:\n            result[threshold] = (x_values, y_values)\n    return result","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0cea3fc8a02ea64a660e4ba7885662d6002d66cc"},"cell_type":"markdown","source":"**weighted_quantile (code)**"},{"metadata":{"trusted":true,"_uuid":"4b383c8ac75dd80da3f4417baa302b71b3f92e02","_kg_hide-input":true},"cell_type":"code","source":"def weighted_quantile(array, quantiles, sample_weight=None, array_sorted=False, old_style=False):\n    \"\"\"Computing quantiles of array. Unlike the numpy.percentile, this function supports weights,\n    but it is inefficient and performs complete sorting.\n    :param array: distribution, array of shape [n_samples]\n    :param quantiles: floats from range [0, 1] with quantiles of shape [n_quantiles]\n    :param sample_weight: optional weights of samples, array of shape [n_samples]\n    :param array_sorted: if True, the sorting step will be skipped\n    :param old_style: if True, will correct output to be consistent with numpy.percentile.\n    :return: array of shape [n_quantiles]\n    Example:\n    >>> weighted_quantile([1, 2, 3, 4, 5], [0.5])\n    Out: array([ 3.])\n    >>> weighted_quantile([1, 2, 3, 4, 5], [0.5], sample_weight=[3, 1, 1, 1, 1])\n    Out: array([ 2.])\n    \"\"\"\n    array = numpy.array(array)\n    quantiles = numpy.array(quantiles)\n    sample_weight = check_sample_weight(array, sample_weight)\n    assert numpy.all(quantiles >= 0) and numpy.all(quantiles <= 1), 'Percentiles should be in [0, 1]'\n\n    if not array_sorted:\n        array, sample_weight = reorder_by_first(array, sample_weight)\n\n    weighted_quantiles = numpy.cumsum(sample_weight) - 0.5 * sample_weight\n    if old_style:\n        # To be convenient with numpy.percentile\n        weighted_quantiles -= weighted_quantiles[0]\n        weighted_quantiles /= weighted_quantiles[-1]\n    else:\n        weighted_quantiles /= numpy.sum(sample_weight)\n    return numpy.interp(quantiles, weighted_quantiles, array)\n\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"29d93eaa98eba62ec9b10c3c501a420e0ed5d990"},"cell_type":"markdown","source":"**check_sample_weight (code)**"},{"metadata":{"trusted":true,"_uuid":"4546d540926c23b0bf41f88b67bfe581b666fdda","_kg_hide-input":true},"cell_type":"code","source":"def check_sample_weight(y_true, sample_weight):\n    \"\"\"Checks the weights, if None, returns array.\n    :param y_true: labels (or any array of length [n_samples])\n    :param sample_weight: None or array of length [n_samples]\n    :return: numpy.array of shape [n_samples]\n    \"\"\"\n    if sample_weight is None:\n        return numpy.ones(len(y_true), dtype=numpy.float)\n    else:\n        sample_weight = numpy.array(sample_weight, dtype=numpy.float)\n        assert len(y_true) == len(sample_weight), \\\n            \"The length of weights is different: not {0}, but {1}\".format(len(y_true), len(sample_weight))\n        return sample_weight\n\n\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d3f7ca9f92e3a2f1ad399e777db9f0f4f2176789"},"cell_type":"markdown","source":"**reorder_by_first (code)**"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"97212a15fd9fe3656f4dc6a1c90d210450a65a12"},"cell_type":"code","source":"def reorder_by_first(*arrays):\n    \"\"\"\n    Applies the same permutation to all passed arrays,\n    permutation sorts the first passed array\n    \"\"\"\n    arrays = check_arrays(*arrays)\n    order = numpy.argsort(arrays[0])\n    return [arr[order] for arr in arrays]\n\nclass Binner(object):\n    def __init__(self, values, bins_number):\n        \"\"\"\n        Binner is a class that helps to split the values into several bins.\n        Initially an array of values is given, which is then splitted into 'bins_number' equal parts,\n        and thus we are computing limits (boundaries of bins).\n        \"\"\"\n        percentiles = [i * 100.0 / bins_number for i in range(1, bins_number)]\n        self.limits = numpy.percentile(values, percentiles)\n\n    def get_bins(self, values):\n        \"\"\"Given the values of feature, compute the index of bin\n        :param values: array of shape [n_samples]\n        :return: array of shape [n_samples]\n        \"\"\"\n        return numpy.searchsorted(self.limits, values)\n\n    def set_limits(self, limits):\n        \"\"\"Change the thresholds inside bins.\"\"\"\n        self.limits = limits\n\n    @property\n    def bins_number(self):\n        \"\"\":return: number of bins\"\"\"\n        return len(self.limits) + 1\n\n    def split_into_bins(self, *arrays):\n        \"\"\"\n        :param arrays: data to be splitted, the first array corresponds\n        :return: sequence of length [n_bins] with values corresponding to each bin.\n        \"\"\"\n        values = arrays[0]\n        for array in arrays:\n            assert len(array) == len(values), \"passed arrays have different length\"\n        bins = self.get_bins(values)\n        result = []\n        for bin in range(len(self.limits) + 1):\n            indices = bins == bin\n            result.append([numpy.array(array)[indices] for array in arrays])\n        return result\nfrom collections import OrderedDict\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"22ef1c47fd35e2f77b3907e59e2f1cbce22802d5"},"cell_type":"code","source":"if DO_MASS_PLOT:\n    eff = get_efficiencies(y_pred, df_corr_check.mass, thresholds=[0.5]) #, thresholds=[0.2, 0.4, 0.5, 0.6, 0.8])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f23d0f3f99acca2e40a37e7d9af79f9772d4dd05"},"cell_type":"code","source":"if DO_MASS_PLOT:\n    eff.keys()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a7b2ea4955a818a9ea39bcde1261d589baef59d5"},"cell_type":"code","source":"if DO_MASS_PLOT:\n    for label_name, eff_data in eff.items():\n        pyplot.plot(eff_data[0], eff_data[1], label=\"global eff  %.1f\" % label_name)\n    pyplot.xlabel('mass')\n    pyplot.ylabel('Efficiency')\n    pyplot.legend();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1eadbd732f8c4845715c913e18d53e65d251a5e7"},"cell_type":"code","source":"if DO_MASS_CORR:\n    corr_metric = check_correlation(y_pred, df_corr_check['mass'])\n    print (corr_metric)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6a6137c7e2187fa9b8ba6fc76e837170b5848853"},"cell_type":"markdown","source":"## MonteCarlo vs Real difference"},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"85d40610857c89e36000d2738e46b991f7404563"},"cell_type":"code","source":"df_agreement = add_features(df_agreement)\ndf_agreement['lines'] = GP(df_agreement).values\ndf_agreement = add_lines(df_agreement)\ndf_agreement.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d5d647abc5fe0e6cc9cddf92dbb9f67be5c8609a"},"cell_type":"code","source":"df_agreement.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"82688a6d84b7726653ceace1b2c50aa7bad89d57"},"cell_type":"code","source":"df_agreement[features].head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c9b17ca8c09f674c572130467af05d94842b5283"},"cell_type":"code","source":"agreement_probs = ugbc.predict_proba(df_agreement[features])[:, 1]\n\nks = compute_ks(\n    agreement_probs[df_agreement['signal'].values == 0],\n    agreement_probs[df_agreement['signal'].values == 1],\n    df_agreement[df_agreement['signal'] == 0]['weight'].values,\n    df_agreement[df_agreement['signal'] == 1]['weight'].values)\nprint ('UGBC KS metric:', ks, \"is OK:\", ks < 0.09)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8161450b4de7dcae0396a6b4de6d5460c8491316"},"cell_type":"markdown","source":"**plot_ks (code)**"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"72f54e1aa4530cc4437ec00450bff66d1fd7a03a"},"cell_type":"code","source":"def plot_ks(X_agreement, y_pred):\n    sig_ind = X_agreement[X_agreement['signal'] == 1].index\n    bck_ind = X_agreement[X_agreement['signal'] == 0].index\n\n    mc_prob = y_pred[sig_ind]\n    mc_weight = numpy.array(X_agreement.loc[sig_ind]['weight'])\n    data_prob = y_pred[bck_ind]\n    data_weight = numpy.array(X_agreement.loc[bck_ind]['weight'])\n    inds = data_weight < 0\n    mc_weight = numpy.array(list(mc_weight) + list(-data_weight[inds]))\n    mc_prob = numpy.array(list(mc_prob) + list(data_prob[inds]))\n    data_prob = data_prob[data_weight >= 0]\n    data_weight = data_weight[data_weight >= 0]\n    hist(data_prob, weights=data_weight, color='r', histtype='step', normed=True, bins=60, label='data')\n    hist(mc_prob, weights=mc_weight, color='b', histtype='step', normed=True, bins=60, label='mc')\n    xlabel(\"prediction\")\n    legend(loc=2)\n    show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"66e4f0f45c0664d76713298dd33459793d87d932"},"cell_type":"code","source":"plot_ks(df_agreement, agreement_probs)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"851af1e5b383c1e0a3d6ddf0ea81536008d1cc11"},"cell_type":"markdown","source":"## Let's see if adding some noise can improve the agreement"},{"metadata":{"_uuid":"def5656de589ac9d7d219def76026fab57a0859f"},"cell_type":"markdown","source":"**add_noise (code)**"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"7e7e263b7cd45279790b311cb436a8ac08e5fef4"},"cell_type":"code","source":"def add_noise(array, level=0.40, random_seed=34):\n    numpy.random.seed(random_seed)\n    return level * numpy.random.random(size=array.size) + (1 - level) * array","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d5249267c20d93a80f2da85bcb06f0b2eec62d64"},"cell_type":"code","source":"if DO_NOISE:\n    agreement_probs_noise = add_noise(ugbc.predict_proba(df_agreement[features])[:, 1])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"cec8b234d282fb19d9aa481acc10ad29dff9d15e"},"cell_type":"code","source":"if DO_NOISE:\n    ks_noise = compute_ks(\n        agreement_probs_noise[df_agreement['signal'].values == 0],\n        agreement_probs_noise[df_agreement['signal'].values == 1],\n        df_agreement[df_agreement['signal'] == 0]['weight'].values,\n        df_agreement[df_agreement['signal'] == 1]['weight'].values)\n    print ('KS metric:', ks_noise, \"is OK:\", ks_noise < 0.09)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"ee7cc568e0172d481bb07477b76c1f3c19f73b8d"},"cell_type":"code","source":"if DO_NOISE:\n    plot_ks(df_agreement, agreement_probs_noise)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ca6081b96dabcfc1d3a3377db02106f7c2a767fa"},"cell_type":"markdown","source":"## Check ROC with noise"},{"metadata":{"trusted":true,"_uuid":"02610f077ca24a766c86f7aecfbe124e5d2b323f"},"cell_type":"code","source":"if DO_NOISE:\n    test.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"11696a18576fa8beec1390562be4a25e026c1405"},"cell_type":"code","source":"if DO_NOISE:\n    y_pred = add_noise(ugbc.predict_proba(test[features])[:, 1])\n\n    plot_metrics(test['signal'], y_pred)\n    test.shape, y_pred.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c29cfaeef5403a6affd02a0f4e4a76ad9fae75a7"},"cell_type":"markdown","source":"Feature selection based on [Beware Default Random Forest Importances](http://explained.ai/rf-importance/index.html)"},{"metadata":{"_uuid":"1b4eb3ae5c91d1cf8542639655e0639690b326c4"},"cell_type":"markdown","source":"If you have your own model(s) that you are ensembling you should implement the followng method yourself. I am trying to keep the rest of the following code non-dependent on specifics of the ensembled models. But for weighting on the Gramolin, this seems necessary?"},{"metadata":{"trusted":true,"_uuid":"b2d8718fd6dc4812fd68567216c8cae4164cf397"},"cell_type":"code","source":"def ensemble_preds(preds):\n    \"\"\"\n    Take the mean or apply whatever weighting you want to the passed predictions.\n    Returns the single ensembled prediction set\n    \"\"\"\n    p_weight = 0.94\n    return 0.5*(p_weight*preds[0,] + (1 - p_weight)*preds[1,])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3394b8ed0ce2957a7c350a6b0e9c2c4d4364e74f"},"cell_type":"code","source":"def ensemble_metric(models, X_train, y_train, features):\n    preds = []\n    for model in models:\n        preds.append(model.predict_proba(X_train[features])[:, 1])\n    ens_preds = ensemble_preds(preds)\n    roc_auc = roc_auc_truncated(y_train, ens_preds)\n    return roc_auc","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1d5042098e7e8b34382935e23ed49617f23ca285"},"cell_type":"code","source":"def metric(model, X_train, y_train, features):\n    y_pred = model.predict_proba(X_train[features])[:, 1]\n    y_agree = model.predict_proba(df_agreement[features])[:, 1]\n    y_corr = model.predict_proba(df_corr_check[features])[:, 1]\n#    roc_auc = roc_auc_score(y_train, y_pred) \n# use truncated since thats what evaluation is actually done on\n    roc_auc = roc_auc_truncated(y_train, y_pred)\n    ks_imp = compute_ks(\n        y_agree[df_agreement['signal'].values == 0],\n        y_agree[df_agreement['signal'].values == 1],\n        df_agreement[df_agreement['signal'] == 0]['weight'].values,\n        df_agreement[df_agreement['signal'] == 1]['weight'].values)\n    corr_metric = check_correlation(y_corr, df_corr_check['mass'])\n    return roc_auc, ks_imp, corr_metric","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0d0551d4204b02257373dc979685be8b5b06a46a"},"cell_type":"code","source":"def perm_ens_imp(models, X_train, y_train, feature_sets):\n    base_auc, base_agree, base_corr = ensemble_metric(model, X_train, y_train, feature_sets)\n    print(\"Baseline = \",base_auc,\"KS metric:\",base_agree,\"MC metric:\",base_corr)\n    imp_auc = [], imp_agree = [], imp_corr = []\n    for features in feature_sets:\n        for col in features:\n            save = X_train[col].copy()\n            X_train[col] = np.random.permutation(X_train[col])\n            m_auc, m_agree, m_corr = metric(model, X_train, y_train, features)\n            X_train[col] = save\n            imp_auc.append(base_auc - m_auc)\n            imp_agree.append(base_agree - m_agree)\n            imp_corr.append(base_corr - m_corr)\n    return np.array([imp_auc, imp_agree, imp_corr])    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"dcd40fbdd99d6fe590f74057ab3ee1a9611ee64d"},"cell_type":"code","source":"def permutation_importances(model, X_train, y_train, features):\n    base_auc, base_agree, base_corr = metric(model, X_train, y_train, features)\n    print(\"Baseline = \",base_auc,\"KS metric:\",base_agree,\"MC metric:\",base_corr)\n    imp_auc = [], imp_agree = [], imp_corr = []\n    for col in features:\n        save = X_train[col].copy()\n        X_train[col] = np.random.permutation(X_train[col])\n        m_auc, m_agree, m_corr = metric(model, X_train, y_train, features)\n        X_train[col] = save\n        imp_auc.append(base_auc - m_auc)\n        imp_agree.append(base_agree - m_agree)\n        imp_corr.append(base_corr - m_corr)\n    return np.array([imp_auc, imp_agree, imp_corr])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"fbf1ab35ebf36356a0102fca3f9035b823b551b3"},"cell_type":"code","source":"if DO_IMP:\n    print(\"Importance preparation\")\n    filter_imp = ['id', 'min_ANNmuon', 'production', 'mass', 'signal',\n                 'iso', 'flight_dist_sig2', 'isolatione', 'isolationa', 'iso_min',\n                 'isolationf']\n    features_imp = list(f for f in train.columns if f not in filter_imp)\n    ugbc_imp = UGradientBoostingClassifier(loss=loss, n_estimators=550,\n                                 max_depth=6,\n                                 learning_rate=0.15,\n                                 train_features=features_imp,\n                                 subsample=0.7,\n                                 random_state=123)\n    ugbc_imp.fit(train[features_imp + ['mass']], train['signal'])\n    print(\"Importance preparation complete...\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"94491007e3949fa9699fa159633b3d8a281b8f35","scrolled":false},"cell_type":"code","source":"if DO_IMP:\n    imp = permutation_importances(ugbc_imp, train, train['signal'], features_imp)\n    imp_df=pd.DataFrame(data=imp,index=features_imp).sort_values(0,ascending=False)\n    print(imp_df.tail(10))\n    roc_auc_truncated(train['signal'],ugbc_imp.predict_proba(train[features_imp])[:, 1])\n    features = features_imp","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"29c1e2f473855c3b12c828fa4a68b9db22861b48"},"cell_type":"code","source":"# models and feature sets are assumed to be from prior gramolin run\nif DO_GRAMOLIN_IMP:\n    imp = perm_ens_imp([gramolin1,gramolin2], train['signal'], [list1,list2])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"608a6264f481d6998d8d34c62affcd83500396de"},"cell_type":"markdown","source":"## Train the model using the whole training sample"},{"metadata":{"trusted":true,"_uuid":"f832c13b04fbdad0ebd2f8db5f574d20cf95350b"},"cell_type":"code","source":"%time ugbc.fit(train_ugbc[features+['mass']], train_ugbc['signal'])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"541d477d005cd2f29cf03d212b9fbbf393befb6e"},"cell_type":"markdown","source":"Compute prediction (noise is not added as this model doesn't need it, as Coursera one did, and ROC get's worse)"},{"metadata":{"trusted":true,"_uuid":"cdd28d5e9eb0d174eae426666c8b3bc19404157b"},"cell_type":"code","source":"#--------------------  prediction ---------------------#\nprint ('----------------------------------------------')\nprint(\"Make predictions on the test set\")\ntest_probs = ugbc.predict_proba(test_ugbc[features])[:,1]\nif DO_5_ENS or DO_GRAMOLIN:\n    test_predictions['ugbc_pred'] = test_probs\nsubmission = pd.DataFrame({\"id\": test_ugbc[\"id\"], \"prediction\": test_probs})\nsubmission.to_csv(\"ugbc_features.csv\", index=False)\nprint(\"UGBC Predictions done...\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"44127ce9455a2fedc6e208b4b250fd459b9a9805"},"cell_type":"markdown","source":"**Average Our Predictions**"},{"metadata":{"trusted":true,"_uuid":"8f41956388ac98d34abac5ec23cfe32d8b678061"},"cell_type":"code","source":"if DO_5_ENS or DO_GRAMOLIN:\n#    test_predictions['avg_preds'] = test_predictions.mean(axis=1)\n    g_weight = .92      # gramolin weight\n    u_weight = .08      # ugbc weight\n    # test_predictions['prediction'] = add_noise(0.5*((g_weight*test_predictions['gramolin']) + (u_weight*test_predictions['ugbc_pred'])),.2)\n    test_predictions['prediction'] = 0.5*((g_weight*test_predictions['gramolin']) + (u_weight*test_predictions['ugbc_pred']))\n    test_predictions['id'] = test_ugbc['id']\n    test_predictions[['id', 'prediction']].to_csv(\"ensembled.csv\", index=False, header=[\"id\", \"prediction\"])\n    if (DO_5_ENS):\n        g_weight = .49      # gramolin weight\n        u_weight = .01      # ugbc weight\n        f_weight = .50      # five lines weight\n        # test_predictions['prediction'] = add_noise((g_weight*test_predictions['gramolin'] + (u_weight*test_predictions['ugbc_pred']) + f_weight*test_preds_1)/3)\n        test_predictions['prediction'] = (g_weight*test_predictions['gramolin'] + (u_weight*test_predictions['ugbc_pred']) + f_weight*test_preds_1)/3\n        test_predictions[['id', 'prediction']].to_csv(\"ensemble_all.csv\", index=False, header=[\"id\", \"prediction\"])\nprint(\"Ensemble Predictions done...\")","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}