Initial commit

.gitignore

+log.*
+/data/
+/pictures
+*.png
+*.dot
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+pip-wheel-metadata/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+.python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/

README.md

+# Build the population tree from D statistics.
+
+## Prepare
+
+python==3.6.8
+
+Install graphviz using apt-get or conda.
+```
+conda install graphviz
+```
+Install python packages.
+```
+pip install -r requirements.txt
+``` 
+
+## Run
+
+Configure data file, save path, threshold, and populations.
+
+```
+python main.py
+```

data_structures.py

+from itertools import combinations
+from collections import defaultdict
+import numpy as np
+import copy
+
+class Graph: 
+    """
+    get graph from a list of edges.
+    note that if a fully-connected graph is NOT a DAG,
+    then the smallest graph (a triangle) in it must also be a DCG.
+    """
+    def __init__(self, vertices): 
+        #self.graph = defaultdict(list) 
+        self.graph = {} 
+        self.start_nodes = set() 
+        #-------- for equal chain --------#
+        self.eq_map_in = {}
+        self.eq_map_out = {}
+        #--------------------------------#        
+        #-------- for SCC --------#
+        self.V = vertices
+        self.Time = 0
+        #-------------------------#
+  
+    def addEdge_eq_as_new_nodes(self,u, v, equal=False):
+        if equal:
+            if u not in self.graph:
+                self.graph.update({u: {'in': set(), 'out': set(), 'eq': set()}})
+            if v not in self.graph:
+                self.graph.update({v: {'in': set(), 'out':set(), 'eq': set()}})
+            self.graph[u]['eq'].add(v)
+            self.graph[v]['eq'].add(u)
+        else:
+            if u not in self.graph:
+                self.graph.update({u: {'in': set(), 'out': set(), 'eq': set()}})
+                self.start_nodes.add(u)
+            if v not in self.graph:
+                self.graph.update({v: {'in': set(), 'out':set(), 'eq': set()}})
+            self.graph[u]['out'].add(v)
+            self.graph[v]['in'].add(u)
+            self.start_nodes.discard(v)
+
+    def get_adj(self):
+        pass
+
+    def addEdge(self,u, v, equal=False, value=None):
+        u = self.eq_map_out[u] if u in self.eq_map_out else u
+        v = self.eq_map_out[v] if v in self.eq_map_out else v
+        if u == v: return
+        if equal:
+            # things must to know about equal:
+            # 1. equal is NOT transitive, this is caused by the unreliable threshold.
+            # 2. equal is symmetric.
+            # 3. for a fully-connected graph, there is at most one edge between two edges,
+            #    which results in the egde-graph a 2-hop DAG.
+            #    proof: for any two nonadjacent edges AB CD, we can connect AC BD, or AD BC two make them connected.
+            if u not in self.graph:
+                self.graph.update({u: {'in': set(), 'out': set(), 'eq': set()}})
+                self.start_nodes.add(u)
+            if v not in self.graph:
+                self.graph.update({v: {'in': set(), 'out':set(), 'eq': set()}})
+                self.start_nodes.add(v)
+            # self.graph[u]['eq'].add(v)
+            # self.graph[v]['eq'].add(u)
+
+            #------------- eq as bi-directed edge ----------------#
+            self.graph[u]['in'].add(v) 
+            self.graph[u]['out'].add(v)
+            self.graph[v]['in'].add(u)
+            self.graph[v]['out'].add(u)
+
+            if value < 0:
+                self.graph[u]['eq'].add((v, -1*value))
+            else:
+                self.graph[v]['eq'].add((u, value))
+            #-----------------------------------------------------#
+            #------------- eq with ['in'] pointing to ----------------#
+            #self.addEqualWithInPointingTo_eq(u, v)
+            #---------------------------------------------------------#
+            #-------------- add equal nodes with an equal chain ---------------#
+            #self.addEqualWithEqualChain(u, v)
+            #------------------------------------------------------------------#
+
+        else:
+            if u not in self.graph:
+                self.graph.update({u: {'in': set(), 'out': set(), 'eq': set()}})
+                self.start_nodes.add(u)
+            if v not in self.graph:
+                self.graph.update({v: {'in': set(), 'out':set(), 'eq': set()}})
+            self.graph[u]['out'].add(v)
+            self.graph[v]['in'].add(u)
+            self.start_nodes.discard(v)
+        
+
+            #------------- eq with ['in'] pointing to ----------------#
+            #self.addEqualWithInPointingTo_in(u, v)
+            #---------------------------------------------------------#
+            #------------- eq ----------------#
+            #----------------------------------#
+
+    def addEqualWithInPointingTo_eq(self, u, v):
+        self.graph[u]['in'].discard(v)
+        self.graph[u]['out'].discard(v)
+        self.graph[v]['in'].discard(u)
+        self.graph[v]['out'].discard(u)
+
+        for u_in in self.graph[u]['in']:
+            if u_in not in self.graph[v]['eq'] and u_in not in self.graph[v]['out']:
+                self.graph[u_in]['out'].add(v)
+                self.graph[v]['in'].add(u_in)
+                self.start_nodes.discard(v)
+        for v_in in self.graph[v]['in']:
+            if v_in not in self.graph[u]['eq'] and v_in not in self.graph[u]['out']:
+                self.graph[v_in]['out'].add(u)
+                self.graph[u]['in'].add(v_in)
+                self.start_nodes.discard(u)
+
+    def addEqualWithInPointingTo_in(self, u, v):
+        self.graph[u]['in'].discard(v)
+        self.graph[v]['out'].discard(u)
+        for v_eq in self.graph[v]['eq']:
+            if v_eq not in self.graph[u]['eq']:
+                self.graph[u]['out'].add(v_eq)
+                self.graph[v_eq]['in'].add(u)
+                self.start_nodes.discard(v_eq)
+
+    def addEqualWithEqualChain(self, u, v):
+        # max node -> min node
+        min_node = min(u, v)
+        max_node = max(u, v)
+        if min_node not in self.graph:
+            if len(self.graph) == 0:
+                self.start_nodes.add(min_node)
+            self.graph.update({min_node: {'in': set(), 'out': set(), 'eq': set()}})
+        self.graph[min_node]['eq'].add(max_node)
+        # replace the start node
+        if max_node in self.start_nodes:
+            self.start_nodes.remove(max_node)
+            self.start_nodes.add(min_node)
+        # the in-out edges and eq nodes of max node are added to min node
+        if max_node in self.graph:
+            self.graph[min_node]['in'] |= self.graph[max_node]['in']
+            self.graph[min_node]['out'] |= self.graph[max_node]['out']
+            self.graph[min_node]['eq'] |= self.graph[max_node]['eq']
+            # can be deleted if eq node inserted first
+            for node in self.graph[max_node]['in']:
+                self.graph[node]['out'].add(min_node)
+                self.graph[node]['out'].discard(max_node)
+            for node in self.graph[max_node]['out']:
+                self.graph[node]['in'].add(min_node)
+                self.graph[node]['in'].discard(max_node)
+
+            del self.graph[max_node]
+        # max node -> min node
+        if min_node not in self.eq_map_in:
+            self.eq_map_in.update({min_node: set()})
+        self.eq_map_in[min_node].add(max_node)
+        self.eq_map_out.update({max_node: min_node})
+        # nodes to max node -> nodes to min node
+        if max_node in self.eq_map_in:
+            self.eq_map_in[min_node] |= self.eq_map_in[max_node]
+            for node in self.eq_map_in[max_node]:
+                self.eq_map_out[node] = min_node
+            del self.eq_map_in[max_node]
+    
+    def copy(self, vertices=None):
+        if vertices is None:
+            return copy.deepcopy(self)
+        else:
+            new_graph = Graph(max(vertices) + 1) 
+            for v in vertices:
+                new_graph.graph[v] = copy.deepcopy(self.graph[v]) 
+            return new_graph
+    
+    def isReachable(self, s, d):
+        if s not in self.graph or d not in self.graph:
+            return False
+        # Mark all the vertices as not visited
+        visited =[False]*(self.V)
+  
+        # Create a queue for BFS
+        queue=[]
+  
+        # Mark the source node as visited and enqueue it
+        queue.append(s)
+        visited[s] = True
+  
+        while queue:
+ 
+            #Dequeue a vertex from queue
+            n = queue.pop(0)
+             
+            # If this adjacent node is the destination node,
+            # then return true
+            if n == d:
+                   return True
+ 
+            #  Else, continue to do BFS
+            for i in self.graph[n]['out']:
+                if i in self.graph and visited[i] == False:
+                    queue.append(i)
+                    visited[i] = True
+         # If BFS is complete without visited d
+        return False
+
+    def isReachableWithoutEq(self, s, d, eq_edges):
+        path = [s]
+
+        # Create a queue for BFS
+        queue = []
+  
+        # Mark the source node as visited and enqueue it
+        queue.append(path)
+
+        self.visited_edges = 0
+  
+        while queue:
+ 
+            #Dequeue a vertex from queue
+            tmp_path = queue.pop(0)
+            n = tmp_path[-1]
+             
+            # If this adjacent node is the destination node,
+            # then return true
+            if n == d:
+                # print(len(tmp_path))
+                # print(tmp_path)
+                if not self.onlyEqEdges(tmp_path, eq_edges):
+                    return True
+ 
+            #  Else, continue to do BFS
+            for i in self.graph[n]['out']:
+                self.visited_edges += 1
+                if i in self.graph and i not in tmp_path:
+                    queue.append(tmp_path + [i])
+         # If BFS is complete without visited d
+        return False
+    
+    def haveEdge(self, n1, n2):
+        if n1 not in self.graph or n2 not in self.graph:
+            return False
+        if n2 in self.graph[n1]['in'] or n2 in self.graph[n1]['out']:
+            return True
+        else:
+            return False
+            
+
+    def onlyEqEdges(self, path, eq_edges):
+        edges = list(combinations(path, 2))
+        for n1, n2 in edges:
+            if self.haveEdge(n1, n2) and (n1, n2) not in eq_edges:
+                return False
+        return True
+
+
+        
+    def SCCUtil(self, u, low, disc, stackMember, st, scc, vertices=None, size=2):
+        # Initialize discovery time and low value
+        disc[u] = self.Time
+        low[u] = self.Time
+        self.Time += 1
+        stackMember[u] = True
+        st.append(u)
+ 
+        # Go through all vertices adjacent to this
+        for v in self.graph[u]['out']:
+ 
+            if vertices is not None:
+                if v not in vertices:
+                    continue
+            # If v is not visited yet, then recur for it
+            if disc[v] == -1:
+ 
+                self.SCCUtil(v, low, disc, stackMember, st, scc, vertices, size)
+ 
+                # Check if the subtree rooted with v has a connection to
+                # one of the ancestors of u
+                # Case 1 (per above discussion on Disc and Low value)
+                low[u] = min(low[u], low[v])
+ 
+            elif stackMember[v] == True:
+ 
+                '''Update low value of 'u' only if 'v' is still in stack
+                (i.e. it's a back edge, not cross edge).
+                Case 2 (per above discussion on Disc and Low value) '''
+                low[u] = min(low[u], disc[v])
+ 
+        # head node found, pop the stack and print an SCC
+        w = -1  # To store stack extracted vertices
+        if low[u] == disc[u]:
+            tmp = []
+            while w != u:
+                w = st.pop()
+                # print(w, end=" ")
+                tmp.append(w)
+                stackMember[w] = False
+            if len(tmp) > size:
+                scc.append(tmp)
+            # print()
+ 
+    # The function to do DFS traversal.
+    # It uses recursive SCCUtil()
+ 
+    def SCC(self, vertices=None, size=2):
+        # Mark all the vertices as not visited
+        # and Initialize parent and visited,
+        # and ap(articulation point) arrays
+        self.Time = 0
+        disc = [-1] * (self.V)
+        low = [-1] * (self.V)
+        stackMember = [False] * (self.V)
+        st = []
+        scc = []
+        # Call the recursive helper function
+        # to find articulation points
+        # in DFS tree rooted with vertex 'i'
+        if vertices is None:
+            for i in range(self.V):
+                if i not in self.graph:
+                    continue
+                if disc[i] == -1:
+                    self.SCCUtil(i, low, disc, stackMember, st, scc, size=size)
+        else:
+            for i in vertices:
+                if disc[i] == -1:
+                    self.SCCUtil(i, low, disc, stackMember, st, scc, vertices, size=size)
+        return scc
+
+#    def topologicalSortUtil(self,v,visited,stack): 
+#  
+#        visited[v] = True
+#  
+#        for i in self.graph[v]: 
+#            if visited[i] == False: 
+#                self.topologicalSortUtil(i,visited,stack) 
+#  
+#        stack.insert(0,v) 
+#  
+#    def topologicalSort(self): 
+#        visited = [False]*self.V 
+#        stack =[] 
+#  
+#        for i in range(self.V): 
+#            if visited[i] == False: 
+#                self.topologicalSortUtil(i,visited,stack) 
+
+
+class UnionFindSet():
+    def __init__(self, data_list):
+        self.father_dict = {}
+        self.size_dict = {}
+        self.different_sets = {}
+        self.data_list = data_list
+        nodes = np.array(data_list).flatten().tolist()
+        for node in nodes:
+            self.father_dict[node] = node
+            self.size_dict[node] = 1
+            self.different_sets.update({node: set([node])})
+ 
+    def find(self, node):
+        father = self.father_dict[node]
+        if(node != father):
+            if father != self.father_dict[father]:
+                self.size_dict[father] -= 1
+            father = self.find(father)
+        self.father_dict[node] = father
+        return father
+ 
+    def is_same_set(self, node_a, node_b):
+        return self.find(node_a) == self.find(node_b)
+ 
+    def union(self, node_a, node_b):
+        if node_a is None or node_b is None:
+            return
+        a_root = self.find(node_a)
+        b_root = self.find(node_b)
+ 
+        if(a_root != b_root):
+            a_set_size = self.size_dict[a_root]
+            b_set_size = self.size_dict[b_root]
+            if(a_set_size >= b_set_size):
+                self.father_dict[b_root] = a_root
+                self.size_dict[a_root] = a_set_size + b_set_size
+                self.different_sets[a_root] |= self.different_sets[b_root]
+                del self.different_sets[b_root]
+            else:
+                self.father_dict[a_root] = b_root
+                self.size_dict[b_root] = a_set_size + b_set_size
+                self.different_sets[b_root] |= self.different_sets[a_root]
+                del self.different_sets[a_root]
+        
+    def init_union(self):
+        for i in range(len(self.data_list)):
+            node1, node2 = self.data_list[i][0], self.data_list[i][1]
+            self.union(node1, node2)
+
+
+if __name__ == '__main__':
+    g = Graph(6)
+    g.addEdge(0, 1)
+    g.addEdge(2, 3)
+    g.addEdge(2, 4)
+    g.addEdge(3, 4)
+    g.addEdge(0, 5)
+    print(g.graph)
+    print(g.isReachableWithoutEq(2,4,[]))
\ No newline at end of file

geneflow_P.py

+
+# D1 = d of D1(S2, Tar; S1, Og)[4]
+# Z1 = Z of D1(S2, Tar; S1, Og)[5]
+# D2 = d of D2(S2, S1'; S1, Og)[4]
+# Z2 = Z of D2(S2, S1'; S1, Og)[5]
+# Mix = mixture from S1 to Tar
+# M_PD = the Probability distribution of the mixture
+
+def Probability_distribution(D1,D2,Z1,Z2):
+    import matplotlib.pyplot as plt
+    import numpy as np
+    #x=D2(S1'):
+    u2 = D2
+    se2 = D2/Z2 
+    #y=D1(Tar):
+    u1 = D1
+    se1 = D1/Z1
+    Mix=D1/D2
+    import scipy.integrate
+    from numpy import exp
+    from scipy.stats import norm
+    if Mix >=0 and Mix <=1:
+        f = lambda y,x: (norm.pdf(x, loc=u2, scale=se2))*(norm.pdf(y, loc=u1, scale=se1))
+        g = 0
+        h = 1
+        i1_1 = scipy.integrate.dblquad(f, 0, 1, g, h)[0]
+        g = 0
+        h = 1
+        i1_2 = scipy.integrate.dblquad(f, -1, 0, g, h)[0]
+        g = -1
+        h = 0
+        i1_3 = scipy.integrate.dblquad(f, -1, 0, g, h)[0]
+        g = -1
+        h = 0
+        i1_4 = scipy.integrate.dblquad(f, 0, 1, g, h)[0]
+        #i1_1-i1_4:1st-4th quadrant
+        g = -1
+        h = 1
+        i1_5 = scipy.integrate.dblquad(f, -1, 1, g, h)[0]
+        #i1_5:x[-1:1],y[-1,1]
+        g = lambda x: x
+        h = 0
+        i2 = scipy.integrate.dblquad(f, -1, 0, g, h)[0]
+        g = 0
+        h = lambda x: x
+        i3 = scipy.integrate.dblquad(f, 0, 1, g, h)[0]
+        i=i2+i3
+        if i1_5 <= 0 or i1_5 <= i:
+            i1=i1_1 + i1_2 + i1_3 + i1_4
+        if i1_5 > 0 and i1_5 >= i:
+            i1=i1_5
+        M_PD=i/i1
+        if M_PD > 1:
+            M_PD=1
+    else:
+       M_PD="Warning: No calculation required!" 
+    return(Mix,M_PD)

main.py

+from treebuilder import TreeBuilder
+from itertools import combinations
+from utils import read_data, equal_vertices_judgement, simple_cycles
+import sys
+
+# integer programming
+import networkx as nx
+from cvxopt import matrix
+from cvxopt.glpk import ilp
+import numpy as np
+import os
+
+DEBUG = False 
+args = {
+    'change_greenstar': True,
+    'threshold1': 3, # red star threshold
+    'threshold2': 0, # gene flow display threshold
+    'threshold3': 3, # gene flow merge threshold 
+    'bidirectional': True, # calculate the reversed gene flow
+    'calculate_pd': False, # calculate the gene flow confidence
+    'data_path': 'data/test.log',
+    'save_path': 'pictures/tree2.png',
+    'draw_gene_flow': True,
+}
+
+populations = ['Kostenki14',
+               'Tianyuan',
+               'AR19K',
+               'AR14K',
+               'Zongri5.1k',
+               'Boshan',
+               'Qihe3',
+               'Donghulin',
+               ]
+
+#args = {
+#    'change_greenstar': True,
+#    'threshold1': 3, # red star threshold
+#    'threshold2': 0, # gene flow display threshold
+#    'threshold3': 3, # gene flow merge threshold 
+#    'bidirectional': True, # calculate the reversed gene flow
+#    'calculate_pd': False, # calculate the gene flow confidence
+#    'data_path': 'data/sub.f4.GX.log',
+#    'save_path': 'pictures/tree.png',
+#    'draw_gene_flow': True,
+#}
+#
+#populations = ['Kostenki14',
+#               'Tianyuan',
+#               'Boshan',
+#               'Qihe',
+#               'Longlin',
+#               'Dushan',
+#               'Baojianshan',
+#               'La368'
+#               ]
+
+
+#data_path = 'data/f4_20230111.log'
+#populations = [
+#    'Kostenki14',
+#    'Tianyuan',
+#    'Yana',
+#    'Liangdao2',
+#    'Qihe',
+#    'Yumin',
+#    'Boshan',
+#    'Xiaogao',
+#    'Kolyma',
+#    'Shamanka_EN',
+#    'Longlin',
+#    'Dushan',
+#    'Baojianshan',
+#    'La368'
+#]
+#
+#data_path = 'data/f4.YN.log'
+#populations = [
+#    'Mbuti',
+#    'Tianyuan',
+#    'La368',
+#    'Longlin',
+#    'Xingyi_EN',
+#    'Xingyi_LN',
+#    'Qihe3',
+#    'Boshan',
+#    'Baiyangcun'
+#]
+
+change_greenstar = args.get('change_greenstar', True) 
+threshold1 = args.get('threshold1', 3)
+threshold2 = args.get('threshold2', 0)
+threshold3 = args.get('threshold3', 3)
+data_path = args.get('data_path')
+save_path = args.get('save_path')
+print_gene_flow = args.get('draw_gene_flow', True)
+
+threshold = threshold1
+gene_flow_threshold = threshold3
+
+
+os.makedirs('/'.join(save_path.split('/')[:-1]), exist_ok=True)
+relations = list(combinations(sorted(populations), 2))
+
+#--------------- check scc with threshold = 0 ------------------#
+print()
+print('Analyzing the minimum Z...')
+# get all edges
+bottom_threshold = 0
+D_statistics, adjusted_D_values, triangles, redstar_indices, min_redstar, s1_record \
+    = read_data(data_path, bottom_threshold, populations)
+builder = TreeBuilder(populations, triangles, s1_record)
+builder.D2graph(D_statistics)
+g = builder.g
+scc0 = g.SCC()
+num_scc0 = len(scc0)
+if DEBUG:
+    print('num_scc0', num_scc0)
+    print('size of s0:', [len(s) for s in scc0])
+
+unchangeable_edges = []
+for s in scc0:
+    for idx in s:
+        # the inversed edge is unchangeable
+        for node_idx in g.graph[idx]['in']:
+            if node_idx in s:
+                unchangeable_edges.append((idx, node_idx))
+for e in unchangeable_edges:
+    e1 = set(builder.relations[e[0]])
+    e2 = set(builder.relations[e[1]])
+    common_pop = (e1 & e2).pop()
+    key = tuple(sorted(list(e1 | e2)))
+    value = triangles[key]['D'][key.index(common_pop)]
+    bottom_threshold = max(bottom_threshold, abs(value))
+    if key.index(common_pop) != redstar_indices[key]:
+        print('Warning: changing the value which is not the red star!')
+        if not change_greenstar:
+            raise ValueError('Error: this graph cannot be plotted without changing the green stars, please set green=True!')
+        
+        for i in range(key.index(common_pop), key.index(common_pop)+3):
+            print(key[(i+1)%3], key[(i+2)%3], key[i%3], triangles[key]['D'][i%3])
+if DEBUG:
+    print('bottom_threshold', bottom_threshold)
+    print('unchangeable_edges', len(unchangeable_edges))
+print()
+print('Minimum Z:', bottom_threshold)
+
+if threshold < bottom_threshold:
+    ValueError('the given Z is smaller than the minimum Z!')
+
+#---------------------------------------------------------------#
+D_statistics, adjusted_D_values, triangles, redstar_indices, min_redstar, s1_record \
+    = read_data(data_path, threshold, populations)
+if DEBUG:
+    if min_redstar:
+        print('min red star:', min(min_redstar))
+    print(adjusted_D_values)
+builder = TreeBuilder(populations, triangles, s1_record)
+builder.D2graph(D_statistics)
+g = builder.g
+
+
+print()
+print('Adjusting values...')
+
+eq_scc0 = {}
+for scc_idx, s in enumerate(scc0):
+    # if len(s) <= 2: continue
+    if DEBUG:
+        print('-------------- scc ----------------')
+        print(s)
+    represent = s[-1]
+    for idx in s[:-1]:
+        for node in g.graph[idx]['in']:
+            if node in s: continue
+            g.graph[represent]['in'].add(node)
+            g.graph[node]['out'].add(represent)
+            g.graph[node]['out'].remove(idx)
+            # note that, the eq set of the node may contain duplicated nodes with different values!!!!!!!
+            # tmp_node_record is used to duplicate them.
+#            tmp_node_record = []
+            new_eq_set = set()
+            for eq_node, value in g.graph[node]['eq']:
+#                if eq_node in tmp_node_record: continue
+                if eq_node == idx:
+                    eq_node = represent
+                new_eq_set.add((eq_node, value))
+#                tmp_node_record.append(eq_node) 
+            g.graph[node]['eq'] = new_eq_set
+        for node in g.graph[idx]['out']:
+            if node in s: continue
+            g.graph[represent]['out'].add(node)
+            g.graph[node]['in'].add(represent)
+            g.graph[node]['in'].remove(idx)
+        for node, value in g.graph[idx]['eq']:
+            if node in s: continue
+            g.graph[represent]['eq'].add((node, value))
+        del g.graph[idx]
+        g.graph[represent]['in'].discard(idx)
+        g.graph[represent]['out'].discard(idx)
+        g.start_nodes.discard(idx) 
+    # note that, the eq set of the node may contain duplicated nodes with different values!!!!!!!
+    # tmp_node_record is used to duplicate them.
+#    tmp_node_record = []
+    new_eq_set = set()
+    for node, value in g.graph[represent]['eq']:
+        if node in s: continue #or node in tmp_node_record: continue
+        new_eq_set.add((node, value))
+#        tmp_node_record.append(node)
+#    del tmp_node_record
+    g.graph[represent]['eq'] = new_eq_set
+    if DEBUG:
+        print(g.graph[represent]['in'])
+        print(g.graph[represent]['out'])
+        print(g.graph[represent]['eq'])
+    eq_scc0[represent] = s[:-1]
+
+if DEBUG:
+    print('start_nodes', g.start_nodes)
+scc = g.SCC()
+num_scc = len(scc)
+if DEBUG:
+    print('num_scc', num_scc)
+    print('size of s:', [len(s) for s in scc])
+
+num_adjusted = []
+total_edges_adjusted = []
+eq_edges_map = {}
+for scc_idx, s in enumerate(scc):
+    # if scc_idx != 1: continue
+    print()
+    print('Adjusting %d/%d scc...' % (scc_idx + 1, num_scc))
+    if DEBUG:
+        print(s)
+    edges_to_adjust = []
+    edges_to_adjust_with_value = []
+    eq_edges = []
+    eq_nodes = set()
+    nodes_to_adjust = set()
+    edges_to_adjust_weights = []
+    for idx in s:
+        tmp_in = set()
+        tmp_out = set()
+        tmp_eq = set() 
+        for n in g.graph[idx]['in']:
+            if n in s:
+                tmp_in.add(n)
+        for n in g.graph[idx]['out']:
+            if n in s:
+                tmp_out.add(n)
+        for n, v in g.graph[idx]['eq']:
+            if n in s:
+                tmp_eq.add(n)
+#        print('----------------')
+#        print(tmp_out)
+#        print(tmp_eq)
+        if tmp_out != tmp_in:
+                eq_nodes.add(idx)
+        for node_idx, value in g.graph[idx]['eq']:
+            if node_idx in s:
+                eq_edges.append((idx, node_idx))
+                eq_edges.append((node_idx, idx))
+            if node_idx in s and (idx, node_idx) not in unchangeable_edges:
+                # print('eq:', idx, node_idx)
+                if (idx, node_idx) not in edges_to_adjust:
+                    edges_to_adjust_weights.append(1)
+                    edges_to_adjust.append((idx, node_idx))
+                    nodes_to_adjust.add(idx)
+                    nodes_to_adjust.add(node_idx)
+                    edges_to_adjust_with_value.append((idx, node_idx, value))
+                    eq_edges_map[(idx, node_idx)] = [(idx, node_idx)]
+                else:
+                    edges_to_adjust_weights[edges_to_adjust.index((idx, node_idx))] += 1
+                    eq_edges_map[(idx, node_idx)].append((idx, node_idx))
+    if DEBUG: 
+        print('calculating adjusted...')
+    if DEBUG:
+        print(len(s), len(edges_to_adjust_with_value))
+        print(len(eq_edges))
+    adjusted = []
+    for e_idx, (node_o, node_i, value) in enumerate(edges_to_adjust_with_value):
+        g_tmp = g.copy(s)
+        g_tmp.graph[node_o]['out'].discard(node_i)
+        g_tmp.graph[node_i]['out'].discard(node_o)
+        g_tmp.graph[node_i]['in'].discard(node_o)
+        g_tmp.graph[node_o]['in'].discard(node_i)
+        g_tmp.graph[node_o]['eq'].discard((node_i, value))
+        for e_idx2, (node_o2, node_i2, value2) in enumerate(edges_to_adjust_with_value):
+            if e_idx2 == e_idx: continue 
+            g_tmp.graph[node_o2]['out'].discard(node_i2)
+            g_tmp.graph[node_i2]['in'].discard(node_o2)
+            g_tmp.graph[node_o2]['eq'].discard((node_i2, value2))
+        if g_tmp.isReachableWithoutEq(node_i, node_o, eq_edges):
+            adjusted.append((node_o, node_i, value))
+    if DEBUG:
+        print('adjusted', len(adjusted))
+        print(adjusted)
+    num_adjusted.append(len(adjusted))
+    for e_idx, (node_o, node_i, value) in enumerate(adjusted):
+        g.graph[node_o]['out'].remove(node_i)
+        g.graph[node_i]['in'].remove(node_o)
+        edges_to_adjust.remove((node_o, node_i))
+        edges_to_adjust_with_value.remove((node_o, node_i, value))
+
+        # count adjusted edges
+        for e in eq_edges_map[(node_o, node_i)]:
+            n1, n2 = e
+            r1 = set(builder.relations[n1])
+            r2 = set(builder.relations[n2])
+            if not (r1 & r2): continue
+            common_pop = (r1 & r2).pop()
+            key = tuple(sorted(list(r1 | r2)))
+            k_idx = key.index(common_pop)
+            value = triangles[key]['D'][k_idx]
+            total_edges_adjusted.append((key[(k_idx+1)%3], key[(k_idx+2)%3], key[k_idx], value))
+    # print('adjusted', adjusted)
+    if DEBUG:
+        print(equal_vertices_judgement(g, s))
+        print(g.SCC(s))
+        print(s)
+
+    G = nx.DiGraph()
+    for node_idx in s:
+        G.add_node(node_idx)
+        for node_idx2 in g.graph[node_idx]['in']:
+            if node_idx2 in s:
+                G.add_edge(node_idx2, node_idx)
+    if DEBUG:
+        print('here')
+        print('num of nodes:', len(G.nodes))
+        print('num of edges:', len(G.edges))
+        print('num of adjust nodes:', len(nodes_to_adjust))
+        print('num of adjust edges:', len(edges_to_adjust))
+        print(edges_to_adjust)
+        # print('eq_edges', eq_edges)
+        print('eq_edges', len(eq_edges))
+        print('eq_nodes', len(eq_nodes))
+    cycles = list(simple_cycles(G, eq_edges, edges_to_adjust))
+
+    num_edges = len(edges_to_adjust)
+    num_cycles = len(cycles)
+    if DEBUG:
+        print(num_edges)
+        print(num_cycles)
+        print(edges_to_adjust_weights)
+        print(edges_to_adjust)
+    if num_cycles == 0:
+        continue
+
+    c = matrix(np.array(edges_to_adjust_weights))
+    G2 = matrix(-1*np.array(cycles, dtype=float))
+    h = matrix(-1 * np.ones(num_cycles))
+    I = set()
+    B = set(range(num_edges))
+
+    print()
+    print('------------------------------- Solving ILP -------------------------------')
+    (status,x) = ilp(c,G2,h,B=B)
+    print('---------------------------------------------------------------------------')
+    if DEBUG:
+        print(x)
+    x = np.array(x)
+    adjusted_edge_idx = np.where(x==1)[0]
+    num_adjusted[-1] += len(adjusted_edge_idx)
+    for e_idx in adjusted_edge_idx:
+        node_o, node_i = edges_to_adjust[e_idx]
+        g.graph[node_o]['out'].remove(node_i)
+        g.graph[node_i]['in'].remove(node_o)
+
+        # count adjusted edges
+        for e in eq_edges_map[(node_o, node_i)]:
+            n1, n2 = e
+            r1 = set(builder.relations[n1])
+            r2 = set(builder.relations[n2])
+            if not (r1 & r2): continue
+            common_pop = (r1 & r2).pop()
+            key = tuple(sorted(list(r1 | r2)))
+            k_idx = key.index(common_pop)
+            value = triangles[key]['D'][k_idx]
+            total_edges_adjusted.append((key[(k_idx+1)%3], key[(k_idx+2)%3], key[k_idx], value))
+
+if DEBUG:
+    print('num adjusted', num_adjusted)
+print()
+print('Adjusted values', total_edges_adjusted)
+print('Total adjusted:', len(total_edges_adjusted))
+
+if DEBUG:
+    print('------------------------')
+scc = g.SCC(size=1)
+num_scc = len(scc)
+eq_scc = eq_scc0 
+for scc_idx, s in enumerate(scc):
+    # if len(s) <= 2: continue
+    if not equal_vertices_judgement(g, s):
+        print('This scc are not eq:', scc_idx)
+        sys.exit(0)
+    if DEBUG:
+        print('--------- scc ----------')
+    represent = s[-1]
+    for idx in s:
+        if idx in eq_scc:
+            represent = idx
+            break
+    s.remove(represent)
+    for idx in s:
+        if DEBUG:
+            print(idx, builder.relations[idx])
+        for node in g.graph[idx]['in']:
+            g.graph[represent]['in'].add(node)
+            g.graph[node]['out'].add(represent)
+            g.graph[node]['out'].remove(idx)
+        for node in g.graph[idx]['out']:
+            g.graph[represent]['out'].add(node)
+            g.graph[node]['in'].add(represent)
+            g.graph[node]['in'].remove(idx)
+        del g.graph[idx]
+        g.start_nodes.discard(idx)
+    g.graph[represent]['in'].remove(represent)
+    g.graph[represent]['out'].remove(represent)
+    if represent in eq_scc:
+        eq_scc[represent].extend(s)
+    else:
+        eq_scc[represent] = s
+    if DEBUG:
+        print(g.graph[represent])
+    
+#    print(edges_to_adjust)
+if DEBUG:
+    print('------------------------')
+    for node in g.graph:
+        if len(g.graph[node]['in']) == 0:
+            print(node)
+    print(g.start_nodes)
+
+print('eq_scc', eq_scc)
+print()
+print('Building the tree...')
+builder.build_tree(eq_scc=eq_scc, D=None, **args)
+print()
+print('Plotting...')
+builder.export_graph(save_path, print_gene_flow=True)
+print()
+print('Saved to %s' % save_path)
+
+print()
+print('Done.')
\ No newline at end of file

preprocess.py

+DEBUG = False
+triangle_templates = [
+#    [ 0, 0, 0],
+#    [ 1,-1, 1],
+#    [ 1,-1, 0],
+#    [ 1,-1,-1],
+#    [-1,-1, 1],
+#    [-1, 0, 1],
+#    [-1, 1, 1],
+#    [ 1, 1,-1],
+#    [ 0, 1,-1],
+#    [-1, 1,-1],
+#    [ 1, 0,-1],
+#    [-1, 1, 0],
+#    [ 0,-1, 1],
+
+#    [ 0, 0, 0],
+
+#    [ 10,-10, 10],
+#    [ 10,-10, 0],
+#    [ 10,-10,-10],
+#    [-10,-10, 10],
+#    [-10, 0,  10],
+#    [-10, 10, 10],
+#    [ 10, 10,-10],
+#    [ 0,  10,-10],
+#    [-10, 10,-10],
+
+    [ 10, 10,-10],
+    [ 10, 0 ,-10],
+    [ 10,-10,-10],
+    [-10, 10,-10],
+    [-10, 10, 0 ],
+    [-10, 10, 10],
+    [ 10,-10, 10],
+    [ 0 ,-10, 10],
+    [-10,-10, 10],
+]
+redstar_index = [2, 2, 2, 1, 1, 1, 0, 0, 0]
+
+def adjust_D_bak(triangle, threshold):
+    adjusted_data = []
+    to_adjust = []
+    for i, value in enumerate(triangle):
+        if abs(value) < threshold:
+            adjusted_data.append(0)
+            to_adjust.append((i, value))
+        elif value >= threshold:
+            adjusted_data.append(10)
+        elif value <= -threshold:
+            adjusted_data.append(-10)
+    to_adjust = sorted(to_adjust, key=lambda x: abs(x[1]), reverse=True)
+    adjusted_value = []
+    adjust_threshold = None 
+    while adjusted_data not in triangle_templates:
+        if len(to_adjust) == 0:
+            raise Exception('Cannot adjust triangle: %s' % triangle)
+        adjust_idx, adjust_threshold = to_adjust.pop(0)
+        adjusted_data[adjust_idx] = 10 if adjust_threshold > 0 else -10
+        adjusted_value.append((adjust_idx, adjust_threshold))
+    return adjusted_value, adjusted_data 
+
+def adjust_D(key, triangle, threshold, min_redstar, s1_record):
+    adjusted_data = []
+    to_adjust = []
+    for i, value in enumerate(triangle):
+#        if abs(value) < threshold:
+#            adjusted_data.append(0)
+#            to_adjust.append((i, value))
+        if value > 0:
+            adjusted_data.append(10)
+        elif value < 0:
+            adjusted_data.append(-10)
+        else:
+            raise ValueError('value is 0')
+    if adjusted_data not in triangle_templates:
+        raise Exception('Cannot adjust triangle: %s, %s' % (key, triangle))
+    index = triangle_templates.index(adjusted_data)
+    redstar_idx = redstar_index[index]
+
+
+    def add_s1_record(s2, tar, s1):
+        if (key[s2], key[tar]) not in s1_record:
+            s1_record[(key[s2], key[tar])] = []
+        s1_record[(key[s2], key[tar])].append(key[s1])
+    T = 0
+    A = 1
+    B = 2
+    if redstar_idx == 2:
+        add_s1_record(T, A, B)
+        add_s1_record(B, A, T)
+    elif redstar_idx == 1:
+        add_s1_record(T, B, A)
+        add_s1_record(A, B, T)
+    elif redstar_idx == 0:
+        add_s1_record(A, T, B)
+        add_s1_record(B, T, A)
+    if DEBUG:
+        print('------------------------------')
+        print(key)
+        print(triangle)
+
+    min_redstar.append(abs(triangle[redstar_idx]))
+    if abs(triangle[redstar_idx]) < threshold:
+        adjusted_data[redstar_idx] = 0
+    return None, adjusted_data, redstar_idx 
+
+def get_adjusted_D_statistics(D_statistics, adjusted_D_values, triangles, A, B, C, f4, value, threshold, min_redstar, redstar_indices, key_visited, s1_record):
+    key = tuple(sorted([A, B, C]))
+    if key in key_visited: return
+    pos = key.index(C)
+    sign = 1 if (key.index(A) - pos) % 3 < (key.index(B) - pos) % 3 else -1
+    value *= sign
+    f4 *= sign
+    # print(key, value)
+    if key not in triangles:
+        triangles[key] = {'D':[None, None, None], 'f4':[None, None, None]}
+    triangles[key]['f4'][pos] = f4
+    triangles[key]['D'][pos] = value
+
+    if None not in triangles[key]['D']:
+        key_visited.append(key)
+        adjusted_value, adjusted_data, redstar_index = adjust_D(key, triangles[key]['D'], threshold, min_redstar, s1_record)
+        redstar_indices[key] = redstar_index
+#        if len(adjusted_value) > 0:
+#            for idx, value in adjusted_value:
+#                adjusted_D_values.append((key[idx], key[(idx+1)%3], key[(idx+2)%3], value))
+        for i, value in enumerate(adjusted_data):
+            # value == 0: record the original number
+            if value == 0:
+                D_statistics.append((key[(i+1)%3], key[(i+2)%3], key[i], value, triangles[key]['D'][i]))
+            else:
+                D_statistics.append((key[(i+1)%3], key[(i+2)%3], key[i], value, triangles[key]['D'][i]))
+        print_key = key[:] + key[:1]
+        # print(print_key, triangles[key], adjusted_data)
+#        if key == ('Dushan', 'Longlin', 'Qihe'):
+#            print(D_statistics)
+        # del triangles[key]
\ No newline at end of file

propotion.py

+import matplotlib.pyplot as plt
+from sklearn.linear_model import LinearRegression
+import numpy as np
+
+
+
+populations = ['Kostenki14',
+               'Tianyuan',
+               'Boshan',
+               'Qihe',
+               'Longlin',
+               'Dushan',
+               'Baojianshan',
+               'La368']
+
+all_statistics = dict()
+all_populations = set()
+with open('data/UPA34_4.f4.AI.log') as f:
+    for l in f.readlines():
+        if not l.startswith('result'): continue
+        A, B, C = l.split()[1:4]
+        all_populations.add(A)
+        all_populations.add(B)
+        all_populations.add(C)
+        value = float(l.split()[5])
+        all_statistics[(A, B, C)] = value
+
+test = 'Dushan'
+R1 = 'Baojianshan'
+x = []
+y = [] 
+print(all_populations)
+for R2 in all_populations:
+    if R2 != 'Qihe': continue
+    for O1 in all_populations:
+        if O1 in populations: continue
+        if R2 == O1 or R1 == test or R2 == test or R1 == O1: continue
+        y.append(all_statistics[(test, R1, O1)])
+        x.append(all_statistics[(test, R2, O1)])
+    
+print(len(x))
+
+plt.scatter(x, y)
+plt.savefig('test.png')
+
+x = np.array(x).reshape(-1, 1)
+y = np.array(y).reshape(-1, 1)
+model = LinearRegression()
+model.fit(x, y)
+
+print("k", model.coef_)
+print("b", model.intercept_)
+print("R2", model.score(x, y))
\ No newline at end of file

requirements.txt

+anytree==2.8.0
+cvxopt==1.3.0
+graphviz==0.19.1
+matplotlib==3.3.4
+networkx==2.5.1
+numpy==1.19.5
+Pillow==8.4.0
+pydot==1.4.2
+pyparsing==3.0.9
+scikit-learn==0.24.2
+six==1.16.0

test_cases.py

+
+class TestCase(object):
+    def __init__(self, populations=['A', 'B', 'T'], D_key=[('A', 'B', 'T'), ('A', 'T', 'B'), ('B', 'T', 'A')], D_value=[0, 0, 0]):
+        self.populations = populations
+        self.D_key = D_key
+        self.D_value = D_value
+        self.D_statistics = []
+        for k, v in zip(self.D_key, self.D_value):
+            self.D_statistics.append((*k, v))
+
+D_key = [('A', 'B', 'T'), ('A', 'T', 'B'), ('B', 'T', 'A')]
+populations = ['A', 'B', 'T']
+three_population_cases = []
+three_population_cases.append(TestCase(['A', 'B', 'T'], D_key, [ 0, 0, 0]))
+three_population_cases.append(TestCase(['A', 'B', 'T'], D_key, [ 1, 1, 1]))
+three_population_cases.append(TestCase(['A', 'B', 'T'], D_key, [ 1, 1, 0]))
+three_population_cases.append(TestCase(['A', 'B', 'T'], D_key, [ 1, 1,-1]))
+three_population_cases.append(TestCase(['A', 'B', 'T'], D_key, [-1, 1, 1]))
+three_population_cases.append(TestCase(['A', 'B', 'T'], D_key, [-1, 0, 1]))
+three_population_cases.append(TestCase(['A', 'B', 'T'], D_key, [-1,-1, 1]))
+three_population_cases.append(TestCase(['A', 'B', 'T'], D_key, [ 1,-1,-1]))
+three_population_cases.append(TestCase(['A', 'B', 'T'], D_key, [ 0,-1,-1]))
+three_population_cases.append(TestCase(['A', 'B', 'T'], D_key, [-1,-1,-1]))
+three_population_cases.append(TestCase(['A', 'B', 'T'], D_key, [ 1, 0,-1]))
+three_population_cases.append(TestCase(['A', 'B', 'T'], D_key, [-1,-1, 0]))
+three_population_cases.append(TestCase(['A', 'B', 'T'], D_key, [ 0, 1, 1]))
\ No newline at end of file

treebuilder.py

+from itertools import combinations
+from data_structures import UnionFindSet, Graph 
+from anytree import Node, PreOrderIter
+from anytree.exporter import DotExporter
+from collections import defaultdict
+# from sklearn.linear_model import LinearRegression
+import numpy as np
+import sys
+from geneflow_P import Probability_distribution
+
+DEBUG = False
+class TreeBuilder(object):
+    def __init__(self, populations, triangles, s1_record):
+        self.triangles = triangles
+        self.populations = sorted(populations)
+        self.num_populations = len(self.populations)
+        self.relations = list(combinations(self.populations, 2))
+        self.s1_record = s1_record
+        self.root = None 
+        self.gene_flow = None 
+
+    def D2graph(self, D):
+        g = Graph(len(self.relations))
+        for i, rule in enumerate(D):
+            node1 = tuple(sorted([rule[0], rule[2]]))
+            node2 = tuple(sorted([rule[1], rule[2]]))
+            relation = rule[3]
+            if relation == 10:
+                g.addEdge(self.relations.index(node1), self.relations.index(node2))
+            elif relation == -10:
+                g.addEdge(self.relations.index(node2), self.relations.index(node1))
+            elif -10 < relation and relation < 10:
+                g.addEdge(self.relations.index(node1), self.relations.index(node2), equal=True, value=rule[4])
+            else:
+                raise ValueError('relation must be 1, -1 or 0, but got {}'.format(relation))
+        self.g = g
+        return g 
+    
+    def get_order_from_graph_back(self, D):
+        # can be speed up by combining with build_tree
+        def helper(queue, order):
+            if len(queue) == 0: return
+            node = queue.pop()
+            order.append([node] + list(g.graph[node]['eq']))
+            for out_node in g.graph[node]['out']:
+                g.graph[out_node]['in'].remove(node)
+                if len(g.graph[out_node]['in']) == 0:
+                    queue.add(out_node)
+            del g.graph[node]
+            helper(queue, order)
+    
+        g = self.D2graph(D)
+        queue = g.start_nodes
+        order = []
+        helper(queue, order)
+        if len(g.graph) > 0:
+            raise ValueError('Graph is not a DAG')
+        return order
+
+    def get_order_from_graph(self, eq_scc=None, D=None):
+        # can be speed up by combining with build_tree
+        def helper(queue, order):
+            if len(queue) == 0: return
+            nodes = queue.pop(0)
+            order.append(nodes)
+            next_nodes = []
+            for node in nodes:
+                for out_node in g.graph[node]['out']:
+                    g.graph[out_node]['in'].remove(node)
+                    tmp_eq = set()
+                    for node_eq, _ in g.graph[out_node]['eq']:
+                        tmp_eq.add(node_eq)
+                    if len(g.graph[out_node]['in']) == 0:
+                        next_nodes.append(out_node)
+                    elif (g.graph[out_node]['in'] | tmp_eq) == tmp_eq:
+                        next_nodes.append(out_node)
+                        for node_i in g.graph[out_node]['in']:
+                            g.graph[node_i]['out'].remove(out_node)
+                del g.graph[node]
+            if len(next_nodes) > 0:
+                queue.append(next_nodes)
+            helper(queue, order)
+        if self.g is None: 
+            self.D2graph(D)
+
+        g = self.g
+        queue = [list(g.start_nodes)]
+        order = []
+        helper(queue, order)
+        if len(g.graph) > 0:
+            print('graph', g.graph)
+            raise ValueError('Graph is not a DAG')
+        for o_id, o in enumerate(order):
+            for node in o:
+                if node in eq_scc: 
+                    order[o_id].extend(eq_scc[node])
+
+        if DEBUG:
+            print('order', order)
+            print([[self.relations[idx] for idx in o] for o in order])
+        return order
+
+    def get_order_from_degree(self, D):
+        count = [0] * len(self.relations)
+        for i, rule in enumerate(D):
+            node1 = tuple(sorted([rule[0], rule[2]]))
+            node2 = tuple(sorted([rule[1], rule[2]]))
+            relation = rule[3]
+            if relation == 1:
+                count[self.relations.index(node2)] += 1
+            elif relation == -1:
+                count[self.relations.index(node1)] += 1
+            elif relation == 0:
+                continue
+            else:
+                raise ValueError('relation must be 1, -1 or 0, but got {}'.format(relation))
+        reverse_index_count = [[] for _ in range(self.num_populations)]
+        for i, c in enumerate(count):
+            reverse_index_count[c].append(i)
+        return reverse_index_count
+    
+    def record_gene_flow(self, visited_leaves, l1, l2, order_id):
+        if visited_leaves[l1]['id'] > visited_leaves[l2]['id']:
+            if visited_leaves[l1]['last_gene_flow_id'] < order_id:
+                visited_leaves[l1]['last_gene_flow_id'] = order_id
+            visited_leaves[l1]['gene_flow'].append(l2)
+        elif visited_leaves[l1]['id'] < visited_leaves[l2]['id']:
+            if visited_leaves[l2]['last_gene_flow_id'] < order_id:
+                visited_leaves[l2]['last_gene_flow_id'] = order_id
+            visited_leaves[l2]['gene_flow'].append(l1)
+        else:
+            if visited_leaves[l1]['last_gene_flow_id'] < order_id:
+                visited_leaves[l1]['last_gene_flow_id'] = order_id
+            if visited_leaves[l2]['last_gene_flow_id'] < order_id:
+                visited_leaves[l2]['last_gene_flow_id'] = order_id
+            visited_leaves[l1]['gene_flow'].append(l2)
+            # visited_leaves[l2]['gene_flow'].append(l1)
+
+    def add_gene_flow(self, gene_flow, visited_leaves, l1, l2, order_id, gene_flow_propotion_population, threshold=0):
+        node_l1 = visited_leaves[l1]['node']
+        node_l2 = visited_leaves[l2]['node']
+        if node_l1 is None or node_l2 is None:
+            return
+        # print('gene flow in {} and {}'.format(l1, l2))
+        # print(gene_flow)
+        # print(node_l1.depth)
+        # print(node_l2.depth)
+        # print(visited_leaves[l1]['gene_flow'])
+#        print('---------- propotion ----------')
+        if node_l1.depth > node_l2.depth and visited_leaves[l2]['last_gene_flow_id'] < order_id:
+            for leaf in visited_leaves[l2]['gene_flow']:
+                if abs(self.get_value(l1, leaf, l2)) > threshold:
+                    gene_flow.append((l2, leaf))
+#                print(l2, leaf, l1)
+                gene_flow_propotion_population.append((l2, leaf, l1))
+            visited_leaves[l2]['gene_flow'] = []
+        elif node_l1.depth < node_l2.depth and visited_leaves[l1]['last_gene_flow_id'] < order_id:
+            for leaf in visited_leaves[l1]['gene_flow']:
+                if abs(self.get_value(l2, leaf, l1)) > threshold:
+                    gene_flow.append((l1, leaf))
+#                print(l1, leaf, l2)
+                gene_flow_propotion_population.append((l1, leaf, l2))
+            visited_leaves[l1]['gene_flow'] = []
+        elif node_l1.depth == node_l2.depth and visited_leaves[l1]['last_gene_flow_id'] < order_id and visited_leaves[l2]['last_gene_flow_id'] < order_id:
+            for leaf in visited_leaves[l1]['gene_flow']:
+                if abs(self.get_value(l2, leaf, l1)) > threshold:
+                    gene_flow.append((l1, leaf))
+#                print(l1, leaf, l2)
+                gene_flow_propotion_population.append((l1, leaf, l2))
+            for leaf in visited_leaves[l2]['gene_flow']:
+                if abs(self.get_value(l1, leaf, l2)) > threshold:
+                    gene_flow.append((l2, leaf))
+#                print(l2, leaf, l1)
+                gene_flow_propotion_population.append((l2, leaf, l1))
+            visited_leaves[l1]['gene_flow'] = []
+            visited_leaves[l2]['gene_flow'] = []
+
+    def get_value(self, A, B, C, typ='D'):
+        key = tuple(sorted([A, B, C]))
+        pos = key.index(C)
+        sign = 1 if (key.index(A) - pos) % 3 < (key.index(B) - pos) % 3 else -1
+        value = sign * self.triangles[key][typ][pos]
+        return value
+    
+    def merge_gene_flow(self, gene_flow, visited_leaves, threshold=100):
+        def merge(merged_gene_flow, parents, threshold):
+            for p, nodes in parents.items():
+                merge_flag = True
+                for node1, node2 in combinations(nodes, 2):
+                    # check D statistics
+                    A, B, C = old_l1, node1.name, node2.name
+                    key = tuple(sorted([A, B, C]))
+                    value = self.triangles[key]['D'][key.index(A)]
+                    if abs(value) < threshold:
+                        merge_flag = False
+
+                if merge_flag and len(nodes) > 1 and visited_leaves[old_l1]['node'].parent != p:
+                    if len(nodes) == len(p.children):
+                        new_node = p
+                    else:
+                        new_node = Node(int(p.name) - 0.5, parent=p)
+                        if DEBUG:
+                            print('%s --1--> %s'% (p.name, str(p.name - 0.5)))
+                    merged_gene_flow.append((old_l1, new_node.name))
+                    for n in nodes:
+                        n.parent = new_node
+                        if DEBUG:
+                            print('%s --2--> %s'% (new_node.name, n.name))
+                else:
+                    for n in nodes:
+                        merged_gene_flow.append((old_l1, n.name))
+
+        if DEBUG:
+            print('=============================================================')
+        merged_gene_flow = []
+        gene_flow = sorted(gene_flow)
+        if DEBUG:
+            print(gene_flow)
+        parents = defaultdict(set)
+        old_l1 = None
+        for i, (l1, l2) in enumerate(gene_flow):
+            if DEBUG:
+                print(l1, l2)
+            if l1 != old_l1:
+                merge(merged_gene_flow, parents, threshold)
+                parents = defaultdict(set)
+                node_l2 = visited_leaves[l2]['node']
+                parents[node_l2.parent].add(node_l2)
+                if DEBUG:
+                    print('parents', parents)
+                old_l1 = l1
+            else:
+                node_l2 = visited_leaves[l2]['node']
+                parents[node_l2.parent].add(node_l2)
+                if DEBUG:
+                    print('parents', parents)
+        merge(merged_gene_flow, parents, threshold)
+        print()
+        print('merged_gene_flow', merged_gene_flow)
+        return merged_gene_flow
+
+    def build_tree(self, eq_scc=None, D=None, **kwargs):
+        threshold2 = kwargs.get('threshold2', 0)
+        gene_flow_threshold = kwargs.get('gene_flow_threshold', 100)
+
+        order = self.get_order_from_graph(eq_scc, D)
+        if DEBUG:
+            print('order', order)
+        visited_leaves = {} # {leaf: {'node': Node class, 'i100d': order_id, 'gene_flow': gene_flow}}
+        glue_id = 0
+        # pic_id = 0
+        cyclic_flag = False 
+        gene_flow = [] 
+        gene_flow_propotion_population = []
+        tmp_gene_flow = []
+        for order_id, edges in enumerate(order):
+            # if there are abruptions in the order, then the tree is cyclic
+            if len(edges) != 0 and cyclic_flag:
+                raise ValueError('cyclic graph')
+            elif edges == []:
+                cyclic_flag = True 
+
+            nodes_to_be_clustered = []
+            for e in edges:
+                l1, l2 = self.relations[e]
+                nodes_to_be_clustered.append([l1, l2])
+#                if e == 10:
+#                   print('out of add gene flow', l1, l2)
+                if l1 in visited_leaves and l2 in visited_leaves:
+                    self.add_gene_flow(tmp_gene_flow, visited_leaves, l1, l2, order_id, gene_flow_propotion_population, threshold2)
+                if l1 not in visited_leaves:
+                    visited_leaves.update({l1: {'node': None, 'id': order_id, 'gene_flow': [], 'last_gene_flow_id': 0}})
+                if l2 not in visited_leaves:
+                    visited_leaves.update({l2: {'node': None, 'id': order_id, 'gene_flow': [], 'last_gene_flow_id': 0}})
+                # gene flow record
+                self.record_gene_flow(visited_leaves, l1, l2, order_id)
+            # cluster the nodes with unionfindset, the clustered nodes will have the same parent (i.e. n-ary tree)
+            ufs = UnionFindSet(nodes_to_be_clustered)
+            # print(nodes_to_be_clustered)
+            ufs.init_union()
+
+            for cluster in ufs.different_sets.values():
+                # print('cluster', cluster)
+                glue_node = Node(glue_id)
+                glue_id += 1
+                for leaf in cluster:
+                    if visited_leaves[leaf]['node'] is None:
+                        new_node = Node(leaf, parent=glue_node)
+                        visited_leaves[leaf]['node'] = new_node
+                    else:
+                        tmp_node = visited_leaves[leaf]['node']
+                        if tmp_node.root != glue_node:
+                            tmp_node.root.parent = glue_node
+                if len(glue_node.children) <= 1:
+                    glue_node.children[0].parent = None
+                    glue_id -= 1
+
+            # root = visited_leaves[populations[0]]['node'].root
+            # RenderTreeGraph(root).to_picture("tree_%d.png" % pic_id)
+            # pic_id += 1
+        print()
+        print('all gene flows', tmp_gene_flow)
+        if len(tmp_gene_flow) > 0:
+            if DEBUG:
+                print('tmp_gene_flow', tmp_gene_flow)
+            tmp_gene_flow = self.calculate_propotion(tmp_gene_flow, order, visited_leaves, **kwargs)
+            merged_gene_flow = self.merge_gene_flow(tmp_gene_flow, visited_leaves, gene_flow_threshold)
+            gene_flow = merged_gene_flow
+        root = visited_leaves[self.populations[0]]['node'].root
+        self.root = root
+        self.gene_flow = gene_flow
+        if DEBUG:
+            print(gene_flow)
+        # self.gene_flow = self.merge_gene_flow(gene_flow, self.root)
+        # print('gene_flow:\n', self.gene_flow)
+        # self.calculate_propotion(gene_flow_propotion_population)
+        return root
+    
+    def calculate_propotion(self, gene_flow, order, visited_leaves, **kwargs):
+        def find_s1_(s1, tar, s2):
+            def find_first_common_ancestor(s1, tar):
+                s1_node = visited_leaves[s1]['node']
+                tar_node = visited_leaves[tar]['node']
+                s1_ancestors = []
+                tar_ancestors = []
+                while s1_node.parent is not None:
+                    s1_ancestors.append(s1_node.parent)
+                    s1_node = s1_node.parent
+                while tar_node.parent is not None:
+                    tar_ancestors.append(tar_node.parent)
+                    tar_node = tar_node.parent
+                common_ancestor = None
+                for node_idx, s1_ancestor in enumerate(s1_ancestors):
+                    if s1_ancestor in tar_ancestors:
+                        common_ancestor = s1_ancestor
+                        break
+                s1__block_nodes = [tar_ancestors[tar_ancestors.index(common_ancestor) - 1]]
+                # s1__block_nodes = tar_ancestors[tar_ancestors.index(common_ancestor)-1:tar_ancestors.index(common_ancestor)]
+                return common_ancestor, s1__block_nodes
+            def bfs_s1_(queue, block_nodes, visited):
+                while queue:
+                    if DEBUG:
+                        print('--------------------')
+                        print(queue)
+                    node = queue.pop(0)
+                    visited.append(node)
+                    if node in block_nodes:
+                        continue
+                    #if node.is_leaf and \
+                    #    (((node.name, s1) not in gene_flow and (s1, node.name) not in gene_flow) or \
+                    #    ((node.name, tar) not in gene_flow and (tar, node.name) not in gene_flow)):
+                    if node.is_leaf:
+                        return node.name
+                    if node.parent and node.parent not in visited:
+                        queue.append(node.parent)
+                    for child in node.children:
+                        if child not in visited:
+                            queue.append(child)
+                return None
+            s1_node = visited_leaves[s1]['node']
+            tar_node = visited_leaves[tar]['node']
+            s2_node = visited_leaves[s2]['node']
+            common_ancestor, s1__block_nodes = find_first_common_ancestor(s1, tar)
+            s1__block_nodes.extend([s1_node, tar_node, s2_node])
+            if DEBUG:
+                print(common_ancestor)
+                print(s1__block_nodes)
+            queue = [tar_node.parent]
+            queue = [s1_node.parent]
+            visited = []
+            s1_ = bfs_s1_(queue, s1__block_nodes, visited)
+            return s1_
+
+        def calculate_m1(s1, tar, s2, pd=False):
+            if DEBUG:
+                print(s1, tar)
+            s1_ = find_s1_(s1, tar, s2)
+            if DEBUG:
+                print('s1: %s, tar: %s, s2: %s, s1_: %s' % (s1, tar, s2, s1_))
+            if s2 is None or s1_ is None:
+                return None, None
+            key = tuple(sorted((s1, s2, tar)))
+            pos = key.index(s1)
+            sign = 1 if (key.index(s2) - pos) % 3 < (key.index(tar) - pos) % 3 else -1
+            D1 = sign * self.triangles[key]['f4'][pos]
+            Z1 = sign * self.triangles[key]['D'][pos]
+
+
+            key = tuple(sorted((s1, s2, s1_)))
+            pos = key.index(s1)
+            sign = 1 if (key.index(s2) - pos) % 3 < (key.index(s1_) - pos) % 3 else -1
+            D2 = sign * self.triangles[key]['f4'][pos]
+            Z2 = sign * self.triangles[key]['D'][pos]
+            # m1 = D1 / D2
+
+            if pd:
+                mix, m_pd = Probability_distribution(D1,D2,Z1,Z2)
+            else:
+                mix = D1 / D2
+                m_pd = None
+
+            if DEBUG:
+                print('D1: %f, D2: %f' % (D1, D2))
+            return mix, m_pd
+
+        bidirectional = kwargs.get('bidirectional', True)
+        calculate_pd = kwargs.get('calculate_pd', False)
+        # calculate_m1('Dushan', 'Baojianshan')
+        # calculate_m1('Baojianshan', 'Dushan')
+        # calculate_m1('Boshan', 'Dushan')
+        # print(self.s1_record)
+        propotion = {}
+        print('Calculating gene flows...')
+        gf_cnt = 0
+        for s2, tar in gene_flow:
+            propotion[(s2, tar)] = []
+            if (s2, tar) not in self.s1_record:
+                continue
+            for s1 in self.s1_record[(s2, tar)]:
+                if DEBUG:
+                    print('-----------------------------------------------') 
+                # pd = False to make it fast
+                mix_bc, m_pd_bc = calculate_m1(s1, tar, s2)
+                if DEBUG:
+                    print('mix:', mix_bc)
+                # mix_a1b, mix_pd_a1b = calculate_m1(s2, tar, s1)
+                # print('mix:', mix_a1b)
+                # if mix_bc is None or mix_a1b is None: continue
+                if mix_bc is None: continue
+                # if mix_bc > 1 or mix_bc < 0 or mix_a1b > 1 or mix_a1b < 0: continue
+                if mix_bc > 1 or mix_bc < 0: continue
+                print('Calculating gene flow %d' % gf_cnt)
+                if calculate_pd:
+                    mix_bc, m_pd_bc = calculate_m1(s1, tar, s2, pd=True)
+                else:
+                    m_pd_bc = float('nan')
+                propotion[(s2, tar)].append((s1, mix_bc, m_pd_bc))
+                gf_cnt += 1
+
+        print()
+        gf_cnt = 0
+        if bidirectional:
+            print('Calculating reversed gene flows...')
+            for tar, s2 in gene_flow:
+                propotion[(s2, tar)] = []
+                if (s2, tar) not in self.s1_record:
+                    continue
+                for s1 in self.s1_record[(s2, tar)]:
+                    mix_bc, m_pd_bc = calculate_m1(s1, tar, s2)
+                    # mix_a1b, mix_pd_a1b = calculate_m1(s2, tar, s1)
+                    # print('mix_bc, mix_a1b', mix_bc, mix_a1b)
+                    # if mix_bc is None or mix_a1b is None: continue
+                    if mix_bc is None: continue
+                    if mix_bc > 1 or mix_bc < 0: continue
+                    # if mix_bc > 1 or mix_bc < 0 or mix_a1b > 1 or mix_a1b < 0: continue
+                    print('Calculating gene flow %d' % gf_cnt)
+                    if calculate_pd:
+                        mix_bc, m_pd_bc = calculate_m1(s1, tar, s2, pd=True)
+                    else:
+                        m_pd_bc = float('nan')
+                    propotion[(s2, tar)].append((s1, mix_bc, m_pd_bc))
+                    gf_cnt += 1
+
+        print()
+        print('---------------------------------------------------------------') 
+        # mat = "{:30}\t{:1}\t{:1}"
+        mat = "{:10}\t{:10}\t{:10}\t{:1}\t{:1}"
+        print(mat.format('s2', 'tar', 's1', 'mix', 'm_pd'))
+        mat = "{:10}\t{:10}\t{:10}\t{:.3f}\t{:.3f}"
+        
+        gene_flows_checked_by_propotion = []
+        for gene_flow, m1s in propotion.items():
+            s2, tar = gene_flow
+            for m1 in m1s:
+                s1, mix, m_pd = m1
+                # print(s2, tar, s1, mix, m_pd)
+                print(mat.format(s2, tar, s1, mix, m_pd))
+                if (s2, tar) not in gene_flows_checked_by_propotion:
+                    gene_flows_checked_by_propotion.append((s2, tar))
+        return gene_flows_checked_by_propotion
+
+    def calculate_propotion_bak2(self, gene_flow, order, visited_leaves, bidirectional=True):
+        def find_s2_s1_(s1, tar):
+            def find_first_common_ancestor(s1, tar):
+                s1_node = visited_leaves[s1]['node']
+                tar_node = visited_leaves[tar]['node']
+                s1_ancestors = []
+                tar_ancestors = []
+                while s1_node.parent is not None:
+                    s1_ancestors.append(s1_node.parent)
+                    s1_node = s1_node.parent
+                while tar_node.parent is not None:
+                    tar_ancestors.append(tar_node.parent)
+                    tar_node = tar_node.parent
+                common_ancestor = None
+                for node_idx, s1_ancestor in enumerate(s1_ancestors):
+                    if s1_ancestor in tar_ancestors:
+                        common_ancestor = s1_ancestor
+                        break
+                s2_block_nodes = s1_ancestors[node_idx-1:node_idx]
+                s1__block_nodes = [common_ancestor]
+                # s1__block_nodes = tar_ancestors[tar_ancestors.index(common_ancestor)-1:tar_ancestors.index(common_ancestor)]
+                return common_ancestor, s2_block_nodes, s1__block_nodes
+            def bfs_s1_(queue, block_nodes):
+                while queue:
+                    if DEBUG:
+                        print('--------------------')
+                        print(queue)
+                    node = queue.pop(0)
+                    if node in block_nodes:
+                        continue
+                    #if node.is_leaf and \
+                    #    (((node.name, s1) not in gene_flow and (s1, node.name) not in gene_flow) or \
+                    #    ((node.name, tar) not in gene_flow and (tar, node.name) not in gene_flow)):
+                    if node.is_leaf:
+                        return node.name
+                    if node.parent:
+                        queue.append(node.parent)
+                    for child in node.children:
+                        queue.append(child)
+                return None
+            def bfs_s2(queue, block_nodes):
+                while queue:
+                    if DEBUG:
+                        print('--------------------')
+                        print(queue)
+                    node = queue.pop(0)
+                    if node in block_nodes:
+                        continue
+                    if node.is_leaf:
+                        return node.name
+                    if node.parent:
+                        queue.append(node.parent)
+                    for child in node.children:
+                        queue.append(child)
+                return None 
+            s1_node = visited_leaves[s1]['node']
+            tar_node = visited_leaves[tar]['node']
+            common_ancestor, s2_block_nodes, s1__block_nodes = find_first_common_ancestor(s1, tar)
+            s2_block_nodes.extend([s1_node, tar_node])
+            s1__block_nodes.extend([s1_node, tar_node])
+            if DEBUG:
+                print(common_ancestor)
+                print(s2_block_nodes)
+                print(s1__block_nodes)
+            queue = [tar_node.parent]
+            s2 = bfs_s2(queue, s2_block_nodes)
+            queue = [s1_node.parent]
+            s1_ = bfs_s1_(queue, s1__block_nodes)
+            return s2, s1_
+
+        def calculate_m1(s1, tar):
+            print('-----------------------------------------------') 
+            if DEBUG:
+                print(s1, tar)
+            s2, s1_ = find_s2_s1_(s1, tar)
+            print('s1: %s, tar: %s, s2: %s, s1_: %s' % (s1, tar, s2, s1_))
+            if s2 is None or s1_ is None:
+                return float('nan')
+            key = tuple(sorted((s1, s2, tar)))
+            pos = key.index(s1)
+            sign = 1 if (key.index(s2) - pos) % 3 < (key.index(tar) - pos) % 3 else -1
+            D1 = sign * self.triangles[key]['f4'][pos]
+
+            key = tuple(sorted((s1, s2, s1_)))
+            pos = key.index(s1)
+            sign = 1 if (key.index(s2) - pos) % 3 < (key.index(s1_) - pos) % 3 else -1
+            D2 = sign * self.triangles[key]['f4'][pos]
+            m1 = D1 / D2
+
+            print('D1: %f, D2: %f' % (D1, D2))
+            return m1
+
+        # calculate_m1('Dushan', 'Baojianshan')
+        # calculate_m1('Baojianshan', 'Dushan')
+        # calculate_m1('Boshan', 'Dushan')
+        propotion = {}
+        for s1, tar in gene_flow:
+            # s2 == s1 is impossible, otherwise there mustn't be a gene flow from s1 to tar
+            m1 = calculate_m1(s1, tar)
+            propotion[(s1, tar)] = {'forward': m1, 'backward': float('nan')}
+
+        if bidirectional:
+            for tar, s1 in gene_flow:
+                m1 = calculate_m1(s1, tar)
+                propotion[(tar, s1)]['backward'] = m1
+        print('-----------------------------------------------') 
+        mat = "{:30}\t{:1}\t{:1}"
+        print(mat.format('gene flow', 'forward', 'backward'))
+        mat = "{:30}\t{:.3f}\t{:.3f}"
+        
+        gene_flows_checked_by_propotion = []
+        for gene_flow, m1 in propotion.items():
+            print(mat.format(str(gene_flow), m1['forward'], m1['backward']))
+            if m1['forward'] and m1['forward'] > 0 and m1['forward'] < 1:
+                gene_flows_checked_by_propotion.append(gene_flow)
+            if m1['backward'] and m1['backward'] > 0 and m1['backward'] < 1:
+                gene_flows_checked_by_propotion.append((gene_flow[1], gene_flow[0]))
+        return gene_flows_checked_by_propotion
+
+
+    def calculate_propotion_bak(self, gene_flow, order, bidirectional=True):
+        print(gene_flow)
+        def find_nearest_population(tar, block_nodes, order):
+            for nodes in order:
+                for n in nodes:
+                    if tar in self.relations[n]:
+                        is_block_node = False
+                        for block_node in block_nodes:
+                            if block_node in self.relations[n]:
+                                is_block_node = True
+                                break
+                        if is_block_node: continue
+                        return self.relations[n][0] if self.relations[n][0] != tar else self.relations[n][1]
+
+        def calculate_m1(s1, tar):
+            print('------------------------------')
+            print(s1, tar)
+            s2 = find_nearest_population(tar, [s1], order)
+            key = tuple(sorted((s1, s2, tar)))
+            pos = key.index(s1)
+            sign = 1 if (key.index(s2) - pos) % 3 < (key.index(tar) - pos) % 3 else -1
+            D1 = sign * self.triangles[key]['f4'][pos]
+
+            s1_ = find_nearest_population(s1, [tar, s2], order)
+            key = tuple(sorted((s1, s2, s1_)))
+            pos = key.index(s1)
+            sign = 1 if (key.index(s2) - pos) % 3 < (key.index(s1_) - pos) % 3 else -1
+            D2 = sign * self.triangles[key]['f4'][pos]
+            m1 = D1 / D2
+
+            print('s1: %s, tar: %s, s2: %s, s1_: %s' % (s1, tar, s2, s1_))
+            print('D1: %f, D2: %f' % (D1, D2))
+            return m1
+
+        propotion = {}
+        for s1, tar in gene_flow:
+            # s2 == s1 is impossible, otherwise there mustn't be a gene flow from s1 to tar
+            m1 = calculate_m1(s1, tar)
+            propotion[(s1, tar)] = {'forward': m1, 'backward': None}
+
+        if bidirectional:
+            for tar, s1 in gene_flow:
+                m1 = calculate_m1(s1, tar)
+                propotion[(tar, s1)]['backward'] = m1
+        print('---------------------') 
+        mat = "{:30}\t{:1}\t{:1}"
+        print(mat.format('gene flow', 'forward', 'backward'))
+        mat = "{:30}\t{:.3f}\t{:.3f}"
+        for gene_flow, m1 in propotion.items():
+            print(mat.format(str(gene_flow), m1['forward'], m1['backward']))
+
+
+#    def calculate_propotion_regression(self, gene_flow_propotion_population):
+#        print('calculate propotion')
+#        out_group = set()
+#        all_statistics = dict()
+#        with open('data/D_ai.log') as f:
+#            for l in f.readlines():
+#                if not l.startswith('result'): continue
+#                A, B, C, D = l.split()[1:5]
+#                for population in [A, B, C, D]: 
+#                    if population not in self.populations:
+#                        out_group.add(population)
+#                key = tuple(sorted((A, B, C, D)))
+#                f4 = float(l.split()[5])
+#                all_statistics[key] = f4 
+#        out_groups = list(combinations(out_group, 2))
+#        for test, R1, R2 in gene_flow_propotion_population:
+#            print(test, R1, R2)
+#            x = []
+#            y = []
+#            for O1, O2 in out_groups:
+#                # print(test, R1, O1, O2)
+#                y.append(all_statistics[tuple(sorted((test, R1, O1, O2)))])
+#                x.append(all_statistics[tuple(sorted((test, R2, O1, O2)))])
+#            x = np.array(x).reshape(-1, 1)
+#            y = np.array(y).reshape(-1, 1)
+#            model = LinearRegression()
+#            model.fit(x, y)
+#
+#            k = model.coef_[0][0]
+#            b = model.intercept_[0]
+#            R2 = model.score(x, y)
+#
+#            print("k", k)
+#            print("a", 1 / (k))
+#            print("b", b)
+#            print("R2", R2)
+
+    def get_dotexporter(self, print_gene_flow=True):
+        options = [] 
+        if print_gene_flow:
+            for l1, l2 in self.gene_flow:
+                options.append("%s -> %s [color=red]" % (l1, l2))
+        return DotExporter(self.root, options=options)
+    
+    def export_dot(self, filename, print_gene_flow=True):
+        self.get_dotexporter(print_gene_flow).to_dotfile(filename)
+
+    def export_graph(self, filename, print_gene_flow=True):
+        self.get_dotexporter(print_gene_flow).to_picture(filename)

utils.py

+from preprocess import get_adjusted_D_statistics
+from itertools import combinations
+import copy
+
+def read_data(data_path, threshold, populations):
+    D_statistics = []
+    adjusted_D_values = []
+    triangles = {}
+    redstar_indices = {}
+    key_visited = []
+    eq_scc = None
+    min_redstar = []
+    s1_record = {}
+    with open(data_path) as f:
+        for l in f.readlines():
+            if not l.startswith('result'): continue
+            A, B, C = l.split()[1:4]
+            if A not in populations or B not in populations or C not in populations: continue
+            f4 = float(l.split()[5])
+            value = float(l.split()[6])
+            get_adjusted_D_statistics(D_statistics, adjusted_D_values, triangles, A, B, C, f4, value, threshold, min_redstar, redstar_indices, key_visited, s1_record)
+            if len(D_statistics) >= len(populations) * (len(populations) - 1) * (len(populations) - 2) / 2:
+                break
+    return D_statistics, adjusted_D_values, triangles, redstar_indices, min_redstar, s1_record
+
+def equal_vertices_judgement(g, vertices):
+    edges = list(combinations(vertices, 2))
+    for v1, v2 in edges:
+#        print(v1, v2)
+#        print(g.graph[v1]['eq'])
+#        print(g.graph[v2]['eq'])
+        if v2 in [eq[0] for eq in g.graph[v1]['eq']] or v1 in [eq[0] for eq in g.graph[v2]['eq']]: continue
+        if v2 in g.graph[v1]['out'] or v2 in g.graph[v1]['in']: return False
+    return True
+
+
+from collections import defaultdict
+
+class Graph:
+    def __init__(self, vertices):
+        self.V = vertices
+        self.adj = defaultdict(list)
+
+    def addEdge(self, u, v):
+        self.adj[u].append(v)
+
+    def DFS(self, v, visited, stack):
+        visited[v] = True
+        stack.append(v)
+
+        for neighbor in self.adj[v]:
+            if not visited[neighbor]:
+                if self.DFS(neighbor, visited, stack):
+                    return True
+            elif neighbor in stack:
+                cycle = []
+                index = stack.index(neighbor)
+                for i in range(index, len(stack)):
+                    cycle.append(stack[i])
+                return cycle
+
+        stack.pop()
+        return False
+
+    def findCycles(self):
+        visited = [False] * self.V
+        stack = []
+        for i in range(self.V):
+            if not visited[i]:
+                result = self.DFS(i, visited, stack)
+                if result != False:
+                    return result
+
+import networkx as nx 
+
+def simple_cycles(G, eq_edges=[], edges_to_adjust=[]):
+    def _unblock(thisnode, blocked, B):
+        stack = {thisnode}
+        while stack:
+            node = stack.pop()
+            if node in blocked:
+                blocked.remove(node)
+                stack.update(B[node])
+                B[node].clear()
+
+    # Johnson's algorithm requires some ordering of the nodes.
+    # We assign the arbitrary ordering given by the strongly connected comps
+    # There is no need to track the ordering as each node removed as processed.
+    # Also we save the actual graph so we can mutate it. We only take the
+    # edges because we do not want to copy edge and node attributes here.
+    subG = type(G)(G.edges())
+    sccs = [scc for scc in nx.strongly_connected_components(subG) if len(scc) > 2]
+
+    # Johnson's algorithm exclude self cycle edges like (v, v)
+    # To be backward compatible, we record those cycles in advance
+    # and then remove from subG
+    for v in subG:
+        if subG.has_edge(v, v):
+            # yield [v]
+            subG.remove_edge(v, v)
+    
+    cycle_ids = set()
+
+    visited_startnodes = set()
+    cycles = []
+    while sccs:
+        scc = sccs.pop()
+        sccG = subG.subgraph(scc)
+        # order of scc determines ordering of nodes
+        # print('scc', scc)
+        for startnode in scc:
+            if startnode in visited_startnodes:
+                continue
+            visited_startnodes.add(startnode)
+            # startnode = scc.pop()
+
+            # Processing node runs "circuit" routine from recursive version
+            path = [startnode]
+            blocked = set()  # vertex: blocked from search?
+            closed = set()  # nodes involved in a cycle
+            blocked.add(startnode)
+            B = defaultdict(set)  # graph portions that yield no elementary circuit
+
+            startnbrs = list(sccG[startnode])
+            for node in copy.deepcopy(startnbrs):
+                if (startnode, node) in eq_edges:
+                    startnbrs.remove(node)
+            if not startnbrs:
+                continue
+            
+            # print(startnode, startnbrs)
+            stack = [(startnode, copy.deepcopy(startnbrs))]  # sccG gives comp nbrs
+
+            while stack:
+                thisnode, nbrs = stack[-1]
+                if nbrs:
+                    nextnode = nbrs.pop()
+                    #print('nextnode', nextnode)
+                    if nextnode == startnode:
+                        # paths.add(''.join(path[:]))
+                        # print(path)
+                        edge_id = [0] * len(edges_to_adjust)
+                        for node_idx in range(len(path)):
+                            edge = (path[node_idx], path[(node_idx + 1) % len(path)])
+                            if edge in edges_to_adjust:
+                                edge_id[edges_to_adjust.index(edge)] = 1
+                        cycle_ids.add(tuple(edge_id))
+                        #yield path[:]
+                        closed.update(path)
+                        # print('path', path)
+                        #print('closed', closed)
+                    #                        #print "Found a cycle", path, closed
+                    elif nextnode not in blocked:
+                        path.append(nextnode)
+                        stack.append((nextnode, list(sccG[nextnode])))
+                        closed.discard(nextnode)
+                        blocked.add(nextnode)
+                        #print('path', path)
+                        continue
+                # done with nextnode... look for more neighbors
+                if not nbrs:  # no more nbrs
+                    if thisnode in closed:
+                        _unblock(thisnode, blocked, B)
+                    else:
+                        for nbr in sccG[thisnode]:
+                            if thisnode not in B[nbr]:
+                                B[nbr].add(thisnode)
+                    stack.pop()
+                    #                assert path[-1] == thisnode
+                    path.pop()
+            # done processing this node
+            # print('paths', paths)
+            H = subG.subgraph(scc)  # make smaller to avoid work in SCC routine
+            H = type(G)(H.edges())
+            # print('remove', startnode, startnbrs)
+            # for node in startnbrs:
+            for node in startnbrs:
+                # print('remove', startnode, node)
+                H.remove_edge(startnode, node) 
+                subG.remove_edge(startnode, node) 
+            # for new_scc in nx.strongly_connected_components(H):
+                # print('new_scc', new_scc)
+            sccs.extend(new_scc for new_scc in nx.strongly_connected_components(H) if len(new_scc) > 2)
+            del H
+            break
+            #print(startnode, 'finished')
+    
+    return cycle_ids 
+
+
+if __name__ == '__main__':
+#    g = Graph(4)
+#    g.addEdge(0, 1)
+#    g.addEdge(0, 2)
+#    g.addEdge(1, 2)
+#    g.addEdge(2, 0)
+#    g.addEdge(2, 3)
+#    g.addEdge(3, 3)
+#
+#    print(g.findCycles())
+
+    import networkx as nx
+
+    # Create Directed Graph
+    G=nx.DiGraph()
+
+    # Add a list of nodes:
+    G.add_nodes_from(["a","b","c","d","e"])
+
+    # Add a list of edges:
+    # G.add_edges_from([("a","b"),("b","c"), ("c","a"), ("b","d"), ("d","e"), ("e","a")])
+    G.add_edges_from([('b','a'), ('c','b'), ("a","b"),("b","c"), ("c","a"), ("b","d"), ("d","a")])
+    # G.add_edges_from([('a','d'), ('a','c'), ('b','a'), ('c','b'), ("a","b"),("b","c"), ("c","a")])
+
+    #Return a list of cycles described as a list o nodes
+    print(list(simple_cycles(G, [('a', 'b'), ('b', 'c'), ('b', 'a'), ('c', 'b')], [('c', 'a')])))
+
+#    G=nx.DiGraph()
+#
+#    # Add a list of nodes:
+#    G.add_nodes_from([1,2,3,4,5])
+#
+#    # Add a list of edges:
+#    G.add_edges_from([(1,2), (2,1)])
+
+    #Return a list of cycles described as a list o nodes
+#    cycles = list(nx.simple_cycles(G))
+#    print(cycles)
+#    new_cycles = []
+#    for cycle in cycles:
+#        edges = []
+#        for i in range(len(cycle)):
+#            edges.append((cycle[i], cycle[(i+1) % len(cycle)]))
+#        new_cycles.append(edges)
+#    cycles = new_cycles
+#
+#    import cvxopt
+#    from cvxopt import matrix
+#    from cvxopt.glpk import ilp
+#    import numpy as np
+#
+#    num_edges = 200
+#    num_cycles = 100
+#    edges = list(G.edges)
+#    c = matrix(np.ones(num_edges, dtype=float))
+#    G2 = matrix(np.zeros((num_cycles, num_edges), dtype=float))
+#    print(cycles)
+#    print(edges[0])
+#    for i in range(num_cycles):
+#        for j in range(num_edges):
+#            if edges[j % 6] in cycles[i % 2]:
+#                G2[i,j] = -1
+#    h = matrix(-1 * np.ones(num_cycles, dtype=float))
+#    I = set()
+#    B = set(range(num_edges))
+#    (status,x)=ilp(c,G2,h,B=B)
+#    print(status)
+#    print(x)
+#
\ No newline at end of file