def _get_random_test_setup(nstates):
    """

    Returns
    -------
    T : undirected networkx graph
        Edges are annotated with transition matrix P.
    root : integer
        Root node.
    root_distn : dict
        Probability distribution at the root.
    node_to_allowed_states : dict
        Map from node to set of allowed states.

    """
    # Sample a random tree.
    branching_distn = [0.7, 0.1, 0.1, 0.1]
    T = get_random_branching_tree(branching_distn, maxnodes=6)
    root = 0

    # For each edge on the tree,
    # sample a random sparse state transition matrix.
    for na, nb in nx.bfs_edges(T, root):
        T[na][nb]['P'] = _get_random_nx_transition_matrix(nstates)

    # Sample a root distribution.
    # It should be a little bit sparse, for testing.
    weights = np.random.exponential(size=nstates)
    imissing = np.random.randint(nstates)
    pairs = [(i, w) for i, w in enumerate(weights) if i != imissing]
    weights[imissing] = 0
    total_weight = np.sum(weights)
    root_distn = dict((i, w / total_weight) for i, w in pairs)

    # Sample allowed states at each node.
    # Disallow a random state at each node.
    states = range(nstates)
    node_to_allowed_states = dict((n, set(states)) for n in T)
    for n in T:
        imissing = np.random.randint(nstates)
        node_to_allowed_states[n].remove(imissing)

    # Final check on transition matrices on edges of T.
    for na, nb in nx.bfs_edges(T, root):
        edge_object = T[na][nb]
        P = edge_object.get('P', None)
        if P is None:
            raise Exception('internal error')

    # Return the random info for testing.
    return T, root, root_distn, node_to_allowed_states

def purge(graph, seeds):
    top = top_nodes(graph)

    dead = set(seeds)
    alive = top.difference(dead)

    dead_tree = nx.DiGraph()
    for n in dead:
        dead_tree.add_edges_from(nx.bfs_edges(graph, n))
    alive_tree = nx.DiGraph()
    for n in alive:
        alive_tree.add_edges_from(nx.bfs_edges(graph, n))

    return set(dead_tree.nodes()).difference(alive_tree.nodes())

def _sample_states_preprocessed(T, edge_to_P, root,
        v_to_subtree_partial_likelihoods):
    """
    Jointly sample states on a tree.

    This variant requires subtree partial likelihoods.

    """
    root_partial_likelihoods = v_to_subtree_partial_likelihoods[root]
    if not root_partial_likelihoods:
        return None
    v_to_sampled_state = {}
    v_to_sampled_state[root] = dict_random_choice(root_partial_likelihoods)
    for edge in nx.bfs_edges(T, root):
        va, vb = edge
        P = edge_to_P[edge]

        # For the relevant parent state,
        # compute an unnormalized distribution over child states.
        sa = v_to_sampled_state[va]

        # Construct conditional transition probabilities.
        fset = set(P[sa]) & set(v_to_subtree_partial_likelihoods[vb])
        sb_weights = {}
        for sb in fset:
            a = P[sa][sb]['weight']
            b = v_to_subtree_partial_likelihoods[vb][sb]
            sb_weights[sb] = a * b

        # Sample the state using the unnormalized dictionary of weights.
        v_to_sampled_state[vb] = dict_random_choice(sb_weights)

    return v_to_sampled_state

def get_expm_augmented_tree(T, root, Q_default=None):
    """
    Add transition probability matrices to edges.

    Construct the augmented tree by annotating each edge
    with the appropriate state transition probability matrix.

    Parameters
    ----------
    T : weighted undirected networkx graph
        This tree is possibly annotated with edge-specific
        rate matrices Q.
    root : integer
        Root node.
    Q_default : 2d ndarray, optional
        Default rate matrix.

    Returns
    -------
    T_aug : weighted undirected networkx graph
        Tree annotated with transition probability matrices P.

    """
    T_aug = nx.Graph()
    for na, nb in nx.bfs_edges(T, root):
        edge = T[na][nb]
        weight = edge['weight']
        Q = edge.get('Q', Q_default)
        _density.check_square_dense(Q)
        P = custom_expm(Q, weight)
        T_aug.add_edge(na, nb, weight=weight, P=P)
    return T_aug

def _sample_states_preprocessed(T, edge_to_P, root,
        v_to_subtree_partial_likelihoods):
    """
    Jointly sample states on a tree.

    This variant requires subtree partial likelihoods.

    """
    root_partial_likelihoods = v_to_subtree_partial_likelihoods[root]
    n = root_partial_likelihoods.shape[0]
    if not root_partial_likelihoods.any():
        return None
    distn1d = normalized(root_partial_likelihoods)
    root_state = weighted_choice(n, p=distn1d)
    v_to_sampled_state = {root : root_state}
    for edge in nx.bfs_edges(T, root):
        va, vb = edge
        P = edge_to_P[edge]

        # For the relevant parent state,
        # compute an unnormalized distribution over child states.
        sa = v_to_sampled_state[va]

        # Construct conditional transition probabilities.
        sb_weights = P[sa] * v_to_subtree_partial_likelihoods[vb]

        # Sample the state.
        distn1d = normalized(sb_weights)
        v_to_sampled_state[vb] = weighted_choice(n, p=distn1d)

    return v_to_sampled_state

 def compute_overlap(self):
     self.middle_modules = None
     self.mean_overlap = None
     if self.size() == 0:
         return
     self.middle_modules = list()
     for j in range(self.parameters.nodes_in):
         if not j in self:
             continue
         middle = set()
         for (u, v) in nx.bfs_edges(self, j):
             if self.parameters.is_middle(u):
                 middle.add(u)
             if self.parameters.is_middle(v):
                 middle.add(v)
         self.middle_modules.append(middle)
     num_in = len(self.middle_modules)
     if num_in == 0:
         return
     # compute overlap
     combinations = (num_in * (num_in - 1)) / 2
     middle_overlap = numpy.zeros(combinations)
     i = 0
     for j in range(num_in - 1):
         for k in range(j + 1, num_in):
             nominator = len(self.middle_modules[j].intersection(self.middle_modules[k]))
             denominator = float(len(self.middle_modules[j].union(self.middle_modules[k])))
             if denominator == 0.0:
                 middle_overlap[i] = numpy.nan
             else:
                 middle_overlap[i] = nominator / denominator
             i += 1
     middle_overlap = numpy.ma.masked_invalid(middle_overlap)
     self.mean_overlap = numpy.mean(middle_overlap[~middle_overlap.mask])

def _forward(T, edge_to_P, root, v_to_subtree_partial_likelihoods):
    """
    Forward pass.

    Return a map from node to posterior state distribution.

    """
    root_partial_likelihoods = v_to_subtree_partial_likelihoods[root]
    v_to_posterior_distn = {}
    v_to_posterior_distn[root] = dict_distn(root_partial_likelihoods)
    for edge in nx.bfs_edges(T, root):
        va, vb = edge
        P = edge_to_P[edge]

        # For each parent state, compute the distribution over child states.
        distn = defaultdict(float)
        parent_distn = v_to_posterior_distn[va]
        for sa, pa in parent_distn.items():

            # Construct conditional transition probabilities.
            fset = set(P[sa]) & set(v_to_subtree_partial_likelihoods[vb])
            sb_weights = {}
            for sb in fset:
                a = P[sa][sb]['weight']
                b = v_to_subtree_partial_likelihoods[vb][sb]
                sb_weights[sb] = a * b
            sb_distn = dict_distn(sb_weights)

            # Add to the marginal distribution.
            for sb, pb in sb_distn.items():
                distn[sb] += pa * pb

        v_to_posterior_distn[vb] = dict(distn)

    return v_to_posterior_distn

    def _color_by_mfpt(self, min1, T=1.):
        print "coloring by the mean first passage time to get to minimum", min1._id
        # get a list of transition states in the same cluster as min1
        edges = nx.bfs_edges(self.graph, min1)
        transition_states = [ self.graph.get_edge_data(u, v)["ts"] for u, v in edges ]
        if not check_thermodynamic_info(transition_states):
            raise Exception("The thermodynamic information is not yet computed")

        
        # get an arbitrary second minimum2
        for ts in transition_states:
            if ts.minimum2 != min1:
                min2 = ts.minimum2
                break
        A = [min1]
        B = [min2]
        rcalc = RatesLinalg(transition_states, A, B, T=T)
        rcalc.compute_rates()
        mfptimes = rcalc.get_mfptimes()
        tmax = max(mfptimes.itervalues())
        def get_mfpt(m):
            try:
                return mfptimes[m]
            except KeyError:
                return tmax
        self.dg.color_by_value(get_mfpt)
        self.redraw_disconnectivity_graph()

def get_edge_to_fvec2d(*args):
    """
    For each edge, get the joint feasibility of states at edge endpoints.

    Parameters
    ----------
    {params}

    Returns
    -------
    edge_to_fvec2d : map from directed edge to networkx DiGraph
        For each directed edge in the rooted tree report the networkx DiGraph
        among states, for which presence/absence of an edge defines the
        posterior feasibility of the corresponding state transition
        along the edge.

    """
    args = validated_params(*args)
    T, edge_to_A, root, root_prior_fvec1d, node_to_data_fvec1d = args

    v_to_fvec1d = get_node_to_fvec1d(*args)
    edge_to_fvec2d = {}
    for edge in nx.bfs_edges(T, root):
        va, vb = edge
        A = edge_to_A[edge]
        fa = v_to_fvec1d[va]
        fb = v_to_fvec1d[vb]
        edge_to_fvec2d[edge] = A & np.outer(fa, fb)
    return edge_to_fvec2d

def get_expm_augmented_tree(T, root, Q_default=None):
    """
    Add transition probability matrices to edges.

    Construct the augmented tree by annotating each edge
    with the appropriate state transition probability matrix.

    Parameters
    ----------
    T : weighted undirected networkx graph
        This tree is possibly annotated with edge-specific
        rate matrices Q.
    root : integer
        Root node.
    Q_default : weighted directed networkx graph, optional
        Sparse rate matrix.

    Returns
    -------
    T_aug : weighted undirected networkx graph
        Tree annotated with transition probability matrices P.

    """
    T_aug = nx.Graph()
    for na, nb in nx.bfs_edges(T, root):
        edge = T[na][nb]
        weight = edge['weight']
        Q = edge.get('Q', Q_default)
        if Q is None:
            raise ValueError('no rate matrix is available for this edge')
        P = sparse_expm(Q, weight)
        T_aug.add_edge(na, nb, weight=weight, P=P)
    return T_aug

def km_random(g,k=5,m=3,start=None):
    """ k nodes of breath first sequence; m add and del number."""
    if start==None:
        start=g.nodes().pop()
    bfList=list(nx.bfs_edges(g,start))
    bfList.reverse()
    bfList.append((start,start))
    tempk=[]
    try:
        while bfList:
            for each in range(k):
                tempk.append(bfList.pop()[1])
            
            tg=nx.subgraph(g,tempk)
            e=del_edge(tg,m)
            g.remove_edges_from(e)

            tg=nx.subgraph(g,tempk)
            e=add_edge(tg,m)
            g.add_edges_from(e)

            tempk=[]

    except IndexError:
        print "pop finishing"

    def ordering_graph(self):
        """Ordering graph

        t1 --> t2 in the ordering graph indicates that t1 happens before t2.
        A missing edge simply means that it is not clear yet.

        """

        g = nx.DiGraph()

        # add times
        for t in self.nodes_iter():
            g.add_node(t)

        # add existing edges
        for t1, t2 in self.edges_iter():
            g.add_edge(t1, t2)

        # connect every pair of anchored times
        anchored = sorted(self.anchored())
        for t1, t2 in itertools.combinations(anchored, 2):
            g.add_edge(t1, t2)

        # connect every time with its sucessors
        _g = g.copy()
        for t1 in _g:
            for t2 in set([target for (_, target) in nx.bfs_edges(_g, t1)]):
                g.add_edge(t1, t2)

        return g

def firstCommonAncestors(dag, synListsSet):
	"""Return the list of first common ancestors"""
	i = 0
	nodesAndDistance = dict()
	
	for synList in synListsSet:
		for a, b in nx.bfs_edges(dag, synList, reverse=True):
				dag[b][i] = True
		i += 1
	
	for node in dag.nodes():
		test = True
		for key in xrange(i):
			if not key in dag[node]:
				test = False
		if test:
			nodesAndDistance[node] = len(nx.shortest_path(dag, target=node))
	
	maxDistance = max(nodesAndDistance.values())
	
	for e in nodesAndDistance.keys():
		if nodesAndDistance[e] < maxDistance:
			del nodesAndDistance[e]
	
	for key in xrange(i):
		removeKey(dag, key)
	
	return nodesAndDistance.keys()

 def information_diffusion(self, num_nodes, beta):
     # 		patients_zero = [random.randint(0,num_nodes) for r in xrange(beta)]
     # 		for i in patients_zero:
     # 			self.network.node[i]['color'] = 1
     # 		print patients_zero
     # 		for i in patients_zero:
     # 			for j in self.network.neighbors(i):
     root_node = random.randint(0, num_nodes - 1)
     self.network.node[root_node]["color"] = 1
     ordered_edges = list(nx.bfs_edges(self.network, root_node))
     print ordered_edges
     t_name = "plots/file_name"
     count = 0
     for i in ordered_edges:
         count = count + 1
         # 			print self.network.node[i[0]]['color']==1,  self.network.node[i[1]]['color']==0
         if self.network.node[i[0]]["color"] == 1 and self.network.node[i[1]]["color"] == 0:
             # 				probability =100* self.network.node[i[1]]['mew_final']*self.network.edge[i[0]][i[1]]['gossip']
             probability = random.random()
             print i, probability
             if probability > beta:
                 # 					print "hello from other side"
                 self.network.node[i[1]]["color"] = 1
         if count % 100 == 0:
             name = t_name + str(count) + ".gml"
             nx.write_gml(self.network, name)

    def chooseplan(self, costfunc=None):
        """Return a join sequence object based on the join graph.  This one is
        simple -- it just adds the joins according to a breadth first search"""
        # choose first node in the original insertion order (networkx does not
        # guarantee even a deterministic order)
        firstnode = [n for n in self.joingraph.nodes() if n.originalorder == 0][0]

        # get a BFS ordering of the edges.  Ignores costs.
        edgesequence = [x for x in nx.bfs_edges(self.joingraph, firstnode)]

        LOG.debug("BFS: edgesequence: %s", edgesequence)

        # Make it deterministic but still in BFS order
        deterministic_edge_sequence = []
        while len(edgesequence) > 0:
            # Consider all edges that have the same first node -- these are all
            # "ties" in BFS order.
            firstx = edgesequence[0][0]
            new_edges = [(x, y) for (x, y) in edgesequence if x == firstx]
            # Sort edges on the originalorder of the source and destination
            deterministic_edge_sequence.extend(
                sorted(new_edges, key=lambda (x, y): (x.originalorder, y.originalorder))
            )  # noqa
            # Remove all those edges from edgesequence
            edgesequence = [(x, y) for (x, y) in edgesequence if x != firstx]

        LOG.debug("BFS: deterministic edge seq: %s", deterministic_edge_sequence)

        # Generate a concrete sequence of terms with conditions properly
        # adjusted
        joinsequence = self.toJoinSequence(deterministic_edge_sequence)
        LOG.debug("BFS: joinsequence: %s", joinsequence)
        return joinsequence

	def __second_order_bad_neighbour(self):  
		
		for node in self.OG.nodes():	
			#calculate the weighted sum of "bad neighbour in" score
			bfs_edge_list = list(nx.bfs_edges(self.OG, node))
			sumedgein = self.__get_sum_of_edge_in(bfs_edge_list)
			sumedgeout = self.__get_sum_of_edge_out(bfs_edge_list)
			total_in_pressure = 0.0
			total_out_pressure = 0.0

			for edge in bfs_edge_list:
				# weight * bni-score 
				# ignore if from node
				try:
					# put a test here??
					bni_score = self.OG.node[edge[1]]['bad_neighbour_in']
					in_weight = self.OG[edge[1]][node]['weight'] / sumedgein
					# get correct sumedge (bfs tree again)
					bni_score = self.__clean_bni(bni_score, node, edge[1])
					total_in_pressure += in_weight * bni_score
					

					bno_score = self.OG.node[edge[1]]['bad_neighbour_out']
					out_weight = self.OG[node][edge[1]]['weight'] / sumedgeout
					total_out_pressure += out_weight * bno_score
				except:
					total_out_pressure += 0
				

			self.OG.node[node]['bad_neighbour_in2'] = total_in_pressure
			self.OG.node[node]['bad_neighbour_out2'] = total_out_pressure

def resolve(relations, entities):
    """Resolve the coreference chains."""
    mentions = graphize(relations, entities)
    for start_entity, _ in mentions.selfloop_edges():  # moet start er nog bij?
        chain = set(entity for edge in nx.bfs_edges(mentions, start_entity)
                    for entity in edge)
        yield chain if chain else set([start_entity])

def get_edge_to_nxfset(T, edge_to_adjacency, root,
        root_prior_fset, node_to_data_fset):
    """
    For each edge, get the joint feasibility of states at edge endpoints.

    Parameters
    ----------
    {params}

    Returns
    -------
    edge_to_nxfset : map from directed edge to networkx DiGraph
        For each directed edge in the rooted tree report the networkx DiGraph
        among states, for which presence/absence of an edge defines the
        posterior feasibility of the corresponding state transition
        along the edge.

    """
    v_to_fset = get_node_to_fset(T, edge_to_adjacency, root,
            root_prior_fset, node_to_data_fset)
    edge_to_nxfset = {}
    for edge in nx.bfs_edges(T, root):
        A = edge_to_adjacency[edge]
        J = nx.DiGraph()
        va, vb = edge
        for sa in v_to_fset[va]:
            sbs = set(A[sa]) & v_to_fset[vb]
            J.add_edges_from((sa, sb) for sb in sbs)
        edge_to_nxfset[edge] = J
    return edge_to_nxfset

	def choose_initial_labeling(self):
		''' Choose a labeling of the local graph that respects the n-hop decision rule'''
		source_neighbors = self.neighbors(self.source)

		# Choose the next virtual source
		weights = []
		for neighbor in source_neighbors:
			subtree = dfs_tree(self, neighbor)
			leaves = [i for i in subtree.nodes() if (subtree.out_degree(i)==0 and \
													 not nx.has_path(subtree,self.source,i))]
			weights.append(len(leaves))
		weights = [float(i)/sum(weights) for i in weights]
		virtual_source = np.random.choice(source_neighbors, p=weights)
		
		# Relabel the nodes so that only the infected nodes have positive values
		label = 1
		mapping = {virtual_source: label}
		# Generate a directed tree emanating from the virtual source away from the source
		directed_self = dfs_tree(self, self.source)
		infected_subgraph = dfs_tree(directed_self, virtual_source)
		bfs_edges = list(nx.bfs_edges(infected_subgraph, virtual_source))
		# In this first timestep, we only need the direct neighbors of V.S. to be infected
		bfs_edges = [item[1] for item in bfs_edges if (item[1] != self.source and item[0] == virtual_source)]
		for edge in bfs_edges:
			label += 1
			mapping[edge] = label
		self = nx.relabel_nodes(self, mapping, copy=False)

def fix_face_winding(mesh):
    '''
    Traverse and change mesh faces in-place to make sure winding is coherent, 
    or that edges on adjacent faces are in opposite directions
    '''
    # we create the face adjacency graph: 
    # every node in g is an index of mesh.faces
    # every edge in g represents two faces which are connected
    graph_all = nx.from_edgelist(mesh.face_adjacency)
    flipped   = 0
    # we are going to traverse the graph using BFS, so we have to start
    # a traversal for every connected component
    for graph in nx.connected_component_subgraphs(graph_all):
        start = graph.nodes()[0]
        # we traverse every pair of faces in the graph
        # we modify mesh.faces and mesh.face_normals in place 
        for face_pair in nx.bfs_edges(graph, start):
            # for each pair of faces, we convert them into edges,
            # find the edge that both faces share, and then see if the edges
            # are reversed in order as you would expect in a well constructed mesh
            pair    = mesh.faces[[face_pair]]            
            edges   = faces_to_edges(pair)
            overlap = group_rows(np.sort(edges,axis=1), require_count=2)
            if len(overlap) == 0:
                # only happens on non-watertight meshes
                continue
            edge_pair = edges[[overlap[0]]]
            if edge_pair[0][0] == edge_pair[1][0]:
                # if the edges aren't reversed, invert the order of one of the faces
                flipped += 1
                mesh.faces[face_pair[1]] = mesh.faces[face_pair[1]][::-1]
    log.info('Flipped %d/%d edges', flipped, len(mesh.faces)*3)