Add make_supercell

docs/src/utilities.md

@@ -26,12 +26,14 @@ far from the reference.
 Several functions transform a periodic graph into another isomorphic to the input, by
 renumbering the vertices ([`vertex_permutation`](@ref PeriodicGraphs.vertex_permutation))
 or the axes ([`swap_axes!`](@ref)), or by offsetting the chosen representatives for each
-vertex ([`offset_representatives!`](@ref)).
+vertex ([`offset_representatives!`](@ref)). It is also possible to make an isomorphic graph
+with more vertices per unit cell by using a supercell ([`make_supercell`](@ref)).
 
 ```@docs
 PeriodicGraphs.vertex_permutation
 swap_axes!
 offset_representatives!
+make_supercell
 ```
 
 ## Dimension reduction

src/algorithms/other.jl

@@ -2,6 +2,7 @@
 
 export offset_representatives!,
        swap_axes!,
+       make_supercell,
        truncated_graph,
        quotient_graph,
        slice_graph
@@ -64,6 +65,64 @@ function swap_axes!(g::PeriodicGraph{N}, t) where N
 end
 @noinline __throw_invalid_axesswap() = throw(DimensionMismatch("The number of axes must match the dimension of the graph"))
 
+struct MetaClock{T}
+    maxs::T
+    factors::Vector{Int}
+    len::Int
+end
+function MetaClock(t)
+    t isa AbstractVector || (t = collect(t))
+    factors = Vector{Int}(undef, length(t))
+    cumprod!(factors, t)
+    len = pop!(factors)
+    pushfirst!(factors, 1)
+    return MetaClock(t, factors, len)
+end
+Base.length(x::MetaClock) = x.len
+Base.eltype(::Type{MetaClock{T}}) where {T} = Vector{Int}
+function Base.iterate(x::MetaClock, state::Vector{Int}=[-1; zeros(Int, length(x.maxs)-1)])
+    maxs = x.maxs
+    for i in 1:length(maxs)
+        y = state[i] + 1
+        if y == maxs[i]
+            state[i] = 0
+        else
+            state[i] = y
+            return (state, state)
+        end
+    end
+    return nothing
+end
+
+"""
+    make_supercell(g::PeriodicGraph, t)
+
+Return a graph isomorphic to the input `g` whose its unit cell is a repetition of that of
+`g`, each dimension `i` being repeated `t[i]` times.
+It follows that the number of vertices of `make_supercell(g, t)` is `prod(t)*nv(g)`
+
+`t` must be an interator over positive integers.
+"""
+function make_supercell(g::PeriodicGraph{N}, t::S) where {N,S}
+    length(t) == N || __throw_invalid_axesswap()
+    N == 0 && return g
+    newedges = PeriodicEdge{N}[]
+    __check_nonpositive_axe(minimum(t))
+    n = nv(g)
+    clock = MetaClock(t)
+    for e in edges(g)
+        (src, (dst, ofs)) = e
+        for pos in clock
+            factor = n*sum(prod, zip(clock.factors, pos); init=0)
+            newofs, newpos = eachrow(reinterpret(reshape, Int, fldmod.(ofs .+ pos, t)))
+            newfactor = n*sum(prod, zip(clock.factors, newpos); init=0)
+            push!(newedges, PeriodicEdge{N}(factor+src, newfactor+dst, newofs))
+        end
+    end
+    return PeriodicGraph{N}(length(clock)*n, newedges)
+end
+__check_nonpositive_axe(i) = i > 0 || throw(DomainError(i, "All supercell dimensions must be strictly positive"))
+
 """
     truncated_graph(g::PeriodicGraph)
 

src/precompile.jl

@@ -176,8 +176,12 @@ function _precompile_()
         # Other utils
         @enforce Base.precompile(Tuple{typeof(offset_representatives!),PeriodicGraph{i},Vector{Int}})
         @enforce Base.precompile(Tuple{typeof(offset_representatives!),PeriodicGraph{i},Vector{NTuple{i,Int}}})
-        @enforce Base.precompile(Tuple{typeof(swap_axes!),PeriodicGraph{i},Vector{Int}})
         @enforce Base.precompile(Tuple{typeof(swap_axes!),PeriodicGraph{i},NTuple{i,Int}})
+        @enforce Base.precompile(Tuple{typeof(swap_axes!),PeriodicGraph{i},SVector{i,Int}})
+        @enforce Base.precompile(Tuple{typeof(swap_axes!),PeriodicGraph{i},Vector{Int}})
+        @enforce Base.precompile(Tuple{typeof(make_supercell),PeriodicGraph{i},NTuple{i,Int}})
+        @enforce Base.precompile(Tuple{typeof(make_supercell),PeriodicGraph{i},SVector{i,Int}})
+        @enforce Base.precompile(Tuple{typeof(make_supercell),PeriodicGraph{i},Vector{Int}})
         @enforce Base.precompile(Tuple{typeof(truncated_graph),PeriodicGraph{i}})
         @enforce Base.precompile(Tuple{typeof(quotient_graph),PeriodicGraph{i}})
         for N in 1:i

test/runtests.jl

@@ -158,6 +158,9 @@ end
     @test ne(gg) == nv(gg) == 0 && isempty(edges(gg))
     @test ne(g) == nv(g) == 0 && isempty(edges(g))
     @test g == gg
+    @test_throws DimensionMismatch make_supercell(g, (2, 1, 1, 1))
+    @test_throws DomainError make_supercell(g, (2, 0, 1))
+    @test g == make_supercell(g, (2, 1, 3))
     @test isempty(connected_components(g))
     @test edgetype(g) == PeriodicEdge3D
     @test eltype(g) == PeriodicVertex3D
@@ -439,6 +442,19 @@ end
     @test edges(gg) < edges(gcopy) # lexicographical ordering on the list of edges
 end
 
+@testset "Supercell" begin
+    function check_make_supercell(g, t)
+        x = make_supercell(g, t)
+        @test nv(x) == prod(t)*nv(g)
+        x
+    end
+    g = PeriodicGraph("3 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")
+    @test make_supercell(g, (1 for _ in 1:3)) == g
+    @test check_make_supercell(g, [3, 2, 4]) == check_make_supercell(check_make_supercell(check_make_supercell(g, [3, 1, 1]), [1, 2, 1]), [1, 1, 4])
+    g312 = check_make_supercell(g, SVector{3,Int}(3, 1, 2))
+    @test truncated_graph(g) == truncated_graph(g312[1:4]) == truncated_graph(g312[5:8]) == truncated_graph(g312[9:12])
+end
+
 const sqc7399 = PeriodicGraph3D("3 1 6 0 -1 -1 1 10 -1 -1 0 1 10 -1 0 0 2 5 0 0 0 2 7 -1 0 0 2 7 0 0 0 3 4 0 1 0 3 7 -1 1 0 3 7 0 0 -1 4 10 -1 0 1 4 10 0 -1 0 5 10 -1 0 0 5 10 0 0 0 6 7 0 0 0 6 7 0 1 0 7 7 0 0 1 7 7 0 1 0 7 7 1 -1 -1 7 7 1 0 0 7 8 0 0 0 7 8 0 0 1 7 10 -1 0 1 7 10 0 -1 0 7 10 0 0 0 7 10 0 0 1 8 9 0 -1 -1 9 10 0 0 0 9 10 0 0 1 10 10 0 0 1 10 10 0 1 0 10 10 1 -1 -1 10 10 1 0 0");
 const sqc10673 = PeriodicGraph3D("3 1 3 0 0 0 1 25 -1 -1 -1 1 26 -1 -1 -1 2 5 0 0 0 2 23 -1 0 0 2 24 -1 0 0 3 4 0 0 -1 3 6 0 -1 0 3 7 0 0 0 4 5 0 0 0 4 16 -1 0 0 4 19 0 0 1 5 6 0 0 0 5 8 0 0 0 6 16 0 1 0 6 19 -1 0 0 7 9 0 0 0 7 10 0 0 0 8 11 0 0 0 8 12 0 0 0 9 13 0 0 0 9 17 0 0 0 10 14 0 0 0 10 18 0 0 0 11 13 0 0 0 11 20 0 0 1 12 14 0 0 0 12 15 0 1 0 13 14 0 0 0 13 21 0 0 0 14 22 0 0 0 15 16 0 0 0 15 26 0 -1 0 16 17 0 0 0 17 23 0 0 0 18 19 0 0 0 18 24 0 0 0 19 20 0 0 0 20 25 0 0 -1 21 22 0 0 0 21 23 0 0 0 21 25 0 0 0 22 24 0 0 0 22 26 0 0 0")