{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.\n\nimport matplotlib.pyplot as plt","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"4ccfe1037f352e6fcc6c48a9f2c8ef4446c45f79","_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# from https://www.kaggle.com/khahuras/0-53x-clustering-using-hough-features-basic\n\n\"\"\"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":{"trusted":true,"collapsed":true,"_uuid":"76b23b4782632a668a41c95f568b4f4ae53622d2"},"cell_type":"code","source":"# from https://www.kaggle.com/khahuras/0-53x-clustering-using-hough-features-basic\n\npath_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":{"trusted":true,"_uuid":"740ae18e11deb8e39e4c2251f34f09d9a52e079d"},"cell_type":"code","source":"_ = plt.hist(particles.vz[particles.vz <= 5][particles.vz >= -5], bins=11)","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}