{"cells":[{"metadata":{"_uuid":"ca7afb91490ad9322a6904d146c26fea8b83d482"},"cell_type":"markdown","source":"Many thanks to Luis Andre Dutra e Silva's [kernel](https://www.kaggle.com/mindcool/unrolling-of-helices-outliers-removal)."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-output":true,"_kg_hide-input":true,"collapsed":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport tensorflow as tf\nfrom tqdm import tqdm\nimport glob\nimport os\n\nfrom trackml.dataset import load_event, load_dataset\nfrom trackml.score import score_event\nfrom IPython.display import display\n\nfrom keras.models import Sequential\nfrom keras.layers import Dense, Dropout, BatchNormalization\nfrom keras.layers.embeddings import Embedding","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"trusted":true,"_uuid":"800999c86031aa047322b979505538b775f3f66a","collapsed":true},"cell_type":"code","source":"train = np.unique([p.split('-')[0] for p in sorted(glob.glob('../input/train_*/**'))])\ntest = np.unique([p.split('-')[0] for p in sorted(glob.glob('../input/test/**'))])\n# det = pd.read_csv('../input/detectors.csv')\n# sub = pd.read_csv('../input/sample_submission.csv')\n\ndef read_files(path):\n    return [p.split('-')[0] for p in sorted(glob.glob(path))]\n\ntrain_file_names = np.unique(read_files('../input/train_*/**.csv'))\ntest_file_names = np.unique(read_files('../input/test/**.csv'))\ndetector = pd.read_csv('../input/detectors.csv')\nsubmission = pd.read_csv('../input/sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e1b2f7f4a809498a920e58d9f5ec3a5dffaf24d9"},"cell_type":"markdown","source":"# Preview of the Data"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-output":false,"_kg_hide-input":true,"scrolled":true,"collapsed":true},"cell_type":"code","source":"hits, cells, particles, truth = load_event(train[0])\nprint(len(hits), len(cells), len(particles), len(truth))\nn=10\nprint(\"hits:\")\ndisplay(hits.head(n))\nprint(\"cells:\")\ndisplay(cells.head(n))\nprint(\"particles:\")\ndisplay(particles.head(n))\nprint(\"truth:\")\ndisplay(truth.head(n))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d922a795dfbd5b6e51d445ed3296402f7c6e0100"},"cell_type":"markdown","source":"# Denoise"},{"metadata":{"_uuid":"34c0c5a09c26bf995d89349583825e5acdd26dfa"},"cell_type":"markdown","source":"### Create Noise Removal Model\n\nThis is just a simple NN, could be made more complex later."},{"metadata":{"trusted":true,"_uuid":"221f78b1a17bf08f8ccd81c8b9f8257547449f88","collapsed":true,"_kg_hide-output":true},"cell_type":"code","source":"denoise_model = Sequential()\ndenoise_model.add(Dense(256, activation='relu', input_shape=[25, ]))\ndenoise_model.add(Dense(256, activation='relu'))\ndenoise_model.add(Dense(256, activation='relu'))\ndenoise_model.add(Dropout(0.5))\ndenoise_model.add(Dense(256, activation='relu'))\ndenoise_model.add(Dense(256, activation='relu'))\ndenoise_model.add(Dense(256, activation='relu'))\ndenoise_model.add(BatchNormalization())\ndenoise_model.add(Dense(1, activation='sigmoid'))\n\ndenoise_model.compile(loss='binary_crossentropy',\n                      optimizer='adam',\n                      metrics=['accuracy'])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ef623f3c621ccbdbdd461ff0e068ea953dc05335"},"cell_type":"markdown","source":"### Train Noise Removal Model"},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true,"collapsed":true,"_uuid":"e60be9d83e7f16e6a77d8a72598bd7998c9ebe65"},"cell_type":"code","source":"def get_train_data_for_denoise(e):\n    # https://www.kaggle.com/meaninglesslives/classifier-hdbscan-helixfitting\n    hits, cells, truth = load_event(e, parts=['hits', 'cells', 'truth'])\n    hits['event_id'] = int(e[-9:])\n    cells = cells.groupby(by=['hit_id'])['ch0', 'ch1', 'value'].agg(\n        ['mean', 'median', 'max', 'min', 'sum']).reset_index()\n    cells.columns = ['hit_id'] + ['-'.join([c2, c1]) for c1 in [\n        'mean', 'median', 'max', 'min', 'sum']\n        for c2 in ['ch0', 'ch1', 'value']]\n    hits = pd.merge(hits, cells, how='left', on='hit_id')\n    tcols = list(truth.columns)\n    hits = pd.merge(hits, truth, how='left', on='hit_id')\n    hits = norm_points(hits)\n    cols = [c for c in hits.columns if c not in [\n        'event_id', 'hit_id', 'particle_id'] + tcols]\n\n    # noise marking\n    # noise -> target==0\n    # genuine data -> target==1\n    hits['target'] = hits['particle_id'].map(lambda x: 0 if x == 0 else 1)\n\n    # return x,y\n    return hits[cols], hits['target'].values\n\n\n# https://www.kaggle.com/mikhailhushchyn/dbscan-benchmark\n# https://www.kaggle.com/mikhailhushchyn/hough-transform\ndef norm_points(df, g=0.7):\n    x = df.x.values\n    y = df.y.values\n    z = df.z.values\n    r = np.sqrt(x**2 + y**2 + z**2)\n    df['x2'] = x / r\n    df['y2'] = y / r\n    r = np.sqrt(x**2 + y**2)\n    df['z2'] = (z / r) * g\n\n    df['r'] = r\n    # df['r2'] = np.sqrt(x**2 + y**2)\n    # df['r3'] = np.sqrt(x**2 + z**2)\n    # df['r4'] = np.sqrt(y**2 + z**2)\n    # df['phi'] = np.arctan2(y, x)\n    # df['phi2'] = np.arctan2(y, z)\n    # df['phi3'] = np.arctan2(x, z)\n    # df['hm'] = (2. * np.cos(df['phi'] - g) / df['r2']).values\n    return df\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"933e61777e2db7929b9cb0da1f1913246b5e7d7b","_kg_hide-output":false,"collapsed":true},"cell_type":"code","source":"# denoise_train_batch_size = 4096\n# denoise_train_epochs = 2\n\n# for f in tqdm(train_file_names[0:2]):\n#     X, y = get_train_data_for_denoise(f)\n#     denoise_model.fit(X, y, batch_size=denoise_train_batch_size, \n#                       epochs=denoise_train_epochs)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"aad33633ecfdfae2662a5206af24fc65a2a22689"},"cell_type":"markdown","source":"### Evaluate Denoise Model"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"b5709a77c106da0b546ea0b27b4753a5be03b424"},"cell_type":"code","source":"# use sklearn's train_test_split","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d770950dd126cac6572ddaab0cf51ce034bd458e"},"cell_type":"markdown","source":"### Create New Dataset With Noise Removed"},{"metadata":{"trusted":true,"_uuid":"521ee886e9143c28ce6cbb8d41107cc21927d453","collapsed":true},"cell_type":"code","source":"# proper_train_data = []\n# noise_train_data = []\n\n# for d in train_data:\n# #     noise==0, proper==1\n#     if denoise_model.predict(d) >= 0.5:\n#         proper_train_data.append(d)\n#     else:\n#         noise_train_data.append(d)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"338a763b12a55633088cd38ff14ebcdb540bee57"},"cell_type":"markdown","source":"# Clustering"},{"metadata":{"trusted":true,"_uuid":"cd1eebf913deab7075afb5e17103047f03bbca6c","_kg_hide-output":true,"collapsed":true},"cell_type":"code","source":"import xgboost as xgb\nfrom sklearn import preprocessing, model_selection\nimport hdbscan\n\nscl = preprocessing.StandardScaler()\n#dbscan = cluster.DBSCAN(eps=0.007555, min_samples=1, algorithm='kd_tree', n_jobs=-1)\ndbscan = hdbscan.HDBSCAN(min_samples=3, \n                         min_cluster_size=5, \n                         cluster_selection_method='leaf', \n                         prediction_data=False, \n                         metric='braycurtis')","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true,"collapsed":true,"_uuid":"2144ceff43b725ce43cd01eeed45085b9f6d0f68"},"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler\nimport hdbscan\nfrom scipy import stats\nfrom tqdm import tqdm\nfrom sklearn.cluster import DBSCAN\n\nclass Clusterer(object):\n    def __init__(self,rz_scales=[0.65, 0.965, 1.528]):                        \n        self.rz_scales=rz_scales\n    \n    def _eliminate_outliers(self,labels,M):\n        norms=np.zeros((len(labels)),np.float32)\n        indices=np.zeros((len(labels)),np.float32)\n        for i, cluster in tqdm(enumerate(labels),total=len(labels)):\n            if cluster == 0:\n                continue\n            index = np.argwhere(self.clusters==cluster)\n            index = np.reshape(index,(index.shape[0]))\n            indices[i] = len(index)\n            x = M[index]\n            norms[i] = self._test_quadric(x)\n        threshold1 = np.percentile(norms,90)*5\n        threshold2 = 25\n        threshold3 = 6\n        for i, cluster in enumerate(labels):\n            if norms[i] > threshold1 or indices[i] > threshold2 or indices[i] < threshold3:\n                self.clusters[self.clusters==cluster]=0   \n    def _test_quadric(self,x):\n        if x.size == 0 or len(x.shape)<2:\n            return 0\n        Z = np.zeros((x.shape[0],10), np.float32)\n        Z[:,0] = x[:,0]**2\n        Z[:,1] = 2*x[:,0]*x[:,1]\n        Z[:,2] = 2*x[:,0]*x[:,2]\n        Z[:,3] = 2*x[:,0]\n        Z[:,4] = x[:,1]**2\n        Z[:,5] = 2*x[:,1]*x[:,2]\n        Z[:,6] = 2*x[:,1]\n        Z[:,7] = x[:,2]**2\n        Z[:,8] = 2*x[:,2]\n        Z[:,9] = 1\n        v, s, t = np.linalg.svd(Z,full_matrices=False)        \n        smallest_index = np.argmin(np.array(s))\n        T = np.array(t)\n        T = T[smallest_index,:]        \n        norm = np.linalg.norm(np.dot(Z,T), ord=2)**2\n        return norm\n\n    def _preprocess(self, hits):\n        \n        x = hits.x.values\n        y = hits.y.values\n        z = hits.z.values\n\n        r = np.sqrt(x**2 + y**2 + z**2)\n        hits['x2'] = x/r\n        hits['y2'] = y/r\n\n        r = np.sqrt(x**2 + y**2)\n        hits['z2'] = z/r\n\n        ss = StandardScaler()\n        X = ss.fit_transform(hits[['x2', 'y2', 'z2']].values)\n        for i, rz_scale in enumerate(self.rz_scales):\n            X[:,i] = X[:,i] * rz_scale\n       \n        return X\n    \n    def _init(self,dfh):\n        dfh['r'] = np.sqrt(dfh['x'].values**2+dfh['y'].values**2+dfh['z'].values**2)\n        dfh['rt'] = np.sqrt(dfh['x'].values**2+dfh['y'].values**2)\n        dfh['a0'] = np.arctan2(dfh['y'].values,dfh['x'].values)\n        dfh['z1'] = dfh['z'].values/dfh['rt'].values\n        dfh['x2'] = 1/dfh['z1'].values\n        dz0 = -0.00070\n        stepdz = 0.00001\n        stepeps = 0.000005\n        mm = 1\n        for ii in tqdm(range(100)):\n            mm = mm*(-1)\n            dz = mm*(dz0+ii*stepdz)\n            dfh['a1'] = dfh['a0'].values+dz*dfh['z'].values*np.sign(dfh['z'].values)\n            dfh['sina1'] = np.sin(dfh['a1'].values)\n            dfh['cosa1'] = np.cos(dfh['a1'].values)\n            dfh['x1'] = dfh['a1'].values/dfh['z1'].values\n            ss = StandardScaler()\n            dfs = ss.fit_transform(dfh[['sina1','cosa1','z1','x1','x2']].values)\n            cx = np.array([1, 1, 0.75, 0.5, 0.5])\n            for k in range(5):\n                dfs[:,k] *= cx[k]\n            clusters=DBSCAN(eps=0.0035+ii*stepeps,min_samples=1,metric='euclidean',n_jobs=4).fit(dfs).labels_            \n            if ii==0:\n                dfh['s1'] = clusters\n                dfh['N1'] = dfh.groupby('s1')['s1'].transform('count')\n            else:\n                dfh['s2'] = clusters\n                dfh['N2'] = dfh.groupby('s2')['s2'].transform('count')\n                maxs1 = dfh['s1'].max()\n                cond = np.where((dfh['N2'].values>dfh['N1'].values) & (dfh['N2'].values<20))\n                s1 = dfh['s1'].values\n                s1[cond] = dfh['s2'].values[cond]+maxs1\n                dfh['s1'] = s1\n                dfh['s1'] = dfh['s1'].astype('int64')\n                dfh['N1'] = dfh.groupby('s1')['s1'].transform('count')\n        return dfh['s1'].values    \n    def predict(self, hits):         \n        self.clusters = self._init(hits) \n        X = self._preprocess(hits) \n        cl = hdbscan.HDBSCAN(min_samples=1,min_cluster_size=7,\n                             metric='braycurtis',cluster_selection_method='leaf',algorithm='best', leaf_size=50)\n        labels = np.unique(self.clusters)\n        self._eliminate_outliers(labels,X)\n        n_labels = 0\n        while n_labels < len(labels):\n            n_labels = len(labels)            \n            max_len = np.max(self.clusters)\n            mask = self.clusters == 0\n            self.clusters[mask] = cl.fit_predict(X[mask])+max_len\n        return self.clusters","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"33c395a8bafb76bbef06f95c47b28749f9945986","collapsed":true},"cell_type":"code","source":"model = Clusterer()\nlabels = model.predict(hits)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"335e0a72ae66f5daed4dfdf6e69e9346f0c62c79"},"cell_type":"markdown","source":"# Test Out The Model"},{"metadata":{"_uuid":"bf9161826d82f38a64fcc813d867a62a61110bae"},"cell_type":"markdown","source":"# Create Submission File"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"016fafab67ea9bc5994b0b49ed5092e6b8afb9b5"},"cell_type":"code","source":"def create_one_event_submission(event_id, hits, labels):\n    sub_data = np.column_stack(([event_id]*len(hits), hits.hit_id.values, labels))\n    submission = pd.DataFrame(data=sub_data, columns=[\"event_id\", \"hit_id\", \"track_id\"]).astype(int)\n    return submission","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b8c127672975f3fd3625a20cc34432c31bd8043f","collapsed":true},"cell_type":"code","source":"submission = create_one_event_submission(0, hits, labels)\nscore = score_event(truth, submission)\nprint(\"Your score: \", score)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8449dad5ca248eb30970f2af7d0708c977a60ef9","collapsed":true},"cell_type":"code","source":"dataset_submissions = []\ndataset_scores = []\npath_to_train = \"../input/train_1\"\n\nfor event_id, hits, cells, particles, truth in load_dataset(path_to_train, skip=0, nevents=5):\n    # Track pattern recognition\n    model = Clusterer()\n    labels = model.predict(hits)\n\n    # Prepare submission for an event\n    one_submission = create_one_event_submission(event_id, hits, labels)\n    dataset_submissions.append(one_submission)\n\n    # Score for the event\n    score = score_event(truth, one_submission)\n    dataset_scores.append(score)\n\n    print(\"Score for event %d: %.8f\" % (event_id, score))\nprint('Mean score: %.8f' % (np.mean(dataset_scores)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f75bf2aa85a7d6a083d3c5a4475f10e6b8e4b9a7","collapsed":true},"cell_type":"code","source":"path_to_test = \"../input/test\"\ntest_dataset_submissions = []\n\ncreate_submission = True # True for submission \nif create_submission:\n    for event_id, hits, cells in load_dataset(path_to_test, parts=['hits', 'cells']):\n\n        # Track pattern recognition \n        model = Clusterer()\n        labels = model.predict(hits)\n\n        # Prepare submission for an event\n        one_submission = create_one_event_submission(event_id, hits, labels)\n        test_dataset_submissions.append(one_submission)\n        \n        print('Event ID: ', event_id)\n\n    # Create submission file\n    submission = pd.concat(test_dataset_submissions, axis=0)\n    submission.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.5","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}