SeparableFunctions.jl
Documentation | Build Status | Code Coverage |
---|---|---|
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:
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. Yet, these lazy versions are currently NOT recommended, since they are consistently slower that the separabel implementations. This is why the specific versions are currently also not exported.
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 tooffset
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().\div2 .+ 1
is chosen. -
scale
: multiplies the axes with these factors. This can be interpreted as apixelsize
. Some functions have additional named arguments. E.g.gaussian
has the additional named argumentsigma
specifying the width of the gaussian, even though this can also be controlled viascale
. Note that this nomenclature is in large parts identical to the packageIndexFunArrays.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)
.