Extract bin subnetworks from the main network

    Identify bins, extract subnets, draws the adjacency matrices,
    saves it all in a specified output directory.

    Parameters
    ----------
    partition_file : file, str or pathlib.Path
        The file containing, for each chunk, the communities it was
        assigned to at each iteration.
    network_file : file, str or pathlib.Path
        The file containing the network in sparse (edge list) format
    output_dir : str or pathlib.Path
        The output directory to write the subnetworks into.
    max_cores : int, optional
        The maximum number of bins to extract. Default is 100.
    max_size_matrix : int, optional
        When rendering contact maps for each bin, the maximum size for the
        matrix. Default is 2000.
    saturation_threshold : float, optional
        When rendering contact maps for each bin, the percentile value over
        which the color map should be saturated. Default is 80.

def extract_subnetworks(
    partition_file,
    network_file,
    output_dir,
    max_cores=DEFAULT_MAX_CORES,
    max_size_matrix=DEFAULT_MAX_SIZE_MATRIX,
    saturation_threshold=DEFAULT_SATURATION_THRESHOLD,
):
    """Extract bin subnetworks from the main network

    Identify bins, extract subnets, draws the adjacency matrices,
    saves it all in a specified output directory.

    Parameters
    ----------
    partition_file : file, str or pathlib.Path
        The file containing, for each chunk, the communities it was
        assigned to at each iteration.
    network_file : file, str or pathlib.Path
        The file containing the network in sparse (edge list) format
    output_dir : str or pathlib.Path
        The output directory to write the subnetworks into.
    max_cores : int, optional
        The maximum number of bins to extract. Default is 100.
    max_size_matrix : int, optional
        When rendering contact maps for each bin, the maximum size for the
        matrix. Default is 2000.
    saturation_threshold : float, optional
        When rendering contact maps for each bin, the percentile value over
        which the color map should be saturated. Default is 80.
    """

    logger.info("Loading partition...")
    data_chunks = np.loadtxt(partition_file, usecols=(1,), dtype=np.int32)

    logger.info("Loading network...")
    network = np.loadtxt(network_file, dtype=np.int32)
    cores = data_chunks

    core_network = np.copy(network)

    core_network[:, 0] = cores[network[:, 0]]
    core_network[:, 1] = cores[network[:, 1]]

    n = np.amax(cores) + 1

    def extract(network_to_keep, filename):

        subnetwork = np.copy(network[network_to_keep])
        subnetwork[:, 0] -= 1
        subnetwork[:, 1] -= 1

        np.savetxt(filename, subnetwork, fmt="%i")

        return subnetwork

    def draw(subnetwork, filename):

        try:
            # Numpy array format
            row = subnetwork[:, 0]
            col = subnetwork[:, 1]
            data = subnetwork[:, 2]
        except TypeError:
            # Scipy sparse format
            row = subnetwork.row
            col = subnetwork.col
            data = subnetwork.data

        row_indices = stats.rankdata(
            np.concatenate((row, col)), method="dense"
        )
        col_indices = stats.rankdata(
            np.concatenate((col, row)), method="dense"
        )
        data = np.concatenate((data, data))

        # print("Row length: {}, col length: {}, data length: {}"
        #       "".format(len(row_indices), len(col_indices), len(data)))

        unique_row = np.unique(row)
        unique_col = np.unique(col)

        # print("Network shape: {},{}".format(len(unique_row),
        #                                     len(unique_col)))

        size = len(np.unique(np.concatenate((unique_row, unique_col)))) + 1
        # print("Size of matrix to draw: {}".format(size))

        try:
            sparse_subnet = sparse.coo_matrix(
                (data, (row_indices, col_indices)), shape=(size, size)
            )
            binning_factor = (size // max_size_matrix) + 1
            binned_subnet = hcs.bin_sparse(
                sparse_subnet, subsampling_factor=binning_factor
            )
            dense_subnet = binned_subnet.todense()

            diagonal = np.diag(np.diag(dense_subnet))
            normed_subnet = hcs.normalize_dense(dense_subnet - diagonal)

            vmax = np.percentile(normed_subnet, saturation_threshold)

            spaceless_pdf_plot_maker(normed_subnet, filename, vmax=vmax)

        except MemoryError:
            logger.warning(
                "Warning, couldn't save matrix due to memory issues"
            )

    def extract_and_draw(network_to_keep, filename_text, filename_image):

        subnetwork = extract(network_to_keep, filename=filename_text)
        draw(subnetwork, filename=filename_image)

    #   Extract and draw subnetworks for chosen cores and draw 2D arrays
    global_network_indices_list = []

    for i in range(1, n):

        if i > max_cores:
            break

        # print("Bin {}:".format(i))
        network_to_keep_1 = core_network[:, 0] == i
        network_to_keep_2 = core_network[:, 1] == i

        network_to_keep = network_to_keep_1 * network_to_keep_2

        nonzero_indices, = np.nonzero(network_to_keep)
        global_network_indices_list += nonzero_indices.tolist()

        subnetwork_file = os.path.join(
            output_dir, "subnetwork_core_{}.dat".format(i)
        )

        image_name = os.path.join(output_dir, "core_{}.eps".format(i))

        extract_and_draw(
            network_to_keep=network_to_keep,
            filename_text=subnetwork_file,
            filename_image=image_name,
        )