Build and optimise MRI sequences

This is a copy from the "jump-sequences" branch in `MCMRSimulator`.

src/MRIBuilder.jl

+"""
+Builds and optimises NMR/MRI sequences.
+"""
 module MRIBuilder
 
-# Write your package code here.
+include("building_blocks.jl")
+include("sequence_builders.jl")
+include("wait.jl")
+include("gradients/gradients.jl")
+include("pulses/pulses.jl")
+include("readouts/readouts.jl")
+
+import .BuildingBlocks: BuildingBlock, scanner_constraints!
+export BuildingBlock, scanner_constraints!
+
+import .Wait: WaitBlock
+export WaitBlock
+
+import .SequenceBuilders: SequenceBuilder, start_time, end_time, duration, TR
+export SequenceBuilder, start_time, end_time, duration, TR
+
+import .Gradients: PulsedGradient, InstantGradientBlock, qval, rise_time, flat_time, slew_rate, gradient_strength, bval
+export PulsedGradient, InstantGradientBlock, qval, rise_time, flat_time, slew_rate, gradient_strength, bval
+
+import .Pulses: InstantRFPulseBlock, ConstantPulse, SincPulse, flip_angle, phase, frequency, bandwidth, N_left, N_right
+export InstantRFPulseBlock, ConstantPulse, SincPulse, flip_angle, phase, frequency, bandwidth, N_left, N_right
+
+import .Readouts: InstantReadout
+export InstantReadout
 
 end

src/building_blocks.jl

+module BuildingBlocks
+import JuMP: has_values, GenericVariableRef, value, Model, @constraint, @objective, owner_model, objective_function
+import Printf: @sprintf
+import ...Sequences: RFPulse, InstantRFPulse, MRGradients, InstantGradient
+import ...Scanners: Scanner
+
+"""
+Parent type for all individual components out of which a sequence can be built.
+
+Required methods:
+- [`duration`](@ref)(block, parameters): returns block duration in ms
+- [`to_mcmr_components`](@ref)(block, parameters): converts the block into components recognised by the MCMR simulator
+- [`properties`](@ref): A list of all functions that are used to compute properties of the building block. Any of these can be used in constraints or objective functions.
+"""
+abstract type BuildingBlock end
+
+
+"""
+    duration(building_block)
+
+The duration of the building block in ms.
+"""
+function duration end
+
+"""
+    to_mcmr_components(building_block, parameters)
+
+Converts the building block into components recognised by the MCMR simulator. These components are:
+- [`RFPulse`](@ref)
+- [`InstantRFPulse`](@ref)
+- [`MRGradients`](@ref)
+- [`InstantGradient`](@ref)
+"""
+function to_mcmr_components end
+
+
+"""
+    gradient_strength(gradient)
+
+Maximum gradient strength in kHz/um.
+
+If a [`Scanner`](@ref) is provided, this will be constrained to be lower than the maximum scanner gradient strength.
+"""
+function gradient_strength end
+
+"""
+    slew_rate(gradient)
+
+Maximum rate of increase (and decrease) of the gradient strength in kHz/um/ms.
+
+If a [`Scanner`](@ref) is provided, this will be constrained to be lower than the maximum scanner slew rate.
+"""
+function slew_rate end
+
+"""
+    scanner_constraints!([model, ]building_block, scanner)
+
+Adds any constraints from a specific scanner to a [`BuildingBlock`]{@ref}.
+If a [`Scanner`](@ref) is provided to the `SequenceBuilder`, this will be called under the hood on every [`BuildingBlock`]@(ref).
+"""
+function scanner_constraints!(building_block::BuildingBlock, scanner::Scanner)
+    scanner_constraints!(owner_model(building_block), building_block, scanner)
+end
+
+function scanner_constraints!(model::Model, building_block::BuildingBlock, scanner::Scanner)
+    for (func, property_name) in [
+        (gradient_strength, :gradient)
+        (slew_rate, :slew_rate)
+    ]
+        if (func in properties(building_block)) && isfinite(getproperty(scanner, property_name))
+            @constraint model func(building_block) <= getproperty(scanner, property_name)
+        end
+    end
+end
+
+"""
+    properties(building_block)
+
+Returns a list of function that can be called to constrain the `building_block`.
+"""
+properties(bb::BuildingBlock) = properties(typeof(bb))
+
+struct _BuildingBlockPrinter
+    bb :: BuildingBlock
+    number :: Integer
+end
+
+function Base.show(io::IO, block::BuildingBlock)
+    print(io, string(typeof(block)), "(")
+    if has_values(block)
+        if iszero(value(duration(block)))
+            print(io, "time=", @sprintf("%.3f", value(start_time(block))), ", ")
+        else
+            print(
+                io, "time=", 
+                @sprintf("%.3f", value(start_time(block))), " - ",
+                @sprintf("%.3f", value(end_time(block))), ", "
+            )
+        end
+    end
+    for name in propertynames(block)
+        value = getproperty(block, name)
+        if value isa GenericVariableRef || name == :builder || string(name)[1] == '_'
+            continue
+        end
+        print(io, name, "=", repr(value), ", ")
+    end
+
+    if has_values(block)
+        for fn in properties(block)
+            if fn == duration
+                continue
+            end
+            numeric_value = value(fn(block))
+            printed_value = @sprintf("%.3g", numeric_value)
+            print(io, "$(nameof(fn))=$(printed_value), ")
+        end
+    end
+    print(io, ")")
+end
+
+# The `start_time` and `end_time` functions are properly defined in "sequence_builders.jl"
+function start_time end
+function end_time end
+
+
+"""
+    set_simple_constraints!(model, block, kwargs)
+
+Add any constraints or objective functions to the properties of a [`BuildingBlock`](@ref) in the JuMP `model`.
+
+Each keyword argument has to match one of the functions in [`properties`](@ref)(block).
+If set to a numeric value, a constraint will be added to fix the function value to that numeric value.
+If set to `:min` or `:max`, minimising or maximising this function will be added to the cost function.
+"""
+function set_simple_constraints!(model::Model, block::BuildingBlock, kwargs)
+    to_funcs = Dict(nameof(fn) => fn for fn in properties(block))
+    for (key, value) in kwargs
+        apply_simple_constraint!(model, to_funcs[key](block), value)
+    end
+    nothing
+end
+
+"""
+    apply_simple_constraint!(model, variable, value)
+
+Add a single constraint or objective to the JuMP `model`.
+This is an internal function used by [`set_simple_constraints`](@ref).
+"""
+apply_simple_constraint!(model::Model, variable, ::Nothing) = nothing
+apply_simple_constraint!(model::Model, variable, value::Symbol) = apply_simple_constraint!(model, variable, Val(value))
+apply_simple_constraint!(model::Model, variable, ::Val{:min}) = @objective model Min objective_function(model) + variable
+apply_simple_constraint!(model::Model, variable, ::Val{:max}) = @objective model Min objective_function(model) - variable
+apply_simple_constraint!(model::Model, variable, value::Number) = @constraint model variable == value
+
+
+"""
+    match_blocks!(block1, block2[, property_list])
+
+Matches the listed properties between two [`BuildingBlock`](@ref) objects.
+By default all shared properties (i.e., those with the same name) are matched.
+"""
+function match_blocks!(block1::BuildingBlock, block2::BuildingBlock, property_list)
+    model = owner_model(block1)
+    @assert model == owner_model(block2)
+    for fn in property_list
+        @constraint model fn(block1) == fn(block2)
+    end
+end
+
+function match_blocks!(block1::BuildingBlock, block2::BuildingBlock)
+    property_list = intersect(properties(block1), properties(block2))
+    match_blocks!(block1, block2, property_list)
+end
+
+
+"""
+Stores the parameters passed on a [`BuildingBlock`](@ref) constructor (of type `T`).
+
+The parameters are temporarily stored in this format, until they can be added to a `SequenceBuilder`.
+
+For example, the following
+```julia
+pg = PulsedGradient(qval=2)
+```
+will return a `BuildingBlockPlaceholder` object rather than a `PulsedGradient` object.
+Only when this object is added to a `SequenceBuilder`, is the `PulsedGradient` actually initialised using the JuMP model of that `SequenceBuilder`:
+```julia
+sb = SequenceBuilder([pg])
+```
+You can access the initialised `PulsedGradient` through the `BuildingBlockPlaceholder` (`pg.concrete[]`) or directly through the `SequenceBuilder` (`sb[1]`)
+
+Each Placeholder can be added to only a single `SequenceBuilder`, but it can be added multiple times to the same `SequenceBuilder`.
+If added multiple times to the same `SequenceBuilder`, all variables will be matched between them.
+"""
+struct BuildingBlockPlaceholder{T<:BuildingBlock}
+    args
+    kwargs
+    concrete :: Ref{T}
+    BuildingBlockPlaceholder{T}(args...; kwargs...) where {T<:BuildingBlock} = new{T}(args, kwargs, Ref{T}())
+end
+
+function Base.show(io::IO, placeholder::BuildingBlockPlaceholder{T}) where {T}
+    if isassigned(placeholder.concrete)
+        print(io, "Assigned BuildingBlockPlaceholder for $(placeholder.concrete[])")
+    else
+        args = join(placeholder.args, ", ")
+        kwargs = join(["$key=$value" for (key, value) in pairs(placeholder.kwargs)], ", ")
+        print(io, "Unassigned BuildingBlockPlaceholder{$T}($args; $kwargs)")
+    end
+end
+
+
+end
\ No newline at end of file

src/gradients/gradients.jl

+module Gradients
+include("integrate_gradients.jl")
+include("pulsed_gradients.jl")
+include("instant_gradients.jl")
+
+import .IntegrateGradients: qval, bval
+import .PulsedGradients: PulsedGradient, rise_time, flat_time, slew_rate, gradient_strength
+import .InstantGradients: InstantGradientBlock
+end
\ No newline at end of file

src/gradients/instant_gradients.jl

+module InstantGradients
+import JuMP: @constraint, @variable, VariableRef
+import ....Sequences: InstantGradient
+import ...BuildingBlocks: BuildingBlock, properties, BuildingBlockPlaceholder, set_simple_constraints!, duration, to_mcmr_components
+import ...SequenceBuilders: SequenceBuilder, owner_model, start_time
+import ..IntegrateGradients: qval, bval
+
+"""
+    InstantGradientBlock(; orientation=:bvec, qval=nothing)
+
+Defines an instantaneous gradient.
+
+This is a [`BuildingBlock`](@ref) for the [`SequenceBuilder`](@ref).
+
+## Parameters
+- `orientation` sets the gradient orienation. Can be set to a vector for a fixed orientation. Alternatively, can be set to :bvec (default) to rotate with the user-provided `bvecs` or to :neg_bvec to always be the reverse of the `bvecs`.
+
+## Variables
+- [`qval`](@ref): Spatial scale on which spins will be dephased due to this pulsed gradient in rad/um.
+"""
+struct InstantGradientBlock <: BuildingBlock
+    builder::SequenceBuilder
+    orientation :: Any
+    qval :: VariableRef
+end
+
+InstantGradientBlock(; kwargs...) = BuildingBlockPlaceholder{InstantGradientBlock}(; kwargs...)
+
+function InstantGradientBlock(builder::SequenceBuilder; orientation=:bvec, kwargs...)
+    model = owner_model(builder)
+    res = InstantGradientBlock(
+        builder,
+        orientation,
+        @variable(model)
+    )
+    set_simple_constraints!(model, res, kwargs)
+    return res
+end
+
+
+qval(instant::InstantGradientBlock) = instant.qval
+bval(instant::InstantGradientBlock) = 0.
+duration(instant::InstantGradientBlock) = 0.
+properties(::Type{<:InstantGradientBlock}) = [qval]
+
+
+function to_mcmr_components(block::InstantGradientBlock)
+    if block.orientation == :bvec
+        rotate = true
+        qvec = [value(qval(block)), 0., 0.]
+    elseif block.orientation == :neg_bvec
+        rotate = true
+        qvec = [-value(qval(block)), 0., 0.]
+    elseif block.orientation isa AbstractVector && size(block.orientation) == (3, )
+        rotate = false
+        qvec = block.orientation .* (value(qval(block)) / norm(block.orientation))
+    else
+        error("Gradient orientation should be :bvec, :neg_bvec or a length-3 vector, not $(block.orienation)")
+    end
+    return InstantGradient(qvec, zeros(3), value(start_time(block)), rotate)
+end
+
+
+end
\ No newline at end of file

src/gradients/integrate_gradients.jl

+module IntegrateGradients
+import ...BuildingBlocks: BuildingBlock
+import ...SequenceBuilders: SequenceBuilder, TR, duration, builder, index
+
+
+"""
+    qval(blocks)
+    qval(builder::SequenceBuilder, indices)
+
+Computes the total q-value summed over multiple gradient [`BuildingBlock`](@ref) objects.
+
+The `blocks` should contain the actual [`BuildingBlock`](@ref) objects, possibly interspersed with one of:
+- `:TR` wait for one TR (has no effect, but included for consistency with [`bval`](@ref)).
+- `:flip` flips the current `qval` value (e.g., if a refocus pulse has happened).
+
+When a `builder` is provided, the `indices` are expected to be the integer indices of the individual blocks rather than the actual [`BuildingBlock`](@ref) objects.
+`:TR` and `:flip` can also still be provided.
+In addition, in this interface one can provide negative indices to indicate that a specific gradient should have a negative contribution to the q-value (this is an alternative for using :flip).
+
+The integral can occur over multiple repetition times by including [`BuildingBlock`](@ref) objects out of order (or by using :TR).
+"""
+qval(builder::SequenceBuilder, indices::AbstractVector) = full_integral(builder, indices)[1]
+qval(indices::AbstractVector) = full_integral(indices)[1]
+
+"""
+    bval(blocks)
+    bval(builder::SequenceBuilder, indices)
+
+Computes the total b-value combined over multiple gradient [`BuildingBlock`](@ref) objects.
+
+The `blocks` should contain the actual [`BuildingBlock`](@ref) objects, possibly interspersed with one of:
+- `:TR` wait for one TR (has no effect, but included for consistency with [`bval`](@ref)).
+- `:flip` flips the current `qval` value (e.g., if a refocus pulse has happened).
+
+When a `builder` is provided, the `indices` are expected to be the integer indices of the individual blocks rather than the actual [`BuildingBlock`](@ref) objects.
+`:TR` and `:flip` can also still be provided.
+In addition, in this interface one can provide negative indices to indicate that a specific gradient should have a negative contribution to the q-value (this is an alternative for using :flip).
+
+The integral can occur over multiple repetition times by including [`BuildingBlock`](@ref) objects out of order (or by using :TR).
+"""
+bval(builder::SequenceBuilder, indices::AbstractVector) = full_integral(builder, indices)[2]
+bval(indices::AbstractVector) = full_integral(indices)[2]
+
+function full_integral(blocks::AbstractVector)
+    actual_blocks = filter(b -> b isa BuildingBlock)
+    if length(actual_blocks) == 0
+        return (0., 0.)
+    end
+    sb = builder(blocks[1])
+
+    return full_integral(sb, map(b -> b isa BuildingBlock ? index(sb, b) : b))
+end
+
+function full_integral(builder::SequenceBuilder, indices::AbstractVector)
+    qval_current = 0.
+    current_index = 0
+    bval_current = 0.
+    for index in indices
+        if index == :flip
+            qval_current = -qval_current
+        elseif index == :TR
+            bval_current = bval_current + qval_current^2 * TR(builder)
+        elseif index isa Integer
+            next_gradient = index < 0 ? -index : index
+            wait_time = duration(builder, current_index + 1, next_gradient - 1)
+
+            bval_current = (
+                bval_current + 
+                qval_current^2 * wait_time +
+                bval(builder[next_gradient], index < 0 ? -qval_current : qval_current)
+            )
+            qval_current = (
+                qval_current +
+                sign(index) * qval(builder[next_gradient])
+            )
+        else
+            error("qval/bval indices should refer to gradient blocks or :flip/:TR, not $index")
+        end
+    end
+    return (qval_current, bval_current)
+end
+
+end
\ No newline at end of file

src/gradients/pulsed_gradients.jl

+"""
+Defines a set of different options for MRI gradients.
+"""
+module PulsedGradients
+
+import JuMP: @constraint, @variable, Model, VariableRef, owner_model, value
+import StaticArrays: SVector
+import ....Sequences: MRGradients
+import ...BuildingBlocks: BuildingBlock, duration, properties, set_simple_constraints!, BuildingBlockPlaceholder, gradient_strength, slew_rate, to_mcmr_components
+import ...SequenceBuilders: SequenceBuilder, start_time
+import ..IntegrateGradients: qval, bval
+
+
+"""
+    PulsedGradient(; orientation=:bvec, variables...)
+
+Defines a trapezoidal pulsed gradient
+
+This is a [`BuildingBlock`](@ref) for a [`SequenceBuilder`](@ref).
+
+## Parameters
+- `orientation` sets the gradient orienation. Can be set to a vector for a fixed orientation. Alternatively, can be set to :bvec (default) to rotate with the user-provided `bvecs` or to :neg_bvec to always be the reverse of the `bvecs`.
+
+## Variables
+Variables can be set during construction or afterwards as an attribute.
+If not set, they will be determined during the sequence optimisation.
+### Timing variables
+- [`rise_time`](@ref): Time of the gradient to reach from 0 to maximum in ms. If explicitly set to 0, the scanner slew rate will be ignored.
+- [`flat_time`](@ref): Time that the gradient stays at maximum strength in ms.
+- [`δ`](@ref): effective pulse duration (`rise_time` + `flat_time`) in ms.
+- [`duration`](@ref): total pulse duration (2 * `rise_time` + `flat_time`) in ms.
+### Gradient variables
+- [`gradient_strength`](@ref): Maximum gradient strength achieved during the pulse in kHz/um
+- [`qval`](@ref): Spatial scale on which spins will be dephased due to this pulsed gradient in rad/um (given by `δ` * `gradient_strength`).
+
+The [`bvalue`](@ref) can be constrained for multiple gradient pulses.
+"""
+mutable struct PulsedGradient <: BuildingBlock
+    builder::SequenceBuilder
+    orientation :: Any
+    slew_rate :: VariableRef
+    rise_time :: VariableRef
+    flat_time :: VariableRef
+end
+
+function PulsedGradient(; kwargs...)
+    return BuildingBlockPlaceholder{PulsedGradient}(; kwargs...)
+end
+
+function PulsedGradient(builder::SequenceBuilder; orientation=:bvec, kwargs...)
+    model = owner_model(builder)
+    res = PulsedGradient(
+        builder,
+        orientation,
+        @variable(model),
+        @variable(model),
+        @variable(model)
+    )
+
+    set_simple_constraints!(model, res, kwargs)
+    @constraint model flat_time(res) >= 0
+    @constraint model rise_time(res) >= 0
+    @constraint model slew_rate(res) >= 0
+    return res
+end
+
+"""
+    rise_time(pulsed_gradient)
+
+The time from 0 till the maximum gradient strength in ms.
+"""
+rise_time(pg::PulsedGradient) = pg.rise_time
+
+"""
+    flat_time(pulsed_gradient)
+
+The time spent at the maximum gradient strength in ms.
+"""
+flat_time(pg::PulsedGradient) = pg.flat_time
+
+gradient_strength(g::PulsedGradient) = rise_time(g) * slew_rate(g)
+
+slew_rate(g::PulsedGradient) = g.slew_rate
+
+"""
+    δ(pulsed_gradient)
+
+Pulse gradient duration (`rise_time + `flat_time`).  This is the effective duration of the gradient. The real duration is longer (and given by [`duration`](@ref)).
+"""
+δ(g::PulsedGradient) = rise_time(g) + flat_time(g)
+
+duration(g::PulsedGradient) = 2 * rise_time(g) + flat_time(g)
+
+"""
+    qval(gradient)
+
+q-value at the end of the gradient (rad/um).
+"""
+qval(g::PulsedGradient) = (g.orientation == :neg_bvec ? -1 : 1) * gradient_strength(g) * δ(g)
+
+
+function bval(g::PulsedGradient, qstart=0.)
+    tr = rise_time(g)
+    td = δ(g)
+    grad = gradient_strength(g)
+    return (
+        # b-value due to just the gradient
+        grad * (1//60 * tr^3 - 1//12 * tr^2 * td + 1//2 * tr * td^2 + 1//3 * td^3) + 
+        # b-value due to cross-term
+        qstart * grad * (td * (td + tr)) +
+        # b-value due to just `qstart`
+        (td + tr) * qstart^2
+    )
+end
+
+properties(::Type{<:PulsedGradient}) = [qval, δ, gradient_strength, duration, rise_time, flat_time, slew_rate]
+
+
+function to_mcmr_components(block::PulsedGradient)
+    if block.orientation == :bvec
+        rotate = true
+        qvec = [value(qval(block)), 0., 0.]
+    elseif block.orientation == :neg_bvec
+        rotate = true
+        qvec = [-value(qval(block)), 0., 0.]
+    elseif block.orientation isa AbstractVector && size(block.orientation) == (3, )
+        rotate = false
+        qvec = block.orientation .* (value(qval(block)) / norm(block.orientation))
+    else
+        error("Gradient orientation should be :bvec, :neg_bvec or a length-3 vector, not $(block.orienation)")
+    end
+    t_start = value(start_time(block))
+    t_rise = value(rise_time(block))
+    t_d = value(δ(block))
+    return MRGradients([
+        (t_start, zeros(3)),
+        (t_start + t_rise, qvec),
+        (t_start + t_d, qvec),
+        (t_start + t_d + t_rise, zeros(3)),
+    ])
+end
+
+
+end
\ No newline at end of file

src/pulses/constant_pulses.jl

+module ConstantPulses
+import JuMP: VariableRef, @constraint, @variable, value
+import ....Sequences: RFPulse
+import ...BuildingBlocks: BuildingBlock, properties, BuildingBlockPlaceholder, set_simple_constraints!, duration, to_mcmr_components
+import ...SequenceBuilders: SequenceBuilder, owner_model, start_time, end_time
+import ..Properties: flip_angle, phase, amplitude, frequency, bandwidth
+
+"""
+    ConstantPulse(; variables...)
+
+Represents an radio-frequency pulse with a constant amplitude and frequency (i.e., a rectangular function).
+
+## Properties
+- [`flip_angle`](@ref): rotation expected for on-resonance spins in degrees.
+- [`duration`](@ref): duration of the RF pulse in ms.
+- [`amplitude`](@ref): amplitude of the RF pulse in kHz.
+- [`phase`](@ref): phase at the start of the RF pulse in degrees.
+- [`frequency`](@ref): frequency of the RF pulse relative to the Larmor frequency (in kHz).
+"""
+struct ConstantPulse <: BuildingBlock
+    builder :: SequenceBuilder
+    amplitude :: VariableRef
+    duration :: VariableRef
+    phase :: VariableRef
+    frequency :: VariableRef
+end
+
+ConstantPulse(; kwargs...) = BuildingBlockPlaceholder{ConstantPulse}(; kwargs...)
+function ConstantPulse(builder::SequenceBuilder; kwargs...) 
+    model = owner_model(builder)
+    res = ConstantPulse(
+        builder,
+        @variable(model),
+        @variable(model),
+        @variable(model),
+        @variable(model)
+    )
+    @constraint model amplitude(res) >= 0
+    set_simple_constraints!(model, res, kwargs)
+    return res
+end
+
+amplitude(pulse::ConstantPulse) = pulse.amplitude
+duration(pulse::ConstantPulse) = pulse.duration
+phase(pulse::ConstantPulse) = pulse.phase
+frequency(pulse::ConstantPulse) = pulse.frequency
+flip_angle(pulse::ConstantPulse) = amplitude(pulse) * duration(pulse) * 360
+bandwidth(pulse::ConstantPulse) = 3.79098854 / duration(pulse)
+
+properties(::Type{<:ConstantPulse}) = [amplitude, duration, phase, frequency, flip_angle, bandwidth]
+
+function to_mcmr_components(block::ConstantPulse)
+    final_phase = phase(block) + duration(block) * frequency(block) * 360
+    return RFPulse(
+        value.([start_time(block), end_time(block)]),
+        value.([amplitude(block), amplitude(block)]),
+        value.([phase(block), final_phase])
+    )
+end
+
+
+end
\ No newline at end of file

src/pulses/instant_pulses.jl

+module InstantPulses
+import JuMP: @constraint, @variable, VariableRef, value
+import ....Sequences: InstantRFPulse
+import ...BuildingBlocks: BuildingBlock, properties, BuildingBlockPlaceholder, set_simple_constraints!, duration, to_mcmr_components
+import ...SequenceBuilders: SequenceBuilder, owner_model, start_time
+import ..Properties: flip_angle, phase
+
+struct InstantRFPulseBlock <: BuildingBlock
+    builder :: SequenceBuilder
+    flip_angle :: VariableRef
+    phase :: VariableRef
+end
+
+InstantRFPulseBlock(; kwargs...) = BuildingBlockPlaceholder{InstantRFPulseBlock}(; kwargs...)
+function InstantRFPulseBlock(builder::SequenceBuilder; kwargs...) 
+    model = owner_model(builder)
+    res = InstantRFPulseBlock(
+        builder,
+        @variable(model),
+        @variable(model)
+    )
+    @constraint model flip_angle(res) >= 0
+    set_simple_constraints!(model, res, kwargs)
+    return res
+end
+
+flip_angle(instant::InstantRFPulseBlock) = instant.flip_angle
+phase(instant::InstantRFPulseBlock) = instant.phase
+duration(instant::InstantRFPulseBlock) = 0.
+properties(::Type{<:InstantRFPulseBlock}) = [flip_angle, phase]
+
+
+function to_mcmr_components(block::InstantRFPulseBlock)
+    return InstantRFPulse(
+        time=value(start_time(block)),
+        flip_angle=value(flip_angle(block)),
+        phase=value(phase(block)),
+    )
+end
+
+end
\ No newline at end of file

src/pulses/properties.jl

+module Properties
+
+"""
+    flip_angle(pulse_block)
+
+The flip angle of the RF pulse in a [`BuildingBlock`](@ref) in degrees.
+"""
+function flip_angle end
+
+"""
+    phase(pulse_block)
+
+The angle of the phase at the start of the RF pulse in a [`BuildingBlock`](@ref) in degrees.
+"""
+function phase end
+
+"""
+    amplitude(pulse_block)
+
+The maximum amplitude during the RF pulse in a [`BuildingBlock`](@ref) in kHz.
+"""
+function amplitude end
+
+"""
+    frequency(pulse_block)
+
+The maximum frequency during the RF pulse in a [`BuildingBlock`](@ref) relative to the Larmor frequency in kHz.
+"""
+function frequency end
+
+
+"""
+    bandwidth(pulse_block)
+
+FWHM of the frequency content of the RF pulse in kHz.
+"""
+function bandwidth end
+
+end
\ No newline at end of file

src/pulses/pulses.jl

+module Pulses
+include("properties.jl")
+include("instant_pulses.jl")
+include("constant_pulses.jl")
+include("sinc_pulses.jl")
+
+import .Properties: flip_angle, phase, amplitude, frequency, bandwidth
+import .InstantPulses: InstantRFPulseBlock
+import .ConstantPulses: ConstantPulse
+import .SincPulses: SincPulse, N_left, N_right
+
+end
\ No newline at end of file

src/pulses/sinc_pulses.jl

+module SincPulses
+
+import JuMP: VariableRef, @constraint, @variable, value
+import QuadGK: quadgk
+import Polynomials: fit, Polynomial
+import ....Sequences: RFPulse
+import ...BuildingBlocks: BuildingBlock, properties, BuildingBlockPlaceholder, set_simple_constraints!, duration, to_mcmr_components
+import ...SequenceBuilders: SequenceBuilder, owner_model, start_time, end_time
+import ..Properties: flip_angle, phase, amplitude, frequency, bandwidth
+
+"""
+    SincPulse(; symmetric=true, max_Nlobes=nothing, apodise=true, variables...)
+
+Represents an radio-frequency pulse with a constant amplitude and frequency.
+
+## Parameters
+- `symmetric`: by default the sinc pulse will be symmetric (i.e., `N_left` == `N_right`). Set `symmetric=false` to independently control `N_left` and `N_right`.
+- `max_Nlobes`: by default the computed [`flip_angle`](@ref) is only approximated as `amplitude` * `lobe_duration`. By setting `max_Nlobes` the flip_angle will be given by the actual integral of the sinc function, which will be more accurate for small number of lobes. However, the number of lobes will not be able to exceed this `max_Nlobes`.
+- `apodise`: if true (default) applies a Hanning apodising window to the sinc pulse.
+
+## Variables
+- [`N_left`](@ref): number of zero-crossings of the sinc function before the main peak (minimum of 1).
+- [`N_right`](@ref): number of zero-crossings of the sinc function after the main peak (minimum of 1).
+- [`flip_angle`](@ref): rotation expected for on-resonance spins in degrees (only approximate if `max_Nlobes` is not set).
+- [`duration`](@ref): duration of the RF pulse in ms.
+- [`amplitude`](@ref): amplitude of the RF pulse in kHz.
+- [`phase`](@ref): phase at the start of the RF pulse in degrees.
+- [`frequency`](@ref): frequency of the RF pulse relative to the Larmor frequency (in kHz).
+- [`bandwidth`](@ref): width of the rectangular function in frequency space (in kHz). If the `duration` is short (compared with 1/`bandwidth`), this bandwidth will only be approximate.
+"""
+struct SincPulse <: BuildingBlock
+    builder :: SequenceBuilder
+    symmetric :: Bool
+    apodise :: Bool
+    nlobe_integral :: Polynomial
+    N_left :: VariableRef
+    N_right :: VariableRef
+    amplitude :: VariableRef
+    phase :: VariableRef
+    frequency :: VariableRef
+    lobe_duration :: VariableRef
+end
+
+SincPulse(; kwargs...) = BuildingBlockPlaceholder{SincPulse}(; kwargs...)
+function SincPulse(builder::SequenceBuilder; symmetric=true, max_Nlobes=nothing, apodise=true, kwargs...) 
+    model = owner_model(builder)
+    if symmetric
+        N_lobes = @variable(model, integer=true)
+        N_left_var = N_right_var = N_lobes
+    else
+        N_left_var = @variable(model, integer=true)
+        N_right_var = @variable(model, integer=true)
+    end
+    res = SincPulse(
+        builder,
+        symmetric,
+        apodise,
+        nlobe_integral_params(max_Nlobes, apodise),
+        N_left_var,
+        N_right_var,
+        @variable(model),
+        @variable(model),
+        @variable(model),
+        @variable(model)
+    )
+    @constraint model amplitude(res) >= 0
+    @constraint model N_left(res) >= 1
+    if !symmetric
+        @constraint model N_right(res) >= 1
+    end
+    set_simple_constraints!(model, res, kwargs)
+    return res
+end
+
+function normalised_function(x; apodise=false)
+    if apodise
+        return (0.54 + 0.46 * cos(π * x)) * sin(π * x) / (π * x)
+    else
+        return sin(π * x) / (π * x)
+    end
+end
+
+function nlobe_integral_params(Nlobe_max, apodise=false)
+    if isnothing(Nlobe_max)
+        return fit([1], [0.5], 0)
+    end
+    f = x -> normalised_function(x; apodise=apodise)
+    integral_values = [quadgk(f, 0, i)[1] for i in 1:Nlobe_max]
+    return fit(1:Nlobe_max, integral_values, Nlobe_max - 1)
+end
+
+amplitude(pulse::SincPulse) = pulse.amplitude
+N_left(pulse::SincPulse) = pulse.N_left
+N_right(pulse::SincPulse) = pulse.N_right
+duration(pulse::SincPulse) = (N_left(pulse) + N_right(pulse)) * lobe_duration(pulse)
+phase(pulse::SincPulse) = pulse.phase
+frequency(pulse::SincPulse) = pulse.frequency
+flip_angle(pulse::SincPulse) = (pulse.nlobe_integral(N_left(pulse)) + pulse.nlobe_integral(N_right(pulse))) * amplitude(pulse) * lobe_duration(pulse) * 360
+lobe_duration(pulse::SincPulse) = pulse.lobe_duration
+bandwidth(pulse::SincPulse) = 1 / lobe_duration(pulse)
+properties(::Type{<:SincPulse}) = [amplitude, N_left, N_right, duration, phase, frequency, flip_angle, lobe_duration, bandwidth]
+
+function to_mcmr_components(block::SincPulse)
+    normed_times = -value(N_left(block)):0.1:value(N_right(block)) + 1e-5
+    times = ((normed_times .+ value(N_left(block))) .* value(lobe_duration(block))) .+ value(start_time(block))
+    amplitudes = value(amplitude(block)) .* (normalised_function.(normed_times; apodise=block.apodise))
+    phases = (value(frequency(block)) .* value(lobe_duration(block))) .* normed_times * 360
+    return RFPulse(times, amplitudes, phases)
+end
+
+end
\ No newline at end of file

src/readouts/instant_readouts.jl

+module InstantReadouts
+import ...BuildingBlocks: BuildingBlock, BuildingBlockPlaceholder, duration, properties, to_mcmr_components
+import ...SequenceBuilders: SequenceBuilder, start_time, to_block
+
+"""
+    InstantReadout()
+
+Represents an instantaneous `Readout` of the signal.
+
+It has no parameters or properties to set.
+"""
+struct InstantReadout <: BuildingBlock
+    builder::SequenceBuilder
+end
+
+InstantReadout() = BuildingBlockPlaceholder{InstantReadout}()
+
+duration(::InstantReadout) = 0.
+properties(::Type{<:InstantReadout}) = []
+to_block(builder::SequenceBuilder, cls::Type{<:InstantReadout}) = cls(builder)
+
+to_mcmr_components(readout::InstantReadout) = Readout(value(start_time(readout)))
+end
\ No newline at end of file

src/readouts/readouts.jl

+module Readouts
+include("instant_readouts.jl")
+
+import .InstantReadouts: InstantReadout
+end
\ No newline at end of file

src/sequence_builders.jl

+module SequenceBuilders
+import JuMP: Model, owner_model, index, VariableRef, @constraint, @variable, has_values, optimize!, value, optimizer_with_attributes
+import Juniper
+import Ipopt
+import ..BuildingBlocks: BuildingBlock, BuildingBlockPlaceholder, match_blocks!, duration, apply_simple_constraint!, scanner_constraints!, to_mcmr_components, start_time, end_time
+import ...Sequences: Sequence
+import ...Scanners: Scanner
+
+"""
+    SequenceBuilder(blocks...)
+
+Defines a sequence as a series of [`BuildingBlock`](@ref) objects.
+
+After defining the blocks, the user can add one or more constraints and an objective function to the properties of the [`BuildingBlock`](@ref) objects.
+A sequence matching these constraints will be produced by calling [`solve`](@ref)(builder) or [`Sequence`](@ref)(builder).
+"""
+struct SequenceBuilder
+    model :: Model
+    scanner :: Scanner
+    blocks :: Vector{<:BuildingBlock}
+    TR :: VariableRef
+    function SequenceBuilder(model::Model, blocks...; scanner=3., TR=nothing) 
+        if scanner isa Number
+            scanner = Scanner(B0=scanner)
+        end
+        builder = new(model, scanner, BuildingBlock[], @variable(model))
+        for block in blocks
+            push!(builder.blocks, to_block(builder, block))
+        end
+        @constraint model builder.TR >= duration(builder)
+        apply_simple_constraint!(model, builder.TR, TR)
+        for block in builder.blocks 
+            scanner_constraints!(block, scanner)
+        end
+        return builder
+    end
+end
+
+function to_block(model::SequenceBuilder, placeholder::BuildingBlockPlaceholder{T}) where {T}
+    block = T(model, placeholder.args...; placeholder.kwargs...)
+    if isassigned(placeholder.concrete)
+        match_blocks!(placeholder.concrete[], block)
+    else
+        placeholder.concrete[] = block
+    end
+    return block
+end
+
+Base.getindex(model::SequenceBuilder, i::Integer)  = model.blocks[i]
+
+function SequenceBuilder(blocks...; kwargs...)
+    ipopt_opt = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
+    juniper_opt = optimizer_with_attributes(Juniper.Optimizer, "nl_solver" => ipopt_opt)
+    model = Model(juniper_opt)
+    SequenceBuilder(model, blocks...; kwargs...)
+end
+
+
+function Base.show(io::IO, builder::SequenceBuilder)
+    if has_values(builder)
+        print(io, "Solved SequenceBuilder:\n")
+        for block in builder.blocks
+            print(io, "- ", block, "\n")
+        end
+    else
+        print(io, "SequenceBuilder($(builder.blocks)) being solved by ")
+        show(io, builder.model)
+    end
+end
+
+
+Base.length(sb::SequenceBuilder) = length(sb.blocks)
+builder(bb::BuildingBlock) = bb.builder
+owner_model(bb::BuildingBlock) = owner_model(builder(bb))
+owner_model(sb::SequenceBuilder) = sb.model
+has_values(object::Union{BuildingBlock, SequenceBuilder}) = has_values(owner_model(object))
+optimize!(sb::SequenceBuilder) = optimize!(owner_model(sb))
+
+"""
+    Sequence(builder::SequenceBuilder)
+
+Creates a sequence that can be run in MCMR from the [`SequenceBuilder`](@ref).
+If the `builder` has not been optimised yet, `optimize!(builder)` will be called first.
+"""
+function Sequence(sb::SequenceBuilder)
+    if !has_values(sb)
+        optimize!(sb)
+    end
+    components = []
+    for block in sb.blocks
+        new_components = to_mcmr_components(block)
+        if new_components isa AbstractVector
+            append!(components, new_components)
+        else
+            push!(components, new_components)
+        end
+    end
+    return Sequence(scanner=sb.scanner, components=components, TR=value(TR(sb)))
+end
+
+"""
+    TR(sequence::SequenceBuilder)
+
+Return the repetition time of the sequence (in ms).
+"""
+TR(sb::SequenceBuilder) = sb.TR
+
+"""
+    index(bb::BuildingBlock)
+
+Returns the index of the [`BuildingBlock`](@ref) in the parent [`SequenceBuilder`](@ref).
+"""
+index(sb::SequenceBuilder, bb::BuildingBlock) = findfirst(isequal(bb), sb.blocks)
+index(bb::BuildingBlock) = index(builder(bb), bb)
+
+
+"""
+    duration(bb1::BuildingBlock, bb2::BuildingBlock)
+    duration(sb::SequenceBuilder, bb1::Integer, bb2::Integer)
+
+The duration of the sequence from the start of [`BuildingBlock`](@ref) `bb1` till the end of `bb2`.
+Returns an error if they are not part of the same sequence.
+If `bb1` plays after `bb2` then this will calculate the time between `bb2` and the `bb1` in the next repetition time
+(i.e., `TR(sb) + end_time(bb2) - start_time(bb1)`)
+"""
+function duration(bb1::BuildingBlock, bb2::BuildingBlock)
+    if builder(bb1) != builder(bb2)
+        error("Cannot compute duration between blocks that are part of different sequences!")
+    end
+    sb = builder(bb1)
+    return duration(sb, index(sb, bb1), index(sb, bb2))
+end
+
+function duration(sb::SequenceBuilder, index1::Integer, index2::Integer)
+    if index2 == index1
+        return duration(sb[index1])
+    elseif index2 == index1 - 1
+        return TR(sb)
+    elseif index2 > index1
+        return sum(duration(sb, i) for i in index1:index2)
+    else
+        return TR(sb) - duration(sb, index2+1, index1-1)
+    end
+end
+
+duration(sb::SequenceBuilder) = duration(sb, 1, length(sb))
+duration(sb::SequenceBuilder, index::Integer) = duration(sb[index])
+
+
+"""
+    start_time(building_block::BuildingBlock)
+    start_time(seq::SequenceBuilder, building_block::Integer)
+
+Return the start time of the given [`BuildingBlock`](@ref) within the [`SequenceBuilder`](@ref).
+You can pass on the actual [`BuildingBlock`](@ref) object or the [`SequenceBuilder`](@ref) together with the index of the `building_block`
+"""
+start_time(bb::BuildingBlock) = start_time(builder(bb), bb)
+
+function start_time(sb::SequenceBuilder, bb::BuildingBlock)
+    if builder(bb) != sb
+        error("No start time of $bb within $sb, because this BuildingBlock is not part of that sequence.")
+    end
+    return start_time(sb, index(bb))
+end
+
+start_time(sb::SequenceBuilder, index::Integer) = index == 1 ? 0. : duration(sb, 1, index - 1)
+
+"""
+    end_time(building_block::BuildingBlock)
+    end_time(seq::SequenceBuilder, building_block::Integer)
+
+Return the end time of the given [`BuildingBlock`](@ref) within the [`SequenceBuilder`](@ref).
+You can pass on the actual [`BuildingBlock`](@ref) object or the [`SequenceBuilder`](@ref) together with the index of the `building_block`
+"""
+end_time(bb::BuildingBlock) = end_time(builder(bb), bb)
+
+function end_time(sb::SequenceBuilder, bb::BuildingBlock)
+    if builder(bb) != sb
+        error("No end time of $bb within $sb, because this BuildingBlock is not part of that sequence.")
+    end
+    return end_time(sb, index(bb))
+end
+
+end_time(sb::SequenceBuilder, index::Integer) = duration(sb, 1, index)
+
+
+end

src/wait.jl

+module Wait
+import JuMP: Model, @constraint, @variable, VariableRef, owner_model
+import ..BuildingBlocks: BuildingBlock, duration, properties, apply_simple_constraint!, BuildingBlockPlaceholder, to_mcmr_components
+import ..SequenceBuilders: SequenceBuilder, to_block
+import ...Scanners: Scanner
+
+"""
+    WaitBlock(duration)
+
+An empty [`BuildingBlock`](@ref) of given `duration` (in ms).
+
+Duration can be set to one of:
+- numeric value to fix it
+- `:min` to minimise its value given any external constraints
+- `:max` to maximise its value given any external constraints
+- `nothing` to make it fully determined by external constraints and objectives
+"""
+struct WaitBlock <: BuildingBlock
+    builder :: SequenceBuilder
+    duration :: VariableRef
+end
+
+function WaitBlock(builder::SequenceBuilder, duration_constraint=nothing)
+    model = owner_model(builder)
+    res = WaitBlock(builder, @variable(model))
+    @constraint model duration(res) >= 0
+    if !isnothing(duration_constraint)
+        apply_simple_constraint!(model, duration(res), duration_constraint)
+    end
+    return res
+end
+
+WaitBlock(duration_constraint=nothing) = BuildingBlockPlaceholder{WaitBlock}(duration_constraint)
+
+to_block(builder::SequenceBuilder, time::Union{Number, Symbol, Nothing, Val{:min}, Val{:max}}) = WaitBlock(builder, time)
+
+
+properties(::Type{WaitBlock}) = [duration]
+
+duration(wb::WaitBlock) = wb.duration
+
+scanner_constraints!(::Model, ::WaitBlock, ::Scanner) = nothing
+
+to_mcmr_components(::WaitBlock) = []
+
+end
\ No newline at end of file