{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"collapsed":true},"cell_type":"markdown","source":"This notebook was based on someone else's code which I downloaded and edited, but I can't find the original kernel.  I also took a lot from Grzegorz Sionkowski's Bayesian Optimization kernel, and used a modified version of that kernel to tune the parameters."},{"metadata":{"trusted":true,"_uuid":"fd70413fdd13e0f32f8647a24b1ccf25a098d219","collapsed":true},"cell_type":"code","source":"%matplotlib inline\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport os\nimport math\n\nfrom trackml.dataset import load_event, load_dataset\nfrom trackml.score import score_event\n\nfrom multiprocessing import Pool\n\ndef 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":{"trusted":true,"collapsed":true,"_uuid":"e80042d686334164dccf5e36c5d8c505cfa0d586"},"cell_type":"code","source":"# Change this according to your directory preferred setting\npath_to_train = \"../input/train_1/\"\npath_to_test = \"../input/test/\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"24ef0fc220257022427cb3020d62bd2d494ccb8f"},"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler\nimport hdbscan\nfrom scipy import stats\nfrom tqdm import tqdm_notebook as tqdm\nfrom sklearn.cluster import DBSCAN, AgglomerativeClustering\n\nclass Clusterer(object):\n    def __init__(self,rz_scales=[0.65, 0.965, 1.528], eps=0.0035, dz0 = -0.00070, stepdz = 0.00001, stepeps = 0.000005, num_loops=100):                        \n        self.rz_scales=rz_scales\n        self.epsilon = eps\n        self.dz0 = dz0\n        self.stepdz = stepdz\n        self.stepeps = stepeps\n        self.num_loops = num_loops\n        \n    # remove outliers\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    \n    # not sure what this function does?\n    def _test_quadric(self,x):\n        if x.size == 0 or len(x.shape)<2:\n            return 0\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    # standard scale our data\n    def _preprocess(self, hits):\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    \n    def _init(self,dfh):\n        dfh['s1'] =dfh.hit_id\n        dfh['N1'] =1\n        dfh['r'] = 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']/dfh['r']\n        dfh['z3'] = 1/dfh['z1'].values\n        mm = 1\n\n        for ii in range(self.num_loops):\n            mm = mm*(-1)                               \n            dfh['a1'] = dfh['a0'].values+mm*(dfh['rt']+self.stepeps*dfh['rt']**2)/1000*(ii/2)/180*math.pi\n            dfh['z'].values*np.sign(dfh['z'].values)\n            dfh['x1'] = dfh['a1'].values/dfh['z1'].values\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', 'z3', 'x1']].values)\n            cx = np.array([1.08, 1.08, .38, .25, 0.009, 0.0001])\n            dfs = np.multiply(dfs, cx)\n\n            clusters=DBSCAN(eps=self.epsilon,min_samples=1,metric='euclidean',n_jobs=1).fit(dfs).labels_ \n\n            if ii==0:\n                dfh['s1'] = clusters\n                dfh['N1'] = dfh.groupby('s1')['s1'].transform('count')\n\n            # else update our hits conditionally, if it's a better fit\n            else:\n                # put our new clusters to another feature\n                dfh['s2'] = clusters\n\n                # get the count of those clusters\n                dfh['N2'] = dfh.groupby('s2')['s2'].transform('count')\n                maxs1 = dfh['s1'].max()\n\n                # if our new clusters are bigger, but less than 20, use the new ones instead\n                cond = np.where((dfh['N2'].values>dfh['N1'].values) & (dfh['N2'].values<20))\n                s1 = dfh['s1'].values\n                s1[cond] = dfh['s2'].values[cond]+maxs1\n\n                # write the new clusters back to our dataframe\n                dfh['s1'] = s1\n                dfh['s1'] = dfh['s1'].astype('int64')\n                dfh['N1'] = dfh.groupby('s1')['s1'].transform('count')\n        \n        # return our clusters\n        return dfh['s1'].values    \n    \n    def predict(self, hits):    \n        # init our clusters\n        self.clusters = self._init(hits) \n        \n        # preprocess our data\n        X = self._preprocess(hits) \n        \n        # create our last clusterer\n        cl = hdbscan.HDBSCAN(min_samples=1,min_cluster_size=7, metric='braycurtis',cluster_selection_method='leaf',algorithm='best', leaf_size=50)\n            \n        # labels = unique clusters\n        labels = np.unique(self.clusters)\n        \n        # remove outliers\n        self._eliminate_outliers(labels,X)\n        \n        # init n_labels\n        n_labels = 0\n        \n        # now we loop through the points that haven't been assigned to a cluster and assign them with\n        # HDBSCAN\n        while n_labels < len(labels):\n            n_labels = len(labels)\n            max_len = np.max(self.clusters)\n            mask = self.clusters == 0\n            self.clusters[mask] = cl.fit_predict(X[mask])+max_len\n            \n        return self.clusters","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"ba5892bbd1ad85398d4ce930571ca0caa8e2ef92"},"cell_type":"code","source":"eps = 0.0035\ndz0 = -0.000720\nstepdz = 3e-6\nstepeps = 2e-16\n\ndef one_loop(event_id):\n    hits  = pd.read_csv(path_to_test + '/event%s-hits.csv'%event_id)\n    cells = pd.read_csv(path_to_test + '/event%s-cells.csv'%event_id)\n    print('Event ID: ', event_id)\n                \n    # Track pattern recognition \n    model = Clusterer(eps=eps, dz0=dz0, stepdz=stepdz, stepeps=stepeps)\n    labels = model.predict(hits)\n\n    # Prepare submission for an event\n    one_submission = create_one_event_submission(event_id, hits, labels)\n            \n    return one_submission\n\ndef create_test_submissions(path_to_test = \"../input/test\", start=0, end=15):\n    event_ids = [ '%09d'%i for i in range(start,end) ]\n\n    pool = Pool(processes=4)\n    results = pool.map(one_loop, event_ids)\n    pool.close()\n    \n    return results","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"53cb2342d1c75ba3e39eb8ca0e0f78ffe14dcc0c","collapsed":true},"cell_type":"code","source":"results = create_test_submissions(start=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"b3c14101091b3378299546f7a289a5e6b00e2a10"},"cell_type":"code","source":"submission = pd.concat(results, axis=0)\nsubmission.to_csv('20180707_dbscan_outliers_agg.csv.gz', index=False, compression='gzip')\nprint(len(submission))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"80342906e12271f3013a3fd68814634ba9a1186f"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"7af2a8b8aacea1c3e3d41baf58d4c00f30b4f5df"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"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}