Custom array with padding through a view

General Usage
1

I would like to pad my vector/array for doing FFT stuff, specifically mirroring - and since my data is pretty big, would be great if there was no need for extra allocations.

My thought was to have a custom array type, that would use clever indexing to have e.g. view of N first samples in the real array in reversed order as the first N samples of the custom array, then all the “real” samples in normal order and N reversed samples from the end similarly.
Following examples from packages that do similar things (like PaddedArray or ShiftedArray) I made custom properties (like size) and modified the getindex to do the clever-thing.
When I just use array[i] where 0 < i < 3*array_size it gives back the proper values, but Julia is not “tricked” and sees the original size.

Maybe I have missed a package that does this or there is something wrong with this approach?
Would be grateful for any hint or pointing it any direction.

2

Can you explain what you mean by this? In what way is Julia not “tricked?”

Let’s set up a type and a convenient constructor:

julia> struct MyVec{T} v::Vector{T} end

julia> Base.getindex(::Type{MyVec}, x...) = MyVec([x...])

julia> xs = MyVec[1,2,3]
MyVec{Int64}([1, 2, 3])

Now set up some behaviors you expect vectors to have:

julia> Base.length(v::MyVec) = 3length(v.v)

julia> Base.getindex(v::MyVec, i) = begin
           @assert 1 ≤ i ≤ length(v) "Please be considerate"
           if isodd((i-1)÷length(v.v))  v.v[mod1(i, length(v.v))]
           else  v.v[length(v.v) - mod1(i, length(v.v))+1]
           end
       end

julia> length(xs)
9

julia> [xs[i] for i=1:length(xs)]
9-element Vector{Int64}:
 3
 2
 1
 1
 2
 3
 3
 2
 1

You can also set up iteration utilities like so:

julia> Base.iterate(v::MyVec, n=1) = n ≤ length(v) ? (v[n], n+1) : nothing

julia> (xs...,)
(3, 2, 1, 1, 2, 3, 3, 2, 1)

And of course, for size:

julia> Base.size(v::MyVec) = (length(v),)

julia> size(xs)
(9,)

Sky’s the limit! (just wait until you hit OffsetArrays :wink:)

3

Hey, thanks for the response!

I actually deleted my non-working code, so can’t reproduce the exact problems - I just copied code for some other custom array type and tried to modify it to my needs.
Julia was not tricked in a sense, that doing MyArray([1,2,3]) would show a 3-element array object and not 9. Also iteration did not work.

However, bouncing of your MWE I think I got the idea (reading docs more closely also helped, as always).
I think I need only one mirror of the original array, so I have this:

struct MirrorVector{T} <: AbstractVector{T}
    source::Vector{T}
end

Base.length(vec::MirrorVector) = 2*length(vec.source)
Base.size(vec::MirrorVector) = (length(vec),)
Base.iterate(vec::MirrorVector, n=1) = (vec[n]. n+1)
Base.IndexStyle(vec::MirrorVector) = IndexLinear()

function Base.getindex(vec::MirrorVector, i)
    @boundscheck checkbounds(vec, i...)
    x = Int(length(vec.source)+0.5 - abs(length(vec.source)+0.5 - i))
    vec.source[x]
end

Which seems to work great. I have not tested it against FFTW, but I hope being a subtype of AbstractVector will allow it to work ok.

One thing I noticed is the penalty of using custom getindex logic.
I tested in on sum() with normal and mirrored vector with the same elements.

julia> @btime sum($NormalVector)
  14.400 μs (0 allocations: 0 bytes)

julia> @btime sum($MirrorVector)
  240.800 μs (0 allocations: 0 bytes)

And on one hand, it does not seem that bad. On the other, all of the difference is the extecution of the getindex code (which I gues is specially super-optimized for normal vector and does not appear anywhere in the profiler data). So I am wondering is there anything to be done here.

4

I don’t know specifically about FFTW but in general wrapped C libraries are not so amenable to Julia indexing tricks. They expect a certain memory layout, which is usually not more exotic than a strided array. The wrapper might copy your AbstractArray behind your back into a regular Array before sending it to the library, which makes things work, but then you don’t really gain much from your custom array.

5

You are definitely right, I just have not gotten that far earlier to face this problem. :wink:

And yes, FFTW requires a strided array (which I assume is the passed to the underlying C code that uses the pointers from the required methods?). Jerry-rigging it neverthless, gives gibberish results, as documentation warns. No options to translate the index-tricks to pointer-tricks in C?

However, this should work fine for a native Julia FFT, if it ever comes to be, right?

To be honest, this is probably solving a non-existing problem, since allocation of a doubling buffer that would be then reused, is only a small fraction of the execution time. But I am trying to push my limits in understanding Julia and that usually leads me for hours down the rabbit hole of trying different things (since I have little intuitions of what will work or not).

6

It should! Although to really capitalize on it, you’d want a custom FFT function too which is optimized for memory efficiency.

With all the brainiacs here, I don’t doubt that someday we’ll have a full-Julia FFT package. It could get really fun.

I wrote a toy FFT in fourteen lines (ten excluding comments) as a demo for a syntax proposal here and did some rough benchmarking. IIRC, on initial runs it was about 20% slower than FFTW with absolutely zero optimization, no compiler hints, no nothing. Not bad!

But then, on subsequent runs FFTW runs like ten times faster :sweat_smile: it turns out that most of the time is spent computing the twiddle factors (values of e^{-i\theta} where 0\le \theta<\pi). For a given dimensionality of data, the twiddle factors are constant so FFTW will memoize the values to allow later runs to be faster.

I was just toying with syntax, so I couldn’t be bothered to optimize it. Also, FFTW works with non-radix-2 data and multi-dimensional arrays… humanity has spent a fair bit of time developing the package so it’s darn near perfect.

But who knows! Maybe someday somebody will be motivated to port it to Julia (maybe to get automatic differentiation or something?), and then the real fun can begin.