Determines the radius 'eps' for DBSCAN clustering of 'data' in an adaptive, data-dependent way.

    Parameters
    ----------
    hdf5_file_name : file object or string
        The handle or name of an HDF5 data structure where any array needed for DBSCAN
        and too large to fit into memory is to be stored.

    data : array of shape (n_samples, n_features)
        An array of features retained from the data-set to be analysed. 
        Subsamples of this curated data-set can also be analysed by a call to DBSCAN by providing an appropriate 
        list of selected samples labels, stored in 'subsamples_matrix' (see below).

    subsamples_matrix : array of shape (n_runs, n_subsamples), optional (default = None)
        Each row of this matrix contains a set of indices identifying the samples selected from the whole data-set
        for each of 'n_runs' independent rounds of DBSCAN clusterings.

    minPts : int
        The number of points within an epsilon-radius hypershpere for the said region to qualify as dense.

    eps : float, optional (default = None)
        Sets the maximum distance separating two data-points for those data-points to be considered 
        as part of the same neighborhood.

    quantile : int, optional (default = 50)
        If 'eps' is not provided by the user, it will be determined as the 'quantile' of the distribution 
        of the k-nearest distances to each sample, with k set to 'minPts'.

    samples_weights : array of shape (n_runs, n_samples), optional (default = None)
        Holds the weights of each sample. A sample with weight greater than 'minPts' is guaranteed to be
        a core sample; a sample with negative weight tends to prevent its 'eps'-neighbors from being core. 
        Weights are absolute and default to 1.

    metric : string or callable, optional (default = 'euclidean')
        The metric to use for computing the pairwise distances between samples 
        (each sample corresponds to a row in 'data'). If metric is a string or callable, it must be compatible 
        with metrics.pairwise.pairwise_distances.

    p : float, optional (default = 2)
        If a Minkowski metric is used, 'p' determines its power.

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

    Returns
    -------
    eps : float
        The parameter of DBSCAN clustering specifying if points are density-reachable. 
        This is either a copy of the value provided at input or, if the user did not specify a value of 'eps' at input, 
        the return value if the one determined from k-distance graphs from the data-set.

    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 load(hdf5_file_name, data, minPts, eps = None, quantile = 50, subsamples_matrix = None, samples_weights = None, 
metric = 'minkowski', p = 2, verbose = True):
    """Determines the radius 'eps' for DBSCAN clustering of 'data' in an adaptive, data-dependent way.

    Parameters
    ----------
    hdf5_file_name : file object or string
        The handle or name of an HDF5 data structure where any array needed for DBSCAN
        and too large to fit into memory is to be stored.

    data : array of shape (n_samples, n_features)
        An array of features retained from the data-set to be analysed. 
        Subsamples of this curated data-set can also be analysed by a call to DBSCAN by providing an appropriate 
        list of selected samples labels, stored in 'subsamples_matrix' (see below).

    subsamples_matrix : array of shape (n_runs, n_subsamples), optional (default = None)
        Each row of this matrix contains a set of indices identifying the samples selected from the whole data-set
        for each of 'n_runs' independent rounds of DBSCAN clusterings.

    minPts : int
        The number of points within an epsilon-radius hypershpere for the said region to qualify as dense.

    eps : float, optional (default = None)
        Sets the maximum distance separating two data-points for those data-points to be considered 
        as part of the same neighborhood.

    quantile : int, optional (default = 50)
        If 'eps' is not provided by the user, it will be determined as the 'quantile' of the distribution 
        of the k-nearest distances to each sample, with k set to 'minPts'.

    samples_weights : array of shape (n_runs, n_samples), optional (default = None)
        Holds the weights of each sample. A sample with weight greater than 'minPts' is guaranteed to be
        a core sample; a sample with negative weight tends to prevent its 'eps'-neighbors from being core. 
        Weights are absolute and default to 1.

    metric : string or callable, optional (default = 'euclidean')
        The metric to use for computing the pairwise distances between samples 
        (each sample corresponds to a row in 'data'). If metric is a string or callable, it must be compatible 
        with metrics.pairwise.pairwise_distances.

    p : float, optional (default = 2)
        If a Minkowski metric is used, 'p' determines its power.

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

    Returns
    -------
    eps : float
        The parameter of DBSCAN clustering specifying if points are density-reachable. 
        This is either a copy of the value provided at input or, if the user did not specify a value of 'eps' at input, 
        the return value if the one determined from k-distance graphs from the data-set.

    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
    """
    
    data = np.array(data, copy = False)
    if data.ndim > 2:
        raise ValueError("\nERROR: DBSCAN_multiplex @ load:\n" 
                         "the data array is of dimension %d. Please provide a two-dimensional "
                         "array instead.\n" % data.ndim)

    if subsamples_matrix is None:
        subsamples_matrix = np.arange(data.shape[0], dtype = int)
        subsamples_matrix = subsamples_matrix.reshape(1, -1)
 
    else:
        subsamples_matrix = np.array(subsamples_matrix, copy = False)

    if subsamples_matrix.ndim > 2:
        raise ValueError("\nERROR: DBSCAN_multiplex @ load:\n"
                         "the array of subsampled indices is of dimension %d. "
                         "Please provide a two-dimensional array instead.\n" % subsamples_matrix.ndim)

    if (data.dtype.char in np.typecodes['AllFloat'] and not np.isfinite(data.sum()) and not np.all(np.isfinite(data))):
        raise ValueError('\nERROR: DBSCAN_multiplex @ load:\n'
                         'the data vector contains at least one infinite or NaN entry.\n')

    if (subsamples_matrix.dtype.type is np.int_ and not np.isfinite(subsamples_matrix.sum()) and not np.all(np.isfinite(subsamples_matrix))):
        raise ValueError('\nERROR: DBSCAN_multiplex @ load:\n' 
                         'the array of subsampled indices contains at least one infinite or NaN entry.\n')

    if not np.all(subsamples_matrix >= 0):
        raise ValueError('\nERROR: DBSCAN_multiplex @ load:\n'
                         'the sampled indices should all be positive integers.\n') 

    N_samples = data.shape[0]
    N_runs, N_subsamples = subsamples_matrix.shape

    if N_subsamples > N_samples:
        raise ValueError('\nERROR: DBSCAN_multiplex @ load:\n'
                         'the number of sampled indices cannot exceed the total number of samples in the whole data-set.\n')

    for i in range(N_runs):
        subsamples_matrix[i] = np.unique(subsamples_matrix[i])
 
    if not isinstance(minPts, int):
        raise TypeError("\nERROR: DBSCAN_multiplex @ load:\n"
                        "the parameter 'minPts' must be an integer.\n")

    if minPts < 2:
        raise ValueError("\nERROR: DBSCAN_multiplex @ load:\n"
                         "the value of 'minPts' must be larger than 1.\n")        

    if eps is None:
        # Determine the parameter 'eps' as the median of the distribution
        # of the maximum of the minPts-nearest neighbors distances for each sample.
        if verbose:
            print(("INFO: DBSCAN_multiplex @ load:\n"
                  "starting the determination of an appropriate value of 'eps' for this data-set"
                  " and for the other parameter of the DBSCAN algorithm set to {minPts}.\n"
                  "This might take a while.".format(**locals())))

        beg_eps = time.time()

        quantile = np.rint(quantile)
        quantile = np.clip(quantile, 0, 100)

        k_distances = kneighbors_graph(data, minPts, mode = 'distance', metric = metric, p = p).data
 
        radii = np.zeros(N_samples, dtype = float)
        for i in range(0, minPts):
            radii = np.maximum(radii, k_distances[i::minPts]) 
             
        if quantile == 50:     
            eps = round(np.median(radii, overwrite_input = True), 4)
        else:
            eps = round(np.percentile(radii, quantile), 4)

        end_eps = time.time()

        if verbose:
            print(("\nINFO: DBSCAN_multiplex @ load:\n"
                  "done with evaluating parameter 'eps' from the data-set provided."
                  " This took {} seconds. Value of epsilon: {}.".format(round(end_eps - beg_eps, 4), eps)))

    else:
        if not (isinstance(eps, float) or isinstance(eps, int)):
            raise ValueError("\nERROR: DBSCAN_multiplex @ load:\n"
                             "please provide a numeric value for the radius 'eps'.\n")

        if not eps > 0.0:
            raise ValueError("\nERROR: DBSCAN_multiplex @ load:\n"
                             "the radius 'eps' must be positive.\n")

        eps = round(eps, 4)

    # For all samples with a large enough neighborhood, 'neighborhoods_indices' 
    # and 'neighborhoods_indptr' help us find the neighbors to every sample. Note
    # that this definition of neighbors leaves the original point in,
    # which will be considered later.
    if verbose:
       print(("\nINFO: DBSCAN_multiplex @ load:\n"
             "identifying the neighbors within an hypersphere of radius {eps} around each sample,"
             " while at the same time evaluating the number of epsilon-neighbors for each sample.\n"
             "This might take a fair amount of time.".format(**locals())))

    beg_neigh = time.time()

    fileh = tables.open_file(hdf5_file_name, mode = 'r+')
    DBSCAN_group = fileh.create_group(fileh.root, 'DBSCAN_group')

    neighborhoods_indices = fileh.create_earray(DBSCAN_group, 'neighborhoods_indices', tables.Int32Atom(), (0,), 
                                                'Indices array for sparse matrix of neighborhoods', 
                                                expectedrows = int((N_samples ** 2) / 50))

    # 'neighborhoods_indptr' is such that for each of row i of the data-matrix 
    # neighborhoods_indices[neighborhoods_indptr[i]:neighborhoods_indptr[i+1]]
    # contains the column indices of row i from the array of 
    # 'eps'-neighborhoods.
    neighborhoods_indptr = np.zeros(1, dtype = np.int64)

    # For each sample, 'neighbors_counts' will keep a tally of the number 
    # of its  neighbors within a hypersphere of radius 'eps'. 
    # Note that the sample itself is counted as part of this neighborhood.
    neighbors_counts = fileh.create_carray(DBSCAN_group, 'neighbors_counts', tables.Int32Atom(), (N_runs, N_samples), 
                                           'Array of the number of neighbors around each sample of a set of subsampled points', 
                                           filters = None)   

    chunks_size = get_chunk_size(N_samples, 3)
    for i in range(0, N_samples, chunks_size):
        chunk = data[i:min(i + chunks_size, N_samples)]

        D = pairwise_distances(chunk, data, metric = metric, p = p, n_jobs = 1)
            
        D = (D <= eps)

        if samples_weights is None:
            for run in range(N_runs):
                x = subsamples_matrix[run]
                M = np.take(D, x, axis = 1)

                legit_rows = np.intersect1d(i + np.arange(min(chunks_size, N_samples - i)), x, assume_unique = True)
                M = np.take(M, legit_rows - i, axis = 0)
                
                neighbors_counts[run, legit_rows] = M.sum(axis = 1)

                del M
        else:
            for run in range(N_runs):
                x = subsamples_matrix[run]

                M = np.take(D, x, axis = 1)

                legit_rows = np.intersect1d(i + np.arange(min(chunks_size, N_samples - i)), x, assume_unique = True)
                M = np.take(M, legit_rows - i, axis = 0)

                neighbors_counts[run, legit_rows] = np.array([np.sum(samples_weights[x[row]]) for row in M])

                del M

        candidates = np.where(D == True)

        del D

        neighborhoods_indices.append(candidates[1])

        _, nbr = np.unique(candidates[0], return_counts = True)
        counts = np.cumsum(nbr) + neighborhoods_indptr[-1]

        del candidates

        neighborhoods_indptr = np.append(neighborhoods_indptr, counts)

    fileh.create_carray(DBSCAN_group, 'neighborhoods_indptr', tables.Int64Atom(), (N_samples + 1,), 
                        'Array of cumulative number of column indices for each row', filters = None)
    fileh.root.DBSCAN_group.neighborhoods_indptr[:] = neighborhoods_indptr[:]

    fileh.create_carray(DBSCAN_group, 'subsamples_matrix', tables.Int32Atom(), (N_runs, N_subsamples), 
                        'Array of subsamples indices', filters = None)
    fileh.root.DBSCAN_group.subsamples_matrix[:] = subsamples_matrix[:]

    fileh.close()

    end_neigh = time.time()

    if verbose:
        print(("\nINFO: DBSCAN_multiplex @ load:\n"
              "done with the neighborhoods. This step took {} seconds.".format(round(end_neigh - beg_neigh, 4))))

    gc.collect()

    return eps