{"cells":[{"metadata":{"_uuid":"dfeb2ea91f47c09a39df8aa0b45df062f7b74dd9"},"cell_type":"markdown","source":"____\n\nIn this kernel I've made a small simulation framework to test your features (based on https://www.kaggle.com/c/trackml-particle-identification/discussion/62804#latest-368156) and provide my own attempts to find best features.\n\n___\n\nI have been struggling to solve this clustering competition. I understand that we need to find the proper features to get beyond 0.6, but without public kernels I can't get more that 0.2. I believe that I'm missing something very obivous.\n\nI've read almost all forum topics, I think, these are the most valuable posts (may be this list will help someone):\n\nhttps://www.kaggle.com/c/trackml-particle-identification/discussion/62804#latest-368156\nhttps://www.kaggle.com/khahuras/0-53x-clustering-using-hough-features-basic/notebook\nhttps://www.kaggle.com/c/trackml-particle-identification/discussion/61590#latest-368105\nhttps://www.kaggle.com/c/trackml-particle-identification/discussion/61081#latest-367422\n\n_____\n\nI'm afraid I won't have time to make and calculate some top-score solution. So let's try to win at least a kernel medal :)"},{"metadata":{"_uuid":"e906a332eb2a2e3527bcb2057aaf1be9041cd84f"},"cell_type":"markdown","source":"Let's assume that helix has this formula:\n\n$$x = -R \\sin{(\\omega \\cdot t + \\theta_0)} + (R - D)\\sin{(\\theta_0)}$$\n$$y = -R \\cos{(\\omega \\cdot t + \\theta_0)} + (R - D)\\cos{(\\theta_0)}$$\n$$z = z_0 + v \\cdot t $$\n\nAs we can see, helix has 4 variables: $x,y,z,t$ and 6 parameters: $\\omega, v, R, \\theta_0, D, z_0$.\n\nLet's fix 3 of them: $R, \\theta_0, z_0 $. We will iterate over their values from a certain range to reduce the number of variables in the equation.\n\nPseudocode:\n``` python\nfinal_clusters = \nfor theta0 in np.linspace(theta0_from, theta0_to, theta0_iterations):\n    for z_0 in np.linspace(z_0_from, z_0_to, z_0_iterations):\n        for R in np.linspace(R_from, R_to, R_iterations):\n            features = calculate_features(hits_data, theta0, z0, R)\n            clusters = clisterize(features)\n            print(score(clusters))\n            final_clusters = merge(clusters=clusters, to=final_clusters)\n            \nprint(\"Final score: {}\".format(score(clusters)))\n```\n\nNow we have three equations with only three unknowns: $w, v, D$. I'll show you that we can combine $w$ and $v$ to $\\frac{w}{v}$, so we will have only 2 unknowns.\n\nThese variables are constants for each helix. Let's find their equations:\n\n$$D = f(x,y,z,t,R, \\theta_0, z_0)$$\n$$\\frac{w}{v} = g(x,y,z,t,R, \\theta_0, z_0)$$\n\n_____"},{"metadata":{"_uuid":"370f4404a41d64cbd50192b596bf087775502324"},"cell_type":"markdown","source":"# The equation for $D$\n\n$$R \\sin{(\\omega \\cdot t + \\theta_0)} =  (R - D)\\sin{(\\theta_0)} - x$$\n$$R \\cos{(\\omega \\cdot t + \\theta_0)} = (R - D)\\cos{(\\theta_0)} - y$$\n\n\n$$R^2 \\sin^2{(\\omega \\cdot t + \\theta_0)} =  (R - D)^2(\\sin{(\\theta_0)} - x)^2 $$ \n$$R^2 \\cos^2{(\\omega \\cdot t + \\theta_0)} = (R - D)^2(\\cos{(\\theta_0)} - y)^2 $$\n\nSum:\n\n$$R^2 = (R-D)^2 + x^2 + y^2 - 2(R-D)\\cdot(x\\sin\\theta_0 + y\\cos\\theta_0)$$\n\nLet's denote the expression\n\n$$\\gamma = x\\sin\\theta_0  +  y\\cos\\theta_0$$\n\nThen\n\n$$D^2 - 2RD + x^2 + y^2 - 2R\\gamma + 2D\\gamma = 0$$\n\n$$D^2 - 2D(R-\\gamma) + x^2 +y^2 - 2R\\gamma = 0$$\n\nSolve the quadratic equation:\n\n$$D = R - \\gamma ± \\sqrt{R^2 + \\gamma^2 - x^2 -y^2} $$\n\nOr using some trigonometric identities\n\n$$D = R - (x\\sin\\theta_0  +  y\\cos\\theta_0) ± \\sqrt{R^2 - (x\\cos\\theta_0  +  y\\sin\\theta_0)^2} $$.\n\nNotice that there are two possible solutions.\n_____"},{"metadata":{"_uuid":"1f29712e729446b42f5dc654854a131eea2a45dc"},"cell_type":"markdown","source":"# The equation for $\\frac{\\omega}{v}$\nLet's find the equation for \\frac{\\omega}{v}  too.\n\nLet's get rid of $D$ using a simple proportion:\n$$\\frac{x + R\\sin(\\omega t + \\theta_0)}{\\sin(\\theta_0)} = R - D = \\frac{y + R\\cos(\\omega t + \\theta_0)}{\\cos(\\theta_0)}$$\n\n$$x\\cos(\\theta_0) - y\\sin(\\theta_0) = R[\\sin(\\theta_0)\\cdot \\cos(\\omega t + \\theta_0) - \ncos(\\theta_0) \\cdot \\sin(\\omega t + \\theta_0)] = R \\sin(-\\omega t)$$\n\nExpress the time from the equation for the z coordinate\n\n$$ t = \\frac{z-z0}{v}$$\n\nAnd finally\n\n\n$$\\frac{\\omega}{v} = \\frac{1}{z-z_0} \\cdot \\arcsin \\frac{y\\sin( \\theta_0) - x\\cos( \\theta_0)}{R}$$\n\n______"},{"metadata":{"_uuid":"5d5f90566c7b92768755b960ee7a8b595c3304c8"},"cell_type":"markdown","source":"#### Let's test this features on the perfect helixes"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"collapsed":true},"cell_type":"code","source":"import os\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"a80a8d8c83306308c4ee75c11f2d8bdddeb37391"},"cell_type":"code","source":"def get_helix(R=2000, w=np.pi * 3, v=3000, theta0=np.pi / 3, D=500, z0=100):\n    # Generates perfect helix data for feature testing\n    t = np.random.rand(1000)\n    df = pd.DataFrame(index=t)\n    df[\"z\"] = z0 + v * t\n    df[\"x\"] = -R * np.sin(w * t + theta0)+(R - D) * np.sin(theta0)\n    df[\"y\"] = -R * np.cos(w * t + theta0)+(R - D) * np.cos(theta0)\n    return df\n\ndef plot_helix(df):\n    # 3D plot\n    from mpl_toolkits.mplot3d import Axes3D\n    import matplotlib.pyplot as plt\n    fig = plt.figure()\n    ax = fig.add_subplot(111, projection='3d')\n    ax.scatter(df.x, df.y, df.z)\n    ax.set_xlabel('X Label')\n    ax.set_ylabel('Y Label')\n    ax.set_zlabel('Z Label')\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ce4459493e6841eec545214d0876827d02f20e7e"},"cell_type":"markdown","source":"Let's visualize our generated helixes."},{"metadata":{"trusted":true,"_uuid":"0900122380f30875512df2acb7cfc1dd0d3e4cc3"},"cell_type":"code","source":"df = get_helix()\ndf.plot.scatter(x=\"x\", y=\"y\")\nplot_helix(df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e8048afe2f03e60992e3ead9fc7825d5f877bfcf","collapsed":true},"cell_type":"code","source":"def add_features(df, theta0=0, R=2000, z0=0):\n    x,y,z = df.x, df.y, df.z\n    \n    x_sin_y_cos = x * np.sin(theta0) + y * np.cos(theta0)\n    x_cos_y_sin = x * np.cos(theta0) + y * np.sin(theta0)\n    df[\"d1\"] = (R - x_sin_y_cos + np.sqrt(R ** 2 - (x_cos_y_sin) ** 2))\n    df[\"d2\"] = (R - x_sin_y_cos - np.sqrt(R ** 2 - (x_cos_y_sin) ** 2))\n    df[\"w_div_v\"] = np.arcsin((y * np.sin(theta0) - x * np.cos(theta0)) / R) / (z - z0)\n    return df","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a33ed38e2f91f12f52741853bdb4a6000562c15a"},"cell_type":"markdown","source":"Let's generate several helixes with different parameters to check our features."},{"metadata":{"trusted":true,"_uuid":"21e124c141be94712eb597bbeafdce8511e6dfbf","collapsed":true},"cell_type":"code","source":"h1 = get_helix(theta0=0, R=2000, z0=0)\nh2 = get_helix(theta0=-0.6, R=1000, z0=250)\nh3 = get_helix(theta0=0.6, R=1500, z0=50)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6b93c93c50a2de4ce9bd5d852dd9e0e8ac9bad7b"},"cell_type":"markdown","source":"Let's assume that we find the params for one helix and calculate the features based on these parameters."},{"metadata":{"trusted":true,"_uuid":"c77130723af0298dcb466c2c0f0a8242ce6939ec"},"cell_type":"code","source":"theta0 = -0.6\nR = 1000\nz0 = 250\n\nh1 = add_features(h1, theta0=theta0, R=R, z0=z0)\nh2 = add_features(h2, theta0=theta0, R=R, z0=z0)\nh3 = add_features(h3, theta0=theta0, R=R, z0=z0)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ebcd1c4904bc4cefede734158e895c4c6731ec01"},"cell_type":"markdown","source":"Ooops! We saw some errors inside the arcsin. That means the $(x,y,z)$ point is definately not on the helix with params $(\\theta_0, R, z_0)$."},{"metadata":{"_uuid":"5d45c0df8b8aa705201f4476932eedcf29c208f8"},"cell_type":"markdown","source":"Let's plot the features for every helix. Remember that the green helix is the right one."},{"metadata":{"trusted":true,"_uuid":"511896bb8ed7a811b6bb6a3b854d11cbc8c78300"},"cell_type":"code","source":"ax = h1.plot.scatter(x=\"x\", y=\"d1\", c=\"r\")\nax = h2.plot.scatter(x=\"x\", y=\"d1\", c=\"g\", ax=ax)\nax = h3.plot.scatter(x=\"x\", y=\"d1\", c=\"b\", ax=ax, title=\"Feature: D1\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"271a6ad04bd4319e9385e3793fa3c9e3da07f7bc"},"cell_type":"code","source":"ax = h1.plot.scatter(x=\"x\", y=\"d2\", c=\"r\")\nax = h2.plot.scatter(x=\"x\", y=\"d2\", c=\"g\", ax=ax)\nax = h3.plot.scatter(x=\"x\", y=\"d2\", c=\"b\", ax=ax, title=\"Feature: D2\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"70b8d586f494283ce53da29acd6b7a94ec79bd89"},"cell_type":"code","source":"ax = h1.plot.scatter(x=\"x\", y=\"w_div_v\", c=\"r\")\nax = h2.plot.scatter(x=\"x\", y=\"w_div_v\", c=\"g\", ax=ax)\nax = h3.plot.scatter(x=\"x\", y=\"w_div_v\", c=\"b\", ax=ax, title=\"Feature: w/v\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"eae69f346d54e76fdf634b4238bf8c7c804cee77"},"cell_type":"markdown","source":"The D1 and D2 features failed. May be we should try some combination. Does it mean that the $D$ of a helix is not constant?"},{"metadata":{"_uuid":"f13c4ee35a10f824d1fb746de297b86cf9834aa7"},"cell_type":"markdown","source":"Let's wee the `value_counts` of the last feature - it seems to be the best:"},{"metadata":{"trusted":true,"_uuid":"1f7ef364930749169055d2b707faf4820de13388"},"cell_type":"code","source":"h2.w_div_v.value_counts().head(20)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"75c384ac6911acf617e9635ba1b83c1f0607868b"},"cell_type":"markdown","source":"We can see that a lot of values are the same. That means I have chosen at least one right feature. Right, @yuval?"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"e026f93ea0e2b8e39afe402bcd909f6ffea29046"},"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}