{"cells":[{"metadata":{"_uuid":"d99098fcbfa6a9ec3ca72e4f681233c323c255c0"},"cell_type":"markdown","source":"In this notebook, I would like to use the Hough transform to cluster hits. This notebook is therefore, get some materials from the past published kernels,\n\nIn the previous notebook, we see a function relating phi and r like (where $r = \\sqrt{x^2 + y^2}$, $\\phi = arctan2(y/x)$):\n$$ \\phi_{new} = \\phi + i(ar + br^2),$$\nwhere $i$ is increased incrementally from 0 (straight tracks) to some number (curve tracks).\n\n\nHowever, the above equation is not exact to relate those two features. Instead, one might want to use the Hough transform:\n$$  \\frac{r}{2r_0} =  \\cos(\\phi - \\theta) $$\n\nIn the above equation, $\\phi$ and $r$ are the original $\\phi$ and $r$ of each hit, while $r_0$ and $\\theta$ are the $r$ and $\\phi$ of a specific point in the XY plane, that is the origin of a circle in XY plane. That circle passes through the inspected hit. \n\nThen, our clustering problem can be stated this way:\n- For each $\\frac{1}{2r_0}$, starting from 0 (corresponding to straight tracks), to an appropriate stopping point, we calculate $\\theta = \\phi - \\arccos(\\frac{r}{2r_0})$\n- Group all hits with the near $\\theta$ and some other features to a detected track by DBSCAN. Since $\\theta$ can take very large or small values, using $\\sin(\\theta)$ and $\\cos(\\theta)$ is better.\n"},{"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\nfrom sklearn.preprocessing import StandardScaler\nimport hdbscan\nfrom scipy import stats\nfrom tqdm import tqdm\nfrom sklearn.cluster import DBSCAN\n#from trackml.dataset import load_event, load_dataset\n#from trackml.score import score_event\nfrom IPython.display import clear_output","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a4837535df93ea5af70db4e8525711b31ed1d82a"},"cell_type":"markdown","source":"I don't know why I can't use the trackml library in this notebook, so I copied all the trackml library contents needed in the below block. You can skip reading the below block (I hide it by the way)."},{"metadata":{"_uuid":"062a395857325d7e018ec6045019bf3f5f92e8ed","_kg_hide-input":true,"_kg_hide-output":true,"trusted":true,"collapsed":true},"cell_type":"code","source":"\"\"\"TrackML dataset loading\"\"\"\n\n__authors__ = ['Moritz Kiehn', 'Sabrina Amrouche']\n\nimport glob\nimport os\nimport os.path as op\nimport re\nimport zipfile\n\nimport pandas\n\nCELLS_DTYPES = dict([\n    ('hit_id', 'i4'),\n    ('ch0', 'i4'),\n    ('ch1', 'i4'),\n    ('value', 'f4'),\n])\nHITS_DTYPES = dict([\n    ('hit_id', 'i4'),\n    ('x', 'f4'),\n    ('y', 'f4'),\n    ('z','f4'),\n    ('volume_id', 'i4'),\n    ('layer_id', 'i4'),\n    ('module_id', 'i4'),\n])\nPARTICLES_DTYPES = dict([\n    ('particle_id', 'i8'),\n    ('vx', 'f4'),\n    ('vy', 'f4'),\n    ('vz', 'f4'),\n    ('px', 'f4'),\n    ('py', 'f4'),\n    ('pz', 'f4'),\n    ('q', 'i4'),\n    ('nhits', 'i4'),\n])\nTRUTH_DTYPES = dict([\n    ('hit_id', 'i4'),\n    ('particle_id', 'i8'),\n    ('tx', 'f4'),\n    ('ty', 'f4'),\n    ('tz', 'f4'),\n    ('tpx', 'f4'),\n    ('tpy', 'f4'),\n    ('tpz', 'f4'),\n    ('weight', 'f4'),\n])\nDTYPES = {\n    'cells': CELLS_DTYPES,\n    'hits': HITS_DTYPES,\n    'particles': PARTICLES_DTYPES,\n    'truth': TRUTH_DTYPES,\n}\nDEFAULT_PARTS = ['hits', 'cells', 'particles', 'truth']\n\ndef _load_event_data(prefix, name):\n    \"\"\"Load per-event data for one single type, e.g. hits, or particles.\n    \"\"\"\n    expr = '{!s}-{}.csv*'.format(prefix, name)\n    files = glob.glob(expr)\n    dtype = DTYPES[name]\n    if len(files) == 1:\n        return pandas.read_csv(files[0], header=0, index_col=False, dtype=dtype)\n    elif len(files) == 0:\n        raise Exception('No file matches \\'{}\\''.format(expr))\n    else:\n        raise Exception('More than one file matches \\'{}\\''.format(expr))\n\ndef load_event_hits(prefix):\n    \"\"\"Load the hits information for a single event with the given prefix.\n    \"\"\"\n    return _load_event_data(prefix, 'hits')\n\ndef load_event_cells(prefix):\n    \"\"\"Load the hit cells information for a single event with the given prefix.\n    \"\"\"\n    return _load_event_data(prefix, 'cells')\n\ndef load_event_particles(prefix):\n    \"\"\"Load the particles information for a single event with the given prefix.\n    \"\"\"\n    return _load_event_data(prefix, 'particles')\n\ndef load_event_truth(prefix):\n    \"\"\"Load only the truth information for a single event with the given prefix.\n    \"\"\"\n    return _load_event_data(prefix, 'truth')\n\ndef load_event(prefix, parts=DEFAULT_PARTS):\n    \"\"\"Load data for a single event with the given prefix.\n    Parameters\n    ----------\n    prefix : str or pathlib.Path\n        The common prefix name for the event files, i.e. without `-hits.csv`).\n    parts : List[{'hits', 'cells', 'particles', 'truth'}], optional\n        Which parts of the event files to load.\n    Returns\n    -------\n    tuple\n        Contains a `pandas.DataFrame` for each element of `parts`. Each\n        element has field names identical to the CSV column names with\n        appropriate types.\n    \"\"\"\n    return tuple(_load_event_data(prefix, name) for name in parts)\n\ndef load_dataset(path, skip=None, nevents=None, parts=DEFAULT_PARTS):\n    \"\"\"Provide an iterator over (all) events in a dataset.\n    Parameters\n    ----------\n    path : str or pathlib.Path\n        Path to a directory or a zip file containing event files.\n    skip : int, optional\n        Skip the first `skip` events.\n    nevents : int, optional\n        Only load a maximum of `nevents` events.\n    parts : List[{'hits', 'cells', 'particles', 'truth'}], optional\n        Which parts of each event files to load.\n    Yields\n    ------\n    event_id : int\n        The event identifier.\n    *data\n        Event data element as specified in `parts`.\n    \"\"\"\n    # extract a sorted list of event file prefixes.\n    def list_prefixes(files):\n        regex = re.compile('^event\\d{9}-[a-zA-Z]+.csv')\n        files = filter(regex.match, files)\n        prefixes = set(op.basename(_).split('-', 1)[0] for _ in files)\n        prefixes = sorted(prefixes)\n        if skip is not None:\n            prefixes = prefixes[skip:]\n        if nevents is not None:\n            prefixes = prefixes[:nevents]\n        return prefixes\n\n    # TODO use yield from when we increase the python requirement\n    if op.isdir(path):\n        for x in _iter_dataset_dir(path, list_prefixes(os.listdir(path)), parts):\n            yield x\n    else:\n        with zipfile.ZipFile(path, mode='r') as z:\n            for x in _iter_dataset_zip(z, list_prefixes(z.namelist()), parts):\n                yield x\n\ndef _extract_event_id(prefix):\n    \"\"\"Extract event_id from prefix, e.g. event_id=1 from `event000000001`.\n    \"\"\"\n    return int(prefix[5:])\n\ndef _iter_dataset_dir(directory, prefixes, parts):\n    \"\"\"Iterate over selected events files inside a directory.\n    \"\"\"\n    for p in prefixes:\n        yield (_extract_event_id(p),) + load_event(op.join(directory, p), parts)\n\ndef _iter_dataset_zip(zipfile, prefixes, parts):\n    \"\"\"\"Iterate over selected event files inside a zip archive.\n    \"\"\"\n    for p in prefixes:\n        files = [zipfile.open('{}-{}.csv'.format(p, _), mode='r') for _ in parts]\n        dtypes = [DTYPES[_] for _ in parts]\n        data = tuple(pandas.read_csv(f, header=0, index_col=False, dtype=d)\n                                     for f, d in zip(files, dtypes))\n        yield (_extract_event_id(p),) + data\n        \n\"\"\"TrackML scoring metric\"\"\"\n\n__authors__ = ['Sabrina Amrouche', 'David Rousseau', 'Moritz Kiehn',\n               'Ilija Vukotic']\n\nimport numpy\nimport pandas\n\ndef _analyze_tracks(truth, submission):\n    \"\"\"Compute the majority particle, hit counts, and weight for each track.\n    Parameters\n    ----------\n    truth : pandas.DataFrame\n        Truth information. Must have hit_id, particle_id, and weight columns.\n    submission : pandas.DataFrame\n        Proposed hit/track association. Must have hit_id and track_id columns.\n    Returns\n    -------\n    pandas.DataFrame\n        Contains track_id, nhits, major_particle_id, major_particle_nhits,\n        major_nhits, and major_weight columns.\n    \"\"\"\n    # true number of hits for each particle_id\n    particles_nhits = truth['particle_id'].value_counts(sort=False)\n    total_weight = truth['weight'].sum()\n    # combined event with minimal reconstructed and truth information\n    event = pandas.merge(truth[['hit_id', 'particle_id', 'weight']],\n                         submission[['hit_id', 'track_id']],\n                         on=['hit_id'], how='left', validate='one_to_one')\n    event.drop('hit_id', axis=1, inplace=True)\n    event.sort_values(by=['track_id', 'particle_id'], inplace=True)\n\n    # ASSUMPTIONs: 0 <= track_id, 0 <= particle_id\n\n    tracks = []\n    # running sum for the reconstructed track we are currently in\n    rec_track_id = -1\n    rec_nhits = 0\n    # running sum for the particle we are currently in (in this track_id)\n    cur_particle_id = -1\n    cur_nhits = 0\n    cur_weight = 0\n    # majority particle with most hits up to now (in this track_id)\n    maj_particle_id = -1\n    maj_nhits = 0\n    maj_weight = 0\n\n    for hit in event.itertuples(index=False):\n        # we reached the next track so we need to finish the current one\n        if (rec_track_id != -1) and (rec_track_id != hit.track_id):\n            # could be that the current particle is the majority one\n            if maj_nhits < cur_nhits:\n                maj_particle_id = cur_particle_id\n                maj_nhits = cur_nhits\n                maj_weight = cur_weight\n            # store values for this track\n            tracks.append((rec_track_id, rec_nhits, maj_particle_id,\n                particles_nhits[maj_particle_id], maj_nhits,\n                maj_weight / total_weight))\n\n        # setup running values for next track (or first)\n        if rec_track_id != hit.track_id:\n            rec_track_id = hit.track_id\n            rec_nhits = 1\n            cur_particle_id = hit.particle_id\n            cur_nhits = 1\n            cur_weight = hit.weight\n            maj_particle_id = -1\n            maj_nhits = 0\n            maj_weights = 0\n            continue\n\n        # hit is part of the current reconstructed track\n        rec_nhits += 1\n\n        # reached new particle within the same reconstructed track\n        if cur_particle_id != hit.particle_id:\n            # check if last particle has more hits than the majority one\n            # if yes, set the last particle as the new majority particle\n            if maj_nhits < cur_nhits:\n                maj_particle_id = cur_particle_id\n                maj_nhits = cur_nhits\n                maj_weight = cur_weight\n            # reset runnig values for current particle\n            cur_particle_id = hit.particle_id\n            cur_nhits = 1\n            cur_weight = hit.weight\n        # hit belongs to the same particle within the same reconstructed track\n        else:\n            cur_nhits += 1\n            cur_weight += hit.weight\n\n    # last track is not handled inside the loop\n    if maj_nhits < cur_nhits:\n        maj_particle_id = cur_particle_id\n        maj_nhits = cur_nhits\n        maj_weight = cur_weight\n    # store values for the last track\n    tracks.append((rec_track_id, rec_nhits, maj_particle_id,\n        particles_nhits[maj_particle_id], maj_nhits, maj_weight / total_weight))\n\n    cols = ['track_id', 'nhits',\n            'major_particle_id', 'major_particle_nhits',\n            'major_nhits', 'major_weight']\n    return pandas.DataFrame.from_records(tracks, columns=cols)\n\ndef score_event(truth, submission):\n    \"\"\"Compute the TrackML event score for a single event.\n    Parameters\n    ----------\n    truth : pandas.DataFrame\n        Truth information. Must have hit_id, particle_id, and weight columns.\n    submission : pandas.DataFrame\n        Proposed hit/track association. Must have hit_id and track_id columns.\n    \"\"\"\n    tracks = _analyze_tracks(truth, submission)\n    purity_rec = numpy.divide(tracks['major_nhits'], tracks['nhits'])\n    purity_maj = numpy.divide(tracks['major_nhits'], tracks['major_particle_nhits'])\n    good_track_idx = (0.5 < purity_rec) & (0.5 < purity_maj)\n\t#good_track_df = tracks[good_track_idx]\n    return tracks['major_weight'][good_track_idx].sum()\n\t\n\t\ndef score_event_kha(truth, submission):\n    \"\"\"Compute the TrackML event score for a single event.\n    Parameters\n    ----------\n    truth : pandas.DataFrame\n        Truth information. Must have hit_id, particle_id, and weight columns.\n    submission : pandas.DataFrame\n        Proposed hit/track association. Must have hit_id and track_id columns.\n    \"\"\"\n\n    tracks = _analyze_tracks(truth, submission)\n    purity_rec = numpy.divide(tracks['major_nhits'], tracks['nhits'])\n    purity_maj = numpy.divide(tracks['major_nhits'], tracks['major_particle_nhits'])\n    good_track = (0.5 < purity_rec) & (0.5 < purity_maj)\n\t#sensentivity = \n    return tracks['major_weight'][good_track].sum()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"572fcbb6-8c7b-4a09-8916-8ec76689130f","_uuid":"63414de98667e95f60407c9155899a25a321cffc","trusted":true,"collapsed":true},"cell_type":"code","source":"path_to_train = \"../input/train_1\"\nevent_prefix = \"event000001000\"\nhits, 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":"def merge(cl1, cl2): # merge cluster 2 to cluster 1\n    d = pd.DataFrame(data={'s1':cl1,'s2':cl2})\n    d['N1'] = d.groupby('s1')['s1'].transform('count')\n    d['N2'] = d.groupby('s2')['s2'].transform('count')\n    maxs1 = d['s1'].max()\n    cond = np.where((d['N2'].values>d['N1'].values) & (d['N2'].values<25)) # tìm vị trí hit với nhit của cluster mới > nhits cluster cũ\n    s1 = d['s1'].values \n    s1[cond] = d['s2'].values[cond]+maxs1 # gán tất cả các hit đó thuộc về track mới (+maxs1 để tăng label cho track để nó khác ban đầu)\n    return s1\n\ndef extract_good_hits(truth, submission):\n    tru = truth[['hit_id', 'particle_id', 'weight']].merge(submission, how='left', on='hit_id')\n    tru['count_both'] = tru.groupby(['track_id', 'particle_id']).hit_id.transform('count')    \n    tru['count_particle'] = tru.groupby(['particle_id']).hit_id.transform('count')\n    tru['count_track'] = tru.groupby(['track_id']).hit_id.transform('count')\n    return tru[(tru.count_both > 0.5*tru.count_particle) & (tru.count_both > 0.5*tru.count_track)]\n\ndef fast_score(good_hits_df):\n    return good_hits_df.weight.sum()\n\n\ndef analyze_truth_perspective(truth, submission):\n    tru = truth[['hit_id', 'particle_id', 'weight']].merge(submission, how='left', on='hit_id')\n    tru['count_both'] = tru.groupby(['track_id', 'particle_id']).hit_id.transform('count')    \n    tru['count_particle'] = tru.groupby(['particle_id']).hit_id.transform('count')\n    tru['count_track'] = tru.groupby(['track_id']).hit_id.transform('count')\n    good_hits = tru[(tru.count_both > 0.5*tru.count_particle) & (tru.count_both > 0.5*tru.count_track)]\n    score = good_hits.weight.sum()\n    \n    anatru = tru.particle_id.value_counts().value_counts().sort_index().to_frame().rename({'particle_id':'true_particle_counts'}, axis=1)\n    #anatru['true_particle_ratio'] = anatru['true_particle_counts'].values*100/np.sum(anatru['true_particle_counts'])\n\n    anatru['good_tracks_counts'] = np.zeros(len(anatru)).astype(int)\n    anatru['good_tracks_intersect_nhits_avg'] = np.zeros(len(anatru))\n    anatru['best_detect_intersect_nhits_avg'] = np.zeros(len(anatru))\n    for nhit in tqdm(range(4,20)):\n        particle_list  = tru[(tru.count_particle==nhit)].particle_id.unique()\n        intersect_count = 0\n        good_tracks_count = 0\n        good_tracks_intersect = 0\n        for p in particle_list:\n            nhit_intersect = tru[tru.particle_id==p].count_both.max()\n            intersect_count += nhit_intersect\n            corresponding_track = tru.loc[tru[tru.particle_id==p].count_both.idxmax()].track_id\n            leng_corresponding_track = len(tru[tru.track_id == corresponding_track])\n            \n            if (nhit_intersect >= nhit/2) and (nhit_intersect >= leng_corresponding_track/2):\n                good_tracks_count += 1\n                good_tracks_intersect += nhit_intersect\n        intersect_count = intersect_count/len(particle_list)\n        anatru.at[nhit,'best_detect_intersect_nhits_avg'] = intersect_count\n        anatru.at[nhit,'good_tracks_counts'] = good_tracks_count\n        if good_tracks_count > 0:\n            anatru.at[nhit,'good_tracks_intersect_nhits_avg'] = good_tracks_intersect/good_tracks_count\n    \n    return score, anatru, good_hits\n\ndef precision(truth, submission,min_hits):\n    tru = truth[['hit_id', 'particle_id', 'weight']].merge(submission, how='left', on='hit_id')\n    tru['count_both'] = tru.groupby(['track_id', 'particle_id']).hit_id.transform('count')    \n    tru['count_particle'] = tru.groupby(['particle_id']).hit_id.transform('count')\n    tru['count_track'] = tru.groupby(['track_id']).hit_id.transform('count')\n    #print('Analyzing predictions...')\n    predicted_list  = tru[(tru.count_track>=min_hits)].track_id.unique()\n    good_tracks_count = 0\n    ghost_tracks_count = 0\n    fp_weights = 0\n    tp_weights = 0\n    for t in predicted_list:\n        nhit_track = tru[tru.track_id==t].count_track.iloc[0]\n        nhit_intersect = tru[tru.track_id==t].count_both.max()\n        corresponding_particle = tru.loc[tru[tru.track_id==t].count_both.idxmax()].particle_id\n        leng_corresponding_particle = len(tru[tru.particle_id == corresponding_particle])\n        if (nhit_intersect >= nhit_track/2) and (nhit_intersect >= leng_corresponding_particle/2): #if the predicted track is good\n            good_tracks_count += 1\n            tp_weights += tru[(tru.track_id==t)&(tru.particle_id==corresponding_particle)].weight.sum()\n            fp_weights += tru[(tru.track_id==t)&(tru.particle_id!=corresponding_particle)].weight.sum()\n        else: # if the predicted track is bad\n                ghost_tracks_count += 1\n                fp_weights += tru[(tru.track_id==t)].weight.sum()\n    all_weights = tru[(tru.count_track>=min_hits)].weight.sum()\n    precision = tp_weights/all_weights*100\n    print('Precision: ',precision,', good tracks:', good_tracks_count,', total tracks:',len(predicted_list),\n           ', loss:', fp_weights, ', reco:', tp_weights, 'reco/loss', tp_weights/fp_weights)\n    return precision\n\n\nclass Clusterer(object):\n    def __init__(self):                        \n        self.abc = []\n          \n    def initialize(self,dfhits):\n        self.cluster = range(len(dfhits))\n        \n    def Hough_clustering(self,dfh,coef,epsilon,min_samples=1,n_loop=180,verbose=True): # [phi_coef,phi_coef,zdivrt_coef,zdivr_coef,xdivr_coef,ydivr_coef]\n        merged_cluster = self.cluster\n        mm = 1\n        stepii = 0.000005\n        count_ii = 0\n        adaptive_eps_coefficient = 1\n        for ii in np.arange(0, n_loop*stepii, stepii):\n            count_ii += 1\n            for jj in range(2):\n                mm = mm*(-1)\n                eps_new = epsilon + count_ii*adaptive_eps_coefficient*10**(-5)\n                dfh['a1'] = dfh['a0'].values - np.nan_to_num(np.arccos(mm*ii*dfh['rt'].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','zdivrt','zdivr','xdivr','ydivr']].values) \n                #dfs = scale_ignore_nan(dfh[['sina1','cosa1','zdivrt','zdivr','xdivr','ydivr']])\n                dfs = np.multiply(dfs, coef)\n                new_cluster=DBSCAN(eps=eps_new,min_samples=min_samples,metric='euclidean',n_jobs=4).fit(dfs).labels_\n                merged_cluster = merge(merged_cluster, new_cluster)\n                if verbose == True:\n                    sub = create_one_event_submission(0, hits, merged_cluster)\n                    good_hits = extract_good_hits(truth, sub)\n                    score_1 = fast_score(good_hits)\n                    print('2r0_inverse:', ii*mm ,'. Score:', score_1)\n                    #clear_output(wait=True)\n        self.cluster = merged_cluster\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  \n\ndef preprocess_hits(h,dz):\n    h['z'] =  h['z'].values + dz\n    h['r'] = np.sqrt(h['x'].values**2+h['y'].values**2+h['z'].values**2)\n    h['rt'] = np.sqrt(h['x'].values**2+h['y'].values**2)\n    h['a0'] = np.arctan2(h['y'].values,h['x'].values)\n    h['zdivrt'] = h['z'].values/h['rt'].values\n    h['zdivr'] = h['z'].values/h['r'].values\n    h['xdivr'] = h['x'].values / h['r'].values\n    h['ydivr'] = h['y'].values / h['r'].values\n    return h","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"0a333cd4-351e-4274-aa7c-4cf8ab7fca1a","_uuid":"70ce31d93086e022159d6227f35c6488bf80eb22","trusted":true,"collapsed":true},"cell_type":"code","source":"# Clustering by varying \n#model = Clusterer()\n#model.initialize(hits) \nc = [1.5,1.5,0.73,0.17,0.027,0.027] #[phi_coef,phi_coef,zdivrt_coef,zdivr_coef,xdivr_coef,ydivr_coef]\nmin_samples_in_cluster = 1\n\nmodel = Clusterer()\nmodel.initialize(hits) \nhits_with_dz = preprocess_hits(hits, 0)\nmodel.Hough_clustering(hits_with_dz,coef=c,epsilon=0.0048,min_samples=min_samples_in_cluster,\n                       n_loop=300,verbose=True)\n\nsubmission = create_one_event_submission(0, hits, model.cluster)\nprint('\\n') ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2df8584057590d42adef24a1286a92599f886df3"},"cell_type":"markdown","source":"Now, let us see some analysis on the clustering result:"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"d496bcbe21a377925f64c69de5a6c811c80d58a7"},"cell_type":"code","source":"pr = precision(truth,submission,min_hits=10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"fc582df7cf36f13b9ee0331572db305521e6336e"},"cell_type":"code","source":"pr = precision(truth,submission,min_hits=7)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"8d22eb6747f01fa2af6bd2086cde122825e323a3"},"cell_type":"code","source":"pr = precision(truth,submission,min_hits=4)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"71aeb104416c23c2e57c45f7569939ff975780de"},"cell_type":"code","source":"pr = precision(truth,submission,min_hits=1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"184b5b3541d8f8a82bdce122b3dc9b041af0d984"},"cell_type":"markdown","source":"As one can see, long tracks have high precision, low lost weights. On the other hand, there are too many ghost short tracks. Then, we can use multi-stage clustering, using min_hits in DBSCAN for each stage (i.e., cluster long tracks first, then cluster short tracks with different parameters without touching the long tracks...)."},{"metadata":{"_uuid":"efd3174b9193b3180af1648d45d43db042e1fc55"},"cell_type":"markdown","source":"Some other notes:\n\n+ Use too many loops can decrease the performance, as one can see from the log result above.\n\n+ No z-shifting is performed  (dz = 0), although the function preprocess already offer it. Some may want to use z-shifting right away just by change dz from 0 to any number between [-5.5, 5.5]\n\n+ Features are not optimized. Honestly, I am also stuck at searching for good features (and good weights). It would be very nice if someone secretly tell me those magic features :-).\n\n+ When $r/(2r_0) > 1$ or $< 1$, arccos is undefined, hence a warning appears (if running on local notebook). The problem, more importantly, is not about the warning. It is a technical issue: all hits with $r/(2r_0) > 1$ or $< 1$ MUST BE EXCLUDED from DBSCAN, because there will be NO track with that parameter pass through those hits. This can be done by some indexing techniques that I do not provide here. (DBSCAN uses a raw matrix to cluster, then we must be careful when exclude hits from the original full hit dataframe).\n\nKV."},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"b8a951a9c4989186d8716a25e087fde36bcc30e2"},"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.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}