From 0702e164d9964d05108fd7ce73a6120056570bc6 Mon Sep 17 00:00:00 2001
From: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
Date: Wed, 31 Jan 2024 12:55:54 +0000
Subject: [PATCH] No longer distinguish between variable and fixed gradients

---
 src/gradients/changing_gradient_blocks.jl |  48 ++++-----
 src/gradients/constant_gradient_blocks.jl |  38 +++----
 src/gradients/fixed_gradients.jl          | 117 ----------------------
 src/gradients/instant_gradients.jl        |   6 +-
 src/gradients/integrate_gradients.jl      |  70 -------------
 src/gradients/slice_selects.jl            |  27 ++---
 6 files changed, 52 insertions(+), 254 deletions(-)
 delete mode 100644 src/gradients/fixed_gradients.jl
 delete mode 100644 src/gradients/integrate_gradients.jl

diff --git a/src/gradients/changing_gradient_blocks.jl b/src/gradients/changing_gradient_blocks.jl
index 9086770..6635335 100644
--- a/src/gradients/changing_gradient_blocks.jl
+++ b/src/gradients/changing_gradient_blocks.jl
@@ -6,41 +6,35 @@ import JuMP: @constraint, @variable, Model, owner_model
 import ...Variables: qvec, bmat_gradient, gradient_strength, slew_rate, duration, variables, VariableType
 import ...BuildingBlocks: GradientBlock, fixed, RFPulseBlock
 
-struct AbstractChangingGradientBlock{T<:VariableType} <: GradientBlock
-    gradient_strength_start :: SVector{3, T}
-    slew_rate :: SVector{3, T}
-    duration :: T
-    rotate :: Union{Nothing, Symbol}
-    scale :: Union{Nothing, Symbol}
-end
-
 """
-    ChangingGradientBlock(grad1, grad2, duration, rotate, scale)
+    ChangingGradientBlock(grad1, slew_rate, duration, rotate, scale)
 
 Underlying type for any linearly changing part in a gradient waveform.
 
-The corresponding type in a fixed gradient is [`FixedChangingGradientBlock`](@ref)
-
-Do not create this directly, it will be created by the gradient waveforms.
-"""
-const ChangingGradientBlock = AbstractChangingGradientBlock{VariableType}
+Usually, you do not want to create this object directly, use a gradient waveform instead.
 
+## Parameters
+- `rotate`: with which user-set parameter will this gradient be rotated (e.g., :bvec). Default is no rotation.
+- `scale`: with which user-set parameter will this gradient be scaled (e.g., :bval). Default is no scaling.
 """
-    FixedChangingGradientBlock(grad1, grad2, duration, rotate, scale)
+struct ChangingGradientBlock <: GradientBlock
+    gradient_strength_start :: SVector{3, <:VariableType}
+    slew_rate :: SVector{3, <:VariableType}
+    duration :: VariableType
+    rotate :: Union{Nothing, Symbol}
+    scale :: Union{Nothing, Symbol}
+end
 
-Fixed equivalent of [`ChangingGradientBlock`](@ref).
-"""
-const FixedChangingGradientBlock = AbstractChangingGradientBlock{Float64}
 
-duration(cgb::AbstractChangingGradientBlock) = cgb.duration
+duration(cgb::ChangingGradientBlock) = cgb.duration
 
-grad_start(cgb::AbstractChangingGradientBlock) = cgb.gradient_strength_start
-slew_rate(cgb::AbstractChangingGradientBlock) = cgb.slew_rate
-grad_end(cgb::AbstractChangingGradientBlock) = grad_start(cgb) .+ slew_rate(cgb) .* duration(cgb)
-gradient_strength(cgb::AbstractChangingGradientBlock) = max.(grad_start(cgb), grad_end(cgb))
-qvec(cgb::AbstractChangingGradientBlock) = (grad_start(cgb) .+ grad_end(cgb)) .* (duration(cgb) * π)
+grad_start(cgb::ChangingGradientBlock) = cgb.gradient_strength_start
+slew_rate(cgb::ChangingGradientBlock) = cgb.slew_rate
+grad_end(cgb::ChangingGradientBlock) = grad_start(cgb) .+ slew_rate(cgb) .* duration(cgb)
+gradient_strength(cgb::ChangingGradientBlock) = max.(grad_start(cgb), grad_end(cgb))
+qvec(cgb::ChangingGradientBlock) = (grad_start(cgb) .+ grad_end(cgb)) .* (duration(cgb) * π)
 
-function bmat_gradient(cgb::AbstractChangingGradientBlock, qstart)
+function bmat_gradient(cgb::ChangingGradientBlock, qstart)
     # grad = (g1 * (duration - t) + g2 * t) / duration
     #      = g1 + (g2 - g1) * t / duration
     # q = qstart + g1 * t + (g2 - g1) * t^2 / (2 * duration)
@@ -60,7 +54,7 @@ function bmat_gradient(cgb::AbstractChangingGradientBlock, qstart)
     )
 end
 
-function bmat_gradient(cgb::AbstractChangingGradientBlock)
+function bmat_gradient(cgb::ChangingGradientBlock)
     diff = slew_rate(cgb) .* duration(cgb)
     return (2Ï€)^2 .* (
         grad_start(cgb) .* permutedims(grad_start(cgb)) ./ 3 .+
@@ -70,6 +64,6 @@ function bmat_gradient(cgb::AbstractChangingGradientBlock)
 end
 
 
-variables(::Type{<:AbstractChangingGradientBlock}) = [duration, slew_rate]
+variables(::Type{<:ChangingGradientBlock}) = [duration, slew_rate, gradient_strength, qvec]
 
 end
diff --git a/src/gradients/constant_gradient_blocks.jl b/src/gradients/constant_gradient_blocks.jl
index 1f1a8dd..5ca15ed 100644
--- a/src/gradients/constant_gradient_blocks.jl
+++ b/src/gradients/constant_gradient_blocks.jl
@@ -6,43 +6,45 @@ import JuMP: @constraint, @variable, Model, owner_model
 import ...Variables: qvec, bmat_gradient, gradient_strength, slew_rate, duration, variables, VariableType
 import ...BuildingBlocks: GradientBlock, fixed, RFPulseBlock
 
-struct AbstractConstantGradientBlock{T<:VariableType} <: GradientBlock
-    gradient_strength :: SVector{3, T}
-    duration :: T
-    rotate :: Union{Nothing, Symbol}
-    scale :: Union{Nothing, Symbol}
-end
-
 """
     ConstantGradientBlock(gradient_strength, duration, rotate, scale)
 
 Underlying type for any flat part in a gradient waveform.
 It can overlap with an [`RFPulseBlock`], in which case it will be a `SliceSelectPulse`
 
-The corresponding type in a fixed gradient is [`FixedConstantGradientBlock`](@ref)
+Usually, you do not want to create this object directly, use a gradient waveform instead.
 
-Do not create this directly, it will be created by the gradient waveforms.
+## Parameters
+- `rotate`: with which user-set parameter will this gradient be rotated (e.g., :bvec). Default is no rotation.
+- `scale`: with which user-set parameter will this gradient be scaled (e.g., :bval). Default is no scaling.
 """
-const ConstantGradientBlock = AbstractConstantGradientBlock{VariableType}
+struct ConstantGradientBlock{T<:VariableType} <: GradientBlock
+    gradient_strength :: SVector{3, T}
+    duration :: T
+    rotate :: Union{Nothing, Symbol}
+    scale :: Union{Nothing, Symbol}
+end
+
+const ConstantGradientBlock = ConstantGradientBlock{VariableType}
 
 """
     FixedConstantGradientBlock(gradient_strength, duration, rotate, scale)
 
 Fixed equivalent of [`ConstantGradientBlock`](@ref).
 """
-const FixedConstantGradientBlock = AbstractConstantGradientBlock{Float64}
+const FixedConstantGradientBlock = ConstantGradientBlock{Float64}
 
-duration(cgb::AbstractConstantGradientBlock) = cgb.duration
-gradient_strength(cgb::AbstractConstantGradientBlock) = cgb.gradient_strength
-slew_rate(::AbstractConstantGradientBlock) = zero(SVector{3, Float64})
-qvec(cgb::AbstractConstantGradientBlock) = duration(cgb) .* gradient_strength(cgb) .* 2Ï€
+duration(cgb::ConstantGradientBlock) = cgb.duration
+gradient_strength(cgb::ConstantGradientBlock) = cgb.gradient_strength
+slew_rate(::ConstantGradientBlock) = zero(SVector{3, Float64})
+qvec(cgb::ConstantGradientBlock) = duration(cgb) .* gradient_strength(cgb) .* 2Ï€
 
-function bmat_gradient(cgb::AbstractConstantGradientBlock)
+function bmat_gradient(cgb::ConstantGradientBlock)
     grad = 2Ï€ .* gradient_strength(cgb)
     return (grad .* permutedims(grad)) .* duration(cgb)^3 ./3
 end
 
-function bmat_gradient(cgb::AbstractConstantGradientBlock, qstart)
+function bmat_gradient(cgb::ConstantGradientBlock, qstart)
     # \int dt (qstart + t * grad)^2 = 
     #   qstart^2 * duration +
     #   qstart * grad * duration^2 +
@@ -55,6 +57,6 @@ function bmat_gradient(cgb::AbstractConstantGradientBlock, qstart)
     )
 end
 
-variables(::Type{<:AbstractConstantGradientBlock}) = [duration, gradient_strength]
+variables(::Type{<:ConstantGradientBlock}) = [duration, gradient_strength, qvec]
 
 end
diff --git a/src/gradients/fixed_gradients.jl b/src/gradients/fixed_gradients.jl
deleted file mode 100644
index 4fc3e3d..0000000
--- a/src/gradients/fixed_gradients.jl
+++ /dev/null
@@ -1,117 +0,0 @@
-module FixedGradients
-
-import Printf: @sprintf
-import LinearAlgebra: norm
-import StaticArrays: SVector
-import ...BuildingBlocks: ContainerBlock, fixed, BuildingBlock, BuildingBlockPrinter, get_children_blocks, GradientBlock
-import ...Variables: variables, duration, qvec, gradient_strength, slew_rate, start_time
-import ..ChangingGradientBlocks: FixedChangingGradientBlock
-
-
-"""
-    FixedGradient(time, Gx, Gy, Gz)
-    FixedGradient(time, Gx)
-    FixedGradient(time, gradients)
-
-Create a Gradient profile that has been fully defined by N control point.
-All arguments should be arrays of the same length N defining these control points.
-
-- `time`: time since the start of this [`BuildingBlock`](@ref) in ms.
-- `Gx`: gradient in the x-direction at every timepoint.
-- `Gy`: gradient in the y-direction at every timepoint (assumed zero if not provided).
-- `Gz`: gradient in the z-direction at every timepoint (assumed zero if not provided).
-- `gradients`: length-3 array at every timepoint defining the x-, y-, and z- direction of the gradient.
-"""
-struct FixedGradient <: ContainerBlock
-    time :: Vector{Float64}
-    gradient_strength :: Vector{SVector{3, Float64}}
-    rotate :: Union{Symbol, Nothing}
-    scale :: Union{Symbol, Nothing}
-    function FixedGradient(time::AbstractVector{<:Number}, Gx::AbstractVector{<:Number}, Gy::AbstractVector{<:Number}, Gz::AbstractVector{<:Number}; rotate=nothing, scale=nothing)
-        @assert length(time) == length(Gx)
-        @assert length(time) == length(Gy)
-        @assert length(time) == length(Gz)
-        grad = [SVector{3, Float64}(Gx[i], Gy[i], Gz[i]) for i in eachindex(Gx)]
-        new(Float64.(time), grad, rotate, scale)
-    end
-end
-
-get_children_indices(fg::FixedGradient) = 1:(length(fg.time)-1)
-Base.getindex(fg::FixedGradient, index::Int) = FixedChangingGradientBlock(
-    fg.gradient_strength[index],
-    fg.gradient_strength[index + 1],
-    fg.times[index + 1] - fg.times[index],
-    rotate,
-    scale,
-)
-start_time(fg::FixedGradient, index::Int) = fg.times[index]
-
-function FixedGradient(time::AbstractVector{<:Number}, arr::AbstractVector{<:AbstractVector{<:Number}}; kwargs...) 
-    @assert all(length.(arr) .== 3)
-    FixedGradient(
-        time,
-        [a[1] for a in arr],
-        [a[2] for a in arr],
-        [a[3] for a in arr];
-        kwargs...
-    )
-end
-
-FixedGradient(time::AbstractVector{<:Number}, Gx::AbstractVector{<:Number}; kwargs...) = FixedGradient(time, Gx, zeros(length(time)), zeros(length(time)); kwargs...)
-
-variables(::Type{<:FixedGradient}) = []
-
-duration(fg::FixedGradient) = maximum(fg.time)
-
-function gradient_strength(fg::FixedGradient)
-    if isnothing(fg.rotate)
-        return maximum(map(g -> max(abs.(g)...), fg.gradient_strength))
-    else
-        return maximum(map(norm, fg.gradient_strength))
-    end
-end
-
-function slew_rate(fg::FixedGradient)
-    diff = (fg.gradient_strength[2:end] .- fg.gradient_strength[1:end-1]) ./ (fg.time[2:end] - fg.time[1:end-1])
-    if isnothing(fg.rotate)
-        return maximum(map(d -> max(abs.(d)...), diff))
-    else
-        return maximum(map(norm, diff))
-    end
-end
-
-
-function Base.show(io::IO, printer::BuildingBlockPrinter{<:FixedGradient})
-    fp = printer.bb
-    pp(value::Number) = @sprintf("%.3g", value)
-
-    print(io, "FixedGradient(")
-    if isnothing(printer.start_time)
-        print(io, "duration=", pp(duration(fp)))
-    else
-        print(io, "t=", pp(printer.start_time), "-", pp(printer.start_time + duration(fp)))
-    end
-    for fn in (qval, gradient_strength, slew_rate)
-        print(io, ", $(nameof(fn))=", pp(fn(fp)))
-    end
-    print(io, ")")
-end
-
-"""
-    FixedInstantGradient(orientation, qval)
-
-Instantaneous MR gradient with no free variables.
-"""
-struct FixedInstantGradient <: GradientBlock
-    qvec :: SVector{3, Number}
-    rotate :: Union{Nothing, Symbol}
-    scale :: Union{Nothing, Symbol}
-end
-
-duraction(instant::FixedInstantGradient) = 0.
-qvec(instant::FixedInstantGradient) = instant.qvec
-
-fixed(f::Union{FixedGradient, FixedInstantGradient}) = f
-
-
-end
\ No newline at end of file
diff --git a/src/gradients/instant_gradients.jl b/src/gradients/instant_gradients.jl
index 2869cf6..9b1222b 100644
--- a/src/gradients/instant_gradients.jl
+++ b/src/gradients/instant_gradients.jl
@@ -80,13 +80,9 @@ end
 
 
 qvec(instant::InstantGradientBlock) = instant.qvec
-bmat_gradient(instant::InstantGradientBlock, qstart=nothing) = zeros(3, 3)
+bmat_gradient(::InstantGradientBlock, qstart=nothing) = zeros(3, 3)
 duration(instant::InstantGradientBlock) = 0.
 variables(::Type{<:InstantGradientBlock}) = [qvec, qval]
 
-function fixed(block::InstantGradientBlock)
-    return FixedInstantGradient(qvec(block), block.rotate, block.scale)
-end
-
 
 end
\ No newline at end of file
diff --git a/src/gradients/integrate_gradients.jl b/src/gradients/integrate_gradients.jl
deleted file mode 100644
index 8b28be4..0000000
--- a/src/gradients/integrate_gradients.jl
+++ /dev/null
@@ -1,70 +0,0 @@
-module IntegrateGradients
-import ...Variables: qval, bval
-import ...BuildingBlocks: BuildingBlock
-import ...Containers: Sequence
-
-
-"""
-    qval(container, 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::Sequence, indices::AbstractVector) = full_integral(builder, indices)[1]
-
-"""
-    bval(container, 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::Sequence, indices::AbstractVector) = full_integral(builder, indices)[2]
-
-function full_integral(builder::Sequence, 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
diff --git a/src/gradients/slice_selects.jl b/src/gradients/slice_selects.jl
index f68853d..c6ab2b4 100644
--- a/src/gradients/slice_selects.jl
+++ b/src/gradients/slice_selects.jl
@@ -6,21 +6,6 @@ import ...Variables: VariableType, slice_thickness, get_free_variable, effective
 import ...Variables: flip_angle, amplitude, phase, frequency, bandwidth, inverse_bandwidth, N_left, N_right
 import ..ConstantGradientBlocks: ConstantGradientBlock
 
-struct AbstractSliceSelect{T<:VariableType}
-    model :: Model
-    time_before :: T
-    flat_before :: ConstantGradientBlock
-    pulse :: RFPulseBlock
-    time_after :: T
-    flat_after :: ConstantGradientBlock
-    flat_total :: ConstantGradientBlock
-end
-
-"""
-The fixed equivalent of [`SliceSelect`](@ref).
-"""
-const FixedSliceSelect = AbstractSliceSelect{Float64}
-
 """
     SliceSelect(gradient_strength, time_before_pulse, pulse, time_after_pulse, rotate, scale)
 
@@ -31,9 +16,17 @@ Do not create it directly.
 
 The fixed equivalent is [`FixedSliceSelect`](@ref).
 """
-const SliceSelect = AbstractSliceSelect{VariableType}
+struct SliceSelect
+    model :: Model
+    time_before :: <:VariableType
+    flat_before :: ConstantGradientBlock
+    pulse :: RFPulseBlock
+    time_after :: <:VariableType
+    flat_after :: ConstantGradientBlock
+    flat_total :: ConstantGradientBlock
+end
 
-function AbstractSliceSelect{VariableType}(model::Model, gradient_strength::SVector{3, VariableType}, time_before_pulse, pulse::RFPulseBlock, time_after_pulse, rotate::Symbol, scale::Symbol)
+function SliceSelect(model::Model, gradient_strength::SVector{3, VariableType}, time_before_pulse, pulse::RFPulseBlock, time_after_pulse, rotate::Symbol, scale::Symbol)
     time_before_pulse = get_free_variable(model, time_before_pulse),
     time_after_pulse = get_free_variable(model, time_after_pulse),
     return SliceSelect(
-- 
GitLab