{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Outline"},{"metadata":{},"cell_type":"markdown","source":"**Objective** - To search Physics beyond Standard Model. <br><br>\nStandard Model was complete with the discovery of Higgs Boson. But that model is still incomplete, as it is unable to explain certain phenomenon like Dark Matter, gravity etc.. <br><br>\n\n**Lepton Flavor** (e<sup>-</sup>, u<sup>-</sup>, t<sup>-</sup>, and the corresponding neutrinos) is a conserved quantity in Standard Model. But the rare decay addressed in this particular problem t->3u violates it.<br><br>\n\nThis problem is focussed on finding this rare decay, which could pave the way for Physics beyond Standard Model"},{"metadata":{},"cell_type":"markdown","source":"The dataset is divided among following files-\n\n* **training.csv** - It is the dataset consisting of real data and simulated data\n* **test.csv** - Unseen data which is to be used to check our results\n* **check_agreement.csv** -     data with a control channel, Ds → φ(2u)π, which has a similar topology as of t(tau)->3u(muons). This test is necessary to avoid bias of classifier on the simulated data\n* **check_correlation.csv** -   this is to check whether classifier is not too correlated with the mass, so as to avoid false peaks. Also as shown later, that mass will separate background and signal almost perfectly"},{"metadata":{},"cell_type":"markdown","source":"# Importing Libraries"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy\nimport pandas as pd\nfrom sklearn.model_selection import train_test_split\nimport xgboost as xgb\nfrom xgboost import plot_importance\nimport matplotlib.pyplot as plt\nfrom hep_ml.gradientboosting import UGradientBoostingClassifier\nfrom hep_ml.losses import BinFlatnessLossFunction","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Basic Utility functions**"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"from sklearn.metrics import roc_curve, auc\n\n\ndef __rolling_window(data, window_size):\n    \"\"\"\n    Rolling window: take window with definite size through the array\n    :param data: array-like\n    :param window_size: size\n    :return: the sequence of windows\n    Example: data = array(1, 2, 3, 4, 5, 6), window_size = 4\n        Then this function return array(array(1, 2, 3, 4), array(2, 3, 4, 5), array(3, 4, 5, 6))\n    \"\"\"\n    shape = data.shape[:-1] + (data.shape[-1] - window_size + 1, window_size)\n    strides = data.strides + (data.strides[-1],)\n    return numpy.lib.stride_tricks.as_strided(data, shape=shape, strides=strides)\n\n\ndef __cvm(subindices, total_events):\n    \"\"\"\n    Compute Cramer-von Mises metric.\n    Compared two distributions, where first is subset of second one.\n    Assuming that second is ordered by ascending\n    :param subindices: indices of events which will be associated with the first distribution\n    :param total_events: count of events in the second distribution\n    :return: cvm metric\n    \"\"\"\n    target_distribution = numpy.arange(1, total_events + 1, dtype='float') / total_events\n    subarray_distribution = numpy.cumsum(numpy.bincount(subindices, minlength=total_events), dtype='float')\n    subarray_distribution /= 1.0 * subarray_distribution[-1]\n    return numpy.mean((target_distribution - subarray_distribution) ** 2)\n\n\ndef compute_cvm(predictions, masses, n_neighbours=200, step=50):\n    \"\"\"\n    Computing Cramer-von Mises (cvm) metric on background events: take average of cvms calculated for each mass bin.\n    In each mass bin global prediction's cdf is compared to prediction's cdf in mass bin.\n    :param predictions: array-like, predictions\n    :param masses: array-like, in case of Kaggle tau23mu this is reconstructed mass\n    :param n_neighbours: count of neighbours for event to define mass bin\n    :param step: step through sorted mass-array to define next center of bin\n    :return: average cvm value\n    \"\"\"\n    predictions = numpy.array(predictions)\n    masses = numpy.array(masses)\n    assert len(predictions) == len(masses)\n\n    # First, reorder by masses\n    predictions = predictions[numpy.argsort(masses)]\n\n    # Second, replace probabilities with order of probability among other events\n    predictions = numpy.argsort(numpy.argsort(predictions, kind='mergesort'), kind='mergesort')\n\n    # Now, each window forms a group, and we can compute contribution of each group to CvM\n    cvms = []\n    for window in __rolling_window(predictions, window_size=n_neighbours)[::step]:\n        cvms.append(__cvm(subindices=window, total_events=len(predictions)))\n    return numpy.mean(cvms)\n\n\ndef __roc_curve_splitted(data_zero, data_one, sample_weights_zero, sample_weights_one):\n    \"\"\"\n    Compute roc curve\n    :param data_zero: 0-labeled data\n    :param data_one:  1-labeled data\n    :param sample_weights_zero: weights for 0-labeled data\n    :param sample_weights_one:  weights for 1-labeled data\n    :return: roc curve\n    \"\"\"\n    labels = [0] * len(data_zero) + [1] * len(data_one)\n    weights = numpy.concatenate([sample_weights_zero, sample_weights_one])\n    data_all = numpy.concatenate([data_zero, data_one])\n    fpr, tpr, _ = roc_curve(labels, data_all, sample_weight=weights)\n    return fpr, tpr\n\n\ndef compute_ks(data_prediction, mc_prediction, weights_data, weights_mc):\n    \"\"\"\n    Compute Kolmogorov-Smirnov (ks) distance between real data predictions cdf and Monte Carlo one.\n    :param data_prediction: array-like, real data predictions\n    :param mc_prediction: array-like, Monte Carlo data predictions\n    :param weights_data: array-like, real data weights\n    :param weights_mc: array-like, Monte Carlo weights\n    :return: ks value\n    \"\"\"\n    assert len(data_prediction) == len(weights_data), 'Data length and weight one must be the same'\n    assert len(mc_prediction) == len(weights_mc), 'Data length and weight one must be the same'\n\n    data_prediction, mc_prediction = numpy.array(data_prediction), numpy.array(mc_prediction)\n    weights_data, weights_mc = numpy.array(weights_data), numpy.array(weights_mc)\n\n    assert numpy.all(data_prediction >= 0.) and numpy.all(data_prediction <= 1.), 'Data predictions are out of range [0, 1]'\n    assert numpy.all(mc_prediction >= 0.) and numpy.all(mc_prediction <= 1.), 'MC predictions are out of range [0, 1]'\n\n    weights_data /= numpy.sum(weights_data)\n    weights_mc /= numpy.sum(weights_mc)\n\n    fpr, tpr = __roc_curve_splitted(data_prediction, mc_prediction, weights_data, weights_mc)\n\n    Dnm = numpy.max(numpy.abs(fpr - tpr))\n    return Dnm\n\n\ndef roc_auc_truncated(labels, predictions, tpr_thresholds=(0.2, 0.4, 0.6, 0.8),\n                      roc_weights=(4, 3, 2, 1, 0)):\n    \"\"\"\n    Compute weighted area under ROC curve.\n    :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":{},"cell_type":"markdown","source":"# M3Body Event Selection"},{"metadata":{},"cell_type":"markdown","source":"**Event Selection** <br><br>\n\nTau particles are the heavier brothers of electron, with mass of aroung 1.78GeV/c<sup>2</sup>. M3Body selection defines the geometric and kinematic properties of the 3 track system which helps in identifying the rare decay. Some good tips regarding event selection are mentioned in the [paper](https://arxiv.org/pdf/1409.8548.pdf)\n\n1. Tau leptons usually originate from a heavy quark decay, which are formed via proton-proton collision.\n2. Candidates are selected via 3 tracks, which form a vertex with particle mass (around 1780Mev/c<sup>2</sup> aka mass of tau particle), and also which is displaced siginificatnly from PV\n\nThe point where proton-proton collision occurs is called Primary Vertex(PV) and the point where tau particle will decay into 3 particles is called secondary vertex. Since the 3 particle decay happens at Secondary vertex, there must be a significant distance between tracks and the PV.\n\n3. Due to low Q-values the angle between momentum vector and collision-decay vertex is small, and is almost collinear. So dira will be having quite high value\n\nThe paper uses 2 classifiers M3Body(using the geometric properties to form a vertex of 3 tracks), and MPID(particle identification of Muon). But here in kaggle we are focussed on M3Body classification only. Instead of using MPID, we are given feature 'minANN_muon' which provides us with probability of particle being a Muon\n\n\n"},{"metadata":{},"cell_type":"markdown","source":"# Data Analysis"},{"metadata":{"trusted":true},"cell_type":"code","source":"train= pd.read_csv(\"../input/flavours-of-physics-kernels-only/training.csv\")\nsignal = train[train.signal==1]\nbackground = train[train.signal==0]\nplt.rcParams[\"figure.figsize\"] = (8, 5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Dira is the angle between momentum vector of 3 particle system, and the line joining PV(p-p collision vertex) and DV(Decay Vertex). Its value is almost collinear due to low loss during decay"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(signal.dira, range=(0.997,1), bins=100, label='Signal', alpha=0.7,normed=True)\nplt.hist(background.dira, range=(0.997,1), bins=100, label='Background', alpha=0.7,normed=True)\nplt.xlabel('Dira')\nplt.legend()\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"VertexChiSquare signifies the vertex fit quality from the three particle decay body. Lower the value, better the fit. Since t->3u decay has a better vertex fit, it is a high-valued feature in our classifier"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(signal.VertexChi2, range=(0,16), bins=100, label='Signal', alpha=0.7,normed=True)\nplt.hist(background.VertexChi2,range=(0,16),bins=100, label='Background',alpha=0.7,normed=True)\nplt.xlabel(r'Vertex $\\chi^2$')\nplt.ylabel('Normalized Fraction')\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"IPSig seems to be a highly valued feature. It signifies the closest distance between PV and the track"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(signal.IPSig, range=(0,25), bins=100, label='Signal', alpha=0.7,normed=True)\nplt.hist(background.IPSig,range=(0,25),bins=100, label='Background',alpha=0.7,normed=True)\nplt.xlabel('IPSig')\nplt.ylabel('Normalized Fraction')\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Mass is a feature which can almost perfectly distinguish between signal and background. Signal is a peak concenterated around the mass of tau particle(1780Mev/c<sup>2</sup>).<br><br>\nAnd background is constructed by removing the actual data around the mass region. These observations are clearly visible in the below mentioned graphs"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.rcParams[\"figure.figsize\"] = (10, 3)\nfig, (ax1,ax2)= plt.subplots(1,2,sharey=False)\nax1.legend(\"Background\")\ntrain['mass'][train['signal']==0].hist(bins=100,ax=ax1)\nax1.plot()\nax2.legend(\"Signal\")\ntrain['mass'][train['signal']==1].hist(bins=100,ax=ax2)\nax2.plot()\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"*min_ANNmuon* This feature is the result of using ANN for particle identification. Here it is completely visible that high probability of this feature correclty coreesponds to the high signal. <br><br>\nThis feature is not to be used in the classifier training"},{"metadata":{"trusted":true},"cell_type":"code","source":"lowMuonSignal= train[train['min_ANNmuon']<=0.4]['signal'].sum() / train[train['min_ANNmuon']<=0.4].shape[0]\nhighMuonSignal= train[train['min_ANNmuon']>0.4]['signal'].sum() / train[train['min_ANNmuon']>0.4].shape[0]\n\nprint(f\"The signal-to-background ratio having muon probability less than 0.4 is {lowMuonSignal*100}%\")\nprint(f\"The signal-to-background ratio having muon probability higher than 0.4 is {highMuonSignal*100}%\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Validation Set"},{"metadata":{},"cell_type":"markdown","source":"As recommended, *min_ANNmuon>0.4* is to be used for evaluating score. So we have created validation set satisfying the test set criteria "},{"metadata":{"trusted":true},"cell_type":"code","source":"def getTrnValPair(train_data):\n    tr,vl=train_test_split( train_data[train_data['min_ANNmuon']>0.4], test_size=0.10, random_state=42)\n    training_set= pd.concat([ train_data[train_data['min_ANNmuon']<=0.4],tr ])\n    return training_set,vl","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Kinematic Features"},{"metadata":{"trusted":true},"cell_type":"code","source":"def addKinematicFeatures(df):\n    p0_pz = (df.p0_p**2 - df.p0_pt**2)**0.5\n    p1_pz = (df.p1_p**2 - df.p1_pt**2)**0.5\n    p2_pz = (df.p2_p**2 - df.p2_pt**2)**0.5\n    df['pz'] = p0_pz + p1_pz + p2_pz\n    df['p'] = (df.pt**2 + df.pz**2)**0.5\n    df['NEW_FD_SUMP']=df['FlightDistance']/(df['p0_p']+df['p1_p']+df['p2_p'])\n    df['vel']= (df['FlightDistance'])/df['LifeTime']\n    return df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Geometric Features"},{"metadata":{"trusted":true},"cell_type":"code","source":"def addGemetricFeatures(df):\n    df['Flight_Dist_Sig']= df['FlightDistance']/df['FlightDistanceError']\n    df['iso_diff']=((df['isolationa']+df['isolationb']+df['isolationc'])/3+\n                 (df['isolationd']+df['isolatione']+df['isolationf'])/3)/2-df['iso']\n    df['DOC_max']= df.loc[:,['DOCAone','DOCAtwo','DOCAthree']].max(axis=1)\n    df['iso_bdt_p_min']= df.loc[:,['p0_IsoBDT','p1_IsoBDT','p2_IsoBDT']].min(axis=1)\n    df['p_track_Chi2Dof_MAX'] = df.loc[:, ['p0_track_Chi2Dof', 'p1_track_Chi2Dof', 'p2_track_Chi2Dof']].max(axis=1)\n    df['CDF_sum'] = df['CDF1']+df['CDF2']+df['CDF3']\n    \n    df['NEW5_lt']=df['LifeTime']*(df['p0_IP']+df['p1_IP']+df['p2_IP'])/3\n    return df\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"drop_cols= ['min_ANNmuon','production','signal','mass','SPDhits','FlightDistance','id',\n            'isolationa','isolationb','isolationc','isolationd','isolatione','isolationf',\n           'DOCAone', 'DOCAtwo', 'DOCAthree','p0_IsoBDT','p1_IsoBDT','p2_IsoBDT',\n           'p0_track_Chi2Dof', 'p1_track_Chi2Dof', 'p2_track_Chi2Dof','FlightDistanceError',\n           'CDF2','CDF3','CDF1','p0_eta', 'p1_eta', 'p2_eta','LifeTime']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# XGBoost"},{"metadata":{"trusted":true},"cell_type":"code","source":"def fitXGBClassifier(training_set):\n    clfXGB = xgb.XGBClassifier(learning_rate =0.1, n_estimators=200, max_depth=3,random_state = 42,n_jobs=-1)\n    clfXGB.fit(training_set.drop(drop_cols,axis=1), training_set['signal'])\n    return clfXGB\n\ndef xgbImportance(clfXGB):\n    plt.rcParams[\"figure.figsize\"] = (15, 30)\n    plot_importance(clfXGB)\n    plt.show()\n\ndef getXGBValScore(clfXGB,validation_set):\n    train_probs= clfXGB.predict_proba(validation_set.drop(drop_cols,axis=1))[:, 1]\n    AUC = roc_auc_truncated(validation_set['signal'], train_probs)\n    return AUC","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ntrn,val= getTrnValPair(train)\nclfXGB= fitXGBClassifier(trn)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"xgbImportance(clfXGB)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Uniform boost"},{"metadata":{},"cell_type":"markdown","source":"This is a boost algorithm based on [paper](https://arxiv.org/pdf/1410.4140.pdf) where we have customized our loss function, which maintains a trade-off between classifier accuracy and mass correlation<br><br>\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"removeUBoostFeatures= []\ndef fitUBoostClassifier(training_set):\n    tempTrn= training_set.drop(drop_cols+removeUBoostFeatures,axis=1)\n    loss = BinFlatnessLossFunction(['mass'], n_bins=15, fl_coefficient=0.1, uniform_label=0)\n    uBoost = UGradientBoostingClassifier(loss=loss, n_estimators=40, subsample=0.1, \n                                           max_depth=7, min_samples_leaf=10,\n                                           learning_rate=0.1, train_features=tempTrn.columns, random_state=11)\n    uBoost.fit(training_set, training_set['signal'])\n    return uBoost\n\ndef plotUBoostFeatureImportance(baseline,cols):\n    plt.rcParams[\"figure.figsize\"] = (15, 30)\n    features= baseline.feature_importances_\n    indices= numpy.argsort(features)\n    # Create plot\n    plt.figure()\n\n    # Create plot title\n    plt.title(\"Feature Importance\")\n\n    # Add bars\n    plt.bar(range(cols.shape[0]), features[indices])\n\n    # Add feature names as x-axis labels\n    plt.xticks(range(cols.shape[0]), cols[indices], rotation=90)\n\n    # Show plot\n    plt.show()\n    \ndef getValUBoostScore(uBoost,validation_set):\n    train_probs= uBoost.predict_proba(validation_set.drop(drop_cols+removeUBoostFeatures,axis=1))[:, 1]\n    AUC = roc_auc_truncated(validation_set['signal'], train_probs)\n    return AUC","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"%%time\ntrain= addGemetricFeatures(train);\ntrain= addKinematicFeatures(train);\ntrn,val= getTrnValPair(train);\nuBoost= fitUBoostClassifier(trn);\ngetValUBoostScore(uBoost,val)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plotUBoostFeatureImportance(uBoost,trn.drop(drop_cols+removeUBoostFeatures,axis=1).columns);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Agreement Test"},{"metadata":{"trusted":true},"cell_type":"code","source":"ks_drop_cols=[d for d in drop_cols if d not in ['min_ANNmuon','production','mass','id']]\ncheck_agreement = pd.read_csv('../input/flavours-of-physics-kernels-only/check_agreement.csv', index_col='id')\nagreement_probs = uBoost.predict_proba(check_agreement.drop(ks_drop_cols,axis=1))[:, 1]\n\nks = compute_ks(\n    agreement_probs[check_agreement['signal'].values == 0],\n    agreement_probs[check_agreement['signal'].values == 1],\n    check_agreement[check_agreement['signal'] == 0]['weight'].values,\n    check_agreement[check_agreement['signal'] == 1]['weight'].values)\nprint(ks,ks < 0.09)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Correlation Check"},{"metadata":{"trusted":true},"cell_type":"code","source":"cvm_drop_cols=[d for d in drop_cols if d not in ['min_ANNmuon','production','signal','id','mass']]\ncheck_correlation = pd.read_csv('../input/flavours-of-physics-kernels-only/check_correlation.csv', index_col='id')\ncorrelation_probs = uBoost.predict_proba(check_correlation.drop(cvm_drop_cols,axis=1))[:, 1]\ncvm = compute_cvm(correlation_probs, check_correlation['mass'])\nprint('CvM metric', cvm, cvm < 0.002)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Test Submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"test= pd.read_csv('../input/flavours-of-physics-kernels-only/test.csv')\ntest_probs= uBoost.predict_proba(test)[:, 1]\ntest['prediction']= test_probs\ntest[['id','prediction']].to_csv(\"submission.csv\",index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}