def generate_test_graph(original_graph, state_type_node_count, chosen_state, state_wise_node_list,
                        no_of_types):
    global available_type, node_latitude, node_longitude, node_state, node_type
    test_graph_nodes = list()
    node_queue = list()
    print state_type_node_count
    for i in range(len(chosen_state)):
        for j in range(no_of_types):
            print 'Running for state', chosen_state[i]
            print state_type_node_count[i, j]
            while not state_type_node_count[i, j] == 0:
                node_index = random.randint(1, len(state_wise_node_list[chosen_state[i]]))
                node = state_wise_node_list[chosen_state[i]][node_index-1]
                state_wise_node_list[chosen_state[i]].remove(node)
                if not node_type[node] in available_type[:no_of_types] or node in node_queue or node in test_graph_nodes:
                    continue
                node_queue.append(node)
                state_type_node_count[i, j] -= 1

                expand_graph(node_queue, test_graph_nodes, original_graph, chosen_state, no_of_types, state_type_node_count)

    test_graph = networkx.DiGraph()
    test_graph = original_graph.subgraph(test_graph_nodes)


    components = networkx.weakly_connected_component_subgraphs(test_graph)
    i = 1
    print 'Components Before:'
    print '******************'
    for component in components:
        print 'Component: ' + str(i) + '- ' + str(networkx.number_of_nodes(component))
        i += 1


    # Check connectivity
    if i > 1:
        resolve_connectivity_issue(test_graph)
        components = networkx.weakly_connected_component_subgraphs(test_graph)
        i = 1
        print 'Components After:'
        print '******************'
        for component in components:
            print 'Component: ' + str(i) + '- ' + str(networkx.number_of_nodes(component))
            i += 1


    node_mapping = dict()
    test_graph_node_state_assgn = list()
    for i in range(len(test_graph_nodes)):
        node_mapping[i] = test_graph_nodes[i]
        test_graph_node_state_assgn.append(node_state[test_graph_nodes[i]])
    print 'No of nodes: ', len(test_graph_nodes)

    return test_graph, node_mapping, test_graph_node_state_assgn

def decompose(paths, args):
    """ runs decomposition
    Parameters
    ----------
    paths.bundle_file       : file
    paths.tmp1_file         : file
    paths.tmp2_file         : file
    paths.decomp_file       : file
    args.msize              : integer
    """

    # load the bundle graph.
    logging.info("loading info")
    BG = nx.read_gpickle(paths.bundle_file)
    #BG = test_bi()
    #BG = test_tri()

    # run decomposition until satisfied.
    BG.graph['redo'] = False
    while 1 == 1:

        # decomposition.
        DC = decomp0(BG, paths.tmp1_file, paths.tmp2_file, msize=args.msize)

        # check if only once.
        if args.msize == None or BG.graph['redo'] == False:
            break
        elif BG.graph['redo'] == True:
            BG.graph['redo'] = False

        # remove temp files.
        if os.path.isfile(paths.tmp1_file) == True:
            subprocess.call(["rm","-f",paths.tmp1_file])
        if os.path.isfile(paths.tmp2_file) == True:
            subprocess.call(["rm","-f",paths.tmp2_file])

    # compact decomposition.
    _compact_outter(DC)
    for subcc in nx.weakly_connected_component_subgraphs(DC):

        # call recursive compaction.
        _compact_inner(DC)

    # verify decomposition.
    for subcc in nx.weakly_connected_component_subgraphs(DC):

        # check its consistency.
        _validate_comp(subcc)

    # write to disk.
    nx.write_gpickle(DC, paths.decomp_file)
    nx.write_gpickle(BG, paths.bundle_file)

    def weakly_connected_subgraphs(self):
        """
        Yields weakly connected subgraphs and their topolgical sort.

        """
        for subgraph in nx.weakly_connected_component_subgraphs(self.G):
            yield (subgraph, nx.topological_sort(subgraph))

def longestPath(g, dict_seq):
    paths = []
    if g.number_of_nodes() == 0:
        return paths
    if g.number_of_nodes() == 1:
        paths.append(g.nodes())
        return paths
    if is_linear_graph(g)[0]:
        p = get_path_linear_graph(g)
        return [p]
    for c in nx.weakly_connected_component_subgraphs(g):
        if c.number_of_nodes() == 1:
            paths.append(c.nodes())
            continue
        dist = {}
        for node in nx.topological_sort(c):
            pairs = [(dist[v][0]+len(dict_seq[node])- g[v][node]['weight'], v) for v in c.pred[node]]
            if pairs:
                dist[node] = max(pairs)
            else:
                dist[node] = (len(dict_seq[node]), node)
        node, (length, _) = max(dist.items(), key=lambda x:x[1])
        path = []
        while length > len(dict_seq[node]):
            path.append(node)
            length, node = dist[node]
        paths.append(list(reversed(path)))
    return paths

 def split(self):
     '''splits into weakly connected component subgraphs'''
     # get connected components of graph which represent independent genes
     # unconnected components are considered different genes
     Gsubs = list(nx.weakly_connected_component_subgraphs(self.G))
     if len(Gsubs) == 1:
         yield self
         return
     # map nodes to components
     node_subgraph_map = {}
     subgraph_transfrag_map = collections.defaultdict(list)
     for i, Gsub in enumerate(Gsubs):
         for n_id in Gsub:
             n = self.get_node_interval(n_id)
             node_subgraph_map[n] = i
     # assign transfrags to components
     for t in self.itertransfrags():
         for n in split_transfrag(t, self.node_bounds):
             subgraph_id = node_subgraph_map[n]
             subgraph_transfrag_map[subgraph_id].append(t)
             break
     # create new graphs using the separate components
     for subgraph_transfrags in subgraph_transfrag_map.itervalues():
         yield SpliceGraph.create(subgraph_transfrags,
                                  guided_ends=self.guided_ends,
                                  guided_assembly=self.guided_assembly)

def layer_layout(g, level_attribute = "t"):
    '''Lay out a directed graph by layer
    
    g - a NetworkX directed graph with the layer defined as the node's "t"
        attribute. The graph must be acyclic - a restriction that's guaranteed
        by TrackObjects since edges are always going forward in time.
        
    level_attribute - the attribute in the node attribute dictionary that
        specifies the level of the node
        
    on exit, each node will have a y attribute that can be used to place
    the node vertically on a display. "t" can be used for the horizontal
    display.
    
    The algorithm is a partial implementation of 
    Sugiyama, Kozo, Tagawa, Shojiro; Toda, Mitsuhiko (1981), 
    "Methods for visual understanding of hierarchical system structures", 
    IEEE Transactions on Systems, Man, and Cybernetics SMC-11 (2):109-125,
    doi:10.1109/TSMC.1981.4308636
	
    as described by sydney.edu.au/engineering/it/~visual/comp4048/slides03.ppt
    '''
    
    subgraphs = nx.weakly_connected_component_subgraphs(g)
    y = 0
    for subgraph in subgraphs:
        y = layer_layout_subgraph(g, subgraph, y, level_attribute)

def simplify_graph(g):
    for e in g.selfloop_edges():
        g.remove_edge(e[0], e[1])
    for node in g.nodes():
        neighbors = list(nx.all_neighbors(g, node))
        edges = g.in_edges(node, data=True)
        edges.extend(g.out_edges(node, data=True))
        plus = []
        minus = []
        for e in edges:
            if e[2][node] == '+':
                plus.append(e)
            else:
                minus.append(e)
            if not plus or not minus:
                continue
            if len(plus) >= len(minus):
                for e in minus:
                    if g.has_edge(e[0], e[1]):
                        g.remove_edge(e[0], e[1])
            if len(plus) <= len(minus):
                for e in plus:
                    if g.has_edge(e[0], e[1]):
                        g.remove_edge(e[0], e[1])
    remove_out_tips(g)
    remove_in_tips(g)
    for c in nx.weakly_connected_component_subgraphs(g):
        if c.number_of_nodes() <= 2:
            continue
        isLinear, ends, source, sink= is_linear_graph(c)
        if isLinear:
            if sink==1 and source == 1:
                continue
            adjust_edge_di(g, c, ends[0], ends[1])

def prune_transcript_graph(G, strand, transcript_map,
                           min_trim_length=0, 
                           trim_utr_fraction=0.0,
                           trim_intron_fraction=0.0):
    '''
    trim_utr_fraction: float specifying the fraction of the average UTR
    coverage below which the ends of the UTR will be trimmed

    trim_intron_fraction: float specifying the fraction of the average 
    intron coverage below which intronic nodes will be removed
    '''
    # trim utrs and intron retentions
    trim_nodes = trim_graph(G, strand, min_trim_length, 
                            trim_utr_fraction, 
                            trim_intron_fraction)
    G.remove_nodes_from(trim_nodes)
    # collapse consecutive nodes in graph
    H = collapse_strand_specific_graph(G, transcript_map, introns=True)
    # get connected components of graph which represent independent genes
    # unconnected components are considered different genes
    Gsubs = nx.weakly_connected_component_subgraphs(H)
    for Gsub in Gsubs:
        # get partial path data supporting graph
        transcript_node_map = get_transcript_node_map(Gsub)
        path_score_dict = collections.defaultdict(lambda: 0)
        for t_id, nodes in transcript_node_map.iteritems():
            # reverse path for negative strand transcripts
            if strand == NEG_STRAND:
                nodes.reverse()
            # get transcript scores
            t = transcript_map[t_id]
            path_score_dict[tuple(nodes)] += t.score
        yield Gsub, strand, path_score_dict.items()

def draw_graphs(G, folder_name):
    domain_name = G.graph['domain']
    dir = folder_name + '/' + domain_name
    if not os.path.exists(dir):
        os.makedirs(dir)
    subgraphs = nx.weakly_connected_component_subgraphs(G)
    add_root_to_subgraphs(subgraphs)
    subgraphs.sort(key=lambda subgraph: subgraph.number_of_nodes())
    for i in xrange(len(subgraphs)): 
        subgraph = subgraphs[i]
        pos = nx.spring_layout(subgraph)
        node_labels = get_node_labels(subgraph)
        positive_nodes = node_labels['positive'].keys()
        negative_nodes = node_labels['negative'].keys()
        labels = dict(node_labels['positive'], **(node_labels['negative']))
        edge_labels = get_edge_labels(subgraph)
        pl.figure(figsize=(16, 12))
        nx.draw_networkx_nodes(subgraph, pos, positive_nodes, alpha=0.5, node_color='w')
        nx.draw_networkx_nodes(subgraph, pos, negative_nodes, alpha=0.5, node_color='b')
        nx.draw_networkx_nodes(subgraph, pos, ['root'], node_color='g')
        nx.draw_networkx_edges(subgraph, pos, color='k')
        nx.draw_networkx_labels(subgraph, pos, labels, font_size=20)
        nx.draw_networkx_edge_labels(subgraph, pos, edge_labels, font_size=20)
        pl.axis('off')
        pl.savefig('%s/%s_subgraph_%d.png' % (dir, domain_name, i+1))

    def compute_dependent_cohorts(self, objects, deletion):
        model_map = defaultdict(list)
        n = len(objects)
        r = range(n)
        indexed_objects = zip(r, objects)

        mG = self.model_dependency_graph[deletion]

        oG = DiGraph()

        for i in r:
            oG.add_node(i)

        for v0, v1 in mG.edges():
            try:
                for i0 in range(n):
                   for i1 in range(n):
                       if i0 != i1:
                            if not deletion and self.concrete_path_exists(
                                    objects[i0], objects[i1]):
                                oG.add_edge(i0, i1)
                            elif deletion and self.concrete_path_exists(objects[i1], objects[i0]):
                                oG.add_edge(i0, i1)
            except KeyError:
                pass

        components = weakly_connected_component_subgraphs(oG)
        cohort_indexes = [reversed(topological_sort(g)) for g in components]
        cohorts = [[objects[i] for i in cohort_index]
                   for cohort_index in cohort_indexes]

        return cohorts

def get_alternative_paths(subg,path):
	paths = []
	subg1 = subg.copy()
	for node in path:
		subg1.remove_node(node)

	for comp in nx.weakly_connected_component_subgraphs(subg1):
		if len(comp.nodes()) == 1:
			paths.append(comp.nodes())
		else:
			p = []
			for node in comp.nodes():
				if comp.out_degree(node) == 1 and comp.in_degree(node) == 0:
					p.append(node)
			for node in comp.nodes():
				if comp.out_degree(node) == 0 and comp.in_degree(node) == 1:
					p.append(node)

			if len(p) == 2:
				try:
					paths.append(nx.shortest_path(comp,p[0],p[1]))
				except:
					continue

	return paths

def find_reach_topsort(dags, c2n):
    node_reach = dict()
    cluster_reach = dict()

    wccs = nx.weakly_connected_component_subgraphs(dags)

    for hub in wccs:
        # treat hubs of size 1 and 2 specially
        if len(hub) == 1:
            cluster = hub.nodes()[0]
            cluster_reach[cluster] = c2n[cluster]

            node_reach.update(dict(zip(c2n[cluster], [len(c2n[cluster])]*len(c2n[cluster]))))
        elif len(hub) == 2:
            cluster1, cluster2 = hub.edges()[0]

            cluster_reach[cluster2] = c2n[cluster2]
            cluster_reach[cluster1] = c2n[cluster1] + c2n[cluster2]

            node_reach.update(dict(zip(c2n[cluster1], [len(cluster_reach[cluster1])]*len(c2n[cluster1]))))
            node_reach.update(dict(zip(c2n[cluster2], [len(cluster_reach[cluster2])]*len(c2n[cluster2]))))
        else:
            hub_ts = nx.topological_sort(hub, reverse=True)
            for cluster in hub_ts:
                reach = set()
                for _, out_cluster in dags.out_edges(cluster):
                    reach.update(cluster_reach[out_cluster])
                reach.update(c2n[cluster])
                cluster_reach[cluster] = reach

                node_reach.update(dict(zip(c2n[cluster], [len(reach)]*len(c2n[cluster]))))
    return node_reach

def main():
    file_path = sys.argv[1]
    global user_graph

    # Constructs the graph based on the dataset
    make_graph(file_path)

    # Get the weakly connected graph components. HITS is to be run on the largest of such components.
    weakly_connected_graph_components = nx.weakly_connected_component_subgraphs(user_graph)

    # Get the largest weekly connected graph component
    largest_weakly_connected_graph = weakly_connected_graph_components[0]

    (hub_score_counter, authority_score_counter) = run_hits_algorithm(largest_weakly_connected_graph)

    # Sort the lists
    sorted_hub_score_list = sorted(hub_score_counter.items(), key = lambda item: item[1], reverse = True)
    sorted_authority_score_list = sorted(authority_score_counter.items(), key = lambda item: item[1], reverse = True)

    # Print top 20 hubs
    print "Top 20 Hubs"
    print "==========="
    for i in range(0, 20):
        if sorted_hub_score_list[i] != None:
            print sorted_hub_score_list[i][0]

    print ""

    # Print top 20 authorities
    print "Top 20 Authorities"
    print "=================="
    for i in range(0, 20):
        if sorted_authority_score_list[i] != None:
            print sorted_authority_score_list[i][0]

def __cut(graph):
    ''' param: 
            graph:a nx.DiGraph obj
	    return:
		    cs : edge cut set of the graph
		    g1 , g2 : subgraphs induced by cs
	
    '''
    assert isinstance(graph, nx.DiGraph), "graph class: %s " % graph.__class__
    assert graph.number_of_nodes() > 1,   "Number of nodes: %d" % graph.number_of_nodes()
    unigraph = nx.Graph( graph )          
    cs = nx.minimum_edge_cut( unigraph ) 
    if not cs:
        raise Exception,"Cut Set of this graph is Empty"

    #CS中的边，可能不存在于原来的有向图中，所以需要将这种边的方向颠倒
    #将所有real edge,存到RCS中
    rcs = []
    for eachEdge in cs:
        if not graph.has_edge( eachEdge[0], eachEdge[1] ):
            eachEdge = (eachEdge[1], eachEdge[0]) #调换方向
        rcs.append(eachEdge)
    graph.remove_edges_from(rcs)
    glist = []
    for eachCntComp in nx.weakly_connected_component_subgraphs(graph, copy = False):
        glist.append(eachCntComp)
    assert len(glist) == 2
    return rcs, glist[0], glist[1]

def connect_digraph(D):
	""" Take a DiGraph with weakly connected components, and coalesce into one component."""

	s = nx.weakly_connected_component_subgraphs(D)
	#s is sorted by the size of the subgraph

	if len(s) > 1:
		largest = s[0]
		remaining = s[1:]

		largest_edges = largest.edges()

		#Let's filter out the one degree edges (otherwise we'll disconnect 
		#the graph when we swap edges around).
		candidates = []
		for u,v in largest_edges:
			if D.degree(u) > 1 and D.degree(v) > 1:
				candidates.append((u,v))

		if len(candidates) < len(remaining):
			raise Exception("There are not enough candidates for swapping.")

		#Connect the largest subgraph to the remaining.
		for G in remaining:
			u,v = random.choice(candidates)
			x,y = random.choice(G.edges())

			D.remove_edge(u, v)
			D.remove_edge(x, y)
			D.add_edge(u, y)
			D.add_edge(x, v)

			largest_edges.remove((u,v))

def main():
    G, karmas = read_data("karma.txt")
    cs = nx.weakly_connected_component_subgraphs(G)
    cs.sort(key=lambda c: c.number_of_nodes(), reverse=True)

    plt.clf()
    draw(cs[126], karmas)
    plt.show()

 def __init__(self, scaffold_graph):
     print "Entering PathFinder module:", str(datetime.now())
     self.G = scaffold_graph.copy()
     #Build strandless list of sequences
     sequences = set([n for n in self.G.nodes() if n > 0])
     #Define weakly connected components
     print "1... Defining weakly connected components"
     component_graphs = set([g for g in nx.weakly_connected_component_subgraphs(self.G)])
     single_node_graphs = set([g for g in component_graphs if len(g.nodes()) == 1])
     multi_node_graphs = set([g for g in component_graphs if len(g.nodes()) > 1])
     print "Number of single-node components:", len(single_node_graphs)
     print "Number of multi-node components:", len(multi_node_graphs)
     #Consolidate unscaffolded nodes, discard reverse strand
     print "2... Consolidating single-node components"
     unscaffolded = set([g.nodes()[0] for g in single_node_graphs])
     discard_nodes = set([n for n in unscaffolded if n < 0])
     for g in iter(single_node_graphs.copy()):
         if g.nodes()[0] in discard_nodes:
             single_node_graphs.discard(g)
     print "Number of unscaffolded sequences:", len(single_node_graphs)
     #Classify multi-node graphs
     print "3... Classifying multi-node components"
     DAG = set([])
     Euler = set([])
     for g in multi_node_graphs:
         if nx.is_directed_acyclic_graph(g):
             DAG.add(g)
         elif nx.is_eulerian(g):
             Euler.add(g)
         else:
             sys.exit("FATAL ERROR: Unknown multi-node graph type!")
     print "Number of directed acyclic graphs:",  len(DAG)
     print "Number of Eulerian graphs:", len(Euler)
     #Build scaffolds from DAGs
     print "4... Building scaffolds from directed acyclic graphs"
     self.scaffolds = set([])
     for g in DAG:
         self.build_dag_scaffold(g)
     #Consolidating complementary scaffolds, keep first found
     print "5... Consolidating complementary scaffolds"
     consolidated_scaff = set([])
     for seq in iter(self.scaffolds):
         comp = self.revc(seq)
         if comp in self.scaffolds:
             if comp not in consolidated_scaff:
                 consolidated_scaff.add(seq)
         else:
             print "WARNING: non-complemented scaffold"
     self.scaffolds = consolidated_scaff
     print "Number of scaffolds assembled:", len(self.scaffolds)
     #Build scaffolds from Eulerian graphs
     
     #Add unscaffolded seqs to scaffolds list
     print "6... Adding unscaffolded sequences to output"
     for g in single_node_graphs:
         seq = self.G.node[g.nodes()[0]]['seq']
         self.scaffolds.add(seq)
     print "Leaving PathFinder module:", str(datetime.now())

    def find_largest_component(self):
	G = self.graph
	list_Graphs = nx.weakly_connected_component_subgraphs(G)
	max_component = list_Graphs[0]
	for g in list_Graphs:
    	    if nx.number_of_nodes(g) > nx.number_of_nodes(max_component):
                max_component = g

        return  max_component

    def keep_weakly_connected(self):
        '''This method filters out exons (nodes) not involved in AS events'''
        # find weakly connected subgraphs
        weakly_connected_list = nx.weakly_connected_component_subgraphs(self.sub_graph)

        # iterate to find which subgraph has the target exon
        for subgraph in weakly_connected_list:
            if self.target in subgraph.nodes():
                self.sub_graph = subgraph  # assign subgraph that actually connects to target exon

def check_connected_balanced(graph):
    """
    :type graph: nx.DiGraph
    """
    for v in graph.nodes():
        assert graph.in_degree(v) == graph.out_degree(v)
    sub_graph_list = nx.weakly_connected_component_subgraphs(graph, True)
    for sub_graph in sub_graph_list:
        print 'connected component:', sub_graph.edges()
    print