def compare_graphs(graph):
    n = nx.number_of_nodes(graph)
    m = nx.number_of_edges(graph)
    k = np.mean(list(nx.degree(graph).values()))
    erdos = nx.erdos_renyi_graph(n, p=m/float(n*(n-1)/2))
    barabasi = nx.barabasi_albert_graph(n, m=int(k)-7)
    small_world = nx.watts_strogatz_graph(n, int(k), p=0.04)
    print(' ')
    print('Compare the number of edges')
    print(' ')
    print('My network: ' + str(nx.number_of_edges(graph)))
    print('Erdos: ' + str(nx.number_of_edges(erdos)))
    print('Barabasi: ' + str(nx.number_of_edges(barabasi)))
    print('SW: ' + str(nx.number_of_edges(small_world)))
    print(' ')
    print('Compare average clustering coefficients')
    print(' ')
    print('My network: ' + str(nx.average_clustering(graph)))
    print('Erdos: ' + str(nx.average_clustering(erdos)))
    print('Barabasi: ' + str(nx.average_clustering(barabasi)))
    print('SW: ' + str(nx.average_clustering(small_world)))
    print(' ')
    print('Compare average path length')
    print(' ')
    print('My network: ' + str(nx.average_shortest_path_length(graph)))
    print('Erdos: ' + str(nx.average_shortest_path_length(erdos)))
    print('Barabasi: ' + str(nx.average_shortest_path_length(barabasi)))
    print('SW: ' + str(nx.average_shortest_path_length(small_world)))
    print(' ')
    print('Compare graph diameter')
    print(' ')
    print('My network: ' + str(nx.diameter(graph)))
    print('Erdos: ' + str(nx.diameter(erdos)))
    print('Barabasi: ' + str(nx.diameter(barabasi)))
    print('SW: ' + str(nx.diameter(small_world)))

def algorithm(w1,w2,w3,w4,G1,G2,G3,G4):
	try:
		cc=np.array([nx.average_clustering(G1,weight='weight'),nx.average_clustering(G2,weight='weight'),nx.average_clustering(G3,weight='weight'),nx.average_clustering(G4,weight='weight')])
		spl=np.array([nx.average_shortest_path_length(G1,weight='weight'),nx.average_shortest_path_length(G2,weight='weight'),nx.average_shortest_path_length(G3,weight='weight'),nx.average_shortest_path_length(G4,weight='weight')])
		nds=np.array([nx.number_of_nodes(G1),nx.number_of_nodes(G2),nx.number_of_nodes(G3),nx.number_of_nodes(G4)])
		edgs= np.array([nx.number_of_edges(G1),nx.number_of_edges(G2),nx.number_of_edges(G3),nx.number_of_edges(G4)])
		if valid(cc):
			cc=stats.zscore(cc)
		else:
			cc=np.array([.1,.1,.1,.1])
		cc= cc-min(cc)+.1
		if valid(spl):
			spl=stats.zscore(spl)
		else:
			spl=np.array([.1,.1,.1,.1])
		spl= spl-min(spl)+.1
		if valid(nds):
			nds=stats.zscore(nds)
		else:
			nds=np.array([.1,.1,.1,.1])
		nds = nds-min(nds)+.1
		if valid(edgs):
			edgs=stats.zscore(edgs)
		else:
			edgs=np.array([.1,.1,.1,.1])
		edgs=edgs-min(edgs)+.1
		r1=(w1*cc[0]+w2*spl[0]+w3*nds[0]+w4*edgs[0])*1000
		r2=(w1*cc[1]+w2*spl[1]+w3*nds[1]+w4*edgs[1])*1000
		r3=(w1*cc[2]+w2*spl[2]+w3*nds[2]+w4*edgs[2])*1000
		r4=(w1*cc[3]+w2*spl[3]+w3*nds[3]+w4*edgs[3])*1000
		d={'Player 1:': r1, 'Player 2:': r2,'Player 3:': r3, 'Player 4:': r4}
		rank = sorted(d.items(), key=lambda x: x[1], reverse=True)
		return ["USAU RANKINGS",str(rank[0][0])+ " " + str(int(rank[0][1])),str(rank[1][0])+" "+ str(int(rank[1][1])),str(rank[2][0])+" "+ str(int(rank[2][1])),str(rank[3][0])+" "+str(int(rank[3][1]))]
	except:
		return ["Unable to compute rankings!  Need data","Player 1","Player 2","Player 3","Player 4"]

def run_main():
    file = str(sys.argv[1])
    f = open(file, 'r')
    print "\nReading inputfile:", file, "..."
    
    edgelist = []
    for line in f.readlines():
        edgelist.append((int(line.split()[0]), int(line.split()[1])))
    
    
    Directed_G = nx.DiGraph(edgelist)
    Undirected_G = Directed_G.to_undirected()
    #plt.figure(figsize=(8,8))
    #nx.draw(Directed_G,pos=nx.spring_layout(Directed_G))
    #plt.draw()
    #time.sleep(0.1)

    # compute other things
    print "Number of nodes involved in network:", nx.number_of_nodes(Undirected_G)
    print "Number of edges:", nx.number_of_edges(Undirected_G)
    print "Average degree:", nx.number_of_edges(Undirected_G) / float(nx.number_of_nodes(Undirected_G))
    t0 = time.clock()
    print "Average clustering coefficient:", compute_clustering_coefficient(Directed_G, Undirected_G)
    print "Took:", time.clock() - t0, "seconds"
    t1 = time.clock()
    print "Average path length:", average_shortest_path(Directed_G, Undirected_G)
    print "Took:", time.clock() - t1, "seconds"
    print "Total time:", time.clock() - t0, "seconds"
           
    report_final_stats()
    counter += 1
    second_counter += 1

def reduceGraph(read_g, write_g, minEdgeWeight, minNodeDegree, Lp, Sp):
    """
    Simplify the undirected graph and then update the 3 undirected weight properties.
    :param read_g: is the graph pickle to read
    :param write_g: is the updated graph pickle to write
    :param minEdgeWeight: the original weight of each edge should be >= minEdgeWeight
    :param minNodeDegree: the degree of each node should be >= minNodeDegree. the degree here is G.degree(node), NOT G.degree(node,weight='weight)
    :return: None
    """
    G=nx.read_gpickle(read_g)
    print 'number of original nodes: ', nx.number_of_nodes(G)
    print 'number of original edges: ', nx.number_of_edges(G)

    for (u,v,w) in G.edges(data='weight'):
        if w < minEdgeWeight:
            G.remove_edge(u,v)

    for n in G.nodes():
        if G.degree(n)<minNodeDegree:
            G.remove_node(n)

    print 'number of new nodes: ', nx.number_of_nodes(G)
    print 'number of new edges: ', nx.number_of_edges(G)

    for (a, b, w) in G.edges_iter(data='weight'):
        unweight_allocation(G, a, b, w,Lp,Sp)

    print 'update weight ok'
    nx.write_gpickle(G, write_g)

    return

def mcmc_subgraph_sample (master, eg):
	#import pdb; pdb.set_trace()
	
	new_graph = eg.copy()
	new_num_edges = curr_num_edges = num_eg_edges = nx.number_of_edges (eg)	
	iterations = 2 * nx.number_of_nodes (eg)    # play around w/ this number
	step_size = 1								# play around w/ this number
	sample_set = set(master.nodes_iter()) - set(eg.nodes_iter())
	coeff = 1
	for i in range(iterations):
		new_nodes = random.sample(sample_set, step_size) 
		old_nodes = random.sample(new_graph.nodes(), step_size)
		new_graph = replace_in_context(new_nodes, old_nodes, new_graph, master)
		new_num_edges = nx.number_of_edges(new_graph)
		new_stat = abs(new_num_edges - num_eg_edges) 
		old_stat = abs(curr_num_edges - num_eg_edges) 
		rval = random.random()
		if (new_stat <= old_stat) or (rval < math.exp(coeff * (old_stat - new_stat))):
			# use new graph
			curr_num_edges = new_num_edges
			sample_set = (sample_set | set(new_nodes)) - set(old_nodes) 
		else:
			# swap back old graph
			new_graph = replace_in_context(old_nodes, new_nodes, new_graph, master)
		if new_stat > num_eg_edges:
			coeff *= 1
	return new_graph.nodes_iter()

def test_union_all_and_compose_all():
    K3=nx.complete_graph(3)
    P3=nx.path_graph(3)

    G1=nx.DiGraph()
    G1.add_edge('A','B')
    G1.add_edge('A','C')
    G1.add_edge('A','D')
    G2=nx.DiGraph()
    G2.add_edge('1','2')
    G2.add_edge('1','3')
    G2.add_edge('1','4')

    G=nx.union_all([G1,G2])
    H=nx.compose_all([G1,G2])
    assert_edges_equal(G.edges(),H.edges())
    assert_false(G.has_edge('A','1'))
    assert_raises(nx.NetworkXError, nx.union, K3, P3)
    H1=nx.union_all([H,G1],rename=('H','G1'))
    assert_equal(sorted(H1.nodes()),
        ['G1A', 'G1B', 'G1C', 'G1D',
         'H1', 'H2', 'H3', 'H4', 'HA', 'HB', 'HC', 'HD'])

    H2=nx.union_all([H,G2],rename=("H",""))
    assert_equal(sorted(H2.nodes()),
        ['1', '2', '3', '4',
         'H1', 'H2', 'H3', 'H4', 'HA', 'HB', 'HC', 'HD'])

    assert_false(H1.has_edge('NB','NA'))

    G=nx.compose_all([G,G])
    assert_edges_equal(G.edges(),H.edges())

    G2=nx.union_all([G2,G2],rename=('','copy'))
    assert_equal(sorted(G2.nodes()),
        ['1', '2', '3', '4', 'copy1', 'copy2', 'copy3', 'copy4'])

    assert_equal(G2.neighbors('copy4'),[])
    assert_equal(sorted(G2.neighbors('copy1')),['copy2', 'copy3', 'copy4'])
    assert_equal(len(G),8)
    assert_equal(nx.number_of_edges(G),6)

    E=nx.disjoint_union_all([G,G])
    assert_equal(len(E),16)
    assert_equal(nx.number_of_edges(E),12)

    E=nx.disjoint_union_all([G1,G2])
    assert_equal(sorted(E.nodes()),[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])

    G1=nx.DiGraph()
    G1.add_edge('A','B')
    G2=nx.DiGraph()
    G2.add_edge(1,2)
    G3=nx.DiGraph()
    G3.add_edge(11,22)
    G4=nx.union_all([G1,G2,G3],rename=("G1","G2","G3"))
    assert_equal(sorted(G4.nodes()),
        ['G1A', 'G1B', 'G21', 'G22',
         'G311', 'G322'])

def validate_constituency_parse(tokenization):
    """
    Args:
      tokenization (concrete.structure.ttypes.Tokenization)

    Returns:
      bool: True if tokenization's constituency parse is valid, False otherwise
    """
    valid = True

    if tokenization.parse:
        total_constituents = len(tokenization.parse.constituentList)
        logging.debug(ilm(6, "tokenization '%s' has %d constituents" % (tokenization.uuid, total_constituents)))

        total_uuid_mismatches = 0
        constituent_id_set = set()
        constituent_parse_tree = nx.DiGraph()

        for constituent in tokenization.parse.constituentList:
            # Add nodes to parse tree
            constituent_parse_tree.add_node(constituent.id)

            if constituent.id not in constituent_id_set:
                constituent_id_set.add(constituent.id)
            else:
                valid = False
                logging.error(ilm(7, "constituent ID %d has already been used in this sentence's tokenization" % constituent.id))

            # Per the Concrete 'structure.thrift' file, tokenSequence may not be defined:
            #   "Typically, this field will only be defined for leaf constituents (i.e., constituents with no children)."
            if constituent.tokenSequence and constituent.tokenSequence.tokenizationId != tokenization.uuid:
                total_uuid_mismatches += 1

        if total_uuid_mismatches > 0:
            valid = False
            logging.error(ilm(6, "tokenization '%s' has UUID mismatch for %d/%d constituents" %
                              (tokenization.uuid, total_uuid_mismatches, total_constituents)))

        # Add edges to constituent parse tree
        for constituent in tokenization.parse.constituentList:
            if constituent.childList:
                for child_id in constituent.childList:
                    constituent_parse_tree.add_edge(constituent.id, child_id)

        # Check if constituent parse "tree" is actually a tree
        undirected_graph = constituent_parse_tree.to_undirected()
        if not nx.is_connected(undirected_graph):
            valid = False
            logging.error(ilm(6, "The constituent parse \"tree\" is not a fully connected graph - the graph has %d components" %
                len(nx.connected_components(undirected_graph))))
        if nx.number_of_nodes(constituent_parse_tree) != nx.number_of_edges(constituent_parse_tree) + 1:
            valid = False
            logging.error(ilm(6, "The constituent parse \"tree\" is not a tree.  |V| != |E|+1  (|V|=%d, |E|=%d)" %
                (nx.number_of_nodes(constituent_parse_tree), nx.number_of_edges(constituent_parse_tree))))

    return valid

def get_number_of_edges(filename):
  import networkx as nx
  threshold = 0
  f = open(filename[:-4]+'_edges.dat','w')
  for i in range(0,101):
    threshold = float(i)/100
    G = get_threshold_matrix(filename, threshold)
    print 'number of edges:', nx.number_of_edges(G)
    max_number_edges = nx.number_of_nodes(G) * (nx.number_of_nodes(G) - 1.) / 2
    f.write("%f\t%d\t%f\n" % (threshold, nx.number_of_edges(G), nx.number_of_edges(G)/max_number_edges))
  f.close()

def eval_proximity_vertices(network,graph_xml) :
    '''returns the proximity of proportion of vertices between synthetic network(test) and real network (goal)'''
    number_of_nodes_test = float(nx.number_of_nodes(network))
    if network.isDirected() : 
        proportion_edges_test = nx.number_of_edges(network)/(number_of_nodes_test*(number_of_nodes_test-1))
    else :
        proportion_edges_test = 2.*nx.number_of_edges(network)/(number_of_nodes_test*(number_of_nodes_test-1))
    
    proportion_edges_goal = eval(graph_xml.find('vertices').get('value'))
    proximity = proximity_numbers(proportion_edges_goal,proportion_edges_test )
    return proximity

	def reformat2Igraph(self,graph):
		print nx.number_of_nodes(graph)
		print nx.number_of_edges(graph)
		G=Graph(0)
		for i in range(nx.number_of_nodes(graph)):
			G.add_vertices(1)
		for i in graph.edge:
			for j in graph.edge[i]:
				if i<=j:
					G.add_edges([(i,j)])
		return G

def modularity(subgs, G):
    Q = 0
    total_edges = float(nx.number_of_edges(G))

    for g in subgs:
        nodes = g.node.keys()
        degree_sum = sum(nx.degree(G, nodes).values())
        edges_num = nx.number_of_edges(g)

        Q += edges_num / total_edges - (degree_sum / (2 * total_edges))**2

    return Q

def test_union_all_and_compose_all():
    K3 = nx.complete_graph(3)
    P3 = nx.path_graph(3)

    G1 = nx.DiGraph()
    G1.add_edge("A", "B")
    G1.add_edge("A", "C")
    G1.add_edge("A", "D")
    G2 = nx.DiGraph()
    G2.add_edge("1", "2")
    G2.add_edge("1", "3")
    G2.add_edge("1", "4")

    G = nx.union_all([G1, G2])
    H = nx.compose_all([G1, G2])
    assert_edges_equal(G.edges(), H.edges())
    assert_false(G.has_edge("A", "1"))
    assert_raises(nx.NetworkXError, nx.union, K3, P3)
    H1 = nx.union_all([H, G1], rename=("H", "G1"))
    assert_equal(sorted(H1.nodes()), ["G1A", "G1B", "G1C", "G1D", "H1", "H2", "H3", "H4", "HA", "HB", "HC", "HD"])

    H2 = nx.union_all([H, G2], rename=("H", ""))
    assert_equal(sorted(H2.nodes()), ["1", "2", "3", "4", "H1", "H2", "H3", "H4", "HA", "HB", "HC", "HD"])

    assert_false(H1.has_edge("NB", "NA"))

    G = nx.compose_all([G, G])
    assert_edges_equal(G.edges(), H.edges())

    G2 = nx.union_all([G2, G2], rename=("", "copy"))
    assert_equal(sorted(G2.nodes()), ["1", "2", "3", "4", "copy1", "copy2", "copy3", "copy4"])

    assert_equal(G2.neighbors("copy4"), [])
    assert_equal(sorted(G2.neighbors("copy1")), ["copy2", "copy3", "copy4"])
    assert_equal(len(G), 8)
    assert_equal(nx.number_of_edges(G), 6)

    E = nx.disjoint_union_all([G, G])
    assert_equal(len(E), 16)
    assert_equal(nx.number_of_edges(E), 12)

    E = nx.disjoint_union_all([G1, G2])
    assert_equal(sorted(E.nodes()), [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])

    G1 = nx.DiGraph()
    G1.add_edge("A", "B")
    G2 = nx.DiGraph()
    G2.add_edge(1, 2)
    G3 = nx.DiGraph()
    G3.add_edge(11, 22)
    G4 = nx.union_all([G1, G2, G3], rename=("G1", "G2", "G3"))
    assert_equal(sorted(G4.nodes()), ["G1A", "G1B", "G21", "G22", "G311", "G322"])

    def test_equivalence_transform(self, ch2, ch3, methane):
        ch2_atoms = list(ch2.particles())
        methane_atoms = list(methane.particles())
        equivalence_transform(ch2, ch2_atoms[0], methane_atoms[0], add_bond=False)
        assert (ch2_atoms[0].pos == methane_atoms[0].pos).all()
        equivalence_transform(ch2, ch2['up'], ch3['up'])
        assert ch2.n_bonds == 2

        assert nx.number_of_edges(ch2.root.bond_graph) == 3
        assert nx.number_of_edges(ch3.root.bond_graph) == 4

        ethyl = mb.Compound([ch2, ch3])
        assert ethyl.n_bonds == 6

def main(DAG, cite, out):
    DAG = create(cite, DAG)
    print nx.number_of_edges(DAG)
   # largest_component = component(DAG) #uses component function to select the largest subgraph in the data
   # if not largest_component.is_directed():
   #     out.write('\n' + 'yup not directed')
   # out.write('\n' + 'There are %d nodes in the largest component' %nx.number_of_nodes(largest_component))
   # largest_component_DAG = redirect(DAG,largest_component)
   # if not largest_component_DAG.is_directed():
   #     out.write('\n' + 'Not directed')
   # else:
   #     out.write('\n' + 'Directed!')
   # out.write('\n' + 'There are now %d nodes in the largest component' %nx.number_of_nodes(largest_component_DAG))
    #check if it's connected etc!
   # out.write(str(DAG.number_of_edges()))
   # if not nx.is_connected(largest_component_DAG.to_undirected()):
   #     out.write('\n'+'The network is disconnected')
   #     out.write('\n' + "There are %d connected components" %nx.number_connected_components(largest_component_DAG.to_undirected()))
    #out.write(str(nx.average_shortest_path_length(largest_component_DAG)))
    print age('0001001', '9401139')
    print age('0103030', '0204161')
    print simple_age('33', '100')
    print simple_age('45', '2')
    paths = find_all_paths(DAG,'8', '0')
 #   print paths
    empty = []
    listofpaths = make_list_of_paths(paths, empty)
    pathset = []
    for path in listofpaths:
        if path not in pathset:
            pathset.append(path)
    for path in pathset:
        out.write('\n' + str(path))
    out.write('\n' + 'The longest path is ' + str(max(pathset, key=len)))
#    setofpaths = set(listofpaths)

        #Uses method 2 of assigning a number to each node reflecting the path length between that node and the youngest node
        print 'Testing using numbering started'  
        paths3 = distance_to_all_nodes(box, extremal_points[0] , n) #creates list containing lists of lists of lists...containing a list of the nodes in path
        empty3 = [] #creates the list which will contain all of the paths
        listofpaths3 = make_list_of_paths(paths3, empty3) #looks through 'paths' to extract proper paths from the various levels of sublists
        pathset3 = [] #creates the list which will contain all unique paths
        for path in listofpaths3:
            if path not in pathset3:
                pathset3.append(path) #only adds unique paths to pathset
        out.write('\n' + 'Path testing using number...%d unique paths found' %(len(pathset3)))
        for path in pathset3:
            out.write('\n' + str(path)) #prints all the unique paths
#        longestpath3 = max(pathset3, key=len) #identifies the longest path in the set of paths
#        out.write('\n' + 'The longest path from numbering is %s which is %d nodes long' %(str(longestpath3), len(longestpath3)))
        print 'Testing using numbering completed'  

def compareNumberOfEdges(masterGraph,wordGraph,worksheet,row):
    numberOfEdgesMasterGraph = nx.number_of_edges(masterGraph)
    numberOfEdgesWordGraph = nx.number_of_edges(wordGraph)
#    worksheet.write(row,1,numberOfEdgesMasterGraph)
#    worksheet.write(row,2,numberOfEdgesWordGraph)
    result = False
    if(numberOfEdgesMasterGraph >= numberOfEdgesWordGraph):
        result = True
#    worksheet.write(row,3,result)

    if result == True:
        return 1
    else:
        return -1

def split_edges(G):
    Gsmaller = nx.DiGraph()
    Glarger = nx.DiGraph()
    for node in G.nodes():
        Gsmaller.add_node(node)
        Glarger.add_node(node)
    for (u,v) in G.edges():
        if G.node[u]['index'] < G.node[v]['index']:
            Gsmaller.add_edge(u,v)
        elif G.node[u]['index'] > G.node[v]['index']:
            Glarger.add_edge(u,v)
    if nx.number_of_edges(Gsmaller) > nx.number_of_edges(Glarger):
        return topological_sort(Gsmaller)
    else:
        return topological_sort(Glarger)

        def generate_smtcode_old(self): #convert_mpnf
            st = time.time()
            self.mdg = get_weighted_graph(self.graph, self.pair_mp_prop, self.opts)

            self.info.extend([('g_nodes', nx.number_of_nodes(self.graph)),
                                    ('g_edges', nx.number_of_edges(self.graph)),
                                    ('mdg_nodes', nx.number_of_nodes(self.mdg)),
                                    ('mdg_edges', nx.number_of_edges(self.mdg))])

            sccs = self._gen_sccs(self.mdg)

            self._gen_smt_from_graph(sccs)

            en = time.time()
            self.time.append(('total gen_smt: graph->mdg->mdg(scc)', en - st))

def get_SMTCode(mdg, mp_prop, opts):
    # input mdg
    # output SMTCode list, result per 
    results = {'SCCs':[], 'AccSCCs':[]}
    if opts.tgba:
        all_acceptlist = mdg.graph['acccond']
    SCCs = nx.strongly_connected_component_subgraphs(mdg)
    for scc in SCCs:
        results['SCCs'].append({'edges':nx.number_of_edges(scc), 'nodes':nx.number_of_nodes(scc)})
        if scc.edges() == []:
            pass
        else:
            # nx.draw_circular(scc)
            accept = 0
            if opts.tgba:
                a = reduce(lambda x, y : x | y , [set(edge[2]['acc']) for edge in scc.edges(data=True)])
                print a
                acceptlist = list(a)
                print acceptlist
                """for edge in scc.edges(data=True):
                    for acc in edge[2]['acc']:
                        if not acc in acceptlist:
                            acceptlist.append(acc)
                    if opts.debug:
                        print acceptlist
                        """
                if all_acceptlist.sort() == acceptlist.sort():
                    if opts.pdebug:
                        print acceptlist,' ==? ',all_acceptlist
                    print '  TGBA-Accept'
                    accept = 1
            else:
                for node in scc.nodes():
                    if 'accept' in node:
                        accept = 1
                        break

            if accept == 1:
                if opts.debug:
                    print '\t found BA/TGBA-Accept SCC'
                code_time = runsmt.gen_smtcodelist_time(scc, mp_prop, opts)
                results['AccSCCs'].append({'SMTCode':code_time['code'],
                                           'edges':nx.number_of_edges(scc), 'nodes':nx.number_of_nodes(scc),
                                           'gen_time':code_time['time'], 'sat_time':0})
            else:
                if opts.pdebug:
                    print ' non accept SCC'
    return results

def calculate(g, voltage):
    edges_num = nx.number_of_edges(g)
    # sort nodes in edges
    edges = [edge if edge[0] < edge[1] else (edge[1], edge[0])
             for edge in nx.edges(g)]

    a = np.zeros((edges_num, edges_num))
    b = np.zeros((edges_num, 1))
    i = 0
    # first law
    for node in [node for node in nx.nodes(g) if node != 0]:
        for neighbor in nx.all_neighbors(g, node):
            edge = tuple(sorted((node, neighbor)))
            a[i][edges.index(edge)] = 1 if neighbor < node else -1
        i += 1
    # second law
    cycles = nx.cycle_basis(g, 0)
    for cycle in cycles:
        for j in range(0, len(cycle)):
            node = cycle[j]
            next_node = cycle[(j + 1) % len(cycle)]
            edge = tuple(sorted((node, next_node)))
            resistance = g[node][next_node]['weight']
            a[i][edges.index(edge)] = resistance if node < next_node else -resistance
        if 0 in cycle:
            b[i] = voltage
        i += 1
    # solve
    x = np.linalg.solve(a, b)
    for (x1, x2), res in zip(edges, x):
        g[x1][x2]['current'] = res[0]

def modularity(graph, membership):
    try:
        noOfEdges = nx.number_of_edges(graph)
        types = max(membership) + 1
        print types, len(membership)
        a = {}
        e = {}
        for i in range(types):
            a[i] = 0
            e[i] = 0
        for edge in graph.edges():
            fromNode = int(edge[0])
            toNode = int(edge[1])
            c1 = membership[fromNode]
            c2 = membership[toNode]
            if c1 == c2:
                e[c1] += 2
            # print "c1= ", c1, "c2=", c2
            a[c1] += 1
            a[c2] += 1
            modularity = 0.0
            if noOfEdges > 0:
                for i in range(types):
                    tmp = a[i] / 2 / noOfEdges
                    modularity += e[i] / 2 / noOfEdges
                    modularity -= tmp * tmp
    except Exception as e:
        exc_type, exc_obj, exc_tb = sys.exc_info()
        fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
        print (exc_type, fname, exc_tb.tb_lineno)
        raise e
    return modularity