Infer arraytype for plan_nfft
from k
#142
Open
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
This PR (hopefully) closes #100. The goal of this PR is to seamlessly create either CPU or GPU plans depending on the given argument(s). At the bottom of this comment I'll list some limitations of the current situation of NFFT.jl in regards to this switching.
The basic problem of correctly inferring the type is that we need the "base" type of a given array. For example if given a
Array{ComplexF32, 2}
we wantArray
and if givenCuArray{Float32, 2, ...}
we wantCuArray
. We require this base type, since it is needed forAdapt.jl
to correctly move temporary plan fields to the GPU. At the moment there is no stable API in Julia to achieve this, so we have to use some workaround.A first option would be to introduce a new function. Since there is no suitable base function we could extend for this case (as far as I know) we have to define our own, for example
plan_type(arr::AbstractArray) = Array
. Then with package extensions we can implementplan_type(::CuArray{...} = CuArray
and so on. Any array type we don't define, would have to be defined by a user and/or a new PR needs to be made.The second option is to try to determine the "base" type of a given array using a (not stable) API as described in here. The benefit is that this also works without us specifying types for the different GPU backends, but it can break in future Julia versions. Potentially, we would need to support different variants of this function for different Julia versions.
Both approaches are just "best-effort" approaches and there will likely be cases where they fall. For example if
k
is not a dense array but some lazy-wrapped array, for exampleadjoint(k)
. In those edge-cases a user still needs to be able to create a fitting plan. For the first option a user could implement theirplan_type
, but since theplan_nfft(Array, k, ...; ...)
interface is still available, I feel like this is superfluous and opted to implement the second option.For convenience I also implemented the inferrence for
adjoint
andtranspose
. If there are any more common wrappers we want to support we can add them as package extensions.With this PR we can chose the correct plan_type based on k. However, at the moment the
initParams
functions in whichk
is used only accepts dense CPU matrices. This means that for the GPU backends, we first infer the type correctly fromk
, then movek
to the CPU and construct our GPU plan. I think this limitation exceeds the scope of this PR, but it is a little bit counter-intuitive from a user perspective.Similarly neither our
initParams
nor the inferrence will work with a matrix-free construct like from LinearOperators.jl or LinearMappings.jl.