text:
const noalloc_doc = """This method does its computation in place, performing no dynamic memory allocation."""

"""
$(SIGNATURES)

Return the mass of a subtree of a `Mechanism`, rooted at `base` (including
the mass of `base`).
"""
function subtree_mass(base::RigidBody{T}, mechanism::Mechanism{T}) where {T}
    result = isroot(base, mechanism) ? zero(T) : spatial_inertia(base).mass
    for joint in joints_to_children(base, mechanism)
        result += subtree_mass(successor(joint, mechanism), mechanism)
    end
    result
end

"""
$(SIGNATURES)

Return the total mass of the `Mechanism`.
"""
mass(m::Mechanism) = subtree_mass(root_body(m), m)

"""
$(SIGNATURES)

Compute the center of mass of an iterable subset of a `Mechanism`'s bodies in
the given state. Ignores the root body of the mechanism.
"""
function center_of_mass(state::MechanismState, itr)
    update_transforms!(state)
    T = cache_eltype(state)
    mechanism = state.mechanism
    frame = root_frame(mechanism)
    com = Point3D(frame, zero(SVector{3, T}))
    mass = zero(T)
    for body in itr
        if body.inertia !== nothing
            inertia = spatial_inertia(body)
            if inertia.mass > 0
                bodycom = transform_to_root(state, body, false) * center_of_mass(inertia)
                com += inertia.mass * FreeVector3D(bodycom)
                mass += inertia.mass
            end
        end
    end
    com /= mass
    com
end

"""
$(SIGNATURES)

Compute the center of mass of the whole `Mechanism` in the given state.
"""
center_of_mass(state::MechanismState) = center_of_mass(state, bodies(state.mechanism))

const geometric_jacobian_doc = """Compute a geometric Jacobian (also known as a
basic, or spatial Jacobian) associated with a directed path in the `Mechanism`'s
spanning tree, (a collection of `Joint`s and traversal directions) in the given state.

A geometric Jacobian maps the `Mechanism`'s joint velocity vector ``v``
to the twist of the target of the joint path (the body succeeding the last joint
in the path) with respect to the source of the joint path (the body preceding the
first joint in the path).

See also [`path`](@ref), [`GeometricJacobian`](@ref), [`Twist`](@ref).
"""

"""
$(SIGNATURES)

$geometric_jacobian_doc

`transformfun` is a callable that may be used to transform the individual motion
subspaces of each of the joints to the frame in which `out` is expressed.

$noalloc_doc
"""
function geometric_jacobian!(jac::GeometricJacobian, state::MechanismState, path::TreePath, transformfun)
    @boundscheck num_velocities(state) == size(jac, 2) || throw(DimensionMismatch())
    @framecheck jac.body default_frame(target(path))
    @framecheck jac.base default_frame(source(path))
    update_motion_subspaces!(state)
    fill!(jac.angular, 0)
    fill!(jac.linear, 0)
    for i in eachindex(path.edges) # TODO: just use @inbounds here; currently messes with frame check in set_col!
        @inbounds joint = path.edges[i]
        vrange = velocity_range(state, joint)
        @inbounds direction = directions(path)[i]
        for col in eachindex(vrange)
            @inbounds vindex = vrange[col]
            @inbounds Scol = transformfun(state.motion_subspaces.data[vindex])
            direction == PathDirections.up && (Scol = -Scol)
            set_col!(jac, vindex, Scol)
        end
    end
    jac
end

"""
$(SIGNATURES)

$geometric_jacobian_doc

`root_to_desired` is the transform from the `Mechanism`'s root frame to the frame
in which `out` is expressed.

$noalloc_doc
"""
function geometric_jacobian!(out::GeometricJacobian, state::MechanismState, path::TreePath, root_to_desired::Transform3D)
    geometric_jacobian!(out, state, path, S -> transform(S, root_to_desired))
end

"""
$(SIGNATURES)

$geometric_jacobian_doc

See [`geometric_jacobian!(out, state, path, root_to_desired)`](@ref). Uses `state`
to compute the transform from the `Mechanism`'s root frame to the frame in which
`out` is expressed.

$noalloc_doc
"""
function geometric_jacobian!(out::GeometricJacobian, state::MechanismState, path::TreePath)
    if out.frame == root_frame(state.mechanism)
        geometric_jacobian!(out, state, path, identity)
    else
        geometric_jacobian!(out, state, path, inv(transform_to_root(state, out.frame)))
    end
end

"""
$(SIGNATURES)

$geometric_jacobian_doc

The Jacobian is computed in the `Mechanism`'s root frame.

See [`geometric_jacobian!(out, state, path)`](@ref).
"""
function geometric_jacobian(state::MechanismState{X, M, C}, path::TreePath{RigidBody{M}, Joint{M}}) where {X, M, C}
    nv = num_velocities(state)
    angular = Matrix{C}(undef, 3, nv)
    linear = Matrix{C}(undef, 3, nv)
    bodyframe = default_frame(target(path))
    baseframe = default_frame(source(path))
    jac = GeometricJacobian(bodyframe, baseframe, root_frame(state.mechanism), angular, linear)
    geometric_jacobian!(jac, state, path, identity)
end

const point_jacobian_doc = """
Compute the Jacobian mapping the `Mechanism`'s joint velocity vector ``v`` to
the velocity of a point fixed to the target of the joint path (the body
succeeding the last joint in the path) with respect to the source of the joint
path (the body preceding the first joint in the path).
"""

"""
$(SIGNATURES)

$point_jacobian_doc

$noalloc_doc
"""
function _point_jacobian!(Jp::PointJacobian, state::MechanismState, path::TreePath,
                         point::Point3D, transformfun)
    @framecheck Jp.frame point.frame
    update_motion_subspaces!(state)
    fill!(Jp.J, 0)
    p̂ = Spatial.hat(point.v)
    @inbounds for i in eachindex(path.edges)
        joint = path.edges[i]
        vrange = velocity_range(state, joint)
        direction = path.directions[i]
        for col in eachindex(vrange)
            vindex = vrange[col]
            Scol = transformfun(state.motion_subspaces.data[vindex])
            direction == PathDirections.up && (Scol = -Scol)
            newcol =  -p̂ * angular(Scol) + linear(Scol)
            for row in 1:3
                Jp.J[row, vindex] = newcol[row]
            end
        end
    end
    Jp
end

"""
$(SIGNATURES)

$point_jacobian_doc

Uses `state` to compute the transform from the `Mechanism`'s root frame to the
frame in which `out` is expressed if necessary.

$noalloc_doc
"""
function point_jacobian!(out::PointJacobian, state::MechanismState, path::TreePath, point::Point3D)
    if out.frame == root_frame(state.mechanism)
        _point_jacobian!(out, state, path, point, identity)
    else
        root_to_desired = inv(transform_to_root(state, out.frame))
        let root_to_desired = root_to_desired
            _point_jacobian!(out, state, path, point, S -> transform(S, root_to_desired))
        end
    end
end

"""
$(SIGNATURES)

$point_jacobian_doc
"""
function point_jacobian(state::MechanismState{X, M, C},
                        path::TreePath{RigidBody{M}, Joint{M}},
                        point::Point3D) where {X, M, C}
    Jp = PointJacobian(point.frame, Matrix{C}(undef, 3, num_velocities(state)))
    point_jacobian!(Jp, state, path, point)
    Jp
end

const mass_matrix_doc = """Compute the joint-space mass matrix
(also known as the inertia matrix) of the `Mechanism` in the given state, i.e.,
the matrix ``M(q)`` in the unconstrained joint-space equations of motion
```math
M(q) \\dot{v} + c(q, v, w_\\text{ext}) = \\tau
```

This method implements the composite rigid body algorithm.
"""

"""
$(SIGNATURES)

$mass_matrix_doc

$noalloc_doc

The `out` argument must be an ``n_v \\times n_v`` lower triangular
`Symmetric` matrix, where ``n_v`` is the dimension of the `Mechanism`'s joint
velocity vector ``v``.
"""
function mass_matrix!(M::Symmetric, state::MechanismState)
    nv = num_velocities(state)
    @boundscheck size(M, 1) == nv || throw(DimensionMismatch("mass matrix has wrong size"))
    @boundscheck M.uplo == 'L' || throw(ArgumentError(("expected a lower triangular symmetric matrix")))
    update_motion_subspaces!(state)
    update_crb_inertias!(state)
    motion_subspaces = state.motion_subspaces.data
    @inbounds for i in Base.OneTo(nv)
        jointi = velocity_index_to_joint_id(state, i)
        bodyi = successorid(jointi, state)
        Ici = crb_inertia(state, bodyi, false)
        Si = motion_subspaces[i]
        Fi = Ici * Si
        for j in Base.OneTo(i)
            jointj = velocity_index_to_joint_id(state, j)
            M.data[i, j] = if supports(jointj, bodyi, state)
                Sj = motion_subspaces[j]
                (transpose(Fi) * Sj)[1]
            else
                zero(eltype(M))
            end
        end
    end
    M
end

mass_matrix!(result::DynamicsResult, state::MechanismState) = mass_matrix!(result.massmatrix, state)

"""
$(SIGNATURES)

$mass_matrix_doc
"""
function mass_matrix(state::MechanismState{X, M, C}) where {X, M, C}
    nv = num_velocities(state)
    ret = Symmetric(Matrix{C}(undef, nv, nv), :L)
    mass_matrix!(ret, state)
    ret
end

const momentum_matrix_doc = """Compute the momentum matrix ``A(q)`` of the
`Mechanism` in the given state.

The momentum matrix maps the `Mechanism`'s joint velocity vector ``v`` to
its total momentum.

See also [`MomentumMatrix`](@ref).
"""

const momentum_matrix!_doc = """$momentum_matrix_doc
The `out` argument must be a mutable `MomentumMatrix` with as many columns as
the dimension of the `Mechanism`'s joint velocity vector ``v``.
"""

"""
$(SIGNATURES)

$momentum_matrix!_doc

`transformfun` is a callable that may be used to transform the individual
momentum matrix blocks associated with each of the joints to the frame in which
`out` is expressed.

$noalloc_doc
"""
function momentum_matrix!(mat::MomentumMatrix, state::MechanismState, transformfun)
    nv = num_velocities(state)
    @boundscheck size(mat, 2) == nv || throw(DimensionMismatch())
    update_motion_subspaces!(state)
    update_crb_inertias!(state)
    @inbounds for i in 1 : nv
        jointi = velocity_index_to_joint_id(state, i)
        bodyi = successorid(jointi, state)
        Ici = crb_inertia(state, bodyi)
        Si = state.motion_subspaces.data[i]
        Fi = transformfun(Ici * Si)
        set_col!(mat, i, Fi)
    end
    mat
end

"""
$(SIGNATURES)

$momentum_matrix!_doc

`root_to_desired` is the transform from the `Mechanism`'s root frame to the frame
in which `out` is expressed.

$noalloc_doc
"""
function momentum_matrix!(mat::MomentumMatrix, state::MechanismState, root_to_desired::Transform3D)
    momentum_matrix!(mat, state, F -> transform(F, root_to_desired))
end

"""
$(SIGNATURES)

$momentum_matrix!_doc

See [`momentum_matrix!(out, state, root_to_desired)`](@ref). Uses `state`
to compute the transform from the `Mechanism`'s root frame to the frame in which
`out` is expressed.

$noalloc_doc
"""
function momentum_matrix!(out::MomentumMatrix, state::MechanismState)
    if out.frame == root_frame(state.mechanism)
        momentum_matrix!(out, state, identity)
    else
        momentum_matrix!(out, state, inv(transform_to_root(state, out.frame)))
    end
end

"""
$(SIGNATURES)

$momentum_matrix_doc

See [`momentum_matrix!(out, state)`](@ref).
"""
function momentum_matrix(state::MechanismState)
    ncols = num_velocities(state)
    T = cache_eltype(state)
    ret = MomentumMatrix(root_frame(state.mechanism), Matrix{T}(undef, 3, ncols), Matrix{T}(undef, 3, ncols))
    momentum_matrix!(ret, state, identity)
    ret
end

function bias_accelerations!(out::AbstractDict{BodyID, SpatialAcceleration{T}}, state::MechanismState{X, M}) where {T, X, M}
    update_bias_accelerations_wrt_world!(state)
    gravitybias = convert(SpatialAcceleration{T}, -gravitational_spatial_acceleration(state.mechanism))
    for jointid in state.treejointids
        bodyid = successorid(jointid, state)
        out[bodyid] = gravitybias + bias_acceleration(state, bodyid, false)
    end
    nothing
end

function spatial_accelerations!(out::AbstractDict{BodyID, SpatialAcceleration{T}}, state::MechanismState, v̇::SegmentedVector{JointID}) where T
    update_transforms!(state)
    update_twists_wrt_world!(state)

    # Compute joint accelerations
    joints = state.treejoints
    qs = values(segments(state.q))
    vs = values(segments(state.v))
    v̇s = values(segments(v̇))
    foreach_with_extra_args(state, out, joints, qs, vs, v̇s) do state, accels, joint, qjoint, vjoint, v̇joint
        bodyid = successorid(JointID(joint), state)
        accels[bodyid] = joint_spatial_acceleration(joint, qjoint, vjoint, v̇joint)
    end

    # Recursive propagation
    # TODO: manual frame changes. Find a way to not avoid the frame checks here.
    mechanism = state.mechanism
    root = root_body(mechanism)
    out[root] = convert(SpatialAcceleration{T}, -gravitational_spatial_acceleration(mechanism))
    @inbounds for jointid in state.treejointids
        parentbodyid, bodyid = predsucc(jointid, state)
        toroot = transform_to_root(state, bodyid, false)
        parenttwist = twist_wrt_world(state, parentbodyid, false)
        bodytwist = twist_wrt_world(state, bodyid, false)
        jointaccel = out[bodyid]
        jointaccel = SpatialAcceleration(jointaccel.body, parenttwist.body, toroot.to,
            Spatial.transform_spatial_motion(jointaccel.angular, jointaccel.linear, rotation(toroot), translation(toroot))...)
        out[bodyid] = out[parentbodyid] + (-bodytwist) × parenttwist + jointaccel
    end
    nothing
end

spatial_accelerations!(result::DynamicsResult, state::MechanismState) = spatial_accelerations!(result.accelerations, state, result.v̇)

function relative_acceleration(accels::AbstractDict{BodyID, SpatialAcceleration{T}}, body::RigidBody{M}, base::RigidBody{M}) where {T, M}
    -accels[BodyID(base)] + accels[BodyID(body)]
end

# TODO: ensure that accelerations are up-to-date
relative_acceleration(result::DynamicsResult, body::RigidBody, base::RigidBody) = relative_acceleration(result.accelerations, body, base)

function newton_euler!(
        out::AbstractDict{BodyID, Wrench{T}}, state::MechanismState{X, M},
        accelerations::AbstractDict{BodyID, SpatialAcceleration{T}},
        externalwrenches::AbstractDict{BodyID, Wrench{W}}) where {T, X, M, W}
    update_twists_wrt_world!(state)
    update_spatial_inertias!(state)
    for jointid in state.treejointids
        bodyid = successorid(jointid, state)
        wrench = newton_euler(state, bodyid, accelerations[bodyid], false)
        out[bodyid] = haskey(externalwrenches, bodyid) ? wrench - externalwrenches[bodyid] : wrench
    end
end

function joint_wrenches_and_torques!(
        torquesout::SegmentedVector{JointID},
        net_wrenches_in_joint_wrenches_out::AbstractDict{BodyID, <:Wrench}, # TODO: consider having a separate Associative{Joint{M}, Wrench{T}} for joint wrenches
        state::MechanismState)
    update_motion_subspaces!(state)
    # Note: pass in net wrenches as wrenches argument. wrenches argument is modified to be joint wrenches
    @boundscheck length(torquesout) == num_velocities(state) || error("length of torque vector is wrong")
    wrenches = net_wrenches_in_joint_wrenches_out
    for jointid in reverse(state.treejointids)
        parentbodyid, bodyid = predsucc(jointid, state)
        jointwrench = wrenches[bodyid]
        wrenches[parentbodyid] += jointwrench # action = -reaction
        for vindex in velocity_range(state, jointid)
            Scol = state.motion_subspaces.data[vindex]
            torquesout[vindex] = (transpose(Scol) * jointwrench)[1]
        end
    end
end

const dynamics_bias_doc = """Compute the 'dynamics bias term', i.e. the term
```math
c(q, v, w_\\text{ext})
```
in the unconstrained joint-space equations of motion
```math
M(q) \\dot{v} + c(q, v, w_\\text{ext}) = \\tau
```
given joint configuration vector ``q``, joint velocity vector ``v``,
joint acceleration vector ``\\dot{v}`` and (optionally) external
wrenches ``w_\\text{ext}``.

The `externalwrenches` argument can be used to specify additional
wrenches that act on the `Mechanism`'s bodies.
"""

"""
$(SIGNATURES)

$dynamics_bias_doc

$noalloc_doc
"""
function dynamics_bias!(
        torques::SegmentedVector{JointID},
        biasaccelerations::AbstractDict{BodyID, <:SpatialAcceleration},
        wrenches::AbstractDict{BodyID, <:Wrench},
        state::MechanismState{X},
        externalwrenches::AbstractDict{BodyID, <:Wrench} = NullDict{BodyID, Wrench{X}}()) where X
    bias_accelerations!(biasaccelerations, state)
    newton_euler!(wrenches, state, biasaccelerations, externalwrenches)
    joint_wrenches_and_torques!(torques, wrenches, state)
    torques
end

function dynamics_bias!(result::DynamicsResult, state::MechanismState)
    dynamics_bias!(result.dynamicsbias, result.accelerations, result.jointwrenches, state, result.totalwrenches)
end

"""
$(SIGNATURES)

$dynamics_bias_doc
"""
function dynamics_bias(
        state::MechanismState{X, M},
        externalwrenches::AbstractDict{BodyID, Wrench{W}} = NullDict{BodyID, Wrench{X}}()) where {X, M, W}
    T = promote_type(X, M, W)
    mechanism = state.mechanism
    torques = similar(velocity(state), T)
    rootframe = root_frame(mechanism)
    jointwrenches = BodyDict(BodyID(b) => zero(Wrench{T}, rootframe) for b in bodies(mechanism))
    accelerations = BodyDict(BodyID(b) => zero(SpatialAcceleration{T}, rootframe, rootframe, rootframe) for b in bodies(mechanism))
    dynamics_bias!(torques, accelerations, jointwrenches, state, externalwrenches)
    torques
end

const inverse_dynamics_doc = """Do inverse dynamics, i.e. compute ``\\tau``
in the unconstrained joint-space equations of motion
```math
M(q) \\dot{v} + c(q, v, w_\\text{ext}) = \\tau
```
given joint configuration vector ``q``, joint velocity vector ``v``,
joint acceleration vector ``\\dot{v}`` and (optionally) external
wrenches ``w_\\text{ext}``.

The `externalwrenches` argument can be used to specify additional
wrenches that act on the `Mechanism`'s bodies.

This method implements the recursive Newton-Euler algorithm.

Currently doesn't support `Mechanism`s with cycles.
"""

"""
$(SIGNATURES)

$inverse_dynamics_doc

$noalloc_doc
"""
function inverse_dynamics!(
        torquesout::SegmentedVector{JointID},
        jointwrenchesout::AbstractDict{BodyID, Wrench{T}},
        accelerations::AbstractDict{BodyID, SpatialAcceleration{T}},
        state::MechanismState,
        v̇::SegmentedVector{JointID},
        externalwrenches::AbstractDict{BodyID, <:Wrench} = NullDict{BodyID, Wrench{T}}()) where T
    @boundscheck length(tree_joints(state.mechanism)) == length(joints(state.mechanism)) || error("This method can currently only handle tree Mechanisms.")
    spatial_accelerations!(accelerations, state, v̇)
    newton_euler!(jointwrenchesout, state, accelerations, externalwrenches)
    joint_wrenches_and_torques!(torquesout, jointwrenchesout, state)
end

"""
$(SIGNATURES)

$inverse_dynamics_doc
"""
function inverse_dynamics(
        state::MechanismState{X, M},
        v̇::SegmentedVector{JointID, V},
        externalwrenches::AbstractDict{BodyID, Wrench{W}} = NullDict{BodyID, Wrench{X}}()) where {X, M, V, W}
    T = promote_type(X, M, V, W)
    mechanism = state.mechanism
    torques = similar(velocity(state), T)
    rootframe = root_frame(mechanism)
    jointwrenches = BodyDict(BodyID(b) => zero(Wrench{T}, rootframe) for b in bodies(mechanism))
    accelerations = BodyDict(BodyID(b) => zero(SpatialAcceleration{T}, rootframe, rootframe, rootframe) for b in bodies(mechanism))
    inverse_dynamics!(torques, jointwrenches, accelerations, state, v̇, externalwrenches)
    torques
end

function constraint_jacobian!(jac::AbstractMatrix, rowranges, state::MechanismState)
    # TODO: traversing jac in the wrong order
    update_motion_subspaces!(state)
    update_constraint_wrench_subspaces!(state)
    fill!(jac, 0)
    @inbounds for nontreejointid in state.nontreejointids
        path = state.constraint_jacobian_structure[nontreejointid]
        for cindex in constraint_range(state, nontreejointid)
            Tcol = state.constraint_wrench_subspaces.data[cindex]
            for i in eachindex(path.edges)
                treejoint = path.edges[i]
                vrange = velocity_range(state, treejoint)
                direction = directions(path)[i]
                sign = ifelse(direction == PathDirections.up, -1, 1)
                for col in eachindex(vrange)
                    vindex = vrange[col]
                    Scol = state.motion_subspaces.data[vindex]
                    jacelement = flipsign((transpose(Tcol) * Scol)[1], sign)
                    jac[cindex, vindex] = jacelement
                end
            end
        end
    end
    jac
end

function constraint_jacobian!(result::DynamicsResult, state::MechanismState)
    constraint_jacobian!(result.constraintjacobian, result.constraintrowranges, state)
end

"""
$(SIGNATURES)

Return the default Baumgarte constraint stabilization gains. These gains result in
critical damping, and correspond to ``T_{stab} = 0.1`` in Featherstone (2008), section 8.3.
"""
function default_constraint_stabilization_gains(scalar_type::Type{T}) where T
    ConstDict{JointID}(SE3PDGains(PDGains(T(100), T(20)), PDGains(T(100), T(20))))
end

const stabilization_gains_doc = """
The `stabilization_gains` keyword argument can be used to set PD gains for Baumgarte
stabilization, which can be used to prevent separation of non-tree (loop) joints.
See Featherstone (2008), section 8.3 for more information. There are several options
for specifying gains:

* `nothing` can be used to completely disable Baumgarte stabilization.
* Gains can be specifed on a per-joint basis using any `AbstractDict{JointID, <:RigidBodyDynamics.PDControl.SE3PDGains}`,
  which maps the `JointID` for the non-tree joints of the mechanism to the gains for that joint.
* As a special case of the second option, the same gains can be used for all joints
  by passing in a `RigidBodyDynamics.CustomCollections.ConstDict{JointID}`.

The [`default_constraint_stabilization_gains`](@ref) function is called to produce the default gains,
which use the last option.
"""

function constraint_bias!(bias::SegmentedVector, state::MechanismState{X};
        stabilization_gains::Union{AbstractDict{JointID, <:SE3PDGains}, Nothing}=default_constraint_stabilization_gains(X)) where X
    update_transforms!(state)
    update_twists_wrt_world!(state)
    update_bias_accelerations_wrt_world!(state)
    update_constraint_wrench_subspaces!(state)
    constraint_wrench_subspaces = state.constraint_wrench_subspaces.data
    for nontreejointid in state.nontreejointids
        predid, succid = predsucc(nontreejointid, state)

        predtwist = twist_wrt_world(state, predid, false)
        succtwist = twist_wrt_world(state, succid, false)
        crossterm = succtwist × predtwist

        succbias = bias_acceleration(state, succid, false)
        predbias = bias_acceleration(state, predid, false)
        jointbias = -predbias + succbias

        biasaccel = crossterm + jointbias # what's inside parentheses in 8.47 in Featherstone

        if stabilization_gains !== nothing
            # TODO: make this nicer (less manual frame juggling and no manual transformations)
            jointtransform = joint_transform(state, nontreejointid, false)
            jointtwist = -predtwist + succtwist
            jointtwist = Twist(jointtransform.from, jointtransform.to, jointtwist.frame, jointtwist.angular, jointtwist.linear) # make frames line up
            successor_to_root = transform_to_root(state, jointtransform.from) # TODO: expensive
            jointtwist = transform(jointtwist, inv(successor_to_root)) # twist needs to be in successor frame for pd method
            stabaccel = pd(stabilization_gains[nontreejointid], jointtransform, jointtwist, SE3PDMethod{:Linearized}()) # in successor frame
            stabaccel = SpatialAcceleration(stabaccel.body, stabaccel.base, biasaccel.frame,
                Spatial.transform_spatial_motion(stabaccel.angular, stabaccel.linear,
                rotation(successor_to_root), translation(successor_to_root))...) # back to world frame. TODO: ugly way to do this
            stabaccel = SpatialAcceleration(biasaccel.body, biasaccel.body, stabaccel.frame, stabaccel.angular, stabaccel.linear) # make frames line up
            @inbounds biasaccel = biasaccel + -stabaccel
        end

        for cindex in constraint_range(state, nontreejointid)
            Tcol = constraint_wrench_subspaces[cindex]
            # TODO: make nicer:
            @framecheck Tcol.frame biasaccel.frame
            bias[cindex] = (transpose(Tcol.angular) * biasaccel.angular + transpose(Tcol.linear) * biasaccel.linear)[1]
        end
    end
    bias
end

function constraint_bias!(result::DynamicsResult, state::MechanismState{X};
        stabilization_gains::Union{AbstractDict{JointID, <:SE3PDGains}, Nothing}=default_constraint_stabilization_gains(X)) where X
    constraint_bias!(result.constraintbias, state; stabilization_gains=stabilization_gains)
end

function contact_dynamics!(result::DynamicsResult{T, M}, state::MechanismState{X, M, C}) where {X, M, C, T}
    update_transforms!(state)
    update_twists_wrt_world!(state)
    mechanism = state.mechanism
    root = root_body(mechanism)
    frame = default_frame(root)
    for body in bodies(mechanism)
        bodyid = BodyID(body)
        wrench = zero(Wrench{T}, frame)
        points = contact_points(body)
        if !isempty(points)
            # TODO: AABB
            body_to_root = transform_to_root(state, bodyid, false)
            twist = twist_wrt_world(state, bodyid, false)
            states_for_body = contact_states(state, bodyid)
            state_derivs_for_body = contact_state_derivatives(result, bodyid)
            for i = 1 : length(points)
                @inbounds c = points[i]
                point = body_to_root * location(c)
                velocity = point_velocity(twist, point)
                states_for_point = states_for_body[i]
                state_derivs_for_point = state_derivs_for_body[i]
                for j = 1 : length(mechanism.environment.halfspaces)
                    primitive = mechanism.environment.halfspaces[j]
                    contact_state = states_for_point[j]
                    contact_state_deriv = state_derivs_for_point[j]
                    model = contact_model(c)
                    # TODO: would be good to move this to Contact module
                    # arguments: model, state, state_deriv, point, velocity, primitive
                    if point_inside(primitive, point)
                        separation, normal = Contact.detect_contact(primitive, point)
                        force = Contact.contact_dynamics!(contact_state_deriv, contact_state,
                            model, -separation, velocity, normal)
                        wrench += Wrench(point, force)
                    else
                        Contact.reset!(contact_state)
                        Contact.zero!(contact_state_deriv)
                    end
                end
            end
        end
        set_contact_wrench!(result, body, wrench)
    end
end

function dynamics_solve!(result::DynamicsResult, τ::AbstractVector)
    # version for general scalar types
    # TODO: make more efficient
    M = result.massmatrix
    c = parent(result.dynamicsbias)
    v̇ = parent(result.v̇)

    K = result.constraintjacobian
    k = result.constraintbias
    λ = result.λ

    nv = size(M, 1)
    nl = size(K, 1)
    G = [Matrix(M) K'; # TODO: Matrix because of https://github.com/JuliaLang/julia/issues/21332
         K zeros(nl, nl)]
    r = [τ - c; -k]
    v̇λ = G \ r
    v̇ .= view(v̇λ, 1 : nv)
    λ .= view(v̇λ, nv + 1 : nv + nl)
    nothing
end

function dynamics_solve!(result::DynamicsResult{T, S}, τ::AbstractVector{T}) where {S, T<:LinearAlgebra.BlasReal}
    # optimized version for BLAS floats
    M = result.massmatrix
    c = parent(result.dynamicsbias)
    v̇ = parent(result.v̇)

    K = result.constraintjacobian
    k = result.constraintbias
    λ = result.λ

    L = result.L
    A = result.A
    z = result.z
    Y = result.Y

    L .= M.data
    uplo = M.uplo
    LinearAlgebra.LAPACK.potrf!(uplo, L) # L <- Cholesky decomposition of M; M == L Lᵀ (note: Featherstone, page 151 uses M == Lᵀ L instead)
    τbiased = v̇
    τbiased .= τ .- c

    if size(K, 1) > 0
        # Loops.
        # Basic math:
        # DAE:
        #   [M Kᵀ] [v̇] = [τ - c]
        #   [K 0 ] [λ]  = [-k]
        # Solve for v̇:
        #   v̇ = M⁻¹ (τ - c - Kᵀ λ)
        # Plug into λ equation:
        #   K M⁻¹ (τ - c - Kᵀ λ) = -k
        #   K M⁻¹ Kᵀ λ = K M⁻¹(τ - c) + k
        # which can be solved for λ, after which λ can be used to solve for v̇.

        # Implementation loosely follows page 151 of Featherstone, Rigid Body Dynamics Algorithms.
        # It is slightly different because Featherstone uses M = Lᵀ L, whereas BLAS uses M = L Lᵀ.
        # M = L Lᵀ (Cholesky decomposition)
        # Y = K L⁻ᵀ, so that Y Yᵀ = K L⁻ᵀ L⁻¹ Kᵀ = K M⁻¹ Kᵀ = A
        # z = L⁻¹ (τ - c)
        # b = Y z + k
        # solve A λ = b for λ
        # solve M v̇  = τ - c for v̇

        # Compute Y = K L⁻ᵀ
        Y .= K
        LinearAlgebra.BLAS.trsm!('R', uplo, 'T', 'N', one(T), L, Y)

        # Compute z = L⁻¹ (τ - c)
        z .= τbiased
        LinearAlgebra.BLAS.trsv!(uplo, 'N', 'N', L, z) # z <- L⁻¹ (τ - c)

        # Compute A = Y Yᵀ == K * M⁻¹ * Kᵀ
        LinearAlgebra.BLAS.gemm!('N', 'T', one(T), Y, Y, zero(T), A) # A <- K * M⁻¹ * Kᵀ

        # Compute b = Y z + k
        b = λ
        b .= k
        LinearAlgebra.BLAS.gemv!('N', one(T), Y, z, one(T), b) # b <- Y z + k

        # Compute λ = A⁻¹ b == (K * M⁻¹ * Kᵀ)⁻¹ * (K * M⁻¹ * (τ - c) + k)
        # LinearAlgebra.LAPACK.posv!(uplo, A, b) # NOTE: doesn't work in general because A is only guaranteed to be positive semidefinite
        singular_value_zero_tolerance = 1e-10 # TODO: more principled choice
        # TODO: https://github.com/JuliaLang/julia/issues/22242
        b[:], _ = LinearAlgebra.LAPACK.gelsy!(A, b, singular_value_zero_tolerance) # b == λ <- (K * M⁻¹ * Kᵀ)⁻¹ * (K * M⁻¹ * (τ - c) + k)

        # Update τbiased: subtract Kᵀ * λ
        LinearAlgebra.BLAS.gemv!('T', -one(T), K, λ, one(T), τbiased) # τbiased <- τ - c - Kᵀ * λ

        # Solve for v̇ = M⁻¹ * (τ - c - Kᵀ * λ)
        LinearAlgebra.LAPACK.potrs!(uplo, L, τbiased) # τbiased ==v̇ <- M⁻¹ * (τ - c - Kᵀ * λ)
    else
        # No loops.
        LinearAlgebra.LAPACK.potrs!(uplo, L, τbiased) # τbiased == v̇ <- M⁻¹ * (τ - c)
    end
    nothing
end

"""
$(SIGNATURES)

Compute the joint acceleration vector ``\\dot{v}`` and Lagrange multipliers
``\\lambda`` that satisfy the joint-space equations of motion
```math
M(q) \\dot{v} + c(q, v, w_\\text{ext}) = \\tau - K(q)^{T} \\lambda
```
and the constraint equations
```math
K(q) \\dot{v} = -k
```
given joint configuration vector ``q``, joint velocity vector ``v``, and
(optionally) joint torques ``\\tau`` and external wrenches ``w_\\text{ext}``.

The `externalwrenches` argument can be used to specify additional
wrenches that act on the `Mechanism`'s bodies.

$stabilization_gains_doc
"""
function dynamics!(result::DynamicsResult, state::MechanismState{X},
        torques::AbstractVector = ConstVector(zero(X), num_velocities(state)),
        externalwrenches::AbstractDict{BodyID, <:Wrench} = NullDict{BodyID, Wrench{X}}();
        stabilization_gains=default_constraint_stabilization_gains(X)) where X
    configuration_derivative!(result.q̇, state)
    contact_dynamics!(result, state)
    for jointid in state.treejointids
        bodyid = successorid(jointid, state)
        contactwrench = result.contactwrenches[bodyid]
        result.totalwrenches[bodyid] = haskey(externalwrenches, bodyid) ? externalwrenches[bodyid] + contactwrench : contactwrench
    end
    dynamics_bias!(result, state)
    mass_matrix!(result, state)
    if has_loops(state.mechanism)
        constraint_jacobian!(result, state)
        constraint_bias!(result, state; stabilization_gains=stabilization_gains)
    end
    dynamics_solve!(result, parent(torques))
    nothing
end

"""
$(SIGNATURES)

Convenience function for use with standard ODE integrators that takes a
`Vector` argument
```math
x = \\left(\\begin{array}{c}
q\\\\
v
\\end{array}\\right)
```
and returns a `Vector` ``\\dot{x}``.
"""
function dynamics!(ẋ::StridedVector{X},
        result::DynamicsResult, state::MechanismState{X}, x::AbstractVector{X},
        torques::AbstractVector = ConstVector(zero(X), num_velocities(state)),
        externalwrenches::AbstractDict{BodyID, <:Wrench} = NullDict{BodyID, Wrench{X}}();
        stabilization_gains=default_constraint_stabilization_gains(X)) where X
    copyto!(state, x)
    dynamics!(result, state, torques, externalwrenches; stabilization_gains=stabilization_gains)
    copyto!(ẋ, result)
    ẋ
end


text:
# Items are assembled by the rss_generator in a global feed and sub-feeds
# for each of the tag. So each item is a tuple with the string of the item
# and
struct RSSItem
    item::String
    date::Date
    tags::Vector{String}
end

const RSS_ITEMS = Vector{RSSItem}()

"""
$SIGNATURES

If there's an RSS feed to generate, check that the template files are there, if not
then get them from FranklinTemplates. Also check that a few key variables
are defined.
"""
function prepare_for_rss()::Nothing
    # check for key variables
    key_vars = (:website_title, :website_url, :website_description)
    flag = any(isempty, globvar(key_var)::String for key_var in key_vars)

    # XXX [April 2020, 0.10.35] if the flag is false and generate_rss is false
    # switch it to true for backward compatibility; otherwise return
    gen_rss = globvar(:generate_rss)::Bool
    if flag
        gen_rss || return nothing
        print_warning("""
            RSS is set to be generated but the RSS feed is improperly described:
            at least one of the following variables have not been defined in
            your 'config.md': 'website_title', 'website_url', 'website_description'.
            The feed will not be (re)generated.
            \nRelevant pointer:
            $POINTER_PV
            """)
        return nothing
    elseif !gen_rss
        set_var!(GLOBAL_VARS, "generate_rss", true)
    end
    # check if there's an _rss folder, if there isn't generate one from
    # template (see FranklinTemplates)
    isdir(path(:rss)) || mkdir(path(:rss))
    for template in ("head.xml", "item.xml")
        dst = joinpath(path(:rss), template)
        if !isfile(dst)
            src = joinpath(
                dirname(pathof(FranklinTemplates)),
                "templates", "common", "_rss", template
            )
            cp(src, dst)
            # ensure the file is not read-only
            chmod(dst, 0o644)
        end
    end
    return nothing
end

"""
$SIGNATURES

Extract the entries from RSS_ITEMS and assemble the RSS. If RSS_ITEMS is empty, nothing
is generated. If tags are associated with items, also create one feed per tag.

Note: if a file `head_rtag.xml` is available, where `rtag` is the refstring of a tag, it
will take precedence for the corresponding tag feed; otherwise `head.xml` will be used.
"""
function rss_generator()::Nothing
    # is there anything to go in the RSS feed?
    isempty(RSS_ITEMS) && return nothing

    # sort items by reverse chronological order
    sorted_items = sort(RSS_ITEMS, rev=true, by=x->x.date)

    feed_name = globvar(:rss_file)::String * ".xml"
    global_feed_path = joinpath(PATHS[:site], feed_name)
    global_feed_head = replace(
        read(joinpath(path(:rss), "head.xml"), String),
        r"<!--(.|\n)*?-->" => ""
    )
    if FD_ENV[:FINAL_PASS]
        pp = globvar(:prepath)::String
        isempty(pp) || (pp = "/$pp")
        global_feed_head = replace(
            global_feed_head,
            r"(\/[a-zA-Z0-9\_-]+\.xsl)" => SubstitutionString("$pp\\1")
        )
    end

    open(global_feed_path, "w") do io
        write(io, convert_html(global_feed_head))
        for item in sorted_items
            write(io, item.item)
        end
        write(io, "</channel></rss>")
    end

    # Tag specific feed; filter items by tag
    # > Collect all raw tags
    tags = Set{String}()
    foreach(x -> union!(tags, x.tags), RSS_ITEMS)
    # > write a feed per tag, use head_tag if it exists
    for tag in tags
        rtag = refstring(tag)
        tag_dir = joinpath(path(:tag), refstring(tag))
        isdir(tag_dir) || mkpath(tag_dir)
        tag_feed_path = joinpath(tag_dir, feed_name)
        open(tag_feed_path, "w") do io
            tag_feed_head = global_feed_head
            tag_feed_head_cand = joinpath(path(:rss), "head_$rtag.xml")
            if isfile(tag_feed_head_cand)
                tag_feed_head = read(tag_feed_head_cand)
            end
            write(io, convert_html(tag_feed_head))
            for item in sorted_items
                tag ∈ item.tags || continue
                write(io, item.item)
            end
            write(io, "</channel></rss>")
        end
    end
    return nothing
end

"""
$SIGNATURES

This is called in `convert_and_write` to add an RSS item.
"""
function add_rss_item()
    item_template = replace(
        read(joinpath(path(:rss), "item.xml"), String),
        r"<!--(.|\n)*?-->" => ""
    )
    # if the rss_title is not given, infer from title
    if isempty(locvar(:rss_title)::String)
        set_var!(LOCAL_VARS, "rss_title", locvar(:title)::String)
    end
    item = replace(convert_html(item_template), r"\n\n" => "")
    item = replace(item, Regex(raw"""<a\shref=(?:"|')?\#.*?>(.*?)</a>""") => s"\1")

    rss_item = RSSItem(
        item,
        locvar(:rss_pubdate)::Date,
        locvar(:tags)::Vector{String}
    )
    push!(RSS_ITEMS, rss_item)
end

ext: jl


text:
function BracketingNthOrderBrentSolver()
    return BracketingNthOrderBrentSolver(())
end

function BracketingNthOrderBrentSolver(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jint)
    return BracketingNthOrderBrentSolver((jdouble, jdouble, jdouble, jint), arg0, arg1, arg2, arg3)
end

function BracketingNthOrderBrentSolver(arg0::jdouble, arg1::jdouble, arg2::jint)
    return BracketingNthOrderBrentSolver((jdouble, jdouble, jint), arg0, arg1, arg2)
end

function BracketingNthOrderBrentSolver(arg0::jdouble, arg1::jint)
    return BracketingNthOrderBrentSolver((jdouble, jint), arg0, arg1)
end

function get_maximal_order(obj::BracketingNthOrderBrentSolver)
    return jcall(obj, "getMaximalOrder", jint, ())
end

function solve(obj::BracketingNthOrderBrentSolver, arg0::jint, arg1::UnivariateFunction, arg2::jdouble, arg3::jdouble, arg4::AllowedSolution)
    return jcall(obj, "solve", jdouble, (jint, UnivariateFunction, jdouble, jdouble, AllowedSolution), arg0, arg1, arg2, arg3, arg4)
end

function solve(obj::BracketingNthOrderBrentSolver, arg0::jint, arg1::UnivariateFunction, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::AllowedSolution)
    return jcall(obj, "solve", jdouble, (jint, UnivariateFunction, jdouble, jdouble, jdouble, AllowedSolution), arg0, arg1, arg2, arg3, arg4, arg5)
end

function solve_interval(obj::BracketingNthOrderBrentSolver, arg0::jint, arg1::UnivariateFunction, arg2::jdouble, arg3::jdouble, arg4::jdouble)
    return jcall(obj, "solveInterval", BracketedUnivariateSolver_Interval, (jint, UnivariateFunction, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4)
end


text:
# Purpose of this file:
# Command-line driven test case for octsock5 high-speed data interface
# Usually, this will be run as two independent processes for server / client interaction.
# Command line arguments "client", "server" define, which side is covered by the process (both is possible).
# Testcases are deterministic (srand), therefore both processes will run in sync
using octsock5;

function usererror(msg::String)
    error("\n\n#####" * msg * "#####\n");
end

function measureRoundtripTime(iOs)
    t::UInt64 = time_ns();
    nRuns::UInt64 = 300000;
    for ix = 1 : nRuns
        octsock5_write(iOs, ix);
        octsock5_read(iOs);
    end

    t = time_ns() - t;
    print("Average roundtrip time ", @sprintf("%1.3f", Float64(t)/1e6/Float64(nRuns)), " ms\n");
end

function measureThroughput(iOs::octsock5_cl)
    nRuns::UInt64 = 3000;
    m = Array{Float64}(rand(60, 100, 40));

    # === make first call without timing ===
    t::UInt64 = time_ns();
    octsock5_write(iOs, m);
    obj = octsock5_read(iOs);

    for ix = 1 : nRuns
        octsock5_write(iOs, m);
        octsock5_read(iOs);
    end

    t_s = Float64(time_ns() - t)/1e9;
    print("Average throughput ", @sprintf("%1.1f", Float64(nRuns*sizeof(m)/2^20)/t_s), " MBytes (2^20) per second, round-trip (one-way: about 2x)\n");
end

function testSpecials(iOs::octsock5_cl)
    arg::Tuple = (Inf, -Inf, NaN);
    octsock5_write(iOs, arg);
    obj::Tuple = octsock5_read(iOs);
    if (obj !== arg)
        error("verify fail in Inf/-Inf/NaN");
    end
    print("OK\n");
end

function testLarge(iOs::octsock5_cl)
    a1::Array{Float64} = rand(1000000, 10);
    octsock5_write(iOs, a1);
    a2::Array{Float64} = octsock5_read(iOs);
    if (a1 != a2)
        error("verify fail in large test");
    end
    print("OK\n");
end

function testAllTypes(iOs::octsock5_cl, nRuns::Int, profiling::Bool)
    # === build data ===
    arg0 = sin.(1.1:10.1);
    arg0f = Array{Float32}(arg0);
    v_UInt8 = Array{UInt8}(0:127);
    v_Int8 = Array{Int8}(0:127);
    v_UInt16 = Array{UInt16}(0:127);
    v_Int16 = Array{Int16}(0:127);
    v_UInt32 = Array{UInt32}(0:127);
    v_Int32 = Array{Int32}(0:127);
    v_UInt64 = Array{UInt64}(0:127);
    v_Int64 = Array{Int64}(0:127);
    m_2 = Array{Float64}(rand(3, 3));
    m_3 = Array{Float64}(rand(3, 3, 3));
    d = Dict();
    d[1] = "two";
    d["three"] = 4;

    # === run tests ===
    arg = Void;
    for v = 0 : nRuns-1
        vv = v+13;

        vec = [mod(v, 127)+im*mod(v, 126),
               mod(v, 125)+im*mod(v, 124),
               mod(v, 123)+im*mod(v, 122), ];
        if arg == Void
            arg = (
                # UInt, Int, float scalar
                UInt8(mod(v, 256)),
                UInt16(mod(v, 65536)),
                UInt32(v),
                UInt64(v),
                Int8(mod(v, 256)-128),
                Int16(mod(v, 65536)-32768),
                Int32(v),
                Int64(v),
                Float32(1.01*v),
                Float64(1.01*v),

                # UInt, Int, float scalar complex
                UInt8(mod(v, 256))+im*UInt8(mod(vv, 256)),
                UInt16(mod(v, 65536))+im*UInt16(mod(vv, 65536)),
                UInt32(v)+im*UInt32(vv),
                UInt64(v)+im*UInt64(vv),
                Int8(mod(v, 256)-128)+im*Int8(mod(vv, 256)-128),
                Int16(mod(v, 65536)-32768)+im*Int16(mod(vv, 65536)-32768),
                Int32(v)+im*Int32(v),
                Int64(v)+im*Int64(v),
                Float32(1.01*v)+im*Float32(1.02*vv),
                Float64(1.01*v)+im*Float64(1.02*vv),

                # UInt, int, float Array
                v_UInt8, v_Int8,
                v_UInt16, v_Int16,
                v_UInt32, v_Int32,
                v_UInt64, v_Int64,
                arg0, arg0f,

                # complex vector
                Array{Complex{UInt8}}(vec),
                Array{Complex{UInt16}}(vec),
                Array{Complex{UInt32}}(vec),
                Array{Complex{UInt64}}(vec),
                Array{Complex{Int8}}(vec),
                Array{Complex{Int16}}(vec),
                Array{Complex{Int32}}(vec),
                Array{Complex{Int64}}(vec),

                # n-d matrix
                m_2, m_3,

                # complex matrix
                m_2-2.5*im*m_2, m_3-2.5*im*m_3,

                # string
                "Hello world. The quick brown fox jumps over the lazy dog. This sentence no verb. Lorem Ipsum!",

                # Dict
                d,
            );

            #arg = Array{Complex{UInt64}}(1:3)
            #arg = Array{Complex{Float64}}(rand(3)+1im*rand(3))
            #arg = Array{Complex{UInt64}}(floor.((2.^50)*rand(10000)) + 1im*floor.((2.^50)*rand(10000)));
        end

        octsock5_write(iOs, arg);
        obj = octsock5_read(iOs);
        if (profiling == false)
            if (obj != arg)
                print("reference:", arg, "\n");
                print("received:", obj, "\n");
                error("verify fail");
            end
        end
        # if profiling == false && (mod(v, 1000) == 0) info(v, "\n"); end
        if profiling == false
            arg = Void;
        end
    end
    if (!profiling)
        print("no errors detected");
    end
end

function main()
    srand(0);

    args::Dict{String, Bool} = Dict(
        "server" => false,
        "client" => false,
        "roundtrip" => false,
        "throughput" => false,
        "profiling" => false,
        "alltypes" => false,
        "tcpip" => false,
        "specials" => false,
        "helloworld" => false,
        "large" => false
    );

    for arg::String in ARGS
        if haskey(args, arg)
            args[arg] = true;
        else
            usererror("invalid command line argument '" * arg * "'");
        end
    end

    # negative port number: Use Windows named pipes
    # positive port number: Use TCP/IP port 2000
    portNum::Int64 = args["tcpip"] ? 2000 : -12345

    if (args["server"] && args["client"])
        usererror("must use server or client argument"); # avoid lockup
    elseif (args["server"])
        for arg in ARGS
            if (!(arg == "server" || arg == "tcpip"))
                error("argument unsupported in server mode: " * arg);
            end
        end
        # === open link ===
        iOs::octsock5_cl = octsock5_new(isServer=true, portNum=portNum);
        # Agreed arbitrary token via STDOUT to guarantee that the server is up when the client is started
        # From the command line, start the client process only when this line has appeared
        # An automated startup would e.g. use redirected STDOUT or simply retry with timeout
        print("SERVER_READY\n");

        # Execution blocks here until connected
        octsock5_accept(iOs);
        print("got connection\n");

        # === run loopback server ===
        while (true)
            tmp = octsock5_read(iOs);
            octsock5_write(iOs, tmp);
            # agreed arbitrary token to stop the loopback server
            if tmp == "end loopback and have a nice day" break; end
        end
        octsock5_delete(iOs);
        print("SERVER_EXIT\n");
        return;
    elseif (args["client"])

        # === open link ===
        iOs = octsock5_new(isServer=false, portNum=portNum);

        # === run tests ===
        if (args["helloworld"])
            octsock5_write(iOs, "Hello World");
            res::String = octsock5_read(iOs); assert(res == "Hello World");
        end

        if (args["roundtrip"])
            measureRoundtripTime(iOs);
        end

        if (args["throughput"])
            measureThroughput(iOs);
        end

        if (args["alltypes"])
            testAllTypes(iOs, 10000, args["profiling"]);
        end

        if (args["specials"])
            testSpecials(iOs);
        end

        if (args["large"])
            testLarge(iOs);
        end

        # agreed arbitrary token to stop loopback server
        octsock5_write(iOs, "end loopback and have a nice day");
        tmpStr::String = octsock5_read(iOs);
        assert(tmpStr == "end loopback and have a nice day");
        octsock5_delete(iOs);
    end # elseif client
end
main()


text:
module FastLegendre

# From
# O(1) Computation of Legendre polynomials and Gauss-Legendre nodes and weights for parallel computing
# Bogaert, Ignace and Michiels, Bart and Fostier, Jan
# SIAM JOURNAL ON SCIENTIFIC COMPUTING
# 34, 2012
# http://dx.doi.org/10.1137/110855442
#
# as mentioned https://github.com/JuliaMath/SpecialFunctions.jl/issues/124
#
# Only Float64
#
# Gives a speed up over Clenshaw, which is expected, as this is O(1), not O(n):
#
# n= 1000 Clenshaw might be faster
# n= 10_000 10x faster
# n= 100_000 10x faster
# n= 1_000_000 500x faster
# n= 10_000_000 10,000x faster
# (measurement of basis(P,n)(x) compared to Basis(P,n)(x)

using FastGaussQuadrature
using SpecialFunctions

# 4 regions of evaluation
# specialized for Float64 values
function  fastlegendre(n, x)
    n == 0 && return one(x)
    n == 1 && return x
    n < 100 &&  return  direct_32(n,x)
    x < 0  &&  return  (iseven(n) ? 1 : -1) * fastlegendre(n,-x)
    θ = acos(x)
    (n+1)*sin(θ) > 25 && return asy_32(n,x)
    asy_33(n,x)
end

## ---- n < 100 case

# precompute first 99
const xᵢs = [FastGaussQuadrature.gausslegendre(n)[1][1:(n÷2)] for n in 1:99]
const xᵢs⁻¹ = [inv.(xᵢ) for xᵢ in xᵢs]
const xᵢs2⁻¹ = [inv.(1 .- xᵢ.^2) for xᵢ in xᵢs]

# precompute
const Cn = zeros(Float64, 99)
Cn[1] = 1.0
for i in 1:49
    n = 2i
    Cn[n]  = sqrt(pi) / gamma(i + 1)  / gamma((1-n)/2)
    Cn[n+1] = (n+1) * Cn[n]
end

# use cosine or sine
function  direct_32(n,x)
    if abs(x) > 0.707
        if iseven(n)
            direct_sine(Val(true), n,x)
        else
            direct_sine(Val(false), n,x)
        end
    else
        if iseven(n)
            direct_cosine(Val(true), n,x)
        else
            direct_cosine(Val(false), n,x)
        end
    end
end

#x = cos(theta) 1 - x^2 = sintheta0
function direct_sine(::Val{true}, n,x)
    tot = one(Float64)
    xs = xᵢs2⁻¹[n]
    for aᵢ in xs
        tot *= 1 - (1-x^2)*aᵢ
    end
    tot
end

direct_sine(::Val{false}, n,x) =  x * direct_sine(Val(true),n,x)

function direct_cosine(p::Val{true}, n,x)
    tot =  Cn[n]
    xs = xᵢs⁻¹[n]
    for aᵢ in  xs
        tot *= (1 - (x*aᵢ)^2)
    end
    tot
end
direct_cosine(p::Val{false}, n,x) = x * direct_cosine(Val(true), n, x)

## ----
const Cl0s = [gamma(n+1)/gamma(n + 3/2) for n in 1:9]
τ(x) = evalpoly(1/x, (1/1, 0, -1/64, 0, 21/8192, 0, -671/524288, 0, 180323/80323134, 0,
                    -20898423/8589934592, - 7426362705/1099511627776))

function Cl0(l)
    1 <= l < 10 && return Cl0s[l]
    1/sqrt(l+3/4)*τ(l+3/4)
end

function asy_32(n,x)
    M = 17
    θ = acos(x)
    val = sqrt((2/pi)*(1/sin(θ)))
    tot = 0.0
    sθ = sin(θ)
    λ = Cl0(n)
    for m in 0:M-1
        α = (n+m+1/2)*θ - (m+1/2)*pi/2
        tot += λ * cos(α)
        λ *= (m+1/2)^2/2/(m+1)/(n+m+3/2) / sθ
    end
    val * tot
end

## ---
h(n,y) = y^n * besselj(n,y)

const cns = [[1/8, -1/12],
      [11/384, -7/160, 1/160],
      [173/15360, -101/3584, 671/80640, -61/120960],
      [22931/3440640,  -90497/3870720, 217/20480, -1261/967680, 1261/29030400],
      [1319183/247726080, -10918993/454164480, +1676287/113541120, -7034857/2554675200,
       +1501/8110080, -79/20275200],
      [233526463/43599790080, -1396004969/47233105920, +2323237523/101213798400,
       -72836747/12651724800, +3135577/5367398400, -1532789/61993451520, +66643/185980354560]
      ]

function fn(i,y)
    n = 2i
    cs = cns[i]
    tot = 0.0
    Δ = i-1
    for j in i:2i
        tot += cs[j-Δ]*h(j,y)
    end
    tot
end

function asy_33(n,x)
    v = n + 1/2
    θ = acos(x)
    y = v * θ

    # Eq (3.19)
    val = besselj(0,y)

    λ = 1/v^2
    for i in 1:6
        val += fn(i,y)*λ
        λ /= v^2
    end
    val

end

end


text:
import ControlSystems: feedback

"""
    FRD(w,r)

Represents frequency-response data. `w` holds the frequency vector and `r` the response. Methods defined on this type include
- `+-*`
- `length, vec, sqrt`
- `plot`
- `feedback`
"""
struct FRD{WT<:AbstractVector,RT<:AbstractVector} <: LTISystem
    w::WT
    r::RT
end

struct Hz <: Number
    i::Any
end
Base.:*(i, ::Type{Hz}) = Hz(i)
struct rad <: Number
    i::Any
end
Base.:*(i, ::Type{rad}) = rad(i)
(::Colon)(start::Union{Hz,rad}, stop::Union{Hz,rad}) = (start, stop)

import Base: +, -, *, /, length, sqrt, getindex
FRD(w, s::LTISystem) = FRD(w, freqresp(s, w)[:, 1, 1])
Base.vec(f::FRD) = f.r
*(f::FRD, f2)    = FRD(f.w, f.r .* vec(f2))
+(f::FRD, f2)    = FRD(f.w, f.r .+ vec(f2))
-(f::FRD, f2)    = FRD(f.w, f.r .- vec(f2))
/(f::FRD, f2)    = FRD(f.w, f.r ./ vec(f2))

# Below are required for ambiguities
*(f::FRD, f2::FRD) = FRD(f.w, f.r .* f2.r)
+(f::FRD, f2::FRD) = FRD(f.w, f.r .+ f2.r)
-(f::FRD, f2::FRD) = FRD(f.w, f.r .- f2.r)
/(f::FRD, f2::FRD) = FRD(f.w, f.r ./ f2.r)

*(f::FRD, f2::LTISystem) = *(f, FRD(f.w, f2))
+(f::FRD, f2::LTISystem) = +(f, FRD(f.w, f2))
-(f::FRD, f2::LTISystem) = -(f, FRD(f.w, f2))

-(f::FRD) = FRD(f.w, -f.r)
length(f::FRD) = length(f.w)
Base.size(f::FRD) = (1, 1) # Size in the ControlSystems sense
Base.lastindex(f::FRD) = length(f)
function Base.getproperty(f::FRD, s::Symbol)
    s === :Ts && return 1 / ((f.w[2] - f.w[2]) / (2π))
    getfield(f, s)
end
Base.propertynames(f::FRD, private::Bool = false) = (fieldnames(typeof(f))..., :Ts)
ControlSystems.noutputs(f::FRD) = 1
ControlSystems.ninputs(f::FRD) = 1

sqrt(f::FRD) = FRD(f.w, sqrt.(f.r))
getindex(f::FRD, i) = FRD(f.w[i], f.r[i])
getindex(f::FRD, i::Int) = f.r[i]
getindex(f::FRD, i::Int, j::Int) = (@assert(i == 1 && j == 1); f)
function getindex(f::FRD, r::Tuple{Hz,Hz})
    s = findfirst(2pi * r[1].i .< f.w)
    e = findlast(f.w .< 2pi * r[2].i)
    f[s:e]
end
function getindex(f::FRD, r::Tuple{rad,rad})
    s = findfirst(r[1].i .< f.w)
    e = findlast(f.w .< r[2].i)
    f[s:e]
end
Base.isapprox(f1::FRD, f2::FRD; kwargs...) =
    (f1.w == f2.w) && isapprox(f1.r, f2.r; kwargs...)
Base.:(==)(f1::FRD, f2::FRD) = (f1.w == f2.w) && ==(f1.r, f2.r)

sensitivity(P::FRD, K) = FRD(P.w, 1.0 ./ (1.0 .+ vec(P) .* vec(K)))
feedback(P::FRD, K) = FRD(P.w, vec(P) ./ (1.0 .+ vec(P) .* vec(K)))
feedback(P, K::FRD) = FRD(K.w, vec(P) ./ (1.0 .+ vec(P) .* vec(K)))
feedback(P::FRD, K::FRD) = feedback(P, vec(K))

feedback(P::FRD, K::LTISystem) = feedback(P, freqresp(K, P.w)[:, 1, 1])
feedback(P::LTISystem, K::FRD) = feedback(freqresp(P, K.w)[:, 1, 1], K)

@recipe function plot_frd(frd::FRD; hz = false, plotphase=false)
    yscale --> :log10
    xscale --> :log10
    xguide --> (hz ? "Frequency [Hz]" : "Frequency [rad/s]")
    yguide --> "Magnitude"
    title --> "Bode Plot"
    legend --> false
    layout --> (plotphase ? 2 : 1)
    @series begin
        inds = findall(x -> x == 0, frd.w)
        subplot --> 1
        useinds = setdiff(1:length(frd.w), inds)
        (hz ? 1 / (2π) : 1) .* frd.w[useinds], abs.(frd.r[useinds])
    end
    if plotphase
        @series begin
            inds = findall(x -> x == 0, frd.w)
            subplot --> 2
            useinds = setdiff(1:length(frd.w), inds)
            (hz ? 1 / (2π) : 1) .* frd.w[useinds], unwrap(angle.(frd.r[useinds]))
        end
    end
    nothing
end

freqvec(h, k) = LinRange(0, π / h, length(k))

"""
    H, N = tfest(data, σ = 0.05)

Estimate a transfer function model using the Correlogram approach.
    Both `H` and `N` are of type `FRD` (frequency-response data).

`σ` determines the width of the Gaussian window applied to the estimated correlation functions before FFT. A larger `σ` implies less smoothing.
- `H` = Syu/Suu             Process transfer function
- `N` = Sy - |Syu|²/Suu     Noise PSD
"""
function tfest(d, σ::Real = 0.05)
    y, u, h = time1(output(d)), time1(input(d)), sampletime(d)
    Syy, Suu, Syu = fft_corr(y, u, σ)
    w = freqvec(h, Syu)
    H = FRD(w, Syu ./ Suu)
    N = FRD(w, @.(Syy - abs2(Syu) / Suu) ./ length(y))
    return H, N
end

function fft_corr(y, u, σ = 0.05)
    n = length(y)
    w = gaussian(2n - 1, σ)

    Syu = rfft(ifftshift(w .* xcorr(y, u)))
    Syy = rfft(ifftshift(w .* xcorr(y, y)))
    Suu = rfft(ifftshift(w .* xcorr(u, u)))
    Syy, Suu, Syu
end

"""
    κ = coherence(d; n = length(d)÷10, noverlap = n÷2, window=hamming)

Calculates the magnitude-squared coherence Function. κ close to 1 indicates a good explainability of energy in the output signal by energy in the input signal. κ << 1 indicates that either the system is nonlinear, or a strong noise contributes to the output energy.
κ: Coherence function (not squared)
N: Noise model
"""
function coherence(d; n = length(d) ÷ 10, noverlap = n ÷ 2, window = hamming)
    y, u, h = time1(output(d)), time1(input(d)), sampletime(d)
    Syy, Suu, Syu = wcfft(y, u, n = n, noverlap = noverlap, window = window)
    k = (abs2.(Syu) ./ (Suu .* Syy))#[end÷2+1:end]
    Sch = FRD(freqvec(h, k), k)
    return Sch
end

function wcfft(y, u; n = length(y) ÷ 10, noverlap = n ÷ 2, window = hamming)
    win, norm2 = DSP.Periodograms.compute_window(window, n)
    uw = arraysplit(u, n, noverlap, nextfastfft(n), win)
    yw = arraysplit(y, n, noverlap, nextfastfft(n), win)
    Syy = zeros(length(uw[1]) ÷ 2 + 1)
    Suu = zeros(length(uw[1]) ÷ 2 + 1)
    Syu = zeros(ComplexF64, length(uw[1]) ÷ 2 + 1)
    for i in eachindex(uw)
        xy = rfft(yw[i])
        xu = rfft(uw[i])
        Syu .+= xy .* conj.(xu)
        Syy .+= abs2.(xy)
        Suu .+= abs2.(xu)
    end
    Syy, Suu, Syu
end

@userplot Coherenceplot
"""
    coherenceplot(d, [(;n=..., noverlap=...); hz=false)

Calculates and plots the (squared) coherence Function κ. κ close to 1 indicates a good explainability of energy in the output signal by energy in the input signal. κ << 1 indicates that either the system is nonlinear, or a strong noise contributes to the output energy.

`hz` indicates Hertz instead of rad/s

Keyword arguments to `coherence` are supplied as a named tuple as a second positional argument .
"""
coherenceplot

@recipe function coherenceplot(p::Coherenceplot; hz = false)
    ae =
        ArgumentError("Call like this: coherenceplot(iddata; hz=false)")
    d = p.args[1]
    d isa AbstractIdData || throw(ae)
    if length(p.args) >= 2
        kwargs = p.args[2]
    else
        kwargs = NamedTuple()
    end
    y, u, h = output(d), input(d), sampletime(d)
    yscale --> :identity
    xscale --> :log10
    ylims --> (0, 1)
    xguide --> (hz ? "Frequency [Hz]" : "Frequency [rad/s]")
    title --> "Coherence"
    legend --> false
    frd = coherence(d; kwargs...)
    @series begin
        inds = findall(x -> x == 0, frd.w)
        useinds = setdiff(1:length(frd.w), inds)
        (hz ? 1 / (2π) : 1) .* frd.w[useinds], abs.(frd.r[useinds])
    end
    nothing
end

"""
    ir, t, Σ = impulseest(d::AbstractIdData, n; λ=0, estimator=ls)

Estimates the system impulse response by fitting an `n`:th order FIR model. Returns impulse-response estimate, time vector and covariance matrix.
See also `impulseestplot`
"""
function impulseest(d::AbstractIdData, n; λ = 0, estimator = ls)
    h, y, u = d.Ts, time1(output(d)), time1(input(d))
    N = min(length(u), length(y))
    @views yt, A = getARXregressor(y[1:N], u[1:N], 0, n)
    ir = estimator(A, yt, λ)
    t = range(h, length = n, step = h)
    Σ = parameter_covariance(yt, A, ir, λ)
    ir ./ h, t, Σ
end

@userplot Impulseestplot

"""
    impulseestplot(data,n)

Estimates the system impulse response by fitting an `n`:th order FIR model and plots the result with a 95% confidence band.
See also `impulseestplot`
"""
impulseestplot
@recipe function impulseestplot(p::Impulseestplot; λ = 0)
    d = p.args[1]
    n = length(p.args) >= 2 ? p.args[2] : 25
    ir, t, Σ = impulseest(d, n; λ = λ)
    title --> "Estimated Impulse Response"
    xguide --> "Time [s]"

    @series begin
        label --> ""
        t, ir
    end
    linestyle := :dash
    color := :black
    label := ""
    seriestype := :hline
    @series begin
        t, 2 .* sqrt.(diag(Σ))
    end
    @series begin
        t, -2 .* sqrt.(diag(Σ))
    end
end

@userplot Crosscorplot

"""
    crosscorplot(data, [lags])

Plot the cross correlation betweein input and output for `lags` that default to 10% of the length of the dataset on the negative side and 50% on the positive side but no more than 100 on each side.
"""
crosscorplot

@recipe function crosscorplot(p::Crosscorplot)
    d = p.args[1]
    N = length(d)
    lags = length(p.args) >= 2 ? p.args[2] : -max(N ÷ 10, 100):max(N ÷ 2, 100)
    xc = crosscor(time1(d.u), time1(d.y), lags, demean = true)
    title --> "Input-Output cross correlation"
    xguide --> "Lag [s]"

    @series begin
        seriestype --> :sticks
        label --> ""
        lags .* d.Ts, xc
    end
    linestyle := :dash

    seriescolor := :black
    label := ""
    primary := false
    # Ni = N .- abs.(lags)
    @series begin
        seriestype := :hline
        # lags.*d.Ts, 2 .*sqrt.(1 ./ Ni) # The denominator in crosscorr already takes care of this
        [2 .* sqrt.(1 ./ N)]
    end
    @series begin
        seriestype := :hline
        # lags.*d.Ts, -2 .*sqrt.(1 ./ Ni)
        [-2 .* sqrt.(1 ./ N)]
    end
end

function ControlSystems.gangoffour(P::FRD, C::FRD, ω = nothing)
    ω === nothing || ω == P.ω || error("Incosistent frequency vectors")
    S = sensitivity(P, C)
    D = (P * S)
    N = (C * S)
    T = (P * N)
    return S, D, N, T
end


text:
module MCMCTempering

import AbstractMCMC
import Distributions
import Random

include("adaptation.jl")
include("tempered.jl")
include("ladders.jl")
include("stepping.jl")
include("model.jl")
include("swapping.jl")
include("plotting.jl")

export tempered, TemperedSampler, plot_swaps, plot_ladders, make_tempered_model, get_tempered_loglikelihoods_and_params, make_tempered_loglikelihood, get_params

function AbstractMCMC.bundle_samples(
    ts::Vector,
    model::AbstractMCMC.AbstractModel,
    sampler::TemperedSampler,
    state::TemperedState,
    chain_type::Type;
    kwargs...
)
    AbstractMCMC.bundle_samples(ts, model, sampler.internal_sampler, state, chain_type; kwargs...)
end

end


text:
function dirichlet!(rng::AbstractRNG, probabilities::AbstractVector{T}, α::NTuple{N,T}) where {N,T}
    cumulative_γ = zero(eltype(α))
    @inbounds for i ∈ eachindex(α)
        γ = randgamma(rng, α[i])
        cumulative_γ += γ
        probabilities[i] = γ
    end
    inv_cumulative_γ = 1 / cumulative_γ
    @inbounds for i ∈ eachindex(α)
        probabilities[i] *= inv_cumulative_γ
    end
end

@inline update_probabilities!(probabilities, α) = update_probabilities!(Random.GLOBAL_RNG, probabilities, α)

"""
This function sets the priors of each successive group.
Each Wishart3x3's scale is set to the identity matrix, with 5 degrees of freedom.
The prior probability of group membership is set to 1.5^(number of groups - group number).
"""
@generated function set_priors!(Λs::AbstractVector{InverseWishart{T}}, ::Val{NG}) where {NG, T}
    quote
        $(Expr(:meta, :inline))
        @inbounds @nexprs $NG g ->
            Λs[g] = InverseWishart(one(T), zero(T), zero(T), one(T), zero(T), one(T), T(3), T($(1/NG))*T(1.5)^($NG-g))
        nothing
    end
end

"""
Calculates the Wishart posteriors, conditional on group membership assignments.
It does this via first resetting them to the priors, summing up based on group membership,
and then finally also calculating the corresponding Cholesky factor and inverse.
"""
function calc_Wisharts!(revcholwisharts::AbstractVector{RevCholWishart{T}},
                        cholinvwisharts::AbstractVector{CholInvWishart{T}},
                        invwisharts::AbstractVector{InverseWishart{T}},
                        groups::Groups{NG}, rank1covs::Vector{InverseWishart{T}}) where {T,NG}
    set_priors!(invwisharts, Val(NG))

    for (i,g) ∈ enumerate(groups)
        @inbounds invwisharts[g] += rank1covs[i]
    end

    inv_and_cholesky!(revcholwisharts, cholinvwisharts, invwisharts)
end

# @generated function update_individual_probs!(probabilities::AbstractMatrix{T},
#                                 baseπ::AbstractVector{T}, LiV::AbstractVector{RevCholWishart{T}},
#                                 ν::NTuple{NG,T}, x::AbstractMatrix{T}) where {T,NG}
#     quote
#         @inbounds for g ∈ 1:NG
#             # L is scaled too large by a factor of √ν; cancels out 1/ν factor on quadratic form in t-pdf
#             exponent = T(-0.5) * ν[g] + T(-1.5)
#             Li = LiV[g]
#             Li11             = Li[1]
#             Li21, Li22       = Li[2], Li[4]
#             Li31, Li32, Li33 = Li[3], Li[5], Li[6]
#
#             base = log(baseπ[g]) + log(Li11) + log(Li22) + log(Li33) +
#                     lgamma(-exponent) - lgamma(T(0.5)*ν[g]) - T(1.5)*log(ν[g])
#
#             @vectorize $T for i ∈ 1:size(x, 1)
#                 lx₁ = Li11*x[i,1]
#                 lx₂ = Li21*x[i,1] + Li22*x[i,2]
#                 lx₃ = Li31*x[i,1] + Li32*x[i,2] + Li33*x[i,3]
#                 probabilities[i,g] = exp(base + exponent * log(one(T) + lx₁*lx₁ + lx₂*lx₂ + lx₃*lx₃))
#             end
#         end
#     end
# end
# @generated function update_individual_probs!(probabilities::AbstractMatrix{T},
#                                 baseπ::AbstractVector{T}, LiV::AbstractVector{RevCholWishart{T}},
#                                 ν::NTuple{NG,T}, x::AbstractMatrix{T}) where {T,NG}
#     quote
#         @inbounds for g ∈ 1:NG
#             # L is scaled too large by a factor of √ν; cancels out 1/ν factor on quadratic form in t-pdf
#             exponent = T(-0.5) * ν[g] + T(-1.5)
#             Li = LiV[g]
#             Li11             = Li[1]
#             Li21, Li22       = Li[2], Li[4]
#             Li31, Li32, Li33 = Li[3], Li[5], Li[6]
#
#             base = log(baseπ[g]) + log(Li11) + log(Li22) + log(Li33) +
#                     lgamma(-exponent) - lgamma(T(0.5)*ν[g]) - T(1.5)*log(ν[g])
#
#             @vectorize $T for i ∈ 1:size(x, 1)
#                 lx₁ = Li11*x[i,1]
#                 lx₂ = Li21*x[i,1] + Li22*x[i,2]
#                 lx₃ = Li31*x[i,1] + Li32*x[i,2] + Li33*x[i,3]
#                 probabilities[i,g] = log(one(T) + lx₁*lx₁ + lx₂*lx₂ + lx₃*lx₃)
#             end
#             @vectorize $T for i ∈ 1:size(x, 1)
#                 probabilities[i,g] = exp(base + exponent * probabilities[i,g])
#             end
#         end
#     end
# end

@generated function update_individual_probs!(probabilities::AbstractMatrix{T},
                                baseπ::AbstractVector{T}, LiV::AbstractVector{RevCholWishart{T}},
                                ν::NTuple{NG,T}, x::AbstractMatrix{T}) where {T,NG}

    # unroll is 3 for AVX512f, 1 otherwise
    unroll = max(1, ( VectorizationBase.REGISTER_COUNT - 4) ÷ 9)
    if unroll > 1
        iters = cld(NG, unroll)
         # more evenly distributes iters. Eg, if unroll was 3 and NG 4, now unroll = 2.
        unroll = cld(NG, iters)
    end
    rem = NG % unroll

    q = quote
        base = MVector{$NG, $T}(undef)
        exponent = MVector{$NG, $T}(undef)
        @inbounds for rep ∈ 0:$unroll:$(NG-unroll)
                # L is scaled too large by a factor of √ν; cancels out 1/ν factor on quadratic form in t-pdf
            Base.Cartesian.@nexprs $unroll k -> begin
                exponent_k = T(-0.5) * ν[rep+k] + T(-1.5)
                Li_k = LiV[rep+k]
                Li11_k                 = Li_k[1]
                Li21_k, Li22_k         = Li_k[2], Li_k[4]
                Li31_k, Li32_k, Li33_k = Li_k[3], Li_k[5], Li_k[6]

                exponent[rep+k] = exponent_k
                base[rep+k] = log(baseπ[rep+k]) + log(Li11_k) + log(Li22_k) + log(Li33_k) +
                        lgamma(-exponent_k) - lgamma(T(0.5)*ν[rep+k]) - T(1.5)*log(ν[rep+k])
            end
            @vectorize $T for i ∈ 1:size(x, 1)
                Base.Cartesian.@nexprs $unroll k -> begin
                    lx₁_k = Li11_k*x[i,1]
                    lx₂_k = Li21_k*x[i,1] + Li22_k*x[i,2]
                    lx₃_k = Li31_k*x[i,1] + Li32_k*x[i,2] + Li33_k*x[i,3]
                    probabilities[i,rep+k] = one(T) + lx₁_k*lx₁_k + lx₂_k*lx₂_k + lx₃_k*lx₃_k
                end
            end
        end
    end
    if rem > 0
        NGmrep = NG - rem
        push!(q.args, quote
                    Base.Cartesian.@nexprs $rem k -> begin
                        exponent_k = T(-0.5) * ν[$NGmrep+k] + T(-1.5)
                        Li_k = LiV[$NGmrep+k]
                        Li11_k                 = Li_k[1]
                        Li21_k, Li22_k         = Li_k[2], Li_k[4]
                        Li31_k, Li32_k, Li33_k = Li_k[3], Li_k[5], Li_k[6]

                        exponent[$NGmrep+k] = exponent_k
                        base[$NGmrep+k] = log(baseπ[$NGmrep+k]) + log(Li11_k) + log(Li22_k) + log(Li33_k) +
                                lgamma(-exponent_k) - lgamma(T(0.5)*ν[$NGmrep+k]) - T(1.5)*log(ν[$NGmrep+k])
                    end
                    @vectorize $T for i ∈ 1:size(x, 1)
                        Base.Cartesian.@nexprs $rem k -> begin
                            lx₁_k = Li11_k*x[i,1]
                            lx₂_k = Li21_k*x[i,1] + Li22_k*x[i,2]
                            lx₃_k = Li31_k*x[i,1] + Li32_k*x[i,2] + Li33_k*x[i,3]
                            probabilities[i,$NGmrep+k] = one(T) + lx₁_k*lx₁_k + lx₂_k*lx₂_k + lx₃_k*lx₃_k
                        end
                    end
                end)

    end
    push!(q.args, quote
        @vectorize $T for i ∈ 1:length(probabilities)
            probabilities[i] = log(probabilities[i])
        end
        # @inbounds for g ∈ 1:$NG
        #     @vectorize $T for i ∈ 1:size(x, 1)
        #         probabilities[i,g] = log(probabilities[i,g])
        #     end
        # end
        @inbounds for g ∈ 1:$NG
            expg = exponent[g]
            baseg = base[g]
            @vectorize $T for i ∈ 1:size(x, 1)
                probabilities[i,g] = exp(baseg + expg * probabilities[i,g])
            end
        end
    end)
    q
end

function update_individual_probs!(probabilities::AbstractMatrix{T}, baseπ, groups::Groups{NG},
                            Li::AbstractMatrix{T}, ν, x::AbstractMatrix{T}) where {T,NG}
    update_individual_probs!(Random.GLOBAL_RNG, probabilities, baseπ, groups, Li, ν, x)
end

function run_sample!(rng::PCG_Scalar_and_Vector, mcmcres::MCMCResult, workingdata::WorkingData{NG},
                    X::AbstractMatrix{T}, rank1covs::AbstractVector{InverseWishart{T}},
                    baseπ::AbstractVector{T}, chain::Integer = 1, iter = 10000, warmup = 4000) where {T,NG}

    N = length(rank1covs)
    resize!(workingdata, N)

    invwisharts, individual_probs, groups =
        workingdata.inverse_wisharts, workingdata.individual_probs, workingdata.groups

    probs, revcholwisharts, cholinvwisharts =
                    mcmcres.Probs, mcmcres.RevCholWisharts, mcmcres.CholInvWisharts

    # chain jump
    CJ = (chain - 1) * iter + 1

    rand!(rng.vector, groups, baseπ)
    # returning mcmcres keeps revcholwisharts, cholinvwisharts
    # therefore no need for GC.@preserve
    @uviews probs begin
        @inbounds @views begin
            calc_Wisharts!(revcholwisharts[:,CJ], cholinvwisharts[:,CJ], invwisharts, groups, rank1covs)
            for i ∈ 1:warmup
                # @show extract_α(invwisharts, Val(NG))
                randdirichlet!(rng.scalar, probs[:,CJ], extract_α(invwisharts, Val(NG)))
                update_individual_probs!(individual_probs, probs[:,CJ],
                        revcholwisharts[:,CJ], extract_ν(invwisharts, Val(NG)), X)
                rand!(rng.vector, groups, individual_probs)
                calc_Wisharts!(revcholwisharts[:,CJ], cholinvwisharts[:,CJ], invwisharts, groups, rank1covs)
            end
            randdirichlet!(rng.scalar, probs[:,CJ], extract_α(invwisharts, Val(NG)))
            for i ∈ 1:iter-1
                update_individual_probs!(individual_probs, probs[:,CJ+i-1],
                        revcholwisharts[:,CJ+i-1], extract_ν(invwisharts, Val(NG)), X)
                rand!(rng.vector, groups, individual_probs)
                calc_Wisharts!(revcholwisharts[:,CJ+i], cholinvwisharts[:,CJ+i], invwisharts, groups, rank1covs)
                randdirichlet!(rng.scalar, probs[:,CJ+i], extract_α(invwisharts, Val(NG)))
            end
        end
    end
    mcmcres
end
function run_sample!(rng::AbstractRNG, mcmcres::MCMCResult, workingdata::WorkingDataUnifCache{NG},
                    X::AbstractMatrix{T}, rank1covs::AbstractVector{InverseWishart{T}},
                    baseπ::AbstractVector{T}, chain::Integer = 1, iter = 10000, warmup = 4000) where {T,NG}

    N = length(rank1covs)
    resize!(workingdata, N)

    invwisharts, individual_probs, uniform_probs, groups =
        workingdata.inverse_wisharts, workingdata.individual_probs,
        workingdata.uniform_probs, workingdata.groups

    probs, revcholwisharts, cholinvwisharts =
                    mcmcres.Probs, mcmcres.RevCholWisharts, mcmcres.CholInvWisharts

    # chain jump
    CJ = (chain - 1) * iter + 1

    rand!(rng, groups, uniform_probs, baseπ)
    # returning mcmcres keeps revcholwisharts, cholinvwisharts
    # therefore no need for GC.@preserve
    @uviews probs begin
        @inbounds @views begin
            calc_Wisharts!(revcholwisharts[:,CJ], cholinvwisharts[:,CJ], invwisharts, groups, rank1covs)
            for i ∈ 1:warmup
                randdirichlet!(rng, probs[:,CJ], extract_α(invwisharts, Val(NG)))
                update_individual_probs!(individual_probs, probs[:,CJ],
                        revcholwisharts[:,CJ], extract_ν(invwisharts, Val(NG)), X)
                rand!(rng, groups, uniform_probs, individual_probs)
                calc_Wisharts!(revcholwisharts[:,CJ], cholinvwisharts[:,CJ], invwisharts, groups, rank1covs)
            end
            randdirichlet!(rng, probs[:,CJ], extract_α(invwisharts, Val(NG)))
            for i ∈ 1:iter-1
                update_individual_probs!(individual_probs, probs[:,CJ+i-1],
                        revcholwisharts[:,CJ+i-1], extract_ν(invwisharts, Val(NG)), X)
                rand!(rng, groups, uniform_probs, individual_probs)
                calc_Wisharts!(revcholwisharts[:,CJ+i], cholinvwisharts[:,CJ+i], invwisharts, groups, rank1covs)
                randdirichlet!(rng, probs[:,CJ+i], extract_α(invwisharts, Val(NG)))
            end
        end
    end
    mcmcres
end

function singlechain_sample!(mcmcres, X, rank1covs, workingdata, BPP, baseπ, iter = size(mcmcres.Probs, 2), warmup = 4000)

    process_big_prop_points!(X, BPP)
    generate_rank1covariances!(rank1covs, X)

    run_sample!(GLOBAL_PCG, mcmcres, workingdata, X, rank1covs, baseπ, 1, iter, warmup)

end
function singlechain_sample!(rng::AbstractRNG, mcmcres, X, rank1covs, workingdata, BPP, baseπ, iter = size(mcmcres.Probs, 2), warmup = 4000)

    process_big_prop_points!(X, BPP)
    generate_rank1covariances!(rank1covs, X)

    run_sample!(rng, mcmcres, workingdatachains, X, rank1covs, baseπ, 1, iter, warmup)

end

function thread_sample!(mcmcres, X, rank1covs, workingdatachains, BPP, baseπ, iter,
                                        warmup = 4000, chains = length(workingdatachains))

    process_big_prop_points!(X, BPP)
    generate_rank1covariances!(rank1covs, X)

    @threads for chain ∈ 1:chains
        run_sample!(TRNG[chain], mcmcres, workingdatachains[chain], X, rank1covs, baseπ, chain, iter, warmup)
    end

end


text:
export NeverObservable, never

import Base: show

"""
    NeverObservable{D}()

An Observable that emits no items to the Observer and never completes.

# Type parameters
- `D`: Type of Observable data

See also: [`Subscribable`](@ref), [`never`](@ref)
"""
@subscribable struct NeverObservable{D} <: Subscribable{D} end

function on_subscribe!(observable::NeverObservable, actor)
    return voidTeardown
end

"""
    never(T = Any)

Creation operator for the `NeverObservable` that emits neither values nor errors nor the completion notification.
It can be used for testing purposes or for composing with other Observables.
Please note that by never emitting a complete notification, this Observable keeps
the subscription from being disposed automatically. Subscriptions need to be manually
disposed.

# Arguments
- `T`: Type of Observable data, optional, `Any` is the default

# Examples

```jldoctest
using Rocket

source = never()
subscribe!(source, logger())
;

# output

```

See also: [`NeverObservable`](@ref), [`subscribe!`](@ref), [`logger`](@ref)
"""
never(T = Any) = NeverObservable{T}()

Base.show(io::IO, ::NeverObservable{D}) where D = print(io, "NeverObservable($D)")


text:
function optimize(f::Function, g::Function, h::Function, x::AbstractVector{T};
        alg::MultivariateAlgorithm,
        abstol=1e-12,
        maxiter::Integer=1_000_000,
        savetrace::Bool=false,
        callback=()->nothing) where {T<:AbstractFloat}
    reset!(alg)
    x = copy(x)
    metadata = Dict{Symbol, Any}()
    curr_state = state(order(alg), f, g, h, x)
    trace = typeof(curr_state)[]
    converged = false
    iter = 0
    while iter <= maxiter
        callback()
        if savetrace
            push!(trace, curr_state)
            x = copy(x)
        end
        if norm(curr_state.∇f) < abstol
            converged = true
            break
        end
        (iter == maxiter || isinf(curr_state.f)) && break
        iter += 1
        step!(x, alg, f, curr_state)
        curr_state = state(order(alg), f, g, h, x)
    end
    savetrace && (metadata[:trace] = trace)
    return Solution(converged, iter, curr_state.x, curr_state.f, metadata)
end

function optimize(f::Function, g::Function, x; kwargs...)
    h = x -> Zygote.hessian(f, x)
    return optimize(f, g, h, x; kwargs...)
end

function optimize(f::Function, x; kwargs...)
    g = x -> Zygote.gradient(f, x)[1]
    h = x -> Zygote.hessian(f, x)
    return optimize(f, g, h, x; kwargs...)
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 DisallowedConfiguration <: SwaggerModel
    vmDiskType::Any # spec type: Union{ Nothing, String } # spec name: vmDiskType

    function DisallowedConfiguration(;vmDiskType=nothing)
        o = new()
        validate_property(DisallowedConfiguration, Symbol("vmDiskType"), vmDiskType)
        setfield!(o, Symbol("vmDiskType"), vmDiskType)
        o
    end
end # type DisallowedConfiguration

const _property_map_DisallowedConfiguration = Dict{Symbol,Symbol}(Symbol("vmDiskType")=>Symbol("vmDiskType"))
const _property_types_DisallowedConfiguration = Dict{Symbol,String}(Symbol("vmDiskType")=>"String")
Base.propertynames(::Type{ DisallowedConfiguration }) = collect(keys(_property_map_DisallowedConfiguration))
Swagger.property_type(::Type{ DisallowedConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DisallowedConfiguration[name]))}
Swagger.field_name(::Type{ DisallowedConfiguration }, property_name::Symbol) =  _property_map_DisallowedConfiguration[property_name]

const _allowed_DisallowedConfiguration_vmDiskType = ["None", "Unmanaged"]

function check_required(o::DisallowedConfiguration)
    true
end

function validate_property(::Type{ DisallowedConfiguration }, name::Symbol, val)
    if name === Symbol("vmDiskType")
        Swagger.validate_param(name, "DisallowedConfiguration", :enum, val, _allowed_DisallowedConfiguration_vmDiskType)
    end
end


text:
using Documenter, Sounds, Unitful
makedocs(
  modules = [Sounds],
  format = :html,
  sitename = "Sounds.jl",
  html_prettyurls = true,
  pages = [
    "Introduction" => "index.md",
    "Manual" => "manual.md",
    "Reference" => "reference.md"
  ]
)
deploydocs(
  repo = "github.com/haberdashPI/Sounds.jl.git",
  julia = "0.6",
  osname = "osx",
  deps = nothing,
  make = nothing,
  target = "build"
)


text:
#Pkg.add("JuliaBerry")
#Pkg.add("PyCall")
#Run with `sudo` -- needed for keyboard.

using JuliaBerry
using JuliaBerry.ExplorerHat
using PyCall

@pyimport keyboard

function movef(x)
    forward(ExplorerHat.motor[1], 50)
    forward(ExplorerHat.motor[2], 50)
    sleep(x)
    stop(ExplorerHat.motor[1])
    stop(ExplorerHat.motor[2])
end

function moveb(x)
    on(ExplorerHat.led)
    backward(ExplorerHat.motor[1], 50)
    backward(ExplorerHat.motor[2], 50)
    sleep(x)
    stop(ExplorerHat.motor[1])
    stop(ExplorerHat.motor[2])
    on(ExplorerHat.led)
end

function movel(x)
    forward(ExplorerHat.motor[2], 75)
    forward(ExplorerHat.motor[1], 5)
    sleep(x)
    stop(ExplorerHat.motor[1])
    stop(ExplorerHat.motor[2])
end

function mover(x)
    forward(ExplorerHat.motor[2], 5)
    forward(ExplorerHat.motor[1], 75)
    sleep(x)
    stop(ExplorerHat.motor[1])
    stop(ExplorerHat.motor[2])
end

function movebl(x)
    explorerhat.light.on()
    backward(ExplorerHat.motor[2], 5)
    backward(ExplorerHat.motor[1], 75)
    sleep(x)
    stop(ExplorerHat.motor[1])
    stop(ExplorerHat.motor[2])
    explorerhat.light.off()
end

function movebr(x)
    explorerhat.light.on()
    backward(ExplorerHat.motor[2], 75)
    backward(ExplorerHat.motor[1], 5)
    sleep(x)
    stop(ExplorerHat.motor[1])
    stop(ExplorerHat.motor[2])
    explorerhat.light.off()

end

function square()
    for i in 1:4
        movef(3)
        mover(2)
    end
end

function eight()
    for i in 1:3
        movef(3)
        mover(2)
    end
    movef(6)
    for i in 1:3
        movel(2)
        movef(3)
    end
end

function disco(x)
    y = x/8
    forward(ExplorerHat.motor[1], 50)
    backward(ExplorerHat.motor[2], 50)
    for i in 1:8
        on(ExplorerHat.led)
        sleep(y)
        off(ExplorerHat.led)
        sleep(y)
    end
    stop(ExplorerHat.motor[1])
    stop(ExplorerHat.motor[2])
end

println("Starting main loop")
while true
    if keyboard.is_pressed("w")
        movef(0.01)
    elseif keyboard.is_pressed("a")
        movel(0.01)
    elseif keyboard.is_pressed("d")
        mover(0.01)
    elseif keyboard.is_pressed("s")
        moveb(0.01)
    elseif keyboard.is_pressed("z")
        movebl(0.01)
    elseif keyboard.is_pressed("c")
        movebr(0.01)
    elseif keyboard.is_pressed("4")
        square()
    elseif keyboard.is_pressed("8")
        eight()
    elseif keyboard.is_pressed("l")
        disco(7)
    elseif keyboard.is_pressed(" ")
        if is_on(ExplorerHat.output[1])
            off(ExplorerHat.output[1])
        else
            on(ExplorerHat.output[1])
        end
    end
end


text:
export PrioritizedDQNLearner

"""
    PrioritizedDQNLearner(;kwargs...)

See paper: [Prioritized Experience Replay](https://arxiv.org/abs/1511.05952)
And also https://danieltakeshi.github.io/2019/07/14/per/

# Keywords

- `approximator`::[`AbstractApproximator`](@ref): used to get Q-values of a state.
- `target_approximator`::[`AbstractApproximator`](@ref): similar to `approximator`, but used to estimate the target (the next state).
- `loss_func`: the loss function.
- `γ::Float32=0.99f0`: discount rate.
- `batch_size::Int=32`
- `update_horizon::Int=1`: length of update ('n' in n-step update).
- `min_replay_history::Int=32`: number of transitions that should be experienced before updating the `approximator`.
- `update_freq::Int=4`: the frequency of updating the `approximator`.
- `target_update_freq::Int=100`: the frequency of syncing `target_approximator`.
- `stack_size::Union{Int, Nothing}=4`: use the recent `stack_size` frames to form a stacked state.
- `default_priority::Float64=100.`: the default priority for newly added transitions.
- `rng = Random.GLOBAL_RNG`
"""
mutable struct PrioritizedDQNLearner{
    Tq<:AbstractApproximator,
    Tt<:AbstractApproximator,
    Tf,
    R<:AbstractRNG,
} <: AbstractLearner
    approximator::Tq
    target_approximator::Tt
    loss_func::Tf
    sampler::NStepBatchSampler
    min_replay_history::Int
    update_freq::Int
    target_update_freq::Int
    update_step::Int
    default_priority::Float32
    β_priority::Float32
    rng::R
    # for logging
    loss::Float32
end

function PrioritizedDQNLearner(;
    approximator::Tq,
    target_approximator::Tt,
    loss_func::Tf,
    stack_size::Union{Int,Nothing} = nothing,
    γ::Float32 = 0.99f0,
    batch_size::Int = 32,
    update_horizon::Int = 1,
    min_replay_history::Int = 32,
    update_freq::Int = 1,
    target_update_freq::Int = 100,
    update_step::Int = 0,
    default_priority::Float32 = 100.0f0,
    β_priority::Float32 = 0.5f0,
    traces = SARTS,
    rng = Random.GLOBAL_RNG,
) where {Tq,Tt,Tf}
    copyto!(approximator, target_approximator)
    sampler = NStepBatchSampler{traces}(;γ=γ, n=update_horizon,stack_size=stack_size,batch_size=batch_size)
    PrioritizedDQNLearner(
        approximator,
        target_approximator,
        loss_func,
        sampler,
        min_replay_history,
        update_freq,
        target_update_freq,
        update_step,
        default_priority,
        β_priority,
        rng,
        0.0f0,
    )
end

Flux.functor(x::PrioritizedDQNLearner) = (Q = x.approximator, Qₜ = x.target_approximator),
y -> begin
    x = @set x.approximator = y.Q
    x = @set x.target_approximator = y.Qₜ
    x
end

"""

!!! note
    The state of the observation is assumed to have been stacked,
    if `!isnothing(stack_size)`.
"""
function (learner::PrioritizedDQNLearner)(env)
    env |>
    get_state |>
    x -> Flux.unsqueeze(x, ndims(x) + 1) |>
    x -> send_to_device(device(learner), x) |>
    learner.approximator |>
    vec |>
    send_to_host
end

function RLBase.update!(learner::PrioritizedDQNLearner, batch::NamedTuple)
    Q = learner.approximator
    Qₜ = learner.target_approximator
    γ = learner.sampler.γ
    β = learner.β_priority
    loss_func = learner.loss_func
    n = learner.sampler.n
    batch_size = learner.sampler.batch_size

    D = device(Q)
    s, a, r, t, s′ = (send_to_device(D,batch[x]) for x in SARTS)
    a = CartesianIndex.(a, 1:batch_size)

    updated_priorities = Vector{Float32}(undef, batch_size)
    w = 1.0f0 ./ ((batch.priority .+ 1f-10) .^ β)
    w ./= maximum(w)
    w = send_to_device(D, w)

    target_q = Qₜ(s′)
    if haskey(batch, :next_legal_actions_mask)
        l′ = send_to_device(D, batch[:next_legal_actions_mask])
        target_q .+= ifelse.(l′, 0.f0, typemin(Float32))
    end

    q′ = dropdims(maximum(target_q; dims = 1), dims = 1)
    G = r .+ γ^n .* (1 .- t) .* q′

    gs = gradient(params(Q)) do
        q = Q(s)[a]
        batch_losses = loss_func(G, q)
        loss = dot(vec(w), vec(batch_losses)) * 1 // batch_size
        ignore() do
            updated_priorities .= send_to_host(vec((batch_losses .+ 1f-10) .^ β))
            learner.loss = loss
        end
        loss
    end

    update!(Q, gs)
    updated_priorities
end


text:
###############################################################################
#
# Capacitance for plants
#
###############################################################################
"""
    $(TYPEDEF)

Hierachy of AbstractCapacity
- [`PVCurveLinear`](@ref)
- [`PVCurveSegmented`](@ref)
"""
abstract type AbstractCapacity{FT<:AbstractFloat} end

"""
    $(TYPEDEF)

Struct that contains information for linear PV curve

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct PVCurveLinear{FT<:AbstractFloat} <: AbstractCapacity{FT}
    "Slope of the linear PV curve (relative to maximum) `[MPa⁻¹]`"
    slope   ::FT = 0.2
    "Conductance for refilling (relative to maximum) `[MPa⁻¹ s⁻¹]`"
    k_refill::FT = 0.0001
end

"""
    $(TYPEDEF)

"""
Base.@kwdef mutable struct PVCurveSegmented{FT} <: AbstractCapacity{FT}
    "n_o / maximum V `[mol m⁻³]`"
    c_all   ::FT = 300.0
    "Bulk modulus of elasticity `[MPa]`"
    ϵ_bulk  ::FT = 20.0
    "Apoplastic water content relative to maximum water volume"
    RWC_apo ::FT = 0.2
    "Relative water content at turgor loss point"
    RWC_TLP ::FT = 0.8
    "Conductance for refilling (relative to maximum) `[MPa⁻¹ s⁻¹]`"
    k_refill::FT = 0.0001
end

###############################################################################
#
# Vulnerability curves for xylem
#
###############################################################################
"""
    $(TYPEDEF)

Hierachy of AbstractXylemVC
- [`WeibullSingle`](@ref)
- [`WeibullDual`](@ref)
"""
abstract type AbstractXylemVC{FT<:AbstractFloat} end

"""
    $(TYPEDEF)

Struct that contains dual Weibull function parameters.

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct WeibullDual{FT} <: AbstractXylemVC{FT}
    "B of first part `[MPa]`"
    b1::FT = 1
    "C of first part"
    c1::FT = 5
    "F of first part"
    f1::FT = 0.5
    "B of second part `[MPa]`"
    b2::FT = 2
    "C of second part"
    c2::FT = 5
    "F of second part"
    f2::FT = 1 - f1
end

"""
    $(TYPEDEF)

Struct that contains single Weibull function parameters.

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct WeibullSingle{FT} <: AbstractXylemVC{FT}
    "B of first part `[MPa]`"
    b::FT = 2
    "C of first part"
    c::FT = 5
end

###############################################################################
#
# Vulnerability curves for soil
#
###############################################################################
"""
    abstract type AbstractSoilVC{FT}

Hierachy of AbstractSoilVC:
- [`BrooksCorey`](@ref)
- [`VanGenuchten`](@ref)
"""
abstract type AbstractSoilVC{FT<:AbstractFloat} end

"""
    mutable struct BrooksCorey{FT<:AbstractFloat}

Brooks Corey soil parameters

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct BrooksCorey{FT} <: AbstractSoilVC{FT}
    "Soil type"
    stype::String = "TBD"
    "Soil b"
    b ::FT = 8.52
    "ϕ at saturation `[MPa]`"
    ϕs::FT = 0.63 * ρg_MPa(FT)
    "Soil water content (Θ) at saturation"
    Θs::FT = 0.476
    "Residual soil water content"
    Θr::FT = 0.1
end

"""
    mutable struct VanGenuchten{FT<:AbstractFloat}

Van Gunechten soil parameters

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct VanGenuchten{FT} <: AbstractSoilVC{FT}
    "Soil type"
    stype::String = "Sandy Clay Loam"
    "Soil α is related to the inverse of the air entry suction, α > 0"
    α ::FT = 602.0419
    "Soil n is Measure of the pore-size distribution"
    n ::FT = 1.48
    "Soil m = 1 - 1/n"
    m ::FT = 1 - 1/n
    "Soil water content (Θ) at saturation"
    Θs::FT = 0.39
    "Residual soil water content"
    Θr::FT = 0.1
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land