diff --git a/Project.toml b/Project.toml index 53ecdb9..abbdebb 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PeriodicArrays" uuid = "343d6138-6384-4525-8bee-38906309ab36" authors = ["Andreas Feuerpfeil "] -version = "1.0.0" +version = "1.1.0" [compat] julia = "1.10" diff --git a/README.md b/README.md index e4d2fa0..9fa2e59 100644 --- a/README.md +++ b/README.md @@ -48,10 +48,10 @@ This package is compatible with [`OffsetArrays.jl`](https://github.com/JuliaArra ## Installation -The package is not yet registered in the Julia general registry. It can be installed trough the package manager with the following command: +The package is registered in the Julia general registry. It can be installed trough the package manager with the following command: ```julia-repl -pkg> add git@github.com:ManyBodyLab/PeriodicArrays.jl.git +pkg> add PeriodicArrays ``` ## Code Samples diff --git a/docs/files/README.jl b/docs/files/README.jl index f197aa0..7ab4808 100644 --- a/docs/files/README.jl +++ b/docs/files/README.jl @@ -19,10 +19,10 @@ # ## Installation -# The package is not yet registered in the Julia general registry. It can be installed trough the package manager with the following command: +# The package is registered in the Julia general registry. It can be installed trough the package manager with the following command: # ```julia-repl -# pkg> add git@github.com:ManyBodyLab/PeriodicArrays.jl.git +# pkg> add PeriodicArrays # ``` # ## Code Samples diff --git a/src/PeriodicArrays.jl b/src/PeriodicArrays.jl index f976908..14891d0 100644 --- a/src/PeriodicArrays.jl +++ b/src/PeriodicArrays.jl @@ -229,4 +229,75 @@ function Base.insert!(a::PeriodicVector, i::Integer, item) return a end +function Base.repeat(A::PeriodicArray{T, N}; inner = nothing, outer = nothing) where {T, N} + map = A.map + # If no outer repetition is requested, just repeat the parent array as usual + A_new = repeat(parent(A); inner = inner) + + if !isnothing(outer) + # allow passing a single integer or a tuple/ntuple for per-dimension repeats + if isa(outer, Number) + outer = ntuple(i -> Int(outer), N) + else + outer = ntuple(i -> Int(outer[i]), N) + end + + # if `inner` was provided, A_new already contains the repeated parent + base = A_new + axs = axes(base) + ps = size(base) + newsize = ntuple(i -> ps[i] * outer[i], N) + + # create a tiled parent filled with translated values from `map` + A_tiled = similar(base, newsize) + tile_ranges = ntuple(i -> 0:(outer[i] - 1), N) + for tile in CartesianIndices(tile_ranges) + shifts = Tuple(Int(tile[i]) for i in 1:N) + for pos in CartesianIndices(base) + tgt = ntuple(i -> tile[i] * ps[i] + (pos[i] - firstindex(axs[i]) + 1), N) + @inbounds A_tiled[tgt...] = map(base[pos], shifts...) + end + end + + @inline function map_new(x::T, shift::Vararg{Int, N}) + # shifts passed to this map refer to super-cell shifts; amplify + # by `outer` to convert them to original unit-cell shifts. + amplified = ntuple(i -> shift[i] * outer[i], N) + return map(x, amplified...) + end + + return PeriodicArray(A_tiled, map_new) + end + + return PeriodicArray(A_new, map) +end + +function Base.reverse(arr::PeriodicArray{T, N, A, F}; dims = :) where {T, N, A, F} + dims == Colon() && return _reverse(arr) + return _reverse(arr, dims) +end + +function _reverse(arr::PeriodicArray{T, N, A, F}) where {T, N, A, F} + base = reverse(parent(arr)) + + @inline function map_rev(x::T, shifts::Vararg{Int, N}) + neg = ntuple(i -> -shifts[i], N) + return arr.map(x, neg...) + end + + return PeriodicArray(base, map_rev) +end + +function _reverse(arr::PeriodicArray{T, N, A, F}, dims...) where {T, N, A, F} + base = reverse(parent(arr); dims = dims) + dimsset = Set(dims) + + @inline function map_rev(x::T, shifts::Vararg{Int, N}) + adj = ntuple(i -> (i in dimsset) ? -shifts[i] : shifts[i], N) + return arr.map(x, adj...) + end + + return PeriodicArray(base, map_rev) +end + end diff --git a/test/test_nontrivial_boundary.jl b/test/test_nontrivial_boundary.jl index a49918e..c0b5bc7 100644 --- a/test/test_nontrivial_boundary.jl +++ b/test/test_nontrivial_boundary.jl @@ -325,4 +325,78 @@ for f in translation_functions # TODO: Figure out how to fix indexing for non-trivial f #@test a[a .> 4] == 5:9 end + + @testset "repeat 1D" begin + a = PeriodicArray([1, 2, 3], f) + # outer as scalar + ar = repeat(a; outer = 2) + base = parent(a) + # expected tiled parent: tiles over shifts 0,1 + expected = vcat([f(base[i], 0) for i in eachindex(base)]..., [f(base[i], 1) for i in eachindex(base)]...) + @test parent(ar) == expected + @test size(parent(ar)) == (length(base) * 2,) + + val = 5 + @test ar.map(val, 1) == a.map(val, 2) + @test ar.map(val, 3) == a.map(val, 6) + + # inner repetition + ai = repeat(a; inner = 2) + @test parent(ai) == repeat(parent(a); inner = 2) + @test size(parent(ai)) == (length(base) * 2,) + + # combined inner+outer + aio = repeat(a; inner = 2, outer = 3) + @test size(parent(aio)) == (length(base) * 2 * 3,) + + # outer as tuple (1D tuple) + ar2 = repeat(a; outer = (2,)) + @test parent(ar2) == expected + end + + @testset "repeat 2D" begin + b = PeriodicArray(reshape(1:6, 3, 2), f) + o1, o2 = 2, 3 + br = repeat(b; outer = (o1, o2)) + base = parent(b) + expected2 = similar(base, size(base, 1) * o1, size(base, 2) * o2) + for t2 in 0:(o2 - 1), t1 in 0:(o1 - 1) + for pos in CartesianIndices(base) + tgt = ( + t1 * size(base, 1) + (pos[1] - firstindex(axes(base, 1)) + 1), + t2 * size(base, 2) + (pos[2] - firstindex(axes(base, 2)) + 1), + ) + @inbounds expected2[tgt...] = f(base[pos], t1, t2) + end + end + @test parent(br) == expected2 + @test size(parent(br)) == (size(base, 1) * o1, size(base, 2) * o2) + + @test br.map(1, 1, 2) == b.map(1, 2, 6) + @test br.map(7, 0, 1) == b.map(7, 0, 3) + + # inner repetition in 2D + bi = repeat(b; inner = (2, 1)) + @test parent(bi) == repeat(parent(b); inner = (2, 1)) + @test size(parent(bi)) == (size(base, 1) * 2, size(base, 2) * 1) + end + + @testset "reverse" begin + a = PeriodicArray([1, 2, 3], f) + ra = reverse(a) + @test parent(ra) == reverse(parent(a)) + @test ra[1] == a[3] + @test ra[2] == a[2] + + # 2D reverse across first dimension + b = PeriodicArray(reshape(1:6, 3, 2), f) + rb1 = reverse(b; dims = 1) + @test parent(rb1) == reverse(parent(b); dims = 1) + @test all(rb1[i, j] == b[4 - i, j] for i in -10:10, j in -10:10) + + # full reverse + rb = reverse(b) + @test parent(rb) == reverse(parent(b)) + @test all(rb[i, j] == b[4 - i, 3 - j] for i in -10:10, j in -10:10) + end end