Add key time points in DiscreteSequence when evaluating periodic and flow-related motions

[pvillacorta]

Try with:

add https://github.com/JuliaHealth/KomaMRI.jl#key-time-points:KomaMRIBase
add https://github.com/JuliaHealth/KomaMRI.jl#key-time-points:KomaMRICore
add https://github.com/JuliaHealth/KomaMRI.jl#key-time-points:KomaMRIFiles

From #620 (comment):

We have the following implementation of the add_jump_times! function for FlowPath actions:

function add_jump_times!(t, a::FlowPath, tc::TimeCurve) 
    jump_times = (tc.t_end - tc.t_start)/(size(a.spin_reset)[2]-1) * (getindex.(findall(a.spin_reset .== 1), 2) .- 1) .- 1e-6
    append!(t, jump_times)
end

which adds to the simulation the time points at which a spin reset is produced. The substraction of 1us (.- 1e-6) to these time points ensures that we are evaluating within the time interval in which interpolation will return 1 (related #511). This way, we do not miss any spin reset.

However, this add_jump_times! function certainly presents some problems. I will put them as checkboxes so I can address them in a PR:

• Until now, I hadn't realized that, if periodic == true, no time points are added for next periods. I must extend the function with push!(jump_times, jump_times .+ n*period) or something similar.
• Since a big "jump" or spatial discontinuity may be produced when spins re-appear at their starting point at the beginning of each period, we might well ensure that no phase is accumulated during a simulation time step (Δt in the drawing) due to this position change. To achieve this, we must include in jump_times two more time points (as depicted in blue in the drawing): tf and tf+eps.

This is the PR intended to solve the above. I have created and renamed some functions to obtain the following structure:

• add_key_time_points!: this function is called from get_variable_times and adds all the necessary time points to the sequence when motion is present. These key points can be divided into two categories: period-related and flow-related, which are obtained through two different and consecutive function calls from add_key_time_points!:

  • add_period_times!: adds to the time vector all the key points related to pseudo-periodic motions (i.e., motions where the periods field is a vector). This is necessary for all motion types, not only FlowPath actions.
  • add_reset_times!: this function is equivalent to the old add_jump_times!, but it has been improved to support pseudo-periodicity.

  Finally, when periodic == true, all key time points obtained from both add_period_times! and add_reset_times! are periodically extended.

This PR does not fully resolve issue #620 but provides a partial fix. Based on the figure in this comment, it addresses the point marked with the blue bullet.

This change might slightly impact performance in dynamic simulations, since we now check periodicity and add key time points for all motion types, not just for FlowPath motions. I realized that any motion (not only FlowPath) can lead to discontinuities when extended periodically, which is why these checks are now always performed (except, of course, in the NoMotion case).

[pvillacorta]

This should be ready.
Buildkite tests might be re-run, since the following error:

is not related with the code, and might be solved by simply re-running.
Documentation is failing due to a shutdown signal during building. Maybe a timeout error?

[cncastillo]

I am re-trying the failed jobs

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 91.22%. Comparing base (aadb879) to head (bb6042e).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #638      +/-   ##
==========================================
+ Coverage   90.14%   91.22%   +1.07%     
==========================================
  Files          57       57              
  Lines        3228     3248      +20     
==========================================
+ Hits         2910     2963      +53     
+ Misses        318      285      -33     

Flag	Coverage Δ
base	91.48% <100.00%> (+2.40%)	⬆️
core	89.73% <100.00%> (ø)
files	94.98% <100.00%> (+0.92%)	⬆️
komamri	88.13% <100.00%> (-0.08%)	⬇️
plots	90.90% <ø> (ø)

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

Files with missing lines	Coverage Δ
KomaMRIBase/src/KomaMRIBase.jl	100.00% <ø> (ø)
KomaMRIBase/src/datatypes/Phantom.jl	86.01% <100.00%> (+0.14%)	⬆️
KomaMRIBase/src/motion/Action.jl	100.00% <100.00%> (+25.00%)	⬆️
KomaMRIBase/src/motion/Motion.jl	100.00% <100.00%> (+7.69%)	⬆️
KomaMRIBase/src/motion/MotionList.jl	100.00% <100.00%> (+7.00%)	⬆️
KomaMRIBase/src/motion/NoMotion.jl	96.00% <100.00%> (ø)
KomaMRIBase/src/motion/TimeCurve.jl	95.83% <100.00%> (+35.83%)	⬆️
KomaMRIBase/src/motion/actions/ArbitraryAction.jl	75.00% <ø> (+8.33%)	⬆️
...se/src/motion/actions/arbitraryactions/FlowPath.jl	100.00% <100.00%> (+75.00%)	⬆️
...MRIBase/src/motion/actions/simpleactions/Rotate.jl	100.00% <100.00%> (+7.69%)	⬆️
... and 4 more

🚀 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.

[pvillacorta]

Ok, I have to review code coverage

[pvillacorta]

The way I tested the included functions is by manually determining key points from a specified "TimeCurve":

t        = [0.0, 0.1, 0.3]
t_unit   = [0.0, 0.4, 1.0]
periods  = [1.0, 0.5, 2.0]
periodic = true
dx = dy  = [0.0 0.0 0.0 0.0]
dz       = [3.0 4.0 -4. -3.]
reset    = [false false true false]
Δ        = 1e-6

fpth = flowpath(dx, dy, dz, reset, TimeCurve(t, t_unit, periodic, periods), AllSpins())

period_times = reduce(vcat, [t .+ [-Δ, Δ] for t in [0.0, 0.3, 0.45, 1.05, 1.35, 1.5]])
reset_times  = [0.2, 0.4, 0.85, 1.25, 1.45] .- Δ

[github-actions]

KomaMRI Benchmarks

Benchmark suite	Current: e8d5ca2	Previous: ee3bc2f	Ratio
MRI Lab/Bloch/CPU/2 thread(s)	339465202.5 ns	341921413 ns	0.99
MRI Lab/Bloch/CPU/4 thread(s)	272873331 ns	275323882 ns	0.99
MRI Lab/Bloch/CPU/8 thread(s)	205654119 ns	285513733 ns	0.72
MRI Lab/Bloch/CPU/1 thread(s)	556512232.5 ns	552505595.5 ns	1.01
MRI Lab/Bloch/GPU/CUDA	21615042.5 ns	21193443 ns	1.02
MRI Lab/Bloch/GPU/oneAPI	76004686 ns	76632822 ns	0.99
MRI Lab/Bloch/GPU/Metal	96513083 ns	96112375 ns	1.00
MRI Lab/Bloch/GPU/AMDGPU	24666705 ns	28621139 ns	0.86
Slice Selection 3D/Bloch/CPU/2 thread(s)	1598294292 ns	1622705599 ns	0.98
Slice Selection 3D/Bloch/CPU/4 thread(s)	885748828 ns	891208860 ns	0.99
Slice Selection 3D/Bloch/CPU/8 thread(s)	556034306 ns	571183580.5 ns	0.97
Slice Selection 3D/Bloch/CPU/1 thread(s)	3040144004 ns	3037417759.5 ns	1.00
Slice Selection 3D/Bloch/GPU/CUDA	32496243 ns	32085949 ns	1.01
Slice Selection 3D/Bloch/GPU/oneAPI	121113069.5 ns	121193856.5 ns	1.00
Slice Selection 3D/Bloch/GPU/Metal	110020292 ns	112065187.5 ns	0.98
Slice Selection 3D/Bloch/GPU/AMDGPU	32605152 ns	36595885 ns	0.89

This comment was automatically generated by workflow using github-action-benchmark.

[cncastillo]

Pleases remove all of this "sampling_params["Δt_rise"]" and use MIN_RISE_TIME if you need a constant for a small time.

[pvillacorta]

I tried to address all the requested changes. @cncastillo feel free to resolve the comments if you agree.

[cncastillo]

Re-running again the failed job, let me know if it passes to bump and register after merge.

[pvillacorta]

@cncastillo documentation is still failing, but only in CI/CD servers. When run locally, docs are correctly built. It seems as a timeout error or something similar. Is this happening in other PRs?

[cncastillo]

That is weird, I run it again with debug logs on, to see what is going on.

[pvillacorta]

I still don't get why docs are failing :-/

[cncastillo]

It says that we need to re-run it. https://github.com/JuliaHealth/KomaMRI.jl/actions/runs/19511907808/job/55858775999

[cncastillo]

I can't figure out it maybe add a timeout, like this one:

https://github.com/JuliaHealth/KomaMRI.jl/blob/master/.github%2Fworkflows%2FCI.yml#L15

[pvillacorta]

It turns out this line in add_key_time_points! (only executed with periodic motions):

for n in 1:floor(t_max/period) append!(aux, aux .+ n*period) end

was absurdly inefficient. When the input t got large enough, the function would actually break (which is why the docs failed for this case). Just looking at the line and the append! call, you can see that on every iteration aux kept getting bigger and bigger, making each successive append! more expensive than the previous one (not to mention that we were adding a huge amount of duplicated points). In short, that line was a complete mess.

I fixed this in 8695493, using method dispatch, the sizehint! function, and fixing the append! call:

sizehint!(aux, initial_size * (n_periods + 1))
for n in 1:n_periods
    append!(aux, aux[1:initial_size] .+ n*period)
end

[pvillacorta]

This should be ready to merge! It might be good to check that this improvement actually meets its goal of providing accurate simulations when periodic flow motions are present. @mcgrathcm do you have any experiment or scenario where this is easy to verify? It would be helpful to at least see the differences between the results before and after this PR (even if we don’t include it in the tests).

[cncastillo]

It looks good to me, can you point to the issue this PR is solving so @mcgrathcm can remember what the problem was and potentially test it?

[pvillacorta]

It is #620

[cncastillo]

@mcgrathcm is finishing testing, in the meantime, can you bump versions and tell me wish packages to register?

[mcgrathcm]

I ran these change with one of the encodes of my 4DFlow sequence, and see no noticeable difference, which is a good thing in this case.
For this simulation, magnetization is numerically spoiled (Mxy set to zero) prior to each RF pulse. As such we do not expect to see long spin histories. Additionally, for this case the reseeding can only really affect the inlet region magnitude, but this would only be a small change in signal.

Overall, I cannot say whether this change fixes what it was set out to fix (as I don't have an example demonstrating the problem), but I can say it did not break anything that already worked!

[pvillacorta]

Very nice example and simulations :) And thank you for testing this!
I think that, even merging this PR, we might not want to close #620 yet, since there is still a pending change for the future:

Ideally, for flow simulation, there should be a step that remaps magnetization, such that if I took a snapshot of all particle before and after the end of the period, the overall magnetization distribution should looks the same, but the actual particle positions should reflect either the last or the first timepoint.

(PD: I have just noticed that these test images are gifs :O)

[cncastillo]

It seems to me that some changes in base are breaking for core and files (like the use of CenterOfMass and removal of add_jump_times!), so i bumped the compat versions. As this is not breaking for the user API it is fine as a minor bump. Do you agree?

[cncastillo]

Just check that the versions of the files are above master and the compats make sense as I think I registered a minor version for komafiles.

[pvillacorta]

Just check that the versions of the files are above master and the compats make sense as I think I registered a minor version for komafiles.

KomaFiles version is set to 0.9.7 in master. This PR sets it to 0.9.8