diff --git a/.gitignore b/.gitignore index e2c98f8..f88b391 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ docs/build/ docs/site/ .vscode/ Manifest.toml + diff --git a/Project.toml b/Project.toml index 57979f5..db9408e 100644 --- a/Project.toml +++ b/Project.toml @@ -8,9 +8,11 @@ Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" EmbeddedGraphs = "9c1af47c-29f5-11e9-0b47-9f5334384f20" LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d" +MetaGraphs = "626554b9-1ddb-594c-aa3c-2596fe9399a5" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SpatialIndexing = "d4ead438-fe20-5cc5-a293-4fd39a41b74c" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] diff --git a/src/NodeTypes.jl b/src/NodeTypes.jl new file mode 100644 index 0000000..6e90b25 --- /dev/null +++ b/src/NodeTypes.jl @@ -0,0 +1,21 @@ +"""For the incorporation of different node types to SyntheticNetworks package.""" + +struct NodeType + name::String + probability::Float32 + neighbor_probability::Dict # Dict(Symbol(NodeType_i.name) => prob) + method_parameter +end + +NodeType() = NodeType("Node", 1., Dict(:Node => 1.), 0) + +""" From a list of probabilities of drawing a node type, determines one randomly and + returns the index of node type""" +function draw_type(node_types::Array{NodeType,1}) + prob = map(x -> x.probability, node_types) + t0 = copy(prob) + pushfirst!(t0 , 0.) + type_interval = [sum(t0[1:i+1]) for i in 1:length(t0)-1] + t_n = findfirst(type_interval .>= rand()) + return node_types[t_n] +end diff --git a/src/SyntheticNetworks.jl b/src/SyntheticNetworks.jl index 02b5168..28994ba 100755 --- a/src/SyntheticNetworks.jl +++ b/src/SyntheticNetworks.jl @@ -1,10 +1,15 @@ +__precompile__(false) + module SyntheticNetworks +import StatsBase: countmap, sample using Random using Parameters using LightGraphs using EmbeddedGraphs using Distances +using MetaGraphs +include("NodeTypes.jl") function __init__() @warn "The parameter `u` is currently not implemented." @@ -19,7 +24,9 @@ end f(i,j,G) = (d_G(i,j) + 1) ^ r / (dist_spatial(x_i,x_j)) """ abstract type SyntheticNetwork end -export SyntheticNetwork, RandomPowerGrid, initialise, generate_graph, grow! +export SyntheticNetwork, NodeType, RandomPowerGrid, initialise, generate_graph, grow! + +max_iter = 200 """ RandomPowerGrid(num_layers, n, n0, p, q, r, s, u, sampling, α, β, γ, debug) @@ -63,7 +70,6 @@ export SyntheticNetwork, RandomPowerGrid, initialise, generate_graph, grow! # export RandomPowerGrid # - struct RandomPowerGrid n::Int n0::Int @@ -72,60 +78,71 @@ struct RandomPowerGrid r::Float32 s::Float32 u::Float32 + types::Array{NodeType,1} + heuristic::Symbol function RandomPowerGrid(n, n0, p, q, r, s, u) new(n, n0, p, q, r, s, u) end end -RandomPowerGrid(n, n0) = RandomPowerGrid(n, n0, rand(5)...) +RandomPowerGrid(n, n0) = RandomPowerGrid(n, n0, rand(5)...,[NodeType()], :standard) +RandomPowerGrid(n,n0,p,q,r,s,u) = RandomPowerGrid(n,n0,p,q,r,s,u,[NodeType()], :standard) function generate_graph(RPG) + @assert RPG.heuristic in [:standard, :maxdegree, :neighbor_probability, :degree_cdf, :neighbor_probability_simple] # (n, n0, p, q, r, s, u) = (RPG.n, RPG.n0, RPG.p, RPG.q, RPG.r, RPG.s, RPG.s) - eg = initialise(RPG.n0, RPG.p, RPG.q, RPG.r, RPG.s, RPG.u) - grow!(eg, RPG.n, RPG.n0, RPG.p, RPG.q, RPG.r, RPG.s, RPG.u) - eg + eg, t_list = initialise(RPG.n0, RPG.p, RPG.q, RPG.r, RPG.s, RPG.u, + RPG.types, RPG.heuristic) + grow!(eg, t_list, RPG.n, RPG.n0, RPG.p, RPG.q, RPG.r, RPG.s, RPG.u, + RPG.types, RPG.heuristic) + mg = Embedded_to_MetaGraph(eg, t_list) + set_props!(mg, Dict(:pqrs => [RPG.p, RPG.q, RPG.r, RPG.s], :heuristic => RPG.heuristic)) + mg + return mg end -function initialise( - n0::Int, - p::Real, - q::Real, - r::Real, - s::Real, - u::Real; - vertex_density_prob::Function = rand_uniform_2D, -) - # we can skip this when only a single node is given - if n0 == 1 - positions = [vertex_density_prob(n0), ] - graph = EmbeddedGraph(complete_graph(n0), positions) - else - # STEP I1 - """If the locations x_1...x_N are not given, draw them independently at - random from ρ.""" - positions = [vertex_density_prob(i) for i = 1:n0] - graph = EmbeddedGraph(SimpleGraph(n0), positions) - - # STEP I2 - """Initialize G to be a minimum spanning tree (MST) for x_1...x_N w.r.t. - the distance function dist_spatial(x, y) (using Kruskal’s simple or - Prim’s more efficient algorithm). """ - mst_graph = EmbeddedGraph(complete_graph(n0), positions) - edges = prim_mst(mst_graph.graph, weights(mst_graph, dense = true)) - for edge in edges - add_edge!(graph, edge) - end - # STEP I3 - """With probability q, draw a node i ∈ {1,...,N} uniformly at - random, find that node l ∈ {1,...,N} which is not yet linked to i and - for which f (i,l,G) is maximal, and add the link i–l to G.""" - m = Int(round(n0 * (1 - s) * (p + q), RoundDown)) - for dummy = 1:m - i = rand(1:nv(graph)) - dist_spatial = - map(j -> euclidean(graph.vertexpos[i], graph.vertexpos[j]), 1:nv(graph)) - l_edge = Step_G34(graph, i, dist_spatial, r) +function Embedded_to_MetaGraph(eg::EmbeddedGraph, t_arr) + mg = MetaGraph(eg.graph,1.) + for (i,j) in enumerate(eg.vertexpos) + set_props!(mg,i,Dict(:pos => j, :type => t_arr[i].name)) + end + return mg +end + + +function initialise(n0::Int, p::Real, q::Real, r::Real, s::Real, u::Real, + node_types::Array{NodeType,1}, heuristic::Symbol; + vertex_density_prob::Function=rand_uniform_2D) + # STEP I1 + """If the locations x_1...x_N are not given, draw them independently at + random from ρ.""" + positions = [vertex_density_prob(i) for i=1:n0 ] + types = [draw_type(node_types) for i = 1:n0] + graph = EmbeddedGraph(SimpleGraph(n0), positions) + + # STEP I2 + """Initialize G to be a minimum spanning tree (MST) for x_1...x_N w.r.t. + the distance function dist_spatial(x, y) (using Kruskal’s simple or + Prim’s more efficient algorithm). """ + mst_graph = EmbeddedGraph(complete_graph(n0), positions) + edges = prim_mst(mst_graph.graph, weights(mst_graph, dense=true)) + for edge in edges + add_edge!(graph, edge) + end + + # STEP I3 + """With probability q, draw a node i ∈ {1,...,N} uniformly at + random, find that node l ∈ {1,...,N} which is not yet linked to i and + for which f (i,l,G) is maximal, and add the link i–l to G.""" + m = Int(round(n0*(1-s)*(p+q), RoundDown)) + for dummy in 1:m + i = rand(1:nv(graph)) + dist_spatial = map(j -> euclidean(graph.vertexpos[i], + graph.vertexpos[j]), 1:nv(graph)) + # + l_edge = Step_G34(graph, i, dist_spatial, r, types, heuristic) + if l_edge != 0 add_edge!(graph, l_edge, i) end end @@ -136,25 +153,19 @@ function initialise( # for which f(i,j,G) [eqn. 1] is maximal, and add the link i–j to G.""" # Step_I3!(graph, r, m) # OLD LOGIC - graph + return graph, types end -function grow!( - graph::EmbeddedGraph, - n::Int, - n0::Int, - p, - q, - r, - s, - u; - vertex_density_prob::Function = rand_uniform_2D, -) - for n_actual = n0+1:n +function grow!(graph::EmbeddedGraph, types, n::Int, n0::Int, p, q, r, s, u, + node_types::Array{NodeType}, heuristic::Symbol; + vertex_density_prob::Function=rand_uniform_2D) + for n_actual in n0+1:n # STEP G0 + t = draw_type(node_types) + push!(types,t) """With probabilities 1−s and s, perform either steps G1–G4 or step G5, respectively.""" - if (rand() >= s) | isempty(edges(graph)) + if rand() >= s || ne(graph) == 0 # STEP G1 """If x i is not given, draw it at random from ρ.""" pos = vertex_density_prob(n_actual) @@ -163,20 +174,48 @@ function grow!( # STEP G2 """ Find that node j ∈ {1,...,N} for which dist_spatial(x_i,x_j) is minimal and add the link i–j to G.""" - dist_spatial = map( - i -> euclidean(graph.vertexpos[nv(graph)], graph.vertexpos[i]), - 1:nv(graph), - ) - dist_spatial[nv(graph)] = 100000.0 #Inf - min_dist_vertex = argmin(dist_spatial) + dist_spatial = map(i -> euclidean(graph.vertexpos[nv(graph)], graph.vertexpos[i]), 1:nv(graph)) + dist_spatial[nv(graph)] = 100000. #Inf + if heuristic == :standard + min_dist_vertex = argmin(dist_spatial) + elseif heuristic == :neighbor_probability_simple + min_dist_vertex = argmin(dist_spatial) + found = types[min_dist_vertex] !== types[nv(graph)] + iter1 = 0 + while !found && iter1 <= max_iter + if rand() < 1/2 + dist_spatial[min_dist_vertex] = 100000. + min_dist_vertex = argmin(dist_spatial) + found = types[min_dist_vertex] !== types[nv(graph)] + else + found = true + end + iter1 += 1 + end + elseif heuristic == :maxdegree + d = degree(graph) + candidates = [d[i] < types[i].method_parameter for i in 1:nv(graph)] + dist_spatial[.!candidates] .= 100000. + min_dist_vertex = argmin(dist_spatial) + elseif heuristic == :degree_cdf + d = degree(graph) + candidates = [rand() <= probs(d[j],types[j].method_parameter) for j in 1:nv(graph)] + dist_spatial[.!candidates] .= 100000. + min_dist_vertex = argmin(dist_spatial) + end add_edge!(graph, min_dist_vertex, nv(graph)) + + # STEP G3 """ With probability p, find that node l ∈ {1,...,N} ⍀ {j} for which f(i,l,G) is maximal, and add the link i–l to G.""" if rand() <= p - l_edge = Step_G34(graph, nv(graph), dist_spatial, r) - add_edge!(graph, l_edge, nv(graph)) + l_edge = Step_G34(graph, nv(graph), dist_spatial, r, types, heuristic) + if l_edge !== 0 + add_edge!(graph, l_edge, nv(graph)) + end + end # STEP G4 """ With probability q, draw a node i' ∈ {1,...,N} uniformly at @@ -184,11 +223,28 @@ function grow!( i' and for which f(i',l',G) is maximal, and add the link i'–l' to G.""" if rand() <= q - i = rand(1:nv(graph)) - dist_spatial = - map(j -> euclidean(graph.vertexpos[i], graph.vertexpos[j]), 1:nv(graph)) - l_edge = Step_G34(graph, i, dist_spatial, r) - add_edge!(graph, l_edge, i) + if heuristic == :degree_cdf + d = degree(graph) + candidates = [rand() <= probs(d[j],types[j].method_parameter) for j in 1:nv(graph)] + iter2 = 0 + while isempty(findall(candidates)) && iter2 <= max_iter + candidates = [rand() <= probs(d[j],types[j].method_parameter) for j in 1:nv(graph)] + iter2 += 1 + end + i = rand(findall(candidates)) + elseif heuristic == :maxdegree + d = degree(graph) + candidates = [d[i] < types[i].method_parameter for i in 1:nv(graph)] + i = rand(findall(candidates)) + else + i = rand(1:nv(graph)) + end + dist_spatial = map(j -> euclidean(graph.vertexpos[i], graph.vertexpos[j]), 1:nv(graph)) + l_edge = Step_G34(graph, i, dist_spatial, r, types, heuristic) + if l_edge !== 0 + add_edge!(graph, l_edge, i) + end + end else @@ -199,10 +255,7 @@ function grow!( New logic splits the nodes somewhere, not in the middle.""" edge = rand(edges(graph)) splitval = rand() - newpos = ( - graph.vertexpos[src(edge)] * splitval + - graph.vertexpos[dst(edge)] * (1 - splitval) - ) + newpos = (graph.vertexpos[src(edge)] * splitval + graph.vertexpos[dst(edge)] * (1-splitval)) add_vertex!(graph, newpos) add_edge!(graph, src(edge), nv(graph)) add_edge!(graph, dst(edge), nv(graph)) @@ -213,26 +266,117 @@ end function Step_I3!(g::EmbeddedGraph, r::Real, m::Int) - for dummy = 1:m - spatial = weights(g, dense = true) + for dummy in 1:m + spatial = weights(g, dense=true) A = floyd_warshall_shortest_paths(g, weights(g)).dists A = ((A .+ spatial) .^ r) ./ spatial - map(i -> A[i, i] = 0, 1:size(A)[1]) - add_edge!(g, Tuple(argmax(A))...) + map(i -> A[i,i] = 0, 1:size(A)[1]) + add_edge!(g,Tuple(argmax(A))...) + end +end + +function Step_G34(g::EmbeddedGraph, i::Int, dist_spatial, r, types::Array{NodeType}, heuristic::Symbol) + if heuristic == :standard + return Step_G34_standard(g, i, dist_spatial, r, types) + elseif heuristic == :maxdegree + return Step_G34_maxdegree(g, i, dist_spatial, r, types) + elseif heuristic == :neighbor_probability + return Step_G34_neighbours(g, i, dist_spatial, r, types) + elseif heuristic == :degree_cdf + return Step_G34_degree_cdf(g, i, dist_spatial, r, types) + elseif heuristic == :neighbor_probability_simple + return Step_G34_neighbours_simple(g, i, dist_spatial, r, types) + end +end + + +function Step_G34_standard(g::EmbeddedGraph, i::Int, dist_spatial, r, types::Array{NodeType}) + V = dijkstra_shortest_paths(g, i).dists + V = ((V .+ dist_spatial) .^ r) ./ dist_spatial + V[i] = 0 + V[neighbors(g, i)] .= 0 + m = argmax(V) + return maximum(V) > 0. ? m : 0 +end + + +function Step_G34_maxdegree(g::EmbeddedGraph, i::Int, dist_spatial, r, types::Array{NodeType}) + d = degree_centrality(g, normalize = false) + candidates = [d[i] < types[i].method_parameter for i in 1:nv(g)] + if true in candidates + V = dijkstra_shortest_paths(g, i).dists + V = ((V .+ dist_spatial) .^ r) ./ dist_spatial + V[i] = 0 + V[neighbors(g, i)] .= 0 + V[.!candidates] .= 0 + m = argmax(V) + return maximum(V) > 0 ? m : 0 + else + return 0 end end +function Step_G34_neighbours(g::EmbeddedGraph, i::Int, dist_spatial, r, types::Array{NodeType}) + n_prob = map(x -> x.neighbor_probability[Symbol(types[i].name)], types) + candidates = n_prob .> 0. + if true in candidates + V = dijkstra_shortest_paths(g, i).dists + V = ((V .+ dist_spatial) .^ r) ./ dist_spatial + V[i] = 0 + V[neighbors(g, i)] .= 0 + V[.!candidates] .= 0 + m = argmax(V) + return maximum(V) > 0 ? m : 0 + else + return 0 + end +end -function Step_G34(g::EmbeddedGraph, i::Int, dist_spatial, r) +function Step_G34_neighbours_simple(g::EmbeddedGraph, i::Int, dist_spatial, r, types::Array{NodeType}) V = dijkstra_shortest_paths(g, i).dists V = ((V .+ dist_spatial) .^ r) ./ dist_spatial V[i] = 0 V[neighbors(g, i)] .= 0 - argmax(V) + if maximum(V) > 0 + j = argmax(V) + same = types[i] == types[j] + iter3 = 0 + while same && iter3 <= max_iter + if rand() < 1/2 + V[j] = 0 + j = argmax(V) + same = types[i] == types[j] + else + same = false + end + iter3 += 1 + end + return j + else + return 0 + end +end + + +probs(k,p) = exp(p[1] + p[2] * k) +function Step_G34_degree_cdf(g::EmbeddedGraph, i::Int, dist_spatial, r, types::Array{NodeType}) + d = degree(g) + candidates = [rand() <= probs(d[j],types[j].method_parameter) for j in 1:nv(g)] + if true in candidates + V = dijkstra_shortest_paths(g, i).dists + V = ((V .+ dist_spatial) .^ r) ./ dist_spatial + V[i] = 0 + V[neighbors(g, i)] .= 0 + V[.!candidates] .= 0. + m = argmax(V) + return maximum(V) > 0 ? m : 0 + else + return 0 + end end -rand_uniform_2D(i) = [rand_uniform(i), rand_uniform(i)] -rand_uniform(i) = 2 * (0.5 - rand()) +rand_uniform_2D(i) = [rand_uniform(i),rand_uniform(i)] +rand_uniform(i) = 2*(0.5-rand()) end # module diff --git a/test/runtests.jl b/test/runtests.jl index ffeb7e1..a8fc8eb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,13 +3,20 @@ using EmbeddedGraphs: EmbeddedGraph, nv using Random Random.seed!(42); -const RPG = RandomPowerGrid(100, 1) -const g = generate_graph(RPG) +const p = 0.3 +const q = 0.2 +const r = 1//3 +const s = 0.1 +const u = 0. + +RPG = RandomPowerGrid(100, 1, p, q, r, s, u, " ", [1], [(g::EmbeddedGraph, i::Int) -> true]) +eg, t_list = initialise(RPG.n0, RPG.p, RPG.q, RPG.r, RPG.s, RPG.u, RPG.t_name, RPG.t_prob, RPG.t_method) +grow!(eg, t_list, RPG.n, RPG.n0, RPG.p, RPG.q, RPG.r, RPG.s, RPG.u, RPG.t_name, RPG.t_prob, RPG.t_method) @testset "RPG" begin @test RPG isa RandomPowerGrid - @test g isa EmbeddedGraph - @test nv(g) == RPG.n + @test eg isa EmbeddedGraph + @test t_list isa Array{Int} + @test nv(eg) == RPG.n end -