def mean_geodesic(pg, debug=0):
    """
    mean_geodesic() calculates the mean geodesic (shortest) distance
    between two vertices in a network.
    """
    length_sum = 0
    if networkx.is_directed_acyclic_graph(pg):
        n_pairs_with_paths = 0
    else:
        n_pairs_with_paths = ( pg.order() * ( pg.order() + 1 ) ) / 2
    tg = networkx.subgraph(pg, pg.nodes())
    for u in pg.nodes_iter():
        tg.delete_node(u)
        for v in tg.nodes_iter():
            try:
                length = networkx.shortest_path_length(pg,u,v)
                if length > 0:
                    length_sum = length_sum + length
                    if networkx.is_directed_acyclic_graph(pg):
                        n_pairs_with_paths = n_pairs_with_paths + 1
            except networkx.exception.NetworkXError:
                pass
    try:
        geodesic = float(length_sum) / float(n_pairs_with_paths)
    except:
        geodesic = -999.
    if debug:
        print 'length_sum:\t', length_sum
        print 'n_pairs_with_paths:\t', n_pairs_with_paths
    return geodesic

def comb_fas( graph):
    '''@param: graph, a nx.DiGraph obj
    '''
    assert isinstance( graph, nx.DiGraph)
    origin_weight = nx.get_edge_attributes( graph, 'weight')
    weight = origin_weight.copy()

    assert len(weight) == graph.number_of_edges(), "Some edge doesnot has a weight attr."
    fas = []
    while( not nx.is_directed_acyclic_graph(graph) ):
        c = list( nx.simple_cycles(graph) )[0]
        mini_weight = min( [ weight[edge] for edge in get_edges(c)] )

        cycle_edges_weight = {edge:weight[edge] for edge in get_edges(c) }
        for eachEdge in cycle_edges_weight.keys():
            cycle_edges_weight[eachEdge] -= mini_weight
            weight[eachEdge ] -= mini_weight
            if cycle_edges_weight[eachEdge] == 0:
                fas.append( eachEdge )
                graph.remove_edge( eachEdge[0], eachEdge[1] )

    for eachEdge in copy.copy(fas):
        graph.add_edge( eachEdge[0], eachEdge[1], {'weight' : origin_weight[eachEdge]} )
        if nx.is_directed_acyclic_graph( graph):
            fas.remove(eachEdge)
            continue
        else:
            graph.remove_edge( eachEdge[0], eachEdge[1] )

    return fas

    def test_topological_sort2(self):
        DG = nx.DiGraph({1: [2], 2: [3], 3: [4],
                         4: [5], 5: [1], 11: [12],
                         12: [13], 13: [14], 14: [15]})
        assert_raises(nx.NetworkXUnfeasible, consume, nx.topological_sort(DG))

        assert_false(nx.is_directed_acyclic_graph(DG))

        DG.remove_edge(1, 2)
        consume(nx.topological_sort(DG))
        assert_true(nx.is_directed_acyclic_graph(DG))

def minkowski_causality(D,N,show_plot=False):
    """
        Instantiates "event_in_minkowski" to return a list of N points.
        """
    def points_in_minkowski(D,N):
        
        points_in_minkowski = []
        
        n=1
        while n<=N:
            point_n = event_in_minkowski(D)
            coords_n = point_n.coord_value_point_n(n)
            points_in_minkowski.append(coords_n)
            n+=1
        return points_in_minkowski
    
    good_points = points_in_minkowski(D,N)
    
    #List --> Dict as nx needs hashable object to add nodes/edges from.
    dict_of_points = {}
    for i in range(len(good_points)):
        good_points[i] = tuple(good_points[i])
        dict_of_points[i] = good_points[i]
    
    #Add nodes to empty nx graph object
    G=nx.DiGraph()
    for point in dict_of_points:
        G.add_node(point)
    print nx.is_directed_acyclic_graph(G)

    #Add edge (from i to j) to empty nx graph object if node j falls within the future light cone of i
    for i in range(len(dict_of_points)):
        for j in range(len(dict_of_points)):
            if i==j:
                continue
            t_separation = dict_of_points[j][0] - dict_of_points[i][0]
            space_separation=0
            for d in range(1,D):
                space_separation += (dict_of_points[i][d] - dict_of_points[j][d])
            if t_separation>=abs(space_separation):
                G.add_edge(i,j)
            else:
                pass

    #Check G is a DAG, print model info

    if nx.is_directed_acyclic_graph(G):
        print "This is a DAG of causal relations between randomly placed events in ",D,"D Minkowski space-time."

    #Show plot
    if show_plot==True:
        draw_in_minkowski(G,dict_of_points)

    return G

    def test_topological_sort2(self):
        DG = nx.DiGraph({1: [2], 2: [3], 3: [4], 4: [5], 5: [1], 11: [12], 12: [13], 13: [14], 14: [15]})
        assert_raises(nx.NetworkXUnfeasible, nx.topological_sort, DG)
        assert_raises(nx.NetworkXUnfeasible, nx.topological_sort_recursive, DG)

        assert_false(nx.is_directed_acyclic_graph(DG))

        DG.remove_edge(1, 2)
        assert_equal(nx.topological_sort_recursive(DG), [11, 12, 13, 14, 15, 2, 3, 4, 5, 1])
        assert_equal(nx.topological_sort(DG), [11, 12, 13, 14, 15, 2, 3, 4, 5, 1])
        assert_true(nx.is_directed_acyclic_graph(DG))

def make_acyclic(G):
	G_copy = G.copy()
	F = []
	original_G = G.copy()
	while not nx.is_directed_acyclic_graph(G_copy):
		#iterate through cycles in G
		
		for cycle in nx.simple_cycles(G_copy):
			min_weight = 100000
			min_u = 0
			min_v = 0
			#Find minimum weight edge in the cycle, weight
			#here is bundle size
			#TODO: start with smallest cycle by sorting
			#print G.edges(data=True)
			for i in xrange(0,len(cycle)-1):
				u = cycle[i]
				v = cycle[i+1]
				if G[u][v]['bsize'] < min_weight:	
					min_weight = G[u][v]['bsize']
					min_u = u
					min_v = v
			if G[cycle[- 1]][cycle[0]]['bsize'] < min_weight:
				min_weight = G[cycle[-1]][cycle[0]]['bsize']
				min_u = cycle[-1]
				min_v = cycle[0]

			#reduce the edge weights by min_weight and remove the edge if its weight is 0
			if min_weight != 100000:
				for i in xrange(0,len(cycle)-1):
					u = cycle[i]
					v = cycle[i+1]
					G[u][v]['bsize'] -= min_weight
				
				G[cycle[-1]][cycle[0]]['bsize'] -= min_weight
				G.remove_edge(min_u,min_v)
				F.append((min_u,min_v,original_G.get_edge_data(min_u,min_v)))
				G_copy = G.copy()
				break

	#Now try adding edges from F to G, TODO do in non-increasing order

		if len(G.edges()) == 0:
			continue
		# if len(G.nodes()) == 0:
		# 	continue
		for edge in F:
			u = edge[0]
			v = edge[1]
			G.add_edge(u,v,edge[2])
			if not nx.is_directed_acyclic_graph(G):
				G.remove_edge(u,v)

	return G

def make_dag(g):
    if nx.is_directed_acyclic_graph(g):
        return
    p = nx.periphery(g)
    for c in nx.weakly_connected_component_subgraphs(g):
        if nx.is_directed_acyclic_graph(g):
            continue
        cycles = nx.simple_cycles(c)
        for c in cycles:
            edges = zip(c[:-1], c[1:])
            edges.append((c[-1], c[0]))
            for e in edges:
                data = g.edges(e[0], e[1])[0][2]
                c.remove_edge(e[0], e[1])

def longest_subsequence_dag(a, sign):
    '''Return a longest increasing (if sign=1) or decreasing (if sign=-1) sub-sequence in the
    permutation a of the first n natural integers. Time and storage are O(n). If multiple longest
    sub-sequences exist, arbitrarily returns one of them.'''

    # Dan Cook's idea: use symmetry to solve the decreasing case in terms of the increasing case 
    if sign < 0:
        return list(reversed(longest_subsequence_dag(list(reversed(a)), 1)))
    
    G = build_dag(np.array(a))  # Construct a DAG whose edges represent all candidate pairs of consecutive elements of the longest subsequence
    assert nx.is_directed_acyclic_graph(G)
    # print 'Edges', G.edges()

    depth = longest_path_length(G)  # For each node, calculate the longest path length       
    # print 'depth', depth 
    
    # Back-track from a node of maximum depth to its ancestors to reconstruct the longest path
    x = np.argmax(depth)
    seq = [x]
    # print 'x', x, 'depth', depth[x]
    while G.in_degree(x) > 0:
        # To find the maximum path, choose a parent of minimum depth
        parents = G.predecessors(x)
        # print 'parents', parents
        x = parents[np.argmax(depth[parents])]
        # print 'x', x, 'depth', depth[x]
        seq.append(x)
        # print 'seq', seq
    # print 'final seq', list(reversed(seq))
    return list(reversed(seq))

	def generate(self):
		"""Workhorse factory method for producing all valid DAGs for this schema and set
		of constraints."""
		graphs = []

		sys.stderr.write("gen:\n" + self.dumpEdgePossibleSettings())

		edgesPossible = self.getAllVarPairs()		
		edgeCombos = factorialDict(self.edgePossibleSettings)
		for edgeCombo in edgeCombos: # edgeCombo is a dict (s, t) -> 1
			graph = nx.DiGraph()
			for i, ev in enumerate(self.entVars):
				lat = True if (ev in self.latents) else False
				det = self.determines.get(ev, None)
				graph.add_node(ev, latent=lat, determines=det, order=i) # order they were passed in is preserved

			graphSig = ""	
			for s, t in edgesPossible:
				setting = edgeCombo.get((s, t), 0)
				if (setting == 1):
					graph.add_edge(s, t)
				elif (setting == 2):
					graph.add_edge(t, s)
				if (s in self.indexSet) and (t in self.indexSet):
					graphSig += str(setting)
			graph.graph['index'] = int(graphSig, 3)
			
			if (not self.dagsOnly) or (nx.is_directed_acyclic_graph(graph)):
				graphs.append(graph)
		if (len(graphs) < len(edgeCombos)):
			sys.stderr.write("eliminated %d cyclic graphs\n" % (len(edgeCombos) - len(graphs)))
		return sorted(graphs, key=lambda x: x.graph['index'])

    def __init__(self, tasks_reqs):
        """Construct a PipelineFramework based on the given Tasks and their requirements.

        A PipelineFramework is the structure of the pipeline, it contains no patient data.

        :param tasks_reqs: the Tasks and their requirements
        :type tasks_reqs: iterable of tuples, each with a Task and its list of required UIDs
        :raises: ValueError
        """
        self.dag = DiGraph()
        task_dict = {}
        for task, _ in tasks_reqs:
            if task_dict.get(task._uid) is not None:
                raise ValueError("Pipeline contains duplicate Task {}".format(task._uid))
            self.dag.add_node(task, done=False)
            task_dict[task._uid] = task

        for task, reqs in tasks_reqs:
            for req_uid in reqs:
                uid = task_dict.get(req_uid)
                if uid is None:
                    raise KeyError("Unknown UID {} set as requirement for {}".format(req_uid, task._uid))
                self.dag.add_edge(uid, task)

        if not is_directed_acyclic_graph(self.dag):
            raise ValueError("Pipeline contains a cycle.")

def balance(graph):
    '''parame: graph, a DAG its.__class__ == nx.DiGraph
       return: r,     removed edges set so makr the input graph a b-structure
    '''
    # 只处理整数形式的图，每一个整数对应的节点可以在后面查到
    # 输入进来的图应该是连通的，如果存在非连通图，minimum_edge_cut就会产生问题
    assert nx.is_directed_acyclic_graph(graph),\
        "The target graph you want to banlance is not a DAG"
    r = [] # removed set
    if check(graph):
        return r
    #非B-Stucture时，一直循环下去
    # BUGY: 如果cs为空呢，那么不可能有两个图返回来，这时候怎么办
    print "\nCutting Graph"
    cs, g1, g2 = cut(graph) 
    r = balance(g1) + balance(g2) + cs
    csl = []
    for eachEdge in cs:
        under_check_graph = graph.copy()
        under_check_graph.remove_edges_from(r)
        under_check_graph.add_edges_from(csl)
        under_check_graph.add_edge(eachEdge[0],eachEdge[1])
        if check(under_check_graph):
            print "Edge: %s added back" % str(eachEdge)
            csl.append(eachEdge)
            graph.add_edge(eachEdge[0],eachEdge[1])
    for eachEdge in csl:
        r.remove(eachEdge)
    print "Removed Edge Set: %s" % str(r)
    return r

 def set_params(self, node, delta, eta, marginal):
     self.clear_memory()
     assert node in self.graph.nodes()
     self.graph.node[node]['delta'] = delta
     self.graph.node[node]['eta'] = eta
     self.graph.node[node]['marginal'] = marginal
     assert nx.is_directed_acyclic_graph(self.graph)

    def build_graph(self):
        """Build graph of relationships between hills 

        Each hill is a list of things that can be used for that hill.

        Each of these may have inputs (names of other hills).

        A graph is build to show the input relations.

        Checks in case graph is cyclic.

        Does a topological sort on the hills to give and order in
        which they should be processed.
        """
        graph = nx.DiGraph()

        for hill, data in self.hills.items():

            for item in data:
                for link in item.inputs:
                    graph.add_edge(link, hill)

        # check if graph is acyclic
        is_dag = nx.is_directed_acyclic_graph(graph)

        if not is_dag:
            raise ValueError("hills must be acyclic")

        self.hill_order = nx.topological_sort(
            graph)

def _validate(G):
    '''
    Validates dependency graph to ensure it has no missing or cyclic dependencies
    '''
    for name in G.nodes():
        if 'value' not in G.node[name] and 'template' not in G.node[name]:
            msg = 'Dependency unsatisfied in variable "%s"' % name
            raise ParamException(msg)

    if not nx.is_directed_acyclic_graph(G):
        graph_cycles = nx.simple_cycles(G)

        variable_names = []
        for cycle in graph_cycles:
            try:
                variable_name = cycle[0]
            except IndexError:
                continue

            variable_names.append(variable_name)

        variable_names = ', '.join(sorted(variable_names))
        msg = ('Cyclic dependency found in the following variables: %s. Likely the variable is '
               'referencing itself' % (variable_names))
        raise ParamException(msg)

    def set_indices(self):
        if not nx.is_directed_acyclic_graph(self):
            raise ValueError('The graph is not DAG')
        if not nx.is_connected(self.to_undirected()):
            raise ValueError('The graph is not connected')

        self.base_digraph = nx.DiGraph(self)
        self.ordered_nodes = nx.topological_sort(self)
        for idx, node in enumerate(self.ordered_nodes):
            self.node[node]['index'] = idx

        self.ordered_edges = OrderedDict({})

        index = 0
        for tail in self.ordered_nodes:
            for head in sorted(self[tail]):
                self.ordered_edges[(tail, head)] = index
                self.base_digraph[tail][head]['capacity'] = len(self[tail][head])
                for idx in sorted(self[tail][head]):
                    self[tail][head][idx]['index'] = index
                    index = index + 1

        # reset data structures
        self.coding_matrix = None
        self.dst_evolution_rec = None
        self.alignment_nodes = []

def count_common_subgraphs(graph1, graph2, n1, n2,
                     node_attrib='label', edge_attrib='label'):
    """
    Counts the number of common (dependency parse) subgraphs rooted at n1 and
    n2. This is an implementation of Cm(n1, n2) for dependency structures from
    Collins and Duffy (2001). Parsing with a Single Neuron.
    """
    for graph in (graph1, graph2):
        assert nx.is_directed_acyclic_graph(graph)
    
    if graph1.node[n1][node_attrib] != graph2.node[n2][node_attrib]:
        return 0

    n1_children = dependency_children(graph1, n1, edge_attrib=edge_attrib)
    n2_children = dependency_children(graph2, n2, edge_attrib=edge_attrib)

    if not n1_children or not n2_children:
        return 0
    else:
        result = 1  # neutral element of multiplication
        for n1_target, n2_target in common_dependency_targets(graph1, graph2, n1, n2,
                                                        node_attrib=node_attrib):
            result *= (count_common_subgraphs(graph1, graph2,
                                        n1_target, n2_target,
                                        node_attrib='label',
                                        edge_attrib='label') + 2)
        return result - 1

def mean_degree_centrality(pg, normalize=0):
    """
    mean_degree_centrality(pg) calculates mean in- and out-degree
    centralities for directed graphs and simple degree-centralities
    for undirected graphs. If the normalize flag is set, each node's
    centralities are weighted by the number of edges in the (di)graph.
    """
    centrality = {}
    try:
        if networkx.is_directed_acyclic_graph(pg):
            cent_sum_in, cent_sum_out = 0, 0
            for n in pg.nodes():
                n_cent_in = pg.in_degree(n)
                n_cent_out = pg.out_degree(n)
                if normalize:
                    n_cent_in = float(n_cent_in) / float(pg.size()-1)
                    n_cent_out = float(n_cent_out) / float(pg.size()-1)
                cent_sum_in = cent_sum_in + n_cent_in
                cent_sum_out = cent_sum_out + n_cent_out
            centrality['in'] = cent_sum_in / float(pg.order())
            centrality['out'] = cent_sum_out / float(pg.order())
        else:
            cent_sum = 0
            for n in pg.nodes():
                if not normalize:
                    n_cent = pg.degree(n)
                else:
                    n_cent = networkx.degree_centrality(pg,n)
                cent_sum = cent_sum + n_cent
            centrality['all'] = cent_sum / float(pg.order())
    except:
        logging.error('pyp_network.mean_degree_centrality() failed!')
    return centrality

def graph(dataframe=None):

    G = nx.DiGraph()
    
    nrow = dataframe.shape[0]

    for i in xrange(nrow):
        
        source = dataframe['module_id'][i]
        G.add_node(source)
        
        if pd.isnull(dataframe['children'][i]) is not True:
            try:
                targets = dataframe['children'][i].split()
            except:
                raise ValueError('Data type is not correct:', i, 
                        dataframe.loc[i,], type(dataframe['children'][i]))
            
            for key in targets:
                G.add_edge(source, key)

    # Sanity check
    selfLoop = G.selfloop_edges()
    assert len(selfLoop) == 0, ValueError('self loop:', selfLoop)
    assert nx.is_directed_acyclic_graph(G), valueError('loop exists!')

    return G

def dyad_census(pg, debug=0, debuglog=0):
    """
    dyad_census() calculates the number of null, asymmetric, and
    mutual edges between all pairs of nodes in a directed graph.
    """
    if not networkx.is_directed_acyclic_graph(pg):
        logging.error('pyp_network.dyad_census() requires a directed graph as input!')
        return 0
    else:
        census = {}
        census['null'] = 0
        census['asymmetric'] = 0
        census['mutual'] = 0
        tg = networkx.subgraph(pg, pg.nodes())
        for u in pg.nodes_iter():
            tg.delete_node(u)
            for v in tg.nodes_iter():
                if not pg.has_neighbor(u,v):
                    census['null'] = census['null'] + 1
                elif u in pg.predecessors(v) and v in pg.successors(u):
                    census['mutual'] = census['mutual'] + 1
                    if debug:
                        print 'Nodes %s and %s link to one another!' % ( u, v )
                    if debuglog:
                        logging.error('Nodes %s and %s link to one another!',u, v)
                elif u in pg.predecessors(v) and v not in pg.successors(u):
                    census['asymmetric'] = census['asymmetric'] + 1
                elif u not in pg.predecessors(v) and v in pg.successors(u):
                    census['asymmetric'] = census['asymmetric'] + 1
                else:
                    pass
        del(tg)
        return census

def read_pedigree_from_test_file(file_name, genotyped_id_file=None):
    '''Load a pedigree from a PLINK TFAM file.'''
    data = np.genfromtxt(file_name, np.dtype(int))
    p = io_pedigree.read(file_name, genotyped_id_file=genotyped_id_file)
    assert_equal(p._graph.number_of_nodes(), data.shape[0], 'Incorrect number of nodes')
    assert nx.is_directed_acyclic_graph(p._graph), 'Pedigree is not a DAG'
    return p