def METHOD_NAME(am, ploidy, single_linkage=False):
    """
    Based on the connectivity of the reads, we want to divide the phasing input, as non- or poorly
    connected regions can be phased independently. This is done based on how pairs of variants are
    connected. There are two modes how to decide whether two variants are connected:
    single_linkage=True -- If there exists a read in the allele matrix which covers both variants,
                           they count as connected
    single_linkage=False -- In order to connect two variants, we need at least reads from ploidy-1
                            different haplotypes. Two variants count as connected, if there
                            sufficiently many reads covering both variants, with "sufficient"
                            meaning, that the connecting reads have a chance of at least 98% that
                            they cover at least ploidy-1 haplotypes.
    First, only consecutive pairs are inspected. Then, this connectivity is made transitive, i.e.
    if the pair (A,C) is connected, as well as the pair (B,C), then (A,B) is also connected. If the
    special case occurs, that for three variants (in this order) A and C are connected, but neither
    is connected to variant B in between them, then the variants are still divided as A|B|C, even
    though A and C are actually connected. This is because the following steps require the splits
    to be intervals with no "holes" inside them.
    """
    num_vars = am.getNumPositions()
    # special case
    if num_vars == 0:
        return []
    # cut threshold
    if ploidy == 2 or single_linkage:
        cut_threshold = 1
    else:
        cut_threshold = ploidy * ploidy
        for i in range(ploidy - 1, ploidy * ploidy):
            # chance to cover at most ploidy-2 haplotypes with i reads
            cut_threshold = i
            if ploidy * pow((ploidy - 2) / ploidy, i) < 0.02:
                cut_threshold = i
                break
    logger.debug(f"Cut position threshold: coverage >= {cut_threshold}")
    # start by looking at neighbouring
    link_to_next = [0 for i in range(num_vars)]
    for read in am:
        pos_list = [pos for (pos, allele) in read]
        for i in range(len(pos_list) - 1):
            if pos_list[i] + 1 == pos_list[i + 1]:
                link_to_next[pos_list[i]] += 1
    pos_clust = [0 for i in range(num_vars)]
    for i in range(1, num_vars):
        if link_to_next[i - 1] >= cut_threshold:
            pos_clust[i] = pos_clust[i - 1]
        else:
            pos_clust[i] = pos_clust[i - 1] + 1
    num_clust = pos_clust[-1] + 1
    # find linkage between clusters
    link_coverage = [defaultdict(int) for i in range(num_clust)]
    for i, read in enumerate(am):
        covered_pos_clusts = {pos_clust[pos] for (pos, allele) in read}
        for p1 in covered_pos_clusts:
            for p2 in covered_pos_clusts:
                link_coverage[p1][p2] += 1
    # merge clusters
    merged_clust = [-1 for i in range(num_clust)]
    new_num_clust = 0
    for i in range(num_clust):
        if merged_clust[i] >= 0:
            continue
        q = Queue()
        q.put(i)
        merged_clust[i] = new_num_clust
        while not q.empty():
            cur = q.get()
            for linked in link_coverage[cur]:
                if merged_clust[linked] < 0 and link_coverage[cur][linked] >= cut_threshold:
                    q.put(linked)
                    merged_clust[linked] = new_num_clust
        new_num_clust += 1
    # determine cut positions
    cuts = [0]
    for i in range(1, num_vars):
        if merged_clust[pos_clust[i]] != merged_clust[pos_clust[i - 1]]:
            cuts.append(i)
    return cuts