Skip to content
Snippets Groups Projects
Unverified Commit e980ac80 authored by Michiel Cottaar's avatar Michiel Cottaar
Browse files

Set up framework for new variable types

parent 0bc897c6
No related branches found
No related tags found
1 merge request!2Define variables through new @defvar macro
...@@ -10,6 +10,7 @@ Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" ...@@ -10,6 +10,7 @@ Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572" JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84" Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b" MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
......
...@@ -6,15 +6,15 @@ module MRIBuilder ...@@ -6,15 +6,15 @@ module MRIBuilder
include("scanners.jl") include("scanners.jl")
include("build_sequences.jl") include("build_sequences.jl")
include("variables.jl") include("variables.jl")
include("components/components.jl") #include("components/components.jl")
include("containers/containers.jl") #include("containers/containers.jl")
include("pathways.jl") #include("pathways.jl")
include("parts/parts.jl") #include("parts/parts.jl")
include("post_hoc.jl") #include("post_hoc.jl")
include("sequences/sequences.jl") #include("sequences/sequences.jl")
include("printing.jl") #include("printing.jl")
include("sequence_io/sequence_io.jl") #include("sequence_io/sequence_io.jl")
include("plot.jl") #include("plot.jl")
import .BuildSequences: build_sequence, global_model, global_scanner, fixed import .BuildSequences: build_sequence, global_model, global_scanner, fixed
export build_sequence, global_model, global_scanner, fixed export build_sequence, global_model, global_scanner, fixed
...@@ -25,31 +25,31 @@ export Scanner, B0, Siemens_Connectom, Siemens_Prisma, Siemens_Terra, Default_Sc ...@@ -25,31 +25,31 @@ export Scanner, B0, Siemens_Connectom, Siemens_Prisma, Siemens_Terra, Default_Sc
import .Variables: variables, duration, effective_time, flip_angle, amplitude, phase, frequency, bandwidth, N_left, N_right, qval, δ, rise_time, flat_time, slew_rate, gradient_strength, qvec, qval_square, slice_thickness, inverse_slice_thickness, fov, inverse_fov, voxel_size, inverse_voxel_size, resolution, nsamples, oversample, dwell_time, ramp_overlap, spoiler_scale, repetition_time, TR, Δ, get_gradient, get_pulse, get_readout, TE, echo_time, diffusion_time, make_generic, slew_rate3, gradient_strength3, qval3 import .Variables: variables, duration, effective_time, flip_angle, amplitude, phase, frequency, bandwidth, N_left, N_right, qval, δ, rise_time, flat_time, slew_rate, gradient_strength, qvec, qval_square, slice_thickness, inverse_slice_thickness, fov, inverse_fov, voxel_size, inverse_voxel_size, resolution, nsamples, oversample, dwell_time, ramp_overlap, spoiler_scale, repetition_time, TR, Δ, get_gradient, get_pulse, get_readout, TE, echo_time, diffusion_time, make_generic, slew_rate3, gradient_strength3, qval3
export variables, duration, effective_time, flip_angle, amplitude, phase, frequency, bandwidth, N_left, N_right, qval, δ, rise_time, flat_time, slew_rate, gradient_strength, qvec, qval_square, slice_thickness, inversne_slice_thickness, fov, inverse_fov, voxel_size, inverse_voxel_size, resolution, nsamples, oversample, dwell_time, ramp_overlap, spoiler_scale, repetition_time, TR, Δ, get_gradient, get_pulse, get_readout, TE, echo_time, diffusion_time, make_generic, slew_rate3, gradient_strength3, qval3 export variables, duration, effective_time, flip_angle, amplitude, phase, frequency, bandwidth, N_left, N_right, qval, δ, rise_time, flat_time, slew_rate, gradient_strength, qvec, qval_square, slice_thickness, inversne_slice_thickness, fov, inverse_fov, voxel_size, inverse_voxel_size, resolution, nsamples, oversample, dwell_time, ramp_overlap, spoiler_scale, repetition_time, TR, Δ, get_gradient, get_pulse, get_readout, TE, echo_time, diffusion_time, make_generic, slew_rate3, gradient_strength3, qval3
import .Components: InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC, CompositePulse, edge_times #import .Components: InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC, CompositePulse, edge_times
export InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC, CompositePulse, edge_times #export InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC, CompositePulse, edge_times
import .Containers: ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR, readout_times, iter_blocks, iter_instant_gradients, iter_instant_pulses #import .Containers: ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR, readout_times, iter_blocks, iter_instant_gradients, iter_instant_pulses
export ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR, readout_times, iter_blocks, iter_instant_gradients, iter_instant_pulses #export ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR, readout_times, iter_blocks, iter_instant_gradients, iter_instant_pulses
import .Pathways: Pathway, duration_transverse, duration_dephase, bval, bmat, get_pathway #import .Pathways: Pathway, duration_transverse, duration_dephase, bval, bmat, get_pathway
export Pathway, duration_transverse, duration_dephase, bval, bmat, get_pathway #export Pathway, duration_transverse, duration_dephase, bval, bmat, get_pathway
import .Parts: dwi_gradients, readout_event, excitation_pulse, refocus_pulse, Trapezoid, SliceSelect, LineReadout, opposite_kspace_lines, SpoiltSliceSelect, SliceSelectRephase, EPIReadout, interpret_image_size #import .Parts: dwi_gradients, readout_event, excitation_pulse, refocus_pulse, Trapezoid, SliceSelect, LineReadout, opposite_kspace_lines, SpoiltSliceSelect, SliceSelectRephase, EPIReadout, interpret_image_size
export dwi_gradients, readout_event, excitation_pulse, refocus_pulse, Trapezoid, SliceSelect, LineReadout, opposite_kspace_lines, SpoiltSliceSelect, SliceSelectRephase, EPIReadout, interpret_image_size #export dwi_gradients, readout_event, excitation_pulse, refocus_pulse, Trapezoid, SliceSelect, LineReadout, opposite_kspace_lines, SpoiltSliceSelect, SliceSelectRephase, EPIReadout, interpret_image_size
import .PostHoc: adjust, merge_sequences #import .PostHoc: adjust, merge_sequences
export adjust, merge_sequences #export adjust, merge_sequences
import .Sequences: GradientEcho, SpinEcho, DiffusionSpinEcho, DW_SE, DWI #import .Sequences: GradientEcho, SpinEcho, DiffusionSpinEcho, DW_SE, DWI
export GradientEcho, SpinEcho, DiffusionSpinEcho, DW_SE, DWI #export GradientEcho, SpinEcho, DiffusionSpinEcho, DW_SE, DWI
import .SequenceIO: read_sequence, write_sequence #import .SequenceIO: read_sequence, write_sequence
export read_sequence, write_sequence #export read_sequence, write_sequence
import .Plot: plot_sequence #import .Plot: plot_sequence
export plot_sequence #export plot_sequence
import JuMP: @constraint, @objective, objective_function, value, Model #import JuMP: @constraint, @objective, objective_function, value, Model
export @constraint, @objective, objective_function, value, Model #export @constraint, @objective, objective_function, value, Model
end end
...@@ -14,6 +14,7 @@ In addition this defines: ...@@ -14,6 +14,7 @@ In addition this defines:
module Variables module Variables
import JuMP: @constraint, @variable, Model, @objective, objective_function, AbstractJuMPScalar import JuMP: @constraint, @variable, Model, @objective, objective_function, AbstractJuMPScalar
import StaticArrays: SVector import StaticArrays: SVector
import MacroTools
import ..Scanners: gradient_strength, slew_rate, Scanner import ..Scanners: gradient_strength, slew_rate, Scanner
import ..BuildSequences: global_model, global_scanner, fixed import ..BuildSequences: global_model, global_scanner, fixed
...@@ -53,98 +54,107 @@ Can return one of: ...@@ -53,98 +54,107 @@ Can return one of:
adjustable(::AbstractBlock) = :false adjustable(::AbstractBlock) = :false
abstract type AnyVariable end
all_variables_symbols = [ """
:block => [ A sequence property that can be constrained and/or optimised.
:duration => "duration of the building block in ms.",
], It acts as a function, so you can call it on a sequence or building block to get the actual values (e.g., `v(sequence)`).
:sequence => [ It can return one of the following:
:TR => "Time on which an MRI sequence repeats itself in ms. Defaults to the result of [`repetition_time`](@ref).", - a number
:repetition_time => "Time on which an MRI sequence repeats itself in ms.", - a vector of number
:TE => "Echo time of the sequence in ms. Defaults to the result of [`echo_time`](@ref).", - a NamedTuple with the values for individual sequence components
:echo_time => "Echo time of the sequence in ms.", """
:diffusion_time => "Diffusion time in ms (i.e., time between start of the diffusion-weighted gradients).", struct Variable <: AnyVariable
:Δ => "Diffusion time in ms (i.e., time between start of the diffusion-weighted gradients). Defaults to the result of [`diffusion_time`](@ref).", name :: Symbol
:qvec => "Net dephasing due to gradients in rad/um.", f :: Function
:area_under_curve => "Net dephasing due to gradients in rad/um (same as [`qvec`](@ref)).", end
:bmat => "Full 3x3 diffusion-weighting matrix in ms/um^2.",
:bval => "Size of diffusion-weighting in ms/um^2 (trace of [`bmat`](@ref)).",
:duration_dephase => "Net T2' dephasing experienced over a sequence (in ms).",
:duration_transverse => "Net T2 signal loss experienced over a sequence (in ms). This is the total duration spins spent in the transverse plane.",
:delay => "Delay between readout and spin echo in an asymmetric spin echo in ms.",
],
:pulse => [
:flip_angle => "The flip angle of the RF pulse in degrees",
:amplitude => "The maximum amplitude of an RF pulse in kHz",
:phase => "The angle of the phase of an RF pulse in KHz",
:frequency => "The off-resonance frequency of an RF pulse (relative to the Larmor frequency of water) in KHz",
:bandwidth => "Bandwidth of the RF pulse in kHz. To set constraints it is often better to use [`inverse_bandwidth`](@ref).",
:inverse_bandwidth => "Inverse of bandwidth of the RF pulse in 1/kHz. Also see [`bandwidth`](@ref).",
:N_left => "The number of zero crossings of the RF pulse before the main peak",
:N_right => "The number of zero crossings of the RF pulse after the main peak",
:slice_thickness => "Slice thickness of an RF pulse that is active during a gradient in mm. To set constraints it is often better to use [`inverse_slice_thickness`](@ref).",
:inverse_slice_thickness => "Inverse of slice thickness of an RF pulse that is active during a gradient in 1/mm. Also, see [`slice_thickness`](@ref).",
],
:gradient => [
:qval3 => "The spatial range with orientation on which the displacements can be detected due to this gradient in rad/um (also see [`qval`](@ref)).",
:qval => "The spatial range on which the displacements can be detected due to this gradient in rad/um. This will be a scalar if the orientation is fixed and a scalar otherwise. If you always want a vector, use [`qvec`](@ref).",
:qval_square => "Square of [`qval`](@ref) in rad^2/um^2.",
:δ => "Effective duration of a gradient pulse ([`rise_time`](@ref) + [`flat_time`](@ref)) in ms.",
:rise_time => "Time for gradient pulse to reach its maximum value in ms.",
:flat_time => "Time of gradient pulse at maximum value in ms.",
:gradient_strength => "Vector with maximum strength of a gradient along each dimension (kHz/um)",
:slew_rate => "Vector with maximum slew rate of a gradient along each dimension (kHz/um/ms)",
:spoiler_scale => "Length-scale on which spins will be dephased by exactly 2π in mm.",
],
:readout => [
:dwell_time => "Time between two samples in an `ADC` in ms.",
:nsamples => "Number of samples during a readout. During the optimisation this might produce non-integer values.",
:fov => "Size of the field of view in mm. To set constraints it is often better to use [`inverse_fov`](@ref).",
:inverse_fov => "Inverse of size of the field of view in 1/mm. Also see [`fov`](@ref).",
:voxel_size => "Size of each voxel in mm. To set constraints it is often better to use [`inverse_voxel_size`](@ref).",
:inverse_voxel_size => "Inverse of voxel size in 1/mm. Also see [`voxel_size`](@ref).",
:resolution => "Number of voxels in the final readout. During the optimisation this might produce non-integer values, but this will be fixed after optimsation.",
:oversample => "How much to oversample with ([`nsamples`](@ref) / [`resolution`](@ref))",
:ramp_overlap => "Fraction of overlap between ADC event and underlying gradient pulse ramp (between 0 and 1)."
],
]
""" """
Collection of all functions that return variables that can be used to query or constrain their values. variable_defined_for(var, Val(type))
Check whether variable is defined for a specific sub-type.
""" """
variables = Dict{Symbol, Function}() variable_defined_for(var::Variable, ::Val{T}) where {T <: AbstractBlock} = hasmethod(var.f, (T, ))
all_variables = Dict{Symbol, AnyVariable}()
macro defvar(func_def)
func_names = []
function adjust_function(ex)
try
fn_def = MacroTools.splitdef(ex)
push!(func_names, fn_def[:name])
fn_def[:name] = Expr(:., fn_def[:name], QuoteNode(:f))
return MacroTools.combinedef(fn_def)
catch e
if e isa AssertionError
return ex
end
rethrow()
end
end
new_func_def = MacroTools.prewalk(adjust_function, func_def)
function fix_function_name(ex)
if ex in func_names
return esc(ex)
else
return ex
end
end
new_func_def = MacroTools.postwalk(fix_function_name, new_func_def)
for (block_symbol, all_functions) in all_variables_symbols expressions = Expr[]
for (func_symbol, description) in all_functions for func_name in func_names
as_string = " $func_symbol($block_symbol)\n\n$description\n\nThis represents a variable within the sequence. Variables can be set during the construction of a [`AbstractBlock`](@ref) or used to create constraints after the fact." push!(expressions, quote
new_func = @eval begin $(esc(func_name)) = if $(QuoteNode(func_name)) in keys(all_variables)
function $func_symbol end all_variables[$(QuoteNode(func_name))]
@doc $as_string $func_symbol else
$func_symbol function $(func_name) end
all_variables[$(QuoteNode(func_name))] = Variable($(QuoteNode(func_name)), $(func_name))
end
end end
variables[func_symbol] = new_func )
end end
return Expr(
:block,
expressions...,
new_func_def
)
end end
@defvar function duration end
@doc """ duration(block)
TE(ab::AbstractBlock) = echo_time(ab) Duration of the sequence or building block in ms.
TR(ab::AbstractBlock) = repetition_time(ab) """ duration
Δ(ab::AbstractBlock) = diffusion_time(ab)
"""
Dictionary with alternative versions of specific function. struct AlternateVariable <: AnyVariable
name :: Symbol
Setting constraints on these alternative functions can be helpful as it avoids some operations, which the optimiser might struggle with. other_var :: Symbol
""" to_other :: Function
alternative_variables = Dict( from_other :: Function
qval => (qval_square, n->n^2, sqrt, false), inverse :: Bool
slice_thickness => (inverse_slice_thickness, inv, inv, true), end
spoiler_scale => (qval, q->1e-3 * 2π/q, l->1e-3 * 2π/l, true),
bandwidth => (inverse_bandwidth, inv, inv, true), all_variables[:qval] = AlternateVariable(:qval, :qval_square, q->q^2, sqrt, false)
fov => (inverse_fov, inv, inv, true), all_variables[:spoiler_scale] = AlternateVariable(:spoiler_scale, :spoiler_scale, q->1e-3 * 2π/q, l->1e-3 * 2π/l, true)
voxel_size => (inverse_voxel_size, inv, inv, true),
) for name in [:slice_thickness, :bandwidth, :fov, :voxel_size]
inv_name = Symbol("inverse_" * string(name))
all_variables[name] = AlternateVariable(name, inv_name, inv, inv, true)
end
for (name, alt_name) in [
(:TE, :echo_time),
(:TR, :repetition_time),
(:Δ, :diffusion_time),
]
all_variables[name] = AlternateVariable(name, alt_name, identity, identity, false)
end
""" """
...@@ -245,6 +255,8 @@ mutable struct VariableNotAvailable <: Exception ...@@ -245,6 +255,8 @@ mutable struct VariableNotAvailable <: Exception
end end
VariableNotAvailable(bb::Type{<:AbstractBlock}, variable::Function) = VariableNotAvailable(bb, variable, nothing) VariableNotAvailable(bb::Type{<:AbstractBlock}, variable::Function) = VariableNotAvailable(bb, variable, nothing)
function Base.showerror(io::IO, e::VariableNotAvailable) function Base.showerror(io::IO, e::VariableNotAvailable)
if isnothing(e.alt_variable) if isnothing(e.alt_variable)
print(io, e.variable, " is not available for block of type ", e.bb, ".") print(io, e.variable, " is not available for block of type ", e.bb, ".")
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment