def draw_degree_rank_plot(orig_g, mG):
    ori_degree_seq = sorted(nx.degree(orig_g).values(), reverse=True)  # degree sequence
    deg_seqs = []
    for newg in mG:
        deg_seqs.append(sorted(nx.degree(newg).values(), reverse=True))  # degree sequence
    df = pd.DataFrame(deg_seqs)

    plt.xscale("log")
    plt.yscale("log")
    plt.fill_between(df.columns, df.mean() - df.sem(), df.mean() + df.sem(), color="blue", alpha=0.2, label="se")
    h, = plt.plot(df.mean(), color="blue", aa=True, linewidth=4, ls="--", label="H*")
    orig, = plt.plot(ori_degree_seq, color="black", linewidth=4, ls="-", label="H")

    plt.title("Degree Distribution")
    plt.ylabel("Degree")
    plt.ylabel("Ordered Vertices")

    plt.tick_params(
        axis="x",  # changes apply to the x-axis
        which="both",  # both major and minor ticks are affected
        bottom="off",  # ticks along the bottom edge are off
        top="off",  # ticks along the top edge are off
        labelbottom="off",
    )  # labels along the bottom edge are off

    plt.legend([orig, h], ["$H$", "HRG $H^*$"], loc=3)
    # fig = plt.gcf()
    # fig.set_size_inches(5, 4, forward=True)
    plt.show()

def plot_co_x(cox, start, end, size = (20,20), title = '', weighted=False, weight_threshold=10):

        """ Plotting function for keyword graphs

        Parameters
        --------------------
        cox: the coword networkx graph; assumes that nodes have attribute 'topic'
        start: start year
        end: end year
        """

        plt.figure(figsize=size)
        plt.title(title +' %s - %s'%(start,end), fontsize=18)
        if weighted:
            elarge=[(u,v) for (u,v,d) in cox.edges(data=True) if d['weight'] >weight_threshold]
            esmall=[(u,v) for (u,v,d) in cox.edges(data=True) if d['weight'] <=weight_threshold]
            pos=nx.graphviz_layout(cox) # positions for all nodes
            nx.draw_networkx_nodes(cox,pos,
                node_color= [s*4500 for s in nx.eigenvector_centrality(cox).values()],
                node_size = [s*6+20  for s in nx.degree(cox).values()],
                alpha=0.7)
            # edges
            nx.draw_networkx_edges(cox,pos,edgelist=elarge,
                                width=1, alpha=0.5, edge_color='black') #, edge_cmap=plt.cm.Blues
            nx.draw_networkx_edges(cox,pos,edgelist=esmall,
                                width=0.3,alpha=0.5,edge_color='yellow',style='dotted')
            # labels
            nx.draw_networkx_labels(cox,pos,font_size=10,font_family='sans-serif')
            plt.axis('off')
        else:
            nx.draw_graphviz(cox, with_labels=True,
                         alpha = 0.8, width=0.1,
                         fontsize=9,
                         node_color = [s*4 for s in nx.eigenvector_centrality(cox).values()],
                         node_size = [s*6+20 for s in nx.degree(cox).values()])

def compute_eps_multi(G, G2, k_arr):
    max_deg_G2 = max(nx.degree(G2).itervalues())
    print "max_deg_G2 =", max_deg_G2
    deg_list = nx.degree(G)     # G not G2
    deg_list_G2 = nx.degree(G2)     
    
    X = compute_X(G2, G2.number_of_nodes(), max_deg_G2, deg_list_G2)
    Y = compute_Y(X, G2.number_of_nodes(), max_deg_G2, deg_list_G2)
    # check X, Y
    check_X_and_Y(X, Y, G2.number_of_nodes(), max_deg_G2, deg_list_G2)
    
    ent = compute_entropy(Y, G2.number_of_nodes(), max_deg_G2, deg_list_G2)
    print "len(ent) =", len(ent)
#    print "entropy =", ent
    
    eps_arr = []
    for k in k_arr:
        num_violated = 0
        LOG2K = math.log(k,2)
#        print "LOG2K =", LOG2K
        for (v, deg) in deg_list.iteritems():   # check the original graph
            if deg <= max_deg_G2:               # in the case of max_deg_G2 < max_deg_G
                if ent[deg] > 0.0 and ent[deg] < LOG2K:   # do not check zero-column of ent
                    num_violated += 1
        # check and update eps_min (if ok)
        eps2 = float(num_violated)/G2.number_of_nodes()
    
        eps_arr.append(eps2)
    #
    return eps_arr

def nodeMaxDegree(G):
    degree=0
    for n in G.nodes():
        if nx.degree(G,n)> degree :
            degree=nx.degree(G,n)
            node=n
    return node

def print_degree_distributions(dataset, context):
    """
    Extracts degree distribution values from networks, and print them to
    cvs-file.

    **warning** overwrites if file exists.
    """
    print '> Reading data..', dataset
    corpus_path = '../data/'+dataset+'_text'
    (documents, labels) = data.read_files(corpus_path)

    degsfile = open('output/properties/cooccurrence/degrees_docs_'+dataset.replace('/','.'), 'w')

    giant = nx.DiGraph()
    print '> Building networks..'
    for i, text in enumerate(documents):
        if i%10==0: print '   ',str(i)+'/'+str(len(documents))
        g = graph_representation.construct_cooccurrence_network(text,context=context)
        giant.add_edges_from(g.edges())
        degs = nx.degree(g).values()
        degs = [str(d) for d in degs]
        degsfile.write(','.join(degs)+'\n')
    degsfile.close()

    print '> Writing giant\'s distribution'
    with open('output/properties/cooccurrence/degrees_giant_'+dataset.replace('/','.'), 'w') as f:
        ds = nx.degree(giant).values()
        ds = [str(d) for d in ds]
        f.write(','.join(ds))

def cand_gen(G, k, a, c, cores=None):
    '''
    generate all the branch given a k-plex a
    a: current k-plex
    c: the block
    '''
    b = neighbors(G, a)
    b.difference_update(c)

    # the strict nodes
    subg = G.subgraph(a)
    strict_nodes = {node for node in a if nx.degree(subg, node) == len(a)-k }
    for node in strict_nodes:
        b.intersection_update(G.neighbors(node))


    # always reshape by optimal
    if cores is None:
        b = {node for node in list(b) if nx.degree(G, node) >= len(optimal)-k}
    else:
        b = {node for node in list(b) if cores[node]>=len(optimal)-k}

    # calculate the valid candidates
    b = {node for node in b if len(set(G.neighbors(node)).intersection(a)) >= len(a)+1-k }

    # sort the candidate list
    b = list(b)
    # b.sort(key = lambda x: len(set(G.neighbors(x)).intersection(a)), reverse=True)

    return b

def compute_eps(G2, k, u, S1, S2):
    max_deg_G2 = max(nx.degree(G2).itervalues())
    print "max_deg_G2 =", max_deg_G2
    deg_list = nx.degree(G2)
    
#    start = time.clock()
    X = compute_X(G2, max_deg_G2)
#    print "compute_X: done"
#    print "Elapsed ", (time.clock() - start
    Y = compute_Y(X, G2.number_of_nodes(), max_deg_G2)
#    print "compute_Y: done"
    ent = compute_entropy(Y, G2.number_of_nodes(), max_deg_G2)
#    print "len(ent) =", len(ent)
    
    num_violated = 0
    LOG2K = math.log(k,2)
    print "LOG2K =", LOG2K
    for (v, deg) in deg_list.iteritems():   # check the original graph
        if (v == u or v in S1) and ent[deg] < LOG2K:
            num_violated += 1
    # check and update eps_min (if ok)
    eps2 = float(num_violated)/G2.number_of_nodes()
    
    #
    return eps2

def PartitionGraph(graph,kc_nodes,anchor=EMPTY_SET,k=2):
    G = graph.copy()
    R_cand = set()
    S_cand = set()
    
    G_ccs = nx.connected_component_subgraphs(G)
    
    for g_cc in G_ccs:
        cc_nodes = set(g_cc.nodes())
        kc_overlap = cc_nodes.intersection(kc_nodes)
        if len(kc_overlap) > 0:
            root = kc_overlap.pop()
            R_nodes = set()
            for n in cc_nodes:
                d = nx.degree(G,n) 
                if n not in anchor and n not in kc_nodes and d > 0 and d < k:
                    R_nodes.add(n)
                R_cand = R_cand.union(set((u,root) for u in R_nodes))                
        else:
            S_nodes = set()
            for n in cc_nodes:
                d = nx.degree(G,n) 
                if n not in anchor and d > 0 and d < k:
                    S_nodes.add(n)
            S_cand = S_cand.union(set((u,v) for u,v in combinations(S_nodes,2)))    

    return R_cand,S_cand

    def __generate_paths(self):
        self.__paths = []
        graph = copy.deepcopy(self.__graph)
        paths = []
        if not graph.nodes():
            return
        start_nodes = [n for n in graph.nodes() if (n.style() != "curve" or (nx.degree(graph, n) == 1))]
        start_nodes = sorted(start_nodes, key=lambda n: nx.degree(graph, n))
        if start_nodes:
            path = [start_nodes[0]]
        else:
            path = [graph.nodes()[0]]

        while path:
            neighbors = nx.neighbors(graph, path[-1])
            if neighbors:
                node = neighbors[0]
                graph.remove_edge(path[-1], node)
                path.append(node)
            else:
                paths.append(copy.copy(path))
                while path and not graph.neighbors(path[-1]):
                    path.pop()

        for path in paths:
            self.__paths.append(reduce_path(path))

        return

def leaf_removal(g, verbose=False):
    G = g.copy()
    stop = 0;
    potential_mis = [];
    isolated = [x for x in g.nodes() if nx.degree(g,x)==0];
    potential_mis.extend(isolated);
    G.remove_nodes_from(isolated);
    while stop==0:
        deg = G.degree();
        if 1 in deg.values():
            for n in G.nodes_iter():
                if deg[n]==1:
                    L = n;
                    break;
            nn = nx.neighbors(G,L)[0]
            G.remove_node(L);
            G.remove_node(nn);
            potential_mis.append(L);
            isolated = [x for x in G.nodes() if nx.degree(G,x)==0];
            potential_mis.extend(isolated);
            G.remove_nodes_from(isolated);
        else:
            stop=1;
    core_mis = [];
    if G.number_of_nodes()>=1:
        core_mis = nx.maximal_independent_set(G);
        if verbose==True:
            print len(potential_mis), len(core_mis), N;
        potential_mis.extend(core_mis);
    else:
        if verbose==True:
            print len(potential_mis), len(core_mis), N;
    return potential_mis, core_mis;

def compute_eps_deterministic_multi(G, G2, k_arr):
    max_deg_G2 = max(nx.degree(G2).itervalues())
    print "max_deg =", max_deg_G2
    deg_list = nx.degree(G) 
    deg_list_G2 = nx.degree(G2)  
    
    deg_count_G2 = [0 for j in range(max_deg_G2+1)]
    for deg in deg_list_G2.itervalues():
        deg_count_G2[deg] += 1
    
    ent = [0.0 for j in range(max_deg_G2+1)]
    for j in range(max_deg_G2+1):
        if deg_count_G2[j] > 0:
            ent[j] = math.log(deg_count_G2[j],2)           # simple
    
    print "len(ent) =", len(ent)
#    print "entropy =", ent
    
    eps_arr = []
    for k in k_arr:
        num_violated = 0
        LOG2K = math.log(k,2)
        for (v, deg) in deg_list.iteritems():   # check the original graph
            if deg <= max_deg_G2:               # in the case of max_deg_G2 < max_deg_G
                if ent[deg] > 0.0 and ent[deg] < LOG2K:   # do not check zero-column of ent
                    num_violated += 1
        # check and update eps_min (if ok)
        eps2 = float(num_violated)/G2.number_of_nodes()
    
        eps_arr.append(eps2)
    #
    return eps_arr

def RW_Size_col(G,r = 30000):
    sampled = []
    now_node = random.choice(G.nodes())
    sampled.append(now_node)
    sumA = 0.0
    sumB = 0.0
    sumA += nx.degree(G,now_node)
    sumB += 1.0/nx.degree(G,now_node)
    count = 0
    while True:
        next_node = random.choice(nx.neighbors(G,now_node))
        now_node = next_node
        sumA += nx.degree(G,now_node)
        sumB += 1.0/nx.degree(G,now_node)
        sampled.append(now_node)
        count += 1
        if count >= r:
            break
    count2 = 0
    for i in range(0,len(sampled)-1):
        for j in range(i+1,len(sampled)):
            if(sampled[i] == sampled[j]):
                count2 += 1

    return sumA*sumB/(2*count2)

def RW_Size(G,r = 1000,m=100):
    sampled = []
    now_node = random.choice(G.nodes())
    sampled.append(now_node)
    while True:
        next_node = random.choice(nx.neighbors(G,now_node))
        now_node = next_node
        sampled.append(now_node)
        if len(sampled) >= r:
            break
    print(1)
    lst = []
    for i in range(0,r-m):
        if i+m <= r-1:
            for j in range(i+m,r):
               # l1 = set(nx.neighbors(G,sampled[i]))
               # l2 = set(nx.neighbors(G,sampled[j]))
               # if len(list(l1 & l2)) >= 1:
                lst.append((sampled[i],sampled[j]))
                lst.append((sampled[j],sampled[i]))
    sumA = 0.0
    sumB = 0.0
    print(len(lst))
    for nodes in lst:
        sumA += float(nx.degree(G,nodes[0]))/nx.degree(G,nodes[1])
        l1 = set(nx.neighbors(G,nodes[0]))
        l2 = set(nx.neighbors(G,nodes[1]))
        count = len(list(l1&l2))
        sumB += count/(float(nx.degree(G,nodes[0]))*nx.degree(G,nodes[1]))
    return sumA/sumB

def _tester():
    # think of this like a networkx scratchpad
    G = nx.Graph()  # this is an undirected graph
    G.add_edge(1, 2)
    G.add_edge(2, 3)
    G.add_node(4)
    print nx.degree(G)
    print nx.info(G)

def mindeg_GSK(BG, variables_index=0, verbose=False):
    Vprime1 = [];
    Vprime2 = [];
    
    layer = nx.get_node_attributes(BG,'bipartite');
    var = [x for x in BG.nodes() if layer[x] == variables_index]
    fac = [x for x in BG.nodes() if layer[x] != variables_index]
    
    if verbose==True:
        print 'Initial variable nodes:', var;
        print 'Initial factor nodes:', fac;

    isolated_variables = [x for x in BG.nodes() if nx.degree(BG,x)==0 and layer[x]==variables_index];
    [var.remove(x) for x in isolated_variables]
    
    G = BG.copy();
    Vprime1.extend(isolated_variables);
    G.remove_nodes_from(isolated_variables)
    
    isolated_factors = [x for x in G.nodes() if nx.degree(BG,x)==0 and layer[x]!=variables_index];
    [fac.remove(x) for x in isolated_factors]
    G.remove_nodes_from(isolated_factors);

    while len(var)>0:
        if verbose==True:
            print '#var:',len(var),'#fac:', len(fac), '#nodes in depleted graph:', G.number_of_nodes(),'#original BG:',BG.number_of_nodes();

        pendant = return_mindeg_pendant(G,layer,variables_index);
        if len(pendant)==0:
            ## if not, choose randomly and do the game. 
            if verbose==True:
                print var
            m = G.number_of_nodes()*2;
            degs = G.degree();
            for e in G.edges():
                if degs[e[0]] + degs[e[1]] < m:
                    m = degs[e[0]] + degs[e[1]];
                    v = e;
            if e[0] in var:
                v = e[0];
            else:
                v = e[1];
            pendant = []
            pendant.append(v);
            pendant.extend(nx.neighbors(G,v));
            Vprime2.append(pendant[0]);
        else:
            Vprime1.append(pendant[0]);
        augmented_pendant = []
        augmented_pendant.extend(pendant);
        for n in pendant[1:]:
            augmented_pendant.extend(nx.neighbors(G,n));
        augmented_pendant = list(set(augmented_pendant));
        G.remove_nodes_from(augmented_pendant);        
        [var.remove(x) for x in augmented_pendant if x in var];
        [fac.remove(x) for x in augmented_pendant if x in fac];

    return Vprime1,Vprime2;

 def k_dependency_feats(self):
     wordindex = self.index + 1
     headindex = dep_head_of(self.deptree,wordindex)
     D = {}
     D["k_dist_to_root"] = len(dep_pathtoroot(self.deptree,wordindex))
     D["k_deprel"] = self.deptree[headindex][wordindex]["deprel"]
     D["k_headdist"] = abs(headindex - wordindex) # maybe do 0 for root?
     D["k_head_degree"] = nx.degree(self.deptree,headindex)
     D["k_child_degree"] = nx.degree(self.deptree,wordindex)
     return D

def sub_by_degree(min_degree,net):
    #return a subset of net by biger than min_degree
    remove_node_list = [item for item in nx.degree(net) \
            if nx.degree(net)[item] < min_degree]
    new_net = nx.Graph()
    new_net.add_nodes_from(net.nodes())
    new_net.add_edges_from(net.edges())
    for name in remove_node_list:
        new_net.remove_node(name)
    return new_net

def main():
    LOG = True

    #if (len(sys.argv) != 3):
     #       print "ERROR: genRandomGeorml <nodes> <raio>"
      #      sys.exit(1)

    NMAX = int(sys.argv[1])
    RAIO = float(sys.argv[2])
    #NMAX=40
    #RAIO=0.1
    ALCANCE=250

    G=nx.random_geometric_graph(NMAX,RAIO,2)

    while not nx.is_connected(G):
         RAIO=RAIO+.005
         G=nx.random_geometric_graph(NMAX,RAIO,2)
         if LOG: print "Graph is not full connected"

    pos=nx.get_node_attributes(G,'pos')
    network(G,pos,1)

    #Remove vizinhos que estejam demasiado perto
    while nodeNear(G)<1000 :
        G.remove_node(nodeNear(G))

    if nx.is_connected(G):
        pos=nx.get_node_attributes(G,'pos')
        network(G,pos,2)

        #Remove no que tem mais vizinhos
        T=G
        if not nodeSolo(T,nodeMaxDegree(T)): T.remove_node(nodeMaxDegree(T))
        if nx.is_connected(T):
                G=T

        pos=nx.get_node_attributes(G,'pos')
        network(G,pos,3)



        for n in G.neighbors(nodeMaxDegree(G)):
            if nx.degree(G,n)== 2 :
                degree=nx.degree(G,n)
                node=n
                print "node=",n
                if not nodeSolo(G,n): G.remove_node(n)
                break
        
        pos=nx.get_node_attributes(G,'pos')
        network(G,pos,4)
    else:
        if LOG: print "SubGraph is not full connected"

	def apply_sifi_surcharge(self):
		degree_sum = 0
		for bank in self.network.contracts:
			degree_sum += float(nx.degree(self.network.contracts)[bank])
		average_degree = float(degree_sum / len(self.network.contracts.nodes()))
		
		for bank in self.network.contracts:
			# the sifi surcharge is the product of the sifiSurchargeFactor and the connectedness as measured
			# by degree/average_degree 
			# the maximum ensures that no bank has to hold less than 1.0 times their banking capital
			sifiSurcharge = max(self.get_state(0).sifiSurchargeFactor*( float(nx.degree(self.network.contracts)[bank]) / average_degree), 1.0)
			bank.apply_sifi_surcharge(sifiSurcharge)

def compute_measures(bigDict):
    """ Computes the measures for each network
    
    Measures to compute:
    
    nr_of_nodes
    nr_of_edges
    
    max_edge_value
    min_edge_value
    
    is_connected
    number_connected_components
    
    average_unweighted_node_degree
    average_weighted_node_degree
    
    average_clustering_coefficient
    average_weighted_shortest_path_length
    average_unweighted_shortest_path_length
    
    To be added:
     single node values, e.g. node degree of brainstem etc.
    
    Non-scalar return values: (not used yet)
     degree_distribution
     edge_weight_distribution
    
    """
    
    returnMeasures = {}

    for key, netw in bigDict.items():

        outm = {}
        
        outm['nr_of_nodes'] = netw.number_of_nodes()
        outm['nr_of_edges'] = netw.number_of_edges()
        
        outm['max_edge_value'] = np.max([d['weight']for f,t,d in netw.edges(data=True)])
        outm['min_edge_value'] = np.min([d['weight']for f,t,d in netw.edges(data=True)])
        
        outm['is_connected'] = nx.is_connected(netw)
        outm['number_connected_components'] = nx.number_connected_components(netw)
        outm['average_unweighted_node_degree'] =  np.mean(nx.degree(netw, weighted = False).values())
        outm['average_weighted_node_degree'] = np.mean(nx.degree(netw, weighted = True).values())
        outm['average_clustering_coefficient'] = nx.average_clustering(netw)
        outm['average_weighted_shortest_path_length'] = nx.average_shortest_path_length(netw, weighted = True)
        outm['average_unweighted_shortest_path_length'] = nx.average_shortest_path_length(netw, weighted = False)
        
        returnMeasures[key] = outm
        
    return returnMeasures