import warnings warnings.filterwarnings("ignore", category=DeprecationWarning) import pandas as pd import numpy as np import os from sklearn.ensemble import AdaBoostClassifier from sklearn.metrics import roc_curve, roc_auc_score, auc from sklearn.model_selection import train_test_split from sklearn.utils.validation import column_or_1d from hep_ml.gradientboosting import UGradientBoostingClassifier, LogLossFunction from hep_ml.uboost import uBoostClassifier, uBoostBDT from hep_ml.losses import BinFlatnessLossFunction from hep_ml.losses import KnnAdaLossFunction from hep_ml.losses import KnnFlatnessLossFunction from hep_ml.losses import CompositeLossFunction import random from hep_ml.nnet import MLPClassifier import matplotlib.pyplot as plt from hep_ml import metrics import time import math from scipy.optimize import minimize #%% # -------------- loading data files -------------- # DO_ON_Serv = True if DO_ON_Serv: print("Load the train/test/eval data using pandas") data_path = '../input' train = pd.read_csv("".join((data_path, '//training.csv'))) train = train[train['min_ANNmuon'] > 0.4] test = pd.read_csv("".join((data_path, '//test.csv'))) check_agreement = pd.read_csv("".join((data_path, '//check_agreement.csv')), index_col='id') corr_check = pd.read_csv("".join((data_path, '//check_correlation.csv'))) signal = train.signal trainids = train.index.values else: print("Load the train/test/eval data using pandas") data_path = os.path.join("D:/", "Timur", "Python", "Data", "LHC") train = pd.read_csv("".join((data_path, '\\training.csv'))) train = train[train['min_ANNmuon'] > 0.4] test = pd.read_csv("".join((data_path, '\\test.csv'))) check_agreement = pd.read_csv("".join((data_path, '\\check_agreement.csv')), index_col='id') corr_check = pd.read_csv("".join((data_path, '\\check_correlation.csv'))) signal = train.signal trainids = train.index.values #%% #--------------- metrics ------------------ # # maybe these could be tuned to better match up training and leaderboard results? # these have seemed off to me sometimes def __rolling_window(data, window_size): """ Rolling window: take window with definite size through the array :param data: array-like :param window_size: size :return: the sequence of windows Example: data = array(1, 2, 3, 4, 5, 6), window_size = 4 Then this function return array(array(1, 2, 3, 4), array(2, 3, 4, 5), array(3, 4, 5, 6)) """ shape = data.shape[:-1] + (data.shape[-1] - window_size + 1, window_size) strides = data.strides + (data.strides[-1],) return np.lib.stride_tricks.as_strided(data, shape=shape, strides=strides) def __cvm(subindices, total_events): """ Compute Cramer-von Mises metric. Compared two distributions, where first is subset of second one. Assuming that second is ordered by ascending :param subindices: indices of events which will be associated with the first distribution :param total_events: count of events in the second distribution :return: cvm metric """ target_distribution = np.arange(1, total_events + 1, dtype='float') / total_events subarray_distribution = np.cumsum(np.bincount(subindices, minlength=total_events), dtype='float') subarray_distribution /= 1.0 * subarray_distribution[-1] return np.mean((target_distribution - subarray_distribution) ** 2) def check_correlation(probabilities, mass): probabilities, mass = map(column_or_1d, [probabilities, mass]) y_pred = np.zeros(shape=(len(probabilities), 2)) y_pred[:, 1] = probabilities y_pred[:, 0] = 1 - probabilities y_true = [0] * len(probabilities) df_mass = pd.DataFrame({'mass': mass}) cvm = metrics.BinBasedCvM(uniform_features=['mass'], uniform_label=0, n_bins=15) # mjh overrode n_bins number cvm.fit(df_mass, y_true) return cvm(y_true, y_pred, sample_weight=None) # knn based version of CvM def check_correlation_knn(probabilities, mass): probabilities, mass = map(column_or_1d, [probabilities, mass]) y_pred = np.zeros(shape=(len(probabilities), 2)) y_pred[:, 1] = probabilities y_pred[:, 0] = 1 - probabilities y_true = [0] * len(probabilities) df_mass = pd.DataFrame({'mass': mass}) cvm = metrics.KnnBasedCvM(uniform_features=['mass'], uniform_label=0, n_neighbours=50) cvm.fit(df_mass, y_true) return cvm(y_true, y_pred, sample_weight=None) def roc_auc_truncated(labels, predictions, tpr_thresholds=(0.2, 0.4, 0.6, 0.8), roc_weights=(4, 3, 2, 1, 0)): """ Compute weighted area under ROC curve. :param labels: array-like, true labels :param predictions: array-like, predictions :param tpr_thresholds: array-like, true positive rate thresholds delimiting the ROC segments :param roc_weights: array-like, weights for true positive rate segments :return: weighted AUC """ assert np.all(predictions >= 0.) and np.all(predictions <= 1.), 'Data predictions are out of range [0, 1]' assert len(tpr_thresholds) + 1 == len(roc_weights), 'Incompatible lengths of thresholds and weights' fpr, tpr, _ = roc_curve(labels, predictions) area = 0. tpr_thresholds = [0.] + list(tpr_thresholds) + [1.] for index in range(1, len(tpr_thresholds)): tpr_cut = np.minimum(tpr, tpr_thresholds[index]) tpr_previous = np.minimum(tpr, tpr_thresholds[index - 1]) area += roc_weights[index - 1] * (auc(fpr, tpr_cut, reorder=True) - auc(fpr, tpr_previous, reorder=True)) tpr_thresholds = np.array(tpr_thresholds) # roc auc normalization to be 1 for an ideal classifier area /= np.sum((tpr_thresholds[1:] - tpr_thresholds[:-1]) * np.array(roc_weights)) return area def __roc_curve_splitted(data_zero, data_one, sample_weights_zero, sample_weights_one): """ Compute roc curve :param data_zero: 0-labeled data :param data_one: 1-labeled data :param sample_weights_zero: weights for 0-labeled data :param sample_weights_one: weights for 1-labeled data :return: roc curve """ labels = [0] * len(data_zero) + [1] * len(data_one) weights = np.concatenate([sample_weights_zero, sample_weights_one]) data_all = np.concatenate([data_zero, data_one]) fpr, tpr, _ = roc_curve(labels, data_all, sample_weight=weights) return fpr, tpr def compute_cvm(predictions, masses, n_neighbours=200, step=50): """ Computing Cramer-von Mises (cvm) metric on background events: take average of cvms calculated for each mass bin. In each mass bin global prediction's cdf is compared to prediction's cdf in mass bin. :param predictions: array-like, predictions :param masses: array-like, in case of Kaggle tau23mu this is reconstructed mass :param n_neighbours: count of neighbours for event to define mass bin :param step: step through sorted mass-array to define next center of bin :return: average cvm value """ predictions = np.array(predictions) masses = np.array(masses) assert len(predictions) == len(masses) # First, reorder by masses predictions = predictions[np.argsort(masses)] # Second, replace probabilities with order of probability among other events predictions = np.argsort(np.argsort(predictions, kind='mergesort'), kind='mergesort') # Now, each window forms a group, and we can compute contribution of each group to CvM cvms = [] for window in __rolling_window(predictions, window_size=n_neighbours)[::step]: cvms.append(__cvm(subindices=window, total_events=len(predictions))) return np.mean(cvms) def compute_ks(data_prediction, mc_prediction, weights_data, weights_mc): """ Compute Kolmogorov-Smirnov (ks) distance between real data predictions cdf and Monte Carlo one. :param data_prediction: array-like, real data predictions :param mc_prediction: array-like, Monte Carlo data predictions :param weights_data: array-like, real data weights :param weights_mc: array-like, Monte Carlo weights :return: ks value """ assert len(data_prediction) == len(weights_data), 'Data length and weight one must be the same' assert len(mc_prediction) == len(weights_mc), 'Data length and weight one must be the same' data_prediction, mc_prediction = np.array(data_prediction), np.array(mc_prediction) weights_data, weights_mc = np.array(weights_data), np.array(weights_mc) assert np.all(data_prediction >= 0.) and np.all(data_prediction <= 1.), 'Data predictions are out of range [0, 1]' assert np.all(mc_prediction >= 0.) and np.all(mc_prediction <= 1.), 'MC predictions are out of range [0, 1]' weights_data /= np.sum(weights_data) weights_mc /= np.sum(weights_mc) fpr, tpr = __roc_curve_splitted(data_prediction, mc_prediction, weights_data, weights_mc) Dnm = np.max(np.abs(fpr - tpr)) return Dnm #--------------- metric importances ----------------# def mc_metric(model, corr_check, features): y_corr = model.predict_proba(corr_check[features])[:, 1] corr_metric = check_correlation(y_corr, corr_check['mass']) return corr_metric def mc_importances(model, corr_check, features): base_mc = mc_metric(model, corr_check, features) print("MC metric:",base_mc) imp_mc = [] for col in features: save = corr_check[col].copy() corr_check[col] = np.random.permutation(corr_check[col]) mc = mc_metric(model, corr_check, features) corr_check[col] = save imp_mc.append(mc - base_mc) return imp_mc def auc_metric(model, X_train, y_train, features): y_pred = model.predict_proba(X_train[features])[:, 1] roc_auc = roc_auc_truncated(y_train, y_pred) return roc_auc def auc_importances(model, X_train, y_train, features): base_auc = auc_metric(model, X_train, y_train, features) print("Baseline = ",base_auc) imp_auc = [] for col in features: save = X_train[col].copy() X_train[col] = np.random.permutation(X_train[col]) m_auc = auc_metric(model, X_train, y_train, features) X_train[col] = save imp_auc.append(base_auc - m_auc) return imp_auc #--------------- feature engineering -------------- # def Output(p): return 1/(1.+np.exp(-p)) def GP(data): 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"])))))) + 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"])))))))))) + 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))))))) + 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"]))) + 1.0*np.tanh((((((((((11.05855369567871094)) / (((((11.05855751037597656)) + (((((data["SPDhits"]) / 2.0)) / 2.0)))/2.0)))) - (data["IP"]))) - (((data["SPDhits"]) / (data["p1_pt"]))))) * 2.0)) + 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)) + 1.0*np.tanh((((-1.0*((((data["DOCAthree"]) / (data["CDF2"])))))) + (np.minimum(((((data["p2_pt"]) / (data["p0_p"])))), ((np.minimum(((data["CDF2"])), ((((np.sin((1.570796))) / 2.0)))))))))) + 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"])))))) / (0.318310)))))) + 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)) + 1.0*np.tanh(((((np.cos(((((data["ISO_SumBDT"]) + (0.318310))/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))))))) p1 = 11.05855369567871094 p2 = 0.318310 p3 = 1.570796 #%% ###### GRAMOLIN # Physical constants: c = 299.792458 # Speed of light m_mu = 105.6583715 # Muon mass (in MeV) m_tau = 1776.82 # Tau mass (in MeV) # List of the features for the first booster: list1 = [ # Original features: 'FlightDistance', 'FlightDistanceError', 'LifeTime', 'IP', 'IPSig', 'VertexChi2', 'dira', 'pt', 'DOCAone', 'DOCAtwo', 'DOCAthree', 'IP_p0p2', 'IP_p1p2', 'isolationa', 'isolationb', 'isolationc', 'isolationd', 'isolatione', 'isolationf', 'iso', 'CDF1', 'CDF2', 'CDF3', 'ISO_SumBDT', 'p0_IsoBDT', 'p1_IsoBDT', 'p2_IsoBDT', 'p0_track_Chi2Dof', 'p1_track_Chi2Dof', 'p2_track_Chi2Dof', 'p0_IP', 'p0_IPSig', 'p1_IP', 'p1_IPSig', 'p2_IP', 'p2_IPSig', # Extra features: 'E', 'FlightDistanceSig', 'DOCA_sum', 'isolation_sum', 'IsoBDT_sum', 'track_Chi2Dof', 'IP_sum', 'IPSig_sum', 'CDF_sum' ] # List of the features for the second booster: list2 = [ # Original features: 'dira', 'pt', 'p0_pt', 'p0_p', 'p0_eta', 'p1_pt', 'p1_p', 'p1_eta', 'p2_pt', 'p2_p', 'p2_eta', # Extra features: 'E', 'pz', 'beta', 'gamma', 'beta_gamma', 'Delta_E', 'Delta_M', 'flag_M', 'E0', 'E1', 'E2', 'E0_ratio', 'E1_ratio', 'E2_ratio', 'p0_pt_ratio', 'p1_pt_ratio', 'p2_pt_ratio', 'eta_01', 'eta_02', 'eta_12', 't_coll' ] # Function to add extra features: def add_features(df): df['NEW_FD_SUMP']=df['FlightDistance']/(df['p0_p']+df['p1_p']+df['p2_p']) df['NEW5_lt']=df['LifeTime']*(df['p0_IP']+df['p1_IP']+df['p2_IP'])/3 df['p_track_Chi2Dof_MAX'] = df.loc[:, ['p0_track_Chi2Dof', 'p1_track_Chi2Dof', 'p2_track_Chi2Dof']].max(axis=1) #df['flight_dist_sig'] = df['FlightDistance']/df['FlightDistanceError'] # modified to: df['flight_dist_sig2'] = (df['FlightDistance']/df['FlightDistanceError'])**2 # features from phunter df['flight_dist_sig'] = df['FlightDistance']/df['FlightDistanceError'] df['NEW_IP_dira'] = df['IP']*df['dira'] df['p0p2_ip_ratio']=df['IP']/df['IP_p0p2'] df['p1p2_ip_ratio']=df['IP']/df['IP_p1p2'] df['DCA_MAX'] = df.loc[:, ['DOCAone', 'DOCAtwo', 'DOCAthree']].max(axis=1) df['iso_bdt_min'] = df.loc[:, ['p0_IsoBDT', 'p1_IsoBDT', 'p2_IsoBDT']].min(axis=1) df['iso_min'] = df.loc[:, ['isolationa', 'isolationb', 'isolationc','isolationd', 'isolatione', 'isolationf']].min(axis=1) # My: # new combined features just to minimize their number; # their physical sense doesn't matter df['NEW_iso_abc'] = df['isolationa']*df['isolationb']*df['isolationc'] df['NEW_iso_def'] = df['isolationd']*df['isolatione']*df['isolationf'] df['NEW_pN_IP'] = df['p0_IP']+df['p1_IP']+df['p2_IP'] df['NEW_pN_p'] = df['p0_p']+df['p1_p']+df['p2_p'] df['NEW_IP_pNpN'] = df['IP_p0p2']*df['IP_p1p2'] df['NEW_pN_IPSig'] = df['p0_IPSig']+df['p1_IPSig']+df['p2_IPSig'] #My: # "super" feature changing the result from 0.988641 to 0.991099 df['NEW_FD_LT']=df['FlightDistance']/df['LifeTime'] #### Gramoline # Number of events: N = len(df) # Internal arrays: p012_p = np.zeros(3) p012_pt = np.zeros(3) p012_z = np.zeros(3) p012_eta = np.zeros(3) # p012_IsoBDT = np.zeros(3) # p012_track_Chi2Dof = np.zeros(3) # p012_IP = np.zeros(3) # p012_IPSig = np.zeros(3) # CDF123 = np.zeros(3) # isolation = np.zeros(6) # Kinematic features related to the mother particle: E = np.zeros(N) pz = np.zeros(N) beta = np.zeros(N) gamma = np.zeros(N) beta_gamma = np.zeros(N) M_lt = np.zeros(N) M_inv = np.zeros(N) # Delta_E = np.zeros(N) # Delta_M = np.zeros(N) flag_M = np.zeros(N) # Kinematic features related to the final-state particles p0, p1, and p2: E012 = np.zeros((N,3)) # E012_ratio = np.zeros((N,3)) # p012_pt_ratio = np.zeros((N,3)) # eta_01 = np.zeros(N) # eta_02 = np.zeros(N) # eta_12 = np.zeros(N) # t_coll = np.zeros(N) # Other extra features: # FlightDistanceSig = np.zeros(N) # DOCA_sum = np.zeros(N) # isolation_sum = np.zeros(N) # IsoBDT_sum = np.zeros(N) # track_Chi2Dof = np.zeros(N) # IP_sum = np.zeros(N) # IPSig_sum = np.zeros(N) # CDF_sum = np.zeros(N) for i in range(N): # Read some of the original features: pt = df['pt'].values[i] LifeTime = df['LifeTime'].values[i] FlightDistance = df['FlightDistance'].values[i] # FlightDistanceError = df['FlightDistanceError'].values[i] # DOCAone = df['DOCAone'].values[i] # DOCAtwo = df['DOCAtwo'].values[i] # DOCAthree = df['DOCAthree'].values[i] # isolation[0] = df['isolationa'].values[i] # isolation[1] = df['isolationb'].values[i] # isolation[2] = df['isolationc'].values[i] # isolation[3] = df['isolationd'].values[i] # isolation[4] = df['isolatione'].values[i] # isolation[5] = df['isolationf'].values[i] for j in range(3): p012_p[j] = df['p'+str(j)+'_p'].values[i] p012_pt[j] = df['p'+str(j)+'_pt'].values[i] p012_eta[j] = df['p'+str(j)+'_eta'].values[i] # p012_IsoBDT[j] = df['p'+str(j)+'_IsoBDT'].values[i] # p012_track_Chi2Dof[j] = df['p'+str(j)+'_track_Chi2Dof'].values[i] # p012_IP[j] = df['p'+str(j)+'_IP'].values[i] # p012_IPSig[j] = df['p'+str(j)+'_IPSig'].values[i] # CDF123[j] = df['CDF'+str(j+1)].values[i] # Differences between pseudorapidities of the final-state particles: # eta_01[i] = p012_eta[0] - p012_eta[1] # eta_02[i] = p012_eta[0] - p012_eta[2] # eta_12[i] = p012_eta[1] - p012_eta[2] # Transverse collinearity of the final-state particles (equals to 1 if they are collinear): # t_coll[i] = sum(p012_pt[:])/pt # Longitudinal momenta of the final-state particles: p012_z[:] = p012_pt[:]*np.sinh(p012_eta[:]) # Energies of the final-state particles: E012[i,:] = np.sqrt(np.square(m_mu) + np.square(p012_p[:])) # Energy and momenta of the mother particle: E[i] = sum(E012[i,:]) pz[i] = sum(p012_z[:]) p = np.sqrt(np.square(pt) + np.square(pz[i])) # Energies and momenta of the final-state particles relative to those of the mother particle: # E012_ratio[i,:] = E012[i,:]/E[i] # p012_pt_ratio[i,:] = p012_pt[:]/pt # Mass of the mother particle calculated from FlightDistance and LifeTime: beta_gamma[i] = FlightDistance/(LifeTime*c) M_lt[i] = p/beta_gamma[i] # If M_lt is around the tau mass then flag_M = 1 (otherwise 0): if np.fabs(M_lt[i] - m_tau - 1.44) < 17: flag_M[i] = 1 # Invariant mass of the mother particle calculated from its energy and momentum: M_inv[i] = np.sqrt(np.square(E[i]) - np.square(p)) if (np.isnan(M_inv[i])): # mjh for about 11 records this is true M_inv[i] = 0 gamma[i] = 0 beta[i] = 0 else: # Relativistic gamma and beta of the mother particle: gamma[i] = E[i]/M_inv[i] beta[i] = np.sqrt(np.square(gamma[i]) - 1.)/gamma[i] # Difference between M_lt and M_inv: # Delta_M[i] = M_lt[i] - M_inv[i] # Difference between energies of the mother particle calculated in two different ways: # Delta_E[i] = np.sqrt(np.square(M_lt[i]) + np.square(p)) - E[i] # Other extra features: # FlightDistanceSig[i] = FlightDistance/FlightDistanceError # DOCA_sum[i] = DOCAone + DOCAtwo + DOCAthree # isolation_sum[i] = sum(isolation[:]) # IsoBDT_sum[i] = sum(p012_IsoBDT[:]) # track_Chi2Dof[i] = np.sqrt(sum(np.square(p012_track_Chi2Dof[:] - 1.))) # IP_sum[i] = sum(p012_IP[:]) # IPSig_sum[i] = sum(p012_IPSig[:]) # CDF_sum[i] = sum(CDF123[:]) # Kinematic features related to the mother particle: df['E'] = E df['pz'] = pz df['beta'] = beta df['gamma'] = gamma df['beta_gamma'] = beta_gamma df['M_lt'] = M_lt df['M_inv'] = M_inv # df['Delta_E'] = Delta_E # df['Delta_M'] = Delta_M df['flag_M'] = flag_M # Kinematic features related to the final-state particles: df['E0'] = E012[:,0] df['E1'] = E012[:,1] df['E2'] = E012[:,2] # df['E0_ratio'] = E012_ratio[:,0] # df['E1_ratio'] = E012_ratio[:,1] # df['E2_ratio'] = E012_ratio[:,2] # df['p0_pt_ratio'] = p012_pt_ratio[:,0] # df['p1_pt_ratio'] = p012_pt_ratio[:,1] # df['p2_pt_ratio'] = p012_pt_ratio[:,2] # df['eta_01'] = eta_01 # df['eta_02'] = eta_02 # df['eta_12'] = eta_12 # df['t_coll'] = t_coll # Other extra features: # df['FlightDistanceSig'] = FlightDistanceSig # df['DOCA_sum'] = DOCA_sum # df['isolation_sum'] = isolation_sum # df['IsoBDT_sum'] = IsoBDT_sum # df['track_Chi2Dof'] = track_Chi2Dof # df['IP_sum'] = IP_sum # df['IPSig_sum'] = IPSig_sum # df['CDF_sum'] = CDF_sum df['GP'] = GP(df) # df['line1'] = 1.0*np.tanh(((((((((df["IPSig"]) + (df["ISO_SumBDT"]))) - (np.minimum(((-2.0)), ((df["ISO_SumBDT"])))))) / (df["ISO_SumBDT"]))) / (np.minimum((((-1.0*((df["ISO_SumBDT"]))))), ((df["IPSig"])))))) # df['line2'] = 1.0*np.tanh((-1.0*((((df["iso"]) + (((((((((((((df["VertexChi2"]) + ((3.0)))) / (df["ISO_SumBDT"]))) * (df["IP"]))) * 2.0)) / (df["ISO_SumBDT"]))) * (((((((((((df["VertexChi2"]) + ((3.0)))) / (df["ISO_SumBDT"]))) * (df["IP"]))) * 2.0)) / (df["ISO_SumBDT"])))))))))) # df['line3'] = 1.0*np.tanh((-1.0*(((((((((df["IPSig"]) * ((((df["iso"]) + (((df["IP"]) * 2.0)))/2.0)))) + (np.tanh((df["p0_IsoBDT"]))))/2.0)) * ((((df["p0_IsoBDT"]) + (df["IPSig"]))/2.0))))))) # df['line4'] = 1.0*np.tanh(((np.minimum(((np.cos((((np.cos((((df["p0_track_Chi2Dof"]) * (np.cos((df["p0_track_Chi2Dof"]))))))) * (np.log((df["IP_p0p2"])))))))), ((np.cos((df["p0_track_Chi2Dof"])))))) * (df["p0_track_Chi2Dof"]))) # df['line5'] = 1.0*np.tanh((((((((((p1)) / (((((p1)) + (((((df["SPDhits"]) / 2.0)) / 2.0)))/2.0)))) - (df["IP"]))) - (((df["SPDhits"]) / (df["p1_pt"]))))) * 2.0)) # df['line6'] = 1.0*np.tanh((((((((((((((df["CDF3"]) / (df["dira"]))) > (df["CDF3"]))*1.)) > (df["CDF3"]))*1.)) / 2.0)) + ((-1.0*((((((df["CDF3"]) * (df["p2_track_Chi2Dof"]))) * (((df["CDF3"]) * (df["p2_track_Chi2Dof"])))))))))/2.0)) # df['line7'] = 1.0*np.tanh((((-1.0*((((df["DOCAthree"]) / (df["CDF2"])))))) + (np.minimum(((((df["p2_pt"]) / (df["p0_p"])))), ((np.minimum(((df["CDF2"])), ((((np.sin((p3))) / 2.0)))))))))) # df['line8'] = 1.0*np.tanh(np.minimum((((-1.0*(((((((df["FlightDistance"]) < (df["IPSig"]))*1.)) / 2.0)))))), ((((np.minimum(((np.cos((np.log((df["p0_pt"])))))), ((np.cos((df["p1_track_Chi2Dof"])))))) / (p2)))))) # df['line9'] = 1.0*np.tanh(((np.sin((np.where(df["iso"]>0, ((((df["iso"]) - ((-1.0*((((df["IPSig"]) / 2.0))))))) / 2.0), ((((3.0) * (df["IP"]))) * 2.0) )))) / 2.0)) # df['line10'] = 1.0*np.tanh(((((np.cos(((((df["ISO_SumBDT"]) + (p2))/2.0)))) - (np.sin((np.log((df["p1_eta"]))))))) - ((((((df["ISO_SumBDT"]) + (np.cos((df["p2_IsoBDT"]))))/2.0)) * ((((df["ISO_SumBDT"]) + (np.cos((df["p2_IsoBDT"]))))/2.0)))))) return df #%% print("Add features") train = add_features(train) test = add_features(test) corr_check = add_features(corr_check) check_agreement = add_features(check_agreement) check_agreement = check_agreement.reset_index(drop=True) #%% print("Eliminate features") filter_out = ['id', 'min_ANNmuon', 'production', 'mass', 'signal' ,'isolationa', 'isolationb', 'isolationc', 'isolationd', 'isolatione', 'isolationf' ,'p0_IP', 'p1_IP', 'p2_IP' ,'p0_p', 'p1_p', 'p2_p' ,'IP_p0p2', 'IP_p1p2' ,'p0_IPSig', 'p1_IPSig', 'p2_IPSig' ,'SPDhits','CDF1', 'CDF2', 'CDF3' ,'p0_pt', 'p1_pt', 'p2_pt' ,'p0_eta', 'p1_eta', 'p2_eta' ,'p0_IsoBDT', 'p1_IsoBDT', 'p2_IsoBDT' ,'p0_track_Chi2Dof', 'p1_track_Chi2Dof', 'p2_track_Chi2Dof' ,'DOCAone', 'DOCAtwo', 'DOCAthree' # ## Kinematic features related to the mother particle: # ,'E', 'pz' # , 'beta', 'gamma', 'beta_gamma' , 'M_lt', 'M_inv' , 'Delta_E', 'Delta_M' ,'flag_M' # ## Kinematic features related to the final-state particles: # ,'E0', 'E1', 'E2', 'E0_ratio', 'E1_ratio', 'E2_ratio', 'p0_pt_ratio', 'p1_pt_ratio', 'p2_pt_ratio', 'eta_01', 'eta_02', 'eta_12', 't_coll' # ## Other extra features: # ,'FlightDistanceSig', 'DOCA_sum', 'isolation_sum', 'IsoBDT_sum', 'track_Chi2Dof', 'IP_sum', 'IPSig_sum', 'CDF_sum' ] #features = list(f for f in train.columns if f not in filter_out) features = ['LifeTime', 'dira', 'FlightDistance', 'FlightDistanceError', 'IP', 'IPSig', 'VertexChi2', 'pt', 'iso', 'ISO_SumBDT', 'E', 'pz', 'beta', 'gamma', 'beta_gamma'] features_mlpc = ['LifeTime', 'dira', 'FlightDistance', 'FlightDistanceError', 'IP' , 'IPSig', 'VertexChi2', 'pt', 'iso', 'ISO_SumBDT', 'NEW_FD_SUMP' , 'NEW5_lt', 'p_track_Chi2Dof_MAX', 'flight_dist_sig2', 'flight_dist_sig' , 'NEW_IP_dira', 'p0p2_ip_ratio', 'p1p2_ip_ratio', 'DCA_MAX' , 'iso_bdt_min', 'iso_min', 'NEW_iso_abc', 'NEW_iso_def', 'NEW_pN_IP' , 'NEW_pN_p', 'NEW_IP_pNpN', 'NEW_pN_IPSig', 'NEW_FD_LT'] # mjh filter on metric importances, improve auc or reduce mass correlation filter_imp = ['NEW_iso_def','NEW_iso_abc','NEW_pN_IP','iso_min','NEW_IP_pNpN', 'FlightDistance','flight_dist_sig'] #features = list(f for f in features if f not in filter_imp) # get validation data train_val, test_val, y_train, y_test, train_id, test_id = train_test_split(train, signal, trainids, random_state=100, test_size=0.15, shuffle=True) print('Number of features = ', len(features)) print('Features: ', features) print('...') #%% print("Train a Basic") start_time = time.time() # loss = BinFlatnessLossFunction(['mass'], n_bins=15, uniform_label=1 , fl_coefficient=10000, power=2) #loss = LogLossFunction() loss_bin = BinFlatnessLossFunction(['mass'], n_bins=15, uniform_label=1 , fl_coefficient=100, power=2) #n_estimators=500, ugbc_bin = UGradientBoostingClassifier(loss=loss_bin, n_estimators=5, max_depth=6, learning_rate=0.15, train_features=features, subsample=0.7, random_state=123) ugbc_bin.fit(train_val[features + ['mass']], train_val['signal']) pred_train_basic = ugbc_bin.predict_proba(train[features])[:, 1] pred_test_basic = ugbc_bin.predict_proba(test[features])[:, 1] pred_test_val_basic = ugbc_bin.predict_proba(test_val[features])[:, 1] pred_corr_check_basic = ugbc_bin.predict_proba(corr_check[features])[:, 1] agreement_probs_basic = ugbc_bin.predict_proba(check_agreement[features])[:, 1] ### number of signals n_sig_train=train['signal'].value_counts()[1] ### cutting 20% according to metric alpha = 0.8 beta = 1 n_sig_train = math.ceil(alpha*n_sig_train) ### threshold probability of signal pred_train_sort = np.flip(np.sort(pred_train_basic), axis=0) train_threshold = pred_train_sort[n_sig_train] train_threshold_cut = 1 - beta*(1-train_threshold) pred_train_shuffle_basic = pred_train_basic.copy() for ind, element in enumerate(pred_train_shuffle_basic.tolist()): if element=train_threshold_cut: # pred_test_shuffle[ind] = 1-(1-train_threshold_cut)*random.random() pred_test_shuffle_basic = pred_test_basic.copy() for ind, element in enumerate(pred_test_shuffle_basic.tolist()): if element=train_threshold_cut: # pred_test_shuffle[ind] = 1-(1-train_threshold_cut)*random.random() pred_test_val_shuffle_basic = pred_test_val_basic.copy() for ind, element in enumerate(pred_test_val_shuffle_basic.tolist()): if element=train_threshold_cut: # pred_test_val_shuffle[ind] = 1-(1-train_threshold_cut)*random.random() pred_corr_check_shuffle_basic = pred_corr_check_basic.copy() for ind, element in enumerate(pred_corr_check_shuffle_basic.tolist()): if element=train_threshold_cut: # pred_corr_check_shuffle[ind] = 1-(1-train_threshold_cut)*random.random() agreement_probs_shuffle_basic = agreement_probs_basic.copy() for ind, element in enumerate(agreement_probs_shuffle_basic.tolist()): if element=train_threshold_cut: # agreement_probs_shuffle[ind] = 1-(1-train_threshold_cut)*random.random() #%% features = ['GP'] print("Train a GP") start_time = time.time() # loss = BinFlatnessLossFunction(['mass'], n_bins=15, uniform_label=1 , fl_coefficient=10000, power=2) #loss = LogLossFunction() loss=CompositeLossFunction() #n_estimators=1000, ugbc_bin = UGradientBoostingClassifier(loss=loss, n_estimators=10, max_depth=6, learning_rate=0.15, train_features=features, subsample=0.7, random_state=123) ugbc_bin.fit(train_val[features], train_val['signal']) pred_train_GP = ugbc_bin.predict_proba(train[features])[:, 1] pred_test_GP = ugbc_bin.predict_proba(test[features])[:, 1] pred_test_val_GP = ugbc_bin.predict_proba(test_val[features])[:, 1] pred_corr_check_GP = ugbc_bin.predict_proba(corr_check[features])[:, 1] agreement_probs_GP = ugbc_bin.predict_proba(check_agreement[features])[:, 1] ### number of signals n_sig_train=train['signal'].value_counts()[1] ### cutting 20% according to metric alpha = 0.8 beta = 1 n_sig_train = math.ceil(alpha*n_sig_train) ### threshold probability of signal pred_train_sort = np.flip(np.sort(pred_train_GP), axis=0) train_threshold = pred_train_sort[n_sig_train] train_threshold_cut = 1 - beta*(1-train_threshold) pred_train_shuffle_GP = pred_train_GP.copy() for ind, element in enumerate(pred_train_shuffle_GP.tolist()): if element=train_threshold_cut: # pred_test_shuffle[ind] = 1-(1-train_threshold_cut)*random.random() pred_test_shuffle_GP = pred_test_GP.copy() for ind, element in enumerate(pred_test_shuffle_GP.tolist()): if element=train_threshold_cut: # pred_test_shuffle[ind] = 1-(1-train_threshold_cut)*random.random() pred_test_val_shuffle_GP = pred_test_val_GP.copy() for ind, element in enumerate(pred_test_val_shuffle_GP.tolist()): if element=train_threshold_cut: # pred_test_val_shuffle[ind] = 1-(1-train_threshold_cut)*random.random() pred_corr_check_shuffle_GP = pred_corr_check_GP.copy() for ind, element in enumerate(pred_corr_check_shuffle_GP.tolist()): if element=train_threshold_cut: # pred_corr_check_shuffle[ind] = 1-(1-train_threshold_cut)*random.random() agreement_probs_shuffle_GP = agreement_probs_GP.copy() for ind, element in enumerate(agreement_probs_shuffle_GP.tolist()): if element=train_threshold_cut: # agreement_probs_shuffle[ind] = 1-(1-train_threshold_cut)*random.random() #%% features = ['flag_M'] print("Train a FLAG_M") start_time = time.time() # loss = BinFlatnessLossFunction(['mass'], n_bins=15, uniform_label=1 , fl_coefficient=10000, power=2) #loss = LogLossFunction() loss=CompositeLossFunction() # n_estimators=1000, ugbc_bin = UGradientBoostingClassifier(loss=loss, n_estimators=10, max_depth=6, learning_rate=0.15, train_features=features, subsample=0.7, random_state=123) ugbc_bin.fit(train_val[features], train_val['signal']) pred_train_Mlt = ugbc_bin.predict_proba(train[features])[:, 1] pred_test_Mlt = ugbc_bin.predict_proba(test[features])[:, 1] pred_test_val_Mlt = ugbc_bin.predict_proba(test_val[features])[:, 1] pred_corr_check_Mlt = ugbc_bin.predict_proba(corr_check[features])[:, 1] agreement_probs_Mlt = ugbc_bin.predict_proba(check_agreement[features])[:, 1] ### number of signals n_sig_train=train['signal'].value_counts()[1] ### cutting 20% according to metric alpha = 0.8 beta = 1 n_sig_train = math.ceil(alpha*n_sig_train) ### threshold probability of signal pred_train_sort = np.flip(np.sort(pred_train_Mlt), axis=0) train_threshold = pred_train_sort[n_sig_train] train_threshold_cut = 1 - beta*(1-train_threshold) pred_train_shuffle_Mlt = pred_train_Mlt.copy() for ind, element in enumerate(pred_train_shuffle_Mlt.tolist()): if element=train_threshold_cut: # pred_test_shuffle[ind] = 1-(1-train_threshold_cut)*random.random() pred_test_shuffle_Mlt = pred_test_Mlt.copy() for ind, element in enumerate(pred_test_shuffle_Mlt.tolist()): if element=train_threshold_cut: # pred_test_shuffle[ind] = 1-(1-train_threshold_cut)*random.random() pred_test_val_shuffle_Mlt = pred_test_val_Mlt.copy() for ind, element in enumerate(pred_test_val_shuffle_Mlt.tolist()): if element=train_threshold_cut: # pred_test_val_shuffle[ind] = 1-(1-train_threshold_cut)*random.random() pred_corr_check_shuffle_Mlt = pred_corr_check_Mlt.copy() for ind, element in enumerate(pred_corr_check_shuffle_Mlt.tolist()): if element=train_threshold_cut: # pred_corr_check_shuffle[ind] = 1-(1-train_threshold_cut)*random.random() agreement_probs_shuffle_Mlt = agreement_probs_Mlt.copy() for ind, element in enumerate(agreement_probs_shuffle_Mlt.tolist()): if element=train_threshold_cut: # agreement_probs_shuffle[ind] = 1-(1-train_threshold_cut)*random.random() #%% print("Train a MPLC") start_time = time.time() clf = MLPClassifier(layers=[12,12], epochs=600) clf.fit(train[features_mlpc], train['signal']) pred_train_MLPC = clf.predict_proba(train[features_mlpc])[:, 1] pred_test_MLPC = clf.predict_proba(test[features_mlpc])[:, 1] pred_test_val_MLPC = clf.predict_proba(test_val[features_mlpc])[:, 1] pred_corr_check_MLPC = clf.predict_proba(corr_check[features_mlpc])[:, 1] agreement_probs_MLPC = clf.predict_proba(check_agreement[features_mlpc])[:, 1] ### number of signals n_sig_train=train['signal'].value_counts()[1] ### cutting 20% according to metric alpha = 0.8 beta = 1 n_sig_train = math.ceil(alpha*n_sig_train) ### threshold probability of signal pred_train_sort = np.flip(np.sort(pred_train_MLPC), axis=0) train_threshold = pred_train_sort[n_sig_train] train_threshold_cut = 1 - beta*(1-train_threshold) pred_train_shuffle_MLPC = pred_train_MLPC.copy() for ind, element in enumerate(pred_train_shuffle_MLPC.tolist()): if element