Make qvec always a vector and qval dependent on whether orientation is fixed

src/MRIBuilder.jl

@@ -21,8 +21,11 @@ export build_sequence, global_model, global_scanner
 import .Scanners: Scanner, B0, Siemens_Connectom, Siemens_Prisma, Siemens_Terra
 export Scanner, B0, Siemens_Connectom, Siemens_Prisma, Siemens_Terra
 
-import .Variables: variables, duration, start_time, end_time, 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
-export variables, duration, start_time, end_time, 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
+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
+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
+
+import .ContainerBlocks: ContainerBlock, start_time, end_time
+export start_time, end_time
 
 import .Components: InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC
 export InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC

src/all_building_blocks/base_building_blocks.jl

@@ -138,21 +138,21 @@ end
 
 
 """
-    qvec(overlapping[, first_event, last_event])
+    qval(overlapping[, first_event, last_event])
 
 Computes the area under the curve for the gradient waveform in [`BaseBuildingBlock`](@ref).
 
 If `first_event` is set to something else than `nothing`, only the gradient waveform after this RF pulse/Readout will be considered.
 Similarly, if `last_event` is set to something else than `nothing`, only the gradient waveform up to this RF pulse/Readout will be considered.
 """
-function qvec(bb::BaseBuildingBlock, index1::Union{Nothing, Integer}, index2::Union{Nothing, Integer})
+function qval(bb::BaseBuildingBlock, index1::Union{Nothing, Integer}, index2::Union{Nothing, Integer})
     @assert isnothing(index1) || isnothing(index2) || index2 >= index1
     if (index1 isa Number) && (index1 == index2)
         return zeros(3)
     end
-    sum(qvec.(waveform_sequence(bb, index1, index2)); init=zero(SVector{3, Float64}))
+    sum(qval.(waveform_sequence(bb, index1, index2)); init=zero(SVector{3, Float64}))
 end
-qvec(bb::BaseBuildingBlock) = qvec(bb, nothing, nothing)
+qval(bb::BaseBuildingBlock) = qval(bb, nothing, nothing)
 
 function bmat_gradient(bb::BaseBuildingBlock, qstart, index1::Union{Nothing, Integer}, index2::Union{Nothing, Integer})
     @assert isnothing(index1) || isnothing(index2) || index2 >= index1

src/all_building_blocks/spoilt_slice_selects.jl

@@ -4,7 +4,7 @@ import LinearAlgebra: norm
 import StaticArrays: SVector
 import JuMP: @constraint, @objective, objective_function
 import ...BuildSequences: global_model, global_scanner
-import ...Variables: VariableType, duration, rise_time, flat_time, effective_time, qvec, gradient_strength, slew_rate, inverse_slice_thickness, get_free_variable, get_pulse
+import ...Variables: VariableType, duration, rise_time, flat_time, effective_time, qval, gradient_strength, slew_rate, inverse_slice_thickness, get_free_variable, get_pulse
 import ...Components: ChangingGradient, ConstantGradient, RFPulseComponent
 import ..BaseBuildingBlocks: BaseBuildingBlock
 
@@ -54,14 +54,14 @@ function SpoiltSliceSelect(pulse::RFPulseComponent; orientation=[0, 0, 1], group
     @constraint model res.diff_time <= res.rise_time1
     @constraint model res.diff_time <= res.fall_time2
 
-    @constraint model qvec(res, nothing, 1) == qvec(res, 1, nothing)
+    @constraint model qval(res, nothing, 1) == qval(res, 1, nothing)
     if !isnothing(spoiler_scale)
-        if spoiler_scale == :min # note that spoiler_scale is inverse of qvec
-            @objective model Min objective_function(model) - qvec(res, nothing, 1)
+        if spoiler_scale == :min # note that spoiler_scale is inverse of qval
+            @objective model Min objective_function(model) - qval(res, nothing, 1)
         elseif spoiler_scale == :max
-            @objective model Min objective_function(model) + qvec(res, nothing, 1)
+            @objective model Min objective_function(model) + qval(res, nothing, 1)
         else
-            @constraint model qvec(res, nothing, 1) >= 2π * 1e-3 / spoiler_scale
+            @constraint model qval(res, nothing, 1) >= 2π * 1e-3 / spoiler_scale
         end
     end
     set_simple_constraints!(res, kwargs)

src/all_building_blocks/trapezoids.jl

@@ -6,7 +6,7 @@ module Trapezoids
 import JuMP: @constraint
 import StaticArrays: SVector
 import LinearAlgebra: norm
-import ...Variables: qvec, rise_time, flat_time, slew_rate, gradient_strength, variables, duration, δ, get_free_variable, VariableType, inverse_bandwidth, effective_time, qval_square, duration, set_simple_constraints!, scanner_constraints!, inverse_slice_thickness
+import ...Variables: qval, rise_time, flat_time, slew_rate, gradient_strength, variables, duration, δ, get_free_variable, VariableType, inverse_bandwidth, effective_time, qval_square, duration, set_simple_constraints!, scanner_constraints!, inverse_slice_thickness
 import ...Variables: Variables, all_variables_symbols, dwell_time, inverse_fov, inverse_voxel_size, fov, voxel_size, get_gradient, get_pulse, get_readout
 import ...BuildSequences: global_model
 import ...Components: ChangingGradient, ConstantGradient, RFPulseComponent, ADC
@@ -44,7 +44,7 @@ If not set, they will be determined during the sequence optimisation.
 - [`duration`](@ref): total pulse duration (2 * `rise_time` + `flat_time`) in ms.
 ### Gradient variables
 - [`gradient_strength`](@ref): Maximum gradient strength achieved during the pulse in kHz/um
-- [`qval`](@ref)/[`qvec`](@ref): Spatial scale on which spins will be dephased due to this pulsed gradient in rad/um (given by `δ` * `gradient_strength`).
+- [`qval`](@ref): Spatial scale on which spins will be dephased due to this pulsed gradient in rad/um (given by `δ` * `gradient_strength`).
 
 The `bvalue` can be constrained for multiple gradient pulses by creating a `Pathway`.
 """
@@ -108,7 +108,7 @@ slew_rate(g::Trapezoid) = g.slew_rate
 δ(g::Trapezoid) = rise_time(g) + flat_time(g)
 duration(g::Trapezoid) = 2 * rise_time(g) + flat_time(g)
 
-qvec(g::BaseTrapezoid, ::Nothing, ::Nothing) = δ(g) .* gradient_strength(g) .* 2π
+qval(g::BaseTrapezoid, ::Nothing, ::Nothing) = δ(g) .* gradient_strength(g) .* 2π
 
 """
     SliceSelect(pulse; orientation=nothing, group=nothing, variables...)

src/components/gradient_waveforms/changing_gradient_blocks.jl

 module ChangingGradientBlocks
 import StaticArrays: SVector
-import ....Variables: VariableType, duration, qvec, bmat_gradient, gradient_strength, slew_rate, get_free_variable
+import ....Variables: VariableType, duration, qval, bmat_gradient, gradient_strength, slew_rate, get_free_variable
 import ...AbstractTypes: GradientWaveform
 
 
@@ -34,7 +34,7 @@ grad_start(cgb::ChangingGradient) = cgb.gradient_strength_start
 slew_rate(cgb::ChangingGradient) = cgb.slew_rate
 grad_end(cgb::ChangingGradient) = grad_start(cgb) .+ slew_rate(cgb) .* duration(cgb)
 gradient_strength(cgb::ChangingGradient) = max.(grad_start(cgb), grad_end(cgb))
-qvec(cgb::ChangingGradient) = (grad_start(cgb) .+ grad_end(cgb)) .* (duration(cgb) * π)
+qval(cgb::ChangingGradient) = (grad_start(cgb) .+ grad_end(cgb)) .* (duration(cgb) * π)
 
 _mult(g1::VariableType, g2::VariableType) = g1 * g2
 _mult(g1::AbstractVector, g2::AbstractVector) = g1 .* permutedims(g2)

src/components/gradient_waveforms/constant_gradient_blocks.jl

 module ConstantGradientBlocks
 import StaticArrays: SVector
-import ....Variables: VariableType, duration, qvec, bmat_gradient, gradient_strength, slew_rate, get_free_variable
+import ....Variables: VariableType, duration, qval, bmat_gradient, gradient_strength, slew_rate, get_free_variable
 import ...AbstractTypes: GradientWaveform
 
 """
@@ -28,8 +28,8 @@ duration(cgb::ConstantGradient) = cgb.duration
 gradient_strength(cgb::ConstantGradient) = cgb.gradient_strength
 slew_rate(::ConstantGradient1D) = 0.
 slew_rate(::ConstantGradient3D) = zero(SVector{3, Float64})
-qvec(cgb::ConstantGradient1D) = duration(cgb) * gradient_strength(cgb) * 2π
-qvec(cgb::ConstantGradient3D) = @. duration(cgb) * gradient_strength(cgb) * 2π
+qval(cgb::ConstantGradient1D) = duration(cgb) * gradient_strength(cgb) * 2π
+qval(cgb::ConstantGradient3D) = @. duration(cgb) * gradient_strength(cgb) * 2π
 
 _mult(g1::VariableType, g2::VariableType) = g1 * g2
 _mult(g1::AbstractVector, g2::AbstractVector) = g1 .* permutedims(g2)

src/components/instant_gradients.jl

 module InstantGradients
 import StaticArrays: SVector, SMatrix
-import ...Variables: VariableType, duration, qvec, bmat_gradient, get_free_variable, set_simple_constraints!, qval, effective_time, make_generic
+import ...Variables: VariableType, duration, qval, bmat_gradient, get_free_variable, set_simple_constraints!, qval, effective_time, make_generic
 import ..AbstractTypes: EventComponent, GradientWaveform
 
 """
-Parent type for [`InstantGradient1D`](@ref) and [`InstantGradient3D`].
+    InstantGradient1D(; orientation=[1, 0, 0], group=nothing, variables...)
+
+Defines an instantaneous gradient with given fixed `orientation` (default: x-direction).
+
+To have a variable `orientation`, see [`InstantGradient3D`](@ref).
 
-This is an [`EventComponent`](@ref) rather than part of the [`GradientWaveform`](@ref),
-because it is more easily thought of as interrupting a finite gradient waveform rather than being part of it.
+## Parameters
+- `orientation` sets the gradient orienation as a length-3 vector. If not set, the gradient can be in any direction.
+- `group`: name of the group to which this gradient belongs (used for scaling and rotating).
+
+## Variables
+- [`qval`](@ref): Spatial frequency on which spins will be dephased due to this pulsed gradient in rad/um (scalar if `orientation` is set and vector otherwise).
+- [`spoiler_scale`](@ref): Length-scale on which spins will be dephased by exactly 2π in mm.
 """
 abstract type InstantGradient <: EventComponent end
 
@@ -39,7 +48,6 @@ function InstantGradient1D(; orientation=[1, 0, 0], group=nothing, qval=nothing,
 end
 
 qval(ig::InstantGradient1D) = ig.qval
-qvec(ig::InstantGradient1D) = qval(ig) .* ig.orientation
 
 
 """
@@ -53,8 +61,7 @@ To have a fixed `orientation`, see [`InstantGradient1D`](@ref).
 - `group`: name of the group to which this gradient belongs (used for scaling and rotating).
 
 ## Variables
-- [`qvec`](@ref): Vector with spatial frequency on which spins will be dephased due to this pulsed gradient in rad/um.
-- [`qval`](@ref): Norm of spatial frequency on which spins will be dephased due to this pulsed gradient in rad/um.
+- [`qval`](@ref): Vector of spatial frequency on which spins will be dephased due to this pulsed gradient in rad/um.
 - [`spoiler_scale`](@ref): Vector with length-scale on which spins will be dephased by exactly 2π in mm.
 """
 struct InstantGradient3D <: InstantGradient
@@ -62,16 +69,16 @@ struct InstantGradient3D <: InstantGradient
     group :: Union{Nothing, Symbol}
 end
 
-function InstantGradient3D(; qvec=[nothing, nothing, nothing], group=nothing, variables...)
-    if isnothing(qvec)
-        qvec = [nothing, nothing, nothing]
+function InstantGradient3D(; qval=[nothing, nothing, nothing], group=nothing, variables...)
+    if isnothing(qval)
+        qval = [nothing, nothing, nothing]
     end
-    res = InstantGradient3D(get_free_variable.(qvec), group)
+    res = InstantGradient3D(get_free_variable.(qval), group)
     set_simple_constraints!(res, variables)
     return res
 end
 
-qvec(ig::InstantGradient3D) = ig.qvec
+qval(ig::InstantGradient3D) = ig.qvec
 
 
 duration(::InstantGradient) = 0.

src/pathways.jl

@@ -4,7 +4,7 @@ import StaticArrays: SVector, SMatrix
 import ..Components: NoGradient, RFPulseComponent, ReadoutComponent, InstantGradient, GradientWaveform
 import ..AllSequences: BaseSequence, Sequence
 import ..AllBuildingBlocks: BaseBuildingBlock, waveform, events, waveform_sequence
-import ..Variables: qvec, qval, bmat_gradient, VariableType, effective_time, duration, qval_square, TR
+import ..Variables: qvec, qval, bmat_gradient, VariableType, effective_time, duration, TR
 import ..ContainerBlocks: start_time
 import ..Alternatives: AlternativeBlocks
 
@@ -127,12 +127,6 @@ You can set `scale` and/or `rotate` to specific symbols to only consider gradien
 """
 qvec(pathway::Pathway; scale=nothing, rotate=nothing) = get(pathway.qvec, (scale, rotate), zero(SVector{3, Float64}))
 
-function qval_square(pathway::Pathway; kwargs...)
-    vec = qvec(pathway; kwargs...)
-    return vec[1]^2 + vec[2]^2 + vec[3]^2
-end
-qval(pathway::Pathway; kwargs...) = sqrt(qval_square(pathway; kwargs...))
-
 
 """
     area_under_curve(pathway::Pathway; scale=nothing, rotate=nothing)

src/variables.jl

@@ -9,6 +9,7 @@ In addition this defines:
 - [`set_simple_constraints`](@ref): call [`apply_simple_constraint`](@ref) for each keyword argument.
 - [`apply_simple_constraint`](@ref): set a simple equality constraint.
 - [`get_pulse`](@ref)/[`get_gradient`](@ref)/[`get_readout`](@ref): Used to get the pulse/gradient/readout part of a building block
+- [`gradient_orientation`](@ref): returns the gradient orientation of a waveform if fixed.
 """
 module Variables
 import JuMP: @constraint, @variable, Model, @objective, objective_function, value, AbstractJuMPScalar
@@ -40,8 +41,8 @@ all_variables_symbols = [
         :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 => [
-        :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)). To set constraints it is often better to use [`qvec`](@ref) or [`qval_square`](@ref).",
+        :qvec => "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.",
@@ -172,6 +173,20 @@ Instead, the `bmat` and `bval` should be constrained for specific `Pathways`
 function bmat_gradient end
 
 
+"""
+    gradient_orientation(building_block)
+
+Returns the gradient orientation.
+"""
+function gradient_orientation(bb::AbstractBlock)
+    if hasproperty(bb, :orientation)
+        return bb.orientation
+    else
+        return gradient_orientation(get_gradient(bb))
+    end
+end
+
+
 function effective_time end
 
 
@@ -231,12 +246,15 @@ for (target_name, all_vars) in all_variables_symbols
     end
 end
 
-function Variables.qval_square(bb::AbstractBlock, args...; kwargs...)
-    vec = Variables.qvec(bb, args...; kwargs...)
-    return vec[1]^2 + vec[2]^2 + vec[3]^2
-end
 
-Variables.qval(bb::AbstractBlock, args...; kwargs...) = sqrt(Variables.qval_square(bb, args...; kwargs...))
+function qvec(bb::AbstractBlock, args...; kwargs...)
+    value = qval(bb, args...; kwargs...)
+    if value isa AbstractVector
+        return value
+    else
+        return gradient_orientation(bb) .* value
+    end
+end
 
 
 """