[![Build Status](https://github.com/Miha Fuderer/BLAKJac.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/Miha Fuderer/BLAKJac.jl/actions/workflows/CI.yml?query=branch%3Amain)

BLAKJac is implemented in the Julia language. It relates to the MR-STAT functionality and aims to predict or minimize the noise levels in the reconstructed quantitative maps.

Two of its main functions are

- Analyzer: given a sequence of encoding values and of RF pulse angles, it predicts the resulting noise level.
- Optimizer: given a sequence of encoding values, it generates a sequence of RF pulse angles that is presumed to be optimal in terms of noise level. The Optimizer makes use of the Analyzer.

In essence, the analyzer has as its input

- The RF sequence: a vector of
$N_{TR}$ complex values, with e.g.$N_{TR}=5\cdot 224=1120$ . The magnitude of the complex values expresses the flip angle in degrees. - The
*trajectorySet*. See dedicated section below. - A large bunch of options, in the form of a dictionary. For details, see BLAKJac_analysis and BLAKJac_defaults

There is one additional parameter:

- [Actually an output] Optionally, a dictionary of
$I^{-1}$ matrices, where$I$ is the diagonalized Fisher Information Matrix. Such a dictionary may be useful if the optimizer is run a multiplicity of times.

As an *Output*, the Analyzer returns

- A vector of noise levels of length
$p$ , with$p$ being the number of properties. For the time being,$p$ is always 3. - An estimate of the
*information content*that can be expected from the scan; for the time being, a pretty immature metric. -
*b1factorRMS*: an indication on how sensitive the sequence is to deviations of B1+.

Mostly, only the first output is used, the vector of noise levels. These have a somewhat quantitative meaning: these are scaled such that if

The interface, in full:

```
function BLAKJac_analysis!(
RFdeg::Vector{ComplexF64},
trajectorySet::Vector{<:Vector{<:TrajectoryElement}},
options::Dict,
saved_H::Dict=Dict())
```

In brief, the optimizer has a trajectorySet as input and it spits out an RF sequence. And it needs a lot of options.

The input:

- A
*trajectorySet*. See dedicated section below. - A large bunch of options, in the form of a dictionary.
- (Optionally) a random-value seed. Relevant if the option opt_initialize (see below) is set to cRandom30. By default, every run of the optimizer will always return the same pattern. If different realizations are needed, each call to the optimizer must be done with a different seed (Int).

Output

- It returns an RF sequence: a vector of
$N_{TR}$ complex values, with e.g.$N_{TR}=5\cdot 224=1120$ . The magnitude of the complex values expresses the flip angle in degrees.

The interface, in full:

function BLAKJac_optimize(trajectorySet, options::Dict, mersenneTwisterSeed=1)

The *trajectorySet*. This is a vector of *TrajectoryElement*. And a TrajectoryElement is a tuple, mostly used for

Note that in the 2D case,

The reason of making the trajectorySet a vector of vectors is for non-Cartesian: then, it can be e.g. a vector of 20 spirals of 1000 samples each and “$(k_y,k_z)$” has to be interpreted as

The function TrajectorySet(kyVector, kzVector) may be useful before invoking the Analyzer: it converts an

There are lots of option parameters, some of them are deprecated and some have defaults that are very rarely adapted. For that purpose, the following function is available:

BLAKJac_defaults!(trajectorySet::Vector{<:Vector{<:TrajectoryElement}}, options::Dict)

It does not return anything, but the “!” indicates that options will be modified by it. Note: it only assigns uninitialized options to a default value. Options that are already filled in will be left unmodified. To set all options to default, use

recon_options = Dict() # erase all existing settings

BLAKJac.BLAKJac_defaults!(trajectorySet, recon_options)

The *trajectorySet* is a parameter, since some defaults are set accordingly to the encoding sequence.

The list of options:

Name | Intended type | Default | Explanation |
---|---|---|---|

rfFile | string | (empty) | Functionally disused. Its only purpose is to tag figures (which may be generated inside the Analyzer) with the selected string. |

rflabel | string | (empty) | Functionally disused. Its only purpose is to tag figures (which may be generated inside the Analyzer) with the selected string. |

nTR | Int | Number of trajectories | Number of excitation pulses |

TR |
Float | 0.01 | TR in seconds |

T1ref | Float | 0.67 | The reference T1 value - by now, mainly used for noise scaling and for "relative" noise deviations |

T2ref | Float | 0.076 | The reference T2 value - by now, mainly used for noise scaling and for "relative" noise deviations |

startstate |
Float | 1.0 | Although it is formally a float, it acts like a Boolean: it indicates the starting state of the z-magnetisation; +1 for no prepulse, -1 for inversion |

maxstate | Int | 64 | Length of history taken into account in simulating magnetisation from sequence |

nky | Int | Derived from range of ky values in trajectorySet | Length of range of ky encoding values. |

nkz | Int | Derived from range of kz values in trajectorySet | Length of range of kz encoding values. So nTR/nky/nkz is number of samples per encoding. |

maxMeas | Int | Filled in as 4nTR/(nkynkz) |
(R1) (see remarks below) |

handleB1 | string | “no” | Available for testing the B1-sensitivity of sequences. Three valid values are "no", "sensitivity", "co-reconstruct". If “sensitivity” is selected, then the value of b1factorRMS is a meaningful output of the Analyzer: a value of 0 indicates that the T1 or T2 maps will not be affected by a deviation of B1; an output of e.g. -1.5 for e.g. T2 (which is typical) indicates that 1% of error in B1 will cause a 1.5% bias (underestimation) of T2. With “co-reconstruct”, B1 is just considered as a 4th parameter to be reconstructed (alongside ρ, T1 and T2) |

considerCyclic | Bool | false | If true, it is assumed that the RF sequence repeats itself endlessly without pause. Only compatible with startstate==1 |

useSymmetry | Bool | false | (Oops. Blooper in default setting: should have been made “true”); Assume real-ness of rho, T1 and T2 and therefor symmetry in k-space. |

invregval | Vector{Float} | [200.0^(-2), 200.0^(-2), 200.0^(-2)] | (R2) |

useSurrogate | Bool | false | (Relic. Don’t touch) |

sigma_ref | Float | 0.2 | (Irrelevant unless studying the information-content criterion, yet immature). Reference normalized noise level for information-content estimation. See logbook around 2021-06-03 |

sar_limit | Float | 40^2/0.01 | SAR limit. Initially set to an RMS level of 40 degrees at TR=10ms |

T1T2set | Array{(Float,Float)} | (R3) | A set of (e.g. 7) points in (T1,T2)-space around which the Analyzer will evaluate the performance. The output is the average over the performances at each point. The current default has been set to span a clinically relevant range of relaxation values. It has been set with 1.5T and 3T in mind. For very different fields (e.g. 7T), adaptation of the T1 values may be indicated. |

sizeSteps | Vector{Int} | [10] | (R4) |

portion_size | Int | 50 | (R4) |

portion_step | Int | 30 | (R4) |

opt_iterations_limit | Int | (R4) | (R4) |

optpars | Optim.Options | (R5) | (R5) |

opt_method | NelderMead | Optimizer algorithm. Do not change: most other choices require the explicit knowledge of the gradient of the criterion, or require a prohibitive number of iterations. | |

opt_expand_type | string | “spline” | Deprecated. Do not change |

opt_keep_positive | Bool | False | If true, resulting RF angles are not allowed to be negative. |

opt_initialize | string | “cRandom30” | Optimizer option: defines the starting pattern for the optimization process. Valid values are: "cRandom30", "ones", "ernst", "init_angle", "quadraticPhase30", "RF_shape". With cRandom30, a random pattern of angles, with standard deviation of 30 degrees, is chosen. “ones” and “ernst” initialize to a flat initialization. With the choice “RF_shape”, it is initialized by the function rfFunction (see below) |

rfFunction | function | (undef) | A function with the parameters nTR and nky that outputs a vector of RF pulse values. |

opt_complex | Bool | True | Optimizer results in a complex output. (This option is actually deprecated, since the following one gives much better results; so the current default is strange: has to be changed to false) |

opt_slow_phase | Bool | false | The phase of the pulses becomes an integral of an integral of a pattern that is subject to optimization alongside the amplitude (or, strictly speaking: the real part) of the RF pulses. |

opt_imposed_2nd_derivative_of_phase | Vector{Float} | [] | Allows to overrule aforementioned optimized pattern. If empty, no overruling takes place. So it can be used for a “hand-drawn” 2nd derivative of the phase, provided in rad. The length of the vector must equal nTR. |

opt_focus | String | “mean” | Valid values are: "rho","T1","T2","mean","weighted", "max". If e.g. selecting “T2”, then the sequence will optimize on the noisiness of the T2 map. Despite its default value, currently, the “max” is en vogue. This selects the criterion |

emphasize_low_freq | Bool | false | if switched on, then the optimization will steer more on reducing the noise factor on low spatial frequencies. |

opt_criterion | string | "sar-limited noise" | (R6) |

account_SAR | Bool | false | SAR-limit has to be taken into account |

lambda_B1 | Float | 10.0 | penalty for B1 sensitivity |

lambda_CSF | Float | 0.0 | penalty for CSF sensitivity |

opt_emergeCriterion | Int | 1000 | Every opt_emergeCriterion iterations, the Optimizer will display an intermediate result. |

plottypes | Vector of string | [] | Controls (debugging) plots. Can be a collection of the following values: "first", "trajectories", "weights", "angles", "noisebars", "noiseSpectrum", "infocon" |

optcount | Int | 0 | (For internal use – do not change) |

stage_text | string | “” | (For internal use – do not change) |

plotfuncs | Dict{functions} | (For internal use – do not change) |

Remarks:

(R1)
: maxMeas: Maximum expected number of measurements per phase-encoding value. A ToDo item: this should be automatic, but for the time being, it has to be provided by the user. If there is no variable density, then the default value has to be set to 4 times the number of k-spaces, i.e. 4*nTR/nky in 2D sequences and to 8*nTR/(nky*nkz) if two phase-encoding directions are applied (which points to a flaw in the default, in the case of 3D).

Why 4 or 8 and not 1? One factor of 2 is needed because of assumed positivity of T1 and T2, for which the signal at

(R2)
: invregval: Although the actual reconstruction does not apply explicit regularization, the Analyzer assumes that it does. With the assumed default of

(R3) : The current default of the T1T2set is set to [(0.8183, 0.0509), (0.67, 0.2066), (1.2208, 0.1253), (0.2465, 0.0461), (2.2245, 0.3082), (0.3677, 0.0461), (0.3677, 0.1253)]

(R4)
: The options sizeStep, portion_size, portion_step and opt_iterations_limit are Opimizer options. Following approach is used: in a first stage the Vector sizeStep is relevant; let it be *individual* RF pulses. Since the optimization process cannot reasonably handle a large number (over 50 or so) of parameters simultaneously, it only addresses a portion of size portion_size at a time, taking all other RF pulses as fixed. Once such a portion is considered optimized, the portion advances by portion_step pulses. In some cases, this may result in sequences with pulses around 180 degrees, but not if a sar_limit has been imposed to e.g. 40 degrees RMS.
The total number of internal runs should therefore be

length(sizeStep)+ceil(nTR/portion_step).

This can be reduced by the option parameter opt_iterations_limit. When omitting the second stage (which is advised), the default opt_iterations_limit=length(sizeStep) is appropriate.

(R5)
: Internal Optimizer options. Default is
Optim.Options(time_limit = 3000.0, iterations = 100000, f_tol=1.0e-2, g_tol = 1.0e-3)
The first two options are useful to stop the prevent “the optimizer being stuck, either by non-convergence or due to a parameter choice that would require ages to converge.
On f_tol and g_tol, the latest insights are that these should be set stricter in order to avoid grossly-suboptimal false minima. A value of around

(R6) : Valid values are:

- "noise level" (not accounting for SAR or high flip angles),
- "low-flip corrected noise" (penalty on the maximum flip angle in the sequence, on the assumption that, on a fixed TR, it will thereby reduce the sampling-time and consequently increase the noise by a certain factor),
- "information content", "information content penalized" (immature – do not use),
- (Legacy )"sar-limited noise" (a heavy penalty term is added if the root-mean-square of the RF pulses exceeds a limit). Deprecated. Instead, combine "noise level" with account_SAR=true