{"cells":[{"metadata":{"_uuid":"d99098fcbfa6a9ec3ca72e4f681233c323c255c0"},"cell_type":"markdown","source":"This python notebook inspired in **Grzegorz Sionkowsk** algorithm to initialize the clusters with DBSCAN.\nAfter the initialization there is a post processing of clusters using  some ideas of the paper **Fitting helices to data by total least squares** by **Yves Nievergelt** and thresholds of minimum and maximum clusters sizes.\n\nThe method to determine if a set of points in space fits in a helix is shown below"},{"metadata":{"collapsed":true,"trusted":true,"_uuid":"2ba21c9c61d5bc3dc4fd588997b6b549fa62aec7"},"cell_type":"code","source":"import numpy as np","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"dd406b0b10810500e10921a58c59c2887509b2db"},"cell_type":"markdown","source":"### Data points\nThe data considered here consists of p points $$\\vec{x_1}, ...,\\vec{x_p}$$"},{"metadata":{"collapsed":true,"trusted":true,"_uuid":"cebc02749a7e94349d73c20e77bb3fccf3633f5e"},"cell_type":"code","source":"x = np.array([(62,397,103),(82,347,107),(93,288,120),\n     (94,266,128),(65,163,169),(12,102,198),\n     (48,138,180),(77,187,157),(85,209,149),(89,316,112)])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"699d84b5c58a7e6ff4a9d0abeeb29420bbe08bc5"},"cell_type":"markdown","source":"### 1. Estimating the degree of colinearity of data"},{"metadata":{"_uuid":"68a07c3607e134e06dd962e526c0a351b66cee94"},"cell_type":"markdown","source":"#### 1.1 Compute the average\n$$\\bar{\\vec{x}}:= (1/p)\\sum_{j=1}^{p}\\vec{x_j}$$"},{"metadata":{"trusted":true,"_uuid":"16b08cbe2fac9158bab17e91fe2942db21427ea3","collapsed":true},"cell_type":"code","source":"xm = np.mean(x,axis=0)\nprint(xm.shape)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e9c377195c450dd14dbc2a790e6da878d3ae8fee"},"cell_type":"markdown","source":"#### 1.2 Form the matrix\n$$X\\in\\mathbb{M}_{p\\times{3}}$$"},{"metadata":{"collapsed":true,"trusted":true,"_uuid":"12bf91cfac1dba26e4be2c692f527e0286b3652a"},"cell_type":"code","source":"x = x - xm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"173b568bf1fbe0864f32d7c213a8f35e66c68e44","collapsed":true},"cell_type":"code","source":"print(x.shape)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"bc1f099250a7129bc42d1b6e503ce0745be2d736"},"cell_type":"markdown","source":"#### 1.3 Compute the singular values of X\n$$\\sigma_1\\geqslant\\sigma_2\\geqslant\\sigma_3\\geqslant0$$\nand the corresponding orthonormal vectors\n$$\\vec{v_1},\\vec{v_2},\\vec{v_3}\\in\\mathbb{R}^3$$"},{"metadata":{"collapsed":true,"trusted":true,"_uuid":"82bbe8c18b4cb863953514e1993e0487ae966382"},"cell_type":"code","source":"v, s, t = np.linalg.svd(x,full_matrices=True)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"trusted":true,"_uuid":"fb2aedd240dca5f1c00e32371f22756d9a60b0f0"},"cell_type":"code","source":"sigma1 = s[0]\nsigma2 = s[1]\nsigma3 = s[2]\nv1 = t[0]\nv2 = t[1]\nv3 = t[2]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9813a5c5bb41ca333b438741328901b3eeadcf7d"},"cell_type":"markdown","source":"##### If $$\\sigma_2>\\sigma_3\\geqslant0$$ \nthe plane of total least squares satisfies the equation $$\\langle\\vec{x}-\\bar{\\vec{x}},\\vec{v_3}\\rangle=0$$\nin particular if $$\\sigma_3=0$$\nall data lie in that plane "},{"metadata":{"_uuid":"394189f6bbcf968e92615db002fb8f7c269b835e"},"cell_type":"markdown","source":"#####  if\n$$\\sigma_1>\\sigma_2=0=\\sigma_3$$\nthen all data lie in a straight line"},{"metadata":{"_uuid":"058ae54af4515c554454fef13671fc014573f55f"},"cell_type":"markdown","source":"##### Similarly, if\n$$\\sigma_1=\\sigma_2=\\sigma_3=0$$\nall the data coalesce at one point"},{"metadata":{"_uuid":"badc42fc468b99d8c1d5bdae06c333793d8866dc"},"cell_type":"markdown","source":"### 2. Fitting the axis and the radius of the helix"},{"metadata":{"_uuid":"60ebbd37f502a2f42eeee09a5c6a7b26034ba087"},"cell_type":"markdown","source":"#### 2.1 Fitting a quadric surface to the data"},{"metadata":{"_uuid":"0c708bf0dc71a3a754eb78f9566485a330f4f15e"},"cell_type":"markdown","source":"Each affine quadric surface satisfies the equation $$F(S;\\vec{x})=0$$ defined by a quadratic\nform $$F(S;\\vec{x})=(\\bar{x}^T,1).S.(\\bar{x}^T,1)^T=\\begin{bmatrix}x_1 & x_2 & x_3 & 1\\end{bmatrix}\\begin{bmatrix}\n    s_{11} & s_{12} & s_{13} & s_{1} \\\\\n    s_{12} & s_{22} & s_{23} & s_{2} \\\\\n    s_{13} & s_{23} & s_{33} & s_{3} \\\\\n    s_{1} & s_{2} & s_{3} & s_{44}\n\\end{bmatrix}\\begin{bmatrix} x_1 \\\\ x_2 \\\\ x_3 \\\\ 1\\end{bmatrix}$$"},{"metadata":{"_uuid":"5b5b4426741b2db135aedb1dca3a3afcdaf36eb9"},"cell_type":"markdown","source":"This algorithm will determine the matrix or matrices S minimizing the total least-squares objective\n$$G(S):=\\sum_{j=1}^{p}[F(S;\\vec{x_j})]^2$$\nsubject to the constraint that\n$$\\sum_{k=1}^{4}\\sum_{\\ell=k}^4|{S_{k,\\ell}}|^2=1$$\nComputationally, arrange matrix S in one-dimensional vector in lexicographic order:\n$$S=\\vec{S}:=(s_{11}, s_{12}, s_{13}, s_{14}; s_{22}, s_{23}, s_{24}; s_{33}, s_{34}; s_{44})$$\nSimilarly, for each point x, form the vector z:\n$$\\vec{z}:=(x_{1}^2,2x_{1}x_{2},2x_{1}x_{3},2x_{1};x_{2}^2,2x_{2}x_{3},2x_{2};x_{3}^2,2x_{3}; 1)$$"},{"metadata":{"_uuid":"3f87e4f0bfcd0901909b0230b38fda90c4c35a47"},"cell_type":"markdown","source":"#### 2.1.1 Form the matrix Z"},{"metadata":{"collapsed":true,"trusted":true,"_uuid":"eef568c8d3cf544c4932802da6d20ad4116001e5"},"cell_type":"code","source":"Z = np.zeros((x.shape[0],10), np.float32)\nZ[:,0] = x[:,0]**2\nZ[:,1] = 2*x[:,0]*x[:,1]\nZ[:,2] = 2*x[:,0]*x[:,2]\nZ[:,3] = 2*x[:,0]\nZ[:,4] = x[:,1]**2\nZ[:,5] = 2*x[:,1]*x[:,2]\nZ[:,6] = 2*x[:,1]\nZ[:,7] = x[:,2]**2\nZ[:,8] = 2*x[:,2]\nZ[:,9] = 1","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"465bb028dba6098820d2c2447670aed21ad64902"},"cell_type":"markdown","source":"#### 2.1.2 Compute the smallest singular value and the corresponding right-singular vectors of the matrix Z"},{"metadata":{"trusted":true,"_uuid":"a6b1ba314c5cb892e4189f7d4201febb0f1c95f5","collapsed":true},"cell_type":"code","source":"v, s, t = np.linalg.svd(Z,full_matrices=True)\nsmallest_value = np.min(np.array(s))\nsmallest_index = np.argmin(np.array(s))\nT = np.array(t)\nT = T[smallest_index,:]\nS = np.zeros((4,4),np.float32)\nS[0,0] = T[0]\nS[0,1] = S[1,0] = T[1]\nS[0,2] = S[2,0] = T[2]\nS[0,3] = S[3,0] = T[3]\nS[1,1] = T[4]\nS[1,2] = S[2,1] = T[5]\nS[1,3] = S[3,1] = T[6]\nS[2,2] = T[7]\nS[2,3] = S[3,2] = T[8]\nS[3,3] = T[9]\nnorm = np.linalg.norm(np.dot(Z,T), ord=2)**2\nprint(norm)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"705d1703b65918fdb1ab0ba3f5e3c3f1ffb1435a"},"cell_type":"markdown","source":"##### The norm value near zero shows that x points above are fitted in a quadric surface of a cylinder or helix."},{"metadata":{"_cell_guid":"e081740e-8169-4481-b1df-f5dd5488314f","_uuid":"0bee86255243664f24e4bcf48af2228a3100a8b7","trusted":true,"collapsed":true},"cell_type":"code","source":"%matplotlib inline\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport os\n\nfrom trackml.dataset import load_event, load_dataset\nfrom trackml.score import score_event","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"572fcbb6-8c7b-4a09-8916-8ec76689130f","_uuid":"63414de98667e95f60407c9155899a25a321cffc","collapsed":true,"trusted":true},"cell_type":"code","source":"# Change this according to your directory preferred setting\npath_to_train = \"../input/train_1\"","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"d0f5916b-8270-4ff5-af17-6fbbb8d00553","_uuid":"3e45554ab05c14faf63a2c423f69ebbe7c108541","collapsed":true,"trusted":true},"cell_type":"code","source":"# This event is in Train_1\nevent_prefix = \"event000001000\"","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"595f01a7-fa03-4398-abb5-2354ca359fa6","_uuid":"0ace6a8761680565b177f0a1b12f85949fecb599","trusted":true,"collapsed":true},"cell_type":"code","source":"hits, cells, particles, truth = load_event(os.path.join(path_to_train, event_prefix))","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"e06d1ed7-5091-4d67-abb4-5984b137e2e6","_uuid":"c2f70ae63abffcc09a534bb17fb89df8ffddb722","scrolled":true,"trusted":true,"collapsed":true},"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler\nimport hdbscan\nfrom scipy import stats\nfrom tqdm import tqdm\nfrom sklearn.cluster import DBSCAN\n\nclass Clusterer(object):\n    def __init__(self,rz_scales=[0.65, 0.965, 1.528]):                        \n        self.rz_scales=rz_scales\n    \n    def _eliminate_outliers(self,labels,M):\n        norms=np.zeros((len(labels)),np.float32)\n        indices=np.zeros((len(labels)),np.float32)\n        for i, cluster in tqdm(enumerate(labels),total=len(labels)):\n            if cluster == 0:\n                continue\n            index = np.argwhere(self.clusters==cluster)\n            index = np.reshape(index,(index.shape[0]))\n            indices[i] = len(index)\n            x = M[index]\n            norms[i] = self._test_quadric(x)\n        threshold1 = np.percentile(norms,90)*5\n        threshold2 = 25\n        threshold3 = 6\n        for i, cluster in enumerate(labels):\n            if norms[i] > threshold1 or indices[i] > threshold2 or indices[i] < threshold3:\n                self.clusters[self.clusters==cluster]=0   \n    def _test_quadric(self,x):\n        if x.size == 0 or len(x.shape)<2:\n            return 0\n        xm = np.mean(x,axis=0)\n        x = x - xm\n        Z = np.zeros((x.shape[0],10), np.float32)\n        Z[:,0] = x[:,0]**2\n        Z[:,1] = 2*x[:,0]*x[:,1]\n        Z[:,2] = 2*x[:,0]*x[:,2]\n        Z[:,3] = 2*x[:,0]\n        Z[:,4] = x[:,1]**2\n        Z[:,5] = 2*x[:,1]*x[:,2]\n        Z[:,6] = 2*x[:,1]\n        Z[:,7] = x[:,2]**2\n        Z[:,8] = 2*x[:,2]\n        Z[:,9] = 1\n        v, s, t = np.linalg.svd(Z,full_matrices=False)        \n        smallest_index = np.argmin(np.array(s))\n        T = np.array(t)\n        T = T[smallest_index,:]        \n        norm = np.linalg.norm(np.dot(Z,T), ord=2)**2\n        return norm\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        for i, rz_scale in enumerate(self.rz_scales):\n            X[:,i] = X[:,i] * rz_scale\n       \n        return X\n    def _init(self, dfh, w1, w2, w3, w4, w5, w6, w7, epsilon, Niter):\n        dfh['r'] = np.sqrt(dfh['x'].values ** 2 + dfh['y'].values ** 2 + dfh['z'].values ** 2)\n        dfh['rt'] = np.sqrt(dfh['x'].values ** 2 + dfh['y'].values ** 2)\n        dfh['a0'] = np.arctan2(dfh['y'].values, dfh['x'].values)\n        dfh['z1'] = dfh['z'].values / dfh['rt'].values\n        dfh['z2'] = dfh['z'].values / dfh['r'].values\n        dfh['s1'] = dfh['hit_id']\n        dfh['N1'] = 1\n        dfh['z1'] = dfh['z'].values / dfh['rt'].values\n        dfh['z2'] = dfh['z'].values / dfh['r'].values\n        dfh['x1'] = dfh['x'].values / dfh['y'].values\n        dfh['x2'] = dfh['x'].values / dfh['r'].values\n        dfh['x3'] = dfh['y'].values / dfh['r'].values\n        dfh['x4'] = dfh['rt'].values / dfh['r'].values\n        mm = 1\n        for ii in tqdm(range(int(Niter))):\n            mm = mm * (-1)\n            dfh['a1'] = dfh['a0'].values + mm * (dfh['rt'].values + 0.000005\n                                                 * dfh['rt'].values ** 2) / 1000 * (ii / 2) / 180 * np.pi\n            dfh['sina1'] = np.sin(dfh['a1'].values)\n            dfh['cosa1'] = np.cos(dfh['a1'].values)\n            ss = StandardScaler()\n            dfs = ss.fit_transform(dfh[['sina1', 'cosa1', 'z1', 'z2','x1','x2','x3','x4']].values)\n            cx = np.array([w1, w1, w2, w3, w4, w5, w6, w7])\n            dfs = np.multiply(dfs, cx)\n            clusters = DBSCAN(eps=epsilon, min_samples=1, metric=\"euclidean\", n_jobs=32).fit(dfs).labels_\n            dfh['s2'] = clusters\n            dfh['N2'] = dfh.groupby('s2')['s2'].transform('count')\n            maxs1 = dfh['s1'].max()\n            dfh.loc[(dfh['N2'] > dfh['N1']) & (dfh['N2'] < 20),'s1'] = dfh['s2'] + maxs1\n            dfh['N1'] = dfh.groupby('s1')['s1'].transform('count')\n        return dfh['s1'].values\n    def predict(self, hits):         \n        self.clusters = self._init(hits,2.7474448671796874,1.3649721713529086,0.7034918842926337,\n                                        0.0005549122352940002,0.023096034747190672,0.04619756315527515,\n                                        0.2437077420144654,0.009750302717746615,338)\n        X = self._preprocess(hits) \n        cl = hdbscan.HDBSCAN(min_samples=1,min_cluster_size=7,\n                             metric='braycurtis',cluster_selection_method='leaf',algorithm='best', leaf_size=50)\n        labels = np.unique(self.clusters)\n        self._eliminate_outliers(labels,X)          \n        max_len = np.max(self.clusters)\n        mask = self.clusters == 0\n        self.clusters[mask] = cl.fit_predict(X[mask])+max_len\n        return self.clusters","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"70612062632493a78bec5bd5c69c5d4d523b8b83","collapsed":true},"cell_type":"code","source":"model = Clusterer()\nlabels = model.predict(hits)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"3d87954eeb3413d98fab9732172c5ef56624602c"},"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\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c0bcf488d3b05ba63ad0b15b13db62b445ddbe3b","collapsed":true},"cell_type":"code","source":"submission = create_one_event_submission(0, hits, labels)\nscore = score_event(truth, submission)\nprint(\"Your score: \", score)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"2b20cf0c-2754-48dd-ab1e-4f489c2aa05c","_uuid":"7f8de52b9022581bf10aa813d2db005b842f0be7","scrolled":true,"trusted":true,"collapsed":true},"cell_type":"code","source":"path_to_test = \"../input/test\"\ntest_dataset_submissions = []\n\ncreate_submission = False # True for submission \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()\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    submission = pd.concat(test_dataset_submissions, axis=0)\n    submission.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"0a333cd4-351e-4274-aa7c-4cf8ab7fca1a","_uuid":"70ce31d93086e022159d6227f35c6488bf80eb22","collapsed":true,"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"anaconda-cloud":{},"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}