Split `split_times` from `linearise`

src/MRIBuilder.jl

@@ -26,8 +26,8 @@ export variables, duration, effective_time, flip_angle, amplitude, phase, freque
 import .Components: InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC
 export InstantPulse, ConstantPulse, SincPulse, GenericPulse, InstantGradient, SingleReadout, ADC
 
-import .Containers: ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR, readout_times, edge_times, linearise
-export ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR, readout_times, edge_times, linearise
+import .Containers: ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR, readout_times, edge_times, SequencePart, split_times, split_timestep, linearise
+export ContainerBlock, start_time, end_time, waveform, waveform_sequence, events, BaseBuildingBlock, BuildingBlock, Wait, BaseSequence, nrepeat, Sequence, AlternativeBlocks, match_blocks!, get_index_single_TR, readout_times, edge_times, SequencePart, split_times, split_timestep, linearise
 
 import .Pathways: Pathway, duration_transverse, duration_dephase, bval, bmat, get_pathway
 export Pathway, duration_transverse, duration_dephase, bval, bmat, get_pathway

src/containers/containers.jl

@@ -9,6 +9,6 @@ import .Abstract: ContainerBlock, start_time, end_time, effective_time, readout_
 import .BuildingBlocks: BaseBuildingBlock, BuildingBlock, Wait, waveform, waveform_sequence, events, ndim_grad
 import .BaseSequences: BaseSequence, Sequence, nrepeat, get_index_single_TR
 import .Alternatives: AlternativeBlocks, match_blocks!
-import .Linearise: linearise
+import .Linearise: SequencePart, split_times, split_timestep, linearise
 
 end
\ No newline at end of file

src/containers/linearise.jl

@@ -2,7 +2,7 @@ module Linearise
 import StaticArrays: SVector
 import ...Components: GradientWaveform, split_timestep
 import ...Variables: amplitude, phase, gradient_strength3, duration
-import ..Abstract: edge_times, start_time, end_time
+import ..Abstract: edge_times, start_time, end_time, ContainerBlock
 import ..BaseSequences: BaseSequence, Sequence
 import ..BuildingBlocks: BaseBuildingBlock
 
@@ -14,7 +14,13 @@ end
 
 
 """
-Represents a small part of a larger [`Sequence`](@ref), where the gradient, RF amplitude, and RF phase can all be modeled as changing linearly.
+    SequencePart(sequence, time1, time2)
+
+Represents the time between `time1` and `time2` of a larger [`Sequence`](@ref)
+
+The gradient, RF amplitude, and RF phase are all be modeled as changing linearly during this time.
+
+See [`linearise`](@ref) to split a sequence into such linear parts.
 """
 struct SequencePart
     gradient :: LinearPart{SVector{3, Float64}}
@@ -24,22 +30,67 @@ struct SequencePart
     duration :: Float64
 end
 
-duration(sp::SequencePart) = sp.duration
 
+function SequencePart(sequence::BaseSequence{N}, time1::Number, time2::Number) where {N}
+    if !(0 <= time1 <= time2 <= duration(sequence))
+        error("Sequence timings are out of bound")
+    end
+    for key in 1:N
+        if (end_time(sequence, key) > time1)
+            return SequencePart(sequence[key], time1 - start_time(sequence, key), time2 - start_time(sequence, key))
+        end
+    end
+end
+
+function SequencePart(bb::BaseBuildingBlock, time1::Number, time2::Number)
+    if !(0 <= time1 <= time2 <= duration(bb))
+        error("Sequence timings are out of bound")
+    end
+
+    tmean = (time1 + time2) / 2
 
-function linearise(seq::BaseSequence; kwargs...)
-    parts = SequencePart[]
-    for block in seq
-        append!(parts, linearise(block; kwargs...))
+    function fit_linear_part(fn)
+        start = fn(bb, time1)
+        final = fn(bb, time2)
+        middle = fn(bb, tmean)
+        mid_deviation = @. middle - (start + final) / 2
+        return LinearPart(start .+ mid_deviation ./ 3, final .+ mid_deviation ./ 3)
     end
-    return parts
+
+    return SequencePart(
+        fit_linear_part(gradient_strength3),
+        zero(SVector{3, Float64}),
+        fit_linear_part(amplitude),
+        fit_linear_part(phase),
+        time2 - time1
+    )
 end
 
-function linearise(bb::BaseBuildingBlock; precision=0.01)
-    parts = SequencePart[]
-    if duration(bb) <= 0.
-        return parts
+duration(sp::SequencePart) = sp.duration
+
+"""
+    split_times(sequence; precision=0.01)
+
+Suggests at what times to split a sequence into linear parts (see [`linearise`](@ref)).
+
+The split times will include any time when:
+- the starting/end points of building blocks, gradients, RF pulses, or ADC readouts.
+- a gradient or RF pulse is discontinuous in the first derivative
+- the time of any instantaneous gradients, RF pulses, or readouts.
+
+Continuous gradient waveforms or RF pulses might be split up further to ensure the linear approximations meet the required `precision` (see [`split_timestep`](@ref)).
+"""
+function split_times(sequence::BaseSequence; precision=0.01)
+    splits = Float64[]
+    for block in sequence
+        append!(splits, split_times(block; precision=precision))
     end
+    return sort(unique(splits))
+end
+
+
+function split_times(bb::BaseBuildingBlock; precision=0.01)
+    splits = Float64[]
     edges = edge_times(bb)
     for (t1, t2) in zip(edges[1:end-1], edges[2:end])
         tmean = (t1 + t2) / 2
@@ -50,36 +101,28 @@ function linearise(bb::BaseBuildingBlock; precision=0.01)
             if !(start_time(bb, key) < tmean < end_time(bb, key))
                 continue
             end
-            new_timestep = split_timestep(block, precision)
+            new_timestep = split_timestep(bb[key], precision)
             if new_timestep < timestep
                 timestep = new_timestep
             end
         end
-        nsteps = div(t2 - t1, timestep, RoundUp)
-        tsplits = range(t1, t2, length=nsteps)
-
-        for (t1b, t2b) in zip(tsplits[1:end-1], tsplits[2:end])
-            tmeanb = (t1b + t2b) / 2
-
-            function fit_linear_part(fn)
-                start = fn(bb, t1b)
-                final = fn(bb, t2b)
-                middle = fn(bb, tmeanb)
-                mid_deviation = @. middle - (start + final) / 2
-                return LinearPart(start .+ mid_deviation ./ 3, final .+ mid_deviation ./ 3)
-            end
-
-            push!(parts, SequencePart(
-                fit_linear_part(gradient_strength3),
-                zero(SVector{3, Float64}),
-                fit_linear_part(amplitude),
-                fit_linear_part(phase),
-                t2b - t1b
-            ))
-        end
+        nsteps = (isinf(timestep) ? 1 : Int(div(t2 - t1, timestep, RoundUp))) + 1
+        append!(splits, range(t1, t2, length=nsteps))
     end
-    return parts
+    return sort(unique(splits))
 end
 
 
+
+"""
+    linearise(sequence[, times])
+
+Splits any [`BaseSequence`](@ref) or [`BaseBuildingBlock`](@ref) into a series of [`SequencePart`](@ref) objects where the gradients/pulses are approximated to be linear.
+
+If the `times` are not explicitly set they will be obtained from [`split_times`](@ref).
+"""
+linearise(container::ContainerBlock) = linearise(container, split_times(container))
+linearise(container::ContainerBlock, times::AbstractVector{<:Number}) = [SequencePart(container, t1, t2) for (t1, t2) in zip(times[1:end-1], times[2:end])]
+
+
 end
\ No newline at end of file

src/pathways.jl

@@ -311,10 +311,6 @@ function walk_pathway!(block::BaseBuildingBlock, walker::PathwayWalker, pulse_ef
     current_time = block_start_time
     for (index_inter, interruption) in events(block)
 
-        if interruption isa Tuple
-            interruption = interruption[2]
-        end
-
         # determine if action should be taken
         if interruption isa RFPulseComponent
             if iszero(length(pulse_effects))

src/variables.jl

@@ -13,6 +13,7 @@ In addition this defines:
 """
 module Variables
 import JuMP: @constraint, @variable, Model, @objective, objective_function, AbstractJuMPScalar
+import StaticArrays: SVector
 import ..Scanners: gradient_strength, slew_rate, Scanner
 import ..BuildSequences: global_model, global_scanner, fixed
 
@@ -290,7 +291,7 @@ for base_fn in [:qval, :gradient_strength, :slew_rate]
     @eval function $fn3(bb::AbstractBlock, args...; kwargs...)
         value = $base_fn(bb, args...; kwargs...)
         if value isa Number && iszero(value)
-            return value
+            return zero(SVector{3, Float64})
         elseif value isa AbstractVector
             return value
         else