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
 - `mass_scaling`: Scaling for state variables corresponding to masses/inertias
 - `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`: `2/dt*u + udot` for each beam element from the previous time step
 - `θdot`: `2/dt*θ + θdot` for each beam element from the previous time step
 - `Pdot`: `2/dt*C'*Cab*P + C'*Cab*Pdot` for each beam element from the previous time step
 - `Hdot`: `2/dt*C'*Cab*H + C'*Cab*Hdot` for each beam element from the previous time step
 - `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
    mass_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}}
    Pdot::Vector{SVector{3,TF}}
    Hdot::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.
 - `mass_scaling`: Factor used to scale system mass/inertia internally.  If not
    specified, a suitable default will be chosen based on the entries of the
    mass 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),
    mass_scaling = default_mass_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]
    Pdot = [@SVector zeros(TF, 3) for i = 1:nelem]
    Hdot = [@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, mass_scaling,
        irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem, udot, θdot,
        Pdot, Hdot, t)
end

function default_force_scaling(assembly)

    TF = eltype(assembly)

    # create vector of all compliance entries
    compliance_entries = vcat([vcat(elem.C11..., elem.C12..., elem.C22...) 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(nonzero_compliance_entries)/100)
    end

    return force_scaling
end

function default_mass_scaling(assembly)

    TF = eltype(assembly)

    # create vector of all inverse mass matrix entries
    minv_entries = vcat([vcat(elem.minv11..., elem.minv12..., elem.minv22...) for elem in assembly.elements]...)

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

    # set mass scaling based on nonzero inverse mass matrix entries
    if isempty(minv_nonzero_indices)
        mass_scaling = 1.0
    else
        nonzero_minv_entries = minv_entries[minv_nonzero_indices]
        mass_scaling = nextpow(2.0, length(nonzero_minv_entries)/sum(nonzero_minv_entries))
    end

    return mass_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 global frame) for each beam element.
 - `theta_e`: Vector containing rotation variables (in the global frame) for each beam element.
 - `F_e`: Vector containing resultant forces (in the beam coordinate frame) for each beam element.
 - `M_e`: Vector containing resultant moments (in the beam coordinate frame) for each beam element.
 - `P_e`: Vector containing linear momenta (in the beam coordinate frame) for each beam element.
 - `H_e`: Vector containing angular momenta (in the beam coordinate frame) for each beam element.
 - `u_p`: Vector containing initial displacements (in the global frame) for each point.
 - `theta_p`: Vector containing rotation variables (in the global frame) for each point.
 - `F_p`: Vector containing externally applied forces (in the global frame) for each point.
 - `M_p`: Vector containing externally applied moments (in the global 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, mass_scaling,
    prescribed_conditions; u_e = nothing, theta_e = nothing,
    F_e = nothing, M_e = nothing, P_e = nothing, H_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_deflections!(x, icol, u_e)
        end
        if !isnothing(theta_e)
            set_element_rotations!(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(P_e)
            set_element_linear_momenta!(x, icol, P_e, mass_scaling)
        end
        if !isnothing(H_e)
            set_element_angular_momenta!(x, icol, H_e, mass_scaling)
        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_deflections!([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_deflections!(system::System, u_e, ielem)
    set_element_deflections!(system.x, system, u_e, ielem)
    return system
end

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

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

"""
    set_element_rotations!([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_rotations!(system::System, θ_e, ielem)
    set_element_rotations!(system.x, system, θ_e, ielem)
    return system
end

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

function set_element_rotations!(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_momenta!([x,] system, u_e, ielem)

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

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

function set_element_linear_momenta!(x, icol, P_e, mass_scaling)
    x[icol+12:icol+14] .= P_e ./ mass_scaling
    return x
end

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

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

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

function set_element_angular_momenta!(x, icol, H_e, mass_scaling)
    x[icol+15:icol+17] .= H_e ./ mass_scaling
    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

"""
    system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem)
    system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem, x0, v0, ω0)
    system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem, x0, v0, ω0, u, θ, udot, θdot)
    system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem, x0, v0, ω0, udot_init, θdot_init, CtCabPdot_init,
        CtCabHdot_init, dt)

Populate the residual vector `resid` with the results of the residual equations
for the system.

There are four implementations corresponding to the following analysis types:
 - Static
 - Dynamic - Steady State
 - Dynamic - Initial Step (for initializing time domain simulations)
 - Dynamic - Time Marching

 See "GEBT: A general-purpose nonlinear analysis tool for composite beams" by
 Wenbin Yu and "Geometrically nonlinear analysis of composite beams using
 Wiener-Milenković parameters" by Qi Wang and Wenbin Yu.

# Arguments
 - `resid`: System residual vector
 - `x`: Current state variables of the system
 - `assembly`: Assembly of rigidly connected nonlinear beam elements
 - `prescribed_conditions`: Dictionary of prescribed conditions at all time steps
 - `distributed_loads`: Dictionary of distributed loads at all time steps
 - `force_scaling`: Scaling parameter for forces/moments
 - `mass_scaling`: Scaling parameter for masses/inertias
 - `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

# Additional Arguments for Dynamic Analyses
 - `x0`: Global frame origin (for the current time step)
 - `v0`: Global frame linear velocity (for the current time step)
 - `ω0`: Global frame angular velocity (for the current time step)

# Additional Arguments for Initial Step Analyses
 - `u`: deflection variables for each beam element
 - `θ`: rotation variables for each beam element
 - `udot`: time derivative of u for each beam element
 - `θdot`: time derivative of θ for each beam element

# Additional Arguments for Time Marching Analyses
 - `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
 - `CtCabPdot_init`: `-2/dt*C'*Cab*P - C'*Cab*Pdot` for each beam element from the previous time step
 - `CtCabHdot_init`: `-2/dt*C'*Cab*H - C'*Cab*Hdot` for each beam element from the previous time step
 - `dt`: Time step size. If set to `nothing`, `udot_init`,
"""
system_residual!

# static
function system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
    force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

    return static_system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)
end

# dynamic - steady state
function system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
    force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0)

    return steady_state_system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
        x0, v0, ω0)
end

# dynamic - initial step
function system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
    force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0, u, θ, udot, θdot)

    return initial_step_system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
    force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0, u, θ, udot, θdot)
end

# dynamic - newmark scheme time marching
function system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
    force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0, udot_init, θdot_init, CtCabPdot_init, CtCabHdot_init, dt)

    # dynamic - newmark scheme time marching
    return newmark_system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
        x0, v0, ω0, udot_init, θdot_init, CtCabPdot_init, CtCabHdot_init, dt)
end

# static
function static_system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
    force_scaling, mass_scaling, irow_point, irow_elem, 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_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]]

        resid = static_element_residual!(resid, x, ielem, assembly.elements[ielem],
             distributed_loads, force_scaling, mass_scaling, icol, irow_e, 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]

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

    return resid
end

# dynamic - steady state
function steady_state_system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
    force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω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]]

        resid = steady_state_element_residual!(resid, x, ielem, assembly.elements[ielem],
             distributed_loads, force_scaling, mass_scaling, icol, irow_e, irow_e1,
            irow_p1, irow_e2, irow_p2, x0, v0, ω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]

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

    return resid
end

# dynamic - initial step
function initial_step_system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
    force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω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]]

        resid = initial_step_element_residual!(resid, x, ielem, assembly.elements[ielem],
             distributed_loads, force_scaling, mass_scaling, icol, irow_e, irow_e1,
            irow_p1, irow_e2, irow_p2, x0, v0, ω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]

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

    return resid
end

# dynamic - newmark scheme time marching
function newmark_system_residual!(resid, x, assembly, prescribed_conditions, distributed_loads,
    force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0, udot_init, θdot_init, CtCabPdot_init, CtCabHdot_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]]

        resid = newmark_element_residual!(resid, x, ielem, assembly.elements[ielem],
             distributed_loads, force_scaling, mass_scaling, icol, irow_e, irow_e1, irow_p1, irow_e2, irow_p2,
             x0, v0, ω0,
             udot_init[ielem], θdot_init[ielem],
             CtCabPdot_init[ielem], CtCabHdot_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]

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

    return resid
end

# dynamic - general
function dynamic_system_residual!(resid, x, dx, assembly, prescribed_conditions, distributed_loads,
    force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω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])
        Pdot = SVector(dx[icol+12], dx[icol+13], dx[icol+14]) .* mass_scaling
        Hdot = SVector(dx[icol+15], dx[icol+16], dx[icol+17]) .* mass_scaling

        resid = dynamic_element_residual!(resid, x, ielem, assembly.elements[ielem],
             distributed_loads, force_scaling, mass_scaling, icol, irow_e, irow_e1, irow_p1, irow_e2, irow_p2,
             x0, v0, ω0, udot, θdot, Pdot, Hdot)

    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]

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

    return resid
end

"""
    system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem)
    system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem, x0, v0, ω0)
    system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem, x0, v0, ω0, u, θ, udot, θdot)
    system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2,
        icol_point, icol_elem, x0, v0, ω0, udot_init, θdot_init, CtCabPdot_init,
        CtCabHdot_init, dt)

Populate the jacobian matrix `jacob` with the jacobian of the residual vector
with respect to the state variables.

There are four implementations corresponding to the following analysis types:
 - Static
 - Dynamic - Steady State
 - Dynamic - Initial Step (for initializing time domain simulations)
 - Dynamic - Time Marching

 See "GEBT: A general-purpose nonlinear analysis tool for composite beams" by
 Wenbin Yu and "Geometrically nonlinear analysis of composite beams using
 Wiener-Milenković parameters" by Qi Wang and Wenbin Yu.

# Arguments
 - `jacob`: Jacobian matrix
 - `x`: Vector containing current state variables of the system
 - `assembly`: Assembly of rigidly connected nonlinear beam elements
 - `prescribed_conditions`: Dictionary of prescribed conditions at all time steps
 - `distributed_loads`: Dictionary of distributed loads at all time steps
 - `force_scaling`: Scaling parameter for forces/moments
 - `mass_scaling`: Scaling parameter for masses/inertia
 - `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

# Additional Arguments for Dynamic Analyses
 - `x0`: Global frame origin (for the current time step)
 - `v0`: Global frame linear velocity (for the current time step)
 - `ω0`: Global frame angular velocity (for the current time step)

# Additional Arguments for Initial Step Analyses
 - `u`: deflection variables for each beam element
 - `θ`: rotation variables for each beam element
 - `udot`: time derivative of u for each beam element
 - `θdot`: time derivative of θ for each beam element

# Additional Arguments for Time Marching Analyses
 - `udot`: `-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
 - `CtCabPdot`: `-2/dt*C'*Cab*P - C'*Cab*Pdot` for each beam element from the previous time step
 - `CtCabHdot`: `-2/dt*C'*Cab*H - C'*Cab*Hdot` for each beam element from the previous time step
 - `dt`: time step size
"""
system_jacobian!

# static
@inline function system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
    force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

    return static_system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)
end

# dynamic - steady state
@inline function system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
    force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0)

    return steady_state_system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
        x0, v0, ω0)
end

# dynamic - initial step
@inline function system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
    force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0, u, θ, udot, θdot)

    return initial_step_system_jacobian!(jacob, x, assembly, prescribed_conditions, distributed_loads,
        force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
        x0, v0, ω0, u, θ, udot, θdot)
end

# dynamic - newmark scheme time marching
@inline function system_jacobian!(jacob, x, assembly, prescribed_conditions,
    distributed_loads, force_scaling, mass_scaling,
    irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0, udot_init, θdot_init, CtCabPdot_init, CtCabHdot_init, dt)

    return newmark_system_jacobian!(jacob, x, assembly, prescribed_conditions,
        distributed_loads, force_scaling, mass_scaling,
        irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
        x0, v0, ω0, udot_init, θdot_init, CtCabPdot_init, CtCabHdot_init, dt)
end

# static
@inline function static_system_jacobian!(jacob, x, assembly,
    prescribed_conditions, distributed_loads, force_scaling, mass_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 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]]

        jacob = element_jacobian!(jacob, x, ielem, assembly.elements[ielem],
            distributed_loads, force_scaling, mass_scaling, icol, irow_e, 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

# dynamic - steady state
@inline function steady_state_system_jacobian!(jacob, x, assembly,
    prescribed_conditions, distributed_loads, force_scaling, mass_scaling,
    irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω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]]

        jacob = element_jacobian!(jacob, x, ielem, assembly.elements[ielem],
            distributed_loads, force_scaling, mass_scaling, icol, irow_e, irow_e1,
            irow_p1, irow_e2, irow_p2, x0, v0, ω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]

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

    return jacob
end

# dynamic - initial step
@inline function initial_step_system_jacobian!(jacob, x, assembly,
    prescribed_conditions, distributed_loads, force_scaling, mass_scaling,
    irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω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]]

        jacob = element_jacobian!(jacob, x, ielem, assembly.elements[ielem],
            distributed_loads, force_scaling, mass_scaling, icol, irow_e, irow_e1,
            irow_p1, irow_e2, irow_p2, x0, v0, ω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]

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

    return jacob
end

# dynamic - newmark scheme time marching
@inline function newmark_system_jacobian!(jacob, x, assembly,
    prescribed_conditions, distributed_loads, force_scaling, mass_scaling,
    irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω0, udot_init, θdot_init, CtCabPdot_init, CtCabHdot_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]]

        jacob = element_jacobian!(jacob, x, ielem, assembly.elements[ielem],
            distributed_loads, force_scaling, mass_scaling, icol, irow_e, irow_e1,
            irow_p1, irow_e2, irow_p2, x0, v0, ω0,
            udot_init[ielem], θdot_init[ielem],
            CtCabPdot_init[ielem], CtCabHdot_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]

        jacob = 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 - general
@inline function dynamic_system_jacobian!(jacob, x, dx, assembly,
    prescribed_conditions, distributed_loads, force_scaling, mass_scaling,
    irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
    x0, v0, ω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])
        Pdot = SVector(dx[icol+12], dx[icol+13], dx[icol+14]) .* mass_scaling
        Hdot = SVector(dx[icol+15], dx[icol+16], dx[icol+17]) .* mass_scaling

        jacob = dynamic_element_jacobian!(jacob, x, ielem, assembly.elements[ielem],
            distributed_loads, force_scaling, mass_scaling, icol, irow_e, irow_e1,
            irow_p1, irow_e2, irow_p2, x0, v0, ω0,
            udot, θdot, Pdot, Hdot)
    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

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

    return jacob
end

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

Populate the system "mass matrix", the jacobian of the residual vector
with respect to the time derivatives of the state variables.

See "GEBT: A general-purpose nonlinear analysis tool for composite beams" by
Wenbin Yu and "Geometrically nonlinear analysis of composite beams using
Wiener-Milenković parameters" by Qi Wang and Wenbin Yu.

# Arguments
 - `jacob`: Jacobian matrix
 - `x`: Vector containing current state variables of the system
 - `assembly`: Assembly of rigidly connected nonlinear beam elements
 - `force_scaling`: Scaling parameter for forces
 - `mass_scaling`: Scaling parameter for masses
 - `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 system_mass_matrix!(jacob, x, assembly, force_scaling, mass_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, assembly.elements[ielem], force_scaling,
            mass_scaling, icol, irow_e, irow_e1, irow_p1, irow_e2, 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, force_scaling, mass_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, force_scaling, mass_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, assembly.elements[ielem],
            force_scaling, mass_scaling, icol, irow_e, irow_e1, irow_p1, irow_e2, irow_p2)
    end

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

    return jacob
end


text:
@testset "dist DSL" begin
  @dist f(x) = exp(normal(x, 0.001))
  @test isapprox(1, f(0); atol = 5)

  @dist labeled_cat(labels, probs) = labels[categorical(probs)]
  @test labeled_cat([:a, :b], [0., 1.]) == :b

  dict = Dict(1 => :a, 2 => :b)
  @dist dict_cat(probs) = dict[categorical(probs)]
  @test dict_cat([0., 1.]) == :b

  @enum Fruit apple orange
  @dist enum_cat(probs) = Fruit(categorical(probs) - 1)
  @test enum_cat([0., 1.]) == orange
end

ext: jl

max_stars_repo_name: cameronraysmith/Gen

max_issues_repo_name: cameronraysmith/Gen

max_forks_repo_name: cameronraysmith/Gen


text:
module JUDI4Flux

    using JUDI: judiJacobian, judiPDEfull, judiVector, judiModeling
    using Zygote, Flux, JOLI
    using Zygote: @adjoint
    import Base.*

    export ExtendedQForward, ExtendedQAdjoint, Forward

    function my_norm(x; dt=1, p=2)
        x = dt * sum(abs.(vec(x)).^p)
        return x^(1.f0/p)
    end

    function convert_to_cell(w)
        nsrc = size(w, 4)
        w_cell = Array{Any}(undef, nsrc)
        for j=1:nsrc
            w_cell[j] = w[:,:,1,j]
        end
        return w_cell
    end

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

    # Linearized Born scattering
    @adjoint *(J::judiJacobian, x::AbstractVecOrMat) = *(J, x), Δ -> (nothing, transpose(J) * Δ)

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

    # Extended source forward modeling
    struct ExtendedQForward{T1}
        F::T1   # forward modeling operator
        m_grad::Bool
    end

    ExtendedQForward(F; m_grad=true) = ExtendedQForward(F, m_grad)

    # Extended source forward modeling: forward mode
    function (EQF::ExtendedQForward)(w::AbstractArray, m::AbstractArray)
        Flocal = deepcopy(EQF.F)
        Flocal.model.m .= m[:, :, 1,1]
        out = Flocal * vec(w)
        nt = Flocal.recGeometry.nt[1]
        nrec = length(Flocal.recGeometry.xloc[1])
        return reshape(out, nt, nrec, 1, Flocal.info.nsrc)
    end

    function grad_w(EQF::ExtendedQForward, m, Δd)
        Flocal = deepcopy(EQF.F)
        Flocal.model.m .= m[:,:,1,1]
        Δw = adjoint(Flocal) * vec(Δd)
        return reshape(Δw, EQF.F.model.n[1], EQF.F.model.n[2], 1, EQF.F.info.nsrc)
    end

    function grad_m(EQF::ExtendedQForward, w, m, Δd)
        Flocal = deepcopy(EQF.F)
        Flocal.model.m .= m[:,:,1,1]
        J = judiJacobian(Flocal, w[:,:,1,:])
        Δm = adjoint(J) * vec(Δd)
        return reshape(Δm, EQF.F.model.n[1], EQF.F.model.n[2], 1, 1)
    end

    @adjoint function (EQF::ExtendedQForward)(w::AbstractArray, m::AbstractArray)
        if EQF.m_grad
            return EQF(w, m), Δ -> (nothing, grad_w(EQF, m, Δ), grad_m(EQF, w, m, Δ))
        else
            return EQF(w, m), Δ -> (nothing, grad_w(EQF, m, Δ), nothing)
        end
    end

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

    # Extended source adjoint modeling
    struct ExtendedQAdjoint{T1}
        F::T1   # adjoint modeling operator
        m_grad::Bool
    end

    ExtendedQAdjoint(F; m_grad=true) = ExtendedQAdjoint(F, m_grad)

    # Extended source adjoint modeling: forward mode
    function (EQT::ExtendedQAdjoint)(d::AbstractArray, m::AbstractArray)
        Flocal = deepcopy(EQT.F)
        Flocal.model.m .= m[:,:,1,1]
        out = adjoint(Flocal)*vec(d)
        out = reshape(out, Flocal.model.n[1], Flocal.model.n[2], 1, Flocal.info.nsrc)
        return out
    end

    function grad_d(EQT::ExtendedQAdjoint, m, Δw)
        Flocal = deepcopy(EQT.F)
        Flocal.model.m .= m[:,:,1,1]
        Δd = Flocal * vec(Δw)
        nt = Flocal.recGeometry.nt[1]
        nrec = length(Flocal.recGeometry.xloc[1])
        return reshape(Δd, nt, nrec, 1, EQT.F.info.nsrc)
    end

    function grad_m(EQT::ExtendedQAdjoint, d, m, Δw)
        Flocal = deepcopy(EQT.F)
        Flocal.model.m .= m[:,:,1,1]
        J = judiJacobian(Flocal, Δw[:,:,1,:])
        Δm = adjoint(J) * vec(d)
        return reshape(Δm, EQT.F.model.n[1], EQT.F.model.n[2], 1, 1)
    end

    @adjoint function (EQT::ExtendedQAdjoint)(d::AbstractArray, m::AbstractArray)
        if EQT.m_grad
            return EQT(d, m), Δ -> (nothing, grad_d(EQT, m, Δ), grad_m(EQT, d, m, Δ))
        else
            return EQT(d, m), Δ -> (nothing, grad_d(EQT, m, Δ), nothing)
        end
    end

    # Fixed source forward modeling
    struct Forward{T1}
        F::T1   # forward modeling operator
        q::judiVector
    end

    # Fixed source forward modeling: forward mode
    (FWD::Forward)(m::AbstractArray) = FWD.F(;m=m)* FWD.q
    @adjoint (FWD::Forward)(m::AbstractArray) = FWD(m), Δ -> (nothing, reshape(transpose(judiJacobian(FWD.F(;m=m), FWD.q)) * vec(Δ), size(m)))
    @adjoint *(F::judiModeling, x::judiVector) = *(F, x), Δ -> (nothing, transpose(F) * vec(Δ))
end


text:
module TextPlots
    export plot, logplot

    include("plot.jl")
end


text:
#Constants and setup
using OrdinaryDiffEq
initial = [0.01, 0.01, 0.01, 0.01]
tspan = (0.,100.)

#Define the problem
function double_pendulum_hamiltonian(udot,u,p,t)
    α  = u[1]
    lα = u[2]
    β  = u[3]
    lβ = u[4]
    udot .=
    [2(lα-(1+cos(β))lβ)/(3-cos(2β)),
    -2sin(α) - sin(α+β),
    2(-(1+cos(β))lα + (3+2cos(β))lβ)/(3-cos(2β)),
    -sin(α+β) - 2sin(β)*(((lα-lβ)lβ)/(3-cos(2β))) + 2sin(2β)*((lα^2 - 2(1+cos(β))lα*lβ + (3+2cos(β))lβ^2)/(3-cos(2β))^2)]
end

#Pass to solvers
poincare = ODEProblem(double_pendulum_hamiltonian, initial, tspan)

sol = solve(poincare, Tsit5())

using Plots; gr()
plot(sol, vars=[(0,3),(0,4)], leg=false, plotdensity=10000)

plot(sol, vars=(3,4), leg=false)

f = (t) -> sol(t,idxs=4)

using Optim
opt = optimize(f,18.0,22.0)

println(opt.minimizer)
println(opt.minimum)

f = (t) -> -sol(first(t),idxs=4)
opt2 = optimize(f,0.0,22.0)

plot(sol, vars=(0,4), plotdensity=10000)
scatter!([opt.minimizer],[opt.minimum],label="Local Min")
scatter!([opt2.minimizer],[-opt2.minimum],label="Local Max")

f = (t) -> -sol(first(t),idxs=4)
opt = optimize(f,[20.0],BFGS())

import NLopt, ForwardDiff

count = 0 # keep track of # function evaluations

function g(t::Vector, grad::Vector)
  if length(grad) > 0
    #use ForwardDiff for the gradients
    grad[1] = ForwardDiff.derivative((t)->sol(first(t),idxs=4),t)
  end
  sol(first(t),idxs=4)
end
opt = NLopt.Opt(:GN_ORIG_DIRECT_L, 1)
NLopt.lower_bounds!(opt, [0.0])
NLopt.upper_bounds!(opt, [40.0])
NLopt.xtol_rel!(opt,1e-8)
NLopt.min_objective!(opt, g)
(minf,minx,ret) = NLopt.optimize(opt,[20.0])
println(minf," ",minx," ",ret)
NLopt.max_objective!(opt, g)
(maxf,maxx,ret) = NLopt.optimize(opt,[20.0])
println(maxf," ",maxx," ",ret)

plot(sol, vars=(0,4), plotdensity=10000)
scatter!([minx],[minf],label="Global Min")
scatter!([maxx],[maxf],label="Global Max")

using DiffEqTutorials
DiffEqTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
"""
    sentence_tokenize([lang,] s)

Splits string `s` into sentences using `WordTokenizers.split_sentences`
function to perform the tokenization. If a language `lang` is provided,
it ignores it ;)
"""
sentence_tokenize(s::T) where T<:AbstractString =
    WordTokenizers.split_sentences(s)
sentence_tokenize(lang::S, s::T) where {S<:Language, T<:AbstractString} =
    WordTokenizers.split_sentences(s)

"""
    tokenize_default([lang,] s)

Splits string `s` into tokens on whitespace using `WordTokenizers.tokenize`
function to perform the tokenization. If a language `lang` is provided,
it ignores it ;)
"""
tokenize_default(s::T) where T<:AbstractString = DEFAULT_WT_TOKENIZER(s)
tokenize_default(lang::S, s::T) where {S<:Language, T<:AbstractString} = DEFAULT_WT_TOKENIZER(s)

"""
    tokenize_stringanalysis(doc [;splitter])

Function that quickly tokenizes `doc` based on the splitting
pattern specified by `splitter::RegEx`.
Supported types for `doc` are: `AbstractString`, `Vector{AbstractString}`,
`StringDocument` and `NGramDocument`.
"""
tokenize_stringanalysis(doc::S; splitter::Regex=DEFAULT_TOKENIZATION_REGEX
                       ) where S<:AbstractString =
    strip.(split(doc, splitter, keepempty=false))

tokenize_stringanalysis(doc::Vector{S}; splitter::Regex=DEFAULT_TOKENIZATION_REGEX
                       ) where S<:AbstractString =
    return vcat((tokenize_stringanalysis(words, splitter=splitter) for words in doc)...)

tokenize_stringanalysis(doc::NGramDocument; splitter::Regex=DEFAULT_TOKENIZATION_REGEX) =
    unique!(tokenize_stringanalysis(collect(keys(doc.ngrams)), splitter=splitter))

tokenize_stringanalysis(doc::StringDocument; splitter::Regex=DEFAULT_TOKENIZATION_REGEX) =
    tokenize_stringanalysis(doc.text, splitter=splitter)

""""
    tokenize(doc [;method, splitter])

Tokenizes the document `doc` based on the `mehtod` (default `:default`, i.e.
a `WordTokenizers.jl` tokenizer) and the `splitter`, which is a `Regex` used
if `method=:stringanalysis`.
"""
function tokenize(doc;
                  method::Symbol=DEFAULT_TOKENIZER,
                  splitter::Regex=DEFAULT_TOKENIZATION_REGEX)
    if method == :stringanalysis
        return tokenize_stringanalysis(doc, splitter=splitter)
    else
        return tokenize_default(doc)
    end
end


text:
# Will be included explicitly
using Base: @deprecate

# Now the batch dimension is the last axis!
@deprecate RunSR(X, y; kw...) EquationSearch(copy(transpose(X)), y; kw...)
@deprecate sqrtm(x) sqrt_abs(x)
@deprecate logm(x) log_abs(x)
@deprecate logm2(x) log2_abs(x)
@deprecate logm10(x) log10_abs(x)


text:
#does any data processign that is needed
#generally called from ComputeNCCS since it is cacheable
#TODO: Make sure the ratios are logged
function prepsensitivityregressions!(data::NCCSData;
    toosmalltolog::Float64 = TOO_SMALL_TO_LOG)

  #logs if possible or returns missing
  plog(f::MFloat64)::MFloat64 = ismissing(f) || f ≤ toosmalltolog ? missing : log(f)

  #cleans up a common-size ratio
  cleanratio(f::MFloat64) = (!ismissing(f)) && isfinite(f) && f ≤ 1.0  && f > 0.0 ? f : missing

  #log our focal variables
  data.df[:lcontributions] = (plog).(data.df[:contributions])
  data.df[:lprogramexpense] = (plog).(data.df[:programexpense])
  data.df[:lprogramexpense_avg] = (plog).(data.df[:programexpense_avg])

  #taking the average for control purposes only with valid contribution fields
  #data.df[:lregnetassets] = ((c::MFloat64, a::MFloat64)->
  #  ismissing(c) ? missing : a).(data.df[:contributions], data.df[:lnetassets])
  lagwithin!(data, :netassets)

  #scale everythign to lagged net assets
  data.df[:pcontributions] = data.df[:contributions] ./ data.df[:Lnetassets]
  data.df[:pcontributions] .= (cleanratio).(data.df[:pcontributions])
  data.df[:lpcontributions] = (plog).(data.df[:pcontributions])

  data.df[:pprogramexpense] = data.df[:programexpense] ./ data.df[:Lnetassets]
  data.df[:pprogramexpense] .= (cleanratio).(data.df[:pprogramexpense])
  data.df[:lpprogramexpense] = (plog).(data.df[:pprogramexpense])

  #get the average expense
  aggdf::DataFrame = aggregatefielddf(data.df[[:ein, :lpprogramexpense,
    :pprogramexpense]], rejoin=false)
  unique!(aggdf)
  data.df = join(data.df, aggdf, on=[:ein], kind=:left)

  #lag our focal variables
  lagwithin!(data, :lcontributions)
  lagwithin!(data, :pcontributions)
  lagwithin!(data, :lpcontributions)
  lagwithin!(data, :publicsupport)
  lagwithin!(data, :lprogramexpense)
  lagwithin!(data, :pprogramexpense)
  lagwithin!(data, :lpprogramexpense)

  println(describe(data.df[[:ein, :lpprogramexpense, :lpprogramexpense_avg, :Llpprogramexpense]]))

  #generate the interactions
  data.df[:LpublicsupportXLlcont] = data.df[:Lpublicsupport] .* data.df[:Llcontributions]
  data.df[:LpublicsupportXLpcont] = data.df[:Lpublicsupport] .* data.df[:Lpcontributions]
  data.df[:publicsupport_avgXLlcont] = data.df[:publicsupport_avg] .* data.df[:Llcontributions]
  data.df[:publicsupport_avgXLpcont] = data.df[:publicsupport_avg] .* data.df[:Lpcontributions]

  data.df[:LlprogramexpenseXLlcont] = data.df[:Llprogramexpense] .* data.df[:Llcontributions]
  data.df[:lprogramexpense_avgXLlcont] = data.df[:lprogramexpense_avg] .* data.df[:Llcontributions]

  data.df[:LlpprogramexpenseXLlcont] = data.df[:Llpprogramexpense] .* data.df[:Llcontributions]
  data.df[:LpprogramexpenseXLpcont] = data.df[:Lpprogramexpense] .* data.df[:Lpcontributions]
  data.df[:lpprogramexpense_avgXLlcont] = data.df[:lpprogramexpense_avg] .* data.df[:Llcontributions]
  data.df[:pprogramexpense_avgXLpcont] = data.df[:pprogramexpense_avg] .* data.df[:Lpcontributions]

  #generate the first differences
  #differencewithin!(data, :lcontributions)
  #differencewithin!(data, :publicsupport)
  #differencewithin!(data, :lprogramexpense)
  #differencewithin!(data, :lprogramexpense)

  #I think regressions work better with symbols instead of strings
  rename!(data.df, :ein=>:ein_old)
  data.df[:ein] = (Symbol).(data.df[:ein_old])
  deletecols!(data.df, :ein_old)

  return nothing
end

#runs the basic panel form regressions
function olssensitivitypanelregPE(df::AbstractDataFrame)
  specs::FMSpecs = FMSpecs(6)

  #set up the results
  push!(specs,
    xspec=Meta.parse("pprogramexpense_avgXLpcont + Lpcontributions"),
    xnames = [:intercept, :LagContribAvgExp, :LagContrib],
    clusterspec = :ein,
    withinspec = :ein,
    yspec=:lreturn)

  push!(specs,
    xspec=Meta.parse("LpprogramexpenseXLpcont + Lpcontributions + Lpprogramexpense"),
    xnames = [:intercept, :LagContribExp, :LagContrib],
    clusterspec = :ein,
    withinspec = :ein,
    yspec=:lreturn)

  push!(specs,
    xspec=Meta.parse("pprogramexpense_avgXLpcont + Lpcontributions"),
    xnames = [:intercept, :LagContribAvgExp, :LagContrib],
    clusterspec = :ein,
    withinspec = :ein,
    yspec=:lreturnnetbmsp5001yr)

  push!(specs,
    xspec=Meta.parse("LpprogramexpenseXLpcont + Lpcontributions + Lpprogramexpense"),
    xnames = [:intercept, :LagContribExp, :LagContrib],
    clusterspec = :ein,
    withinspec = :ein,
    yspec=:lreturnnetbmsp5001yr)

  push!(specs,
    xspec=Meta.parse("pprogramexpense_avgXLpcont + Lpcontributions"),
    xnames = [:intercept, :LagContribAvgExp, :LagContrib],
    clusterspec = :ein,
    withinspec = :ein,
    yspec=:lreturnnetbmff31yr)

  push!(specs,
    xspec=Meta.parse("LpprogramexpenseXLpcont + Lpcontributions + Lpprogramexpense"),
    xnames = [:intercept, :LagContribExp, :LagContrib],
    clusterspec = :ein,
    withinspec = :ein,
    yspec=:lreturnnetbmff31yr)

  #run the regressions
  computeFMLMresults!(df, specs)

  return specs
end

#runs the basic panel form regressions
function olssensitivitypanelregPS(df::AbstractDataFrame)
  specs::FMSpecs = FMSpecs(6)

  #set up the results
  #NOTE: delete the l before public support on next pass
  push!(specs,
    xspec=Meta.parse("publicsupport_avgXLpcont + Lpcontributions"),
    xnames = [:intercept, :LagContribAvgSup, :LagContrib],
    withinspec = :ein,
    clusterspec = :ein,
    yspec=:lreturn)

  push!(specs,
    xspec=Meta.parse("LpublicsupportXLpcont + Lpcontributions + Lpublicsupport"),
    xnames = [:intercept, :LagContribSup, :LagContrib],
    clusterspec = :ein,
    withinspec = :ein,
    yspec=:lreturn)

  push!(specs,
    xspec=Meta.parse("publicsupport_avgXLpcont + Lpcontributions"),
    xnames = [:intercept, :LagContribAvgSup, :LagContrib],
    clusterspec = :ein,
    withinspec = :ein,
    yspec=:lreturnnetbmsp5001yr)

  push!(specs,
    xspec=Meta.parse("LpublicsupportXLpcont + Lpcontributions + Lpublicsupport"),
    xnames = [:intercept, :LagContribSup, :LagContrib],
    clusterspec = :ein,
    withinspec = :ein,
    yspec=:lreturnnetbmsp5001yr)

  push!(specs,
    xspec=Meta.parse("publicsupport_avgXLpcont + Lpcontributions"),
    xnames = [:intercept, :LagContribAvgSup, :LagContrib],
    clusterspec = :ein,
    withinspec = :ein,
    yspec=:lreturnnetbmff31yr)

  push!(specs,
    xspec=Meta.parse("LpublicsupportXLpcont + Lpcontributions + Lpublicsupport"),
    xnames = [:intercept, :LagContribSup, :LagContrib],
    clusterspec = :ein,
    withinspec = :ein,
    yspec=:lreturnnetbmff31yr)

  #run the regressions
  computeFMLMresults!(df, specs)

  return specs
end

#runs the basic panel form regressions
function olssensitivity2pass(df::DataFrame;
  rhs2ndpassfields_avg::Vector{Symbol} = [:pprogramexpense_avg],#, :publicsupport_avg=#],
  rhs2ndpassfields::Vector{Symbol} = [:Lpprogramexpense#=, :Lpublicsupport_avg=#],
  winsorizeprop::Float64 = WINSORIZE_PROP)

  specs::FMSpecs = FMSpecs(6)

  local βdf::DataFrame
  local βdfcross::DataFrame
  local focalfields::Vector{Symbol}
  #run the first pass regressions to get the betas
  #println("q expense complete: ",
  #  sum(completecases(df[[:lreturn, :Llcontributions, :lprogramexpense_avg]])))

  regressbynonprofit!(df,
    YField = :lreturn,
    XFields = [:Lpcontributions],
    βNames = [:b0temp1, :b1lRetLCon],
    groupbyfield = :ein,
    predictwithresults = false,
    checkunique=false)

  regressbynonprofit!(df,
    YField = :lreturnnetbmsp5001yr,
    XFields = [:Lpcontributions],
    βNames = [:b0temp2, :b1lRetNetSP500LCon],
    groupbyfield = :ein,
    predictwithresults = false,
    checkunique=false)

  regressbynonprofit!(df,
    YField = :lreturnnetbmff31yr,
    XFields = [:Lpcontributions],
    βNames = [:b0temp3, :b1lRetNetFF3LCon],
    groupbyfield = :ein,
    predictwithresults = false,
    checkunique=false)

  deletecols!(df, [:b0temp1, :b0temp2, :b0temp3])

  βfields::Vector{Symbol} = [:b1lRetLCon, :b1lRetNetSP500LCon,
    :b1lRetNetFF3LCon]

  (f->winsorize!(df[f], prop=winsorizeprop)).(βfields)
  focalfields = [:ein; βfields]

  #now fold in the 2nd pass regression terms through repeated inner joins and uniqueness filters
  βdf = unique(df[focalfields])
  for f ∈ rhs2ndpassfields_avg
    subdf::DataFrame = view(df, completecases(df[:,[:ein,f]]), [:ein, f])
    βdf = join(βdf, unique(subdf), on=:ein, kind=:left)
  end

  CSV.write("$WORKING_PATH\\betatest.csv", βdf)

  #make sure the joins went as planned
  @assert size(unique(βdf),1) == size(βdf,1)
  focalfields = [:ein; :fisyr; βfields]

  βdfcross = df[focalfields]

  for f ∈ rhs2ndpassfields
    subdf::DataFrame = view(df, completecases(df[:,[:ein,:fisyr,f]]), [:ein,:fisyr,f])
    βdfcross = join(βdfcross, unique(subdf), on=[:ein, :fisyr], kind=:left)
  end

  #make sure all of the betas within each entity are unique
  @assert size(unique(βdfcross),1) == size(βdfcross,1)

  #make the fisyrs into symbols
  rename!(βdfcross, [:fisyr=>:fisyrold])
  βdfcross[:fisyr] = (y::Int->Symbol("y", y)).(βdfcross[:fisyrold])
  deletecols!(βdfcross, :fisyrold)

  #Regress the betas on the focal variables
  push!(specs,
    xspec=Meta.parse("pprogramexpense_avg"),
    xnames = [:intercept, :AvgExp],
    yspec=:b1lRetLCon)

  push!(specs,
    xspec=Meta.parse("Lpprogramexpense"),
    xnames = [:intercept, :LagExp],
    withinspec = :fisyr,
    yspec=:b1lRetLCon)

  push!(specs,
    xspec=Meta.parse("pprogramexpense_avg"),
    xnames = [:intercept, :AvgExp],
    yspec=:b1lRetNetSP500LCon)

  push!(specs,
    xspec=Meta.parse("Lpprogramexpense"),
    xnames = [:intercept, :LagExp],
    withinspec = :fisyr,
    yspec=:b1lRetNetSP500LCon)

  push!(specs,
    xspec=Meta.parse("pprogramexpense_avg"),
    xnames = [:intercept, :AvgExp],
    yspec=:b1lRetNetFF3LCon)

  push!(specs,
    xspec=Meta.parse("Lpprogramexpense"),
    xnames = [:intercept, :LagExp],
    withinspec = :fisyr,
    yspec=:b1lRetNetFF3LCon)

  println(describe(βdf))
  println(describe(βdfcross))

  #run the regressions
  βdfs::Vector{DataFrame} = [βdf, βdfcross, βdf, βdfcross, βdf, βdfcross]
  computeFMLMresults!(βdfs, specs)

  return specs
end

#primarily generates the gregeression table
function olssensitivitypanel(df::AbstractDataFrame;
  decimals::Int=DECIMALS,
  suffix::String = "",
  displayedvars::Union{Vector{Symbol},Nothing} = nothing,
  contentRowNames::Union{Vector{String},Nothing}=nothing,
  displayedvarsscaling::Union{Vector{Float64},Nothing} = nothing,
  regressionfunction::Function = olssensitivitypanelreg,
  titleCaption::String = "TBD",
  outputpath::String = TABLE_PATH,
  errorfunction::Function = getClustered!,
  scalingvalue::Float64 = 10000.0)

  #first run the regressions
  local specs::FMSpecs = regressionfunction(df)

  #apply the defaults based on the sepcs
  isnothing(displayedvars) &&  (displayedvars = setdiff(unique([specs.xnames...;]), [:intercept]))
  isnothing(contentRowNames) &&  (contentRowNames = (string).(displayedvars))
  isnothing(displayedvarsscaling) && (displayedvarsscaling = fill(scalingvalue, length(displayedvars)))

  local descRowNames::Vector{String} = [
    "Nonprofit(NP) FE",
    "Year FE",
    "NP Clustered",
    "\$R^{2}\$ (within)",
    "N"]

  #allocate space for the descriptive rows
  local descContent::Vector{Vector{String}} =
    ((i::Int)->Vector{String}(undef, specs.N[])).(1:length(descRowNames))

    #this builds the descriptive rows. There is Probably a better way to do this,
  #but its fairly specific to the project.
  for i ∈ 1:specs.N[]
    r::Int = 0
    descContent[r+=1][i] = "$(specs.withinspecs[i]==:ein ||
        occursin("ein", string(specs.xspecs[i])) ? "X" : "")"
    descContent[r+=1][i] = "$(specs.withinspecs[i]==:fisyr ||
        occursin("fisyr", string(specs.xspecs[i])) ? "X" : "")"
    descContent[r+=1][i] = "$(specs.clusterspecs[i]==:ein ? "X" : "")"
    descContent[r+=1][i] =
      "$(num2Str(getR²(specs.results[i], adjusted=false),  decimals))"
    descContent[r+=1][i] = "$(specs.results[i].N)"
  end

  #NOTE delete this if the console gets too cluttered
  display(descContent)

  local colNames::Vector{Vector{String}} = [
    ["return", "SP500 Adj", "FF3 Adj"],
    specs.specnames]

  local widthColNames::Vector{Vector{Int}} = [[2,2,2], ones(Int, specs.N[])]

  local alignmentColNames::Vector{Vector{String}} =
    [fill("c",length(widthColNames[1])), fill("r", specs.N[])]

  #checks based on the dimensions of the table parameters
  @assert (minimum((sum).(widthColNames)) == specs.N[] &&
    maximum((sum).(widthColNames)) == specs.N[] &&
    length(contentRowNames) == length(displayedvars)
    )

  local tabletext::String = texTable(specs.results,
    errorfunction#=getNeweyWestFunc(5) (lm::CTLM)->getNeweyWestSlow(lm, 5)=#,
    displayedvars,
    titleCaption = titleCaption,
    colNames = colNames,
    widthColNames = widthColNames,
    contentRowNames = contentRowNames,
    descRowNames = descRowNames,
    descContent = descContent,
    decimalDigits = 1,
    scaling= displayedvarsscaling,
    stars=true,
    starStrings = OVERRIDE_STAR_STRINGS,
    #clearMem = USE_AGGRESSIVE_GC,
    caption = "to be written",
    nakedTable=true)

  writeNakedTable(tabletext, path=outputpath, outName="olspanel$suffix.tex")
  println("Regressions for olspanel$suffix.tex complete and recorded.")

  return tabletext
end

function prefilterforolssensitivity(df::DataFrame;
  minpointspluskforβ::Int = MIN_POINTS_PLUS_K_FOR_BETA,
  essentialfields::Vector{Symbol}=[:Llcontributions, :Lpcontributions])::DataFrame

  local subdf::SubDataFrame
  local ssubdf::SubDataFrame

  #make sure we ahve enough unique points for each organization
  df[:useforregressions] = trues(size(df,1))

  checkfields::Vector{Symbol} = [:lreturn; essentialfields; :useforregressions]

  #loop over all unique eins
  for subdf ∈ groupby(df, :ein)

    ssubdf = view(subdf, completecases(subdf[checkfields]), checkfields)
    subN = minimum((f->size(unique(ssubdf[f]),1)).(essentialfields))

    #make sure we have enough points for the regression
    if subN < minpointspluskforβ + 2 #add two for the existing dof
      ssubdf[:useforregressions] .= false
    end
  end

  return df[df[:useforregressions],:]
end

function makeolspaneltables(data::NCCSData)

  #this gets repeated a lot so we can save time doing it now
  regdf::DataFrame = prefilterforolssensitivity(data.df)

  olssensitivitypanel(regdf, suffix = "programexpense",
    #displayedvars = nothing,
    #contentRowNames = nothing,
    regressionfunction = olssensitivitypanelregPE,
    titleCaption = "Program Expense vs Contribution Return Sensitivity")

  #NOTE: Failed proxy
  olssensitivitypanel(regdf, suffix = "publicsupport",
    #displayedvars = nothing,
    #contentRowNames = nothing,
    regressionfunction = olssensitivitypanelregPS,
    titleCaption = "PublicSupport vs Contribution Return Sensitivity")

  olssensitivitypanel(regdf, suffix = "2pass",
    #displayedvars = nothing,
    #contentRowNames = nothing,
    regressionfunction = olssensitivity2pass,
    errorfunction = getWhiteΣ!,
    titleCaption = "Effect of public support and program expense on sensitivity")

end

function olssensitivitytables(data::NCCSData)
  makeolspaneltables(data)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Golgi
using Test

@testset "Golgi.jl" begin
    # Write your tests here.
end

using Enzyme
using ChainRules
Enzyme.pullback(sin, 2.0)(1.0)
@time Enzyme.pullback(tanh, 2.0)(1.0)
ChainRules.rrule(sin, 2.0)[2](1.0)

f(x) = cos(sin(x))

Golgi.rrule_call(f, 2.0)
Golgi.rrule_call(sin, 2.0)

foo(x) = x
ChainRules.rrule(foo, 2.0)

using Enzyme

struct Foo
    x::Float64
end

function (f::Foo)(x::Float64)
    return f.x + x
end

Enzyme.pullback(Foo(1.0), 2.0)(1.0)

@generated function goo(xs...)
    quote
        return $(xs)
    end
end

@generated function foo(x)
    quote
        goo(x)
    end
end

using Zygote

Zygote.pullback(sin, 2.0)
goo(1, 2, 3.0)

function rrule(::Type{f}, args...) where {f <: Operation}
end

function rrule(::Type{typeof(apply!)}, r::AbstractRegister, op::Operation)
    apply!(r, op), function pullback(Δ)
        reverse_r, grad = Δ
        apply!(reverse_r, adjoint(op))
        # ...
        # NOTE: grad_op is of type ∇Operation
        # which is a struct contains same field
        # of Operation, but are gradients of the fields
        return NoTangent(), grad_r, grad_op
    end
end

function rrule(f::Type{F}, args...) where {F <: Operation}
    f(args...), function pullback(Δ::∇Operation)
        return unpack(Δ) # unpack the struct back to a Tuple in argument order
    end
end

# this only works when return value is an expectation or measure
# error when return value is not a Real
# how do we extract expectation loops from function body?
# canonicalize measure accumulation like for _ in 1:1000::UnitRange; measure end? but how
function frule((∇self, ∇r, ∇op), ::Type{typeof(apply!)}, r::AbstractRegister, op::Operation)
    apply!(r, op)
    # ∇r can be NoTangent() if compile target is a quantum device
    # how do we add/sub pi/2 in the IR?
end


text:
function DEMO(; kwargs...)
  cur_dict = Dict{Symbol, AbstractFloat}(
    :N_G => 1.2,
    :kappa => 1.59,
    :epsilon => 0.323,
    :beta_N => 0.0259,
    :Q => 39.9,
    :f_D => 0.775,
    :P_F => 2037,
    :B_max => 12.3,
    :B_max_hat => 12.9,
    :sigma_max => 660e6,
    :sigma_max_hat => 660e6,
    :tau_P => 2,
    :b => 1.63,
    :nu_n => 0.27,
    :nu_T => 1.094,
    :eta_T => 0.4,
    :eta_RF => 0.75
  )

  merge!(cur_dict, kwargs)

  Reactor(; cur_dict...)
end


text:
using Documenter, JET

makedocs(; modules = [JET],
           sitename="JET.jl",
           pages = [
             "README" => "index.md",
             "Usages" => "usages.md",
             "Configurations" => "config.md",
             "Internals" => "internals.md",
           ],
           format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
           )

deploydocs(; repo = "github.com/aviatesk/JET.jl.git",
             push_preview = true,
             )


text:
# The file test the uniform inflow condition

include("const.jl")
include("mathfunctions.jl")
include("flapreponse.jl")
include("uniforminflow.jl")
include("aoaget.jl")
include("clcdget.jl")
include("rotorforce.jl")
include("trim_windtu.jl")

tic() # 程序起始标志
print("===计算初始化===")

ψ = 0.0
beta_lat = 0.0
beta_lon = 0.0
θcp = zeros(Float64,NR)
θlat = zeros(Float64,NR)
θlon = zeros(Float64,NR)
θ0 = Array{Float64}(NR,Nbe)
for k in 1:NR
  θcp[k] = θ7
  θlat[k] = thelat
  θlon[k] = thelon
  for i in 1:Nbe
    θ0[k,i] = θcp[k]+θtw*(rb[k,i]/R-0.7)
  end
end
index = 1
betanow = βang0
betaxnow = dβ0
rot = 0.0
blat = 0.0
blon = 0.0
uitmp = uniforminflow(ψ,θcp,θlat,θlon,betanow,betaxnow)

while index<=60*npsi # 如果计算了100圈还没收敛就结束计算
  vind_r = uitmp[1]
  vall_r = uitmp[2]
  beta = uitmp[4]
  # if index%npsi==0
  #   print("===当前挥舞角$(beta/π*180)°===\n")
  # end
  dbeta = uitmp[5]
  blat = blat+uitmp[4]*sin(ψ)
  blon = blon+uitmp[4]*cos(ψ)

  # This file is used for get angle-of-attack
  aoatmp = aoaget(vall_r,ψ,θ0,θlat,θlon,beta)
  α_aero = aoatmp[1]
  θ = aoatmp[2]

  # This file is to get the Cl and Cd of Blade Element
  clcdtmp = clcdget(α_aero,vall_r)
  Cl = clcdtmp[1]
  Cd = clcdtmp[2]

  # This file calculate the whole aerodynamic force of the rotor in the hub
  rftmp = rotorforce(ψ,vall_r,α_aero,θ,Cl,Cd)
  rot = rot+rftmp[3]
  Mbeta_aero = rftmp[5][1]
  # if index%npsi==0
  #   print("===当前挥舞力矩：$(Mbeta_aero)N·m===\n")
  # end
  # print("$(θ/π*180)\n")
  # print("$(α_aero/π*180)\n")
  # print("\n")

  # 此处要输入转过一周进行配平的条件，今天来不及了明天完成，作此标志 12/2/2017
  if index%npsi==0
    rot = rot/npsi
    print("===当前拉力：$(rot)N===\n")
    beta_lat = blat/npsi*2
    beta_lon = blon/npsi*2
    print("===当前横向挥舞：$(beta_lat*180/π)deg===\n")
    print("===当前纵向挥舞：$(beta_lon*180/π)deg===\n")

    # rot = 0.0
    # blat = 0.0
    # blon = 0.0

    ψ = 0.0

    # 此处开始进行配平
    trimtmp = trimwt(uitmp,rot,beta_lat,beta_lon,θcp,θlat,θlon)
    if trimtmp[1]
      # print(trimtmp)
      print("配平总距：$(trimtmp[5]*180/π)\n")
      print("配平横向变距：$(trimtmp[3]*180/π)\n")
      print("配平纵向变距：$(trimtmp[4]*180/π)\n")
      break
    else
      θcp = trimtmp[5]
      θ0 = trimtmp[2]
      print("配平总距：$(trimtmp[5]*180/π)\n")
      print("配平横向变距：$(trimtmp[3]*180/π)\n")
      print("配平纵向变距：$(trimtmp[4]*180/π)\n")
      θlat = trimtmp[3]
      θlon = trimtmp[4]
      rot = 0.0
      blat = 0.0
      blon = 0.0
    end
  end

  # The file calculate the uniform induced velocity
  uitmp = uniforminflow(ψ,θcp,θlat,θlon,beta,dbeta,Mbeta_aero)
  ψ = ψ+dψ
  index = index+1
end

print("===计算收敛===\n")
print("===总计算步数$(index-1)===\n")
toc() # 程序结束

max_stars_repo_name: dingzhw/RotorTrim

max_issues_repo_name: helicopterliu/RotorTrim

max_forks_repo_name: helicopterliu/RotorTrim


text:
const Connectivity{T} =
    Union{AbstractVector{T}, NTuple{N,T} where N} where {T <: Integer}

"""
    MeshCell

Single cell element in unstructured or polygonal grid.

It is characterised by a cell type (for instance, `VTKCellType.TRIANGLE` or
`PolyData.Strips`) and by a connectivity vector determining the points on the
grid defining this cell.

---

    MeshCell(cell_type, connectivity::AbstractVector)

Define a single cell element of an unstructured grid.

The `cell_type` argument characterises the type of cell (e.g. vertex, triangle,
hexaedron, ...):

- cell types for unstructured datasets are defined in the [`VTKCellTypes`](@ref)
module;
- cell types for polygonal datasets are defined in the [`PolyData`](@ref) module.

The `connectivity` argument is a vector containing the indices of the points
passed to [`vtk_grid`](@ref) which define this cell.

# Example

Define a triangular cell passing by points with indices `[3, 5, 42]`.

```jldoctest
julia> cell = MeshCell(VTKCellTypes.VTK_TRIANGLE, [3, 5, 42])
MeshCell{VTKCellType,Array{Int64,1}}(VTKCellType("VTK_TRIANGLE", 0x05, 3), [3, 5, 42])
```
"""
struct MeshCell{CellType, V <: Connectivity}
    ctype::CellType  # cell type identifier (see VTKCellTypes.jl)
    connectivity::V  # indices of points (one-based, following the convention in Julia)
    function MeshCell(ctype, conn)
        if nodes(ctype) ∉ (length(conn), -1)
            error("Wrong number of nodes in connectivity vector.")
        end
        C = typeof(ctype)
        V = typeof(conn)
        new{C,V}(ctype, conn)
    end
end

Base.eltype(::Type{<:MeshCell}) = VTKCellTypes.VTKCellType
cell_type(cell::MeshCell) = cell.ctype

# By default, MeshCell are attached to unstructured grids.
grid_type(::Type{<:MeshCell}) = VTKUnstructuredGrid()

function add_cells!(vtk, xml_piece, number_attr, xml_name, cells;
                    with_types::Val=Val(true))
    Cell = eltype(cells)
    Ncls = length(cells)
    write_types = with_types === Val(true)

    # Create data arrays.
    offsets = Array{Int32}(undef, Ncls)

    # Write `types` array? This must be true for unstructured grids,
    # and false for polydata.
    if write_types
        types = Array{UInt8}(undef, Ncls)
    end

    Nconn = 0     # length of the connectivity array
    if Ncls >= 1  # it IS possible to have no cells
        offsets[1] = length(cells[1].connectivity)
    end

    for (n, c) in enumerate(cells)
        Npts_cell = length(c.connectivity)
        Nconn += Npts_cell
        if write_types
            types[n] = cell_type(c).vtk_id
        end
        if n >= 2
            offsets[n] = offsets[n-1] + Npts_cell
        end
    end

    # Create connectivity array.
    conn = Array{Int32}(undef, Nconn)
    n = 1
    for c in cells, i in c.connectivity
        # We transform to zero-based indexing, required by VTK.
        conn[n] = i - 1
        n += 1
    end

    # Add arrays to the XML file (DataArray nodes).
    set_attribute(xml_piece, number_attr, Ncls)

    xnode = new_child(xml_piece, xml_name)
    data_to_xml(vtk, xnode, conn, "connectivity")
    data_to_xml(vtk, xnode, offsets, "offsets")
    if write_types
        data_to_xml(vtk, xnode, types, "types")
    end

    vtk
end

const CellVector = AbstractVector{<:MeshCell}

# Possible specifications for unstructured points.
const UnstructuredCoords = Union{
    AbstractMatrix,  # array with dimensions (3, Np)
    Tuple{Vararg{AbstractVector{T},3}} where T,  # tuple of 3 vectors of length Np
}

function num_points(::UnstructuredVTKDataset, x::AbstractMatrix)
    if size(x, 1) != 3
        throw(DimensionMismatch("first dimension must have length 3"))
    end
    size(x, 2)
end

function num_points(::UnstructuredVTKDataset, x::Tuple)
    Ns = map(length, x)
    N = first(Ns)
    if !all(Ns .== N)
        throw(DimensionMismatch("tuples (x, y, z) must have the same length"))
    end
    N
end

function vtk_grid(dtype::VTKUnstructuredGrid, filename::AbstractString,
                  points::UnstructuredCoords, cells::CellVector;
                  kwargs...)
    Npts = num_points(dtype, points)
    Ncls = length(cells)

    xvtk = XMLDocument()
    vtk = DatasetFile(dtype, xvtk, filename, Npts, Ncls; kwargs...)

    xroot = vtk_xml_write_header(vtk)
    xGrid = new_child(xroot, vtk.grid_type)

    xPiece = new_child(xGrid, "Piece")
    set_attribute(xPiece, "NumberOfPoints", vtk.Npts)

    xPoints = new_child(xPiece, "Points")
    data_to_xml(vtk, xPoints, points, "Points", 3)

    add_cells!(vtk, xPiece, "NumberOfCells", "Cells", cells)

    vtk
end

# Variant of vtk_grid with 2-D array "points".
#   size(points) = (dim, num_points), with dim ∈ {1, 2, 3}
function vtk_grid(filename::AbstractString, points::AbstractArray{T,2},
                  cells::CellVector, args...; kwargs...) where T
    dim, Npts = size(points)
    gtype = grid_type(eltype(cells))

    if dim == 3
        return vtk_grid(gtype, filename, points, cells, args...; kwargs...)
    end

    xyz = if dim == 2
        (
            view(points, 1, :),
            view(points, 2, :),
            Zeros{T}(Npts),
        )
    elseif dim == 1
        (
            view(points, 1, :),
            Zeros{T}(Npts),
            Zeros{T}(Npts),
        )
    else
        throw(DimensionMismatch(
            "`points` array must be of size (dim, num_points) with dim ∈ 1:3"
        ))
    end

    vtk_grid(gtype, filename, xyz, cells, args...; kwargs...)
end

# Variant of vtk_grid with 1-D arrays x, y, z.
# Size of each array: (num_points)
function vtk_grid(filename::AbstractString, x::AbstractVector{T},
                  y::AbstractVector{T}, z::AbstractVector{T},
                  cells::CellVector, args...; kwargs...) where {T}
    if !(length(x) == length(y) == length(z))
        throw(DimensionMismatch("length of x, y and z arrays must be the same."))
    end
    points = (x, y, z)
    gtype = grid_type(eltype(cells))
    vtk_grid(gtype, filename, points, cells, args...; kwargs...)
end

# 2D version
function vtk_grid(filename::AbstractString, x::AbstractVector{T},
                  y::AbstractVector{T}, cells::CellVector, args...;
                  kwargs...) where {T}
    N = length(x)
    vtk_grid(filename, x, y, Zeros{T}(N), cells, args...; kwargs...)
end

# 1D version
function vtk_grid(filename::AbstractString, x::AbstractVector{T},
                  cells::CellVector, args...; kwargs...) where {T}
    N = length(x)
    vtk_grid(filename, x, Zeros{T}(N), Zeros{T}(N), cells, args...; kwargs...)
end

# Variant with 4-D Array (for "pseudo-unstructured" datasets, i.e., those that
# actually have a 3D structure) -- maybe this variant should be removed...
function vtk_grid(filename::AbstractString, points::AbstractArray{T,4},
                  cells::CellVector, args...; kwargs...) where T
    dim, Ni, Nj, Nk = size(points)
    points_r = reshape(points, dim, Ni * Nj * Nk)
    gtype = grid_type(eltype(cells))
    vtk_grid(gtype, filename, points_r, cells, args...; kwargs...)
end


text:
function mh_accept(rng::AbstractRNG, H::AbstractFloat, H_new::AbstractFloat)
    logα = min(0, H - H_new)
    return log(rand(rng)) < logα, exp(logα)
end
mh_accept(H::AbstractFloat, H_new::AbstractFloat) = mh_accept(GLOBAL_RNG, logα)

sample(h::Hamiltonian, prop::AbstractProposal, θ::AbstractVector{T}, n_samples::Int; verbose::Bool=true) where {T<:Real} =
    sample(GLOBAL_RNG, h, prop, θ, n_samples; verbose=verbose)

function sample(rng::AbstractRNG, h::Hamiltonian, prop::AbstractProposal, θ::AbstractVector{T}, n_samples::Int; verbose::Bool=true) where {T<:Real}
    θs = Vector{Vector{T}}(undef, n_samples)
    Hs = Vector{T}(undef, n_samples)
    αs = Vector{T}(undef, n_samples)
    time = @elapsed for i = 1:n_samples
        θs[i], Hs[i], αs[i] = step(rng, h, prop, i == 1 ? θ : θs[i-1])
    end
    verbose && @info "Finished sampling with $time (s)" typeof(h.metric) typeof(prop) EBFMI(Hs) mean(αs)
    return θs
end

sample(h::Hamiltonian, prop::AbstractProposal, θ::AbstractVector{T}, n_samples::Int, adaptor::Adaptation.AbstractAdaptor,
       n_adapts::Int=min(div(n_samples, 10), 1_000); verbose::Bool=true) where {T<:Real} =
       sample(GLOBAL_RNG, h, prop, θ, n_samples, adaptor, n_adapts; verbose=verbose)

function sample(rng::AbstractRNG, h::Hamiltonian, prop::AbstractProposal, θ::AbstractVector{T}, n_samples::Int, adaptor::Adaptation.AbstractAdaptor,
                n_adapts::Int=min(div(n_samples, 10), 1_000); verbose::Bool=true) where {T<:Real}
    θs = Vector{Vector{T}}(undef, n_samples)
    Hs = Vector{T}(undef, n_samples)
    αs = Vector{T}(undef, n_samples)
    time = @elapsed for i = 1:n_samples
        θs[i], Hs[i], αs[i] = step(rng, h, prop, i == 1 ? θ : θs[i-1])
        if i <= n_adapts
            adapt!(adaptor, θs[i], αs[i])
            h, prop = update(h, prop, adaptor)
            if verbose
                if i == n_adapts
                    @info "Finished $n_adapts adapation steps" typeof(adaptor) prop.integrator.ϵ h.metric
                elseif i % Int(n_adapts / 10) == 0
                    @info "Adapting $i of $n_adapts steps" typeof(adaptor) prop.integrator.ϵ h.metric
                end
            end
        end
    end
    verbose && @info "Finished $n_samples sampling steps in $time (s)" typeof(h.metric) typeof(prop) EBFMI(Hs) mean(αs)
    return θs
end

function step(rng::AbstractRNG, h::Hamiltonian, prop::TakeLastProposal{I}, θ::AbstractVector{T}) where {T<:Real,I<:AbstractIntegrator}
    r = rand_momentum(rng, h)
    H = hamiltonian_energy(h, θ, r)
    θ_new, r_new = transition(prop, h, θ, r)
    H_new = hamiltonian_energy(h, θ_new, r_new)
    # Accept via MH criteria
    is_accept, α = mh_accept(rng, H, H_new)
    if !is_accept
        return θ, H, α
    end
    return θ_new, H_new, α
end

function step(rng::AbstractRNG, h::Hamiltonian, prop::NUTS{I}, θ::AbstractVector{T}) where {T<:Real,I<:AbstractIntegrator}
    r = rand_momentum(rng, h)
    θ_new, r_new, α = transition(rng, prop, h, θ, r)
    H_new = hamiltonian_energy(h, θ_new, r_new)
    # We always accept in NUTS
    return θ_new, H_new, α
end

step(h::Hamiltonian, p::AbstractProposal, θ::AbstractVector{T}) where {T<:Real} = step(GLOBAL_RNG, h, p, θ)


text:
println(VERSION)
println(VERSION.major)
println(VERSION.minor)
println(VERSION.patch)

VERSION < v"0.7.0" && exit()

using Libdl
using Pkg

println(Base.Sys.BINDIR)
println(Libdl.dlpath(string("lib", splitext(Base.julia_exename())[1])))
println(unsafe_string(Base.JLOptions().image_file))

pkg = Base.PkgId(Base.UUID(0x438e738f_606a_5dbb_bf0a_cddfbfd45ab0), "PyCall")
modpath = Base.locate_package(pkg)
if modpath !== nothing
    PyCall_depsfile = joinpath(dirname(modpath),"..","deps","deps.jl")
    if isfile(PyCall_depsfile)
        include(PyCall_depsfile)
        println(pyprogramname)
        println(libpython)
    end
end

max_stars_repo_name: darrencl/pyjulia

max_issues_repo_name: darrencl/pyjulia

max_forks_repo_name: darrencl/pyjulia


text:
# This file was generated, do not modify it. # hide
train, test = partition(collect(eachindex(targets)), 0.70, shuffle=true,
                        rng=StableRNG(52))

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials