def test_random_seed(self):
        """Tests that each call with the same random seed generates the
        same graph.

        """
        deg_seq = [3] * 12
        G1 = nx.configuration_model(deg_seq, seed=1000)
        G2 = nx.configuration_model(deg_seq, seed=1000)
        assert_true(nx.is_isomorphic(G1, G2))
        G1 = nx.configuration_model(deg_seq, seed=10)
        G2 = nx.configuration_model(deg_seq, seed=10)
        assert_true(nx.is_isomorphic(G1, G2))

def generateNetwork():
    ## use a networkx function to create a degree sequence that follows a power law
    degreeSequence=nx.utils.create_degree_sequence(numberOfNodes,powerlaw_sequence, 100)
    ## use aforementioned degree sequence to configure a pseudograph that contains self-loops & hyper-edges
    pseudoGraph=nx.configuration_model(degreeSequence)
    ## remove hyper (parallel) edges
    Graph = nx.Graph(pseudoGraph)
    ## remove self edges
    Graph.remove_edges_from(Graph.selfloop_edges())
    ## loop through all nodes and set capacity equal to degree
    for bankID in range(0, len(Graph.node)):
        Graph.node[bankID]['bankID'] = bankID
        Graph.node[bankID]['capacity'] = Graph.degree(bankID)
        Graph.node[bankID]['solventNeighbors'] = Graph.degree(bankID)
        ## right now capacity = degree
        Graph.node[bankID]['cumulativeShock'] = 0
        ## solvent = normal
        ## exposed = cumulative shock is less than capacity
        ## fail = recently failed, about to spread
        ## dead = can no longer spread or receive shocks
        Graph.node[bankID]['status'] = 'solvent'
        ## here we set the timestep that the bank becomes insolvent to a big number
        Graph.node[bankID]['insolventTimestep'] = 50
        ## here we set the size of the shock to be propagated (zero at sim start)
        Graph.node[bankID]['shockToPropagate'] = 0
    return Graph

    def createGraph(self):
        """ Calculate the basic data for the graph """
        """ Generate everything for the scale free network """

        # Create a graph with degrees following a power law distribution

        s = []

        count = 0

        while len(s) < self.N:
            nextval = int(nx.utils.powerlaw_sequence(int(self.k), self.e)[0])
            
            if nextval != 0:
                count += nextval
                s.append(nextval)
                
        # s scaled and rounded such that the average degree equals k
        s = s / np.mean(s) * self.k
        s = np.around(s).astype(int)

        # Sum of degrees must be even. I added one edge to the first node to fix this
        if sum(s) % 2:
            s[0] += 1
            
        G = nx.configuration_model(s)
        G = nx.Graph(G)
           
        # Remove self-loops
        G.remove_edges_from(G.selfloop_edges())
            
        self.G = G
        
        self.generateInfected()

def get_corr_rand_set(G,disease_seeds,num_reps=5,alpha=.5,num_its=20,conserve_heat=True):
    '''
    
    Calculate the dot-product of heat propagated on N disease gene sets (disease_seeds: dict with keys disease names and values lists of disease genes), on an edge-shuffled, degree-preserving random matrix, with number of repetitions = num_reps, alpha=alpha, num_its = num_its.
    
    Return the mean of the dot-product averaged over num_reps, and the standard deviation over num_reps, over all pairs of gene sets in disease_seeds.  This way we only have to create one random matrix for each pair, which will speed up processing time a bit.
    
    '''
    
    num_Ds = len(disease_seeds)
    
    dnames = disease_seeds.keys()
    
    dname_pairs = list(itertools.combinations(dnames, 2))
        
    
    dot_rand=dict()
    dot_rand_mean = dict()
    dot_rand_std = dict()
    for d in dname_pairs:
        # initialize dictionaries
        dot_rand[d] = []
        dot_rand_mean[d] = []
        dot_rand_std[d] = []
    
    for r in range(num_reps):
        G_temp = nx.configuration_model(G.degree().values())
        G_rand = nx.Graph()  # switch from multigraph to digraph
        G_rand.add_edges_from(G_temp.edges())
        # remove self-loops
        #G_rand.remove_edges_from(G_rand.selfloop_edges())
        G_rand = nx.relabel_nodes(G_rand,dict(zip(range(len(G_rand.nodes())),G.degree().keys())))
        Wprime_rand = normalized_adj_matrix(G_rand,conserve_heat=conserve_heat)
        
        
        for i in range(len(dname_pairs)):
            seeds_D1 = disease_seeds[dname_pairs[i][0]]
            seeds_D2 = disease_seeds[dname_pairs[i][1]]
        
            Fnew_D1 = network_propagation(G_rand,Wprime_rand,seeds_D1,alpha=alpha,num_its=num_its)
            Fnew_D1_norm = Fnew_D1/np.linalg.norm(Fnew_D1)
        
            rand_seeds = seeds_D2 #set(random.sample(G.nodes(),size_rand_set))
            Fnew_D2 = network_propagation(G_rand,Wprime_rand,seeds_D2,alpha=alpha,num_its=num_its)
            Fnew_D2_norm = Fnew_D2/np.linalg.norm(Fnew_D2)

            idx_g0 = list(Fnew_D1[(Fnew_D1>0)&(Fnew_D2>0)].index)
            idx_ND1ND2 = list(np.setdiff1d(list(Fnew_D1.index),np.union1d(seeds_D1,seeds_D2)))
            
            dot_D1_D1 = np.dot(Fnew_D1_norm,Fnew_D2_norm)
                
            dot_rand[dname_pairs[i]].append(dot_D1_D1)
        
    
    for d in dname_pairs:
        dot_rand_mean[d] = np.mean(dot_rand[d])
        dot_rand_std[d] = np.std(dot_rand[d])

    
    return dot_rand_mean,dot_rand_std

def scale_free_network(n=100, gamma=2.5, avrdeg=8):
    """
    Generates a scale free network by configuration model,
    which is shown in :cite:`Wu2011`

    :param int   n:      The number of nodes
    :param float gamma:  Exponent of degree distribution
    :param float avrdeg: Average degree
    :return: Scale-free network
    :rtype: networkx.Graph

    Generates degree sequence from calling
    :func:`~.powerlaw_sequence` and
    :func:`networkx.utils.random_sequence.create_degree_sequence`
    then, build network by configuration model
    (see :func:`networkx.generators.degree_seq.configuration_model`)
    finally, remove self-loop and parallel edge to simplify.

    .. Note:: Actual average degree is smaller than `avrdeg`
    """
    seq = create_degree_sequence(
        n, powerlaw_sequence, gamma=gamma, avrdeg=avrdeg)
    G = configuration_model(seq)
    # multigraph -> graph
    G.remove_edges_from(G.selfloop_edges())
    return nx.Graph(G)

def generate_graph(type = 'PL', n = 100, seed = 1.0, parameter = 2.1):
    if type == 'ER':
        G = nx.erdos_renyi_graph(n, p=parameter, seed=seed, directed=True)
        G = nx.DiGraph(G)
        G.remove_edges_from(G.selfloop_edges())
    elif type == 'PL':
        z = nx.utils.create_degree_sequence(n, nx.utils.powerlaw_sequence, exponent = parameter)
        while not nx.is_valid_degree_sequence(z):
            z = nx.utils.create_degree_sequence(n, nx.utils.powerlaw_sequence, exponent = parameter)
        G = nx.configuration_model(z)
        G = nx.DiGraph(G)
        G.remove_edges_from(G.selfloop_edges())
    elif type == 'BA':
        G = nx.barabasi_albert_graph(n, 3, seed=None)
        G = nx.DiGraph(G)
    elif type == 'grid':
        G = nx.grid_2d_graph(int(np.ceil(np.sqrt(n))), int(np.ceil(np.sqrt(n))))
        G = nx.DiGraph(G)
    elif type in ['facebook', 'enron', 'twitter', 'students', 'tumblr', 'facebookBig']:
        #print 'start reading'
        #_, G, _, _ = readRealGraph(os.path.join("..","..","Data", type+".txt"))
        _, G, _, _ = readRealGraph(os.path.join("..","Data", type+".txt"))
        print 'size of graph', G.number_of_nodes()
        #Gcc = sorted(nx.connected_component_subgraphs(G.to_undirected()), key = len, reverse=True)
        #print Gcc[0].number_of_nodes()
        #print 'num of connected components', len(sorted(nx.connected_component_subgraphs(G.to_undirected()), key = len, reverse=True))
        #exit()
        if G.number_of_nodes() > n:
            G = getSubgraph(G, n)
        #G = getSubgraphSimulation(G, n, infP = 0.3)
    #nx.draw(G)
    #plt.show()
    return G

def normalise(G, func, n=500, retNorm=True, inVal=None, printLCC=False):
    """
    This function normalises the function G by generating a series of n random
    graphs and averaging the results. If retNorm is specified, the normalised 
    value is returned, else a list of n values for a random graph are returned.
    """
    vals = []
    for i in range(n):
        rand = configuration_model(G.degree().values())
        if nx.components.number_connected_components(rand)>1:
            graphList = ccs(rand)
            rand = graphList.next()
            if printLCC:
                print "Selecting largest connected components of "+str(len(rand.nodes()))+" nodes"
        try:
            vals.append(func(rand)) # collect values using the function
        except KeyError: # exception raised if the spatial information is missing
            nx.set_node_attributes(rand, 'xyz', {rn:G.node[v]['xyz'] for rn,v in enumerate(G.nodes())})
            vals.append(func(rand)) # collect values using the function

    if retNorm: # return the normalised values
        if not inVal:
            inVal = func(G)
        return(inVal/np.mean(vals))
        
    else: # return a list of the values from the random graph
        return(vals)

	def initScaleFreeGraph(self,N,k,e):
		s = []
		while len(s) < N:
		    nextval = int(nx.utils.powerlaw_sequence(k, e)[0])
		    if nextval != 0:
		        s.append(nextval)

		#TODO: s must be scaled
		s = s/np.mean(s) * k
		s = np.around(s)
		s = s.astype(int)

		# Sum of degrees must be even. One way to solve this is to add a single node with degree 1
		if sum(s) % 2:
			s[len(s)-1] = s[len(s)-1] + 1

		print "length s:", len(s)
		print "Sum:", sum(s)

		G = nx.configuration_model(s)
		G = nx.Graph(G)
		
		# Remove self-loops
		G.remove_edges_from(G.selfloop_edges())
		return G

def normaliseNodeWise(G, func, n=500, retNorm=True, inVal=None):
    """
    This function normalises the function G by generating a series of n random
    graphs and averaging the results. If retNorm is specified, the normalised 
    value is returned, else a list of n values for a random graph are returned.
    """
    nodesDict = {v:[] for v in G.nodes()}
    for i in range(n):
        rand = configuration_model(G.degree().values())
        rand = nx.Graph(rand) # convert to simple graph from multigraph
        nodeList = [v for v in rand.nodes() if rand.degree(v)==0]
        rand.remove_nodes_from(nodeList)
        
        try:
            res = func(rand) # collect values using the function
        except KeyError: # exception raised if the spatial information is missing
            print "Adding spatial info"
            nx.set_node_attributes(rand, 'xyz', {rn:G.node[v]['xyz'] for rn,v in enumerate(G.nodes())}) # copy across spatial information
            res = func(rand) # collect values using the function
            
        for x,node in enumerate(nodesDict):
            nodesDict[node].append(res[x])

    if retNorm: # return the normalised values
        if not inVal:
            inVal = func(G)
        for node in nodesDict:
            nodesDict[node] = inVal[node]/np.mean(nodesDict[node])
        return(nodesDict)
        
    else: # return a list of the values from the random graph
        return(nodesDict)

def test_configuration():
    seeds = [2718183590, 2470619828, 1694705158, 3001036531, 2401251497]
    for seed in seeds:
        deg_seq = nx.random_powerlaw_tree_sequence(20, seed=seed, tries=5000)
        G = nx.Graph(nx.configuration_model(deg_seq, seed=seed))
        G.remove_edges_from(nx.selfloop_edges(G))
        _check_augmentations(G)

def powerlaw_degree_sequence(n, a):
    """
    Create a graph without self-loops or parallel edges; having power law degree
    distribution with an exponent 'around' a.
    
    Parameters
    ----------

    n : int
       Number of nodes in graph.
       
    a : float
       Ideal exponent.

    Returns
    -------
    
    G : networkx.Graph
       The constructed graph.

    """
    dsq = nx.create_degree_sequence(n, nx.utils.powerlaw_sequence, exponent = a)
    G = nx.Graph(nx.configuration_model(dsq))
    G.remove_edges_from(G.selfloop_edges())
    
    # Check for a disconnected graph just in case...
    if not nx.is_connected(G):
        emsg = "The generated power-law graph is not connected!"
        logger.error(emsg)
        raise NepidemiXBaseException(emsg)

    return G

def randomly_clustering(g, tries = 10):
	"""
		Comparing the average clustering coefficient of g with other graphs h
		which share identical degree sequence. This function returns the comparison ratio.

		Parameters:
		-----------
			g: NetworkX Graph, NetworkX DiGraph
			tries: int, optional, (default = 10)
				number of tries (compared graphs)
		See also:
		---------
			mean_clustering
		Returns:
		--------
			float, the ratio of avg clustering coefficient, avg_cc(g) / mean(avg_cc(h))
	"""
	from scipy import average
	g = to_undirected(g)
	d = g.degree().values()
	c = mean_clustering(g, normalized = False)
	p = list()
	for t in xrange(tries):
		ng = nx.configuration_model(d, create_using = nx.Graph())
		p.append(mean_clustering(ng))
		del ng
	return(c / average(p))

    def test_degree_zero(self):
        """Tests that a degree sequence of all zeros yields the empty
        graph.

        """
        G = nx.configuration_model([0, 0, 0])
        assert_equal(len(G), 3)
        assert_equal(G.number_of_edges(), 0)

def test_configuration_directed():
    # seeds = [671221681, 2403749451, 124433910, 672335939, 1193127215]
    seeds = [67]
    for seed in seeds:
        deg_seq = nx.random_powerlaw_tree_sequence(20, seed=seed, tries=5000)
        G = nx.DiGraph(nx.configuration_model(deg_seq, seed=seed))
        G.remove_edges_from(nx.selfloop_edges(G))
        _check_edge_connectivity(G)

def measure_connectivity(G,focal_nodes,method='jaccard',num_reps=10):
    avg_focal_degree = np.mean(G.degree(focal_nodes).values())
    var_focal_degree = np.std(G.degree(focal_nodes).values())
    
    
    
    if method=='jaccard':
        
        #focal_sim = jaccard_sim_array(G,focal_nodes)
        
        # select random node set
        rand_sim=[]
        focal_sim=[]
        for r in range(num_reps):
            
            # bootstrap a sample from focal_nodes
            focal_subset = np.random.choice(list(focal_nodes),size=len(focal_nodes),replace=True)
            
            focal_sim.append(jaccard_sim_array(G,focal_subset))
            
            print('calculating random set ' + str(r) + ' out of ' + str(num_reps))
            G_temp = nx.configuration_model(G.degree().values())
            G_rand = nx.Graph()  # switch from multigraph to digraph
            G_rand.add_edges_from(G_temp.edges())
            # remove self-loops
            #G_rand.remove_edges_from(G_rand.selfloop_edges())
            G_rand = nx.relabel_nodes(G_rand,dict(zip(range(len(G_rand.nodes())),G.degree().keys())))
            
            rand_sim.append(jaccard_sim_array(G_rand,focal_subset)) # measure the jaccard similarity of degree-preserving edge shuffled network
            
    elif method=='edge_overlap':
        
        focal_sim = num_shared_neighbors(G,focal_nodes)
        # select random node set
        rand_sim=[]
        for r in range(num_reps):
            print('calculating random set ' + str(r) + ' out of ' + str(num_reps))
            G_rand = nx.configuration_model(G.degree().values())
            G_rand = nx.relabel_nodes(G_rand,dict(zip(range(len(G_rand.nodes())),G.degree().keys())))
            
            rand_sim.append(num_shared_neighbors(G_rand,focal_nodes)) # rand_sim is array of length num_reps
        
        
    
    return focal_sim,rand_sim

def generate_random_graph_powerlaw(n):
	m = 1
	p = 0.7
	#H = powerlaw_cluster_graph(n, m, p)

	gamma = 2.5
	powerlaw_gamma = lambda N: nx.utils.powerlaw_sequence(N, exponent=gamma)
	z = nx.utils.create_degree_sequence(n, powerlaw_gamma, max_tries=1000)
	G = nx.configuration_model(z)

	while not nx.is_connected(G):
		powerlaw_gamma = lambda N: nx.utils.powerlaw_sequence(N, exponent=gamma)
		z = nx.utils.create_degree_sequence(n, powerlaw_gamma, max_tries=100)
		G = nx.configuration_model(z)

	G = [n for n in nx.connected_component_subgraphs(G)][0]
	G.remove_edges_from(G.selfloop_edges())

	return G

def calc_localization(Gint,genes_focal,write_file_name='localization_results',num_reps=5,num_genes=20,
                     conserve_heat=True, replace=True,subsample=True, savefile=True):
    
    seed_FOCAL = list(np.intersect1d(list(genes_focal),Gint.nodes()))
    
    if subsample:
        num_genes_S=num_genes
    else:
        num_genes_S=len(seed_FOCAL)

    Wprime = normalized_adj_matrix(Gint,conserve_heat=conserve_heat)
    
    kurt_FOCAL =[]
    kurt_Srand=[]
    var_FOCAL, var_Srand=[],[]
    sumTop_FOCAL= []
    sumTop_Srand= []
    for r in range(num_reps):
        print(r)
        
        subset_FOCAL = np.random.choice(seed_FOCAL,size=num_genes_S,replace=replace)

        Fnew_FOCAL = network_propagation(Gint,Wprime,subset_FOCAL,alpha=.5,num_its=20)
        Fnew_FOCAL.sort()
        kurt_FOCAL.append(scipy.stats.kurtosis(Fnew_FOCAL))
        var_FOCAL.append(np.var(Fnew_FOCAL))
        sumTop_FOCAL.append(np.sum(Fnew_FOCAL.head(1000)))

        G_temp = nx.configuration_model(Gint.degree().values())
        G_rand = nx.Graph()  # switch from multigraph to digraph
        G_rand.add_edges_from(G_temp.edges())
        # remove self-loops
        #G_rand.remove_edges_from(G_rand.selfloop_edges())
        G_rand = nx.relabel_nodes(G_rand,dict(zip(range(len(G_rand.nodes())),Gint.degree().keys())))
        Wprime_rand = normalized_adj_matrix(G_rand,conserve_heat=conserve_heat)

        Fnew_Srand = network_propagation(G_rand,Wprime_rand,subset_FOCAL,alpha=.5,num_its=20)
        Fnew_Srand.sort()
        kurt_Srand.append(scipy.stats.kurtosis(Fnew_Srand))
        var_Srand.append(np.var(Fnew_Srand))
        sumTop_Srand.append(np.sum(Fnew_Srand.head(1000)))

        print(var_FOCAL[-1])
        print(var_Srand[-1])
    
    results_dict = {'kurtosis':kurt_FOCAL,'kurt_rand':kurt_Srand,
                   'var':var_FOCAL,'var_rand':var_Srand,
                   'sumTop':sumTop_FOCAL, 'sumTop_rand':sumTop_Srand,
                   'num_reps':num_reps, 'conserve_heat':conserve_heat,
                   'replace':replace,'subsample':subsample,'num_genes':num_genes}
    
    if savefile:
        json.dump(results_dict,open(write_file_name,'w'))
        
    return results_dict

def generate_graph(gamma, kmin, kmax, dk):
    # map probability to power law:
    p = np.random.random(graph_size)
    C = (-gamma+1) / ((kmax+dk)**(-gamma+1) - (kmin+dk)**(-gamma+1))
    k = ((-gamma+1)/C * p + (kmin+dk)**(-gamma + 1)) ** (1/(-gamma+1)) - dk
    sequence = map(int, sorted(np.round(k)))                    # map sorted intervalls to int
    if np.mod(np.sum(sequence), 2) != 0:                        # if odd, get rid of last element
        sequence[-1]-=1                                 
    G = nx.configuration_model(sequence)                
    G = nx.Graph(G)                                             # remove multiple connections
    return G

    def test_degree_sequence(self):
        """Tests that the degree sequence of the generated graph matches
        the input degree sequence.

        """
        deg_seq = [5, 3, 3, 3, 3, 2, 2, 2, 1, 1, 1]
        G = nx.configuration_model(deg_seq, seed=12345678)
        assert_equal(sorted((d for n, d in G.degree()), reverse=True),
                     [5, 3, 3, 3, 3, 2, 2, 2, 1, 1, 1])
        assert_equal(sorted((d for n, d in G.degree(range(len(deg_seq)))),
                            reverse=True),
                     [5, 3, 3, 3, 3, 2, 2, 2, 1, 1, 1])

 def connect_configuration(self,N,degseq):
     """
     Uses the configuration model to wire up a pulse oscillator network with a given
     input degree sequence. The length of degseq must equal N.  The resulting network
     is pruned of both self loops and parallel (multi) edges.
     """
     assert len(degseq) == N, "ERROR. Each node needs an input degree."
     self.connect_empty(N)
     G = nx.configuration_model(degseq)
     G = nx.Graph(G)
     G.remove_edges_from(G.selfloop_edges())
     self.add_edges_from(G.edges())