text:
"""
    System{TF, TV<:AbstractVector{TF}, TM<:AbstractMatrix{TF}}

Contains the system state, residual vector, and jacobian matrices as well as
pointers to be able to access their contents.  Also contains additional storage
needed for time domain simulations.

# Fields:
 - `static`: Flag indicating whether system matrices are only valid for static analyses
 - `x`: State vector
 - `r`: Residual vector
 - `K`: System jacobian matrix with respect to the state variables
 - `M`: System jacobian matrix with respect to the time derivative of the state variables
 - `force_scaling`: Scaling for state variables corresponding to forces/moments
 - `irow_point`: Row index of first equilibrium equation for each point
 - `irow_elem`: Row index of first equation for just this beam element
 - `irow_elem1`: Row index of first equation for the left side of each beam
 - `irow_elem2`: Row index of first equation for the right side of each beam
 - `icol_point`: Row/Column index of first state variable for each point
 - `icol_elem`: Row/Column index of first state variable for each beam element
 - `udot`: Time derivative of state variable `u` for each beam element
 - `θdot`: Time derivative of state variable `θ` for each beam element
 - `Vdot`: Time derivative of state variable `V` for each beam element
 - `Ωdot`: Time derivative of state variable `Ω` for each beam element
 - `t`: Current system time
"""
mutable struct System{TF, TV<:AbstractVector{TF}, TM<:AbstractMatrix{TF}}
    static::Bool
    x::TV
    r::TV
    K::TM
    M::TM
    force_scaling::TF
    irow_point::Vector{Int}
    irow_elem::Vector{Int}
    irow_elem1::Vector{Int}
    irow_elem2::Vector{Int}
    icol_point::Vector{Int}
    icol_elem::Vector{Int}
    udot::Vector{SVector{3,TF}}
    θdot::Vector{SVector{3,TF}}
    Vdot::Vector{SVector{3,TF}}
    Ωdot::Vector{SVector{3,TF}}
    t::TF
end
Base.eltype(::System{TF, TV, TM}) where {TF, TV, TM} = TF

"""
    System([TF=eltype(assembly),] assembly, static; kwargs...)

Initialize an object of type `System` which stores the system state.

# Arguments:
 - `TF:`(optional) Floating point type, defaults to the floating point type of `assembly`
 - `assembly`: Assembly of rigidly connected nonlinear beam elements
 - `static`: Flag indicating whether the system corresponds to a static system.

# Keyword Arguments
 - `prescribed_points`: Point indices corresponding to the points whose equations
    and state variables should be included in the system of equations.  By default,
    all point indices are included in the system of equations.
 - `force_scaling`: Factor used to scale system forces/moments internally.  If
    not specified, a suitable default will be chosen based on the entries of the
    compliance matrix.

Note that points with prescribed conditions must be included in the system of
equations.
"""
function System(assembly, static; kwargs...)

    return System(eltype(assembly), assembly, static; kwargs...)
end

function System(TF, assembly, static;
    prescribed_points = 1:length(assembly.points),
    force_scaling = default_force_scaling(assembly)
    )

    # system dimensions
    npoint = length(assembly.points)
    nelem = length(assembly.elements)

    # initialize system pointers
    N, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem =
        system_indices(assembly.start, assembly.stop, static; prescribed_points)

    # initialize system matrices
    x = zeros(TF, N)
    r = zeros(TF, N)
    K = spzeros(TF, N, N)
    M = spzeros(TF, N, N)

    # initialize storage for time domain simulations
    udot = [@SVector zeros(TF, 3) for i = 1:nelem]
    θdot = [@SVector zeros(TF, 3) for i = 1:nelem]
    Vdot = [@SVector zeros(TF, 3) for i = 1:nelem]
    Ωdot = [@SVector zeros(TF, 3) for i = 1:nelem]

    # initialize current time
    t = 0.0

    # set system types
    TV = promote_type(typeof(x), typeof(r))
    TM = promote_type(typeof(K), typeof(M))

    return System{TF, TV, TM}(static, x, r, K, M, force_scaling,
        irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem, udot, θdot,
        Vdot, Ωdot, t)
end

function default_force_scaling(assembly)

    TF = eltype(assembly)

    # create vector of all compliance entries
    compliance_entries = vcat([elem.compliance for elem in assembly.elements]...)

    # get all nonzero entries
    compliance_nonzero_indices = findall(xi -> abs(xi) > eps(TF), compliance_entries)

    # set force scaling based on nonzero compliance matrix entries
    if isempty(compliance_nonzero_indices)
        force_scaling = 1.0
    else
        nonzero_compliance_entries = compliance_entries[compliance_nonzero_indices]
        force_scaling = nextpow(2.0, length(nonzero_compliance_entries)/sum(abs, nonzero_compliance_entries)/100)
    end

    return force_scaling
end

"""
    system_indices(start, stop, static; kwargs...)

Return indices for accessing the equations and state variables associated with
each point and beam element in a system given its connectivity.

# Arguments:
 - `start`: Vector containing the point indices where each beam element starts
 - `stop`: Vector containing the point indices where each beam element stops
 - `static`: Flag indicating whether the analysis is static (rather than dynamic).
    Defaults to `false`.

# Keyword Arguments:
 - `prescribed_points`: Point indices corresponding to the points whose equations
    and state variables should be included in the system of equations.  By default,
    all point indices are included in the system of equations.

# Return Arguments:
 - `N`: total number of equations and unknowns in the system
 - `irow_point`: Row index of the first equation corresponding to each point
 - `irow_elem`: Row index of the first equation corresponding to each beam element
 - `irow_elem1`: Row index of the first equation corresponding to the start of
        each beam element
 - `irow_elem2`: Row index of the first equation corresponding to the end of
        each beam element
 - `icol_point`: Column index of the first state variable corresponding to each
        point
 - `icol_elem`: Column index of the first state variable corresponding to each
        beam element

Negative indices indicate that the equations and/or state variables associated
with the point/beam element have been omitted from the system of equations.
"""
function system_indices(start, stop, static;
    prescribed_points = 1:max(maximum(start), maximum(stop)))

    npoint = max(maximum(start), maximum(stop))
    nelem = length(start)

    keep = [i in prescribed_points for i = 1:npoint]

    add_necessary_points!(keep, start, stop)

    # indicates whether a point has associated equations and state variables
    assigned = fill(false, npoint)

    # initialize pointers for equations
    irow_point = Vector{Int}(undef, npoint)
    irow_elem = Vector{Int}(undef, nelem)
    irow_elem1 = Vector{Int}(undef, nelem)
    irow_elem2 = Vector{Int}(undef, nelem)

    # initialize pointers for state variables
    icol_point = Vector{Int}(undef, npoint)
    icol_elem = Vector{Int}(undef, nelem)

    irow = 1
    icol = 1
    for ielem = 1:nelem
        # add state variables/equations for the start of the beam element
        ipt = start[ielem]

        # check if the point has associated equations and state variables
        if !assigned[ipt]
            # add equations and state variables for this point
            assigned[ipt] = true

            # add 6 equilibrium equations + 6 compatability equations
            irow_point[ipt] = irow
            irow_elem1[ielem] = irow
            irow += 12

            # add 6 state variables for this point (if necessary)
            if keep[ipt]
                icol_point[ipt] = icol
                icol += 6
            else
                icol_point[ipt] = -1
            end
        else
            # add additional equations for this point

            # add compatibility equations for this point (if necessary)
            if keep[ipt]
                irow_elem1[ielem] = irow
                irow += 6
            else
                irow_elem1[ielem] = -1
            end
        end

        # add state variables/equations for the beam element

        # add 12 state variables for this element
        icol_elem[ielem] = icol
        icol += 12

        # add an additional 6 state variables and equations for this element (if necessary)
        if static
            irow_elem[ielem] = -1
        else
            irow_elem[ielem] = irow
            irow += 6
            icol += 6
        end

        # add state variables/equations for the end of the beam element
        ipt = stop[ielem]

        # check if the point has associated equations and state variables
        if !assigned[ipt]
            # add equations and state variables for this point
            assigned[ipt] = true

            # add 6 equilibrium equations + 6 compatability equations
            irow_point[ipt] = irow
            irow_elem2[ielem] = irow
            irow += 12

            # add 6 state variables for this point (if necessary)
            if keep[ipt]
                icol_point[ipt] = icol
                icol += 6
            else
                icol_point[ipt] = -1
            end
        else
            # add additional equations for this point

            # add compatibility equations for this point (if necessary)
            if keep[ipt]
                irow_elem2[ielem] = irow
                irow += 6
            else
                irow_elem2[ielem] = -1
            end
        end
    end

    # number of state variables/equations
    N = irow - 1

    return N, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem
end

function add_necessary_points!(keep, start, stop)

    npoint = length(keep)

    for ipt = 1:npoint
        # calculate number of connections to this point
        ncon = count(x -> x == ipt, start) + count(x -> x == ipt, stop)

        if ncon != 2
            # this point must be preserved in the system of equations
            if keep[ipt] != true
                # only mutate if necessary
                keep[ipt] = true
            end
        end
    end

    return keep
end

"""
    system_state(system)

Return a vector containing the state variables of `system`.
"""
system_state(system) = system.x

"""
    set_state!([x,] system, prescribed_conditions; kwargs...)

Set the state variables in `system` (or in the vector `x`) to the provided values.
If values are not provided for a given keyword argument, then the state variables
corresponding to the keyword argument are not updated.

# Keyword Arguments
 - `u_e`: Vector containing initial displacements (in the body frame) for each beam element.
 - `theta_e`: Vector containing rotation variables (in the body frame) for each beam element.
 - `F_e`: Vector containing resultant forces (in the deformed local beam element frame) for each beam element.
 - `M_e`: Vector containing resultant moments (in the deformed local beam element frame) for each beam element.
 - `V_e`: Vector containing linear velocity (in the deformed local beam element frame) for each beam element.
 - `Ω_e`: Vector containing angular velocity (in the deformed local beam element frame) for each beam element.
 - `u_p`: Vector containing initial displacements (in the body frame) for each point.
 - `theta_p`: Vector containing rotation variables (in the body frame) for each point.
 - `F_p`: Vector containing externally applied forces (in the body frame) for each point.
 - `M_p`: Vector containing externally applied moments (in the body frame) for each point.
"""
set_state!

function set_state!(system, prescribed_conditions; kwargs...)
    x = set_state!(system.x, system, prescribed_conditions; kwargs...)
    return system
end

function set_state!(x, system, prescribed_conditions; kwargs...)
    return set_state!(x, system.icol_elem, system.icol_point; kwargs...)
end

function set_state!(x, icol_elem, icol_point, force_scaling,
    prescribed_conditions; u_e = nothing, theta_e = nothing,
    F_e = nothing, M_e = nothing, V_e = nothing, Ω_e = nothing,
    u_p = nothing, theta_p = nothing, F_p = nothing, M_p = nothing)

    nelem = length(icol_elem)
    npoint = length(icol_point)

    for ielem = 1:nelem
        icol = icol_elem[ielem]
        if !isnothing(u_e)
            set_element_deflection!(x, icol, u_e)
        end
        if !isnothing(theta_e)
            set_element_rotation!(x, icol, theta_e)
        end
        if !isnothing(F_e)
            set_element_forces!(x, icol, F_e, force_scaling)
        end
        if !isnothing(M_e)
            set_element_moments!(x, icol, M_e, force_scaling)
        end
        if !isnothing(V_e)
            set_element_linear_velocity!(x, icol, V_e)
        end
        if !isnothing(Ω_e)
            set_element_angular_velocity!(x, icol, Ω_e)
        end
    end

    for ipoint = 1:npoint
        icol = icol_point[ipoint]
        if !isnothing(u_p)
            set_point_deflections!(x, icol, u_p, prescribed_conditions)
        end
        if !isnothing(theta_p)
            set_point_rotations!(x, icol, theta_p, prescribed_conditions)
        end
        if !isnothing(F_p)
            set_point_forces!(x, icol, F_p, force_scaling, prescribed_conditions)
        end
        if !isnothing(M_p)
            set_point_moments!(x, icol, M_p, force_scaling, prescribed_conditions)
        end
    end

    return x
end

"""
    set_element_deflection!([x,] system, u_e, ielem)

Set the state variables in `system` (or in the vector `x`) corresponding to the
deflections of element `ielem` to the provided values.
"""
function set_element_deflection!(system::System, u_e, ielem)
    set_element_deflection!(system.x, system, u_e, ielem)
    return system
end

function set_element_deflection!(x, system::System, u_e, ielem)
    icol = system.icol_elem[ielem]
    set_element_deflection!(x, icol, u_e)
    return x
end

function set_element_deflection!(x, icol, u_e)
    x[icol:icol+2] .= u_e
    return x
end

"""
    set_element_rotation!([x,] system, θ_b, ielem)

Set the state variables in `system` (or in the vector `x`) corresponding to the
rotations of element `ielem` to the provided values.
"""
function set_element_rotation!(system::System, θ_e, ielem)
    set_element_rotation!(system.x, system, θ_e, ielem)
    return system
end

function set_element_rotation!(x, system::System, θ_e, ielem)
    icol = system.icol_elem[ielem]
    set_element_rotation!(x, icol, θ_e)
    return x
end

function set_element_rotation!(x, icol, θ_e)
    x[icol+3:icol+5] .= θ_e
    return x
end

"""
    set_element_forces!([x,] system, F_e, ielem)

Set the state variables in `system` (or in the vector `x`) corresponding to the
resultant forces of element `ielem` to the provided values.
"""
function set_element_forces!(system::System, F_e, ielem)
    set_element_forces!(system.x, system, F_e, ielem)
    return system
end

function set_element_forces!(x, system::System, F_e, ielem)
    icol = system.icol_elem[ielem]
    force_scaling = system.force_scaling
    set_element_forces!(x, icol, F_e, force_scaling)
    return x
end

function set_element_forces!(x, icol, F_e, force_scaling)
    x[icol+6:icol+8] .= F_e ./ force_scaling
    return x
end

"""
    set_element_moments!([x,] system, u_e, ielem)

Set the state variables in `system` (or in the vector `x`) corresponding to the
resultant moments of element `ielem` to the provided values.
"""
function set_element_moments!(system::System, M_e, ielem)
    set_element_moments!(system.x, system, M_e, ielem)
    return system
end

function set_element_moments!(x, system::System, M_e, ielem)
    icol = system.icol_elem[ielem]
    force_scaling = system.force_scaling
    set_element_moments!(x, icol, M_e, force_scaling)
    return x
end

function set_element_moments!(x, icol, M_e, force_scaling)
    x[icol+9:icol+11] .= M_e ./ force_scaling
    return x
end

"""
    set_element_linear_velocity!([x,] system, V_e, ielem)

Set the state variables in `system` (or in the vector `x`) corresponding to the
linear velocity of element `ielem` to the provided values.
"""
function set_element_linear_velocity!(system::System, V_e, ielem)
    set_element_linear_velocity!(x, system, V_e, ielem)
    return system
end

function set_element_linear_velocity!(x, system::System, V_e, ielem)
    @assert !system.static
    icol = system.icol_elem[ielem]
    set_element_linear_velocity!(x, icol, V_e)
    return x
end

function set_element_linear_velocity!(x, icol, V_e)
    x[icol+12:icol+14] .= V_e
    return x
end

"""
   set_element_angular_velocity!([x,] system, Ω_e, ielem)

Set the state variables in `system` (or in the vector `x`) corresponding to the
angular velocity of element `ielem` to the provided values.
"""
function set_element_angular_velocity!(system::System, Ω_e, ielem)
    set_element_angular_velocity!(x, system, Ω_e, ielem)
    return system
end

function set_element_angular_velocity!(x, system::System, Ω_e, ielem)
    @assert !system.static
    icol = system.icol_elem[ielem]
    set_element_angular_velocity!(x, icol, Ω_e)
    return x
end

function set_element_angular_velocity!(x, icol, Ω_e)
    x[icol+15:icol+17] .= Ω_e
    return x
end

"""
    set_point_deflections!([x,] system, u_e, ipoint, prescribed_conditions)

Set the state variables in `system` (or in the vector `x`) corresponding to the
deflections (or externally applied forces) at point `ipoint` to the provided values.
"""
set_point_deflections!

function set_point_deflections!(system::System, u_p, ipoint, prescribed_conditions)
    x = set_point_deflections!(system.x, system, u_p, ipoint, prescribed_conditions)
    return system
end

function set_point_deflections!(x, system::System, u_p, ipoint, prescribed_conditions)

    icol = system.icol_point[ipoint]

    if ipoint in keys(prescribed_conditions)
        prescribed_forces = prescribed_conditions[ipoint].force
    else
        prescribed_forces = @SVector ones(Bool, 6)
    end

    set_point_deflections!(x, icol, u_p, prescribed_forces)

    return x
end

function set_point_deflections!(x, icol, u_p, prescribed_forces)

    if icol <= 0
        # point variables are not state variables
        return x
    end

    # point variables are state variables
    for k = 1:3
        if prescribed_forces[k]
            # applied force is prescribed, deflection is a state variable
            x[icol+k-1] = u_p[k]
        end
    end

    return x
end

"""
    set_point_rotations!([x,] system, θ_e, ipoint, prescribed_conditions)

Set the state variables in `system` (or in the vector `x`) corresponding to the
rotations (or externally applied moments) at point `ipoint` to the provided values.
"""
set_point_rotations!

function set_point_rotations!(system::System, θ_p, ipoint, prescribed_conditions)
    set_point_rotations!(system.x, system, θ_p, ipoint, prescribed_conditions)
    return system
end

function set_point_rotations!(x, system::System, θ_p, ipoint, prescribed_conditions)

    icol = system.icol_point[ipoint]

    if ipoint in keys(prescribed_conditions)
        prescribed_forces = prescribed_conditions[ipoint].force
    else
        prescribed_forces = @SVector ones(Bool, 6)
    end

    set_point_rotations!(x, icol, θ_p, prescribed_forces)

    return x
end

function set_point_rotations!(x, icol, θ_p, prescribed_forces)

    if icol <= 0
        # point variables are not state variables
        return x
    end

    # point variables are state variables
    for k = 1:3
        if prescribed_forces[3+k]
            # applied moment is prescribed, rotation is a state variable
            x[icol+k+2] = θ_p[k]
        end
    end

    return x
end

"""
    set_point_forces!([x,] system, F_p, ipoint, prescribed_conditions)

Set the state variables in `system` (or in the vector `x`) corresponding to the
external forces applied on point `ipoint` to the provided values.
"""
set_point_forces!

function set_point_forces!(system::System, F_p, ipoint, prescribed_conditions)
    set_point_forces!(system.x, system, F_p, ipoint, prescribed_conditions)
    return system
end

function set_point_forces!(x, system::System, F_p, ipoint, prescribed_conditions)

    icol = icol_point[ipoint]

    if ipoint in keys(prescribed_conditions)
        prescribed_forces = prescribed_conditions[ipoint].force
    else
        prescribed_forces = @SVector ones(Bool, 6)
    end

    set_point_forces!(x, icol, F_p, prescribed_forces, force_scaling)

    return x
end

function set_point_forces!(x, icol, F_p, prescribed_forces, force_scaling)

    if icol <= 0
        # point variables are not state variables
        return x
    end

    # point variables are state variables
    for k = 1:3
        # prescribed conditions exist for this point
        if !prescribed_forces[k]
            # deflections are prescribed, applied force is a state variable
            x[icol+k-1] = F_p[k] / force_scaling
        end
    end

    return x
end

"""
    set_point_moments!([x,] system, M_p, ipoint, prescribed_conditions)

Set the state variables in `system` (or in the vector `x`) corresponding to the
external moments applied on point `ipoint` to the provided values.
"""
function set_point_moments!(system::System, M_p, ipoint, prescribed_conditions)
    set_point_moments!(system.x, system, M_p, ipoint, prescribed_conditions)
    return system
end

function set_point_moments!(x, system::System, M_p, ipoint, prescribed_conditions)

    icol = icol_point[ipoint]

    if ipoint in keys(prescribed_conditions)
        prescribed_forces = prescribed_conditions[ipoint].force
    else
        prescribed_forces = @SVector ones(Bool, 6)
    end

    set_point_moments!(x, icol, M_p, prescribed_forces, force_scaling)

    return x
end

function set_point_moments!(x, icol, M_p, prescribed_forces, force_scaling)

    if icol <= 0
        # point variables are not state variables
        return x
    end

    # point variables are state variables
    for k = 1:3
        # prescribed conditions exist for this point
        if !prescribed_forces[3+k]
            # rotations are prescribed, applied moment is a state variable
            x[icol+k+2] = M_p[k] / force_scaling
        end
    end

    return x
end

"""
    reset_state!(system)

Reset the state variables in `system` (stored in `system.x`).
"""
function reset_state!(system)
    system.x .= 0
    return system
end

"""
    get_sparsity(system, assembly)

Return a matrix indicating the sparsity structure of the jacobian matrix.
"""
function get_sparsity(system, assembly)

    N = length(system.x)
    irow_point = system.irow_point
    irow_elem = system.irow_elem
    irow_elem1 = system.irow_elem1
    irow_elem2 = system.irow_elem2
    icol_point = system.icol_point
    icol_elem = system.icol_elem

    sparsity = spzeros(Bool, N, N)

    # --- add constributions from beam element state variables --- #

    nelem = length(icol_elem)

    for ielem = 1:nelem
        # indices for this beam
        icol = icol_elem[ielem]
        irow_e = irow_elem[ielem]
        irow_e1 = irow_elem1[ielem]
        irow_p1 = irow_point[assembly.start[ielem]]
        irow_e2 = irow_elem2[ielem]
        irow_p2 = irow_point[assembly.stop[ielem]]

        # --- static, beam element jacobian contributions --- #

        # left endpoint equilibrium equations
        sparsity[irow_p1:irow_p1+2, icol:icol+8] .= true
        sparsity[irow_p1+3:irow_p1+5, icol+3:icol+11] .= true
        # left point compatability equations
        irow = ifelse(irow_e1 == irow_p1 || irow_e1 <= 0, irow_p1+6, irow_e1)
        sparsity[irow:irow+2, icol:icol+11] .= true
        sparsity[irow+3:irow+5, icol+3:icol+11] .= true
        # right endpoint equilibrium equations
        sparsity[irow_p2:irow_p2+2, icol:icol+8] .= true
        sparsity[irow_p2+3:irow_p2+5, icol+3:icol+11] .= true
        # right point compatability equations
        irow = ifelse(irow_e2 == irow_p2 || irow_e2 <= 0, irow_p2 + 6, irow_e2)
        sparsity[irow:irow+2, icol:icol+11] .= true
        sparsity[irow+3:irow+5, icol+3:icol+11] .= true

        # --- dynamic, beam element jacobian contributions --- #
        if !system.static
            # left endpoint equilibrium equations
            sparsity[irow_p1:irow_p1+2, icol+12:icol+14] .= true
            sparsity[irow_p1+3:irow_p1+5, icol+12:icol+17] .= true
            # right endpoint equilibrium equations
            sparsity[irow_p2:irow_p2+2, icol+12:icol+14] .= true
            sparsity[irow_p2+3:irow_p2+5, icol+12:icol+17] .= true
            # element residual equations
            sparsity[irow_e:irow_e+2, icol:icol+5] .= true
            sparsity[irow_e:irow_e+2, icol+12:icol+17] .= true
            sparsity[irow_e+3:irow_e+5, icol+3:icol+5] .= true
            sparsity[irow_e+3:irow_e+5, icol+12:icol+17] .= true
        end
    end

    # --- add constributions from point state variables --- #

    npoint = length(icol_point)

    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_point[ipoint] <= 0
            continue
        end

        icol = icol_point[ipoint]
        irow_p = irow_point[ipoint]

        for ielem = 1:nelem
            # check left side of beam
            if ipoint == assembly.start[ielem]
                # add jacobian entries for the beam endpoint
                irow_e = irow_elem1[ielem]
                if irow_e == irow_p
                    for i = 1:6
                        for j = 4:6
                            sparsity[irow_e+i-1, icol+j-1] = true
                        end
                        sparsity[irow_e+i+5, icol+i-1] = true
                        sparsity[irow_e+i-1, icol+i-1] = true
                    end
                else
                    for i = 1:6
                        sparsity[irow_e+i-1, icol+i-1] = true
                    end
                end
            end
            # check right side of beam
            if ipoint == assembly.stop[ielem]
                # add jacobian entries for the beam endpoint
                irow_e = irow_elem2[ielem]
                if irow_e == irow_p
                    for i = 1:6
                        for j = 4:6
                            sparsity[irow_e+i-1, icol+j-1] = true
                        end
                        sparsity[irow_e+i+5, icol+i-1] = true
                        sparsity[irow_e+i-1, icol+i-1] = true
                    end
                else
                    for i = 1:6
                        sparsity[irow_e+i-1, icol+i-1] = true
                    end
                end
            end
        end
    end

    return sparsity
end

"""
    static_system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads, point_masses, gvec,
        force_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

Populate the system residual vector `resid` for a static analysis

# Arguments
 - `resid`: system residual vector
 - `x`: current states of the system
 - `assembly`: assembly of nonlinear beam elements
 - `prescribed_conditions`: dictionary of prescribed conditions
 - `distributed_loads`: dictionary of distributed loads
 - `point_masses`: dictionary of point masses
 - `gvec`: gravity vector
 - `force_scaling`: scaling parameter for forces/moments
 - `irow_point`: row index of first equilibrium equation for each point
 - `irow_elem`: row index of first equation for just this beam element
 - `irow_elem1`: row index of first equation for the left side of each beam
 - `irow_elem2`: row index of first equation for the right side of each beam
 - `icol_point`: column index of first state variable for each point
 - `icol_elem`: column index of first state variable for each beam element
"""
function static_system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads, point_masses, gvec,
    force_scaling, irow_point, irow_elem1, irow_elem2, icol_point, icol_elem)

    npoint = length(assembly.points)
    nelem = length(assembly.elements)

    # add contributions to residual equations from the elements
    for ielem = 1:nelem

        # get pointers for element
        icol = icol_elem[ielem]
        irow_e1 = irow_elem1[ielem]
        irow_p1 = irow_point[assembly.start[ielem]]
        irow_e2 = irow_elem2[ielem]
        irow_p2 = irow_point[assembly.stop[ielem]]

        static_element_residual!(resid, x, ielem, assembly.elements[ielem],
            distributed_loads, point_masses, gvec, force_scaling, icol, irow_e1,
            irow_p1, irow_e2, irow_p2)
    end

    # add contributions to the residual equations from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_point[ipoint] <= 0
            continue
        end

        icol = icol_point[ipoint]
        irow_p = irow_point[ipoint]

        point_residual!(resid, x, ipoint, assembly, prescribed_conditions,
            force_scaling, icol, irow_p, irow_elem1, irow_elem2)
    end

    return resid
end

"""
    steady_state_system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads, point_masses, gvec,
        force_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
        x0, v0, ω0, a0, α0)

Populate the system residual vector `resid` for a steady state analysis.

# Arguments
 - `resid`: system residual vector
 - `x`: current states of the system
 - `assembly`: assembly of nonlinear beam elements
 - `prescribed_conditions`: dictionary of prescribed conditions
 - `distributed_loads`: dictionary of distributed loads
 - `point_masses`: dictionary of point masses
 - `gvec`: gravity vector
 - `force_scaling`: scaling parameter for forces/moments
 - `irow_point`: row index of first equilibrium equation for each point
 - `irow_elem`: row index of first equation for just this beam element
 - `irow_elem1`: row index of first equation for the left side of each beam
 - `irow_elem2`: row index of first equation for the right side of each beam
 - `icol_point`: column index of first state variable for each point
 - `icol_elem`: column index of first state variable for each beam element
 - `x0`: body frame origin
 - `v0`: body frame linear velocity
 - `ω0`: body frame angular velocity
 - `a0`: body frame linear acceleration
 - `α0`: body frame angular acceleration
"""
function steady_state_system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads, point_masses, gvec,
    force_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0, a0, α0)

    npoint = length(assembly.points)
    nelem = length(assembly.elements)

    # add contributions to residual equations from the elements
    for ielem = 1:nelem

        # get pointers for element
        icol = icol_elem[ielem]
        irow_e = irow_elem[ielem]
        irow_e1 = irow_elem1[ielem]
        irow_p1 = irow_point[assembly.start[ielem]]
        irow_e2 = irow_elem2[ielem]
        irow_p2 = irow_point[assembly.stop[ielem]]

        steady_state_element_residual!(resid, x, ielem, assembly.elements[ielem],
             distributed_loads, point_masses, gvec, force_scaling, icol, irow_e, irow_e1,
            irow_p1, irow_e2, irow_p2, x0, v0, ω0, a0, α0)
    end

    # add contributions to the residual equations from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_point[ipoint] <= 0
            continue
        end

        icol = icol_point[ipoint]
        irow_p = irow_point[ipoint]

        point_residual!(resid, x, ipoint, assembly, prescribed_conditions,
            force_scaling, icol, irow_p, irow_elem1, irow_elem2)
    end

    return resid
end

"""
    initial_condition_system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads, point_masses, gvec,
        force_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
        x0, v0, ω0, a0, α0, u, θ, udot, θdot)

Populate the system residual vector `resid` for an initial conditions analysis.

 # Arguments
 - `resid`: system residual vector
 - `x`: current states of the system
 - `assembly`: assembly of nonlinear beam elements
 - `prescribed_conditions`: dictionary of prescribed conditions
 - `distributed_loads`: dictionary of distributed loads
 - `point_masses`: dictionary of point masses
 - `gvec`: gravity vector
 - `force_scaling`: scaling parameter for forces/moments
 - `irow_point`: row index of first equilibrium equation for each point
 - `irow_elem`: row index of first equation for just this beam element
 - `irow_elem1`: row index of first equation for the left side of each beam
 - `irow_elem2`: row index of first equation for the right side of each beam
 - `icol_point`: column index of first state variable for each point
 - `icol_elem`: column index of first state variable for each beam element
 - `x0`: body frame origin
 - `v0`: body frame linear velocity
 - `ω0`: body frame angular velocity
 - `a0`: body frame linear acceleration
 - `α0`: body frame angular acceleration
 - `u`: initial linear deflections for each beam element
 - `θ`: initial angular deflections for each beam element
 - `udot`: initial linear deflection rates for each beam element
 - `θdot`: initial angular deflection rates for each beam element
"""
function initial_condition_system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads, point_masses, gvec,
    force_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0, a0, α0, u, θ, udot, θdot)

    npoint = length(assembly.points)
    nelem = length(assembly.elements)

    # add contributions to residual equations from the elements
    for ielem = 1:nelem

        # get pointers for element
        icol = icol_elem[ielem]
        irow_e = irow_elem[ielem]
        irow_e1 = irow_elem1[ielem]
        irow_p1 = irow_point[assembly.start[ielem]]
        irow_e2 = irow_elem2[ielem]
        irow_p2 = irow_point[assembly.stop[ielem]]

        initial_condition_element_residual!(resid, x, ielem, assembly.elements[ielem],
            distributed_loads, point_masses, gvec, force_scaling, icol, irow_e, irow_e1,
            irow_p1, irow_e2, irow_p2, x0, v0, ω0, a0, α0,
            u[ielem], θ[ielem], udot[ielem], θdot[ielem])
    end

    # add contributions to the residual equations from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_point[ipoint] <= 0
            continue
        end

        icol = icol_point[ipoint]
        irow_p = irow_point[ipoint]

        point_residual!(resid, x, ipoint, assembly, prescribed_conditions,
            force_scaling, icol, irow_p, irow_elem1, irow_elem2)
    end

    return resid
end

"""
    newmark_system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads, point_masses, gvec,
        force_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
        x0, v0, ω0, a0, α0, udot_init, θdot_init, Vdot_init, Ωdot_init, dt)

Populate the system residual vector `resid` for a Newmark scheme time-marching analysis.

# Arguments
 - `resid`: system residual vector
 - `x`: current states of the system
 - `assembly`: assembly of nonlinear beam elements
 - `prescribed_conditions`: dictionary of prescribed conditions
 - `distributed_loads`: dictionary of distributed loads
 - `point_masses`: dictionary of point masses
 - `gvec`: gravity vector
 - `force_scaling`: scaling parameter for forces/moments
 - `irow_point`: row index of first equilibrium equation for each point
 - `irow_elem`: row index of first equation for just this beam element
 - `irow_elem1`: row index of first equation for the left side of each beam
 - `irow_elem2`: row index of first equation for the right side of each beam
 - `icol_point`: column index of first state variable for each point
 - `icol_elem`: column index of first state variable for each beam element
 - `x0`: body frame origin
 - `v0`: body frame linear velocity
 - `ω0`: body frame angular velocity
 - `a0`: body frame linear acceleration
 - `α0`: body frame angular acceleration
 - `udot_init`: `2/dt*u + udot` for each beam element from the previous time step
 - `θdot_init`: `2/dt*θ + θdot` for each beam element from the previous time step
 - `Vdot_init`: `2/dt*V + Vdot` for each beam element from the previous time step
 - `Ωdot_init`: `2/dt*Ω + Ωdot` for each beam element from the previous time step
 - `dt`: time step size
"""
function newmark_system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads, point_masses, gvec,
    force_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0, a0, α0, udot_init, θdot_init, Vdot_init, Ωdot_init, dt)

    nelem = length(assembly.elements)
    npoint = length(assembly.points)

    # add contributions to residual equations from the beam elements
    for ielem = 1:nelem

        # get pointers for element
        icol = icol_elem[ielem]
        irow_e = irow_elem[ielem]
        irow_e1 = irow_elem1[ielem]
        irow_p1 = irow_point[assembly.start[ielem]]
        irow_e2 = irow_elem2[ielem]
        irow_p2 = irow_point[assembly.stop[ielem]]

        newmark_element_residual!(resid, x, ielem, assembly.elements[ielem],
            distributed_loads, point_masses, gvec, force_scaling, icol, irow_e, irow_e1, irow_p1, irow_e2, irow_p2,
            x0, v0, ω0, a0, α0,
            udot_init[ielem], θdot_init[ielem],
            Vdot_init[ielem], Ωdot_init[ielem], dt)
    end

    # add contributions to the residual equations from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_point[ipoint] <= 0
            continue
        end

        icol = icol_point[ipoint]
        irow_p = irow_point[ipoint]

        point_residual!(resid, x, ipoint, assembly, prescribed_conditions,
            force_scaling, icol, irow_p, irow_elem1, irow_elem2)
    end

    return resid
end

"""
    dynamic_system_residual!(resid, dx, x, assembly, prescribed_conditions, distributed_loads, point_masses, gvec,
        force_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
        x0, v0, ω0, a0, α0)

Populate the system residual vector `resid` for a general dynamic system analysis.

# Arguments
 - `resid`: system residual vector
 - `dx`: current state rates of the system
 - `x`: current states of the system
 - `assembly`: assembly of nonlinear beam elements
 - `prescribed_conditions`: dictionary of prescribed conditions
 - `distributed_loads`: dictionary of distributed loads
 - `point_masses`: dictionary of point masses
 - `gvec`: gravity vector
 - `force_scaling`: scaling parameter for forces/moments
 - `irow_point`: row index of first equilibrium equation for each point
 - `irow_elem`: row index of first equation for just this beam element
 - `irow_elem1`: row index of first equation for the left side of each beam
 - `irow_elem2`: row index of first equation for the right side of each beam
 - `icol_point`: column index of first state variable for each point
 - `icol_elem`: column index of first state variable for each beam element
 - `x0`: body frame origin
 - `v0`: body frame linear velocity
 - `ω0`: body frame angular velocity
 - `a0`: body frame linear acceleration
 - `α0`: body frame angular acceleration
"""
function dynamic_system_residual!(resid, dx, x, assembly, prescribed_conditions, distributed_loads, point_masses, gvec,
    force_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0, a0, α0)

    nelem = length(assembly.elements)
    npoint = length(assembly.points)

    # add contributions to residual equations from the beam elements
    for ielem = 1:nelem

        # get pointers for element
        icol = icol_elem[ielem]
        irow_e = irow_elem[ielem]
        irow_e1 = irow_elem1[ielem]
        irow_p1 = irow_point[assembly.start[ielem]]
        irow_e2 = irow_elem2[ielem]
        irow_p2 = irow_point[assembly.stop[ielem]]

        # set state rates for element
        udot = SVector(dx[icol], dx[icol+1], dx[icol+2])
        θdot = SVector(dx[icol+3], dx[icol+4], dx[icol+5])
        Vdot = SVector(dx[icol+12], dx[icol+13], dx[icol+14])
        Ωdot = SVector(dx[icol+15], dx[icol+16], dx[icol+17])

        dynamic_element_residual!(resid, x, ielem, assembly.elements[ielem],
             distributed_loads, point_masses, gvec, force_scaling, icol, irow_e, irow_e1, irow_p1, irow_e2, irow_p2,
             x0, v0, ω0, a0, α0, udot, θdot, Vdot, Ωdot)

    end

    # add contributions to the residual equations from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_point[ipoint] <= 0
            continue
        end

        icol = icol_point[ipoint]
        irow_p = irow_point[ipoint]

        point_residual!(resid, x, ipoint, assembly, prescribed_conditions,
            force_scaling, icol, irow_p, irow_elem1, irow_elem2)
    end

    return resid
end

"""
    static_system_jacobian!(jacob, x, assembly,
        prescribed_conditions, distributed_loads, point_masses, gvec, force_scaling,
        irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

Populate the system jacobian matrix `jacob` for a static analysis.

# Arguments
 - `jacob`: system jacobian matrix
 - `x`: current states of the system
 - `assembly`: assembly of nonlinear beam elements
 - `prescribed_conditions`: dictionary of prescribed conditions
 - `distributed_loads`: dictionary of distributed loads
 - `point_masses`: dictionary of point masses
 - `gvec`: gravity vector
 - `force_scaling`: scaling parameter for forces/moments
 - `irow_point`: row index of first equilibrium equation for each point
 - `irow_elem1`: row index of first equation for the left side of each beam
 - `irow_elem2`: row index of first equation for the right side of each beam
 - `icol_point`: column index of first state variable for each point
 - `icol_elem`: column index of first state variable for each beam element
"""
@inline function static_system_jacobian!(jacob, x, assembly,
    prescribed_conditions, distributed_loads, point_masses, gvec, force_scaling,
    irow_point, irow_elem1, irow_elem2, icol_point, icol_elem)

    jacob .= 0

    npoint = length(assembly.points)
    nelem = length(assembly.elements)

    # add contributions to residual equations from the elements
    for ielem = 1:nelem

        icol = icol_elem[ielem]
        irow_e1 = irow_elem1[ielem]
        irow_p1 = irow_point[assembly.start[ielem]]
        irow_e2 = irow_elem2[ielem]
        irow_p2 = irow_point[assembly.stop[ielem]]

        jacob = static_element_jacobian!(jacob, x, ielem, assembly.elements[ielem],
            distributed_loads, point_masses, gvec, force_scaling, icol, irow_e1,
            irow_p1, irow_e2, irow_p2)
    end

    # add contributions to the system jacobian matrix from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_point[ipoint] <= 0
            continue
        end

        icol = icol_point[ipoint]
        irow_p = irow_point[ipoint]

        jacob = point_jacobian!(jacob, x, ipoint, assembly, prescribed_conditions,
            force_scaling, icol, irow_p, irow_elem1, irow_elem2)

    end

    return jacob
end

"""
    steady_state_system_jacobian!(jacob, x, assembly,
        prescribed_conditions, distributed_loads, point_masses, gvec, force_scaling,
        irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
        x0, v0, ω0, a0, α0)

Populate the system jacobian matrix `jacob` for a steady state analysis.

# Arguments
 - `jacob`: system jacobian matrix
 - `x`: current states of the system
 - `assembly`: assembly of nonlinear beam elements
 - `prescribed_conditions`: dictionary of prescribed conditions
 - `distributed_loads`: dictionary of distributed loads
 - `point_masses`: dictionary of point masses
 - `gvec`: gravity vector
 - `force_scaling`: scaling parameter for forces/moments
 - `irow_point`: row index of first equilibrium equation for each point
 - `irow_elem`: row index of first equation for just this beam element
 - `irow_elem1`: row index of first equation for the left side of each beam
 - `irow_elem2`: row index of first equation for the right side of each beam
 - `icol_point`: column index of first state variable for each point
 - `icol_elem`: column index of first state variable for each beam element
 - `x0`: body frame origin
 - `v0`: body frame linear velocity
 - `ω0`: body frame angular velocity
 - `a0`: body frame linear acceleration
 - `α0`: body frame angular acceleration
"""
@inline function steady_state_system_jacobian!(jacob, x, assembly,
    prescribed_conditions, distributed_loads, point_masses, gvec, force_scaling,
    irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0, a0, α0)

    jacob .= 0

    npoint = length(assembly.points)
    nelem = length(assembly.elements)

    # add contributions to residual equations from the elements
    for ielem = 1:nelem

        icol = icol_elem[ielem]
        irow_e = irow_elem[ielem]
        irow_e1 = irow_elem1[ielem]
        irow_p1 = irow_point[assembly.start[ielem]]
        irow_e2 = irow_elem2[ielem]
        irow_p2 = irow_point[assembly.stop[ielem]]

        steady_state_element_jacobian!(jacob, x, ielem, assembly.elements[ielem],
            distributed_loads, point_masses, gvec, force_scaling, icol, irow_e, irow_e1,
            irow_p1, irow_e2, irow_p2, x0, v0, ω0, a0, α0)
    end

    # add contributions to the system jacobian matrix from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_point[ipoint] <= 0
            continue
        end

        icol = icol_point[ipoint]
        irow_p = irow_point[ipoint]

        point_jacobian!(jacob, x, ipoint, assembly, prescribed_conditions,
             force_scaling, icol, irow_p, irow_elem1, irow_elem2)
    end

    return jacob
end

"""
    initial_condition_system_jacobian!(jacob, x, assembly,
        prescribed_conditions, distributed_loads, point_masses, gvec, force_scaling,
        irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
        x0, v0, ω0, a0, α0, u, θ, udot, θdot)

Populate the system jacobian matrix `jacob` for an initial conditions analysis.

# Arguments
 - `jacob`: system jacobian matrix
 - `x`: current states of the system
 - `assembly`: assembly of nonlinear beam elements
 - `prescribed_conditions`: dictionary of prescribed conditions
 - `distributed_loads`: dictionary of distributed loads
 - `point_masses`: dictionary of point masses
 - `gvec`: gravity vector
 - `force_scaling`: scaling parameter for forces/moments
 - `irow_point`: row index of first equilibrium equation for each point
 - `irow_elem`: row index of first equation for just this beam element
 - `irow_elem1`: row index of first equation for the left side of each beam
 - `irow_elem2`: row index of first equation for the right side of each beam
 - `icol_point`: column index of first state variable for each point
 - `icol_elem`: column index of first state variable for each beam element
 - `x0`: body frame origin
 - `v0`: body frame linear velocity
 - `ω0`: body frame angular velocity
 - `a0`: body frame linear acceleration
 - `α0`: body frame angular acceleration
 - `u`: initial linear deflections for each beam element
 - `θ`: initial angular deflections for each beam element
 - `udot`: initial linear deflection rates for each beam element
 - `θdot`: initial angular deflection rates for each beam element
"""
@inline function initial_condition_system_jacobian!(jacob, x, assembly,
    prescribed_conditions, distributed_loads, point_masses, gvec, force_scaling,
    irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0, a0, α0, u, θ, udot, θdot)

    jacob .= 0

    npoint = length(assembly.points)
    nelem = length(assembly.elements)

    # add contributions to residual equations from the elements
    for ielem = 1:nelem

        icol = icol_elem[ielem]
        irow_e = irow_elem[ielem]
        irow_e1 = irow_elem1[ielem]
        irow_p1 = irow_point[assembly.start[ielem]]
        irow_e2 = irow_elem2[ielem]
        irow_p2 = irow_point[assembly.stop[ielem]]

        initial_condition_element_jacobian!(jacob, x, ielem, assembly.elements[ielem],
            distributed_loads, point_masses, gvec, force_scaling, icol, irow_e, irow_e1,
            irow_p1, irow_e2, irow_p2, x0, v0, ω0, a0, α0,
            u[ielem], θ[ielem], udot[ielem], θdot[ielem])
    end

    # add contributions to the system jacobian matrix from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_point[ipoint] <= 0
            continue
        end

        icol = icol_point[ipoint]
        irow_p = irow_point[ipoint]

        point_jacobian!(jacob, x, ipoint, assembly, prescribed_conditions,
            force_scaling, icol, irow_p, irow_elem1, irow_elem2)
    end

    return jacob
end

"""
    newmark_system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads, point_masses, gvec,
        force_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem, x0, v0, ω0, a0, α0, udot_init, θdot_init, Vdot_init,
        Ωdot_init, dt)

Populate the system jacobian matrix `jacob` for a Newmark scheme time-marching analysis.

# Arguments
 - `jacob`: system jacobian matrix
 - `x`: current states of the system
 - `assembly`: assembly of nonlinear beam elements
 - `prescribed_conditions`: dictionary of prescribed conditions
 - `distributed_loads`: dictionary of distributed loads
 - `point_masses`: dictionary of point masses
 - `gvec`: gravity vector
 - `force_scaling`: scaling parameter for forces/moments
 - `irow_point`: row index of first equilibrium equation for each point
 - `irow_elem`: row index of first equation for just this beam element
 - `irow_elem1`: row index of first equation for the left side of each beam
 - `irow_elem2`: row index of first equation for the right side of each beam
 - `icol_point`: column index of first state variable for each point
 - `icol_elem`: column index of first state variable for each beam element
 - `x0`: body frame origin
 - `v0`: body frame linear velocity
 - `ω0`: body frame angular velocity
 - `a0`: body frame linear acceleration
 - `α0`: body frame angular acceleration
 - `udot_init`: `2/dt*u + udot` for each beam element from the previous time step
 - `θdot_init`: `2/dt*θ + θdot` for each beam element from the previous time step
 - `Vdot_init`: `2/dt*V + Vdot` for each beam element from the previous time step
 - `Ωdot_init`: `2/dt*Ω + Ωdot` for each beam element from the previous time step
 - `dt`: time step size
"""
@inline function newmark_system_jacobian!(jacob, x, assembly,
    prescribed_conditions, distributed_loads, point_masses, gvec, force_scaling,
    irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0, a0, α0, udot_init, θdot_init, Vdot_init, Ωdot_init, dt)

    jacob .= 0

    npoint = length(assembly.points)
    nelem = length(assembly.elements)

    # add contributions to residual equations from the elements
    for ielem = 1:nelem

        icol = icol_elem[ielem]
        irow_e = irow_elem[ielem]
        irow_e1 = irow_elem1[ielem]
        irow_p1 = irow_point[assembly.start[ielem]]
        irow_e2 = irow_elem2[ielem]
        irow_p2 = irow_point[assembly.stop[ielem]]

        newmark_element_jacobian!(jacob, x, ielem, assembly.elements[ielem],
            distributed_loads, point_masses, gvec, force_scaling, icol, irow_e, irow_e1,
            irow_p1, irow_e2, irow_p2, x0, v0, ω0, a0, α0,
            udot_init[ielem], θdot_init[ielem],
            Vdot_init[ielem], Ωdot_init[ielem], dt)
    end

    # add contributions to the system jacobian matrix from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_point[ipoint] <= 0
            continue
        end

        icol = icol_point[ipoint]
        irow_p = irow_point[ipoint]

        point_jacobian!(jacob, x, ipoint, assembly, prescribed_conditions,
            force_scaling, icol, irow_p, irow_elem1, irow_elem2)
    end

    # # zero out near-zero values ( < eps() )
    # jacob = droptol!(jacob, eps(eltype(jacob)))

    return jacob
end

"""
    dynamic_system_jacobian!(jacob, dx, x, assembly,
        prescribed_conditions, distributed_loads, point_masses, gvec, force_scaling,
        irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
        x0, v0, ω0, a0, α0)

Populate the jacobian matrix `jacob` for a general dynamic analysis.

# Arguments
 - `jacob`: system jacobian matrix
 - `dx`: current state rates of the system
 - `x`: current states of the system
 - `assembly`: assembly of nonlinear beam elements
 - `prescribed_conditions`: dictionary of prescribed conditions
 - `distributed_loads`: dictionary of distributed loads
 - `point_masses`: dictionary of point masses
 - `gvec`: gravity vector
 - `force_scaling`: scaling parameter for forces/moments
 - `irow_point`: row index of first equilibrium equation for each point
 - `irow_elem`: row index of first equation for just this beam element
 - `irow_elem1`: row index of first equation for the left side of each beam
 - `irow_elem2`: row index of first equation for the right side of each beam
 - `icol_point`: column index of first state variable for each point
 - `icol_elem`: column index of first state variable for each beam element
 - `x0`: body frame origin
 - `v0`: body frame linear velocity
 - `ω0`: body frame angular velocity
 - `a0`: body frame linear acceleration
 - `α0`: body frame angular acceleration
"""
@inline function dynamic_system_jacobian!(jacob, dx, x, assembly,
    prescribed_conditions, distributed_loads, point_masses, gvec, force_scaling,
    irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0, a0, α0)

    jacob .= 0

    npoint = length(assembly.points)
    nelem = length(assembly.elements)

    # add contributions to residual equations from the elements
    for ielem = 1:nelem

        icol = icol_elem[ielem]
        irow_e = irow_elem[ielem]
        irow_e1 = irow_elem1[ielem]
        irow_p1 = irow_point[assembly.start[ielem]]
        irow_e2 = irow_elem2[ielem]
        irow_p2 = irow_point[assembly.stop[ielem]]

        # set state rates for element
        udot = SVector(dx[icol], dx[icol+1], dx[icol+2])
        θdot = SVector(dx[icol+3], dx[icol+4], dx[icol+5])
        Vdot = SVector(dx[icol+12], dx[icol+13], dx[icol+14])
        Ωdot = SVector(dx[icol+15], dx[icol+16], dx[icol+17])

        dynamic_element_jacobian!(jacob, x, ielem, assembly.elements[ielem],
            distributed_loads, point_masses, gvec, force_scaling, icol, irow_e, irow_e1,
            irow_p1, irow_e2, irow_p2, x0, v0, ω0, a0, α0,
            udot, θdot, Vdot, Ωdot)
    end

    # add contributions to the system jacobian matrix from the prescribed point conditions
    for ipoint = 1:npoint

        # skip if the unknowns have been eliminated from the system of equations
        if icol_point[ipoint] <= 0
            continue
        end

        icol = icol_point[ipoint]
        irow_p = irow_point[ipoint]

        point_jacobian!(jacob, x, ipoint, assembly, prescribed_conditions,
            force_scaling, icol, irow_p, irow_elem1, irow_elem2)
    end

    # # zero out near-zero values ( < eps() )
    # jacob = droptol!(jacob, eps(eltype(jacob)))

    return jacob
end

"""
    system_mass_matrix!(jacob, x, assembly, point_masses, force_scaling,
        irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

Populate the system mass matrix for a general dynamic analysis

# Arguments
 - `jacob`: jacobian of the residuals with respect to the state rates
 - `x`: current states of the system
 - `assembly`: assembly of nonlinear beam elements
 - `point_masses`: dictionary of point masses
 - `force_scaling`: scaling parameter for forces/moments
 - `irow_point`: row index of the first equilibrium equation for each point
 - `irow_elem`: row index of the first linear/angular velocity residual for each element
 - `irow_elem1`: row index of the first equation for the left side of each beam
 - `irow_elem2`: row index of the first equation for the right side of each beam
 - `icol_point`: column index of the first state variable for each point
 - `icol_elem`: column index of the first state variable for each beam element
"""
function system_mass_matrix!(jacob, x, assembly, point_masses, force_scaling,
    irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

    jacob .= 0

    npoint = length(assembly.points)
    nelem = length(assembly.elements)

    # add contributions to "mass matrix" from the beam elements
    for ielem = 1:nelem

        icol = icol_elem[ielem]
        irow_e = irow_elem[ielem]
        irow_e1 = irow_elem1[ielem]
        irow_p1 = irow_point[assembly.start[ielem]]
        irow_e2 = irow_elem2[ielem]
        irow_p2 = irow_point[assembly.stop[ielem]]

        element_mass_matrix!(jacob, x, ielem, assembly.elements[ielem], point_masses,
            force_scaling, icol, irow_e, irow_p1, irow_p2)
    end

    # no contributions to "mass matrix" from point state variables

    # # zero out near-zero values ( < eps() )
    # jacob = droptol!(jacob, eps(eltype(jacob)))

    return jacob
end

"""
    system_mass_matrix!(jacob, gamma, x, dx, assembly, point_masses, force_scaling,
        irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

Add the system mass matrix to `jacob`, scaled by the scaling parameter `gamma`.
"""
function system_mass_matrix!(jacob, gamma, x, assembly, point_masses, force_scaling,
    irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

    npoint = length(assembly.points)
    nelem = length(assembly.elements)

    # add contributions to "mass matrix" from the beam elements
    for ielem = 1:nelem

        icol = icol_elem[ielem]
        irow_e = irow_elem[ielem]
        irow_e1 = irow_elem1[ielem]
        irow_p1 = irow_point[assembly.start[ielem]]
        irow_e2 = irow_elem2[ielem]
        irow_p2 = irow_point[assembly.stop[ielem]]

        element_mass_matrix!(jacob, gamma, x, ielem, assembly.elements[ielem], point_masses,
            force_scaling, icol, irow_e, irow_p1, irow_p2)
    end

    # no contributions to "mass matrix" from point state variables

    return jacob
end


text:
const BATCH_SPAWN_REGEX = r"Spawning (array )?job: (?<id>[0-9a-f\-]+)(?(1) \(n=(?<n>\d+)\))"

# Test inner constructor
@testset "AWSBatchManager" begin
    @testset "constructors" begin
        @testset "inner" begin
            mgr = AWSBatchManager(
                1,
                2,
                "job-definition",
                "job-name",
                "job-queue",
                1000,
                "us-east-2",
                Minute(10),
                ip"1.0.0.0",
                ip"2.0.0.0",
            )

            # Validate that no additional fields were added without the tests being updated
            @test fieldcount(AWSBatchManager) == 10

            @test mgr.min_workers == 1
            @test mgr.max_workers == 2
            @test mgr.job_definition == "job-definition"
            @test mgr.job_name == "job-name"
            @test mgr.job_queue == "job-queue"
            @test mgr.job_memory == 1000
            @test mgr.region == "us-east-2"
            @test mgr.timeout == Minute(10)
            @test mgr.min_ip == ip"1.0.0.0"
            @test mgr.max_ip == ip"2.0.0.0"

            @test launch_timeout(mgr) == Minute(10)
            @test desired_workers(mgr) == (1, 2)
        end

        @testset "keywords" begin
            mgr = AWSBatchManager(
                3,
                4;
                definition="keyword-def",
                name="keyword-name",
                queue="keyword-queue",
                memory=1000,
                region="us-west-1",
                timeout=Second(5),
                min_ip=ip"3.0.0.0",
                max_ip=ip"4.0.0.0",
            )

            @test mgr.min_workers == 3
            @test mgr.max_workers == 4
            @test mgr.job_definition == "keyword-def"
            @test mgr.job_name == "keyword-name"
            @test mgr.job_queue == "keyword-queue"
            @test mgr.job_memory == 1000
            @test mgr.region == "us-west-1"
            @test mgr.timeout == Second(5)
            @test mgr.min_ip == ip"3.0.0.0"
            @test mgr.max_ip == ip"4.0.0.0"
        end

        @testset "defaults" begin
            mgr = AWSBatchManager(0)

            @test mgr.min_workers == 0
            @test mgr.max_workers == 0
            @test isempty(mgr.job_definition)
            @test isempty(mgr.job_name)
            @test isempty(mgr.job_queue)
            @test mgr.job_memory == -1
            @test mgr.region == "us-east-1"
            @test mgr.timeout == AWSClusterManagers.BATCH_TIMEOUT
            @test mgr.min_ip == ip"0.0.0.0"
            @test mgr.max_ip == ip"255.255.255.255"

            @test launch_timeout(mgr) == AWSClusterManagers.BATCH_TIMEOUT
            @test desired_workers(mgr) == (0, 0)
        end

        @testset "num workers" begin
            @test_throws ArgumentError AWSBatchManager(-1)
            @test_throws ArgumentError AWSBatchManager(2, 1)
            @test desired_workers(AWSBatchManager(0, 0)) == (0, 0)
            @test desired_workers(AWSBatchManager(2)) == (2, 2)
            @test desired_workers(AWSBatchManager(3:4)) == (3, 4)
            @test_throws MethodError AWSBatchManager(3:1:4)
            @test_throws MethodError AWSBatchManager(3:2:4)
        end

        @testset "queue env variable" begin
            # Mock being run on an AWS batch job
            withenv("WORKER_JOB_QUEUE" => "worker") do
                # Fall back to using the WORKER_JOB_QUEUE environmental variable
                mgr = AWSBatchManager(3)
                @test mgr.job_queue == "worker"

                # Use the queue passed in
                mgr = AWSBatchManager(3; queue="special")
                @test mgr.job_queue == "special"
            end
        end
    end

    @testset "equality" begin
        @test AWSBatchManager(3) == AWSBatchManager(3)
    end

    @testset "addprocs" begin
        mock_queue_arn = "arn:aws:batch:us-east-1:000000000000:job-queue/queue"

        # Note: due to the `addprocs` running our code with @async it can be difficult to
        # debug failures in these tests. If a failure does occur it is recommended you run
        # the code with `launch(AWSBatchManager(...), Dict(), [], Condition())` to get a
        # useful stacktrace.

        @testset "success" begin
            patches = [
                @patch AWSBatch.JobQueue(queue::AbstractString) = JobQueue(mock_queue_arn)
                @patch AWSBatch.max_vcpus(::JobQueue) = 1
                @patch function AWSBatch.run_batch(; kwargs...)
                    @async run(kwargs[:cmd])
                    return BatchJob("00000000-0000-0000-0000-000000000001")
                end
            ]

            apply(patches) do
                # Get an initial list of processes
                init_procs = procs()
                # Add a single AWSBatchManager worker
                added_procs = @test_log LOGGER "notice" BATCH_SPAWN_REGEX begin
                    addprocs(AWSBatchManager(1))
                end
                # Check that the workers are available
                @test length(added_procs) == 1
                @test procs() == vcat(init_procs, added_procs)
                # Remove the added workers
                rmprocs(added_procs; waitfor=5.0)
                # Double check that rmprocs worked
                @test init_procs == procs()
            end
        end

        @testset "worker timeout" begin
            patches = [
                @patch AWSBatch.JobQueue(queue::AbstractString) = JobQueue(mock_queue_arn)
                @patch AWSBatch.max_vcpus(::JobQueue) = 1
                @patch function AWSBatch.run_batch(; kwargs...)
                    # Avoiding spawning a worker process
                    return BatchJob("00000000-0000-0000-0000-000000000002")
                end
            ]

            @test_throws TaskFailedException apply(patches) do
                @test_log LOGGER "notice" BATCH_SPAWN_REGEX begin
                    addprocs(AWSBatchManager(1; timeout=Second(1)))
                end
            end
        end

        @testset "VCPU limit" begin
            @testset "minimum exceeds" begin
                patches = [
                    @patch function AWSBatch.JobQueue(queue::AbstractString)
                        return JobQueue(mock_queue_arn)
                    end
                    @patch AWSBatch.max_vcpus(::JobQueue) = 3
                ]

                apply(patches) do
                    @test_throws TaskFailedException addprocs(
                        AWSBatchManager(4; timeout=Second(5))
                    )
                    @test nprocs() == 1
                end
            end

            @testset "maximum exceeds" begin
                patches = [
                    @patch function AWSBatch.JobQueue(queue::AbstractString)
                        return JobQueue(mock_queue_arn)
                    end
                    @patch AWSBatch.max_vcpus(::JobQueue) = 1
                    @patch function AWSBatch.run_batch(; kwargs...)
                        for _ in 1:kwargs[:num_jobs]
                            @async run(kwargs[:cmd])
                        end
                        return BatchJob("00000000-0000-0000-0000-000000000004")
                    end
                ]
                msg = string(
                    "Due to the max VCPU limit (1) most likely only a partial amount ",
                    "of the requested workers (2) will be spawned.",
                )
                apply(patches) do
                    added_procs = @test_log LOGGER "warn" msg begin
                        addprocs(AWSBatchManager(0:2; timeout=Second(5)))
                    end
                    @test length(added_procs) > 0
                    rmprocs(added_procs; waitfor=5.0)
                end
            end
        end
    end
end

ext: jl


text:
# Matrix setup
#==========================================================================================#

"""
    influence_coefficient(horseshoe_1 :: Horseshoe, horseshoe_2 :: Horseshoe)

Compute the influence coefficient of the first `Horseshoe` at the collocation point of the second `Horseshoe`.
"""
influence_coefficient(horseshoe :: Horseshoe, horseshoe_j) = dot(velocity(collocation_point(horseshoe_j), horseshoe, 1.), horseshoe_normal(horseshoe_j))
influence_coefficient(horseshoe :: Horseshoe, horseshoe_j, V_hat) = dot(velocity(collocation_point(horseshoe_j), horseshoe, 1., V_hat), horseshoe_normal(horseshoe_j))

"""
    influence_matrix(horseshoes)
    influence_matrix(horseshoes, u_hat)

Assemble the Aerodynamic Influence Coefficient (AIC) matrix given an array of `Horseshoes` and the freestream direction ``̂u``.
"""
influence_matrix(horseshoes) = [ influence_coefficient(horsie_j, horsie_i) for horsie_i ∈ horseshoes, horsie_j ∈ horseshoes ]
influence_matrix(horseshoes, V_hat) = [ influence_coefficient(horsie_j, horsie_i, V_hat) for horsie_i ∈ horseshoes, horsie_j ∈ horseshoes ]

"""
    boundary_condition(horseshoes, U, Ω)

Assemble the boundary condition vector given an array of `Horseshoes`, the freestream velocity ``U``, and a quasi-steady rotation vector  ``Ω``.
"""
boundary_condition(horseshoes, U, Ω) = map(hs -> dot(U + Ω × collocation_point(hs), horseshoe_normal(hs)), horseshoes)

# Matrix-free setup for nonlinear analyses
#==========================================================================================#

induced_velocity(r, horseshoes, Γs, U_hat) = @views sum(x -> velocity(r, x[1], x[2], U_hat), zip(horseshoes, Γs))

# In-place version
function induced_velocity!(vel, r, horseshoes, Γs, U_hat)
    for i in eachindex(horseshoes)
        vel += @views velocity(r, horseshoes[i], Γs[i], U_hat)
    end

    vel
end

induced_trailing_velocity(r, horseshoes, Γs, U_hat) = @views sum(x -> trailing_velocity(r, x[1], x[2], U_hat), zip(horseshoes, Γs))

# In-place version
function induced_trailing_velocity!(vel, r, horseshoes, Γs, U_hat)
    for i in eachindex(horseshoes)
        vel += @views trailing_velocity(r, horseshoes[i], Γs[i], U_hat)
    end

    vel
end

function induced_velocity(r, hs, Γs, U, Ω)
    # vel = zeros(eltype(r), 3)
    induced_velocity(r, hs, Γs, -normalize(U)) - (U + Ω × r)
end

function induced_trailing_velocity(r, horseshoes, Γs, U, Ω)
    # vel = zeros(eltype(r), 3)
    induced_trailing_velocity(r, horseshoes, Γs, -normalize(U)) - (U + Ω × r)
end

# Residual computations
function residual(r, n, hs, Γs, U, Ω)
   dot(induced_velocity(r, hs, Γs, U, Ω), n)
end

solve_nonlinear(horseshoes, Γs, U_hat, Ω_hat) = map(hs -> residual(collocation_point(hs), horseshoe_normal(hs), horseshoes, Γs, U_hat, Ω_hat), horseshoes)

solve_nonlinear!(R, horseshoes, Γs, U_hat, Ω_hat) = map!(hs -> residual(collocation_point(hs), horseshoe_normal(hs), horseshoes, Γs, U_hat, Ω_hat), R, horseshoes)


text:
const DISCOURSE_DOMAIN = "https://discourse.julialang.org"

"""
    DiscoursePost

A Discourse post is unique identified by an `id` and a `slug`.
"""
struct DiscoursePost
    id::Int
    slug::String
end

"""
    url(p::DiscoursePost)

Return the URL for the specified post.
"""
url(p::DiscoursePost) = "$DISCOURSE_DOMAIN/t/$(p.slug)/$(p.id)"

struct DiscourseData
    posts::Vector{DiscoursePost}
    index::Int
end

function Base.show(io::IO, d::DiscourseData)
    print(io, "DiscourseData[index=", d.index, ", posts=", length(d.posts), "]")
end

"""
    current(r::DiscourseData)

Return the current post.
"""
current(r::DiscourseData) = r.posts[r.index]

Base.length(r::DiscourseData) = length(r.posts)

function next(r::DiscourseData)
    new_index = r.index + 1
    return DiscourseData(r.posts, new_index > length(r) ? 1 : new_index)
end

function previous(r::DiscourseData)
    new_index = r.index - 1
    return DiscourseData(r.posts, new_index == 0 ? length(r) : new_index)
end

"""
    message(r::DiscourseData)

Return a formatted string with the current post.
"""
function message(r::DiscourseData)
    link = url(current(r))
    total = length(r)
    return "$(r.index)/$(total): $link"
end

# JSON3 struct bindings
StructTypes.StructType(::Type{HoJBot.DiscourseData}) = StructTypes.Struct()
StructTypes.StructType(::Type{HoJBot.DiscoursePost}) = StructTypes.Struct()

"""
Construct a new `DiscourseData` based upon JSON3 data returned from a search.
"""
function DiscourseData(json::JSON3.Array, make_post::Function)
    list = unique(make_post(t) for t in json)
    return DiscourseData(list, 1)
end


text:
# This file includes code that was formerly a part of Julia.
# Further modifications and additions: Scott P. Jones
# License is MIT: LICENSE.md

using ModuleInterfaceTools

#ModuleInterfaceTools.debug[] = true

@api test ChrBase

for C in (ASCIIChr, LatinChr, UCS2Chr, UTF32Chr, Char)
@testset "$C" begin

    maxch = typemax(C)

    V6_COMPAT || @testset "Ranges" begin
        # This caused JuliaLang/JSON.jl#82
        rng = C('\x00'):C('\x7f')
        @test first(rng) === C('\x00')
        @test last(rng) === C('\x7f')
    end

    C != Char && @testset "Casefold character" begin
        for c = 0:UInt(maxch)
            is_valid(C, c) || continue
            ch = C(c)
            cj = Char(c)
            uj = uppercase(cj)
            if uj <= maxch
                uc = uppercase(ch)
                uc == uj || println(" $c: $maxch $uc $uj")
                @test uc == uj
            end
            @test lowercase(ch) == lowercase(cj)
        end
    end

    @testset "Edge conditions" begin
        for (val, pass) in (
            (0, true), (0xd7ff, true),
            (0xd800, false), (0xdfff, false),
            (0xe000, true), (0xffff, true),
            (0x10000, true), (0x10ffff, true),
            (0x110000, false))
            pass &= Char(val) <= maxch
            @test is_valid(C, val) == pass
        end
    end

    V6_COMPAT || @testset "Char in invalid string" for (val, pass) in (
            (String(b"\x00"), true),
            (String(b"\x7f"), true),
            (String(b"\x80"), false),
            (String(b"\xbf"), false),
            (String(b"\xc0"), false),
            (String(b"\xff"), false),
            (String(b"\xc0\x80"), false),
            (String(b"\xc1\x80"), false),
            (String(b"\xc2\x80"), Char(0x80) <= maxch),
            (String(b"\xc2\xc0"), false),
            (String(b"\xed\x9f\xbf"), Char(0xd7ff) <= maxch),
            (String(b"\xed\xa0\x80"), false),
            (String(b"\xed\xbf\xbf"), false),
            (String(b"\xee\x80\x80"), Char(0xe000) <= maxch),
            (String(b"\xef\xbf\xbf"), Char(0xffff) <= maxch),
            (String(b"\xf0\x80\x80\x80"), false),
            (String(b"\xf0\x90\x80\x80"), Char(0x10000) <= maxch),
            (String(b"\xf4\x8f\xbf\xbf"), Char(0x10ffff) <= maxch),
            (String(b"\xf4\x90\x80\x80"), false),
            (String(b"\xf5\x80\x80\x80"), false),
            (String(b"\ud800\udc00"), false),
            (String(b"\udbff\udfff"), false),
            (String(b"\ud800\u0100"), false),
            (String(b"\udc00\u0100"), false),
            (String(b"\udc00\ud800"), false))
            @test is_valid(C, val[1]) == pass
    end

    @testset "Invalid Chars" begin
        @test  is_valid(C, 'a')
        @test  is_valid(C, 'ÿ')    == ('ÿ' <= maxch)
        @test  is_valid(C, '柒')   == ('柒' <= maxch)
        @test  is_valid(C, 0xd7ff) == (Char(0xd7ff) <= maxch)
        @test  is_valid(C, 0xe000) == (Char(0xe000) <= maxch)
        @test !is_valid(C, Char(0xd800))
        @test !is_valid(C, Char(0xdfff))
    end
end
end


text:
using JLD2
@load "data/groundref_b8acd9_wynglas.jld2"
ground_ref = deepcopy(d)
varIds = collect(keys(ground_ref))
sort!(varIds, lt=DFG.natural_lt)
refx = []
refy = []
refθ = []
for idx in varIds
  push!(refx, ground_ref[idx].val[1])
  push!(refy, ground_ref[idx].val[2])
  push!(refθ, ground_ref[idx].val[3])
end

cd(@__DIR__)
using DistributedFactorGraphs
using RoME
## load FG
folder = "results/skuif_X_l7_8_9_14/mm/"
folder = "results/skuif_X_l7_8_init/mm/"
filels =  readdir(pwd()*"/"*folder)
sort!(filels, lt=natural_lt)
filels = filels[2:2:end]
n = length(filels)
big_ppi_mat_x = zeros(68,n)#Matrix{Float64}(undef, 68, 20)
big_ppi_mat_y = zeros(68,n)#Matrix{Float64}(undef, 68, 20)
big_ppi_mat_θ = zeros(68,n)#Matrix{Float64}(undef, 68, 20)
for (i,filein) in enumerate(filels)
  fg = initfg()
  @info "loading: $(filename*filein)"
  loadDFG!(fg, folder*filein)  ##
  varIds = ls(fg)
  sort!(varIds, lt=DFG.natural_lt)
    for (k, idx) in enumerate(varIds)
    # segppe = getVariableSolverData(fg, idx, :parametric).val
    # segppe = getPPE(fg, idx) |> getMeanPPE
    segppe = getPPE(fg, idx) |> getMaxPPE
    big_ppi_mat_x[k,i] = segppe[1]
    big_ppi_mat_y[k,i] = segppe[2]
    big_ppi_mat_θ[k,i] = segppe[3]
  end
end##using Plotserr_mat_x = Matrix{Float64}(undef, 68, n)
err_mat_y = Matrix{Float64}(undef, 68, n)
err_mat_θ = Matrix{Float64}(undef, 68, n)
Xrefx =  refx[9:end]
Xrefy =  refy[9:end]
Xrefθ =  refθ[9:end]
p = scatter(refx[1:8],refy[1:8])
scatter!(p, big_ppi_mat_x[1:2:16,1], big_ppi_mat_y[1:2:16,1])
scatter!(p, big_ppi_mat_x[2:2:16,1], big_ppi_mat_y[2:2:16,1])
Xbig_ppi_mat_x =  big_ppi_mat_x[17:end,:]
Xbig_ppi_mat_y =  big_ppi_mat_y[17:end,:]
Xbig_ppi_mat_θ =  big_ppi_mat_θ[17:end,:]
err_mat_x = Matrix{Float64}(undef, 52, n)
err_mat_y = Matrix{Float64}(undef, 52, n)
err_mat_θ = Matrix{Float64}(undef, 52, n)
for i = 1:52
  err_mat_x[i,:] = Xrefx[i] .- Xbig_ppi_mat_x[i,:]
  err_mat_y[i,:] = Xrefy[i] .- Xbig_ppi_mat_y[i,:]
  err_mat_θ[i,:] = Xrefθ[i] .- Xbig_ppi_mat_θ[i,:]
end
for i = 1:68
  err_mat_x[i,:] = big_ppi_mat_x[i,1] .- big_ppi_mat_x[i,:]
  err_mat_y[i,:] = big_ppi_mat_y[i,1] .- big_ppi_mat_y[i,:]
  err_mat_θ[i,:] = big_ppi_mat_θ[i,1] .- big_ppi_mat_θ[i,:]
end
Plots.plot(err_mat_x[9:end,:]')
Plots.plot(err_mat_x[1:9,:]')
Plots.plot(big_ppi_mat_x)
Plots.plot(big_ppi_mat_y)
Plots.plot(big_ppi_mat_θ)
noises = collect(noiselevels)
noises = collect(multihypos)
# mx = Plots.plot(noises, sqrt.(mean(err_mat_x[9:end,:].^2; dims=1))', labels = "Mean Error x")
mx = Plots.plot(noises, sqrt.(mean(err_mat_x[9:end,:].^2; dims=1))', labels = "mm-x")
# mx = Plots.plot(noises, mean(err_mat_x[9:end,:]; dims=1)', labels = "")
# mc = Plots.plot!(mx, noises, sqrt.(mean(err_mat_y[9:end,:].^2; dims=1))', labels="Mean Error y")
mc = Plots.plot!(mx, noises, sqrt.(mean(err_mat_y[9:end,:].^2; dims=1))', labels="par-x")
xlabel!("Offset on landmarks l7 and l8 [m]")
ylabel!("RMS perturbation error over all poses [m]")
ylims!(0,0.15)
xlims!(0,0.5)
savefig("plots/skuif_X_l7_8_init/mm.pdf")
Plots.plot!(mc, )
Plots.plot(err_mat_θ)
ground_ref[:l10].val


text:
function LinearOptics(physic::Union{Scalar, Vectorial}, atoms, laser)
    @debug "start: LinearOptics - $( typeof(physic) )"

    kernelFunction = get_kernelFunction(physic, atoms)
    spectrum = Dict(:λ=>zeros(ComplexF64, atoms.N), :ψ=>zeros(ComplexF64, atoms.N, atoms.N))
    data = Dict()

    @debug "end  : LinearOptics - $( typeof(physic) )"
    return LinearOptics(physic, atoms, laser, kernelFunction, spectrum, data)
end

# get_kernelFunction(physics::Scalar,    dimension::Atom{<:TwoD}) = sin
# get_kernelFunction(physics::Vectorial, dimension::Atom{<:TwoD}) = cos

get_kernelFunction(physics::Scalar,    dimension::Atom{<:ThreeD}) = green_scalar!
get_kernelFunction(physics::Vectorial, dimension::Atom{<:ThreeD}) = green_vectorial!


text:
function Triangle{T}(domain::Triangle{T}, parent::ArbField) where {T<:Union{arb,fmpq}}
    Triangle(domain.angles[1], domain.angles[2], parent)
end

function Base.show(io::IO, domain::Triangle{fmpq})
    angles_string = ["$(numerator(a))π/$(denominator(a))" for a in domain.angles]
    print(
        io,
        "Triangle with $(domain.parent.prec) bits of precision and angles ($(angles_string[1]), $(angles_string[2]), $(angles_string[3]))",
    )
end

function Base.show(io::IO, domain::Triangle{arb})
    angles_string = ["$a" for a in domain.angles]
    print(
        io,
        "Triangle with $(domain.parent.prec) bits of precision and angles ($(angles_string[1]), $(angles_string[2]), $(angles_string[3]))",
    )
end

boundaries(::Triangle) = 1:3

"""
    angle(domain::Triangle, i::Integer)

Return the angle for vertex `i` of the triangle.
"""
function angle(domain::Triangle{fmpq}, i::Integer)
    if i ∈ boundaries(domain)
        return domain.parent(π) * domain.angles[i]
    else
        throw(ArgumentError("attempt to get vertex $i from a $(typeof(domain))"))
    end
end

function angle(domain::Triangle{arb}, i::Integer)
    if i ∈ boundaries(domain)
        return domain.angles[i]
    else
        throw(ArgumentError("attempt to get vertex $i from a $(typeof(domain))"))
    end
end

"""
    angledivπ(domain::Triangle, i::Integer)

Return the angle for vertex `i` of the triangle divided by `π`.
"""
function angledivπ(domain::Triangle{fmpq}, i::Integer)
    if i ∈ boundaries(domain)
        return domain.angles[i]
    else
        throw(ArgumentError("attempt to get vertex $i from a $(typeof(domain))"))
    end
end

function angledivπ(domain::Triangle{arb}, i::Integer)
    if i ∈ boundaries(domain)
        return domain.angles[i] / domain.parent(π)
    else
        throw(ArgumentError("attempt to get vertex $i from a $(typeof(domain))"))
    end
end

"""
    angles(domain::Triangle)

Return the angles of the triangle.
"""
angles(domain::Triangle) = tuple((angle(domain, i) for i in boundaries(domain))...)

"""
    anglesdivπ(domain::Triangle)

Return the angles of the triangle divided by `π`.
"""
anglesdivπ(domain::Triangle) = tuple((angledivπ(domain, i) for i in boundaries(domain))...)

"""
    vertex(domain::Triangle, i::Integer)

Return the Cartesian coordinates of vertex `i` of the triangle.
"""
function vertex(domain::Triangle{T}, i::Integer) where {T}
    if i == 1
        return SVector(domain.parent(0), domain.parent(0))
    elseif i == 2
        return SVector(domain.parent(1), domain.parent(0))
    elseif i == 3
        if T == arb
            x = sin(angle(domain, 2)) / sin(angle(domain, 3))
            s, c = sincos(angle(domain, 1))
        else
            x =
                sinpi(angledivπ(domain, 2), domain.parent) /
                sinpi(angledivπ(domain, 3), domain.parent)
            s, c = sincospi(angledivπ(domain, 1), domain.parent)
        end
        return SVector(c * x, s * x)
    else
        throw(ArgumentError("attempt to get vertex $i from a $(typeof(domain))"))
    end
end

"""
    vertices(domain::Triangle)

Return all three vertices of the triangle.
"""
vertices(domain::Triangle) = tuple((vertex(domain, i) for i in boundaries(domain))...)

area(domain::Triangle) = vertex(domain, 3)[2] / 2

"""
    center(domain::Triangle)

Return the centroid of the triangle.
"""
center(domain::Triangle) = sum(vertices(domain)) / 3

function Base.in(xy, domain::Triangle)
    x, y = xy
    return (
        y >= 0 &&
        atan(y, x) <= angle(domain, 1) &&
        domain.parent(π) - atan(y, x - 1) <= angle(domain, 2) # To the left of the right edge
    )
end

function boundary_parameterization(t, domain::Triangle, i::Integer)
    i ∈ boundaries(domain) ||
        throw(ArgumentError("attempt to get vertex $i from a $(typeof(domain))"))
    v = vertex(domain, mod1(i + 1, 3))
    w = vertex(domain, mod1(i + 2, 3))
    return v .+ t .* (w - v)
end

function boundary_points(
    domain::Triangle,
    i::Integer,
    n::Integer;
    distribution = :chebyshev,
)
    i ∈ boundaries(domain) ||
        throw(ArgumentError("attempt to get vertex $i from a $(typeof(domain))"))

    v = vertex(domain, mod1(i + 1, 3))
    w = vertex(domain, mod1(i + 2, 3))

    if distribution == :linear
        f = j -> domain.parent(j // (n + 1))
    elseif distribution == :chebyshev
        f = j -> (cospi(domain.parent((2j - 1) // 2n)) + 1) / 2
    elseif distribution == :exponential
        m = (n + 1) // 2
        f = j -> begin
            d = exp(-domain.parent(m - j))
            if j < m
                return d / 2
            elseif j > m
                return 1 - inv(d) / 2
            else
                return domain.parent(0.5)
            end
        end
    elseif distribution == :root_exponential
        m = (n + 1) // 2
        f = j -> begin
            d = exp(-domain.parent(m - j) / sqrt(domain.parent(n)))
            if j < m
                return d / 2
            elseif j > m
                return 1 - inv(d) / 2
            else
                return domain.parent(0.5)
            end
        end
    else
        throw(ArgumentError("no distribution named $distribution"))
    end

    points = Vector{SVector{2,arb}}(undef, n)
    for j = 1:n
        t = f(j)
        points[j] = v .+ t .* (w - v)
    end

    return points, fill(i, n)
end

function interior_points(domain::Triangle, n::Integer; rng = MersenneTwister(42))
    v = vertex(domain, 2)
    w = vertex(domain, 3)

    points = Vector{SVector{2,arb}}(undef, n)
    for i = 1:n
        u1, u2 = rand(rng), rand(rng)
        if u1 + u2 > 1
            u1 = 1 - u1
            u2 = 1 - u2
        end
        points[i] = u1 .* v + u2 .* w
    end

    return points
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule lazygit_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("lazygit")
JLLWrappers.@generate_main_file("lazygit", UUID("01a6f0b2-7079-58b6-a815-2223d49c44bc"))
end  # module lazygit_jll


text:
struct FourCornersVariant <: AbstractPolicy
    high_prob::Float64
    low_prob::Float64
    state_to_weights_map::Dict{Array{Int64, 1}, Weights}
    function FourCornersVariant(high_prob, low_prob)
        map = Dict{Array{Int64, 1}, Weights}(
            [3,2]=>Weights([high_prob, low_prob, low_prob, high_prob]),
            [3,9]=>Weights([high_prob, high_prob, low_prob, low_prob]),
            [10,2]=>Weights([low_prob, low_prob, high_prob, high_prob]),
            [10,9]=>Weights([low_prob, high_prob, high_prob, low_prob])
        )
        new(high_prob, low_prob, map)
    end
end

FourCornersVariant() = FourCornersVariant(0.45, 0.05)

function Base.get(π::FourCornersVariant, state_t, action_t, state_tp1, action_tp1, preds_tp1)
    prob = -1.0
    if state ∈ keys(π.state_to_weights_map)
        prob = π.state_to_weights_map[state][action_t]
    else
        prob = 0.25
    end
    return prob
end

function StatsBase.sample(rng::Random.AbstractRNG, π::FourCornersVariant, state)
    act = -1
    if state ∈ keys(π.state_to_weights_map)
        act = StatsBase.sample(rng, 1:4, π.state_to_weights_map[state])
    else
        act = rand(rng, 1:4)
    end

    return act
end

function find_same(states)
    s=0
    for p in collect(Iterators.product(states,states))
        s += if (p[1] == p[2]) 1 else 0 end
    end
    return (s - length(states))/2
end

function gen_states(seed, num)
    mt = Random.MersenneTwister(seed)
    states = []
    while length(states) < num
        new_state = rand(mt, 1:11, 2)
        if !(new_state in states)
            append!(states, [new_state])
        end
    end
    return states
end

struct RandomStateVariant <: AbstractPolicy
    states::Array{Array{Int64, 1}, 1}
    weights::Weights
end

RandomStateVariant(action, prob, num_states, seed) =
    RandomStateVariant(gen_states(seed, num_states),
                       Weights([
                           act == action ? prob : (1.0-prob)/3
                           for act in 1:4
                       ]))

RandomStateVariant() =
    RandomStateVariant(FourRoomsParams.DOWN, 0.05, 25, 1143139448)

function Base.get(π::RandomStateVariant, state_t::Array{Int64, 1}, action_t, state_tp1, action_tp1, preds_tp1)
    prob = -1.0
    if state_t ∈ π.states
        prob = π.weights[action_t]
    else
        prob = 0.25
    end
    return prob
end

function Base.get(π::RandomStateVariant, state_t::CartesianIndex{2}, action_t, state_tp1, action_tp1, preds_tp1)
    return StatsBase.get(π, [state_t[1], state_t[2]], action_t, state_tp1, action_tp1, preds_tp1)
end

function StatsBase.sample(rng::Random.AbstractRNG, π::RandomStateVariant, state::Array{Int64, 1})
    act = -1
    if state ∈ π.states
        act = StatsBase.sample(rng, 1:4, π.weights)
    else
        act = rand(rng, 1:4)
    end
    return act
end

function StatsBase.sample(rng::Random.AbstractRNG, π::RandomStateVariant, state::CartesianIndex{2})
    return StatsBase.sample(rng, π, [state[1], state[2]])
end

# # Changing behaviour policy
# mutable struct changing_behavior
#     freq
#     step
# end

# function next_action_cp(CB::changing_behavior, state_t; rng=Random.GLOBAL_RNG)

# end

# function probability_cp(CB::changing_behavior, state_t, action_t; rng=Random.GLOBAL_RNG)

# end

# function create_changing_behavior_policy(freq)
#     changing_behavior(freq, 0)

# end


text:
module SparseWrappers

using LinearAlgebra
using SparseArrays

import LinearAlgebra:   mul!
import SparseArrays:    SparseMatrixCSCUnion
import Base.Order: Forward
import LinearAlgebra: AbstractTriangular, UnitLowerTriangular, UnitUpperTriangular

import SparseArrays: nzrange, rowvals, nonzeros, nnz, nzvalview, SparseMatrixCSCView

if VERSION >= v"1.1.0-DEV"
import SparseArrays: LowerTriangularPlain, UpperTriangularPlain
else
const UnitDiagonalTriangular = Union{UnitUpperTriangular,UnitLowerTriangular}

const LowerTriangularPlain{T} = Union{
            LowerTriangular{T,<:SparseMatrixCSCUnion{T}},
            UnitLowerTriangular{T,<:SparseMatrixCSCUnion{T}}}

const LowerTriangularWrapped{T} = Union{
            Adjoint{T,<:UpperTriangular{T,<:SparseMatrixCSCUnion{T}}},
            Adjoint{T,<:UnitUpperTriangular{T,<:SparseMatrixCSCUnion{T}}},
            Transpose{T,<:UpperTriangular{T,<:SparseMatrixCSCUnion{T}}},
            Transpose{T,<:UnitUpperTriangular{T,<:SparseMatrixCSCUnion{T}}}} where T

const UpperTriangularPlain{T} = Union{
            UpperTriangular{T,<:SparseMatrixCSCUnion{T}},
            UnitUpperTriangular{T,<:SparseMatrixCSCUnion{T}}}

const UpperTriangularWrapped{T} = Union{
            Adjoint{T,<:LowerTriangular{T,<:SparseMatrixCSCUnion{T}}},
            Adjoint{T,<:UnitLowerTriangular{T,<:SparseMatrixCSCUnion{T}}},
            Transpose{T,<:LowerTriangular{T,<:SparseMatrixCSCUnion{T}}},
            Transpose{T,<:UnitLowerTriangular{T,<:SparseMatrixCSCUnion{T}}}} where T

const UpperTriangularSparse{T} = Union{
            UpperTriangularWrapped{T}, UpperTriangularPlain{T}} where T

const LowerTriangularSparse{T} = Union{
            LowerTriangularWrapped{T}, LowerTriangularPlain{T}} where T

const TriangularSparse{T} = Union{
            LowerTriangularSparse{T}, UpperTriangularSparse{T}} where T

end

"""
    Conjugate(A::AbstractMatrix)

Represent the elementwise conjugate of a given matrix.
"""
struct Conjugate{T,S} <:AbstractMatrix{T}
    parent::S
    Conjugate(A::AbstractMatrix{T}) where T = new{T,typeof(A)}(A)
end
Base.size(A::Conjugate) = size(A.parent)
Base.getindex(A::Conjugate,ind...) = conj.(getindex(A.parent, ind...))
Base.parent(A::Conjugate) = A.parent

"""
    HermiteTridiagonal(diagonal::Vector,upper::Vector)

Represent a Hermitian tridiagonal matrix.
"""
struct HermiteTridiagonal{T,V<:AbstractVector{T}} <: AbstractMatrix{T}
    dv::V
    ev::V
end
Base.size(A::HermiteTridiagonal) = (length(dv), length(dv))

struct XSubArray{T,Ti,N,S<:SparseMatrixCSC{T,Ti},I,B,R<:AbstractVector{Ti}} <: AbstractSparseArray{T,Ti,N}
    sub::SubArray{T,N,S,I,B}
    rowval::R
end
function XSubArray(s::SubArray{T,N,SparseMatrixCSC{T,Ti},I,B}) where {T,N,Ti,I,B}
    rv = rowvals(s)
    XSubArray(s, rv)
end

const SparseUnion{Tv,Ti} = Union{SparseMatrixCSC{Tv,Ti},
                                 SubArray{Tv,2,SparseMatrixCSC{Tv,Ti},<:Tuple,false},
                                 XSubArray{Tv,Ti,2,SparseMatrixCSC{Tv,Ti},<:Tuple,false}}

struct ArrayViewWrapper{T,S<:AbstractMatrix,W} <: AbstractMatrix{T}
    data::W
end
ArrayViewWrapper(X::ArrayViewWrapper) = X
function ArrayViewWrapper(X::S) where {T,S<:AbstractArray{T}}
    ArrayViewWrapper{T,array_storage(S),S}(X)
end
Base.parent(X::ArrayViewWrapper) = X.data
Base.size(X::ArrayViewWrapper) = size(X.data)
Base.getindex(X::ArrayViewWrapper, I...) = getindex(X.data, I...)

"""
    SparseMatrixCSCSymmHerm

`Symmetric` or `Hermitian` of a `SparseMatrixCSC` or `SparseMatrixCSCView`.
"""
const SparseMatrixCSCSymmHerm{Tv,Ti} = Union{Symmetric{Tv,<:SparseUnion{Tv,Ti}},
                                            Hermitian{Tv,<:SparseUnion{Tv,Ti}}}

export iswrsparse, indextype, isupper, inflate, depth, sparsecsc, sparseaccess, sparsecopy
export unwrap
export array_storage
export Conjugate, HermiteTridigonal

up(B::Union{Symmetric,Hermitian,Bidiagonal}) = B.uplo == 'U' ? :U : :L
toggle(s::Symbol) = s == :U ? :L : :U

#export nziterator
#export Conjugate, HermitianTriangular

#include("nziterators.jl")
#include("operations.jl")
#include("sparsify.jl")
#include("combine.jl")
#include("semigroup.jl")
#include("linalg.jl")
#include("fallback.jl")
#include("views.jl")
#include("timings.jl")
#include("mmatrix.jl")
#include("universal.jl")
include("wrappers.jl")

end # module