{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"#%matplotlib inline\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport os\nfrom trackml.dataset import load_event, load_dataset\nfrom trackml.score import score_event\n\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.neighbors import KNeighborsClassifier\n\npath_to_train = \"../input/train_1\"\nevent_prefix = \"event000001000\"\n\nhits, cells, particles, truth = load_event(os.path.join(path_to_train, event_prefix))\n\n\ndef get_training_sample(path_to_data, event_names):\n    events = []\n    track_id = 0\n\n    for name in event_names:\n        # Read an event\n        hits, cells, particles, truth = load_event(os.path.join(path_to_data, name))\n\n        # Generate new vector of particle id\n        particle_ids = truth.particle_id.values\n        particle2track = {}\n        for pid in np.unique(particle_ids):\n            particle2track[pid] = track_id\n            track_id += 1\n        hits['particle_id'] = [particle2track[pid] for pid in particle_ids]\n\n        # Collect hits\n        events.append(hits)\n\n    # Put all hits into one sample with unique tracj ids\n    data = pd.concat(events, axis=0)\n    return data\n\n\nstart_event_id = 1000\nn_train_samples = 5\ntrain_event_names = [\"event0000{:05d}\".format(i) for i in range(start_event_id, start_event_id+n_train_samples)]\ntrain_data = get_training_sample(path_to_train, train_event_names)\n\n\nclass KNNScaledClusterer(object):\n    def __init__(self):\n        self.classifier = None\n\n    def _preprocess(self, hits):\n        x = hits.x.values\n        y = hits.y.values\n        z = hits.z.values\n        r = np.sqrt(x**2 + y**2 + z**2)\n        hits['x2'] = x/r\n        hits['y2'] = y/r\n        hits['z2'] = z/r\n\n        ss = StandardScaler()\n        X = ss.fit_transform(hits[['x2', 'y2', 'z2']].values)\n        return X\n\n    def fit(self, hits):\n        X = self._preprocess(hits)\n        y = hits.particle_id.values\n        self.classifier = KNeighborsClassifier(n_neighbors=1, n_jobs=-1)\n        self.classifier.fit(X, y)\n\n    def predict(self, hits):\n        X = self._preprocess(hits)\n        labels = self.classifier.predict(X)\n        return labels\n\n\nmodel = KNNScaledClusterer()\nmodel.fit(train_data)\n\npath_to_event = os.path.join(path_to_train, \"event0000{:05d}\".format(start_event_id + n_train_samples + 1))\nhits, cells, particles, truth = load_event(path_to_event)\n\nlabels = model.predict(hits)\n\ndef 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\n\nsubmission = create_one_event_submission(0, hits, labels)\nscore = score_event(truth, submission)\n\nprint(\"Std KNN score: \", score)\n","execution_count":1,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"class KNNPolarClusterer(object):\n    def __init__(self):\n        self.classifier = None\n\n    def _preprocess(self, hits):\n        x = hits.x.values\n        y = hits.y.values\n        z = hits.z.values\n        r = np.sqrt(x**2 + y**2)\n        t = np.arctan2(y, x)\n        hits['r'] = r\n        hits['t'] = t\n\n        ss = StandardScaler()\n        X = ss.fit_transform(hits[['r', 't', 'z']].values)\n        return X\n\n    def fit(self, hits):\n        X = self._preprocess(hits)\n        y = hits.particle_id.values\n        self.classifier = KNeighborsClassifier(n_neighbors=1, n_jobs=-1)\n        self.classifier.fit(X, y)\n\n    def predict(self, hits):\n        X = self._preprocess(hits)\n        labels = self.classifier.predict(X)\n        return labels\n    \nmodel = KNNPolarClusterer()\nmodel.fit(train_data)\n\nhits, cells, particles, truth = load_event(path_to_event)\nlabels = model.predict(hits)\n\nsubmission = create_one_event_submission(0, hits, labels)\nscore = score_event(truth, submission)\nprint(\"Polar KNN score: \", score)","execution_count":3,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f961a956819069a859b23d8d90bdeaad4dc940f9"},"cell_type":"code","source":"\nfrom mpl_toolkits.mplot3d import Axes3D\n\ndef cart2spherical(cart):\n    r = np.linalg.norm(cart, axis=0)\n    theta = np.degrees(np.arccos(cart[2] / r))\n    phi = np.degrees(np.arctan2(cart[1], cart[0]))\n    return np.vstack((r, theta, phi))\n\n\n\nNUM_PARTICLES = 100\ntruth_dedup = truth.drop_duplicates('particle_id')\ntruth_sort = truth_dedup.sort_values('weight', ascending=False)\ntruth_head = truth_sort.head(NUM_PARTICLES)\n\n# Get points where the same particle intersected subsequent layers of the observation material\np_traj_list = []\nfor _, tr in truth_head.iterrows():\n    p_traj = truth[truth.particle_id == tr.particle_id][['tx', 'ty', 'tz']]\n    # Add initial position.\n    #p_traj = (p_traj\n    #          .append({'tx': particle.vx, 'ty': particle.vy, 'tz': particle.vz}, ignore_index=True)\n    #          .sort_values(by='tz'))\n    p_traj_list.append(p_traj)\n    \n# Convert to spherical coordinate.\nrtp_list = []\nfor p_traj in p_traj_list:\n    xyz = p_traj.loc[:, ['tx', 'ty', 'tz']].values.transpose()\n    rtp = cart2spherical(xyz).transpose()\n    rtp_df = pd.DataFrame(rtp, columns=('r', 'theta', 'phi'))\n    rtp_list.append(rtp_df)\n\n# Plot with Cartesian coordinates.\nfig = plt.figure(figsize=(10, 10))\nax = fig.add_subplot(111, projection='3d')\nfor p_traj in p_traj_list:\n    ax.plot(\n        xs=p_traj.tx,\n        ys=p_traj.ty,\n        zs=p_traj.tz,\n        marker='o')\nax.set_xlabel('X (mm)')\nax.set_ylabel('Y (mm)')\nax.set_zlabel('Z (mm) -- Detection layers')\nplt.title('Trajectories of top weights particles in Cartesian coordinates.')\n\n# Plot with spherical coordinates.\nfig2 = plt.figure(figsize=(10, 10))\nax = fig2.add_subplot(111, projection='3d')\nfor rtp_df in rtp_list:\n    ax.plot(\n        xs=rtp_df.theta,\n        ys=rtp_df.phi,\n        zs=rtp_df.r,\n        marker='o')\nax.set_xlabel('Theta (deg)')\nax.set_ylabel('Phi (deg)')\nax.set_zlabel('R  (mm) -- Detection layers')\nplt.title('Trajectories of top weights particles in spherical coordinates.')\nplt.show()","execution_count":4,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"75cbe75156e10b16e4885bc351d8ba23bfc5752c"},"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler\nfrom sklearn.cluster import DBSCAN\n\nclass DBSCANClusterer(object):\n    \n    def __init__(self, eps):\n        self.eps = eps\n        \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        \n        return X\n    \n    \n    def predict(self, hits):\n        \n        X = self._preprocess(hits)\n        \n        cl = DBSCAN(eps=self.eps, min_samples=1, algorithm='kd_tree')\n        labels = cl.fit_predict(X)\n        \n        return labels\n    \nmodel = DBSCANClusterer(eps=0.008)\nlabels = model.predict(hits)\nsubmission = create_one_event_submission(0, hits, labels)\nscore = score_event(truth, submission)\nprint(\"Your score: \", score)","execution_count":5,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b0689903bd79df26947e24ebafaa2b89700adb39"},"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler\nfrom sklearn.cluster import DBSCAN\n\nclass DBSCANPolarClusterer(object):\n    \n    def __init__(self, eps):\n        self.eps = eps\n        \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        t = np.arctan2(y, x)\n        p = np.arccos(y / r)\n        hits['r'] = r\n        hits['t'] = t\n        hits['p'] = p\n\n        ss = StandardScaler()\n        X = ss.fit_transform(hits[['r', 't', 'p']].values)\n        X = ss.fit_transform(hits[['t', 'p']].values)\n        \n        return X\n    \n    \n    def predict(self, hits):\n        \n        X = self._preprocess(hits)\n        \n        cl = DBSCAN(eps=self.eps, min_samples=1, algorithm='kd_tree')\n        labels = cl.fit_predict(X)\n        \n        return labels\n    \nmodel = DBSCANPolarClusterer(eps=0.008)\nlabels = model.predict(hits)\nsubmission = create_one_event_submission(0, hits, labels)\nscore = score_event(truth, submission)\nprint(\"Your score: \", score)","execution_count":6,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"24ede8d43e4967624f2d740eddd57854b0a6d876"},"cell_type":"code","source":"path_to_train = \"../input/train_1\"\nevent_prefix = \"event000001000\"\nhits, cells, particles, truth = load_event(os.path.join(path_to_train, event_prefix))\n\ndata = pd.merge(hits, truth[['hit_id', 'particle_id']], on='hit_id')\n\ndata.head()\ndata['r'] = np.sqrt(data.x*data.x + data.y*data.y + data.z*data.z)\n#data.groupby('particle_id').transform(lambda x: print(x))\nrefs = []\nfor pid, dat in data.groupby('particle_id'):\n    #print(dat[(dat.r - dat.r.median()) == (dat.r - dat.r.median()).min()])\n    ref = dat[(dat.r - dat.r.median()) == (dat.r - dat.r.median()).min()].iloc[0]\n    refs.append({'particle_id': pid, 'ref_x': ref.x, 'ref_y': ref.y, 'ref_z': ref.z})\n    \nrefs = pd.DataFrame(refs)\ndata = pd.merge(data, refs, on='particle_id')\ndata[data.x == data.ref_x].head()","execution_count":7,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d720abe2c66ff229755e2b114f212279a47a30e8"},"cell_type":"code","source":"data.shape","execution_count":15,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"941e8bf26edce0810627682876fda001aab30226"},"cell_type":"code","source":"# want to ignore particle_id == 0 because they aren't actually a single track and will fuck up the regression\ndatanonzero = data[data.particle_id != 0]\ndatanonzero.shape","execution_count":44,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9de4cdf0e1597486579ae1363f3a7d92475ec861"},"cell_type":"code","source":"# https://github.com/jcjohnson/pytorch-examples/blob/master/nn/two_layer_net_nn.py\n\nimport torch\nimport torch.optim as optim\n\ndevice = torch.device('cpu')\n# device = torch.device('cuda') # Uncomment this to run on GPU\n\n# N is batch size; D_in is input dimension;\n# H is hidden dimension; D_out is output dimension.\nN, D_in, H, D_out = 64, 3, 100, 3\n\n# Create random Tensors to hold inputs and outputs\nx = np.asarray(datanonzero[['x', 'y', 'z']])\ny = np.asarray(datanonzero[['ref_x', 'ref_y', 'ref_z']])\nfrom sklearn.preprocessing import StandardScaler\nxs = StandardScaler()\nys = StandardScaler()\nx = xs.fit_transform(x)\ny = ys.fit_transform(y)\nx = torch.autograd.Variable(torch.from_numpy(x)).float()\ny = torch.autograd.Variable(torch.from_numpy(y)).float()\n\n# Use the nn package to define our model as a sequence of layers. nn.Sequential\n# is a Module which contains other Modules, and applies them in sequence to\n# produce its output. Each Linear Module computes output from input using a\n# linear function, and holds internal Tensors for its weight and bias.\n# After constructing the model we use the .to() method to move it to the\n# desired device.\nmodel = torch.nn.Sequential(\n          torch.nn.Linear(D_in, H),\n          torch.nn.ReLU(),\n          torch.nn.Dropout(),\n          torch.nn.Linear(H, D_out),\n        ).to(device)\n#model = model.train()\n\n# The nn package also contains definitions of popular loss functions; in this\n# case we will use Mean Squared Error (MSE) as our loss function.\nloss_fn = torch.nn.MSELoss(size_average=False)\nlearning_rate = 1e-4\nlearning_rate = 1e-6\nweight_decay = 0.01\nopt = optim.SGD(model.parameters(), lr=learning_rate, weight_decay=weight_decay)\n\nfor t in range(500):\n    model.zero_grad()\n    opt.zero_grad()\n    y_pred = model(x)\n    loss = loss_fn(y_pred, y)\n    print(t, loss.item())\n    loss.backward()\n    opt.step()\n\n    # Update the weights using gradient descent. Each parameter is a Tensor, so\n    # we can access its data and gradients like we did before.\n    #with torch.no_grad():\n    #for param in model.parameters():\n    #        param.data -= learning_rate * param.grad\n    #        print(param.grad)","execution_count":64,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"45226165aaa385b30c51c5046527716a5efd6d10"},"cell_type":"code","source":"scaled = xs.transform(data[['x', 'y', 'z']])\ndata['sx'], data['sy'], data['sz'] = scaled.T\n\n# polar\ndata['polar_r'] = np.linalg.norm(data[['x', 'y', 'z']], axis=1)\ndata['polar_p'] = np.degrees(np.arccos(data['z'] / data['polar_r']))\ndata['polar_t'] = np.degrees(np.arctan2(data['y'], data['x']))\n\nss_polar = StandardScaler()\npolar_scaled = ss_polar.fit_transform(data[['polar_r', 'polar_p', 'polar_t']])\ndata['ss_polar_r'], data['ss_polar_p'], data['ss_polar_t'] = polar_scaled.T\n\n# cylindrical\ndata['cylindrical_r'] = np.linalg.norm(data[['x', 'y']], axis=1)\ndata['cylindrical_t'] = np.degrees(np.arctan2(data['y'], data['x']))\ndata['cylindrical_z'] = data['z']\n\nss_cylindrical = StandardScaler()\ncylindrical_scaled = ss_cylindrical.fit_transform(data[['cylindrical_r', 'cylindrical_t', 'cylindrical_z']])\ndata['ss_cylindrical_r'], data['ss_cylindrical_t'], data['ss_cylindrical_z'] = cylindrical_scaled.T\n\n# radial norm\ndata['radn_x'] = data['x'] / np.linalg.norm(data[['x', 'y', 'z']], axis=1)\ndata['radn_y'] = data['y'] / np.linalg.norm(data[['x', 'y', 'z']], axis=1)\ndata['radn_z'] = data['z'] / np.linalg.norm(data[['x', 'y']], axis=1)\n\nss_radn = StandardScaler()\nradn_scaled = ss_radn.fit_transform(data[['radn_x', 'radn_y', 'radn_z']])\ndata['ss_radn_x'], data['ss_radn_y'], data['ss_radn_z'] = radn_scaled.T\n\ndata['nn_med_x'], data['nn_med_y'], data['nn_med_z'] = model(torch.from_numpy(scaled).float()).data.numpy().T","execution_count":97,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9cbe64a6d8f61ffbe226a33fbe64db54905c87e2"},"cell_type":"code","source":"fig = plt.figure(figsize=(10, 10))\nax = fig.add_subplot(111, projection='3d')\nfor (pid, dat), _ in zip(data.groupby('particle_id'), range(50)):\n    if pid == 0: continue\n    ax.plot(\n        xs=dat.sx,\n        ys=dat.sy,\n        zs=dat.sz,\n        marker='o')\nax.set_xlabel('x')\nax.set_ylabel('y')\nax.set_zlabel('z')\nplt.title('scaled space')\nplt.show()","execution_count":89,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e6523c1e2aaf1da991a2bb69acbe49edd1333f5c"},"cell_type":"code","source":"fig = plt.figure(figsize=(10, 10))\nax = fig.add_subplot(111, projection='3d')\nfor (pid, dat), _ in zip(data.groupby('particle_id'), range(50)):\n    if pid == 0: continue\n    ax.plot(\n        xs=dat.ss_polar_t,\n        ys=dat.ss_polar_p,\n        zs=dat.ss_polar_r,\n        marker='o')\nax.set_xlabel('t')\nax.set_ylabel('p')\nax.set_zlabel('r')\nplt.title('scaled polar space')\nplt.show()","execution_count":90,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7b06913b17caf1889cef9b9083ae4bdddf141a55"},"cell_type":"code","source":"fig = plt.figure(figsize=(10, 10))\nax = fig.add_subplot(111, projection='3d')\nfor (pid, dat), _ in zip(data.groupby('particle_id'), range(50)):\n    if pid == 0: continue\n    ax.plot(\n        xs=dat.ss_cylindrical_t,\n        ys=dat.ss_cylindrical_z,\n        zs=dat.ss_cylindrical_r,\n        marker='o')\nax.set_xlabel('t')\nax.set_ylabel('z')\nax.set_zlabel('r')\nplt.title('scaled cylindrical space')\nplt.show()","execution_count":96,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"97d201fe2e019cc2abdf272ee6d9cd4aba503041"},"cell_type":"code","source":"fig = plt.figure(figsize=(10, 10))\nax = fig.add_subplot(111, projection='3d')\nfor (pid, dat), _ in zip(data.groupby('particle_id'), range(50)):\n    if pid == 0: continue\n    ax.plot(\n        xs=dat.ss_radn_x,\n        ys=dat.ss_radn_z,\n        zs=dat.ss_radn_y,\n        marker='o')\nax.set_xlabel('x')\nax.set_ylabel('z')\nax.set_zlabel('y')\nplt.title('scaled radial norm space')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"97d22cf50c35fb33bf1652482e14b164dc72b488"},"cell_type":"code","source":"scaled = xs.transform(data[['x', 'y', 'z']])\ndata['sx'], data['sy'], data['sz'] = scaled.T\n\nfig = plt.figure(figsize=(10, 10))\nax = fig.add_subplot(111, projection='3d')\nfor (pid, dat), _ in zip(data.groupby('particle_id'), range(50)):\n    if pid == 0: continue\n    ax.plot(\n        xs=dat.nn_med_x,\n        ys=dat.nn_med_y,\n        zs=dat.nn_med_z,\n        marker='o')\nax.set_xlabel('x')\nax.set_ylabel('y')\nax.set_zlabel('z')\nplt.title('Trajectories of top weights particles in spherical coordinates.')\nplt.show()","execution_count":91,"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}