From 1b66a519d88f5afd2c2fcb614c2cb1a6b441d49c Mon Sep 17 00:00:00 2001 From: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk> Date: Mon, 5 Feb 2024 17:00:04 +0000 Subject: [PATCH] Support new variables interface --- src/building_blocks.jl | 28 ++++----- src/overlapping/abstract.jl | 6 +- .../gradient_readouts/single_lines.jl | 6 +- src/readouts/ADCs.jl | 9 +-- src/variables.jl | 57 +++---------------- 5 files changed, 34 insertions(+), 72 deletions(-) diff --git a/src/building_blocks.jl b/src/building_blocks.jl index 252b4b7..653502c 100644 --- a/src/building_blocks.jl +++ b/src/building_blocks.jl @@ -129,7 +129,7 @@ for variable_func in keys(variables) value = alt_var(bb) catch e if e isa VariableNotAvailable - throw(VariableNotAvailable(typeof(bb), variable_func)) + throw(VariableNotAvailable(typeof(bb), Variables.$variable_func)) end rethrow() end @@ -138,9 +138,9 @@ for variable_func in keys(variables) elseif value isa AbstractArray{<:Number} return backward.(value) end - throw(VariableNotAvailable(typeof(bb), variable_func, alt_var)) + throw(VariableNotAvailable(typeof(bb), Variables.$variable_func, alt_var)) end - throw(VariableNotAvailable(typeof(bb), variable_func)) + throw(VariableNotAvailable(typeof(bb), Variables.$variable_func)) end end @@ -208,21 +208,21 @@ function Base.show(io::IO, printer::BuildingBlockPrinter) end try numeric_value = _robust_value(fn(block)) + if isnothing(numeric_value) + continue + end + if numeric_value isa AbstractVector + printed_value = "[" * join(map(v -> @sprintf("%.3g", v), numeric_value), ", ") * "]" + else + printed_value = @sprintf("%.3g", numeric_value) + end + print(io, "$(nameof(fn))=$(printed_value), ") catch e if e isa VariableNotAvailable continue end rethrow() end - if isnothing(numeric_value) - continue - end - if numeric_value isa AbstractVector - printed_value = "[" * join(map(v -> @sprintf("%.3g", v), numeric_value), ", ") * "]" - else - printed_value = @sprintf("%.3g", numeric_value) - end - print(io, "$(nameof(fn))=$(printed_value), ") end print(io, ")") end @@ -239,8 +239,8 @@ If set to `:min` or `:max`, minimising or maximising this function will be added """ function set_simple_constraints!(block::BuildingBlock, kwargs) for (key, value) in kwargs - if key in keys(alternative_variables) - alt_var, forward, backward, to_invert = alternative_variables[Variables.$variable_func] + if variables[key] in keys(alternative_variables) + alt_var, forward, backward, to_invert = alternative_variables[variables[key]] invert_value(value::VariableType) = forward(value) invert_value(value::Symbol) = invert_value(Val(value)) invert_value(::Val{:min}) = to_invert ? Val(:max) : Val(:min) diff --git a/src/overlapping/abstract.jl b/src/overlapping/abstract.jl index f2c9d76..4936fe1 100644 --- a/src/overlapping/abstract.jl +++ b/src/overlapping/abstract.jl @@ -1,6 +1,6 @@ module Abstract -import ...BuildingBlocks: BuildingBlock, ContainerBlock, RFPulseBlock, get_children_indices +import ...BuildingBlocks: BuildingBlock, ContainerBlock, RFPulseBlock, get_children_indices, VariableNotAvailable import ...Wait: WaitBlock import ...Readouts: InstantReadout import ...Gradients: GradientBlock, split_gradient @@ -59,9 +59,9 @@ for (func_symbol, _) in [Dict(all_variables_symbols)[:pulse]..., :effective_time @eval function $func_symbol(ao::AbstractOverlapping, args...; kwargs...) single_pulse = pulses(ao) if iszero(length(single_pulse)) - error("Building block does not contain any RF pulse, so cannot compute $func_symbol.") + throw(VariableNotAvailable(typeof(ao), $func_symbol)) elseif length(single_pulse) > 1 - error("Building block has multipl pulses, so cannot compute $func_symbol.") + error("Building block has multipl pulses, so cannot compute ", $func_symbol, ".") end $func_symbol(single_pulse[1], args...; kwargs...) end diff --git a/src/overlapping/gradient_readouts/single_lines.jl b/src/overlapping/gradient_readouts/single_lines.jl index 4a174d1..6228e2a 100644 --- a/src/overlapping/gradient_readouts/single_lines.jl +++ b/src/overlapping/gradient_readouts/single_lines.jl @@ -1,7 +1,7 @@ module SingleLines import ....BuildingBlocks: set_simple_constraints! import ....Readouts: ADC -import ....Variables: dwell_time, fov, resolution, nsamples, effective_time, gradient_strength, slew_rate, voxel_size, rise_time, duration, variables +import ....Variables: dwell_time, inverse_fov, resolution, nsamples, effective_time, gradient_strength, slew_rate, inverse_voxel_size, rise_time, duration, variables import ...Abstract: AbstractOverlapping, waveform, interruptions import ...GradientPulses: TrapezoidGradient @@ -27,8 +27,8 @@ effective_time(sl::SingleLine) = rise_time(sl.grad) + effective_time(sl.adc) dwell_time(sl::SingleLine) = dwell_time(sl.adc) slew_rate(sl::SingleLine) = slew_rate(sl.grad)[1] gradient_strenth(sl::SingleLine) = slew_rate(sl) * rise_time(sl.grad) -fov_inverse(sl::SingleLine) = 1e3 * dwell_time(sl) * gradient_strenth(sl) -voxel_size_inverse(sl::SingleLine) = 1e3 * duration(sl.adc) * gradient_strenth(sl) +inverse_fov(sl::SingleLine) = 1e3 * dwell_time(sl) * gradient_strenth(sl) +inverse_voxel_size(sl::SingleLine) = 1e3 * duration(sl.adc) * gradient_strenth(sl) nsamples(sl::SingleLine) = nsamples(sl.adc) resolution(sl::SingleLine) = resolution(sl.adc) duration(sl::SingleLine) = duration(sl.grad) diff --git a/src/readouts/ADCs.jl b/src/readouts/ADCs.jl index d40e38b..1f3b30c 100644 --- a/src/readouts/ADCs.jl +++ b/src/readouts/ADCs.jl @@ -1,5 +1,5 @@ module ADCs -import JuMP: @constraint +import JuMP: @constraint, value import ...Variables: variables, dwell_time, duration, effective_time, get_free_variable, VariableType, nsamples, resolution, oversample import ...BuildingBlocks: BuildingBlock, apply_simple_constraint!, set_simple_constraints!, fixed import ...BuildSequences: global_model @@ -53,9 +53,10 @@ resolution(adc::ADC) = adc.resolution function fixed(adc::ADC) # round nsamples during fixing - nsamples = Int(round(adc.nsamples)) - dwell_time = duration(adc) / nsamples - return ADC(nsamples, dwell_time, time_to_center(adc)) + r = Int(round(value(resolution(adc)))) + n = Int(round(value(nsamples(adc)))) + oversample = n // r + return ADC(r, value(adc.dwell_time), value(adc.time_to_center), oversample) end end \ No newline at end of file diff --git a/src/variables.jl b/src/variables.jl index 3e43d15..9d65d52 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -16,16 +16,19 @@ all_variables_symbols = [ :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)", + :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).", ], # gradients :gradient => [ :qvec => "The spatial range and orientation on which the displacements can be detected due to this gradient in rad/um.", - :qval => "The spatial range on which the displacements can be detected due to this gradient in rad/um (i.e., norm of [`qvec`](@ref)).", + :qval => "The spatial range on which the displacements can be detected due to this gradient in rad/um (i.e., norm of [`qvec`](@ref)). To set constraints it is often better to use [`qvec`](@ref) or [`qval_square`](@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.", @@ -35,8 +38,10 @@ all_variables_symbols = [ :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.", - :voxel_size => "Size of each voxel in mm.", + :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))", ] @@ -60,54 +65,10 @@ end # helper functions -""" - inverse_slice_thickness(pulse) - -Inverse of [`slice_thickness`](@ref) in ms. - -It is defined separately from `slice_thickness` to avoid unnecessary divisions in any constraint. -""" -function inverse_slice_thickness end - -""" - inverse_bandwidth(pulse) - -Inverse of [`bandwidth`](@ref) in 1/mm. - -It is defined separately from `bandwidth` to avoid unnecessary divisions in any constraint. -""" -function inverse_bandwidth end - -""" - inverse_fov(readout) - -Inverse of [`fov`](@ref) in 1/mm. - -It is defined separately from `fov` to avoid unnecessary divisions in any constraint. -""" -function inverse_fov end - -""" - inverse_voxel_size(readout) - -Inverse of [`voxel_size`](@ref) in 1/mm. - -It is defined separately from `voxel_size` to avoid unnecessary divisions in any constraint. -""" -function inverse_voxel_size end - -""" - qval_square(gradient) - -Square of [`qval`](@ref) in rad/um. - -It is defined separately from `qval` to avoid unnecessary square root in any constraint. -""" function qval_square(bb, args...; kwargs...) vec = qvec(bb, args...; kwargs...) return vec[1]^2 + vec[2]^2 + vec[3]^2 end -qval(bb; kwargs...) = sqrt(qval_square(bb)) alternative_variables = Dict( qval => (qval_square, n->n^2, sqrt, false), -- GitLab