def _calc_counts(G, min_size, gap, scaf_only=False):
    sizes = list()
    for comp in nx.weakly_connected_components(G):

        # skip non scaffolds.
        if scaf_only == True:
            if len(comp) < 2:
                continue

        # add contig size.
        size = 0
        for n in comp:
            size += G.node[n]['width']

        # add gap size.
        if gap != False:
            for p,q in G.edges(comp):
                size += G[p][q]['gap']

        # skip this.
        if min_size != False and size < min_size:
            continue

        # save the size.
        sizes.append(len(comp))

    return sizes

	def comps(self):
		''' generator for connected components '''
		
		# loop over components.
		for comp in nx.weakly_connected_components(self):
			
			# create subgraph.
			subg = nx.DiGraph()
			
			# build node list.
			nlist = []
			for n in comp:
				nlist.append( (n, self.node[n]) )
				
			# add nodes.
			subg.add_nodes_from(nlist)
			
			# build edge list.
			elist = []
			for e in self.edges(comp):
				elist.append( (e[0], e[1], self[e[0]][e[1]]) )
				
			# add edges.
			subg.add_edges_from(elist)
			
			# yield the subgraph.
			yield subg

def lesion_met_largest_weak_component(G, orig_order=None):
    """
    Get largest weak component size of a graph.

    Parameters
    ----------
    G : directed networkx graph
        Graph to compute largest component for
    orig_order : int
        Define orig_order if you'd like the largest component proportion

    Returns
    -------
    largest weak component size : int
        Proportion of largest remaning component size if orig_order
        is defined. Otherwise, return number of nodes in largest component.
    """

    components = sorted(nx.weakly_connected_components(G), key=len,
                        reverse=True)
    if len(components) > 0:
        largest_component = len(components[0])
    else:
        largest_component = 0.

    # Check if original component size is defined
    if orig_order is not None:
        return largest_component / float(orig_order)
    else:
        return largest_component

def mark_vpn(graph, vpn_macs):
    components = map(frozenset, nx.weakly_connected_components(graph))
    components = filter(vpn_macs.intersection, components)
    nodes = reduce(lambda a, b: a | b, components, set())

    for node in nodes:
        for k, v in graph[node].items():
            v['vpn'] = True

def report_stats(G, params):
    print 'Nodes: %d.  Edges: %d'%(G.number_of_nodes(), G.number_of_edges())
    sccs = nx.strongly_connected_components(G)
    wccs = nx.weakly_connected_components(G)
    print 'Strongly ccs: %d, Weakly ccs: %d'%(len(sccs), len(wccs))

    sizes_sccs, sizes_wccs = ([len(c) for c in sccs], [len(c) for c in wccs])
    print 'Singletons. Strongly: %d, Weakly: %d'%(sum(np.array(sizes_sccs)==1), sum(np.array(sizes_wccs)==1))
    print [len(c) for c in sccs[:10]]

    def singletons(self):
        """A singleton is a weakly connected component that has only one node.

        Returns:
            A list with the singleton nodes.
        """
        components = networkx.weakly_connected_components(self.nxgraph)
        return [component[0] for component in components
                if len(component) == 1]

def test_weakly_connected_components(testgraph):
    """
    Test strongly connected components
    """

    comps0 = nx.weakly_connected_components(testgraph[0])
    comps1 = sg.components.weak(testgraph[1])

    assert_components_equal(comps0, comps1)

def find_len_2_control_kernals(deterministic_transition_graph, nodes_list, attractor_ID):
    """uses the deterministic_transition_graph and specified attractor to find 
    all pairs of control nodes if they exist

    note: these aren't "strict" control kernels because they specify the states needed to be in the
    main attractor. controlling them doesn't necessarily change what attractor you'll be in.
    """


    subgraphs = [g for g in nx.weakly_connected_components(deterministic_transition_graph)]
    # index_of_largest_subgraph = max(enumerate(all_subgraph_sets), key = lambda tup: len(tup[1]))[0]

    all_subgraph_sets = []
    all_but_attractor_subgraph = []
    for sg in subgraphs:
        if attractor_ID not in sg:
            all_state_sets = []
            for nID in sg:
                all_state_sets.append(set(time_evol.decimal_to_binary(nodes_list,nID).items()))
            all_subgraph_sets.append(all_state_sets)
            all_but_attractor_subgraph.append(all_state_sets)
        else:
            all_state_sets = []
            for nID in sg:
                all_state_sets.append(set(time_evol.decimal_to_binary(nodes_list,nID).items()))
            all_subgraph_sets.append(all_state_sets)

    state_0_list = [0 for i in range(len(nodes_list))]
    state_1_list = [1 for i in range(len(nodes_list))]

    possible_states = zip(nodes_list,state_0_list) + zip(nodes_list,state_1_list) # list of (node,state)

    possible_pairs = [combo for combo in combinations(possible_states, 2)]

    possible_pairs_pared = copy.deepcopy(possible_pairs)

    # remove pairs where both keys are the same (eg. gF and gF)
    for pair in possible_pairs:

        if pair[0][0] == pair[1][0] and (pair in possible_pairs_pared):

            possible_pairs_pared.remove(pair)

    # remove pairs when any of the networks in the non-attractor subgraph contain that pair
    # (ie. that pair can not possibly be a control kernel because it is present in the wrong attractor)
    for sg in all_but_attractor_subgraph:

        for state_set in sg:

            for pair in possible_pairs:

                if (pair[0] in state_set) and (pair[1] in state_set) and (pair in possible_pairs_pared):

                    possible_pairs_pared.remove(pair)

    return possible_pairs_pared # return a list ((node,state),(node,state)) pairs that are control kernels

def main():
    if len(sys.argv) != 2:
        print "Error!\nCorrect usage is:\n\t"
        print "python visualize_osm_test_region.py [osm_test_region_for_draw.dat]"
        return
   
    G = nx.read_gpickle(sys.argv[1])

    components = nx.weakly_connected_components(G)
    print "There are %d connected components."%len(components)

    H = G.subgraph(components[0])
    G = H

    osm_for_drawing = OSM_DRAW(G)
    easting, northing = osm_for_drawing.node_list()
    edge_list = osm_for_drawing.edge_list()

    #map_decomposition = MapDecomposition()
    #map_decomposition.primitive_decomposition(G, 10.0)

    #return

    fig = plt.figure(figsize=const.figsize)
    ax = fig.add_subplot(111, aspect='equal')
    #print edge_list
    arrow_params = {'length_includes_head':True, 'shape':'full', 'head_starts_at_zero':False}
    for segment in edge_list:
        u = segment[1][0] - segment[0][0]
        v = segment[1][1] - segment[0][1]
        ax.arrow(segment[0][0], segment[0][1], u, v, width=0.5, head_width=5,\
                    head_length=10, overhang=0.5, **arrow_params)
    #edge_collection = LineCollection(edge_list, colors='gray', linewidths=2)
    #ax.add_collection(edge_collection)

    # Junction nodes
    for node in G.nodes():
        if G.degree(node) > 2:
            ax.plot(G.node[node]['data'].easting,
                    G.node[node]['data'].northing,
                    'ro')

    # Connected components
    #for index in range(0, len(components)):
    #    color = const.colors[index%7]
    #    print len(components[index])
    #    for node in components[index]:
    #        ax.plot(G.node[node]['data'].easting,
    #                G.node[node]['data'].northing,
    #                'o', color=color)
    #    break

    ax.set_xlim([const.RANGE_SW[0], const.RANGE_NE[0]])
    ax.set_ylim([const.RANGE_SW[1], const.RANGE_NE[1]])

    plt.show()

def main():
    parser = OptionParser()
    parser.add_option("-m", "--osm", dest="osm_data", help="Input open street map data (typically in gpickle format)", metavar="OSM_DATA", type="string")
    parser.add_option("-t", "--track_data", dest="track_data", help="Input GPS tracks", metavar="TRACK_DATA", type="string")
    parser.add_option("-o", "--output_osm", dest="output_osm", help="Output file name (suggested extention: gpickle)", metavar="OUTPUT", type="string")
    parser.add_option("--test_case", dest="test_case", type="int", help="Test cases: 0: region-0; 1: region-1; 2: SF-region.", default=0)
    (options, args) = parser.parse_args()
    
    if not options.osm_data:
        parser.error("Input osm_data not found!")
    if not options.track_data:
        parser.error("Input track_data not found!")
    if not options.output_osm:
        parser.error("Output image not specified!")

    R = const.R 
    if options.test_case == 0:
        LOC = const.Region_0_LOC
    elif options.test_case == 1:
        LOC = const.Region_1_LOC
    elif options.test_case == 2:
        LOC = const.SF_LOC
    else:
        parser.error("Test case indexed %d not supported!"%options.test_case)

    G = nx.read_gpickle(options.osm_data)
    components = nx.weakly_connected_components(G)
    H = G.subgraph(components[0])
    G = H

    tracks = 

    osm_for_drawing = OSM_DRAW(G)
    edge_lists = osm_for_drawing.edge_list()
    line_strings = []
    for edge_list in edge_lists:
        line = LineString(edge_list)
        line_strings.append(line)


    fig = plt.figure(figsize=(10, 10))
    ax = plt.Axes(fig, [0., 0., 1., 1.], aspect='equal')
    ax.set_axis_off()
    fig.add_axes(ax)

    ROAD_WIDTH = 7 # in meters
    for line_string in line_strings:
        polygon = line_string.buffer(ROAD_WIDTH)
        patch = PolygonPatch(polygon, facecolor='k', edgecolor='k')
        ax.add_patch(patch)
    ax.set_xlim([LOC[0]-R, LOC[0]+R])
    ax.set_ylim([LOC[1]-R, LOC[1]+R])
    fig.savefig(options.output_img, dpi=100)
    plt.close()
    
    return

def trn_stats(genes, trn, t_factors, version):
    LOGGER.info("Computing TRN statistics")
    nodes = sorted(trn.nodes_iter())
    node2id = {n: i for (i, n) in enumerate(nodes)}
    id2node = {i: n for (i, n) in enumerate(nodes)}
    (grn, node2id) = to_simple(trn.to_grn(), return_map=True)
    nodes = sorted(grn.nodes_iter())
    regulating = {node for (node, deg) in grn.out_degree_iter() if deg > 0}
    regulated = set(nodes) - regulating
    components = sorted(nx.weakly_connected_components(grn), key=len,
            reverse=True)
    data = dict()
    for (a, b) in itertools.product(("in", "out"), repeat=2):
        data["{a}_{b}_ass".format(a=a, b=b)] = nx.degree_assortativity_coefficient(grn, x=a, y=b)
    census = triadic_census(grn)
    forward = census["030T"]
    feedback = census["030C"]
    num_cycles = sum(1 for cyc in nx.simple_cycles(grn) if len(cyc) > 2)
    in_deg = [grn.in_degree(node) for node in regulated]
    out_deg = [grn.out_degree(node) for node in regulating]
    data["version"] = version,
    data["release"] = pd.to_datetime(RELEASE[version]),
    data["num_genes"] = len(genes),
    data["num_tf"] = len(t_factors),
    data["num_nodes"] = len(nodes),
    data["num_regulating"] = len(regulating),
    data["num_regulated"] = len(regulated),
    data["num_links"] = grn.size(),
    data["density"] = nx.density(grn),
    data["num_components"] = len(components),
    data["largest_component"] = len(components[0]),
    data["feed_forward"] = forward,
    data["feedback"] = feedback,
    data["fis_out"] = trn.out_degree(TranscriptionFactor[FIS_ID, version]),
    data["hns_out"] = trn.out_degree(TranscriptionFactor[HNS_ID, version]),
    data["cycles"] = num_cycles,
    data["regulated_in_deg"] = mean(in_deg),
    data["regulating_out_deg"] = mean(out_deg),
    data["hub_out_deg"] = max(out_deg)
    stats = pd.DataFrame(data, index=[1])
    in_deg = [grn.in_degree(node) for node in nodes]
    out_deg = [grn.out_degree(node) for node in nodes]
    bc = nx.betweenness_centrality(grn)
    bc = [bc[node] for node in nodes]
    dists = pd.DataFrame({
            "version": version,
            "release": [pd.to_datetime(RELEASE[version])] * len(nodes),
            "node": [id2node[node].unique_id for node in nodes],
            "regulated_in_degree": in_deg,
            "regulating_out_degree": out_deg,
            "betweenness": bc
        })
    return (stats, dists)

    def rooted_core_interface_pairs(self, root, thickness=None,  for_base=False,
                                        hash_bitmask=None,
                                      radius_list=[],
                                      thickness_list=None,
                                      node_filter=lambda x, y: True):
        """

        Parameters
        ----------
        root:
        thickness:
        args:

        Returns
        -------

        """

        ciplist = super(self.__class__, self).rooted_core_interface_pairs(root, thickness, for_base=for_base,
                                        hash_bitmask=hash_bitmask,
                                      radius_list=radius_list,
                                      thickness_list=thickness_list,
                                      node_filter=node_filter)



        # numbering shards if cip graphs not connected
        for cip in ciplist:
            if not nx.is_weakly_connected(cip.graph):
                comps = [list(node_list) for node_list in nx.weakly_connected_components(cip.graph)]
                comps.sort()

                for i, nodes in enumerate(comps):

                    for node in nodes:
                        cip.graph.node[node]['shard'] = i

        '''
        solve problem of single-ede-nodes in the core
        this may replace the need for fix_structure thing
        this is a little hard.. may fix later

        it isnt hard if i write this code in merge_core in ubergraphlearn

        for cip in ciplist:
            for n,d in cip.graph.nodes(data=True):
                if 'edge' in d and 'interface' not in d:
                    if 'interface' in cip.graph.node[ cip.graph.successors(n)[0]]:
                        #problem found
        '''

        return ciplist

 def save(self, filename):
     # output the biggest weakly connected components
     vertices = list(sorted(nx.weakly_connected_components(self.graph),
                            key=len, reverse=True))[0]
     output_graph = nx.DiGraph()
     output_graph.add_nodes_from([(node, data)
                                  for node, data
                                  in self.graph.nodes(data=True)
                                  if node in vertices])
     output_graph.add_edges_from([(src, dest, data)
                                  for src, dest, data
                                  in self.graph.edges(data=True)
                                  if src in vertices and dest in vertices])
     nx.drawing.nx_pydot.write_dot(output_graph, filename)
     return self

 def _connected_component(self):
     if self.graph.is_directed():
         edges_in_wcc = set()
         nodes_in_wcc = set()
         wcc_counter = 0
         for wcc in nx.weakly_connected_components(self.graph):
             wcc_counter += 1
             print(type(wcc))
             nodes_in_wcc.add(wcc.nodes())
             edges_in_wcc.add(wcc.edges())
     else:  # undirected
         scc_counter = 0
         for scc in nx.strongly_connected_components(self.graph):
             scc_counter += 0
             print(type(scc))

    def test_zero_d_to_molecule_graph(self):
        comp_graphs = [self.mol_structure.graph.subgraph(c) for c in
                       nx.weakly_connected_components(self.mol_structure.graph)]

        mol_graph = zero_d_graph_to_molecule_graph(self.mol_structure,
                                                   comp_graphs[0])

        self.assertEqual(mol_graph.get_connected_sites(0)[0].index, 1)
        self.assertEqual(mol_graph.get_connected_sites(1)[1].index, 2)
        self.assertEqual(mol_graph.molecule.num_sites, 3)

        # test catching non zero dimensionality graphs
        comp_graphs = [self.graphite.graph.subgraph(c) for c in
                       nx.weakly_connected_components(self.graphite.graph)]
        self.assertRaises(ValueError, zero_d_graph_to_molecule_graph,
                          self.graphite, comp_graphs[0])

        # test for a troublesome structure
        s = loadfn(os.path.join(test_dir, "PH7CN3O3F.json.gz"))
        bs = CrystalNN().get_bonded_structure(s)
        comp_graphs = [bs.graph.subgraph(c) for c in
                       nx.weakly_connected_components(bs.graph)]
        mol_graph = zero_d_graph_to_molecule_graph(bs, comp_graphs[0])
        self.assertEqual(mol_graph.molecule.num_sites, 12)

def get_largest_component(G, strongly=False):
    """
    Return a subgraph of the largest weakly or strongly connected component
    from a directed graph.

    Parameters
    ----------
    G : networkx multidigraph
    strongly : bool
        if True, return the largest strongly instead of weakly connected
        component

    Returns
    -------
    G : networkx multidigraph
        the largest connected component subgraph from the original graph
    """

    start_time = time.time()
    original_len = len(list(G.nodes()))

    if strongly:
        # if the graph is not connected retain only the largest strongly connected component
        if not nx.is_strongly_connected(G):
            
            # get all the strongly connected components in graph then identify the largest
            sccs = nx.strongly_connected_components(G)
            largest_scc = max(sccs, key=len)
            G = induce_subgraph(G, largest_scc)
            
            msg = ('Graph was not connected, retained only the largest strongly '
                   'connected component ({:,} of {:,} total nodes) in {:.2f} seconds')
            log(msg.format(len(list(G.nodes())), original_len, time.time()-start_time))
    else:
        # if the graph is not connected retain only the largest weakly connected component
        if not nx.is_weakly_connected(G):
            
            # get all the weakly connected components in graph then identify the largest
            wccs = nx.weakly_connected_components(G)
            largest_wcc = max(wccs, key=len)
            G = induce_subgraph(G, largest_wcc)
            
            msg = ('Graph was not connected, retained only the largest weakly '
                   'connected component ({:,} of {:,} total nodes) in {:.2f} seconds')
            log(msg.format(len(list(G.nodes())), original_len, time.time()-start_time))

    return G

def create_node_conncomp_graph(G,layer1,layer2,layer3):
    # print layer1
    npartition = list(nx.weakly_connected_components(G))
    print len(npartition)
    # G=nx.Graph(G)
    layers={'layer1':layer1,'layer2':layer2,'layer3':layer3}
    broken_partition={}
    for i,v in enumerate(npartition):
        vs=set(v)
        for ii,vv in layers.items():
            papa=vs.intersection(set(vv))
            if len(papa)==len(v):
                broken_partition['a_%i_%s_s' %(i,ii)]=v
            elif len(papa)>0:
                broken_partition['b_%i_%s' %(i,ii)]=list(papa)
                vs=vs-set(vv)
    # print rbroken_partition
    broken_graph=nx.Graph()
    rbroken_partition=dict()
    colors=[name for name,hex in matplotlib.colors.cnames.iteritems()]
    colors=list(set(colors)-set(['red','blue','green']))
    cl=dict()
    for i,v in broken_partition.items():
        name=i.split('_')
        for ii in v:
            # print ii
            rbroken_partition[ii]=i
        if name[-1]=='s':
            cl[name[1]]=colors.pop()
        elif name[0]=='b' and not cl.has_key(name[1]):
            cl[name[1]]=colors.pop()
    # print rbroken_partition
    for i,v in rbroken_partition.items():
        name=v.split('_')
        broken_graph.add_node(v,color=cl[name[1]])
        edg=G[i]
        # print edg
        for j in edg:
            if j not in broken_partition[v]:
                # print j
                if not broken_graph.has_edge(v,rbroken_partition[j]):
                    broken_graph.add_edge(v,rbroken_partition[j])
    
    return broken_graph,broken_partition,npartition

def resolve_connectivity_issue(test_graph):
    weakly_connected_components = networkx.weakly_connected_components(test_graph)
    wcc_list = list()

    # get components sorted  based on the length
    for component in weakly_connected_components:
        index = 0
        for wcc in wcc_list:
            if len(component) < len(wcc):
                break
            index += 1
        wcc_list.insert(index, component)

    print "sorted list"
    for wcc in wcc_list:
        print 'component length' + str(len(wcc))

    for i in range(len(wcc_list)-1):
        for j in range(i+1, len(wcc_list)):
            connect_weakly_connected_component(wcc_list[i], wcc_list[j], test_graph)

    def rooted_core_interface_pairs(self, root, thickness=None, **args):
        '''

        Args:
            root: int
            thickness:  
            **args:

        Returns:

        '''

        ciplist = super(self.__class__, self).rooted_core_interface_pairs(root, thickness, **args)


        #numbering shards if cip graphs not connected
        for cip in ciplist:
            if not nx.is_weakly_connected(cip.graph):
                comps = [list(node_list) for node_list in nx.weakly_connected_components(cip.graph)]
                comps.sort()

                for i, nodes in enumerate(comps):

                    for node in nodes:
                        cip.graph.node[node]['shard'] = i

        '''
        solve problem of single-ede-nodes in the core
        this may replace the need for fix_structure thing
        this is a little hard.. may fix later

        it isnt hard if i write this code in merge_core in ubergraphlearn

        for cip in ciplist:
            for n,d in cip.graph.nodes(data=True):
                if 'edge' in d and 'interface' not in d:
                    if 'interface' in cip.graph.node[ cip.graph.successors(n)[0]]:
                        #problem found
        '''

        return ciplist

def main():
    print "time_evol module is the main code."
    edge_file = '../data/inputs/edges-init.dat' 
    # edge_file = '../data/inputs/edges_1-init.dat' 
    node_file = '../data/inputs/ant-nodes-init.dat'
    
    G = inet.read_network_from_file(edge_file, node_file)
    nodes_list = G.nodes() # print graph nodes without states. Could also do nx.nodes(G)

    # print G.nodes()
    initial_state = nx.get_node_attributes(G,'state')
    # print initial_state #OrderedDict(sorted(initial_state.items(), key=lambda t: t[0]))

    #####

    # timeseriesdata = time_series_all(G, nodes_list, 10)#, timesteps=20)

    # network_ID = 1 #this is the ID of the init state
    # biStates = decimal_to_binary(nodes_list, network_ID)
    # print 'initial state', biStates
    # for node in G.nodes():
    #     print node, timeseriesdata[node][1]

    G_transition_graph = net_state_transition(G, nodes_list) 
    
    ####################################################
    #### calculate the max number of steps to get to any given attractor in my system
    # target_from_source_dict = nx.single_source_shortest_path(G_transition_graph.to_undirected(), 43)
    # print max([len(target_from_source_dict[i]) for i in target_from_source_dict])
    # exit()
    ####################################################
    # nx.draw(G_transition_graph)
    # plt.show()

    attractors = find_attractor(G_transition_graph) # why does this change if I update G?

    print attractors.keys()
    # print nx.number_weakly_connected_components(G_transition_graph)
    # print [i for i in nx.weakly_connected_components(G_transition_graph)]
    print [len(i) for i in nx.weakly_connected_components(G_transition_graph)]