Simulation

MR-zero provides two simulations:

• a purposefully simple isochromat based simulation to check for correctness
• a state-of-the-art PDG simulation, built for speed and accuracy.

Both simulations are differentiable and capable to run on any device supported by PyTorch (all CPUs and NVIDIA GPUs). While the isochromat simulation is rather slow and only there for confirmation purposes, PDG is one of the fastest currently available simulations. You can find more information about it in the Paper.

It is strongly adviced to use the PDG simulation - isochromat_sim only exists for testing purposes.

Phase Distribution Graph (PDG) simulation

Phase Distribution Graphs split the magnetization into multiple states, forming a Graph over the duration of the sequence. More details can be found in the PDG Paper, once published.

The simulation runs in two passes:

• The pre-pass (compute_graph) simulates a large number of states, but only roughly estimates the signal. This information is used to assess which states are important for the signal, and which can be skipped because they don’t contribute anything.
• The main-pass (execute_graph) uses the exact signal equation and takes all voxels into account, but simulates as few states as possible. It does so by using the information generated by the pre-pass.

The precision / speed tradeoff for both passes can be chosen as follows:

# Default parameters: max 200 states per repetition
# with at leat 0.0001 of maximum possible magnetization
graph = mr0.compute_graph(seq, data, 200, 1e-4)
# Minimum emitted and latent signal of 0.01 (normalized)
signal = mr0.execute_graph(graph, seq, data, 0.01, 0.01)

The pre-pass limits the state count by both the maximum number (higher = more states / slower + more precise simulaiton) of states and the minimum magnetization (lower = more states). In practice, it turns out that the second parameter is less surprising in how it influences the quality of the simulation, which means that the first parameter should be set to a high number and the second parameter should be increased if the pre-pass is too slow (usually the main-pass is more likely to be the bottleneck).

The main-pass has two parameters which are used as threshold to determine which states are used in the signal equation and which states are simulated. Usually they are both set to the same number. Lower numbers mean that more states are above the thresold and will be simulated. It is not useful to set min_emitted_signal to a value lower than min_latent_signal, as this is trying to include to states into the signal equation that are not simulated.

Examples:

# Very precise simulation
graph = mr0.compute_graph(seq, data, 10000, 1e-7)
signal = mr0.execute_graph(graph, seq, data, 1e-7, 1e-7)

# Average precision
graph = mr0.compute_graph(seq, data, 10000, 1e-4)
signal = mr0.execute_graph(graph, seq, data, 1e-3, 1e-4)

# FID only
graph = mr0.compute_graph(seq, data, 10000, 1e-3)
signal = mr0.execute_graph(graph, seq, data, 1, 0.1)

Optimal settings can always depend on the sequence itself. Some sequences are more demanding and require a very precise simulation, some are already decently described by only a few states.

compute_graph

This utility function computes the parameters for compute_graph_ext from the passed data. Returns the computed Graph.

compute_graph_ext

Compute the Phase Distribution Graph for the given sequence and settings. Runs a simulation of seq for a single voxel with the given properties. Discards states with magnetization below min_state_mag and the states with the lowest magnetization if surpassing max_state_count. Returns the computed Graph.

execute_graph

Takes the Sequence, the SimData and the Graph and calculates the measured ADC signal. Because of the work done by the pre-pass, only the minimal work needed in order to achieve the desired precision is executed. This precision can be tuned by the min_emitted_signal and min_latent_signal thresholds. Higher values lead to less states being simulated, which improves speed and reduces accuracy. A value of 1 will mean that 2 states will be simulated (z0 and one +), resulting in the FID signal. A value of 0 means that everything will be simulated that somehow contributes to the signal, even if vanishingly small.

Graph

A wrapper for the retuned graph, which provides two additional methods:

Graph.plot

Creates a matplotlib scatter plot of all PDG states. Does not create or show the figure, which allows to modify it, for example:

graph = compute_graph(seq, data)

plt.figure()
graph.plot(dephasing="k_x")
plt.grid()
plt.show()

Graph.plot_data

Returns the data used in Graph.plot.

isochromat_sim

Isochromat based Bloch simulation. The code is implemented for readability, not speed. This allows to verify its correctness, bugs are not expected.

Warning

This simulation

• is slow
• does not support diffusion
• always uses "box" - voxels, which is usually a bad choice and will introduce slight blurring. See Voxel Shape for more information. This choice was made for simplicity.