def get_graph(filename, with_root=False):
    DG = nx.DiGraph()
    f = open(filename, 'r')
    line = None
    edges = []
    coordinates = []
    terms = []
    if with_root:
        root = None
    while line != 'EOF':
        line = f.readline().strip()
        toks = line.split(' ')
        if toks[0] == 'A':
            t = tuple(int(x) for x in toks[1:])
            edges.append(t)
        if toks[0] == 'T':
            terms.append(int(toks[1]))
        if toks[0] == 'Root':
            if with_root:
                root = int(toks[1])
        if toks[0] == 'DD':
            t = tuple(int(x) for x in toks[1:])
            coordinates.append(t)
    for coord in coordinates:
        DG.add_node(coord[0], pos=(coord[1], coord[2]))
    terms.sort()
    DG.add_weighted_edges_from(edges)
    # print_graph(DG)
    # nx.draw(DG, node_size=50)
    # plt.show()
    # f.close()
    if with_root:
        return DG, terms, root
    else:
        print_graph(DG)
        max_len = 0
        max_node = None
        for node in nx.nodes(DG):
            # print(node, tr_cl.out_edges(node))
            descs = nx.descendants(DG, node)
            # desc_numb = len(descs)
            if len(set(terms) & set(descs)) == len(descs):
                # max_len = desc_numb
                max_node = node
        if max_len == len(nx.nodes(DG)):
            return DG, terms, max_node
        else:
            reachable = set(nx.descendants(DG, max_node)) | {max_node}
            unreachable = set(nx.nodes(DG)) - reachable
            for node in unreachable:
                DG.remove_node(node)
        terms = list(set(terms) & reachable)
        print('terms =', len(terms))
        return DG, terms, max_node

	def descendants(self, nbunch):
		self.validate_input_nodes(nbunch)
		if not self.acceptable_iterable(nbunch):	#single input node
			return nx.descendants(self, nbunch)
		else:
			if len(nbunch) == 1:	#still a single node
				return nx.descendants(self, nbunch[0])
			else:	#multiple input nodes
				DG = self.copy()
				s = DG.add_node_unique()
				for node in nbunch:
					DG.add_edge(s, node) # this automatically adds s to DG too
				return nx.descendants(DG, s) - set(nbunch) # returns a SET

    def forward_reachable(self, state):
        """Return states reachable from given state.

        Iterated post(), a wrapper of networkx.descendants.
        """
        descendants = nx.descendants(self, state)
        return descendants

def print_impacting_modules(single_node=None, json_out=None):
    """
    For each module, print a list of modules that the module is depending on,
    i.e. modules whose change can potentially impact the module. The function
    shows all levels of dependency, not just the immediately imported
    modules.  If the json_out argument is not None, then the output will be
    recorded there instead of on stdout.
    :return:
    """
    if json_out is None:
        print('\n===Impacting Modules===')
    else:
        json_out['impacting_modules'] = {}
    for node_name in G.nodes_iter():
        if single_node and (node_name!=single_node):
            continue
        descendants = nx.descendants(G, node_name)
        if json_out is None:
            print(augment_format_string(node_name, '\n%s:') % node_name)
        else:
            json_out['impacting_modules'][node_name] = []
        for d in descendants:
            if json_out is None:
                print(augment_format_string(d, '    %s') % d)
            else:
                json_out['impacting_modules'][node_name].append(d)

    def _resolve_update_list(self, changed_properties):
        """
        Returns a list of all plasma models which are affected by the
        changed_modules due to there dependency in the
        the plasma_graph.

        Parameters
        ----------

        changed_modules: ~list
            all modules changed in the plasma

        Returns
        -------

            : ~list
            all affected modules.
        """

        descendants_ob = []

        for plasma_property in changed_properties:
            node_name = self.outputs_dict[plasma_property].name
            descendants_ob += nx.descendants(self.graph, node_name)

        descendants_ob = list(set(descendants_ob))
        sort_order = nx.topological_sort(self.graph)

        descendants_ob.sort(key=lambda val: sort_order.index(val))

        logger.debug("Updating modules in the following order:".format("->".join(descendants_ob)))

        return descendants_ob

def modify_downstream_edges_faster(G,source,modified_edges,time_to_solve,og_delay):
#  downstream_nodes = nx.descendants(G,source)
#  for node in downstream_nodes:
#  
#    #Getting incoming edges to this node.
#    in_edges = G.in_edges(node)
#    #Get the weights of  in_edges.
#    weights = [z for x,y,z in in_edges]
#    #The maximum weight (which is when this downstream node is ready to solve)
#    ready_to_solve = copy(max(weights))
#    
#    for u,v in in_edges:
#      if (u == source or u in downstream_nodes):
#        if not modified_edges:
#          G[u][v]['weight'] += delay
#          modified_edges.append((u,v))
#        elif (u,v) not in modified_edges:
#          G[u][v]['weight'] += delay
#          modified_edges.append((u,v))
  
  
  downstream_nodes = list(nx.descendants(G,source))
  #Add the source node to the downstream nodes.
  downstream_nodes = [source] + downstream_nodes
  num_downstream_nodes = len(downstream_nodes)
  #We get when each downstream node is ready to solve.
  ready_to_solve_all = {}
  for n in range(0,num_downstream_nodes):
    current_node = downstream_nodes[n]
    #Get incoming edge with the maximum weight to this node.
    ready_to_solve_all[current_node] = get_max_incoming_weight(G,current_node)
  
  
  #Sorting the downstream nodes in order of when they solve.
  ready_to_solve_all = dict(sorted(ready_to_solve_all.items(),key=lambda x:x[1]))
  
  
  for k,val in ready_to_solve_all.items():
    #The current node.
    node = k
    #When the current node is ready to solve.
    ready_to_solve = val
    #Get outgoing edges of this node.
    out_edges = G.out_edges(node)
    for u,v in out_edges:
      if (v in downstream_nodes):
        if not modified_edges:
          delay = time_to_solve[node] + ready_to_solve - G[u][v]['weight']
          if delay > 0.0:
            G[u][v]['weight'] += delay
            ready_to_solve_all[v] = get_max_incoming_weight(G,v)
            #modified_edges.append((u,v))
        elif (u,v) not in modified_edges:
          delay = time_to_solve[node] + ready_to_solve - G[u][v]['weight']
          if delay > 0.0:
            G[u][v]['weight'] += delay
            ready_to_solve_all[v] = get_max_incoming_weight(G,v)
            #modified_edges.append((u,v))
  
  return G

def out_component(G, source):
    '''rather than following the pseudocode in figure 6.15 of Kiss, Miller & Simon, this uses a built in networkx command.  I plan to improve this algorithm.

    finds the set of nodes (including source) which are reachable from nodes in source.

    Parameters
    ----------
    G : NetworkX Graph
        The network the disease will transmit through.
    source : either a node or an iterable of nodes (set, list, tuple)
        The nodes from which the infections start.

    Returns
    -------
    reachable_nodes : set
        the set of nodes reachable from source (including source).
    '''
    try:
        #testing whether this is an iterable
        iterator = iter(source)
    except TypeError:
        #It's not an iterable.  It "must" be a node.
        if G.has_node(source):
            source_nodes = set([source])
    else:
        #it's an iterable.  
        source_nodes = set(source)
    reachable_nodes = set([])
    for node in source_nodes:
        reachable_nodes = reachable_nodes.union(set(nx.descendants(G, node)))
    return reachable_nodes

    def as_dependency_list(self, limit_to=None):
        """returns a list of list of nodes, eg. [[0,1], [2], [4,5,6]]. Each element contains nodes whose
        dependenices are subsumed by the union of all lists before it. In this way, all nodes in list `i`
        can be run simultaneously assuming that all lists before list `i` have been completed"""

        if limit_to is None:
            graph_nodes = set(self.graph.nodes())
        else:
            graph_nodes = set()
            for node in limit_to:
                graph_nodes.add(node)
                if node in self.graph:
                    graph_nodes.update(nx.descendants(self.graph, node))
                else:
                    raise RuntimeError("Couldn't find model '{}' -- does it exist or is it diabled?".format(node))

        depth_nodes = defaultdict(list)

        for node in graph_nodes:
            num_ancestors = len(nx.ancestors(self.graph, node))
            depth_nodes[num_ancestors].append(node)

        dependency_list = []
        for depth in sorted(depth_nodes.keys()):
            dependency_list.append(depth_nodes[depth])

        return dependency_list

    def _calculate_scores(self):
        """Calculate the 'value' of each node in the graph based on how many
        blocking descendants it has. We use this score for the internal
        priority queue's ordering, so the quality of this metric is important.

        The score is stored as a negative number because the internal
        PriorityQueue picks lowest values first.

        We could do this in one pass over the graph instead of len(self.graph)
        passes but this is easy. For large graphs this may hurt performance.

        This operates on the graph, so it would require a lock if called from
        outside __init__.

        :return Dict[str, int]: The score dict, mapping unique IDs to integer
            scores. Lower scores are higher priority.
        """
        scores = {}
        for node in self.graph.nodes():
            score = -1 * len([
                d for d in nx.descendants(self.graph, node)
                if self._include_in_cost(d)
            ])
            scores[node] = score
        return scores

 def filter_graph(graph):
     from_s = nx.descendants(graph, start_id)
     from_s.add(start_id)
     to_e = nx.ancestors(graph, end_id)
     to_e.add(end_id)
     del_cross = (from_s | to_e) - (from_s & to_e)
     graph.remove_nodes_from(del_cross)

 def reset(self, cell):
     if cell.value is None:
         return
     cell.value = None
     for descendant in descendants(self.graph, cell):
         if isinstance(descendant, CellRange) or descendant.formula:
             descendant.value = None

def _determine_t_death(tree, target):

    # find the time of parent and the distance from it
    parent = tree.predecessors(target)[0]

    start_dist = tree.edge[parent][target]['distance']
    start_time = tree.node[parent]['t_death']

    # build list of descendants within the same species
    descendants = [n for n in nx.descendants(tree, target)
                   if tree.node[n]['S'] == tree.node[target]['S']]

    # find the most distant descendant with 't_death' label
    distances_times = []
    for node in descendants:
        distance_time = (
            nx.shortest_path_length(tree, source=target, target=node, weight='distance'),
            tree.node[node].get('t_death', None)
        )
        distances_times.append(distance_time)

    # max_dist = max(distances)
    distances_times.sort(key=lambda x: x[0])

    end_dist, end_time = distances_times[-1]

    # t_death for node is between that of parent and descendant
    # proportionate to the distance to each
    t_death = start_time + (end_time - start_time) * (start_dist / (start_dist + end_dist))

    tree.node[target]['t_death'] = t_death

    def subgraph_needed_for(self, start_at, end_at):
        """Find the subgraph of all dependencies to run these tasks. Returns a
        new graph.
        """
        assert start_at or end_at, "one of {start_at,end_at} must be a task id"
        start, end = map(self.task_dict.get, [start_at, end_at])
        if None in [start, end]:
            graph = self.get_networkx_graph()
            if start:
                task_subset = nx.descendants(graph, start)
                task_subset.add(start)
            elif end:
                task_subset = nx.ancestors(graph, end)
                task_subset.add(end)
        elif start == end:
            task_subset = set([start])
        else:
            graph = self.get_networkx_graph()
            task_subset = set()
            for path in nx.all_simple_paths(graph, start, end):
                task_subset.update(path)

        # make sure the tasks are added to the subgraph in the same
        # order as the original configuration file
        tasks_kwargs_list = [task.yaml_data for task in self.task_list
                             if task in task_subset]
        subgraph = TaskGraph(self.config_path, tasks_kwargs_list)
        return subgraph

def subtree(G, node):
    GS = G.copy()
    GS.remove_node(node)
    sd = nx.descendants(G, node)
    sd.add(node)
    s = set(sd)
    S = G.subgraph(s).copy()

    for n in sd:
        if n == node:
            continue
        ns = nx.ancestors(GS, n)
        if not ns.issubset(sd):
            S.remove_node(n)
            s.discard(n)

    pn = set(G.predecessors_iter(node))

    gs = set(
        itertools.chain.from_iterable(
            nx.shortest_path(G, "REPO", n) for n in pn
        ))

    GS = G.subgraph(gs.union(s)).copy()

    for n in pn.difference(s):
        GS.node[n]["fontcolor"] = "#FF0000"

    for n in s:
        GS.node[n]["fontcolor"] = "#006600"

    GS.remove_node("REPO")

    return S, GS

    def parents(self, gid):
        """Return direct asscendants in the hierarchy for this GeoName ID.

        If the location has not parents in the hierarchy it will attempt to
        find them nonetheless using the following algorithm:

        1. Find all descendants
        2. Find the 1000 nearest locations, if any of them has the same name
           or has more population and it's not a descendant then it's the
           new parent.

        The descendants check is to avoid loops in the hierarchy.
        """
        try:
            p = self._hierarchy.predecessors(gid)
        except nx.NetworkXError:
            p = []

        if not p and gid not in self._root:
            name = self.name(gid)
            population = self.population(gid)
            try:
                descendants = nx.descendants(self._hierarchy, gid)
            except nx.NetworkXError:
                descendants = set()
            for neighbor in self.nearest(gid, 1000):
                match_name = self.name(neighbor) == name
                bigger = (population > 0) and (self.population(neighbor) > population)
                if ((match_name or bigger) and (neighbor not in descendants)):
                    p.append(neighbor)
                    self._hierarchy.add_edge(neighbor, gid)
                    break
            if not p:
                self._root.add(gid)
        return p

    def OnClick(self, node_id):
        self.color_nodes()
        self._current_node_id = node_id
        node_ea = self[node_id]

        self._remove_target_handler.unregister()
        self._disable_source_handler.unregister()
        self._enable_source_handler.unregister()

        if node_ea in self._targets:
            self._remove_target_handler.register()
            self._attach_to_popup(self._remove_target_handler.get_name())

            for ea in nx.ancestors(self._lca_graph, node_ea):
                if ea not in self._targets and ea not in self._sources:
                    self._set_node_bg_color(self._node_ids[ea], COLOR_PATH)

        if node_ea in self._sources:
            if node_ea in self._disabled_sources:
                self._enable_source_handler.register()
                self._attach_to_popup(self._enable_source_handler.get_name())
            else:
                self._disable_source_handler.register()
                self._attach_to_popup(self._disable_source_handler.get_name())

                for ea in nx.descendants(self._lca_graph, node_ea):
                    if ea not in self._targets and ea not in self._sources:
                        self._set_node_bg_color(self._node_ids[ea], COLOR_PATH)

        return False

def mark_reachable_nodes(ea, source_color=COLOR_SOURCE, other_color=COLOR_REACHABLE):
    graph = get_nx_graph(ea)
    block_ea = get_block_start(ea)
    for descendant in nx.descendants(graph, block_ea):
        CodeBlock(descendant).color = other_color

    CodeBlock(ea).color = source_color

def mark_reaching_nodes(ea, source_color=COLOR_SOURCE, other_color=COLOR_REACHING):
    graph = get_nx_graph(ea)
    graph = graph.reverse()
    block_ea = get_block_start(ea)
    for descendant in nx.descendants(graph, block_ea):
        CodeBlock(descendant).color = other_color

    CodeBlock(ea).color = source_color

def getReachability(G):
	numNodes = len(G.nodes())
	# Obtains all the reability values for the nodes in the graph
	reaches = list(imap(lambda node: len(nx.descendants(G, 
		node))/numNodes, G.nodes()))
	avgReach = mean(reaches)
	indReache = dict(izip(G.nodes(), reaches))
	return avgReach

def getcentral(g1):
    # get different centrality
    return pd.DataFrame({
        u'anc': {x: len(nx.ancestors(g1, x)) for x in g1.nodes()},
        u'des': {x: len(nx.descendants(g1, x)) for x in g1.nodes()},
        u'indeg': g1.in_degree(),
        u'outdeg': g1.out_degree()
        })