stev47 / StaticKernels.jl

Julia-native non-allocating kernel operations on arrays
Other
25 stars 0 forks source link

Feature Request: More granular control over boundary conditions #3

Open AlfTetzlaff opened 4 years ago

AlfTetzlaff commented 4 years ago

Hi, I don't know how one could implement this efficiently, but it would be cool to have more granular control over boundary conditions. Imagine for example a cube with periodic boundary conditions in two directions and two different constant values e.g. at top and bottom.

stev47 commented 4 years ago

It is currently still undocumented, but you have full control over the boundary conditions already. A quick way to do that is defining your own extension type (here mostly copied from ExtensionCircular in src/extension.jl):

using StaticKernels
struct ExtensionCustom{T} <: StaticKernels.Extension
    a::T
    b::T
end
Base.@propagate_inbounds @inline function StaticKernels.getindex_extension(w, wi, ext::ExtensionCustom)
    wi[1] < 0 && return ext.a
    wi[1] > 0 && return ext.b
    pi = position(w) + wi
    pimod = CartesianIndex(mod1.(Tuple(pi), size(parent(w))))
    return parent(w)[pimod]
end
a = rand(3, 3, 3)
k = Kernel{(-1:1,-1:1,-1:1)}(w -> sum(Tuple(w)), ExtensionCustom(5, -5))
map(k, a)

Note that there is currently high compilation time for higher dimensional kernels with boundaries.

This implementation checks the index wi[1]. You could check on axes_inner(w) instead, since that is known at compile-time to get better performance but you probably won't notice much of a difference since the interior dominates very quickly and there the extension is not called anyway.

AlfTetzlaff commented 4 years ago

Thanks for the quick reply! I'll need some time to think into it and will give it a try. I leave the decision to close this to you, as open issues have better visibility if someone else looks for the same thing (at least if there are only a few). Btw: naming a variable "pi" looks strange to me.