{"cells":[{"metadata":{"_uuid":"e81259b6db616bb60121dd1f439fd1bad34b144a"},"cell_type":"markdown","source":"This notebook contains ideas borrowed from [Luis Andre Dutra e Silva](https://www.kaggle.com/mindcool/nievergelt-helix-fitting)  and [the1own](http://https://www.kaggle.com/the1owl/the-martian).\n\nThe  main ideas are:\n1.  Train a xgboost classifier to predict event labels as 0 and 1.\n2. Use HDBSCAN to cluster only events with label 1.\n3. After clustering, eliminate hits that don't lie on a quadric surface."},{"metadata":{"_cell_guid":"02a77365-f727-4acf-9025-ed67eb5763e6","_uuid":"6e637231aa068b22b8d0ce699526d6c7631a1eff"},"cell_type":"markdown","source":"Import Libraries and Events\n================"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true,"collapsed":true},"cell_type":"code","source":"from trackml.dataset import load_event\nfrom trackml.score import score_event\nfrom IPython.display import display\nimport tensorflow as tf\nimport pandas as pd\nimport numpy as np\nimport glob\n\ntrain = np.unique([p.split('-')[0] for p in sorted(glob.glob('../input/train_1/**'))])\ntest = np.unique([p.split('-')[0] for p in sorted(glob.glob('../input/test/**'))])\ndet = pd.read_csv('../input/detectors.csv')\nsub = pd.read_csv('../input/sample_submission.csv')\nprint(len(train), len(test), len(det), len(sub))","execution_count":21,"outputs":[]},{"metadata":{"_cell_guid":"a0b6cf0a-9d57-4469-96cb-a11ba4e4372c","_uuid":"5bc545e7f2449d73fadf516de336cca8716632e7"},"cell_type":"markdown","source":"Review Train Events\n========================"},{"metadata":{"_cell_guid":"80011cc9-5c91-4a76-8fa2-38615acb7e7f","_uuid":"28a4bbca57dbbee46e6b5e8ba1443b6bc1fb7b01","trusted":true,"collapsed":true},"cell_type":"code","source":"for e in train:\n    hits, cells, particles, truth = load_event(e)\n    print(len(hits), len(cells), len(particles), len(truth))\n    display(hits.head(2))\n    display(cells.head(2))\n    display(particles.head(2))\n    display(truth.head(2))\n    break","execution_count":2,"outputs":[]},{"metadata":{"_cell_guid":"626eec4c-3528-46e2-bec2-746f9f919423","_uuid":"1c18b6445cadea37521f16ed3d04b8599a6cccc7"},"cell_type":"markdown","source":"Review Test Events\n=========================="},{"metadata":{"_cell_guid":"5b79bd65-43d8-4024-bae4-a63b523c81ea","_uuid":"9f9729968f977e93cc1b145d42a39cc5811dc11c","trusted":true,"collapsed":true},"cell_type":"code","source":"for e in test:\n    hits, cells = load_event(e, parts=['hits', 'cells'])\n    print(len(hits), len(cells))\n    display(hits.head(2))\n    display(cells.head(2))\n    break","execution_count":3,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"a7813b86e86d1805ae7d350323246ec9e9651ecb"},"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler\nimport hdbscan\nfrom sklearn.neighbors import NearestNeighbors\nfrom scipy import stats\n\"\"\"\nupdated - added self.rz_scale\n\"\"\"\nclass Clusterer(object):\n    \n    def __init__(self,rz_scales=[0.65, 0.965, 1.428]):                        \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.int32)\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            x = M[index]\n            indices[i] = len(index)\n            local_mask = np.ones((len(index)),np.bool)\n            norms[i] = self._test_quadric(x)\n        threshold = np.percentile(norms,90)\n        for i, cluster in tqdm(enumerate(labels),total=len(labels)):\n            if norms[i] > threshold:\n                self.clusters[self.clusters==cluster]=0    \n                \n    def _test_quadric(self,x):\n        if len(x.shape)==3:\n            x = np.reshape(x,(x.shape[0],3))\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, rz_scales):\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(rz_scales):\n            X[:,i] = X[:,i] * rz_scale\n        \n        return X   \n    def predict(self, hits, rz_scales):        \n        volumes = np.unique(hits['volume_id'].values)\n        X = self._preprocess(hits, rz_scales)\n        self.clusters = np.zeros((len(X),1),np.int32)\n        max_len = 1\n        cl = hdbscan.HDBSCAN(min_samples=1,min_cluster_size=7,\n                             metric='braycurtis',cluster_selection_method='leaf',algorithm='best', leaf_size=50)\n        self.clusters = cl.fit_predict(X)+1\n        labels = np.unique(self.clusters)\n        n_labels = 0\n        while n_labels < len(labels):\n            n_labels = len(labels)\n            self._eliminate_outliers(labels,X)\n            max_len = np.max(self.clusters)\n            self.clusters[self.clusters==0] = cl.fit_predict(X[self.clusters==0])+max_len\n            labels = np.unique(self.clusters)\n        return self.clusters","execution_count":4,"outputs":[]},{"metadata":{"_cell_guid":"3db0d5bd-4fb7-4056-8252-693a9ab86fc5","_uuid":"c47cdf4dec2f522a077f85d7553ed3b622dc2bdc"},"cell_type":"markdown","source":"Create simple train and test datasets for further exploration\n==================================="},{"metadata":{"_cell_guid":"ead78633-64ad-4128-92e3-520916b6bc6d","_uuid":"0a26f517a46397c9af7a35031095043b3d68f52a","trusted":true,"collapsed":true},"cell_type":"code","source":"import xgboost as xgb\nfrom sklearn import *\nimport hdbscan\nfrom tqdm import tqdm\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=2, min_cluster_size=7, cluster_selection_method='leaf', prediction_data=False, metric='braycurtis')","execution_count":5,"outputs":[]},{"metadata":{"_cell_guid":"377a397f-b8bd-496b-afdb-1eef4ef5f916","_uuid":"104c8aab114bf51eafd6b7739ce2ade9765bdc7b"},"cell_type":"markdown","source":"Visualize\n==============\n* https://grechka.family/dmitry/sandbox/trackML_event_viewer/\n* https://github.com/dgrechka/TrackML_EventViewer/"},{"metadata":{"_cell_guid":"30dcb8cb-2b9d-40a1-a506-3248dfa33e20","_uuid":"46e9c79031ac483262bb024366e9036c47809061","collapsed":true,"trusted":true},"cell_type":"code","source":"#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","execution_count":6,"outputs":[]},{"metadata":{"_cell_guid":"2326db39-094e-4cf5-b2fe-3435db76cf14","_uuid":"96443569cf92a20b474f0ca6aee9fbbc68633892","trusted":true,"collapsed":true},"cell_type":"code","source":"limit = 0\ndf_train = []\ndf_val = []\nfor e in train:\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(['mean','median','max','min','sum']).reset_index()\n    cells.columns = ['hit_id'] + ['-'.join([c2,c1]) for c1 in ['mean','median','max','min','sum'] 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 ['event_id', 'hit_id', 'particle_id']+tcols]\n\n    #Classification\n    hits['target'] = hits['particle_id'].map(lambda x: 0 if x==0 else 1)\n    #tree = neighbors.KDTree(hits[['x2', 'y2', 'z2']].values)\n    #hits['dist_kd'] = hits.apply(lambda r: tree.query([[r['x2'], r['y2'], r['z2']]], k=1)[0][0][0], axis=1)\n    #hits['n_in_r_kd'] = hits.apply(lambda r: tree.query_radius([[r['x2'], r['y2'], r['z2']]], r=0.3, count_only=True)[0], axis=1)\n    #tree = neighbors.BallTree(hits[['x2', 'y2', 'z2']].values)\n    #hits['dist_ball'] = hits.apply(lambda r: tree.query([[r['x2'], r['y2'], r['z2']]], k=1)[0][0][0], axis=1)\n    #hits['n_in_r_ball'] = hits.apply(lambda r: tree.query_radius([[r['x2'], r['y2'], r['z2']]], r=0.3, count_only=True)[0], axis=1)\n    x1, x2, y1, y2 = model_selection.train_test_split(hits[cols], hits['target'].values, test_size=0.20, random_state=19)\n\n    train_ = x1.copy()\n    train_['target'] = y1\n    df_train.append(train_.copy())\n    \n    val = x2.copy()\n    val['target'] = y2\n    df_val.append(val.copy())\n    \n    params = {'eta': 0.2, 'max_depth': 7, 'objective':'binary:logistic', 'seed': 19, 'silent': True}\n    \n    if limit == 0:\n        watchlist = [(xgb.DMatrix(x1, y1), 'train'), (xgb.DMatrix(x2, y2), 'valid')]\n        model = xgb.train(params, xgb.DMatrix(x1, y1), 500,  watchlist, verbose_eval=10, early_stopping_rounds=20)\n    else:\n        train_ = pd.concat(df_train, ignore_index=True).reset_index(drop=True)\n        val = pd.concat(df_val, ignore_index=True).reset_index(drop=True)\n        watchlist = [(xgb.DMatrix(train_[cols], train_['target'].values), 'train'), (xgb.DMatrix(val[cols], val['target'].values), 'valid')]\n        #model.load_model('tml_xgb_model.model')\n        #model = xgb.train(params, xgb.DMatrix(train_[cols], train_['target'].values), 500,  watchlist, verbose_eval=10, early_stopping_rounds=20, xgb_model=model)\n        model = xgb.train(params, xgb.DMatrix(train_[cols], train_['target'].values), 500,  watchlist, verbose_eval=10, early_stopping_rounds=20)\n    #model.save_model('tml_xgb_model.model')\n    \n    preds = model.predict(xgb.DMatrix(hits[cols]), ntree_limit=model.best_ntree_limit)\n    hits['track_id'] = [1 if p>0.2 else 0 for p in preds]\n\n    #Add TensorFlow reinforcement learning here for N cluster optimizer...\n    \n    if len(hits[hits['track_id']==0])>0:\n        hits0 = hits[hits['track_id']==0].reset_index(drop=True)\n        hits1 = hits[hits['track_id']==1].reset_index(drop=True)\n        print(len(hits0), len(hits1))\n        hits1['track_id'] = dbscan.fit_predict(scl.fit_transform(hits1[['x2', 'y2', 'z2']].values)) + 1\n        hits = pd.concat([hits0, hits1], ignore_index=True).fillna(0).reset_index(drop=True)\n    else:\n        hits['track_id'] = dbscan.fit_predict(scl.fit_transform(hits[['x2', 'y2', 'z2']].values)) + 1\n\n    score = score_event(hits[tcols], hits[['event_id','hit_id','track_id']])\n    print(e, len(hits), len(truth['particle_id'].unique()), len(hits['track_id'].unique()), score)\n    if limit > 5:\n        break\n    limit += 1","execution_count":7,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"c5a1174e1730a60d8c18a0fb9fda3245eef70e37"},"cell_type":"code","source":"RZ_SCALE = [0.65, 0.965, 1.428]\nLEAF_SIZE = 50","execution_count":18,"outputs":[]},{"metadata":{"_cell_guid":"0b548826-4fe4-466a-92f7-2101753733b0","_uuid":"87f3d68cbab22fda4c99cce9e9071113e4b0b980","trusted":true,"collapsed":true},"cell_type":"code","source":"df_test = []\nfor e in test:\n    hits, cells = load_event(e, parts=['hits', 'cells'])\n    hits['event_id'] = int(e[-9:])\n    cells = cells.groupby(by=['hit_id'])['ch0', 'ch1', 'value'].agg(['mean','median','max','min','sum']).reset_index()\n    cells.columns = ['hit_id'] + ['-'.join([c2,c1]) for c1 in ['mean','median','max','min','sum'] for c2 in ['ch0', 'ch1', 'value']]\n    hits = pd.merge(hits, cells, how='left', on='hit_id')\n    col = [c for c in hits.columns if c not in ['event_id', 'hit_id', 'particle_id']]\n    hits = norm_points(hits)\n\n    #Classification\n    #tree = neighbors.KDTree(hits[['x2', 'y2', 'z2']].values)\n    #hits['dist_kd'] = hits.apply(lambda r: tree.query([[r['x2'], r['y2'], r['z2']]], k=1)[0][0][0], axis=1)\n    #hits['n_in_r_kd'] = hits.apply(lambda r: tree.query_radius([[r['x2'], r['y2'], r['z2']]], r=0.3, count_only=True)[0], axis=1)\n    #tree = neighbors.BallTree(hits[['x2', 'y2', 'z2']].values)\n    #hits['dist_ball'] = hits.apply(lambda r: tree.query([[r['x2'], r['y2'], r['z2']]], k=1)[0][0][0], axis=1)\n    #hits['n_in_r_ball'] = hits.apply(lambda r: tree.query_radius([[r['x2'], r['y2'], r['z2']]], r=0.3, count_only=True)[0], axis=1)\n    preds = model.predict(xgb.DMatrix(hits[cols]), ntree_limit=model.best_ntree_limit)\n    hits['particle_id'] = [1 if p>0.2 else 0 for p in preds]\n\n    if len(hits[hits['particle_id']==0])>0:\n        hits0 = hits[hits['particle_id']==0].reset_index(drop=True)\n        hits1 = hits[hits['particle_id']==1].reset_index(drop=True)\n        # Track pattern recognition\n        try:\n            model_cluster = Clusterer()\n            labels = model_cluster.predict(hits1, RZ_SCALE)\n            hits1['particle_id'] = labels + 1\n        except:\n            pass\n#         hits1['particle_id'] = dbscan.fit_predict(scl.fit_transform(hits1[['x2', 'y2', 'z2']].values)) + 1\n        hits = pd.concat([hits0, hits1], ignore_index=True).fillna(0).reset_index(drop=True)\n    else:\n        try:\n            model_cluster = Clusterer()\n            labels = model_cluster.predict(hits, RZ_SCALE)\n            hits['particle_id'] = labels + 1\n        except:\n            pass\n#         hits1['particle_id'] = labels + 1  #seems fishy !!!\n          \n#         hits['particle_id'] = dbscan.fit_predict(scl.fit_transform(hits[['x2', 'y2', 'z2']].values)) + 1\n    df_test.append(hits[['event_id','hit_id','particle_id']].copy())\n    print(e, len(hits['particle_id'].unique()))\n\ndf_test = pd.concat(df_test, ignore_index=True)\ndf_test.head()","execution_count":20,"outputs":[]},{"metadata":{"_cell_guid":"9532ec4b-a752-41f2-8d78-e0664888316c","_uuid":"15b6afef0647da30162dcf05bbe57386f6b0ec19","collapsed":true,"trusted":true},"cell_type":"code","source":"sub = pd.merge(sub, df_test, how='left', on=['event_id','hit_id'])\nsub['track_id'] = sub['particle_id'].astype(int)\nsub[['event_id','hit_id','track_id']].to_csv('submission-001.csv', index=False)\n#sub.to_csv('submission-001.csv.gzip',index=False, compression='gzip')\n#!kaggle competitions submit -c trackml-particle-identification -f submission-001.csv.gzip -m \"Happy Kaggling :)\"","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}