From f93bebeb7056963af114ae8e0dc54b2014064d5e Mon Sep 17 00:00:00 2001
From: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
Date: Sun, 28 Jan 2024 21:17:41 +0000
Subject: [PATCH] Switch to bmat_gradient
---
src/gradients/changing_gradient_blocks.jl | 10 +++++-----
src/gradients/constant_gradient_blocks.jl | 10 +++++-----
src/gradients/instant_gradients.jl | 6 +++---
src/variables.jl | 19 ++++++++++++++++---
4 files changed, 29 insertions(+), 16 deletions(-)
diff --git a/src/gradients/changing_gradient_blocks.jl b/src/gradients/changing_gradient_blocks.jl
index b424871..a759eac 100644
--- a/src/gradients/changing_gradient_blocks.jl
+++ b/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
diff --git a/src/gradients/constant_gradient_blocks.jl b/src/gradients/constant_gradient_blocks.jl
index cbcae10..dce6fb3 100644
--- a/src/gradients/constant_gradient_blocks.jl
+++ b/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
diff --git a/src/gradients/instant_gradients.jl b/src/gradients/instant_gradients.jl
index 58c96bb..bb23766 100644
--- a/src/gradients/instant_gradients.jl
+++ b/src/gradients/instant_gradients.jl
@@ -1,6 +1,6 @@
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)))
diff --git a/src/variables.jl b/src/variables.jl
index 8630ade..e8d0180 100644
--- a/src/variables.jl
+++ b/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
--
GitLab