Do not compute scalar version of gradient_strength and slew_rate

src/building_blocks.jl

@@ -122,11 +122,32 @@ end
 
 function scanner_constraints!(model::Model, building_block::BuildingBlock, scanner::Scanner)
     for func in [gradient_strength, slew_rate]
-        if (func in variables(building_block)) && isfinite(func(scanner))
-            @constraint model func(building_block) <= func(scanner)
+        if isfinite(func(scanner))
+            scanner_constraints!(model, building_block, scanner, func)
         end
     end
-    if building_block isa ContainerBlock
+end
+
+function scanner_constraints!(model::Model, building_block::BuildingBlock, scanner::Scanner, func::Function)
+    if func in variables(building_block)
+        # apply constraint at this level
+        res_bb = func(building_block)
+        if res_bb isa AbstractVector
+            if isnothing(getproperty(building_block, :rotate, true))
+                # no rotation; apply constraint to each dimension independently
+                for expr in res_bb
+                    @constraint model expr <= func(scanner)
+                end
+            else
+                # with rotation: apply constraint to total squared
+                total_squared = sum(map(n->n^2, res_bb))
+                @constraint model total_squared <= func(scanner)^2
+            end
+        else
+            @constraint model res_bb <= func(scanner)
+        end
+    elseif building_block isa ContainerBlock
+        # apply constraints at lower level
         for (_, child_block) in get_children_blocks(building_block)
             scanner_constraints!(model, child_block, scanner)
         end

src/gradients/changing_gradient_blocks.jl

@@ -35,6 +35,7 @@ const FixedChangingGradientBlock = AbstractChangingGradientBlock{Float64}
 duration(cgb::AbstractChangingGradientBlock) = cgb.duration
 
 slew_rate(cgb::AbstractChangingGradientBlock) = (cgb.gradient_strength_end .- cgb.gradient_strength_start) ./ duration(cgb)
+gradient_strength(cgb::AbstractChangingGradientBlock) = max.(cgb.gradient_strength_start, cgb.gradient_strength_start)
 qvec(cgb::AbstractChangingGradientBlock) = (cgb.gradient_strength_start .+ cgb.gradient_strength_end) .* (duration(cgb)/2)
 
 function bmat_gradient(cgb::AbstractChangingGradientBlock, qstart)

src/gradients/pulsed_gradients.jl

@@ -37,13 +37,15 @@ If not set, they will be determined during the sequence optimisation.
 
 The [`bvalue`](@ref) can be constrained for multiple gradient pulses.
 """
-mutable struct PulsedGradient <: ContainerBlock
+struct PulsedGradient <: ContainerBlock
     model :: Model
     rise :: ChangingGradientBlock
     flat :: ConstantGradientBlock
     fall :: ChangingGradientBlock
-    slew_rate_vec :: SVector{3, VariableType}
+    slew_rate :: SVector{3, VariableType}
     _scaling :: Union{Nothing, VariableType}
+    rotate :: Union{Nothing, Symbol}
+    scale :: Union{Nothing, Symbol}
 end
 
 @global_model_constructor PulsedGradient
@@ -79,15 +81,17 @@ function PulsedGradient(model::Model; orientation=nothing, rise_time=nothing, fl
     end
     rise_time = get_free_variable(model, rise_time)
     flat_time = get_free_variable(model, flat_time)
-    grad_vec = slew_rate .* rise_time
+    grad = slew_rate .* rise_time
 
     res = PulsedGradient(
         model,
-        ChangingGradientBlock(zeros(3), grad_vec, rise_time, rotate, scale),
-        ConstantGradientBlock(grad_vec, flat_time, rotate, scale),
-        ChangingGradientBlock(grad_vec, zeros(3), rise_time, rotate, scale),
+        ChangingGradientBlock(zeros(3), grad, rise_time, rotate, scale),
+        ConstantGradientBlock(grad, flat_time, rotate, scale),
+        ChangingGradientBlock(grad, zeros(3), rise_time, rotate, scale),
         slew_rate,
         rate_1d,
+        rotate,
+        scale
     )
 
     set_simple_constraints!(model, res, kwargs)
@@ -98,14 +102,11 @@ end
 
 rise_time(pg::PulsedGradient) = duration(pg.rise)
 flat_time(pg::PulsedGradient) = duration(pg.flat)
-gradient_strength_vec(g::PulsedGradient) = rise_time(g) * slew_rate_vec(g)
-gradient_strength(g::PulsedGradient) = isnothing(g._scaling) ? maximum(gradient_strength_vec(g)) : (res._scaling * rise_time(g))
-slew_rate_vec(g::PulsedGradient) = g.slew_rate_vec
-slew_rate(g::PulsedGradient) = isnothing(g._scaling) ? maximum(abs.(slew_rate_vec(g))) : res._scaling
+gradient_strength(g::PulsedGradient) = gradient_strength(g.flat)
+slew_rate(g::PulsedGradient) = g.slew_rate
 δ(g::PulsedGradient) = rise_time(g) + flat_time(g)
 duration(g::PulsedGradient) = 2 * rise_time(g) + flat_time(g)
-qvec(g::PulsedGradient) = δ(g) .* gradient_strength_vec(g)
-qval(g::PulsedGradient) = δ(g) * gradient_strength(g)
+qvec(g::PulsedGradient) = δ(g) .* gradient_strength(g)
 
 get_children_indices(::PulsedGradient) = (:rise, :flat, :fall)
 Base.getindex(pg::PulsedGradient, symbol::Symbol) = pg[Val(symbol)]
@@ -119,16 +120,16 @@ start_time(pg::PulsedGradient, ::Val{:flat}) = rise_time(pg)
 start_time(pg::PulsedGradient, ::Val{:fall}) = δ(pg)
 
 
-variables(::Type{<:PulsedGradient}) = [qval, δ, gradient_strength_vec, duration, rise_time, flat_time]
+variables(::Type{<:PulsedGradient}) = [qval, δ, gradient_strength, duration, rise_time, flat_time]
 
 
 function fixed(block::PulsedGradient)
-    grad_vec = value.(gradient_strength_vec(block))
+    grad = value.(gradient_strength(block))
     t_rise = value(rise_time(block))
     t_d = value(δ(block))
     return FixedGradient(
         [0., t_rise, t_d, td + t_rise],
-        [zeros(3), grad_vec, grad_vec, zeros(3)]; 
+        [zeros(3), grad, grad, zeros(3)]; 
         rotate=rotate
     )
 end

src/variables.jl

@@ -19,14 +19,13 @@ all_variables_symbols = [
     :slice_thickness => (:pulse, "Slice thickness of an RF pulse that is active during a gradient."),
 
     # gradients
-    :qvec => (:gradient, "The spatial range and orientation on which the displacements can be detected due to this gradient in 1/um (equivalent to [`area_under_curve_vec`](@ref))"),
-    :qval => (:gradient, "The spatial range on which the displacements can be detected due to this gradient in 1/um (equivalent to [`area_under_curve`](@ref))"),
-    :qval_sqr => (:gradient, "Square of qval in 1/um^2. This can be more efficient to calculate and hence to contrain than `qval` (as it avoids computing a square root)."),
+    :qvec => (:gradient, "The spatial range and orientation on which the displacements can be detected due to this gradient in 1/um."),
+    :qval => (:gradient, "The spatial range on which the displacements can be detected due to this gradient in 1/um."),
     :δ => (:gradient, "Effective duration of a gradient pulse ([`rise_time`](@ref) + [`flat_time`](@ref)) in ms."),
     :rise_time => (:gradient, "Time for gradient pulse to reach its maximum value in ms."),
     :flat_time => (:gradient, "Time of gradient pulse at maximum value in ms."),
-    :slew_rate => (:gradient, "maximum increase of a gradient (kHz/um/ms)"),
-    :gradient_strength => (:gradient, "maximum strength of a gradient (kHz/um)"),
+    :gradient_strength => (:gradient, "vector with maximum strength of a gradient along each dimension (kHz/um)"),
+    :slew_rate => (:gradient, "vector with maximum slew rate of a gradient along each dimension (kHz/um)"),
 ]
 
 symbol_to_func = Dict{Symbol, Function}()
@@ -54,11 +53,8 @@ variables() = [values(symbol_to_func)...]
 
 
 # Some universal truths
-qval(bb) = sqrt(qval_sqr(bb))
-function qval_sqr(bb)
-    vec = qvec(bb)
-    return vec[1] * vec[1] + vec[2] * vec[2] + vec[3] * vec[3]
-end
+qval(bb) = norm(qvec(bb))
+
 
 
 # These functions are more fully defined in building_blocks.jl