text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    This struct keeps track of all uses of some mutable struct allocated
    in the current function. `uses` are all instances of `getfield` on the
    struct. `defs` are all instances of `setfield!` on the struct. The terminology
    refers to the uses/defs of the ``slot bundle'' that the mutable struct represents.

    In addition we keep track of all instances of a foreigncall preserve of this mutable
    struct. Somewhat counterintuitively, we don't actually need to make sure that the
    struct itself is live (or even allocated) at a ccall site. If there are no other places
    where the struct escapes (and thus e.g. where its address is taken), it need not be
    allocated. We do however, need to make sure to preserve any elements of this struct.
"""
struct SSADefUse
    uses::Vector{Int}
    defs::Vector{Int}
    ccall_preserve_uses::Vector{Int}
end
SSADefUse() = SSADefUse(Int[], Int[], Int[])

function try_compute_fieldidx_expr(@nospecialize(typ), @nospecialize(use_expr))
    field = use_expr.args[3]
    isa(field, QuoteNode) && (field = field.value)
    isa(field, Union{Int, Symbol}) || return nothing
    return try_compute_fieldidx(typ, field)
end

function lift_defuse(cfg::CFG, ssa::SSADefUse)
    # We remove from `uses` any block where all uses are dominated
    # by a def. This prevents insertion of dead phi nodes at the top
    # of such a block if that block happens to be in a loop
    ordered = Tuple{Int, Int, Bool}[(x, block_for_inst(cfg, x), true) for x in ssa.uses]
    for x in ssa.defs
        push!(ordered, (x, block_for_inst(cfg, x), false))
    end
    ordered = sort(ordered, by=x->x[1])
    bb_defs = Int[]
    bb_uses = Int[]
    last_bb = last_def_bb = 0
    for (_, bb, is_use) in ordered
        if bb != last_bb && is_use
            push!(bb_uses, bb)
        end
        last_bb = bb
        if last_def_bb != bb && !is_use
            push!(bb_defs, bb)
            last_def_bb = bb
        end
    end
    SSADefUse(bb_uses, bb_defs, Int[])
end

function find_curblock(domtree::DomTree, allblocks::Vector{Int}, curblock::Int)
    # TODO: This can be much faster by looking at current level and only
    # searching for those blocks in a sorted order
    while !(curblock in allblocks)
        curblock = domtree.idoms_bb[curblock]
    end
    return curblock
end

function val_for_def_expr(ir::IRCode, def::Int, fidx::Int)
    if isexpr(ir[SSAValue(def)], :new)
        return ir[SSAValue(def)].args[1+fidx]
    else
        # The use is whatever the setfield was
        return ir[SSAValue(def)].args[4]
    end
end

function compute_value_for_block(ir::IRCode, domtree::DomTree, allblocks::Vector{Int}, du::SSADefUse, phinodes::IdDict{Int, SSAValue}, fidx::Int, curblock::Int)
    curblock = find_curblock(domtree, allblocks, curblock)
    def = 0
    for stmt in du.defs
        if block_for_inst(ir.cfg, stmt) == curblock
            def = max(def, stmt)
        end
    end
    def == 0 ? phinodes[curblock] : val_for_def_expr(ir, def, fidx)
end

function compute_value_for_use(ir::IRCode, domtree::DomTree, allblocks::Vector{Int}, du::SSADefUse, phinodes::IdDict{Int, SSAValue}, fidx::Int, use_idx::Int)
    # Find the first dominating def
    curblock = stmtblock = block_for_inst(ir.cfg, use_idx)
    curblock = find_curblock(domtree, allblocks, curblock)
    defblockdefs = let curblock = curblock
        Int[stmt for stmt in du.defs if block_for_inst(ir.cfg, stmt) == curblock]
    end
    def = 0
    if !isempty(defblockdefs)
        if curblock != stmtblock
            # Find the last def in this block
            def = 0
            for x in defblockdefs
                def = max(def, x)
            end
        else
            # Find the last def before our use
            def = 0
            for x in defblockdefs
                def = max(def, x >= use_idx ? 0 : x)
            end
        end
    end
    if def == 0
        if !haskey(phinodes, curblock)
            # If this happens, we need to search the predecessors for defs. Which
            # one doesn't matter - if it did, we'd have had a phinode
            return compute_value_for_block(ir, domtree, allblocks, du, phinodes, fidx, first(ir.cfg.blocks[stmtblock].preds))
        end
        # The use is the phinode
        return phinodes[curblock]
    else
        return val_for_def_expr(ir, def, fidx)
    end
end

function simple_walk(compact::IncrementalCompact, @nospecialize(defssa#=::AnySSAValue=#), pi_callback=(pi, idx)->false)
    while true
        if isa(defssa, OldSSAValue)
            if already_inserted(compact, defssa)
                rename = compact.ssa_rename[defssa.id]
                if isa(rename, AnySSAValue)
                    defssa = rename
                    continue
                end
                return rename
            end
        end
        def = compact[defssa]
        if isa(def, PiNode)
            if pi_callback(def, defssa)
                return defssa
            end
            def = def.val
            if isa(def, SSAValue)
                is_old(compact, defssa) && (def = OldSSAValue(def.id))
            else
                return def
            end
            defssa = def
        elseif isa(def, AnySSAValue)
            pi_callback(def, defssa)
            if isa(def, SSAValue)
                is_old(compact, defssa) && (def = OldSSAValue(def.id))
            end
            defssa = def
        elseif isa(def, Union{PhiNode, PhiCNode, Expr, GlobalRef})
            return defssa
        else
            return def
        end
    end
end

function simple_walk_constraint(compact::IncrementalCompact, @nospecialize(defidx), @nospecialize(typeconstraint) = types(compact)[defidx])
    callback = function (@nospecialize(pi), @nospecialize(idx))
        isa(pi, PiNode) && (typeconstraint = typeintersect(typeconstraint, widenconst(pi.typ)))
        return false
    end
    def = simple_walk(compact, defidx, callback)
    return Pair{Any, Any}(def, typeconstraint)
end

"""
    walk_to_defs(compact, val, intermediaries)

Starting at `val` walk use-def chains to get all the leaves feeding into
this val (pruning those leaves rules out by path conditions).
"""
function walk_to_defs(compact::IncrementalCompact, @nospecialize(defssa), @nospecialize(typeconstraint), visited_phinodes::Vector{Any}=Any[])
    if !isa(defssa, AnySSAValue) || !isa(compact[defssa], PhiNode)
        return Any[defssa]
    end
    # Step 2: Figure out what the struct is defined as
    def = compact[defssa]
    ## Track definitions through PiNode/PhiNode
    found_def = false
    ## Track which PhiNodes, SSAValue intermediaries
    ## we forwarded through.
    visited = IdDict{Any, Any}()
    worklist_defs = Any[]
    worklist_constraints = Any[]
    leaves = Any[]
    push!(worklist_defs, defssa)
    push!(worklist_constraints, typeconstraint)
    while !isempty(worklist_defs)
        defssa = pop!(worklist_defs)
        typeconstraint = pop!(worklist_constraints)
        visited[defssa] = typeconstraint
        def = compact[defssa]
        if isa(def, PhiNode)
            push!(visited_phinodes, defssa)
            possible_predecessors = Int[]
            for n in 1:length(def.edges)
                isassigned(def.values, n) || continue
                val = def.values[n]
                if is_old(compact, defssa) && isa(val, SSAValue)
                    val = OldSSAValue(val.id)
                end
                edge_typ = widenconst(compact_exprtype(compact, val))
                typeintersect(edge_typ, typeconstraint) === Union{} && continue
                push!(possible_predecessors, n)
            end
            for n in possible_predecessors
                pred = def.edges[n]
                val = def.values[n]
                if is_old(compact, defssa) && isa(val, SSAValue)
                    val = OldSSAValue(val.id)
                end
                if isa(val, AnySSAValue)
                    new_def, new_constraint = simple_walk_constraint(compact, val, typeconstraint)
                    if isa(new_def, AnySSAValue)
                        if !haskey(visited, new_def)
                            push!(worklist_defs, new_def)
                            push!(worklist_constraints, new_constraint)
                        elseif !(new_constraint <: visited[new_def])
                            # We have reached the same definition via a different
                            # path, with a different type constraint. We may have
                            # to redo some work here with the wider typeconstraint
                            push!(worklist_defs, new_def)
                            push!(worklist_constraints, tmerge(new_constraint, visited[new_def]))
                        end
                        continue
                    end
                    val = new_def
                end
                if def === val
                    # This shouldn't really ever happen, but
                    # patterns like this can occur in dead code,
                    # so bail out.
                    break
                else
                    push!(leaves, val)
                end
                continue
            end
        else
            push!(leaves, defssa)
        end
    end
    leaves
end

function process_immutable_preserve(new_preserves::Vector{Any}, compact::IncrementalCompact, def::Expr)
    for arg in (isexpr(def, :new) ? def.args : def.args[2:end])
        if !isbitstype(widenconst(compact_exprtype(compact, arg)))
            push!(new_preserves, arg)
        end
    end
end

function already_inserted(compact::IncrementalCompact, old::OldSSAValue)
    id = old.id
    if id < length(compact.ir.stmts)
        return id < compact.idx
    end
    id -= length(compact.ir.stmts)
    if id < length(compact.ir.new_nodes)
        error("")
    end
    id -= length(compact.ir.new_nodes)
    @assert id <= length(compact.pending_nodes)
    return !(id in compact.pending_perm)
end

function is_pending(compact::IncrementalCompact, old::OldSSAValue)
    return old.id > length(compact.ir.stmts) + length(compact.ir.new_nodes)
end

function is_getfield_captures(@nospecialize(def), compact::IncrementalCompact)
    isa(def, Expr) || return false
    length(def.args) >= 3 || return false
    is_known_call(def, getfield, compact) || return false
    which = compact_exprtype(compact, def.args[3])
    isa(which, Const) || return false
    which.val === :captures || return false
    oc = compact_exprtype(compact, def.args[2])
    return oc ⊑ Core.OpaqueClosure
end

function lift_leaves(compact::IncrementalCompact, @nospecialize(stmt),
        @nospecialize(result_t), field::Int, leaves::Vector{Any})
    # For every leaf, the lifted value
    lifted_leaves = IdDict{Any, Any}()
    maybe_undef = false
    for leaf in leaves
        leaf_key = leaf
        if isa(leaf, AnySSAValue)
            function lift_arg(ref::Core.Compiler.UseRef)
                lifted = ref[]
                if is_old(compact, leaf) && isa(lifted, SSAValue)
                    lifted = OldSSAValue(lifted.id)
                end
                if isa(lifted, GlobalRef) || isa(lifted, Expr)
                    lifted = insert_node!(compact, leaf, effect_free(NewInstruction(lifted, compact_exprtype(compact, lifted))))
                    ref[] = lifted
                    (isa(leaf, SSAValue) && (leaf.id < compact.result_idx)) && push!(compact.late_fixup, leaf.id)
                end
                lifted_leaves[leaf_key] = RefValue{Any}(lifted)
                nothing
            end
            function walk_leaf(@nospecialize(leaf))
                if isa(leaf, OldSSAValue) && already_inserted(compact, leaf)
                    leaf = compact.ssa_rename[leaf.id]
                    if isa(leaf, AnySSAValue)
                        leaf = simple_walk(compact, leaf)
                    end
                    if isa(leaf, AnySSAValue)
                        def = compact[leaf]
                    else
                        def = leaf
                    end
                elseif isa(leaf, AnySSAValue)
                    def = compact[leaf]
                else
                    def = leaf
                end
                return Pair{Any, Any}(def, leaf)
            end
            (def, leaf) = walk_leaf(leaf)
            if is_tuple_call(compact, def) && 1 <= field < length(def.args)
                lift_arg(UseRef(def, 1 + field))
                continue
            elseif isexpr(def, :new)
                typ = widenconst(types(compact)[leaf])
                if isa(typ, UnionAll)
                    typ = unwrap_unionall(typ)
                end
                (isa(typ, DataType) && !isabstracttype(typ)) || return nothing
                @assert !ismutabletype(typ)
                if length(def.args) < 1 + field
                    if field > fieldcount(typ)
                        return nothing
                    end
                    ftyp = fieldtype(typ, field)
                    if !isbitstype(ftyp)
                        # On this branch, this will be a guaranteed UndefRefError.
                        # We use the regular undef mechanic to lift this to a boolean slot
                        maybe_undef = true
                        lifted_leaves[leaf_key] = nothing
                        continue
                    end
                    return nothing
                    # Expand the Expr(:new) to include it's element Expr(:new) nodes up until the one we want
                    compact[leaf] = nothing
                    for i = (length(def.args) + 1):(1+field)
                        ftyp = fieldtype(typ, i - 1)
                        isbitstype(ftyp) || return nothing
                        ninst = effect_free(NewInstruction(Expr(:new, ftyp), result_t))
                        push!(def.args, insert_node!(compact, leaf, ninst))
                    end
                    compact[leaf] = def
                end
                lifted = def.args[1+field]
                if is_old(compact, leaf) && isa(lifted, SSAValue)
                    lifted = OldSSAValue(lifted.id)
                end
                if isa(lifted, GlobalRef) || isa(lifted, Expr)
                    lifted = insert_node!(compact, leaf, effect_free(NewInstruction(lifted, compact_exprtype(compact, lifted))))
                    def.args[1+field] = lifted
                    (isa(leaf, SSAValue) && (leaf.id < compact.result_idx)) && push!(compact.late_fixup, leaf.id)
                end
                lifted_leaves[leaf_key] = RefValue{Any}(lifted)
                continue
            elseif is_getfield_captures(def, compact)
                # Walk to new_opaque_closure
                ocleaf = def.args[2]
                if isa(ocleaf, AnySSAValue)
                    ocleaf = simple_walk(compact, ocleaf)
                end
                ocdef, _ = walk_leaf(ocleaf)
                if isexpr(ocdef, :new_opaque_closure) && isa(field, Int) && 1 <= field <= length(ocdef.args)-5
                    lift_arg(UseRef(ocdef, 5 + field))
                    continue
                end
                return nothing
            else
                typ = compact_exprtype(compact, leaf)
                if !isa(typ, Const)
                    # TODO: (disabled since #27126)
                    # If the leaf is an old ssa value, insert a getfield here
                    # We will revisit this getfield later when compaction gets
                    # to the appropriate point.
                    # N.B.: This can be a bit dangerous because it can lead to
                    # infinite loops if we accidentally insert a node just ahead
                    # of where we are
                    return nothing
                end
                leaf = typ.val
                # Fall through to below
            end
        elseif isa(leaf, QuoteNode)
            leaf = leaf.value
        elseif isa(leaf, GlobalRef)
            mod, name = leaf.mod, leaf.name
            if isdefined(mod, name) && isconst(mod, name)
                leaf = getfield(mod, name)
            else
                return nothing
            end
        elseif isa(leaf, Union{Argument, Expr})
            return nothing
        end
        !ismutable(leaf) || return nothing
        isdefined(leaf, field) || return nothing
        val = getfield(leaf, field)
        is_inlineable_constant(val) || return nothing
        lifted_leaves[leaf_key] = RefValue{Any}(quoted(val))
    end
    lifted_leaves, maybe_undef
end

make_MaybeUndef(@nospecialize(typ)) = isa(typ, MaybeUndef) ? typ : MaybeUndef(typ)

function lift_comparison!(compact::IncrementalCompact, idx::Int,
        @nospecialize(c1), @nospecialize(c2), stmt::Expr,
        lifting_cache::IdDict{Pair{AnySSAValue, Any}, AnySSAValue})
    if isa(c1, Const)
        cmp = c1
        typeconstraint = widenconst(c2)
        val = stmt.args[3]
    else
        cmp = c2
        typeconstraint = widenconst(c1)
        val = stmt.args[2]
    end

    is_type_only = isdefined(typeof(cmp), :instance)

    if isa(val, Union{OldSSAValue, SSAValue})
        val, typeconstraint = simple_walk_constraint(compact, val, typeconstraint)
    end

    visited_phinodes = Any[]
    leaves = walk_to_defs(compact, val, typeconstraint, visited_phinodes)

    # Let's check if we evaluate the comparison for each one of the leaves
    lifted_leaves = IdDict{Any, Any}()
    for leaf in leaves
        r = egal_tfunc(compact_exprtype(compact, leaf), cmp)
        if isa(r, Const)
            lifted_leaves[leaf] = RefValue{Any}(r.val)
        else
            # TODO: In some cases it might be profitable to hoist the ===
            # here.
            return
        end
    end

    lifted_val = perform_lifting!(compact, visited_phinodes, cmp, lifting_cache, Bool, lifted_leaves, val)
    @assert lifted_val !== nothing

    #global assertion_counter
    #assertion_counter::Int += 1
    #insert_node_here!(compact, Expr(:assert_egal, Symbol(string("assert_egal_", assertion_counter)), SSAValue(idx), lifted_val), nothing, 0, true)
    #return
    compact[idx] = lifted_val.x
end

struct LiftedPhi
    ssa::AnySSAValue
    node::Any
    need_argupdate::Bool
end

function is_old(compact, @nospecialize(old_node_ssa))
    isa(old_node_ssa, OldSSAValue) &&
        !is_pending(compact, old_node_ssa) &&
        !already_inserted(compact, old_node_ssa)
end

function perform_lifting!(compact::IncrementalCompact,
        visited_phinodes::Vector{Any}, @nospecialize(cache_key),
        lifting_cache::IdDict{Pair{AnySSAValue, Any}, AnySSAValue},
        @nospecialize(result_t), lifted_leaves::IdDict{Any, Any}, @nospecialize(stmt_val))
    reverse_mapping = IdDict{Any, Any}(ssa => id for (id, ssa) in enumerate(visited_phinodes))

    # Insert PhiNodes
    lifted_phis = LiftedPhi[]
    for item in visited_phinodes
        if (item, cache_key) in keys(lifting_cache)
            ssa = lifting_cache[Pair{AnySSAValue, Any}(item, cache_key)]
            push!(lifted_phis, LiftedPhi(ssa, compact[ssa], false))
            continue
        end
        n = PhiNode()
        ssa = insert_node!(compact, item, effect_free(NewInstruction(n, result_t)))
        lifting_cache[Pair{AnySSAValue, Any}(item, cache_key)] = ssa
        push!(lifted_phis, LiftedPhi(ssa, n, true))
    end

    # Fix up arguments
    for (old_node_ssa, lf) in zip(visited_phinodes, lifted_phis)
        old_node = compact[old_node_ssa]
        new_node = lf.node
        lf.need_argupdate || continue
        for i = 1:length(old_node.edges)
            edge = old_node.edges[i]
            isassigned(old_node.values, i) || continue
            val = old_node.values[i]
            orig_val = val
            if is_old(compact, old_node_ssa) && isa(val, SSAValue)
                val = OldSSAValue(val.id)
            end
            if isa(val, Union{NewSSAValue, SSAValue, OldSSAValue})
                val = simple_walk(compact, val)
            end
            if val in keys(lifted_leaves)
                push!(new_node.edges, edge)
                lifted_val = lifted_leaves[val]
                if lifted_val === nothing
                    resize!(new_node.values, length(new_node.values)+1)
                    continue
                end
                lifted_val = lifted_val.x
                if isa(lifted_val, Union{NewSSAValue, SSAValue, OldSSAValue})
                    lifted_val = simple_walk(compact, lifted_val, (pi, idx)->true)
                end
                push!(new_node.values, lifted_val)
            elseif isa(val, Union{NewSSAValue, SSAValue, OldSSAValue}) && val in keys(reverse_mapping)
                push!(new_node.edges, edge)
                push!(new_node.values, lifted_phis[reverse_mapping[val]].ssa)
            else
                # Probably ignored by path condition, skip this
            end
        end
    end

    for lf in lifted_phis
        count_added_node!(compact, lf.node)
    end

    # Fixup the stmt itself
    if isa(stmt_val, Union{SSAValue, OldSSAValue})
        stmt_val = simple_walk(compact, stmt_val)
    end

    if stmt_val in keys(lifted_leaves)
        stmt_val = lifted_leaves[stmt_val]
    elseif isa(stmt_val, Union{NewSSAValue, SSAValue, OldSSAValue}) && stmt_val in keys(reverse_mapping)
        stmt_val = RefValue{Any}(lifted_phis[reverse_mapping[stmt_val]].ssa)
    end

    return stmt_val
end

assertion_counter = 0
function getfield_elim_pass!(ir::IRCode)
    compact = IncrementalCompact(ir)
    insertions = Vector{Any}()
    defuses = IdDict{Int, Tuple{IdSet{Int}, SSADefUse}}()
    lifting_cache = IdDict{Pair{AnySSAValue, Any}, AnySSAValue}()
    revisit_worklist = Int[]
    #ndone, nmax = 0, 200
    for ((_, idx), stmt) in compact
        isa(stmt, Expr) || continue
        #ndone >= nmax && continue
        #ndone += 1
        result_t = compact_exprtype(compact, SSAValue(idx))
        is_getfield = is_setfield = false
        field_ordering = :unspecified
        is_ccall = false
        # Step 1: Check whether the statement we're looking at is a getfield/setfield!
        if is_known_call(stmt, setfield!, compact)
            is_setfield = true
            4 <= length(stmt.args) <= 5 || continue
            if length(stmt.args) == 5
                field_ordering = compact_exprtype(compact, stmt.args[5])
            end
        elseif is_known_call(stmt, getfield, compact)
            is_getfield = true
            3 <= length(stmt.args) <= 5 || continue
            if length(stmt.args) == 5
                field_ordering = compact_exprtype(compact, stmt.args[5])
            elseif length(stmt.args) == 4
                field_ordering = compact_exprtype(compact, stmt.args[4])
                widenconst(field_ordering) === Bool && (field_ordering = :unspecified)
            end
        elseif is_known_call(stmt, isa, compact)
            # TODO
            continue
        elseif is_known_call(stmt, typeassert, compact) && length(stmt.args) == 3
            # Canonicalize
            #   X = typeassert(Y, T)::S
            # into
            #   typeassert(Y, T)
            #   X = PiNode(Y, S)
            # N.B.: Inference may have a more precise type for `S`, than
            #       just T, but from here on out, there's no problem with
            #       using just using that.
            # so subsequent analysis only has to deal with the latter
            # form. TODO: This isn't the best place to put this.
            # Also, we should probably have a version of typeassert
            # that's defined not to return its value to make life easier
            # for the backend.
            pi = insert_node_here!(compact,
                NewInstruction(
                    PiNode(stmt.args[2], compact.result[idx][:type]),
                    compact.result[idx][:type],
                    compact.result[idx][:line]), true)
            compact.ssa_rename[compact.idx-1] = pi
            continue
        elseif is_known_call(stmt, (===), compact) && length(stmt.args) == 3
            c1 = compact_exprtype(compact, stmt.args[2])
            c2 = compact_exprtype(compact, stmt.args[3])
            if !(isa(c1, Const) || isa(c2, Const))
                continue
            end
            (isa(c1, Const) && isa(c2, Const)) && continue
            lift_comparison!(compact, idx, c1, c2, stmt, lifting_cache)
            continue
        elseif isexpr(stmt, :foreigncall)
            nccallargs = length(stmt.args[3]::SimpleVector)
            new_preserves = Any[]
            old_preserves = stmt.args[(6+nccallargs):end]
            for (pidx, preserved_arg) in enumerate(old_preserves)
                isa(preserved_arg, SSAValue) || continue
                let intermediaries = IdSet()
                    callback = function(@nospecialize(pi), ssa::AnySSAValue)
                        push!(intermediaries, ssa.id)
                        return false
                    end
                    def = simple_walk(compact, preserved_arg, callback)
                    isa(def, SSAValue) || continue
                    defidx = def.id
                    def = compact[defidx]
                    if is_tuple_call(compact, def)
                        process_immutable_preserve(new_preserves, compact, def)
                        old_preserves[pidx] = nothing
                        continue
                    elseif isexpr(def, :new)
                        typ = widenconst(compact_exprtype(compact, SSAValue(defidx)))
                        if isa(typ, UnionAll)
                            typ = unwrap_unionall(typ)
                        end
                        if typ isa DataType && !ismutabletype(typ)
                            process_immutable_preserve(new_preserves, compact, def)
                            old_preserves[pidx] = nothing
                            continue
                        end
                    else
                        continue
                    end
                    mid, defuse = get!(defuses, defidx, (IdSet{Int}(), SSADefUse()))
                    push!(defuse.ccall_preserve_uses, idx)
                    union!(mid, intermediaries)
                end
                continue
            end
            if !isempty(new_preserves)
                old_preserves = filter(ssa->ssa !== nothing, old_preserves)
                new_expr = Expr(:foreigncall, stmt.args[1:(6+nccallargs-1)]...,
                    old_preserves..., new_preserves...)
                compact[idx] = new_expr
            end
            continue
        else
            continue
        end
        ## Normalize the field argument to getfield/setfield
        field = stmt.args[3]
        isa(field, QuoteNode) && (field = field.value)
        isa(field, Union{Int, Symbol}) || continue

        struct_typ = unwrap_unionall(widenconst(compact_exprtype(compact, stmt.args[2])))
        if isa(struct_typ, Union) && struct_typ <: Tuple
            struct_typ = unswitchtupleunion(struct_typ)
        end
        isa(struct_typ, DataType) || continue

        struct_typ.name.atomicfields == C_NULL || continue # TODO: handle more
        if !(field_ordering === :unspecified || (field_ordering isa Const && field_ordering.val === :not_atomic))
            continue
        end

        def, typeconstraint = stmt.args[2], struct_typ

        if ismutabletype(struct_typ)
            isa(def, SSAValue) || continue
            let intermediaries = IdSet()
                callback = function(@nospecialize(pi), ssa::AnySSAValue)
                    push!(intermediaries, ssa.id)
                    return false
                end
                def = simple_walk(compact, def, callback)
                # Mutable stuff here
                isa(def, SSAValue) || continue
                mid, defuse = get!(defuses, def.id, (IdSet{Int}(), SSADefUse()))
                if is_setfield
                    push!(defuse.defs, idx)
                else
                    push!(defuse.uses, idx)
                end
                union!(mid, intermediaries)
            end
            continue
        elseif is_setfield
            continue
        end

        if isa(def, Union{OldSSAValue, SSAValue})
            def, typeconstraint = simple_walk_constraint(compact, def, typeconstraint)
        end

        visited_phinodes = Any[]
        leaves = walk_to_defs(compact, def, typeconstraint, visited_phinodes)

        isempty(leaves) && continue

        field = try_compute_fieldidx(struct_typ, field)
        field === nothing && continue

        r = lift_leaves(compact, stmt, result_t, field, leaves)
        r === nothing && continue
        lifted_leaves, any_undef = r

        if any_undef
            result_t = make_MaybeUndef(result_t)
        end

#        @Base.show result_t
#        @Base.show stmt
#        for (k,v) in lifted_leaves
#            @Base.show (k, v)
#            if isa(k, AnySSAValue)
#                @Base.show compact[k]
#            end
#            if isa(v, RefValue) && isa(v.x, AnySSAValue)
#                @Base.show compact[v.x]
#            end
#        end
        val = perform_lifting!(compact, visited_phinodes, field, lifting_cache, result_t, lifted_leaves, stmt.args[2])

        # Insert the undef check if necessary
        if any_undef
            if val === nothing
                insert_node!(compact, SSAValue(idx),
                    non_effect_free(NewInstruction(Expr(:throw_undef_if_not, Symbol("##getfield##"), false), Nothing)))
            else
                # val must be defined
            end
        else
            @assert val !== nothing
        end

        global assertion_counter
        assertion_counter::Int += 1
        #insert_node_here!(compact, Expr(:assert_egal, Symbol(string("assert_egal_", assertion_counter)), SSAValue(idx), val), nothing, 0, true)
        #continue
        compact[idx] = val === nothing ? nothing : val.x
    end

    non_dce_finish!(compact)
    # Copy the use count, `simple_dce!` may modify it and for our predicate
    # below we need it consistent with the state of the IR here (after tracking
    # phi node arguments, but before dce).
    used_ssas = copy(compact.used_ssas)
    simple_dce!(compact)
    ir = complete(compact)

    # Compute domtree, needed below, now that we have finished compacting the
    # IR. This needs to be after we iterate through the IR with
    # `IncrementalCompact` because removing dead blocks can invalidate the
    # domtree.
    @timeit "domtree 2" domtree = construct_domtree(ir.cfg.blocks)

    # Now go through any mutable structs and see which ones we can eliminate
    for (idx, (intermediaries, defuse)) in defuses
        intermediaries = collect(intermediaries)
        # Check if there are any uses we did not account for. If so, the variable
        # escapes and we cannot eliminate the allocation. This works, because we're guaranteed
        # not to include any intermediaries that have dead uses. As a result, missing uses will only ever
        # show up in the nuses_total count.
        nleaves = length(defuse.uses) + length(defuse.defs) + length(defuse.ccall_preserve_uses)
        nuses = 0
        for idx in intermediaries
            nuses += used_ssas[idx]
        end
        nuses_total = used_ssas[idx] + nuses - length(intermediaries)
        nleaves == nuses_total || continue
        # Find the type for this allocation
        defexpr = ir[SSAValue(idx)]
        isexpr(defexpr, :new) || continue
        typ = ir.stmts[idx][:type]
        if isa(typ, UnionAll)
            typ = unwrap_unionall(typ)
        end
        # Could still end up here if we tried to setfield! and immutable, which would
        # error at runtime, but is not illegal to have in the IR.
        ismutabletype(typ) || continue
        # Partition defuses by field
        fielddefuse = SSADefUse[SSADefUse() for _ = 1:fieldcount(typ)]
        ok = true
        for use in defuse.uses
            stmt = ir[SSAValue(use)]
            # We may have discovered above that this use is dead
            # after the getfield elim of immutables. In that case,
            # it would have been deleted. That's fine, just ignore
            # the use in that case.
            stmt === nothing && continue
            field = try_compute_fieldidx_expr(typ, stmt)
            field === nothing && (ok = false; break)
            push!(fielddefuse[field].uses, use)
        end
        ok || continue
        for use in defuse.defs
            field = try_compute_fieldidx_expr(typ, ir[SSAValue(use)])
            field === nothing && (ok = false; break)
            push!(fielddefuse[field].defs, use)
        end
        ok || continue
        # Check that the defexpr has defined values for all the fields
        # we're accessing. In the future, we may want to relax this,
        # but we should come up with semantics for well defined semantics
        # for uninitialized fields first.
        for (fidx, du) in pairs(fielddefuse)
            isempty(du.uses) && continue
            if fidx + 1 > length(defexpr.args)
                ok = false
                break
            end
        end
        ok || continue
        preserve_uses = IdDict{Int, Vector{Any}}((idx=>Any[] for idx in IdSet{Int}(defuse.ccall_preserve_uses)))
        # Everything accounted for. Go field by field and perform idf
        for (fidx, du) in pairs(fielddefuse)
            ftyp = fieldtype(typ, fidx)
            if !isempty(du.uses)
                push!(du.defs, idx)
                ldu = compute_live_ins(ir.cfg, du)
                phiblocks = Int[]
                if !isempty(ldu.live_in_bbs)
                    phiblocks = idf(ir.cfg, ldu, domtree)
                end
                phinodes = IdDict{Int, SSAValue}()
                for b in phiblocks
                    n = PhiNode()
                    phinodes[b] = insert_node!(ir, first(ir.cfg.blocks[b].stmts),
                        NewInstruction(n, ftyp))
                end
                # Now go through all uses and rewrite them
                allblocks = sort(vcat(phiblocks, ldu.def_bbs))
                for stmt in du.uses
                    ir[SSAValue(stmt)] = compute_value_for_use(ir, domtree, allblocks, du, phinodes, fidx, stmt)
                end
                if !isbitstype(ftyp)
                    for (use, list) in preserve_uses
                        push!(list, compute_value_for_use(ir, domtree, allblocks, du, phinodes, fidx, use))
                    end
                end
                for b in phiblocks
                    for p in ir.cfg.blocks[b].preds
                        n = ir[phinodes[b]]
                        push!(n.edges, p)
                        push!(n.values, compute_value_for_block(ir, domtree,
                            allblocks, du, phinodes, fidx, p))
                    end
                end
            end
            for stmt in du.defs
                stmt == idx && continue
                ir[SSAValue(stmt)] = nothing
            end
            continue
        end
        isempty(defuse.ccall_preserve_uses) && continue
        push!(intermediaries, idx)
        # Insert the new preserves
        for (use, new_preserves) in preserve_uses
            useexpr = ir[SSAValue(use)]
            nccallargs = length(useexpr.args[3]::SimpleVector)
            old_preserves = let intermediaries = intermediaries
                filter(ssa->!isa(ssa, SSAValue) || !(ssa.id in intermediaries), useexpr.args[(6+nccallargs):end])
            end
            new_expr = Expr(:foreigncall, useexpr.args[1:(6+nccallargs-1)]...,
                old_preserves..., new_preserves...)
            ir[SSAValue(use)] = new_expr
        end
    end
    ir
end

function adce_erase!(phi_uses::Vector{Int}, extra_worklist::Vector{Int}, compact::IncrementalCompact, idx::Int)
    # return whether this made a change
    if isa(compact.result[idx][:inst], PhiNode)
        return maybe_erase_unused!(extra_worklist, compact, idx, val -> phi_uses[val.id] -= 1)
    else
        return maybe_erase_unused!(extra_worklist, compact, idx)
    end
end

function count_uses(@nospecialize(stmt), uses::Vector{Int})
    for ur in userefs(stmt)
        use = ur[]
        if isa(use, SSAValue)
            uses[use.id] += 1
        end
    end
end

function mark_phi_cycles(compact::IncrementalCompact, safe_phis::BitSet, phi::Int)
    worklist = Int[]
    push!(worklist, phi)
    while !isempty(worklist)
        phi = pop!(worklist)
        push!(safe_phis, phi)
        for ur in userefs(compact.result[phi][:inst])
            val = ur[]
            isa(val, SSAValue) || continue
            isa(compact[val], PhiNode) || continue
            (val.id in safe_phis) && continue
            push!(worklist, val.id)
        end
    end
end

function adce_pass!(ir::IRCode)
    phi_uses = fill(0, length(ir.stmts) + length(ir.new_nodes))
    all_phis = Int[]
    compact = IncrementalCompact(ir)
    for ((_, idx), stmt) in compact
        if isa(stmt, PhiNode)
            push!(all_phis, idx)
        end
    end
    non_dce_finish!(compact)
    for phi in all_phis
        count_uses(compact.result[phi][:inst]::PhiNode, phi_uses)
    end
    # Perform simple DCE for unused values
    extra_worklist = Int[]
    for (idx, nused) in Iterators.enumerate(compact.used_ssas)
        idx >= compact.result_idx && break
        nused == 0 || continue
        adce_erase!(phi_uses, extra_worklist, compact, idx)
    end
    while !isempty(extra_worklist)
        adce_erase!(phi_uses, extra_worklist, compact, pop!(extra_worklist))
    end
    # Go back and erase any phi cycles
    changed = true
    while changed
        changed = false
        safe_phis = BitSet()
        for phi in all_phis
            # Save any phi cycles that have non-phi uses
            if compact.used_ssas[phi] - phi_uses[phi] != 0
                mark_phi_cycles(compact, safe_phis, phi)
            end
        end
        for phi in all_phis
            if !(phi in safe_phis)
                push!(extra_worklist, phi)
            end
        end
        while !isempty(extra_worklist)
            if adce_erase!(phi_uses, extra_worklist, compact, pop!(extra_worklist))
                changed = true
            end
        end
    end
    return complete(compact)
end

function type_lift_pass!(ir::IRCode)
    type_ctx_uses = Vector{Vector{Int}}[]
    has_non_type_ctx_uses = IdSet{Int}()
    lifted_undef = IdDict{Int, Any}()
    insts = ir.stmts
    for idx in 1:length(insts)
        stmt = insts[idx][:inst]
        stmt isa Expr || continue
        if (stmt.head === :isdefined || stmt.head === :undefcheck)
            # after optimization, undef can only show up by being introduced in
            # a phi node (or an UpsilonNode() argument to a PhiC node), so lift
            # all these nodes that have maybe undef values
            val = stmt.args[(stmt.head === :isdefined) ? 1 : 2]
            if stmt.head === :isdefined && (val isa Slot || val isa GlobalRef ||
                    isexpr(val, :static_parameter) || val isa Argument || val isa Symbol)
                # this is a legal node, so assume it was not introduced by
                # slot2ssa (at worst, we might leave in a runtime check that
                # shouldn't have been there)
                continue
            end
            # otherwise, we definitely have a corrupt node from slot2ssa, and
            # must fix or delete that now
            processed = IdDict{Int, Union{SSAValue, Bool}}()
            def = val
            while true
                # peek through PiNodes
                isa(val, SSAValue) || break
                def = insts[val.id][:inst]
                isa(def, PiNode) || break
                val = def.val
            end
            if !isa(val, SSAValue) || (!isa(def, PhiNode) && !isa(def, PhiCNode))
                # in most cases, reaching this statement implies we had a value
                if stmt.head === :undefcheck
                    insts[idx][:inst] = nothing
                else
                    insts[idx][:inst] = true
                end
                continue
            end
            stmt_id = val.id
            worklist = Tuple{Int, Int, SSAValue, Int}[(stmt_id, 0, SSAValue(0), 0)]
            if !haskey(lifted_undef, stmt_id)
                first = true
                while !isempty(worklist)
                    item, w_up_id, which, use = pop!(worklist)
                    def = insts[item][:inst]
                    if isa(def, PhiNode)
                        edges = copy(def.edges)
                        values = Vector{Any}(undef, length(edges))
                        new_phi = if length(values) == 0
                            false
                        else
                            insert_node!(ir, item, NewInstruction(PhiNode(edges, values), Bool))
                        end
                    else
                        values = Vector{Any}(undef, length(def.values))
                        new_phi = if length(values) == 0
                            false
                        else
                            insert_node!(ir, item, NewInstruction(PhiCNode(values), Bool))
                        end
                    end
                    processed[item] = new_phi
                    if first
                        lifted_undef[stmt_id] = new_phi
                        first = false
                    end
                    local id::Int = 0
                    for i = 1:length(values)
                        if !isassigned(def.values, i)
                            val = false
                        elseif !isa(def.values[i], SSAValue)
                            val = true
                        else
                            up_id = id = def.values[i].id
                            @label restart
                            if !isa(ir.stmts[id][:type], MaybeUndef)
                                val = true
                            else
                                node = insts[id][:inst]
                                if isa(node, UpsilonNode)
                                    if !isdefined(node, :val)
                                        val = false
                                    elseif !isa(node.val, SSAValue)
                                        val = true
                                    else
                                        id = node.val.id
                                        @goto restart
                                    end
                                else
                                    while isa(node, PiNode)
                                        id = node.val.id
                                        node = insts[id][:inst]
                                    end
                                    if isa(node, Union{PhiNode, PhiCNode})
                                        if haskey(processed, id)
                                            val = processed[id]
                                        else
                                            push!(worklist, (id, up_id, new_phi, i))
                                            continue
                                        end
                                    else
                                        val = true
                                    end
                                end
                            end
                        end
                        if isa(def, PhiNode)
                            values[i] = val
                        else
                            values[i] = insert_node!(ir, up_id, NewInstruction(UpsilonNode(val), Bool))
                        end
                    end
                    if which !== SSAValue(0)
                        phi = ir[which]
                        if isa(phi, PhiNode)
                            phi.values[use] = new_phi
                        else
                            phi = phi::PhiCNode
                            phi.values[use] = insert_node!(ir, w_up_id, NewInstruction(UpsilonNode(new_phi), Bool))
                        end
                    end
                end
            end
            inst = lifted_undef[stmt_id]
            if stmt.head === :undefcheck
                inst = Expr(:throw_undef_if_not, stmt.args[1], inst)
            end
            insts[idx][:inst] = inst
        end
    end
    ir
end

function cfg_simplify!(ir::IRCode)
    bbs = ir.cfg.blocks
    merge_into = zeros(Int, length(bbs))
    merged_succ = zeros(Int, length(bbs))
    function follow_merge_into(idx::Int)
        while merge_into[idx] != 0
            idx = merge_into[idx]
        end
        return idx
    end
    function follow_merged_succ(idx::Int)
        while merged_succ[idx] != 0
            idx = merged_succ[idx]
        end
        return idx
    end

    # Walk the CFG from the entry block and aggressively combine blocks
    for (idx, bb) in enumerate(bbs)
        if length(bb.succs) == 1
            succ = bb.succs[1]
            if length(bbs[succ].preds) == 1
                # Prevent cycles by making sure we don't end up back at `idx`
                # by following what is to be merged into `succ`
                if follow_merged_succ(succ) != idx
                    merge_into[succ] = idx
                    merged_succ[idx] = succ
                end
            end
        end
    end

    # Assign new BB numbers
    max_bb_num = 1
    bb_rename_succ = zeros(Int, length(bbs))
    for i = 1:length(bbs)
        # Drop blocks that will be merged away
        if merge_into[i] != 0
            bb_rename_succ[i] = -1
        end
        # Drop blocks with no predecessors
        if i != 1 && length(ir.cfg.blocks[i].preds) == 0
            bb_rename_succ[i] = -1
        end

        bb_rename_succ[i] != 0 && continue

        curr = i
        while true
            bb_rename_succ[curr] = max_bb_num
            max_bb_num += 1
            # Now walk the chain of blocks we merged.
            # If we end in something that may fall through,
            # we have to schedule that block next
            curr = follow_merged_succ(curr)
            terminator = ir.stmts[ir.cfg.blocks[curr].stmts[end]][:inst]
            if isa(terminator, GotoNode) || isa(terminator, ReturnNode)
                break
            end
            curr += 1
        end
    end

    # Figure out how predecessors should be renamed
    bb_rename_pred = zeros(Int, length(bbs))
    for i = 1:length(bbs)
        if merged_succ[i] != 0
            # Block `i` should no longer be a predecessor (before renaming)
            # because it is being merged with its sole successor
            bb_rename_pred[i] = -1
            continue
        end
        bbnum = follow_merge_into(i)
        bb_rename_pred[i] = bb_rename_succ[bbnum]
    end

    # Compute map from new to old blocks
    result_bbs = Int[findfirst(j->i==j, bb_rename_succ) for i = 1:max_bb_num-1]

    # Compute new block lengths
    result_bbs_lengths = zeros(Int, max_bb_num-1)
    for (idx, orig_bb) in enumerate(result_bbs)
        ms = orig_bb
        while ms != 0
            result_bbs_lengths[idx] += length(bbs[ms].stmts)
            ms = merged_succ[ms]
        end
    end

    # Compute statement indices the new blocks start at
    bb_starts = Vector{Int}(undef, 1+length(result_bbs_lengths))
    bb_starts[1] = 1
    for i = 1:length(result_bbs_lengths)
        bb_starts[i+1] = bb_starts[i] + result_bbs_lengths[i]
    end

    cresult_bbs = let result_bbs = result_bbs,
                      merged_succ = merged_succ,
                      merge_into = merge_into,
                      bbs = bbs,
                      bb_rename_succ = bb_rename_succ

        # Compute (renamed) successors and predecessors given (renamed) block
        function compute_succs(i)
            orig_bb = follow_merged_succ(result_bbs[i])
            return map(i -> bb_rename_succ[i], bbs[orig_bb].succs)
        end
        function compute_preds(i)
            orig_bb = result_bbs[i]
            preds = bbs[orig_bb].preds
            return map(pred -> bb_rename_pred[pred], preds)
        end

        BasicBlock[
            BasicBlock(StmtRange(bb_starts[i],
                                 i+1 > length(bb_starts) ?
                                    length(compact.result) : bb_starts[i+1]-1),
                       compute_preds(i),
                       compute_succs(i))
            for i = 1:length(result_bbs)]
    end

    compact = IncrementalCompact(ir, true)
    # Run instruction compaction to produce the result,
    # but we're messing with the CFG
    # so we don't want compaction to do so independently
    compact.fold_constant_branches = false
    compact.bb_rename_succ = bb_rename_succ
    compact.bb_rename_pred = bb_rename_pred
    compact.result_bbs = cresult_bbs
    result_idx = 1
    for (idx, orig_bb) in enumerate(result_bbs)
        ms = orig_bb
        while ms != 0
            for i in bbs[ms].stmts
                node = ir.stmts[i]
                compact.result[compact.result_idx] = node
                if isa(node[:inst], GotoNode) && merged_succ[ms] != 0
                    # If we merged a basic block, we need remove the trailing GotoNode (if any)
                    compact.result[compact.result_idx][:inst] = nothing
                else
                    process_node!(compact, compact.result_idx, node, i, i, ms, true)
                end
                # We always increase the result index to ensure a predicatable
                # placement of the resulting nodes.
                compact.result_idx += 1
            end
            ms = merged_succ[ms]
        end
    end
    compact.active_result_bb = length(bb_starts)
    return finish(compact)
end


text:
####
#### Helper functions
####

# Element matrix for refocusing pulse with angle α (in degrees); acts on the magnetization state vector (MPSV)
@inline element_flipmat(α::T) where {T} = SA{Complex{T}}[
        cosd(α/2)^2    sind(α/2)^2 -im*sind(α);
        sind(α/2)^2    cosd(α/2)^2  im*sind(α);
    -im*sind(α)/2   im*sind(α)/2       cosd(α)]

####
####
####

####
# EPG parameterization interface
# For each parameterization, define:
#   - struct MyParameterization{T,ETL} <: FieldVector{N,T}
#   - restructure(::MyParameterization, ::NTuple{N,T})
#   - parameter getters

struct EPGOptions{T,ETL} <: FieldVector{5,T}
    α::T
    TE::T
    T2::T
    T1::T
    β::T
end
Base.NamedTuple(θ::EPGOptions) = NamedTuple{(:α, :TE, :T2, :T1, :β)}(Tuple(θ))
@inline EPGOptions(xs::NamedTuple{(:α, :TE, :T2, :T1, :β)}, ::Val{ETL}, ::Type{T} = floattype(xs)) where {T,ETL} = EPGOptions{T,ETL}(Tuple(xs))
@inline restructure(::EPGOptions{<:Any, ETL}, xs::NTuple{5,T}) where {T,ETL} = EPGOptions{T,ETL}(xs)

@inline Base.eltype(::EPGOptions{T}) where {T} = T
@inline echotrainlength(::EPGOptions{T,ETL}) where {T,ETL} = ETL
@inline B1correction(θ::EPGOptions{T}) where {T} = T(θ.α / 180) # Multiplicative FA correction: A = α/180
@inline flipangle(θ::EPGOptions{T}, i::Int) where {T} = ifelse(i == 0, T(90), ifelse(i == 1, T(180), θ.β)) # Pulse sequence: 90, 180, β, β, ...
@inline echotime(θ::EPGOptions{T}) where {T} = θ.TE
@inline T2time(θ::EPGOptions{T}) where {T} = θ.T2
@inline T1time(θ::EPGOptions{T}) where {T} = θ.T1

struct EPGIncreasingFlipAnglesOptions{T,ETL} <: FieldVector{6,T}
    α::T
    α1::T
    α2::T
    TE::T
    T2::T
    T1::T
end
Base.NamedTuple(θ::EPGIncreasingFlipAnglesOptions) = NamedTuple{(:α, :α1, :α2, :TE, :T2, :T1)}(Tuple(θ))
@inline EPGIncreasingFlipAnglesOptions(xs::NamedTuple{(:α, :α1, :α2, :TE, :T2, :T1)}, ::Val{ETL}, ::Type{T} = floattype(xs)) where {T,ETL} = EPGIncreasingFlipAnglesOptions{T,ETL}(Tuple(xs))
@inline restructure(::EPGIncreasingFlipAnglesOptions{<:Any, ETL}, xs::NTuple{6,T}) where {T,ETL} = EPGIncreasingFlipAnglesOptions{T,ETL}(xs)

@inline Base.eltype(::EPGIncreasingFlipAnglesOptions{T}) where {T} = T
@inline echotrainlength(::EPGIncreasingFlipAnglesOptions{T,ETL}) where {T,ETL} = ETL
@inline B1correction(θ::EPGIncreasingFlipAnglesOptions{T}) where {T} = T(θ.α / 180) # Multiplicative FA correction: A = α/180
@inline flipangle(θ::EPGIncreasingFlipAnglesOptions{T}, i::Int) where {T} = ifelse(i == 0, T(90), ifelse(i == 1, θ.α1, ifelse(i == 2, θ.α2, T(180)))) # Pulse sequence: 90, α1, α2, 180, 180, ...
@inline echotime(θ::EPGIncreasingFlipAnglesOptions{T}) where {T} = θ.TE
@inline T2time(θ::EPGIncreasingFlipAnglesOptions{T}) where {T} = θ.T2
@inline T1time(θ::EPGIncreasingFlipAnglesOptions{T}) where {T} = θ.T1

const EPGParameterization{T,ETL} = Union{
    EPGOptions{T,ETL},
    EPGIncreasingFlipAnglesOptions{T,ETL},
}

#### Destructuring/restructuring to/from vectors

@generated function destructure(θ, ::Val{Fs}) where {Fs}
    vals = [:(getproperty(θ, $(QuoteNode(F)))) for F in Fs]
    :(Base.@_inline_meta; SVector{$(length(Fs)), $(eltype(θ))}(tuple($(vals...))))
end
destructure(θ, Fs::NTuple{N,Symbol}) where {N} = SVector{N, eltype(θ)}(map(F -> getproperty(θ, F), Fs))

@generated function restructure(θ, x::AbstractVector{T}, ::Val{Fs}) where {T, Fs}
    idxmap = NamedTuple{Fs}(ntuple(i -> i, length(Fs)))
    vals   = [F ∈ Fs ? :(@inbounds(x[$(getproperty(idxmap, F))])) : :($T(getproperty(θ, $(QuoteNode(F))))) for F in fieldsof(θ)]
    :(Base.@_inline_meta; restructure(θ, tuple($(vals...))))
end
function restructure(θ, x::AbstractVector{T}, Fs::NTuple{N,Symbol}) where {T, N}
    fields = fieldsof(typeof(θ))
    vals   = ntuple(length(θ)) do i
        @inbounds for j in 1:N
            (Fs[j] == fields[i]) && return x[j]
        end
        @inbounds T(getproperty(θ, fields[i]))
    end
    return restructure(θ, vals)
end
restructure(θ, x::NamedTuple{Fs}) where {Fs} = restructure(θ, SVector(Tuple(x)), Val(Fs))

####
#### Abstract Interface
####

abstract type AbstractEPGWorkspace{T,ETL} end

@inline Base.eltype(::AbstractEPGWorkspace{T}) where {T} = T
@inline echotrainlength(::AbstractEPGWorkspace{T,ETL}) where {T,ETL} = ETL
@inline mpsv(work::AbstractEPGWorkspace) = work.MPSV
@inline decaycurve(work::AbstractEPGWorkspace) = work.dc

@inline EPGdecaycurve_work(::EPGParameterization{T,ETL}) where {T,ETL} = EPGdecaycurve_work(T, ETL)
@inline EPGdecaycurve_work(::Type{T}, ETL::Int) where {T} = EPGWork_ReIm_DualMVector_Split(T, ETL) # fallback
@inline EPGdecaycurve_work(::Type{T}, ETL::Int) where {T <: FloatingTypes} = EPGWork_ReIm_DualMVector_Split(T, ETL) # default for T <: SIMD.FloatingTypes

"""
    EPGdecaycurve(ETL::Int, α::Real, TE::Real, T2::Real, T1::Real, β::Real)

Computes the normalized echo decay curve for a multi spin echo sequence
using the extended phase graph algorithm using the given input parameters.

The sequence of flip angles used is slight generalization of the standard
90 degree excitation pulse followed by 180 degree pulse train.
Here, the sequence used is `A*90, A*180, A*β, A*β, ...` where `A = α/180`
accounts for B1 inhomogeneities. Equivalently, the pulse sequence can
be written as `α/2, α, α * (β/180), α * (β/180), ...`.
Note that if `α = β = 180`, we recover the standard `90, 180, 180, ...`
pulse sequence.

# Arguments
- `ETL::Int`:   echo train length, i.e. number of echos
- `α::Real`:    angle of refocusing pulses (Units: degrees)
- `TE::Real`:   inter-echo time (Units: seconds)
- `T2::Real`:   transverse relaxation time (Units: seconds)
- `T1::Real`:   longitudinal relaxation time (Units: seconds)
- `β::Real`:    value of Refocusing Pulse Control Angle (Units: degrees)

# Outputs
- `decay_curve::AbstractVector`: normalized echo decay curve with length `ETL`
"""
@inline EPGdecaycurve(ETL::Int, α::Real, TE::Real, T2::Real, T1::Real, β::Real) = EPGdecaycurve(EPGOptions{floattype((α,TE,T2,T1,β)),ETL}((α,TE,T2,T1,β)))
@inline EPGdecaycurve(θ::EPGParameterization{T,ETL}) where {T,ETL} = EPGdecaycurve!(EPGdecaycurve_work(θ), θ)
@inline EPGdecaycurve!(work::AbstractEPGWorkspace{T,ETL}, args::Real...) where {T,ETL} = EPGdecaycurve!(decaycurve(work), work, EPGParameterization{T,ETL}(args...))
@inline EPGdecaycurve!(work::AbstractEPGWorkspace{T,ETL}, θ::EPGParameterization{T,ETL}) where {T,ETL} = EPGdecaycurve!(decaycurve(work), work, θ)
@inline EPGdecaycurve!(dc::AbstractVector{T}, work::AbstractEPGWorkspace{T,ETL}, θ::EPGParameterization{T,ETL}) where {T,ETL} = epg_decay_curve!(dc, work, θ)

####
#### Jacobian utilities (currently hardcoded for `EPGWork_ReIm_DualMVector_Split`)
####

struct EPGWorkCacheDict{ETL} <: AbstractDict{DataType, Any}
    dict::Dict{DataType, Any}
    EPGWorkCacheDict{ETL}() where {ETL} = new{ETL}(Dict{DataType, Any}())
end
@inline Base.keys(caches::EPGWorkCacheDict) = Base.keys(caches.dict)
@inline Base.values(caches::EPGWorkCacheDict) = Base.values(caches.dict)
@inline Base.length(caches::EPGWorkCacheDict) = Base.length(caches.dict)
@inline Base.iterate(caches::EPGWorkCacheDict, state...) = Base.iterate(caches.dict, state...)

@inline function Base.getindex(caches::EPGWorkCacheDict{ETL}, ::Type{T}) where {T,ETL}
    R = cachetype(caches, T)
    get!(caches.dict, T) do
        EPGWork_ReIm_DualMVector_Split(T, ETL)::R
    end::R
end
@inline cachetype(::EPGWorkCacheDict{ETL}, ::Type{T}) where {T,ETL} = EPGWork_ReIm_DualMVector_Split{T, ETL, MVector{ETL, SVector{3, T}}, MVector{ETL, T}}

#=
struct EPGWorkDualCache{T, D}
    work::T
    dual_work::D
end

function EPGWorkDualCache(::Type{T}, ::Val{ETL}, ::Val{chunk_size}) where {T, ETL, chunk_size}
    D = ForwardDiff.Dual{nothing, T, chunk_size}
    work = EPGWork_ReIm_DualVector_Split(T, ETL)
    dual_work = EPGWork_ReIm_DualVector_Split(D, ETL)
    EPGWorkDualCache(work, dual_work)
end

getindex(c::EPGWorkDualCache, ::Type{T}) where {T} = c.work
getindex(c::EPGWorkDualCache, ::Type{D}) where {D <: ForwardDiff.Dual} = remake(c.dual_work, D)

remake(x::Array, ::Type{T}) where {T} = x
remake(x::Array, ::Type{D}) where {D <: ForwardDiff.Dual} = reinterpret(D, x)

@inline function remake(work::EPGWork_ReIm_DualVector_Split{<:Any, ETL}, ::Type{D}) where {D <: ForwardDiff.Dual, ETL}
    MPSV₁ = reinterpret(SVector{3,D}, work.MPSV₁.data)
    MPSV₁ = SizedVector{ETL,SVector{3,D},typeof(MPSV₁)}(MPSV₁)
    MPSV₂ = reinterpret(SVector{3,D}, work.MPSV₂.data)
    MPSV₂ = SizedVector{ETL,SVector{3,D},typeof(MPSV₂)}(MPSV₂)
    dc    = reinterpret(D, work.dc.data)
    dc    = SizedVector{ETL,D,typeof(dc)}(dc)
    EPGWork_ReIm_DualVector_Split{D,ETL,typeof(MPSV₁),typeof(dc)}(MPSV₁, MPSV₂, dc)
end
=#

struct EPGFunctor{T, ETL, Fs, Tθ <: EPGParameterization{T,ETL}}
    caches::EPGWorkCacheDict{ETL}
    θ::Tθ
end
EPGFunctor(θ::EPGParameterization{T,ETL}, Fs::NTuple{N,Symbol}) where {T,ETL,N} = EPGFunctor{T,ETL,Fs,typeof(θ)}(EPGWorkCacheDict{ETL}(), θ)
EPGFunctor(f!::EPGFunctor{T,ETL,Fs}, θ::EPGParameterization{T,ETL}) where {T,ETL,Fs} = EPGFunctor{T,ETL,Fs,typeof(θ)}(f!.caches, θ)

@inline parameters(f!::EPGFunctor) = f!.θ
@inline optfields(::EPGFunctor{T,ETL,Fs}) where {T,ETL,Fs} = Val(Fs)

function (f!::EPGFunctor)(y::AbstractVector{D}, epg_work::AbstractEPGWorkspace{D,ETL}, x::AbstractVector{D}) where {D,ETL}
    θ = restructure(parameters(f!), x, optfields(f!))
    DECAES.EPGdecaycurve!(y, epg_work, θ)
end
(f!::EPGFunctor)(x::AbstractVector{D}) where {D} = f!(decaycurve(f!.caches[D]), f!.caches[D], x)
(f!::EPGFunctor)(y::AbstractVector{D}, x::AbstractVector{D}) where {D} = f!(y, f!.caches[D], x)

struct EPGJacobianFunctor{T, ETL, Fs, F <: EPGFunctor{T,ETL,Fs}, R <: DiffResults.DiffResult, C <: ForwardDiff.JacobianConfig}
    f!::F
    res::R
    cfg::C
end
function EPGJacobianFunctor(θ::EPGParameterization{T,ETL}, Fs::NTuple{N,Symbol}) where {T,ETL,N}
    f!  = EPGFunctor(θ, Fs)
    res = DiffResults.JacobianResult(zeros(T, ETL), zeros(T, N))
    cfg = ForwardDiff.JacobianConfig(f!, zeros(T, ETL), zeros(T, N), ForwardDiff.Chunk(N))
    return EPGJacobianFunctor(f!, res, cfg)
end

@inline parameters(j!::EPGJacobianFunctor) = parameters(j!.f!)
@inline optfields(j!::EPGJacobianFunctor) = optfields(j!.f!)

function (j!::EPGJacobianFunctor{T,ETL})(J::Union{AbstractMatrix, DiffResults.DiffResult}, y::AbstractVector{T}, θ::EPGParameterization{T,ETL}) where {T,ETL}
    @unpack f!, cfg = j!
    f! = EPGFunctor(f!, θ)
    x  = destructure(parameters(f!), optfields(f!))
    ForwardDiff.jacobian!(J, f!, y, x, cfg)
    return J isa AbstractMatrix ? J : DiffResults.jacobian(J)
end
(j!::EPGJacobianFunctor{T,ETL})(y::AbstractVector{T}, θ::EPGParameterization{T,ETL}) where {T,ETL} = j!(j!.res, y, θ)

####
#### EPGWork_Basic_Cplx
####

struct EPGWork_Basic_Cplx{T, ETL, MPSVType <: AbstractVector{SVector{3,Complex{T}}}, DCType <: AbstractVector{T}} <: AbstractEPGWorkspace{T,ETL}
    MPSV::MPSVType
    dc::DCType
end
function EPGWork_Basic_Cplx(T, ETL::Int)
    MSPV = SizedVector{ETL,SVector{3,Complex{T}}}(undef)
    dc   = SizedVector{ETL,T}(undef)
    EPGWork_Basic_Cplx{T,ETL,typeof(MSPV),typeof(dc)}(MSPV, dc)
end

# Compute a basis function under the extended phase graph algorithm. The magnetization phase state vector (MPSV) is
# successively modified by applying relaxation for TE/2, then a refocusing pulse as described by Hennig (1988),
# then transitioning phase states as given by Hennig (1988) but corrected by Jones (1997), and a finally relaxing for TE/2.
# See the appendix in Prasloski (2012) for details:
#    https://doi.org/10.1002/mrm.23157

function epg_decay_curve!(dc::AbstractVector{T}, work::EPGWork_Basic_Cplx{T,ETL}, θ::EPGOptions{T,ETL}) where {T,ETL}
    # Unpack workspace
    @unpack MPSV = work
    A   = B1correction(θ)
    αₑₓ = A * 90
    α₁  = A * 180
    αᵢ  = A * θ.β
    TE  = echotime(θ)
    T2  = T2time(θ)
    T1  = T1time(θ)
    V   = SA{Complex{T}} # alias

    # Precompute compute element flip matrices and other intermediate variables
    E1, E2 = exp(-(TE/2)/T1), exp(-(TE/2)/T2)
    E      = SA{T}[E2, E2, E1]
    R₁     = element_flipmat(α₁)
    Rᵢ     = element_flipmat(αᵢ)

    # Initialize magnetization phase state vector (MPSV)
    @inbounds for j in 1:ETL
        MPSV[j] = V[0, 0, 0]
    end
    @inbounds MPSV[1] = V[sind(αₑₓ), 0, 0] # initial magnetization in F1 state

    @inbounds for i in 1:ETL
        # Relaxation for TE/2, followed by flip matrix
        R = i == 1 ? R₁ : Rᵢ
        for j in 1:ETL
            MPSV[j] = R * (E .* MPSV[j])
        end

        # Transition between phase states
        Mᵢ, Mᵢ₊₁ = MPSV[1], MPSV[2]
        MPSV[1]  = V[Mᵢ[2], Mᵢ₊₁[2], Mᵢ[3]] # (F₁, F₁*, Z₁)⁺ = (F₁*, F₂*, Z₁)
        for j in 2:ETL-1
            Mᵢ₋₁, Mᵢ, Mᵢ₊₁ = Mᵢ, Mᵢ₊₁, MPSV[j+1]
            MPSV[j]        = V[Mᵢ₋₁[1], Mᵢ₊₁[2], Mᵢ[3]] # (Fᵢ, Fᵢ*, Zᵢ)⁺ = (Fᵢ₋₁, Fᵢ₊₁*, Zᵢ)
        end
        Mᵢ₋₁, Mᵢ  = Mᵢ, Mᵢ₊₁
        MPSV[ETL] = V[Mᵢ₋₁[1], 0, Mᵢ[3]] # (Fₙ, Fₙ*, Zₙ)⁺ = (Fₙ₋₁, 0, Zₙ)

        # Relaxation for TE/2
        for j in 1:ETL
            MPSV[j] = E .* MPSV[j] # Relaxation for TE/2
        end
        dc[i] = abs(MPSV[1][1]) # first echo amplitude
    end

    return dc
end

function epg_decay_curve!(dc::AbstractVector{T}, work::EPGWork_Basic_Cplx{T,ETL}, θ::EPGIncreasingFlipAnglesOptions{T,ETL}) where {T,ETL}
    # Unpack workspace
    @unpack MPSV = work
    A   = B1correction(θ)
    αₑₓ = A * flipangle(θ, 0)
    TE  = echotime(θ)
    T2  = T2time(θ)
    T1  = T1time(θ)
    V   = SA{Complex{T}} # alias

    # Precompute compute element flip matrices and other intermediate variables
    E1, E2 = exp(-(TE/2)/T1), exp(-(TE/2)/T2)
    E      = SA{T}[E2, E2, E1]

    # Initialize magnetization phase state vector (MPSV)
    @inbounds for j in 1:ETL
        MPSV[j] = V[0, 0, 0]
    end
    @inbounds MPSV[1] = V[sind(αₑₓ), 0, 0] # initial magnetization in F1 state

    @inbounds for i in 1:ETL
        # Relaxation for TE/2, followed by flip matrix
        R = element_flipmat(A * flipangle(θ, i))
        for j in 1:ETL
            MPSV[j] = R * (E .* MPSV[j])
        end

        # Transition between phase states
        Mᵢ, Mᵢ₊₁ = MPSV[1], MPSV[2]
        MPSV[1]  = V[Mᵢ[2], Mᵢ₊₁[2], Mᵢ[3]] # (F₁, F₁*, Z₁)⁺ = (F₁*, F₂*, Z₁)
        for j in 2:ETL-1
            Mᵢ₋₁, Mᵢ, Mᵢ₊₁ = Mᵢ, Mᵢ₊₁, MPSV[j+1]
            MPSV[j]        = V[Mᵢ₋₁[1], Mᵢ₊₁[2], Mᵢ[3]] # (Fᵢ, Fᵢ*, Zᵢ)⁺ = (Fᵢ₋₁, Fᵢ₊₁*, Zᵢ)
        end
        Mᵢ₋₁, Mᵢ  = Mᵢ, Mᵢ₊₁
        MPSV[ETL] = V[Mᵢ₋₁[1], 0, Mᵢ[3]] # (Fₙ, Fₙ*, Zₙ)⁺ = (Fₙ₋₁, 0, Zₙ)

        # Relaxation for TE/2
        for j in 1:ETL
            MPSV[j] = E .* MPSV[j] # Relaxation for TE/2
        end
        dc[i] = abs(MPSV[1][1]) # first echo amplitude
    end

    return dc
end

####
#### EPGWork_ReIm
####

struct EPGWork_ReIm{T, ETL, MPSVType <: AbstractVector{SVector{3,T}}, DCType <: AbstractVector{T}} <: AbstractEPGWorkspace{T,ETL}
    MPSV::MPSVType
    dc::DCType
end
function EPGWork_ReIm(T, ETL::Int)
    MSPV = SizedVector{ETL,SVector{3,T}}(undef)
    dc   = SizedVector{ETL,T}(undef)
    EPGWork_ReIm{T,ETL,typeof(MSPV),typeof(dc)}(MSPV, dc)
end

function epg_decay_curve!(dc::AbstractVector{T}, work::EPGWork_ReIm{T,ETL}, θ::EPGOptions{T,ETL}) where {T,ETL}
    # Unpack workspace
    @unpack MPSV = work
    A  = B1correction(θ)
    α₁ = deg2rad(A * 180)
    αᵢ = deg2rad(A * θ.β)
    TE = echotime(θ)
    T2 = T2time(θ)
    T1 = T1time(θ)
    V  = SA{T} # alias

    # Precompute intermediate variables
    E₁, E₂           = exp(-(TE/2)/T1), exp(-(TE/2)/T2)
    sin½α₁, cos½α₁   = sincos(α₁/2)
    sin²½α₁, cos²½α₁ = sin½α₁^2, cos½α₁^2
    sinα₁            = 2*sin½α₁*cos½α₁
    sinαᵢ, cosαᵢ     = sincos(αᵢ)
    cos²½αᵢ          = (1+cosαᵢ)/2
    sin²½αᵢ          = 1-cos²½αᵢ
    a₁, b₁, c₁       = E₂^2*cos²½α₁, E₂^2*sin²½α₁, E₁*E₂*sinα₁
    aᵢ, bᵢ, cᵢ, dᵢ   = E₂^2*cos²½αᵢ, E₂^2*sin²½αᵢ, E₁*E₂*sinαᵢ, E₁^2*cosαᵢ
    F, F̄, Z          = V[aᵢ, bᵢ, cᵢ], V[bᵢ, aᵢ, -cᵢ], V[-cᵢ/2, cᵢ/2, dᵢ]

    # Initialize magnetization phase state vector (MPSV), pulling i=1 iteration out of loop
    @inbounds begin
        m₀      = sin½α₁ # since αₑₓ = ½α₁
        Mᵢ⁺     = V[b₁*m₀, 0, -c₁*m₀/2]
        dc[1]   = abs(Mᵢ⁺[1])
        MPSV[1] = Mᵢ⁺
        MPSV[2] = V[a₁*m₀, 0, 0]
        MPSV[3] = V[0, 0, 0]
    end

    @inbounds for i in 2:ETL-1
        # j = 1, initialize and update `dc`
        Mᵢ, Mᵢ₊₁ = MPSV[1], MPSV[2]
        Mᵢ⁺      = V[F̄⋅Mᵢ, F̄⋅Mᵢ₊₁, Z⋅Mᵢ]
        dc[i]    = abs(Mᵢ⁺[1])
        MPSV[1]  = Mᵢ⁺

        # inner loop
        jup = min(i, ETL-i)
        @simd for j in 2:jup
            Mᵢ₋₁, Mᵢ, Mᵢ₊₁ = Mᵢ, Mᵢ₊₁, MPSV[j+1]
            MPSV[j]        = V[F⋅Mᵢ₋₁, F̄⋅Mᵢ₊₁, Z⋅Mᵢ]
        end

        # cleanup for next iteration
        if i == jup
            Mᵢ₋₁      = Mᵢ
            MPSV[i+1] = V[F⋅Mᵢ₋₁, 0, 0]
            MPSV[i+2] = V[0, 0, 0]
        end
    end

    @inbounds dc[ETL] = abs(F̄⋅MPSV[1])

    return dc
end

####
#### EPGWork_ReIm_Generated
####

struct EPGWork_ReIm_Generated{T, ETL, MPSVType <: AbstractVector{SVector{3,T}}, DCType <: AbstractVector{T}} <: AbstractEPGWorkspace{T,ETL}
    MPSV::MPSVType
    dc::DCType
end
function EPGWork_ReIm_Generated(T, ETL::Int)
    MSPV = SizedVector{ETL,SVector{3,T}}(undef)
    dc   = SizedVector{ETL,T}(undef)
    EPGWork_ReIm_Generated{T,ETL,typeof(MSPV),typeof(dc)}(MSPV, dc)
end

function epg_decay_curve_impl!(dc::Type{A}, work::Type{W}, θ::Type{O}) where {T, ETL, A<:AbstractVector{T}, W<:EPGWork_ReIm_Generated{T,ETL}, O<:EPGOptions{T,ETL}}
    MPSV(i::Int) = Symbol(:MPSV, i)
    quote
        # Unpack workspace
        A = B1correction(θ)
        α₁  = deg2rad(A * 180)
        αᵢ  = deg2rad(A * θ.β)
        TE  = echotime(θ)
        T2  = T2time(θ)
        T1  = T1time(θ)

        # Precompute intermediate variables
        # α                = deg2rad(α)
        # α₁, αᵢ           = α, α*β/180
        V                = SA{$T} # alias
        E₁, E₂           = exp(-(TE/2)/T1), exp(-(TE/2)/T2)
        sin½α₁,  cos½α₁  = sincos(α₁/2)
        sin²½α₁, cos²½α₁ = sin½α₁^2, cos½α₁^2
        sinα₁            = 2*sin½α₁*cos½α₁
        sinαᵢ, cosαᵢ     = sincos(αᵢ)
        cos²½αᵢ          = (1+cosαᵢ)/2
        sin²½αᵢ          = 1-cos²½αᵢ
        a₁, b₁, c₁       = E₂^2*cos²½α₁, E₂^2*sin²½α₁, E₁*E₂*sinα₁
        aᵢ, bᵢ, cᵢ, dᵢ   = E₂^2*cos²½αᵢ, E₂^2*sin²½αᵢ, E₁*E₂*sinαᵢ, E₁^2*cosαᵢ
        F, F̄, Z          = V[aᵢ, bᵢ, cᵢ], V[bᵢ, aᵢ, -cᵢ], V[-cᵢ/2, cᵢ/2, dᵢ]

        # Initialize MPSV vector elements
        $([
            :($(MPSV(i)) = zero(SVector{3,$T}))
            for i in 1:ETL
        ]...)

        # Initialize magnetization phase state vector (MPSV), pulling i=1 iteration out of loop
        @inbounds begin
            m₀         = sin½α₁ # since αₑₓ = ½α₁
            Mᵢ⁺        = V[b₁*m₀, 0, -c₁*m₀/2]
            dc[1]      = abs(Mᵢ⁺[1])
            $(MPSV(1)) = Mᵢ⁺
            $(MPSV(2)) = V[a₁*m₀, 0, 0]
        end

        # Main loop
        $([
            quote
                # Initialize and update `dc` (j = 1)
                @inbounds begin
                    Mᵢ, Mᵢ₊₁   = $(MPSV(1)), $(MPSV(2))
                    Mᵢ⁺        = V[F̄⋅Mᵢ, F̄⋅Mᵢ₊₁, Z⋅Mᵢ]
                    dc[$i]     = abs(Mᵢ⁺[1])
                    $(MPSV(1)) = Mᵢ⁺
                end

                # Inner loop
                $([
                    quote
                        (Mᵢ₋₁, Mᵢ, Mᵢ₊₁) = (Mᵢ, Mᵢ₊₁, $(MPSV(j+1)))
                        $(MPSV(j))       = V[F⋅Mᵢ₋₁, F̄⋅Mᵢ₊₁, Z⋅Mᵢ]
                    end
                    for j in 2:min(i, ETL-i)+1
                ]...)
            end
            for i in 2:ETL-1
        ]...)

        # Last echo
        @inbounds dc[$ETL] = abs(F̄⋅$(MPSV(1)))

        return dc
    end
end

@generated function epg_decay_curve!(dc::AbstractVector{T}, work::EPGWork_ReIm_Generated{T,ETL}, θ::EPGOptions{T,ETL}) where {T,ETL}
    return epg_decay_curve_impl!(dc, work, θ)
end

####
#### EPGWork_ReIm_DualVector
####

struct EPGWork_ReIm_DualVector{T, ETL, MPSVType <: AbstractVector{SVector{3,T}}, DCType <: AbstractVector{T}} <: AbstractEPGWorkspace{T,ETL}
    MPSV₁::MPSVType
    MPSV₂::MPSVType
    dc::DCType
end
function EPGWork_ReIm_DualVector(T, ETL::Int)
    MPSV₁ = SizedVector{ETL,SVector{3,T}}(undef)
    MPSV₂ = SizedVector{ETL,SVector{3,T}}(undef)
    dc    = SizedVector{ETL,T}(undef)
    EPGWork_ReIm_DualVector{T,ETL,typeof(MPSV₁),typeof(dc)}(MPSV₁, MPSV₂, dc)
end

function epg_decay_curve!(dc::AbstractVector{T}, work::EPGWork_ReIm_DualVector{T,ETL}, θ::EPGOptions{T,ETL}) where {T,ETL}
    # Unpack workspace
    @unpack MPSV₁, MPSV₂ = work
    A  = B1correction(θ)
    α₁ = deg2rad(A * 180)
    αᵢ = deg2rad(A * θ.β)
    TE = echotime(θ)
    T2 = T2time(θ)
    T1 = T1time(θ)
    V  = SA{T} # alias

    # Precompute intermediate variables
    # α                = deg2rad(α)
    # α₁, αᵢ           = α, α*β/180
    E₁, E₂           = exp(-(TE/2)/T1), exp(-(TE/2)/T2)
    sin½α₁, cos½α₁   = sincos(α₁/2)
    sin²½α₁, cos²½α₁ = sin½α₁^2, cos½α₁^2
    sinα₁            = 2*sin½α₁*cos½α₁
    sinαᵢ, cosαᵢ     = sincos(αᵢ)
    cos²½αᵢ          = (1 + cosαᵢ)/2
    sin²½αᵢ          = 1-cos²½αᵢ
    a₁, b₁, c₁       = E₂^2*cos²½α₁, E₂^2*sin²½α₁, E₁*E₂*sinα₁
    aᵢ, bᵢ, cᵢ, dᵢ   = E₂^2*cos²½αᵢ, E₂^2*sin²½αᵢ, E₁*E₂*sinαᵢ, E₁^2*cosαᵢ
    F, F̄, Z          = V[aᵢ, bᵢ, cᵢ], V[bᵢ, aᵢ, -cᵢ], V[-cᵢ/2, cᵢ/2, dᵢ]

    # Initialize magnetization phase state vector (MPSV), pulling i=1 iteration out of loop
    @inbounds begin
        m₀           = sin½α₁ # since αₑₓ = ½α₁
        Mᵢ⁺          = V[b₁*m₀, 0, -c₁*m₀/2]
        dc[1]        = abs(Mᵢ⁺[1])
        MPSV₁[1]     = Mᵢ⁺
        MPSV₁[2]     = V[a₁*m₀, 0, 0]
        MPSV₁[3]     = V[0, 0, 0]
        MPSV₁, MPSV₂ = MPSV₂, MPSV₁
    end

    @inbounds for i in 2:ETL-1
        # j = 1, initialize and update `dc`
        Mᵢ, Mᵢ₊₁ = MPSV₂[1], MPSV₂[2]
        Mᵢ⁺      = V[F̄⋅Mᵢ, F̄⋅Mᵢ₊₁, Z⋅Mᵢ]
        dc[i]    = abs(Mᵢ⁺[1])
        MPSV₁[1] = Mᵢ⁺

        # inner loop
        jup = min(i, ETL-i)
        @simd for j in 2:jup
            Mᵢ₋₁, Mᵢ, Mᵢ₊₁ = Mᵢ, Mᵢ₊₁, MPSV₂[j+1]
            MPSV₁[j]       = V[F⋅Mᵢ₋₁, F̄⋅Mᵢ₊₁, Z⋅Mᵢ]
        end

        # cleanup for next iteration
        if i == jup
            Mᵢ₋₁       = Mᵢ
            MPSV₁[i+1] = V[F⋅Mᵢ₋₁, 0, 0]
            MPSV₁[i+2] = V[0, 0, 0]
        end
        MPSV₁, MPSV₂ = MPSV₂, MPSV₁
    end

    @inbounds dc[ETL] = abs(F̄⋅MPSV₂[1])

    return dc
end

####
#### EPGWork_ReIm_DualVector_Split
####

struct EPGWork_ReIm_DualVector_Split{T, ETL, MPSVType <: AbstractVector{SVector{3,T}}, DCType <: AbstractVector{T}} <: AbstractEPGWorkspace{T,ETL}
    MPSV₁::MPSVType
    MPSV₂::MPSVType
    dc::DCType
end
function EPGWork_ReIm_DualVector_Split(T, ETL::Int)
    MPSV₁ = SizedVector{ETL,SVector{3,T}}(undef)
    MPSV₂ = SizedVector{ETL,SVector{3,T}}(undef)
    dc    = SizedVector{ETL,T}(undef)
    EPGWork_ReIm_DualVector_Split{T,ETL,typeof(MPSV₁),typeof(dc)}(MPSV₁, MPSV₂, dc)
end

function epg_decay_curve!(dc::AbstractVector{T}, work::EPGWork_ReIm_DualVector_Split{T,ETL}, θ::EPGOptions{T,ETL}) where {T,ETL}
    # Unpack workspace
    @unpack MPSV₁, MPSV₂ = work
    A  = B1correction(θ)
    α₁ = deg2rad(A * 180)
    αᵢ = deg2rad(A * θ.β)
    TE = echotime(θ)
    T2 = T2time(θ)
    T1 = T1time(θ)
    V  = SA{T} # alias

    # Precompute intermediate variables
    # α                = deg2rad(α)
    # α₁, αᵢ           = α, α*β/180
    E₁, E₂           = exp(-(TE/2)/T1), exp(-(TE/2)/T2)
    sin½α₁, cos½α₁   = sincos(α₁/2)
    sin²½α₁, cos²½α₁ = sin½α₁^2, cos½α₁^2
    sinα₁            = 2*sin½α₁*cos½α₁
    sinαᵢ, cosαᵢ     = sincos(αᵢ)
    cos²½αᵢ          = (1+cosαᵢ)/2
    sin²½αᵢ          = 1-cos²½αᵢ
    a₁, b₁, c₁       = E₂^2*cos²½α₁, E₂^2*sin²½α₁, E₁*E₂*sinα₁
    aᵢ, bᵢ, cᵢ, dᵢ   = E₂^2*cos²½αᵢ, E₂^2*sin²½αᵢ, E₁*E₂*sinαᵢ, E₁^2*cosαᵢ
    F, F̄, Z          = V[aᵢ, bᵢ, cᵢ], V[bᵢ, aᵢ, -cᵢ], V[-cᵢ/2, cᵢ/2, dᵢ]

    # Initialize magnetization phase state vector (MPSV), pulling i=1 iteration out of loop
    @inbounds begin
        m₀           = sin½α₁ # since αₑₓ = ½α₁
        Mᵢ⁺          = V[b₁*m₀, 0, -c₁*m₀/2]
        dc[1]        = abs(Mᵢ⁺[1])
        MPSV₁[1]     = Mᵢ⁺
        MPSV₁[2]     = V[a₁*m₀, 0, 0]
        MPSV₁, MPSV₂ = MPSV₂, MPSV₁
    end

    @inbounds for i in 2:ETL÷2
        Mᵢ, Mᵢ₊₁ = MPSV₂[1], MPSV₂[2] # j = 1, initialize and update `dc`
        Mᵢ⁺      = V[F̄⋅Mᵢ, F̄⋅Mᵢ₊₁, Z⋅Mᵢ]
        dc[i]    = abs(Mᵢ⁺[1])
        MPSV₁[1] = Mᵢ⁺
        @simd for j in 2:i-1
            Mᵢ₋₁, Mᵢ, Mᵢ₊₁ = Mᵢ, Mᵢ₊₁, MPSV₂[j+1]
            MPSV₁[j]       = V[F⋅Mᵢ₋₁, F̄⋅Mᵢ₊₁, Z⋅Mᵢ]
        end
        MPSV₁[i]     = V[F⋅Mᵢ, 0, Z⋅Mᵢ₊₁]
        MPSV₁[i+1]   = V[F⋅Mᵢ₊₁, 0, 0]
        MPSV₁, MPSV₂ = MPSV₂, MPSV₁
    end

    @inbounds for i in ETL÷2+1:ETL-1
        Mᵢ, Mᵢ₊₁ = MPSV₂[1], MPSV₂[2] # j = 1, initialize and update `dc`
        Mᵢ⁺      = V[F̄⋅Mᵢ, F̄⋅Mᵢ₊₁, Z⋅Mᵢ]
        dc[i]    = abs(Mᵢ⁺[1])
        MPSV₁[1] = Mᵢ⁺
        @simd for j in 2:ETL-i
            Mᵢ₋₁, Mᵢ, Mᵢ₊₁ = Mᵢ, Mᵢ₊₁, MPSV₂[j+1]
            MPSV₁[j]       = V[F⋅Mᵢ₋₁, F̄⋅Mᵢ₊₁, Z⋅Mᵢ]
        end
        MPSV₁, MPSV₂ = MPSV₂, MPSV₁
    end

    @inbounds dc[ETL] = abs(F̄⋅MPSV₂[1])

    return dc
end

####
#### EPGWork_ReIm_DualMVector_Split
####

struct EPGWork_ReIm_DualMVector_Split{T, ETL, MPSVType <: AbstractVector{SVector{3,T}}, DCType <: AbstractVector{T}} <: AbstractEPGWorkspace{T,ETL}
    MPSV₁::MPSVType
    MPSV₂::MPSVType
    dc::DCType
end
function EPGWork_ReIm_DualMVector_Split(T, ETL::Int)
    MPSV₁ = MVector{ETL,SVector{3,T}}(undef)
    MPSV₂ = MVector{ETL,SVector{3,T}}(undef)
    dc    = MVector{ETL,T}(undef)
    EPGWork_ReIm_DualMVector_Split{T,ETL,typeof(MPSV₁),typeof(dc)}(MPSV₁, MPSV₂, dc)
end

function epg_decay_curve!(dc::AbstractVector{T}, work::EPGWork_ReIm_DualMVector_Split{T,ETL}, θ::EPGOptions{T,ETL}) where {T,ETL}
    # Unpack workspace
    @unpack MPSV₁, MPSV₂ = work
    A  = B1correction(θ)
    α₁ = deg2rad(A * 180)
    αᵢ = deg2rad(A * θ.β)
    TE = echotime(θ)
    T2 = T2time(θ)
    T1 = T1time(θ)
    V  = SA{T} # alias

    # Precompute intermediate variables
    # α                = deg2rad(α)
    # α₁, αᵢ           = α, α*β/180
    E₁, E₂           = exp(-(TE/2)/T1), exp(-(TE/2)/T2)
    sin½α₁, cos½α₁   = sincos(α₁/2)
    sin²½α₁, cos²½α₁ = sin½α₁^2, cos½α₁^2
    sinα₁            = 2*sin½α₁*cos½α₁
    sinαᵢ, cosαᵢ     = sincos(αᵢ)
    cos²½αᵢ          = (1+cosαᵢ)/2
    sin²½αᵢ          = 1-cos²½αᵢ
    a₁, b₁, c₁       = E₂^2*cos²½α₁, E₂^2*sin²½α₁, E₁*E₂*sinα₁
    aᵢ, bᵢ, cᵢ, dᵢ   = E₂^2*cos²½αᵢ, E₂^2*sin²½αᵢ, E₁*E₂*sinαᵢ, E₁^2*cosαᵢ
    F, F̄, Z          = V[aᵢ, bᵢ, cᵢ], V[bᵢ, aᵢ, -cᵢ], V[-cᵢ/2, cᵢ/2, dᵢ]

    # Initialize magnetization phase state vector (MPSV), pulling i=1 iteration out of loop
    @inbounds begin
        m₀           = sin½α₁ # since αₑₓ = ½α₁
        Mᵢ⁺          = V[b₁*m₀, 0, -c₁*m₀/2]
        dc[1]        = abs(Mᵢ⁺[1])
        MPSV₁[1]     = Mᵢ⁺
        MPSV₁[2]     = V[a₁*m₀, 0, 0]
        MPSV₁, MPSV₂ = MPSV₂, MPSV₁
    end

    @inbounds for i in 2:ETL÷2
        Mᵢ, Mᵢ₊₁ = MPSV₂[1], MPSV₂[2] # j = 1, initialize and update `dc`
        Mᵢ⁺      = V[F̄⋅Mᵢ, F̄⋅Mᵢ₊₁, Z⋅Mᵢ]
        dc[i]    = abs(Mᵢ⁺[1])
        MPSV₁[1] = Mᵢ⁺
        @simd for j in 2:i-1
            Mᵢ₋₁, Mᵢ, Mᵢ₊₁ = Mᵢ, Mᵢ₊₁, MPSV₂[j+1]
            MPSV₁[j]       = V[F⋅Mᵢ₋₁, F̄⋅Mᵢ₊₁, Z⋅Mᵢ]
        end
        MPSV₁[i]     = V[F⋅Mᵢ, 0, Z⋅Mᵢ₊₁]
        MPSV₁[i+1]   = V[F⋅Mᵢ₊₁, 0, 0]
        MPSV₁, MPSV₂ = MPSV₂, MPSV₁
    end

    @inbounds for i in ETL÷2+1:ETL-1
        Mᵢ, Mᵢ₊₁ = MPSV₂[1], MPSV₂[2] # j = 1, initialize and update `dc`
        Mᵢ⁺      = V[F̄⋅Mᵢ, F̄⋅Mᵢ₊₁, Z⋅Mᵢ]
        dc[i]    = abs(Mᵢ⁺[1])
        MPSV₁[1] = Mᵢ⁺
        @simd for j in 2:ETL-i
            Mᵢ₋₁, Mᵢ, Mᵢ₊₁ = Mᵢ, Mᵢ₊₁, MPSV₂[j+1]
            MPSV₁[j]       = V[F⋅Mᵢ₋₁, F̄⋅Mᵢ₊₁, Z⋅Mᵢ]
        end
        MPSV₁, MPSV₂ = MPSV₂, MPSV₁
    end

    @inbounds dc[ETL] = abs(F̄⋅MPSV₂[1])

    return dc
end

####
#### EPGWork_ReIm_DualPaddedMVector_Vec_Split
####

struct EPGWork_ReIm_DualPaddedMVector_Vec_Split{T, ETL, MPSVType <: AbstractVector{Vec{4,T}}, DCType <: AbstractVector{T}} <: AbstractEPGWorkspace{T,ETL}
    MPSV₁::MPSVType
    MPSV₂::MPSVType
    dc::DCType
end
function EPGWork_ReIm_DualPaddedMVector_Vec_Split(T, ETL::Int)
    MPSV₁ = MVector{ETL,Vec{4,T}}(undef)
    MPSV₂ = MVector{ETL,Vec{4,T}}(undef)
    dc    = MVector{ETL,T}(undef)
    EPGWork_ReIm_DualPaddedMVector_Vec_Split{T,ETL,typeof(MPSV₁),typeof(dc)}(MPSV₁, MPSV₂, dc)
end

function epg_decay_curve!(dc::AbstractVector{T}, work::EPGWork_ReIm_DualPaddedMVector_Vec_Split{T,ETL}, θ::EPGOptions{T,ETL}) where {T,ETL}
    # Unpack workspace
    @unpack MPSV₁, MPSV₂ = work
    A  = B1correction(θ)
    α₁ = deg2rad(A * 180)
    αᵢ = deg2rad(A * θ.β)
    TE = echotime(θ)
    T2 = T2time(θ)
    T1 = T1time(θ)
    V  = Vec{4,T} # alias

    # Precompute intermediate variables
    # α                = deg2rad(α)
    # α₁, αᵢ           = α, α*β/180
    E₁, E₂           = exp(-(TE/2)/T1), exp(-(TE/2)/T2)
    sin½α₁, cos½α₁   = sincos(α₁/2)
    sin²½α₁, cos²½α₁ = sin½α₁^2, cos½α₁^2
    sinα₁            = 2*sin½α₁*cos½α₁
    sinαᵢ, cosαᵢ     = sincos(αᵢ)
    cos²½αᵢ          = (1+cosαᵢ)/2
    sin²½αᵢ          = 1-cos²½αᵢ
    a₁, b₁, c₁       = E₂^2*cos²½α₁, E₂^2*sin²½α₁, E₁*E₂*sinα₁
    aᵢ, bᵢ, cᵢ, dᵢ   = E₂^2*cos²½αᵢ, E₂^2*sin²½αᵢ, E₁*E₂*sinαᵢ, E₁^2*cosαᵢ
    F, F̄, Z          = V((aᵢ, bᵢ, cᵢ, 0)), V((bᵢ, aᵢ, -cᵢ, 0)), V((-cᵢ/2, cᵢ/2, dᵢ, 0))

    # Initialize magnetization phase state vector (MPSV), pulling i=1 iteration out of loop
    @inbounds begin
        m₀           = sin½α₁ # since αₑₓ = ½α₁
        Mᵢ⁺          = V((b₁*m₀, 0, -c₁*m₀/2, 0))
        dc[1]        = abs(Mᵢ⁺[1])
        MPSV₁[1]     = Mᵢ⁺
        MPSV₁[2]     = V((a₁*m₀, 0, 0, 0))
        MPSV₁, MPSV₂ = MPSV₂, MPSV₁
    end

    @inbounds for i in 2:ETL÷2
        Mᵢ, Mᵢ₊₁ = MPSV₂[1], MPSV₂[2] # j = 1, initialize and update `dc`
        Mᵢ⁺      = V((sum(F̄*Mᵢ), sum(F̄*Mᵢ₊₁), sum(Z*Mᵢ), 0))
        dc[i]    = abs(Mᵢ⁺[1])
        MPSV₁[1] = Mᵢ⁺
        @simd for j in 2:i-1
            Mᵢ₋₁, Mᵢ, Mᵢ₊₁ = Mᵢ, Mᵢ₊₁, MPSV₂[j+1]
            MPSV₁[j]       = V((sum(F*Mᵢ₋₁), sum(F̄*Mᵢ₊₁), sum(Z*Mᵢ), 0))
        end
        MPSV₁[i]     = V((sum(F*Mᵢ), 0, sum(Z*Mᵢ₊₁), 0))
        MPSV₁[i+1]   = V((sum(F*Mᵢ₊₁), 0, 0, 0))
        MPSV₁, MPSV₂ = MPSV₂, MPSV₁
    end

    @inbounds for i in ETL÷2+1:ETL-1
        Mᵢ, Mᵢ₊₁ = MPSV₂[1], MPSV₂[2] # j = 1, initialize and update `dc`
        Mᵢ⁺      = V((sum(F̄*Mᵢ), sum(F̄*Mᵢ₊₁), sum(Z*Mᵢ), 0))
        dc[i]    = abs(Mᵢ⁺[1])
        MPSV₁[1] = Mᵢ⁺
        @simd for j in 2:ETL-i
            Mᵢ₋₁, Mᵢ, Mᵢ₊₁ = Mᵢ, Mᵢ₊₁, MPSV₂[j+1]
            MPSV₁[j]       = V((sum(F*Mᵢ₋₁), sum(F̄*Mᵢ₊₁), sum(Z*Mᵢ), 0))
        end
        MPSV₁, MPSV₂ = MPSV₂, MPSV₁
    end

    @inbounds dc[ETL] = abs(sum(F̄*MPSV₂[1]))

    return dc
end

####
#### EPGWork_ReIm_DualPaddedVector_Split
####

struct EPGWork_ReIm_DualPaddedVector_Split{T, ETL, MPSVType <: AbstractVector{SVector{4,T}}, DCType <: AbstractVector{T}} <: AbstractEPGWorkspace{T,ETL}
    MPSV₁::MPSVType
    MPSV₂::MPSVType
    dc::DCType
end
function EPGWork_ReIm_DualPaddedVector_Split(T, ETL::Int)
    MPSV₁ = SizedVector{ETL,SVector{4,T}}(undef)
    MPSV₂ = SizedVector{ETL,SVector{4,T}}(undef)
    dc    = SizedVector{ETL,T}(undef)
    EPGWork_ReIm_DualPaddedVector_Split{T,ETL,typeof(MPSV₁),typeof(dc)}(MPSV₁, MPSV₂, dc)
end

function epg_decay_curve!(dc::AbstractVector{T}, work::EPGWork_ReIm_DualPaddedVector_Split{T,ETL}, θ::EPGOptions{T,ETL}) where {T,ETL}
    # Unpack workspace
    @unpack MPSV₁, MPSV₂ = work
    A  = B1correction(θ)
    α₁ = deg2rad(A * 180)
    αᵢ = deg2rad(A * θ.β)
    TE = echotime(θ)
    T2 = T2time(θ)
    T1 = T1time(θ)
    V  = SA{T} # alias

    # Precompute intermediate variables
    # α                = deg2rad(α)
    # α₁, αᵢ           = α, α*β/180
    E₁, E₂           = exp(-(TE/2)/T1), exp(-(TE/2)/T2)
    sin½α₁, cos½α₁   = sincos(α₁/2)
    sin²½α₁, cos²½α₁ = sin½α₁^2, cos½α₁^2
    sinα₁            = 2*sin½α₁*cos½α₁
    sinαᵢ, cosαᵢ     = sincos(αᵢ)
    cos²½αᵢ          = (1+cosαᵢ)/2
    sin²½αᵢ          = 1-cos²½αᵢ
    a₁, b₁, c₁       = E₂^2*cos²½α₁, E₂^2*sin²½α₁, E₁*E₂*sinα₁
    aᵢ, bᵢ, cᵢ, dᵢ   = E₂^2*cos²½αᵢ, E₂^2*sin²½αᵢ, E₁*E₂*sinαᵢ, E₁^2*cosαᵢ
    F, F̄, Z          = V[aᵢ, bᵢ, cᵢ, 0], V[bᵢ, aᵢ, -cᵢ, 0], V[-cᵢ/2, cᵢ/2, dᵢ, 0]

    # Initialize magnetization phase state vector (MPSV), pulling i=1 iteration out of loop
    @inbounds begin
        m₀           = sin½α₁ # since αₑₓ = ½α₁
        Mᵢ⁺          = V[b₁*m₀, 0, -c₁*m₀/2, 0]
        dc[1]        = abs(Mᵢ⁺[1])
        MPSV₁[1]     = Mᵢ⁺
        MPSV₁[2]     = V[a₁*m₀, 0, 0, 0]
        MPSV₁, MPSV₂ = MPSV₂, MPSV₁
    end

    @inbounds for i in 2:ETL÷2
        Mᵢ, Mᵢ₊₁ = MPSV₂[1], MPSV₂[2] # j = 1, initialize and update `dc`
        Mᵢ⁺      = V[F̄⋅Mᵢ, F̄⋅Mᵢ₊₁, Z⋅Mᵢ, 0]
        dc[i]    = abs(Mᵢ⁺[1])
        MPSV₁[1] = Mᵢ⁺
        @simd for j in 2:i-1
            Mᵢ₋₁, Mᵢ, Mᵢ₊₁ = Mᵢ, Mᵢ₊₁, MPSV₂[j+1]
            MPSV₁[j]       = V[F⋅Mᵢ₋₁, F̄⋅Mᵢ₊₁, Z⋅Mᵢ, 0]
        end
        MPSV₁[i]     = V[F⋅Mᵢ, 0, Z⋅Mᵢ₊₁, 0]
        MPSV₁[i+1]   = V[F⋅Mᵢ₊₁, 0, 0, 0]
        MPSV₁, MPSV₂ = MPSV₂, MPSV₁
    end

    @inbounds for i in ETL÷2+1:ETL-1
        Mᵢ, Mᵢ₊₁ = MPSV₂[1], MPSV₂[2] # j = 1, initialize and update `dc`
        Mᵢ⁺      = V[F̄⋅Mᵢ, F̄⋅Mᵢ₊₁, Z⋅Mᵢ, 0]
        dc[i]    = abs(Mᵢ⁺[1])
        MPSV₁[1] = Mᵢ⁺
        @simd for j in 2:ETL-i
            Mᵢ₋₁, Mᵢ, Mᵢ₊₁ = Mᵢ, Mᵢ₊₁, MPSV₂[j+1]
            MPSV₁[j]       = V[F⋅Mᵢ₋₁, F̄⋅Mᵢ₊₁, Z⋅Mᵢ, 0]
        end
        MPSV₁, MPSV₂ = MPSV₂, MPSV₁
    end

    @inbounds dc[ETL] = abs(F̄⋅MPSV₂[1])

    return dc
end

####
#### EPGWork_Vec
####

# Flip matrix and relaxation matrix steps are combined into one loop, and SIMD.jl `Vec` types are used instead of `Complex`.
# As this function is called many times during T2mapSEcorr, the micro-optimizations are worth the loss of code readability.
# See `EPGWork_Basic_Cplx` for a more readable, mathematically identicaly implementation.

struct EPGWork_Vec{T, ETL, MPSVType <: AbstractVector{Vec{2,T}}, DCType <: AbstractVector{T}} <: AbstractEPGWorkspace{T,ETL}
    MPSV::MPSVType
    dc::DCType
end
function EPGWork_Vec(T, ETL::Int)
    MSPV = SizedVector{3*ETL,Vec{2,T}}(undef)
    dc = SizedVector{ETL,T}(undef)
    EPGWork_Vec{T,ETL,typeof(MSPV),typeof(dc)}(MSPV, dc)
end

function epg_decay_curve!(dc::AbstractVector{T}, work::EPGWork_Vec{T,ETL}, θ::EPGOptions{T,ETL}) where {T,ETL}
    ###########################
    # Setup
    @unpack MPSV = work
    A  = B1correction(θ)
    α₁ = deg2rad(A * 180)
    αᵢ = deg2rad(A * θ.β)
    TE = echotime(θ)
    T2 = T2time(θ)
    T1 = T1time(θ)

    @inbounds begin
        # Initialize magnetization phase state vector (MPSV)
        E2, E1  = exp(-TE/T2), exp(-TE/T1)
        E2_half = exp(-(TE/2)/T2)
        m₀      = E2_half * sin(α₁/2) # initial population; since αₑₓ = α₁/2
        M1x     =  m₀ * cos(α₁/2)^2   # M1x, M1y, M1z are elements resulting from first refocusing pulse applied to [m₀, 0, 0]
        M1y     =  m₀ - M1x           # M1y = m₀ * sin(α₁/2)^2 = m₀ - m₀ * cos(α₁/2)^2 = m₀ - M1x
        M1z     = -m₀ * sin(α₁)/2     # Note: this is the imaginary part
        dc[1]   = E2_half * abs(M1y)  # first echo amplitude

        # Apply first relaxation matrix iteration on non-zero states
        MPSV[1] = Vec((E2 * M1y, zero(T)))
        MPSV[2] = zero(Vec{2,T})
        MPSV[3] = Vec((zero(T), E1 * M1z))
        MPSV[4] = Vec((E2 * M1x, zero(T)))

        # Extract matrix elements + initialize temporaries
        a1, a2, a3, a4, a5 = sin(αᵢ), cos(αᵢ), sin(αᵢ/2)^2, cos(αᵢ/2)^2, sin(αᵢ)/2 # independent elements of T2mat
        b1, b2, b3, b4, b5 = E2*a1, E1*a2, E2*a3, E2*a4, E1*a5
        c1, c3, c4         = E2_half*a1, E2_half*a3, E2_half*a4
        b1F, b5F, c1F      = Vec((-b1, b1)), Vec((-b5, b5)), Vec((-c1, c1))
        Mz3                = MPSV[3]
    end

    @inbounds for i = 2:ETL-1
        ###########################
        # Unroll first flipmat/relaxmat iteration
        Vx, Vy  = MPSV[1], MPSV[2]
        c1z     = shufflevector(c1F * Mz3, Val((1,0)))
        Mz2     = muladd(c3, Vx, muladd(c4, Vy, -c1z)) # flipmat: 2 -> dc
        Mz4     = muladd(b4, Vx, muladd(b3, Vy, E2_half * c1z)) # relaxmat: 1 -> 4, save in buffer
        dc[i]   = sqrt(sum(Mz2 * Mz2)) # decay curve coefficient
        MPSV[1] = E2_half * Mz2 # relaxmat: 2 -> 1
        b5xy    = shufflevector(b5F * (Vx - Vy), Val((1,0)))
        Mz3     = muladd(b2, Mz3, b5xy) # relaxmat: 3 -> 3, save in buffer

        ###########################
        # flipmat + relaxmat loop
        @inbounds for j in 4:3:3*min(i-1, ETL)
            Vx, Vy, Vz = MPSV[j], MPSV[j+1], MPSV[j+2]
            b1z        = shufflevector(b1F * Vz, Val((1,0)))
            MPSV[j  ]  = Mz4 # relaxmat: assign forward, j -> j+3
            Mz4        = muladd(b4, Vx, muladd(b3, Vy,  b1z))
            MPSV[j-2]  = muladd(b3, Vx, muladd(b4, Vy, -b1z)) # relaxmat: assign backwards, j+1 -> j+1-3
            b5xy       = shufflevector(b5F * (Vx - Vy), Val((1,0)))
            MPSV[j+2]  = muladd(b2, Vz, b5xy) # relaxmat: j+2 -> j+2
        end

        ###########################
        # cleanup + zero next elements
        j         = 3i-2
        Vx        = MPSV[j]
        MPSV[j  ] = Mz4 # relaxmat: assign forward, j -> j+3
        MPSV[j-2] = b3 * Vx # relaxmat: assign backwards, j+1 -> j+1-3
        MPSV[j+2] = shufflevector(b5F * Vx, Val((1,0))) # relaxmat: j+2 -> j+2
        MPSV[j+3] = b4 * Vx # relaxmat: assign forward, j -> j+3
        MPSV[j+1] = Vec((zero(T), zero(T))) # relaxmat: assign backwards, j+1 -> j+1-3
        MPSV[j+5] = Vec((zero(T), zero(T))) # relaxmat: j+2 -> j+2
    end

    ###########################
    # decay curve coefficient
    @inbounds begin
        c1z     = shufflevector(c1F * Mz3, Val((1,0)))
        Mz2     = muladd(c3, MPSV[1], muladd(c4, MPSV[2], -c1z)) # last iteration of flipmat unrolled
        dc[end] = sqrt(sum(Mz2 * Mz2))
    end

    return dc
end

####
#### Algorithm list
####

const EPG_Algorithms = Any[
    EPGWork_Basic_Cplx,
    EPGWork_Vec,
    EPGWork_ReIm,
    EPGWork_ReIm_DualVector,
    EPGWork_ReIm_DualVector_Split,
    EPGWork_ReIm_DualMVector_Split,
    EPGWork_ReIm_DualPaddedMVector_Vec_Split,
    EPGWork_ReIm_DualPaddedVector_Split,
    EPGWork_ReIm_Generated,
]

ext: jl


text:
function simulate(data,p,v,m,nsteps,dt,iprint)
  f = [ zeros(2) for i in 1:data.N ]
  for i in 1:nsteps
    u, f = energyandforce!(p,f,data)
    for ip in 1:length(p)
      @. p[i] = p[i] + v[i]*dt + f[i]*dt^2/(2*m)
      @. v[i] = v[i] + f[i]*dt/m
    end
    if i%iprint == 0
      println(" Step $i, Energy = $u ")
    end
  end
end


text:
import PotentialFlow.Utils: @get, MappedVector

@testset "Utils" begin
    @testset "get" begin
        z = rand(ComplexF64)

        @get z (re, im) (r, i)
        @test i == imag(z)
        @test r == real(z)

    end

    @testset "MappedVector" begin
        x = [π, 0.0, π]
        y = MappedVector(cos, x, 1)
        @test y[0] == -1.0
        @test y[1] == 1.0
        @test y[2] == -1.0

        buff = IOBuffer()
        show(buff, y)
        @test String(take!(buff)) == "Vector{Float64} → typeof(cos) (0:2)"
        #@test String(take!(buff)) == "Array{Float64,1} → typeof(cos) (0:2)"

    end
end


text:
using KissThreading
using Test
using Statistics
using Random
using Random: GLOBAL_RNG

include("test_tmap.jl")
include("test_tmapreduce.jl")
include("bootstrap.jl")
include("bubblesort.jl")
include("sort_batch.jl")
include("summation.jl")
include("mapreduce.jl")

println("========")


text:
using BinaryBuilder

try
    @show readdir("products32")
    @show readdir("products32/bin")
end

file, hash = package(Prefix("products32"), "MySQL"; platform=Windows(:i686), verbose=true, force=true)
file2, hash2 = package(Prefix("products64"), "MySQL"; platform=Windows(:x86_64), verbose=true, force=true)

@show readdir()

println("""
Windows(:i686) => (\"\$bin_prefix/$file\", \"$hash\"),
Windows(:x86_64) => (\"\$bin_prefix/$file2\", \"$hash2\")
""")


text:
# ---
# title: 1117. Building H2O
# id: problem1117
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/building-h2o/description/>
# hidden: true
# ---
#
# There are two kinds of threads, `oxygen` and `hydrogen`. Your goal is to group
# these threads to form water molecules. There is a barrier where each thread
# has to wait until a complete molecule can be formed. Hydrogen and oxygen
# threads will be given `releaseHydrogen` and `releaseOxygen` methods
# respectively, which will allow them to pass the barrier. These threads should
# pass the barrier in groups of three, and they must be able to immediately bond
# with each other to form a water molecule. You must guarantee that all the
# threads from one molecule bond _before_ any other threads from the next
# molecule do.
#
# In other words:
#
#   * If an oxygen thread arrives at the barrier when no hydrogen threads are present, it has to wait for two hydrogen threads.
#   * If a hydrogen thread arrives at the barrier when no other threads are present, it has to wait for an oxygen thread and another hydrogen thread.
#
# We don't have to worry about matching the threads up explicitly; that is, the
# threads do not necessarily know which other threads they are paired up with.
# The key is just that threads pass the barrier in complete sets; thus, if we
# examine the sequence of threads that bond and divide them into groups of
# three, each group should contain one oxygen and two hydrogen threads.
#
# Write synchronization code for oxygen and hydrogen molecules that enforces
# these constraints.
#
#
#
# **Example 1:**
#
#
#
#     Input: "HOH"
#     Output: "HHO"
#     Explanation: "HOH" and "OHH" are also valid answers.
#
#
# **Example 2:**
#
#
#
#     Input: "OOHHHH"
#     Output: "HHOHHO"
#     Explanation: "HOHHHO", "OHHHHO", "HHOHOH", "HOHHOH", "OHHHOH", "HHOOHH", "HOHOHH" and "OHHOHH" are also valid answers.
#
#
#
#
# **Constraints:**
#
#   * Total length of input string will be 3 _n_ , where 1 ≤  _n_  ≤ 20.
#   * Total number of `H` will be 2 _n_  in the input string.
#   * Total number of `O` will be _n_  in the input string.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Astroshaper
using Documenter

makedocs(sitename="Astroshaper.jl")

max_stars_repo_name: MasanoriKanamaru/Astroshaper

max_issues_repo_name: MasanoriKanamaru/Astroshaper

max_forks_repo_name: MasanoriKanamaru/Astroshaper


text:
function test_volume(npoints = 100, dim = 3)
    x = rand(dim, npoints)
    @time v, p = voronoi_random(x, 10000)
    @time a = adjacency(v)
    @time av = [boundary_area_verts(a, v, p) for a in a]
    @time ae = [boundary_area_edges(a, p) for a in a]
    delta = filter(isfinite, ae - av)
    quota = sum(delta .< 1e-6) / length(delta)
    @assert quota > 2/3
    @assert length(delta) / length(av) > 2/3
end

function test_random(d=6, n=200)
    x = rand(d, n)
    @time v,p = voronoi_random(rand(d,n), 10_000_000; maxstuck=100_000);
    @show length(v)
    a = adjacency(v)
    println("avg. no. of neighbors: ", length(a)/n)
    println("avg. no. of vertices per face: ", mean(length.(values(a))))
    @time c = connectivity_matrix(v, p)
    v,p,c
end

function test_grid(n=5, iter=10000)
    plot(legend=false);
    x = hcat(hexgrid(n)...)
    x .+= randn(2,n*n) .* 0.01
    @time v, P  = voronoi_random(x, iter)

    v = Dict(filter(collect(v)) do (k,v)
        norm(v) < 10
        end)

    #c = extractconn(v)
    @time A, Vs = connectivity_matrix(v, P)

    AA = map(x->x>.0, A)
    plot_connectivity!(AA .* 2, P)
    scatter!(eachrow(hcat(values(v)...))...)
    xlims!(1,6); ylims!(0,5)
end

#=
function benchmark(n=100, d=6, iter=100, particles=10)
    x = rand(d, n)
    @benchmark voronoi($x, $iter, $particles)
end
=#

function hexgrid(n)
    P = []
    for i in 1:n
        for j in 1:n
            x = i + j/2
            y = j * sqrt(3) / 2
            push!(P, [x,y])
        end
    end
    P
end

function tests()
    test_volume()
    test_random()
    test_grid()
end


text:
using Gnuplot
t = 0:0.001:1
@gp t sin.(10π*t) "w l tit 'sin' lc 'black'"
#save("plt2_ex2.gp") # hide
#save(term="pdf size 5,3", output="plt2_ex2.pdf") # hide
save(term="pngcairo size 600,400", output="./code/plt2_ex2.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
#Gingrich stocking chart calculated using the original upland oak parameters
#for more information see:http://oak.snr.missouri.edu/silviculture/tools/gingrich.html

using RecipesBase
## eqs for gingrich chart
amd_convert(qmd) = -0.259 + 0.973qmd
get_tpa(j, i, a) = (1*j*10)/(a[1]+a[2]*amd_convert(i)+a[3]*i^2)

function stklines(stk, dia, a)
  tpa = get_tpa.(stk', dia, Ref(a))
  ba = pi * (dia ./ 24).^2 .* tpa
  tpa, ba
end

@userplot gingrich_chart
@recipe function f(dat::gingrich_chart)

  #define parameters
  a = [-0.0507, 0.1698, 0.0317]
  b = [0.175, 0.205, 0.06]
  #QMD lines
  dia = [7, 8, 10, 12, 14, 16, 18, 20, 22]
  #stocking percent horizontal lines
  stk_percent= 20:10:110

  # plot attributes
  xlabel --> "Trees Per Acre"
  ylabel --> "Basal Area (sq ft./acre)"
  xticks --> 0:50:450
  yticks --> 0:20:200
  legend := false
  primary := false

  tpa, ba = stklines(stk_percent, dia, a)

  # horizontal lines
  @series begin
    linecolor := :black
    tpa, ba
  end

  # vertical lines
  @series begin
    linecolor := :black
    tpa', ba'
  end

  # upper red line
  @series begin
    linecolor := :red
    stklines([100], dia, a)
  end

  # lower red line
  @series begin
    linecolor := :red
    stklines([100], dia, b)
  end

  # the point and the annotations
  @series begin
    primary := true
    seriestype := :scatter
    annotations := (100:38:442, 20:10:110, ["$i%" for i in 20:10:110])

    x,y = dat.args
    [x], [y]

  end
end

function gingrich_chart(s::Stand)
    gingrich_chart(s.trees_per_area,s.basal_area)
end

##Reineke SDI chart
@userplot reineke_chart
@recipe function f(dat::reineke_chart ; maxsdi=450)
diarng = 1:1:50
#tparng = [1,1200] #unusued for now, TODO: figure out how to autoscale plot?
maxline() = maxsdi ./ (diarng ./ 10.0).^1.605

  # plot attributes
  xlabel --> "Trees Per Acre"
  ylabel --> "Quadratic Mean Diameter"

  xticks --> [1,5,10,50,100,500,1000]
  yticks --> [1,2,5,10,20]

  xscale --> :log10
  yscale --> :log10

  legend := false
  primary := false
  formatter --> x -> string(convert(Int,round(x,digits=2)))

  title --> "Maximum SDI of $maxsdi"

#max SDI line
@series begin
    linecolor := :black
    maxline(), diarng
end

#competition induced mortality SDI line
@series begin
    linecolor := :orange
    maxline() * 0.55, diarng
end

#crown closure SDI line
@series begin
    linecolor := :red
    maxline() * 0.35, diarng
end

 # the point and the annotations
 @series begin
    primary := true
    seriestype := :scatter
    #annotations := (100:38:442, 20:10:110, ["$i%" for i in 20:10:110])
    x,y = dat.args
    [x], [y]
 end
end

function reineke_chart(s::Stand)
    reineke_chart(s.trees_per_area, s.qmd)
end


text:
export iauC2ixys
"""
Form the celestial to intermediate-frame-of-date matrix given the CIP
X,Y and the CIO locator s.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   x,y      double         Celestial Intermediate Pole (Note 1)
   s        double         the CIO locator s (Note 2)

Returned:
   rc2i     double[3][3]   celestial-to-intermediate matrix (Note 3)

Notes:

   1. The Celestial Intermediate Pole coordinates are the x,y
   components of the unit vector in the Geocentric Celestial
   Reference System.

   2. The CIO locator s (in radians) positions the Celestial
   Intermediate Origin on the equator of the CIP.

   3. The matrix rc2i is the first stage in the transformation from
   celestial to terrestrial coordinates:

      [TRS] = RPOM * R_3(ERA) * rc2i * [CRS]

            = RC2T * [CRS]

   where [CRS] is a vector in the Geocentric Celestial Reference
   System and [TRS] is a vector in the International Terrestrial
   Reference System (see IERS Conventions 2003), ERA is the Earth
   Rotation Angle and RPOM is the polar motion matrix.

Called:
   iauIr        initialize r-matrix to identity
   iauRz        rotate around Z-axis
   iauRy        rotate around Y-axis

Reference:

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

This revision:  2014 November 7

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauC2ixys(double x, double y, double s, double rc2i[3][3])

function iauC2ixys(x::Real, y::Real, s::Real)

   # Allocate return value
   rc2i = zeros(Float64, 3, 3)

   ccall((:iauC2ixys, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Cdouble,
         Ptr{Cdouble}),
         convert(Float64, x), convert(Float64, y), convert(Float64, s),
         rc2i)

   return SMatrix{3,3}(rc2i')
end