{"cells":[{"metadata":{"_uuid":"22ced56347fedd0d1919718e4a822c941c583655"},"cell_type":"markdown","source":"# TrackML Challenge Kernel for Kaggle Course "},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true,"_uuid":"32c7d658276d38cdff24979ac7ce4134d9bf5f5a"},"cell_type":"code","source":"from IPython.core.display import display, HTML\ndisplay(HTML(\"\"))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fb08a57477ec149f3e138ea43fa06631786c47a1"},"cell_type":"markdown","source":"This kernel outlines the process taken towards the Kaggle Course offered by AI Academy at Aalto University."},{"metadata":{"_uuid":"40bd10f6f4e9dc7ac753015ec94c222050a3f502"},"cell_type":"markdown","source":"**Due to the nature of the TrackML Challenge, this notebook will deviate from the requirements to clearly present the problem and associated soltuions that were tested. Please read the rest of the document and evaluate it accordingly.**"},{"metadata":{"_uuid":"ae3a3ef505f5be0adfce772e6deeff5661d64cf3"},"cell_type":"markdown","source":"I will also be extensively using content from public kernels that explain the data and problem much better than I ever can hope to."},{"metadata":{"_uuid":"7a38f2d39b78c6e2b3eb8bcddbdf4220326bc20f"},"cell_type":"markdown","source":"## Problem Overview\nThe TrackML Challenge was put forth by CERN to find solutions to the well known particle path reconstruction problem. The problem will be explained in stages. "},{"metadata":{"_uuid":"92cc9a05e6429c9a76a0ab89db8695f161ee7ce2"},"cell_type":"markdown","source":"Firstly an into to CERN, they have large particle colliders that speed up elementary particles to high speeds and then collide them and find inferences from the data generated from the collision. \n\nEach collision results in a myriad of particles that are fleeting in nature and then decompose into even more particles. The detectors in CERN are **silicon slabs** that are able to detect small changes that correspond to when a particle hits the detector. The detectors are arranged in a manner as shown in the figure below.
\n\n
\nThis gives us an view of how the detectors are arranged in 3 dimensions, which will come in handy when doing visualizations.A very simple analogy to understand the problem:
\n\n\n**Consider that we shoot a bullet through multiple vertical cardboard cutouts that are separated horizonatally. Now after the bullet has passed through them we are left with the holes in the cardboard, now our job is to reconstruct the path of the bullet given the holes detected in the cardboard**. \n\nThe following image may provide an intuition behind the problem, where the red dots are the holes left behind and the arc is the possible part of the bullet.\n
\n"},{"metadata":{"_uuid":"52c109e9d605c283b05383a46d5c262c73911954"},"cell_type":"markdown","source":"## Data\nUpon getting an understanding of the problem, it may be easy to ask, **what is so hard about this?**. A valid question, this brings us to the complexity of the problem and the difficulty in directly throwing any well known ML algorithm without thinking. This is due to the enormous amount of data that is generated and given to us.
\nConsider each collision as an **event**. Each **event** creates roughly **10,000** particles. Each particle passes through an average of 10 silicon detectors. On average each **event** has **100,000** data points. The whole dataset has **8,850** events leading to a grand total of **70GB** of training data. \n\n* __Hits__ are given as x, y, z positions and the detector groups they interacted with: volume, layer, module in increasing granularity (volume big, module small). \n* For each hit, we're also also given the __cells__ position and charge that did the detecting. These are the smallest unit of detector group, even finer than module. \n* For each event we're also given __particle__ information x, y, z momentums, and charge, and number of detectors hit from that particle. \n* Finally, we're given __truth__, telling us for certain which hits belong to which particles. Maybe this is kind of like our \"y\", and the other features are kind of like our \"X\".\n\nThis is all explained in more detail below, later.\n\nFirst let's see some visualizations"},{"metadata":{"_uuid":"cb3b7bb3de31d0f6ec72dc46068fc494da41174c","trusted":true},"cell_type":"code","source":"from trackml.dataset import load_event\nimport os\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport scipy.stats as stats\nimport numpy as np\nimport pandas as pd\nfrom IPython.display import Image,display\nfrom trackml.dataset import load_dataset\nfrom trackml.score import score_event","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b496553a9ba19d683da118386541744cb54ec2c1","trusted":true},"cell_type":"code","source":"import matplotlib as mpl\nfrom mpl_toolkits.mplot3d import Axes3D\n\nmpl.rcParams['figure.dpi'] = 300","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8ad69c52552b4fa266b46baa28092bd922ad0fd8"},"cell_type":"markdown","source":"## Load one event to create visualizations"},{"metadata":{"_uuid":"2faff7c0d145f06e38302621024efb2af972bd05","trusted":true},"cell_type":"code","source":"event_prefix = 'event000001000'\npath_to_train = \"../input/trackml-particle-identification/train_1\"\npath_to_test = \"../input/trackml-particle-identification/test/\"\nhits, cells, particles, truth = load_event(os.path.join('../input/trackml-particle-identification/train_1/', event_prefix))\nhits.describe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"76f9c74d9b6277164a7e1e2bd68f7611fd3bade6"},"cell_type":"markdown","source":"__Credits to : @Joshua Bonatt
\n__hits__ \nThe hits file contains the following values for each hit/entry: \n* __hit_id__: numerical identifier of the hit inside the event.\n* __x, y, z__: measured x, y, z position (in millimeter) of the hit in global coordinates.\n* __volume_id__: numerical identifier of the detector group.\n* __layer_id__: numerical identifier of the detector layer inside the group.\n* __module_id__: numerical identifier of the detector module inside the layer. \n\nThe volume/layer/module id could in principle be deduced from x, y, z. They are given here to simplify detector-specific data handling.\n\n__cells__ \nThe cells file contains the constituent active detector cells that comprise each hit. The cells can be used to refine the hit to track association. A cell is the smallest granularity inside each detector module, much like a pixel on a screen, except that depending on the volume_id a cell can be a square or a long rectangle. It is identified by two channel identifiers that are unique within each detector module and encode the position, much like column/row numbers of a matrix. A cell can provide signal information that the detector module has recorded in addition to the position. Depending on the detector type only one of the channel identifiers is valid, e.g. for the strip detectors, and the value might have different resolution. \n* __hit_id__: numerical identifier of the hit as defined in the hits file.\n* __ch0, ch1__: channel identifier/coordinates unique within one module.\n* __value__: signal value information, e.g. how much charge a particle has deposited.\n\n__particles__ \nThe particles files contains the following values for each particle/entry: \n* __particle_id__: numerical identifier of the particle inside the event.\n* __vx, vy, vz__: initial position or vertex (in millimeters) in global coordinates.\n* __px, py, pz__: initial momentum (in GeV/c) along each global axis.\n* __q__: particle charge (as multiple of the absolute electron charge).\n* __nhits__: number of hits generated by this particle.\n\n__truth__ \nThe truth file contains the mapping between hits and generating particles and the true particle state at each measured hit. Each entry maps one hit to one particle. \n* __hit_id__: numerical identifier of the hit as defined in the hits file.\n* __particle_id__: numerical identifier of the generating particle as defined in the particles file. A value of 0 means that the hit did not originate from a reconstructible particle, but e.g. from detector noise.\n* __tx, ty, tz__ true intersection point in global coordinates (in millimeters) between the particle trajectory and the sensitive surface.\n* __tpx, tpy, tpz__ true particle momentum (in GeV/c) in the global coordinate system at the intersection point. The corresponding vector is tangent to the particle trajectory at the intersection point.\n* __weight__ per-hit weight used for the scoring metric; total sum of weights within one event equals to one.\n\n"},{"metadata":{"_uuid":"7ff0b016cc5c5d9ed75fd67b1c947614cc0376db"},"cell_type":"markdown","source":"## Hits"},{"metadata":{"_uuid":"c0a39df7022de4ec060f457f2176ad90413338b8","trusted":true},"cell_type":"code","source":"hits.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"be92619c7198df23b0c56a92ae65d276436bbac8","trusted":true},"cell_type":"code","source":"hits.tail()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0c5f109cac97d330b9e5ae3b9df86b965e74eb56","trusted":true},"cell_type":"code","source":"g = sns.jointplot(hits.x, hits.y, s=1, height=12)\ng.ax_joint.cla()\nplt.sca(g.ax_joint)\n\nvolumes = hits.volume_id.unique()\nfor volume in volumes:\n v = hits[hits.volume_id == volume]\n plt.scatter(v.x, v.y, s=3, label='volume {}'.format(volume))\n\nplt.xlabel('X (mm)')\nplt.ylabel('Y (mm)')\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e22f3a52fdf85e6ef133f93932c0103c357e8cf4"},"cell_type":"markdown","source":"**The solid green core in the middle is more concentric tigtly packed detectors. Let's zoom into them and see them.**"},{"metadata":{"_uuid":"72dd9f5ee05e288bd0f89ece517fcec64008c2e6","trusted":true},"cell_type":"code","source":"g = sns.jointplot(hits.x, hits.y, s=1, height=12)\ng.ax_joint.cla()\nplt.sca(g.ax_joint)\n\nvolumes = hits.volume_id.unique()\nfor volume in volumes:\n v = hits[hits.volume_id == volume]\n plt.scatter(v.x, v.y, s=3, label='volume {}'.format(volume))\n\nplt.xlim([-200,200])\nplt.ylim([-200,200])\nplt.xlabel('X (mm)')\nplt.ylabel('Y (mm)')\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"14c88ead5099fb1c0d34651d62742d64195e14c8"},"cell_type":"markdown","source":"That scattering we see (the blue and green dots) in between the orange rings (detectors) are just events detected in the vertical caps. Upon Removing them we see only the detections on the X and Y planes. We set a value of **|z|<=300**, these are the locations of the first set of vertical detectors, hence we remove all blue and green events.\n\nUpon even closer observation, each orange slab is actually a detector that has been saturated with hits to become visible"},{"metadata":{"_uuid":"a1ea176ea8c0c74dc58ff331ef596134243960d9","trusted":true},"cell_type":"code","source":"g = sns.jointplot(hits.x, hits.y, s=1, height=12)\ng.ax_joint.cla()\nplt.sca(g.ax_joint)\n\nvolumes = hits.volume_id.unique()\nfor volume in volumes:\n v = hits[hits.volume_id == volume]\n plt.scatter(v[v.z.abs()<300].x, v[v.z.abs()<300].y, s=3, label='volume {}'.format(volume))\n\nplt.xlim([-200,200])\nplt.ylim([-200,200])\nplt.xlabel('X (mm)')\nplt.ylabel('Y (mm)')\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a1e6450d7828598d0a6180d8cecc2e1a1fde16f3"},"cell_type":"markdown","source":"## A view along Y and Z axes\n\nNow we can see the orange detectors from before as the concentric rings from above. We can now also see the position of the vertical detectors along the Z axis in blue and green near the origin. "},{"metadata":{"_uuid":"d18d1c19a9c53f3ec41b1b8ffb89927c5941d4ee","trusted":true},"cell_type":"code","source":"volumes = hits.volume_id.unique()\nplt.figure(figsize=(24,8))\n\nfor volume in volumes:\n v = hits[hits.volume_id == volume]\n plt.scatter(v.z, v.y, s=1, label='volume {}'.format(volume))\n\nplt.xlabel('Z (mm)')\nplt.ylabel('Y (mm)')\nplt.ylim([-200,200])\nplt.legend()\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8a5279ef847b202372fd08c741fcd88e7c7a0f4a"},"cell_type":"markdown","source":"## Undertanding Concentric arrangement\n\nIf we take apart the concentric orange rings and then plot them according to the volume id that denotes the location of the detectors. \n\n__Volume 7__ is to the __left__, __Volume 8__ is the __middle__ (concentric in X axis) and __Volume 9__ in the __right__. Similarly the next 3 volume ids also follow the same order.\n\nWe see that __Volume 8,13,17__ are the ones that have concentric rings of detectors along the X axis"},{"metadata":{"_uuid":"e862c1d76bc0ae71600798482a225a1fa7138d64","trusted":true},"cell_type":"code","source":"g = sns.FacetGrid(hits, col=\"volume_id\",aspect=0.5,sharey=True,sharex=True)\ng=(g.map(plt.scatter,\"x\", \"y\",s=0.01))\ng.set_axis_labels('')\nfor i in range(len(hits.volume_id.unique())):\n g.facet_axis(0,i).axis('off')\nplt.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"858bd5db710a11589a05f247c678bb6ee3828bd9"},"cell_type":"markdown","source":"### A 3D plot of volume ids of the whole detector"},{"metadata":{"_uuid":"a4ab2b3a5a42c123458c42674cfae16df300b4cb","trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(12, 12))\nax = fig.add_subplot(111, projection='3d')\nsample = hits.sample(10000)\nfor volume in volumes:\n v = sample[sample.volume_id == volume]\n ax.scatter(v.z, v.x, v.y, s=1, label='volume {}'.format(volume), alpha=0.5)\nax.set_title('Hit Locations')\nax.set_xlabel('Z (millimeters)')\nax.set_ylabel('X (millimeters)')\nax.set_zlabel('Y (millimeters)')\n#ax.view_init(elev=10,azim=10)\nax.legend()\n\nax.scatter(3000,3000,3000, s=0)\nax.scatter(-3000,-3000,-3000, s=0)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c923067645d77af65b866f0e67d0f70dd566bdc7"},"cell_type":"markdown","source":"### A 3D plot for Layer IDs \n\nLayer IDs are different from volume and as seen in the plot, the increasing Layer ID indicates how close to the centre they were."},{"metadata":{"_uuid":"4858d831110541a9309d88a0a010badc1ec8c825","trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(12, 12))\nax = fig.add_subplot(111, projection='3d')\nlayers = hits.layer_id.unique()\n\nfor layer in layers:\n v = sample[sample.layer_id == layer]\n ax.scatter(v.z, v.x, v.y, s=1, label='layer {}'.format(layer), alpha=0.5)\nax.set_title('Hit Locations')\nax.set_xlabel('Z (millimeters)')\nax.set_ylabel('X (millimeters)')\nax.set_zlabel('Y (millimeters)')\nax.scatter(3000,3000,3000, s=0)\nax.scatter(-3000,-3000,-3000, s=0)\nax.legend()\n#ax.view_init(elev=10,azim=10)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"bf506875ecd729d44195569861a6dbbed1e775e1"},"cell_type":"markdown","source":"# Particles\nThe particles files contains the following values for each particle/entry: \n* __particle_id__: numerical identifier of the particle inside the event.\n* __vx, vy, vz__: initial position or vertex (in millimeters) in global coordinates.\n* __px, py, pz__: initial momentum (in GeV/c) along each global axis.\n* __q__: particle charge (as multiple of the absolute electron charge).\n* __nhits__: number of hits generated by this particle."},{"metadata":{"_uuid":"ea5ea1a54c3b69f3ad94284f2dbd0cc837cf1337","trusted":true},"cell_type":"code","source":"particles.describe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8fec95568cafdbe9925c0536eb399ba3424c5243"},"cell_type":"markdown","source":"# Truth\nThe truth file contains the mapping between hits and generating particles and the true particle state at each measured hit. Each entry maps one hit to one particle. \n* __hit_id__: numerical identifier of the hit as defined in the hits file.\n* __particle_id__: numerical identifier of the generating particle as defined in the particles file. A value of 0 means that the hit did not originate from a reconstructible particle, but e.g. from detector noise.\n* __tx, ty, tz__ true intersection point in global coordinates (in millimeters) between the particle trajectory and the sensitive surface.\n* __tpx, tpy, tpz__ true particle momentum (in GeV/c) in the global coordinate system at the intersection point. The corresponding vector is tangent to the particle trajectory at the intersection point.\n* __weight__ per-hit weight used for the scoring metric; total sum of weights within one event equals to one."},{"metadata":{"_uuid":"4c97ab6d4753755f934079d84137912829065006","trusted":true},"cell_type":"code","source":"truth.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"565eb27cf3c12c0daab24634f8925f03901cc324"},"cell_type":"markdown","source":"# True Particle Trajectories\n\nThis plot shows what we wish to achieve, the true paths of the particles given the hits and particels information. This is only a plot of 100th particle track"},{"metadata":{"_uuid":"aeac438f4569d67381a82ad01ae664f3fc16d6d7","trusted":true},"cell_type":"code","source":"# Get every 100th particle\n\ntracks = truth.particle_id.unique()[1::100]\n\nplt.figure(figsize=(15,15))\nax = plt.axes(projection='3d')\nfor track in tracks:\n t = truth[truth.particle_id == track]\n ax.plot3D(t.tz, t.tx, t.ty)\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()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"385984c682c84aa0befc4279b32b106aeffba39e"},"cell_type":"markdown","source":"## Some Physics (Unfortunately) for Feature Engineering and Selection\n\nBy utilizing some particle physics we can come up with ways to convert the hit locations into coordinates that may be easier to group. This rests on the idea that a charged particle moving in a magnetic field will assume a helical orbit as seen in the image below. \n
\n\n
\n\nAs the slide from the High Energy Physis (HEP) domain suggests by converting our cartesian coordinates to other coordinates we can achieve a better grouping of the __hits__ that are detected.\n\nThe following involve some partice physics calculations that show how we arive at the transformation that is the easiest to implement. "},{"metadata":{"_uuid":"4b2d3f400c96f100fa9f337b5b429173a79b4e3d"},"cell_type":"markdown","source":"## Coordinate Transformations\n**Credits: @Mark JD Hamilton**
\n**Can Be Skipped**
\nSome notations to keep in common \nWe fix the following notation:\n1. Particle rest mass $m_0$ and electric charge $q$.\n2. $m=\\gamma m_0$ with relativistic factor $\\gamma$ (considered constant).\n3. Coordinate system $(x,y,z)$.\n4. Magnetic field $\\vec{B}=Be_z$.\n5. Initial coordinates $x_0=y_0=z_0=0$ (for simplicity, see comment below).\n6. Initial velocities $v_{x0}, v_{y0}, v_{z0}$ and initial momenta $p_{x0}, p_{y0}, p_{z0}$ with $p_{i0}=m v_{i0}$ for $i=x,y,z$.\n7. Longitudinal momentum $p_\\parallel = p_{z0}$, transversal momentum \n\\begin{equation*}\np_\\perp = \\sqrt{p_{x0}^2+p_{y0}^2},\n\\end{equation*}\nand total momentum \n\\begin{equation*}\np=\\sqrt{p_\\parallel^2+p_\\perp^2}.\n\\end{equation*}\n\nWe define:\n\\begin{align*}\n\\omega &= \\frac{qB}{m}\\\\\nR&= \\frac{p_\\perp}{qB}\n\\end{align*}\nand an angle $\\phi_0$ by\n\\begin{align*}\n\\cos\\phi_0&=-\\frac{p_{y0}}{p_\\perp}\\\\\n\\sin\\phi_0&=\\frac{p_{x0}}{p_\\perp}.\n\\end{align*}\nThe system of differential equations\n\\begin{align*}\nm\\ddot{x} &= qB\\dot{y}\\\\\nm\\ddot{y} &= -qB\\dot{x}\\\\\nm\\ddot{z} &= 0\n\\end{align*}\nhas the helix solution\n\\begin{align*}\nx(t) &= R(\\cos(\\phi_0-\\omega t)-\\cos(\\phi_0))\\\\\ny(t) &= R(\\sin(\\phi_0-\\omega t)-\\sin(\\phi_0))\\\\\nz(t) &= R\\frac{p_\\parallel}{p_\\perp}\\omega t.\n\\end{align*}\nIf the initial coordinates are non-zero, we just add $x_0$, $y_0$ and $z_0$ to these expressions.\n\nProjected onto the $x$-$y$-plane, the helix forms a circle of radius $|R|$ with center at the point \n\\begin{equation*}\n-R(\\cos\\phi_0,\\sin\\phi_0).\n\\end{equation*}\nThe motion in $z$-direction is of constant velocity."},{"metadata":{"_uuid":"4dde932eb875a982f7bc5129efb7b64ac44196d5"},"cell_type":"markdown","source":"Using\n\\begin{align*}\n\\cos^2(\\alpha)+\\sin^2(\\alpha)&=1\\\\\n\\cos(\\alpha)\\cos(\\beta)+\\sin(\\alpha)\\sin(\\beta)&=\\cos(\\alpha-\\beta),\n\\end{align*}\nwe have\n\\begin{align*}\nx^2+y^2 &= 2R^2(1-\\cos(\\omega t))\\\\\n&\\approx R^2\\left(\\omega^2t^2-\\frac{1}{12}\\omega^4t^4 \\right),\n\\end{align*}\nwhere the approximation in the second line follows from Taylor's formula and is valid for small $\\omega t$.\n\nWe define a new coordinate\n\\begin{equation*}\nz_2=\\frac{z}{\\sqrt{x^2+y^2}}=\\pm\\frac{p_\\parallel}{p_\\perp}\\frac{\\omega t}{\\sqrt{2(1-\\cos(\\omega t))}},\n\\end{equation*}\nwhich is independent of $R$ and $\\phi_0$ and the sign $\\pm$ is the sign of the charge $q$.\n\nFor small $\\omega t$ we get\n\\begin{align*}\nz_2&\\approx \\frac{p_\\parallel}{p_\\perp}\\frac{1}{\\sqrt{1-\\frac{1}{12}\\omega^2t^2}}\\\\\n&\\approx\\frac{p_\\parallel}{p_\\perp}\\left(1+\\frac{1}{24}\\omega^2t^2\\right).\n\\end{align*}\nHence up to first order in $\\omega t$ the coordinate $z_2$ is constant for each helix and given by\n\\begin{equation*}\nz_2\\approx \\frac{p_\\parallel}{p_\\perp}.\n\\end{equation*}"},{"metadata":{"_uuid":"85cab2dc0cc115397443b3e6537ee5ef2d2c6d3c"},"cell_type":"markdown","source":"We set\n\\begin{align*}\nx_2 &= \\frac{x}{\\sqrt{x^2+y^2+z^2}}\\\\\ny_2 &= \\frac{y}{\\sqrt{x^2+y^2+z^2}}.\n\\end{align*}\nThen the expression\n\\begin{align*}\nx_2^2+y_2^2 &= \\frac{x^2+y^2}{x^2+y^2+z^2}\\\\\n&=\\frac{1}{1+z_2^2}\n\\end{align*}\nis independent of $R$ and $\\phi_0$.\n\nFor small $\\omega t$, where $z_2\\approx \\frac{p_\\parallel}{p_\\perp}$, the point $(x_2,y_2)$ lies approximately on a circle of radius\n\\begin{equation*}\n\\frac{1}{\\sqrt{1+\\frac{p_\\parallel^2}{p_\\perp^2}}} = \\frac{p_\\perp}{\\sqrt{p_\\parallel^2+p_\\perp^2}} = \\frac{p_\\perp}{p},\n\\end{equation*}\nwhere $p$ is the total momentum.\n\nUsing Taylor expansion for small $\\omega t$ we can calculate\n\\begin{align*}\nx(t)&\\approx R\\sin(\\phi_0)\\omega t\\\\\ny(t)&\\approx -R\\cos(\\phi_0)\\omega t.\n\\end{align*}\nThis implies\n\\begin{align*}\n\\lim_{t\\rightarrow 0}x_2(t)&=\\sin(\\phi_0)\\frac{p_\\perp}{p}\\\\\n\\lim_{t\\rightarrow 0}y_2(t)&=-\\cos(\\phi_0)\\frac{p_\\perp}{p}.\n\\end{align*}\nWe see that the (idealized) initial point $(x_2,y_2)$ on the circle of radius $\\frac{p_\\perp}{p}$ depends on the angle $\\phi_0$, even though $x_0=y_0=0$. Hence in the coordinates $x_2,y_2$ the different helices get separated, depending on the angle $\\phi_0$."},{"metadata":{"_uuid":"3355daf826554f4e2488644e7691e88707db30f0"},"cell_type":"markdown","source":"## The Angular Track Transformation \n\nThe final transformation that is essential is as follows :\n$$ \\{x,y,z\\} \\rightarrow \\{x_2,y_2,z_2\\}`$$\nWhere \n$$ \\{x_2,y_2,z_2\\} \\rightarrow \\left\\{ \\frac{x}{\\sqrt{x^2+y^2+z^2}},\\frac{y}{\\sqrt{x^2+y^2+z^2}},\\frac{z}{\\sqrt{x^2+y^2}} \\right\\} $$\n\nThe slide below also mentions the $\\{x,y\\}$ transformation and also lists it's benifits. \n
\n\n
"},{"metadata":{"_uuid":"77313b543ed01ab1231192e66c397129bd6cb2eb"},"cell_type":"markdown","source":"## Spherical coordinates\n\nTo compare the usefulness of the above __Angulr Track Transform__,I also compared the effect of tranforming the Cartesian coordintes to the well known __Spherical Coordinate__ system. The maths behind that is as follows \n\nSpherical coordinates $(r,\\phi,\\theta)$, with $r\\geq 0, \\phi\\in[0,2\\pi], \\theta\\in[0,\\pi]$, are given by\n\\begin{align*}\nx &= r\\sin\\theta\\cos\\phi\\\\\ny &= r\\sin\\theta\\sin\\phi\\\\\nz &= r\\cos\\theta.\n\\end{align*}\nUsing our calculations above we can describe the helix by\n\\begin{align*}\nr(t)&=|R|\\sqrt{\\left(\\frac{p_\\parallel}{p_\\perp}\\right)^2\\omega^2t^2+4\\sin^2\\left(\\frac{\\omega t}{2}\\right)}\\\\\n\\phi(t)&=\\phi_0-\\frac{1}{2}\\omega t-\\frac{\\pi}{2}\\\\\n\\theta(t)&=\\arctan\\left(\\frac{p_\\perp}{p_\\parallel}\\frac{2\\sin\\left(\\frac{\\omega t}{2}\\right)}{\\omega t}\\right).\n\\end{align*}\nUp to terms of second order in $\\omega t$ we get\n\\begin{align*}\nr(t)&\\approx \\frac{p}{p_\\perp}|R\\omega t|\\\\\n\\phi(t)&=\\phi_0-\\frac{1}{2}\\omega t-\\frac{\\pi}{2}\\\\\n\\theta(t)&\\approx\\arctan\\left(\\frac{p_\\perp}{p_\\parallel}\\right).\n\\end{align*}\nIn particular, the angle $\\theta$ is approximately constant.\n"},{"metadata":{"_uuid":"56de53df6c7e9c4c558b6438a7669780a81c505f"},"cell_type":"markdown","source":"## Implementing Feature Engineering and Helix Transformations"},{"metadata":{"_uuid":"9c48ccd75a4b1722439bc9667c5a19cffbfc19de"},"cell_type":"markdown","source":"### The Angular Track Transformation\nAs described it is implemented in the following function. It is called dbscan_preprocess as it will later be used in the clustering step to get the scores."},{"metadata":{"_uuid":"9b4358a53c1dc619ee37df4aac65bd1735d684d2","trusted":true},"cell_type":"code","source":"def dbscan_preprocess(hits):\n \n x = hits[0]\n y = hits[1]\n z = hits[2]\n\n r = np.sqrt(x**2 + y**2 + z**2)\n x2 = x/r\n y2 = y/r\n\n r = np.sqrt(x**2 + y**2)\n z2 = z/r\n \n return np.vstack((x2,y2,z2))\n ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"85246795530bb95c434e442588ada2b92e4a9816"},"cell_type":"markdown","source":"### The Spherical Coordinate Trasnform\nThe below function implements the described spherical transform"},{"metadata":{"_uuid":"4d97665a6be24f0de7bc9332ab0ec880846a2dd4","trusted":true},"cell_type":"code","source":"def 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))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3b66f8be0511742c36605c7f0c4b7fdfe51c3fba"},"cell_type":"markdown","source":"## Plot Cartesian Coordinates Particle Trajectories\n\nThe following image is the paths of the 20 particles that have the highest truth weights, i.e. most likely to be correct. As we see many particles overlapp and some have bizzare paths. \n\n**Important Point:** If we notice clearly we see that not all particles originate from $\\{0,0,0\\}$, this is because particles are generated by the decay of other particles, hence the beginning of the particle path can be anywhere."},{"metadata":{"_uuid":"e3d2aaec312aa9889b36e8dba633fbc159dd4e8c","trusted":true},"cell_type":"code","source":"# Get particle id with highest weights\nNUM_PARTICLES = 20\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 p_traj_list.append(p_traj)\n \n# Convert to spherical coordinate.\nrtp_list = []\nang_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 ang = dbscan_preprocess(xyz).transpose()\n ang_df = pd.DataFrame(ang, columns=('x2', 'y2', 'z2'))\n ang_list.append(ang_df)\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.tz,\n ys=p_traj.tx,\n zs=p_traj.ty,\n marker='o')\nax.set_xlabel('Z (mm) -- Detection layers')\nax.set_ylabel('X (mm)')\nax.set_zlabel('Y (mm) ')\nplt.title('Trajectories of top weights particles in Cartesian coordinates.')\n\nplt.show()\n\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"957bd75e516aca00c6fd75568e8713034cc07317"},"cell_type":"markdown","source":"## Particle Trajectories in Spherical Coordinates $\\{R,\\phi,\\theta\\}$\n\nWe can already see that a Spherical Transformation is able to isolate out the trajectories adn this will be useful if we are to use clustering methods. But there are still erratic paths in this coordinate system"},{"metadata":{"_uuid":"46182436ed85526f17f35e37eaafba5cbbd4c6f5","trusted":true},"cell_type":"code","source":"# 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":null,"outputs":[]},{"metadata":{"_uuid":"73efec1325a7cdd618bd64fce06a132ae1adfa75"},"cell_type":"markdown","source":"## Particle Trajectories in Angular Track Coordinates $\\{x_2.y_2,z_2\\}$\n\nIn the below figure the hits have been transformed into the angular track coordinates that was discussed in the previous sections. We can see that the hits from the same path are quite tightly grouped and much better than the spherical coordinates. This is the major transformation that will be used in the sample submission that will be developed."},{"metadata":{"_uuid":"c8a83987591069c06511a066ca90d8218997f0b7","trusted":true},"cell_type":"code","source":"# Plot with angular track coordinates.\nfig3 = plt.figure(figsize=(10, 10))\nax = fig3.add_subplot(111, projection='3d')\nfor ang_df in ang_list:\n ax.plot(\n xs=ang_df.z2,\n ys=ang_df.x2,\n zs=ang_df.y2,\n marker='o')\nax.set_xlabel('z2')\nax.set_ylabel('x2')\nax.set_zlabel('y2')\nplt.title('Trajectories of top weights particles in anglular track coordinates.')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"21a8a7b184f99bb5412c32dfc53523ee18c5a587"},"cell_type":"markdown","source":"## Other Feature Engineering Options\nThere are plenty of other feature engineering methods and approaches that are there and I didn't try out. Mostly because I wanted to understand what I was doind and not just take stuff from kernels and run for better results. Hence even though this approach only yeilds $\\approx$ __0.2__ on the public leaderboard, I will leave the more advanced methods to when I have time to comprehend the Phsycis and Maths behind them. A non-exhaustive list of the other options are :\n- Hough Transformations\n- Kalman Filters \n- Scaled Angular Track coordinates\n\nNow lets move onto the Algorithm and the implementation"},{"metadata":{"_uuid":"3b4497fe09608c4f3dbdc2bfd92be555ccd86c07"},"cell_type":"markdown","source":"## The Chosen Approach\n\nAs mentioned before this challenge doens't rely on advanced ML algorithms. Given good and reasonable features the approach is typically to use some variant of Clustering to group hits into particles and predict the tracks. The most common kNN clustering is also valid but doens't yield the best results. Typically __DBSCAN__ is the preferred approach for this problem.\n\n### DBSCAN \n**DBSCAN** (Density-Based Spatial Clustering and Application with Noise), is a density-based clusering algorithm, introduced in Ester et al. 1996, which can be used to identify clusters of any shape in a data set containing noise and outliers.\n\nTwo important parameters are required for DBSCAN: \n\n* epsilon (**eps**) : defines the radius of neighborhood around a point x. It’s called called the ϵ-neighborhood of x. \n* minimum points (**MinPts**) : the minimum number of neighbors within “eps” radius.\n\nBy default the public DBSCAN benchmark used __eps=0.008__ and achieves an average score of __0.207__. A major task that can be done using ML is using the existing data to make a good guess on the value of __eps__. This improves the clustering and the overall score.\n\n\nBelow is an implementation of a clustering __Class__ that will preprocess the data according to the Angular Track transformation and then proceed to cluster using __DBSCAN__"},{"metadata":{"_uuid":"561b19cdb791333ffea4deaa1cb5e9c98594ebca","trusted":true},"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler\nfrom sklearn.cluster import DBSCAN\n\nclass Clusterer(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","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e818451d2e4d365e5d44f89c4723103dca6ab133"},"cell_type":"markdown","source":"### Sample clustering of hits of one event "},{"metadata":{"trusted":true,"_uuid":"7b310f2011015be550428fa7c4decc0cb07b56bf"},"cell_type":"code","source":"model = Clusterer(eps=0.00715)\nlabels= model.predict(hits)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"77978290c4c5a472c23c8a98e45828b2b8ecfe5b"},"cell_type":"code","source":"# The cluster labels generated \nprint(labels)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d6fd06096a5e542293703edade397dd3fc7e3481"},"cell_type":"markdown","source":"### Function to create a submission "},{"metadata":{"trusted":true,"_uuid":"7acefb5cb37372be4dbd51b3cb61e012e157081b"},"cell_type":"code","source":"def 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","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"10c1c0de01fdd74ca49b1af1043e007e3dcc8b8b"},"cell_type":"markdown","source":"### Checking the score of the sample clustering above"},{"metadata":{"trusted":true,"_uuid":"8182d4975f2a1080f8b760e82089f28ad0daa595"},"cell_type":"code","source":"submission = create_one_event_submission(0, hits, labels)\nscore = score_event(truth, submission)\nprint(\"The score is : \", score)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"033c304704c1e05ad2f8d4a3a8650164db6b6570"},"cell_type":"markdown","source":"### Cross Validating \nChecking the mean score for 5 events that are occur 1000 events after the sample 1st event.\nWe see that we achieve an average score of $\\approx$ __0.2__"},{"metadata":{"trusted":true,"_uuid":"8c43510b9a1fae952374456736377325b459799a"},"cell_type":"code","source":"dataset_submissions = []\ndataset_scores = []\n\nfor event_id, hits, cells, particles, truth in load_dataset(path_to_train, skip=1000, nevents=5):\n \n # Track pattern recognition\n model = Clusterer(eps=0.00715)\n labels = model.predict(hits)\n \n # Prepare submission for an event\n one_submission = create_one_event_submission(event_id, hits, labels)\n dataset_submissions.append(one_submission)\n \n # Score for the event\n score = score_event(truth, one_submission)\n dataset_scores.append(score)\n \n print(\"Score for event %d: %.3f\" % (event_id, score))\n \nprint('Mean score: %.3f' % (np.mean(dataset_scores)))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4d632d09f9a955a7a70e38084e5b826b31744683"},"cell_type":"markdown","source":"## Create a Submission File\nUsing the test data this time we create a submission file"},{"metadata":{"trusted":true,"_uuid":"62f84223dc621c387db7e1ef2b2373b95ae7e5f9","_kg_hide-output":true},"cell_type":"code","source":"test_dataset_submissions = []\n\ncreate_submission = False # True for submission \n\nif create_submission:\n for event_id, hits, cells in load_dataset(path_to_test, parts=['hits', 'cells']):\n\n # Track pattern recognition\n model = Clusterer(eps=0.00715)\n labels = model.predict(hits)\n\n # Prepare submission for an event\n one_submission = create_one_event_submission(event_id, hits, labels)\n test_dataset_submissions.append(one_submission)\n \n print('Event ID: ', event_id)\n\n # Create submission file\n submussion = pd.concat(test_dataset_submissions, axis=0)\n submussion.to_csv('submission.csv.gz', index=False, compression='gzip')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"af1f51825c42501b6a2ed6e933b9af829e8f9d06"},"cell_type":"markdown","source":"# Results\nAs expected the submission scored __0.20643__ on the public leaderboard and __0.20817__ on the private leaderboard. This is good as we haven't overfir to the public data and still perform equally well on the whole data set. \n\nUnfortunately this is the extent of my Kaggle competition and exploration. There are still many ML methods that I haven't ventured into, because the nature of the competiton makes it either redundant or very computationlly intensive. The are and not limited to:\n- Trying XGBoost on multiple coordinate transformations as features along with particle data\n- Using exploratory kNN and other approaches to optimize the __eps__ value for DBSCAN\n- Utilizing track voting and ensemble-esqe methods to improve final track prediction"},{"metadata":{"_uuid":"06f64546349881c35d735c85c2cb721737d701ef"},"cell_type":"markdown","source":"# Conclusions\n\nThe Kaggle Course was immensely helpful in providing tips and tricks to gain the best model given the limitations/constraints of a Kaggle Competition. The weekly topics were from a wide variety of domains and topics and with expert talks helping me understand the best practices, whether it was cross validation, model selection or ensemble methods.\n\nAlthough in my specific Kaggle Challenge I wasn't able to implement many of the topics that were discussed. I really loved the __TrackML__ competition, it gave a perspective of domain specific problems that cannot be tackled by simple application of well known ML algorithms. It helped me better understand the value of a good data exploration and visualization using seaborn. These are the most valuable learning that I take away from the course and the competition in general.\n\n__Thank You for reading this long and technical notebook__\n\n
\n\n
"},{"metadata":{"trusted":true,"_uuid":"beabe61f6334275032821e03d78fd97d0c202843"},"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.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}