Add periodic FlowPath magnetization remapping

[mcgrathcm]

Summary

Adds GPU magnetization remapping for periodic, non-closed FlowPath motion.

• Accepts a precomputed destination-to-source cycle_map
• Splits simulation blocks at exact cycle boundaries
• Applies the map with a standalone GPU gather kernel and reusable scratch buffer
• Remaps both transverse and longitudinal magnetization
• Preserves existing inflow reset behavior
• Supports phantom serialization with and without a map
• Disallows phantom slicing when a map is present
• Reports a clear error for CPU simulation

Koma does not calculate the mapping; it must be supplied with the phantom.

Testing

• KomaMRIBase: 775 tests passed
• KomaMRICore CPU: 437 tests passed
• KomaMRIFiles: 24,450 tests passed
• GPU behavior validated with a 2-million-spin periodic-flow simulation

Closes #930

[codecov]

Codecov Report

❌ Patch coverage is 80.34188% with 23 lines in your changes missing coverage. Please review.
✅ Project coverage is 88.96%. Comparing base (2cf6c8e) to head (5dd7fff).
⚠️ Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
KomaMRICore/src/simulation/Flow.jl	54.83%	14 Missing ⚠️
KomaMRIBase/src/datatypes/Phantom.jl	63.63%	4 Missing ⚠️
...se/src/motion/actions/arbitraryactions/FlowPath.jl	82.35%	3 Missing ⚠️
KomaMRIBase/src/motion/MotionList.jl	77.77%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #935      +/-   ##
==========================================
- Coverage   91.38%   88.96%   -2.43%     
==========================================
  Files         137      137              
  Lines        7255     7293      +38     
==========================================
- Hits         6630     6488     -142     
- Misses        625      805     +180     

Flag	Coverage Δ
base	90.66% <86.15%> (-0.10%)	⬇️
core	75.09% <70.21%> (-13.84%)	⬇️
files	95.23% <100.00%> (+0.07%)	⬆️
komamri	90.51% <ø> (ø)
plots	91.76% <ø> (ø)

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing lines	Coverage Δ
KomaMRIBase/src/motion/Motion.jl	100.00% <100.00%> (ø)
KomaMRIBase/src/motion/NoMotion.jl	92.30% <100.00%> (+0.30%)	⬆️
KomaMRICore/src/simulation/Functors.jl	72.00% <ø> (-12.00%)	⬇️
KomaMRICore/src/simulation/SimulatorCore.jl	89.13% <100.00%> (-0.26%)	⬇️
KomaMRIFiles/src/Phantom/Phantom.jl	97.50% <100.00%> (+0.91%)	⬆️
KomaMRIBase/src/motion/MotionList.jl	98.16% <77.77%> (-1.84%)	⬇️
...se/src/motion/actions/arbitraryactions/FlowPath.jl	85.00% <82.35%> (-15.00%)	⬇️
KomaMRIBase/src/datatypes/Phantom.jl	84.80% <63.63%> (-1.21%)	⬇️
KomaMRICore/src/simulation/Flow.jl	71.79% <54.83%> (-11.19%)	⬇️

... and 15 files with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.
• 📦 JS Bundle Analysis: Save yourself from yourself by tracking and limiting bundle sizes in JS merges.

[mcgrathcm]

Here is an example of the issue being solved for a spoiled GRE sequence. The artifact presents itself as a flash at the start the of the cycle, most noticeable at the outlet region. NN remapping performs nearest neighbor remapping of magnetization. The magnitude of this issue is very sequence dependent, so it may be more or less present depending on sequence. @cncastillo

[mcgrathcm]

I have come across an issue that needs some investigation here.

Currently this implementation requires block splitting at periodic boundaries. This is fine for precession blocks, however for excitation blocks this requires splitting the RF pulse. I am not entirely sure how this will affect the simulation accuracy, especially for more advanced methods like BlochMagnus.

Unfortunately I believe this is required, as this requires communication of magnetization between spins, which conflicts directly with the kernel architecture which is one spin per thread.

One upside, even if this somewhat breaks the calculation of the RF pulse, it likely will not break it as much as the current implementation where particles can teleport across the entire phantom, resulting in crazy state mixing.

Basically, the end result would mean that sometimes a single RF pulse may split into two excitation blocks with non-zero start/end strength. And in-between them there is a re-mixing of magnetization.

@cncastillo @pvillacorta curious what you think.

[pvillacorta]

Hi! Sorry for the delay in replying. I think this RF-splitting is worthwhile, as long as this behavior in the simulation is isolated to the FlowPath case.

I still need to fully understand this specific issue and see how it fits into some of the other issues we have already addressed. If I understand correctly, in #620 we identified two related problems involving periodic motion and flow:

1. (Blue dot in the diagram in #620) Both when spin remapping occurred and when a new motion period started, the magnetization of the affected spins was reset to their initial magnetization (M0). However, there were cases where, when discretizing the sequence, no key point was created within the reset window (while the spin_reset flag was set to 1). This is why, in Add key time points in DiscreteSequence when evaluating periodic and flow-related motions #638, we forced additional points to be added to the DiscreteSequence near these boundaries.

2. (Add GPU magnetization remapping at periodic motion boundaries #930, yellow dot in the diagram in #620) There are cases where the spin magnetization after being "reset" should not be equal to the initial M0, but should instead start from an intermediate state between 0 and M0 (e.g. M_steady_state). To solve this, we would need to be able to set the spin magnetization, after the reset, to an arbitrary value. We had not solved this problem until now.
   I have a couple of clarifying questions:

   1. What you are proposing is to set this arbitrary value to that of the nearest neighbor when each motion period ends and a new one begins, correct?
   2. Is this only needed when a new period starts, or is it also needed when a spin leaves an outlet and enters an inlet (spin_reset=1)? I find it difficult to distinguish between these two situations. From what I can see in the GIFs, NN remapping is only needed in the first case, right? If there is a reset within a period, could we assume that the initial state of that spin should go back to M0, since it is effectively "coming from outside" the excited slice?

3. There is also this other related issue, Outflow Spin Reset at ADC points #625, which is still unresolved.

I want to reiterate that I think this solution is worthwhile because it simulates the behavior much more accurately.

Just to clarify one more point:

One upside, even if this somewhat breaks the calculation of the RF pulse, it likely will not break it as much as the current implementation where particles can teleport across the entire phantom, resulting in crazy state mixing.

The particles are still going to teleport, right? It's just that their magnetization will be able to acquire more accurate values when the period is reset. So, just to make sure I understand the proposed solution correctly.

Thanks!

[mcgrathcm]

Yes, this directly addresses #620.

1. For adding time points at spin reset events, Add key time points in DiscreteSequence when evaluating periodic and flow-related motions #638 seems to be working well and I do not see any issues. Interestingly for the large scale cases we are using, this ends up adding a simulation time point at every motion point, as there is constant resetting.
2. It is basically what you have stated in ii, this PR is specifically addressing the periodic wrapping of the entire FlowPath, not the intermediate usage of spin_reset. What you are pointing to here is a real issue, but solving this is non-trivial, as we would need to determine the correct steady state signal (including phase if we want to set Mxy). I think this is out of scope for the simulation package. Maybe in the future we can allow the user to set the "steady state" M, but the responsibility to determine the correct value would be on them.
3. I have somewhat forgotten about this, but will check it out again and see if I can incorporate a fix here.

And to your final comment, yes the spins will still teleport, and thanks to your changes in #638 the teleport is over MIN_RISE_TIME so negligible phase will be added.
The flow is:

• Particles reach final step in FlowPath
• Particles take a single step of MIN_RISE_TIME back to starting point, negligible change to magnetization
• Particle magnetization is re-mapped such that the ensemble magnetization structure of the flow is preserved.
• Continue normally

There are a few minor fixes I still need to add, but otherwise this seems to be functional on my end.