From 8c56d924ddf230262ac6acdd74c2c91bdeca35dc Mon Sep 17 00:00:00 2001 From: Lionel Zoubritzky Date: Fri, 10 Jul 2020 23:00:48 +0200 Subject: [PATCH] First commit --- .travis.yml | 3 + Project.toml | 9 +- README.md | 76 ++++ src/PeriodicGraphs.jl | 985 +++++++++++++++++++++++++++++++++++++++++- test/runtests.jl | 343 ++++++++++++++- 5 files changed, 1411 insertions(+), 5 deletions(-) diff --git a/.travis.yml b/.travis.yml index c1945b0..37c12a0 100644 --- a/.travis.yml +++ b/.travis.yml @@ -23,3 +23,6 @@ jobs: exclude: - arch: x86 os: osx +codecov: true +after_success: +- julia -e 'import PeriodicGraphs; cd(joinpath(dirname(pathof(PeriodicGraphs)), "..")); using Pkg; Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())' diff --git a/Project.toml b/Project.toml index 0d7ab53..43d2edc 100644 --- a/Project.toml +++ b/Project.toml @@ -3,11 +3,18 @@ uuid = "18c5b727-b240-4874-878a-f2e242435bab" authors = ["Lionel Zoubritzky lionel.zoubritzky@gmail.com"] version = "0.1.0" +[deps] +LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + [compat] julia = "1" [extras] +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test"] +test = ["Test", "Random"] diff --git a/README.md b/README.md index c3c8d6e..3214856 100644 --- a/README.md +++ b/README.md @@ -2,3 +2,79 @@ [![Build Status](https://travis-ci.com/Liozou/PeriodicGraphs.jl.svg?branch=master)](https://travis-ci.com/Liozou/PeriodicGraphs.jl) [![Build Status](https://ci.appveyor.com/api/projects/status/github/Liozou/PeriodicGraphs.jl?svg=true)](https://ci.appveyor.com/project/Liozou/PeriodicGraphs-jl) +[![codecov](https://codecov.io/gh/Liozou/PeriodicGraphs.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/Liozou/PeriodicGraphs.jl) + + +This module allows to manipulate `N`-dimensional periodic graphs, using the new type `PeriodicGraph{N}`. This is a subtype of `AbstractGraph{Int}` from [LightGraphs.jl](https://github.com/JuliaGraphs/LightGraphs.jl/) and it extends its API. + +The main difference with a simple graph is the notion of offset. Each vertex, of type `PeriodicVertex{N}` is uniquely defined using a numeric identifier (a positive integer, like for any simple graph) associated with the offset of the cell in which the designated vertex is, compared to a fixed reference cell. For instance, all vertices in the reference cell have a zero offset, and may be built like so: +```julia +julia> PeriodicVertex{1}(1, (0)) +PeriodicVertex{1}(1, (0)) + +julia> PeriodicVertex{4}(2) # shorthand for the vertices of the reference cell +PeriodicVertex{4}(2, (0,0,0,0)) +``` + +For example, in 2D (`N = 2`), the vertex which is the representative of vertex `3` in the cell with offset `(0, 1)` (that is, on the same plane as the reference cell along the x-axis, and just next to the reference cell positively along the y-axis) is defined as `PeriodicVertex{2}(3, (0,1))`. + +An edge, of type `PeriodicEdge{N}`, is defined by its representative starting from the reference cell. Hence, it is uniquely defined by the identifier of the source vertex (in the reference cell) and the destination vertex. It can equivalently be defined by the identifiers of both source and destination vertices, and the offset of the destination, the source being in the reference cell. For example, in 3D, the edge between vertex `1` and vertex `4` in cell `(0, 1, 0)` is `PeriodicEdge{3}(1, PeriodicVertex{3}(4, (0,1,0)))`, or, equivalently, `PeriodicEdge{3}(1, 4, (0,1,0))`. More examples: +```julia +julia> PeriodicEdge{4}(2, PeriodicVertex{4}(3, (0,1,0,0))) +PeriodicEdge{4}(2, 3, (0,1,0,0)) + +julia> LightGraphs.src(PeriodicEdge{3}(5, 6, (1,0,2))) +5 + +julia> LightGraphs.dst(PeriodicEdge{3}(5, 6, (1,0,2))) +6 + +julia> PeriodicGraphs.ofs(PeriodicEdge{3}(5, 6, (1,0,2))) +3-element StaticArrays.SVector{3,Int64} with indices SOneTo(3): + 1 + 0 + 2 + +julia> PeriodicEdge{2}(5, 5, (0,0)) +ERROR: LoopException: a loop from vertex 5 to itself in the same unit cell is a forbidden edges. Maybe the offset is wrong? +``` + +Note that loops (that is, edges of the form `(u, u, (0,0,0,...,0))`) are forbidden in the current version of PeriodicGraphs.jl. + +Finally, `N`-periodic graphs, represented by the type `PeriodicGraph{N}`, are defined by the number of vertices in the reference cell and the set of edges starting from the reference cell. When the number of vertices is simply the highest number appearing in the list of edges, it can be omitted. Periodic graphs can be built with several methods: +```julia +julia> PeriodicGraph{5}() # create the empty graph +PeriodicGraph{5}(0, PeriodicEdge{5}[]) + +julia> PeriodicGraph{1}(4) # create a graph with 4 vertices but no edge +PeriodicGraph{1}(4, PeriodicEdge{1}[]) + +julia> PeriodicGraph(2, PeriodicEdge{4}[(1, 1, (0,0,1,1)), (1, 1, (0,1,0,-1))]) # the dimension can be inferred from the list of edges +PeriodicGraph{4}(2, PeriodicEdge{4}[(1, 1, (0,0,1,1)), (1, 1, (0,1,0,-1))]) + +julia> PeriodicGraph{3}(PeriodicEdge{3}[(1, 3, (0,1,0)), (2, 2, (0,0,-1)), (1, 2, (1,0,0)), (2, 3, (0,1,1))]) +PeriodicGraph3D(3, PeriodicEdge3D[(1, 2, (1,0,0)), (1, 3, (0,1,0)), (2, 2, (0,0,1)), (2, 3, (0,1,1))]) + +julia> PeriodicGraph("3 1 2 1 0 0 1 3 0 1 0 2 2 0 0 1 2 3 0 1 1") # compact representation of the previous graph +PeriodicGraph3D(3, PeriodicEdge3D[(1, 2, (1,0,0)), (1, 3, (0,1,0)), (2, 2, (0,0,1)), (2, 3, (0,1,1))]) + +julia> string(ans) # to obtain the compact representation from the graph +"3 1 2 1 0 0 1 3 0 1 0 2 2 0 0 1 2 3 0 1 1" +``` + +Note that `PeriodicGraph`s are undirected: for this reason, any edge in the graph of the form `(u, v, (x1, x2, ..., xN))` has a reverse edge in the graph of the form `(v, u, (-x1, -x2, ..., -xN))`. For this reason, calling `LightGraphs.edges` on a `PeriodicGraph` yields an iterator to the edges of the graph that will only show the edges under the canonical form `(u, v, ofs)` with either `u < v` or `u == v && ofs > zero(ofs)`. + +Neighbors can be obtained using LightGraphs's `neighbors` function, with the same API. It will return the list of neighbors of the given vertex, assuming the vertex is in the reference cell. Each given is given as a `PeriodicVertex`, which contains the offset compared to the reference cell. For instance: +```julia +julia> g = PeriodicGraph("3 1 2 1 0 0 1 3 0 1 0 2 2 0 0 1 2 3 0 1 1") +PeriodicGraph3D(3, PeriodicEdge3D[(1, 2, (1,0,0)), (1, 3, (0,1,0)), (2, 2, (0,0,1)), (2, 3, (0,1,1))]) + +julia> neighbors(g, 2) +4-element Vector{PeriodicVertex3D}: + (1, (-1,0,0)) + (2, (0,0,-1)) + (2, (0,0,1)) + (3, (0,1,1)) +``` + +For convenience, aliases are exported for 2D and 3D (`N = 2` and `N = 3`) under the names `PeriodicGraph2D`, `PeriodicGraph3D`, `PeriodicEdge2D`, `PeriodicEdge3D`, `PeriodicVertex2D` and `PeriodicVertex3D`. diff --git a/src/PeriodicGraphs.jl b/src/PeriodicGraphs.jl index 7ce0b33..773acfd 100644 --- a/src/PeriodicGraphs.jl +++ b/src/PeriodicGraphs.jl @@ -1,5 +1,988 @@ module PeriodicGraphs -# Write your package code here. +using LinearAlgebra +using StaticArrays +using LightGraphs +export PeriodicVertex, PeriodicEdge, PeriodicGraph, + PeriodicVertex2D, PeriodicEdge2D, PeriodicGraph2D, + PeriodicVertex3D, PeriodicEdge3D, PeriodicGraph3D, + coordination_sequence, cellgraph, periodiccellgraph, + offset_representatives!, swap_axes!, find_edges +import Base: (==), isless, convert, show, showerror, eltype, iterate, zero, + length, in, ndims, print, cmp +import Base.Order: Forward, Lt + +""" + PeriodicVertex{N} + +Vertex type for an `N`-periodic graph. + +A vertex is uniquely determined by the identifier of its representative in the +a fixed initial cell, and the offset of the cell containing the vertex compared +to the the initial cell. +Vertex identifiers start at 1. +""" +struct PeriodicVertex{N} + v::Int + ofs::SVector{N,Int} +end +PeriodicVertex{N}(n::Integer) where {N} = PeriodicVertex{N}(n, zero(SVector{N,Int})) +PeriodicVertex(v::Integer, ofs::Union{SVector{N},NTuple{N}}) where {N} = PeriodicVertex{N}(v, ofs) + +function show(io::IO, x::PeriodicVertex{N}) where N + if get(io, :typeinfo, Any) != PeriodicVertex{N} + print(io, PeriodicVertex{N}) + end + print(io, '(', x.v, ", (", join(x.ofs, ','), ')', ')') +end +function convert(::Type{PeriodicVertex{N}}, (dst, offset)::Tuple{Integer,Any}) where N + PeriodicVertex{N}(dst, offset) +end +function convert(::Type{PeriodicVertex}, (dst, offset)::Tuple{Integer,Union{SVector{N},NTuple{N}}}) where N + PeriodicVertex{N}(dst, offset) +end +function cmp(x::PeriodicVertex{N}, y::PeriodicVertex{N}) where N + c = cmp(x.v, y.v) + iszero(c) || return c + return cmp(x.ofs, y.ofs) +end +isless(x::PeriodicVertex{N}, y::PeriodicVertex{N}) where {N} = cmp(x, y) < 0 +==(x::PeriodicVertex{N}, y::PeriodicVertex{M}) where {N,M} = N == M && iszero(cmp(x, y)) + +const PeriodicVertex2D = PeriodicVertex{2} +const PeriodicVertex3D = PeriodicVertex{3} + + +ZtoN(x::Signed) = -(x<0) + 2*abs(x) +function hash_position((x1,x2,x3)::SVector{3,<:Integer}) + x1 = ZtoN(x1); x2 = ZtoN(x2); x3 = ZtoN(x3) + _b = x1 >= x2 + b1 = _b & (x1 >= x3) + b2 = (!_b) & (x2 >= x3) + b3 = (!b1) & (!b2) + return b1*(x2*(x1 + 1) + x3 + x1^3) + + b2*((x2 + 1)*(x1 + x2 + 1) + x3 + x2^3) + + b3*((x3 + 1)*(2x3 + 1) + x3*(x1 + x3^2) + x2) +end + +function hash_position((x1,x2)::SVector{2,<:Integer}) + x1 = ZtoN(x1); x2 = ZtoN(x2) + b = x1 >= x2 + return b*(x2 + x1^2) + (!b)*(x1 + x2*(x2 + 1) + 1) +end + +""" + hash_position(x::PeriodicVertex{N}, n::Integer) where {N} + +Given a `PeriodicVertex{N}` and the number n of vertex identifiers in a graph, compute a +unique hash for the given vertex. + +This hash function is a bijection between the set of all the vertices of the graph and +the set of positive integers. Its value is an integer between `1+n\\*(2d-1)^N` +(or `1` if `d == 0`) and `n*(2d+1)^N`, where `d = max.(abs.(x.ofs))`. +This means that when one unit cell A is further than another B (for the Manhattan +distance), all vertices in A will have a larger hash than all vertices in B. +""" +function hash_position(x::PeriodicVertex, n::Integer) + return x.v + n*hash_position(x.ofs) +end + +""" + LoopException <: Exception + +Error type for constructing an invalid `PeriodicEdge{N}` of the form (u, u, zeros(Int,N)). +Loops are not expected in the algorithms implemented in PeriodicGraphs.jl. If you +still want to construct them, use the `unsafe_edge{N}` constructor instead of +`PeriodicEdge{N}``. +""" +struct LoopException <: Exception + src::Int +end +showerror(io::IO, e::LoopException) = print(io, "LoopException: a loop from vertex $(e.src) to itself in the same unit cell is a forbidden edges. Maybe the offset is wrong?") +@noinline __throw_loopexception(src) = throw(LoopException(src)) + +""" + unsafe_edge{N} + +Internal constructor for `PeriodicEdge{N}` that bypasses the loop check. +This function is not part of the official API, use with caution. +""" +struct unsafe_edge{N} end + +""" + PeriodicEdge{N} <: LightGraphs.SimpleGraphs.AbstractSimpleEdge{Int} + +Edge type for an `N`-periodic graph. + +An edge is uniquely determined by the vertex identifiers of its source and +destination, and the cell offset between the source vertex and the destination vertex. +""" +struct PeriodicEdge{N} <: LightGraphs.SimpleGraphs.AbstractSimpleEdge{Int} + src::Int + dst::PeriodicVertex{N} + function (::Type{unsafe_edge{N}})(src, dst) where N + return new{N}(src, dst) + end +end +unsafe_edge{N}(src, dst, ofs) where {N} = unsafe_edge{N}(src, PeriodicVertex{N}(dst, ofs)) +function PeriodicEdge{N}(src, dst::PeriodicVertex{N}) where N + src == dst.v && all(iszero.(dst.ofs)) && __throw_loopexception(src) + return unsafe_edge{N}(src, dst) +end +PeriodicEdge(src, dst::PeriodicVertex{N}) where {N} = PeriodicEdge{N}(src, dst) +function PeriodicEdge{N}(src, dst, offset) where N + return PeriodicEdge{N}(src, PeriodicVertex{N}(dst, offset)) +end +function PeriodicEdge(src, dst, offset::Union{SVector{N,T},NTuple{N,T}}) where {N,T} + PeriodicEdge{N}(src, dst, offset) +end +function PeriodicEdge{N}((src, dst, offset)::Tuple{Any,Any,Union{SVector{N,T},NTuple{N,T}}}) where {N,T<:Integer} + PeriodicEdge{N}(src, dst, offset) +end +function PeriodicEdge((src, dst, offset)::Tuple{Any,Any,Union{SVector{N,T},NTuple{N,T}}}) where {N,T<:Integer} + PeriodicEdge{N}(src, dst, offset) +end + +function convert(::Type{PeriodicEdge{N}}, (src, dst, offset)::Tuple) where {N} + PeriodicEdge{N}(src, dst, offset) +end +function convert(::Type{PeriodicEdge}, (src, dst, offset)::Tuple{Any,Any,Union{SVector{N,T},NTuple{N,T}}}) where {N,T<:Integer} + PeriodicEdge{N}(src, dst, offset) +end + +function show(io::IO, x::PeriodicEdge{N}) where N + if get(io, :typeinfo, Any) != PeriodicEdge{N} + print(io, PeriodicEdge{N}) + end + print(io, '(', x.src, ", ", x.dst.v, ", (", join(x.dst.ofs, ','), ')', ')') +end + +const PeriodicEdge2D = PeriodicEdge{2} +const PeriodicEdge3D = PeriodicEdge{3} + +function LightGraphs.reverse(e::PeriodicEdge{N}) where N + if e.src == e.dst + return e + end + return unsafe_edge{N}(e.dst.v, e.src, .-e.dst.ofs) +end +LightGraphs.src(e::PeriodicEdge) = e.src +LightGraphs.dst(e::PeriodicEdge) = e.dst.v + +""" + ofs(v::PeriodicVertex) + ofs(e::PeriodicEdge) + +The offset of a vertex or an edge. +""" +ofs(v::PeriodicVertex) = v.ofs +ofs(e::PeriodicEdge) = ofs(e.dst) +isindirectedge(e::PeriodicEdge) = e.src > e.dst.v || (e.src == e.dst.v && e.dst.ofs < zero(e.dst.ofs)) +function cmp(x::PeriodicEdge{N}, y::PeriodicEdge{N}) where N + c = cmp(x.src, y.src) + iszero(c) || return c + return cmp(x.dst, y.dst) +end +isless(x::PeriodicEdge{N}, y::PeriodicEdge{N}) where {N} = cmp(x, y) < 0 +==(x::PeriodicEdge{N}, y::PeriodicEdge{M}) where {N,M} = N == M && iszero(cmp(x, y)) + + +""" + PeriodicGraph{N} <: AbstractGraph{Int} + +Type representing an undirected `N`-periodic graph. +""" +struct PeriodicGraph{N} <: AbstractGraph{Int} + ne::Base.RefValue{Int} + nlist::Vector{Vector{PeriodicVertex{N}}} # For each vertex, the sorted list of its neighbors + directedgestart::Vector{Int} + width::Base.RefValue{Rational{Int}} +end + +const PeriodicGraph2D = PeriodicGraph{2} +const PeriodicGraph3D = PeriodicGraph{3} + +function PeriodicGraph{N}(nv::Integer, t, s) where N + return PeriodicGraph{N}(Ref(nv), t, s, Ref(-1//1)) +end + +""" + PeriodicGraph{N}(n::Integer=0) + +Construct a `PeriodicGraph{N}` with `n` vertices and 0 edges. +""" +function PeriodicGraph{N}(n::Integer = 0) where N + @assert n >= 0 + return PeriodicGraph{N}(0, [PeriodicVertex{N}[] for _ in 1:n], [1 for _ in 1:n]) +end +""" + PeriodicGraph([nv::Integer, ]edge_list::AbstractVector{PeriodicEdge{N}}) + PeriodicGraph{N}([nv::Integer, ]edge_list::AbstractVector{PeriodicEdge{N}}) + +Construct a `PeriodicGraph{N}` from a vector of edges. +If `nv` is unspecified, the number of vertices is the highest that is used +in an edge in `edge_list`. + +### Implementation Notes +This constructor works the fastest when `edge_list` is sorted by the lexical ordering +and does not contain any duplicates. + +## Examples +```julia +julia> el = PeriodicEdge2D[(1, 1, (1,0)), (1, 3, (0,1)), (3, 1, (0,-1))]; + +julia> g = PeriodicGraph(el) +PeriodicGraph2D(3, PeriodicEdge2D[(1, 1, (1, 0)), (1, 3, (0, 1))]) + +julia> ne(g) +2 +``` +""" +function PeriodicGraph{N}(nv::Integer, t::AbstractVector{PeriodicEdge{N}}) where N + sort!(t); unique!(t) + ne = length(t) + g = PeriodicGraph{N}(nv) + for e in t + add_edge!(g, e) + end + return g +end +PeriodicGraph(nv::Integer, t::AbstractVector{PeriodicEdge{N}}) where {N} = PeriodicGraph{N}(nv, t) +function PeriodicGraph{N}(t::AbstractVector{PeriodicEdge{N}}) where N + @static if VERSION < v"1.6-" + isempty(t) && return PeriodicGraph(0, t) + return PeriodicGraph(maximum(max(e.src, e.dst.v) for e in t), t) + else + return PeriodicGraph(maximum(max(e.src, e.dst.v) for e in t; init=0), t) + end +end +PeriodicGraph(t::AbstractVector{PeriodicEdge{N}}) where {N} = PeriodicGraph{N}(t) + + +struct KeyString{T,S<:AbstractString} + x::S + start::Base.RefValue{Int} +end +function KeyString{T}(x) where T + KeyString{T,typeof(x)}(x, Ref(firstindex(x))) +end +iterate(::KeyString, ::Nothing) = nothing +function iterate(k::KeyString{T}, (next_char, idx)) where T + start = idx + state = nothing + while isspace(next_char) + state = iterate(k.x, idx) + state isa Nothing && return nothing + start = idx + next_char, idx = state + end + stop = start + tmp = idx + while !isspace(next_char) + state = iterate(k.x, idx) + state isa Nothing && break + stop = tmp + tmp = idx + next_char, idx = state + end + return (parse(T, SubString(k.x, start:stop)), state) +end +function iterate(k::KeyString) + iterate(k, (' ', k.start[])) +end +function Base.popfirst!(k::KeyString) + next = iterate(k) + next isa Nothing && throw(ArgumentError("collection must be non-empty")) + char, state = next + k.start[] = state isa Nothing ? (ncodeunits(k.x) + 1) : last(state) + return char +end +Base.isempty(k::KeyString) = isempty(SubString(k.x, k.start[])) +Base.IteratorSize(::Type{KeyString}) = SizeUnknown() + +""" + PeriodicGraph(key::AbstractString) + PeriodicGraph{N}(key::AbstractString) + +Construct a `PeriodicGraph{N}` from a `key`, which is a string of whitespace-separated +values of the form +`"N src1 dst1 ofs1_1 ofs1_2 ... ofs1_N src2 dst2 ofs2_1 ofs2_2 ... ofs2_N ... srcm dstm ofsm_1 ofsm_2 ... ofsm_N"` +where `N` is the number of repeating dimensions of the graph, `m` is the number of edges +and for all `i` between `1` and `m`, the number of edges, the `i`-th edge is described as +- `srci`, the vertex identifier of the source vertex, +- `dsti`, the vertex identifier of the destination vertex and +- `(ofsi_1, ofsi_2, ..., ofsi_N)` the offset of the edge. + +This compact representation of a graph can be obtained simply by `print`ing the graph +or with `string`. + +## Examples +```julia +julia> PeriodicGraph("2 1 2 0 0 2 1 1 0 1 1 0 -1") +PeriodicGraph2D(2, PeriodicEdge2D[(1, 1, (0,1)), (1, 2, (-1,0)), (1, 2, (0,0))]) + +julia> PeriodicGraph3D("3 1 1 0 0 1 1 1 0 1 0 1 1 1 0 0") +PeriodicGraph3D(1, PeriodicEdge3D[(1, 1, (0,0,1)), (1, 1, (0,1,0)), (1, 1, (1,0,0))]) + +julia> string(ans) +"3 1 1 0 0 1 1 1 0 1 0 1 1 1 0 0" + +julia> string(PeriodicGraph(ans)) == ans +true +``` +""" +function PeriodicGraph{N}(s::AbstractString) where N + key = KeyString{Int}(s) + M = popfirst!(key) + M != N && throw(DimensionMismatch("Cannot construct a $N-dimensional graph from a $M-dimensional key")) + edgs = PeriodicEdge{N}[] + while !isempty(key) + src = popfirst!(key) + dst = popfirst!(key) + ofs = SVector{N,Int}([popfirst!(key) for _ in 1:N]) + push!(edgs, PeriodicEdge{N}(src, dst, ofs)) + end + return PeriodicGraph{N}(edgs) +end +function PeriodicGraph(s::AbstractString) + N = first(KeyString{Int}(s)) + return PeriodicGraph{N}(s) +end +function PeriodicGraph{N}(g::PeriodicGraph{N}) where N + PeriodicGraph{N}(g.ne[], [copy(x) for x in g.nlist], copy(g.directedgestart)) +end +PeriodicGraph(g::PeriodicGraph{N}) where {N} = PeriodicGraph{N}(g) + + +function show(io::IO, g::PeriodicGraph{N}) where N + if get(io, :typeinfo, Any) != PeriodicGraph{N} + print(io, PeriodicGraph{N}) + end + print(io, '(', nv(g), ',', ' ', collect(edges(g)), ')') +end +function print(io::IO, g::PeriodicGraph{N}) where N + print(io, N) + for e in edges(g) + print(io, ' ', e.src, ' ', e.dst.v, ' ', join(e.dst.ofs, ' ')) + end +end + +""" + ==(g1::PeriodicGraph, g2::PeriodicGraph) + +Test whether two `PeriodicGraph`s have the same representations. Beware that +two `PeriodicGraph`s may be isomorphic but still have different representations. +""" +function ==(g1::PeriodicGraph{N}, g2::PeriodicGraph{M}) where {N,M} + return N == M && nv(g1) == nv(g2) && edges(g1) == edges(g2) +end + +ndims(::PeriodicVertex{N}) where {N} = N +ndims(::PeriodicEdge{N}) where {N} = N +ndims(::PeriodicGraph{N}) where {N} = N + +""" + ndims(::PeriodicVertex) + ndims(::PeriodicEdge) + ndims(::PeriodicGraph) + +Return the number of repeating dimensions. + +## Examples +```julia +julia> ndims(PeriodicVertex{5}(2)) +5 + +julia> ndims(PeriodicEdge3D(1, 3, (0,0,1))) +3 + +julia> ndims(PeriodicGraph2D(7)) +2 +``` +""" +ndims + +LightGraphs.ne(g::PeriodicGraph) = g.ne[] +LightGraphs.nv(g::PeriodicGraph) = length(g.nlist) +LightGraphs.vertices(g::PeriodicGraph) = Base.OneTo(nv(g)) +LightGraphs.edges(g::PeriodicGraph{N}) where {N} = PeriodicEdgeIter{N}(g) +eltype(g::PeriodicGraph{N}) where {N} = PeriodicVertex{N} +LightGraphs.edgetype(::PeriodicGraph{N}) where {N} = PeriodicEdge{N} +function LightGraphs.has_edge(g::PeriodicGraph, s, d) + ((s < 1) | (s > nv(g))) && return false + #=@inbounds=# begin + start = g.directedgestart[s] + lo, hi = s > d ? (1, start-1) : (start, lastindex(g.nlist[s])) + i = searchsortedfirst(g.nlist[s], d, lo, hi, Lt((x,y)->isless(x.v, y))) + return i <= length(g.nlist[s]) && g.nlist[s][i].v == d + end +end + +""" + find_edges(g::PeriodicGraph, s::Int, d::Int) + +Return the set of edges of graph `g` such that the source vertex identifier is `s` +and the destination vertex identifier is `d`. +""" +function find_edges(g::PeriodicGraph, s::Int, d::Int) + ((s < 1) | (s > nv(g))) && return false + #=@inbounds=# begin + start = g.directedgestart[s] + lo, hi = s > d ? (1, start-1) : (start, lastindex(g.nlist[s])) + rng = searchsorted(g.nlist[s], d, lo, hi, Lt((x,y)->isless(x.v, y))) + if s == d.v + rng = (2*first(rng) - last(rng) - 1):last(rng) + end + return g.nlists[rng] + end +end +function LightGraphs.has_edge(g::PeriodicGraph, e::PeriodicEdge) + s, d = e.src, e.dst + ((s < 1) | (s > nv(g))) && return false + #=@inbounds=# begin + start = g.directedgestart[s] + lo, hi = isindirectedge(e) ? (1, start-1) : (start, lastindex(g.nlist[s])) + i = searchsortedfirst(g.nlist[s], d, lo, hi, Forward) + return i <= length(g.nlist[s]) && g.nlist[s][i] == d + end +end +LightGraphs.outneighbors(g::PeriodicGraph, v::Integer) = g.nlist[v] +LightGraphs.inneighbors(g::PeriodicGraph, v::Integer) = outneighbors(g, v) +zero(::Type{PeriodicGraph{N}}) where N = PeriodicGraph{N}(0) +LightGraphs.is_directed(::Type{<:PeriodicGraph}) = false +@static if isdefined(LightGraphs, :has_contiguous_vertices) + @inline LightGraphs.has_contiguous_vertices(::Type{<:PeriodicGraph}) = true +end +LightGraphs.has_vertex(g::PeriodicGraph, v::Integer) = 1 <= v <= nv(g) +function LightGraphs.SimpleGraphs.add_vertices!(g::PeriodicGraph{N}, n::Integer) where N + append!(g.nlist, [PeriodicVertex{N}[] for _ in 1:n]) + append!(g.directedgestart, [1 for _ in 1:n]) + g.width[] = -1 + return n +end +function LightGraphs.SimpleGraphs.add_vertex!(g::PeriodicGraph{N}) where N + push!(g.nlist, PeriodicVertex{N}[]) + push!(g.directedgestart, 1) + g.width[] = -1 + true +end + +function _add_edge!(g::PeriodicGraph, e::PeriodicEdge, ::Val{check}) where check + #=@inbounds=# begin + s, dst = e.src, e.dst + neigh = g.nlist[s] + start = g.directedgestart[s] + indirectedge = isindirectedge(e) + lo, hi = indirectedge ? (1, start-1) : (start, lastindex(neigh)) + i = searchsortedfirst(neigh, dst, lo, hi, Forward) + if check + i <= length(neigh) && neigh[i] == dst && return false + end + g.directedgestart[s] += indirectedge + insert!(neigh, i, dst) + return true + end +end +function LightGraphs.add_edge!(g::PeriodicGraph, e::PeriodicEdge) + (src(e) < 1 || src(e) > nv(g) || dst(e) < 1 || dst(e) > nv(g)) && return false + success = _add_edge!(g, e, Val(true)) && _add_edge!(g, reverse(e), Val(false)) + if success + g.ne[] += 1 + g.width[] = -1 + end + return success +end + +function _rem_edge!(g::PeriodicGraph, e::PeriodicEdge, ::Val{check}) where check + #=@inbounds=# begin + s, dst = e.src, e.dst + neigh = g.nlist[s] + start = g.directedgestart[s] + indirectedge = isindirectedge(e) + lo, hi = indirectedge ? (1, start-1) : (start, lastindex(neigh)) + i = searchsortedfirst(neigh, dst, lo, hi, Forward) + if check + i <= length(neigh) && neigh[i] == dst || return false + end + g.directedgestart[s] -= indirectedge + deleteat!(neigh, i) + return true + end +end +function LightGraphs.rem_edge!(g::PeriodicGraph, e::PeriodicEdge) + (src(e) < 1 || src(e) > nv(g) || dst(e) < 1 || dst(e) > nv(g)) && return false + success = _rem_edge!(g, e, Val(true)) && _rem_edge!(g, reverse(e), Val(false)) + if success + g.ne[] -= 1 + g.width[] = -1 + end + return success +end + + +function LightGraphs.SimpleGraphs.rem_vertices!(g::PeriodicGraph{N}, t::AbstractVector{<:Integer}, keep_order::Bool=false) where N + isempty(t) && return collect(1:nv(g)) + sort!(t) + (first(t) < 1 || last(t) > nv(g)) && throw(ArgumentError("Vertices to be removed must be in the range 1:nv(g).")) + + bt = falses(nv(g)) + bt[t] .= true + + for i in t + g.ne[] -= count((i < x.v || (i == x.v && x.ofs > zero(x.ofs))) for x in g.nlist[i]) + end + + vmap = Int[] + rev_vmap = zeros(Int, nv(g)) + + if keep_order + append!(vmap, collect(1:nv(g))) + deleteat!(vmap, bt) + rev_vmap[vmap] .= 1:length(vmap) + deleteat!(g.nlist, bt) + else + i_next_vertex_to_del = 1 + next_vertex_to_del = t[i_next_vertex_to_del] + t_end = length(t) + g_end = length(g.nlist) + i = 1 + while i <= g_end + if i == next_vertex_to_del + i_next_vertex_to_del += 1 + if i_next_vertex_to_del > t_end + next_vertex_to_del = nv(g) + 1 + else + next_vertex_to_del = t[i_next_vertex_to_del] + end + while t_end >= i_next_vertex_to_del && g_end == t[t_end] + t_end -= 1 + g_end -= 1 + end + if i < g_end + g.nlist[i], g.nlist[g_end] = g.nlist[g_end], g.nlist[i] + push!(vmap, g_end) + rev_vmap[g_end] = length(vmap) + end + g_end -= 1 + else + push!(vmap, i) + rev_vmap[i] = length(vmap) + end + i += 1 + end + end + resize!(g.nlist, g_end) + resize!(g.directedgestart, g_end) + + for i in vertices(g) + neighbors = g.nlist[i] + remove_edges = falses(length(neighbors)) + startoffset = 1 + for (k, x) in enumerate(neighbors) + if bt[x.v] + remove_edges[k] = true + else + neigh = PeriodicVertex{N}(rev_vmap[x.v], x.ofs) + neighbors[k] = neigh + startoffset += isindirectedge(PeriodicEdge{N}(i, neigh)) + end + end + deleteat!(neighbors, remove_edges) + sort!(neighbors) + g.directedgestart[i] = startoffset + end + + g.width[] = -1 + + return vmap +end + +function LightGraphs.SimpleGraphs.rem_vertex!(g::PeriodicGraph, v::Integer) + n = nv(g) + return length(rem_vertices!(g, [v])) == n - 1 +end + + +""" + PeriodicEdgeIter{N} <: AbstractEdgeIter + +Edge iterator type for undirected `N`-periodic graphs. + +The iterator only yields edges in the form `(u, v, ofs)` with either `u < v` or +`u == v && ofs > zero(ofs)`. +This is possible because `PeriodicGraph`s are undirected, hence to each edge +`(u, v, ofs)` in the graph corresponds its reverse edge `(v, u, .-ofs)`. The iterator +thus yields each edge of the graph exactly once. +""" +struct PeriodicEdgeIter{N} <: AbstractEdgeIter + g::PeriodicGraph{N} +end + +eltype(::Type{PeriodicEdgeIter{N}}) where {N} = PeriodicEdge{N} +length(iter::PeriodicEdgeIter) = iter.g.ne[] + +function iterate(iter::PeriodicEdgeIter{N}, (vertex, neigh)) where N + nlists = iter.g.nlist + n = length(nlists) + #=@inbounds=# while vertex <= n + if iszero(neigh) + neigh = iter.g.directedgestart[vertex] + end + neighbors = nlists[vertex] + if neigh > length(neighbors) + vertex += 1 + neigh = 0 + continue + end + return (unsafe_edge{N}(vertex, neighbors[neigh]), (vertex, neigh+1)) + end + return nothing +end + +iterate(iter::PeriodicEdgeIter) = iterate(iter, (1, 0)) + +function in(edge, iter::PeriodicEdgeIter{N}) where N + has_edge(iter.g, edge) +end + +function cmp(it1::PeriodicEdgeIter{N}, it2::PeriodicEdgeIter{N}) where N + n = length(it1) + m = length(it2) + n == m || return cmp(n,m) + n == 0 && return 0 + cmpofs = 0 + e1::PeriodicEdge{N}, st1::Tuple{Int,Int} = iterate(it1) + e2::PeriodicEdge{N}, st2::Tuple{Int,Int} = iterate(it2) + for _ in 1:n-1 + c = cmp((e1.src, e1.dst.v), (e2.src, e2.dst.v)) + if iszero(c) + if iszero(cmpofs) + cmpofs = cmp(e1.dst.ofs, e2.dst.ofs) + end + e1, st1 = iterate(it1, st1) + e2, st2 = iterate(it2, st2) + continue + end + return c + end + c = cmp((e1.src, e1.dst.v), (e2.src, e2.dst.v)) + return iszero(c) ? (iszero(cmpofs) ? cmp(e1.dst.ofs, e2.dst.ofs) : cmpofs) : c +end + +function isless(it1::PeriodicEdgeIter{N}, it2::PeriodicEdgeIter{N}) where N + return cmp(it1, it2) < 0 +end +function ==(it1::PeriodicEdgeIter{N}, it2::PeriodicEdgeIter{N}) where N + return iszero(cmp(it1, it2)) +end + +""" + vertex_permutation(g::PeriodicGraph, vlist) + +Return the `PeriodicGraph` corresponding to `g` with its vertices identifiers +permuted according to `vlist`. `isperm(vlist)` must hold and will not be checked. + +See also `induced_subgraph` for the more general case where `vlist` is not a +permutation. +""" +function vertex_permutation(g::PeriodicGraph{N}, vlist) where N + n = length(vlist) + newvid = Vector{Int}(undef, n) + for i in 1:n + newvid[vlist[i]] = i + end + edges = Vector{Vector{PeriodicVertex{N}}}(undef, n) + startoffsets = [1 for _ in 1:n] + #=@inbounds=# for i in 1:n + startoffset = 1 + neighs = copy(g.nlist[vlist[i]]) + edges[i] = neighs + for j in 1:length(neighs) + dst = neighs[j] + neigh = PeriodicVertex{N}(newvid[dst.v], dst.ofs) + neighs[j] = neigh + startoffsets[i] += isindirectedge(PeriodicEdge{N}(i, neigh)) + end + sort!(neighs) + end + return PeriodicGraph{N}(g.ne[], edges, startoffsets) +end + +function LightGraphs.induced_subgraph(g::PeriodicGraph{N}, vlist::AbstractVector{U}) where {N, U<:Integer} + allunique(vlist) || __throw_unique_vlist() + n = length(vlist) + n == nv(g) && return (vertex_permutation(g, vlist), vlist) + newvid = zeros(Int, nv(g)) + for i in 1:n + newvid[vlist[i]] = i + end + + ne = 0 + edges = Vector{Vector{PeriodicVertex{N}}}(undef, n) + startoffsets = [1 for _ in 1:n] + for i in 1:n + edges[i] = PeriodicVertex{N}[] + startne = ne + for dst in g.nlist[vlist[i]] + v = newvid[dst.v] + iszero(v) && continue + neigh = PeriodicVertex{N}(v, dst.ofs) + push!(edges[i], neigh) + ne += isindirectedge(PeriodicEdge{N}(i, neigh)) + end + startoffsets[i] = 1 + ne - startne + sort!(edges[i]) + end + return (PeriodicGraph{N}(ne, edges, startoffsets), vlist) +end +@noinline __throw_unique_vlist() = throw(ArgumentError("Vertices in subgraph list must be unique")) + +""" + offset_representatives!(g::PeriodicGraph, offsets) + +In-place modifies graph `g` so that the `i`-th vertex of the new initial cell corresponds +to the `i`-th vertex in cell `offsets[i]` compared to the previous initial cell. + +Note that the resulting graph is isomorphic to the initial one, only the representation +has changed. +""" +function offset_representatives!(g::PeriodicGraph{N}, offsets) where N + n = nv(g) + length(offsets) == n || __throw_invalid_offsets() + for i in 1:n + neighs = g.nlist[i] + ofsi = offsets[i] + startoffset = 1 + for j in 1:length(neighs) + x = neighs[j] + neigh = PeriodicVertex{N}(x.v, x.ofs .+ ofsi .- offsets[x.v]) + neighs[j] = neigh + startoffset += isindirectedge(PeriodicEdge{N}(i, neigh)) + end + sort!(neighs) + g.directedgestart[i] = startoffset + g.width[] = -1 + end + g +end +@noinline __throw_invalid_offsets() = throw(ArgumentError("The size of offsets does not match the number of vertices")) + +""" + swap_axes!(g::PeriodicGraph, t::AbstractVector{<:Integer}) + +In-place modifies graph `g` so that the new initial cell corresponds to the previous +one with its axes swapped according to the permutation `t`. + +Note that the resulting graph is isomorphic to the initial one, only the representation +has changed. +""" +function swap_axes!(g::PeriodicGraph{N}, t::AbstractVector{<:Integer}) where N + length(t) == N || __throw_invalid_axesswap() + for i in vertices(g) + neighs = g.nlist[i] + for j in 1:length(neighs) + x = neighs[j] + neigh = PeriodicVertex{N}(x.v, x.ofs[t]) + neighs[j] = neigh + end + sort!(neighs) + # Note: g.width[] is unchanged + end + g +end +@noinline __throw_invalid_axesswap() = throw(ArgumentError("The number of axes must match the dimension of the graph")) + + +function LightGraphs.connected_components(g::PeriodicGraph) + nvg = nv(g) + label = zeros(Int, nvg) + + for u in vertices(g) + label[u] != 0 && continue + label[u] = u + Q = Int[] + push!(Q, u) + #=@inbounds=# while !isempty(Q) + src = popfirst!(Q) + for dst in outneighbors(g, src) + vertex = dst.v + if label[vertex] == 0 + push!(Q, vertex) + label[vertex] = u + end + end + end + end + c, d = LightGraphs.components(label) + return c +end + +""" + graph_width!(g::PeriodicGraph) + +Set the `width` internal field of the graph so that the for all `n` ∈ N\\*, +the `n`-th neighbor of any vertex `v` of the initial cell is in a cell +`(i,j,k)` such that `max(abs.((i,k,k))) ≤ 1 + fld((n - 1), width)`. + +This function is meant for internal use and will be used whenever the `width` field +is required but unset. If you decide to modify the other internal fields of `g`, +it is probably a good idea to do `g.width[] = -1` so that this function gets +automatically called when needed, unless you are sure the width will not be +affected by your change. + +It is never necessary to use this function or to touch the `width` field of the graph +if you only modify the graph using the official API (i.e. if you never directly touch +the fields). +""" +function graph_width!(g::PeriodicGraph{N}) where N + distances = floyd_warshall_shortest_paths(cellgraph(g)).dists + extremalpoints = NTuple{N,NTuple{2,Vector{Tuple{Int,Int}}}}((([],[]),([],[]),([],[]))) + # a, x ∈ extremalpoints[i][j] where i ∈ ⟦1,N⟧ and j ∈ ⟦1,2⟧ means that + # vertex x has a neighbor whose offset is a*(-1)^(j-1) along dimension i + maxa = 1 + for e in edges(g) + iszero(ofs(e)) && continue + offset = ofs(e) + for i in 1:N + iszero(offset[i]) && continue + j = signbit(offset[i]) + 1 + a = abs(offset[i]) + if a > maxa + maxa = a + end + push!(extremalpoints[i][j], (a, src(e))) + push!(extremalpoints[i][3-j], (a, dst(e))) + end + end + + width::Rational{Int} = Rational(nv(g)+1) + for i in 1:N + if any(isempty.(extremalpoints[i])) + # TODO check this part + if width > 1 + width = 1//1 + end + continue + end + for (a1, x1) in extremalpoints[i][1], (a2, x2) in extremalpoints[i][2] + dist = distances[x1, x2] + if dist == typemax(Int) + dist = 1 + end + d = (dist + 1) // (a1 + a2) + if d < width + width = d + end + end + end + g.width[] = width == nv(g)+1 ? Rational(maxa) : width +end + +function LightGraphs._neighborhood(g::Union{PeriodicGraph{2}, PeriodicGraph{3}}, v::Integer, d::Real, distmx::AbstractMatrix{U}, neighborfn::Function) where U <: Real + @assert typeof(neighborfn) === typeof(outneighbors) + N = ndims(g) + Q = Tuple{PeriodicVertex{N}, U}[] + d < zero(U) && return Q + start_vertex = PeriodicVertex{N}(v) + push!(Q, (start_vertex, zero(U),) ) + n = nv(g) + width = g.width[] + if width == -1 + width = graph_width!(g) + end + seen_size = n*(2*(1 + fld(d-1, width)) + 1)^N + seen = falses(seen_size) + seen[hash_position(start_vertex, n)] = true + #=@inbounds=# for (src, currdist) in Q + currdist == d && continue # should be in Q but all its neighbours are too far + for dst in outneighbors(g, src.v) + dst = PeriodicVertex{N}(dst.v, dst.ofs .+ src.ofs) + position = hash_position(dst, n) + if !seen[position] + seen[position] = true + distance = currdist + distmx[src.v, dst.v] + if distance <= d + push!(Q, (dst, distance)) + end + end + end + end + return Q +end + +function LightGraphs._neighborhood(g::PeriodicGraph{N}, v::Integer, d::Real, distmx::AbstractMatrix{U}, neighborfn::Function) where {N,U <: Real} + @assert typeof(neighborfn) === typeof(outneighbors) + Q = Tuple{PeriodicVertex, U}[] + d < zero(U) && return Q + start_vertex = PeriodicVertex{N}(v) + push!(Q, (start_vertex, zero(U),) ) + n = nv(g) + seen = Set{PeriodicVertex{N}}() + push!(seen, start_vertex) + #=@inbounds=# for (src, currdist) in Q + currdist == d && continue # should be in Q but all its neighbours are too far + @simd for dst in outneighbors(g, src.v) + dst = PeriodicVertex(dst.v, dst.ofs .+ src.ofs) + if dst ∉ seen + push!(seen, dst) + distance = currdist + distmx[src.v, dst.v] + if distance <= d + push!(Q, (dst, distance)) + end + end + end + end + return Q +end + +""" + coordination_sequence(g::PeriodicGraph, v::Integer, dmax) + +Compute the list of numbers of `n`-th neighbors of vertex `v` in graph `g`, for +`1 ≤ n ≤ dmax`. +""" +function coordination_sequence(g::PeriodicGraph, v::Integer, dmax) + Q = LightGraphs._neighborhood(g, v, dmax, weights(g), outneighbors) + popfirst!(Q) + ret = zeros(Int, dmax) + for (_, d) in Q + ret[d] += 1 + end + return ret +end + +""" + cellgraph(g::PeriodicGraph) + +Extract a simple graph from `g` by only keeping the edges that are strictly +within the initial cell. + +See also `periodiccellgraph` to keep these edges. +""" +function cellgraph(g::PeriodicGraph) + edgs = [Edge{Int}(src(x), dst(x)) for x in edges(g) if iszero(ofs(x))] + ret = SimpleGraph(edgs) + add_vertices!(ret, nv(g) - nv(ret)) + return ret +end + +""" + periodiccellgraph(g::PeriodicGraph) + +Extract a simple graph from `g` by removing all indications of offset in the edges. +This means that edges that used to cross the boundaries of the initial cell now +bind the source vertex to the representative of the destination vertex that is in +the initial cell. + +Note that these modified edges may turn into loops. + +See also `cellgraph` to remove these edges. +""" +function periodiccellgraph(g::PeriodicGraph) + ret = SimpleGraph([Edge{Int}(i, j.v) for i in vertices(g) for j in outneighbors(g, i)]) + add_vertices!(ret, nv(g) - nv(ret)) + return ret +end + end diff --git a/test/runtests.jl b/test/runtests.jl index 93f51ea..243c6a7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,343 @@ -using PeriodicGraphs using Test +using Random +using LightGraphs +using PeriodicGraphs +import PeriodicGraphs: ofs, vertex_permutation, LoopException +import StaticArrays: SVector + +# using Aqua +# Aqua.test_all(PeriodicGraphs) + +@testset "PeriodicVertex" begin + @test_throws MethodError PeriodicVertex2D() # Always require a vertex identifier + @test iszero(ofs(PeriodicVertex2D(1))) + @test ndims(PeriodicVertex2D(1)) == 2 + @test PeriodicVertex2D(1) != PeriodicVertex3D(1) + @test PeriodicVertex{4}(5) == PeriodicVertex{4}(5) + @test PeriodicVertex{1}(1, SVector{1,Int}(3)) == PeriodicVertex{1}(1, [3]) + @test PeriodicVertex(3, (1,2)) == PeriodicVertex2D(3, (1,2)) == PeriodicVertex(3, SVector{2,Int}(1,2)) + @test PeriodicVertex[(1, (-1,0))] == [PeriodicVertex2D(1, (-1,0))] == PeriodicVertex2D[(1, (-1,0))] + @test_throws MethodError PeriodicVertex(2, [4]) # Performance trap if the dimension is unknown + @test_throws DimensionMismatch PeriodicVertex2D(1, (2,)) + @test_throws DimensionMismatch PeriodicVertex3D(3, (3, 1, 4, 0)) + @test ndims(PeriodicVertex(0, SVector{2,Int}(1,4))) == 2 # Inference from SVector + @test ndims(PeriodicVertex{0}(1)) == 0 # No need of a second argument here + @test PeriodicVertex{0}(1) < PeriodicVertex{0}(2) + @test PeriodicVertex(2, (3,)) > PeriodicVertex(2, (1,)) + @static if VERSION < v"1.6.0-DEV" + @test string(PeriodicVertex3D(12, (1,0,-1))) == "PeriodicVertex{3}(12, (1,0,-1))" + @test string(PeriodicVertex2D[(1, (0,0)), (1, (-1,0))]) == "PeriodicVertex{2}[(1, (0,0)), (1, (-1,0))]" + else + @test string(PeriodicVertex3D(12, (1,0,-1))) == "PeriodicVertex3D(12, (1,0,-1))" + @test string(PeriodicVertex2D[(1, (0,0)), (1, (-1,0))]) == "PeriodicVertex2D[(1, (0,0)), (1, (-1,0))]" + end + @test ofs(PeriodicVertex(1, (2,1))) == SVector{2,Int}(2,1) +end + +@testset "PeriodicEdge" begin + @test_throws LoopException PeriodicEdge2D(1, 1, (0, 0)) + @test_throws LoopException PeriodicEdge(2, PeriodicVertex3D(2)) + @test_throws LoopException PeriodicEdge{0}(3, 3, Int[]) + @test_throws LoopException PeriodicEdge(1, 1, ()) + # There should be no default zero offset to prevent programmers from forgetting the offset + @test_throws MethodError PeriodicEdge{1}(1, 2) + @test PeriodicEdge2D(1, 2, (0,0)) != PeriodicEdge3D(1, 2, (0,0,0)) + @test PeriodicEdge(2, PeriodicVertex2D(1, (1,3))) == PeriodicEdge2D(2, 1, (1,3)) + @test PeriodicEdge(1, 1, SVector{1,Int}(2)) == PeriodicEdge{1}(1, 1, SVector{1,Int}(2)) + @test PeriodicEdge(2, 1, ()) == PeriodicEdge{0}((2, 1, SVector{0,Int}())) + @test PeriodicEdge((1, 2, (1,0,0))) == PeriodicEdge3D((1, 2, (1,0,0))) + @test PeriodicEdge{1}[(1,2,SVector{1,Int}(3))] == PeriodicEdge[(1,2,(3,))] == [PeriodicEdge(1, 2, (3,))] + @static if VERSION < v"1.6.0-DEV" + @test string(PeriodicEdge3D[(1, 2, (1,0,0))]) == "PeriodicEdge{3}[(1, 2, (1,0,0))]" + @test string(PeriodicEdge(3,4, (0,0))) == "PeriodicEdge{2}(3, 4, (0,0))" + else + @test string(PeriodicEdge3D[(1, 2, (1,0,0))]) == "PeriodicEdge3D[(1, 2, (1,0,0))]" + @test string(PeriodicEdge(3,4, (0,0))) == "PeriodicEdge2D(3, 4, (0,0))" + end + @test reverse(reverse(PeriodicEdge(1, 1, (2,3,1)))) == PeriodicEdge3D(1, 1, (2,3,1)) + @test reverse(PeriodicEdge(1, 2, (1,-2))) == PeriodicEdge(2, 1, (-1,2)) + @test src(PeriodicEdge(1, 2, (-1, 0))) == 1 + @test dst(PeriodicEdge(2, 2, (1,))) == 2 + @test ofs(PeriodicEdge(1, 1, (0,1,3))) == SVector{3,Int}(0,1,3) + @test PeriodicEdge2D(1, 1, (2,3)) < PeriodicEdge2D(1, 1, (3,0)) < PeriodicEdge2D(1, 2, (3,0)) < PeriodicEdge2D(2, 2, (3,0)) + @test ndims(PeriodicEdge(1, 2, ())) == 0ndims(PeriodicEdge(1, 2, ())) == 0 + @test ndims(PeriodicEdge(1, 1, (2,0,0,0))) == 4 +end + +@testset "Empty graphs" begin + @test_throws MethodError PeriodicGraph() # Missing dimension + g = PeriodicGraph3D() + @test ne(g) == nv(g) == 0 + @test ndims(g) == 3 + @test !rem_edge!(g, PeriodicEdge3D(1, 2, (0,0,0))) + @test_throws ArgumentError rem_vertex!(g, 0) + @test_throws ArgumentError rem_vertex!(g, 1) + @test isempty(rem_vertices!(g, Int[])) + @test isempty(edges(g)) + @test isempty(vertices(g)) + gg = vertex_permutation(g, []) + @test ne(gg) == nv(gg) == 0 && isempty(edges(gg)) + @test g == gg + gg, vlist = induced_subgraph(g, Int[]) + @test isempty(vlist) + @test ne(gg) == nv(gg) == 0 && isempty(edges(gg)) + @test g == gg + gg = offset_representatives!(g, Int[]) + @test ne(gg) == nv(gg) == 0 && isempty(edges(gg)) + @test ne(g) == nv(g) == 0 && isempty(edges(g)) + @test g == gg + gg = swap_axes!(g, Int[1,2,3]) + @test ne(gg) == nv(gg) == 0 && isempty(edges(gg)) + @test ne(g) == nv(g) == 0 && isempty(edges(g)) + @test g == gg + @test isempty(connected_components(g)) + @test edgetype(g) == PeriodicEdge3D + @test eltype(g) == PeriodicVertex3D + @test edges(g) isa PeriodicGraphs.PeriodicEdgeIter{3} + @test_throws BoundsError coordination_sequence(g, 1, 0) + @test_throws BoundsError coordination_sequence(g, 1, 0) # once width is set +end + +@testset "Basic graph construction and modification" begin + @test_throws MethodError PeriodicGraph(2) # missing dimension + g = PeriodicGraph3D() + + @testset "add_vertex!, rem_vertex!" begin + @test add_vertex!(g) + @test ne(g) == 0 && nv(g) == 1 + @test g == PeriodicGraph3D(1) + @test add_vertex!(g) + @test ne(g) == 0 && nv(g) == 2 + @test g == PeriodicGraph3D(2) + @test_throws ArgumentError rem_vertex!(g, 0) + @test_throws ArgumentError rem_vertex!(g, 3) + @test rem_vertex!(g, 1) + @test ne(g) == 0 && nv(g) == 1 + @test g == PeriodicGraph3D(1) + @test rem_vertex!(g, 1) + @test ne(g) == 0 && nv(g) == 0 + @test g == PeriodicGraph3D(0) + end + + @testset "add_vertices!, rem_vertices!" begin + @test isempty(edges(g)) + @test isempty(vertices(g)) + @test_throws ArgumentError rem_vertices!(g, Int[1]) + @test isempty(rem_vertices!(g, Int[])) + @test add_vertices!(g, 0) == 0 + @test nv(g) == ne(g) == 0 && g == PeriodicGraph3D() + @test add_vertices!(g, 1) == 1 + @test ne(g) == 0 && nv(g) == 1 && g == PeriodicGraph3D(1) + @test only(vertices(g)) == 1 + @test rem_vertices!(g, [1]) == Int[] + @test add_vertices!(g, 5) == 5 + @test ne(g) == 0 && nv(g) == 5 && g == PeriodicGraph3D(5) + @test vertices(g) == 1:5 + @test issetequal(rem_vertices!(g, [2,3]), [1,4,5]) + @test ne(g) == 0 && nv(g) == 3 && g == PeriodicGraph3D(3) + @test vertices(g) == 1:3 + @test rem_vertex!(g, 2) + @test ne(g) == 0 && nv(g) == 2 && g == PeriodicGraph3D(2) + @test vertices(g) == [1,2] + @test rem_vertices!(g, [1,2]) == Int[] + end + + @testset "Single edge" begin + @test isempty(edges(g)) + @test isempty(vertices(g)) + @test !rem_edge!(g, PeriodicEdge(1, 2, (0,0,0))) + @test g == PeriodicGraph3D() + @test !add_edge!(g, PeriodicEdge(1, 1, (-1,0,0))) + @test g == PeriodicGraph3D() + @test !has_edge(g, 1, 1) + @test !has_edge(g, PeriodicEdge3D(1, 2, (0,0,0))) + @test add_vertex!(g) + @test !has_edge(g, 1, 1) + @test !has_edge(g, PeriodicEdge3D(1, 1, (0,1,1))) + @test add_edge!(g, PeriodicEdge(1, 1, (-1,0,0))) + @test !add_edge!(g, PeriodicEdge(1, 1, (-1,0,0))) + @test has_edge(g, 1, 1) + @test has_edge(g, PeriodicEdge3D(1, 1, (1,0,0))) + @test has_edge(g, PeriodicEdge3D(1, 1, (-1,0,0))) + @test g == PeriodicGraph3D(PeriodicEdge3D[(1, 1, (1,0,0))]) + @test g != PeriodicGraph3D(1) + @test rem_edge!(g, PeriodicEdge(1, 1, (-1,0,0))) + @test g == PeriodicGraph3D(1) + @test !has_edge(g, 1, 1) + @test !has_edge(g, PeriodicEdge3D(1, 1, (1,0,0))) + @test !has_edge(g, PeriodicEdge3D(1, 1, (-1,0,0))) + @test add_edge!(g, PeriodicEdge(1, 1, (2,-1,0))) + @test !add_edge!(g, PeriodicEdge(1, 1, (-2,1,0))) + @test g == PeriodicGraph3D(PeriodicEdge3D[(1, 1, (-2,1,0))]) + @test has_edge(g, 1, 1) + @test has_edge(g, PeriodicEdge3D(1, 1, (2,-1,0))) + @test !has_edge(g, PeriodicEdge3D(1, 1, (1,0,0))) + @test rem_edge!(g, PeriodicEdge(1, 1, (2,-1,0))) + @test g == PeriodicGraph3D(1) + @test add_edge!(g, PeriodicEdge(1, 1, (1,0,0))) + @test rem_vertex!(g, 1) + @test g == PeriodicGraph3D(0) + @test add_vertices!(g, 2) == 2 + @test !has_edge(g, 1, 2) + @test g == PeriodicGraph3D(2) + @test add_edge!(g, PeriodicEdge(1, 2, (0,0,0))) + @test has_edge(g, 1, 2) + @test has_edge(g, PeriodicEdge(1, 2, (0,0,0))) + @test !has_edge(g, PeriodicEdge(1, 2, (1,0,0))) + @test g == PeriodicGraph3D(PeriodicEdge3D[(1, 2, (0,0,0))]) + @test g != PeriodicGraph3D(2) + @test rem_edge!(g, PeriodicEdge(1, 2, (0,0,0))) + @test !has_edge(g, 1, 2) + @test g == PeriodicGraph3D(2) + @test add_edge!(g, PeriodicEdge(1, 2, (-1,0,0))) + @test g == PeriodicGraph3D(PeriodicEdge3D[(1, 2, (-1,0,0))]) + @test g == PeriodicGraph3D(PeriodicEdge3D[(2, 1, (1,0,0))]) + @test g != PeriodicGraph3D(PeriodicEdge3D[(1, 2, (1,0,0))]) + @test has_edge(g, 1, 2) + @test has_edge(g, PeriodicEdge(1, 2, (-1,0,0))) + @test has_edge(g, PeriodicEdge(2, 1, (1,0,0))) + @test rem_vertex!(g, 1) + @test g == PeriodicGraph3D(1) + @test !has_edge(g, 1, 2) + @test !has_edge(g, 1, 1) + end +end + +@testset "String graph construction" begin + @test PeriodicGraph("0 1 2 1 3") == PeriodicGraph(PeriodicEdge{0}[(1, 2, ()), (1, 3, ())]) + @test PeriodicGraph("3") == PeriodicGraph3D("3") == PeriodicGraph3D() + @test ne(PeriodicGraph("3 1 1 1 0 0 1 1 -1 0 0 1 1 3 0 0")) == 2 + @test nv(PeriodicGraph("2 11 12 0 0")) == 12 + @test_throws ArgumentError PeriodicGraph2D("") + @test_throws ArgumentError PeriodicGraph("") + @test_throws DimensionMismatch PeriodicGraph2D("4 1 2 0 0 0 1") + @test_throws DimensionMismatch PeriodicGraph3D("2 1 2 0 0 0 1 3 0 1 0") + @test_throws ArgumentError PeriodicGraph("1 2 1") + @test_throws ArgumentError PeriodicGraph{4}("4 1") + @test_throws ArgumentError PeriodicGraph("2 1 1 0 1 0") + @test_throws LoopException PeriodicGraph("3 1 1 0 0 0") + @test string(PeriodicGraph("2 1 2 0 0 2 3 1 0")) == "2 1 2 0 0 2 3 1 0" + @test string(PeriodicGraph("1 23 4 5 10 9 43")) == "1 4 23 -5 9 10 -43" + @test string(PeriodicGraph{5}()) == "5" +end + +@testset "Complex graphs modifications" begin + edgs = PeriodicEdge3D[(1, 1, (1,0,0)), (1, 2, (0,0,0)), (1, 3, (0,0,1)), + (1, 3, (0,1,-1)), (2, 3, (0,0,0)), (2, 4, (0,-1,0))] + g = PeriodicGraph3D(edgs) + @test nv(g) == 4 + @test ne(g) == length(edges(g)) == length(edgs) + @test collect(edges(g)) == edgs + insert!(edgs, 5, PeriodicEdge(1, 4, (-1,-1,0))) + @test add_edge!(g, edgs[5]) + @test ne(g) == length(edges(g)) == length(edgs) + @test collect(edges(g)) == edgs + @test rem_edge!(g, edgs[4]) + deleteat!(edgs, 4) + @test ne(g) == length(edges(g)) == length(edgs) + @test collect(edges(g)) == edgs +end + +@testset "Neighbors" begin + g = PeriodicGraph3D(1) + @test add_vertices!(g, 2) == 2 + @test add_edge!(g, PeriodicEdge3D(2, 3, (0, 0, 1))) + @test only(neighbors(g, 2)) == PeriodicVertex3D(3, (0, 0, 1)) + @test only(neighbors(g, 3)) == PeriodicVertex3D(2, (0, 0, -1)) + @test !has_edge(g, 2, 2) + @test add_edge!(g, PeriodicEdge3D(2, 2, (1, 0, 0))) + @test only(neighbors(g, 3)) == PeriodicVertex3D(2, (0, 0, -1)) + @test length(neighbors(g, 2)) == 3 + @test has_edge(g, 2, 2) + @test neighborhood(g, 1, 1) == pushfirst!(collect(neighbors(g, 1)), PeriodicVertex3D(1, (0, 0, 0))) +end + +@testset "Graph reduction" begin + g = PeriodicGraph3D("3 1 2 0 0 0 1 3 0 1 0 2 3 1 0 -1 + 4 4 0 0 1 4 5 0 -1 0 5 6 0 0 0") + c = cellgraph(g) + p = periodiccellgraph(g) + @test nv(g) == nv(p) == nv(c) + @test issubset(edges(c), edges(p)) + for e in edges(p) + @test has_edge(g, src(e), dst(e)) + end + @test connected_components(p) == connected_components(g) +end + +@testset "Periodic vertex hash" begin + n::Int = 4 + seen = falses((2n-1)^3) + for i in -n+1:n-1, j in -n+1:n-1, k in -n+1:n-1 + x = PeriodicGraphs.hash_position(SVector{3,Int}(i,j,k))+1 + @test !seen[x] + seen[x] = true + end + @test all(seen) + append!(seen, falses((2n+1)^3-(2n-1)^3)) + for i in (-n,n), j in -n:n, k in -n:n + x = PeriodicGraphs.hash_position(SVector{3,Int}(i,j,k))+1 + @test !seen[x] + seen[x] = true + end + for i in -n+1:n-1, j in (-n,n), k in -n:n + x = PeriodicGraphs.hash_position(SVector{3,Int}(i,j,k))+1 + @test !seen[x] + seen[x] = true + end + for i in -n+1:n-1, j in -n+1:n-1, k in (-n,n) + x = PeriodicGraphs.hash_position(SVector{3,Int}(i,j,k))+1 + @test !seen[x] + seen[x] = true + end + + for i in 1:length(seen) + if !seen[i] + @show i + end + end + @test all(seen) + g = PeriodicGraph3D([PeriodicEdge3D(1, 1, (0, 0, 1)), + PeriodicEdge3D(1, 1, (1, 1, 0))]) + @test neighborhood(g, 1, 1) == PeriodicVertex3D[(1, (0, 0, 0)), + (1, (-1, -1, 0)), (1, (0, 0, -1)), (1, (0, 0, 1)), (1, (1, 1, 0))] + @test neighborhood(g, 1, 2) == PeriodicVertex3D[(1, (0, 0, 0)), (1, (-1, -1, 0)), + (1, (0, 0, -1)), (1, (0, 0, 1)), (1, (1, 1, 0)), (1, (-2, -2, 0)), + (1, (-1, -1, -1)), (1, (-1, -1, 1)), (1, (0, 0, -2)), (1, (1, 1, -1)), + (1, (0, 0, 2)), (1, (1, 1, 1)), (1, (2, 2, 0))] +end + +@testset "Edge iteration" begin + g::PeriodicGraph3D = PeriodicGraph3D([PeriodicEdge3D(2, 3, (-1, 0, 0)), + PeriodicEdge3D(2, 3, (0, 1, 0)), + PeriodicEdge3D(4, 4, (-1, 1, 0)), + PeriodicEdge3D(1, 3, (0, 0, 1)), + PeriodicEdge3D(3, 2, (0, 0, 0))]) + expected = PeriodicEdge3D[(1, 3, (0, 0, 1)), (2, 3, (-1, 0, 0)), + (2, 3, (0, 0, 0)), (2, 3, (0, 1, 0)), + (4, 4, (1, -1, 0))] + @test collect(edges(g)) == expected + @test PeriodicEdge3D(4, 4, (1, -1, 0)) ∈ edges(g) + add_vertex!(g) + @test collect(edges(g)) == expected + @test length(edges(g)) == ne(g) +end -@testset "PeriodicGraphs.jl" begin - # Write your tests here. +@testset "[Legacy test] Vertex removal" begin + g::PeriodicGraph3D = PeriodicGraph3D([PeriodicEdge3D(1, 3, (0, 0, 1)), + PeriodicEdge3D(2, 3, (0, 1, 0)), + PeriodicEdge3D(2, 3, (0, -1, 0)), + PeriodicEdge3D(3, 2, (0, 0, 0)), + PeriodicEdge3D(4, 4, (1, 1, 0))]) + gg = deepcopy(g) + @test rem_vertices!(g, Int[]) == [1, 2, 3, 4] + @test g == gg + @test rem_vertices!(g, [2, 4], false) == [1, 3] + @test nv(g) == 2 + @test ne(g) == 1 + @test neighbors(g, 1) == PeriodicVertex3D[(2, (0, 0, 1))] + @test neighbors(g, 2) == PeriodicVertex3D[(1, (0, 0, -1))] + @test rem_vertex!(gg, 2) + @test rem_vertex!(gg, 2) + @test g == gg end