diff --git a/src/MRIBuilder.jl b/src/MRIBuilder.jl
index 06704351c5f189b286671165030e437b1bca1ee2..3da8a36cad50ffaa5a137603a672e0c196ffb160 100644
--- a/src/MRIBuilder.jl
+++ b/src/MRIBuilder.jl
@@ -12,6 +12,7 @@ include("gradients/gradients.jl")
include("pulses/pulses.jl")
include("containers/containers.jl")
include("readouts/readouts.jl")
+include("pathways.jl")
import .BuildSequences: build_sequence
export build_sequence
@@ -40,6 +41,9 @@ export InstantReadout
import .Containers: Sequence, FixedBlock, FixedSequence
export Sequence, FixedBlock, FixedSequence
+import .Pathways: Pathway, duration_transverse, duration_dephase, bval, bmat
+export duration_transverse, duration_dephase, bval, bmat
+
import JuMP: @constraint, @objective, objective_function, optimize!, has_values, value, owner_model, Model
export @constraint, @objective, objective_function, optimize!, has_values, value, owner_model, Model
diff --git a/src/pathways.jl b/src/pathways.jl
index 0255ebcad363c333fd77ca1a7a4d490ab514a274..c9dad816b3ca8224cc521bd931c781722f5647b1 100644
--- a/src/pathways.jl
+++ b/src/pathways.jl
@@ -1,9 +1,10 @@
module Pathways
import LinearAlgebra: norm
-import StaticArrays: SVector, SMatrix
+import StaticArrays: SVector, SMatrix, MVector, MMatrix
+import ..BuildingBlocks: BuildingBlock, GradientBlock, RFPulseBlock
import ..Containers: Sequence
-import ..Variables: qval, bval, area_under_curve
-import ..PathwayWalkers: PathwayWalker
+import ..Variables: qvec, qval, bmat_gradient, VariableType
+
"""
Pathway(sequence::Sequence, pulse_effects::Vector{:Symbol/Number}, readout_index=1)
@@ -48,7 +49,7 @@ You can select which gradients to consider when accessing these values.
struct Pathway
# user provided
sequence :: Sequence
- pulse_effects :: Vector{Union{:Symbol, :Number}}
+ pulse_effects :: Vector{Union{Symbol, Number}}
readout_index :: Integer
# computed
@@ -59,12 +60,11 @@ end
function Pathway(sequence::Sequence, pulse_effects::AbstractVector, readout_index::Integer)
walker = PathwayWalker()
- walk_pathway!(sequence, copy(pulse_effects), walker, Ref(readout_index))
+ walk_pathway!(sequence, interpret_pulse_effects.(pulse_effects), walker, Ref(readout_index))
return Pathway(
sequence,
pulse_effects,
readout_index,
-
Dict(k => SVector{4, VariableType}(v) for (k, v) in pairs(walker.duration_states)),
Dict(k => SVector{3, VariableType}(v) for (k, v) in pairs(walker.qvec)),
Dict(k => SMatrix{3, 3, VariableType, 9}(v) for (k, v) in pairs(walker.bmat)),
@@ -84,13 +84,7 @@ The requested state can be set using `transverse` and `positive` as follows:
- `transverse=true`, `positive=false`: -transverse
"""
function duration_state(pathway::Pathway, transverse, positive)
- index = Dict([
- (false, true) => 1,
- (true, true) => 2,
- (false, false) => 3,
- (true, false) => 4,
- ])[(Bool(transverse), Bool(positive))]
- return pathway.duration_states[index]
+ return pathway.duration_states[duration_state_index(transverse, positive)]
end
"""
@@ -179,4 +173,233 @@ You can set `scale` and/or `rotate` to specific symbols to only consider gradien
"""
bval(pathway::Pathway; scale=nothing, rotate=nothing) = tr(bmat(pathway; scale, rotate))
+
+"""
+ interpret_pulse_effects(number_or_symbol)
+
+Interpret the various numbers and symbols that can be passed on to a `Pathway`.
+
+The result will be one of:
+- :ignore (if input is 0, :ignore, or :skip).
+- :excite (if input is 90 or :excite).
+- :refocus (if input is 180, :refocus, or :excite).
+- :neg_excite (if input is -90, 270, or :negexcite).
+"""
+function interpret_pulse_effects(number::Number)
+ normed = Int(number % 360)
+ if normed == 0
+ return :ignore
+ elseif normed == 90
+ return :excite
+ elseif normed == 180
+ return :refocus
+ elseif normed == 270
+ return :neg_excite
+ else
+ error("The pulse effect along a pathway should be divisible by 90, not $number.")
+ end
+end
+
+function interpret_pulse_effects(sym::Symbol)
+ mapping = Dict{Symbol, Symbol}(
+ :ignore => :ignore,
+ :skip => :ignore,
+ :excite => :excite,
+ :net_excite => :net_excite,
+ :refocus => :refocus,
+ :invert => :refocus,
+ )
+ if sym in keys(mapping)
+ return mapping[sym]
+ else
+ all_symbols = join(keys(mapping), ", ")
+ error("The pulse effect along a pathway should be one of ($all_symbols), not $sym.")
+ end
+end
+
+
+"""
+Helper structure for [`PathwayWalker`](@ref), which is itself a helper for `Pathway`. You are deep down the rabit hole now...
+
+For documentation, see that structure and [`walk_pathway`](@ref) and [`update_walker_gradient!`](@ref).
+"""
+mutable struct GradientTracker
+ last_gradient_time :: VariableType
+ qvec :: MVector{3, VariableType}
+ bmat :: MMatrix{3, 3, VariableType, 9}
+end
+
+GradientTracker() = GradientTracker(0., zeros(3), zeros(3, 3))
+
+
+"""
+Helper structure for `Pathway`.
+
+For documentation, see that structure and [`walk_pathway`](@ref).
+"""
+mutable struct PathwayWalker
+ last_pulse_time :: VariableType
+ is_transverse :: Bool
+ is_positive :: Bool
+ duration_states :: MVector{4, VariableType}
+ gradient_trackers :: Dict{Any, GradientTracker}
+end
+
+PathwayWalker() = PathwayWalker(
+ 0., false, true,
+ zeros(4),
+ Dict{Any, GradientTracker}()
+)
+
+"""
+ walk_pathway!(bb::BuildingBlock, walker::PathwayWalker, pulse_effects::Vector, nreadout::Ref{Int}, start_time)
+
+Computes the effect of a specific [`BuildingBlock`](@ref) (starting at `start_time`) on the [`PathwayWalker`](@ref).
+
+This needs to be defined for any [`ContainerBlock`](@ref) that the walker may encounter.
+
+This might be as simpla as calling [`walk_pathway!`] on all the children (i.e., for a `Sequence`).
+
+For individual pulses and gradients, the following behaviour is implemented:
+- If a pulse is encountered, call [`update_walker_pulse!`](@ref)`(walker, pulse_effects, pulse_effective_time)`
+- If a gradient is encountered, call [`update_walker_gradient!`](@ref)(gradient, walker, gradient_start_time)
+
+For overlapping gradients/pulses, one should first encounter the part of the gradient before the [`effective_time`](@ref) of the pulse,
+then apply the pulse, and then the rest of the gradient.
+This effective time can be passed on to [`update_walker_gradient!`](@ref) to allow part of the gradient waveform to be applied.
+
+The function should return `true` if the `Pathway` has reached its end (i.e., the final readout) and `false` otherwise.
+"""
+function walk_pathway!(grad::GradientBlock, walker::PathwayWalker, pulse_effects::Vector{Symbol}, nreadout::Ref{Int}, start_time=0.::VariableType)
+ update_walker_gradient!(grad, walker, start_time)
+ return false
+end
+
+function walk_pathway!(pulse::RFPulseBlock, walker::PathwayWalker, pulse_effects::Vector{Symbol}, nreadout::Ref{Int}, start_time=0.::VariableType)
+ update_walker_pulse(walker, pulse_effects, nreadout, start_time + effective_time(pulse))
+ return iszero(length(pulse_effects)) && iszero(nreadout[])
+end
+
+
+"""
+ update_walker_pulse!(walker::PathwayWalker, pulse_effects::Vector, pulse_time)
+
+Apply the first element of `pulse_effects` to the `walker` at the given `pulse_time`.
+
+The following steps will be taken if the first `pulse_effect` is not `:ignore`
+- if `walker.transverse` is true before the pulse, increase the `walker.bmat` by the outer product of `walker.qvec` with itself multiplied by the time since the last gradient
+- update `walker.duration_states` with time since last pulse.
+- update `walker.last_pulse_time`
+- update `walker.is_transverse`, and `walker.is_positive` based on the first value in `pulse_effects`.
+- if `walker.is_positive` changed in the previous step than the `walker.qvec` needs to be flipped.
+- remove the first element from `pulse_effects`.
+
+If the first element is `:ignore` the only effect is that the first element is removed from `pulse_effects`.
+"""
+function update_walker_pulse!(walker::PathwayWalker, pulse_effects::AbstractVector{Symbol}, pulse_time::VariableType)
+ if length(pulse_effects) == 0
+ error("Pathway definition is invalid! Another RF pulse was encountered before the number of readouts expected from `nreadout` where detected.")
+ end
+ instruction = popfirst!(pulse_effects)
+ if instruction == :ignore
+ return
+ end
+
+ # update qvec/bmat
+ if walker.is_transverse
+ for gradient_tracker in values(walker.gradient_trackers)
+ gradient_tracker.bmat += (
+ (gradient_tracker.qvec .* gradient_tracker.qvec') .*
+ (pulse_time - gradient_tracker.last_gradient_time)
+ )
+ gradient_tracker.last_gradient_time = pulse_time
+ end
+ end
+ prev_sign = walker.is_positive
+
+ # update durations
+ index = duration_state_index(walker.is_transverse, walker.is_positive)
+ walker.duration_states[index] = walk.duration_state[index] + (pulse_time - walker.last_pulse_time)
+
+ walker.last_pulse_time = pulse_time
+
+ # -transverse, +longitudinal, +transverse, -longitudinal, -transverse, +longitudinal
+ ordering = [(true, false), (false, true), (true, true), (false, false), (true, false)]
+
+ if instruction == :refocus
+ walker.is_positive = !walker.is_positive
+ elseif instruction == :excite
+ index = findfirst(isequal(walker.is_transverse, walker.is_positive), ordering)
+ (walker.is_transverse, walker.is_positive) = ordering[index + 1]
+ elseif instruction == :neg_excite
+ index = findlast(isequal(walker.is_transverse, walker.is_positive), ordering)
+ (walker.is_transverse, walker.is_positive) = ordering[index - 1]
+ else
+ error("Invalid pulse instruction ($instruction); This error should have been caught earlier.")
+ end
+
+ # flip qvec if needed
+ if prev_sign != walker.is_positive
+ for gradient_tracker in values(walker.gradient_trackers)
+ gradient_tracker.qvec = -gradient_tracker.qvec
+ end
+ end
+end
+
+"""
+ update_walker_gradient!(gradient_block::GradientBlock, walker::PathwayWalker, gradient_start_time::VariableType; overlapping_pulses=[], overlapping_readouts=[])
+
+Update the walker's `qvec` and `bmat` based on the given `gradient_block`.
+
+The following steps will be taken:
+- Do nothing if `walker.transverse` is false
+- increase the appropriate `walker.bmat` by the outer product of `walker.qvec` with itself multiplied by the time since the last gradient
+- update the appropriate `walker.qvec` and `walker.bmat` based on the gradient waveform. This will require appropriate `qvec`/`bmat` functions to be defined for the gradient building block.
+- update `walker.last_gradient_time` to the time at the end of the gradient.
+
+This requires [`bmat`](@ref) and [`qvec`](@ref) to be implemented for the [`GradientBlock`](@ref).
+"""
+function update_walker_gradient!(gradient::GradientBlock, walker::PathwayWalker, gradient_start_time::VariableType, internal_start_time, internal_end_time)
+ if walker.transverse
+ return
+ end
+
+ if iszero(internal_start_time) || isnothing(internal_start_time)
+ # only worry about this for the first call
+
+ # make sure the appropriate gradient tracker exists
+ key = (gradient.scale, gradient.rotate)
+ if !(key in keys(walker.gradient_trackers))
+ walker.gradient_trackers[key] = GradientTracker()
+ end
+ tracker = walker.gradient_trackers[key]
+
+ # update bmat till start of gradient
+ tracker.bmat = tracker.bmat .+ (
+ (tracker.qvec .* tracker.qvec') .*
+ (gradient_start_time - tracker.last_gradient_time)
+ )
+ tracker.last_gradient_time = gradient_start_time
+ end
+
+ tracker.bmat = tracker.bmat .+ bmat(gradient, tracker.qvec, internal_start_time, internal_end_time)
+ tracker.qvec = tracker.qvec .+ qvec(gradient, internal_start_time, internal_end_time)
+end
+
+"""
+ duration_state_index(transverse, positive)
+
+Returns the index of a specific state in the `Pathway.duration_state` vector.
+
+This function is used internally to access that vector.
+"""
+function duration_state_index(transverse, positive)
+ return Dict([
+ (false, true) => 1,
+ (true, true) => 2,
+ (false, false) => 3,
+ (true, false) => 4,
+ ])[(Bool(transverse), Bool(positive))]
+end
+
end
\ No newline at end of file