make optimisation optional in `build_sequence`.

src/build_sequences.jl

@@ -14,7 +14,7 @@ Wrapper to build a sequence.
 
 Use as 
 ```julia
-build_sequence([model]) do model
+build_sequence(scanner[, optimiser_constructor];) do
     ...
 end
 ```
@@ -29,33 +29,53 @@ seq = Sequence(
 
 You can also add any arbitrary constraints or objectives using the same syntax as for [`JuMP`](https://jump.dev/JuMP.jl):
 ```
-@constraint model duration(seq) == 30.
+@constraint global_model() duration(seq) == 30.
 ```
 
 As soon as the code block is the optimal sequence matching all your constraints and objectives will be returned.
+
+## Parameters
+- `scanner`: Set to a [`Scanner`](@ref) to limit the gradient strength and slew rate. When this call to `build_sequence` is embedded in another, this parameter can be set to `nothing` to indicate that the same scanner should be used. 
+- `optimiser_constructor`: A `JuMP` solver optimiser as described in the [JuMP documentation](https://jump.dev/JuMP.jl/stable/tutorials/getting_started/getting_started_with_JuMP/#What-is-a-solver?). Defaults to using [Ipopt](https://github.com/jump-dev/Ipopt.jl).
+
+## Variables
+- `optimise`: Whether to optimise and fix the sequence as soon as it is returned. This defaults to `true` if a scanner is provided and `false` if no scanner is provided.
+- `kwargs...`: Other keywords are passed on as attributes to the `optimiser_constructor` (e.g., set `print_level=3` to make the Ipopt optimiser quieter).
 """
-function build_sequence(f::Function, scanner::Scanner, model::Model)
+function build_sequence(f::Function, scanner::Union{Nothing, Scanner}, model::Model, optimise::Bool)
     prev_model = GLOBAL_MODEL[]
     GLOBAL_MODEL[] = model
     prev_scanner = GLOBAL_SCANNER[]
     if !isnothing(scanner)
         GLOBAL_SCANNER[] = scanner
+    elseif !isfinite(gradient_strength(GLOBAL_SCANNER[]))
+        error("Scanner should be explicitly set when creating a new top-level sequence.")
     end
     try
-        sequence = f(model)
-        optimize!(model)
-        return fixed(sequence)
+        sequence = f()
+        if optimise
+            optimize!(model)
+            return fixed(sequence)
+        else
+            return sequence
+        end
     finally
         GLOBAL_MODEL[] = prev_model
         GLOBAL_SCANNER[] = prev_scanner
     end
 end
 
-function build_sequence(f::Function, scanner::Scanner)
-    ipopt_opt = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 3)
-    juniper_opt = optimizer_with_attributes(Juniper.Optimizer, "nl_solver" => ipopt_opt)
-    model = Model(ipopt_opt)
-    build_sequence(f, scanner, model)
+function build_sequence(f::Function, ::Nothing; optimise=false)
+    build_sequence(f, nothing, global_model(), optimise)
+end
+
+function build_sequence(f::Function, scanner::Scanner, optimiser_constructor; optimise=true, kwargs...)
+    model = Model(optimizer_with_attributes(optimiser_constructor, pairs(kwargs)...))
+    build_sequence(f, scanner, model, optimise)
+end
+
+function build_sequence(f::Function, scanner::Scanner; kwargs...)
+    build_sequence(f, scanner, Ipopt.Optimizer; kwargs...)
 end
 
 

src/helper_functions.jl

 module HelperFunctions
 import JuMP: @constraint
-import ..BuildSequences: global_model
+import ..BuildSequences: global_model, build_sequence
 import ..Sequences: Sequence
 import ..Pulses: SincPulse, ConstantPulse, InstantRFPulseBlock
 import ..Overlapping: TrapezoidGradient, SpoiltSliceSelect
@@ -29,6 +29,8 @@ If `min_slice_thickness` is not set or is set to `:min`, then either `bandwidth`
 
 ## Parameters
 For an [`InstantRFPulseBlock`](@ref) (i.e., `shape=:instant`), only the `flip_angle`, `phase`, and `scale` will be used. All other parameters are ignored.
+- `optimise`: set to true to optimise this RF pulse separately from the embedding sequence.
+
 ### Pulse parameters
 - `shape`: The shape of the RF pulse. One of `:sinc` (for [`SincPulse`](@ref)), `:constant`/`:hard` (for [`ConstantPulse`](@ref)), or `:instant` (for [`InstantRFPulseBlock`](@ref)).
 - `flip_angle`: size of the flip due to the RF pulse in degrees (default: 90).
@@ -43,30 +45,33 @@ For an [`InstantRFPulseBlock`](@ref) (i.e., `shape=:instant`), only the `flip_an
 - `slice_thickness`: minimum slice thickness that should be possible without adjusting the sequence timings in um (not mm!) (default: no slice selection). Can be set to `:min` to indicate that this should be minimised given the scanner constraints and user values for `bandwidth` or `duration`.
 - `rephase`: set to false to disable the spin rephasing after the RF pulse.
 - `rotate_grad`: name of the parameter with which the slice selection gradient will be rotated after sequence optimisation (default: `:FOV`).
+- `scanner`: overrides the [`global_scanner`](@ref) for this part of the sequence. Recommended to set only if not part of a larger sequence.
 """
-function excitation_pulse(; flip_angle=90, phase=0., frequency=0., shape=:sinc, slice_thickness=Inf, rephase=true, Nzero=3, scale=:transmit_B1, rotate_grad=:FOV, bandwidth=nothing, duration=nothing)
-    pulse = _get_pulse(shape, flip_angle, phase, frequency, Nzero, scale, bandwidth, duration)
-    if pulse isa InstantRFPulseBlock
-        if !isinf(slice_thickness)
-            error("An instant RF pulse always affects all spins equally, so using `shape=:instant` is incompatible with setting `slice_thickness`.")
+function excitation_pulse(; flip_angle=90, phase=0., frequency=0., shape=:sinc, slice_thickness=Inf, rephase=true, Nzero=3, scale=:transmit_B1, rotate_grad=:FOV, bandwidth=nothing, duration=nothing, scanner=nothing, optimise=false)
+    build_sequence(scanner; optimise=optimise) do
+        pulse = _get_pulse(shape, flip_angle, phase, frequency, Nzero, scale, bandwidth, duration)
+        if pulse isa InstantRFPulseBlock
+            if !isinf(slice_thickness)
+                error("An instant RF pulse always affects all spins equally, so using `shape=:instant` is incompatible with setting `slice_thickness`.")
+            end
+            return pulse
+        end
+        if isinf(slice_thickness)
+            return pulse
         end
-        return pulse
-    end
-    if isinf(slice_thickness)
-        return pulse
-    end
 
-    grad = TrapezoidGradient(pulse=pulse, duration=:min, slice_thickness=[Inf, Inf, slice_thickness], orientation=[0, 0, 1.], rotate=rotate_grad)
-    if !rephase
-        return grad
+        grad = TrapezoidGradient(pulse=pulse, duration=:min, slice_thickness=[Inf, Inf, slice_thickness], orientation=[0, 0, 1.], rotate=rotate_grad)
+        if !rephase
+            return grad
+        end
+        seq = Sequence(
+            grad,
+            TrapezoidGradient(orientation=[0, 0, -1.], rotate=rotate_grad, duration=:min);
+            TR=Inf
+        )
+        @constraint global_model() qvec(grad, 1, nothing)[3] == -qvec(seq[2])[3]
+        return seq
     end
-    seq = Sequence(
-        grad,
-        TrapezoidGradient(orientation=[0, 0, -1.], rotate=rotate_grad, duration=:min);
-        TR=Inf
-    )
-    @constraint global_model() qvec(grad, 1, nothing)[3] == -qvec(seq[2])[3]
-    return seq
 end
 
 
@@ -82,6 +87,8 @@ If `slice_thickness` is not set or is set to `:min`, then either `bandwidth` or
 
 ## Parameters
 For an [`InstantRFPulseBlock`](@ref) (i.e., `shape=:instant`), only the `flip_angle`, `phase`, and `scale` will be used. All other parameters are ignored.
+- `optimise`: set to true to optimise this RF pulse separately from the embedding sequence.
+
 ### Pulse parameters
 - `shape`: The shape of the RF pulse. One of `:sinc` (for [`SincPulse`](@ref)), `:constant`/`:hard` (for [`ConstantPulse`](@ref)), or `:instant` (for [`InstantRFPulseBlock`](@ref)).
 - `flip_angle`: size of the flip due to the RF pulse in degrees (default: 180).
@@ -96,27 +103,30 @@ For an [`InstantRFPulseBlock`](@ref) (i.e., `shape=:instant`), only the `flip_an
 - `slice_thickness`: minimum slice thickness that should be possible without adjusting the sequence timings in um (not mm!) (default: no slice selection). Can be set to `:min` to indicate that this should be minimised given the scanner constraints and user values for `bandwidth` or `duration`.
 - `spoiler`: set to the spatial scale on which the spins should be dephased in mm. For rotating spoilers, this does include the contribution from the slice select gradient as well.
 - `rotate_grad`: name of the parameter with which the slice selection and spoiler gradient will be rotated after sequence optimisation (default: `:FOV`).
+- `scanner`: overrides the [`global_scanner`](@ref) for this part of the sequence. Recommended to set only if not part of a larger sequence.
 """
-function refocus_pulse(; flip_angle=180, phase=0., frequency=0., shape=:sinc, slice_thickness=Inf, Nzero=3, scale=:transmit_B1, rotate_grad=:FOV, bandwidth=nothing, duration=nothing, spoiler=Inf)
-    pulse = _get_pulse(shape, flip_angle, phase, frequency, Nzero, scale, bandwidth, duration)
-    if pulse isa InstantRFPulseBlock && !isinf(slice_thickness)
-        error("An instant RF pulse always affects all spins equally, so using `shape=:instant` is incompatible with setting `slice_thickness`.")
-    end
-    if isinf(spoiler)
-        if isinf(slice_thickness)
-            return pulse
-        else
-            return TrapezoidGradient(pulse=pulse, duration=:min, slice_thickness=[Inf, Inf, slice_thickness], orientation=[0, 0, 1.], rotate=rotate_grad)
+function refocus_pulse(; flip_angle=180, phase=0., frequency=0., shape=:sinc, slice_thickness=Inf, Nzero=3, scale=:transmit_B1, rotate_grad=:FOV, bandwidth=nothing, duration=nothing, spoiler=Inf, scanner=nothing, optimise=false)
+    build_sequence(scanner; optimise=optimise) do
+        pulse = _get_pulse(shape, flip_angle, phase, frequency, Nzero, scale, bandwidth, duration)
+        if pulse isa InstantRFPulseBlock && !isinf(slice_thickness)
+            error("An instant RF pulse always affects all spins equally, so using `shape=:instant` is incompatible with setting `slice_thickness`.")
         end
-    else
-        orientation=isnothing(rotate_grad) ? [1, 1, 1] : [0, 0, 1]
-        if isinf(slice_thickness)
-            grad = TrapezoidGradient(orientation=orientation, duration=:min, rotate=rotate_grad)
-            @constraint global_model() qvec(grad)[3] == 2π * 1e-3 / spoiler
-            return Sequence(grad, pulse, grad; TR=Inf)
+        if isinf(spoiler)
+            if isinf(slice_thickness)
+                return pulse
+            else
+                return TrapezoidGradient(pulse=pulse, duration=:min, slice_thickness=[Inf, Inf, slice_thickness], orientation=[0, 0, 1.], rotate=rotate_grad)
+            end
         else
-            res = SpoiltSliceSelect(pulse; orientation=orientation, duration=:min, rotate=rotate_grad, slice_thickness=slice_thickness, spoiler_scale=spoiler)
-            return res
+            orientation=isnothing(rotate_grad) ? [1, 1, 1] : [0, 0, 1]
+            if isinf(slice_thickness)
+                grad = TrapezoidGradient(orientation=orientation, duration=:min, rotate=rotate_grad)
+                @constraint global_model() qvec(grad)[3] == 2π * 1e-3 / spoiler
+                return Sequence(grad, pulse, grad; TR=Inf)
+            else
+                res = SpoiltSliceSelect(pulse; orientation=orientation, duration=:min, rotate=rotate_grad, slice_thickness=slice_thickness, spoiler_scale=spoiler)
+                return res
+            end
         end
     end
 end