def connect_module_graph(G, outdegree_list):
    """Connect disconnected modules. Note: This function cannot be used to 
    connect the entire modular graph."""
    cc_tot = list(nx.connected_components(G))  # cc returns the connected components of G as lists cc[0], cc[1], etc.
    isolated_comp, outedge_comp, isolated_comp_count, outedge_comp_count = partition_network_components(
        cc_tot, outdegree_list
    )

    while isolated_comp_count > 0:  # while G is not connected, reduce number of components
        # pick a random node in the largest component cc[0] that has degree > 1
        node1 = rnd.choice(isolated_comp[0])
        # pick a node in another component whose degree >1
        node2 = rnd.choice(outedge_comp[rnd.choice([x for x in xrange(outedge_comp_count)])])
        while G.degree(node2) <= 1:
            node2 = rnd.choice(outedge_comp[rnd.choice([x for x in xrange(outedge_comp_count)])])

        # pick neighbors of node1 and node2
        nbr1 = rnd.choice(G.neighbors(node1))
        nbr2 = rnd.choice(G.neighbors(node2))

        # swap connections between node1,nbr1 with connections between node2,nbr2
        #  to attempt to connect the two components
        G.remove_edges_from([(node1, nbr1), (node2, nbr2)])
        G.add_edges_from([(node1, node2), (nbr1, nbr2)])

        cc_tot = list(nx.connected_components(G))
        isolated_comp, outedge_comp, isolated_comp_count, outedge_comp_count = partition_network_components(
            cc_tot, outdegree_list
        )

    def _reduce_graph(self, graph, min0list):
        """determine how much of the graph to include in the disconnectivity graph
        """
        used_nodes = []
        # make sure we include the subgraph containing min0
        if len(min0list) > 0:
            for min0 in min0list:
                nodes = nx.node_connected_component(graph, min0)
                if len(nodes) > 2:
                    used_nodes += nodes
                else:
                    print("dgraph: too few nodes connected to", min0)  
        if len(used_nodes) == 0: 
            # use the biggest connected cluster
            cc = sorted(nx.connected_components(graph), key=len, reverse=True)
            used_nodes += cc[0]  # list is ordered by size of cluster

        if self.subgraph_size is not None:
            node_lists = nx.connected_components(graph)
            for nodes in node_lists:
                if len(nodes) >= self.subgraph_size:
                    used_nodes += nodes

        newgraph = graph.subgraph(used_nodes).copy()
        return newgraph

def runGirvanNewman(G, Orig_deg, m_):
    #let's find the best split of the graph
    BestQ = 0.0
    Q = 0.0
    print "runGirvanNewman"
    while True:    
        CmtyGirvanNewmanStep(G)
        Q = _GirvanNewmanGetModularity(G, Orig_deg, m_);
        print "current modularity: %f" % Q
        if len(nx.connected_components(G)) >= 10 or Q >= 0.5:
            break;
        if Q > BestQ:
            BestQ = Q
            Bestcomps = nx.connected_components(G)    #Best Split
 #           print "comps:"
 #           print Bestcomps
        if G.number_of_edges() == 0:
            break
    if BestQ > 0.0:
#        print "Best Q: %f" % BestQ
        result_data = {};
        result_data['num_clusters'] = len(Bestcomps)
        result_data['list'] = Bestcomps
        return result_data
    else:
#        print "Best Q: %f" % BestQ
        result_data = {};
        result_data['num_clusters'] = len(nx.connected_components(G))
        result_data['list'] = nx.connected_components(G)
        return result_data

def fix_connectivity_of_network(initial_g, pos, tx, threshold=80.0):
    g = initial_g
    x = [(len(c), c) for c in nx.connected_components(g)]
    maximum = max(c[0] for c in x)
    if maximum < len(pos)*threshold/100:
        return False
    while len(x) > 1:
        x.sort(lambda a, b: -1 if a[0] < b[0] else 1)
        sizes = [c[0] for c in x]
        # print len(x), " components with sizes: ", sizes
        idx_c = 0
        c = x[0]
        # print "\thola: ", c[1]
        ttt = [ find_closest(pos, x, idx, idx_c) for idx in c[1] ]
        ttt.sort(lambda a, b: -1 if a[0] < b[0] else 1)
        pair = ttt[0]
        # print "\t\tminimum distance : ", pair[0] - tx, ", to node:", pair[1]
        toward = pair[1]
        x0 = pos[pair[2]][0]
        y0 = pos[pair[2]][1]
        x1 = pos[toward][0]
        y1 = pos[toward][1]
        d = pair[0] - tx + 1
        px = (x1 - x0)/d
        py = (y1 - y0)/d
        # move all nodes in this component
        for idx in c[1]:
            pos[idx][0] += px
            pos[idx][1] += py
        g = build_graph(pos, tx)
        x = [(len(c), c) for c in nx.connected_components(g)]
    return True

	def conn_comps(self, g2):
		# then try to do connectedness
		num_e = len(g2.edges())
		conn_comps = networkx.connected_components(g2)
		conn_trials = 0
		while len(conn_comps) > 1:
			small_comps = conn_comps[1:]
			small_comps.reverse()   # start with smallest components, easiest to fix (?)
			for comp in small_comps:
				edges = []
				for u in comp:
					i = self.state.get_partition_of_node(u)
					nodes_i = self.state.get_nodes(i)
					# sample by edges
					for v in nodes_i:
						if u != v and g2.degree(v) >= 2:
							edges += [ (u, v, w) for w in g2.neighbors(v) \
							if (w != u and not g2.has_edge(u,w)) ]

				if len(edges) == 0:
					continue
				(u, v, w) = random.choice(edges)
				g2.remove_edge(v,w)
				assert u != w
				assert not g2.has_edge(u,w)
				g2.add_edge(u,w)
				assert g2.has_edge(u,w)
				assert num_e == len(g2.edges()), "expected %d vs. actual %d: u=%s v=%s w=%s" % (num_e, len(g2.edges()), str(u), str(v), str(w))
			conn_comps = networkx.connected_components(g2)
			conn_trials += 1
			if conn_trials > 50:
				print "After %d tries, could not connect graph, leaving %d components" % (conn_trials, len(conn_comps))
				break

	def communities(self, nCommunities, weight=None):
		"""
		Compute communities.

		Parameters
		----------
		nCommunities - number of communities to be returned.
			This is added to simplify the process, the original GN algorithm doesn't 
			need predecided number of communities. 
			Other measures like a threshold on betweenness centrality can be used instead.
		
		weight (string) - If None, all edge weights are considered equal. 
			Otherwise holds the name of the edge attribute used as weight. 


		Returns
		--------
		A list of communities where each community is a list of the nodes in the community.	 
		"""
		gr = self.g
		n = nx.number_connected_components(gr)
		components = nx.connected_components(gr)

		while (n < nCommunities):
			gr = self.communitySplits(gr, weight=weight)
			components = nx.connected_components(gr)
			n = nx.number_connected_components(gr)
			if gr.number_of_edges() == 0:
				break
		return components

def select_clusters_for_length(selected_protein_pairs, genome_graph, genome_graph_query):
	# defines the shorteste length of a synteny otholog cluster
	pairs={}
	for n in selected_protein_pairs:	
		pairs[n[0]]=n[1]


	valid_graph=nx.Graph()
	edges=list()
	for edge in genome_graph.edges():


		if edge[0] in pairs.keys() and edge[1] in pairs.keys():
			valid_graph.add_edge(edge[0],edge[1])

	for cluster in nx.connected_components(valid_graph):
		if len(cluster)>4:
			query_nodes={}
			valid_query_graph=nx.Graph()
			for node in cluster:
				
				query_nodes[pairs[node]]=node


			for query_edge in genome_graph_query.edges():
				if query_edge[0] in query_nodes.keys() and query_edge[1] in query_nodes.keys():
					
					valid_query_graph.add_edge(query_edge[0],query_edge[1])
			
			for query_cluster in nx.connected_components(valid_query_graph):
				if len(query_cluster)>4:
					for query_node in sorted(query_cluster):
						print query_nodes[query_node]+"\t"+query_node

def connect_simple_graph(G):
    """check if simple graph G is disconnected and connect if necessary"""

    cc = list(nx.connected_components(G))  # cc returns the connected components of G as lists cc[0], cc[1], etc
    component_count = len(cc)
    while component_count > 1:  # while G is not connected, reduce number of components

        # pick a random node in the largest component cc[0] that has degree > 1
        node1 = rnd.choice(cc[0])
        while G.degree(node1) == 1:
            node1 = rnd.choice(cc[0])

            # pick a node in another component
        node2 = rnd.choice(cc[1])

        # pick neighbors of node1 and node2
        nbr1 = rnd.choice(G.neighbors(node1))
        nbr2 = rnd.choice(G.neighbors(node2))

        # swap connections between node1,nbr1 with connections between node2,nbr2
        #  to attempt to connect the two components
        G.remove_edges_from([(node1, nbr1), (node2, nbr2)])
        G.add_edges_from([(node1, node2), (nbr1, nbr2)])

        cc = list(nx.connected_components(G))
        component_count = len(cc)

def compareGraphs(g1, g2):
    
    """#Compares the quantitative properties of two graph. So I can check the coarse graining. """

    
    #Nodes and edges
    print 'Graph1: #(Nodes, Edges) = (' + str(len(g1.nodes())) + ', ' + str(len(g1.edges())) + ')'
    print 'Graph2: #(Nodes, Edges) = (' + str(len(g2.nodes())) + ', ' + str(len(g2.edges())) + ')'

    #Connected Components
    #print '\n#CCs for graph 1: ' + str(len(nx.connected_components(g1)))
    #print '#CCs for graph 2: ' + str(len(nx.connected_components(g2)))
    
    plt.hist([len(i) for i in nx.connected_components(g1)])
    plt.hist([len(i) for i in nx.connected_components(g2)])
    plt.title('Cluster Size')
    plt.xlabel('Cluster Size')
    plt.ylabel('#Cluster')
    show()
    
    #Degree Distribution
    plt.hist(nx.degree_histogram(g1))
    plt.hist(nx.degree_histogram(g2))
    plt.title('Degree Distribution' )
    plt.xlabel('Degree')
    plt.ylabel('#Nodes')
    show()
    
    #Betweeness --- this is by far the most compuationally demanding.
    plt.hist(nx.betweenness_centrality(g1, normalized = False).values())
    plt.hist(nx.betweenness_centrality(g2, normalized = False).values())
    plt.title('Distribution of Betweenness' )
    plt.xlabel('Betweenness')
    plt.ylabel('#Nodes')
    show()        

def runGirvanNewman(G, Orig_deg, m_):
    #let's find the best split of the graph    
    BestQ = 0.0
    Q = 0.0
    Bestcomps = list(nx.connected_components(G))
    while True:    
        CmtyGirvanNewmanStep(G)
        Q = _GirvanNewmanGetModularity(G, Orig_deg, m_);
        #print "current modularity: %f" % Q
        if Q > BestQ:
            BestQ = Q
            Bestcomps = list(nx.connected_components(G))    #Best Split
            #print "comps:" 
            #print Bestcomps
        if G.number_of_edges() == 0:
            #print "at last use this break"
            break
    """
    if BestQ > 0.0:
        print "Best Q: %f" % BestQ
        print Bestcomps
    else:
        print "Best Q: %f" % BestQ
    """
    return Bestcomps

def betweenness_fracture(infile, outfile, fraction, recalculate = False):
    """
    Removes given fraction of nodes from infile network in reverse order of 
    betweenness centrality (with or without recalculation of centrality values 
    after each node removal) and saves the network in outfile.
    """

    g = networkx.read_gml(infile)
    m = networkx.betweenness_centrality(g)
    l = sorted(m.items(), key = operator.itemgetter(1), reverse = True)
    largest_component = max(networkx.connected_components(g), key = len)
    n = len(g.nodes())
    for i in range(1, n):
        g.remove_node(l.pop(0)[0])
        if recalculate:
            m = networkx.betweenness_centrality(g)
            l = sorted(m.items(), key = operator.itemgetter(1), 
                       reverse = True)
        largest_component = max(networkx.connected_components(g), key = len)
        if i * 1. / n >= fraction:
            break
    components = networkx.connected_components(g)
    component_id = 1
    for component in components:
        for node in component:
            g.node[node]["component"] = component_id
        component_id += 1
    networkx.write_gml(g, outfile)

def closeness(infile, recalculate = False):
    """
    Performs robustness analysis based on closeness centrality,  
    on the network specified by infile using sequential (recalculate = True) 
    or simultaneous (recalculate = False) approach. Returns a list 
    with fraction of nodes removed, a list with the corresponding sizes of 
    the largest component of the network, and the overall vulnerability 
    of the network.
    """

    g = networkx.read_gml(infile)
    m = networkx.closeness_centrality(g)
    l = sorted(m.items(), key = operator.itemgetter(1), reverse = True)
    x = []
    y = []
    largest_component = max(networkx.connected_components(g), key = len)
    n = len(g.nodes())
    x.append(0)
    y.append(len(largest_component) * 1. / n)
    R = 0.0
    for i in range(1, n):
        g.remove_node(l.pop(0)[0])
        if recalculate:
            m = networkx.closeness_centrality(g)
            l = sorted(m.items(), key = operator.itemgetter(1), 
                       reverse = True)
        largest_component = max(networkx.connected_components(g), key = len)
        x.append(i * 1. / n)
        R += len(largest_component) * 1. / n
        y.append(len(largest_component) * 1. / n)
    return x, y, 0.5 - R / n

def eigenvector(g, recalculate=False):
    """
    Performs robustness analysis based on eigenvector centrality,  
    on the network specified by infile using sequential (recalculate = True) 
    or simultaneous (recalculate = False) approach. Returns a list 
    with fraction of nodes removed, a list with the corresponding sizes of 
    the largest component of the network, and the overall vulnerability 
    of the network.
    """

    m = networkx.eigenvector_centrality(g, max_iter=5000)
    l = sorted(m.items(), key=operator.itemgetter(1), reverse=True)
    x = []
    y = []
    largest_component = max(networkx.connected_components(g), key=len)
    n = len(g.nodes())
    x.append(0)
    y.append(len(largest_component) * 1. / n)
    r = 0.0
    for i in range(1, n - 1):
        g.remove_node(l.pop(0)[0])
        if recalculate:

            try:
                m = networkx.eigenvector_centrality(g, max_iter=5000)
            except networkx.NetworkXError:
                break

            l = sorted(m.items(), key=operator.itemgetter(1),
                       reverse=True)
        largest_component = max(networkx.connected_components(g), key=len)
        x.append(i * 1. / n)
        r += len(largest_component) * 1. / n
        y.append(len(largest_component) * 1. / n)
    return x, y, r / n

def rand(infile):
    """
    Performs robustness analysis based on random attack, on the network 
    specified by infile. Returns a list with fraction of nodes removed, a 
    list with the corresponding sizes of the largest component of the 
    network, and the overall vulnerability of the network.
    """

    g = networkx.read_gml(infile)
    l = [(node, 0) for node in g.nodes()]
    random.shuffle(l)
    x = []
    y = []
    largest_component = max(networkx.connected_components(g), key = len)
    n = len(g.nodes())
    x.append(0)
    y.append(len(largest_component) * 1. / n)
    R = 0.0
    for i in range(1, n):
        g.remove_node(l.pop(0)[0])
        largest_component = max(networkx.connected_components(g), key = len)
        x.append(i * 1. / n)
        R += len(largest_component) * 1. / n
        y.append(len(largest_component) * 1. / n)
    return x, y, 0.5 - R / n

def find_communities_modularity(G, max_iter=None):
    '''
    INPUT:
        G: networkx Graph
        max_iter: (optional) if given, maximum number of iterations
    OUTPUT: list of lists of strings (node names)

    Run the Girvan-Newman algorithm on G and find the communities with the
    maximum modularity.
    '''
    degrees = G.degree()
    num_edges = G.number_of_edges()
    G1 = G.copy()
    best_modularity = -1.0
    best_comps = nx.connected_components(G1)
    i = 0
    while G1.number_of_edges() > 0:
        subgraphs = nx.connected_component_subgraphs(G1)
        modularity = get_modularity(subgraphs, degrees, num_edges)
        if modularity > best_modularity:
            best_modularity = modularity
            best_comps = list(nx.connected_components(G1))
        girvan_newman_step(G1)
        i += 1
        if max_iter and i >= max_iter:
            break
    return best_comps

def curation(strange,Z,CV,FV,EV):
	v_,ee_ = TRANS(strange)
	e_ = set(CAT(ee_))
	vv = set(CAT([EV[e] for e in e_]))
	VIEW(STRUCT(AA(POLYLINE)([[EV[e1],EV[e2]] for e1,e2 in ee_])))
	import networkx as nx
	G=nx.Graph()
	G.add_nodes_from(e_)
	G.add_edges_from(ee_)
	nx.connected_components(G)
	edges2remove = []; newEdges = []
	for edges in nx.connected_components(G):
		print edges
		edict = defaultdict(list)
		verts = CAT([EV[e] for e in edges])
		for v in verts: 
			edict[v] += [1]
		newEdge = [v for v in edict.keys() if len(edict[v])==1]
		edges2remove += list(edges)
		newEdges += [newEdge]
	v = 0; 
	W = [None for k in range(len(Z)-len(v_))]
	newindex = [None for k in range(len(Z))]
	for k,vert in enumerate(Z):
		if k not in v_:
			W[v] =  vert
			newindex[k]=v
			v += 1
	CW = [list(set([newindex[v] for v in cv if newindex[v]!=None])) for cv in CV]
	FW = [[newindex[v] for v in fv if newindex[v]!=None] for fv in FV if len(set(fv))>1]
	print "FV = ",FV
	EW = [[newindex[v] for v in ev if newindex[v]!=None] for e,ev in enumerate(EV) if e not in edges2remove]
	EW = [[w1,w2] for w1,w2 in EW if w1!=w2]  + [[newindex[v1],newindex[v2]] for [v1,v2] in newEdges]
	return W,CW,FW,EW

	def conn_comps(self):
		self.components={}
		for n in nx.connected_components(self.G):
			components[len(n)]=components.get(len(n),0)+1
		self.n_comp=len(nx.connected_components(self.G))
		self.comp_mode=mode(components)
		self.comp_average=average(components)

 def pruneBranches(self, G): #Based on high-degree neighbors - use on a pruned MST, before merging chains
     nodesToMerge = set()
     for n in G:
         if G.degree(n) >= 3:
             for nbr in G[n]:
                 if G.degree(nbr) >= 3:
                     nodesToMerge.add(n)
                     nodesToMerge.add(nbr)
                     
     SG = G.subgraph(nodesToMerge)
     print "pruneBranches2 is going to merge the following node clusters:", nx.connected_components(SG)
     #We're going to merge each connected component in this subgraph into a single node
     for C in nx.connected_components(SG):
         if len(C) <= 1: #should always be false
             raise InputError('Single node in pruneBranches2 subgraph')
         else:
             maxDegree = 0
             maxNode = None
             for n in C:
                 if SG.degree(n) > maxNode:
                     maxDegree = SG.degree(n)
                     maxNode = n
             C.remove(maxNode)
             for n in C:
                 self.mergeGraphNodes(G, maxNode, n)

    def detectGroupsNX(self, nodes, links):
        '''
        detect groups using the NetworkX library. It uses a simple algorithm of
        remove "highest between_ness centrality nodes" and then detecting the graph
        split.
        '''
        def make_nx_graph(nodes, links):
            G = nx.Graph()
            G.add_nodes_from(nodes)
            link_tupples = [(node_tupple[0], node_tupple[1], val['count'])
                            for node_tupple, val in links.iteritems()]
            G.add_weighted_edges_from(link_tupples)
            return G

        def calc_betweenness(graph):
            centrality = nx.edge_betweenness_centrality(graph, False)
            centrality = sorted(
                centrality.iteritems(), key=operator.itemgetter(1), reverse=True)
            return centrality

        def calc_average(iter):
            averge = sum(iter) * 1.0 / len(centrality)
            return averge

        def centrality_stddev(centrality):
            # find the standard deviation of centrality
            count = len(centrality)
            average = 0.0
            std_dev = 0.0
            if count > 0:
                total = sum(
                    itertools.imap(operator.itemgetter(1), centrality)) * 1.0
                average = total / count
                variance_sum = sum(map(
                    lambda x: (x - average)**2, itertools.imap(operator.itemgetter(1), centrality)))
                std_dev = math.sqrt(variance_sum / count)
            return average, std_dev

        graph = make_nx_graph(nodes, links)
        groups = nx.connected_components(graph)

        print "number of groups detected %d" % len(groups)
        centrality = calc_betweenness(graph)
        average, censtddev = centrality_stddev(centrality)
        # remove all the edges with centrality > (average+stddev)
        centrality_maxval = average + (censtddev * 1.96)
        edges = [edge_info[0]
                 for edge_info in centrality if edge_info[1] >= centrality_maxval]
        # Store the information about suggested edges to remove
        self.edgesToRemove = edges
        graph.remove_edges_from(edges)
        print "edges removed %d" % len(edges)

        # now extract the groups (or connected components) from the graph.
        groups = nx.connected_components(graph)
        groups = sorted(groups, key=lambda g: len(g), reverse=True)
        print "number of groups detected %d" % len(groups)
        return groups

def linker_graph(seeds, interactome):
    """Creates an expanded graph by iteratively connecting genes of
    interest into a largest connected component (LCC).

    An initial graph is created by finding the largest connected
    component of genes of interest only. Each unincluded GOI is
    included by finding any intermediate interactor to connect
    it to the LCC. That is, if there is path such that
    Unincluded GOI <-> Other Gene <-> Included GOI, the path is
    added to the network. This process is repeated until either
    all GOIs are included, or no more GOIs can be included.
    
    Paramters
    ---------
    seeds : An iterable of genes to generate the network from

    interactome : The consolidate interactome

    Returns
    -------
    An expanded NetworkX Graph
    """

    seeds = set(seeds) & set(interactome)

    nodes = set(max(
        nx.connected_components(
            interactome.subgraph(seeds)),
        key=len))

    orphaned = seeds - nodes

    while orphaned:
        # Find any common interactors for the orphans and the current graph
        orphan_edges = set()
        for orphan in orphaned:
            orphan_edges.update(interactome[orphan].keys())

        linkers = set()
        for node in nodes:
            linkers.update(set(interactome[node].keys()) & orphan_edges)


        # Give up if we can't find any new connections
        if not linkers:
            break

        # Update the graph to include the old graph, new linkers and any
        # newly included orphans
        nodes = set(max(
            nx.connected_components(
                interactome.subgraph(nodes | linkers | orphaned)),
            key=len))

        orphaned -= nodes

    return interactome.subgraph(nodes)