SeparableFunctions.jl

Calculates multidimensional functions faster by exploiting their separability. Often a function involves an operation such as a (complex) exponential which by itself is computationally relatively heavy. Yet a number of multidimenstional functions are separable, which means that they can be written as a product (or sum) of single-dimensional functions. A good example of a separable function is a Gaussian function, which can be written as a product of a purely X-dependend Gaussian with a purely Y-dependend Gaussian.

In this package, multidimensional functions are computed by first calculating their single-dimensional values and then creating the final multidimensional result by an outer product. Since multiplications and the broadcasting mechanism of Julia are fast compared to the evaluation of the function at each multidimensional position, the final result is calculated faster. The typical speedup can be an order of magnitude.

The package offers a general way of calculating separable functions as well as a LazyArrays version of that function which can then be used inside other expressions. The non-lazy version should currently also work with CUDA.jl, however the LazyArrays version does not. To nevertheless use separable expressions in CUDA.jl, you can reside to externally applying the broadcast operator to the separable expression (see the gaussian_sep example below).

The package further offers a number of predifined separable function implementations such as gaussian_col() collects a multidimensional array of a multidimensional Gaussian via a fast seperable implementation, gaussian_lz() yields a lazy representation via LazyArrays and sep = gaussian_sep() yields an iterable of separable pre-oriented vectors which can easily be mutually combined via .*(sep...).

Another noteworthy example is the complex plain wave as represented by the respective function exp_ikx_col(), exp_ikx_lz(), exp_ikx_sep().

All separable functions share a common interface with Standard non-named arguments

• the first optional argument being the type of the result array. Examples are Array{Float64}, CuArray{Float32} and the default depends on the function but uses a 32-bit result type where applicable.
• the next argument is the size of the result data supplied as a tuple.
• optionally a further argument can specify the position of zero of the array with respect to offset as given below. This allows for convenient N-dimensional shifting of the functions.

Named arguments:

• offset: is a optional named argument specifying center of the array. By default the Fourier center, size().÷2 .+ 1 is chosen.
• scale: multiplies the axes with these factors. This can be interpreted as a pixelsize.

Some functions have additional named arguments. E.g. gaussian has the additional named argument sigma specifying the width of the gaussian, even though this can also be controlled via scale. Note that this nomenclature is in large parts identical to the package IndexFunArrays.jl.

In general arguments can be supplied as single scalar values or vectors. If a scalar value is supplied, it will automatically be replicated as a vector. E.g. sigma=10.0 for a 2-dimensional array will be interpreted as sigma=(10.0, 10.0).

calculate_broadcasted([::Type{TA},] fct, sz::NTuple{N, Int}, args...; dims=1:N, pos=zero(real(eltype(DefaultArrType))), offset=sz.÷2 .+1, scale=one(real(eltype(DefaultArrType))), operation = *, kwargs...) where {TA, N}

returns an instantiated broadcasted separable array, which essentially behaves almost like an array yet uses broadcasting. Test revealed maximal speed improvements for this version. Yet, a problem is that reduce operations with specified dimensions cause an error. However this can be avoided by calling collect.

Arguments:

• fct: The separable function, with a number of arguments corresponding to length(args), each corresponding to a vector of separable inputs of length N. The first argument of this function is a Tuple corresponding the centered indices.
• sz: The size of the N-dimensional array to create
• args...: a list of arguments, each being an N-dimensional vector
• dims: a vector [] of valid dimensions. Only these dimension will be calculated but they are oriented in ND.
• offset: position of the center from which the position is measured
• scale: defines the pixel size as vector or scalar. Default: 1.0.
• operation: the separable operation connecting the separable dimensions

julia> fct = (r, sz, pos, sigma)-> exp(-(r-pos)^2/(2*sigma^2))
julia> my_gaussian = calculate_broadcasted(fct, (6,5), (0.1,0.2), (0.5,1.0))
Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}}(*, 
(Float32[4.4963495f-9; 0.00014774836; … ; 0.1978987; 0.0007318024;;], 
Float32[0.088921614 0.48675224 … 0.726149 0.1978987]))
julia> collect(my_gaussian)
6×5 Matrix{Float32}:
 3.99823f-10  2.18861f-9   4.40732f-9   3.26502f-9   8.89822f-10
 1.3138f-5    7.19168f-5   0.000144823  0.000107287  2.92392f-5
 0.00790705   0.0432828    0.0871608    0.0645703    0.0175975
 0.0871608    0.477114     0.960789     0.71177      0.19398
 0.0175975    0.0963276    0.19398      0.143704     0.0391639
 6.50731f-5   0.000356206  0.000717312  0.000531398  0.000144823

calculate_separables([::Type{AT},] fct, sz::NTuple{N, Int}, args...; dims = 1:N, pos=zero(real(eltype(AT))), offset=sz.÷2 .+1, scale=one(real(eltype(AT))), kwargs...) where {AT, N}

creates a list of one-dimensional vectors, which can be combined to yield a separable array. In a way this can be seen as a half-way Lazy operation. The (potentially heavy) work of calculating the one-dimensional functions is done now but the memory-heavy calculation of the array is done later. This function is used in separable_view and separable_create.

#Arguments

• AT: optional type signfying the array result type. You can for example use CuArray{Float32} using CUDA to create the views on the GPU.
• fct: the function to calculate for each axis index (no need for broadcasting!) of this iterable of seperable axes. Note that the first arguments of fct have to be the index of this coordinate and the size of this axis. Any further args and nargs can follow. Often the second argument is not used but it still needs to be present.
• sz: the size of the result array (when appying the one-D axes)
• args: further arguments which are passed over to the function fct.
• dims: a vector [] of valid dimensions. Only these dimension will be calculated but they are oriented in ND.
• pos: a position shifting the indices passed to fct in relationship to the offset.
• offset: specifying the center (zero-position) of the result array in one-based coordinates. The default corresponds to the Fourier-center.
• scale: multiplies the index before passing it to fct

#Example:

julia> fct = (r, sz, sigma)-> exp(-r^2/(2*sigma^2))
julia> gauss_sep = calculate_separables(fct, (6,5), (0.5,1.0), pos = (0.1,0.2))
2-element Vector{Array{Float32}}:
 [4.4963495f-9; 0.00014774836; … ; 0.1978987; 0.0007318024;;]
 [0.088921614 0.48675224 … 0.726149 0.1978987]
 julia> my_gaussian = .*(gauss_sep...) # this is how to broadcast it
 6×5 Matrix{Float32}:
 3.99823f-10  2.18861f-9   4.40732f-9   3.26502f-9   8.89822f-10
 1.3138f-5    7.19168f-5   0.000144823  0.000107287  2.92392f-5
 0.00790705   0.0432828    0.0871608    0.0645703    0.0175975
 0.0871608    0.477114     0.960789     0.71177      0.19398
 0.0175975    0.0963276    0.19398      0.143704     0.0391639
 6.50731f-5   0.000356206  0.000717312  0.000531398  0.000144823

separable_view{N}(fct, sz, args...; pos=zero(real(eltype(AT))), offset =  sz.÷2 .+1, scale = one(real(eltype(AT))), operation = .*)

creates an array view of an N-dimensional separable function. Note that this view consumes much less memory than a full allocation of the collected result. Note also that an N-dimensional calculation expression may be much slower than this view reprentation of a product of N one-dimensional arrays. See the example below.

Arguments:

• fct: The separable function, with a number of arguments corresponding to length(args), each corresponding to a vector of separable inputs of length N. The first argument of this function is a Tuple corresponding the centered indices.
• sz: The size of the N-dimensional array to create
• args...: a list of arguments, each being an N-dimensional vector
• dims: a vector [] of valid dimensions. Only these dimension will be calculated but they are oriented in ND.
• offset: position of the center from which the position is measured
• scale: defines the pixel size as vector or scalar. Default: 1.0.
• operation: the separable operation connecting the separable dimensions

Example:

julia> fct = (r, sz, pos, sigma)-> exp(-(r-pos)^2/(2*sigma^2))
julia> my_gaussian = separable_view(fct, (6,5), (0.1,0.2), (0.5,1.0))
(6×1 Matrix{Float32}) .* (1×5 Matrix{Float32}):
 3.99823e-10  2.18861e-9   4.40732e-9   3.26502e-9   8.89822e-10
 1.3138e-5    7.19168e-5   0.000144823  0.000107287  2.92392e-5
 0.00790705   0.0432828    0.0871609    0.0645703    0.0175975
 0.0871609    0.477114     0.960789     0.71177      0.19398
 0.0175975    0.0963276    0.19398      0.143704     0.0391639
 6.50731e-5   0.000356206  0.000717312  0.000531398  0.000144823

separable_create([::Type{TA},] fct, sz::NTuple{N, Int}, args...; pos=zero(real(eltype(AT))), offset =  sz.÷2 .+1, scale = one(real(eltype(AT))), operation = *, kwargs...) where {TA, N}

creates an array view of an N-dimensional separable function including memory allocation and collection. See the example below.

Arguments:

• fct: The separable function, with a number of arguments corresponding to length(args), each corresponding to a vector of separable inputs of length N. The first argument of this function is a Tuple corresponding the centered indices.
• sz: The size of the N-dimensional array to create
• args...: a list of arguments, each being an N-dimensional vector
• dims: a vector [] of valid dimensions. Only these dimension will be calculated but they are oriented in ND.
• offset: position of the center from which the position is measured
• scale: defines the pixel size as vector or scalar. Default: 1.0.
• operation: the separable operation connecting the separable dimensions

Example:

julia> fct = (r, sz, sigma)-> exp(-r^2/(2*sigma^2))
julia> my_gaussian = separable_create(fct, (6,5), (0.5,1.0); pos=(0.1,0.2))
6×5 Matrix{Float32}:
 3.99823f-10  2.18861f-9   4.40732f-9   3.26502f-9   8.89822f-10
 1.3138f-5    7.19168f-5   0.000144823  0.000107287  2.92392f-5
 0.00790705   0.0432828    0.0871608    0.0645703    0.0175975
 0.0871608    0.477114     0.960789     0.71177      0.19398
 0.0175975    0.0963276    0.19398      0.143704     0.0391639
 6.50731f-5   0.000356206  0.000717312  0.000531398  0.000144823

copy_corners!(arr::AbstractArray{T,N}) where {T,N}

replicates the first N-dimensional quadrant my several mirror operations over the entire array. The overwrites the entire array with the (mirrored) content of the first quadrant!

Arguments

• arr: The array in which the copy operations are performed
• speedup_last_dim=true: if true a one-dimensional assignment trick is used for the last non-singleton dimension.