Tutorial
This tutorial demonstrates how to implement the GrapeMR.jl package by scripting a setup and directly calling its functions.
- Define your physical system: Spins, relaxation values, inhomogeneities, etc.
- Optimization parameters: Scheduler parameters and maximum iterations.
- Grape Parameters: Time steps, cost function and which fields to optimize.
- Generate initial control field: Cubic spline is used by default.
using GrapeMRSaturation contrast without inhomogeneities
Physical system
Define the initial magnetization state, relaxation times (in seconds), and spin labels. Specify target according to the desired outcome (cost function-dependent) and offset for $B_0$ inhomogeneities in Hertz. Finally, create a Spin object with all spins. Each unique inhomogeneity combination is treated as a separate spin.
M0 = [0.0, 0.0, 1.0]
T1 = [0.6, 0.1]
T2 = [0.3, 0.05]
ΔB1 = [1.0]
B0 = 0.0
offset = collect(-B0:1:B0)
label = ["s1", "s2"]
target = ["min", "max"]
spins = Spin(M0, T1, T2, offset, ΔB1, target, label)2-element Vector{Spin}:
Spin([0.0, 0.0, 1.0], 0.6, 0.3, 0.0, 1.0, "min", "s1", 2.0)
Spin([0.0, 0.0, 1.0], 0.1, 0.05, 0.0, 1.0, "max", "s2", 2.0)Parameters
Optimization Parameters
max_iter = get(ENV, "DEV", "false") == "true" ? 1 : 2000 # we set max_iter to 1 if we're in development mode to build the docs faster
poly_start = 0.5
poly_degree = 2
opt_params = OptimizationParams(poly_start, poly_degree, max_iter)OptimizationParams(0.5, 2, 2000)Grape Parameters
N = 2000
cost = GrapeMR.saturation_contrast
fields2opt = Dict("B1x" => true, "B1y" => true, "Bz" => false)
grape_params = GrapeParams(N, cost, fields2opt)GrapeParams{typeof(GrapeMR.saturation_contrast)}(2000, GrapeMR.saturation_contrast, Dict{String, Bool}("B1y" => 1, "B1x" => 1, "Bz" => 0))Parameter struct
params = Parameters(grape_params, opt_params)Parameters{typeof(GrapeMR.saturation_contrast)}(GrapeParams{typeof(GrapeMR.saturation_contrast)}(2000, GrapeMR.saturation_contrast, Dict{String, Bool}("B1y" => 1, "B1x" => 1, "Bz" => 0)), OptimizationParams(0.5, 2, 2000))Initial control field
Generate an initial control field using the spline function. Set the control time control_time and reference amplitude B1ref.
B1ref = 1.0
control_time = 0.5
control_field = spline_RF(N, control_time, B1ref)ControlField{Float64, Matrix{Float64}, Matrix{Float64}}([0.3841906028583226 0.3848705247668154 … 0.9802906446344607 0.9834930691088447], [0.3841906028583226 0.3848705247668154 … 0.9802906446344607 0.9834930691088447], 1.0, [0.0 0.0 … 0.0 0.0], 0.5)Run Optimization
Run the optimization with the contructed Spins, configured Parameters and ControlField.
grape_output = grape(params, control_field, spins);
Final Cost Function Value = 0.018
---------- RF Analysis ----------
Pulse Peak Amplitude = 3.616 [Hz]
Pulse Peak Amplitude = 0.0849 [μT]
Attenuation corresponding to maximum amplitude: 49.3256 [dB]
Maximum power = 1.7088 [W]
Average power = 0.3647 [W]
Pulse energy = 0.1824 [J]Plot
julia> using GrapeMR, Plots; unicodeplots(); # change the backend so that plots go to stdout and can be rendered in CI/headless mode.julia> control_fields = plot_control_fields(grape_output.control_field);julia> display(control_fields)