diff --git a/REQUIRE b/REQUIRE index 5b3c60e..aec8fe9 100644 --- a/REQUIRE +++ b/REQUIRE @@ -6,4 +6,6 @@ Polynomials SimplePartitions Graphs IterTools -Iterators +Iterators +Optim +LightXML diff --git a/src/SimpleGraphs.jl b/src/SimpleGraphs.jl index 2fa3e6b..b1ee571 100644 --- a/src/SimpleGraphs.jl +++ b/src/SimpleGraphs.jl @@ -42,4 +42,6 @@ include("d_simple_matrices.jl") include("d_dist.jl") +include("embedding/GraphEmbeddings.jl") + end # module SimpleGraphs diff --git a/src/embedding/GraphEmbeddings.jl b/src/embedding/GraphEmbeddings.jl new file mode 100644 index 0000000..e5b5bd9 --- /dev/null +++ b/src/embedding/GraphEmbeddings.jl @@ -0,0 +1,616 @@ + +using LinearAlgebra, Statistics + +import Base.show, LinearAlgebra.scale! + +export embed, remove_embedding, has_embedding +export set_fill_color, set_line_color, get_fill_color, get_line_color + + +export GraphEmbedding, show +# export draw, draw_labels, unclip + +export transform, scale, rotate, translate, recenter +export getxy +export set_vertex_size, get_vertex_size +export edge_length + +const DEFAULT_MARKER_SIZE = 10 +MARKER_SIZE = DEFAULT_MARKER_SIZE + + +""" +`GraphEmbedding(G)` creates a drawing of the `SimpleGraph` `G`. It +is given a circular embedding. + +**Note**: Direct usage of `GraphEmbedding` is **deprecated**! + +`GraphEmbedding(G,d)` creates a drawing using the `Dict` `d` to +specify the vertex locations. `d` should map vertices of `G` to +`Vector`s that specify x,y-coordinates. +""" +mutable struct GraphEmbedding + G::SimpleGraph + xy::Dict{Any,Vector{Float64}} + line_color::String + fill_color::String + vertex_size::Int + + + function GraphEmbedding(GG::SimpleGraph) + d = circular_embedding(GG) + new(GG,d,"black","white",DEFAULT_MARKER_SIZE) + end + + function GraphEmbedding(GG::SimpleGraph,d::Dict) + T = vertex_type(GG) + dd = Dict{T,Vector{Float64}}() + + try + for v in GG.V + dd[v] = d[v] + end + catch + error("Dictionary doesn't include all vertices") + end + new(GG,dd,"black","white", DEFAULT_MARKER_SIZE) + end +end + + +""" +`remove_embedding(G)` erases the graph's embedding saved in its cache. +""" +function remove_embedding(G::SimpleGraph) + cache_clear(G,:GraphEmbedding) + nothing +end + +""" +`has_embedding(G)` returns `true` if an embedding has been created +for `G` in its cache. +""" +function has_embedding(G::SimpleGraph) + return cache_check(G,:GraphEmbedding) +end + +""" +`get_embedding_direct(G)` returns a direct reference to a +graph's embedding. This function is not exposed. Users should +use `cache_recall(G,:GraphEmbedding)` instead to retrieve +a copy of the drawing. If there is not drawing, a default +(circular) drawing is created. +""" +function get_embedding_direct(G::SimpleGraph) + if !cache_check(G,:GraphEmbedding) + X = GraphEmbedding(G) + G.cache[:GraphEmbedding] = X + end + return G.cache[:GraphEmbedding] +end + +""" +`set_embedding_direct(G,X)` overwrites (or creates) the drawing +`X` in the graph's cache. This non-exposed function should not be +called by the user. +""" +function set_embedding_direct(G::SimpleGraph, X::GraphEmbedding) + G.cache[:GraphEmbedding] = X + nothing +end + +""" +`set_fill_color(G,color)` sets the color that gets drawn in the +interior of the vertices. (All vertices get the same color.) +""" +function set_fill_color(G::SimpleGraph, color::String="white") + if !has_embedding(G) + embed(G) + end + G.cache[:GraphEmbedding].fill_color=color + nothing +end + +""" +`get_fill_color(G)` returns the color used to fill in vertices. +See `set_fill_color`. +""" +function get_fill_color(G::SimpleGraph) + if !has_embedding(G) + embed(G) + end + return G.cache[:GraphEmbedding].fill_color +end + +""" +`set_line_color(G,color)` sets the color used to draw edges and +the circles around vertices. +""" +function set_line_color(G::SimpleGraph, color::String="black") + if !has_embedding(G) + embed(G) + end + G.cache[:GraphEmbedding].line_color = color + nothing +end + +""" +`get_line_color(G)` returns the color used to draw edges and the +circles around vertices. See `set_line_color`. +""" +function get_line_color(G::SimpleGraph) + if !has_embedding(G) + embed(G) + end + return G.cache[:GraphEmbedding].line_color +end + + + + +""" +`list2dict(list)` takes a list of `(Symbol,Any)` pairs and +converts them to a dictionary mapping the symbols to their +associated values. +""" +function list2dict(list::Vector) + d = Dict{Symbol,Any}() + for it in list + x,y = it + d[x]=y + end + return d +end + + + +""" +`embed(G)` creates a new embedding for `G`. The full call is +``` +embed(G,algorithm;args...) +``` +The `symbol` algorithm indicates the embedding algorithm. +The `args` are a collection of possible arguments to be sent +to the algorithm. + +The `algorithm` defaults to `:circular` and may be one of the following: + +* `:circular`: arrange the vertices evenly in a circle. +* `:random`: arrange the vertices randomly. +* `:spring`: use the `spring` layout from `GraphLayout`. Optional argument: + * `iterations`. +* `:stress`: use the `stress` layout from `GraphLayout`. +* `:distxy`: is my inefficient version of `stress`. Optional arguments: + * `tolerance` -- small positive real number indicating when to stop iterating. + * `verbose`-- a boolean specifying if verbose output is produced. +* `:spectral`: use the `spectral` embedding. Optional arguments: + * `xcol` -- which eigenvalue to use for the `x` coordinate. + * `ycol` -- which eigenvalue to use for the `y` coordinate + +Note that if the graph already has (say) an embedding, that embedding may +be used as the starting point for one of the algorithms. +""" +function embed(G::SimpleGraph, algorithm::Symbol=:circular; args...) + X = get_embedding_direct(G) + arg_dict = list2dict(collect(args)) + + ## FULL LIST OF POSSIBLE ARGUMENTS PRE-PARSED HERE + + verbose::Bool = true + if haskey(arg_dict,:verbose) + verbose = arg_dict[:verbose] + end + + iterations::Int = 0 + if haskey(arg_dict,:iterations) && arg_dict[:iterations]>0 + iterations = arg_dict[:iterations] + end + + tolerance = 0.001 + if haskey(arg_dict,:tolerance) && arg_dict[:tolerance]>0 + tolerance = arg_dict[:tolerance] + end + + xcol::Int = 2 + if haskey(arg_dict,:xcol) && arg_dict[:xcol] > 0 + xcol = arg_dict[:xcol] + end + + ycol::Int = 3 + if haskey(arg_dict,:ycol) && arg_dict[:ycol] > 0 + xcol = arg_dict[:ycol] + end + + if algorithm == :circular + circular!(X) + return nothing + end + + if algorithm == :random + random!(X) + return nothing + end + + if algorithm == :spring + if iterations <= 0 + spring!(X) + else + spring!(X,iterations) + end + return nothing + end + + if algorithm == :stress + stress!(X) + return nothing + end + + if algorithm == :distxy + distxy!(X,tolerance,verbose) + return nothing + end + + if algorithm == :combined + embed(G,:spring,iterations = iterations) + scale(G) + embed(G,:stress) + return nothing + end + + + if algorithm == :spectral + spectral!(X,xcol, ycol) + return nothing + end + + @warn "Unknown embedding algorithm; no action taken" + return nothing +end + + +""" +`embed(G,d)` specifies an embedding of the graph `G` with +a dictionary `d` that maps vertices to coordinates (as two +dimensional vectors `[x,y]`). +""" +function embed(G::SimpleGraph,d::Dict) + X = GraphEmbedding(G,d) + set_embedding_direct(G,X) + nothing +end + + + + + +function circular_embedding(G::SimpleGraph{T}) where T + d = Dict{T, Vector{Float64}}() + n = NV(G) + + vv = vlist(G) # vertices as a list + + s = 2*sin(pi/n) + + t = collect(0:n-1)*2*pi/n + x = map(sin,t)/s + y = map(cos,t)/s + + for i = 1:n + v = vv[i] + d[v] = [x[i],y[i]] + end + return d +end + + +function show(io::IO, X::GraphEmbedding) + print(io,"Embedding of $(X.G)") +end + + +""" +`set_vertex_size(G,sz)` sets the size of the circle used +when drawing vertices. The default is `$DEFAULT_MARKER_SIZE`. +See also `get_vertex_size`. +""" +function set_vertex_size(G::SimpleGraph, m::Int = DEFAULT_MARKER_SIZE) + if m < 1 + error("Vertex size must be nonnegative") + end + if !has_embedding(G) + embed(G) + end + G.cache[:GraphEmbedding].vertex_size = m + nothing +end + +""" +`get_vertex_size(G)` returns the size of the circle used when +drawing vertices. See also `set_vertex_size`. +""" +function get_vertex_size(G::SimpleGraph) + if !has_embedding(G) + embed(G) + end + return G.cache[:GraphEmbedding].vertex_size +end + +""" +`circular!(X)` arranges the vertices of the graph held in the +`GraphEmbedding` around a circle. +""" +function circular!(X::GraphEmbedding) + X.xy = circular_embedding(X.G) + scale!(X) + return +end + +function private_adj(X::GraphEmbedding) + n = NV(X.G) + A = zeros(Int,n,n) + vv = vlist(X.G) + for i=1:n-1 + u = vv[i] + for j=(i+1):n + w = vv[j] + if has(X.G,u,w) + A[i,j] = 1 + A[j,i] = 1 + end + end + end + return A,vv +end + + +function private_dist(X::GraphEmbedding) + n = NV(X.G) + d = dist(X.G) + + A = zeros(n,n) + vv = vlist(X.G) + + for i=1:n + u = vv[i] + for j=1:n + w = vv[j] + A[i,j] = d[u,w] + if A[i,j] < 0 + A[i,j] = n/2 # should be enough to separate the comps + end + end + end + + + return A,vv +end + +# TEMPORARILY OFF LINE DURING TRANSITION TO JULIA 0.7 + +include("my_spring.jl") + +""" +`spring!(X)` gives the graph held in `X` with a spring embedding +(based on code in the `GraphLayout` module). If runs a default number of +iterations (100) of that algorithm. To change the number of +iterations, use `spring!(X,nits)`. +""" +function spring!(X::GraphEmbedding, nits::Int=100) + n = NV(X.G) + A,vv = private_adj(X) + + x,y = layout_spring_adj(A,MAXITER=nits) + + d = Dict{vertex_type(X.G), Vector{Float64}}() + for i = 1:n + v = vv[i] + d[v] = [x[i], y[i]] + end + X.xy = d + return +end + +include("my_stress.jl") + +""" +`stress!(X)` computes a stress major layout using code taken from the +`GraphLayout` package. +""" +function stress!(X::GraphEmbedding) + n = NV(X.G) + A,vv = private_dist(X) + + currentxy = zeros(n,2) + for i=1:n + v = vv[i] + currentxy[i,:] = X.xy[v] + end + + xy = my_layout_stressmajorize_adj(A,2,nothing,currentxy) + + d = Dict{vertex_type(X.G), Vector{Float64}}() + for i = 1:n + v = vv[i] + d[v] = collect(xy[i,:]) + end + X.xy = d + return +end + + +""" +`random!(X)` gives the graph held in `X` a random embedding. +""" +function random!(X::GraphEmbedding) + rootn = sqrt(NV(X.G)) + for v in X.G.V + X.xy[v] = [rand(), rand()]*rootn + end + recenter!(X) +end + + +function transform!(X::GraphEmbedding, + A::Array{S,2}, + b::Vector{T}) where {S<:Real,T<:Real} + + # apply the transformation x |--> Ax+b to all coordinates + for k in keys(X.xy) + X.xy[k] = A*X.xy[k] + b + end + nothing +end + +""" +`transform(G,A,b)` applies an affine transformation to all coordinates +in the graph's drawing. Here `A` is 2-by-2 matrix and `b` is a 2-vector. +Each point `p` is mapped to `A*p+b`. +""" +function transform(G::SimpleGraph, + A::Array{S,2}, + b::Vector{T}) where {S<:Real, T<:Real} + X = get_embedding_direct(G) + transform!(X,A,b) +end + + +function LinearAlgebra.scale!(X::GraphEmbedding, m::T) where T<:Real + A = zeros(T,2,2) + A[1,1] = m + A[2,2] = m + b = zeros(T,2) + transform!(X,A,b) +end + +function LinearAlgebra.scale!(X::GraphEmbedding) + L = mean(edge_length(X)) + if L != 0 + scale!(X,1/L) + else + @warn "Cannot scale by 0. No action taken." + end + nothing +end + + +""" +`scale(G,m)` multiplies all coordinates in the graph's drawing by +`m`. If `m` is omitted, the drawing is rescaled so that the average +length of an edge equals 1. +""" +function scale(G::SimpleGraph, m::T) where T<:Real + X = get_embedding_direct(G) + scale!(X,m) +end + +function scale(G::SimpleGraph) + X = get_embedding_direct(G) + scale!(X) +end + + + +function rotate!(X::GraphEmbedding, theta) + A = zeros(Float64,2,2) + A[1,1] = cos(theta) + A[1,2] = -sin(theta) + A[2,1] = sin(theta) + A[2,2] = cos(theta) + b = zeros(Float64,2) + transform!(X,A,b) + return +end + +""" +`rotate(G,theta)` rotate's the graph's drawing by the angle +`theta`. +""" +rotate(G::SimpleGraph, theta::T) where T<:Real = rotate!(get_embedding_direct(G),theta) + + +function translate!(X::GraphEmbedding, b::Vector{T}) where T + #A = eye(T,2) + A = Matrix{T}(I,2,2) + transform!(X,A,b) + return +end + +""" +`translate(G,b)` translates the graph's drawing by the vector `b`; +that is, every point `p` in the drawing is replaced by `p+b`. +""" +translate(G::SimpleGraph, b::Vector{T}) where T<:Real = translate!(get_embedding_direct(G),b) + + +function recenter!(X::GraphEmbedding) + b = [0;0] + + for v in X.G.V + b += X.xy[v] + end + b /= NV(X.G) + translate!(X,-b) +end + +""" +`recenter(G)` translates the graph's drawing so that the center of mass +of the vertices is at the origin. +""" +recenter(G::SimpleGraph) = recenter!(get_embedding_direct(G)) + + +""" +`draw_labels(G)` will add the vertices names to the drawing +window. + +The results are usually ugly. One can try increasing the size of the +vertices (see `set_vertex_size`). The font size of the labels can be +specified using `draw_labels(G,FontSize)`. The default is 10. +""" +function draw_labels(G::SimpleGraph, FontSize::Int=10) + X = get_embedding_direct(G) + draw_labels(X,FontSize) + nothing +end + + +function draw_labels(X::GraphEmbedding, FontSize::Int) + for v in X.G.V + x,y = X.xy[v] + text(x,y,string(v),fontsize=FontSize) + end +end + + +""" +`getxy(G)` returns a copy of the dictionary mapping vertices to their +x,y-coordinates. This is a way of saving an embedding. +""" +getxy(G::SimpleGraph) = deepcopy(get_embedding_direct(G).xy) + +""" +`edge_length(X::GraphEmbedding,e)` gives the distance between the +embedded endpoints of the edge `e` in the drawing `X` +""" +function edge_length(X::GraphEmbedding, e) + p1 = X.xy[e[1]] + p2 = X.xy[e[2]] + return norm(p1-p2) +end + +""" +`edge_length(X::GraphEmbedding)` returns an array containing the +lengths of the edges in this drawing. +""" +function edge_length(X::GraphEmbedding) + EE = elist(X.G) + + return [ edge_length(X,e) for e in EE ] +end + + +#include("distxy.jl") # distxy! has been replaced by stress! +include("graffle.jl") +include("spectral.jl") +include("geogebra.jl") +include("embedded-graphs.jl") diff --git a/src/embedding/KnightTourDrawing.jl b/src/embedding/KnightTourDrawing.jl new file mode 100644 index 0000000..77cf1f4 --- /dev/null +++ b/src/embedding/KnightTourDrawing.jl @@ -0,0 +1,58 @@ +using SimpleGraphs +using PyPlot + + +""" +`KnightTourDrawing(r,c)` illustrates a Knight's tour of an r-by-c +chessboard returning `true` if the tour exists and `false` if not. +""" +function KnightTourDrawing(r::Int=6, c::Int=6)::Bool + G = Knight(r,c) + println("Searching for a Hamiltonian cycle in an $r-by-$c Knight's move graph") + tic() + h = hamiltonian_cycle(G) + println("Finished") + toc() + + if length(h)==0 + println("Sorry. This graph is not Hamiltonian") + return false + end + + T = vertex_type(G) + H = SimpleGraph{T}() + for v in G.V + add!(H,v) + end + + for k=1:NV(H)-1 + add!(H,h[k],h[k+1]) + end + + add!(H,h[1],h[end]) + + clf() + + xy = getxy(H) + + for v in H.V + xy[v] = collect(v) + end + + embed(H,xy) + + draw(H) + + for a=0:c + plot( [0.5, r+0.5], [a+0.5, a+0.5], color="black", linestyle=":") + end + + for b=0:r + plot( [b+0.5, b+0.5], [0.5, c+0.5], color="black", linestyle=":") + end + + axis([-1,r+1,-1,c+1]) + + return true + +end diff --git a/src/embedding/distxy.jl b/src/embedding/distxy.jl new file mode 100644 index 0000000..64e15e2 --- /dev/null +++ b/src/embedding/distxy.jl @@ -0,0 +1,120 @@ +using Optim +export distxy!, demo_distxy + +# Grab the 2n coordinates as a single vector +function vector_out(X::GraphEmbedding) + vv = vlist(X.G) + n = length(vv) + + x = zeros(2*n) + for k=1:n + v = X.xy[vv[k]] + x[2*k-1] = v[1] + x[2*k] = v[2] + end + + return x +end + +# Inverse of vector_out +function vector_in!(X::GraphEmbedding, x::Vector) + vv = vlist(X.G) + n = length(vv) + + for k=1:n + v = vv[k] + X.xy[v] = x[2*k-1 : 2*k] + end + + return +end + + +function distxy!(X::GraphEmbedding, + nits::Integer=0, verbose::Bool=true) + D = dist_matrix(X.G) + x0 = vector_out(X) + + function score(x::Vector) + nn = length(x) + n = NV(X.G) + s = 0. + + for u=1:(n-1) + pu = x[2*u-1 : 2*u] + for v = (u+1):n + pv = x[2*v-1 : 2*v] + duv = D[u,v] + term = (duv - norm(pu-pv))^2 / duv^(1.75) + s += term + end + end + return s + end + + + if nits > 0 + res = optimize(score,x0,iterations=nits) + else + res = optimize(score,x0) + end + + + x1 = res.minimizer + vector_in!(X,x1) + + if verbose + msg = "$(res.f_calls) function calls\t" * + "score = $(score(x1))" + @info msg + end + + return score(x1) +end + +function distxy!(X::GraphEmbedding, + tol::Float64=0.001, verb::Bool=true) + x0 = Inf + x1 = 0.0 + # if verb + # tic() + # end + while true + x1 = distxy!(X,0,verb) + if abs(x1-x0)/x0 < tol + break + end + x0 = x1 + end + # if verb + # toc() + # end + return x1 +end + +""" +`demo_distxy(G,tol=1e-3)` presents an animation showing the evolving +drawing found by `distxy!`. +""" +function demo_distxy(G::SimpleGraph=BuckyBall(), tol=1e-3) + return demo_distxy(GraphEmbedding(G), tol) +end + +function demo_distxy(X::GraphEmbedding,tol::Real=1e-3) + tic() + x0 = 1.e10 + figure(1) + while true + clf() + draw(X) + title("Embedding with distxy!") + x1 = distxy!(X,0) + if abs(x1-x0)/x0 < tol + break + end + x0 = x1 + end + title("Finished") + toc() + return X +end diff --git a/src/embedding/embedded-graphs.jl b/src/embedding/embedded-graphs.jl new file mode 100644 index 0000000..f859409 --- /dev/null +++ b/src/embedding/embedded-graphs.jl @@ -0,0 +1,43 @@ +export Spindle + +""" +`Spindle()` returns the Moser spindle graph. This is a seven-vertex +unit distance graph with chromatic number equal to 4. +""" +function Spindle() + G = IntGraph(7) + edges = [ + 1 2 + 1 3 + 2 3 + 2 4 + 3 4 + 1 5 + 1 6 + 5 6 + 5 7 + 6 7 + 4 7 ] + add_edges!(G,edges) + + d = Dict{Int,Vector}() + a = sqrt(3)/2 + + pts = [ 0 1/2 -1/2 0 ; 0 a a 2a ] + + theta = acos(5/6)/2 + R = [ cos(theta) -sin(theta); sin(theta) cos(theta) ] + + p1 = R*pts + for k=1:4 + d[k] = p1[:,k] + end + + p2 = R'*pts + for k=5:7 + d[k] = p2[:,k-3] + end + embed(G,d) + SimpleGraphs.name(G,"Moser Spindle") + return G +end diff --git a/src/embedding/geogebra.jl b/src/embedding/geogebra.jl new file mode 100644 index 0000000..3be7cc6 --- /dev/null +++ b/src/embedding/geogebra.jl @@ -0,0 +1,84 @@ +export geogebra + +""" +`geogebra(G,file_name)` takes a `SimpleGraph` and writes out a +script to produce a drawing of this graph in GeoGebra. + +Here is the secret sauce to make this work. + +* Run `geogebra(G,file_name)` to save the script to `file_name`. By + default, the file name is `geogebra.txt`. +* Copy the contents of `file_name` to the clipboard. +* Create a new GeoGebra document. +* In the **Input** zone at the bottom, enter `Button[]` to create a new + button. +* Right click the button and select the **Object Properties...** menu option. +* Go to the **Scripting** tab and paste in the copied commands. +* Save and close the properties window. +* Press the newly created button. + +Some properties of the vertices can be specified in this function with +named parameters as follows: + +* `vertex_labels`: If set to `true`, the vertices in the drawing have + labels. Default is `false` (no labels drawn). Note that the color of + the labels matches the color of the vertices, so if you set the + color to `white` the labels will be invisible. +* `vertex_color`: Use this to specify the fill color of the + vertices. The default is `"black"`. +* `vertex_colors`: This is a dictionary mapping vertices to color + names. Vertices are given the color specified. If, however, a vertex + is missing from this dictionary, then `vertex_color` is used + instead. +* `vertex_size`: Use this to specify the radius of the vertices. The + default is `3`. +""" +function geogebra(G::SimpleGraph, + file_name::String="geogebra.txt"; + vertex_labels::Bool=false, + vertex_color::String="black", + vertex_colors::Dict=Dict(), + vertex_size::Int=3 + ) + X = get_embedding_direct(G) + VV = vlist(X.G) + n = NV(X.G) + F = open(file_name,"w") + + for i=1:n + v = VV[i] + vname = string(v) + (x,y) = X.xy[v] + x = round(x,3) + y = round(y,3) + println(F,"v_{$i} = CopyFreeObject[Point[{$x,$y}]]") + println(F,"SetPointSize[v_{$i}, $(string(vertex_size))]") + + color = "" + try + color = vertex_colors[v] + catch + color = vertex_color + end + println(F,"SetColor[v_{$i}, \"$color\"]") + + + println(F,"SetCaption[v_{$i}, \"$vname\"]") + println(F,"ShowLabel[v_{$i}, $(string(vertex_labels))]") + end + + for i=1:n-1 + u = VV[i] + for j=i+1:n + v = VV[j] + if has(X.G,u,v) + println(F,"e_{$i,$j} = Segment[v_{$i},v_{$j}]") + println(F,"ShowLabel[e_{$i,$j}, false]") + end + end + end + println(F, "ShowAxes[false]") + println(F, "ShowGrid[false]") + println(F, "Execute[{\"Delete[button1]\"}]") + close(F) +end diff --git a/src/embedding/graffle.jl b/src/embedding/graffle.jl new file mode 100644 index 0000000..4afc6ca --- /dev/null +++ b/src/embedding/graffle.jl @@ -0,0 +1,182 @@ +using LightXML + +export graffle + +function bounds(X::GraphEmbedding) + + xmin = Inf + ymin = Inf + xmax = -Inf + ymax = -Inf + + for pt in values(X.xy) + x,y = pt[1],pt[2] + if x < xmin + xmin = x + end + + if x > xmax + xmax = x + end + + if y < ymin + ymin = y + end + + if y > ymax + ymax = y + end + + end + + return (xmin, xmax, ymin, ymax) +end + +function make_scaler(X::GraphEmbedding) + (xmin, xmax, ymin, ymax) = bounds(X) + + f(x,y) = ( round(Int,72*(x-xmin+0.5)), round(Int,72*(ymax-y+0.5))) + + return f +end + + + +function add_key!(dict_node::XMLElement, key::AbstractString) + x = new_child(dict_node,"key") + add_text(x,key) +end + +function add_value!(dict_node::XMLElement, + val_type::AbstractString, val::AbstractString) + x = new_child(dict_node,val_type) + add_text(x,val) +end + + +function add_key_value!(dict_node::XMLElement, + key::AbstractString, + val_type::AbstractString, + val::AbstractString) + add_key!(dict_node,key) + add_value!(dict_node,val_type, val) +end + +function add_dict!(node::XMLElement, key::AbstractString) + add_key!(node, key) + x = new_child(node,"dict") + return x +end + + + + +""" +`graffle(G::SimpleGraph, filename="julia.graffle",rad=9)` creates +an OmniGraffle document of this drawing. + +* `G` is the graph +* `filename` is the name of the OmniGraffle document (be sure to end with `.graffle`) +* `rad` is the radius of the circles representing the vertices +""" +function graffle(G::SimpleGraph, filename="julia.graffle", rad::Int=9) + + X = get_embedding_direct(G) + + # minimal header + xdoc = XMLDocument() + + xtop = create_root(xdoc,"plist") + outer = new_child(xtop, "dict") + + add_key_value!(outer,"GraphDocumentVersion", "integer", "12") + add_key!(outer,"AutoAdjust") + new_child(outer,"true") + + add_key!(outer,"GraphicsList") + + glist = new_child(outer,"array") + + #return xdoc + + # vertices and edges are children of "glist" + + VV = vlist(X.G) + n = NV(X.G) + + lookup = Dict{Any,Int}(zip(VV,1:n)) + + + + f = make_scaler(X) + + for v in VV + k = lookup[v] + vtx = new_child(glist,"dict") + + # Location + pt = X.xy[v] + x,y = f(pt[1],pt[2]) + location = "{{$x,$y},{$(rad),$(rad)}}" + add_key_value!(vtx,"Bounds", "string", location) + + # Class + add_key_value!(vtx, "Class", "string", "ShapedGraphic") + + # ID + add_key_value!(vtx, "ID", "integer", string(k)) + + # Layer + add_key_value!(vtx, "Layer", "integer", "0") + + # Shape + add_key_value!(vtx, "Shape", "string", "Circle") + + # Style + a = add_dict!(vtx,"Style") + b = add_dict!(a,"shadow") + add_key_value!(b,"Draws", "string", "NO") + end + + + EE = elist(X.G) + id = n + + for e in EE + a = lookup[e[1]] + b = lookup[e[2]] + id += 1 + + edge = new_child(glist,"dict") + + # Class + add_key_value!(edge,"Class", "string", "LineGraphic") + + # ID + add_key_value!(edge,"ID", "integer", string(id)) + + # Layer + add_key_value!(edge, "Layer", "integer", string(1)) + + # Head + h = add_dict!(edge,"Head") + add_key_value!(h, "ID", "integer", string(a)) + + # Tail + t = add_dict!(edge,"Tail") + add_key_value!(t, "ID", "integer", string(b)) + + + # Style + a = add_dict!(edge,"Style") + b = add_dict!(a,"shadow") + add_key_value!(b,"Draws", "string", "NO") + + + end + + + save_file(xdoc, filename) + + return filename +end diff --git a/src/embedding/my_spring.jl b/src/embedding/my_spring.jl new file mode 100644 index 0000000..5efa112 --- /dev/null +++ b/src/embedding/my_spring.jl @@ -0,0 +1,90 @@ +# THIS FILE ADAPTED FROM GraphLayout.jl +# see https://github.com/IainNZ/GraphLayout.jl + + +""" + Use the spring/repulsion model of Fruchterman and Reingold (1991): + Attractive force: f_a(d) = d^2 / k + Repulsive force: f_r(d) = -k^2 / d + where d is distance between two vertices and the optimal distance + between vertices k is defined as C * sqrt( area / num_vertices ) + where C is a parameter we can adjust + + Arguments: + adj_matrix Adjacency matrix of some type. Non-zero of the eltype + of the matrix is used to determine if a link exists, + but currently no sense of magnitude + C Constant to fiddle with density of resulting layout + MAXITER Number of iterations we apply the forces + INITTEMP Initial "temperature", controls movement per iteration +""" +function layout_spring_adj(adj_matrix::Array{T,2}; C=2.0, MAXITER=100, INITTEMP=2.0) where T + + size(adj_matrix, 1) != size(adj_matrix, 2) && error("Adj. matrix must be square.") + N = size(adj_matrix, 1) + + # Initial layout is random on the square [-1,+1]^2 + locs_x = 2*rand(N) .- 1.0 + locs_y = 2*rand(N) .- 1.0 + + # The optimal distance bewteen vertices + K = C * sqrt(4.0 / N) + + # Store forces and apply at end of iteration all at once + force_x = zeros(N) + force_y = zeros(N) + + # Iterate MAXITER times + @inbounds for iter = 1:MAXITER + # Calculate forces + for i = 1:N + force_vec_x = 0.0 + force_vec_y = 0.0 + for j = 1:N + i == j && continue + d_x = locs_x[j] - locs_x[i] + d_y = locs_y[j] - locs_y[i] + d = sqrt(d_x^2 + d_y^2) + if adj_matrix[i,j] != zero(eltype(adj_matrix)) || adj_matrix[j,i] != zero(eltype(adj_matrix)) + # F = d^2 / K - K^2 / d + F_d = d / K - K^2 / d^2 + else + # Just repulsive + # F = -K^2 / d^ + F_d = -K^2 / d^2 + end + # d / sin θ = d_y/d = fy/F + # F /| dy fy -> fy = F*d_y/d + # / | cos θ = d_x/d = fx/F + # /--- -> fx = F*d_x/d + # dx fx + force_vec_x += F_d*d_x + force_vec_y += F_d*d_y + end + force_x[i] = force_vec_x + force_y[i] = force_vec_y + end + # Cool down + TEMP = INITTEMP / iter + # Now apply them, but limit to temperature + for i = 1:N + force_mag = sqrt(force_x[i]^2 + force_y[i]^2) + scale = min(force_mag, TEMP)/force_mag + locs_x[i] += force_x[i] * scale + #locs_x[i] = max(-1.0, min(locs_x[i], +1.0)) + locs_y[i] += force_y[i] * scale + #locs_y[i] = max(-1.0, min(locs_y[i], +1.0)) + end + end + + # Scale to unit square + min_x, max_x = minimum(locs_x), maximum(locs_x) + min_y, max_y = minimum(locs_y), maximum(locs_y) + function scaler(z, a, b) + 2.0*((z - a)/(b - a)) - 1.0 + end + locs_x = map(z -> scaler(z, min_x, max_x), locs_x) + locs_y = map(z -> scaler(z, min_y, max_y), locs_y) + + return locs_x,locs_y +end diff --git a/src/embedding/my_stress.jl b/src/embedding/my_stress.jl new file mode 100644 index 0000000..152d4c9 --- /dev/null +++ b/src/embedding/my_stress.jl @@ -0,0 +1,171 @@ +# THIS FILE TAKEN FROM GraphLayout AND EDITED SO IT WORKS + +""" +Compute graph layout using stress majorization + +Inputs: + + δ: Matrix of pairwise distances + p: Dimension of embedding (default: 2) + w: Matrix of weights. If not specified, defaults to + w[i,j] = δ[i,j]^-2 if δ[i,j] is nonzero, or 0 otherwise + X0: Initial guess for the layout. Coordinates are given in rows. + If not specified, default to random matrix of Gaussians + +Additional optional keyword arguments control the convergence of the algorithm +and the additional output as requested: + + maxiter: Maximum number of iterations. Default: 400size(X0, 1)^2 + abstols: Absolute tolerance for convergence of stress. + The iterations terminate if the difference between two + successive stresses is less than abstol. + Default: √(eps(eltype(X0)) + reltols: Relative tolerance for convergence of stress. + The iterations terminate if the difference between two + successive stresses relative to the current stress is less than + reltol. Default: √(eps(eltype(X0)) + abstolx: Absolute tolerance for convergence of layout. + The iterations terminate if the Frobenius norm of two successive + layouts is less than abstolx. Default: √(eps(eltype(X0)) + verbose: If true, prints convergence information at each iteration. + Default: false + returnall: If true, returns all iterates and their associated stresses. + If false (default), returns the last iterate + +Output: + + The final layout X, with coordinates given in rows, unless returnall=true. + +Reference: + + The main equation to solve is (8) of: + + @incollection{ + author = {Emden R Gansner and Yehuda Koren and Stephen North}, + title = {Graph Drawing by Stress Majorization} + year={2005}, + isbn={978-3-540-24528-5}, + booktitle={Graph Drawing}, + seriesvolume={3383}, + series={Lecture Notes in Computer Science}, + editor={Pach, J\'anos}, + doi={10.1007/978-3-540-31843-9_25}, + publisher={Springer Berlin Heidelberg}, + pages={239--250}, + } +""" +function my_layout_stressmajorize_adj(δ, p::Int=2, w=nothing, X0=randn(size(δ, 1), p); + maxiter = 400size(X0, 1)^2, abstols=√(eps(eltype(X0))), + reltols=√(eps(eltype(X0))), abstolx=√(eps(eltype(X0))), + verbose = false, returnall = false) + + @assert size(X0, 2)==p + + if w==nothing + w = δ.^-2 + w[.!isfinite.(w)] .= 0 + end + + @assert size(X0, 1)==size(δ, 1)==size(δ, 2)==size(w, 1)==size(w, 2) + Lw = weightedlaplacian(w) + pinvLw = pinv(Lw) + newstress = stress(X0, δ, w) + Xs = Matrix[X0] + stresses = [newstress] + itcounter = 0 + for iter = 1:maxiter + itcounter = iter + #TODO the faster way is to drop the first row and col from the iteration + X = pinvLw * (LZ(X0, δ, w)*X0) + @assert all(isfinite.(X)) + newstress, oldstress = stress(X, δ, w), newstress + verbose && info("""Iteration $iter + Change in coordinates: $(vecnorm(X - X0)) + Stress: $newstress (change: $(newstress-oldstress)) + """) + push!(Xs, X) + push!(stresses, newstress) + abs(newstress - oldstress) < reltols * newstress && break + abs(newstress - oldstress) < abstols && break + norm(X - X0) < abstolx && break + X0 = X + end + itcounter == maxiter && warn("Maximum number of iterations reached without convergence") + returnall ? (Xs, stresses) : Xs[end] +end + +""" +Stress function to majorize + +Input: + X: A particular layout (coordinates in rows) + d: Matrix of pairwise distances + w: Weights for each pairwise distance + +See (1) of Reference +""" +function stress(X, d=fill(1.0, size(X, 1), size(X, 1)), w=nothing) + s = 0.0 + n = size(X, 1) + if w==nothing + w = d.^-2 + w[.!isfinite.(w)] = 0 + end + @assert n==size(d, 1)==size(d, 2)==size(w, 1)==size(w, 2) + for j=1:n, i=1:j-1 + s += w[i, j] * (norm(X[i,:] - X[j,:]) - d[i,j])^2 + end + @assert isfinite.(s) + s +end + +function _checksquare(M::Matrix) + r,c = size(M) + return r +end +""" +Compute weighted Laplacian given ideal weights w + +Lʷ defined in (4) of the Reference +""" +function weightedlaplacian(w) + n = _checksquare(w) + T = eltype(w) + Lw = zeros(T, n, n) + for i=1:n + D = zero(T) + for j=1:n + i==j && continue + Lw[i, j] = -w[i, j] + D += w[i, j] + end + Lw[i, i]=D + end + Lw +end + +""" +Computes L^Z defined in (5) of the Reference + +Input: Z: current layout (coordinates) + d: Ideal distances (default: all 1) + w: weights (default: d.^-2) +""" +function LZ(Z, d, w) + n = size(Z, 1) + L = zeros(n, n) + for i=1:n + D = 0.0 + for j=1:n + i==j && continue + nrmz = norm(Z[i,:] - Z[j,:]) + nrmz==0 && continue + δ = w[i, j] * d[i, j] + L[i, j] = -δ/nrmz + D -= -δ/nrmz + end + L[i, i] = D + end + @assert all(isfinite.(L)) + L +end diff --git a/src/embedding/spectral.jl b/src/embedding/spectral.jl new file mode 100644 index 0000000..ebeaaf2 --- /dev/null +++ b/src/embedding/spectral.jl @@ -0,0 +1,29 @@ +export spectral! + +""" +`spectral!(X::GraphEmbedding)` gives the graph held in `X` an +embedding based on the eigenvectors of the Laplacian matrix of the +graph. Specifically, the `x`-coordinates come from the eigenvector +associated with the second smallest eigenvalue, and the +`y`-coordinates come from the eigenveector associated with the third +smallest. + +This may also be invoked as `spectral!(X,xcol,ycol)` to choose other +eigenvectors to use for the x and y coordinates of the embedding. +""" +function spectral!(X::GraphEmbedding,xcol::Int=2,ycol::Int=3) + L = laplace(X.G) + EV = eigvecs(L) + x = EV[:,xcol] + y = EV[:,ycol] + + VV = vlist(X.G) + n = length(VV) + + for k=1:n + v = VV[k] + X.xy[v] = [x[k],y[k]] + end + scale!(X) + return X +end