Switch to bmat_gradient

src/gradients/changing_gradient_blocks.jl

@@ -3,7 +3,7 @@ import StaticArrays: SVector
 import ...Variables: VariableType, variables
 import ...BuildingBlocks: GradientBlock
 import JuMP: @constraint, @variable, Model, owner_model
-import ...Variables: qvec, bmat, gradient_strength, slew_rate, duration, variables, VariableType
+import ...Variables: qvec, bmat_gradient, gradient_strength, slew_rate, duration, variables, VariableType
 import ...BuildingBlocks: GradientBlock, fixed, RFPulseBlock
 
 struct AbstractChangingGradientBlock{T<:VariableType} <: GradientBlock
@@ -37,7 +37,7 @@ duration(cgb::AbstractChangingGradientBlock) = cgb.duration
 slew_rate(cgb::AbstractChangingGradientBlock) = (cgb.gradient_strength_end .- cgb.gradient_strength_start) ./ duration(cgb)
 qvec(cgb::AbstractChangingGradientBlock) = (cgb.gradient_strength_start .+ cgb.gradient_strength_end) .* (duration(cgb)/2)
 
-function bmat(cgb::AbstractChangingGradientBlock, qstart)
+function bmat_gradient(cgb::AbstractChangingGradientBlock, qstart)
     # grad = (g1 * (duration - t) + g2 * t) / duration
     #      = g1 + (g2 - g1) * t / duration
     # q = qstart + g1 * t + (g2 - g1) * t^2 / (2 * duration)
@@ -53,11 +53,11 @@ function bmat(cgb::AbstractChangingGradientBlock, qstart)
         duration(cgb)^2 .* qstart .* (
             2 .* cgb.gradient_strength_start' .+
             cgb.gradient_strength_end') ./ 3 .+
-        bmat(cgb)
+        bmat_gradient(cgb)
     )
 end
 
-function bmat(cgb::AbstractChangingGradientBlock)
+function bmat_gradient(cgb::AbstractChangingGradientBlock)
     diff = cgb.gradient_strength_end .- cgb.gradient_strength_start
     return (
         cgb.gradient_strength_start .* cgb.gradient_strength_start' ./ 3 .+
@@ -67,6 +67,6 @@ function bmat(cgb::AbstractChangingGradientBlock)
 end
 
 
-variables(::Type{<:AbstractChangingGradientBlock}) = [duration, qvec, bmat, slew_rate]
+variables(::Type{<:AbstractChangingGradientBlock}) = [duration, slew_rate]
 
 end

src/gradients/constant_gradient_blocks.jl

@@ -3,7 +3,7 @@ import StaticArrays: SVector
 import ...Variables: VariableType, variables
 import ...BuildingBlocks: GradientBlock
 import JuMP: @constraint, @variable, Model, owner_model
-import ...Variables: qvec, bmat, gradient_strength, slew_rate, duration, variables, VariableType
+import ...Variables: qvec, bmat_gradient, gradient_strength, slew_rate, duration, variables, VariableType
 import ...BuildingBlocks: GradientBlock, fixed, RFPulseBlock
 
 struct AbstractConstantGradientBlock{T<:VariableType} <: GradientBlock
@@ -37,11 +37,11 @@ gradient_strength(cgb::AbstractConstantGradientBlock) = cgb.gradient_strength
 slew_rate(::AbstractConstantGradientBlock) = zero(SVector{3, Float64})
 qvec(cgb::AbstractConstantGradientBlock) = duration(cgb) .* gradient_strength(cgb)
 
-function bmat(cgb::AbstractConstantGradientBlock)
+function bmat_gradient(cgb::AbstractConstantGradientBlock)
     return (grad .* grad') .* duration(cgb)^3 ./3
 end
 
-function bmat(cgb::AbstractConstantGradientBlock, qstart)
+function bmat_gradient(cgb::AbstractConstantGradientBlock, qstart)
     # \int dt (qstart + t * grad)^2 = 
     #   qstart^2 * duration +
     #   qstart * grad * duration^2 +
@@ -49,10 +49,10 @@ function bmat(cgb::AbstractConstantGradientBlock, qstart)
     return (
         qstart .* qstart' .* duration(cgb) .+
         qstart .* grad' .* duration(cgb)^2 .+
-        bmat(cgb)
+        bmat_gradient(cgb)
     )
 end
 
-variables(::Type{<:AbstractConstantGradientBlock}) = [duration, gradient_strength, qvec, bmat]
+variables(::Type{<:AbstractConstantGradientBlock}) = [duration, gradient_strength]
 
 end

src/gradients/instant_gradients.jl

 module InstantGradients
 import JuMP: @constraint, @variable, Model, owner_model, AbstractJuMPScalar
-import ...Variables: qvec, bmat, duration, variables, get_free_variable, VariableType
+import ...Variables: qvec, bmat_gradient, duration, variables, get_free_variable, VariableType
 import ...BuildingBlocks: GradientBlock, fixed
 import ...BuildSequences: @global_model_constructor
 import ..FixedGradients: FixedInstantGradient
@@ -70,9 +70,9 @@ end
 
 
 qvec(instant::InstantGradientBlock) = instant.qvec
-bmat(instant::InstantGradientBlock, qstart=nothing) = zeros(3, 3)
+bmat_gradient(instant::InstantGradientBlock, qstart=nothing) = zeros(3, 3)
 duration(instant::InstantGradientBlock) = 0.
-variables(::Type{<:InstantGradientBlock}) = [qval]
+variables(::Type{<:InstantGradientBlock}) = [qvec, qval]
 
 function fixed(block::InstantGradientBlock)
     return FixedInstantGradient(block.orientation, value(qval(block)))

src/variables.jl

@@ -19,8 +19,9 @@ 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))"),
-    :area_under_curve => (:gradient, "Area under the curve of the gradient in 1/um (equivalent to [`qval`](@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)."),
     :δ => (: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."),
@@ -53,7 +54,11 @@ variables() = [values(symbol_to_func)...]
 
 
 # Some universal truths
-area_under_curve(bb) = qval(bb)
+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
 
 
 # These functions are more fully defined in building_blocks.jl
@@ -93,6 +98,14 @@ function get_free_variable(model::Model, ::Val{:max})
     return var
 end
 
-function bval end
+"""
+    bmat_gradient(gradient::GradientBlock, qstart=(0, 0, 0))
+
+Computes the diffusion-weighting matrix due to a single gradient block.
+
+This should be defined for every `GradientBlock`, but not be called directly.
+Instead, the `bmat` and `bval` should be constrained for specific `Pathways`
+"""
+function bmat_gradient end
 
 end
\ No newline at end of file