{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom sklearn.metrics import roc_curve, auc\nimport matplotlib.pyplot as plt\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def __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 np.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 = np.arange(1, total_events + 1, dtype='float') / total_events\n    subarray_distribution = np.cumsum(np.bincount(subindices, minlength=total_events), dtype='float')\n    subarray_distribution /= 1.0 * subarray_distribution[-1]\n    return np.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 = np.array(predictions)\n    masses = np.array(masses)\n    assert len(predictions) == len(masses)\n\n    # First, reorder by masses\n    predictions = predictions[np.argsort(masses)]\n\n    # Second, replace probabilities with order of probability among other events\n    predictions = np.argsort(np.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 np.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 = np.concatenate([sample_weights_zero, sample_weights_one])\n    data_all = np.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 = np.array(data_prediction), np.array(mc_prediction)\n    weights_data, weights_mc = np.array(weights_data), np.array(weights_mc)\n\n    assert np.all(data_prediction >= 0.) and np.all(data_prediction <= 1.), 'Data predictions are out of range [0, 1]'\n    assert np.all(mc_prediction >= 0.) and np.all(mc_prediction <= 1.), 'MC predictions are out of range [0, 1]'\n\n    weights_data /= np.sum(weights_data)\n    weights_mc /= np.sum(weights_mc)\n\n    fpr, tpr = __roc_curve_splitted(data_prediction, mc_prediction, weights_data, weights_mc)\n\n    Dnm = np.max(np.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 np.all(predictions >= 0.) and np.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 = np.minimum(tpr, tpr_thresholds[index])\n        tpr_previous = np.minimum(tpr, tpr_thresholds[index - 1])\n        area += roc_weights[index - 1] * (auc(fpr, tpr_cut) - auc(fpr, tpr_previous))\n    tpr_thresholds = np.array(tpr_thresholds)\n    # roc auc normalization to be 1 for an ideal classifier\n    area /= np.sum((tpr_thresholds[1:] - tpr_thresholds[:-1]) * np.array(roc_weights))\n    return area","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ls ../input/flavours-of-physics-kernels-only","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Load the train/test/eval data using pandas\")\ntrain = pd.read_csv(\"../input/flavours-of-physics-kernels-only/training.csv.zip\")\ntest  = pd.read_csv(\"../input/flavours-of-physics-kernels-only/test.csv.zip\")","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":true},"cell_type":"code","source":"# -------------- loading data files -------------- #\n\n\n\n#--------------- 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\n\nprint(\"Add features\")\ntrain = add_features(train)\ntest = add_features(test)\n\n\nprint(\"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\nfeatures = list(f for f in train.columns if f not in filter_out)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mungedtrain = train[features].astype(complex)\nmungedtest = test[features].astype(complex)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def Output(p):\n    return 1./(1.+np.exp(-(0.476446 +p)))\ndef GPReal(data):\n     return (0.659360*np.tanh(np.real(((np.sin((np.sin(((((((data[\"p_track_Chi2Dof_MAX\"]) * (data[\"p_track_Chi2Dof_MAX\"]))) + (data[\"IPSig\"]))/2.0)))))) - ((((((data[\"p_track_Chi2Dof_MAX\"]) * (np.power((data[\"iso\"]),2.)))) + ((((((data[\"IPSig\"]) + (data[\"VertexChi2\"]))/2.0)) - (((data[\"flight_dist_sig\"]) / (data[\"IPSig\"]))))))/2.0))))) +\n             0.800000*np.tanh(np.real(((np.cos((((np.power((((((data[\"LifeTime\"]) / (((np.tanh((complex(1.86195778846740723)))) * (np.conjugate(data[\"NEW_FD_SUMP\"]))*(complex(0,1)))))) - (((((data[\"dira\"]) + (data[\"LifeTime\"]))) + (data[\"LifeTime\"]))))),2.)) / 2.0)))) - (((data[\"VertexChi2\"]) + (np.power((((np.tanh((data[\"LifeTime\"]))) / (data[\"NEW_FD_SUMP\"]))),3.))))))))\ndef GPComplex(data):\n    return (0.659360*np.tanh(np.imag(((np.sin((np.sin(((((((data[\"p_track_Chi2Dof_MAX\"]) * (data[\"p_track_Chi2Dof_MAX\"]))) + (data[\"IPSig\"]))/2.0)))))) - ((((((data[\"p_track_Chi2Dof_MAX\"]) * (np.power((data[\"iso\"]),2.)))) + ((((((data[\"IPSig\"]) + (data[\"VertexChi2\"]))/2.0)) - (((data[\"flight_dist_sig\"]) / (data[\"IPSig\"]))))))/2.0))))) +\n             0.800000*np.tanh(np.imag(((np.cos((((np.power((((((data[\"LifeTime\"]) / (((np.tanh((complex(1.86195778846740723)))) * (np.conjugate(data[\"NEW_FD_SUMP\"]))*(complex(0,1)))))) - (((((data[\"dira\"]) + (data[\"LifeTime\"]))) + (data[\"LifeTime\"]))))),2.)) / 2.0)))) - (((data[\"VertexChi2\"]) + (np.power((((np.tanh((data[\"LifeTime\"]))) / (data[\"NEW_FD_SUMP\"]))),3.))))))))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15,15))\ncolors=['r','g']\nplt.scatter(GPReal(mungedtrain),GPComplex(mungedtrain),s=1,color=[colors[int(c)] for c in train.signal.values])\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15,15))\nplt.scatter(GPReal(mungedtest[::12]),GPComplex(mungedtest[::12]),s=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"check_agreement = pd.read_csv('../input/flavours-of-physics-kernels-only/check_agreement.csv.zip')\ncheck_correlation = pd.read_csv('../input/flavours-of-physics-kernels-only/check_correlation.csv.zip')\ncheck_agreement.head()\ncheck_agreement = add_features(check_agreement)\ncheck_agreement = check_agreement[features+['weight','signal']]\ncheck_agreement.head()\ncheck_correlation = add_features(check_correlation)\ncheck_correlation = check_correlation[features+['mass']]\ncheck_agreement[check_agreement.columns[:-2]] = check_agreement[check_agreement.columns[:-2]].astype(complex)\nagreement_probs = Output(GPReal(check_agreement))\nprint('Checking agreement...')\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 metric', ks, ks < 0.09)\ncheck_correlation[check_correlation.columns[:-1]] = check_correlation[check_correlation.columns[:-1]].astype(complex)\ncorrelation_probs = Output(GPReal(check_correlation))\nprint ('Checking correlation...')\ncvm = compute_cvm(correlation_probs, check_correlation['mass'])\nprint ('CvM metric', cvm, cvm < 0.002)\n\ntrain_eval_probs = Output(GPReal(mungedtrain))\n\nprint ('Calculating AUC...')\nAUC = roc_auc_truncated(train['signal'], train_eval_probs)\nprint ('AUC', AUC)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.DataFrame({'id':test.id.values,'prediction':Output(GPReal(mungedtest))}).to_csv('submission.csv',index=False)\n","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.3"}},"nbformat":4,"nbformat_minor":1}