def create_sp_graph(n=30, l=0.001, spmax = 0.4):

    #first create a regular ring
    G = nx.watts_strogatz_graph(n, 6, 0, seed=None)

    #then rewire all edges
    nx.double_edge_swap(G, len(G.edges)/2, 100000)
    sp_t = sp_of_graph(G)

    for i in range(20000000):
        orig_G = G.copy()
        sp_orig = sp_t

        if i%100 == 0:
            print('goal:',spmax,' [',i/200000,'%], sp val: ',sp_orig)


        if i%2000 == 0:
            nx.write_adjlist(G, 'z1000_d6_sp'+str(sp_t)+'_iter'+str(i)+'.adjlist')


        edge = random.choice(list(G.edges))

        G.remove_edge(edge[0],edge[1])

        a = random.randint(0,n)
        b = random.randint(0,n)
        while ((a,b) in G.edges()):
            a = random.randint(0,n)
            b = random.randint(0,n)
        G.add_edge(a,b)

        if not nx.has_path(G,edge[0],edge[1]):
            G = orig_G.copy()
            continue

        sp_t = sp_of_graph(G)
        prob = l*(spmax - sp_t)

        if (not(sp_t > sp_orig) or (random.random() < prob)):
            G = orig_G.copy()
            sp_t = sp_orig
            continue

        if (sp_t > spmax):
            print(sp_t)
            return G
    print(sp_t)
    return G

def sp_of_graph(G):
    return np.mean([sp_of_node(G, x) for x in G.nodes])

def sp_of_node(G,n):
    RA = set()
    for node in G.neighbors(n):
        RA.add(node)
        RA.update(G.neighbors(node))
    return np.mean([sp_3(G, n, x) for x in RA])

def kdelta(G,a,b):
    return int((a,b) in G.edges)

def sp_3(G, A, B):
    if A == B:
        return 1
    bneigh = set(nx.neighbors(G,B))
    union = bneigh & set(nx.neighbors(G,A))
    res = 2 * kdelta(G,A,B) + len(union)
    return np.divide(res,len(bneigh) + 1  )