Sharedepistasis_and_posets_public.sagewsOpen in CoCalc
#---------------------------------------------------------------------------------------------------------------#
# prop_extensions_implying_epistasis(poset,partition, N, P ), a function for determining the proportion of linear orders which imply epistasis, and its helper functions.

#given a partition of a poset (or a labelling of a poset with two charcters), in the form of a dictionary, and the set of linear extensions of a poset
#as a list of lists, return all words generated by all linear extensions
def words_from_orders(partition, list_of_total_orders):
    words = []
    for order in list_of_total_orders:
        word = ''
        for element in order:
            word=word+str(partition[element])
        words.append(word)
    return words

#given a string and two characters, determine if string is a Dyck word on the two characters
def is_dyck_word(string, letter1, letter2):
    count1 = 0
    count2 = 0
    for char in string:
        if char == letter1:
            count1 = count1+1
        if char == letter2:
            count2 = count2+1
        if count2 > count1:
            return false
    if (count1 == count2):
        return true
    else: return false

#input: a partial order, a partition, and the negative and positive labels in the partition
#output: a tuple (npositive, nnegative, ninconclusive) where npositive is the number of linear extensions implying pos epistasis, nnegative is the prop. of linear extensions implying negative epistasis, and ninconclusive is the number or linear extensions that do not imply epistasis
def prop_extensions_implying_epistasis(poset,partition, N, P ):
    npositive = 0
    nnegative = 0
    ninconclusive = 0
    extensions = poset.linear_extensions().list()
    words = words_from_orders( partition, extensions)
    for word in words:
        if is_dyck_word(word, N, P):
            npositive = npositive+ 1
        elif is_dyck_word(word, P, N):
            nnegative = nnegative + 1
        else:
            ninconclusive = ninconclusive + 1
    return (npositive, nnegative, ninconclusive)



#---------------------------------------------------------------------------------------------------------------------------------#
#Some functions for working with partial orders that are useful in generating and analyzing examples
# a function for computing the intersection of two partial orders, and a function for finding the linear extensions graph of a poset

#imput: two partial orders
#output: the intersection of the partial orders--that is, the partial order p in which x<y if and only if x<y in p1 and p2
def poset_intersection(poset1, poset2):
    #if poset1.list() != poset2.list():
       # return "error, posets are not on the same set of elements"
    elements = poset1.list()
    list_of_edges = []
    for x in elements:
        for y in elements:
            if poset1.is_less_than(x,y) and poset2.is_less_than(x,y):
                list_of_edges.append((x,y))
    return Poset(DiGraph([elements,list_of_edges], format = 'vertices_and_edges'))

def intersection_gives_digraph(poset1, poset2):
    elements = poset1.list()
    list_of_edges = []
    for x in elements:
        for y in elements:
            if poset1.is_less_than(x,y) and poset2.is_less_than(x,y):
                list_of_edges.append((x,y))
    return DiGraph([elements,list_of_edges], format = 'vertices_and_edges')


#checks if two total orders differ by one transposition or less
def check_adjacency(list1, list2):
    numswaps = 0
    i=0
    while i<(len(list1)-1):
        if list1[i] == list2[i]:
            i = i + 1
        if list1[i] != list2[i]:
            if numswaps == 0 and list1[i]==list2[i+1]:
                if list1[i+1]==list2[i]:
                    numswaps = numswaps+1
                    i = i + 2
                else: return false
            else: return false
    return true


#makes a list of list into a list of strings
def to_string(list_of_lists):
    strings = []
    for list in list_of_lists:
        strings.append(str(list))
    return strings

#input: a poset
#output: a graph whose vertices are the linear extensions of the poset. there is an edge on a pair of vertices if they are connected by a single adjacent transposition
def linear_extensions_graph(poset):
    extensions = poset.linear_extensions().list()
    vertices = to_string(extensions)
    print ('u', vertices)
    edges = []
    for ext1 in extensions:
        for ext2 in extensions:
            if check_adjacency(ext1, ext2):
                edges.append((str(ext1),str(ext2)))
    return DiGraph([vertices, edges])

#input: a list
#output: a finite poset which is the rank order from this list
def rank_order_from_list(list):
    vertices = list
    edges = []
    for i in range (len(list)-1):
        edges.append((list[i], list[i+1]))
    return Poset(DiGraph([vertices, edges]))

#Our efficient implementation for checking whether a partial order implies epistasis, and its helper functions
#input: a directed graph, two labels, and a partition
#output: a bipartite graph whose verties are the elments of the original directed graph. the partition in the graph is between vertices labeled with letter 1 and letter 2. there is an edge between vertex u labeled with letter1 and vertex v labeled with letter2 if there is a directed path from u to v
def one_sided_comparability_bigraph(digraph, letter1, letter2, partition):
    adjacency = {}
    for vertex in digraph.vertices():
        adjacency[vertex]=[]
        if partition[vertex] == letter1:
            reachable = list(digraph.breadth_first_search(vertex))
            neighbors = []
            for neighbor in reachable:
                if partition[neighbor]==letter2:
                    neighbors.append(neighbor)
            adjacency[vertex]=neighbors
    return DiGraph(adjacency)

#input: a graph, a partition, and the partition labels
#output: true if the graph has a perfect matching and false otherwise
def check_matching(graph, partition, letter1, letter2):
    copy = DiGraph(graph)
    size = len(graph.vertices())/2
    for vertex in copy.vertices():
        if partition[vertex]==letter1:
            copy.add_edge('s', vertex)
        if partition[vertex]==letter2:
            copy.add_edge(vertex, 't')
    max_flow = int(copy.flow('s','t'))
    if size == max_flow:
        return true

#input: a directed graph, a partition on the vertices, and the labels used in the partition
#output: "positive epistasis implied", "negative epistasis implied", "no epistasis implied"
def quick_epistasis_check(digraph, partition, letter1, letter2):
    bigraph1 = one_sided_comparability_bigraph(digraph, letter1, letter2, partition)
    if check_matching(bigraph1, partition, letter1, letter2):
        return 1
    bigraph2 = one_sided_comparability_bigraph(digraph, letter2, letter1, partition)
    if check_matching(bigraph2, partition, letter2, letter1):
        return -1
    else:
        return 0

#-------------------------------------------------------------------------------------------------------------------------------------------------#
# Our function for computing the proportion of partial orders on n elements which imply epistasis
# would not reccomend using with n>8
# or n odd

#input: an even number n
#output: a list containing each dictionary on the elements 0, ..., n-1 such that n/2 of the keys have the value 1 and the other n/2 have the value 0
def bipartitions(n):
    partitions_n = []
    for subset in Subsets(n, n/2):
        partition = {}
        for x in range (0, n):
            if x+1 in subset:
                partition[x]=1
            else: partition[x]=0
        partitions_n.append(partition)
    return partitions_n

#input: an even number. quick for n = 2, 4, 6, takes ~2 hours for n=8; with good luck, might finish n=10 this century
#output: a tuple containining the number of posets
def prop_implying_epistasis(n):
    count_pos = 0
    count_neg = 0
    count_none = 0
    posets_n = Posets(n)
    offset = factorial(n/2)^2
    for poset in posets_n:
        hasse = poset.hasse_diagram()
        auto = hasse.automorphism_group().cardinality()
        for bipartition in bipartitions(n):
            if quick_epistasis_check(hasse, bipartition, 0, 1)==1:
                count_pos = count_pos + 1/auto
            elif quick_epistasis_check(hasse, bipartition, 0, 1)==-1:
                count_neg = count_neg+1/auto
            else: count_none = count_none+1/auto
    return (count_pos*offset, count_neg*offset, count_none*offset)