text:
# Variables #
# ========================== #
function JuMP.build_variable(_error::Function, variable::JuMP.ScalarVariable, set::DecisionSet)
    return VariableConstrainedOnCreation(variable, SingleDecisionSet(Decision(variable.info, Float64), set.constraint))
end

function JuMP.build_variable(_error::Function, variables::Vector{<:JuMP.ScalarVariable}, set::DecisionSet)
    return VariablesConstrainedOnCreation(variables, MultipleDecisionSet([Decision(variable.info, Float64) for variable in variables], set.constraint))
end

function JuMP.build_variable(_error::Function, variable::JuMP.ScalarVariable, ::KnownSet)
    return VariableConstrainedOnCreation(variable, SingleKnownSet(KnownDecision(variable.info, Float64)))
end

function JuMP.build_variable(_error::Function, variables::Vector{<:JuMP.ScalarVariable}, ::KnownSet)
    return VariablesConstrainedOnCreation(variables, MultipleKnownSet([KnownDecision(variable.info, Float64) for variable in variables]))
end

function JuMP.add_variable(model::Model, variable::VariableConstrainedOnCreation{<:Union{SingleDecisionSet, SingleKnownSet}}, name::String)
    decisions = get_decisions(model)
    if decisions isa IgnoreDecisions
        # Create a regular JuMP variable if decisions are not handled
        return JuMP.add_variable(model, variable.scalar_variable, name)
    end
    var_index, con_index = MOI.add_constrained_variable(backend(model), variable.set)
    # Map to model decisions after indices are known
    if !has_decision(decisions, var_index)
        # Store decision if is seen for the first time
        set_decision!(decisions, var_index, variable.set)
    else
        # Reuse if decision has been created already
        MOI.set(backend(model), MOI.ConstraintSet(), con_index, reuse(variable.set, decision(decisions, var_index)))
    end
    # Add any given decision constraints
    _moi_constrain_decision(backend(model), var_index, variable.scalar_variable.info, variable.set)
    # Finally, set any given name
    if !isempty(name)
        MOI.set(backend(model), MOI.VariableName(), var_index, name)
    end
    # Return created decision as DecisionRef/KnownRef
    return VariableRef(model, var_index, variable.set)
end

function JuMP.add_variable(model::Model, variable::VariablesConstrainedOnCreation{<:Union{MultipleDecisionSet, MultipleKnownSet}}, names)
    decisions = get_decisions(model)
    if decisions isa IgnoreDecisions
        # Create regular JuMP variables if decisions are not handled
        var_refs = map(zip(variable.scalar_variables, names)) do (scalar_variable, name)
            JuMP.add_variable(model, scalar_variable, name)
        end
        return reshape_vector(var_refs, variable.shape)
    end
    var_indices, con_index = MOI.add_constrained_variables(backend(model), variable.set)
    # Map to model decisions after indices are known
    seen_decisions = Vector{Decision{Float64}}()
    for (i,var_index) in enumerate(var_indices)
        if !has_decision(decisions, var_index)
            # Store decision if is seen for the first time
            set_decision!(decisions, var_index, i, variable.set)
        else
            # Reuse if decision has been created already
            push!(seen_decisions, decision(decisions, var_index))
        end
    end
    if !isempty(seen_decisions)
        # Sanity check
        length(seen_decisions) == length(variable.scalar_variables) || error("Inconsistency in number of seen decisions and created variables.")
        # Update decision set for reuse
        MOI.set(backend(model), MOI.ConstraintSet(), con_index, reuse(variable.set, seen_decisions))
    end
    # Add any given decision constraints
    for (index, scalar_variable) in zip(var_indices, variable.scalar_variables)
        _moi_constrain_decision(backend(model), index, scalar_variable.info, variable.set)
    end
    # Finally, set any given names
    for (var_index, name) in zip(var_indices, JuMP.vectorize(names, variable.shape))
        if !isempty(name)
            MOI.set(backend(model), MOI.VariableName(), var_index, name)
        end
    end
    # Return created decisions as DecisionRefs/KnownRefs
    refs = map(var_indices) do index
        VariableRef(model, index, variable.set)
    end
    return reshape_vector(refs, variable.shape)
end

# Containers #
# ========================== #
const DenseAxisArray = JuMP.Containers.DenseAxisArray
const SparseAxisArray = JuMP.Containers.SparseAxisArray

struct DecisionDenseAxisArray{V <: Union{DecisionRef, KnownRef}, A <: DenseAxisArray, S <: MOI.AbstractSet}
    array::A
    constraint::S

    function DecisionDenseAxisArray{V}(array::DenseAxisArray, constraint::MOI.AbstractSet) where V <: Union{DecisionRef, KnownRef}
        A = typeof(array)
        S = typeof(constraint)
        return new{V, A, S}(array, constraint)
    end
end

struct DecisionSparseAxisArray{V <: Union{DecisionRef, KnownRef}, A <: SparseAxisArray, S <: MOI.AbstractSet}
    array::A
    constraint::S

    function DecisionSparseAxisArray{V}(array::SparseAxisArray, constraint::MOI.AbstractSet) where V <: Union{DecisionRef, KnownRef}
        A = typeof(array)
        S = typeof(constraint)
        return new{V, A, S}(array, constraint)
    end
end

function JuMP.build_variable(_error::Function, variables::DenseAxisArray{<:JuMP.ScalarVariable}, set::DecisionSet)
    return DecisionDenseAxisArray{DecisionRef}(variables, set.constraint)
end

function JuMP.build_variable(_error::Function, variables::DenseAxisArray{<:JuMP.ScalarVariable}, ::KnownSet)
    return DecisionDenseAxisArray{KnownRef}(variables, NoSpecifiedConstraint())
end

function JuMP.build_variable(_error::Function, variables::SparseAxisArray{<:JuMP.ScalarVariable}, set::DecisionSet)
    return DecisionSparseAxisArray{DecisionRef}(variables, set.constraint)
end

function JuMP.build_variable(_error::Function, variables::SparseAxisArray{<:JuMP.ScalarVariable}, ::KnownSet)
    return DecisionSparseAxisArray{KnownRef}(variables, NoSpecifiedConstraint())
end

function JuMP.add_variable(model::Model, variable::DecisionDenseAxisArray{V}, names::DenseAxisArray{String}) where V <: Union{DecisionRef, KnownRef}
    array = variable.array
    refs = DenseAxisArray{V}(undef, array.axes...)
    for idx in eachindex(array)
        var = array[idx]
        refs[idx] = JuMP.add_variable(model,
                                      VariableConstrainedOnCreation(var, set(var, V, variable.constraint)),
                                      names[idx])
    end
    return refs
end

function JuMP.add_variable(model::Model,
                           variable::DecisionSparseAxisArray{V, SparseAxisArray{T,N,K}},
                           names::SparseAxisArray{String}) where {V <: Union{DecisionRef, KnownRef}, T, N, K}
    refs = SparseAxisArray(Dict{K,V}())
    for (idx, var) in variable.array.data
        refs[idx] = JuMP.add_variable(model,
                                      VariableConstrainedOnCreation(var, set(var, V, variable.constraint)),
                                      names[idx])
    end
    return refs
end

# Constraints #
# ========================== #
function JuMP.build_constraint(_error::Function, aff::Union{DecisionAffExpr, DecisionQuadExpr}, set::MOI.AbstractScalarSet)
    offset = constant(aff.variables)
    JuMP.add_to_expression!(aff.variables, -offset)
    shifted_set = MOIU.shift_constant(set, -offset)
    return JuMP.ScalarConstraint(aff, shifted_set)
end

function JuMP.build_constraint(_error::Function, aff::Union{DecisionAffExpr, DecisionQuadExpr}, lb, ub)
    JuMP.build_constraint(_error, aff, MOI.Interval(lb, ub))
end

# Helper function #
# ========================== #
function _moi_constrain_decision(backend::MOI.ModelLike, index, info, set::Union{SingleDecisionSet, MultipleDecisionSet})
    # We don't call the _moi* versions (e.g., _moi_set_lower_bound) because they
    # have extra checks that are not necessary for newly created variables.
    nothing_added = true
    if info.has_lb
        MOI.add_constraint(backend, SingleDecision(index),
                           MOI.GreaterThan{Float64}(info.lower_bound))
        nothing_added &= false
    end
    if info.has_ub
        MOI.add_constraint(backend, SingleDecision(index),
                           MOI.LessThan{Float64}(info.upper_bound))
        nothing_added &= false
    end
    if info.has_fix
        MOI.add_constraint(backend, SingleDecision(index),
                           MOI.EqualTo{Float64}(info.fixed_value))
        nothing_added &= false
    end
    if info.binary
        MOI.add_constraint(backend, SingleDecision(index),
                           MOI.ZeroOne())
        nothing_added &= false
    end
    if info.integer
        MOI.add_constraint(backend, SingleDecision(index), MOI.Integer())
        nothing_added &= false
    end
    if info.has_start
        MOI.set(backend, MOI.VariablePrimalStart(), index,
                Float64(info.start))
        nothing_added &= false
    end
    # Mark decision as free if not bounded or otherwise constrained
    if nothing_added && set.constraint isa NoSpecifiedConstraint
        MOI.add_constraint(backend, SingleDecision(index), FreeDecision())
    end
end

function _moi_constrain_decision(backend::MOI.ModelLike, index, info, ::Union{SingleKnownSet, MultipleKnownSet})
    # Known decisions are not constrained, just return
    return nothing
end


text:
# Backend: a thing upon with other things are drawn.

export finish, absolute_measure

abstract Backend

# Conversion to SimpleMeasure{MillimeterUnit} given the parent canvases
# bounding box in absolute units and its unit box.
function absolute_measure(u::SimpleMeasure{CanvasXUnit},
                          unit_box::BoundingBox, parent_box::BoundingBox)
    ((u.value - unit_box.x0.value) / unit_box.width.value) * parent_box.width
end

function absolute_measure(u::SimpleMeasure{CanvasYUnit},
                          unit_box::BoundingBox, parent_box::BoundingBox)
    ((u.value - unit_box.y0.value) / unit_box.height.value) * parent_box.height
end

function absolute_measure(u::SimpleMeasure{WidthUnit},
                          unit_box::BoundingBox, parent_box::BoundingBox)
    u.value * parent_box.width
end

function absolute_measure(u::SimpleMeasure{HeightUnit},
                          unit_box::BoundingBox, parent_box::BoundingBox)
    u.value * parent_box.height
end

function absolute_measure(u::CompoundMeasure,
                          unit_box::BoundingBox,
                          parent_box::BoundingBox)
    v = 0mm
    for (U, value) in u.values
        v += adsolute_measure(SimpleMeasure{U}(value), unit_box, parent_box)
    end
    v
end

function absolute_measure(point::Point,
                          unit_box::BoundingBox,
                          parent_box::BoundingBox)

    x = parent_box.x0 + absolute_measure(point.x, unit_box, parent_box)
    y = parent_box.y0 + absolute_measure(point.y, unit_box, parent_box)
    Point(x, y)
end

function absolute_measure(a::BoundingBox,
                          unit_box::BoundingBox,
                          parent_box::BoundingBox)
    BoundingBox(absolute_measure(unit_box, parent_box, a.x0),
                absolute_measure(unit_box, parent_box, a.y0),
                absolute_measure(unit_box, parent_box, a.width),
                absolute_measure(unit_box, parent_box, a.height))
end

# Conversion to native units
# The backend is responsible for defining a
# `native_measure{T <: AbsoluteUnit}(u::SimpleMeasure{T}, backend::Backend)`
# function to convert any SimpleMeasure in absolute units to a Measure of its
# particular subtype of NativeUnit. Everything else is handled here.

#function natie_measureu

function native_measure(u::SimpleMeasure{CanvasXUnit},
                        unit_box::BoundingBox,
                        parent_box::NativeBoundingBox,
                        backend::Backend)
    ((u.value - unit_box.x0.value) / unit_box.width.value) * parent_box.width
end

function native_measure(u::SimpleMeasure{CanvasYUnit},
                        unit_box::BoundingBox,
                        parent_box::NativeBoundingBox,
                        backend::Backend)
    ((u.value - unit_box.y0.value) / unit_box.height.value) * parent_box.height
end

function native_measure(u::SimpleMeasure{WidthUnit},
                        unit_box::BoundingBox,
                        parent_box::NativeBoundingBox,
                        backend::Backend)
    u.value * parent_box.width
end

function native_measure(u::SimpleMeasure{HeightUnit},
                        unit_box::BoundingBox,
                        parent_box::NativeBoundingBox,
                        backend::Backend)
    u.value * parent_box.height
end

function native_measure(u::SimpleMeasure,
                        rot::Rotation,
                        unit_box::BoundingBox,
                        parent_box::NativeBoundingBox,
                        backend::Backend)
    native_measure(u, unit_box, parent_box, backenD)
end

function native_measure(u::SimpleMeasure,
                        unit_box::BoundingBox,
                        parent_box::NativeBoundingBox,
                        backend::Backend)
    native_measure(u, backend)
end

function native_measure(u::SimpleMeasure,
                        t::NativeTransform,
                        unit_box::BoundingBox,
                        parent_box::NativeBoundingBox,
                        backend::Backend)
    native_measure(u, unit_box, parent_box, backend)
end

function native_measure(u::CompoundMeasure,
                        t::NativeTransform,
                        unit_box::BoundingBox,
                        parent_box::NativeBoundingBox,
                        backend::Backend)
    v = native_zero(backend)
    for (U, value) in u.values
        v += native_measure(SimpleMeasure{U}(value), t, unit_box,
                            parent_box, backend)
    end
    v
end

function native_measure(point::Point,
                        t::NativeTransform,
                        unit_box::BoundingBox,
                        parent_box::NativeBoundingBox,
                        backend::Backend)

    if typeof(point.x) <: NativeMeasure
        x = point.x
    else
        x = parent_box.x0 + native_measure(point.x, t, unit_box,
                                           parent_box, backend)
    end

    if typeof(point.y) <: NativeMeasure
        y = point.y
    else
        y = parent_box.y0 + native_measure(point.y, t, unit_box,
                                           parent_box, backend)
    end

    xy = [convert(Float64, x), convert(Float64, y), 1.0]
    xyt = t.M * xy

    Point(convert(typeof(x), xyt[1]),
          convert(typeof(y), xyt[2]))
end

function native_measure(rot::Rotation,
                        t::NativeTransform,
                        unit_box::BoundingBox,
                        parent_box::NativeBoundingBox,
                        backend::Backend)
    off = native_measure(rot.offset, t, unit_box, parent_box, backend)
    ct = cos(rot.theta)
    st = sin(rot.theta)
    x0 = off.x - (ct * off.x - st * off.y)
    y0 = off.y - (st * off.x + ct * off.y)
    NativeTransform([ct -st convert(Float64, x0)
                     st  ct convert(Float64, y0)
                     0.0 0.0 1.0])
end

function native_measure(box::BoundingBox,
                        t::NativeTransform,
                        unit_box::BoundingBox,
                        parent_box::NativeBoundingBox,
                        backend::Backend)
    NativeBoundingBox(
          parent_box.x0 + native_measure(box.x0, t, unit_box,
                                         parent_box, backend),
          parent_box.y0 + native_measure(box.y0, t, unit_box,
                                         parent_box, backend),
          native_measure(box.width,  t, unit_box, parent_box, backend),
          native_measure(box.height, t, unit_box, parent_box, backend))
end

ext: jl


text:
module SoundSource

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Sound Source" => Ahorn.EntityPlacement(
        Maple.SoundSource
    )
)

function Ahorn.selection(entity::Maple.SoundSource)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

Ahorn.editingOptions(entity::Maple.SoundSource) = Dict{String, Any}(
    "sound" => EnvironmentSounds.sounds
)

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.SoundSource, room::Maple.Room) = Ahorn.drawImage(ctx, Ahorn.Assets.speaker, -12, -12)

end


text:
@doc raw"""
    Triangle2

An object `t` of the type [`Triangle2`](@ref) is a triangle in the
two-dimensional Euclidean plane ``Ε^2``.

Triangle `t` is oriented, i.e., its boundary has clockwise or
counterclockwise orientation. We call the side to the left of the boundary the
positive side and the side to the right of the boundary the negative side.

The boundary of a triangle splits the plane in two open regions, a bounded one
and an unbounded one.
"""
Triangle2

"""
    Triangle2(p::Point2, q::Point2, r::Point2)

Introduces a triangle `t` with vertices `p`, `q` and `r`.
"""
Triangle2(p::Point2, q::Point2, r::Point2)

"""
    ==(t::Triangle2, t2::Triangle2)

Test for equality: two triangles are equal, iff there exists a cyclic
permutation of the vertices of ``t2``, such that they are equal to the vertices
of `t`.
"""
==(t::Triangle2, t2::Triangle2)

"""
    vertex(t::Triangle2, i::Integer)

Returns the `i`ᵗʰ vertex modulo 3 of `t`.
"""
vertex(t::Triangle2, i::Integer)

"""
    is_degenerate(t::Triangle2)

Triangle `t` is degenerate, if the vertices are collinear.
"""
is_degenerate(t::Triangle2)

"""
    orientation(t::Triangle2)

Returns the orientation of `t`.
"""
orientation(t::Triangle2)

"""
    oriented_side(t::Triangle2, p::Point2)

Returns [`ON_ORIENTED_BOUNDARY`](@ref), or [`POSITIVE_SIDE`](@ref), or the
constant [`ON_NEGATIVE_SIDE`](@ref), determined by the position of point `p`.

!!! info "Precondition"

    `t` is not degenerate.
"""
oriented_side(t::Triangle2, p::Point2)

"""
    bounded_side(t::Triangle2, p::Point2)

Returns [`ON_BOUNDARY`](@ref), [`ON_BOUNDED_SIDE`](@ref), or else
[`ON_UNBOUNDED_SIDE`](@ref), depending on where point `p` is.

!!! info "Precondition"

    `t` is not degenerate.
"""
bounded_side(t::Triangle2, p::Point2)

"""
    has_on_positive_side(t::Triangle2, p::Point2)

Returns `true`, iff `p` lies on the positive side of `t`.
"""
has_on_positive_side(t::Triangle2, p::Point2)

"""
    has_on_negative_side(t::Triangle2, p::Point2)

Returns `true`, iff `p` lies on the negative side of `t`.
"""
has_on_negative_side(t::Triangle2, p::Point2)

"""
    has_on_boundary(t::Triangle2, p::Point2)

Returns `true`, iff `p` lies on the boundary of `t`.
"""
has_on_boundary(t::Triangle2, p::Point2)

"""
    has_on_bounded_side(t::Triangle2, p::Point2)

Returns `true`, iff `p` lies properly inside of `t`.
"""
has_on_bounded_side(t::Triangle2, p::Point2)

"""
    has_on_unbounded_side(t::Triangle2, p::Point2)

Returns `true`, iff `p` lies properly outside of `t`.
"""
has_on_unbounded_side(t::Triangle2, p::Point2)

"""
    opposite(t::Triangle2)

Returns a triangle where the boundary is oriented the other way round (this
flips the positive and negative side, but not the bounded and unbounded side).
"""
opposite(t::Triangle2)

"""
    area(t::Triangle2)

Returns the signed area of `t`.
"""
area(t::Triangle2)

"""
    bbox(t::Triangle2)

Returns a bounding box containing `t`.
"""
bbox(t::Triangle2)

"""
    transform(t::Triangle2, at::AffTransformation2)

Returns the triangle obtained by applying `at` on the three vertices of `t`.
"""
transform(t::Triangle2, at::AffTransformation2)


text:
"""
    PointState

Holds the state variables for a point

# Fields:
 - `u`: Displacement variables for the point (in the global coordinate frame)
 - `theta`: Wiener-Milenkovic rotational displacement variables for the point
 - `F`: Externally applied forces on the point (in the global coordinate frame)
 - `M`: Externally applied moments on the point (in the global coordinate frame)
"""
struct PointState{TF}
    u::SVector{3, TF}
    theta::SVector{3, TF}
    F::SVector{3, TF}
    M::SVector{3, TF}
end

"""
    ElementState

Holds the state variables for an element

# Fields:
 - `u`: Displacement variables for the element (in the global coordinate frame)
 - theta: Wiener-Milenkovic rotational displacement variables for the element
       (in the global coordinate frame)
 - `F`: Resultant forces for the element (in the deformed local coordinate frame)
 - `M`: Resultant moments for the element (in the deformed local coordinate frame)
 - `V`: Linear velocity of the element (in the deformed local coordinate frame)
 - `Omega`: Angular velocity of the element (in the deformed local coordinate frame)
"""
struct ElementState{TF}
    u::SVector{3, TF}
    theta::SVector{3, TF}
    F::SVector{3, TF}
    M::SVector{3, TF}
    V::SVector{3, TF}
    Omega::SVector{3, TF}
end

"""
    AssemblyState{TF, TP<:AbstractVector{PointState{TF}},
        TE<:AbstractVector{ElementState{TF}}}

Struct for storing state variables for the points and elements in an assembly.

# Fields:
 - `points::TP`: Array of `PointState`s for each point in the assembly
 - `elements::TE`: Array of `ElementState`s for each element in the assembly
"""
struct AssemblyState{TF, TP<:AbstractVector{PointState{TF}}, TE<:AbstractVector{ElementState{TF}}}
    points::TP
    elements::TE
end
Base.eltype(::AssemblyState{TF, TP, TE}) where {TF, TP, TE} = TF

"""
    AssemblyState(system, assembly, x = system.x;
        prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}())

Post-process the system state given the solution vector `x`.  Return an object
of type `AssemblyState` that defines the state of the assembly for the time step.

If `prescribed_conditions` is not provided, all point state variables are assumed
to be displacements/rotations, rather than their actual identities as used in the
analysis.
"""
function AssemblyState(system, assembly, x = system.x;
    prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}())

    points = extract_point_states(system, assembly, x; prescribed_conditions)

    elements = extract_element_states(system, x)

    return AssemblyState(points, elements)
end

"""
    extract_point_state(system, assembly, ipoint, x = system.x;
        prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}())

Return the state variables corresponding to point `ipoint` (see [`PointState`](@ref))
given the solution vector `x`.

If `prescribed_conditions` is not provided, all point state variables are assumed
to be displacements/rotations, rather than their actual identities as used in the
analysis.
"""
function extract_point_state(system, assembly, ipoint, x = system.x;
    prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}())

    pc = typeof(prescribed_conditions) <: AbstractDict ? prescribed_conditions :
        prescribed_conditions(system.t)

    force_scaling = system.force_scaling

    icol_p = system.icol_point[ipoint]

    if icol_p <= 0
        # point variables are not state variables, solve for point variables

        # find the first beam connected to the point, checking both sides
        ielem = findfirst(x -> x == ipoint, assembly.start)
        if !isnothing(ielem)
            side = -1
        else
            ielem = findfirst(x -> x == ipoint, assembly.stop)
            side = 1
        end

        elem = assembly.elements[ielem]

        # get beam state variables
        icol_b = system.icol_elem[ielem]
        u_e = SVector(x[icol_b   ], x[icol_b+1 ], x[icol_b+2 ])
        θ_b = SVector(x[icol_b+3 ], x[icol_b+4 ], x[icol_b+5 ])
        F_e = SVector(x[icol_b+6 ], x[icol_b+7 ], x[icol_b+8 ]) .* force_scaling
        M_e = SVector(x[icol_b+9 ], x[icol_b+10], x[icol_b+11]) .* force_scaling

        # get beam element properties
        ΔL_b = elem.L
        Ct_b = get_C(θ_b)'
        Cab_b = elem.Cab
        γ_b = element_strain(elem, F_e, M_e)
        κ_b = element_curvature(elem, F_e, M_e)

        # solve for the displacements/rotations of the point
        u = u_e + side*ΔL_b/2*(Ct_b*Cab_b*(e1 + γ_b) - Cab_b*e1)
        theta = θ_b + side*ΔL_b/2*get_Qinv(θ_b)*Cab_b*κ_b
        F = zero(u)
        M = zero(u)
    else
        # point variables are state variables, extract point state variables
        prescribed = haskey(pc, ipoint)
        if prescribed
            u, theta, F, M = point_variables(x, icol_p, pc[ipoint], force_scaling)
        else
            u, theta, F, M = point_variables(x, icol_p)
        end
    end

    # convert rotation parameter to Wiener-Milenkovic parameters
    scaling = rotation_parameter_scaling(theta)
    theta *= scaling

    # promote all variables to the same type
    u, theta, F, M = promote(u, theta, F, M)

    return PointState(u, theta, F, M)
end

"""
    extract_point_states(system, assembly, x = system.x;
        prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}())

Return the state variables corresponding to each point (see [`PointState`](@ref))
given the solution vector `x`.

If `prescribed_conditions` is not provided, all point state variables are assumed
to be displacements/rotations, rather than their actual identities as used in the
analysis.
"""
function extract_point_states(system, assembly, x = system.x; kwargs...)
    TF = promote_type(eltype(system), eltype(x))
    points = Vector{PointState{TF}}(undef, length(assembly.points))
    return extract_point_states!(points, system, assembly, x; kwargs...)
end

"""
    extract_point_states!(points, system, assembly, x = system.x;
        prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}())

Pre-allocated version of [`extract_point_states`](@ref)
"""
function extract_point_states!(points, system, assembly, x = system.x;
    prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}())

    for ipoint = 1:length(points)
        points[ipoint] = extract_point_state(system, assembly, ipoint, x; prescribed_conditions)
    end

    return points
end

"""
    extract_element_state(system, ielem, x = system.x)

Return the state variables corresponding to element `ielem` (see [`ElementState`](@ref))
given the solution vector `x`.
"""
function extract_element_state(system, ielem, x = system.x)

    force_scaling = system.force_scaling

    icol = system.icol_elem[ielem]
    static = system.irow_elem[ielem] <= 0
    u = SVector(x[icol], x[icol+1], x[icol+2])
    theta = SVector(x[icol+3], x[icol+4], x[icol+5])
    F = SVector(x[icol+6], x[icol+7], x[icol+8]) .* force_scaling
    M = SVector(x[icol+9], x[icol+10], x[icol+11]) .* force_scaling
    V = ifelse(static, zero(u), SVector(x[icol+12], x[icol+13], x[icol+14]))
    Ω = ifelse(static, zero(u), SVector(x[icol+15], x[icol+16], x[icol+17]))

    # convert rotation parameter to Wiener-Milenkovic parameters
    scaling = rotation_parameter_scaling(theta)
    theta *= scaling

    # promote all variables to the same type
    u, theta, F, M, V, Ω = promote(u, theta, F, M, V, Ω)

    return ElementState(u, theta, F, M, V, Ω)
end

"""
    extract_element_states(system, x = system.x)

Return the state variables corresponding to each element (see [`ElementState`](@ref))
given the solution vector `x`.
"""
function extract_element_states(system, x = system.x)
    TF = promote_type(eltype(system), eltype(x))
    elements = Vector{ElementState{TF}}(undef, length(system.icol_elem))
    return extract_element_states!(elements, system, x)
end

"""
    extract_element_states!(elements, system, x = system.x)

Pre-allocated version of [`extract_element_states`](@ref)
"""
function extract_element_states!(elements, system, x = system.x)
    for ielem = 1:length(elements)
        elements[ielem] = extract_element_state(system, ielem, x)
    end
    return elements
end

"""
    rotate(xyz, r, theta)

Rotate the vectors in `xyz` about point `r` using the Wiener-Milenkovic
parameters in `theta`.
"""
rotate(xyz, r, theta) = rotate!(copy(xyz), r, theta)

"""
    rotate!(xyz, r, theta)

Pre-allocated version of [`rotate`](@ref)
"""
function rotate!(xyz, r, theta)
    # reshape cross section points (if necessary)
    rxyz = reshape(xyz, 3, :)
    # convert inputs to static arrays
    r = SVector{3}(r)
    theta = SVector{3}(theta)
    # create rotation matrix
    Ct = wiener_milenkovic(theta)'
    # rotate each point
    for ipt = 1:size(rxyz, 2)
        p = SVector(rxyz[1,ipt], rxyz[2,ipt], rxyz[3,ipt])
        rxyz[:,ipt] .= Ct*(p - r) + r
    end
    # return the result
    return xyz
end

"""
    translate(xyz, u)

Translate the points in `xyz` by the displacements in `u`.
"""
translate(xyz, u) = translate!(copy(xyz), u)

"""
    translate!(xyz, u)

Pre-allocated version of [`translate`](@ref)
"""
function translate!(xyz, u)
    # reshape cross section points (if necessary)
    rxyz = reshape(xyz, 3, :)
    # convert inputs to static arrays
    u = SVector{3}(u)
    # translate each point
    for ipt = 1:size(rxyz, 2)
        p = SVector(rxyz[1,ipt], rxyz[2,ipt], rxyz[3,ipt])
        rxyz[:,ipt] .= p .+ u
    end
    # return the result
    return xyz
end

"""
    deform_cross_section(xyz, r, u, theta)

Rotate the points in `xyz` (of shape (3, :)) about point `r` using
the Wiener-Milenkovic parameters in `theta`, then translate the points by the
displacements in `u`.
"""
deform_cross_section(xyz, r, u, theta) = deform_cross_section!(copy(xyz), r, u, theta)

"""
    deform_cross_section!(xyz, r, u, theta)

Pre-allocated version of [`deform_cross_section`](@ref)
"""
deform_cross_section!(xyz, r, u, theta) = translate!(rotate!(xyz, r, theta), u)

"""
    cross_section_velocities(xyz, r, v, ω)

Calculate the velocities of the points in `xyz` given the linear velocity `v`,
and the angular velocity `ω` about point `r`.
"""
function cross_section_velocities(xyz, r, v, ω)
    TF = promote_type(eltype(xyz), eltype(r), eltype(v), eltype(ω))
    vxyz = similar(xyz, TF)
    return cross_section_velocities!(vxyz, xyz, r, v, ω)
end

"""
    cross_section_velocities!(vxyz, xyz, r, v, ω)

Pre-allocated version of [`cross_section_velocities`](@ref)
"""
function cross_section_velocities!(vxyz, xyz, r, v, ω)
    # reshape cross section points (if necessary)
    rxyz = reshape(xyz, 3, :)
    rvxyz = reshape(vxyz, 3, :)
    # convert inputs to static arrays
    r = SVector{3}(r)
    v = SVector{3}(v)
    ω = SVector{3}(ω)
    # calculate velocities at each point
    for ipt = 1:size(rxyz, 2)
        p = SVector(rxyz[1,ipt], rxyz[2,ipt], rxyz[3,ipt])
        rvxyz[:,ipt] .= v .+ cross(ω, p - r)
    end
    # return the result
    return vxyz
end

"""
    write_vtk(name, assembly::Assembly; kwargs...)
    write_vtk(name, assembly::Assembly, state::AssemblyState; kwargs...)
    write_vtk(name, assembly::Assembly, history::Vector{<:AssemblyState}], dt;
        kwargs...)

Write the deformed geometry (and associated data) to a VTK file for visualization
using ParaView.

The `state` argument may be omitted to write the original geometry to a VTK file
without any associated data.

If the solution time `history` is provided, the time step must also be provided

# Keyword Arguments
 - `sections = nothing`: Cross section geometry corresponding to each point,
    defined in a frame aligned with the global frame but centered around the
    corresponding point. Defined as an array with shape `(3, ncross, np)` where `ncross`
    is the number of points in each cross section and `np` is the number of points.
 - `scaling=1.0`: Parameter to scale the deflections (only valid if state is provided)
 - `metadata=Dict()`: Dictionary of metadata for the file(s)
"""
function write_vtk(name, assembly; sections=nothing, metadata=Dict())

    # get problem dimensions
    npoint = length(assembly.points)
    ncross = isnothing(sections) ? 1 : size(sections, 2)
    nelem = length(assembly.elements)

    if isnothing(sections)
        # extract point locations
        points = Matrix{eltype(assembly)}(undef, 3, npoint)
        for ip = 1:npoint
            for i = 1:3
                points[i,ip] = assembly.points[ip][i]
            end
        end

        # create cells
        cells = [MeshCell(PolyData.Lines(), [assembly.start[i], assembly.stop[i]]) for i = 1:nelem]

    else

        li = LinearIndices((ncross, npoint))

        # extract cell point locations
        points = Matrix{eltype(assembly)}(undef, 3, ncross*npoint)
        if ndims(sections) > 2
            for ip = 1:npoint
                for ic = 1:ncross
                    points[:,li[ic,ip]] = assembly.points[ip] + sections[:,ic,ip]
                end
            end
        else
            for ip = 1:npoint
                for ic = 1:ncross
                    points[:,li[ic,ip]] = assembly.points[ip] + sections[:,ic]
                end
            end
        end

        # construct triangle strip for each beam element
        cells = Vector{MeshCell{VTKCellType, Vector{Int64}}}(undef, nelem)
        for ielem = 1:nelem
            # index of key point corresponding to the start of the beam element
            ipt1 = assembly.start[ielem]
            # index of key point corresponding to the end of the beam element
            ipt2 = assembly.stop[ielem]
            # triangle strip points
            connectivity = Vector{Int}(undef, ncross*2)
            for ic = 1:ncross
                connectivity[2*ic-1] = li[ic, ipt1]
                connectivity[2*ic] = li[ic, ipt2]
            end
            cells[ielem] = MeshCell(VTKCellTypes.VTK_TRIANGLE_STRIP, connectivity)
        end

    end

    # write vtk file
    vtk_grid(name, points, cells) do vtkfile

        # add metadata
        for (key, value) in pairs(metadata)
            vtkfile[string(key)] = value
        end

        # add local axis data
        axis_name = ["x-axis", "y-axis", "z-axis"]
        axis_vector = [e1, e2, e3]
        for i = 1:3
            data = Matrix{eltype(assembly)}(undef, 3, nelem)
            for ielem = 1:nelem
                data[:, ielem] .= assembly.elements[ielem].Cab*axis_vector[i]
            end
            vtkfile[axis_name[i], VTKCellData()] = data
        end
    end

    return nothing
end

function write_vtk(name, assembly, state; sections = nothing, scaling=1.0,
    metadata=Dict())

    # get problem dimensions
    npoint = length(assembly.points)
    ncross = isnothing(sections) ? 1 : size(sections, 2)
    nelem = length(assembly.elements)

    if isnothing(sections)
        # extract point locations
        points = Matrix{eltype(assembly)}(undef, 3, npoint)
        for ip = 1:npoint
            for i = 1:3
                points[i,ip] = assembly.points[ip][i] + scaling*state.points[ip].u[i]
            end
        end

        # create cells
        cells = [MeshCell(PolyData.Lines(), [assembly.start[i], assembly.stop[i]]) for i = 1:nelem]
    else
        li = LinearIndices((ncross, npoint))

        # extract cell point locations
        points = Matrix{eltype(assembly)}(undef, 3, ncross*npoint)
        if ndims(sections) > 2
            for ip = 1:npoint
                for ic = 1:ncross
                    u = scaling*state.points[ip].u
                    c = scaling*state.points[ip].theta
                    Ct = wiener_milenkovic(c)'
                    points[:,li[ic,ip]] = assembly.points[ip] + u + Ct*sections[:,ic,ip]
                end
            end
        else
            for ip = 1:npoint
                for ic = 1:ncross
                    u = scaling*state.points[ip].u
                    c = scaling*state.points[ip].theta
                    Ct = wiener_milenkovic(c)'
                    points[:,li[ic,ip]] = assembly.points[ip] + u + Ct*sections[:,ic]
                end
            end
        end

        # construct triangle strip for each beam element
        cells = Vector{MeshCell{VTKCellType, Vector{Int64}}}(undef, nelem)
        for ielem = 1:nelem
            # index of key point corresponding to the start of the beam element
            ipt1 = assembly.start[ielem]
            # index of key point corresponding to the end of the beam element
            ipt2 = assembly.stop[ielem]
            # triangle strip points
            connectivity = Vector{Int}(undef, ncross*2)
            for ic = 1:ncross
                connectivity[2*ic-1] = li[ic, ipt1]
                connectivity[2*ic] = li[ic, ipt2]
            end
            cells[ielem] = MeshCell(VTKCellTypes.VTK_TRIANGLE_STRIP, connectivity)
        end
    end

    # write vtk file
    vtk_grid(name, points, cells) do vtkfile

        vtkfile["scaling"] = scaling

        # add metadata
        for (key, value) in pairs(metadata)
            vtkfile[string(key)] = value
        end

        # add local axis data
        axis_name = ["x-axis", "y-axis", "z-axis"]
        axis_vector = [e1, e2, e3]
        for i = 1:3
            data = Matrix{eltype(assembly)}(undef, 3, nelem)
            for ielem = 1:nelem
                CtCab = wiener_milenkovic(state.elements[ielem].theta)'*assembly.elements[ielem].Cab
                data[:, ielem] .= CtCab * axis_vector[i]
            end
            vtkfile[axis_name[i], VTKCellData()] = data
        end

        # add point data
        for field in fieldnames(PointState)
            data = Matrix{eltype(assembly)}(undef, 3, npoint*ncross)
            li = LinearIndices((ncross, npoint))
            for ip = 1:npoint
                data[:, li[:,ip]] .= getproperty(state.points[ip], field)
            end
            vtkfile[string(field), VTKPointData()] = data
        end

        # add cell data
        for field in fieldnames(ElementState)
            data = Matrix{eltype(assembly)}(undef, 3, nelem)
            for ielem = 1:nelem
                data[:, ielem] .= getproperty(state.elements[ielem], field)
            end
            vtkfile[string(field), VTKCellData()] = data
        end

    end

    return nothing
end

function write_vtk(name, assembly, history, t; sections=nothing, scaling=1.0,
    metadata=Dict())

    # get problem dimensions
    npoint = length(assembly.points)
    ncross = isnothing(sections) ? 1 : size(sections, 2)
    nelem = length(assembly.elements)

    paraview_collection(name) do pvd

        for (current_step, state) in enumerate(history)

            if isnothing(sections)
                # extract point locations
                points = Matrix{eltype(assembly)}(undef, 3, npoint)
                for ip = 1:npoint
                    for i = 1:3
                        points[i,ip] = assembly.points[ip][i] + scaling*state.points[ip].u[i]
                    end
                end

                # create cells
                cells = [MeshCell(PolyData.Lines(), [assembly.start[i], assembly.stop[i]]) for i = 1:nelem]
            else

                li = LinearIndices((ncross, npoint))

                # extract cell point locations
                points = Matrix{eltype(assembly)}(undef, 3, ncross*npoint)
                if ndims(sections) > 2
                    for ip = 1:npoint
                        for ic = 1:ncross
                            u = scaling*state.points[ip].u
                            c = scaling*state.points[ip].theta
                            Ct = wiener_milenkovic(c)'
                            points[:,li[ic,ip]] = assembly.points[ip] + u + Ct*sections[:,ic,ip]
                        end
                    end
                else
                    for ip = 1:npoint
                        for ic = 1:ncross
                            u = scaling*state.points[ip].u
                            c = scaling*state.points[ip].theta
                            Ct = wiener_milenkovic(c)'
                            points[:,li[ic,ip]] = assembly.points[ip] + u + Ct*sections[:,ic]
                        end
                    end
                end

                # construct triangle strip for each beam element
                cells = Vector{MeshCell{VTKCellType, Vector{Int64}}}(undef, nelem)
                for ielem = 1:nelem
                    # index of key point corresponding to the start of the beam element
                    ipt1 = assembly.start[ielem]
                    # index of key point corresponding to the end of the beam element
                    ipt2 = assembly.stop[ielem]
                    # triangle strip points
                    connectivity = Vector{Int}(undef, ncross*2)
                    for ic = 1:ncross
                        connectivity[2*ic-1] = li[ic, ipt1]
                        connectivity[2*ic] = li[ic, ipt2]
                    end
                    cells[ielem] = MeshCell(VTKCellTypes.VTK_TRIANGLE_STRIP, connectivity)
                end
            end

            # write vtk file
            vtkfile = vtk_grid(name*"-step$current_step", points, cells)

            # add metadata
            vtkfile["scaling"] = scaling
            vtkfile["time"] = t[current_step]
            for (key, value) in pairs(metadata)
                vtkfile[string(key)] = value
            end

            # add local axis data
            axis_name = ["x-axis", "y-axis", "z-axis"]
            axis_vector = [e1, e2, e3]
            for i = 1:3
                data = Matrix{eltype(assembly)}(undef, 3, nelem)
                for ielem = 1:nelem
                    CtCab = wiener_milenkovic(state.elements[ielem].theta)' * assembly.elements[ielem].Cab
                    data[:, ielem] .= CtCab * axis_vector[i]
                end
                vtkfile[axis_name[i], VTKCellData()] = data
            end

            # add point data
            for field in fieldnames(PointState)
                data = Matrix{eltype(assembly)}(undef, 3, npoint*ncross)
                li = LinearIndices((ncross, npoint))
                for ip = 1:npoint
                    data[:, li[:,ip]] .= getproperty(state.points[ip], field)
                end
                vtkfile[string(field), VTKPointData()] = data
            end

            # add cell data
            for field in fieldnames(ElementState)
                data = Matrix{eltype(assembly)}(undef, 3, nelem)
                for ielem = 1:nelem
                    data[:, ielem] .= getproperty(state.elements[ielem], field)
                end
                vtkfile[string(field), VTKCellData()] = data
            end

            pvd[t[current_step]] = vtkfile
        end
    end

    return nothing
end

"""
    write_vtk(name, assembly::Assembly, [state::AssemblyState, ]λ::Number,
        eigenstate::AssemblyState; scaling=1.0, mode_scaling=1.0, cycles=1,
        steps=100)

Write a series of files corresponding to the elastic motion of the `assembly`
about the deformed state encoded in `state` defined by the eigenvalue `λ` and
the eigenvector encoded in `eigenstate` over the time period specified by `time`.

The steady-state deflections can be scaled with `scaling` and the eigenmode
deflections can be scaled using `mode_scaling`.

The current time is encoded in the metadata tag "time"
"""
function write_vtk(name, assembly, state, λ, eigenstate;
    sections = nothing,
    scaling=1.0,
    mode_scaling=1.0,
    cycles = 1,
    steps = 100)

    npoint = length(assembly.points)
    ncross = isnothing(sections) ? 1 : size(sections, 2)
    nelem = length(assembly.elements)

    damping = real(λ)
    period = 2*pi/imag(λ)

    if damping <= 0.0
        start = -period*cycles
        stop = 0.0
    else
        start = 0.0
        stop = period*cycles
    end

    time=range(start, stop, length=steps)

    paraview_collection(name) do pvd

        for (current_step, t) in enumerate(time)

            if isnothing(sections)

                # extract point locations
                points = Matrix{eltype(assembly)}(undef, 3, npoint)
                for ip = 1:npoint
                    for i = 1:3
                        points[i,ip] = assembly.points[ip][i] +
                            scaling*state.points[ip].u[i] +
                            mode_scaling*real(eigenstate.points[ip].u[i]*exp(λ*t))
                    end
                end

                # create cells
                cells = [MeshCell(PolyData.Lines(), [assembly.start[i], assembly.stop[i]]) for i = 1:nelem]

            else

                li = LinearIndices((ncross, npoint))

                # extract cell point locations
                points = Matrix{eltype(assembly)}(undef, 3, ncross*npoint)
                if ndims(sections) > 2
                    for ip = 1:npoint
                        for ic = 1:ncross
                            u = mode_scaling*real.(eigenstate.points[ip].u .* exp(λ*t))
                            c = mode_scaling*real.(eigenstate.points[ip].theta .* exp(λ*t))
                            Ct = wiener_milenkovic(c)'
                            points[:,li[ic,ip]] = assembly.points[ip] + u + Ct*sections[:,ic,ip]
                        end
                    end
                else
                    for ip = 1:npoint
                        for ic = 1:ncross
                            u = mode_scaling*real.(eigenstate.points[ip].u .* exp(λ*t))
                            c = mode_scaling*real.(eigenstate.points[ip].theta .* exp(λ*t))
                            Ct = wiener_milenkovic(c)'
                            points[:,li[ic,ip]] = assembly.points[ip] + u + Ct*sections[:,ic]
                        end
                    end
                end

                # construct triangle strip for each beam element
                cells = Vector{MeshCell{VTKCellType, Vector{Int64}}}(undef, nelem)
                for ielem = 1:nelem
                    # index of key point corresponding to the start of the beam element
                    ipt1 = assembly.start[ielem]
                    # index of key point corresponding to the end of the beam element
                    ipt2 = assembly.stop[ielem]
                    # triangle strip points
                    connectivity = Vector{Int}(undef, ncross*2)
                    for ic = 1:ncross
                        connectivity[2*ic-1] = li[ic, ipt1]
                        connectivity[2*ic] = li[ic, ipt2]
                    end
                    cells[ielem] = MeshCell(VTKCellTypes.VTK_TRIANGLE_STRIP, connectivity)
                end

            end

            # write vtk file
            vtkfile = vtk_grid(name*"-step$(current_step)", points, cells)

            # add metadata
            vtkfile["scaling"] = scaling
            vtkfile["time"] = t
            vtkfile["phase"] = 2*pi*t/period

            # add local axis data
            axis_name = ["x-axis", "y-axis", "z-axis"]
            axis_vector = [e1, e2, e3]
            for i = 1:3
                data = Matrix{eltype(assembly)}(undef, 3, nelem)
                for ielem = 1:nelem
                    CtCab = wiener_milenkovic(real(state.elements[ielem].theta .* exp(λ*t)))' *
                        assembly.elements[ielem].Cab
                    data[:, ielem] .= CtCab * axis_vector[i]
                end
                vtkfile[axis_name[i], VTKCellData()] = data
            end

            # add point data
            for field in fieldnames(PointState)
                data = Matrix{eltype(assembly)}(undef, 3, ncross*npoint)
                li = LinearIndices((ncross, npoint))
                for ip = 1:npoint
                    data[:, li[:,ip]] .= getproperty(state.points[ip], field) +
                         mode_scaling*real(getproperty(eigenstate.points[ip], field) .* exp(λ*t))
                end
                vtkfile[string(field), VTKPointData()] = data
            end

            # add cell data
            for field in fieldnames(ElementState)
                data = Matrix{eltype(assembly)}(undef, 3, nelem)
                for ielem = 1:nelem
                    data[:, ielem] .= getproperty(state.elements[ielem], field) +
                         mode_scaling*real(getproperty(eigenstate.elements[ielem], field).* exp(λ*t))
                end
                vtkfile[string(field), VTKCellData()] = data
            end

            pvd[t] = vtkfile
        end
    end

    return nothing
end

"""
    left_eigenvectors(system, λ, V)
    left_eigenvectors(K, M, λ, V)

Compute the left eigenvector matrix `U` for the `system` using inverse power
iteration given the eigenvalues `λ` and the corresponding right eigenvector
matrix `V`.

The complex conjugate of each left eigenvector is stored in each row of the
matrix `U`

Left and right eigenvectors satisfy the following M-orthogonality condition:
 - u'*M*v = 1 if u and v correspond to the same eigenvalue
 - u'*M*v = 0 if u and v correspond to different eigenvalues
This means that U*M*V = I

This function assumes that `system` has not been modified since the eigenvalues
and right eigenvectors were computed.
"""
left_eigenvectors(system, λ, V) = left_eigenvectors(system.K, system.M, λ, V)

function left_eigenvectors(K, M, λ, V)

    # problem type and dimensions
    TC = eltype(V)
    nx = size(V,1)
    nev = size(V,2)

    # allocate storage
    U = rand(TC, nev, nx)
    u = Vector{TC}(undef, nx)
    tmp = Vector{TC}(undef, nx)

    # # get entries in M
    # iM, jM, valM = findnz(M)

    # compute eigenvectors for each eigenvalue
    for iλ = 1:nev

        # factorize (K + λ*M)'
        KmλMfact = factorize(K' + λ[iλ]'*M')

        # initialize left eigenvector
        for i = 1:nx
            u[i] = U[iλ,i]
        end

        # perform a few iterations to converge the left eigenvector
        for ipass = 1:3
            # get updated u
            mul!(tmp, M, u)
            ldiv!(u, KmλMfact, tmp)
            # normalize u
            unorm = zero(TC)
            for i in axes(M, 1), j in axes(M, 2)
                unorm += conj(u[i])*M[i,j]*V[j,iλ]
            end
            rdiv!(u, conj(unorm))
        end

        # store conjugate of final eigenvector
        for i = 1:nx
            U[iλ,i] = conj(u[i])
        end
    end

    return U
end

function left_eigenvectors(K, M::SparseMatrixCSC, λ, V)

    # problem type and dimensions
    TC = eltype(V)
    nx = size(V,1)
    nev = size(V,2)

    # allocate storage
    U = rand(TC, nev, nx)
    u = Vector{TC}(undef, nx)
    tmp = Vector{TC}(undef, nx)

    # get entries in M
    iM, jM, valM = findnz(M)

    # compute eigenvectors for each eigenvalue
    for iλ = 1:nev

        # factorize (K + λ*M)'
        KmλMfact = factorize(K' + λ[iλ]'*M')

        # initialize left eigenvector
        for i = 1:nx
            u[i] = U[iλ,i]
        end

        # perform a few iterations to converge the left eigenvector
        for ipass = 1:3
            # get updated u
            mul!(tmp, M, u)
            ldiv!(u, KmλMfact, tmp)
            # normalize u
            unorm = zero(TC)
            for k = 1:length(valM)
                unorm += conj(u[iM[k]])*valM[k]*V[jM[k],iλ]
            end
            rdiv!(u, conj(unorm))
        end

        # store conjugate of final eigenvector
        for i = 1:nx
            U[iλ,i] = conj(u[i])
        end
    end

    return U
end

"""
    correlate_eigenmodes(C)

Return the permutation and the associated corruption index vector which associates
eigenmodes from the current iteration with those of the previous iteration given
the correlation matrix `C`.

The correlation matrix can take one of the following forms (in order of preference):
 - `C = U_p*M*V`
 - `C = U*M_p*V_p`
 - `C = V_p'*V`
 - `C = V'*V_p`
where `U` is a matrix of conjugated left eigenvectors, `M` is the system mass
matrix, `V` is a matrix of right eigenvectors, and `()_p` indicates a variable
from the previous iteration.

Note that the following two forms of the correlation matrix seem to be significantly
inferior to their counterparts listed above: `C = U*M*V_p` and `C = U_p*M_p*V`.
This is likely due to the way in which the left eigenvector matrix is calculated.

The corruption index is the largest magnitude in a given row of `C`
that was not chosen divided by the magnitude of the chosen eigenmode.  It is most
meaningful when using one of the forms of the correlation matrix that uses left
eigenvectors since correct eigenmodes will have magnitudes close to 1 and
incorrect eigenmodes will have magnitudes close to 0.

If the new mode number is already assigned, the next highest unassigned mode
number is used.  In this case a corruption index higher than 1 will be returned,
otherwise the values of the corruption index will always be bounded by 0 and 1.

See "New Mode Tracking Methods in Aeroelastic Analysis" by Eldred, Vankayya, and
Anderson.
"""
function correlate_eigenmodes(C)

    # get row permutation that puts maximum values on the diagonals
    nev = size(C, 1)
    tmp = Vector{real(eltype(C))}(undef, nev)
    perm = zeros(Int, nev)
    corruption = Vector{real(eltype(C))}(undef, nev)
    for i = 1:nev

        # rank each value in this row
        for j = 1:nev
            tmp[j] = abs(C[i,j])
        end
        ranked_modes = sortperm(tmp, rev=true)

        # choose the best fit that is not yet assigned
        i1 = ranked_modes[findfirst((x) -> !(x in perm), ranked_modes)]
        i2 = i1 == ranked_modes[1] ? ranked_modes[2] : ranked_modes[1]

        # assign best eigenmode fit, create corruption index
        perm[i] = i1
        corruption[i] = tmp[i2]/tmp[i1]
    end

    return perm, corruption
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct InboundNatRuleListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{InboundNatRule} } # spec name: value
    nextLink::Any # spec type: Union{ Nothing, String } # spec name: nextLink

    function InboundNatRuleListResult(;value=nothing, nextLink=nothing)
        o = new()
        validate_property(InboundNatRuleListResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(InboundNatRuleListResult, Symbol("nextLink"), nextLink)
        setfield!(o, Symbol("nextLink"), nextLink)
        o
    end
end # type InboundNatRuleListResult

const _property_map_InboundNatRuleListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"), Symbol("nextLink")=>Symbol("nextLink"))
const _property_types_InboundNatRuleListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{InboundNatRule}", Symbol("nextLink")=>"String")
Base.propertynames(::Type{ InboundNatRuleListResult }) = collect(keys(_property_map_InboundNatRuleListResult))
Swagger.property_type(::Type{ InboundNatRuleListResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_InboundNatRuleListResult[name]))}
Swagger.field_name(::Type{ InboundNatRuleListResult }, property_name::Symbol) =  _property_map_InboundNatRuleListResult[property_name]

function check_required(o::InboundNatRuleListResult)
    true
end

function validate_property(::Type{ InboundNatRuleListResult }, name::Symbol, val)
end


text:
mutable struct TTestResult<:TestResult
    type::Symbol  # either :one, :paired, :simple, or :welch
    levels::Union{Vector{String}, Nothing}
    n_obs::IntVector
    means::FloatVector
    sds::FloatVector
    ci::FloatMatrix
    statistic::Float64
    dof::Union{Int, Float64}
    pval::Float64
    cohen_d::Union{Float64, Nothing}
end

function show(io::IO, ttr::TTestResult)
    if ttr.type == :one
        println(io, "\n        One Sample t-test\n")
        @printf(io, "Data (Mean ± SD [95%% CI]): %.2f ± %.2f [%.2f, %.2f]\n",
                ttr.means[1], ttr.sds[1], ttr.ci[1,1], ttr.ci[1,2])
        @printf(io, "  difference from %.2f\n",
                ttr.means[2])
    elseif ttr.type == :paired
        println(io, "\n        Paired t-test\n")
        @printf(io, "Difference (Mean ± SD [95%% CI]): %.2f ± %.2f [%.2f, %.2f]\n",
                ttr.means[1], ttr.sds[1], ttr.ci[1,1], ttr.ci[1,2])
    elseif ttr.type == :simple
        println(io, "\n        Simple Two Sample t-test\n")
        println(io, "Data (Mean ± SD [95%% CI])")
        @printf(io, "  %s: %.2f ± %.2f [%.2f, %.2f]\n",
                ttr.levels[1], ttr.means[1], ttr.sds[1],
                ttr.ci[1,1], ttr.ci[1,2])
        @printf(io, "  %s: %.2f ± %.2f [%.2f, %.2f]\n",
                ttr.levels[2], ttr.means[2], ttr.sds[2],
                ttr.ci[2,1], ttr.ci[2,2])
    elseif ttr.type == :welch
        println(io, "\n        Welch Two Sample t-test\n")
        println(io, "Data (Mean ± SD [95%% CI])")
        @printf(io, "  %s: %.2f ± %.2f [%.2f, %.2f]\n",
                ttr.levels[1], ttr.means[1], ttr.sds[1],
                ttr.ci[1,1], ttr.ci[1,2])
        @printf(io, "  %s: %.2f ± %.2f [%.2f, %.2f]\n",
                ttr.levels[2], ttr.means[2], ttr.sds[2],
                ttr.ci[2,1], ttr.ci[2,2])
    end
    println(io, "\nResult:\n", coeftable(ttr))
end

function coeftable(ttr::TTestResult)
    if (ttr.type == :one) || (ttr.type == :paired)
        CoefTable(
            [ttr.statistic, ttr.dof, ttr.pval],
            ["t", "DF", "P-value"],
            [""],
            3, 1
        )
    elseif (ttr.type == :simple) || (ttr.type == :welch)
        CoefTable(
            [ttr.statistic, ttr.dof, ttr.pval, ttr.cohen_d],
            ["t", "DF", "P-value", "Cohen's d"],
            [""],
            3, 1
        )
    else
        error("invalid t-test type $ttr.type")
    end
end

"""
    t_test(x, μ0::AbstractFloat=0.0)

One sample t-test

The argument `x` can be numeric `Vector`.

One sample t-test will be performed to compare between the mean of `x1` and `μ0`.
"""
function t_test(x::NumVector; μ0::AbstractFloat=0.0)::TTestResult
    # check types
    if !(x isa FloatVector)
        @warn "x is converted to Float64"
        x = convert(Vector{Float64}, x)
    end

    # check length
    n = length(x)
    n > 1 || error("size of x is not enough: actual size is $n")

    μ = mean(x)
    σ = std(x)
    t_stat = (μ - μ0) / (σ / sqrt(n))
    ν = n - 1
    distr_t = TDist(ν)
    pval = 2.0 * D.ccdf(distr_t, abs(t_stat))
    ci = μ .+ quantile.(distr_t, [0.025, 0.975]) .* (σ/sqrt(n))

    res = TTestResult(
        :one, nothing, [n], [μ, μ0], [σ], reshape(ci, (1, 2)),
        t_stat, ν, pval, nothing
    )
    return res
end

"""
    t_test(x1, x2, paired::Bool=false, varequal::Bool=false)

Two sample t-test

The arguments `x1` and `x2` can be numeric `Vector`.

Paired t-test will be performed if the keyword argument `paired` is `true`. Simple two sample t-test will be performed if the keyword argument `varequal` is `true`. Otherwise, Welch two sample t-test will be performed.
"""
function t_test(
    x1::NumVector,
    x2::NumVector;
    levels::Union{Vector{String}, Nothing}=nothing,
    paired::Bool=false,
    varequal::Bool=false
)::TTestResult
    # use temporal sample names
    if !paired & isnothing(levels)
        @info "Using the first and second sample names as 'Sample1' and 'Sample2'"
        levels = ["Sample1", "Sample2"]
    end

    # check types
    if !(x1 isa FloatVector)
        @warn "$levels[1] is converted to Float64"
        x1 = convert(Vector{Float64}, x1)
    end
    if !(x2 isa FloatVector)
        @warn "$levels[2] is converted to Float64"
        x2 = convert(Vector{Float64}, x2)
    end

    # check length
    n1 = length(x1)
    n2 = length(x2)
    n1 > 1 || error("size of $levels[1] is not enough: actual size is $n1")
    n2 > 1 || error("size of $levels[2] is not enough: actual size is $n2")

    # paired, simple, or Welch
    if paired
        n1 == n2 || error("size of $levels[1] and $levels[2] must be the same: ($n1, $n2)")
        t_test_paired(x1, x2)
    elseif varequal
        t_test_var_equal(x1, x2, levels)
    else
        t_test_welch(x1, x2, levels)
    end
end

function t_test_paired(
    x1::FloatVector,
    x2::FloatVector,
)::TTestResult
    n = length(x1)
    Δ = x1 .- x2
    μ = mean(Δ)
    σ = std(Δ)
    t_stat = μ / (σ / sqrt(n))
    ν = n - 1
    distr_t = TDist(ν)
    pval = 2.0 * D.ccdf(distr_t, abs(t_stat))
    ci = μ .+ quantile.(distr_t, [0.025, 0.975]) .* (σ/sqrt(n))

    res = TTestResult(
        :paired, nothing, [n], [μ], [σ], reshape(ci, (1,2)),
        t_stat, ν, pval, nothing
    )
    return res
end

function t_test_var_equal(
    x1::FloatVector,
    x2::FloatVector,
    levels::Vector{String}
)::TTestResult
    m1 = mean(x1)
    m2 = mean(x2)
    ms = [m1, m2]
    v1 = var(x1)
    v2 = var(x2)
    sds = sqrt.([v1, v2])
    n1 = length(x1)
    n2 = length(x2)
    ns = [n1, n2]

    Δ = m1 - m2
    ν = n1 + n2 - 2
    σ = sqrt(((n1-1)*v1 + (n2-1)*v2) / ν)
    t_stat = Δ / (σ * sqrt(1/n1 + 1/n2))
    pval = 2.0 * D.ccdf(TDist(ν), abs(t_stat))
    distr_t = TDist.(ns .- 1)
    qm = [quantile(t, p) for (t, p) in Iterators.product(distr_t, [0.025, 0.975])]
    ci = ms .+ qm .* sds ./ sqrt.(ns)
    cohen_d = abs(Δ) / sqrt((v1 + v2) / 2.0)

    res = TTestResult(
        :simple, levels, ns, ms, sds, ci,
        t_stat, ν, pval, cohen_d
    )
    return res
end

function t_test_welch(
    x1::FloatVector,
    x2::FloatVector,
    levels::Vector{String}
)::TTestResult
    m1 = mean(x1)
    m2 = mean(x2)
    ms = [m1, m2]
    v1 = var(x1)
    v2 = var(x2)
    sds = sqrt.([v1, v2])
    n1 = length(x1)
    n2 = length(x2)
    ns = [n1, n2]
    vn1 = v1 / n1
    vn2 = v2 / n2

    Δ = m1 - m2
    σ = sqrt(vn1 + vn2)
    t_stat = Δ / σ
    ν = (vn1 + vn2)^2 / ((vn1^2/(n1-1)) + (vn2^2/(n2-1)))
    pval = 2.0 * D.ccdf(TDist(ν), abs(t_stat))
    distr_t = TDist.(ns .- 1)
    qm = [quantile(t, p) for (t, p) in Iterators.product(distr_t, [0.025, 0.975])]
    ci = ms .+ qm .* sds ./ sqrt.(ns)
    cohen_d = abs(Δ) / sqrt((v1 + v2) / 2.0)

    res = TTestResult(:welch, levels, ns, ms, sds, ci,
                      t_stat, ν, pval, cohen_d)
    return res
end

"""
    TTestModel

A `TestModel` type for t-test

# Members

- `X`: model matrix
- `y`: vector to compare its means
- `type`: either :one, :paired, :simple, or :welch
- `levels`: `String` vector of the unique levels (only 1 or 2 elements are allowed)
"""
struct TTestModel{T<:Real} <: TestModel
    x1::AbstractVector{T}
    x2::AbstractVector{T}
    type::Symbol
    levels::Vector{String}
    function TTestModel(
        X::Matrix{T},
        y::AbstractVector{T},
        type::Symbol,
        levels::Vector{String};
        μ0::T = 0.0
    ) where T
        if type == :one
            new{T}(y, [μ0], :one, levels)
        else
            dummy = X[:,2]
            x1 = y[dummy .== 0]
            x2 = y[dummy .== 1]
            new{T}(x1, x2, type, levels)
        end
    end
end

function fit(obj::TTestModel)
    if obj.type == :one
        t_test(obj.x1, μ0=obj.x2[1])
    elseif obj.type == :paired
        t_test(obj.x1, obj.x2, paired=true, levels=obj.levels)
    else
        t_test(obj.x1, obj.x2, paired=false, varequal=(obj.type == :simple), levels=obj.levels)
    end
end

function fit(
        ::Type{TTestModel},
        X::AbstractMatrix{<:Real},
        y::AbstractVector{<:Real},
        type::Symbol,
        levels::Vector{String};
        μ0::Real = 0.0
    )
    ttest_obj = TTestModel(X, y, type, levels, μ0=μ0)
    fit(ttest_obj)
end
function fit(
        ::Type{TTestModel},
        X::AbstractMatrix{<:Real},
        y::AbstractVector{<:Real},
        type::Symbol
    )
    @info "Using the first and second sample names as 'Sample1' and 'Sample2'"
    fit(TTestModel, X, y, type, ["Sample1", "Sample2"], μ0=0.0)
end

"""
    t_test(X, y, paired::Bool=false, varequal::Bool=false)

An alias for `fit(TTestModel, X, y, levels)`

The arguments `X` and `y` can be a `Matrix` and a `Vector` or a `Formula` and `DataFrame`.

The argument `type` must be either `:one`, `:paired`, `:simple`, or `welch`.

The argument `levels` can be a `Vector{String}`.

The keyword argument `μ0` can be a `Real` to use one sample t-test (default = 0.0). Ignored if the `type` is not `:one`.
"""
function t_test(
    X::AbstractMatrix{<:Real},
    y::AbstractVector{<:Real},
    paired::Bool=false,
    varequal::Bool=false
)
    if paired
        fit(TTestModel, X, y, :paired)
    elseif varequal
        fit(TTestModel, X, y, :simple)
    else
        fit(TTestModel, X, y, :welch)
    end
end
function t_test(
    formula::FormulaTerm{Term,Term},
    df::DataFrame,
    paired::Bool=false,
    varequal::Bool=false
)
    # obtain levels
    f = apply_schema(formula, schema(formula, df))
    levels = DataAPI.levels(df[!,f.rhs.terms[1].sym])
    if eltype(levels) <: Real
        levels = string.(levels)
    end

    # model matrix
    model_frame = ModelFrame(formula, df)
    X = ModelMatrix(model_frame).m
    y = response(model_frame)

    # perform t-test
    if paired
        fit(TTestModel, X, y, :paired, levels)
    elseif varequal
        fit(TTestModel, X, y, :simple, levels)
    else
        fit(TTestModel, X, y, :welch, levels)
    end
end


text:
"""
Power nonlinearity, with energy C ∫ρ^α where ρ is the density
"""
struct PowerNonlinearity
    C::Real
    α::Real
end
(P::PowerNonlinearity)(basis) = TermPowerNonlinearity(basis, P.C, P.α)

struct TermPowerNonlinearity <: TermNonlinear
    basis::PlaneWaveBasis
    C::Real
    α::Real
end

function ene_ops(term::TermPowerNonlinearity, ψ, occ; ρ, kwargs...)
    E = term.C * sum(ρ .^ term.α) * term.basis.dvol
    potential = @. term.C * term.α * ρ^(term.α-1)

    # In the case of collinear spin, the potential is spin-dependent
    ops = [RealSpaceMultiplication(term.basis, kpoint, potential[:, :, :, kpoint.spin])
           for kpoint in term.basis.kpoints]
    (E=E, ops=ops)
end

_pnl_kernel(C, α, ρ) = @. C * α * (α-1) * ρ^(α-2)

function compute_kernel(term::TermPowerNonlinearity; ρ, kwargs...)
    K = Diagonal(vec(_pnl_kernel(term.C, term.α, ρ)))
end

function apply_kernel(term::TermPowerNonlinearity, δρ; ρ, kwargs...)
    _pnl_kernel(term.C, term.α, ρ) .* δρ
end


text:
function bob(stimulus::AbstractString)

end


text:
module Tables

import DataFrames
using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, table, template, register_normal_element

export Column, DataTablePagination, DataTableOptions, DataTable

register_normal_element("q__table", context = @__MODULE__)

const ID = "__id"

#===#

Base.@kwdef mutable struct Column
  name::String
  required::Bool = false
  label::String = name
  align::Symbol = :left
  field::String = name
  sortable::Bool = true
end

function Column(name::String; args...)
  Column(name = name; args...)
end

function Base.Symbol(v::Vector{Column}) :: Vector{Symbol}
  [Symbol(c.name) for c in v]
end

Base.@kwdef mutable struct DataTablePagination
  sort_by::Symbol = :desc
  descending::Bool = false
  page::Int = 1
  rows_per_page::Int = 10
end

Base.@kwdef mutable struct DataTableOptions
  addid::Bool = false
  idcolumn::String = "ID"
  columns::Union{Vector{Column},Nothing} = nothing
end

Base.@kwdef mutable struct DataTable{T<:DataFrames.DataFrame}
  data::T = DataFrames.DataFrame()
  opts::DataTableOptions = DataTableOptions()
end

function DataTable(data::T) where {T<:DataFrames.DataFrame}
  DataTable(data, DataTableOptions())
end

#===#

function active_columns(t::T)::Vector{Column} where {T<:DataTable}
  t.opts.columns !== nothing ? t.opts.columns : [Column(string(name)) for name in DataFrames.names(t.data)]
end

function columns(t::T)::Vector{Column} where {T<:DataTable}
  columns = active_columns(t) |> copy

  if t.opts.addid
    pushfirst!(columns, Column(t.opts.idcolumn, true, t.opts.idcolumn, :right, t.opts.idcolumn, true))
  end

  columns
end

function rows(t::T)::Vector{Dict{String,Any}} where {T<:DataTable}
  rows = []

  for (count, row) in enumerate(DataFrames.eachrow(t.data))
    r = Dict()

    if t.opts.addid
      r[t.opts.idcolumn] = count
    end

    r[ID] = count
    for column in active_columns(t)
      r[column.name] = row[Symbol(column.name)]
    end

    push!(rows, r)
  end

  rows
end

function data(t::T, fieldname::Symbol; datakey = "data_$fieldname", columnskey = "columns_$fieldname")::Dict{String,Any} where {T<:DataTable}
  Dict(
    columnskey  => columns(t),
    datakey     => rows(t)
  )
end

function table( fieldname::Symbol,
                                    args...;
                                    rowkey::String = ID,
                                    title::String = "",
                                    datakey::String = "data_$fieldname",
                                    columnskey::String = "columns_$fieldname",
                                    kwargs...) :: String

  q__table(args...; kw(
    [Symbol(":data") => "$fieldname.$datakey", Symbol(":columns") => "$fieldname.$columnskey", Symbol("row-key") => rowkey,
    :fieldname => fieldname, kwargs...])...)
end

#===#

function Stipple.render(t::T, fieldname::Union{Symbol,Nothing} = nothing) where {T<:DataTable}
  data(t, fieldname)
end

function Stipple.render(dtp::DataTablePagination, fieldname::Union{Symbol,Nothing} = nothing)
  Dict(:sortBy => dtp.sort_by, :descending => dtp.descending, :page => dtp.page, :rowsPerPage => dtp.rows_per_page)
end

#===#

function Stipple.watch(vue_app_name::String, fieldtype::R{T}, fieldname::Symbol, channel::String, model::M)::String where {M<:ReactiveModel,T<:DataTable}
  string(vue_app_name, raw".\$watch('", fieldname, "', function(newVal, oldVal){

  });\n\n")
end

#===#

function Base.parse(::Type{DataTablePagination}, d::Dict{String,Any})
  dtp = DataTablePagination()

  dtp.sort_by = get!(d, "sortBy", "desc") |> Symbol
  dtp.page = get!(d, "page", 1)
  dtp.descending = get!(d, "descending", false)
  dtp.rows_per_page = get!(d, "rowsPerPage", 10)

  dtp
end

function Base.parse(::Type{DataTable}, ::Dict{String,Any})
  # todo: add support
end

function Base.parse(::Type{DataTable{DataFrames.DataFrame}}, ::Dict{String,Any})
  # error("Not implemented") # todo implement parser
end

end


text:
using MCMCChains
using StatsBase
using Test
using Random

# Tests for missing values.
function testdiff(cdf1, cdf2)
    m1 = convert(Array, cdf1)
    m2 = convert(Array, cdf2)
    return all(((x, y),) -> isapprox(x, y; atol=1e-2), zip(m1, m2))
end

@testset "utils" begin
    x = [1, missing, 3, 2]
    @test MCMCChains.cummean(x) == [1., 1., 2., 2.]
end

@testset "diagnostic functions" begin
    Random.seed!(1234)
    nchains = 2
    chn = Chains(randn(5000, 2, nchains))

    # Add missing values.
    chn_m = Chains(cat(chn.value, ones(1, 2, nchains) .* missing, dims = 1))

    # Currently we only check if diagnostic functions throw a `MethodError` if the
    # element type of the chain is not a subtype of `Real`.
    @test_throws MethodError gelmandiag(chn_m)

    gw_1 = gewekediag(chn)
    gw_2 = gewekediag(chn_m)
    hd_1 = heideldiag(chn)
    hd_2 = heideldiag(chn_m)
    rf_1 = rafterydiag(chn)
    rf_2 = rafterydiag(chn_m)

    @testset "diagnostics missing tests" for i in 1:nchains
        @test testdiff(gw_1, gw_2)
        @test testdiff(hd_1, hd_2)
        @test testdiff(rf_1, rf_2)
    end

    @test_throws MethodError discretediag(chn_m)
end


text:
"""
    PriorSampling()

Sampler which will use the prior distribution from the given model to provide samples.
"""
struct PriorSampling <: AbstractBQSampler end

function AbstractMCMC.step(rng::Random.AbstractRNG, model, ::PriorSampling; kwargs...)
    return rand(rng, p_0(model)), nothing
end

function AbstractMCMC.step(
    rng::Random.AbstractRNG, model, ::PriorSampling, ::Any; kwargs...
)
    return rand(rng, p_0(model)), nothing
end


text:
# ---
# title: 1471. The k Strongest Values in an Array
# id: problem1471
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Sort
# link: <https://leetcode.com/problems/the-k-strongest-values-in-an-array/description/>
# hidden: true
# ---
#
# Given an array of integers `arr` and an integer `k`.
#
# A value `arr[i]` is said to be stronger than a value `arr[j]` if `|arr[i] - m|
# > |arr[j] - m|` where `m` is the **median** of the array.
# If `|arr[i] - m| == |arr[j] - m|`, then `arr[i]` is said to be stronger than
# `arr[j]` if `arr[i] > arr[j]`.
#
# Return _a list of the strongest`k`_ values in the array. return the answer
# **in any arbitrary order**.
#
# **Median** is the middle value in an ordered integer list. More formally, if
# the length of the list is n, the median is the element in position `((n - 1) /
# 2)` in the sorted list  **(0-indexed)**.
#
#   * For `arr = [6, -3, 7, 2, 11]`, `n = 5` and the median is obtained by sorting the array `arr = [-3, 2, 6, 7, 11]` and the median is `arr[m]` where `m = ((5 - 1) / 2) = 2`. The median is `6`.
#   * For `arr = [-7, 22, 17, 3]`, `n = 4` and the median is obtained by sorting the array `arr = [-7, 3, 17, 22]` and the median is `arr[m]` where `m = ((4 - 1) / 2) = 1`. The median is `3`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [1,2,3,4,5], k = 2
#     Output: [5,1]
#     Explanation: Median is 3, the elements of the array sorted by the strongest are [5,1,4,2,3]. The strongest 2 elements are [5, 1]. [1, 5] is also **accepted** answer.
#     Please note that although |5 - 3| == |1 - 3| but 5 is stronger than 1 because 5 > 1.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [1,1,3,5,5], k = 2
#     Output: [5,5]
#     Explanation: Median is 3, the elements of the array sorted by the strongest are [5,5,1,1,3]. The strongest 2 elements are [5, 5].
#
#
# **Example 3:**
#
#
#
#     Input: arr = [6,7,11,7,6,8], k = 5
#     Output: [11,8,6,6,7]
#     Explanation: Median is 7, the elements of the array sorted by the strongest are [11,8,6,6,7,7].
#     Any permutation of [11,8,6,6,7] is **accepted**.
#
#
# **Example 4:**
#
#
#
#     Input: arr = [6,-3,7,2,11], k = 3
#     Output: [-3,11,2]
#
#
# **Example 5:**
#
#
#
#     Input: arr = [-7,22,17,3], k = 2
#     Output: [22,17]
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 10^5`
#   * `-10^5 <= arr[i] <= 10^5`
#   * `1 <= k <= arr.length`
#
#
## @lc code=start
using LeetCode

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


text:
function with_args(f, args::Vector{String} = String[])
    old = copy(ARGS)
    empty!(ARGS)
    append!(ARGS, args)
    ret = f()
    empty!(ARGS)
    append!(ARGS, old)
    return ret
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using Random
using Distributions, StatsBase

@testset "BAT.HistogramAsUvDistribution" begin
    Random.seed!(123)
    μ, σ = 1.23, 0.74
    true_dist = Normal(μ, σ)
    h = Histogram(μ-10σ:σ/10:μ+10σ)
    append!(h, rand(true_dist, 10^7))
    d = BAT.HistogramAsUvDistribution(h)
    @test isapprox(μ, d.μ, atol = 0.01)
    @test isapprox(σ, d.σ, atol = 0.01)
    fit_dist = fit(Normal, rand(d, 10^7))
    @test isapprox(μ, fit_dist.μ, atol = 0.01)
    @test isapprox(σ, fit_dist.σ, atol = 0.01)
end


text:
export project_nuclear!

function project_nuclear!(
                          x     ::Array{TF,2},
                          sigma ::TF,
                          mode  #not an option for matrix inputs
                          ) where {TF<:Real}
"""
Project the matrix onto the set of matrices with nuclear norm less then or equal to sigma
outputs a vector. Uses singular value decomposition
"""

  #slow reference, like matlab
  #(U,S,V) = svd(x)
  #x = U[:,1:r]*diagm(S)*(V[:,1:r])'

  #F=svd!(x)
  F = svd(x) #obtain svd
  project_l1_Duchi!(F.S, sigma) # project singular values onto the l1-ball of size sigma

  #copy!(x,F.U * diagm(F.S) * F.Vt) #reconstruct
  x .= F.U * diagm(0 => F.S) * F.Vt #reconstruct
  x = vec(x)
  #return vec(x)
end

function project_nuclear!(
                          x     ::Array{TF,3},
                          sigma ::TF,
                          mode  ::Tuple{String,String}
                          ) where {TF<:Real}
"""
Project each slice of the tensor (x,y or z) onto the set of matrices with nuclear norm less then or equal to sigma
outputs a vector. Uses singular value decomposition
"""
if mode[1] == "slice"
  if mode[2] == "x" #project y-z slices, so one slice per gridpoint in the x direction
    Threads.@threads for i=1:size(x,1)
      F=svd(view(x,i,:,:))
      project_l1_Duchi!(F.S, sigma) # project singular values onto the l1-ball of size sigma
      @inbounds x[i,:,:] .= F.U * diagm(0 => F.S) * F.Vt #reconstruct
    end
  elseif mode[2] == "y"
    Threads.@threads for i=1:size(x,2)
      F=svd(view(x,:,i,:))
      project_l1_Duchi!(F.S, sigma) # project singular values onto the l1-ball of size sigma
      @inbounds x[:,i,:].= F.U * diagm(0 => F.S) * F.Vt
    end
  elseif mode[2] == "z"
    Threads.@threads for i=1:size(x,3)
      F=svd(view(x,:,:,i))
      project_l1_Duchi!(F.S, sigma) # project singular values onto the l1-ball of size sigma
      @inbounds x[:,:,i].= F.U * diagm(0 => F.S) * F.Vt
    end
  end #end if mode[2]
else
 error("mode[1] for nuclear norm projections can only be: slice")
end #if mode[1]

x=vec(x)
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/NodeNumbering.jl/blob/master/LICENSE

"""
    create_RCM_adjacency(adjacency::Dict{Int, Vector{Int}}, finalorder::Dict{Int, Int})

Create an adjacency Dict of the RCM ordered nodes. A part of the Reverse Cuthill-McKee Algorithm.

# Examples

```jldoctest
julia> adjacency = Dict{Int, Vector{Int}}(
                   1 => [3, 8, 9],
                   2 => [3, 8, 7, 6],
                   3 => [1, 2],
                   4 => [8, 9],
                   5 => [7, 8],
                   6 => [2, 7],
                   7 => [5, 2, 6],
                   8 => [1, 2, 4, 5],
                   9 => [1, 4]);

julia> finalorder =  Dict{Int, Int}(
                     7 => 2,
                     9 => 9,
                     4 => 8,
                     2 => 3,
                     3 => 5,
                     5 => 4,
                     8 => 6,
                     6 => 1,
                     1 => 7);

julia> create_RCM_adjacency(adjacency, finalorder)
Dict{Int64,Array{Int64,1}} with 9 entries:
  7 => [5, 6, 9]
  9 => [7, 8]
  4 => [2, 6]
  2 => [4, 3, 1]
  3 => [5, 6, 2, 1]
  8 => [6, 9]
  5 => [7, 3]
  6 => [7, 3, 8, 4]
  1 => [3, 2]

```

# References

* Wikipedia contributors. "Adjacency list". Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia, 7 Jun. 2017. Web. 17 Jul. 2017. https://en.wikipedia.org/wiki/Adjacency_list
"""
function create_RCM_adjacency(adjacency::Dict{Int,Vector{Int}}, finalorder::Dict{Int,Int})
    RCM_adjacency = Dict(finalorder[i] => [finalorder[j] for j in c] for (i,c) in adjacency)
    return RCM_adjacency
end


text:
const primes = [53, 97, 193, 389, 769, 1543, 3079, 6151, 12289, 24593, 49157, 98317, 196613, 393241, 786433, 1572869, 3145739, 6291469, 12582917, 25165843, 50331653, 100663319, 201326611, 402653189, 805306457, 1610612741]

# Layout (for capacity C, value size V)
# -   C bytes for used hash bucket marker
# -   C bytes for used collition bucket marker
# - 16C bytes for 2C key hash values
# -  8C bytes for 2C next entry values
# -  8C bytes for 2C value sizes
# - 2VC bytes for 2C values, each of max V bytes
#
# Note:
# - actual dict capacity is greater than the requested capacity (actual capacity = 2 * next higher value from primes)
# - with this scheme, loadfactor is around 0.5
struct ShmDict
    capacity::UInt32            # requested capacity (actual capacity > 2 * requested capacity)
    maxvalsize::UInt32          # max bytes in value (or value + key if keeping keys)
    maxkeysize::UInt8           # max bytes in key
    keepkeys::Bool              # whether to keep keys along with values (used for exact comparison, take up space from value buffer)
    shmid::Cint                 # shm id
    shmptr::Ptr{Void}           # shm pointer
    lck::NamedSemaphore         # lock for dict operations
    used::Vector{Bool}          # bucket in use marker
    khash::Vector{UInt64}       # hash of key that occupies the bucket
    next::Vector{UInt32}        # the next collision entry for this bucket (0 if none)
    valsize::Vector{UInt32}     # size occupied in val
    val::Vector{UInt8}          # value bytes

    function ShmDict(path::Union{String,Tuple{String,Integer}}, capacity::Integer, maxvalsize::Integer; keepkeys::Bool=true, maxkeysize::Integer=32, create::Bool=false, create_exclusive::Bool=false, permission::Integer=0o660)
        if keepkeys
            maxvalsize = maxvalsize + maxkeysize + 1
        end
        capacity = primes[min(length(primes), searchsortedfirst(primes, capacity))]
        memsz, offsets = _shmsize(capacity, maxvalsize)
        id = 0
        if isa(path, Tuple)
            path,id = path
        end
        lck = NamedSemaphore("/$(hash(path))")
        tok = ftok(path, id)
        shmid = shmget(tok, memsz; create=create, create_exclusive=create_exclusive, permission=permission)
        shmptr = shmat(shmid)

        used = unsafe_wrap(Array, convert(Ptr{Bool}, shmptr), (2*capacity,))
        khash = unsafe_wrap(Array, convert(Ptr{UInt64}, shmptr+shift!(offsets)), (2*capacity,))
        next = unsafe_wrap(Array, convert(Ptr{UInt32}, shmptr+shift!(offsets)), (2*capacity,))
        valsize = unsafe_wrap(Array, convert(Ptr{UInt32}, shmptr+shift!(offsets)), (2*capacity,))
        val = unsafe_wrap(Array, convert(Ptr{UInt8}, shmptr+shift!(offsets)), (2*maxvalsize*capacity,))

        new(capacity, maxvalsize, maxkeysize, keepkeys, shmid, shmptr, lck, used, khash, next, valsize, val)
    end
end

function _shmsize(C, V)
    offsets = Vector{Int}()
    sz = 2*C
    sz = ceil(Int, sz/8) * 8
    push!(offsets, sz)
    sz += 16*C
    sz = ceil(Int, sz/8) * 8
    push!(offsets, sz)
    sz += 8*C
    sz = ceil(Int, sz/8) * 8
    push!(offsets, sz)
    sz += 8*C
    sz = ceil(Int, sz/8) * 8
    push!(offsets, sz)
    sz += 2*V*C
    sz, offsets
end

_byte_repr(val::Vector{UInt8}) = val
_byte_repr(val::String) = convert(Vector{UInt8}, val)
_byte_repr(val::Number) = reinterpret(UInt8, [val])
function _byte_repr(val)
    iob = IOBuffer()
    serialize(iob, val)
    take!(iob)
end

_from_byte_repr(::Type{Vector{UInt8}}, val::Vector{UInt8}) = val
_from_byte_repr(::Type{T}, val::Vector{UInt8}) where {T<:Number} = reinterpret(T, val)[1]
_from_byte_repr(::Type{String}, val::Vector{UInt8}) = String(val)
_from_byte_repr(::Type{T}, val::Vector{UInt8}) where {T} = deserialize(IOBuffer(val))

function _key_matches(D::ShmDict, key, bucket)
    if D.keepkeys
        keystart = Int((D.maxvalsize * (bucket-1)) + 1)
        keysz = D.val[keystart]
        keystart += 1
        keyend = Int(keystart + keysz - 1)
        key == SubArray(D.val, (keystart:keyend,))
    else
        true
    end
end

function _val_bytes(D::ShmDict, bucket)
    if D.keepkeys
        valstart = Int((D.maxvalsize * (bucket-1)) + 1)
        keysz = D.val[valstart]
        keyoffset = 1 + keysz
        valstart += keyoffset
        valend = Int(valstart + D.valsize[bucket] - keyoffset - 1)
    else
        valstart = Int((D.maxvalsize * (bucket-1)) + 1)
        valend = Int(valstart + D.valsize[bucket] - 1)
    end
    SubArray(D.val, (valstart:valend,))
end

haskey(D::ShmDict, key) = haskey(D, _byte_repr(key))
function haskey(D::ShmDict, key::Vector{UInt8})
    khash = hash(key)
    bucket = (khash % D.capacity) + 1
    if D.used[bucket]
        while bucket > 0
            (D.khash[bucket] === khash) && _key_matches(D, key, bucket) && (return true)
            bucket = D.next[bucket]
        end
    end
    return false
end

getindex(D::ShmDict, key) = getindex(D, _byte_repr(key))
function getindex(D::ShmDict, key::Vector{UInt8})
    khash = hash(key)
    bucket = (khash % D.capacity) + 1
    if D.used[bucket]
        while bucket > 0
            (D.khash[bucket] === khash) && _key_matches(D, key, bucket) && (return _val_bytes(D, bucket))
            bucket = D.next[bucket]
        end
    end
    error("key not found")
end

function _setbucket(D::ShmDict, bucket, khash::UInt64, key::Vector{UInt8}, val::Vector{UInt8}, next::UInt32)
    D.used[bucket] = true
    D.khash[bucket] = khash
    D.next[bucket] = next
    D.valsize[bucket] = length(val)
    valstart = (D.maxvalsize * (bucket-1)) + 1
    valend = valstart + D.valsize[bucket] - 1
    if D.keepkeys
        keysz = UInt8(length(key))
        D.valsize[bucket] += (1 + keysz)
        D.val[valstart] = keysz
        copy!(D.val, valstart+1, key, 1, keysz)
        valstart += (1 + keysz)
    end
    copy!(D.val, valstart, val, 1, length(val))
    nothing
end

function _setbucket(D::ShmDict, bucket, khash::UInt64, keyval::Vector{UInt8}, next::UInt32)
    D.used[bucket] = true
    D.khash[bucket] = khash
    D.next[bucket] = next
    D.valsize[bucket] = length(keyval)
    valstart = (D.maxvalsize * (bucket-1)) + 1
    valend = valstart + D.valsize[bucket] - 1
    copy!(D.keyval, valstart, keyval, 1, length(keyval))
    nothing
end

setindex!(D::ShmDict, val, key) = setindex!(D, _byte_repr(val), _byte_repr(key))
function setindex!(D::ShmDict, val::Vector{UInt8}, key::Vector{UInt8})
    lkey = length(key)
    lval = length(val)
    if D.keepkeys
        (lkey > D.maxkeysize) && error("key size $lkey greater than max allowed $(D.maxkeysize)")
        lval += (1 + lkey)
    end
    (lval > D.maxvalsize) && error("value size $(lval-lkey-1) greater than max allowed $(D.maxvalsize-D.maxkeysize-1)")
    khash = hash(key)
    bucket = (khash % D.capacity) + 1
    nextbucket = UInt32(0)
    if D.used[bucket]
        prevbucket = bucket
        while bucket > 0
            # navigate to matching entry or last entry
            (D.khash[bucket] === khash) && _key_matches(D, key, bucket) && break
            prevbucket = bucket
            bucket = D.next[bucket]
        end
        if bucket == 0
            # if no matching entry, find a free collision entry and chain it in
            bucket = findfirst(SubArray(D.used, (Int(D.capacity+1):Int(2*D.capacity),)), false)
            (bucket == 0) && error("no space left in dict")
            bucket += D.capacity
            D.next[prevbucket] = bucket
        else
            nextbucket = D.next[bucket]
        end
    end
    _setbucket(D, bucket, khash, key, val, nextbucket)
    val
end

delete!(D::ShmDict, key) = delete!(D, _byte_repr(key))
function delete!(D::ShmDict, key::Vector{UInt8})
    khash = hash(key)
    bucket = (khash % D.capacity) + 1
    if D.used[bucket]
        prevbucket = bucket
        while bucket > 0
            # navigate to matching entry
            (D.khash[bucket] === khash) && _key_matches(D, key, bucket) && break
            prevbucket = bucket
            bucket = D.next[bucket]
        end
        if bucket > 0
            nextbucket = D.next[bucket]
            if (prevbucket == bucket) && (nextbucket > 0)
                valstart = Int((D.maxvalsize * (nextbucket-1)) + 1)
                valend = Int(valstart + D.valsize[nextbucket] - 1)
                _setbucket(D, bucket, D.khash[nextbucket], SubArray(D.val, (valstart:valend,)), D.next[nextbucket])
                D.used[nextbucket] = false
            else
                D.next[prevbucket] = nextbucket
                D.used[bucket] = false
            end
            return D
        end
    end

    error("key not found")
end

function Base.empty!(D::ShmDict)
    memsz, offsets = _shmsize(D.capacity, D.maxvalsize)
    ccall(:memset, Ptr{Void}, (Ptr{Void}, Cint, Csize_t), D.shmptr, 0, memsz)
    D
end

function Base.delete!(D::ShmDict)
    shmrm(D.shmid)
    delete!(D.lck)
    nothing
end

function Base.close(D::ShmDict)
    shmdt(D.shmptr)
    close(D.lck)
    nothing
end


text:
tOCbd = TimerOutput()
for i in 10:10:500
    global numLin = i
    global grid = Vector{Float64}[]
    for p in product(LinRange(gridMin, gridMax, numLin),LinRange(gridMin, gridMax, numLin))
        push!(grid, collect(p))
    end
    evalBellCbd()
end

max_stars_repo_name: uq-aibe/aim

max_issues_repo_name: uq-aibe/aim

max_forks_repo_name: uq-aibe/aim


text:
@testset "align" begin
    @testset "nw_align" begin
        out1, out2 = nw_align("GTCGATCGACTAGCTA","GTCGATGCGCCTGCTA")
        @test out1 == "GTCGAT-CGACTAGCTA"
        @test out2 == "GTCGATGCGCCT-GCTA"
    end

    @testset "triplet_nw_align" begin
        out1, out2 = triplet_nw_align("GTCGATCGACTAGCT","GTCGATAGTCGACTAGCT")
        @test out1 == "GTCGAT---CGACTAGCT"
        @test out2 == "GTCGATAGTCGACTAGCT"
    end

end


text:
using REPL
include("game.jl")
include("TUI_like.jl")
import .Tetris

function get2c()
    t = REPL.TerminalMenus.terminal
    REPL.TerminalMenus.enableRawMode(t) || error("unable to switch to raw mode")
    c = Char(REPL.TerminalMenus.readKey(t.in_stream))
    REPL.TerminalMenus.disableRawMode(t)

    return c
end

# Keyboard pooling input
function kbd_pooling_input()
    inputBuffer = Channel(100)

    @async begin
        while true
            c = get2c()
            put!(inputBuffer, string(c))
        end
    end
    return inputBuffer
end

function create_pipes(game_pipe, player_pipe)
    if !isfifo(game_pipe)
        run(`mkfifo $game_pipe`)
    end
    if !isfifo(player_pipe)
        run(`mkfifo $player_pipe`)
    end
end

function bind_player_pipe_with_channel(pp::String, pc::Channel)
    @async begin
        while true
            # move = read(run(`cat $pp`; wait=true), String)
            move = read(`cat $pp`, String)
            # move = "a"
            move = strip(move, ['\n', '\r'])|>String
            # println("to lendo $(typeof(move))")
            put!(pc, move)
        end
    end
end

function bind_game_buffer_with_output(gc::Channel, gp::String="")
    @async begin
        while true
            game_info = take!(gc)
            if game_info.status
                print_board(game_info)
                gp!="" && send_data_through_pipe(game_info, gp)
            else
                exit()
            end
        end
    end
end

function send_data_through_pipe(game_info::Tetris.GameStatus, gp::String)
    msg = "\"$(game_info.curr_tetromino.name)\n\""
    run(`./teste/sc.sh $msg $gp`; wait=true)
end

function init()
    game_pipe = "/tmp/game_pipe"
    player_pipe = "/tmp/player_pipe"

    human = nothing

    player_buffer = if length(ARGS) <= 0
        human = true
        kbd_pooling_input()
    else
        create_pipes(game_pipe, player_pipe)
        pc = Channel(1000)
        # println("alo")
        bind_player_pipe_with_channel(player_pipe, pc)
        human = false
        player_program = ARGS[1]

        @async run(`./$player_program $game_pipe $player_pipe`)
        pc
    end

    game_buffer = Tetris.start_game(player_buffer)

    if human
        bind_game_buffer_with_output(game_buffer)
    else
        bind_game_buffer_with_output(game_buffer, game_pipe)
    end

    wait()

end

init()

max_stars_repo_name: GabrielTrettel/tetris

max_issues_repo_name: GabrielTrettel/tetris

max_forks_repo_name: GabrielTrettel/tetris