# Indexing on a pseudo-dimension

**URL:** <https://discourse.julialang.org/t/indexing-on-a-pseudo-dimension/101119>\
**Category:** General Usage\
**Tags:** indexing, arrays\
**Created:** [July 3, 2023, 12:09pm UTC](https://discourse.julialang.org/t/indexing-on-a-pseudo-dimension/101119 "2023-07-03T12:09:10Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![croberts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/croberts/32/9465_2.png) [@croberts](https://discourse.julialang.org/u/croberts)\
**Post date:** [July 3, 2023, 12:09pm UTC](https://discourse.julialang.org/t/indexing-on-a-pseudo-dimension/101119/1 "2023-07-03T12:09:10Z")

</div>

Suppose I have a K x M array, my\_arr, but I want to pretend that it is a K x M x N array, my\_pseudo\_arr , so that indexing it with my\_pseudo\_arr[k, m, n] equals my\_arr[k, m], my\_pseudo\_arr[k, m, :] equals repeat([my\_arr[k,m]], N), and my\_pseudo\_arr[k, m, n+1] throws an error.

Obviously, I could construct a new higher dimensional array where I repeat my\_arr N times. But it seems computationally wasteful to construct this thing when I could solve this problem simply by reinterpretting the getindex function.

Can someone point me in the right direction? Is there existing functionality in Base or some other package that makes this trivial?

Thanks!

---

<div class="post-metadata">

**Author:** ![Henrique\_Becker](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henrique_becker/32/15443_2.png) [@Henrique\_Becker](https://discourse.julialang.org/u/Henrique_Becker)\
**Post date:** [July 3, 2023, 1:27pm UTC](https://discourse.julialang.org/t/indexing-on-a-pseudo-dimension/101119/2 "2023-07-03T13:27:29Z")

</div>

You can create a new type inheriting `AbstractArray` that just wraps an existing array, and then define its `getindex`. I am not sure if there is a package that already does that.

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [July 4, 2023, 9:38am UTC](https://discourse.julialang.org/t/indexing-on-a-pseudo-dimension/101119/3 "2023-07-04T09:38:36Z")

</div>

Also as an exercise for myself, didn’t find a package:

```julia
 struct DimPlusOneArray{T, N, AA} <: AbstractArray{T, N}
     parent::AA
     lengthD::Int
     size::NTuple{N, Int}
     function DimPlusOneArray(A, lengthD)
         sz = ntuple(i -> i ≤ ndims(A) ? size(A, i) : lengthD,
                           Val(ndims(A) + 1))  
         new{eltype(A), ndims(A) + 1, typeof(A)}(A, lengthD, sz) 
     end 
 end
 
 @inline function Base.getindex(M::DimPlusOneArray{T, N, AA},
                                                    I::Vararg{Int, N}) where {T, N, AA} 
     @boundscheck checkbounds(M, I...)
     return Base.getindex(M.parent, I[begin:end-1]...)
 end
 
 Base.size(A::DimPlusOneArray{T, N, AA}) where {T, N, AA} = A.size

```

It works, but there is a single allocation for array access. Not sure where it comes from:

```julia
julia> @time A = DimPlusOneArray(randn((2,2)), 3)
  0.000003 seconds (2 allocations: 144 bytes)
2×2×3 DimPlusOneArray{Float64, 3, Matrix{Float64}}:
[:, :, 1] =
 -0.358432 2.17364
 -0.174398 0.48638

[:, :, 2] =
 -0.358432 2.17364
 -0.174398 0.48638

[:, :, 3] =
 -0.358432 2.17364
 -0.174398 0.48638

julia> @time A[1, 2, 2]
  0.000004 seconds (1 allocation: 16 bytes)
2.1736430112129765

julia> @time A[1, 2, 1]
  0.000005 seconds (1 allocation: 16 bytes)
2.1736430112129765

julia> @time A[1, 1, 1]
  0.000004 seconds (1 allocation: 16 bytes)
-0.35843209612481847

julia> @time A[1, 2, 5]
ERROR: BoundsError: attempt to access 2×2×3 DimPlusOneArray{Float64, 3, Matrix{Float64}} at index [1, 2, 5]
Stacktrace:
 [1] throw_boundserror(A::DimPlusOneArray{Float64, 3, Matrix{Float64}}, I::Tuple{Int64, Int64, Int64})
   @ Base ./abstractarray.jl:744
 [2] checkbounds
   @ ./abstractarray.jl:709 [inlined]
 [3] getindex(::DimPlusOneArray{Float64, 3, Matrix{Float64}}, ::Int64, ::Int64, ::Int64)
   @ Main ~/.julia/dev/RepeatedArrays.jl/RepeatedArrays.jl:16
 [4] top-level scope
   @ ./timing.jl:273 [inlined]
 [5] top-level scope
   @ ./REPL[13]:0

```

---

<div class="post-metadata">

**Author:** ![Henrique\_Becker](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henrique_becker/32/15443_2.png) [@Henrique\_Becker](https://discourse.julialang.org/u/Henrique_Becker)\
**Post date:** [July 4, 2023, 1:27pm UTC](https://discourse.julialang.org/t/indexing-on-a-pseudo-dimension/101119/4 "2023-07-04T13:27:22Z")

</div>

> [@roflmaostc](#):
>
> `I[begin:end-1]...`

I would guess that it is from actually creating an intermediary array from the slice form? Did you try using `@view` (not sure if it will work with Vararg)? Or maybe this needs a generator function so there is guarantee the code will specialize the splat using the `N` information.

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [July 4, 2023, 1:34pm UTC](https://discourse.julialang.org/t/indexing-on-a-pseudo-dimension/101119/5 "2023-07-04T13:34:11Z")

</div>

> [@Henrique\_Becker](#):
>
> creating an intermediary array

That’s a tuple, not an array.

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [July 4, 2023, 1:39pm UTC](https://discourse.julialang.org/t/indexing-on-a-pseudo-dimension/101119/6 "2023-07-04T13:39:17Z")

</div>

There are certainly packages, e.g. mine (which [lists alternatives in its readme](https://github.com/mcabbott/LazyStack.jl#other-stack-like-packages)):

```julia
julia> mat = [1 2 3; 4 5 6];

julia> using LazyStack

julia> lazystack((mat, mat))
2×3×2 lazystack(::Tuple{Matrix{Int64}, Matrix{Int64}}) with eltype Int64:
[:, :, 1] =
 1 2 3
 4 5 6

[:, :, 2] =
 1 2 3
 4 5 6

julia> lazystack(fill(mat, 10))
2×3×10 lazystack(::Vector{Matrix{Int64}}) with eltype Int64:
[:, :, 1] =
 1 2 3
 4 5 6
...

```

However, operations consuming this lazy array will sometimes be much less optimised than the same operations on an Array. What are you doing next?

It is almost always better (IMO) to adjust the next operation to accept the data you really have, rather than wrap the data in something like this.

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [July 4, 2023, 1:43pm UTC](https://discourse.julialang.org/t/indexing-on-a-pseudo-dimension/101119/7 "2023-07-04T13:43:05Z")

</div>

> [@roflmaostc](#):
>
> there is a single allocation for array access

I think that’s a measurement problem, `@btime $A[1, 2, 2]` has zero.

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [July 4, 2023, 1:44pm UTC](https://discourse.julialang.org/t/indexing-on-a-pseudo-dimension/101119/8 "2023-07-04T13:44:55Z")

</div>

Same behaviour for ShiftedArrays

```julia
julia> M2 = ShiftedArrays.circshift(M, (1,2,3))
2×2×2 CircShiftedArray{Float64, 3, Array{Float64, 3}}:
[:, :, 1] =
 -0.140213 0.981398
  0.593498 -1.45221

[:, :, 2] =
 -0.870943 -1.08996
  1.00483 1.72976

julia> using BenchmarkTools

julia> @time M2[1,1,1]
  0.000008 seconds (1 allocation: 16 bytes)
-0.14021348397210173

julia> @time M2[1,1,1]
  0.000017 seconds (1 allocation: 16 bytes)
-0.14021348397210173

julia> @btime $M2[1,1,1]
  3.617 ns (0 allocations: 0 bytes)
-0.14021348397210173

```

---

<div class="post-metadata">

**Author:** ![Henrique\_Becker](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henrique_becker/32/15443_2.png) [@Henrique\_Becker](https://discourse.julialang.org/u/Henrique_Becker)\
**Post date:** [July 4, 2023, 1:49pm UTC](https://discourse.julialang.org/t/indexing-on-a-pseudo-dimension/101119/9 "2023-07-04T13:49:26Z")

</div>

True, but an intermediary tuple should not allocate, so maybe lack of inference led the `getindex`/`...` to be dynamically dispatched? The problem is, `@code_typed` gave me no clue of some problem with inference.

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [July 4, 2023, 2:41pm UTC](https://discourse.julialang.org/t/indexing-on-a-pseudo-dimension/101119/10 "2023-07-04T14:41:52Z")

</div>

> [@roflmaostc](#):
>
> there is a single allocation for array access. Not sure where it comes from

Minimal examples, the difference is in the use of the `parent` function:

1. Allocation

```julia
struct S{T, n, A <: AbstractArray{T}} <: AbstractArray{T, n}
  parent::A
  function S{T,n,A}(a::A) where {T, n, m, A <: AbstractArray{T, m}}
    new{T, n, A}(a)
  end
end

S(a::A) where {T, m, A <: AbstractArray{T, m}} = S{T, m + 1, A}(a)

Base.getindex(a::S{<:Any, n}, I::Vararg{Int, n}) where {n} = (a.parent)[Base.front(I)...]

Base.size(a::S) = (size(a.parent)..., 5)

const arr = S(rand(2))

@allocated arr[1, 2]
@allocated arr[1, 2] # nonzero

```

1. No allocation

```julia
struct S{T, n, A <: AbstractArray{T}} <: AbstractArray{T, n}
  parent::A
  function S{T,n,A}(a::A) where {T, n, m, A <: AbstractArray{T, m}}
    new{T, n, A}(a)
  end
end

S(a::A) where {T, m, A <: AbstractArray{T, m}} = S{T, m + 1, A}(a)

parent(a::S) = a.parent

Base.getindex(a::S{<:Any, n}, I::Vararg{Int, n}) where {n} = parent(a)[Base.front(I)...]

Base.size(a::S) = (size(a.parent)..., 5)

const arr = S(rand(2))

@allocated arr[1, 2]
@allocated arr[1, 2] # zero

```

~~Will report a bug now.~~ EDIT: possibly related bug reports:

> <https://github.com/JuliaLang/julia/issues/49241>
>
> The code below does not allocate in Julia 1.8, but it allocates 16 bytes in Juli…a 1.9 and 48 in Julia 1.10.
> \`\`\`julia
> struct S
> a::Vector{Int}
> end
> 
> const i = S(Int\[\])
> 
> function f(a::S)::Int
> isempty(a.a) && return 1
> return error("we")
> end
> @allocated(f(i))
> 16
> \`\`\`
> 
> \`\`\`
> versioninfo()
> Julia Version 1.9.0-rc2
> Commit 72aec423c2a (2023-04-01 10:41 UTC)
> Platform Info:
> OS: Linux (x86\_64-linux-gnu)
> CPU: 12 × Intel(R) Core(TM) i7-9750H CPU @ 2.60GHz
> WORD\_SIZE: 64
> LIBM: libopenlibm
> LLVM: libLLVM-14.0.6 (ORCJIT, skylake)
> Threads: 1 on 12 virtual core
> \`\`\`

> <https://github.com/JuliaLang/julia/issues/49697>
>
> The execution time is similar for this small example, but there are 10x more all…ocations on 1.9-rc3.
> 
> \`\`\`
> using BenchmarkTools
> using Plots
> const xx = randn(20, 3000)
> @btime plot(xx; color = :blue, label = nothing);
> 
> \# 1.8.5
> \# 1.214 s (392640 allocations: 31.88 MiB)
> 
> \# 1.9-rc3
> \# 1.166 s (3964617 allocations: 86.43 MiB)
> \`\`\`
> 
> From https://github.com/JuliaPlots/Plots.jl/issues/4742

> <https://github.com/JuliaLang/julia/issues/50073>
>
> On upgrading to Julia from 1.8.5. to 1.9.0 I noticed code in one of my packages …slow down significantly. (I've already highlighted a type inference issue that I found. This is an entirely separate issue.) After upgrading, the following code suddenly causes a lot of memory allocations. This is the simplest code which replicates the problem, but most minor changes will suddenly cause all the allocations to disappear:
> 
> \`\`\`julia
> using StaticArrays, ForwardDiff, BenchmarkTools, JET
> 
> function rodrigues(x::T, y::T, z::T) where T
> if x == 0 && y == 0 && z == 0
> # Short cut for derivatives at identity
> return SMatrix{3, 3, T, 9}(T(1), z, -y, -z, T(1), x, y, -x, T(1))
> end
> theta2 = x \* x + y \* y + z \* z
> cosf = T(0.5)
> sinc = T(1)
> if theta2 \> T(2.23e-16)
> theta = sqrt(theta2)
> sinc, cosf = sincos(theta)
> cosf -= 1
> sinc /= theta
> cosf /= -theta2
> end
> a = x \* y \* cosf
> b = sinc \* z
> c = x \* z \* cosf
> d = sinc \* y
> e = y \* z \* cosf
> f = sinc \* x
> return SMatrix{3, 3, T, 9}((x \* x - theta2) \* cosf + 1, a + b, c - d,
> a - b, (y \* y - theta2) \* cosf + 1, e + f,
> c + d, e - f, (z \* z - theta2) \* cosf + 1)
> end
> 
> struct BarrelDistortion{T}
> k1::T
> k2::T
> end
> update(var::BarrelDistortion, updatevec, start=1) = BarrelDistortion(var.k1 + updatevec\[start\], var.k2 + updatevec\[start+1\])
> 
> function ideal2distorted(lens::BarrelDistortion, x)
> z = x' \* x
> z = z \* (lens.k1 + z \* lens.k2) + 1
> return x \* z
> end
> 
> struct SimpleCamera{T}
> f::T
> end
> SimpleCamera(v::T) where T = SimpleCamera{T}(v::T)
> update(var::SimpleCamera, updatevec, start=1) = SimpleCamera(max(var.f, floatmin(var.f)) \* exp(updatevec\[start\]))
> 
> struct Point3D{T}
> v::SVector{3, T}
> end
> Point3D(x, y, z) = Point3D(SVector{3}(x, y, z))
> Point3D() = Point3D(SVector{3}(0., 0., 0.))
> update(var::Point3D, updatevec, start=1) = Point3D(var.v + updatevec\[StaticArrays.SUnitRange(0, 2) .+ start\])
> project(x::Point3D{T}) where T = SVector{2, T}(x.v\[1\], x.v\[2\]) ./ x.v\[3\]
> 
> struct Rotation{T}
> m::SMatrix{3, 3, T, 9}
> end
> Rotation(x, y, z) = Rotation(rodrigues(x, y, z))
> update(var::Rotation, updatevec, start=1) = transform(Rotation(updatevec\[start\], updatevec\[start+1\], updatevec\[start+2\]), var)
> transform(rota::Rotation, rotb::Rotation) = Rotation(rota.m \* rotb.m)
> 
> struct EffPose3D{T}
> rot::Rotation{T}
> camcenter::Point3D{T}
> end
> EffPose3D(rx, ry, rz, cx, cy, cz) = EffPose3D(Rotation(rx, ry, rz), Point3D(cx, cy, cz))
> update(var::EffPose3D, updatevec, start=1) = EffPose3D(update(var.rot, updatevec, start), update(var.camcenter, updatevec, start+3))
> transform(pose::EffPose3D, point::Point3D) = Point3D(pose.rot.m \* (point.v - pose.camcenter.v))
> 
> \# Description of BAL image, and function to transform a landmark from world coordinates to pixel coordinates
> struct Image{T}
> pose::EffPose3D{T}
> sensor::SimpleCamera{T}
> lens::BarrelDistortion{T}
> end
> update(var::Image, updatevec, start=1) = Image(update(var.pose, updatevec, start), update(var.sensor, updatevec, start+6), update(var.lens, updatevec, start+7))
> function Image(rx::T, ry::T, rz::T, tx::T, ty::T, tz::T, f::T, k1::T, k2::T) where T
> R = Rotation(rx, ry, rz)
> return Image{T}(EffPose3D(R, Point3D(R.m' \* -SVector(tx, ty, tz))), SimpleCamera(f), BarrelDistortion(k1, k2))
> end
> transform(im::Image, X::Point3D) = ideal2distorted(im.lens, -project(Point3D(im.pose.rot.m \* (X.v - im.pose.camcenter.v))))
> 
> computeresjac(vars...) = ForwardDiff.jacobian(z -\> transform(update(vars\[1\], z, 1), update(vars\[2\], z, 10)), zeros(SVector{12, Float64}))
> 
> function mytest()
> im = Image(1., 1., 1., 1., 1., 1., 1., 1., 1.)
> lm = Point3D(0., 0., 0.)
> @btime computeresjac($im, $lm)
> show(JET.@report\_opt computeresjac(im, lm))
> end
> 
> mytest()
> \`\`\`

> <https://github.com/JuliaLang/julia/issues/50317>
>
> Here's a small regression that is showing up pretty large in the RAI codebase (w…e have a global lock in our database pager that's way slower now).
> 
> Somehow accessing the fields of a global constant has become type unstable in some cases.
> 
> \`\`\`julia
> struct Wrapper
> lock::ReentrantLock
> end
> 
> Base.lock(f::Function, m::Wrapper) = Base.lock(f, m.lock)
> Base.lock(m::Wrapper) = Base.lock(m.lock)
> Base.unlock(m::Wrapper) = Base.unlock(m.lock)
> 
> const MONITOR = Wrapper(ReentrantLock())
> 
> function foo()
> Base.@lock MONITOR begin
> return 2+2
> end
> end
> foo()
> @time foo() # 0 allocs on 1.8, 1 alloc on 1.9
> \`\`\`
> The allocation only shows up when referencing the global const variable. It disappears if you pass it in as an argument:
> \`\`\`julia
> function bar(x)
> Base.@lock x begin
> return 2+2
> end
> end
> bar(MONITOR)
> @time bar(MONITOR) # 0 allocs on 1.8, 0 allocs on 1.9
> \`\`\`
> 
> This might be related to https://github.com/JuliaLang/julia/issues/50073, but i'm not sure.

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [July 4, 2023, 2:43pm UTC](https://discourse.julialang.org/t/indexing-on-a-pseudo-dimension/101119/11 "2023-07-04T14:43:21Z")

</div>

BTW, regarding the original question, the documentation is here (the first two sections can be skipped):

[https://docs.julialang.org/en/v1/manual/interfaces/](https://docs.julialang.org/en/v1/manual/interfaces/)