When implementing a circular array, doesn’t it make much more sense to use 0-based indexing?
maybe, my toughts were to preserve the 1-based indexing of julia, and preserve the equality circular_x[i,j] == x[i,j]
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.
And there is the mod1
function to avoid all this shifting indices manually:
julia> mod1(5,5)
5
julia> mod1(6,5)
1
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?
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
This is material for a small compact package!
Is there a package with this circular array type available and well-tested already?
@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.
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
The idea is similar to the xgcm python package:
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
.
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.
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!)