text:
"""
    JointConstraint{T} <: Constraint{T}

    constraint restricting translational and rotational degrees of freedom between two Body objects.

    id: a unique identifying number
    name: a unique identifying name
    translational: Translational
    rotational: Rotational
    spring: flag for joint springs on
    damper: flag for joint dampers on
    parent_id: identifying number for parent Body{T}
    child_id: identifying number for child Body{T}
    minimal_index: indices for minimal coordinates
    impulses: joint impulses that maintain constraint between two Body{T} objects
"""
mutable struct JointConstraint{T,N,Nc,TJ,RJ} <: Constraint{T,N}
    # ID
    id::Int64
    name::Symbol

    # joint constraints
    translational::TJ
    rotational::RJ

    # springs and dampers
    spring::Bool
    damper::Bool

    # neighbor IDs
    parent_id::Int
    child_id::Int

    # indices
    minimal_index::SVector{Nc,SVector{2,Int64}} # indices for minimal coordinates, assumes joints # Nc = 2 THIS IS SPECIAL CASED

    # impulses
    impulses::Vector{SVector{N,T}}

    function JointConstraint(data;
        name::Symbol=Symbol("joint_" * randstring(4)))

        @assert data[1][2] == data[2][2] # check parent ids
        @assert data[1][3] == data[2][3] # check child ids

        # joints
        translational = data[1][1]
        rotational = data[2][1]

        # IDs
        parent_id = data[1][2]
        child_id = data[1][3]

        # data dype
        T = typeof(data[1][1]).parameters[1]

        # set springs & dampers off
        spring = false
        damper = false

        minimal_index = Vector{Int64}[]
        N = 0
        for joint_data in data
            joint = joint_data[1]

            # set spring & damper on
            joint.spring != 0 && (spring = true)
            joint.damper != 0 && (damper = true)

            # minimal-coordaintes indices
            Nλ = joint_length(joint)
            Nset = impulses_length(joint)
            if isempty(minimal_index)
                push!(minimal_index, [1;3-Nλ])
            else
                push!(minimal_index, [last(minimal_index)[2]+1; last(minimal_index)[2]+3-Nλ])
            end
            N += Nset
        end

        Nc = 2
        impulses = [zeros(T, N) for i=1:2]

        return new{T,N,Nc,typeof(translational),typeof(rotational)}(getGlobalID(), name, translational, rotational, spring, damper, parent_id, child_id, minimal_index, impulses)
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, constraint::JointConstraint)
    summary(io, constraint)
    println(io, "")
    println(io, "id:            "*string(constraint.id))
    println(io, "name:          "*string(constraint.name))
    println(io, "spring:        "*string(constraint.spring))
    println(io, "damper:        "*string(constraint.damper))
    println(io, "damper:        "*string(constraint.damper))
    println(io, "parent_id:     "*string(constraint.parent_id))
    println(io, "child_id:      "*string(constraint.child_id))
    println(io, "minimal_index: "*string(constraint.minimal_index))
    println(io, "impulses:      "*string(constraint.impulses))
end

# constraints
@generated function constraint(mechanism, joint::JointConstraint)
    pbody = :(get_body(mechanism, joint.parent_id))
    cbody = :(get_body(mechanism, joint.child_id))
    tra = :(constraint(joint.translational,
        $pbody, $cbody,
        joint.impulses[2][joint_impulse_index(joint,1)], mechanism.μ, mechanism.timestep))
    rot = :(constraint(joint.rotational,
        $pbody, $cbody,
        joint.impulses[2][joint_impulse_index(joint,2)], mechanism.μ, mechanism.timestep))
    return :(svcat($tra, $rot))
end

# # constraints Jacobians
# @generated function constraint_jacobian(joint::JointConstraint)
#     tra = :(constraint_jacobian(joint.translational, joint.impulses[2][joint_impulse_index(joint, 1)]))
#     rot = :(constraint_jacobian(joint.rotational, joint.impulses[2][joint_impulse_index(joint, 2)]))
#     return :(cat($tra, $rot, dims=(1,2)))
# end
function constraint_jacobian(joint::JointConstraint)
    tra = constraint_jacobian(joint.translational, joint.impulses[2][joint_impulse_index(joint, 1)])
    rot = constraint_jacobian(joint.rotational, joint.impulses[2][joint_impulse_index(joint, 2)])
    return diagonal_cat(tra, rot)
end

@generated function constraint_jacobian_configuration(mechanism, joint::JointConstraint, body::Body)
    relative = :(body.id == joint.parent_id ? :parent : :child)
    pbody = :(get_body(mechanism, joint.parent_id))
    cbody = :(get_body(mechanism, joint.child_id))
    tra = :(constraint_jacobian_configuration($relative,
        joint.translational,
        $pbody, $cbody,
        joint.impulses[2][joint_impulse_index(joint, 1)], mechanism.timestep))
    rot = :(constraint_jacobian_configuration($relative,
        joint.rotational,
        $pbody, $cbody,
        joint.impulses[2][joint_impulse_index(joint, 2)], mechanism.timestep))
    return :(vcat($tra, $rot))
end

# impulses
function impulses!(mechanism, body::Body, joint::JointConstraint{T,Nλ}) where {T,Nλ}
    (Nλ > 0) && (body.state.d -= impulse_map(mechanism, joint, body) * joint.impulses[2])
    joint.spring && (body.state.d -= spring_impulses(mechanism, joint, body))
    joint.damper && (body.state.d -= damper_impulses(mechanism, joint, body))
    return
end

function impulse_map(mechanism, joint::JointConstraint, body::Body)
    relative = body.id == joint.parent_id ? :parent : :child
    pbody = get_body(mechanism, joint.parent_id)
    cbody = get_body(mechanism, joint.child_id)
    tra = impulse_map(relative, joint.translational,
        pbody, cbody,
        joint.impulses[2][joint_impulse_index(joint, 1)])
    rot = impulse_map(relative, joint.rotational,
        pbody, cbody,
        joint.impulses[2][joint_impulse_index(joint, 2)])
    return hcat(tra, rot)
end

# @generated function impulse_map(mechanism, joint::JointConstraint, body::Body)
#     relative = :(body.id == joint.parent_id ? :parent : :child)
#     pbody = :(get_body(mechanism, joint.parent_id))
#     cbody = :(get_body(mechanism, joint.child_id))
#     tra = :(impulse_map($relative, joint.translational,
#         $pbody, $cbody,
#         joint.impulses[2][joint_impulse_index(joint, 1)]))
#     rot = :(impulse_map($relative, joint.rotational,
#         $pbody, $cbody,
#         joint.impulses[2][joint_impulse_index(joint, 2)]))
#     return :(hcat($tra, $rot))
# end

# impulses Jacobians
function impulses_jacobian_velocity!(mechanism, body::Body, joint::JointConstraint)

    # relative
    relative = (body.id == joint.parent_id ? :parent : (body.id == joint.child_id ? :child : error()))

    # time step
    timestep= mechanism.timestep

    # bodies
    pbody = get_body(mechanism, joint.parent_id)
    cbody = get_body(mechanism, joint.child_id)

    # springs
    joint.spring && (body.state.D -= spring_jacobian_velocity(relative, relative, joint.translational, pbody, cbody, timestep))
    joint.spring && (body.state.D -= spring_jacobian_velocity(relative, relative, joint.rotational, pbody, cbody, timestep))

    # dampers
    joint.damper && (body.state.D -= damper_jacobian_velocity(relative, relative, joint.translational, pbody, cbody, timestep))
    joint.damper && (body.state.D -= damper_jacobian_velocity(relative, relative, joint.rotational, pbody, cbody, timestep))

    return
end

# off-diagonal Jacobians
function off_diagonal_jacobians(mechanism, body::Body, joint::JointConstraint)
    return -impulse_map(mechanism, joint, body), constraint_jacobian_configuration(mechanism, joint, body) * integrator_jacobian_velocity(body, mechanism.timestep)
end

function off_diagonal_jacobians(mechanism, joint::JointConstraint, body::Body)
    return constraint_jacobian_configuration(mechanism, joint, body) * integrator_jacobian_velocity(body, mechanism.timestep), -impulse_map(mechanism, joint, body)
end

function off_diagonal_jacobians(mechanism, pbody::Body, cbody::Body)
    # time step
    timestep = mechanism.timestep

    # dimensions
    Ne = length(mechanism.joints)
    Nb = length(mechanism.bodies)
    Nc = length(mechanism.contacts)

    # Jacobian
    jacobian_parent_child = szeros(6, 6)
    jacobian_child_parent = szeros(6, 6)

    for connectionid in connections(mechanism.system, pbody.id)
        # joints
        if connectionid <= Ne
            joint = get_node(mechanism, connectionid)
            if pbody.id == joint.parent_id
                for element in (joint.translational, joint.rotational)
                    if cbody.id == joint.child_id
                        joint.damper && (jacobian_parent_child -= damper_jacobian_velocity(:parent, :child, element, pbody, cbody, timestep))
                        joint.damper && (jacobian_child_parent -= damper_jacobian_velocity(:child, :parent, element, pbody, cbody, timestep))
                    end
                end
            elseif cbody.id == joint.parent_id
                for element in (joint.translational, joint.rotational)
                    if pbody.id == joint.child_id
                        joint.damper && (jacobian_parent_child -= damper_jacobian_velocity(:parent, :child, element, cbody, pbody, timestep))
                        joint.damper && (jacobian_child_parent -= damper_jacobian_velocity(:child, :parent, element, cbody, pbody, timestep))
                    end
                end
            end
        end

        # contacts
        if connectionid > Ne + Nb
            contact = get_node(mechanism, connectionid)
            if pbody.id == contact.parent_id
                if cbody.id == contact.child_id
                    Jpc = impulse_map_jacobian(:parent, :child, contact.model,
                            pbody,
                            cbody,
                            contact.impulses[2],
                            mechanism.timestep) * integrator_jacobian_velocity(cbody, timestep)

                    Jcp = impulse_map_jacobian(:child, :parent, contact.model,
                            pbody,
                            cbody,
                            contact.impulses[2],
                            mechanism.timestep) * integrator_jacobian_velocity(pbody, timestep)
                    # impulse_map_jacobian_configuration(mechanism, body, contact) * integrator_jacobian_velocity(body, timestep)
                    # impulse_map(mechanism, contact, body) * contact.impulses[2]
                    jacobian_parent_child -= Jpc#damper_jacobian_velocity(:parent, :child, element, pbody, cbody, timestep)
                    jacobian_child_parent -= Jcp#damper_jacobian_velocity(:child, :parent, element, pbody, cbody, timestep)
                end
            elseif cbody.id == contact.parent_id
                if pbody.id == contact.child_id
                    Jpc = impulse_map_jacobian(:parent, :child, contact.model,
                            cbody,
                            pbody,
                            contact.impulses[2],
                            mechanism.timestep) * integrator_jacobian_velocity(pbody, timestep)

                    Jcp = impulse_map_jacobian(:child, :parent, contact.model,
                            cbody,
                            bbody,
                            contact.impulses[2],
                            mechanism.timestep) * integrator_jacobian_velocity(cbody, timestep)

                    jacobian_parent_child -= Jpc #damper_jacobian_velocity(:parent, :child, element, cbody, pbody, timestep)
                    jacobian_child_parent -= Jcp #damper_jacobian_velocity(:child, :parent, element, cbody, pbody, timestep)
                end
            end
        end
    end

    return jacobian_parent_child, jacobian_child_parent
end

# linear system
function set_matrix_vector_entries!(mechanism, matrix_entry::Entry, vector_entry::Entry, joint::JointConstraint)
    matrix_entry.value = constraint_jacobian(joint)
    vector_entry.value = -constraint(mechanism, joint)
end

# springs
function spring_impulses(mechanism, joint::JointConstraint{T}, body::Body{T};
    unitary::Bool=false) where T

    relative = (body.id == joint.parent_id ? :parent : :child)
    impulses = szeros(T,6)

    pbody = get_body(mechanism, joint.parent_id)
    cbody = get_body(mechanism, joint.child_id)

    impulses += spring_impulses(relative, joint.translational,
        pbody,
        cbody,
        mechanism.timestep,
        unitary=unitary)

    impulses += spring_impulses(relative, joint.rotational,
        pbody,
        cbody,
        mechanism.timestep,
        unitary=unitary)

    return impulses
end

# dampers
function damper_impulses(mechanism, joint::JointConstraint{T}, body::Body;
    unitary::Bool=false) where T

    relative = (body.id == joint.parent_id ? :parent : :child)
    impulses = szeros(T,6)

    pbody = get_body(mechanism, joint.parent_id)
    cbody = get_body(mechanism, joint.child_id)

    impulses += damper_impulses(relative, joint.translational,
        pbody,
        cbody,
        mechanism.timestep,
        unitary=unitary)

    impulses += damper_impulses(relative, joint.rotational,
        pbody,
        cbody,
        mechanism.timestep,
        unitary=unitary)

    return impulses
end

# inputs
function set_input!(joint::JointConstraint{T,N,Nc}, input::AbstractVector) where {T,N,Nc}
    @assert length(input) == input_dimension(joint)
    # translational
    r_idx = SUnitRange(joint.minimal_index[1][1], joint.minimal_index[1][2])
    length(r_idx) > 0 && set_input!(joint.translational, input[r_idx])
    # rotational
    r_idx = SUnitRange(joint.minimal_index[2][1], joint.minimal_index[2][2])
    length(r_idx) > 0 && set_input!(joint.rotational, input[r_idx])
    return
end

function add_input!(joint::JointConstraint{T,N,Nc}, input::AbstractVector) where {T,N,Nc}
    @assert length(input) == input_dimension(joint)
    add_input!(joint.translational, input[SUnitRange(joint.minimal_index[1][1], joint.minimal_index[1][2])])
    add_input!(joint.rotational, input[SUnitRange(joint.minimal_index[2][1], joint.minimal_index[2][2])])
    return
end

@generated function input_jacobian_control(mechanism, joint::JointConstraint{T,N,Nc}, body::Body) where {T,N,Nc}
    relative = :(body.id == joint.parent_id ? :parent : :child)
    pbody = :(get_body(mechanism, joint.parent_id))
    cbody = :(get_body(mechanism, joint.child_id))
    rot = :(input_jacobian_control($relative, joint.translational, $pbody, $cbody, mechanism.timestep))
    tra = :(input_jacobian_control($relative, joint.rotational, $pbody, $cbody, mechanism.timestep))
    return :(hcat($rot, $tra))
end

function input_impulse!(joint::JointConstraint{T,N,Nc}, mechanism, clear::Bool=true) where {T,N,Nc}
    pbody = get_body(mechanism, joint.parent_id)
    cbody = get_body(mechanism, joint.child_id)
    input_impulse!(joint.translational, pbody, cbody, mechanism.timestep, clear)
    input_impulse!(joint.rotational, pbody, cbody, mechanism.timestep, clear)
    return
end

# minimal
@generated function minimal_coordinates(mechanism, joint::JointConstraint{T,N,Nc}) where {T,N,Nc}
    pbody = :(get_body(mechanism, joint.parent_id))
    cbody = :(get_body(mechanism, joint.child_id))
    tra = :(minimal_coordinates(joint.translational, $pbody, $cbody))
    rot = :(minimal_coordinates(joint.rotational, $pbody, $cbody))
    return :(svcat($tra, $rot))
end

@generated function minimal_velocities(mechanism, joint::JointConstraint{T,N,Nc}) where {T,N,Nc}
    pbody = :(get_body(mechanism, joint.parent_id))
    cbody = :(get_body(mechanism, joint.child_id))
    tra = :(minimal_velocities(joint.translational, $pbody, $cbody, mechanism.timestep))
    rot = :(minimal_velocities(joint.rotational, $pbody, $cbody, mechanism.timestep))
    return :(svcat($tra, $rot))
end

################################################################################
# Utilities
################################################################################
function get_joint_impulses(joint::JointConstraint{T,N,Nc}, i::Int) where {T,N,Nc}
    n1 = 1
    for j = 1:i-1
        n1 += impulses_length((joint.translational, joint.rotational)[j])
    end
    n2 = n1 - 1 + impulses_length((joint.translational, joint.rotational)[i])

    λi = SVector{n2-n1+1,T}(joint.impulses[2][n1:n2])
    return λi
end

function joint_impulse_index(joint::JointConstraint{T,N,Nc}, i::Int) where {T,N,Nc}
    s = 0
    for j = 1:i-1
        element = (joint.translational, joint.rotational)[j]
        s += impulses_length(element)
    end
    joint_impulse_index((joint.translational, joint.rotational)[i], s)
end

# function reset!(joint::JointConstraint{T,N,Nc};
#     scale::T=1.0) where {T,N,Nc}
#     λ = []
#     for (i, element) in enumerate((joint.translational, joint.rotational))
#         Nλ = joint_length(element)
#         Nb = limits_length(element)
#         push!(λ, [scale * sones(2Nb); szeros(Nλ)])
#     end
#     joint.impulses[1] = vcat(λ...)
#     joint.impulses[2] = vcat(λ...)
#     return
# end

function reset!(joint::JointConstraint{T,N,Nc}; scale::T=1.0) where {T,N,Nc}
    Nλ_tra = joint_length(joint.translational)
    Nb_tra = limits_length(joint.translational)
    Nλ_rot = joint_length(joint.rotational)
    Nb_rot = limits_length(joint.rotational)
    joint.impulses[1] = [scale * sones(2Nb_tra); szeros(Nλ_tra); scale * sones(2Nb_rot); szeros(Nλ_rot)]
    joint.impulses[2] = [scale * sones(2Nb_tra); szeros(Nλ_tra); scale * sones(2Nb_rot); szeros(Nλ_rot)]
    return
end

function input_dimension(joint::JointConstraint{T,N,Nc};
    ignore_floating_base::Bool=false) where {T,N,Nc}
    ignore_floating_base && (N == 0) && return 0
    N̄ = 0
    N̄ = input_dimension(joint.translational) + input_dimension(joint.rotational)
    return N̄
end


text:
#!/usr/bin/env julia

using Match

function process_passport(line::String)
  kv_pairs = split(line, " ")

  found_keys = Set{String}()

  for field in kv_pairs
    key, value = split(field, ':')
    push!(found_keys, key)
  end

  all(["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"]) do k
    k ∈ found_keys
  end
end

function parse_num_range(num, lb, ub)
  if nothing === num
    return false
  end

  lb <= parse(Int, num) <= ub
end

function validate_passport(line::String)
  kv_pairs = split(line, " ")
  d = Dict{String, String}()

  for field in kv_pairs
    key, value = split(field, ':')
    d[key] = value
  end

  validator = Dict([
    ("byr", s -> parse_num_range(s, 1920, 2002)),
    ("iyr", s -> parse_num_range(s, 2010, 2020)),
    ("eyr", s -> parse_num_range(s, 2020, 2030)),
    ("hgt", s -> begin
       m = match(r"(\d+)(cm|in)", s)
       if m === nothing || length(m.captures) < 2
         return false
       end

       hgt, units = m.captures[1], m.captures[2]
       @match units begin
         "cm" => parse_num_range(hgt, 150, 193)
         "in" => parse_num_range(hgt, 59, 76)
         _ => false
       end
     end),
    ("hcl", s -> match(r"^#[0-9a-f]{6}$", s) !== nothing),
    ("ecl", s -> s ∈ ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"]),
    ("pid", s -> match(r"^\d{9}$", s) !== nothing)
  ])

  all(k -> haskey(d, k) && (validator[k])(d[k]), keys(validator))
end

function main()
  res = 0
  passport_lines = String[]
  for line in eachline("./input")
    if match(r"^\W*$", line) !== nothing
      res += validate_passport(join(passport_lines, ' '))
      passport_lines = []
      continue
    end

    push!(passport_lines, line)
  end

  # process last passport
  res += validate_passport(join(passport_lines, ' '))
  println(res)
end

main()

ext: jl


text:
#=
Be careful, these functions are not all coded to be used in a neural network, especially the functions that return 0 or 1.
They return a value and not an array.
However, you can easily modify them and adapt them to your project, as they remain very general.

Coded by Adrien I, 12/27/2020
=#

function sigmoid(z)
    return 1 ./ (1 .+ exp.(-z))
end

function sigmoid_derivative(z)
    return sigmoid(z) .* (1 .- sigmoid(z))
end

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

function tanh(z)
    return 2 ./ ((1 .+ exp.(-2 .* z)) .+ 1 )
end

function tanh_derivative(z)
    return 1 .- tanh(z) .^ 2
end

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

function relu(x)
    if x < 0
        return 0
    end
    if x >= 0
        return x
    end
end

function derivative_relu(x)
    if x < 0
        return 0
    end
    if x >= 0
        return 1
    end
end

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

function prelu(alpha, x)
    if x < 0
        return alpha .* x
    end
    if x >= 0
        return x
    end
end

function derivative_prelu(alpha, x)
    if x < 0
        return alpha
    end
    if x >= 0
        return 1
    end
end

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

function softmax(x)
    return exp.(x) ./ sum(exp.(x))
end

function derivative_softmax(x)
    return (exp.(x[1]) .* (exp.(x[2]) .+ exp.(x[3]))) ./ ((exp.(x[1]) .+ exp.(x[2]) .+ exp.(x[3])) .^ 2 )
end

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

function elu(alpha, x)
    if x < 0
        return alpha .* (exp.(x) .- 1 )
    end
    if x >= 0
        return x
    end
end

function derivative_elu((alpha, x))
    if x < 0
        return elu(alpha, x) .+ alpha
    end
    if x >= 0
        return 1
    end
end

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

function identity(x)
    return x
end

function derivative_identity(x)
    return 1
end

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

function binary(x)
    if x < 0
        return 0
    end
    if x >= 0
        return 1
    end
end

function derivative_binary(x)
    if x != 0
        return 0
    end
    if x == 0
        return nothing
    end
end


text:
import Base.∈

export constrained_dimensions,
       tosimplehrep,
       remove_redundant_constraints,
       remove_redundant_constraints!,
       linear_map

"""
    ∈(x::AbstractVector{N}, P::AbstractPolyhedron{N})::Bool where {N<:Real}

Check whether a given point is contained in a polyhedron.

### Input

- `x` -- point/vector
- `P` -- polyhedron

### Output

`true` iff ``x ∈ P``.

### Algorithm

This implementation checks if the point lies inside each defining half-space.
"""
function ∈(x::AbstractVector{N}, P::AbstractPolyhedron{N})::Bool where {N<:Real}
    @assert length(x) == dim(P) "a $(length(x))-dimensional point cannot be " *
        "an element of a $(dim(P))-dimensional set"

    for c in constraints_list(P)
        if dot(c.a, x) > c.b
            return false
        end
    end
    return true
end

"""
    constrained_dimensions(P::AbstractPolyhedron)::Vector{Int} where {N<:Real}

Return the indices in which a polyhedron is constrained.

### Input

- `P` -- polyhedron

### Output

A vector of ascending indices `i` such that the polyhedron is constrained in
dimension `i`.

### Examples

A 2D polyhedron with constraint ``x1 ≥ 0`` is constrained in dimension 1 only.
"""
function constrained_dimensions(P::AbstractPolyhedron)::Vector{Int}
    zero_indices = zeros(Int, dim(P))
    for constraint in constraints_list(P)
        for i in constrained_dimensions(constraint)
            zero_indices[i] = i
        end
    end
    return filter(x -> x != 0, zero_indices)
end

"""
    tosimplehrep(constraints::AbstractVector{LC})
        where {N<:Real, LC<:LinearConstraint{N}}

Return the simple H-representation ``Ax ≤ b`` from a list of linear constraints.

### Input

- `constraints` -- a list of linear constraints

### Output

The tuple `(A, b)` where `A` is the matrix of normal directions and `b` is the
vector of offsets.
"""
function tosimplehrep(constraints::AbstractVector{LC}
                     ) where {N<:Real, LC<:LinearConstraint{N}}
    n = length(constraints)
    if n == 0
        A = Matrix{N}(undef, 0, 0)
        b = Vector{N}(undef, 0)
        return (A, b)
    end
    A = zeros(N, n, dim(first(constraints)))
    b = zeros(N, n)
    @inbounds begin
        for (i, Pi) in enumerate(constraints)
            A[i, :] = Pi.a
            b[i] = Pi.b
        end
    end
    return (A, b)
end

"""
     remove_redundant_constraints!(constraints::AbstractVector{LC};
         [backend]=GLPKSolverLP())::Bool where {N<:Real,
                                                LC<:LinearConstraint{N}}

Remove the redundant constraints of a given list of linear constraints; the list
is updated in-place.

### Input

- `constraints` -- list of constraints
- `backend`     -- (optional, default: `GLPKSolverLP`) numeric LP solver backend

### Output

`true` if the function was successful and the list of constraints `constraints`
is modified by removing the redundant constraints, and `false` only if the
constraints are infeasible.

### Notes

Note that the result may be `true` even if the constraints are infeasible.
For example, ``x ≤ 0 && x ≥ 1`` will return `true` without removing any
constraint.
To check if the constraints are infeasible, use
`isempty(HPolyhedron(constraints)`.

### Algorithm

If there are `m` constraints in `n` dimensions, this function checks one by one
if each of the `m` constraints is implied by the remaining ones.

To check if the `k`-th constraint is redundant, an LP is formulated using the
constraints that have not yet being removed.
If, at an intermediate step, it is detected that a subgroup of the constraints
is infeasible, this function returns `false`.
If the calculation finished successfully, this function returns `true`.

For details, see [Fukuda's Polyhedra
FAQ](https://www.cs.mcgill.ca/~fukuda/soft/polyfaq/node24.html).
"""
function remove_redundant_constraints!(constraints::AbstractVector{LC};
                                       backend=GLPKSolverLP()
                                      )::Bool where {N<:Real,
                                                     LC<:LinearConstraint{N}}

    A, b = tosimplehrep(constraints)
    m, n = size(A)
    non_redundant_indices = 1:m

    i = 1 # counter over reduced constraints

    for j in 1:m    # loop over original constraints
        α = A[j, :]
        Ar = A[non_redundant_indices, :]
        br = b[non_redundant_indices]
        br[i] = b[j] + one(N)
        lp = linprog(-α, Ar, '<', br, -Inf, Inf, backend)
        if lp.status == :Infeasible
            # the polyhedron is empty
            return false
        elseif lp.status == :Optimal
            objval = -lp.objval
            if objval <= b[j]
                # the constraint is redundant
                non_redundant_indices = setdiff(non_redundant_indices, j)
            else
                # the constraint is not redundant
                i = i+1
            end
        else
            error("LP is not optimal; the status of the LP is $(lp.status)")
        end
    end

    deleteat!(constraints, setdiff(1:m, non_redundant_indices))
    return true
end

"""
    remove_redundant_constraints(constraints::AbstractVector{LC};
        backend=GLPKSolverLP()) where {N<:Real, LC<:LinearConstraint{N}}

Remove the redundant constraints of a given list of linear constraints.

### Input

- `constraints` -- list of constraints
- `backend`     -- (optional, default: `GLPKSolverLP`) numeric LP solver backend

### Output

The list of constraints with the redundant ones removed, or an empty set if the
constraints are infeasible.

### Algorithm

See
[`remove_redundant_constraints!(::AbstractVector{<:LinearConstraint{<:Real}})`](@ref)
for details.
"""
function remove_redundant_constraints(constraints::AbstractVector{LC};
                                      backend=GLPKSolverLP()
                                     ) where {N<:Real, LC<:LinearConstraint{N}}
    constraints_copy = copy(constraints)
    if remove_redundant_constraints!(constraints_copy, backend=backend)
        return constraints_copy
    else  # the constraints are infeasible
        return EmptySet{N}()
    end
end

"""
    linear_map(M::AbstractMatrix{N},
               P::AbstractPolyhedron{N};
               check_invertibility::Bool=true,
               cond_tol::Number=DEFAULT_COND_TOL,
               use_inv::Bool=!issparse(M)
               ) where {N<:Real}

Concrete linear map of a polyhedron in constraint representation.

### Input

- `M` -- matrix
- `P` -- abstract polyhedron
- `check_invertibility` -- (optional, deault: `true`) check if the linear map is
                           invertible, in which case this function uses the matrix
                           inverse; if the invertibility check fails, or if
                           this flag is set to `false`, use the vertex representation
                           to compute the linear map (see below for details)
- `cond_tol` -- (optional) tolerance of matrix condition (used to check whether
                the matrix is invertible)
- `use_inv`  -- (optional, default: `false` if `M` is sparse and `true`
                otherwise) whether to compute the full left division through
                `inv(M)`, or to use the left division for each vector; see below

### Output

The type of the result is "as close as possible" to the the type of `P`.
Let `(m, n)` be the size of `M`, where `m ≠ n` is allowed for rectangular maps.

To fix the type of the output to something different than the default value,
consider post-processing the result of this function with a call to a suitable
`convert` method.

In particular, the output depends on the type of `P`, on `m`, and the algorithm
that was used:

- If the vertex-based approach was used:

    - If `P` is a `VPolygon` and `m = 2` then the output is a `VPolygon`.
    - If `P` is a `VPolytope` then the output is a `VPolytope`.
    - Otherwise, the output is an `Interval` if `m = 1`, a `VPolygon` if `m = 2`
      and a `VPolytope` in other cases.

- If the invertibility criterion was used:

    - The types of `HalfSpace`, `Hyperplane`, `Line` and `AbstractHPolygon` are
      preserved.
    - If `P` is an `AbstractPolytope`, then the output is an `Interval` if `m = 1`,
      an `HPolygon` if `m = 2` and an `HPolytope` in other cases.
    - Otherwise, the output is an `HPolyhedron`.

### Algorithm

This function implements two algorithms for the linear map:

- If the matrix ``M`` is invertible (which we check with a sufficient condition),
  then ``y = M x`` implies ``x = \\text{inv}(M) y`` and we transform the
  constraint system ``A x ≤ b`` to ``A \\text{inv}(M) y ≤ b``.
- Otherwise, we transform the polyhedron to vertex representation and apply the map
  to each vertex, returning a polyhedron in vertex representation.

Note that the vertex representation (second approach) is only available if the
polyhedron is bounded. Hence we check boundedness first.

To switch off the check for invertibility, set the option
`check_invertibility=false`. If `M` is not invertible and the polyhedron
is unbounded, this function returns an exception.

The option `use_inv` lets the user control - in case `M` is invertible - if
the full matrix inverse is computed, or only the left division on the normal
vectors. Note that this helps as a workaround when `M` is a sparse matrix, since
the `inv` function is not available for sparse matrices. In this case, either
use the option `use_inv=false` or convert the type of `M` as in
`linear_map(Matrix(M), P)`.

Internally, this function operates at the level of the `AbstractPolyhedron`
interface, but the actual algorithm uses dispatch on the concrete type of `P`,
depending on the algorithm that is used:

-  `_linear_map_vrep(M, P)` if the vertex approach is used
-  `_linear_map_hrep(M, P, use_inv)` if the invertibility criterion is used

New subtypes of the interface should write their own `_linear_map_vrep`
(resp. `_linear_map_hrep`) for special handling of the linear map; otherwise
the fallback implementation for `AbstractPolyhedron` is used (see below).
"""
function linear_map(M::AbstractMatrix{N},
                    P::AbstractPolyhedron{N};
                    check_invertibility::Bool=true,
                    cond_tol::Number=DEFAULT_COND_TOL,
                    use_inv::Bool=!issparse(M)
                   ) where {N<:Real}
    @assert dim(P) == size(M, 2) "a linear map of size $(size(M)) cannot be " *
        "applied to a set of dimension $(dim(P))"

    if !check_invertibility || !isinvertible(M; cond_tol=cond_tol)
        # vertex representation is enforced or the matrix is not invertible => use vertex approach
        return _linear_map_vrep(M, P)
    else
        # matrix M is invertible => use H-rep approach
        # the normal vectors are vec(c.a' * inv(M))
        return _linear_map_hrep(M, P, use_inv)
    end
end

# handle different numeric types
function linear_map(M::AbstractMatrix{NM},
                    P::AbstractPolyhedron{NP};
                    kwargs...) where {NM<:Real, NP<:Real}
    N = promote_type(NM, NP)
    if N != NP
        error("conversion between numeric types of polyhedra not implemented " *
            "yet (see #1181)")
    else
        return linear_map(N.(M), P; kwargs...)
    end
end

function _linear_map_vrep(M::AbstractMatrix{N}, P::AbstractPolyhedron{N}) where {N<:Real}
    if !isbounded(P)
        throw(ArgumentError("the linear map in vertex representation for an " *
            "unbounded set is not defined"))
    end
    require(:Polyhedra; fun_name="linear_map",
            explanation="of a $(typeof(P)) by a non-invertible matrix")
    # since P is bounded, we pass an HPolytope and then convert it to vertex representation
    P = tovrep(HPolytope(constraints_list(P), check_boundedness=false))
    return _linear_map_vrep(M, P)
end

function _linear_map_hrep(M::AbstractMatrix{N}, P::AbstractPolyhedron{N},
                          use_inv::Bool) where {N<:Real}
    constraints = _linear_map_hrep_helper(M, P, use_inv)
    return HPolyhedron(constraints)
end

function _linear_map_hrep_helper(M::AbstractMatrix{N}, P::AbstractPolyhedron{N},
                                 use_inv::Bool) where {N<:Real}
    constraints_P = constraints_list(P)
    constraints_MP = similar(constraints_P)
    if use_inv
        invM = inv(M)
        @inbounds for (i, c) in enumerate(constraints_P)
            constraints_MP[i] = LinearConstraint(vec(_At_mul_B(c.a, invM)), c.b)
        end
    else
        # take left division for each constraint c, transpose(M) \ c.a
        @inbounds for (i, c) in enumerate(constraints_P)
            constraints_MP[i] = LinearConstraint(_At_ldiv_B(M, c.a), c.b)
        end
    end
    return constraints_MP
end

"""
    plot_recipe(P::AbstractPolyhedron{N}, [ε]::N=zero(N)) where {N<:Real}

Convert a (bounded) polyhedron to a pair `(x, y)` of points for plotting.

### Input

- `P` -- bounded polyhedron
- `ε` -- (optional, default: `0`) ignored, used for dispatch

### Output

A pair `(x, y)` of points that can be plotted.

### Algorithm

We first assert that `P` is bounded (i.e., that `P` is a polytope).

One-dimensional polytopes are converted to an `Interval`.
Three-dimensional or higher-dimensional polytopes are not supported.

For two-dimensional polytopes (i.e., polygons) we compute their set of vertices
using `vertices_list` and then plot the convex hull of these vertices.
"""
function plot_recipe(P::AbstractPolyhedron{N}, ε::N=zero(N)) where {N<:Real}
    @assert dim(P) <= 2 "cannot plot a $(dim(P))-dimensional $(typeof(P))"
    @assert isbounded(P) "cannot plot an unbounded $(typeof(P))"

    if dim(P) == 1
        return plot_recipe(convert(Interval, P), ε)
    else
        vlist = transpose(hcat(convex_hull(vertices_list(P))...))
        if isempty(vlist)
            @warn "received a polyhedron with no vertices during plotting"
            return plot_recipe(EmptySet{N}(), ε)
        end
        x, y = vlist[:, 1], vlist[:, 2]

        if length(x) > 1
            # add first vertex to "close" the polygon
            push!(x, x[1])
            push!(y, y[1])
        end

        return x, y
    end
end


text:
function isweyl(lgir::LGIrrep{3}; timereversal::Bool=true)
    irdim(lgir) ≠ 2 && return false # must be 2D irrep

    H = kdotp(lgir; timereversal)

    # must contain nonzero σ₁, σ₂, σ₃ terms, but no σ₀ terms
    #(length(H.hs) ≥ 3 && I(2) ∉ H.hs) || return false
    (length(H.hs) ≥ 3 && [0 1; 1 0] ∈ H.hs && [0 -1im; 1im 0] ∈ H.hs && [1 0; 0 -1] ∈ H.hs) || return false

    # find the "generalized" Weyl form of the k ⋅ p Hamiltonian, H = (Vk) ⋅ σ, with matrix V
    V = zeros(3,3)
    for a in eachindex(H.cs)
        # **some** prefactor value; any (not-equal & not-too-simply-related) should be OK...
        x = .1*a^(1.25)
        for d in 1:3
            V[:,d] .+= x.*H.cs[a][d][1:3]
        end
    end

    # if V is invertible, the Hamiltonian can be transformed to a coordinate setting where
    # it has the standard form of the Weyl Hamiltonian, i.e. k ⋅ σ
    return rank(V) == 3
end

function is_spinone_weyl(lgir::LGIrrep{3}; timereversal::Bool=true)
    # TODO
    # check whether H is unitarily equivalent to a matrix of the form
    #   H = k ⋅ S
    # with S = [L₁, L₂, L₃] where Lᵢ are any 3×3 matrices that fulfil the commutation
    # relation [Lᵢ,Lⱼ] = iϵᵢⱼₖLₖ. This conceptual classification can probably be generalized
    # to any irrep dimension. There's a quick discussion of the Chern number of the lowest
    # band of a spin-N particle in Girvin & Yang's book, p. 340-341.
end


text:
## Sandboxing for CloudgraphsDFG
using Revise
using Neo4j
using DistributedFactorGraphs
using IncrementalInference
using Test
using RoME

# DANGER: Clear everything from robot+session + Neo4j database
@testset "Setup and clearing the existing graph" begin
    # Create connection
    # Caching: Off
    global cgDFG = CloudGraphsDFG{NoSolverParams}("localhost", 7474, "neo4j", "test",
        "testUser", "testRobot", "testSession",
        nothing,
        nothing,
        IncrementalInference.decodePackedType)

    clearRobot!!(cgDFG)
    @test any(cgDFG.variableCache) == false
    @test any(cgDFG.factorCache) == false
    @test any(cgDFG.labelDict) == false
    @test any(cgDFG.addHistory) == false
end

@testset "Adding/getting/updating variables" begin
    if exists(cgDFG, :x0)
        vExisting = deleteVariable!(cgDFG, :x0)
        @test vExisting.label == :x0
    end
    @test !exists(cgDFG, :x0)
    @test addVariable!(cgDFG, :x0, Pose2) != nothing
    # Shouldn't throw error, should just update
    @test addVariable!(cgDFG, :x0, Pose2) != nothing
    @test exists(cgDFG, :x0)
    @test getVariable(cgDFG, :x0) != nothing
    @test getVariableIds(cgDFG) == [:x0]
    @test getVariableIds(cgDFG, r"x.*") == [:x0]
    @test getVariableIds(cgDFG, r"y.*") == []
    variable = getVariable(cgDFG, :x0)
    push!(variable.tags, :EXTRAEXTRA)
    updateVariable!(cgDFG, variable)
    variableBack = getVariable(cgDFG, :x0,)
    @test variableBack.tags == variable.tags
    @test deleteVariable!(cgDFG, :x0).label == :x0

    # Things that should pass
    helloVariable = addVariable!(cgDFG, :hellothisIsALabel_Something, Pose2)
    @test helloVariable != nothing
    @test getVariableIds(cgDFG, r".*thisIsALabel.*") == [:hellothisIsALabel_Something]
    @test deleteVariable!(cgDFG, helloVariable).label == helloVariable.label
end

# TODO: Comparisons
# Comparisons

# Factors
@testset "Building a decent factor graph" begin
    x1 = addVariable!(cgDFG, :x1, Pose2)
    x2 = addVariable!(cgDFG, :x2, Pose2)
    x3 = addVariable!(cgDFG, :x3, Pose2)
    l1 = addVariable!(cgDFG, :l1, Pose2)
    l2 = addVariable!(cgDFG, :l2, Pose2)

    prior = PriorPose2( MvNormal([10; 10; pi/6.0], Matrix(Diagonal([0.1;0.1;0.05].^2))))
    addFactor!(cgDFG, [:x1], prior )

    pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
    p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))
    x1x2f1 = addFactor!(cgDFG, [:x1, :x2], pp, autoinit=false)
    x2x3f1 = addFactor!(cgDFG, [:x2, :x3], pp, autoinit=false)
    x1l1f1 = addFactor!(cgDFG, [:x1, :l1], p2br, autoinit=false)
    x2l1f1 = addFactor!(cgDFG, [:x2, :l1], p2br, autoinit=false)
    x3l2f1 = addFactor!(cgDFG, [:x3, :l2], p2br, autoinit=false)

    @test getVariable(cgDFG, :x1) != nothing
    @test getFactor(cgDFG, :x1x2f1) != nothing

    # Testing variable orderings
    @test getFactor(cgDFG, :x2x3f1)._variableOrderSymbols ==[:x2, :x3]
end

@testset "ls() and getNeighbors() tests" begin
    x2 = getVariable(cgDFG, :x2)

    # Not fast though, TODO: https://github.com/JuliaRobotics/DistributedFactorGraphs.jl/issues/39
    vars = getVariables(cgDFG)
    vars2 = ls(cgDFG)
    @test map(v->v.label, vars) == vars2
    @test symdiff(map(v->v.label, vars), [:x1, :x2, :x3, :l1, :l2]) == []

    facts = getFactors(cgDFG)
    facts2 = lsf(cgDFG)
    @test map(v->v.label, facts) == facts2
    @test symdiff(map(f->f.label, facts), [:x1x2f1, :x1f1, :x2l1f1, :x2x3f1, :x3l2f1, :x1l1f1]) == []

    @test symdiff(getNeighbors(cgDFG, :x2), [:x2l1f1, :x2x3f1, :x1x2f1]) == []
    @test getNeighbors(cgDFG, :x2) == getNeighbors(cgDFG, x2)
    @test lsf(cgDFG, :x2) == getNeighbors(cgDFG, :x2)

    # TODO: Test ready and backendset filtering
end

# Test update
# Variables
@testset "Update testing" begin
    x1x2f1 = getFactor(cgDFG, :x1x2f1)
    x1 = getVariable(cgDFG, :x1)
    x2 = getVariable(cgDFG, :x2)
    x3 = getVariable(cgDFG, :x3)

    ve = VariableEstimate([0.1, 0.2, 0.3], :testKey, :testKey)
    push!(x1.estimateDict, :testKey => ve)
    updateVariable!(cgDFG, x1)
    x1Ret = getVariable(cgDFG, :x1)
    @test haskey(x1Ret.estimateDict, :testKey)
    @test x1Ret.estimateDict[:testKey].estimate == ve.estimate

    # Factors: Just the body
    updateFactor!(cgDFG, x1x2f1)
    # Factors: And now the links
    updateFactor!(cgDFG, [x1, x2, x3], x1x2f1)
    updateFactor!(cgDFG, [:x1, :x2, :x3], x1x2f1)
    @test symdiff(getNeighbors(cgDFG, x1x2f1.label), [:x1, :x2, :x3]) == []
end

@testset "copySession!() tests" begin
# Should be able to copy
    cgDFGCopy = copySession!(cgDFG)
    @info "Copied session destination = $(cgDFGCopy.sessionId)"
    @test symdiff(map(v->v.label, getVariables(cgDFGCopy)), map(v->v.label, getVariables(cgDFG))) == []
    @test symdiff(map(f->f.label, getFactors(cgDFGCopy)), map(f->f.label, getFactors(cgDFG))) == []

    clearSession!!(cgDFGCopy)
end

@testset "getSubgraphAroundNode() tests" begin
    cgDFGCopy = getSubgraphAroundNode(cgDFG, getVariable(cgDFG, :x2), 2)
    @info "Subgraph session: $(cgDFGCopy.sessionId)"
    # Checking that the extents do have links (was an issue)
    @test symdiff(getVariableIds(cgDFGCopy), [:l1, :x1, :x3, :x2]) == []
    @test symdiff(getFactorIds(cgDFGCopy), [:x1x2f1, :x2x3f1, :x2l1f1]) == []
    @test getNeighbors(cgDFGCopy, :l1) == [:x2l1f1]
    clearSession!!(cgDFGCopy)

    # Right, now trying to update and produce same result
    cgDFGCopy = getSubgraphAroundNode(cgDFG, getVariable(cgDFG, :x2), 1, true)
    @info "Subgraph session: $(cgDFGCopy.sessionId)"
    @test symdiff(getVariableIds(cgDFGCopy), [:x2]) == []
    @test symdiff(getFactorIds(cgDFGCopy), [:x1x2f1, :x2x3f1, :x2l1f1]) == []

    # Increase range to 2 and update existing
    cgDFGCopy = getSubgraphAroundNode(cgDFG, getVariable(cgDFG, :x2), 2, true, cgDFGCopy)
    @info "Subgraph session: $(cgDFGCopy.sessionId)"
    @test symdiff(getVariableIds(cgDFGCopy), [:l1, :x1, :x3, :x2]) == []
    @test symdiff(getFactorIds(cgDFGCopy), [:x1x2f1, :x2x3f1, :x2l1f1]) == []
    # Check that we copied a linked subset of the graph
    @test getNeighbors(cgDFGCopy, :l1) == [:x2l1f1]
    @test symdiff(getNeighbors(cgDFGCopy, :x1x2f1), [:x1, :x2, :x3]) == []

    # Final cleanup
    # clearSession!!(cgDFGCopy)
end

@testset "getAdjacencyMatrix test" begin
    adjMat = getAdjacencyMatrix(cgDFG)
    # Note that by this point, :x3 is bound to :x1x2f1, which is counterintuitive but that's tests for you
    expected =
    [nothing  :l1      :l2      :x1      :x2      :x3;
    :x1f1    nothing  nothing  :x1f1    nothing  nothing;
    :x1l1f1  :x1l1f1  nothing  :x1l1f1  nothing  nothing;
    :x1x2f1  nothing  nothing  :x1x2f1  :x1x2f1  :x1x2f1;
    :x2l1f1  :x2l1f1  nothing  nothing  :x2l1f1  nothing;
    :x2x3f1  nothing  nothing  nothing  :x2x3f1  :x2x3f1;
    :x3l2f1  nothing  :x3l2f1  nothing  nothing  :x3l2f1]

    @test adjMat == expected
end

@testset "isFullyConnected, hasOrphans tests" begin
    @test isFullyConnected(cgDFG) == true
    @test hasOrphans(cgDFG) == false

    # Now let's add an orphan
    lOrphan = addVariable!(cgDFG, :lOrphan, Pose2)
    @test isFullyConnected(cgDFG) == false
    @test hasOrphans(cgDFG) == true
end

# Show it
# DFG.toDotFile(dfg, "/tmp/testRmMarg.dot")


text:
using MIPforQP

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

using JuMP, BARON, Gurobi, LinearAlgebra

const TIME_LIMIT = 5 * 60.0

const NUM_LAYERS = 3
const SOC_LOWER_BOUND = false
const EPIGRAPH = false
const FAMILIES = ("small", "basic", "extended", "extended2")

const RESULT_PATH = joinpath(@__DIR__, "result.csv")
const ROOT_DIR = joinpath(@__DIR__, "box_qp", "data")

function gurobi_direct_factory(cutoff)
    return JuMP.optimizer_with_attributes(
        Gurobi.Optimizer,
        "TimeLimit" => TIME_LIMIT,
        "NonConvex" => 2,
        "Cutoff" => cutoff
    )
end
function gurobi_direct_mipfocus_factory(cutoff, mipfocus::Int)
    return JuMP.optimizer_with_attributes(
        Gurobi.Optimizer,
        "TimeLimit" => TIME_LIMIT,
        "NonConvex" => 2,
        "MIPFocus" => mipfocus,
        "Cutoff" => cutoff
    )
end
function gurobi_mip_factory(cutoff)
    return JuMP.optimizer_with_attributes(
        Gurobi.Optimizer,
        "TimeLimit" => TIME_LIMIT,
        "MIPFocus" => 3,
        "Cutoff" => cutoff
    )
end
function baron_direct_factory(cutoff)
    return JuMP.optimizer_with_attributes(
        BARON.Optimizer,
        # NOTE: To plug in CPLEX as LP/MIP solver, uncomment the following line and fill in the correct path.
        # "CplexLibName" => "/path/to/cplex/library",
        "MaxTime" => TIME_LIMIT,
        "CutOff" => cutoff
    )
end
const primal_factory = JuMP.optimizer_with_attributes(
    Gurobi.Optimizer,
    "TimeLimit" => TIME_LIMIT,
    "NonConvex" => 2,
    "MIPFocus" => 1
)

function nn_bridge(model)
    MIPforQP.formulate_quadratics!(
        model,
        MIPforQP.Reformulation(MIPforQP.NeuralNetRelaxation(NUM_LAYERS), SOC_LOWER_BOUND)
    )
    return
end
function hongbo_bridge(model)
    MIPforQP.formulate_quadratics!(
        model,
        MIPforQP.Reformulation(MIPforQP.HongboRelaxation(NUM_LAYERS), SOC_LOWER_BOUND)
    )
    return
end

const FACTORIES_AND_BRIDGES = Dict(
    ("gurobi", "direct") => (gurobi_direct_factory, identity),
    ("gurobi_mf_3", "direct") => (co -> gurobi_direct_mipfocus_factory(co, 3), identity),
    ("gurobi", "nn")     => (gurobi_mip_factory, nn_bridge),
    ("gurobi", "hongbo") => (gurobi_mip_factory, hongbo_bridge),
    ("baron", "direct")  => (baron_direct_factory, identity),
)

const CONFIGURATIONS = collect(keys(FACTORIES_AND_BRIDGES))

const CONFIGS_TO_ANALYZE = [
    ("gurobi", "direct"),
    ("gurobi_mf_3", "direct"),
    ("gurobi", "nn"),
    ("gurobi", "hongbo"),
    ("baron", "direct"),
]

max_stars_repo_name: joehuchette/compact-mips-for-qp

max_issues_repo_name: joehuchette/compact-mips-for-qp

max_forks_repo_name: joehuchette/compact-mips-for-qp


text:
module Backup

using ArgParse
using DataStructures
using Crayons

const RESET = crayon"reset"

include("config.jl")
include("core.jl")
include("cli.jl")

function get_value(config, args, key; required=true, msg=nothing)
    val = something(args[key], config[Symbol(key)], Some(nothing))
    required && isnothing(val) && error(isnothing(msg) ? "Value required for '$key', but not found in configuration. Please provide it as an argument, or edit your configuration." : msg)
    val
end

function override!(config, args, key)
    val = get_value(config, args, key)
    if !isnothing(val)
        config[Symbol(key)] = val
    end
end

function main()
    args = parse_args(s)
    config = read_config()
    cmd = args["%COMMAND%"]
    cmd_args = args[cmd]
    _get_value(x; kwargs...) = get_value(config, cmd_args, x; kwargs...)
    if cmd == "make"
        config = BackupConfig(
            _get_value("target"),
            _get_value("dst"),
        )
        confirm = !cmd_args["no-confirm"]
        make_backup(config; confirm)
    elseif cmd == "list"
        dst = _get_value("dst")
        backups = existing_backups(dst)
        if isempty(backups)
            println("No backups found.")
            exit(0)
        else
            println("Existing backups:")
        end

        map(backups) do id
            path = backup_path(dst, id)
            println(crayon"cyan", lpad(id, 4), RESET, " ───→ ", crayon"#ffddaa", path, RESET, "    ", String(readchomp(`stat -c %y $path`)))

        end
    elseif cmd == "config"
        override!(config, cmd_args, "dst")
        save_config(config)
    end
end

export
        BackupConfig,
        backup,
        main

end


text:
using Images

# load images from disk
sample_image_path = "sample-images/cats-3061372_640.jpg";
sample_image = nothing

if isfile(sample_image_path)
    sample_image = load(sample_image_path);
else
    info("ERROR: Image not found!")
end

# load images from url
image_url = "https://cdn.pixabay.com/photo/2018/01/04/18/58/cats-3061372_640.jpg?attachment"
downloaded_image_path = download(image_url)
downloaded_image = load(downloaded_image_path)

# download image to a predefined location
image_url = "https://cdn.pixabay.com/photo/2018/01/04/18/58/cats-3061372_640.jpg?attachment"
downloaded_image_path = download(image_url, 'sample-images/cats-3061372_640.jpg')
downloaded_image = load(downloaded_image_path)

# read directory with images
directory_path = "sample-images";
directory_files = readdir(directory_path);
directory_images = filter(x -> ismatch(r"\.(jpg|png|gif){1}$"i, x), directory_files);

for image_name in directory_images
    image_path = joinpath(directory_path, image_name);
    image = load(image_path);
    # other operations
end


text:
#
# State vector descriptions.
#

"""
A supertype for all states in astrodynamics.
"""
abstract type AbstractState{F, LU, TU, AU, T} <: ParameterizedLabelledArray{F, 1, T, LArray{F, 1, MVector{6, F}, T}} end

"""
A supertype for all state representations in astrodynamics.
"""
abstract type StateVector{F, LU, TU, AU, T} <: AbstractState{F, LU, TU, AU, T} end

"""
A supertype for all state representations with local linearizations in astrodynamics.
"""
abstract type StateVectorWithSTM{F, LU, TU, AU, T} <: AbstractState{F, LU, TU, AU, T} end

"""
Returns the lengthunit of the state vector.
"""
Base.@pure lengthunit(::AbstractState{F, LU, TU, AU}) where {F, LU, TU, AU} = LU

"""
Returns the timeunit of the state vector.
"""
Base.@pure timeunit(::AbstractState{F, LU, TU, AU}) where {F, LU, TU, AU} = TU

"""
Returns the angularunit of the state vector.
"""
Base.@pure angularunit(::AbstractState{F, LU, TU, AU}) where {F, LU, TU, AU} = AU

"""
Returns the angularunit of the state vector.
"""
velocityunit(::AbstractState{F, LU, TU, AU}) where {F, LU, TU, AU} = LU/TU

"""
The length of any `StateVector` is 6!
"""
Base.@pure Base.length(::StateVector) = 6

"""
The size of any `StateVector` is (6,)!
"""
Base.@pure Base.size(::StateVector) = (6,)

"""
The length of any `StateVector` is 6!
"""
Base.@pure Base.length(::StateVectorWithSTM) = 42

"""
The size of any `StateVector` is (6,)!
"""
Base.@pure Base.size(::StateVectorWithSTM) = (42,)

"""
A Cartesian state vector with length 6. Internally
uses `MVector` and `LVector` to store data.
Data is accessible via labels, which are
`x, y, z, ẋ, ẏ, ż, r, v`, where `r`
access `x,y,z` and `v` accesses `ẋ, ẏ, ż`.
"""
mutable struct CartesianState{F, LU, TU, AU} <: StateVector{F, LU, TU, AU, (x = 1, y = 2, z = 3, ẋ = 4, ẏ = 5, ż = 6, r = 1:3, v = 4:6)}
    __rawdata::LArray{F, 1, MVector{6, F}, (x = 1, y = 2, z = 3, ẋ = 4, ẏ = 5, ż = 6, r = 1:3, v = 4:6)}

    """
    $(SIGNATURES)

    Constructs a `CartesianState` with types specified.
    """
    CartesianState{F, LU, TU, AU}(data) where {F, LU, TU, AU} = new{F,LU,TU,AU}(LArray{F, 1, MVector{6, F}, (x = 1, y = 2, z = 3, ẋ = 4, ẏ = 5, ż = 6, r = 1:3, v = 4:6)}(data))
end

"""
Constructs a `CartesianState` from provided position and velocity vectors.
"""
function CartesianState(r::AbstractArray, v::AbstractArray;
                        lengthunit=(unit(eltype(r)) isa Unitful.Length ? unit(eltype(r)) : u"km"),
                        timeunit=(unit(eltype(v)) isa Unitful.Velocity ? lengthunit / unit(eltype(v)) : u"s"),
                        angularunit=u"°")
    rr = (eltype(r) isa Unitful.Length ? ustrip.(lengthunit, r) : r)
    vv = (eltype(v) isa Unitful.Velocity ? ustrip.(lengthunit / timeunit, v) : v)
    F = promote_type(eltype(rr), eltype(vv))
    F isa AbstractFloat || (F = Float64)
    return CartesianState(vcat(rr,vv); lengthunit = lengthunit, timeunit = timeunit, angularunit = angularunit)
end

"""
Constructs a `CartesianState`.
"""
function CartesianState(statevector; lengthunit=u"km", timeunit=u"s", angularunit=u"°")
    F = eltype(statevector)
    F isa AbstractFloat || (F = Float64)

    return CartesianState{F, lengthunit, timeunit, angularunit}(
        LArray{F, 1, MVector{6, F}, (x = 1, y = 2, z = 3, ẋ = 4, ẏ = 5, ż = 6, r = 1:3, v = 4:6)}(
            statevector
        )
    )
end

"""
A Cartesian state vector with local linearization
and length 42. Internally
uses `MVector` and `LVector` to store data.
Data is accessible via labels, which are
`x, y, z, ẋ, ẏ, ż, r, v`, where `r`
access `x,y,z` and `v` accesses `ẋ, ẏ, ż`.
"""
mutable struct CartesianStateWithSTM{F, LU, TU, AU} <: StateVectorWithSTM{F, LU, TU, AU, (x = 1, y = 2, z = 3, ẋ = 4, ẏ = 5, ż = 6, r = 1:3, v = 4:6, Φ=(7:42))}
    __rawdata::LArray{F, 1, MVector{42, F}, (x = 1, y = 2, z = 3, ẋ = 4, ẏ = 5, ż = 6, r = 1:3, v = 4:6, Φ=(7:42))}

    """
    $(SIGNATURES)

    Constructs a `CartesianStateWithSTM` with types specified.
    """
    CartesianStateWithSTM{F, LU, TU, AU}(data) where {F, LU, TU, AU} = new{F,LU,TU,AU}(LArray{F, 1, MVector{42, F}, (x = 1, y = 2, z = 3, ẋ = 4, ẏ = 5, ż = 6, r = 1:3, v = 4:6, Φ=(7:42))}(data))
end

"""
Outer constructor for `CartesianStateWithSTM`.
"""
CartesianStateWithSTM(data; lengthunit=u"km", timeunit=u"s", angularunit=u"°") = CartesianStateWithSTM{eltype(data), lengthunit, timeunit, angularunit}(data)

"""
All Cartesian state!
"""
const CartesianStateVector = Union{<:CartesianState, <:CartesianStateWithSTM}

"""
Returns `x`.
"""
get_x(state::CartesianStateVector) = state[1]

"""
Returns `y`.
"""
get_y(state::CartesianStateVector) = state[2]

"""
Returns `z`.
"""
get_z(state::CartesianStateVector) = state[3]

"""
Returns `ẋ`.
"""
get_ẋ(state::CartesianStateVector) = state[4]

"""
Returns `ẏ`.
"""
get_ẏ(state::CartesianStateVector) = state[5]

"""
Returns `ż`.
"""
get_ż(state::CartesianStateVector) = state[6]

"""
Returns `r`.
"""
get_r(state::CartesianStateVector) = state[1:3]

"""
Returns `v`.
"""
get_v(state::CartesianStateVector) = state[4:6]

"""
Returns `ϕ`.
"""
get_ϕ(state::CartesianStateWithSTM) = state[7:42]

"""
Returns the state transition matrix.
"""
get_stm(state::CartesianStateWithSTM) = MMatrix{6,6}(state[7:42])

"""
Returns the whole state vector, without units.
"""
statevector(state::CartesianState) = MVector{6}(get_x(state), get_y(state), get_z(state), get_ẋ(state), get_ẏ(state), get_ż(state))

"""
Returns the whole state vector, without units.
"""
statevector(state::CartesianStateWithSTM) = MVector{42}(get_x(state), get_y(state), get_z(state), get_ẋ(state), get_ẏ(state), get_ż(state), get_ϕ(state)...)

"""
Displays a `CartesianState`.
"""
function Base.show(io::IO, state::CartesianState; showfloats=true, space="")
    println(io, space, "Cartesian State", showfloats ? " with eltype $(eltype(state))" : "", "\n")
    println(io, space, "  ", "x = $(get_x(state) * lengthunit(state))")
    println(io, space, "  ", "y = $(get_y(state) * lengthunit(state))")
    println(io, space, "  ", "z = $(get_z(state) * lengthunit(state))")
    println(io, space, "  ", "ẋ = $(get_ẋ(state) * velocityunit(state))")
    println(io, space, "  ", "ẏ = $(get_ẏ(state) * velocityunit(state))")
    println(io, space, "  ", "ż = $(get_ż(state) * velocityunit(state))")
end
Base.show(io::IO, ::MIME"text/plain", state::CartesianState; showfloats=true, space="") = show(io, state; showfloats=showfloats, space=space)

"""
Displays a `CartesianStateWithSTM`.
"""
function Base.show(io::IO, state::CartesianStateWithSTM; showfloats=true, space="")
    println(io, space, "Cartesian State with local linearization", showfloats ? " and eltype $(eltype(state))" : "", "\n")
    println(io, space, "  ", "x = $(get_x(state) * lengthunit(state))")
    println(io, space, "  ", "y = $(get_y(state) * lengthunit(state))")
    println(io, space, "  ", "z = $(get_z(state) * lengthunit(state))")
    println(io, space, "  ", "ẋ = $(get_ẋ(state) * velocityunit(state))")
    println(io, space, "  ", "ẏ = $(get_ẏ(state) * velocityunit(state))")
    println(io, space, "  ", "ż = $(get_ż(state) * velocityunit(state))")
end
Base.show(io::IO, ::MIME"text/plain", state::CartesianStateWithSTM; showfloats=true, space="") = show(io, state; showfloats=showfloats, space=space)

"""
A Keplerian state vector with length 6. Internally
uses `MVector` and `LVector` to store data.
Data is accessible via labels, which are
`e, a, i, Ω, ω, ν`.
"""
mutable struct KeplerianState{F, LU, TU, AU} <: StateVector{F, LU, TU, AU, (e=1, a=2, i=3, Ω=4, ω=5, ν=6)}
    __rawdata::LArray{F, 1, MVector{6, F}, (e=1, a=2, i=3, Ω=4, ω=5, ν=6)}

    function KeplerianState(statevector; lengthunit=u"km", timeunit=u"s", angularunit=u"°")
        F = eltype(statevector)
        F isa AbstractFloat || (F = Float64)

        return new{F, lengthunit, timeunit, angularunit}(
            statevector
        )
    end
end

"""
Constructs a `KeplerianState` from any `AbstractVector`.
"""
KeplerianState{F, LU, TU, AU}(statevector::AbstractVector{<:Real}) where {F, LU, TU, AU} = KeplerianState(
    convert(LArray{F, 1, MVector{6, F}, (e=1, a=2, i=3, Ω=4, ω=5, ν=6)}, statevector)
)

"""
Constructs a `KeplerianState` from provided position and velocity vectors.
"""
function KeplerianState(e::Real, a, i, Ω, ω, ν;
                        lengthunit=(a isa Unitful.Length ? unit(a) : u"km"),
                        timeunit=u"s",
                        angularunit=(i isa DimensionlessQuantity ? unit(i) : u"°"))
    aa = (a isa Unitful.Length ? ustrip(lengthunit, a) : a)
    ii = (i isa Unitful.DimensionlessQuantity ? ustrip(angularunit, i) : i)
    ΩΩ = (Ω isa Unitful.DimensionlessQuantity ? ustrip(angularunit, Ω) : Ω)
    ωω = (ω isa Unitful.DimensionlessQuantity ? ustrip(angularunit, ω) : ω)
    νν = (ν isa Unitful.DimensionlessQuantity ? ustrip(angularunit, ν) : ν)
    F = promote_type(typeof(e), typeof(aa), typeof(ii), typeof(ΩΩ), typeof(ωω), typeof(νν))
    F isa AbstractFloat || (F = Float64)
    return Keplerian{F, lengthunit, timeunit, angularunit}(@LArray(MVector{6,F}(e,a,i,Ω,ω,ν),(e=1, a=2, i=3, Ω=4, ω=5, ν=6)))
end

"""
Returns `e`.
"""
get_e(state::KeplerianState) = state[1]

"""
Returns `a`.
"""
get_a(state::KeplerianState) = state[2]

"""
Returns `i`.
"""
get_i(state::KeplerianState) = state[3]

"""
Returns `Ω`.
"""
get_Ω(state::KeplerianState) = state[4]

"""
Returns `ω.
"""
get_ω(state::KeplerianState) = state[5]

"""
Returns `ν`.
"""
get_ν(state::KeplerianState) = state[6]

"""
Returns the whole state vector, without units.
"""
statevector(state::KeplerianState) = MVector{6}(get_e(state), get_a(state), get_i(state), get_Ω(state), get_ω(state), get_ν(state))

"""
Displays a `KeplerianState`.
"""
function Base.show(io::IO, state::KeplerianState; showfloats=true, space="")
    println(io, space, "Keplerian State", showfloats ? " with eltype $(eltype(state))" : "", "\n")
    println(io, space, "  ", "e = $(get_e(state))")
    println(io, space, "  ", "a = $(get_a(state) * lengthunit(state))")
    println(io, space, "  ", "Ω = $(get_Ω(state) * angularunit(state))")
    println(io, space, "  ", "ω = $(get_ω(state) * angularunit(state))")
    println(io, space, "  ", "ν = $(get_ν(state) * angularunit(state))")
end
Base.show(io::IO, ::MIME"text/plain", state::KeplerianState; showfloats=true, space="") = show(io, state; showfloats=showfloats, space=space)

"""
Returns `true` if the `StateVector` has `lengthunit` and `timeunit`
parameters of some type `T <: AbstractQuantity`. This is
intended to be used for normalizing `CartesianState` vectors.
"""
isnormalized(state::StateVector) = lengthunit(state) isa Unitful.AbstractQuantity && timeunit(state) isa Unitful.AbstractQuantity

"""
Converts types and units for a `CartesianState`.
"""
function Base.convert(::Type{CartesianState{F, LU, TU, AU}}, state::CartesianState) where {F, LU, TU, AU}
    r = OrbitalStates.get_r(state) .* lengthunit(state) ./ LU .|> upreferred .|> F
    v = OrbitalStates.get_v(state) .* velocityunit(state) ./ (LU/TU) .|> upreferred .|> F
    return CartesianState(vcat(r,v); lengthunit=LU, timeunit=TU, angularunit=AU)
end

"""
Converts types and units for a `CartesianState`.
"""
function Base.convert(::Type{KeplerianState{F, LU, TU, AU}}, state::KeplerianState) where {F, LU, TU, AU}
    e = get_e(state) |> F
    a = ustrip(LU, get_a(state) * lengthunit(state))  |> F
    i = ustrip(AU, get_i(state) * angularunit(state)) |> F
    Ω = ustrip(AU, get_Ω(state) * angularunit(state)) |> F
    ω = ustrip(AU, get_ω(state) * angularunit(state)) |> F
    ν = ustrip(AU, get_ν(state) * angularunit(state)) |> F
    return KeplerianState(MVector{6}(e, a, i, Ω, ω, ν); lengthunit=LU, timeunit=TU, angularunit=AU)
end


text:
module ParallelUtilities

using Distributed
using SplittablesBase

export pmapreduce
export pmapreduce_productsplit
export pmapbatch
export pmapbatch_productsplit
export workerrank

include("productsplit.jl")

include("clusterquery.jl")
using .ClusterQueryUtils: procs_node
using .ClusterQueryUtils: workerpool_nodes
using .ClusterQueryUtils: workers_myhost
using .ClusterQueryUtils: workersactive
using .ClusterQueryUtils: maybetrimmedworkerpool

include("trees.jl")
include("reductionfunctions.jl")
include("mapreduce.jl")

end # module


text:
# TargetInfo
function TargetInfo(opts::TargetOptions, diag::DiagnosticsEngine=DiagnosticsEngine())
    @check_ptrs opts diag
    info = clang_TargetInfo_CreateTargetInfo(diag, opts)
    return TargetInfo(info)
end


text:
# ---
# title: 206. Reverse Linked List
# id: problem206
# author: zhwang
# date: 2022-02-17
# difficulty: Easy
# categories: Linked List
# link: <https://leetcode.com/problems/reverse-linked-list/description/>
# hidden: true
# ---
#
# Reverse a singly linked list.
#
# **Example:**
#
#
#
#     Input: 1->2->3->4->5->NULL
#     Output: 5->4->3->2->1->NULL
#
#
# **Follow up:**
#
# A linked list can be reversed either iteratively or recursively. Could you
# implement both?
#
#
## @lc code=start
using LeetCode

reverse_list(head::Nothing) = head
function reverse_list(head::ListNode)::ListNode
    new_head = pre_node = nothing
    while !isnothing(head)
        new_head = ListNode(head.val)
        pre_node, new_head.next = new_head, pre_node
        head = head.next
    end
    return new_head
end
## @lc code=end


text:
# # Relaxed Two-Way Partitioning Problem
#
# We consider the (nonconvex) two-way partitioning problem from Boyd and Vandenberghe (2004), p.219 \[1\]:
# $$
# \begin{array}{ll} \text{minimize} &  x^\top W x \\
# \text{subject to} &  x_i^2 = 1, \quad i=1,\dots, n,
# \end{array}
# $$
# with $x \in \mathbf{R}^n$ and $W \in \mathbf{S}^n$. The problem can be interpreted as finding a partition of the points $i = 1,\dots, n$ into two sets, where the cost of two points in the same set is $W_{ij}$ and $-W_{ij}$ otherwise. Brute forcing the solution $x^*$ takes $2^n$ attempts and becomes quickly intractable, e.g. for $n \geq 30$. However, a lower-bound for this problem can be computed by solving the convex dual of the problem:
# $$
# \begin{array}{ll} \text{maximize} &  -\sum_i \nu \\
# \text{subject to} &  W + \text{diag}(\nu) \in \mathbf{S}_+^n.
# \end{array}
# $$
# Solving this problem with optimal objective $d^*$ yields a lower bound as least as good as a lower bound based on the minimum eigenvalue $d^* \geq \lambda_{\text{min}}(W)$.
# Let's set up the problem for $n = 20$ and with a specially structured (banded) $W$:

#-
using LinearAlgebra, Random, COSMO, JuMP, IterTools

rng = Random.MersenneTwister(212313);
n = 20;

W = diagm(0 => randn(rng, n));
W += diagm(1 => randn(rng, n-1));
W = Symmetric(W)

# As you can see, the matrix $W$ imposes a structure on the LMI-constraint. Accordingly, COSMO will be able to use chordal decomposition to decompose the LMI constraint into multiple smaller constraints. This can make a significant difference in the performance of the algorithm for large $n$.

#-

model = JuMP.Model(COSMO.Optimizer);
@variable(model, ν[1:n]);
@objective(model, Max, -ones(n)' * ν )
@constraint(model, Symmetric(W + diagm(ν) )  in JuMP.PSDCone());
JuMP.optimize!(model)

# Looking at the solver output you can see how the PSD constraint was decomposed into 19 PSD constraints. Let's look at the lower bound:
JuMP.objective_value(model)

# As $n$ is small, we can verify our result by finding the optimal solution by trying out all possible combinations:

function brute_force_optimisation(W, n)
   opt_obj = Inf
   opt_x = Inf * ones(n)

   for xt in Iterators.product([[1.; -1.] for k = 1:n]...)
      x = [i for i in xt]
      obj_val = x' * W * x
      if obj_val < opt_obj
         opt_obj = obj_val
         opt_x = x
      end
   end
   return opt_obj, opt_x
end
opt_obj, opt_x = brute_force_optimisation(W, n)
# ## References
# [1] Boyd and Vandenberghe - Convex Optimization, Cambridge University Press (2004)


text:
using LinearAlgebra
using Random

function test_collocation_matrix()
    # Test special case of 1×1 matrix
    @testset "1×1 matrix" begin
        B = BandedMatrix{Float64}(undef, (1, 1), (0, 0))
        B[1, 1] = 0
        C = @inferred CollocationMatrix(B)
        @test_throws ZeroPivotException(1) lu!(C)
        C[1, 1] = 3
        F = lu(C)
        @test F.L == [1]'
        @test F.U == [3]'
        y = [6]
        x = F \ y
        @test x == [2]
    end

    @testset "Non-square" begin
        B = BandedMatrix{Float32}(undef, (7, 8), (2, 2))
        C = CollocationMatrix(B)
        @test_throws DimensionMismatch lu(C)  # not supported
    end

    # Special cases of triangular input matrices
    @testset "Upper triangular" begin
        U = CollocationMatrix(brand(Float32, 8, 8, 0, 2))
        F = lu(U)
        @test F.L == I
        @test F.U == U

        U[3, 3] = 0
        @test_throws ZeroPivotException(3) lu!(U)
    end

    @testset "Lower triangular" begin
        L = CollocationMatrix(brand(Float32, 8, 8, 2, 0))
        D = Diagonal(L)
        F = lu(L)
        @test F.U == D
        @test F.L * D ≈ L

        L[4, 4] = 0
        @test_throws ZeroPivotException(4) lu!(L)
    end

    nothing
end

function test_collocation(B::BSplineBasis, ::Type{T} = Float64) where {T}
    xcol = collocation_points(B)
    @inferred collocation_matrix(B, xcol, Matrix{T})
    @inferred collocation_matrix(B, xcol, CollocationMatrix{T})
    @inferred collocation_matrix(B, xcol, SparseMatrixCSC{T})

    C_dense = collocation_matrix(B, xcol, Matrix{T})
    C_banded = collocation_matrix(B, xcol, CollocationMatrix{T})
    C_sparse = collocation_matrix(B, xcol, SparseMatrixCSC{T})

    @test C_dense == C_banded
    @test C_banded == C_sparse

    @testset "Greville sites" begin
        xs = collocation_points(B; method = Collocation.AvgKnots())
        it = Collocation.GrevilleSiteIterator(B)
        @test xs == collect(it)
    end

    @testset "Basis recombination" begin
        C = collocation_matrix(B, xcol)
        N = length(B)

        # Recombined basis for Dirichlet BCs (u = 0)
        @inferred RecombinedBSplineBasis(Derivative(0), B)
        R0 = RecombinedBSplineBasis(Derivative(0), B)
        Nr = length(R0)
        @test Nr == N - 2

        # Collocation points for recombined basis, for implicit boundary
        # conditions (points at the boundaries are removed!).
        x = collocation_points(R0)
        @test length(x) == Nr
        @test x ≈ @view xcol[2:end-1]

        @inferred collocation_matrix(R0, x)
        C0 = collocation_matrix(R0, x)
        @test size(C0) == (Nr, Nr)

        # C0 is simply the interior elements of C.
        let r = 2:(N - 1)
            @test C0 ≈ @view C[r, r]
        end

        # Neumann BCs (du/dx = 0)
        R1 = RecombinedBSplineBasis(Derivative(1), B)
        @test collocation_points(R1) == x  # same as for Dirichlet
        C1 = collocation_matrix(R1, x)
        @test size(C1) == (Nr, Nr)

        let r = 2:(N - 1)
            @test @views C1[:, 1] ≈ C[r, 1] .+ C[r, 2]
            @test @views C1[:, Nr] ≈ C[r, N - 1] .+ C[r, N]
            @test @views C1[:, 2:(Nr - 1)] ≈ C[r, 3:(N - 2)]
        end
    end

    nothing
end

@testset "Collocation matrix" begin
    test_collocation_matrix()
end

@testset "Collocation (k = $k)" for k ∈ (4, 5)
    t = gauss_lobatto_points(10 + k)
    B = BSplineBasis(k, t)
    test_collocation(B)
end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
context("decompose functions") do

context("General") do
    # this fails on travis and I can't reproduce it.. It's not the most important
    # test so I uncommented it for now!
    #@fact_throws ArgumentError decompose(Normal{3, Float32}, Circle{Float32}(Point2f0(0), 1f0))
end

context("HyperRectangles") do
    a = HyperRectangle(Vec(0,0),Vec(1,1))
    pt_expa = Point{2,Int}[(0,0), (1,0), (0,1), (1,1)]
    @fact decompose(Point{2,Int},a) --> pt_expa
    b = HyperRectangle(Vec(1,1,1),Vec(1,1,1))
    pt_expb = Point{3,Int}[(1,1,1),(2,1,1),(1,2,1),(2,2,1),(1,1,2),(2,1,2),(1,2,2),(2,2,2)]
    @fact decompose(Point{3,Int}, b) --> pt_expb
end

context("Faces") do
    @fact decompose(GLTriangle, Face{4, Int, 0}(1,2,3,4)) --> (Face{3,UInt32,-1}(0,1,2), Face{3,UInt32,-1}(0,2,3))
    @fact decompose(Face{3,Int,-1}, Face{4, Int, -1}(1,2,3,4)) --> (Face{3,Int,-1}(1,2,3),Face{3,Int,-1}(1,3,4))
    @fact decompose(Face{3,Int,2}, Face{4, Int, 1}(1,2,3,4)) --> (Face{3,Int,2}(2,3,4),Face{3,Int,2}(2,4,5))
    @fact decompose(Face{2,Int,0}, Face{4, Int, 0}(1,2,3,4)) --> (Face{2,Int,0}(1,2),
                                                                  Face{2,Int,0}(2,3),
                                                                  Face{2,Int,0}(3,4),
                                                                  Face{2,Int,0}(4,1))
end

context("Simplex") do
    s1 = Simplex(:x1,:x2,:x3)
    s2 = Simplex(:x1,:x2,:x3,:x4)
    @fact decompose(Simplex{1}, s1) --> (Simplex(:x1),Simplex(:x2),Simplex(:x3))
    @fact decompose(Simplex{2}, s1) --> (Simplex(:x1,:x2),Simplex(:x2,:x3),Simplex(:x3,:x1))
    @fact decompose(Simplex{2}, s2) --> (Simplex(:x1,:x2),Simplex(:x2,:x3),Simplex(:x3,:x4),Simplex(:x4,:x1))
    @fact decompose(Simplex{3}, s2) --> (Simplex(:x1,:x2,:x3),Simplex(:x1,:x3,:x4))
end

context("SimpleRectangle") do
    r = SimpleRectangle(0,0,1,1)
    pts = decompose(Point, r)
    @fact pts --> Point{2,Int}[
        (0,0),
        (1,0),
        (0,1),
        (1,1)
    ]

    mesh = GLNormalUVMesh(r, (3,3))
    points = decompose(Point{3,Float64}, mesh)
    point_target = Point{3,Float64}[
        (0.0,0.0,0.0)
        (0.5,0.0,0.0)
        (1.0,0.0,0.0)
        (0.0,0.5,0.0)
        (0.5,0.5,0.0)
        (1.0,0.5,0.0)
        (0.0,1.0,0.0)
        (0.5,1.0,0.0)
        (1.0,1.0,0.0)
    ]
    @fact points --> point_target

    faces = decompose(Face{3, Int, 0}, mesh)
    face_target = Face{3, Int, 0}[
        (1,2,5)
        (1,5,4)
        (2,3,6)
        (2,6,5)
        (4,5,8)
        (4,8,7)
        (5,6,9)
        (5,9,8)
    ]
    @fact faces --> face_target

    uvs = decompose(UV{Float64}, mesh)
    uv_target = UV{Float64}[
        (0.0,0.0)
        (0.5,0.0)
        (1.0,0.0)
        (0.0,0.5)
        (0.5,0.5)
        (1.0,0.5)
        (0.0,1.0)
        (0.5,1.0)
        (1.0,1.0)
    ]
    @fact uvs --> uv_target
end

context("Normals") do
    n64 = Normal{3, Float64}[
        (0.0,0.0,-1.0),
        (0.0,0.0,-1.0),
        (0.0,0.0,-1.0),
        (0.0,0.0,-1.0),
        (0.0,0.0,1.0),
        (0.0,0.0,1.0),
        (0.0,0.0,1.0),
        (0.0,0.0,1.0),
        (-1.0,0.0,0.0),
        (-1.0,0.0,0.0),
        (-1.0,0.0,0.0),
        (-1.0,0.0,0.0),
        (1.0,0.0,0.0),
        (1.0,0.0,0.0),
        (1.0,0.0,0.0),
        (1.0,0.0,0.0),
        (0.0,1.0,0.0),
        (0.0,1.0,0.0),
        (0.0,1.0,0.0),
        (0.0,1.0,0.0),
        (0.0,-1.0,0.0),
        (0.0,-1.0,0.0),
        (0.0,-1.0,0.0),
        (0.0,-1.0,0.0),
    ]
    n32 = map(Normal{3,Float32}, n64)
    r = GLPlainMesh(centered(HyperRectangle))
    @fact normals(vertices(r), faces(r), Normal{3, Float32}) --> n32
    @fact normals(vertices(r), faces(r), Normal{3, Float64}) --> n64

    r = PlainMesh{Float64, Face{3, UInt32, 0}}(centered(HyperRectangle))
    @fact normals(vertices(r), faces(r), Normal{3, Float32}) --> n32
    @fact normals(vertices(r), faces(r), Normal{3, Float64}) --> n64

    r = PlainMesh{Float16, Face{3, UInt64, -1}}(centered(HyperRectangle))
    @fact normals(vertices(r), faces(r), Normal{3, Float32}) --> n32
    @fact normals(vertices(r), faces(r), Normal{3, Float64}) --> n64

end

context("HyperSphere") do
    sphere = Sphere{Float32}(Point3f0(0), 1f0)

    points = decompose(Point, sphere, 3)
    point_target = Point3f0[
        (0.0,0.0,1.0)
        (0.86602545,0.0,0.49999997)
        (0.8660254,0.0,-0.50000006)
        (0.0,0.0,1.0)
        (-0.43301278,0.75,0.49999997)
        (-0.43301275,0.75,-0.50000006)
        (0.0,0.0,1.0)
        (-0.43301263,-0.75000006,0.49999997)
        (-0.4330126,-0.75,-0.50000006)
        (0.0,0.0,-1.0)
    ]
    @fact points --> point_target

    faces = decompose(GLTriangle, sphere, 3)
    face_target = GLTriangle[
        (0,3,1)
        (1,3,4)
        (3,6,4)
        (4,6,7)
        (6,0,7)
        (7,0,1)
        (1,4,2)
        (2,4,5)
        (4,7,5)
        (5,7,8)
        (7,1,8)
        (8,1,2)
        (2,5,9)
        (9,5,9)
        (5,8,9)
        (9,8,9)
        (8,2,9)
        (9,2,9)
    ]
    @fact faces --> face_target
end

end


text:
function test_Zalg(alg::InducingPoints.OnIPSA, N::Int=30, D::Int=3; kwargs...)
    @testset "Test Vector{Real}" begin
        for T in [Float64, Float32]
            X = rand(T, N) * 10
            Z = inducingpoints(alg, X; arraytype=Vector{T}, kwargs...)
            @test eltype(Z) == T
            X2 = rand(T, 2 * N) * 10
            Z = updateZ!(Z, alg, X2; kwargs...)
            @test eltype(Z) == T
            Z = updateZ(Z, alg, X2; kwargs...)
            @test eltype(Z) == T
        end
    end
    @testset "Test Vector{Vector}" begin
        for T in [Float64, Float32]
            X = [rand(T, D) * 10 for _ in 1:N]
            Z = inducingpoints(alg, X; arraytype=Vector{T}, kwargs...)
            @test eltype(Z) <: AbstractVector{T}
            X2 = [rand(T, D) * 10 for _ in 1:N]
            Z = updateZ!(Z, alg, X2; kwargs...)
            @test eltype(Z) <: AbstractVector{T}
            Z = updateZ(Z, alg, X2; kwargs...)
            @test eltype(Z) <: AbstractVector{T}
        end
    end
    @testset "Test ColVecs" begin
        for T in [Float64, Float32]
            X = ColVecs(rand(T, D, N) * 10)
            Z = inducingpoints(alg, X; arraytype=Vector{T}, kwargs...)
            @test eltype(Z) <: AbstractVector{T}
            X2 = ColVecs(rand(T, D, 2 * N) * 10)
            X3 = [rand(T, D) * 10 for _ in 1:N]
            Z = updateZ!(Z, alg, X2; kwargs...)
            @test eltype(Z) <: AbstractVector{T}
            Z = updateZ!(Z, alg, X3; kwargs...)
            @test eltype(Z) <: AbstractVector{T}
            Z = updateZ(Z, alg, X2; kwargs...)
            @test eltype(Z) <: AbstractVector{T}
        end
    end
end

function test_Zalg(alg::InducingPoints.OffIPSA, N::Int=30, D::Int=3; kwargs...)
    @testset "Test Vector{Real}" begin
        for T in [Float64, Float32]
            X = rand(T, N) * 10
            Z = inducingpoints(alg, X; kwargs...)
            @test Z isa AbstractVector
            @test eltype(Z) == T
        end
    end
    @testset "Test Vector{Vector}" begin
        for T in [Float64, Float32]
            X = [rand(T, D) * 10 for _ in 1:N]
            Z = inducingpoints(alg, X; kwargs...)
            @test Z isa AbstractVector{<:AbstractVector}
            @test first(Z) isa AbstractVector{<:T}
        end
    end
    @testset "Test ColVecs" begin
        for T in [Float64, Float32]
            X = ColVecs(rand(T, D, N) * 10)
            Z = inducingpoints(alg, X; kwargs...)
            @test Z isa ColVecs{T}
        end
    end
end


text:
function animate_boat(sol, bh2o_0)
    # Points in boat coordinate system
    P_boat = [
        0 xmin0 xmax0 xmax
        0 hb hb 0
    ]
    r_o = [xo, zo]
    # Transform to coordinates associated with water level
    r_po = r_o + [0, -bh2o_0]
    # P - point in boat system, y - solution (x, z, θ)
    R(θ) = begin
        sθ, cθ = sincos(θ)
        [cθ -sθ; sθ cθ]
    end
    T = [1 0; 0 -1]
    # From boat reference system to displaced global + reverse of z by T
    to_water(P, x, z, θ) = T * ((r_po + [x, z]) .+ R(θ) * (P .- r_o))

    P_water = [
        -0.1 xmax+0.1 xmax+0.1 -0.1
        0 0 -hb-0.1 -hb-0.1
    ]

    n_frame = 4
    fps = 1 / sol.t[2] / n_frame
    # println("Animaiton with fps=$fps")
    anim = @animate for i = 1:n_frame:length(sol.t)
        x, z, θ = sol[1, i], sol[2, i], sol[3, i]
        P_boat_i = to_water(P_boat, x, z, θ)
        P_water_i = P_water .+ [x, 0]
        P_cm_i = to_water(r_o, x, z, θ)

        plot(
            xlim = (P_water_i[1, 1], P_water_i[1, 2]),
            ylim = (P_water_i[2, 3], 0.2),
            aspect_ratio = 4,
        )
        plot!(Shape(P_boat_i[1, :], P_boat_i[2, :]), opacity = 0.5, color = :orange)
        plot!(Shape(P_water_i[1, :], P_water_i[2, :]), opacity = 0.25, color = :blue)
        plot!(P_cm_i[1, :], P_cm_i[2, :], seriestype = :scatter)
    end
    # gif(anim, "anim_fps$(fps).gif", fps = fps)
    anim, fps
end


text:
abstract type Animal end

name(x::Animal) = x.name

abstract type Monkey <: Animal end
abstract type Bird <: Animal end

struct HugeMonkey <: Monkey
    name::String
    age::Int
    size::Int
end

struct SmallMonkey <: Monkey
    name::String
end

struct BigBird <: Bird
    name::String
    flyable::Bool
end

struct Swoole <: Bird
end

name(x::Swoole) = "Swoole"

eatable(a::Monkey, b::Bird) = true
eatable(a::Bird, b::Monkey) = false


text:
module TestMath

using Base.Test
using RLEVectors

@testset begin

# mode and implicitly countmap
x = RLEVector([4,5,6,5],cumsum([10,5,3,6]))
@test mode(x) == 5

end # testset

end # module


text:
platform(dict::ExtensionFrame) = ExtensionFramePlatform(platform(basis(dict)), support(dict))

export ExtensionFramePlatform
"""
    struct ExtensionFramePlatform <: FramePlatform

A platform that creates extension frames.
"""
struct ExtensionFramePlatform <: FramePlatform
    basisplatform   ::  Platform
    domain          ::  Domain
end

correctparamformat(platform::ExtensionFramePlatform, param) =
    correctparamformat(platform.basisplatform, param)

support(platform::ExtensionFramePlatform) = platform.domain

SolverStyle(dict::ExtensionFrame, ::OversamplingStyle) = AZStyle()

SamplingStyle(p::ExtensionFramePlatform) = OversamplingStyle()

unsafe_dictionary(p::ExtensionFramePlatform, n) =
    extensionframe(p.domain, unsafe_dictionary(p.basisplatform, n))

measure(platform::ExtensionFramePlatform) = restrict(measure(platform.basisplatform), platform.domain)

dualdictionary(platform::ExtensionFramePlatform, plt_par, measure::Measure; options...) =
   extensionframe(dualdictionary(platform.basisplatform, plt_par, supermeasure(measure); options...), platform.domain)

default_param_path(platform::ExtensionFramePlatform) =
   default_param_path(platform.basisplatform)


text:
function update_env_precip_tendencies(state, k, qt_tendency, θ_liq_ice_tendency, qr_tendency, qs_tendency)

    aux_en = center_aux_environment(state)
    tendencies_pr = center_tendencies_precipitation(state)

    aux_en.qt_tendency_precip_formation[k] = qt_tendency * aux_en.area[k]
    aux_en.θ_liq_ice_tendency_precip_formation[k] = θ_liq_ice_tendency * aux_en.area[k]

    tendencies_pr.q_rai[k] += qr_tendency * aux_en.area[k]
    tendencies_pr.q_sno[k] += qs_tendency * aux_en.area[k]
    return
end

function update_sat_unsat(en_thermo::EnvironmentThermodynamics, state, k, ts)
    q_liq = TD.liquid_specific_humidity(ts)
    q_ice = TD.ice_specific_humidity(ts)
    aux_en = center_aux_environment(state)
    aux_en_sat = aux_en.sat
    aux_en_unsat = aux_en.unsat
    if TD.has_condensate(q_liq + q_ice)
        aux_en.cloud_fraction[k] = 1.0
        aux_en_sat.θ_dry[k] = TD.dry_pottemp(ts)
        aux_en_sat.θ_liq_ice[k] = TD.liquid_ice_pottemp(ts)
        aux_en_sat.T[k] = TD.air_temperature(ts)
        aux_en_sat.q_tot[k] = TD.total_specific_humidity(ts)
        aux_en_sat.q_vap[k] = TD.vapor_specific_humidity(ts)
    else
        aux_en.cloud_fraction[k] = 0.0
        aux_en_unsat.θ_dry[k] = TD.dry_pottemp(ts)
        aux_en_unsat.θ_virt[k] = TD.virtual_pottemp(ts)
        aux_en_unsat.q_tot[k] = TD.total_specific_humidity(ts)
    end
    return
end

function sgs_mean(en_thermo::EnvironmentThermodynamics, grid, state, en, precip, dt, param_set)

    p0_c = center_ref_state(state).p0
    ρ0_c = center_ref_state(state).ρ0
    aux_en = center_aux_environment(state)
    prog_pr = center_prog_precipitation(state)

    @inbounds for k in real_center_indices(grid)
        # condensation
        q_tot_en = aux_en.q_tot[k]
        ts = thermo_state_pθq(param_set, p0_c[k], aux_en.θ_liq_ice[k], q_tot_en)
        # autoconversion and accretion
        mph = precipitation_formation(
            param_set,
            precip.precipitation_model,
            prog_pr.q_rai[k],
            prog_pr.q_sno[k],
            aux_en.area[k],
            ρ0_c[k],
            dt,
            ts,
        )
        update_sat_unsat(en_thermo, state, k, ts)
        update_env_precip_tendencies(
            state,
            k,
            mph.qt_tendency,
            mph.θ_liq_ice_tendency,
            mph.qr_tendency,
            mph.qs_tendency,
        )
    end
    return
end

function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, precip, dt, param_set)
    # TODO: double check this python-> julia translation
    # a, w = np.polynomial.hermite.hermgauss(en_thermo.quadrature_order)
    a, w = FastGaussQuadrature.gausshermite(en_thermo.quadrature_order)
    p0_c = center_ref_state(state).p0
    ρ0_c = center_ref_state(state).ρ0
    aux_en = center_aux_environment(state)
    prog_pr = center_prog_precipitation(state)
    aux_en_unsat = aux_en.unsat
    aux_en_sat = aux_en.sat

    #TODO - remember you output source terms multipierd by dt (bec. of instanteneous autoconcv)
    #TODO - add tendencies for gm H, QT and QR due to rain
    #TODO - if we start using eos_smpl for the updrafts calculations
    #       we can get rid of the two categories for outer and inner quad. points

    abscissas = a
    weights = w
    # arrays for storing quadarature points and ints for labeling items in the arrays
    # a python dict would be nicer, but its 30% slower than this (for python 2.7. It might not be the case for python 3)
    env_len = 8
    src_len = 8

    sqpi_inv = 1.0 / sqrt(π)
    sqrt2 = sqrt(2.0)

    epsilon = 10e-14 #np.finfo(np.float).eps

    # initialize the quadrature points and their labels
    inner_env = zeros(env_len)
    outer_env = zeros(env_len)
    inner_src = zeros(src_len)
    outer_src = zeros(src_len)
    i_ql, i_qi, i_T, i_cf, i_qt_sat, i_qt_unsat, i_T_sat, i_T_unsat = 1:env_len
    i_SH_qt, i_Sqt_H, i_SH_H, i_Sqt_qt, i_Sqt, i_SH, i_Sqr, i_Sqs = 1:src_len

    @inbounds for k in real_center_indices(grid)
        if (
            aux_en.QTvar[k] > epsilon &&
            aux_en.Hvar[k] > epsilon &&
            abs(aux_en.HQTcov[k]) > epsilon &&
            aux_en.q_tot[k] > epsilon &&
            sqrt(aux_en.QTvar[k]) < aux_en.q_tot[k]
        )

            if en_thermo.quadrature_type == "log-normal"
                # Lognormal parameters (mu, sd) from mean and variance
                sd_q = sqrt(log(aux_en.QTvar[k] / aux_en.q_tot[k] / aux_en.q_tot[k] + 1.0))
                sd_h = sqrt(log(aux_en.Hvar[k] / aux_en.θ_liq_ice[k] / aux_en.θ_liq_ice[k] + 1.0))
                # Enforce Schwarz"s inequality
                corr = max(min(aux_en.HQTcov[k] / sqrt(aux_en.Hvar[k] * aux_en.QTvar[k]), 1.0), -1.0)
                sd2_hq =
                    log(corr * sqrt(aux_en.Hvar[k] * aux_en.QTvar[k]) / aux_en.θ_liq_ice[k] / aux_en.q_tot[k] + 1.0)
                sd_cond_h_q = sqrt(max(sd_h * sd_h - sd2_hq * sd2_hq / sd_q / sd_q, 0.0))
                mu_q =
                    log(aux_en.q_tot[k] * aux_en.q_tot[k] / sqrt(aux_en.q_tot[k] * aux_en.q_tot[k] + aux_en.QTvar[k]))
                mu_h = log(
                    aux_en.θ_liq_ice[k] * aux_en.θ_liq_ice[k] /
                    sqrt(aux_en.θ_liq_ice[k] * aux_en.θ_liq_ice[k] + aux_en.Hvar[k]),
                )
            else
                sd_q = sqrt(aux_en.QTvar[k])
                sd_h = sqrt(aux_en.Hvar[k])
                corr = max(min(aux_en.HQTcov[k] / max(sd_h * sd_q, 1e-13), 1.0), -1.0)

                # limit sd_q to prevent negative qt_hat
                sd_q_lim = (1e-10 - aux_en.q_tot[k]) / (sqrt2 * abscissas[1])
                # walking backwards to assure your q_t will not be smaller than 1e-10
                # TODO - check
                # TODO - change 1e-13 and 1e-10 to some epislon
                sd_q = min(sd_q, sd_q_lim)
                qt_var = sd_q * sd_q
                sigma_h_star = sqrt(max(1.0 - corr * corr, 0.0)) * sd_h
            end

            # zero outer quadrature points
            for idx in 1:env_len
                outer_env[idx] = 0.0
            end
            for idx in 1:src_len
                outer_src[idx] = 0.0
            end

            for m_q in 1:(en_thermo.quadrature_order)
                if en_thermo.quadrature_type == "log-normal"
                    qt_hat = exp(mu_q + sqrt2 * sd_q * abscissas[m_q])
                    mu_h_star = mu_h + sd2_hq / sd_q / sd_q * (log(qt_hat) - mu_q)
                else
                    qt_hat = aux_en.q_tot[k] + sqrt2 * sd_q * abscissas[m_q]
                    mu_h_star = aux_en.θ_liq_ice[k] + sqrt2 * corr * sd_h * abscissas[m_q]
                end

                # zero inner quadrature points
                for idx in 1:env_len
                    inner_env[idx] = 0.0
                end
                for idx in 1:src_len
                    inner_src[idx] = 0.0
                end

                for m_h in 1:(en_thermo.quadrature_order)
                    if en_thermo.quadrature_type == "log-normal"
                        h_hat = exp(mu_h_star + sqrt2 * sd_cond_h_q * abscissas[m_h])
                    else
                        h_hat = sqrt2 * sigma_h_star * abscissas[m_h] + mu_h_star
                    end

                    # condensation
                    ts = thermo_state_pθq(param_set, p0_c[k], h_hat, qt_hat)
                    q_liq_en = TD.liquid_specific_humidity(ts)
                    q_ice_en = TD.ice_specific_humidity(ts)
                    T = TD.air_temperature(ts)
                    # autoconversion and accretion
                    mph = precipitation_formation(
                        param_set,
                        precip.precipitation_model,
                        prog_pr.q_rai[k],
                        prog_pr.q_sno[k],
                        aux_en.area[k],
                        ρ0_c[k],
                        dt,
                        ts,
                    )

                    # environmental variables
                    inner_env[i_ql] += q_liq_en * weights[m_h] * sqpi_inv
                    inner_env[i_qi] += q_ice_en * weights[m_h] * sqpi_inv
                    inner_env[i_T] += T * weights[m_h] * sqpi_inv
                    # cloudy/dry categories for buoyancy in TKE
                    if TD.has_condensate(q_liq_en + q_ice_en)
                        inner_env[i_cf] += weights[m_h] * sqpi_inv
                        inner_env[i_qt_sat] += qt_hat * weights[m_h] * sqpi_inv
                        inner_env[i_T_sat] += T * weights[m_h] * sqpi_inv
                    else
                        inner_env[i_qt_unsat] += qt_hat * weights[m_h] * sqpi_inv
                        inner_env[i_T_unsat] += T * weights[m_h] * sqpi_inv
                    end
                    # products for variance and covariance source terms
                    inner_src[i_Sqt] += mph.qt_tendency * weights[m_h] * sqpi_inv
                    inner_src[i_Sqr] += mph.qr_tendency * weights[m_h] * sqpi_inv
                    inner_src[i_Sqs] += mph.qs_tendency * weights[m_h] * sqpi_inv
                    inner_src[i_SH] += mph.θ_liq_ice_tendency * weights[m_h] * sqpi_inv
                    inner_src[i_Sqt_H] += mph.qt_tendency * h_hat * weights[m_h] * sqpi_inv
                    inner_src[i_Sqt_qt] += mph.qt_tendency * qt_hat * weights[m_h] * sqpi_inv
                    inner_src[i_SH_H] += mph.θ_liq_ice_tendency * h_hat * weights[m_h] * sqpi_inv
                    inner_src[i_SH_qt] += mph.θ_liq_ice_tendency * qt_hat * weights[m_h] * sqpi_inv
                end

                for idx in 1:env_len
                    outer_env[idx] += inner_env[idx] * weights[m_q] * sqpi_inv
                end
                for idx in 1:src_len
                    outer_src[idx] += inner_src[idx] * weights[m_q] * sqpi_inv
                end
            end

            # update environmental variables
            update_env_precip_tendencies(
                state,
                k,
                outer_src[i_Sqt],
                outer_src[i_SH],
                outer_src[i_Sqr],
                outer_src[i_Sqs],
            )

            # update cloudy/dry variables for buoyancy in TKE
            aux_en.cloud_fraction[k] = outer_env[i_cf]
            aux_en_unsat.q_tot[k] = outer_env[i_qt_unsat]
            # TD.jl cannot compute θ_unsat when T=0
            aux_en_unsat.θ_dry[k] = if outer_env[i_T_unsat] > 0
                ts_unsat = TD.PhaseEquil_pTq(param_set, p0_c[k], outer_env[i_T_unsat], aux_en_unsat.q_tot[k])
                TD.dry_pottemp(ts_unsat)
            else
                0
            end

            aux_en_sat.T[k] = outer_env[i_T_sat]
            aux_en_sat.q_vap[k] = outer_env[i_qt_sat] - outer_env[i_ql]
            aux_en_sat.q_tot[k] = outer_env[i_qt_sat]
            ts_sat = TD.PhaseEquil_pTq(param_set, p0_c[k], aux_en_sat.T[k], aux_en_sat.q_tot[k])
            aux_en_sat.θ_dry[k] = TD.dry_pottemp(ts_sat)
            aux_en_sat.θ_liq_ice[k] = TD.liquid_ice_pottemp(ts_sat)

            # update var/covar rain sources
            aux_en.Hvar_rain_dt[k] = outer_src[i_SH_H] - outer_src[i_SH] * aux_en.θ_liq_ice[k]
            aux_en.QTvar_rain_dt[k] = outer_src[i_Sqt_qt] - outer_src[i_Sqt] * aux_en.q_tot[k]
            aux_en.HQTcov_rain_dt[k] =
                outer_src[i_SH_qt] - outer_src[i_SH] * aux_en.q_tot[k] + outer_src[i_Sqt_H] -
                outer_src[i_Sqt] * aux_en.θ_liq_ice[k]

        else
            # if variance and covariance are zero do the same as in SA_mean
            ts = thermo_state_pθq(param_set, p0_c[k], aux_en.θ_liq_ice[k], aux_en.q_tot[k])
            mph = precipitation_formation(
                param_set,
                precip.precipitation_model,
                prog_pr.q_rai[k],
                prog_pr.q_sno[k],
                aux_en.area[k],
                ρ0_c[k],
                dt,
                ts,
            )
            update_env_precip_tendencies(
                state,
                k,
                mph.qt_tendency,
                mph.θ_liq_ice_tendency,
                mph.qr_tendency,
                mph.qs_tendency,
            )
            update_sat_unsat(en_thermo, state, k, ts)

            aux_en.Hvar_rain_dt[k] = 0.0
            aux_en.QTvar_rain_dt[k] = 0.0
            aux_en.HQTcov_rain_dt[k] = 0.0
        end
    end

    return
end

function microphysics(en_thermo::EnvironmentThermodynamics, grid, state, en, precip, dt, param_set)

    if en.EnvThermo_scheme == "mean"
        sgs_mean(en_thermo, grid, state, en, precip, dt, param_set)

    elseif en.EnvThermo_scheme == "quadrature"
        sgs_quadrature(en_thermo, grid, state, en, precip, dt, param_set)

    else
        error("EDMF_Environment: Unrecognized EnvThermo_scheme. Possible options: mean, quadrature")
    end

    return
end