{"cells":[{"metadata":{"_uuid":"d99098fcbfa6a9ec3ca72e4f681233c323c255c0"},"cell_type":"markdown","source":"This python notebook uses some ideas of the paper **Fitting helices to data by total least squares** by *Yves Nievergelt*.\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":2,"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":3,"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":{"collapsed":true,"trusted":true,"_uuid":"16b08cbe2fac9158bab17e91fe2942db21427ea3"},"cell_type":"code","source":"xm = np.mean(x)","execution_count":4,"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).T","execution_count":5,"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":6,"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":7,"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}^3,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":8,"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":9,"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":10,"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":11,"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":12,"outputs":[]},{"metadata":{"_cell_guid":"595f01a7-fa03-4398-abb5-2354ca359fa6","_uuid":"0ace6a8761680565b177f0a1b12f85949fecb599","collapsed":true,"trusted":true},"cell_type":"code","source":"hits, cells, particles, truth = load_event(os.path.join(path_to_train, event_prefix))","execution_count":13,"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\n\nclass Clusterer(object):\n    def __init__(self,rz_scales=[0.65, 0.965, 1.428]):                        \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.int32)\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            x = M[index]\n            indices[i] = len(index)\n            local_mask = np.ones((len(index)),np.bool)\n            norms[i] = self._test_quadric(x)\n        threshold = np.percentile(norms,90)*4\n        for i, cluster in tqdm(enumerate(labels),total=len(labels)):\n            if norms[i] > threshold:\n                self.clusters[self.clusters==cluster]=0            \n    def _test_quadric(self,x):\n        if len(x.shape)==3:\n            x = np.reshape(x,(x.shape[0],3))\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        \n        return X\n    \n    def predict(self, hits):        \n        volumes = np.unique(hits['volume_id'].values)\n        X = self._preprocess(hits)\n        self.clusters = np.zeros((len(X),1),np.int32)\n        max_len = 1\n        cl = hdbscan.HDBSCAN(min_samples=1,min_cluster_size=7,\n                             metric='braycurtis',cluster_selection_method='leaf',algorithm='best', leaf_size=50)\n        self.clusters = cl.fit_predict(X)+1\n        labels = np.unique(self.clusters)\n        n_labels = 0\n        while n_labels < len(labels):\n            n_labels = len(labels)\n            self._eliminate_outliers(labels,X)\n            max_len = np.max(self.clusters)\n            self.clusters[self.clusters==0] = cl.fit_predict(X[self.clusters==0])+max_len\n            labels = np.unique(self.clusters)\n        return self.clusters","execution_count":14,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"c0bcf488d3b05ba63ad0b15b13db62b445ddbe3b"},"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":16,"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 = True # 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":17,"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}