def main():
    tempo_dir = "../corpus-local/tempo-txt"
    file_regex = ".*\.txt"

    G = build_graph(tempo_dir, file_regex)
    """
  ccs = nx.clustering(G)
  avg_clust = sum(ccs.values()) / len(ccs)
  """
    print tempo_dir
    print "\tAda " + str(len(G.nodes())) + " node."
    print "\tAda " + str(len(G.edges())) + " edge."
    print "\tClustering coefficient      : " + str(nx.average_clustering(G))
    print "\tAverage shortest path length"
    for g in nx.connected_component_subgraphs(G):
        print "\t\t" + str(nx.average_shortest_path_length(g))

    kompas_dir = "../corpus-local/kompas-txt"
    G = build_graph(kompas_dir, file_regex)
    print kompas_dir
    print "\tAda " + str(len(G.nodes())) + " node."
    print "\tAda " + str(len(G.edges())) + " edge."
    print "\tClustering coefficient      : " + str(nx.average_clustering(G))
    print "\tAverage shortest path length"
    for g in nx.connected_component_subgraphs(G):
        print "\t\t" + str(nx.average_shortest_path_length(g))

def zc(G,list_G1,list_G2,f):#计算z值
    """
    输入参数：原始网络G，不保持连通性置乱网络list_G1,保持连通性置乱网络list_G2，要求网络连通的指标函数名f
    返回：该指标的不保持连通性z值z1,保持连通性z值z2
    """
    list_G_l0 = []; list_G_l1 = []; list_G_l2 = []
    for g in nx.connected_component_subgraphs(G):
        list_G_l0.append(f(g))
    for G1 in list_G1:
        for g1 in nx.connected_component_subgraphs(G1):
            list_G_l1.append(f(g1))#指标值列表
    for G2 in list_G2:
        for g2 in nx.connected_component_subgraphs(G2):
            list_G_l2.append(f(g2))
    #print list_G_l0, list_G_l1, list_G_l2
    G_l0 = np.mean(list_G_l0)
    G_l1 = np.mean(list_G_l1) #求均值 
    G_l2 = np.mean(list_G_l2)
    var_z1 = np.var(list_G_l1) #求方差
    var_z2 = np.var(list_G_l2)
    if var_z1 == 0: #若方差为0，则z值取0
        z1 = 0
    else:
        z1 = (G_l0 - G_l1)/var_z1#z值
    if var_z2 == 0: 
        z2 = 0
    else:
        z2 = (G_l0 - G_l2)/var_z2#z值
    return z1, z2

    def _extract_ccomponents(self, graph, threshold=0, min_size=2):
        # remove all vertices that have a score less then threshold
        cc_list = []

        if self.less_then:
            less_component_graph = graph.copy()
            for v, d in less_component_graph.nodes_iter(data=True):
                if d.get(self.attribute, False):
                    if d[self.attribute] < threshold:
                        less_component_graph.remove_node(v)
            for cc in nx.connected_component_subgraphs(less_component_graph):
                if len(cc) >= min_size:
                    cc_list.append(cc)

        # remove all vertices that have a score more then threshold
        if self.more_than:
            more_component_graph = graph.copy()
            for v, d in more_component_graph.nodes_iter(data=True):
                if d.get(self.attribute, False):
                    if d[self.attribute] >= threshold:
                        more_component_graph.remove_node(v)

            for cc in nx.connected_component_subgraphs(more_component_graph):
                if len(cc) >= min_size:
                    cc_list.append(cc)
        return cc_list

    def _good_k_break(self, old_edges, new_edges):
        """
        Checks that the break does not change chromomsome structure significantly
        """
        MIN_OVLP_SCORE = 0.9
        MAX_K_BREAK = 4
        if len(old_edges) > MAX_K_BREAK:
            return False

        new_adj_graph = self.adj_graph.copy()
        for u, v in old_edges:
            new_adj_graph.remove_edge(u, v)
        for u, v in new_edges:
            new_adj_graph.add_edge(u, v)

        all_nodes = new_adj_graph.nodes()
        old_sets = list(map(lambda g: set(g.nodes()),
                            nx.connected_component_subgraphs(self.adj_graph)))
        new_sets = list(map(lambda g: set(g.nodes()),
                            nx.connected_component_subgraphs(new_adj_graph)))
        if len(old_sets) != len(new_sets):
            return False

        for old_set in old_sets:
            max_overlap = 0
            best_score = 0
            for new_set in new_sets:
                overlap = len(old_set & new_set)
                if overlap > max_overlap:
                    max_overlap = overlap
                    best_score = float(overlap) / len(old_set | new_set)
            if best_score < MIN_OVLP_SCORE:
                return False

        return True

    def _extract_ccomponents(self, graph, threshold=0, min_size=2, max_size=20):
        # remove all vertices that have a score less then threshold
        cc_list = []

        if self.less_then:
            less_component_graph = graph.copy()
            for v, d in less_component_graph.nodes_iter(data=True):
                if self.get_attr_from_noded(d):
                    if self.get_attr_from_noded(d) < threshold:
                        less_component_graph.remove_node(v)
            for cc in nx.connected_component_subgraphs(less_component_graph):
                if len(cc) >= min_size and len(cc) <= max_size:
                    cc_list.append(cc)
                if len(cc) > max_size and self.shrink_graphs:
                    cc_list += list(self.enforce_max_size(cc, min_size, max_size))

        # remove all vertices that have a score more then threshold
        if self.more_than:
            more_component_graph = graph.copy()
            for v, d in more_component_graph.nodes_iter(data=True):
                if self.get_attr_from_noded(d):
                    if self.get_attr_from_noded(d) >= threshold:
                        more_component_graph.remove_node(v)

            for cc in nx.connected_component_subgraphs(more_component_graph):
                if len(cc) >= min_size and len(cc) <= max_size:
                    cc_list.append(cc)

                if len(cc) > max_size and self.shrink_graphs:
                    cc_list += list(self.enforce_max_size(cc, min_size, max_size, choose_cut_node=max))

        return cc_list

def printStats(filename):
	'''
	Converts json adjacency list into networkx to calculate and print the
	graphs's 
	  - average clustering coefficient
	  - overall clustering coefficient
	  - maximum diameter
	  - average diameter
	  - number of paritions using community.best_parition
	  - modularity of community.best_partition
	'''
	g = makeGraphFromJSON(filename)
	
	print "Average Clustering Coefficient: %f" % nx.average_clustering(g)
	print "Overall Clustering Coefficient: %f" % nx.transitivity(g)
	
	connected_subgraphs = list(nx.connected_component_subgraphs(g))
	largest = max(nx.connected_component_subgraphs(g), key=len)
	print "# Connected Components: %d" % len(connected_subgraphs)
	print "    Maximal Diameter: %d" % nx.diameter(largest)
	print "    Average Diameter: %f" % nx.average_shortest_path_length(largest)

	# Find partition that maximizes modularity using Louvain's algorithm
	part = community.best_partition(g)	
	print "# Paritions: %d" % (max(part.values()) + 1)
	print "Louvain Modularity: %f" % community.modularity(part, g)

def get_boundary_for_label(data, classifier, num_label, step):
    # See
    # http://en.wikipedia.org/wiki/Postcodes_in_the_United_Kingdom#Operation_and_application
    # for the various divisions.
    t_start = time.time()
    district = data[data[:,0] == num_label, 1:]

    # Align grid to nearest "step". Also grow border by 25 units to
    # to make sure the marching squares can build a full loop.
    x0, y0 = np.floor(district.min(0) / step - 25) * step
    x1, y1 = np.ceil(district.max(0) / step + 25) * step

    # Use KNN to colour a grid that covers the district
    xx, yy = np.mgrid[x0:x1:step, y0:y1:step]
    prediction = classifier.predict(np.c_[xx.ravel(), yy.ravel()]).reshape(xx.shape)

    # Split predicted labels into inside/outside
    prediction = (prediction == num_label).astype('u1')

    # We transpose to make reasoning about the lookups easier.
    prediction = prediction.transpose()

    # zero-pad predictions to make sure marching squares creates
    # closed outlines.
    tmp = np.zeros((prediction.shape[0] + 2, prediction.shape[1] + 2), dtype='u1')
    tmp[1:-1,1:-1] = prediction
    prediction = tmp

    outline = networkx.Graph()

    h, w = prediction.shape
    
    for i, j in np.ndindex(h - 1, w - 1):
        # We use tostring() as a cheap, hashable lookup type for the
        # marching squared implementation.

        # Dimension 0 ~ y ~ i, dim 1 ~ x ~ j:
        piter = iter(MARCHING_SQUARE_LOOKUP[prediction[i:i+2,j:j+2].tostring()])

        for rel1, rel2 in zip(piter, piter):
            p1 = int(x0 + step * (j + rel1[0])), int(y0 + step * (i + rel1[1]))
            p2 = int(x0 + step * (j + rel2[0])), int(y0 + step * (i + rel2[1]))

            outline.add_node(p1)
            outline.add_node(p2)
            outline.add_edge(p1, p2)

    # Pick the largest subgraph, other graphs are most likely outliers.
    logging.info(
        "%s: Found %s connected graphs in %.2fs",
        num_label,
        len(networkx.connected_component_subgraphs(outline)),
        time.time() - t_start,
    )
    largest = max(
        networkx.connected_component_subgraphs(outline),
        key=lambda x: x.size()
    )
    return list(shapely.ops.polygonize(largest.edges()))[0]

def get_small_worldness(filename):
  import networkx as nx
  threshold = 0
  f = open(filename[:-4]+'_small_worldness.dat','w')
  for i in range(0,101):
    threshold = float(i)/100
    G = get_threshold_matrix(filename, threshold)
    ER_graph = nx.erdos_renyi_graph(nx.number_of_nodes(G), nx.density(G))

    cluster = nx.average_clustering(G)
    ER_cluster = nx.average_clustering(ER_graph)
    
    transi = nx.transitivity(G)
    ER_transi = nx.transitivity(ER_graph)

    print 'threshold: %f, average cluster coefficient: %f, random nw: %f, transitivity: %f, random nw: %f' %(threshold, cluster, ER_cluster, transi, ER_transi)

    f.write("%f\t%f\t%f" % (threshold, cluster, ER_cluster))
    components = nx.connected_component_subgraphs(G)
    ER_components = nx.connected_component_subgraphs(ER_graph)

    values = []
    ER_values = []
    for i in range(len(components)):
      if nx.number_of_nodes(components[i]) > 1:
        values.append(nx.average_shortest_path_length(components[i]))
    for i in range(len(ER_components)):
      if nx.number_of_nodes(ER_components[i]) > 1:
        ER_values.append(nx.average_shortest_path_length(ER_components[i]))
    if len(values) == 0:
      f.write("\t0.")
    else:
      f.write("\t%f" % (sum(values)/len(values)))

    if len(ER_values) == 0:
      f.write("\t0.")
    else:
      f.write("\t%f" % (sum(ER_values)/len(ER_values)))
    
    f.write("\t%f\t%f" % (transi, ER_transi))  
    
    if (ER_cluster*sum(values)*len(values)*sum(ER_values)*len(ER_values)) >0 :
      S_WS = (cluster/ER_cluster) / ((sum(values)/len(values)) / (sum(ER_values)/len(ER_values)))
    else:
      S_WS = 0.
    if (ER_transi*sum(values)*len(values)*sum(ER_values)*len(ER_values)) >0 :
      S_Delta = (transi/ER_transi) / ((sum(values)/len(values)) / (sum(ER_values)/len(ER_values)))
    else:
      S_Delta = 0.
    
    f.write("\t%f\t%f" % (S_WS, S_Delta))  
    f.write("\n")
    
  f.close()  
  print "1:threshold 2:cluster-coefficient 3:random-cluster-coefficient 4:shortest-pathlength 5:random-shortest-pathlength 6:transitivity 7:random-transitivity 8:S-Watts-Strogatz 9:S-transitivity" 

def hemst(G, k):
    nc = 1
    mst = G
    point_set = {}
    while nc != k:
        nc = 1
        mst = nx.minimum_spanning_tree(mst)
        weights = np.array([attrs['weight'] for _,_,attrs in mst.edges(data=True)])
        mean_w = weights.mean()
        std = weights.std()

        for a,b,attrs in mst.edges(data=True):
            w = attrs['weight']
            if w > mean_w + std:
                mst.remove_edge(a,b)
                nc+=1

        if nc < k:
            while nc != k:
                remove_longest_edge(mst)
                nc+=1
            break

        if nc > k:
            sG = nx.connected_component_subgraphs(mst)
            centroid_nodes = []
            for g in sG:
                cl = nx.closeness_centrality(g)
                sorted_set_nodes = sorted(cl.items(), key=lambda a: a[1])
                closest_to_c = sorted_set_nodes[0][0]

                point_set[closest_to_c] = g.nodes()
                for p, _ in sorted_set_nodes[1:]:
                    if p in point_set:
                        point_set[closest_to_c]+= point_set[p]

                centroid_nodes.append(closest_to_c)

            edges=itertools.combinations(centroid_nodes,2)
            mst.clear()
            mst.add_nodes_from(centroid_nodes)
            mst.add_edges_from(edges)
            for u,v in mst.edges():
                weight = G.get_edge_data(u,v)["weight"]
                nx.set_edge_attributes(mst, "weight", {(u,v):weight})
                
    sG = nx.connected_component_subgraphs(mst)
    if point_set:
        for g in sG:
            for node in g.nodes():
                if node in point_set:
                    g.add_nodes_from(point_set[node])

    return sG

def get_small_worldness(G, thr):
	f = open(out_prfx + 'small_worldness.dat', 'a')
	g = open(out_prfx + 'cc_trans_ER.dat', 'a')
	#g.write('r(thre.)\t\cc_A\tcc_ER\ttran_A\ttran_ER\n')
	ER_graph = nx.erdos_renyi_graph(nx.number_of_nodes(G), nx.density(G))
	# erdos-renyi, binomial random graph generator ...(N,D:density)	
	cluster = nx.average_clustering(G)   # clustering coef. of whole network
	ER_cluster = nx.average_clustering(ER_graph)	#cc of random graph
	
	transi = nx.transitivity(G)
	ER_transi = nx.transitivity(ER_graph)

	g.write("%f\t%f\t%f\t%f\t%f\n" % (thr, cluster,ER_cluster,transi,ER_transi ))
	
	f.write("%f\t%f\t%f" % (thr, cluster, ER_cluster))
	components = nx.connected_component_subgraphs(G)
	ER_components = nx.connected_component_subgraphs(ER_graph)

	values = []
	ER_values = []
	for i in range(len(components)):
		if nx.number_of_nodes(components[i]) > 1:
			values.append(nx.average_shortest_path_length(components[i]))
	for i in range(len(ER_components)):
		if nx.number_of_nodes(ER_components[i]) > 1:
			ER_values.append(nx.average_shortest_path_length(ER_components[i]))
	if len(values) == 0:
		f.write("\t0.")
	else:
		f.write("\t%f" % (sum(values)/len(values))) # pathlenght

	if len(ER_values) == 0:
		f.write("\t0.")
	else:
		f.write("\t%f" % (sum(ER_values)/len(ER_values)))

	f.write("\t%f\t%f" % (transi, ER_transi))  

	if (ER_cluster*sum(values)*len(values)*sum(ER_values)*len(ER_values)) >0 :
		S_WS = (cluster/ER_cluster) / ((sum(values)/len(values)) / (sum(ER_values)/len(ER_values)))  
	else:
		S_WS = 0.
	if (ER_transi*sum(values)*len(values)*sum(ER_values)*len(ER_values)) >0 :
		S_Delta = (transi/ER_transi) / ((sum(values)/len(values)) / (sum(ER_values)/len(ER_values)))
	else:
		S_Delta = 0.

	f.write("\t%f\t%f" % (S_WS, S_Delta)) # S_WS ~ small worldness 
	f.write("\n")

	f.close() 
	g.close()	 

def eigenvector_apl(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)
    l = sorted(m.items(), key=operator.itemgetter(1), reverse=True)
    x = []
    y = []

    average_path_length = 0.0
    number_of_components = 0
    n = len(g.nodes())

    for sg in networkx.connected_component_subgraphs(g):
        average_path_length += networkx.average_shortest_path_length(sg)
    number_of_components += 1

    average_path_length /= number_of_components
    initial_apl = average_path_length

    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)
        average_path_length = 0.0
        number_of_components = 0

        for sg in networkx.connected_component_subgraphs(g):
            if len(sg.nodes()) > 1:
                average_path_length += networkx.average_shortest_path_length(sg)
            number_of_components += 1

        average_path_length = average_path_length / number_of_components

        x.append(i * 1. / initial_apl)
        r += average_path_length * 1. / initial_apl
        y.append(average_path_length * 1. / initial_apl)
    return x, y, r / initial_apl

def process_network(G, namespace):
    print 'Nodes:', len(G)
    print 'Edges:', G.number_of_edges()
    if namespace.clustering_coefficient:
        print 'Clustering Coefficient:', nx.average_clustering(G)
    if namespace.components:
        components = nx.connected_component_subgraphs(G)
        print 'Number of Components:', len(components)
        isles = [c for c in components if len(c) == 1]
        print 'Isles:', len(isles)
        print 'Largest Component Size:', len(components[0])
    else: components = None
    if namespace.cpl:
        if namespace.approximate_cpl:
            average_shortest_path_length = approximate_cpl
        else:
            print 'Using full slow CPL'
            average_shortest_path_length = nx.average_shortest_path_length
        if components is None:
            components = nx.connected_component_subgraphs(G)
        for i, g in enumerate(g for g in components if
                float(len(g))/float(len(G)) >
                namespace.component_size):
            print 'CPL %d: (%f)' % (i, float(len(g))/float(len(G)))
            print average_shortest_path_length(g)
    if namespace.assortativity:
        print 'Assortativity: NOT IMPLEMENTED.'
    if namespace.degree_distribution:
        hst = nx.degree_histogram(G)

        plt.subplot(121)
        plt.xscale('log')
        plt.yscale('log')
        plt.title("Degree Distribution")
        plt.ylabel("Occurrencies")
        plt.xlabel("Degree")
        plt.plot(range(len(hst)), hst, marker='+')

        plt.subplot(122)
        ccdf = pynetsym.mathutil.ccdf(hst)
        plt.xscale('log')
        plt.yscale('log')
        plt.title("CCDF Degree Distribution")
        plt.ylabel("$P(X>x)$")
        plt.xlabel("Degree")
        plt.plot(range(len(ccdf)), ccdf, color='red')

        if namespace.degree_distribution_out is None:
            plt.show()
        else:
            plt.save_fig(namespace.degree_distribution_out)

    def one_girvan_newman(self,G):
        def find_best_edge(G0):
            eb = nx.edge_betweenness_centrality(G0)
            eb_il = eb.items()
            eb_il.sort(key=lambda x: x[1], reverse=True)
            return eb_il[0][0]

        num_clusters = len(sorted(nx.connected_component_subgraphs(G),key=len,reverse=True))
        caused_split = False
        while not caused_split:
            G.remove_edge(*find_best_edge(G))
            components = sorted(nx.connected_component_subgraphs(G),key=len,reverse=True)
            if len(components) == num_clusters+1:
                caused_split = True

def rig(x, y, G, labs=None, res=1e-9):
    """ Compute the RIG metric on all components.

    Parameters
    ----------
    x : pd.Series or array_like
       Vector of nodes and their abundance in sample x
    y : pd.Series or array_like
       Vector of nodes and their abundance in sample y
    G : nx.Graph
       A connected graph of weighted edges

    Returns
    -------
    float :
       Distance between sample x and sample y

    Note
    ----
    If x or y is None, then 1 will be added to the total distance.
    If they are both None, then the distance will be zero.

    """
    if labs is not None:
        x = pd.Series(x, index=labs)
        y = pd.Series(y, index=labs)

    cost = 0
    _G = copy.deepcopy(G)
    # This converts all of the weights to integers
    for u, v, d in _G.edges(data=True):
        d["weight"] = int(d["weight"] / res)

    # This calculates the largest edge set to offset the insertion cost.
    weights = []
    for comp in nx.connected_component_subgraphs(_G):
        edges = list(comp.edges(data="weight"))
        if len(edges) > 0:
            weights.append(sum(list(zip(*edges))[2]))
    maxW = max(weights) + 1

    for comp in nx.connected_component_subgraphs(_G):
        nodes = set(comp.nodes())
        subx = x[nodes & set(x.keys())]
        suby = y[nodes & set(y.keys())]

    c = rig_component(comp, subx, suby, maxW)
    cost += c
    return (cost) * res

 def _plot_graphs(self):
     self.f,self.ax = plt.subplots(len(self.transition['all']),4,figsize=(14,10)) # first col motion , second distance
     self.f.suptitle('Scene : '+str(self.scene), fontsize=20)
     for feature in [0,2]:
         # plot the different graphs of motion and distance
         for sub,T in enumerate(self.transition['all']):
             plt.sca(self.ax[sub,feature])
             print 'plotting graph : '+str(sub+1)+' from '+str(len(self.transition['all']))
             if feature == 0: 
                 if T not in self.transition['motion']:
                     for i in self.transition['motion']:
                         if i<T: t=i
                 else: t=T
                 G=self.G_motion[t]['graph']
             elif feature == 2: 
                 if T not in self.transition['touch']:
                     for i in self.transition['touch']:
                         if i<T: t=i
                 else: t=T
                 G=self.G_touch[t]['graph']
             # layout graphs with positions using graphviz neato
             pos=nx.graphviz_layout(G,prog="neato")
             # color nodes the same in each connected subgraph
             C=nx.connected_component_subgraphs(G)
             cK = 0
             for i in C:  cK += 1
             C=nx.connected_component_subgraphs(G)
             colors = np.linspace(.2,.6,cK)
             for count,g in enumerate(C):
                 c=[colors[count]]*nx.number_of_nodes(g) # same color...
                 nx.draw(g,pos,node_size=80,node_color=c,vmin=0.0,vmax=1.0,with_labels=False)
                 #nx.draw_networkx_edges(g,pos, with_labels=False, edge_color=c[0], width=6.0, alpha=0.5)
             nx.draw_networkx_nodes(self.G,pos, node_color='b', node_size=100, alpha=1)
             nx.draw_networkx_nodes(self.G,pos, nodelist=['G'], node_color='r', node_size=100, alpha=1)
             nx.draw_networkx_nodes(self.G,pos, nodelist=[str(self.m_obj)], node_color='c', node_size=100, alpha=1)
             nx.draw_networkx_edges(G,pos, alpha=0.8)
             #nx.draw(G)  # networkx draw()
             self.ax[sub,feature].axis('on')
             self.ax[sub,feature].axis('equal')
             plt.tick_params(axis='x',which='both',bottom='off',top='off',labelbottom='off')
             plt.tick_params(axis='y',which='both',right='off',left='off',labelleft='off')
             if feature == 0:
                 self.ax[sub,feature].set_ylabel('frame : '+str(T))
                 if sub == 0:
                     self.ax[sub,feature].set_title('motion')
             if feature == 2:
                 self.ax[sub,feature].set_ylabel('frame : '+str(T))
                 if sub == 0:
                     self.ax[sub,feature].set_title('connectivity')

    def analyze_graph(self, graph):
        start_time = time.time()
        self.clear_stats()
        self._graph = graph

        self.node_count = nx.number_of_nodes(graph)
        self.edge_count = nx.number_of_edges(graph)

        degree_list = nx.degree(graph).values()

        self.connected_component_count = \
            sum(1 for cx in nx.connected_components(graph))
        if self.connected_component_count is 0:
            return

        self._connected_component_graphs = \
            nx.connected_component_subgraphs(graph)
        self._largest_component_graph = \
            max(nx.connected_component_subgraphs(graph), key=len)

        self.average_degree = sum(degree_list) / float(len(degree_list))
        self._degree_histogram = nx.degree_histogram(graph)
        spc = self.shortest_paths(graph)
        self.shortest_path_count = len(spc)
        self.maximum_shortest_path_length = \
            self.max_shortest_path_length(graph)

        if self.connected_component_count is 1:

            self.diameter = nx.diameter(graph)

            if self.node_count > 1:
                self.average_shortest_path_length = \
                    nx.average_shortest_path_length(graph)
                self.minimum_connectivity = self.min_connectivity(graph)

        if self.node_count > 0:
            self.maximum_degree = max(degree_list)
            self.minimum_degree = min(degree_list)

        if self.node_count > 1:
            dg = nx.degree_centrality(graph)
            self.maximum_degree_centrality = max(list(dg.values()))

        bc = nx.betweenness_centrality(graph)
        self.maximum_between_centrality = max(list(bc.values()))

        self.elapsed_time = time.time() - start_time

def sensi_diameter(G):
    import networkx as nx
    
    """
    Compute graph sensitivity to node removal, in terms of
    the difference in graph diameter on the removal of each
    node in turn.
     
    This uses local function x_diameter(G), which is modified
    from networkx.diamter(G) to work on XGraphs.
    
    DL Urban (9 Feb 2007)
    """
    
    # Starting diameter for full graph:
    
    if nx.is_connected(G):
        d0 = x_diameter(G)
    else:
        G0 = nx.connected_component_subgraphs(G) [0] # the largest subgraph
        d0 = x_diameter(G0)
        nc = nx.number_connected_components(G)	     # how many are there?
    
    sensi = {}
    
    for node in G.nodes():
        ex = G.edges(node) 		# a set of edges adjacent to node; 
        G.delete_edges_from(ex)		# remove all of these,
        G.delete_node(node)		# and then kill the node, too
        if nx.is_connected(G):
            dx = x_diameter(G)
            cuts = 0
        else:
            Gx = nx.connected_component_subgraphs(G) [0]	# the biggest
            ncx = nx.number_connected_components(G)
            if nc == ncx:
                cuts = 0
            else:
                cuts = 1
            dx = x_diameter(Gx)
        delta = d0 - dx
        G.add_node(node)		# put the node and edges back again
        G.add_edges_from(ex)
        sensi[node] = (cuts, delta)
 

    # create and return a tuple (cuts, delta)
    return sensi

def atlas6():
    """ Return the atlas of all connected graphs of 6 nodes or less.
        Attempt to check for isomorphisms and remove.
    """

    Atlas = graph_atlas_g()[0:208]  # 208
    # remove isolated nodes, only connected graphs are left
    U = nx.Graph()  # graph for union of all graphs in atlas
    for G in Atlas:
        zerodegree = [n for n in G if G.degree(n) == 0]
        for n in zerodegree:
            G.remove_node(n)
        U = nx.disjoint_union(U, G)

    # list of graphs of all connected components
    C = nx.connected_component_subgraphs(U)

    UU = nx.Graph()
    # do quick isomorphic-like check, not a true isomorphism checker
    nlist = []  # list of nonisomorphic graphs
    for G in C:
        # check against all nonisomorphic graphs so far
        if not iso(G, nlist):
            nlist.append(G)
            UU = nx.disjoint_union(UU, G)  # union the nonisomorphic graphs
    return UU

def lanl_graph():
    """ Return the lanl internet view graph from lanl.edges
    """
    import networkx as nx
    try:
        fh=open('lanl_routes.edgelist','r')
    except IOError:
        print "lanl.edges not found"
        raise

    G=nx.Graph()

    time={}
    time[0]=0 # assign 0 to center node
    for line in fh.readlines():
        (head,tail,rtt)=line.split()
        G.add_edge(int(head),int(tail))
        time[int(head)]=float(rtt)

    # get largest component and assign ping times to G0time dictionary
    G0=nx.connected_component_subgraphs(G)[0]
    G0.rtt={}
    for n in G0:
        G0.rtt[n]=time[n]

    return G0

def unitigs(args):
    """
    %prog unitigs best.edges

    Reads Celera Assembler's "best.edges" and extract all unitigs.
    """
    p = OptionParser(unitigs.__doc__)
    p.add_option("--maxerr", default=2, type="int", help="Maximum error rate")
    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(not p.print_help())

    bestedges, = args
    G = read_graph(bestedges, maxerr=opts.maxerr, directed=True)
    H = nx.Graph()
    intconv = lambda x: int(x.split("-")[0])
    for k, v in G.iteritems():
        if k == G.get(v, None):
            H.add_edge(intconv(k), intconv(v))

    nunitigs = nreads = 0
    for h in nx.connected_component_subgraphs(H, copy=False):
        st = [x for x in h if h.degree(x) == 1]
        if len(st) != 2:
            continue
        src, target = st
        path = list(nx.all_simple_paths(h, src, target))
        assert len(path) == 1
        path, = path
        print "|".join(str(x) for x in path)
        nunitigs += 1
        nreads += len(path)
    logging.debug("A total of {0} unitigs built from {1} reads.".format(nunitigs, nreads))