The main granular assembling function, that IDENTIFIES assembly points, but does not perform the assembly on its own
    It DOES NOT change the supplied breakpoint graph in any way!!!

def identify_assembly_points(graph, bgtree, target_multicolor, exclude=None, verbose=False, verbose_destination=None):
    """
    The main granular assembling function, that IDENTIFIES assembly points, but does not perform the assembly on its own
    It DOES NOT change the supplied breakpoint graph in any way!!!
    """
    if verbose:
        print(">>Identifying assemblies for target multicolor:",
              [e.name for e in target_multicolor.multicolors.elements()], file=verbose_destination)

    guidance = bgtree.consistent_multicolors[:]
    offset = len(Multicolor.split_colors(target_multicolor, guidance=guidance,
                                         account_for_color_multiplicity_in_guidance=False)) - 1

    threshold = 1 if offset == 0 else 2
    assemblies = []  # the overall result
    if exclude is None:
        exclude = []  # a container with single colors of genomes, that are to be considered fully assembled

    p_t_consistent_multicolors_in_target = Multicolor.split_colors(target_multicolor, guidance=guidance,
                                                                   account_for_color_multiplicity_in_guidance=False)
    t_consistent_multicolors_in_target = []
    for tcmc in p_t_consistent_multicolors_in_target:
        t_consistent_multicolors_in_target.append(sorted(color.name for color in tcmc.colors))
    # we work with each connected component separately, as connected components usually preserve fragmentation points
    ################################################################################################
    #
    # its important that we iterate over connected components making each particular one a deepcopy of an
    #   underlying breakpoint graph connected component
    #
    ################################################################################################
    for i, cc in enumerate(graph.connected_components_subgraphs(copy=True)):

        # we filter current connected component of uninteresting / ambiguous edges and retrieve a list of
        #   connected components that are left in the original connected components after filtration
        irregular_subnets = get_irregular_subnets(cc, target_multicolor, exclude)

        if len(irregular_subnets) > 0 and verbose:
            print(">>Processing", str(i) + "th", "connected component", file=verbose_destination)
            print("\tcontains", len(irregular_subnets), "subnet groups", file=verbose_destination)

        # each subnet can be processed separately
        for subnet in irregular_subnets:

            supporting_edge_scores = get_support_edge_scores(graph, subnet, target_multicolor, bgtree)

            # we create a new dummy graph for the purpose of computing maximum weight matching for support edges in it
            new_graph = nx.Graph()
            if verbose:
                print("\tcontains", len(supporting_edge_scores), "possible assembly points", file=verbose_destination)

            # we'll keep track of possible assembly points for future reference
            support_edge_dict = {}
            for (v1, v2), before, after, ex_data in supporting_edge_scores:
                ex_data["tcmc"] = t_consistent_multicolors_in_target
                ##########################################################################################
                #
                # INSERT YOUR CODE ASSEMBLY SCORE THRESHOLD FILTRATION HERE IF NEED BE
                #
                ##########################################################################################
                if before - after - offset < threshold:
                    continue
                ##########################################################################################
                #
                # by default networkx assumes all edges, that have weight >= 0 are good
                #
                ##########################################################################################
                new_graph.add_edge(v1, v2, weight=before - after - offset)
                support_edge_dict[(v1, v2)] = (before, after + offset, ex_data)
                support_edge_dict[(v2, v1)] = (before, after + offset, ex_data)

            maximal_matching = nx.max_weight_matching(new_graph)

            if verbose:
                print("\t", len(maximal_matching) // 2, "assembly points are identified", file=verbose_destination)

            # as networkx provides a maximum matching in a form of adjacency list, every identified edge
            #   (pair of vertices) is present their twice (i.e. matching[v1]=v2 and matching[v2]=v1)
            #   we need to make sure we only add every edge only once
            visited = set()
            for v1, v2 in maximal_matching.items():
                if v1 in visited or v2 in visited:
                    continue
                visited.add(v1)
                visited.add(v2)
                assemblies.append((v1, v2, support_edge_dict[(v1, v2)]))

    # we return the result as a list of assembly points that were identified for the targeted multicolor
    # as a list of tuples (v1, v2, ("before", "after"))
    # where v1 and v2 correspond to assembly point and "before" and "after" are used to compute the assembly score
    return assemblies