XGBoost.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xgboost as xgb\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from tqdm import tqdm_notebook\n",
    "prefix='../'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "def 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\n",
    "\n",
    "def get_event_sample(event):\n",
    "    hits= pd.read_csv(prefix+'train_sample/%s-hits.csv'%event)\n",
    "    cells= pd.read_csv(prefix+'train_sample/%s-cells.csv'%event)\n",
    "    truth= pd.read_csv(prefix+'train_sample/%s-truth.csv'%event)\n",
    "    particles= pd.read_csv(prefix+'train_sample/%s-particles.csv'%event)\n",
    "    return hits, cells, truth, particles\n",
    "\n",
    "def convert_to_cylindrical_coordinates(hits):\n",
    "    x, y = hits['x'], hits['y']\n",
    "    hits['r'] = np.sqrt(x**2 + y**2)\n",
    "    hits['a0'] = np.arctan2(y, x)\n",
    "    return hits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dtrain = xgb.DMatrix('train.csv?format=csv&label_column=0')\n",
    "dtest = xgb.DMatrix('test.csv?format=csv&label_column=0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "20f6fd1b0d084ea1bf8316375413ce76",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "event000021100 (911412, 17)\n",
      "\n",
      "(3645403, 17)\n"
     ]
    }
   ],
   "source": [
    "train = True\n",
    "if train:\n",
    "    Train = []\n",
    "    for i in tqdm_notebook(range(11,12)):\n",
    "        #event = 'event0000010%02d'%i\n",
    "        event = \"event000021100\"\n",
    "        hits, cells, truth, particles = get_event_sample(event)\n",
    "        hit_cells = cells.groupby(['hit_id']).value.count().values\n",
    "        hit_value = cells.groupby(['hit_id']).value.sum().values\n",
    "        cartesian_features = np.array(hits[['x','y','z']]/1000).reshape(len(hit_cells), 3)\n",
    "        #features = np.array(hits[['x','y','z']]/1000), hit_cells.reshape(len(hit_cells),1)/10,hit_value.reshape(len(hit_cells),1)))\n",
    "        hits = convert_to_cylindrical_coordinates(hits)\n",
    "        features = np.hstack((np.hstack((np.array(hits['z']/1000).reshape(np.array(hits['x']).shape[0], 1), \n",
    "                                         np.array(hits['r']/1000).reshape(np.array(hits['r']).shape[0], 1))), \n",
    "                              np.array(hits['a0']).reshape(len(hits['a0'].iloc[:]), 1)))\n",
    "        hit_signals = np.hstack((hit_cells.reshape(len(hit_cells),1)/10,hit_value.reshape(len(hit_cells),1)))\n",
    "        features = np.hstack((features, hit_signals))\n",
    "        features = np.hstack((features, cartesian_features))\n",
    "        #features = np.hstack((features, np.array(hits[['x','y','z']]/1000).reshape(len(hit_cells), 3)))\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))\n",
    "    del Train0, Train1\n",
    "\n",
    "    np.random.shuffle(Train)\n",
    "    print(Train.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_Train = np.load('train10events.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[07:12:51] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 116 extra nodes, 0 pruned nodes, max_depth=6\n",
      "[07:12:51] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 116 extra nodes, 0 pruned nodes, max_depth=6\n",
      "[07:12:51] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 118 extra nodes, 0 pruned nodes, max_depth=6\n",
      "[07:12:51] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 110 extra nodes, 0 pruned nodes, max_depth=6\n",
      "[07:12:51] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 110 extra nodes, 0 pruned nodes, max_depth=6\n",
      "[07:12:51] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 108 extra nodes, 0 pruned nodes, max_depth=6\n",
      "[07:12:51] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 122 extra nodes, 0 pruned nodes, max_depth=6\n",
      "[07:12:51] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 120 extra nodes, 0 pruned nodes, max_depth=6\n",
      "[07:12:51] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 90 extra nodes, 0 pruned nodes, max_depth=6\n",
      "[07:12:51] src/tree/updater_prune.cc:74: tree pruning end, 1 roots, 124 extra nodes, 0 pruned nodes, max_depth=6\n"
     ]
    }
   ],
   "source": [
    "dtrain = xgb.DMatrix(data = s_Train[:20000,:-1], label = s_Train[:20000,-1:])\n",
    "param = {'booster':'gbtree'}\n",
    "num_round = 10\n",
    "bst = xgb.train(param, dtrain)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'[0]\\teval-rmse:0.353906'"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dtest = xgb.DMatrix(data = s_Train[-200000:,:-1], label = s_Train[-200000:,-1:])\n",
    "bst.eval(dtest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xgboost import XGBClassifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import accuracy_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.6/site-packages/sklearn/preprocessing/label.py:219: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n",
      "/anaconda3/lib/python3.6/site-packages/sklearn/preprocessing/label.py:252: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3h 33min, sys: 1min 3s, total: 3h 34min 3s\n",
      "Wall time: 3h 48min 53s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "model = XGBClassifier(max_depth=25, n_estimators=500)\n",
    "model.fit(Train[:2000000,:-1], Train[:2000000,-1:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy: 92.33%\n",
      "CPU times: user 22.8 s, sys: 113 ms, total: 22.9 s\n",
      "Wall time: 23.1 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "model.nthread = 4\n",
    "y_pred = model.predict(Train[-50000:, :-1])\n",
    "#predictions = [round(value) for value in y_pred]\n",
    "# evaluate predictions\n",
    "accuracy = accuracy_score(Train[-50000:,-1:], y_pred)\n",
    "print(\"Accuracy: %.2f%%\" % (accuracy * 100.0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event = 'event000001001'\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",
    "cartesian_features = np.hstack((hits[['x','y','z']]/1000, hit_cells.reshape(len(hit_cells),1)/10,hit_value.reshape(len(hit_cells),1)))\n",
    "hits = convert_to_cylindrical_coordinates(hits)\n",
    "features = np.hstack((np.hstack((np.array(hits['z']/1000).reshape(np.array(hits['x']).shape[0], 1), \n",
    "                                 np.array(hits['r']/1000).reshape(np.array(hits['r']).shape[0], 1))), \n",
    "                      np.array(hits['a0']).reshape(len(hits['a0'].iloc[:]), 1)))\n",
    "hit_signals = np.hstack((hit_cells.reshape(len(hit_cells),1)/10,hit_value.reshape(len(hit_cells),1)))\n",
    "features = np.hstack((features, hit_signals))\n",
    "features = np.hstack((features, cartesian_features))\n",
    "count = hits.groupby(['volume_id','layer_id','module_id'])['hit_id'].count().values\n",
    "module_id = np.zeros(len(hits), dtype='int32')\n",
    "\n",
    "for i in range(len(count)):\n",
    "    si = np.sum(count[:i])\n",
    "    module_id[si:si+count[i]] = i\n",
    "\n",
    "skip_predict = False\n",
    "\n",
    "if 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)\n",
    "else:\n",
    "    print('load predicts')\n",
    "    preds = np.load('../input/trackml/my_%s.npy'%event)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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",
    "\n",
    "def get_predict2(p):\n",
    "    c = np.zeros(len(preds))\n",
    "    c[preds[p, 0]] = preds[p, 1]          \n",
    "    return c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "skip_reconstruct = False\n",
    "\n",
    "if 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)\n",
    "else:\n",
    "    print('load tracks')\n",
    "    tracks_all = np.load('../input/trackml/my_tracks_all.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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",
    "\n",
    "def 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",
    "\n",
    "def 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",
    "\n",
    "def 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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate track's confidence (about 2 mins)\n",
    "scores = get_track_score(tracks_all, 8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge tracks by confidence and get score\n",
    "idx = np.argsort(scores)[::-1]\n",
    "tracks = np.zeros(len(hits))\n",
    "track_id = 0\n",
    "\n",
    "for 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",
    "\n",
    "evaluate_tracks(tracks, truth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.save_model('XGBoost_500_estimators_25_depth.model')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "hits, cells, particles, truth = get_event_sample('event000021100')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}