Swtich to model where BuildingBlocks do not know their parents

src/MRIBuilder.jl

@@ -3,40 +3,48 @@ Builds and optimises NMR/MRI sequences.
 """
 module MRIBuilder
 
+include("global_model.jl")
 include("scanners.jl")
+include("variables.jl")
 include("building_blocks.jl")
-include("sequence_builders.jl")
 include("concrete_blocks.jl")
 include("wait.jl")
+include("containers/containers.jl")
 include("gradients/gradients.jl")
 include("pulses/pulses.jl")
 include("readouts/readouts.jl")
 
-import .BuildingBlocks: BuildingBlock, scanner_constraints!
-export BuildingBlock, scanner_constraints!
+import .GlobalModel: set_model
+export set_model
 
-import .SequenceBuilders: SequenceBuilder, start_time, end_time, duration, TR
-export SequenceBuilder, start_time, end_time, duration, TR
+import .Scanners: Scanner, B0, Siemens_Connectom, Siemens_Prisma, Siemens_Terra
+export Scanner, B0, Siemens_Connectom, Siemens_Prisma, Siemens_Terra
 
-import .ConcreteBlocks: ConcreteBlock, Sequence
-export ConcreteBlock, Sequence
+import .Variables: variables, duration, start_time, end_time, flip_angle, amplitude, phase, frequency, bandwidth, N_left, N_right, qval, area_under_curve, δ, rise_time, flat_time, slew_rate, gradient_strength
+export variables, duration, start_time, end_time, flip_angle, amplitude, phase, frequency, bandwidth, N_left, N_right, qval, area_under_curve, δ, rise_time, flat_time, slew_rate, gradient_strength
+
+import .BuildingBlocks: BuildingBlock
+export BuildingBlocks
+
+import .ConcreteBlocks: ConcreteBlock, AbstractConcreteBlock
+export ConcreteBlocks, AbstractConcreteBlock
 
 import .Wait: WaitBlock
 export WaitBlock
 
-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 .Containers: Sequence
+export Sequence
 
-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 .Gradients: PulsedGradient, InstantGradientBlock
+export PulsedGradient, InstantGradientBlock
+
+import .Pulses: InstantRFPulseBlock, ConstantPulse, SincPulse
+export InstantRFPulseBlock, ConstantPulse, SincPulse
 
 import .Readouts: InstantReadout
 export InstantReadout
 
-import .Scanners: Scanner, Siemens_Connectom, Siemens_Prisma, Siemens_Terra
-export Scanner, Siemens_Connectom, Siemens_Prisma, Siemens_Terra
-
-using JuMP
+import JuMP: @constraint, @objective, objective_function, optimize!, has_values, value
 export @constraint, @objective, objective_function, optimize!, has_values, value
 
 end

src/building_blocks.jl

 module BuildingBlocks
 import JuMP: has_values, GenericVariableRef, value, Model, @constraint, @objective, owner_model, objective_function
 import Printf: @sprintf
-import ..Scanners: Scanner, gradient_strength, slew_rate
+import ..Scanners: Scanner
+import ..Variables: variables, start_time, duration, end_time, gradient_strength, slew_rate
 
 """
 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_concrete_block`](@ref)(sequence, block): converts the block into a `ConcreteBlock`, which will be part of given `Sequence`.
+- `to_concrete_block`(sequence, block): converts the block into a `ConcreteBlock`, which will be part of given `Sequence`.
 - [`variables`](@ref): A list of all functions that are used to compute variables of the building block. Any of these can be used in constraints or objective functions.
 """
 abstract type BuildingBlock end
 
 
 """
-    duration(building_block)
+    start_time(container, args...)
 
-The duration of the building block in ms.
+Returns the starting time of the specific [`BuildingBlock`](@ref) within the container.
+The [`BuildingBlock`](@ref) is defined by one or more indices as defined below.
 """
-function duration end
+start_time(bb) = 0.
 
-
-"""
-    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
+    end_time(container, args...)
 
+Returns the end time of the specific [`BuildingBlock`](@ref) within the container.
+The [`BuildingBlock`](@ref) is defined by one or more indices as defined below.
 """
-    slew_rate(gradient)
+end_time(bb) = duration(bb)
 
-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)
@@ -65,24 +57,14 @@ Returns a list of function that can be called to constrain the `building_block`.
 """
 variables(bb::BuildingBlock) = variables(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
+    internal_print(io, block)
+    print(io, ")")
+end
+
+function internal_print(io::IO, block::BuildingBlock)
     for name in propertynames(block)
         value = getproperty(block, name)
         if value isa GenericVariableRef || name == :parent || string(name)[1] == '_'
@@ -101,14 +83,12 @@ function Base.show(io::IO, block::BuildingBlock)
             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)
 
@@ -158,42 +138,4 @@ function match_blocks!(block1::BuildingBlock, block2::BuildingBlock)
     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/concrete_blocks.jl

 module ConcreteBlocks
 import JuMP: has_values, optimize!, value
-import ..BuildingBlocks: BuildingBlock, BuildingBlockPlaceholder, properties, duration
-import ..SequenceBuilders: SequenceBuilder, to_block, AbstractSequence, TR, get_blocks
+import ..Variables: variables, duration
+import ..BuildingBlocks: BuildingBlock
 
-abstract type AbstractConcreteBlock <: BuildingBlock end
 
 struct ConcreteRFPulse
     time :: Vector{Number}
@@ -64,13 +63,14 @@ ConcreteGradient(values::Tuple{<:Vector, <:Vector}) = ConcreteGradient(values...
 ConcreteGradient(values::Tuple{<:Vector, <:Vector, <:Vector, <:Vector}) = ConcreteGradient(values...)
 
 
+abstract type AbstractConcreteBlock <: BuildingBlock end
+
 """
     ConcreteBlock(duration; pulse=nothing, gradient=nothing, rotate_bvec=false, readout_times=nothing)
 
 A [`BuildingBlock`](@ref) that is fully defined (i.e., there are no variables to be optimised).
 """
 struct ConcreteBlock <: AbstractConcreteBlock
-    builder :: AbstractSequence
     duration :: Float64
     pulse :: Union{ConcreteRFPulse, Nothing}
     gradient :: Union{ConcreteGradient, Nothing}
@@ -78,31 +78,29 @@ struct ConcreteBlock <: AbstractConcreteBlock
     rotate_gradient :: Bool
 end
 
-ConcreteBlock(args...; kwargs...) = BuildingBlockPlaceholder{ConcreteBlock}(args...; kwargs...)
-
-function ConcreteBlock(builder::AbstractSequence, duration::Number; pulse=nothing, gradient=nothing, readout_times=Number[], rotate_gradient=false)
-    ConcreteBlock(builder, duration, ConcreteRFPulse(pulse), ConcreteGradient(gradient), Float64.(readout_times), rotate_gradient)
+function ConcreteBlock(duration::Number; pulse=nothing, gradient=nothing, readout_times=Number[], rotate_gradient=false)
+    ConcreteBlock(duration, ConcreteRFPulse(pulse), ConcreteGradient(gradient), Float64.(readout_times), rotate_gradient)
 end
 
 
-has_values(c::AbstractConcreteBlock) = has_values(c.builder)
+has_values(c::AbstractConcreteBlock) = true
 duration(c::AbstractConcreteBlock) = 0.
 duration(c::ConcreteBlock) = c.duration
 
 
 """
-    ConcreteBlock(sequence, other_building_block)
+    ConcreteBlock(other_building_block)
 
 Creates a [`ConcreteBlock`](@ref) from another [`BuildingBlock`](@ref).
 
 This will raise an error if the other [`BuildingBlock`](@ref) has not been optimised yet.
 If it has been optimised, then [`to_concrete_block`](@ref) will be called.
 """
-function ConcreteBlock(sequence::AbstractSequence, block::BuildingBlock)
+function ConcreteBlock(block::BuildingBlock)
     if !has_values(block)
         error("Making `BuildingBlock` objects concrete is only possible after optimisation.")
     end
-    return to_concrete_block(sequence, block)
+    return to_concrete_block(block)
 end
 
 
@@ -113,40 +111,11 @@ Internal function used to create [`ConcreteBlock`](@ref) from any [`BuildingBloc
 
 This needs to be defined for every [`BuildingBlock`](@ref)
 """
-function to_concrete_block(sequence::AbstractSequence, cb::ConcreteBlock)
-    return ConcreteBlock(sequence, cb.duration, cb.pulse, cb.gradient, cb.readout_times)
+function to_concrete_block(cb::ConcreteBlock)
+    return ConcreteBlock(cb.duration, cb.pulse, cb.gradient, cb.readout_times, cb.rotate_gradient, value.(cb.start_time))
 end
 
 variables(::Type{<:AbstractConcreteBlock}) = []
 
 
-"""
-    Sequence(builder::SequenceBuilder)
-
-A fully defined sequence with no free variables.
-
-When created from a [`SequenceBuilder`](@ref), all non-fixed variables are optimised given any constraints
-and the resulting sequence is returned.
-"""
-struct Sequence <: AbstractSequence
-    blocks :: Vector{<:AbstractConcreteBlock}
-    TR :: Number
-end
-
-TR(seq::Sequence) = seq.TR
-get_blocks(seq::Sequence) = seq.blocks
-has_values(seq::Sequence) = true
-
-function Sequence(builder::SequenceBuilder)
-    if !has_values(builder)
-        optimize!(builder.model)
-    end
-    seq = Sequence(AbstractConcreteBlock[], value(TR(builder)))
-    for block in builder.blocks
-        @show block
-        push!(seq.blocks, ConcreteBlock(seq, block))
-    end
-    return seq
-end
-
 end
\ No newline at end of file

src/containers/containers.jl

+module Containers
+include("sequences.jl")
+
+import .Sequences: Sequence
+end
\ No newline at end of file

src/containers/sequences.jl

+module Sequences
+import JuMP: Model
+import ...GlobalModel: @global_model_constructor
+import ...Variables: variables, start_time, duration, VariableType
+import ...BuildingBlocks: BuildingBlock
+import ...InsertedBuildingBlocks: InsertedBuildingBlock
+
+"""
+    Sequence(building_blocks...)
+    Sequence([building_blocks])
+
+Represents a series of [`BuildingBlock`](@ref) objects run in turn.
+"""
+struct Sequence <: BuildingBlock
+    model :: Model
+    blocks :: Vector{<:BuildingBlock}
+end
+
+@global_model_constructor Sequence
+
+Sequence(model::Model, blocks...) = Sequence(model, blocks)
+
+Base.length(seq::Sequence) = length(seq)
+Base.getindex(seq::Sequence, index) = seq[index]
+
+"""
+    start_time(sequence::Sequence, index::Integer, args...)
+
+Returns the starting time of the [`BuildingBlock`](@ref) with index `index`.
+Additional `args` can be used to select a sub-block of that [`BuildingBlock`](@ref).
+The starting time is returned with respect to the start of this sequence.
+"""
+start_time(seq::Sequence) = 0.
+start_time(seq::Sequence, index::Integer) = iszero(index) ? start_time(seq) : (start_time(seq, index-1) + duration(seq[index]))
+start_time(seq::Sequence, index::Integer, args...) = start_time(seq, index) + start_time(seq[index], args...)
+
+"""
+    end_time(sequence::Sequence, index::Integer, args...)
+
+Returns the end time of the [`BuildingBlock`](@ref) with index `index`.
+Additional `args` can be used to select a sub-block of that [`BuildingBlock`](@ref).
+The end time is returned with respect to the start of this sequence.
+"""
+end_time(seq::Sequence, index::Integer) = start_time(seq, index) + duration(seq[index])
+end_time(seq::Sequence, index::Integer, args...) = start_time(seq, index) + end_time(seq[index], args...)
+
+duration(seq::Sequence) = end_time(seq, length(seq))
+
+
+
+end
+
+
+
+

src/global_model.jl

+module GlobalModel
+import JuMP: Model
+
+const GLOBAL_MODEL = Ref(Model())
+const IGNORE_MODEL = GLOBAL_MODEL[]
+
+"""
+Sets a global JuMP `Model`.
+
+Use as 
+```julia
+model = set_model([model]) do
+    ...
+end
+```
+Any sequences created within the code block will be assigned the same JuMP `Model`.
+If no model is provided a new one is created (using a Juniper optimizer based on the Ipopt non-linear optimizer).
+
+The function will return the fully formed JuMP `Model`.
+"""
+function set_model(f::Function, model::Model)
+    prev_model = global_model[]
+    global_model[] = model
+    try
+        f()
+    finally
+        global_model[] = prev_model
+    end
+    return model
+end
+
+function set_model(f::Function)
+    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)
+    set_model(f, model)
+end
+
+
+function get_global_model()
+    if GLOBAL_MODEL[] == IGNORE_MODEL
+        error("No global model has been set. Please explicitly set one in the constructor or set a global model using `set_model`.")
+    end
+    return GLOBAL_MODEL[]
+end
+
+
+"""
+    @global_model_constructor BuildingBlockType
+
+Add a onstructor the [`BuildingBlock`](@ref) subtype that fetches the global JuMP model (set by [`set_model`](@ref)) and assigns it to the first argument.
+```
+BuildingBlockType(args...; kwargs...) = BuildingBlockType(global_model::JuMP.Model, args...; kwargs...)
+```
+"""
+macro global_model_constructor(bb)
+    quote
+        $(esc(bb))(args...; kwargs...) = $(esc(bb))(get_global_model(), args...; kwargs...)
+    end
+end
+
+
+end
\ No newline at end of file

src/gradients/gradients.jl

@@ -3,7 +3,6 @@ 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 .PulsedGradients: PulsedGradient
 import .InstantGradients: InstantGradientBlock
 end
\ No newline at end of file

src/gradients/instant_gradients.jl

 module InstantGradients
-import JuMP: @constraint, @variable, VariableRef
-import ...BuildingBlocks: BuildingBlock, properties, BuildingBlockPlaceholder, set_simple_constraints!, duration
-import ...SequenceBuilders: SequenceBuilder, owner_model, start_time
-import ...ConcreteBlocks: to_concrete_block, AbstractConcreteBlock, Sequence, AbstractSequence
-import ..IntegrateGradients: qval, bval
+import JuMP: @constraint, @variable, Model, owner_model
+import ...Variables: qval, bval, start_time, duration, variables, get_free_variable, VariableType
+import ...BuildingBlocks: BuildingBlock
+import ...ConcreteBlocks: to_concrete_block, AbstractConcreteBlock
+import ...GlobalModel: @global_model_constructor
 
 """
     InstantGradientBlock(; orientation=:bvec, qval=nothing)
@@ -19,21 +19,20 @@ This is a [`BuildingBlock`](@ref) for the [`SequenceBuilder`](@ref).
 - [`qval`](@ref): Spatial scale on which spins will be dephased due to this pulsed gradient in rad/um.
 """
 struct InstantGradientBlock <: BuildingBlock
-    builder::SequenceBuilder
+    model::Model
     orientation :: Any
-    qval :: VariableRef
+    qval :: VariableType
 end
 
-InstantGradientBlock(; kwargs...) = BuildingBlockPlaceholder{InstantGradientBlock}(; kwargs...)
+@global_model_constructor InstantGradientBlock
 
-function InstantGradientBlock(builder::SequenceBuilder; orientation=:bvec, kwargs...)
-    model = owner_model(builder)
+function InstantGradientBlock(model::Model; orientation=:bvec, qval=nothing)
     res = InstantGradientBlock(
-        builder,
+        model,
         orientation,
-        @variable(model)
+        get_free_variable(model, qval),
     )
-    set_simple_constraints!(model, res, kwargs)
+    @constraint model model.qval >= 0
     return res
 end
 
@@ -51,13 +50,12 @@ Instantaneous MR gradient with no free variables.
 See [`InstantGradientBlock`](@ref) for a version where [`qval`](@ref) is variable.
 """
 struct ConcreteInstantGradient <: AbstractConcreteBlock
-    builder :: AbstractSequence
     orientation :: Any
     qval :: Number
 end
 
-function to_concrete_block(builder::Sequence, block::InstantGradientBlock)
-    return ConcreteInstantGradient(builder, block.orientation, value(qval(block)))
+function to_concrete_block(block::InstantGradientBlock)
+    return ConcreteInstantGradient(block.orientation, value(qval(block)))
 end
 
 

src/gradients/integrate_gradients.jl

 module IntegrateGradients
+import ...Variables: qval, bval
 import ...BuildingBlocks: BuildingBlock
-import ...SequenceBuilders: SequenceBuilder, TR, duration, builder, index
+import ...Containers: Sequence
 
 
 """
-    qval(blocks)
-    qval(builder::SequenceBuilder, indices)
+    qval(container, indices)
 
 Computes the total q-value summed over multiple gradient [`BuildingBlock`](@ref) objects.
 
@@ -19,12 +19,10 @@ In addition, in this interface one can provide negative indices to indicate that
 
 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]
+qval(builder::Sequence, indices::AbstractVector) = full_integral(builder, indices)[1]
 
 """
-    bval(blocks)
-    bval(builder::SequenceBuilder, indices)
+    bval(container, indices)
 
 Computes the total b-value combined over multiple gradient [`BuildingBlock`](@ref) objects.
 
@@ -38,20 +36,9 @@ In addition, in this interface one can provide negative indices to indicate that
 
 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]
+bval(builder::Sequence, indices::AbstractVector) = full_integral(builder, 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)
+function full_integral(builder::Sequence, indices::AbstractVector)
     qval_current = 0.
     current_index = 0
     bval_current = 0.

src/gradients/pulsed_gradients.jl

@@ -5,10 +5,10 @@ module PulsedGradients
 
 import JuMP: @constraint, @variable, Model, VariableRef, owner_model, value
 import StaticArrays: SVector
-import ...BuildingBlocks: BuildingBlock, duration, properties, set_simple_constraints!, BuildingBlockPlaceholder, gradient_strength, slew_rate
-import ...SequenceBuilders: SequenceBuilder, start_time
-import ...ConcreteBlocks: ConcreteBlock, to_concrete_block, AbstractSequence
-import ..IntegrateGradients: qval, bval
+import ...Variables: qval, bval, rise_time, flat_time, slew_rate, gradient_strength, variables, duration, δ, get_free_variable, VariableType
+import ...BuildingBlocks: BuildingBlock, duration, set_simple_constraints!
+import ...ConcreteBlocks: ConcreteBlock, to_concrete_block
+import ...GlobalModel: @global_model_constructor
 
 
 """
@@ -36,66 +36,36 @@ If not set, they will be determined during the sequence optimisation.
 The [`bvalue`](@ref) can be constrained for multiple gradient pulses.
 """
 mutable struct PulsedGradient <: BuildingBlock
-    builder::SequenceBuilder
+    model :: Model
     orientation :: Any
-    slew_rate :: VariableRef
-    rise_time :: VariableRef
-    flat_time :: VariableRef
+    slew_rate :: VariableType
+    rise_time :: VariableType
+    flat_time :: VariableType
 end
 
-function PulsedGradient(; kwargs...)
-    return BuildingBlockPlaceholder{PulsedGradient}(; kwargs...)
-end
+@global_model_constructor PulsedGradient
 
-function PulsedGradient(builder::SequenceBuilder; orientation=:bvec, kwargs...)
+function PulsedGradient(model::Model; orientation=:bvec, slew_rate=nothing, rise_time=nothing, flat_time=nothing, kwargs...)
     model = owner_model(builder)
     res = PulsedGradient(
         builder,
         orientation,
-        @variable(model),
-        @variable(model),
-        @variable(model)
+        [get_free_variable(model, value) for value in (slew_rate, rise_time, flat_time)]...
     )
 
     set_simple_constraints!(model, res, kwargs)
-    @constraint model flat_time(res) >= 0
-    @constraint model rise_time(res) >= 0
-    @constraint model slew_rate(res) >= 0
+    @constraint model res.flat_time >= 0
+    @constraint model res.rise_time >= 0
+    @constraint model res.slew_rate >= 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)
 
 
@@ -116,7 +86,7 @@ end
 variables(::Type{<:PulsedGradient}) = [qval, δ, gradient_strength, duration, rise_time, flat_time, slew_rate]
 
 
-function to_concrete_block(s::AbstractSequence, block::PulsedGradient)
+function to_concrete_block(block::PulsedGradient)
     if block.orientation == :bvec
         rotate = true
         qvec = [value(qval(block)), 0., 0.]
@@ -131,7 +101,7 @@ function to_concrete_block(s::AbstractSequence, block::PulsedGradient)
     end
     t_rise = value(rise_time(block))
     t_d = value(δ(block))
-    return ConcreteBlock(s, t_d + t_rise, gradient=[
+    return ConcreteBlock(t_d + t_rise, gradient=[
         (0., zeros(3)),
         (t_rise, qvec),
         (t_d, qvec),

src/inserted_building_blocks.jl

+module InsertedBuildingBlocks
+import Printf: @sprintf
+import JuMP: owner_model, has_values, Model
+import ..GlobalModel: @global_model_constructor
+import ..Variables: Variables, duration, start_time, variables, end_time, get_free_variable
+import ..BuildingBlocks: BuildingBlock, internal_print
+
+"""
+    InsertedBuildingBlock(building_block; start_time=nothing)
+    InsertedBuildingBlock(building_block, index)
+
+Represents a specific insertion of the [`BuildingBlock`](@ref) object within a larger sequence.
+
+If `index` is set an existing [`InsertedBuildingBlock`](@ref) is returned. Otherwise, a new one is created with the given `start_time` as constraint.
+"""
+struct InsertedBuildingBlock{T<:BuildingBlock}
+    bb :: T
+    index :: Int
+    function InsertedBuildingBlock(bb, index)
+        if index < 1 || index > length(bb)
+            error("$index is out of range for $bb")
+        end
+        return new{typeof(bb)}(bb, index)
+    end
+end
+
+owner_model(inserted::InsertedBuildingBlock) = owner_model(inserted.bb)
+has_values(inserted::InsertedBuildingBlock) = has_values(owner_model(inserted))
+
+@global_model_constructor InsertedBuildingBlock
+
+function InsertedBuildingBlock(model::Model, bb::BuildingBlock; start_time=nothing)
+    index = length(bb.start_time) + 1
+    push!(bb.start_time, get_free_variable(model, start_time))
+    InsertedBuildingBlock(bb, index)
+end
+
+
+function Base.show(io::IO, inserted::InsertedBuildingBlock)
+    print(io, string(typeof(block)), "(")
+    if has_values(inserted)
+        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
+    internal_print(io, inserted.bb)
+    print(io, ")")
+end
+
+for func in variables()
+    if func in (start_time, end_time)
+        continue
+    end
+    @eval Variables.$(nameof(func))(inserted::InsertedBuildingBlock, args...; kwargs...) = Variables.$(nameof(func))(inserted.bb, args...; kwargs...)
+end
+
+
+end
\ No newline at end of file

src/pulses/constant_pulses.jl

 module ConstantPulses
-import JuMP: VariableRef, @constraint, @variable, value
-import ...BuildingBlocks: BuildingBlock, properties, BuildingBlockPlaceholder, set_simple_constraints!, duration
-import ...SequenceBuilders: SequenceBuilder, owner_model, start_time, end_time, AbstractSequence
+import JuMP: VariableRef, @constraint, @variable, value, Model
+import ...BuildingBlocks: BuildingBlock, set_simple_constraints!
 import ...ConcreteBlocks: ConcreteBlock, to_concrete_block
-import ..Properties: flip_angle, phase, amplitude, frequency, bandwidth
+import ...Variables: variables, get_free_variable, flip_angle, phase, amplitude, frequency, bandwidth, start_time, end_time, VariableType, duration
+import ...GlobalModel: @global_model_constructor
 
 """
     ConstantPulse(; variables...)
@@ -18,24 +18,21 @@ Represents an radio-frequency pulse with a constant amplitude and frequency (i.e
 - [`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
+    model :: Model
+    amplitude :: VariableType
+    duration :: VariableType
+    phase :: VariableType
+    frequency :: VariableType
 end
 
-ConstantPulse(; kwargs...) = BuildingBlockPlaceholder{ConstantPulse}(; kwargs...)
-function ConstantPulse(builder::SequenceBuilder; kwargs...) 
-    model = owner_model(builder)
+@global_model_constructor ConstantPulse
+
+function ConstantPulse(model::Model; amplitude=nothing, duration=nothing, phase=nothing, frequency=nothing, kwargs...) 
     res = ConstantPulse(
-        builder,
-        @variable(model),
-        @variable(model),
-        @variable(model),
-        @variable(model)
+        model,
+        [get_free_variable(model, value) for value in (amplitude, duration, phase, frequency)]...
     )
-    @constraint model amplitude(res) >= 0
+    @constraint model res.amplitude >= 0
     set_simple_constraints!(model, res, kwargs)
     return res
 end
@@ -49,10 +46,10 @@ bandwidth(pulse::ConstantPulse) = 3.79098854 / duration(pulse)
 
 variables(::Type{<:ConstantPulse}) = [amplitude, duration, phase, frequency, flip_angle, bandwidth]
 
-function to_concrete_block(s::AbstractSequence, block::ConstantPulse)
+function to_concrete_block(block::ConstantPulse)
     d = value(duration(block))
     final_phase = value(phase(block)) + d * value(frequency(block)) * 360
-    return ConcreteBlock(s, value(duration(block)), pulse=[
+    return ConcreteBlock(value(duration(block)), pulse=[
         ([0., d]),
         value.([amplitude(block), amplitude(block)]),
         value.([phase(block), final_phase])

src/pulses/instant_pulses.jl

 module InstantPulses
-import JuMP: @constraint, @variable, VariableRef, value
-import ...BuildingBlocks: BuildingBlock, properties, BuildingBlockPlaceholder, set_simple_constraints!, duration
-import ...SequenceBuilders: SequenceBuilder, owner_model, start_time
-import ...ConcreteBlocks: to_concrete_block, AbstractConcreteBlock, Sequence, AbstractSequence
-import ..Properties: flip_angle, phase
+import JuMP: @constraint, @variable, VariableRef, value, Model
+import ...BuildingBlocks: BuildingBlock
+import ...ConcreteBlocks: to_concrete_block, AbstractConcreteBlock
+import ...Variables: flip_angle, phase, start_time, variables, duration, get_free_variable, VariableType
+import ...GlobalModel: @global_model_constructor
 
 struct InstantRFPulseBlock <: BuildingBlock
-    builder :: SequenceBuilder
-    flip_angle :: VariableRef
-    phase :: VariableRef
+    model :: Model
+    flip_angle :: VariableType
+    phase :: VariableType
 end
 
-InstantRFPulseBlock(; kwargs...) = BuildingBlockPlaceholder{InstantRFPulseBlock}(; kwargs...)
-function InstantRFPulseBlock(builder::SequenceBuilder; kwargs...) 
-    model = owner_model(builder)
+@global_model_constructor InstantRFPulseBlock
+
+function InstantRFPulseBlock(model::Model; flip_angle=nothing, phase=nothing) 
     res = InstantRFPulseBlock(
-        builder,
-        @variable(model),
-        @variable(model)
+        model,
+        get_free_variable(model, flip_angle),
+        get_free_variable(model, phase)
     )
-    @constraint model flip_angle(res) >= 0
+    @constraint model res.flip_angle >= 0
     set_simple_constraints!(model, res, kwargs)
     return res
 end
@@ -38,13 +38,12 @@ Instantaneous RF pulse with no free variables.
 See [`InstantRFPulseBlock`](@ref) for a version where [`flip_angle`](@ref) and [`phase`](@ref) are variables.
 """
 struct ConcreteInstantRFPulse <: AbstractConcreteBlock
-    builder :: AbstractSequence
     flip_angle :: Number
     phase :: Number
 end
 
-function to_concrete_block(builder::Sequence, block::InstantRFPulseBlock)
-    return ConcreteInstantRFPulse(builder, value(flip_angle(block)), value(phase(block)))
+function to_concrete_block(block::InstantRFPulseBlock)
+    return ConcreteInstantRFPulse(value(flip_angle(block)), 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
+import .SincPulses: SincPulse
 
 end
\ No newline at end of file

src/pulses/sinc_pulses.jl

 module SincPulses
 
-import JuMP: VariableRef, @constraint, @variable, value
+import JuMP: VariableRef, @constraint, @variable, value, Model
 import QuadGK: quadgk
 import Polynomials: fit, Polynomial
-import ...BuildingBlocks: BuildingBlock, properties, BuildingBlockPlaceholder, set_simple_constraints!, duration
-import ...SequenceBuilders: SequenceBuilder, owner_model, start_time, end_time, AbstractSequence
+import ...BuildingBlocks: BuildingBlock, set_simple_constraints!
 import ...ConcreteBlocks: ConcreteBlock, to_concrete_block
-import ..Properties: flip_angle, phase, amplitude, frequency, bandwidth
+import ...Variables: flip_angle, phase, amplitude, frequency, bandwidth, VariableType, variables, get_free_variable, duration
+import ...GlobalModel: @global_model_constructor
 
 """
     SincPulse(; symmetric=true, max_Nlobes=nothing, apodise=true, variables...)
@@ -29,27 +29,32 @@ Represents an radio-frequency pulse with a constant amplitude and frequency.
 - [`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
+    model :: Model
     symmetric :: Bool
     apodise :: Bool
     nlobe_integral :: Polynomial
-    N_left :: VariableRef
-    N_right :: VariableRef
-    amplitude :: VariableRef
-    phase :: VariableRef
-    frequency :: VariableRef
-    lobe_duration :: VariableRef
+    N_left :: VariableType
+    N_right :: VariableType
+    amplitude :: VariableType
+    phase :: VariableType
+    frequency :: VariableType
+    lobe_duration :: VariableType
 end
 
-SincPulse(; kwargs...) = BuildingBlockPlaceholder{SincPulse}(; kwargs...)
-function SincPulse(builder::SequenceBuilder; symmetric=true, max_Nlobes=nothing, apodise=true, kwargs...) 
-    model = owner_model(builder)
+@global_model_constructor SincPulse
+
+function SincPulse(model::Model; 
+    symmetric=true, max_Nlobes=nothing, apodise=true, N_lobes=nothing, N_left=nothing, N_right=nothing, 
+    amplitude=nothing, phase=nothing, frequency=nothing, lobe_duration=nothing, kwargs...
+) 
     if symmetric
-        N_lobes = @variable(model, integer=true)
+        N_lobes = get_free_variable(model, N_lobes)
+        @assert isnothing(N_left) && isnothing(N_right) "N_left and N_right cannot be set if symmetric=true (default)"
         N_left_var = N_right_var = N_lobes
     else
-        N_left_var = @variable(model, integer=true)
-        N_right_var = @variable(model, integer=true)
+        @assert isnothing(N_lobes) "N_lobes cannot be set if symmetric=true (default)"
+        N_left_var = get_free_variable(model, N_left)
+        N_right_var = get_free_variable(model, N_right)
     end
     res = SincPulse(
         builder,
@@ -58,15 +63,12 @@ function SincPulse(builder::SequenceBuilder; symmetric=true, max_Nlobes=nothing,
         nlobe_integral_params(max_Nlobes, apodise),
         N_left_var,
         N_right_var,
-        @variable(model),
-        @variable(model),
-        @variable(model),
-        @variable(model)
+        [get_free_variable(model, value) for value in (amplitude, phase, frequeny, lobe_duration)]...
     )
-    @constraint model amplitude(res) >= 0
-    @constraint model N_left(res) >= 1
+    @constraint model res.amplitude >= 0
+    @constraint model res.N_left >= 1
     if !symmetric
-        @constraint model N_right(res) >= 1
+        @constraint model res.N_right >= 1
     end
     set_simple_constraints!(model, res, kwargs)
     return res
@@ -100,13 +102,13 @@ lobe_duration(pulse::SincPulse) = pulse.lobe_duration
 bandwidth(pulse::SincPulse) = 1 / lobe_duration(pulse)
 variables(::Type{<:SincPulse}) = [amplitude, N_left, N_right, duration, phase, frequency, flip_angle, lobe_duration, bandwidth]
 
-function to_concrete_block(s::AbstractSequence, block::SincPulse)
+function to_concrete_block(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)))
     amplitudes = value(amplitude(block)) .* (normalised_function.(normed_times; apodise=block.apodise))
     phases = (value(frequency(block)) .* value(lobe_duration(block))) .* normed_times * 360
     return ConcreteBlock(
-        s, value(duration(block)); 
+        value(duration(block)); 
         pulse=(times, amplitudes, phases)
     )
 end

src/readouts/instant_readouts.jl

 module InstantReadouts
-import ...BuildingBlocks: BuildingBlock, BuildingBlockPlaceholder, duration, properties
-import ...SequenceBuilders: SequenceBuilder, start_time, to_block, AbstractSequence
+import ...BuildingBlocks: BuildingBlock
 import ...ConcreteBlocks: AbstractConcreteBlock, to_concrete_block
+import ...Variables: variables
 
 """
     InstantReadout()
@@ -11,13 +11,8 @@ Represents an instantaneous `Readout` of the signal.
 It has no parameters or properties to set.
 """
 struct InstantReadout <: AbstractConcreteBlock
-    builder::AbstractSequence
 end
 
-InstantReadout() = BuildingBlockPlaceholder{InstantReadout}()
-
 variables(::Type{<:InstantReadout}) = []
-to_block(builder::SequenceBuilder, cls::Type{<:InstantReadout}) = cls(builder)
-to_concrete_block(builder::AbstractSequence, ::InstantReadout) = InstantReadout(builder)
-
+to_concrete_block(::InstantReadout) = InstantReadout()
 end
\ No newline at end of file

src/variables.jl

+module Variables
+import JuMP: @variable, Model, @objective, objective_function, owner_model, has_values, value, AbstractJuMPScalar
+import ..Scanners: gradient_strength, slew_rate
+
+all_variables_symbols = [
+    # general
+    :duration => (:block, "duration of the building block in ms."),
+
+    # RF pulse
+    :flip_angle => (:pulse, "The flip angle of the RF pulse in degrees"),
+    :amplitude => (:pulse, "The maximum amplitude of an RF pulse in kHz"),
+    :phase => (:pulse, "The angle of the phase of an RF pulse in KHz"),
+    :frequency => (:pulse, "The off-resonance frequency of an RF pulse (relative to the Larmor frequency of water) in KHz"),
+    :bandwidth => (:pulse, "Bandwidth of the RF pulse in kHz"),
+    :N_left => (:pulse, "The number of zero crossings of the RF pulse before the main peak"),
+    :N_right => (:pulse, "The number of zero crossings of the RF pulse after the main peak"),
+
+    # gradients
+    :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))."),
+    :δ => (: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."),
+    :slew_rate => (:gradient, "maximum increase of a gradient (kHz/um/ms)"),
+    :gradient_strength => (:gradient, "maximum strength of a gradient (kHz/um)"),
+]
+
+symbol_to_func = Dict{Symbol, Function}()
+
+
+
+for (func_symbol, (block_symbol, description)) in all_variables_symbols
+    as_string = "    $func_symbol($block_symbol)\n\n$description\n\nThis represents a variable within the sequence. Variables can be set during the construction of a [`BuildingBlock`](@ref) or used to create constraints after the fact."
+    new_func = @eval begin
+        $as_string
+        function $func_symbol end
+    end
+    symbol_to_func[func_symbol] = new_func
+end
+
+
+"""
+    variables(building_block)
+    variables()
+
+Returns all functions representing properties of a [`BuildingBlock`](@ref) object.
+"""
+variables() = [values(symbol_to_func)...]
+
+
+# Some universal truths
+area_under_curve(bb) = qval(bb)
+
+
+function start_time end
+function end_time end
+
+
+const VariableType = Union{Number, AbstractJuMPScalar}
+
+
+"""
+    get_free_variable(model, value; integer=false)
+
+Get a representation of a given `variable` given a user-defined constraint.
+"""
+get_free_variable(::Model, value::Number; integer=false) = integer ? Int(value) : Float64(value)
+function get_free_variable(model::Model, value::VariableType; integer=false)
+    if owner_model(value) != model
+        if has_values(value)
+            return value(value)
+        end
+        error("Cannot set any constraints between sequences stored in different JuMP models.")
+    end
+    return value
+end
+get_free_variable(model::Model, ::Nothing; integer=false) = @variable(model, integer=integer)
+get_free_variable(model::Model, value::Symbol; integer=false) = integer ? error("Cannot maximise or minimise an integer variable") : get_free_variable(model, Val(value))
+function get_free_variable(model::Model, ::Val{:min})
+    var = get_free_variable(model, nothing)
+    @objective model Min objective_function(model) + var
+    return var
+end
+function get_free_variable(model::Model, ::Val{:max})
+    var = get_free_variable(model, nothing)
+    @objective model Min objective_function(model) - var
+    return var
+end
+
+function bval end
+
+end
\ No newline at end of file

src/wait.jl

 module Wait
 import JuMP: Model, @constraint, @variable, VariableRef, owner_model, value
-import ..BuildingBlocks: BuildingBlock, duration, properties, apply_simple_constraint!, BuildingBlockPlaceholder
-import ..SequenceBuilders: SequenceBuilder, to_block, AbstractSequence
+import ..Variables: VariableType, variables, duration, get_free_variable
+import ..BuildingBlocks: BuildingBlock
 import ..ConcreteBlocks: to_concrete_block, ConcreteBlock
+import ..GlobalModel: @global_model_constructor
 import ...Scanners: Scanner
 
 """
@@ -17,31 +18,29 @@ Duration can be set to one of:
 - `nothing` to make it fully determined by external constraints and objectives
 """
 struct WaitBlock <: BuildingBlock
-    builder :: SequenceBuilder
-    duration :: VariableRef
+    model :: Model
+    duration :: VariableType
 end
 
-function WaitBlock(builder::SequenceBuilder, duration_constraint=nothing)
-    model = owner_model(builder)
-    res = WaitBlock(builder, @variable(model))
+@global_model_constructor WaitBlock
+
+function WaitBlock(model::Model, duration=nothing)
+    res = WaitBlock(
+        model, 
+        get_free_variable(model, duration),
+        VariableType[]
+    )
     @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)
+to_block(time::Union{VariableType, Symbol, Nothing, Val{:min}, Val{:max}}) = WaitBlock(time)
 
 
 variables(::Type{WaitBlock}) = [duration]
 
 duration(wb::WaitBlock) = wb.duration
 
-scanner_constraints!(::Model, ::WaitBlock, ::Scanner) = nothing
-
-to_concrete_block(builder::AbstractSequence, wb::WaitBlock) = ConcreteBlock(builder, value(duration(wb)))
+to_concrete_block(wb::WaitBlock) = ConcreteBlock(value(duration(wb)))
 
 end
\ No newline at end of file