-
Notifications
You must be signed in to change notification settings - Fork 27
Allow convolve to work with Real arrays with arbitrary precision
#159
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Conversation
|
Usual windows dll failure, but otherwise passing. |
nHackel
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
From what I could see the changes didn't lead to any type instabilities and the compiler should be able to generate the same code as before.
If you update the Project.toml to with a version 0.14.2, I can make a release
src/convolution.jl
Outdated
| p.params.m, p.params.σ, scale, j, d, L, p.params.LUTSize) ) | ||
|
|
||
| @nexprs 1 d -> prodWin_{$D} = one(T) | ||
| T = promote_type(Tp, Tg) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Very minor nitpick: T = promote_type(Tp, Tg) is unused in this function
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thank you very much for the quick reivew and for pointing out the unused T.
Actually that review made me realize that I should be more careful with all of the input and output types and that led to further changes including tests to make sure the caller does not try to convolve a Complex input into a Real output.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
BTW, one could argue that Real is too permissive because it would even allow an array of Ints which clearly would not work for the output of convolve!.
I'm willing to modify to be AbstractFloat instead if you prefer.
Given that the previous version used duck typing for some of the arguments, I figured there is an expectation here that any user who tries to convolve with an inappropriate output type like Int will realize their mistake when they get a low-level error from the compiler.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think Real is fine, even if it's just based on 'historical precedent'. Since AbstractNFFTs does not specify an element type, we could argue that restricting it to abstract floats is not a breaking change. As long as we use abstract types (which we want to use) rather than concrete types, users may still encounter a case where a number type will break. Our tests and documentation (should) show the intended and usual cases.
I try to specify types when we want to use them for behaviour/dispatch. For example, you now do this with the promote_type, and in the case of the high-level convolve! function, we want to perform a rough check to ensure that the user isn't passing in any non-numeric elements. But at the point of the high-level function, we don't necessarily know whether any of the lower-level functions will break and one could argue for very relaxed typing here.
I think it's fine the way it is, but maybe we could add some tests showing that the convolve call breaks with incorrect inputs, such as complex input to real output. This would document the intended behaviour. Some of those tests might happen indirectly with your planned changes to the sdc anyway
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Adding tests makes sense. It took me a while to get the tests right because I was using a wrong array size and instead of failing with indexing errors, I was getting garbage numbers. So I also added size tests.
|
BTW, it was challenging to check with Float16 precision because |
src/convolution.jl
Outdated
| throw(DimensionMismatch("size(g)=$(size(g)) ≠ Ñ = $(p.Ñ)")) | ||
| size(fHat) == (p.J,) || | ||
| throw(DimensionMismatch("size(fHat)=$(size(fHat)) ≠ J = $(p.J)")) | ||
| (eltype(g) <: Complex) && (eltype(fHat) <: Real) && |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Tiny nitpick, which you can ignore:
We can determine the wrong eltype via multiple-dispatch by adding a:
function AbstractNFFTs.convolve!( # and similar for convolve_transpose!
::NFFTPlan{<:Real,D,1},
::AbstractArray{<:Complex, D},
::StridedVector{<:Real}) where {D}
throw(ArgumentError("Complex input g requires Complex output fHat"))
endIn a small test, it seemed like the compiler couldn't automagically remove the eltype check even though it does have the necessary information. But even in that case, there shouldn't be any relevant impact on performance
|
Happy new year @JeffFessler, thanks a lot for the added tests! I added a small comment which you can just ignore and if you ping me I'll merge and release this PR. I'm not deep enough into the internals to understand the interpolations issue, but your solution around that looks good. At the moment NFFT.jl is only working for FFTW.jl. If there ever is an application which requires Float16 values and if it's possible to exchange FFT backends like we can do it with AbstractNFFT backends, then we might be able to fully support that via FFTA.jl further down the line. |
|
Happy New Year to you too! I like the multiple dispatch approach to the Complex/Real mismatch tests. |
|
The CI passes now (other than the usual dll issue) so I am merging. Thanks again! |
Currently,
convolve!and its adjoint are (partially) constrained to work withComplexvectors.This constraint is natural for MRI reconstruction, but it is unnatural and inconvenient when designing sampling density compensation functions, because density factors should be
Real(and positive).This PR makes the following changes. The first two are the big ones.
These are all backward-compatible.
Unionto allow input arrays to be eitherRealorComplex.This generality should have no effect in the usual cases where the input array precisions match that of the
Plan, while providing callers more flexibility if needed, with nearly negligible extra coding. I'll note that the current version requires one of the two input arrays to have a precision that matches the plan, while the other input array uses "duck typing" with no precision constraints at all (in fact not even a constraint that the input array haveNumberelements). This PR provides somewhat more equity in the element type freedom of the two input arrays.Edit: the updated version removes the duck typing at the highest level of the interface and adds checks to make sure that a callee does not try to convolve a Complex input into a Real output.
convolve!and its adjoint return the mutated array, which is typical for mutating functions.ininstead of=withforloops; previous code used a mix of both styles.If approved, I will ask for a version bump so I can use the new features to improve
sdc.