Perform DBSCAN clustering with parameters 'minPts' and 'eps'
        (as determined by a prior call to 'load' from this module). 
        If multiple subsamples of the dataset were provided in a preliminary call to 'load', 
        'sample_ID' specifies which one of those subsamples is to undergo DBSCAN clustering. 
         
    Parameters
    ----------
    hdf5_file_name : file object or string
        The handle or name of an HDF5 file where any array needed for DBSCAN and too large to fit into memory 
        is to be stored. Procedure 'shoot' relies on arrays stored in this data structure by a previous
        call to 'load' (see corresponding documentation)

    sample_ID : int, optional (default = 0)
        Identifies the particular set of selected data-points on which to perform DBSCAN.
        If not subsamples were provided in the call to 'load', the whole dataset will be subjected to DBSCAN clustering.

    minPts : int
        The number of points within a 'eps'-radius hypershpere for this region to qualify as dense.

    random_state: np.RandomState, optional (default = None)
        The generator used to reorder the samples. If None at input, will be set to np.random.

    verbose : Boolean, optional (default = True)
        Whether to display messages concerning the status of the computations and the time it took to complete 
        each major stage of the algorithm.

    Returns
    -------
    core_samples : array of shape (n_core_samples, )
        Indices of the core samples.

    labels : array of shape (N_samples, ) 
        Holds the cluster labels of each sample. The points considered as noise have entries -1. 
        The points not initially selected for clustering (i.e. not listed in 'subsampled_indices', if the latter
        has been provided in the call to 'load' from this module) are labelled -2.

    References
    ----------
    Ester, M., H. P. Kriegel, J. Sander and X. Xu, "A Density-Based
    Algorithm for Discovering Clusters in Large Spatial Databases with Noise".
    In: Proceedings of the 2nd International Conference on Knowledge Discovery
    and Data Mining, Portland, OR, AAAI Press, pp. 226-231. 1996

def shoot(hdf5_file_name, minPts, sample_ID = 0, random_state = None, verbose = True): 
    """Perform DBSCAN clustering with parameters 'minPts' and 'eps'
        (as determined by a prior call to 'load' from this module). 
        If multiple subsamples of the dataset were provided in a preliminary call to 'load', 
        'sample_ID' specifies which one of those subsamples is to undergo DBSCAN clustering. 
         
    Parameters
    ----------
    hdf5_file_name : file object or string
        The handle or name of an HDF5 file where any array needed for DBSCAN and too large to fit into memory 
        is to be stored. Procedure 'shoot' relies on arrays stored in this data structure by a previous
        call to 'load' (see corresponding documentation)

    sample_ID : int, optional (default = 0)
        Identifies the particular set of selected data-points on which to perform DBSCAN.
        If not subsamples were provided in the call to 'load', the whole dataset will be subjected to DBSCAN clustering.

    minPts : int
        The number of points within a 'eps'-radius hypershpere for this region to qualify as dense.

    random_state: np.RandomState, optional (default = None)
        The generator used to reorder the samples. If None at input, will be set to np.random.

    verbose : Boolean, optional (default = True)
        Whether to display messages concerning the status of the computations and the time it took to complete 
        each major stage of the algorithm.

    Returns
    -------
    core_samples : array of shape (n_core_samples, )
        Indices of the core samples.

    labels : array of shape (N_samples, ) 
        Holds the cluster labels of each sample. The points considered as noise have entries -1. 
        The points not initially selected for clustering (i.e. not listed in 'subsampled_indices', if the latter
        has been provided in the call to 'load' from this module) are labelled -2.

    References
    ----------
    Ester, M., H. P. Kriegel, J. Sander and X. Xu, "A Density-Based
    Algorithm for Discovering Clusters in Large Spatial Databases with Noise".
    In: Proceedings of the 2nd International Conference on Knowledge Discovery
    and Data Mining, Portland, OR, AAAI Press, pp. 226-231. 1996
    """       
        
    fileh = tables.open_file(hdf5_file_name, mode = 'r+')

    neighborhoods_indices = fileh.root.DBSCAN_group.neighborhoods_indices
    neighborhoods_indptr = fileh.root.DBSCAN_group.neighborhoods_indptr[:]

    neighbors_counts = fileh.root.DBSCAN_group.neighbors_counts[sample_ID]
    subsampled_indices = fileh.root.DBSCAN_group.subsamples_matrix[sample_ID]

    N_samples = neighborhoods_indptr.size - 1
    N_runs, N_subsamples = fileh.root.DBSCAN_group.subsamples_matrix.shape

    if not isinstance(sample_ID, int):
        raise ValueError("\nERROR: DBSCAN_multiplex @ shoot:\n"
                         "'sample_ID' must be an integer identifying the set of subsampled indices "
                         "on which to perform DBSCAN clustering\n")    

    if (sample_ID < 0) or (sample_ID >= N_runs):
        raise ValueError("\nERROR: DBSCAN_multiplex @ shoot:\n"
                 "'sample_ID' must belong to the interval [0; {}].\n".format(N_runs - 1))
      
    # points that have not been sampled are labelled with -2
    labels = np.full(N_samples, -2, dtype = int)
    # among the points selected for clustering, 
    # all are initally characterized as noise
    labels[subsampled_indices] = - 1
    
    random_state = check_random_state(random_state)

    core_samples = np.flatnonzero(neighbors_counts >= minPts)
 
    index_order = np.take(core_samples, random_state.permutation(core_samples.size))

    cluster_ID = 0

    # Look at all the selected samples, see if they qualify as core samples
    # Create a new cluster from those core samples
    for index in index_order:
        if labels[index] not in {-1, -2}:
            continue

        labels[index] = cluster_ID

        candidates = [index]
        while len(candidates) > 0:
            candidate_neighbors = np.zeros(0, dtype = np.int32)
            for k in candidates:
                candidate_neighbors = np.append(candidate_neighbors, 
                                                neighborhoods_indices[neighborhoods_indptr[k]: neighborhoods_indptr[k+1]])
                candidate_neighbors = np.unique(candidate_neighbors)

            candidate_neighbors = np.intersect1d(candidate_neighbors, subsampled_indices, assume_unique = True)
                
            not_noise_anymore = np.compress(np.take(labels, candidate_neighbors) == -1, candidate_neighbors)
            
            labels[not_noise_anymore] = cluster_ID

            # Eliminate as potential candidates the points that have already 
            # been used to expand the current cluster by a trail 
            # of density-reachable points
            candidates = np.intersect1d(not_noise_anymore, core_samples, assume_unique = True) 
     
        cluster_ID += 1
    # Done with building this cluster. 
    # "cluster_ID" is now labelling the next cluster.

    fileh.close()

    gc.collect()

    return core_samples, labels