# Crazy allocations using CartesianIndices

**URL:** <https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262>\
**Category:** General Usage\
**Tags:** question, indexing, memory-allocation\
**Created:** [June 29, 2020, 6:17pm UTC](https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262 "2020-06-29T18:17:08Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [June 29, 2020, 6:17pm UTC](https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262/1 "2020-06-29T18:17:08Z")

</div>

I have been switching over a 2D flux computation method to use CartesianIndices so it will work for arrays of any dimension, ala [Multidimensional algorithms and iteration](https://julialang.org/blog/2016/02/iteration/)

Unfortunately, I’ve been running into a bunch of memory allocation “gotchas” and I don’t understand what I need to do to avoid them. Here’s an example

```julia
function test(a)
    σ = 0.
    R = CartesianIndex(size(a)[1:end-1]); o = oneunit(R)
    @time for I ∈ 2o:R
        σ += a[I,1]-a[I-o,2]^2 # actually, a more complex function...
    end
    return σ
end

```

But when I run this I get slow code with a giant memory issue

```julia
julia> test(rand(256,256,2));
  0.046946 seconds (585.23 k allocations: 13.891 MiB)

```

Please, could someone let me know:

1. What is going on here?
2. How am I supposed to do this elegantly?

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [June 29, 2020, 6:32pm UTC](https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262/2 "2020-06-29T18:32:51Z")

</div>

If you run `@code_warntype test(rand(256, 256, 2))` you’ll see a lot of type-instabilities:

```julia
julia> @code_warntype test(rand(256, 256, 2))
Variables
  #self#::Core.Compiler.Const(test, false)
  a::Array{Float64,3}
  stats::Base.GC_Num
  elapsedtime::UInt64
  val::Nothing
  diff::Base.GC_Diff
  σ::Any
  R::CartesianIndex{_A} where _A
  o::CartesianIndex{_A} where _A
  @_10::Any
  I::Any

Body::Any

```

Those ultimately stem from `size(a)[1:end-1]`, which takes a tuple and then slices it into a vector. The problem here is that the type of `CartesianIndex` depends on the number of dimensions, but the length of that slice of `size(a)` is not known to the compiler (it depends on the run-time value of `end-1`).

Fixing that removes all the allocations:

```julia
julia> function test(a)
           σ = 0.
           R = CartesianIndex((size(a, 1), size(a, 2))); o = oneunit(R)
           @time for I ∈ 2o:R
               σ += a[I,1]-a[I-o,2]^2 # actually, a more complex function...
           end
           return σ
       end

```

```julia
julia> test(rand(256,256,2));
  0.000243 seconds

```

The only issue here is that I’ve implicitly assumed that `a` has 3 dimensions, which your original code did not. We can make a sufficiently general version with a helper function:

```julia
function all_but_last(t::NTuple{N, T}) where {N, T}
  ntuple(i -> t[i], N - 1)
end

```

and use it like so:

```julia
julia> function test(a)
           σ = 0.
           R = CartesianIndex(all_but_last(size(a))); o = oneunit(R)
           @time for I ∈ 2o:R
               σ += a[I,1]-a[I-o,2]^2 # actually, a more complex function...
           end
           return σ
       end
test (generic function with 1 method)

julia> test(rand(256,256,2));
  0.000218 seconds

```

The moral is that the length of a tuple is part of its type, so be wary of operations that might accidentally throw away that information by converting to a `Vector`.

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [June 29, 2020, 6:52pm UTC](https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262/3 "2020-06-29T18:52:32Z")

</div>

That’s a really helpful answer!

As you say, I need the code to work regardless of dimension, but your helper functions works great. However, I’m a little confused since this slicing was used in the blog post

```julia
function expfiltdim(x, dim::Integer, α)
    s = similar(x)
    Rpre = CartesianIndices(size(x)[1:dim-1])
    Rpost = CartesianIndices(size(x)[dim+1:end])
    _expfilt!(s, x, α, Rpre, size(x, dim), Rpost)
end

```

Would helpers like this avoid the issues even though `dim` is a variable?

```julia
function before_dim(dim,t::NTuple{N, T}) where {N, T}
  ntuple(i -> t[i], dim - 1)
end
function after_dim(dim,t::NTuple{N, T}) where {N, T}
  ntuple(i -> t[dim+i], N - dim)
end

```

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [June 29, 2020, 7:11pm UTC](https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262/4 "2020-06-29T19:11:49Z")

</div>

> [@weymouth](#):
>
> However, I’m a little confused since this slicing was used in the blog post

The example there has the same type-instability, but it doesn’t cause a problem because `_expfilt!` is creating a [“function barrier”](https://docs.julialang.org/en/v1/manual/performance-tips/index.html#kernel-functions-1). Even though the type of `Rpre` in `expfiltdim` cannot be inferred, the `_expfilt!` function still receives a value with a concrete type, so within the body of `_expfilt!` there is no type-instability. Check out the section of the manual I linked for more details. You could get the same benefit by moving your loop into its own dedicated function, although I personally think fixing the instability altogether is good practice when possible.

> [@weymouth](#):
>
> Would helpers like this avoid the issues even though `dim` is a variable?

It’s hard to say–it depends on how clever the compiler is with its constant propagation. If you’re dealing with cases like this, then the function barrier approach seems like a good idea.

---

<div class="post-metadata">

**Author:** ![Seif\_Shebl](https://avatars.discourse-cdn.com/v4/letter/s/eada6e/32.png) [@Seif\_Shebl](https://discourse.julialang.org/u/Seif_Shebl)\
**Post date:** [June 30, 2020, 1:56am UTC](https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262/5 "2020-06-30T01:56:13Z")

</div>

By the way, I don’t see any instabilities with Julia 1.5.0-beta1.0, on what version did you test your code:

```julia
julia> @code_warntype test(rand(256, 256, 2))
Variables
  #self#::Core.Compiler.Const(test, false)
  a::Array{Float64,3}
  stats::Base.GC_Num
  elapsedtime::UInt64
  val::Nothing
  diff::Base.GC_Diff
  σ::Float64
  R::CartesianIndex{2}
  o::CartesianIndex{2}
  @_10::Union{Nothing, Tuple{CartesianIndex{2},CartesianIndex{2}}}
  I::CartesianIndex{2}

Body::Float64
. . .

```

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [June 30, 2020, 2:22am UTC](https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262/6 "2020-06-30T02:22:36Z")

</div>

It’s been improved on 1.5, all prior versions fail to infer it. It’s worth noting that there are still some tuple-indexing behaviors that _aren’t_ inferred correctly, but I think it’s fair to say that all the most common ones are.

[https://github.com/JuliaLang/julia/pull/31138](https://github.com/JuliaLang/julia/pull/31138)

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [June 30, 2020, 6:59am UTC](https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262/7 "2020-06-30T06:59:59Z")

</div>

That sounds like it could save me a bunch of memory headaches. Is that release stable?

Since I’m trying to make this a package other people can use with other versions, I am also a bit nervous solving the problem this way.

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [June 30, 2020, 7:04am UTC](https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262/8 "2020-06-30T07:04:42Z")

</div>

If you want to support multiple Julia versions, you’ll have to handle it the “hard way.” Aside from `ntuple`, there’s `Base.tail`, `Base.front`, and `@inline`d splatting. I can’t pull up memory of a blog post on “lispy tuple programming,” anyone have suggestions?

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [June 30, 2020, 8:00am UTC](https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262/9 "2020-06-30T08:00:44Z")

</div>

IMHO Julia is improving fast and for a brand new package like yours, I do not see why you should pay the price with compatibility with old versions… Of course it is totally up to you.

I have spent some time reading your code which is extremely elegant and concise. Coming from C++, I am still amazed that such a dimension agnostic coding style can be obtained at zero cost. If it is the case with Julia 1.5, it should be kept like this.

I also wonder about, the future adaption to GPU. It could be interesting to see if construction like @vchuravy [https://github.com/JuliaGPU/KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl) could GPUify your code efficiently (although I think that this package is not yet compatible with Julia 1.5 nor 1.4).

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 30, 2020, 9:21am UTC](https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262/10 "2020-06-30T09:21:34Z")

</div>

> [@tim.holy](#):
>
> I can’t pull up memory of a blog post on “lispy tuple programming,” anyone have suggestions?

The tricks at the end of

> <https://github.com/mitmath/18S096/blob/master/lectures/lecture6/Types%20and%20Dispatch.ipynb>

cover a lot of cases.

Doing the exercises is recommended, and automatically confers a green belt to the user.

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [June 30, 2020, 9:49am UTC](https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262/11 "2020-06-30T09:49:17Z")

</div>

😀 I used to lead the MIT TKD club, so this is full circle.

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [June 30, 2020, 1:26pm UTC](https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262/12 "2020-06-30T13:26:22Z")

</div>

Ok, I’m going to update to 1.5 (downloaded already) since I don’t enjoy tracking down all these memory thing.

However, I have one more chestnut which I would like to offer up to the hive mind. I still can’t get the function barrier approach to work for applying boundary conditions

```julia
function _apply!(a,dim::Integer,left,right,A)
    for i ∈ left, j ∈ right
        a[i,size(a,dim),j, dim] = A
    end
end
function test(a::Array{T,N},A) where {T,N}
    for dim ∈ 1:N-1
        left = CartesianIndices(size(a)[1:dim-1])
        right = CartesianIndices(size(a)[dim+1:N-1])
        @time _apply!(a,dim,left,right,A)
    end
    return a
end
julia> test(ones(20,20,20,3),0.);
  0.000299 seconds (2.80 k allocations: 93.766 KiB)
  0.000003 seconds (1 allocation: 16 bytes)
  0.000002 seconds (1 allocation: 16 bytes)

```

If I try the `Val{dim}` trick from that ipynb tutorial, I get

```julia
function _apply!(a,::Val{dim},left,right,A) where dim
    for i ∈ left, j ∈ right
        a[i,size(a,dim),j, dim] = A
    end
end
_apply!(a,dim::Integer,left,right,A) = (a,Val{dim}(),left,right,A)

julia> test(ones(20,20,20,3),0.);
  0.000026 seconds (4 allocations: 112 bytes)
  0.000048 seconds (5 allocations: 112 bytes)
  0.000092 seconds (4 allocations: 112 bytes)

```

The @code\_typed is also a mess for a fairly simple idea: `a[:,size(a,2),:,2]=A` etc. Maybe it would be better to do this with code generation?

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [June 30, 2020, 3:00pm UTC](https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262/13 "2020-06-30T15:00:04Z")

</div>

Generic splits are not handled well yet. There is `Base.IteratorsMD.split`, in case it’s useful:

```julia
julia> I = CartesianIndices((5,10,15,20)).indices
(Base.OneTo(5), Base.OneTo(10), Base.OneTo(15), Base.OneTo(20))

julia> pre, post = Base.IteratorsMD.split(I, Val(2));

julia> pre
(Base.OneTo(5), Base.OneTo(10))

julia> post
(Base.OneTo(15), Base.OneTo(20))

```

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [July 2, 2020, 10:23pm UTC](https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262/14 "2020-07-02T22:23:42Z")

</div>

I couldn’t get this to work satisfactorily for the BCs, but the rest of the solver is nicely dimensionally independent. Here’s a slice through 3D Taylor Green Vortex.  
 ![TGVortex](https://global.discourse-cdn.com/julialang/original/3X/1/3/13c683cb8cbfe16ce69f72573ef16543b9b5782b.gif)

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [March 29, 2021, 12:44pm UTC](https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262/15 "2021-03-29T12:44:33Z")

</div>

I know this was long ago and far away, but I looked at it again this weekend and came up with a good solution

```julia
function slice(N::NTuple{n,Int},s::Int,dims::Int)::CartesianIndices{n} where n
    CartesianIndices(ntuple( i-> i==dims ? (s:s) : (1:N[i]), n))
end

```

This let’s me write general multi-dimensional code applied to the boundaries of an array; like setting boundary conditions:

```julia
function BC!(a::Array{T,n}) where {T,n}
    N = size(a)
    for j ∈ 1:n
        @simd for I ∈ slice(N,1,j)
            a[I] = a[I+δ(j,n)] # Neumann
        end
        @simd for I ∈ slice(N,N[j],j)
            a[I] = a[I-δ(j,n)] # Neumann
        end
    end
end
@inline δ(a,d::Int) = CartesianIndex(ntuple(i -> i==a ? 1 : 0, d))

```