Spreadinterponly (CPU and GPU)

[ahbarnett]

This is a proposal for how a convenient access to the spread/interp task is possible via the FINUFFT/cuFINUFFT interface. (since spreadinterp module is not part of our API). It is based on @chaithyagr PR #599 for the new CPU opts field. And the GPU s/i-only is already part of 2.3.1 and master.

Design discussion:

I feel that the GPU spreadinterponly interface has to be tweaked so that upsampfac controls the kernel shape. This would also apply to CPU. Currently GPU forces upsampfac=1.0, but under the hood uses the kernel for the default upsampfac=2, which makes no sense.
Eg, I would like to be able to set upsampfac=infinity, which gives a certain kernel shape needed in Ewald methods.

I will hash this out here. Comments welcome.

Tasks:

• clean up CPU spreadinterponly=1 logic
• add the new opts field to all language interfaces and check
• example which demos it (1D)
• add tester to CI (1D)
• docs CPU
• MATLAB demo CPU
• Revisit gpu_spreadinterponly=1 and edit logic so upsampfac is respected rather than forced to be 1.0
• GPU doc copying CPU doc

[ahbarnett]

@chaithyagr If you look at the example, tester, and docs/opts.rst you will see my proposal for the cpu spreadinterponly=1 implementation. See what you think. If you are not unhappy - let me know - I will apply this to GPU too - the difference being you would not have to set upsampfac=1 to use it, rather would set it as with a NUFFT. This is really the only way I can see it working, without breaking the API.

I still would be curious how you guys access the actual kernel function used, in your DCW for MRI application...

[ahbarnett]

@chaithyagr Still hoping for some feedback on GPU interface tweak, so I can go ahead (early February at this rate) with:

Revisit gpu_spreadinterponly=1 and edit logic so upsampfac is respected rather than forced to be 1.0.

[chaithyagr]

Hey @ahbarnett , I am so sorry for my delay. I think I responded your question regarding ns going to infinity in #564 .
Coming to this update,
I am okay with this new interface, I will have to change some stuff in mind-inria/mri-nufft#195 and mind-inria/mri-nufft#224.
Just some quick questions to be sure I understood correcty:

1. So right now if spreadinterp_only is set to 1, the output is still with upsampfac>1, basically we will need to send arrays of the right sizes on our own. This was not needed when we do NUFFT as the output shape is exactly the same size as what we need.

To be clearer:
If we are doing NUFFT with 1000 points trajectory to a shape of (256, 256) and upsampfac=2,
then the image size in forward NUFFT is (256, 256), but if spreadinterp_only is set, it is (512, 512), right?

Same applies for adjoint operation where we provide a memory location for output image dimensions.

2. Do we have checks on input shapes of the images? I didnt find this originally and I just wanted to make sure that the sizes to API remain the same, which is why I chose the route of making upsampfac=1

Finally, please do let me know when you are ready and I would like to make sure my tests pass and I can make fixes right on my side.

Thank you so much for handeling this and please do let me know when you are ready. I will do my changes on finufft side at mri-nufft today and give you update on whether it works good.

[chaithyagr]

Assuming that we need to provide (512, 512), note that current state of codes lead to following error:

finufft/python/finufft/finufft/_interfaces.py

Lines 467 to 474 in d13dade

	if fshape[-1] != ms:
	raise RuntimeError('FINUFFT f.shape is not consistent with n_modes')
	if dim>1:
	if fshape[-2] != mt:
	raise RuntimeError('FINUFFT f.shape is not consistent with n_modes')
	if dim>2:
	if fshape[-3] != mu:
	raise RuntimeError('FINUFFT f.shape is not consistent with n_modes')

Basically the validity checks are right, except when we are pushing out a larger array when we have spreadinterp_only is true.

[chaithyagr]

@ahbarnett did you get a chance to look at my comments? The main issue of setting spreadinterp_only=1 and upsampfac!=1, implies we need to update the API interface such that it can support having image domain shapes = img_size * upsampfac.