When implementing a circular array, doesn’t it make much more sense to use 0-based indexing?
2 Likes
maybe, my toughts were to preserve the 1-based indexing of julia, and preserve the equality circular_x[i,j] == x[i,j]
3 Likes
There exists an implementation in ShiftedArrays: https://github.com/piever/ShiftedArrays.jl/blob/master/src/circshiftedarray.jl
If you need something more complex but along these lines, hopefully it can serve as a starting point to write your own custom array type.
1 Like
And there is the mod1 function to avoid all this shifting indices manually:
julia> mod1(5,5)
5
julia> mod1(6,5)
1
3 Likes
Thanks again longemen3000 (and others)! This do what I want except this:
julia> x1[:,0]
ERROR: BoundsError: attempt to access 4×4 CircularArray{Int64,2} at index [Base.Slice(Base.OneTo(4)), 0]
Stacktrace:
[1] throw_boundserror(::CircularArray{Int64,2}, ::Tuple{Base.Slice{Base.OneTo{Int64}},Int64}) at ./abstractarray.jl:484
[2] checkbounds at ./abstractarray.jl:449 [inlined]
[3] _getindex(::IndexCartesian, ::CircularArray{Int64,2}, ::Base.Slice{Base.OneTo{Int64}}, ::Int64) at ./multidimensional.jl:641
[4] getindex(::CircularArray{Int64,2}, ::Function, ::Int64) at ./abstractarray.jl:927
[5] top-level scope at none:0
x[1,0] and x[:,1] work so it is only a problem dealing with ranges outside array boundaries. Any ideas how to solve this?
1 Like
That’s just a matter of defining the appropriate checkbounds specialization. In fact, I think you can simply turn off bounds checks entirely:
Base.checkbounds(A::CircularArray, I...) = nothing
3 Likes
This is material for a small compact package!
2 Likes
Is there a package with this circular array type available and well-tested already?
2 Likes
@udistr can you comment on the difference between CyclicArrays.jl and CircularArrays.jl? Aren’t they equivalent?
CyclicArrays.jl is aimed at handling more complex multi-face array topologies. For example, a cubed-sphere grid with 6 faces, in which the faces are interconnected in a non-trivial manner.
1 Like
It is still not clear to me, perhaps a simple example of this sphere use case?
I don’t promise this is right, but it looks like the simplest things you can do are these:
julia> torus = CyclicArray(
[1 2 3; 4 5 6; 7 8 9],
hcat(
reshape([1 1 2 0; 1 1 1 0],1,1,2,4),
reshape([1 2 2 0; 1 2 1 0],1,1,2,4)
))
3×3 CyclicArray{Int64, 2}:
1 2 3
4 5 6
7 8 9
julia> [torus[i,j] for i in -2:6, j in -2:6]
9×9 Matrix{Int64}:
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
7 8 9 7 8 9 7 8 9
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
7 8 9 7 8 9 7 8 9
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
7 8 9 7 8 9 7 8 9
julia> klein = CyclicArray(
[1 2 3; 4 5 6; 7 8 9],
hcat(
reshape([1 1 2 0; 1 1 1 0], 1,1,2,4),
reshape([1 2 2 1; 1 2 1 1], 1,1,2,4)
))
3×3 CyclicArray{Int64, 2}:
1 2 3
4 5 6
7 8 9
julia> [klein[i,j] for i in -2:6, j in -2:6]
9×9 Matrix{Int64}:
3 2 1 3 2 1 3 2 1
6 5 4 6 5 4 6 5 4
9 8 7 9 8 7 9 8 7
1 2 3 1 2 3 1 2 3
4 5 6 4 5 6 4 5 6
7 8 9 7 8 9 7 8 9
3 2 1 3 2 1 3 2 1
6 5 4 6 5 4 6 5 4
9 8 7 9 8 7 9 8 7
3 Likes
The idea is similar to the xgcm python package:
1 Like
Just a note, CyclicArrays define a container with an abstract type. Which may lead to some performance issues. A way around this would be mimic CircularArrays and define an extra type parameter struct CyclicArray{T,N,A}…
This may not be important for some application cases, but defining Base.sum(ca::CyclicArray) = Base.sum(ca.data) should have some inference trouble in @code_warntype.
2 Likes
Here is a quick example to play with:
x0=permutedims(reshape((1:1:96)’,(4,4,6,1)),(3,4,2,1))
x1=CyclicArray(x0,“cubed”)
This example includes 6 faces interconnected in this order:
5
1 2 3 4. # → x direction
6
julia> x1[1,1,1,0]==x1[4,1,1,4]
true
julia> x1[3,1,1,0]==x1[2,1,1,4]
true
The package is still under development. The named grid topology is still not documented but it is active.
for searching purposes, there is now a registered Package that does this:
Note that with this type of approach (as in CircularArrays.jl and CyclicArrays.jl), you pay a computational price for every getindex operations (for a mod in CircularArrays and a much more expensive calculation in CyclicArrays).
In many circumstances, you don’t need to access the array at arbitrary indices, but only at indices in the interior and just outside the boundary, and in such cases you can get both greater efficiency and greater flexibility (in topology etcetera) by using ghost cells.
3 Likes
If you only need one-dimensional arrays, you can use CircularList.jl. It looks to have a really nice interface, where you can have a head pointing to a certain point in the list. (And the developer is great!)