{"cells":[{"metadata":{"_uuid":"1dffb6ddb3f648fbb425c3ef59f52f88eb97b969"},"cell_type":"markdown","source":"####  Transformation\nThe first set of plots shows both Hits and Tracks in 3D & 2D on three different coordinate systems:\n* cartesian (x,y,z)\n* cylindrical (r,phi,z) \n* polar (r,s,c) \n\nThe second set of plots also shows hits and tracks after PCA is applied to reduce dimensionality from 3D to 2D on all the three coordinates. \n\n##### See any major edits or notes at the bottom\n##### Expand code to see the functions (cartesian_to_cylindrical, cartesian_to_3d_polar, PCA, etc).\n\n"},{"metadata":{"scrolled":true,"trusted":true,"_uuid":"ffe24d7282b6cbdc48a33a671013f3053aab6953","_kg_hide-input":true,"collapsed":true},"cell_type":"code","source":"####################################################################### IMPORT\nimport os\nimport numpy as np\nimport pandas as pd\nfrom trackml.dataset import load_event\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\nfrom sklearn.decomposition import PCA \n%matplotlib inline\n####################################################################### TRANSFORMATION FUNCTIONS\n# Convert to Cylindrical \ndef cartesian_to_cylindrical(x, y, z):\n    r = np.sqrt(x**2 + y**2)\n    phi = np.arctan2(y, x)\n    z = z\n    return r, phi, z\n\n# Convert to 3D Polar\ndef cartesian_to_3d_polar(x,y,z):\n    r = np.sqrt(x**2 + y**2)\n    phi = np.arctan2(y, x)\n    s  = np.sin(phi)\n    c  = np.cos(phi)\n    return r, s, c\n\n# transform data to one of the above\ndef transform(data, transformation = \"cylindrical\"):\n    data = data.copy()\n    X = data.x.values\n    Y = data.y.values\n    Z = data.z.values\n    if transformation == \"cylindrical\":\n        data[\"x\"], data[\"y\"], data[\"z\"] = cartesian_to_cylindrical(X,Y,Z)\n        data.rename(columns={'x': 'r', 'y': 'phi', 'z': 'z'}, inplace=True)\n    elif transformation == \"polar\":\n        data[\"x\"], data[\"y\"], data[\"z\"] = cartesian_to_3d_polar(X,Y,Z)\n        data.rename(columns={'x': 'r', 'y': 's', 'z': 'c'}, inplace=True)\n    \n    x, y, z = None, None, None \n    return data\n\n\n####################################################################### Getting Tracks \ndef get_tracks(data, tracks_n=100,\n               include_zero_weights=False,\n               include_zero_ID = False, \n               coordinates = ['x','y','z'] ):\n    # cordinates= normal or polar or cylindrical \n    data = data.copy()\n    # remove zero weight particles (the one that is lower than 3 hits & and random noise)\n    if include_zero_weights == False:\n        data = data[data[\"weight\"] > 0]\n    if include_zero_ID == False:\n        data = data[data[\"particle_id\"] > 0]\n    \n    # get unique particle ids\n    track = truth.particle_id.unique()\n    \n    if tracks_n > 0 and tracks_n < track.size:\n        # select random particle ID (i.e., random tracks)\n        # would work if tracks is a dataframe\n        # tracks = tracks.loc[np.random.choice(tracks.index, size=n, replace = False)] \n        track = track[np.random.choice(track.shape[0], size=tracks_n, replace = False)]\n    \n    # get selected tracks only\n    data = data[data[\"particle_id\"].isin(track)]\n    # replace particle ids by 1,2,3,4,5,6...n ... make iterating through tracks easier  \n    data[\"particle_id\"] = pd.factorize(data[\"particle_id\"])[0]\n    # Change     \n    return data\n    \n\n\n####################################################################### Retrieval Functions \ndef get_hits(train = 1, event = 1000, sample_n = 0):\n    event_prefix = \"event00000\" + str(event)\n    hits, _, _, _ = load_event(os.path.join('../input/train_1', event_prefix))\n    _ = None\n    if sample_n > 0 and sample_n < hits.shape[0]:\n        hits = hits.sample(sample_n)\n    return hits\n\ndef get_truth(train = 1, event = 1000, sample_n = 0):\n    event_prefix = \"event00000\" + str(event)\n    _, _, _, truth = load_event(os.path.join('../input/train_1', event_prefix))\n    _ = None\n    if sample_n > 0 and sample_n < truth.shape[0]:\n        truth = truth.sample(sample_n)\n    return truth\n\n\n\n################## Constants\ntrain_file = 1\nevent = 1000\nsample_n = 80000\ntracks_n = 50\n\n####################################################################### Retrieve Data\n# Get Hits\nhits = get_hits(train = train_file, event = event, sample_n = sample_n)\nhits_c = transform(hits, \"cylindrical\")\nhits_p = transform(hits, \"polar\")\n# Get Truths\ntruth = get_truth(train = train_file, event = event, sample_n = sample_n)\ntruth.rename(columns={'tx': 'x', 'ty': 'y', 'tz': 'z'}, inplace=True)\ntruth_c = transform(truth, \"cylindrical\")\ntruth_p = transform(truth, \"polar\")\n\n\n########################################################################### PCA\n### Apply PCA to Hits (needs cleaning up)\ndef doPCA(data):\n    pca = PCA(n_components=3)\n    pca.fit(data)\n    return pca \n\np = doPCA(hits[[\"x\",\"y\",\"z\"]])\nhits[\"x_\"], hits[\"y_\"] = p.transform(hits[[\"x\",\"y\",\"z\"]])[:,0], p.transform(hits[[\"x\",\"y\",\"z\"]])[:,1]\n#print(\"hits normal score:\" , p.explained_variance_ratio_)\n\np = doPCA(hits_c[[\"r\",\"phi\",\"z\"]])\nhits_c[\"x_\"], hits_c[\"y_\"] = p.transform(hits_c[[\"r\",\"phi\",\"z\"]])[:,0], p.transform(hits_c[[\"r\",\"phi\",\"z\"]])[:,1]\n#print(\"hits cylind score:\" , p.explained_variance_ratio_)\np = doPCA(hits_p[[\"r\",\"s\",\"c\"]])\nhits_p[\"x_\"], hits_p[\"y_\"] = p.transform(hits_p[[\"r\",\"s\",\"c\"]])[:,0], p.transform(hits_p[[\"r\",\"s\",\"c\"]])[:,1]\n#print(\"hits polar score:\" , p.explained_variance_ratio_)\n\n### Apply PCA to Tracks \np = doPCA(truth[[\"x\",\"y\",\"z\"]])\ntruth[\"x_\"], truth[\"y_\"] = p.transform(truth[[\"x\",\"y\",\"z\"]])[:,0], p.transform(truth[[\"x\",\"y\",\"z\"]])[:,1]\n#print(\"truth normal score:\" , p.explained_variance_ratio_)\n\np = doPCA(truth_c[[\"r\",\"phi\",\"z\"]])\ntruth_c[\"x_\"], truth_c[\"y_\"] = p.transform(truth_c[[\"r\",\"phi\",\"z\"]])[:,0], p.transform(truth_c[[\"r\",\"phi\",\"z\"]])[:,1]\n#print(\"truth cylind score:\" , p.explained_variance_ratio_)\n\np = doPCA(truth_p[[\"r\",\"s\",\"c\"]])\ntruth_p[\"x_\"], truth_p[\"y_\"] = p.transform(truth_p[[\"r\",\"s\",\"c\"]])[:,0], p.transform(truth_p[[\"r\",\"s\",\"c\"]])[:,1]\n#print(\"truth polar score:\" , p.explained_variance_ratio_)\n\n\n# Get Tracks (after PCA)\ntracks = get_tracks(truth, tracks_n= tracks_n, coordinates = ['x','y','z'])\ntracks_c = get_tracks(truth_c, tracks_n= tracks_n, coordinates = ['r','phi','z'])\ntracks_p = get_tracks(truth_p, tracks_n= tracks_n, coordinates = ['r','s','c'])\n","execution_count":12,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"6ea2aeb8808650021c124b627889f7cdadbc675f","collapsed":true},"cell_type":"markdown","source":"### SET 1: BEFORE DIMENSIONALITY REDUCTION (PCA)\n**Plotting Hits  (3D)**\n"},{"metadata":{"scrolled":false,"trusted":true,"_uuid":"4ff4b848d6db93bfc276a365b0f2a03022102e03","_kg_hide-input":true},"cell_type":"code","source":"##### PLOTTING Normal data\nplt.figure(1, figsize=(10,10))\n#plt.figure(figsize=(15,15))\nax1 = plt.axes(projection='3d')\nax1.scatter(hits.z, hits.y, hits.x, s=5, alpha=0.5)\nax1.set_xlabel('z (mm)')\nax1.set_ylabel('y (mm)')\nax1.set_zlabel('x (mm)')\nplt.title(\"normal x, y & z\")\nplt.show()\n\nplt.figure(2, figsize=(10,10))\nax2 = plt.axes(projection='3d')\nax2.scatter(hits_c.z, hits_c.phi, hits_c.r, s=5, alpha=0.5)\nax2.set_xlabel('z (mm)')\nax2.set_ylabel('phi')\nax2.set_zlabel('r')\nplt.title(\"Cylindrical z, phi & r\")\nplt.show()\n\n# 3D polar\nplt.figure(3, figsize=(10,10))\nax2 = plt.axes(projection='3d')\nax2.scatter(hits_p.r, hits_p.s, hits_p.c, s=5, alpha=0.5)\nax2.set_xlabel('r')\nax2.set_ylabel('sin(theta)')\nax2.set_zlabel('cos(theta)')\nplt.title(\" 3D Polar r, s & c\")\nplt.show()\n","execution_count":13,"outputs":[]},{"metadata":{"_uuid":"666b3d12400a006b62aa0e5b881b3ee855ecc147"},"cell_type":"markdown","source":"**Plotting Tracks (3D)**\n"},{"metadata":{"trusted":true,"_uuid":"9d6c6d91dc6aa3166296b256bb2d01af563d6147","_kg_hide-input":true,"scrolled":false},"cell_type":"code","source":"# Plotting Tracks \nplt.figure(figsize=(10,10))\nax = plt.axes(projection='3d')\n\nfor i in range(tracks[\"particle_id\"].max()):\n    t = tracks[tracks.particle_id == i]\n    ax.plot3D(t.z, t.x, t.y)\nax.set_xlabel('z (mm)')\nax.set_ylabel('x (mm)')\nax.set_zlabel('y (mm)')\n# These two added to widen the 3D space\nax.scatter(3000,3000,3000, s=0)\nax.scatter(-3000,-3000,-3000, s=0)\nplt.show()\n\nplt.figure(figsize=(10,10))\nax = plt.axes(projection='3d')\n\nfor i in range(tracks_c[\"particle_id\"].max()):\n    t = tracks_c[tracks_c.particle_id == i]\n    ax.plot3D(t.z, t.phi, t.r)\nax.set_xlabel('z (mm)')\nax.set_ylabel('phi')\nax.set_zlabel('r')\n# These two added to widen the 3D space\nplt.show()\n\n\nplt.figure(figsize=(10,10))\nax = plt.axes(projection='3d')\n\nfor i in range(tracks_p[\"particle_id\"].max()):\n    t = tracks_p[tracks_p.particle_id == i]\n    ax.plot3D(t.r, t.s, t.c)\nax.set_xlabel('r')\nax.set_ylabel('s')\nax.set_zlabel('c')\nplt.show()","execution_count":14,"outputs":[]},{"metadata":{"_uuid":"11ac75be65e552e4a6ac681f70174b89a96e0086"},"cell_type":"markdown","source":"**Plotting Hits  (2D)**\n\nNote: Some axes combinations not shown for:\n\n* x & z has the same visiualization as y & z\n* r & z has the same visiualization as phi & r (flipped)\n* r & c has the same visiulaization as r & s (flipped)"},{"metadata":{"scrolled":false,"trusted":true,"_uuid":"16693a7dbc0a10374bea9a6ccdf813d3fa60ec7a","_kg_hide-input":true},"cell_type":"code","source":"##### PLOT 2D x,y & phi,r\n\nplt.figure(4, figsize=(15,15))\nplt.subplot(321)\nplt.scatter(hits.y, hits.x, s=5, alpha=0.5)\nplt.xlabel('y (mm)')\nplt.ylabel('x (mm)')\n\nplt.subplot(322)\nplt.scatter(hits.y, hits.z, s=5, alpha=0.5)\nplt.xlabel('y')\nplt.ylabel('z')\n# x & z looks same as y & z\n\nplt.subplot(323)\nplt.scatter(hits_c.phi, hits_c.r, s=5, alpha=0.5)\nplt.xlabel('phi')\nplt.ylabel('r')\n# r & z looks same as phi & r\n\nplt.subplot(324)\nplt.scatter(hits_c.phi, hits_c.z, s=5, alpha=0.5)\nplt.xlabel('phi')\nplt.ylabel('z')\n\nplt.subplot(325)\nplt.scatter(hits_p.r, hits_p.s, s=5, alpha=0.5)\nplt.xlabel('r')\nplt.ylabel('s')\n# r & c looks same as r & s\n\nplt.subplot(326)\nplt.scatter(hits_p.s, hits_p.c, s=5, alpha=0.5)\nplt.xlabel('s')\nplt.ylabel('c')\nplt.show()","execution_count":15,"outputs":[]},{"metadata":{"_uuid":"7c1a41dd021810aaa3cb2ca1bd08b06083fe61d3"},"cell_type":"markdown","source":"**Plotting Hits with Tracks (2D)**\n\nNote: each unique color is a track\n"},{"metadata":{"trusted":true,"_uuid":"1752a47c2bf32bc05c9b6b446e3f275ec11660ed","_kg_hide-input":true},"cell_type":"code","source":"### Plot Graphs above with Tracks, each unique color is a track\nplt.figure(4, figsize=(20,20))\nplt.subplot(321)\nplt.scatter(hits.y, hits.x, s=5, alpha=0.5)\nfor i in range(tracks[\"particle_id\"].max()):\n    t = tracks[tracks.particle_id == i]\n    plt.scatter(t.y, t.x)\nplt.xlabel('y (mm)')\nplt.ylabel('x (mm)')\n\nplt.subplot(322)\nplt.scatter(hits.y, hits.z, s=5, alpha=0.5)\nfor i in range(tracks[\"particle_id\"].max()):\n    t = tracks[tracks.particle_id == i]\n    plt.scatter(t.y, t.z)\nplt.xlabel('y')\nplt.ylabel('z')\n# x & z looks same as y & z\n\nplt.subplot(323)\nplt.scatter(hits_c.phi, hits_c.r, s=5, alpha=0.5)\nfor i in range(tracks_c[\"particle_id\"].max()):\n    t = tracks_c[tracks_c.particle_id == i]\n    plt.scatter(t.phi, t.r)\nplt.xlabel('phi')\nplt.ylabel('r')\n# r & z looks same as phi & r\n\nplt.subplot(324)\nplt.scatter(hits_c.phi, hits_c.z, s=5, alpha=0.5)\nfor i in range(tracks_c[\"particle_id\"].max()):\n    t = tracks_c[tracks_c.particle_id == i]\n    plt.scatter(t.phi, t.z)\nplt.xlabel('phi')\nplt.ylabel('z')\n\nplt.subplot(325)\nplt.scatter(hits_p.r, hits_p.s, s=5, alpha=0.5)\nfor i in range(tracks_p[\"particle_id\"].max()):\n    t = tracks_p[tracks_p.particle_id == i]\n    plt.scatter(t.r, t.s)\nplt.xlabel('r')\nplt.ylabel('s')\n# r & c looks same as r & s\n\nplt.subplot(326)\nplt.scatter(hits_p.s, hits_p.c, s=5, alpha=0.5)\nfor i in range(tracks_p[\"particle_id\"].max()):\n    t = tracks_p[tracks_p.particle_id == i]\n    plt.scatter(t.s, t.c)\nplt.xlabel('s')\nplt.ylabel('c')\nplt.show()","execution_count":16,"outputs":[]},{"metadata":{"_uuid":"5d091cab48799a951f2a7918f947d99ce4c63e21"},"cell_type":"markdown","source":"### SET 2: AFTER DIMENSIONALITY REDUCTION (PCA)\n**Plotting Hits (3D --> 2D)**"},{"metadata":{"trusted":true,"_uuid":"eee26e48961e448da61686146a3ed0ededb6a210","_kg_hide-input":true,"scrolled":false},"cell_type":"code","source":"## Plot hits after dimensionallity reduction \nplt.figure(5, figsize=(10,10))\nplt.scatter(hits.x_, hits.y_, s=5, alpha=0.5)\nplt.title(\"Cartesian reduced\")\nplt.show()\n\n\nplt.figure(6, figsize=(10,10))\nplt.scatter(hits_c.x_, hits_c.y_, s=5, alpha=0.5)\n\nplt.title(\"Cylindrical Reduced\")\nplt.show()\n\nplt.figure(7, figsize=(10,10))\nplt.scatter(hits_p.x_, hits_p.y_, s=5, alpha=0.5)\n\nplt.title(\"Polar reduced\")\nplt.show()\n","execution_count":17,"outputs":[]},{"metadata":{"_uuid":"aea02cf59eb891d4543cf4411e438ae9f6a37097"},"cell_type":"markdown","source":"**Plotting Hits with Tracks**"},{"metadata":{"scrolled":false,"trusted":true,"_uuid":"c202026ca890ebce4b99ab32592ded40189a4fc5","_kg_hide-input":true},"cell_type":"code","source":"# Plot graphs above with tracks after PCA applied to both\nplt.figure(7, figsize=(10,10))\nplt.scatter(hits.x_, hits.y_, s=5, alpha=0.5)\n\nfor i in range(tracks[\"particle_id\"].max()):\n    t = tracks[tracks.particle_id == i]#.sort_values(\"y2\")\n    plt.plot(t.x_, t.y_)\n\nplt.title(\"Reduced Tracks & Hits (cartesian)\")\nplt.show()\n\n#tracks.sort_values(\"y2\")\nplt.figure(8, figsize=(10,10))\nplt.scatter(hits_c.x_, hits_c.y_, s=5, alpha=0.5)\n\nfor i in range(tracks_c[\"particle_id\"].max()):\n    t = tracks_c[tracks_c.particle_id == i]#.sort_values(\"y2\")\n    plt.scatter(t.x_, t.y_)\n\nplt.title(\"Reduced Tracks & Hits (cylindrical)\")\nplt.show()\n\nplt.figure(9, figsize=(10,10))\nplt.scatter(hits_p.x_, hits_p.y_, s=5, alpha=0.5)\nfor i in range(tracks_p[\"particle_id\"].max()):\n    t = tracks_p[tracks_p.particle_id == i]#.sort_values(\"y2\")\n    plt.scatter(t.x_, t.y_)\nplt.title(\"Reduced Tracks & Hits (polar)\")\nplt.show()\n","execution_count":18,"outputs":[]},{"metadata":{"_uuid":"c6987573e35acaa67b9d322d3b54cc5d355341ae"},"cell_type":"markdown","source":"#### EDITs\n##### what I thought was a problem\nBefore I showed that I had a problem where the tracks were tilted/shifted and did not match the hits data (see example in this discussion: [tilted tracks](https://www.kaggle.com/c/trackml-particle-identification/discussion/57931).  The problem was that I applied PCA after I selected the tracks from the truth dataset (oops!). Now, both transformation and PCA are first applied to the truth dataset, then tracks are selected. As you can see, more accurate results. \n\n##### Notes\nI'll keep updating this kernel in the future. \n\nTo try larger number of tracks, change the \"constants\" in the Functions cell (first code cell). \n\n#### Credits \nSome of the functions and concepts were taken from the following kernels/discussions:  \n\nhttps://www.kaggle.com/c/trackml-particle-identification/discussion/57643 by Heng CherKeng\n\nhttps://www.kaggle.com/mikhailhushchyn/hough-transform by Mikhail Hushchyn \n\nhttps://www.kaggle.com/jbonatt/trackml-eda-etc/notebook by Joshua Bonatt\n"}],"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}