{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"collapsed":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\nfrom keras.models import Model, Sequential, load_model\nfrom keras.layers import Dense, Input\nfrom keras.optimizers import Adam\nfrom tqdm import tqdm_notebook\nprint(os.listdir(\"../input\"))\nprint(os.listdir(\"../input/trackml/\"))\nprefix='../input/trackml-particle-identification/'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"6cfcc4d7929c3744e906026c56c9f8d03a0a46df"},"cell_type":"code","source":"def init_model(fs = 10):\n    model = Sequential()\n    model.add(Dense(800, activation='selu', input_shape=(fs,)))\n    model.add(Dense(400, activation='selu'))\n    model.add(Dense(400, activation='selu'))\n    model.add(Dense(400, activation='selu'))\n    model.add(Dense(200, activation='selu'))\n    model.add(Dense(1, activation='sigmoid'))\n    return model\n\ndef get_event(event):\n    hits= pd.read_csv(prefix+'train_1/%s-hits.csv'%event)\n    cells= pd.read_csv(prefix+'train_1/%s-cells.csv'%event)\n    truth= pd.read_csv(prefix+'train_1/%s-truth.csv'%event)\n    particles= pd.read_csv(prefix+'train_1/%s-particles.csv'%event)\n    return hits, cells, truth, particles","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"22ca78487ca459c532750dd133c7ed401e592e3d"},"cell_type":"markdown","source":"# Step 1 - Prepare training data\n* use 10 events for training\n* input: hit pair\n* output: 1 if two hits are the same particle_id, 0 otherwise.\n* feature size: 10 (5 per hit)"},{"metadata":{"trusted":true,"_uuid":"8a8aa70428d754645a5ab2c74f7bc6c3c99505b7","collapsed":true},"cell_type":"code","source":"# you can jump to step4 for test only.\n\nTrain = []\nfor i in tqdm_notebook(range(10,20)):\n    event = 'event0000010%02d'%i\n    hits, cells, truth, particles = get_event(event)\n    hit_cells = cells.groupby(['hit_id']).value.count().values\n    hit_value = cells.groupby(['hit_id']).value.sum().values\n    features = np.hstack((hits[['x','y','z']]/1000, hit_cells.reshape(len(hit_cells),1)/10,hit_value.reshape(len(hit_cells),1)))\n    particle_ids = truth.particle_id.unique()\n    particle_ids = particle_ids[np.where(particle_ids!=0)[0]]\n    \n    pair = []\n    for particle_id in particle_ids:\n        hit_ids = truth[truth.particle_id == particle_id].hit_id.values-1\n        for i in hit_ids:\n            for j in hit_ids:\n                if i != j:\n                    pair.append([i,j])\n    pair = np.array(pair)   \n    Train1 = np.hstack((features[pair[:,0]], features[pair[:,1]], np.ones((len(pair),1))))\n    \n    if len(Train) == 0:\n        Train = Train1\n    else:\n        Train = np.vstack((Train,Train1))\n    \n    n = len(hits)\n    size = len(Train1)*3\n    p_id = truth.particle_id.values\n    i =np.random.randint(n, size=size)\n    j =np.random.randint(n, size=size)\n    pair = np.hstack((i.reshape(size,1),j.reshape(size,1)))\n    pair = pair[((p_id[i]==0) | (p_id[i]!=p_id[j]))]\n\n    Train0 = np.hstack((features[pair[:,0]], features[pair[:,1]], np.zeros((len(pair),1))))\n\n    print(event, Train1.shape)\n    \n    Train = np.vstack((Train,Train0))\ndel Train0, Train1\n\nnp.random.shuffle(Train)\nprint(Train.shape)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9f38fbd32829f9ca2a18f8b00bf3db26f57f4297"},"cell_type":"markdown","source":"# Step 2 - Train model"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"3bab2db55c8d0abf38119e4f6780269adaaf906e"},"cell_type":"code","source":"model = init_model()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9e96d105f22f0fb6d605332c43994b7608744bcc","collapsed":true},"cell_type":"code","source":"lr=-5\nmodel.compile(loss=['binary_crossentropy'], optimizer=Adam(lr=10**(lr)), metrics=['accuracy'])\nHistory = model.fit(x=Train[:,:-1], y=Train[:,-1], batch_size=8000, epochs=1, verbose=2, validation_split=0.05, shuffle=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"29eb9b3b2db4909cf9993435075bf0561b2cc468","collapsed":true},"cell_type":"code","source":"lr=-4\nmodel.compile(loss=['binary_crossentropy'], optimizer=Adam(lr=10**(lr)), metrics=['accuracy'])\nHistory = model.fit(x=Train[:,:-1], y=Train[:,-1], batch_size=8000, epochs=20, verbose=2, validation_split=0.05, shuffle=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ab41ee83a17a170f907b0789379420cee8da037f","collapsed":true},"cell_type":"code","source":"lr=-5\nmodel.compile(loss=['binary_crossentropy'], optimizer=Adam(lr=10**(lr)), metrics=['accuracy'])\nHistory = model.fit(x=Train[:,:-1], y=Train[:,-1], batch_size=8000, epochs=3, verbose=2, validation_split=0.05, shuffle=True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"185872dc6a27b54f7b810c66a5076796c8c8ec01"},"cell_type":"markdown","source":"# Step 3 - Hard Negative Mining"},{"metadata":{"trusted":true,"_uuid":"20a591651544390c1940aedd8c48a4f1526e1f8f","collapsed":true},"cell_type":"code","source":"# if you skip step2, you still need to run step1 to get training data.\ntry:\n    model\nexcept NameError:\n    print('load model')\n    model = load_model('../input/trackml/my_model.h5')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"64ce296b3f18f0b4dd368ee232bb9aaa7655365f","collapsed":true},"cell_type":"code","source":"Train_hard = []\n\nfor i in tqdm_notebook(range(10,20)):\n\n    event = 'event0000010%02d'%i\n    hits, cells, truth, particles = get_event(event)\n    hit_cells = cells.groupby(['hit_id']).value.count().values\n    hit_value = cells.groupby(['hit_id']).value.sum().values\n    features = np.hstack((hits[['x','y','z']]/1000, hit_cells.reshape(len(hit_cells),1)/10,hit_value.reshape(len(hit_cells),1)))\n\n    size=30000000\n    n = len(truth)\n    i =np.random.randint(n, size=size)\n    j =np.random.randint(n, size=size)\n    p_id = truth.particle_id.values\n    pair = np.hstack((i.reshape(size,1),j.reshape(size,1)))\n    pair = pair[((p_id[i]==0) | (p_id[i]!=p_id[j]))]\n\n    Train0 = np.hstack((features[pair[:,0]], features[pair[:,1]], np.zeros((len(pair),1))))\n    \n    pred = model.predict(Train0[:,:-1], batch_size=20000)\n    s = np.where(pred>0.5)[0]\n    \n    print(event, len(Train0), len(s))\n    \n    if len(Train_hard) == 0:\n        Train_hard = Train0[s]\n    else:\n        Train_hard = np.vstack((Train_hard,Train0[s]))\ndel Train0\nprint(Train_hard.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"72e5d14bd8617a4e3c5175cea7705cdabcaa0f5e"},"cell_type":"code","source":"Train = np.vstack((Train,Train_hard))\nnp.random.shuffle(Train)\nprint(Train.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"ba695b07a4d35af53127564ba11b7884d56fd1a9"},"cell_type":"code","source":"lr=-4\nmodel.compile(loss=['binary_crossentropy'], optimizer=Adam(lr=10**(lr)), metrics=['accuracy'])\nHistory = model.fit(x=Train[:,:-1], y=Train[:,-1], batch_size=8000, epochs=30, verbose=2, validation_split=0.05, shuffle=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"ffecdec926f9645fe4ddac507787980fbf6142a6"},"cell_type":"code","source":"lr=-5\nmodel.compile(loss=['binary_crossentropy'], optimizer=Adam(lr=10**(lr)), metrics=['accuracy'])\nHistory = model.fit(x=Train[:,:-1], y=Train[:,-1], batch_size=8000, epochs=10, verbose=2, validation_split=0.05, shuffle=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"b81f5b8a3e882fae96bc923f79720ef3f130ca5d"},"cell_type":"code","source":"lr=-6\nmodel.compile(loss=['binary_crossentropy'], optimizer=Adam(lr=10**(lr)), metrics=['accuracy'])\nHistory = model.fit(x=Train[:,:-1], y=Train[:,-1], batch_size=8000, epochs=2, verbose=2, validation_split=0.05, shuffle=True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1e5db9b7e283f01e8016e919c0e4d4c4731bf247"},"cell_type":"markdown","source":"# Step 4 - Test event 1001"},{"metadata":{"trusted":true,"_uuid":"ba7f35bbd13b606908207e4dbea86190a681f8cc","collapsed":true},"cell_type":"code","source":"try:\n    model\nexcept NameError:\n    print('load model')\n    model = load_model('../input/trackml/my_model_h.h5')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"b87b9218f8426fc7b86d09156d5d461713993e8b"},"cell_type":"code","source":"event = 'event000001001'\nhits, cells, truth, particles = get_event(event)\nhit_cells = cells.groupby(['hit_id']).value.count().values\nhit_value = cells.groupby(['hit_id']).value.sum().values\nfeatures = np.hstack((hits[['x','y','z']]/1000, hit_cells.reshape(len(hit_cells),1)/10,hit_value.reshape(len(hit_cells),1)))\ncount = hits.groupby(['volume_id','layer_id','module_id'])['hit_id'].count().values\nmodule_id = np.zeros(len(hits), dtype='int32')\n\nfor i in range(len(count)):\n    si = np.sum(count[:i])\n    module_id[si:si+count[i]] = i","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"f882457cf0e251f4e8e8adce0629d381063a3ca2"},"cell_type":"code","source":"def get_path(hit, mask, thr):\n    path = [hit]\n    a = 0\n    while True:\n        c = get_predict(path[-1], thr/2)\n        mask = (c > thr)*mask\n        mask[path[-1]] = 0\n        \n        if 1:\n            cand = np.where(c>thr)[0]\n            if len(cand)>0:\n                mask[cand[np.isin(module_id[cand], module_id[path])]]=0\n                \n        a = (c + a)*mask\n        if a.max() < thr*len(path):\n            break\n        path.append(a.argmax())\n    return path\n\ndef get_predict(hit, thr=0.5):\n    Tx = np.zeros((len(truth),10))\n    Tx[:,5:] = features\n    Tx[:,:5] = np.tile(features[hit], (len(Tx), 1))\n    pred = model.predict(Tx, batch_size=len(Tx))[:,0]\n    # TTA\n    idx = np.where(pred > thr)[0]\n    Tx2 = np.zeros((len(idx),10))\n    Tx2[:,5:] = Tx[idx,:5]\n    Tx2[:,:5] = Tx[idx,5:]    \n    pred1 = model.predict(Tx2, batch_size=len(idx))[:,0]\n    pred[idx] = (pred[idx] + pred1)/2\n    return pred","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"98400d020f1339b0a467047d3acff43925f39080","collapsed":true},"cell_type":"code","source":"# select one hit to construct a track\nfor hit in range(3):\n    path = get_path(hit, np.ones(len(truth)), 0.95)\n    gt = np.where(truth.particle_id==truth.particle_id[hit])[0]\n    print('hit_id = ', hit+1)\n    print('reconstruct :', path)\n    print('ground truth:', gt.tolist())","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b218782577435528ac3e131eac33e7d4b8f64d75"},"cell_type":"markdown","source":"# Step 5 - Predict and Score\n"},{"metadata":{"trusted":true,"_uuid":"6dd697c803477e88b7d0bd98b24e5c3c649898ae","collapsed":true},"cell_type":"code","source":"# Predict all pairs for reconstruct by all hits. (takes 2.5hr but can skip)\nskip_predict = True\n\nif skip_predict == False:\n    TestX = np.zeros((len(features), 10))\n    TestX[:,5:] = features\n\n    # for TTA\n    TestX1 = np.zeros((len(features), 10))\n    TestX1[:,:5] = features\n\n    preds = []\n\n    for i in tqdm_notebook(range(len(features)-1)):\n        TestX[i+1:,:5] = np.tile(features[i], (len(TestX)-i-1, 1))\n\n        pred = model.predict(TestX[i+1:], batch_size=20000)[:,0]                \n        idx = np.where(pred>0.2)[0]\n\n        if len(idx) > 0:\n            TestX1[idx+i+1,5:] = TestX[idx+i+1,:5]\n            pred1 = model.predict(TestX1[idx+i+1], batch_size=20000)[:,0]\n            pred[idx] = (pred[idx]+pred1)/2\n\n        idx = np.where(pred>0.5)[0]\n\n        preds.append([idx+i+1, pred[idx]])\n\n        #if i==0: print(preds[-1])\n\n    preds.append([np.array([], dtype='int64'), np.array([], dtype='float32')])\n\n    # rebuild to NxN\n    for i in range(len(preds)):\n        ii = len(preds)-i-1\n        for j in range(len(preds[ii][0])):\n            jj = preds[ii][0][j]\n            preds[jj][0] = np.insert(preds[jj][0], 0 ,ii)\n            preds[jj][1] = np.insert(preds[jj][1], 0 ,preds[ii][1][j])\n\n    #np.save('my_%s.npy'%event, preds)\nelse:\n    print('load predicts')\n    preds = np.load('../input/trackml/my_%s.npy'%event)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"0d34e0bcdc6ba971b4bc7bab1e8c185e17523956"},"cell_type":"code","source":"def get_path2(hit, mask, thr):\n    path = [hit]\n    a = 0\n    while True:\n        c = get_predict2(path[-1])\n        mask = (c > thr)*mask\n        mask[path[-1]] = 0\n        \n        if 1:\n            cand = np.where(c>thr)[0]\n            if len(cand)>0:\n                mask[cand[np.isin(module_id[cand], module_id[path])]]=0\n                \n        a = (c + a)*mask\n        if a.max() < thr*len(path):\n            break\n        path.append(a.argmax())\n    return path\n\ndef get_predict2(p):\n    c = np.zeros(len(preds))\n    c[preds[p, 0]] = preds[p, 1]          \n    return c","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0a1af92b508896ec7fc6da0d5f89f153246f73b1","collapsed":true},"cell_type":"code","source":"# reconstruct by all hits. (takes 0.6hr but can skip)\nskip_reconstruct = True\n\nif skip_reconstruct == False:\n    tracks_all = []\n    thr = 0.85\n    x4 = True\n    for hit in tqdm_notebook(range(len(preds))):\n        m = np.ones(len(truth))\n        path  = get_path2(hit, m, thr)\n        if x4 and len(path) > 1:\n            m[path[1]]=0\n            path2  = get_path2(hit, m, thr)\n            if len(path) < len(path2):\n                path = path2\n                m[path[1]]=0\n                path2  = get_path2(hit, m, thr)\n                if len(path) < len(path2):\n                    path = path2\n            elif len(path2) > 1:\n                m[path[1]]=1\n                m[path2[1]]=0\n                path2  = get_path2(hit, m, thr)\n                if len(path) < len(path2):\n                    path = path2\n        tracks_all.append(path)\n    #np.save('my_tracks_all', tracks_all)\nelse:\n    print('load tracks')\n    tracks_all = np.load('../input/trackml/my_tracks_all.npy')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"ce9865fb97fd16066e062bc375027cb405cf1f51"},"cell_type":"code","source":"def get_track_score(tracks_all, n=4):\n    scores = np.zeros(len(tracks_all))\n    for i, path in enumerate(tracks_all):\n        count = len(path)\n\n        if count > 1:\n            tp=0\n            fp=0\n            for p in path:\n                tp = tp + np.sum(np.isin(tracks_all[p], path, assume_unique=True))\n                fp = fp + np.sum(np.isin(tracks_all[p], path, assume_unique=True, invert=True))\n            scores[i] = (tp-fp*n-count)/count/(count-1)\n        else:\n            scores[i] = -np.inf\n    return scores\n\ndef score_event_fast(truth, submission):\n    truth = truth[['hit_id', 'particle_id', 'weight']].merge(submission, how='left', on='hit_id')\n    df = truth.groupby(['track_id', 'particle_id']).hit_id.count().to_frame('count_both').reset_index()\n    truth = truth.merge(df, how='left', on=['track_id', 'particle_id'])\n    \n    df1 = df.groupby(['particle_id']).count_both.sum().to_frame('count_particle').reset_index()\n    truth = truth.merge(df1, how='left', on='particle_id')\n    df1 = df.groupby(['track_id']).count_both.sum().to_frame('count_track').reset_index()\n    truth = truth.merge(df1, how='left', on='track_id')\n    truth.count_both *= 2\n    score = truth[(truth.count_both > truth.count_particle) & (truth.count_both > truth.count_track)].weight.sum()\n    particles = truth[(truth.count_both > truth.count_particle) & (truth.count_both > truth.count_track)].particle_id.unique()\n\n    return score, truth[truth.particle_id.isin(particles)].weight.sum(), 1-truth[truth.track_id>0].weight.sum()\n\ndef evaluate_tracks(tracks, truth):\n    submission = pd.DataFrame({'hit_id': truth.hit_id, 'track_id': tracks})\n    score = score_event_fast(truth, submission)[0]\n    track_id = tracks.max()\n    print('%.4f %2.2f %4d %5d %.4f %.4f'%(score, np.sum(tracks>0)/track_id, track_id, np.sum(tracks==0), 1-score-np.sum(truth.weight.values[tracks==0]), np.sum(truth.weight.values[tracks==0])))\n\ndef extend_path(path, mask, thr, last = False):\n    a = 0\n    for p in path[:-1]:\n        c = get_predict2(p)\n        if last == False:\n            mask = (c > thr)*mask\n        mask[p] = 0\n        cand = np.where(c>thr)[0]\n        mask[cand[np.isin(module_id[cand], module_id[path])]]=0\n        a = (c + a)*mask\n\n    while True:\n        c = get_predict2(path[-1])\n        if last == False:\n            mask = (c > thr)*mask\n        mask[path[-1]] = 0\n        cand = np.where(c>thr)[0]\n        mask[cand[np.isin(module_id[cand], module_id[path])]]=0\n        a = (c + a)*mask\n            \n        if a.max() < thr*len(path):\n            break\n\n        path.append(a.argmax())\n        if last: break\n    \n    return path","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"d4670b7668b14f70b2c64727ae48ff319e2c4552"},"cell_type":"code","source":"# calculate track's confidence (about 2 mins)\nscores = get_track_score(tracks_all, 8)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6d5bde61fd95aba5e09536d81360003742857062","collapsed":true},"cell_type":"code","source":"# merge tracks by confidence and get score\nidx = np.argsort(scores)[::-1]\ntracks = np.zeros(len(hits))\ntrack_id = 0\n\nfor hit in idx:\n\n    path = np.array(tracks_all[hit])\n    path = path[np.where(tracks[path]==0)[0]]\n\n    if len(path)>3:\n        track_id = track_id + 1  \n        tracks[path] = track_id\n\nevaluate_tracks(tracks, truth)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"caa0cbbee67c1b4a6eb749d495f865436f3b25eb","collapsed":true},"cell_type":"code","source":"# multistage\nidx = np.argsort(scores)[::-1]\ntracks = np.zeros(len(hits))\ntrack_id = 0\n\nfor hit in idx:\n    path = np.array(tracks_all[hit])\n    path = path[np.where(tracks[path]==0)[0]]\n\n    if len(path)>6:\n        track_id = track_id + 1  \n        tracks[path] = track_id\n\nevaluate_tracks(tracks, truth)\n\nfor track_id in range(1, int(tracks.max())+1):\n    path = np.where(tracks == track_id)[0]\n    path = extend_path(path.tolist(), 1*(tracks==0), 0.6)\n    tracks[path] = track_id\n        \nevaluate_tracks(tracks, truth)\n        \nfor hit in idx:\n    path = np.array(tracks_all[hit])\n    path = path[np.where(tracks[path]==0)[0]]\n\n    if len(path)>3:\n        path = extend_path(path.tolist(), 1*(tracks==0), 0.6)\n        track_id = track_id + 1  \n        tracks[path] = track_id\n        \nevaluate_tracks(tracks, truth)\n\nfor track_id in range(1, int(tracks.max())+1):\n    path = np.where(tracks == track_id)[0]\n    path = extend_path(path.tolist(), 1*(tracks==0), 0.5)\n    tracks[path] = track_id\n        \nevaluate_tracks(tracks, truth)\n\nfor hit in idx:\n    path = np.array(tracks_all[hit])\n    path = path[np.where(tracks[path]==0)[0]]\n\n    if len(path)>1:\n        path = extend_path(path.tolist(), 1*(tracks==0), 0.5)\n    if len(path)>2:\n        track_id = track_id + 1\n        tracks[path] = track_id\n        \nevaluate_tracks(tracks, truth)\n\nfor track_id in range(1, int(tracks.max())+1):\n    path = np.where(tracks== track_id)[0]\n    if len(path)%2 == 0:\n        path = extend_path(path.tolist(), 1*(tracks==0), 0.5, True)\n        tracks[path] = track_id\n        \nevaluate_tracks(tracks, truth)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"64f8a293a76cac80636c80a4adb8fe994fc8cef4"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}