text:
@enum InstructionOp::UInt8 begin
    INSTR_ADD
    INSTR_SUB
    INSTR_NEG
    INSTR_MUL
    INSTR_DIV
    INSTR_MULADD # muladd(x,y,z) = x * y + z
    INSTR_MULSUB # mulsub(x,y,z) = x * y - z
    INSTR_SUBMUL # submul(x,y,z) = z - x * y
    INSTR_SQR
    INSTR_POW
    INSTR_SIN
    INSTR_COS
end

function call_op(op::InstructionOp)
    if op == INSTR_ADD
        :add_fast
    elseif op == INSTR_SUB
        :sub_fast
    elseif op == INSTR_NEG
        :neg_fast
    elseif op == INSTR_MUL
        :mul_fast
    elseif op == INSTR_DIV
        :div_fast
    elseif op == INSTR_MULADD
        :muladd_fast
    elseif op == INSTR_MULSUB
        :mulsub_fast
    elseif op == INSTR_SUBMUL
        :submul_fast
    elseif op == INSTR_SQR
        :sqr_fast
    elseif op == INSTR_POW
        :pow_fast
    elseif op == INSTR_SIN
        :sin
    elseif op == INSTR_COS
        :cos
    else
        error(op)
        :UNKNOWN
    end
end

function taylor_call_op(op::InstructionOp)
    if op == INSTR_ADD
        :taylor_add
    elseif op == INSTR_SUB
        :taylor_sub
    elseif op == INSTR_NEG
        :tayor_neg
    elseif op == INSTR_MUL
        :taylor_mul
    elseif op == INSTR_DIV
        :taylor_div
    elseif op == INSTR_MULADD
        :taylor_muladd
    elseif op == INSTR_MULSUB
        :taylor_mulsub
    elseif op == INSTR_SUBMUL
        :taylor_submul
    elseif op == INSTR_SQR
        :taylor_sqr
    elseif op == INSTR_POW
        :taylor_pow
    elseif op == INSTR_SIN
        :taylor_sin
    elseif op == INSTR_COS
        :taylor_cos
    else
        error(op)
        :UNKNOWN
    end
end

function is_univariate(op::InstructionOp)
    op == INSTR_SQR || op == INSTR_NEG
end

function arity(op::InstructionOp)
    if op == INSTR_ADD
        2
    elseif op == INSTR_SUB
        2
    elseif op == INSTR_NEG
        1
    elseif op == INSTR_MUL
        2
    elseif op == INSTR_DIV
        2
    elseif op == INSTR_MULADD
        3
    elseif op == INSTR_MULSUB
        3
    elseif op == INSTR_SUBMUL
        3
    elseif op == INSTR_SQR
        1
    elseif op == INSTR_POW
        2
    elseif op == INSTR_SIN
        1
    elseif op == INSTR_COS
        1
    else
        error(op)
        :UNKNOWN
    end
end
@inline neg_fast(x) = Base.FastMath.sub_fast(x)
@inline add_fast(x, y) = Base.FastMath.add_fast(x, y)
@inline sub_fast(x, y) = Base.FastMath.sub_fast(x, y)
@inline mul_fast(x, y) = Base.FastMath.mul_fast(x, y)
@inline div_fast(x, y) = Base.FastMath.div_fast(x, y)
@inline pow_fast(x, p::Integer) = Base.power_by_squaring(x, p)
@inline muladd_fast(x, y, z) = @fastmath x * y + z
@inline mulsub_fast(x, y, z) = @fastmath x * y - z
@inline submul_fast(x, y, z) = @fastmath z - x * y
@inline sqr_fast(x) = sqr(x)
@inline sqr(x) = @fastmath x * x
@inline function sqr(z::Complex)
    x, y = reim(z)
    @fastmath Complex((x + y) * (x - y), (x + x) * y)
end
"""
    inv_not_zero(x)

Invert x unless it is 0, then return 0.
"""
@inline inv_not_zero(x) = ifelse(is_zero(x), x, @fastmath inv(x))
@inline inv_fast(x) = @fastmath inv(x)

struct InstructionRef
    i::Int
end
Base.show(io::IO, i::InstructionRef) = print(io, "ι_", i.i)

const InstructionArg = Union{Nothing,Number,Symbol,InstructionRef}

struct Instruction
    op::InstructionOp
    args::Tuple{InstructionArg,InstructionArg,InstructionArg}
end

Instruction(op::InstructionOp, a::InstructionArg) = Instruction(op, (a, nothing, nothing))
Instruction(op, a::InstructionArg, b::InstructionArg) = Instruction(op, (a, b, nothing))
Instruction(op, a::InstructionArg, b::InstructionArg, c::InstructionArg) =
    Instruction(op, (a, b, c))
Base.getindex(I::Instruction, i) = I.args[i]

struct InstructionList
    instructions::Vector{Instruction}
    lifetime::Vector{Union{Nothing,Int}}
end

InstructionList() = InstructionList(Instruction[], Union{Nothing,Int}[])
Base.length(v::InstructionList) = length(v.instructions)
function Base.iterate(I::InstructionList, state = nothing)
    if isnothing(state)
        next_state = iterate(I.instructions)
    else
        next_state = iterate(I.instructions, state)
    end
    if isnothing(next_state)
        return nothing
    end
    next, state = next_state
    (InstructionRef(state - 1), next), state
end

function instruction_list(
    exprs::Vector{Expression};
    subexpressions::Dict{Expression,Expression} = Dict{Expression,Expression}(),
    perform_cse::Bool = true,
)
    list = InstructionList()
    if perform_cse
        exprs, CSE = cse(exprs)
        merge!(CSE, subexpressions)
    else
        CSE = subexpressions
    end
    symb_cse = Dict{Symbol,Expression}()
    for (k, v) in CSE
        push!(symb_cse, Symbol(to_string(k)) => v)
    end
    pse = Dict{Symbol,InstructionRef}()
    assignments = map(ex -> add_instructions!(list, ex, symb_cse, pse), exprs)

    # makesure to not overwrite assignments
    for iref in assignments
        update_lifetime!(list, iref, typemax(Int64))
    end
    list, assignments
end

function Base.show(io::IO, list::InstructionList)
    for (id, instr) in list
        println(
            io,
            :($id = $(Expr(:call, Symbol(instr.op), filter(!isnothing, instr.args)...))),
        )
    end
end

function update_lifetime!(list::InstructionList, x::InstructionRef, n::Int)
    list.lifetime[x.i] = n
    nothing
end
update_lifetime!(list::InstructionList, x::InstructionArg, n::Int) = nothing

function Base.push!(list::InstructionList, x::Instruction)
    push!(list.instructions, x)
    n = length(list)
    # by default does an element have infinite lifetime
    push!(list.lifetime, typemax(Int64))

    (a, b, c) = x.args
    update_lifetime!(list, a, n)
    update_lifetime!(list, b, n)
    update_lifetime!(list, c, n)

    InstructionRef(n)
end

# Helpers for InstructionOps
function add!(list::InstructionList, @nospecialize(a), @nospecialize(b))
    if a === nothing
        b === nothing ? nothing : b
    else
        b === nothing ? a : push!(list, Instruction(INSTR_ADD, a, b))
    end
end
function neg!(list::InstructionList, @nospecialize(a))
    a === nothing ? nothing : push!(list, Instruction(INSTR_NEG, a))
end
function sub!(list::InstructionList, @nospecialize(a), @nospecialize(b))
    if a === nothing
        b === nothing ? nothing : push!(list, Instruction(INSTR_NEG, b))
    else
        b === nothing ? a : push!(list, Instruction(INSTR_SUB, a, b))
    end
end

function mul!(list::InstructionList, @nospecialize(a), @nospecialize(b))
    (a === nothing || b === nothing) && return nothing
    if is_one(a)
        is_one(b) ? a : b
    elseif is_minus_one(a)
        push!(list, Instruction(INSTR_NEG, b))
    elseif is_minus_one(b)
        push!(list, Instruction(INSTR_NEG, a))
    elseif a isa Integer && a == 2
        push!(list, Instruction(INSTR_ADD, b, b))
    else
        is_one(b) ? a : push!(list, Instruction(INSTR_MUL, a, b))
    end
end

function muladd!(
    list::InstructionList,
    @nospecialize(a),
    @nospecialize(b),
    @nospecialize(c),
)
    (a === nothing || b === nothing) && return c
    c === nothing && return mul!(list, a, b)
    is_one(a) && return add!(list, b, c)
    is_one(b) && return add!(list, a, c)
    push!(list, Instruction(INSTR_MULADD, a, b, c))
end
function mulsub!(
    list::InstructionList,
    @nospecialize(a),
    @nospecialize(b),
    @nospecialize(c),
)
    (a === nothing || b === nothing) && return neg!(list, c)
    c === nothing && return mul!(list, a, b)
    is_one(a) && return sub!(list, b, c)
    is_one(b) && return sub!(list, a, c)
    push!(list, Instruction(INSTR_MULSUB, a, b, c))
end
function submul!(
    list::InstructionList,
    @nospecialize(a),
    @nospecialize(b),
    @nospecialize(c),
)
    # c - a * b
    (a === nothing || b === nothing) && return c
    c === nothing && return neg!(list, mul!(list, a, b))
    is_one(a) && return sub!(list, c, b)
    is_one(b) && return sub!(list, c, a)
    push!(list, Instruction(INSTR_SUBMUL, a, b, c))
end

function div!(list::InstructionList, @nospecialize(a), @nospecialize(b))
    a === nothing && return nothing
    is_one(b) ? a : push!(list, Instruction(INSTR_DIV, a, b))
end

function sqr!(list, @nospecialize(a))
    a === nothing && return nothing
    push!(list, Instruction(INSTR_SQR, a))
end
function pow!(list, @nospecialize(a), b::Integer)
    @assert b ≥ 0
    is_one(b) && return a
    b == 2 && return sqr!(list, a)
    is_zero(b) && return Int32(1)
    is_one(a) && return a
    a === nothing ? nothing : push!(list, Instruction(INSTR_POW, a, b))
end

function add_instructions!(
    list::InstructionList,
    ex::Basic,
    # common sub expression
    cse::Dict{Symbol,Expression},
    # processed sub expression
    # mapping of the cse var with the instruction used
    pse::Dict{Symbol,InstructionRef},
)
    t = class(ex)
    if t == :Symbol
        s = Symbol(to_string(ex))
        if haskey(pse, s)
            return pse[s]
        elseif haskey(cse, s)
            val = add_instructions!(list, cse[s], cse, pse)
            if val isa InstructionRef
                push!(pse, s => val)
            end
            val
        else
            s
        end
    elseif t == :Pow
        x, k_expr = args(ex)
        k = convert(Int32, k_expr)
        xinstr = add_instructions!(list, x, cse, pse)
        if k == -1
            div!(list, Int32(1), xinstr)
        elseif k < -1
            div!(list, Int32(1), pow!(list, xinstr, -k))
        else
            pow!(list, xinstr, k)
        end
    elseif t == :Mul
        _add_mul_instructions!(list, args(ex), cse, pse)
    elseif t == :Add
        # For sums we do work to find muladd, mulsub and negmuladd instructions
        vec = args(ex)
        neg = false
        op_arg = nothing
        for i = 1:length(vec)
            vᵢ = vec[i]
            if class(vᵢ) == :Mul
                vᵢ_args = args(vᵢ)
                nᵢ = length(vᵢ_args)
                if is_minus_one(vᵢ_args[1])
                    x = add_instructions!(list, vᵢ_args[2], cse, pse)
                    if length(vᵢ_args) > 2
                        y = _add_mul_instructions!(list, view(vᵢ_args, 3:nᵢ), cse, pse)
                        if isnothing(op_arg)
                            op_arg = mul!(list, x, y)
                            neg = true
                        elseif neg
                            op_arg = muladd!(list, x, y, op_arg)
                            neg = true
                        else
                            op_arg = submul!(list, x, y, op_arg)
                            neg = false
                        end
                    else
                        y = add_instructions!(list, vᵢ_args[2], cse, pse)
                        if isnothing(op_arg)
                            op_arg = add_instructions!(list, vᵢ_args[2], cse, pse)
                            neg = true
                        elseif neg
                            op_arg = add!(list, y, op_arg)
                            neg = true
                        else
                            op_arg = sub!(list, op_arg, y)
                            neg = false
                        end
                    end

                elseif isnothing(op_arg)
                    op_arg = add_instructions!(list, vᵢ, cse, pse)
                else
                    x = add_instructions!(list, vᵢ_args[1], cse, pse)
                    y = _add_mul_instructions!(list, view(vᵢ_args, 2:nᵢ), cse, pse)
                    if neg
                        op_arg = mulsub!(list, x, y, op_arg)
                        neg = false
                    else
                        op_arg = muladd!(list, x, y, op_arg)
                    end
                end
            elseif isnothing(op_arg)
                op_arg = add_instructions!(list, vᵢ, cse, pse)
            elseif neg
                op_arg = sub!(list, add_instructions!(list, vᵢ, cse, pse), op_arg)
                neg = false
            else
                op_arg = add!(list, add_instructions!(list, vᵢ, cse, pse), op_arg)
            end
        end
        neg ? neg!(list, op_arg) : op_arg
    elseif t == :Div
        x, y = ModelKit.args(ex)
        return div!(
            list,
            add_instructions!(list, x, cse, pse),
            add_instructions!(list, y, cse, pse),
        )
    elseif t == :Sin
        x = add_instructions!(list, ModelKit.args(ex)[1], cse, pse)
        return push!(list, Instruction(INSTR_SIN, x))
    elseif t == :Cos
        x = add_instructions!(list, ModelKit.args(ex)[1], cse, pse)
        return push!(list, Instruction(INSTR_COS, x))
    else
        return to_number(ex)
    end
end

function _add_mul_instructions!(list, vec, cse, pse)
    length(vec) == 1 && return add_instructions!(list, vec[1], cse, pse)
    op_arg = nothing
    is_neg = false
    divs = []
    for vec_arg in vec
        if class(vec_arg) == :Pow
            x, k_expr = args(vec_arg)
            k = convert(Int32, k_expr)
            xinstr = add_instructions!(list, x, cse, pse)
            if k == -1
                push!(divs, xinstr)
            elseif k < -1
                push!(divs, pow!(list, xinstr, -k))
            else
                if isnothing(op_arg)
                    op_arg = pow!(list, xinstr, k)
                else
                    op_arg = mul!(list, op_arg, pow!(list, xinstr, k))
                end
            end
        elseif is_minus_one(vec_arg)
            is_neg = !is_neg
        else
            vinstr = add_instructions!(list, vec_arg, cse, pse)
            if isnothing(op_arg)
                op_arg = vinstr
            else
                op_arg = mul!(list, op_arg, vinstr)
            end
        end
    end

    if !isempty(divs)
        denom = first(divs)
        for i = 2:length(divs)
            denom = mul!(list, denom, divs[i])
        end
        if isnothing(op_arg)
            op_arg = div!(list, Int32(1), denom)
        else
            op_arg = div!(list, op_arg, denom)
        end
    end
    is_neg ? neg!(list, op_arg) : op_arg
end

struct DiffMap
    D::Dict{Tuple{Any,Int},Any}
end
DiffMap() = DiffMap(Dict{Tuple{Any,Int},Any}())
#
Base.getindex(D::DiffMap, var, ∂i) = get(D.D, (var, ∂i), nothing)
Base.setindex!(D::DiffMap, val, var, ∂i) = D.D[(var, ∂i)] = val
Base.setindex!(D::DiffMap, ::Nothing, var, ∂i) = nothing

"""
    gradient(list::InstructionList, N::Int, diff_map)

Returns an `InstructionList` list computing the gradients all instructions.
"""
function gradient(list::InstructionList, assignments, variables::Vector{Symbol})
    D = DiffMap()
    for (∂i, var) in enumerate(variables)
        D[var, ∂i] = Int32(1)
    end
    N = length(variables)

    grad_list = InstructionList()
    list_ref_mapping = Dict{InstructionRef,InstructionRef}()
    for (instr_ref, instr) in list
        op = instr.op
        arg₁, arg₂, arg₃ = get.(Ref(list_ref_mapping), instr.args, instr.args)

        # handle pow first here since we change the original instruction
        if op == INSTR_POW
            base = pow!(grad_list, arg₁, arg₂ - 1)
            new_ref = mul!(grad_list, arg₁, base)
            for i = 1:N
                D[new_ref, i] = mul!(grad_list, arg₂, mul!(grad_list, base, D[arg₁, i]))
            end
        else
            new_ref = push!(grad_list, Instruction(op, (arg₁, arg₂, arg₃)))
        end
        list_ref_mapping[instr_ref] = new_ref

        if op == INSTR_ADD
            for i = 1:N
                D[new_ref, i] = add!(grad_list, D[arg₁, i], D[arg₂, i])
            end
        elseif op == INSTR_SUB
            for i = 1:N
                D[new_ref, i] = sub!(grad_list, D[arg₁, i], D[arg₂, i])
            end
        elseif op == INSTR_NEG
            for i = 1:N
                D[new_ref, i] = neg!(grad_list, D[arg₁, i])
            end
        elseif op == INSTR_MUL
            for i = 1:N
                D[new_ref, i] =
                    muladd!(grad_list, arg₁, D[arg₂, i], mul!(grad_list, D[arg₁, i], arg₂))
            end
        elseif op == INSTR_DIV
            denom = sqr!(grad_list, arg₂)
            for i = 1:N
                D[new_ref, i] = div!(
                    grad_list,
                    mulsub!(grad_list, D[arg₁, i], arg₂, mul!(grad_list, arg₁, D[arg₂, i])),
                    denom,
                )
            end
        elseif op == INSTR_MULADD
            for i = 1:N
                D[new_ref, i] = muladd!(
                    grad_list,
                    arg₁,
                    D[arg₂, i],
                    muladd!(grad_list, D[arg₁, i], arg₂, D[arg₃, i]),
                )
            end
        elseif op == INSTR_MULSUB
            for i = 1:N
                D[new_ref, i] = muladd!(
                    grad_list,
                    arg₁,
                    D[arg₂, i],
                    mulsub!(grad_list, D[arg₁, i], arg₂, D[arg₃, i]),
                )
            end
        elseif op == INSTR_SUBMUL
            for i = 1:N
                D[new_ref, i] = submul!(
                    grad_list,
                    arg₁,
                    D[arg₂, i],
                    submul!(grad_list, D[arg₁, i], arg₂, D[arg₃, i]),
                )
            end
        elseif op == INSTR_SQR
            for i = 1:N
                D[new_ref, i] = mul!(grad_list, mul!(grad_list, 2, arg₁), D[arg₁, i])
            end
        elseif op == INSTR_SIN
            for i = 1:N
                if !isnothing(D[arg₁, i])
                    D[new_ref, i] = mul!(
                        grad_list,
                        push!(grad_list, Instruction(INSTR_COS, arg₁)),
                        D[arg₁, i],
                    )
                end
            end
        elseif op == INSTR_COS
            for i = 1:N
                if !isnothing(D[arg₁, i])
                    D[new_ref, i] = neg!(
                        grad_list,
                        mul!(
                            grad_list,
                            push!(grad_list, Instruction(INSTR_SIN, arg₁)),
                            D[arg₁, i],
                        ),
                    )
                end
            end
        elseif op == INSTR_POW
            # handled above
        else
            error(op)
            :UNKNOWN
        end
    end

    new_assignments = [get(list_ref_mapping, instr, instr) for instr in assignments]
    grad_assignments = map([D[instr, ∂i] for instr in new_assignments, ∂i = 1:N]) do x
        isnothing(x) ? 0 : x
    end
    grad_list, new_assignments, grad_assignments
end

"""
    allocation_alias(list::InstructionList)

This performs an analysis of the recorded instructions to figure to decrease the number of
allocations necessary.

## Example

```julia
julia> ex = (a * b * (h - u) - c * d) * (h - u);

julia> list, _ = ModelKit.instruction_list([ex]);

julia> ModelKit.to_expr(list)
quote
    ι1 = sub_fast(h, u)
    ι2 = mul_fast(c, d)
    ι3 = mul_fast(b, ι1)
    ι4 = mulsub_fast(a, ι3, ι2)
    ι5 = mul_fast(ι1, ι4)
end

julia> ModelKit.to_expr(list; allocation_alias = ModelKit.allocation_alias(list))
quote
    ι1 = sub_fast(h, u)
    ι2 = mul_fast(c, d)
    ι3 = mul_fast(b, ι1)
    ι3 = mulsub_fast(a, ι3, ι2)
    ι3 = mul_fast(ι1, ι3)
end

```
"""
function allocation_alias(list::InstructionList, assignments::Set{Int})
    n = length(list.lifetime)
    lifetimes = map((i, x) -> i ∈ assignments ? typemax(Int) : x, 1:n, list.lifetime)
    alias = Int[]
    open_slots = Int[]
    max_slots = 0
    deaths = sort!([(k, i) for (i, k) in enumerate(lifetimes)])
    death_index = 1

    for (k, Δ) in enumerate(lifetimes)
        # check for deaths
        n₁ = length(open_slots)
        # we use here < k and not ≤ k to avoid the potential of accidental overwrites due
        # to shared memory
        while first(deaths[death_index]) < k
            push!(open_slots, alias[last(deaths[death_index])])
            death_index += 1
        end

        if isempty(open_slots)
            max_slots += 1
            push!(alias, max_slots)
        else
            push!(alias, pop!(open_slots))
        end

        n₂ = length(open_slots)
        if n₁ != n₂
            sort!(open_slots, alg = InsertionSort)
        end
    end

    alias
end

to_expr_arg(i::InstructionRef, allocation_alias) = Symbol("ι", allocation_alias[i.i])
to_expr_arg(i::InstructionRef, ::Nothing = nothing) = Symbol("ι", i.i)
to_expr_arg(x, _allocation_alias = nothing) = x

function to_expr(
    f::F,
    list::InstructionList;
    variables = Dict{Symbol,Union{Expr,Symbol}}(),
    assignements = Dict{Symbol,Vector{Expr}}(),
    allocation_alias = nothing,
) where {F}
    block = Expr(:block)
    exprs = block.args
    for (instr_ref, el) in list
        id = to_expr_arg(instr_ref, allocation_alias)
        op, (arg1, arg2, arg3) = el.op, el.args
        a = to_expr_arg(get(variables, arg1, arg1), allocation_alias)
        if arg3 !== nothing && arg2 !== nothing
            c = to_expr_arg(get(variables, arg3, arg3), allocation_alias)
            b = to_expr_arg(get(variables, arg2, arg2), allocation_alias)
            push!(exprs, :($id = $(f(op, (a, b, c)))))
        elseif arg2 !== nothing
            if op == INSTR_POW && arg2 isa Int
                k::Int = arg2
                if k < 0
                    push!(exprs, :($id = inv_fast($(unroll_pow(get(variables, a, a), -k)))))
                else
                    push!(exprs, :($id = $(unroll_pow(get(variables, a, a), k))))
                end
            else
                b = to_expr_arg(get(variables, arg2, arg2), allocation_alias)
                push!(exprs, :($id = $(f(op, (a, b)))))
            end
        else
            push!(exprs, :($id = $(f(op, (a,)))))
        end

        if haskey(assignements, id)
            append!(exprs, assignements[id])
        end
    end
    block
end
to_expr(list::InstructionList; kwargs...) =
    to_expr(list; kwargs...) do op, args
        Expr(:call, call_op(op), args...)
    end

# HIGHER ORDER DIFFERENTIATION
# This follows Chapter 13 of Griewank & Walther - Evaluating derivatives

# Taylor
function untuple(varname, dx)
    :($(Expr(:tuple, (Symbol(varname, k) for k = 0:dx)...)) = $varname)
end
taylor_tuple(ids) = Expr(:tuple, [isnothing(id) ? :(zero(x0)) : id for id in ids]...)

function taylor_impl(f!, K::Int, dx::Int)
    D = DiffMap()
    list = InstructionList()
    for k = 0:dx
        D[:x, k] = Symbol(:x, k)
    end
    ids = to_expr_arg.(f!(list, D), nothing)
    quote
        $(untuple(:x, dx))
        $(to_expr(list; variables = Dict(:x => :x)))
        $(taylor_tuple(ids))
    end
end
function taylor_impl(f!, K::Int, dx::Int, dy::Int)
    D = DiffMap()
    list = InstructionList()
    for k = 0:dx
        D[:x, k] = Symbol(:x, k)
    end
    for k = 0:dy
        D[:y, k] = Symbol(:y, k)
    end
    ids = to_expr_arg.(f!(list, D), nothing)
    quote
        $(untuple(:x, dx))
        $(untuple(:y, dy))
        $(to_expr(list; variables = Dict(:x => :x, :y => :y)))
        $(taylor_tuple(ids))
    end
end
function taylor_impl(f!, K::Int, dx::Int, dy::Int, dz::Int)
    D = DiffMap()
    list = InstructionList()
    for k = 0:dx
        D[:x, k] = Symbol(:x, k)
    end
    for k = 0:dy
        D[:y, k] = Symbol(:y, k)
    end
    for k = 0:dz
        D[:z, k] = Symbol(:z, k)
    end
    ids = to_expr_arg.(f!(list, D), nothing)
    quote
        $(untuple(:x, dx))
        $(untuple(:y, dy))
        $(untuple(:z, dz))
        $(to_expr(list; variables = Dict(:x => :x, :y => :y, :z => :z)))
        $(taylor_tuple(ids))
    end
end

@inline make_ntuple(t::N) where {N<:Number} = (t,)
@inline make_ntuple(t::Tuple{A}) where {A} = t
@inline make_ntuple(t::Tuple{A,B}) where {A,B} = convert(NTuple{2,promote_type(A, B)}, t)
@inline make_ntuple(t::Tuple{A,B,C}) where {A,B,C} =
    convert(NTuple{3,promote_type(A, B, C)}, t)
@inline make_ntuple(t::Tuple{A,B,C,D}) where {A,B,C,D} =
    convert(NTuple{4,promote_type(A, B, C, D)}, t)
@inline make_ntuple(t::Tuple{A,B,C,D,E}) where {A,B,C,D,E} =
    convert(NTuple{5,promote_type(A, B, C, D, E)}, t)

for f in [:taylor_neg, :taylor_sin, :taylor_cos, :taylor_sqr]
    @eval @inline $f(::Val{K}, x::X) where {K,X} = $f(Val(K), make_ntuple(x))
end
for f in [:taylor_add, :taylor_sub, :taylor_mul, :taylor_div]
    @eval $f(::Val{K}, x::X, y::Y) where {K,X,Y} =
        $f(Val(K), make_ntuple(x), make_ntuple(y))
end
taylor_pow(::Val{K}, x::X, r::Integer) where {K,X} = taylor_pow(Val(K), make_ntuple(x), r)
for f in [:taylor_muladd, :taylor_mulsub, :taylor_submul]
    @eval @inline $f(::Val{K}, x::X, y::Y, z::Z) where {K,X,Y,Z} =
        $f(Val(K), make_ntuple(x), make_ntuple(y), make_ntuple(z))
end

@generated function taylor_neg(::Val{K}, x::NTuple{M}) where {K,M}
    taylor_impl(K, M - 1) do list, D
        [neg!(list, D[:x, k]) for k = 0:K]
    end
end

function taylor_sqr_impl(K::Int, dx::Int)
    taylor_impl(K, dx) do list, D
        map(0:K) do k
            # Compute ∑_{j=0}^k u_j U_{k - j} =
            #         2(∑_{j=0}^div(k - 1,2) u_j U_{k - j}) +
            #         iseven(k) * u_{k-1}^2
            k == 0 && return sqr!(list, D[:x, 0])
            w_k = nothing
            for j = 0:div(k - 1, 2)
                w_k = muladd!(list, D[:x, j], D[:x, k-j], w_k)
            end
            if iseven(k)
                muladd!(list, 2, w_k, sqr!(list, D[:x, div(k, 2)]))
            else
                add!(list, w_k, w_k)
            end
        end
    end
end
@generated function taylor_sqr(::Val{K}, x::NTuple{M}) where {K,M}
    taylor_sqr_impl(K, M - 1)
end

@generated function taylor_add(::Val{K}, x::NTuple{M}, y::NTuple{N}) where {K,M,N}
    taylor_impl(K, M - 1, N - 1) do list, D
        [add!(list, D[:x, k], D[:y, k]) for k = 0:K]
    end
end

@generated function taylor_sub(::Val{K}, x::NTuple{M}, y::NTuple{N}) where {K,M,N}
    taylor_impl(K, M - 1, N - 1) do list, D
        [sub!(list, D[:x, k], D[:y, k]) for k = 0:K]
    end
end

@generated function taylor_mul(::Val{K}, x::NTuple{M}, y::NTuple{N}) where {K,M,N}
    taylor_impl(K, M - 1, N - 1) do list, D
        map(0:K) do k
            c_k = nothing
            for j = 0:k
                c_k = muladd!(list, D[:x, j], D[:y, k-j], c_k)
            end
            c_k
        end
    end
end

@generated function taylor_div(::Val{K}, x::NTuple{M}, y::NTuple{N}) where {K,M,N}
    taylor_impl(K, M - 1, N - 1) do list, D
        ids = []
        for k = 0:K
            s = nothing
            for j = 0:(k-1)
                s = muladd!(list, ids[j+1], D[:y, k-j], s)
            end
            push!(ids, div!(list, sub!(list, D[:x, k], s), D[:y, 0]))
        end
        ids
    end
end

function taylor_pow_impl(K, dx)
    D = DiffMap()
    list = InstructionList()
    for k = 0:dx
        D[:x, k] = Symbol(:x, k)
    end
    w = Any[:w₀]
    for k = 1:K
        s = nothing
        for j = 1:k
            ũ_j = mul!(list, j, D[:x, j])
            w_kj = k == j ? :w₀ : w[k-j+1]
            s = muladd!(list, w_kj, ũ_j, s)
        end
        s = mul!(list, :r, s)

        t = nothing
        for j = 1:k-1
            u_kj = D[:x, k-j]
            w̃_j = mul!(list, j, w[j+1])
            t = muladd!(list, u_kj, w̃_j, t)
        end
        w̃_k = mul!(list, :u₀_inv, sub!(list, s, t))
        w_k = div!(list, w̃_k, k)
        push!(w, w_k)
    end
    quote
        $(untuple(:x, dx))
        iszero(x0) && return $(taylor_tuple([nothing for _ = 0:K]))
        w₀ = pow_fast(x0, r)
        $((K > 0 ? (:(u₀_inv = inv_fast(x0)),) : ())...)
        $(to_expr(list; variables = Dict(:x => :x, :u₀_inv => :u₀_inv, :w₀ => :w₀)))
        $(taylor_tuple(to_expr_arg.(w, nothing)))
    end
end
@generated function taylor_pow(::Val{K}, x::NTuple{M,T}, r::Integer) where {K,M,T}
    taylor_pow_impl(K, M - 1)
end

@generated function taylor_muladd(
    ::Val{K},
    x::NTuple{M},
    y::NTuple{N},
    z::NTuple{L},
) where {K,M,N,L}
    taylor_impl(K, M - 1, N - 1, L - 1) do list, D
        map(0:K) do k
            c_k = D[:z, k]
            for j = 0:k
                c_k = muladd!(list, D[:x, j], D[:y, k-j], c_k)
            end
            c_k
        end
    end
end
@generated function taylor_mulsub(
    ::Val{K},
    x::NTuple{M},
    y::NTuple{N},
    z::NTuple{L},
) where {K,M,N,L}
    taylor_impl(K, M - 1, N - 1, L - 1) do list, D
        map(0:K) do k
            c_k = D[:z, k]
            neg = false
            for j = 0:k
                if neg
                    c_k = muladd!(list, D[:x, j], D[:y, k-j], c_k)
                else
                    c_k = mulsub!(list, D[:x, j], D[:y, k-j], c_k)
                    neg = true
                end
            end
            neg ? c_k : neg!(list, c_k)
        end
    end
end

@generated function taylor_submul(
    ::Val{K},
    x::NTuple{M},
    y::NTuple{N},
    z::NTuple{L},
) where {K,M,N,L}
    taylor_impl(K, M - 1, N - 1, L - 1) do list, D
        map(0:K) do k
            c_k = D[:z, k]
            for j = 0:k
                c_k = submul!(list, D[:x, j], D[:y, k-j], c_k)
            end
            c_k
        end
    end
end

function taylor_sin_helper!(list, D, k)
    k == 0 && return :s₀
    s_k = nothing
    for j = 1:k
        s_k =
            muladd!(list, mul!(list, j, D[:x, j]), taylor_cos_helper!(list, D, k - j), s_k)
    end
    div!(list, s_k, k)
end
function taylor_cos_helper!(list, D, k)
    k == 0 && return :c₀
    c_k = nothing
    for j = 1:k
        c_k =
            submul!(list, mul!(list, j, D[:x, j]), taylor_sin_helper!(list, D, k - j), c_k)
    end
    div!(list, c_k, k)
end
function taylor_sin_impl(K, M)
    D = DiffMap()
    list = InstructionList()
    for k = 0:M-1
        D[:x, k] = Symbol(:x, k)
    end

    ids = Any[:s₀]
    for k = 1:K
        push!(ids, taylor_sin_helper!(list, D, k))
    end

    quote
        $(untuple(:x, M - 1))
        $(K > 0 ? :((s₀, c₀) = sincos(x0)) : :(s₀ = sin(x0)))
        $(to_expr(list; variables = Dict(:x => :x, :s₀ => :s₀, :c₀ => :c₀)))
        $(taylor_tuple(to_expr_arg.(ids, nothing)))
    end
end
function taylor_cos_impl(K, M)
    D = DiffMap()
    list = InstructionList()
    for k = 0:M-1
        D[:x, k] = Symbol(:x, k)
    end

    ids = Any[:c₀]
    for k = 1:K
        push!(ids, taylor_cos_helper!(list, D, k))
    end

    quote
        $(untuple(:x, M - 1))
        $(K > 0 ? :((s₀, c₀) = sincos(x0)) : :(c₀ = cos(x0)))
        $(to_expr(list; variables = Dict(:x => :x, :s₀ => :s₀, :c₀ => :c₀)))
        $(taylor_tuple(to_expr_arg.(ids, nothing)))
    end
end
@generated function taylor_sin(::Val{K}, x::NTuple{M}) where {K,M}
    taylor_sin_impl(K, M)
end
@generated function taylor_cos(::Val{K}, x::NTuple{M}) where {K,M}
    taylor_cos_impl(K, M)
end


text:
"""
    JopInterp(dom, rng)

Performs linear interpolation from dom::JetSpace to rng::JetSpace. It is often
used to reduce dimensionality in FWI.  JopInterp is ported from the CVX frequency
domain FWI tools in SeisSpace.

# Notes
* It is required that domain and range have the same dimensionality.
* It is required that domain and range both have more than 2 points per dimension
* It is required that the domain is coarser than the range.
* It is assumed that domain and range have the same boundary (e.g. in 2D: xmin,xmax,zmin,zmax).

# Examples:

## 2D
```julia
dom = JetSpace(Float32, 11, 11)
rng = JetSpace(Float32, 23, 23)
A = JopInterp(dom, rng)
y = A * rand(dom)
x = A' * rand(rng)
```
"""
function JopInterp(dom::JetAbstractSpace, rng::JetAbstractSpace)
    if length(size(dom)) != length(size(rng))
        error("dimensionality of domain ($(length(size(dom)))) and range ($(length(size(rng)))) do not agree!")
    end

    dsize = size(dom)
    rsize = size(rng)

    if rsize[1] < dsize[1] || rsize[2] < dsize[2]
        error("size of domain ($(size(dom))) smaller than range ($(size(rng))) !")
    end

    JopLn(dom=dom, rng=rng, df! = JopInterp_df!, df′! = JopInterp_df′!, s = (dom=dom, rng=rng))
end

export JopInterp

# coarse to dense
function JopInterp_apply_2d(domvec::AbstractArray, rngvec::AbstractArray, dom, rng, isadj::Bool)
    nzdom,nxdom = size(dom)
    nzrng,nxrng = size(rng)
    zmin,zmax = 0.0,1.0
    xmin,xmax = 0.0,1.0
    dzdom,dxdom = (zmax-zmin)/(nzdom-1),(xmax-xmin)/(nxdom-1)
    dzrng,dxrng = (zmax-zmin)/(nzrng-1),(xmax-xmin)/(nxrng-1)
    dxdzdom = dxdom * dzdom
    ϵ = 1000 * eps(eltype(domvec))
    wmin,wmax = 0-ϵ,1+ϵ

    if isadj
        domvec .= 0
    else
        rngvec .= 0
    end

    for kxrng = 1:nxrng
        for kzrng = 1:nzrng
            xx = xmin + dxrng * (kxrng-1)
            zz = zmin + dzrng * (kzrng-1)

            kxdom = clamp(floor(Int, (xx - xmin) / dxdom) + 1, 1, nxdom-1)
            kzdom = clamp(floor(Int, (zz - zmin) / dzdom) + 1, 1, nzdom-1)

            x1 = xmin + dxdom * (kxdom - 1 + 0)
            x2 = xmin + dxdom * (kxdom - 1 + 1)
            z1 = zmin + dzdom * (kzdom - 1 + 0)
            z2 = zmin + dzdom * (kzdom - 1 + 1)

            w1 = (x2 - xx) * (z2 - zz) / dxdzdom
            w2 = (xx - x1) * (z2 - zz) / dxdzdom
            w3 = (x2 - xx) * (zz - z1) / dxdzdom
            w4 = (xx - x1) * (zz - z1) / dxdzdom

            # Sanity checks on weights
            if abs(w1 + w2 + w3 + w4 - 1) > ϵ
                error("sum of weights minus 1 ($(w1+w2+w3+w4-1)) not close to zero")
            end

            if w1 < wmin || w1 > wmax || w2 < wmin || w2 > wmax ||
                w3 < wmin || w3 > wmax || w4 < wmin || w4 > wmax
                error("weights in interpolation ($(w1),$(w2),$(w3),$(w1)) exceed range [$(wmin),$(wmax)]")
            end

            if isadj
                domvec[kzdom+0,kxdom+0] += w1 * rngvec[kzrng,kxrng]
                domvec[kzdom+0,kxdom+1] += w2 * rngvec[kzrng,kxrng]
                domvec[kzdom+1,kxdom+0] += w3 * rngvec[kzrng,kxrng]
                domvec[kzdom+1,kxdom+1] += w4 * rngvec[kzrng,kxrng]
            else
                rngvec[kzrng,kxrng] += w1 * domvec[kzdom+0,kxdom+0]
                rngvec[kzrng,kxrng] += w2 * domvec[kzdom+0,kxdom+1]
                rngvec[kzrng,kxrng] += w3 * domvec[kzdom+1,kxdom+0]
                rngvec[kzrng,kxrng] += w4 * domvec[kzdom+1,kxdom+1]
            end
        end
    end
    nothing
end

JopInterp_df!(d::AbstractArray{T,2}, m::AbstractArray{T,2}; dom, rng, kwargs...) where {T} = begin JopInterp_apply_2d(m, d, dom, rng, false); d end
JopInterp_df′!(m::AbstractArray{T,2}, d::AbstractArray{T,2}; dom, rng, kwargs...) where {T} = begin  JopInterp_apply_2d(m, d, dom, rng, true); m end

ext: jl


text:
#!/usr/bin/env julia

export pe34

function pe34()
    intlength(n::Integer)=length(digits(n))
    fdsum(n::Integer)=sum(factorial.(digits(n)))
    max_digit=2
    while intlength(fdsum(10^max_digit-1))>max_digit
        max_digit +=1
    end
    filter(n->fdsum(n)==n,10:(10^max_digit-1))|>sum
end

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe34()|>print
end

max_stars_repo_path: 34. Digit factorials/34.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_path: 34. Digit factorials/34.jl

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_path: 34. Digit factorials/34.jl

max_forks_repo_name: miRoox/ProjectEuler


text:
# using Revise

using Test
using RoME
using DistributedFactorGraphs
using TensorCast
##

@testset "sanity check on Hex example" begin

##

fg = generateCanonicalFG_Hexagonal(graphinit=false)

getSolverParams(fg).drawtree = false
getSolverParams(fg).showtree = false
getSolverParams(fg).graphinit = false
getSolverParams(fg).treeinit = true
getSolverParams(fg).downsolve = true
getSolverParams(fg).multiproc = false
getSolverParams(fg).async = false

# @test_skip false
getSolverParams(fg).useMsgLikelihoods=true
# @error "Restore useMsgLikelihoods=true"

# direct solve would be
tree = solveTree!(fg)
# tree = solveTree!(fg)

@error "Hex init test degraded quality during useMsgLikelihoods refactor, must restore to 55"
M = getManifold(Pose2())

@cast cs[j,i] := getCoordinates.(Pose2, getPoints(getBelief(fg, :x0)))[i][j]
@test 35 < sum(-3.0 .< cs[1,:] .< 3.0)
@test 35 < sum(-3.0 .< cs[2,:] .< 3.0)
@test 35 < sum(-0.3 .< cs[3,:] .< 0.3)

@cast cs[j,i] := getCoordinates.(Pose2, getPoints(getBelief(fg, :x1)))[i][j]
@test 35 < sum(7.0 .< cs[1,:] .< 13.0)
@test 35 < sum(-3.0 .< cs[2,:] .< 3.0)
@test 35 < sum(0.7 .< cs[3,:] .< 1.3)

@cast cs[j,i] := getCoordinates.(Pose2, getPoints(getBelief(fg, :x2)))[i][j]
@test 35 < sum(12.0 .< cs[1,:] .< 18.0)
@test 35 < sum(6.0 .< cs[2,:] .< 11.0)
@test 35 < sum(1.8 .< cs[3,:] .< 2.4)

@cast cs[j,i] := getCoordinates.(Pose2, getPoints(getBelief(fg, :x3)))[i][j]
@test 35 < sum(7.0 .< cs[1,:] .< 13.0)
@test 35 < sum(15.0 .< cs[2,:] .< 20.0)
# @test 35 < sum(-0.3 .< cs[3,:] .< 0.3)
μ = mean(M, getPoints(getBelief(fg, :x3)))
@test isapprox(μ.parts[1], [11; 17.5], atol=3.0)
@test isapprox(SpecialOrthogonal(2), μ.parts[2], [-1 0; 0 -1], atol=0.5)
Σ = cov(getManifold(Pose2()), getPoints(getBelief(fg, :x3)))
@test all(diag(Σ) .< [5,5,1].^2) #TODO smaller

@cast cs[j,i] := getCoordinates.(Pose2, getPoints(getBelief(fg, :x4)))[i][j]
@test 35 < sum(-5.0 .< cs[1,:] .< 5.0)
@test 35 < sum(13.0 .< cs[2,:] .< 22.0)
@test 35 < sum(-2.8 .< cs[3,:] .< -1.5)

@cast cs[j,i] := getCoordinates.(Pose2, getPoints(getBelief(fg, :x5)))[i][j]
@test 35 < sum(-8.0 .< cs[1,:] .< -2.0)
@test 35 < sum(6.0 .< cs[2,:] .< 11.0)
@test 35 < sum(-1.3 .< cs[3,:] .< -0.7)

@cast cs[j,i] := getCoordinates.(Pose2, getPoints(getBelief(fg, :x6)))[i][j]
@test 35 < sum(-3.0 .< cs[1,:] .< 3.0)
@test 35 < sum(-3.0 .< cs[2,:] .< 3.0)
@test 35 < sum(-0.3 .< cs[3,:] .< 0.3)

@cast cs[j,i] := getCoordinates.(Point2, getPoints(getBelief(fg, :l1)))[i][j]
@test 35 < sum(17.0 .< cs[1,:] .< 23.0)
@test 35 < sum(-5.0 .< cs[2,:] .< 5.0)

##

0
end # testset

#
# #  Do some plotting
# using RoMEPlotting
# #
# Gadfly.set_default_plot_size(35cm,20cm)
# pl = drawPosesLandms(fg, drawhist=true)

## OR manually do clique by clique
#
# # solve by hand, one cliq at a time
# tree = wipeBuildNewTree!(fg)
# # drawTree(tree,show=true)
#
#
# # solve the first cliq
# smt, hist = solveCliq!(fg, tree, :x0, recordcliq=true)
#
# # solve second clq
# smt, hist = solveCliq!(fg, tree, :x2, cliqHistories=hist, recordcliq=true)
# # resetCliqSolve!(fg, tree, :x2)
#
# # solve forth clq
# smt, hist = solveCliq!(fg, tree, :x4, cliqHistories=hist, recordcliq=true)
# # solve forth clq
# smt, hist = solveCliq!(fg, tree, :x6, cliqHistories=hist, recordcliq=true)
# # solve forth clq
# smt, hist = solveCliq!(fg, tree, :x3, cliqHistories=hist, recordcliq=true)


text:
module TestLabeledArray

using FactCheck
using MultidimensionalTables

facts("LabeledArray tests") do
  context("constructor tests") do
    col1 = nalift(rand(10, 50))
    col2 = nalift(reshape(1:500, 10, 50))
    col3 = nalift(reshape(map(i->string("sym_",i), 1:500), 10, 50))
    col4 = nalift(hcat(rand(10,30), fill(:testsym, 10, 20)))
    axis1c1 = DictArray(k1=nalift(collect(101:110)), k2=nalift(collect(201:210)))
    axis1c2 = nalift(map(i->string("a_",i), 1:50))
    d = DictArray(c1=col1, c2=col2, c3=col3, c4=col4)
    larr1 = LabeledArray(d, (axis1c1, axis1c2))
    larr2 = LabeledArray(d, axis1=axis1c1, axis2=axis1c2)
    @fact larr1 --> larr2
    @fact convert(LabeledArray, LabeledArray([LDict(:a=>Nullable(3),:b=>Nullable{Int}()), LDict(:a=>Nullable(5),:b=>Nullable(3))],axis1=darr(k=[:x,:y]))) --> @larr(a=[3,5],b=[NA,3],axis1[darr(k=[:x,:y])])
    @fact convert(LabeledArray, LabeledArray(nalift([LDict(:a=>Nullable(3),:b=>Nullable{Int}()), LDict(:a=>Nullable(5),:b=>Nullable(3))]),axis1=darr(k=[:x,:y]))) --> @larr(a=[3,5],b=[NA,3],axis1[darr(k=[:x,:y])])
    @fact convert(LabeledArray, [LDict(:a=>Nullable(3),:b=>Nullable{Int}()), LDict(:a=>Nullable(5),:b=>Nullable(3))]) --> @larr(a=[3,5],b=[NA,3])
    @fact convert(LabeledArray, nalift([LDict(:a=>Nullable(3),:b=>Nullable{Int}()), LDict(:a=>Nullable(5),:b=>Nullable(3))])) --> @larr(a=[3,5],b=[NA,3])
    @fact @larr(a=[1,2,3],b=1,axis1[k=[:A,:B,:C]]) --> larr(a=[1,2,3],b=[1,1,1],axis1=darr(k=[:A,:B,:C]))
    @fact larr(@larr(a=[1,2,3],b=1,axis1[k=[:A,:B,:C]]), c=[:x,:y,:z], :d=>[1,2,3], axis1=darr(k=["X", "Y", "Z"])) --> larr(axis1=darr(k=["X","Y","Z"]), a=[1,2,3],b=[1,1,1],d=[1,2,3],c=[:x,:y,:z])
    @fact @larr(@larr(a=[1 NA;3 4;NA NA],:b=>[1.0 1.5;:sym 'a';"X" "Y"],c=1,axis1[:U,NA,:W],axis2[r=['m','n']]), c=[NA NA;3 4;5 6], :d=>:X, axis1[k=["g","h","i"]]) --> @larr(a=[1 NA;3 4;NA NA],b=[1.0 1.5;:sym 'a';"X" "Y"],c=[NA NA;3 4;5 6],d=reshape(fill(:X,6),3,2),axis2[r=['m','n']],axis1[k=["g","h","i"]])
  end
  context("array related method tests") do
    larr = LabeledArray(
      DictArray(:a=>nalift(reshape(1:20, 5, 4)),
                2=>nalift(fill(:sym, 5, 4)),
                :third=>nalift(reshape(map(i->string("str_",i), 1:20), 5, 4))),
          axis1=DictArray(k1=nalift(["a","a","b","b","b"]), k2=nalift(collect(101:105))),
          axis2=DictArray(r1=nalift([:alpha,:beta,:gamma,:delta])))
    @fact larr[3,2] --> LDict(:a=>Nullable(8), 2=>Nullable(:sym), :third=>Nullable("str_8"))
    @fact tbltool.getindexpair(larr,3,2).second.second --> LDict(:a=>Nullable(8), 2=>Nullable(:sym), :third=>Nullable("str_8"))
    @fact larr[1:3,2] --> LabeledArray(
                            DictArray(:a=>nalift([6,7,8]),
                                      2=>nalift([:sym,:sym,:sym]),
                                      :third=>nalift(["str_6","str_7","str_8"])),
                            axis1=DictArray(k1=nalift(["a","a","b"]), k2=nalift(collect(101:103))))
    @fact larr[2:2,2:4] --> LabeledArray(
                            DictArray(:a=>nalift([7 12 17]),
                                      2=>nalift([:sym :sym :sym]),
                                      :third=>nalift(["str_7" "str_12" "str_17"])),
                            axis1=DictArray(k1=nalift(["a"]), k2=nalift([102])),
                            axis2=DictArray(r1=nalift([:beta,:gamma,:delta])))
    @fact size(larr) --> (5,4)
    #@fact eltype(larr) --> LDict{Any}
    @fact pick(larr,:a)[1,1].value --> 1
    @fact pick(larr,(:a,))[1][1,1].value --> 1
    @fact endof(larr) --> length(larr)
    @fact transpose(larr) --> LabeledArray(DictArray(mapvalues(transpose, larr.data.data)),
                                           axis1=larr.axes[2],
                                           axis2=larr.axes[1])
    @fact permutedims(larr, (2,1)) --> transpose(larr)
    @fact LabeledArray(DictArray(a=nalift([1 2 3])), axes2=nalift([100,101,102])) ==
          LabeledArray(DictArray(a=nalift([1 2 3])), axes2=nalift([100,101,102])) --> true
    @fact copy(larr) !== larr --> true
    @fact size(cat(1, larr, larr)) --> ntuple(n->n==1 ? size(larr,1)*2 : size(larr,n), ndims(larr))
    @fact size(cat(2, larr, larr)) --> ntuple(n->n==2 ? size(larr,2)*2 : size(larr,n), ndims(larr))
    @fact size(repeat(larr, inner=[5,2], outer=[3,4])) --> (size(larr,1)*5*3, size(larr,2)*2*4)
    @fact @larr(a=[1 2 3;4 5 6], 'x'=[:a :b :c;:x :y :z])[:a] --> nalift([1 2 3;4 5 6])
    @fact @larr(a=[1 2 3;4 5 6], b=[:a :b :c;:x :y :z])[:a,:b] --> Any[nalift([1 2 3;4 5 6]), nalift([:a :b :c;:x :y :z])]
    @fact @larr(a=[1 2 3;4 5 6], b=[:a :b :c;:x :y :z])[[:a,:b]] --> @darr(a=[1 2 3;4 5 6], b=[:a :b :c;:x :y :z])
    @fact @larr(a=[1 2 3;4 5 6], b=[:a :b :c;:x :y :z])[(:a,:b)] --> (nalift([1 2 3;4 5 6]), nalift([:a :b :c;:x :y :z]))
    @fact @larr(@larr(a=[1 2 3;4 5 6],axis1[k=[:x,:y]],axis2[r=['a','b','c']]), axis1[[1,NA]],axis2[[3,NA,1]]) --> @larr(a=[1 2 3;4 5 6],axis1[1,NA],axis2[3,NA,1])
    @fact @larr(@larr(a=[1 2 3;4 5 6],axis1[k=[:x,:y]],axis2[r=['a','b','c']]), axis2[[3,NA,1]]) --> @larr(a=[1 2 3;4 5 6],axis1[k=[:x,:y]],axis2[[3,NA,1]])
    @fact @larr([1 2 3;4 5 6]) --> LabeledArray(nalift([1 2 3;4 5 6]))
  end
  context("additional method tests") do
    d1 = larr(:a=>nalift(reshape(1:20, 5, 4)),
              2=>nalift(fill(:sym, 5, 4)),
              :third=>nalift(reshape(map(i->string("str_",i), 1:20), 5, 4)))
    d2 = larr(:x=>nalift(reshape(1:20, 5, 4)),
              :third=>nalift(fill(:sym, 5, 4)),
              :z=>nalift(reshape(map(i->string("str_",i), 1:20), 5, 4)))
    @fact pick(merge(d1, d2), [:a,2]) --> pick(d1, [:a,2])
    @fact pick(merge(d1, d2), [:x,:third,:z]) --> pick(d2, [:x,:third,:z])
    @fact keys(peel(merge(d1, d2))) --> Any[:a,2,:third,:x,:z]
    @fact delete(merge(d1, d2), :a,2,:third) --> pick(d2, [:x, :z])
    @fact pick(d1, [:a]) --> DictArray(a=pick(d1, :a))
    @fact mapslices(x->LDict(:c=>Nullable(length(x))),@larr(a=[1 2 3;4 5 6],b=["a" "b" "c";"d" "e" "f"],axis1[k=[:x,:y]],axis2[r=[:m,:n,:p]]),[1]) --> @larr(c=[2,2,2], axis1[r=[:m,:n,:p]])
    @fact mapslices(x->LDict(:c=>Nullable(length(x))),@larr(a=[1 2 3;4 5 6],b=["a" "b" "c";"d" "e" "f"],axis1[k=[:x,:y]],axis2[r=[:m,:n,:p]]),[2]) --> @larr(c=[3,3], axis1[k=[:x,:y]])
    @fact mapslices(x->Nullable(length(x)),@larr(a=[1 2 3;4 5 6],b=["a" "b" "c";"d" "e" "f"]),[1]) --> LabeledArray(@nalift([2,2,2]), axes1=@nalift([1,2,3]))
    @fact size(mapslices(identity, @larr(a=[1 2 3;4 5 6]), [])) --> (2,3)
    @fact mapslices(x->LDict(:c1=>MultidimensionalTables.naop_plus(x[:a],x[:b]),:c2=>Nullable(10)), @larr(a=[1 2 3;4 5 6],b=[1.0 2.0 3.0;4.0 5.0 6.0]), []) --> @larr(c1=[2.0 4.0 6.0;8.0 10.0 12.0], c2=@rap reshape(_,(2,3)) fill (10,6)...)
    #@fact typeof(mapslices(identity, @larr(a=[1 2 3;4 5 6]), [])) --> Array{Pair{Nullable{Int64},Pair{Nullable{Int64},MultidimensionalTables.LDict{Symbol,Nullable{Int64}}}},2}
    @fact typeof(mapslices(identity, @larr(a=[1 2 3;4 5 6]), [])) --> MultidimensionalTables.LabeledArray{MultidimensionalTables.LDict{Symbol,Nullable{Int64}},2,Tuple{MultidimensionalTables.DefaultAxis,MultidimensionalTables.DefaultAxis},MultidimensionalTables.DictArray{Symbol,2,MultidimensionalTables.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}},Nullable{Int64}}}
    @fact eltype(mapslices(identity, @larr(a=[1 2 3;4 5 6]), [])) --> MultidimensionalTables.LDict{Symbol,Nullable{Int64}}
    #@fact eltype(mapslices(identity, @larr(a=[1 2 3;4 5 6]), [])) --> Pair{Nullable{Int64},Pair{Nullable{Int64},MultidimensionalTables.LDict{Symbol,Nullable{Int64}}}}
    @fact mapslices(x -> x, @larr(a=[1 2 3;4 5 6]), []) --> @larr(a=[1 2 3;4 5 6])
    @fact map(x->LDict(:c=>x[:a]), @larr(a=[1,2,3],b=[4,5,6],axis1[[:x,:y,:z]])) --> @larr(c=[1,2,3],axis1[[:x,:y,:z]])
    @fact map(x->x[:a], @larr(a=[1,2,3],b=[4,5,6],axis1[[:x,:y,:z]])) --> @larr([1,2,3],axis1[[:x,:y,:z]])
    @fact tbltool.create_dict(@larr(a=[1 NA NA;4 5 6],b=[NA NA 6; 7 8 9],axis2[r=[:x,:y,:z]]))[Nullable(2)][LDict(:r=>Nullable(:z))][:a].value --> 6
    @fact collect(keys(tbltool.create_dict(@larr(a=[1 NA NA;4 5 6],b=[NA NA 6; 7 8 9],axis2[r=[:x,:y,:z]]))[Nullable(1)])) --> [LDict(:r=>Nullable(:x)), LDict(:r=>Nullable(:z))]
    @fact tbltool.create_dict(@larr(a=[1 NA NA;4 5 6],b=[NA NA 6; 7 8 9],axis2[r=[:x,:y,:z]]))[Nullable(1)][LDict(:r=>Nullable(:x))] --> LDict(:a=>Nullable(1), :b=>Nullable{Int}())
    @fact reverse(@larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]]),[1]) --> @larr(a=[4 5 6;1 2 3],axis2[r=[:x,:y,:z]])
    @fact reverse(@larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]]),[2]) --> @larr(a=[3 2 1;6 5 4],axis2[r=[:z,:y,:x]])
    @fact reverse(@larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]]),1:2) --> @larr(a=[6 5 4;3 2 1],axis2[r=[:z,:y,:x]])
    @fact flipdim(@larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]]),2) --> @larr(a=[3 2 1;6 5 4],axis2[r=[:z,:y,:x]])
    @fact flipdim(@larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]]),1,2) --> @larr(a=[6 5 4;3 2 1],axis2[r=[:z,:y,:x]])
    @fact reshape(larr(a=[1 2 3;4 5 6],b=[10 11 12;13 14 15],axis1=darr(k1=[:a,:b],k2=[100,101]),axis2=[:m,:n,:p]),6) --> larr(a=[1,4,2,5,3,6],b=[10,13,11,14,12,15], axis1=darr(k1=repmat([:a,:b],3),k2=repmat([100,101],3),x1=[:m,:m,:n,:n,:p,:p]))
    @fact reshape(larr(a=[1 2 3;4 5 6],b=[10 11 12;13 14 15],axis1=darr(k1=[:a,:b],k2=[100,101]),axis2=[:m,:n,:p]),1,6) --> @rap transpose reshape(_,6,1) larr(a=[1,4,2,5,3,6],b=[10,13,11,14,12,15], axis1=darr(k1=repmat([:a,:b],3),k2=repmat([100,101],3),x1=[:m,:m,:n,:n,:p,:p]))
    @fact reshape(larr([1 2 3;4 5 6],axis1=[10,11]),1,6) --> larr([1 4 2 5 3 6], axis2=repmat([10,11],3))
    @fact reducedim((x,y)->x+y[:a].value, larr(a=reshape(1:24,2,3,4)),[1],0) --> larr([3 15 27 39;7 19 31 43;11 23 35 47])
    @fact reducedim((x,y)->x+y[:a].value, larr(a=reshape(1:24,2,3,4)),[1,2],0) --> larr([21,57,93,129])
    @fact reducedim((x,y)->x+y[:a].value, larr(a=reshape(1:24,2,3,4)),[1,2,3],0).value --> 300
    @fact reducedim((x,y)->x+y[:a].value, larr(a=reshape(1:24,2,3,4),axis1=darr(k1=[:a,:b]),axis3=[:x,:y,:z,:w]),[1],0) --> larr([3 15 27 39;7 19 31 43;11 23 35 47],axis2=[:x,:y,:z,:w])
    @fact reducedim((x,y)->x+y[:a].value, larr(a=reshape(1:24,2,3,4),axis1=darr(k1=[:a,:b]),axis3=[:x,:y,:z,:w]),[1,2],0) --> larr([21,57,93,129],axis1=[:x,:y,:z,:w])
    @fact reducedim((x,y)->x+y[:a].value, larr(a=reshape(1:24,2,3,4),axis1=darr(k1=[:a,:b]),axis3=[:x,:y,:z,:w]),[1,2,3],0).value --> 300
    # runnability test when the result is empty.
    @fact size(@rap transpose @select((@rap transpose larr(r=rand(8,20), axis1=darr(a=rand(8),b=101:108))), :r=_r.*100, where[_r.>1])) --> (0,0)
    @fact reorder(larr('x'=>1:10, 3=>11:20, axis1=101:110), 3) --> larr(3=>11:20, 'x'=>1:10, axis1=101:110)
    @fact reorder(larr(c1=1:10,c2=11:20),:c2,:c1) --> reorder(larr(c1=1:10,c2=11:20),:c2)
    @fact reorder(larr(c1=1:10,c2=11:20),:c2,:c1) --> larr(c2=11:20,c1=1:10)
    @fact rename(larr('x'=>1:10, 3=>11:20, axis1=101:110), 1) --> larr(1=>1:10, 3=>11:20, axis1=101:110)
    @fact rename(larr(c1=1:10,c2=11:20),:a,:b) --> larr(a=1:10,b=11:20)
    @fact cat(1, larr(a=[1 2 3;4 5 6], b=['a' 'b' 'c';'d' 'e' 'f'], axis1=[:u,:v]), larr(b=['x' 'y' 'z'], d=[:m :n :p], axis1=[3])) --> larr(reshape(@larr(a=[1,4,NA,2,5,NA,3,6,NA],b=['a','d','x','b','e','y','c','f','z'],d=[NA,NA,:m,NA,NA,:n,NA,NA,:p]), 3, 3), axis1=[:u,:v,3])
    @fact cat(2, larr(a=[1 2 3;4 5 6], b=['a' 'b' 'c';'d' 'e' 'f'], axis1=[:u,:v]), larr(b=['x','y'], d=[:m,:n], axis1=[:u,:v])) --> larr(reshape(@larr(a=[1,4,2,5,3,6,NA,NA], b=['a','d','b','e','c','f','x','y'], d=[NA,NA,NA,NA,NA,NA,:m,:n]), 2, 4), axis1=[:u,:v])
    @fact cat(2, larr(a=[1 2 3;4 5 6], b=['a' 'b' 'c';'d' 'e' 'f']), larr(b=[10,11], d=[:m,:n])) --> reshape(@larr(a=[1,4,2,5,3,6,NA,NA], b=['a','d','b','e','c','f',10,11], d=[NA,NA,NA,NA,NA,NA,:m,:n]), 2, 4)
    @fact cat(2, larr(a=[1 2 3;4 5 6], b=['a' 'b' 'c';'d' 'e' 'f'], axis2=darr(r=[:x,:y,:z])), larr(b=[10,11], d=[:m,:n])) --> larr(reshape(@darr(a=[1,4,2,5,3,6,NA,NA], b=['a','d','b','e','c','f',10,11], d=[NA,NA,NA,NA,NA,NA,:m,:n]), 2, 4), axis2=@darr(r=[:x,:y,:z,NA]))
    # this will cause a NullException().
    #@fact cat(2, larr(a=[1 2 3;4 5 6], b=['a' 'b' 'c';'d' 'e' 'f'],axis1=[100,200]), larr(b=[10 11]', d=[:m :n]',axis1=[100,200],axis2=darr(r=[:k]))) --> larr(reshape(@darr(a=[1,4,2,5,3,6,NA,NA], b=['a','d','b','e','c','f',10,11], d=[NA,NA,NA,NA,NA,NA,:m,:n]), 2, 4), axis1=[100,200],axis2=[1,2,3,LDict(:r=>Nullable(:k))])
    @fact cat(2, larr(a=[1 2 3;4 5 6], b=['a' 'b' 'c';'d' 'e' 'f'],axis1=[100,200]), larr(b=[10 11]', d=[:m :n]',axis1=[100,200],axis2=[:k])) --> larr(reshape(@darr(a=[1,4,2,5,3,6,NA,NA], b=['a','d','b','e','c','f',10,11], d=[NA,NA,NA,NA,NA,NA,:m,:n]), 2, 4), axis1=[100,200],axis2=[1,2,3,:k])
    @fact cat(2, larr(a=[1 2 3;4 5 6], b=['a' 'b' 'c';'d' 'e' 'f'], axis2=[:x,:y,:z]), larr(b=[10,11], d=[:m,:n])) --> larr(reshape(@darr(a=[1,4,2,5,3,6,NA,NA], b=['a','d','b','e','c','f',10,11], d=[NA,NA,NA,NA,NA,NA,:m,:n]), 2, 4), axis2=[Nullable(:x),Nullable(:y),Nullable(:z),Nullable{Symbol}()])
    @fact merge(larr(a=[1,2,3],b=[:x,:y,:z],axis1=[:a,:b,:c]),darr(c=[4,5,6],b=[:m,:n,:p]),darr(a=["X","Y","Z"])) --> larr(a=["X","Y","Z"],b=[:m,:n,:p],c=[4,5,6],axis1=[:a,:b,:c])
    @fact @larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]])[2:-1:1,2] --> larr(a=[5,2])
    @fact @larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]])[[2,1],2] --> larr(a=[5,2])
    @fact @larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]])[[2,1],2:3] --> larr(a=[5 6;2 3],axis2=darr(r=[:y,:z]))
    @fact sub(@larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]]),[2,1],2) --> getindex(@larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]]),[2,1],2)
    @fact sub(@larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]]),[2,1],2:3) --> getindex(@larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]]),[2,1],2:3)
    @fact slice(@larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]]),[2,1],2) --> larr(a=[5,2])
    @fact slice(@larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]]),1, 2:3) --> larr(a=[2,3], axis1=darr(r=[:y,:z]))
    @fact sub(@larr(a=[1 2 3;4 5 6],axis2[r=[:x,:y,:z]]),1, 2:3) --> larr(a=[2 3], axis2=darr(r=[:y,:z]))
  end
end

end


text:
function instance(;port=12385, log_dir="logs")
    procid = parse(Int,ENV["SLURM_PROCID"])
    nodes = ENV["SLURM_NODELIST"]
    nodes = readlines(`scontrol show hostname $(nodes)`)
    master_node = nodes[1]

    procid == 0 && mkpath(log_dir)

    @sync begin
        procid == 0 && server(port, log_dir)
        sleep(2)
        client(master_node, port, log_dir)
    end
end


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

JLLWrappers.@generate_main_file_header("gdk_pixbuf")
JLLWrappers.@generate_main_file("gdk_pixbuf", UUID("da03df04-f53b-5353-a52f-6a8b0620ced0"))
end  # module gdk_pixbuf_jll


text:
using Utils
using ProblemParser

using Lazy
using MLStyle

struct SNum
    str::AbstractString
end

@data SameOrChanged begin
    Same(SNum)
    Changed(SNum)
end
Changed(inner::String) = Changed(SNum(inner))

macro s_str(str)
    SNum(str)
end

const file = find_input(@__FILE__)
const input = parse(Lines(Map(SNum)), slurp(file))

function getinner(str::AbstractString)
    opened = 0
    comma = -1
    for i in eachindex(str)
        if str[i] == '['
            opened += 1
        elseif str[i] == ']'
            opened -= 1
            opened == 0 && return view(str, 2:i - 1), comma - 1, i - 1
        elseif str[i] == ',' && opened == 1
            comma = i
        end
    end
    error("closing bracket not found")
end

function replacematch1(str::AbstractString, (match, repl)::Pair; reversed=false)
    re = match === nothing ? str : string(
        str[begin:match.offset - 1],
        reversed ? reverse(string(repl)) : repl,
        str[match.offset + length(match.match):end],
    )
    if reversed
        re = reverse(re)
    end
    re
end

anumber = r"\d+"

function explode(num::SNum)::SameOrChanged
    str = num.str
    depth = 0
    for i in eachindex(str)
        c = str[i]
        if c == '['
            @views if depth == 4
                prefix = reverse(str[begin:i - 1])
                lhs, comma, bracket = getinner(str[i:end])
                rhs = lhs[comma + 1:end]
                lhs = lhs[begin:comma - 1]
                suffix = str[i + bracket + 1:end]

                pos = match(anumber, prefix)
                if pos !== nothing
                    number = parse.(Int, [lhs, reverse(pos.match)]) |> sum
                    prefix = replacematch1(prefix, pos => number; reversed=true)
                end

                pos = match(anumber, suffix)
                if pos !== nothing
                    number = parse.(Int, [rhs, pos.match]) |> sum
                    suffix = replacematch1(suffix, pos => number)
                end

                return string(prefix, 0, suffix) |> Changed
            else
                depth += 1
            end
        elseif c == ']'
            depth -= 1
        end
    end
    Same(num)
end

append(lhs::SNum, rhs::SNum) = string('[', lhs.str, ',', rhs.str, ']') |> SNum

append(lhs, rhs) = string('[', lhs, ',', rhs, ']')

function Base.split(num::SNum)::SameOrChanged
    for m in eachmatch(anumber, num.str)
        n = parse(Int, m.match)
        if n > 9
            lr = [n, n + 1] .÷ 2
            return replacematch1(num.str, m => append(lr...)) |> Changed
        end
    end
    Same(num)
end

function Base.reduce(num::SNum)::SNum
    while true
        num = @match explode(num) begin
            Changed(num) => num
            Same(num) => @match split(num) begin
                Changed(num) => num
                _ => return num
            end
        end
    end
end

Base.:+(lhs::SNum, rhs::SNum) = append(lhs, rhs) |> reduce

function magnitude(num::AbstractString)
    if isdigit(num[1])
        m = match(anumber, num)
        parse(Int, m.match)
    else
        lhs, comma = getinner(num)
        rhs = lhs[comma + 1:end]
        lhs = lhs[begin:comma - 1]
        3magnitude(lhs) + 2magnitude(rhs)
    end
end
magnitude(num::SNum) = magnitude(num.str)

part1(input) = @>> input foldl(+) magnitude

assertequal(part1(input), 3725)

part2(input) = maximum(magnitude(a + b) for a in input for b in input)

assertequal(part2(input), 4832)

using Test
@testset "18" begin
    @testset "addition" begin
        @test s"[[[[4,3],4],4],[7,[[8,4],9]]]" + s"[1,1]" == s"[[[[0,7],4],[[7,8],[6,0]]],[8,1]]"
    end
    @testset "replacematch1" begin
        target = "this is a foo string"
        @test replacematch1(target, match(r"fo+", target) => "bar") == "this is a bar string"
        target = reverse(target)
        @test replacematch1(target, match(r"o+f", target) => "bar"; reversed=true) == "this is a bar string"
    end
    @testset "explode" begin
        @test explode(s"[[[[[9,8],1],2],3],4]") == Changed(s"[[[[0,9],2],3],4]")
        @test explode(s"[7,[6,[5,[4,[3,2]]]]]") == Changed(s"[7,[6,[5,[7,0]]]]")
        @test explode(s"[[6,[5,[4,[3,2]]]],1]") == Changed(s"[[6,[5,[7,0]]],3]")
        @test explode(s"[[3,[2,[1,[7,3]]]],[6,[5,[4,[3,2]]]]]") == Changed(s"[[3,[2,[8,0]]],[9,[5,[4,[3,2]]]]]")
        @test explode(s"[[3,[2,[8,0]]],[9,[5,[4,[3,2]]]]]") == Changed(s"[[3,[2,[8,0]]],[9,[5,[7,0]]]]")
        @test explode(s"[[3,[2,[8,0]]],[9,[5,[4,1]]]]") == Same(s"[[3,[2,[8,0]]],[9,[5,[4,1]]]]")
    end
    @testset "split" begin
        @test split(s"[[[[0,7],4],[9,[0,9]]],[1,1]]") == Same(s"[[[[0,7],4],[9,[0,9]]],[1,1]]")
        @test split(s"[[[[0,7],4],[15,[0,13]]],[1,1]]") == Changed(s"[[[[0,7],4],[[7,8],[0,13]]],[1,1]]")
        @test split(s"[[[[0,7],4],[[7,8],[0,13]]],[1,1]]") == Changed(s"[[[[0,7],4],[[7,8],[0,[6,7]]]],[1,1]]")
    end
    @testset "reduce" begin
        @test reduce(s"[[[[[4,3],4],4],[7,[[8,4],9]]],[1,1]]") == s"[[[[0,7],4],[[7,8],[6,0]]],[8,1]]"
    end
    @testset "magnitude" begin
        @test magnitude(s"[[1,2],[[3,4],5]]") == 143
        @test magnitude(s"[[[[0,7],4],[[7,8],[6,0]]],[8,1]]") == 1384
        @test magnitude(s"[[[[1,1],[2,2]],[3,3]],[4,4]]") == 445
        @test magnitude(s"[[[[3,0],[5,3]],[4,4]],[5,5]]") == 791
        @test magnitude(s"[[[[5,0],[7,4]],[5,5]],[6,6]]") == 1137
        @test magnitude(s"[[[[8,7],[7,7]],[[8,6],[7,7]]],[[[0,7],[6,6]],[8,7]]]") == 3488
    end
end

max_stars_repo_name: CmdQ/AoC

max_issues_repo_name: CmdQ/AoC

max_forks_repo_name: CmdQ/AoC


text:
# utility/z-list.jl

include("downsample.jl")
export downsample1
export downsample2
export downsample3
export downsample

include("map_many.jl")
export map_many

include("mask.jl")
export mask_or
export mask_outline
export embed
export maskit

include("ndgrid.jl")
export ndgrid


text:
#_____________________________________________________________________
#                           Charlie's code
#

# Here I use `Gnuplot.jl`, but any other package would work...
#=
using Gnuplot
setverbose(false)
@gp "set size ratio -1" "set grid" "set key bottom right" xr=(-1.5, 2.5) :-
=#

# Methods to plot a Polygon
plotlabel(p::Polygon) = "Polygon (" * string(length(p.x)) * " vert.)"
plotlabel(p::RegularPolygon) = "RegularPolygon (" * string(vertices(p)) * " vert., area=" * string(round(area(p) * 100) / 100) * ")"
plotlabel(p::Named) = name(p) * " (" * string(vertices(p)) * " vert., area=" * string(round(area(p) * 100) / 100) * ")"
function plot(p::AbstractPolygon; dt=1, color="black")
    x = coords_x(p); x = [x; x[1]]
    y = coords_y(p); y = [y; y[1]]
    title = plotlabel(p)
    # @gp :- x y "w l tit '$title' dt $dt lw 2 lc rgb '$color'"
end

# Finally, let's have fun with the shapes!
line = Polygon([0., 1.], [0., 1.])
triangle = RegularPolygon(3, 1)
square = Named{RegularPolygon}("Square", 4, 1)

plot(line, color="black")
plot(triangle, color="blue")
plot(square, color="red")

rotate!(triangle, 90)
move!(triangle, 1, 1)
scale!(triangle, 0.32)
scale!(square, sqrt(2))

plot(triangle, color="blue", dt=2)
plot(square, color="red", dt=2)

# Increase number of vertices
p = Named{RegularPolygon}("Heptagon", 7, 1)
plot(p, color="orange", dt=4)
circle = Named{RegularPolygon}("Circle", 1000, 1)
plot(circle, color="dark-green", dt=4)


text:
###############################################################################
# Make regular snaps
###############################################################################
export @snap

"""
adds a trace entry in happysnaps for each expr in exprstrs with
corresponding val from vals
"""
storesnaps(location, everyN, exprstrs, vals) = begin
    if everyN == 1 || location_counts[location]%everyN == 1
        for (exprstr, val) in zip(exprstrs, vals)
            storesnap(location, exprstr, val, round(Int, location_counts[location]/everyN))
        end
    end
    location_counts[location] += 1
end

storesnap(location, exprstr, val, lcount) = push!(happysnaps, TraceItem(location, exprstr, val, lcount))

macro snapexprs(location, N, exprs)
    res = :(_rtexprstrs = []; _rtvals=[])
    for expr in exprs
        exprstr = string(expr)
        res = quote
            $res
            push!(_rtexprstrs, $exprstr)
            push!(_rtvals,
                try
                    $expr
                catch e
                    typeof(e) != UndefVarError && throw(e)
                    :undefined
                end)
        end
        autowatch && watch_exprstr(exprstr) #called at macro expansion time, not run time
    end
    res = :($res; RickTracy.storesnaps($location, $N, _rtexprstrs, _rtvals))
    res = :($res; RickTracy.happysnaps)
    res |> esc
end

"""
Take a snap/trace of a variable/expression.
#example

    fred = "flintstone"
    barney = 10

    @snap fred barney
    @tracevals fred

outputs:

    1-element Array{String,1}:
    "flintstone"

By default the variable/expression will be added to the watch list,
and logged/snapped on calls to `@snapall` that are parsed/loaded later than
this call.. To disable this behaviour call RickTracy.set_autowatch(false).

A numbered location string will be added to the trace entry to identify
the code location. To specify your own location use:

    @snap loc="decriptive location name" var1 var2
    #or
    @snapat location=@__LINE__ var1 var2

# n.b. `location`, `loc`, or just plain `l` are valid

    for person in ["wilma", "fred", "betty", "barney"]
        @snap N=2 person
    end
    @tracevals person

    returns:

    2-element Array{String,1}:
     "wilma"
     "betty"

"""
macro snap(exprs...)
    kwargs, kwextras, exprs, arginfo = kwargparse(trace_kwargspec, exprs)
    condition = kwargs[:iff]
    quote
        if $condition
            @RickTracy.snapexprs $(kwargs[:location]) $(kwargs[:everyN]) $exprs
        end
    end |> esc
end

max_stars_repo_name: JobJob/RickTracy


text:
using Test, Distributions, StatsBase

@testset "TDist(Inf) is Normal(0,1)" begin
    T = TDist(Inf)
    N = Normal(0, 1)

    x = rand(N)
    z = rand(T, 10000000)

    @test mean(T) == mean(N)
    @test median(T) == median(N)
    @test mode(T) == mode(N)
    @test var(T) == var(N)
    @test skewness(T) == skewness(N)
    @test kurtosis(T) == kurtosis(N)
    @test entropy(T) == entropy(N)
    @test pdf(T, x) ≈ pdf(N, x)
    @test logpdf(T, x) ≈ logpdf(N, x)
    @test gradlogpdf(T, x) ≈ gradlogpdf(N, x)
    @test cf(T, x) ≈ cf(N, x)

    fnecdf = ecdf(z)
    y = [-1.96, -1.644854, -1.281552, -0.6744898, 0, 0.6744898, 1.281552, 1.644854, 1.96]
    @test isapprox(fnecdf(y), [0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.975], atol=1e-3)
end


text:
# Simple Linear Gaussian State Space Model

module LGSSM1Benchmark

using Rocket
using ReactiveMP
using GraphPPL
using Distributions

@model function lgssm(n)
    x = randomvar(n)
    y = datavar(Float64, n)
    c = constvar(1.0)
    P = constvar(1.0)

    x_prior ~ NormalMeanVariance(0.0, 100.0)
    x_prev = x_prior

    for i in 1:n
        x[i] ~ x_prev + c
        y[i] ~ NormalMeanVariance(x[i], P)
        x_prev = x[i]
    end

    return x, y
end

generate_input(rng, n) = collect(1:n) + rand(rng, Normal(0.0, sqrt(1.0)), n);

function benchmark(input)
    n = length(input)

    _, (x, y) = lgssm(n);

    x_buffer  = buffer(Marginal, n)
    marginals = getmarginals(x)

    subscription = subscribe!(marginals, x_buffer)

    update!(y, input)

    unsubscribe!(subscription)

    return getvalues(x_buffer)
end

end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((1, 3, 2, false, true, false, false, false),),
    ((1, 3, 2, false, false, true, false, false),),
    ((1, 3, 2, true, false, false, true, false),),
    ((1, 3, 2, true, false, true, false, false),),
    ((1, 3, 2, true, false, false, true, true),),
    ]
# insts["fast"] = [] # TODO
insts["various"] = [
    ((5, 20, 20, false, true, false, false, false),),
    ((5, 20, 20, false, false, true, false, false),),
    ((5, 20, 20, true, false, false, true, false),),
    ((5, 20, 20, true, false, true, false, false),),
    ((5, 20, 20, true, false, false, true, true),),
    ((10, 30, 20, false, true, false, false, false),),
    ((10, 30, 20, false, false, true, false, false),),
    ((10, 30, 20, true, false, false, true, false),),
    ((10, 30, 20, true, false, true, false, false),),
    ((10, 30, 20, true, false, false, true, true),),
    ((8, 50, 30, false, true, false, false, false),),
    ((8, 50, 30, false, false, true, false, false),),
    ((8, 50, 30, true, false, false, true, false),),
    ((8, 50, 30, true, false, true, false, false),),
    ((8, 50, 30, true, false, false, true, true),),
    ]
return (SparseLMIJuMP, insts)


text:
export
    UnitGaussianNormalizer,
    encode,
    decode,
    get_burgers_data,
    get_darcy_flow_data

struct UnitGaussianNormalizer{T}
    mean::Array{T}
    std::Array{T}
    ϵ::T
end

function UnitGaussianNormalizer(𝐱; ϵ=1f-5)
    dims = 1:ndims(𝐱)-1

    return UnitGaussianNormalizer(mean(𝐱, dims=dims), StatsBase.std(𝐱, dims=dims), ϵ)
end

encode(n::UnitGaussianNormalizer, 𝐱::AbstractArray) = @. (𝐱-n.mean) / (n.std+n.ϵ)
decode(n::UnitGaussianNormalizer, 𝐱::AbstractArray) = @. 𝐱 * (n.std+n.ϵ) + n.mean

function register_burgers()
    register(DataDep(
        "Burgers",
        """
        Burgers' equation dataset from
        [fourier_neural_operator](https://github.com/zongyi-li/fourier_neural_operator)
        """,
        "https://drive.google.com/file/d/17MYsKzxUQVaLMWodzPbffR8hhDHoadPp/view?usp=sharing",
        "9cbbe5070556c777b1ba3bacd49da5c36ea8ed138ba51b6ee76a24b971066ecd",
        fetch_method=gdownload,
        post_fetch_method=unpack
    ))
end

function register_darcy_flow()
    register(DataDep(
        "DarcyFlow",
        """
        Darcy flow dataset from
        [fourier_neural_operator](https://github.com/zongyi-li/fourier_neural_operator)
        """,
        "https://drive.google.com/file/d/1zzVMuGhOG70EnR5L24LWqmX9-Wh_H5Wu/view?usp=sharing",
        "802825de9da7398407296c99ca9ceb2371c752f6a3bdd1801172e02ce19edda4",
        fetch_method=gdownload,
        post_fetch_method=unpack
    ))
end

function register_datasets()
    register_burgers()
    register_darcy_flow()
end

function get_burgers_data(; n=2048, Δsamples=2^3, grid_size=div(2^13, Δsamples), T=Float32)
    file = matopen(joinpath(datadep"Burgers", "burgers_data_R10.mat"))
    x_data = T.(collect(read(file, "a")[1:n, 1:Δsamples:end]'))
    y_data = T.(collect(read(file, "u")[1:n, 1:Δsamples:end]'))
    close(file)

    x_loc_data = Array{T, 3}(undef, 2, grid_size, n)
    x_loc_data[1, :, :] .= reshape(repeat(LinRange(0, 1, grid_size), n), (grid_size, n))
    x_loc_data[2, :, :] .= x_data

    return x_loc_data, y_data
end

function get_darcy_flow_data(; n=1024, Δsamples=5, T=Float32, test_data=false)
    # size(training_data) == size(testing_data) == (1024, 421, 421)
    file = test_data ? "piececonst_r421_N1024_smooth2.mat" : "piececonst_r421_N1024_smooth1.mat"
    file = matopen(joinpath(datadep"DarcyFlow", file))
    x_data = T.(permutedims(read(file, "coeff")[1:n, 1:Δsamples:end, 1:Δsamples:end], (3, 2, 1)))
    y_data = T.(permutedims(read(file, "sol")[1:n, 1:Δsamples:end, 1:Δsamples:end], (3, 2, 1)))
    close(file)

    x_dims = pushfirst!([size(x_data)...], 1)
    y_dims = pushfirst!([size(y_data)...], 1)
    x_data, y_data = reshape(x_data, x_dims...), reshape(y_data, y_dims...)

    x_normalizer, y_normalizer = UnitGaussianNormalizer(x_data), UnitGaussianNormalizer(y_data)

    return encode(x_normalizer, x_data), encode(y_normalizer, y_data), x_normalizer, y_normalizer
end


text:
using Base:NamedTuple
# abstract types -
abstract type VLAbstractGameWorld end
abstract type VLAbstractGameAgent end
abstract type VLAbstractGameStrategy end

# ----------------------------------------------------------------------------- #
# General -
struct VLMinorityGameSimulationResult
    data::Union{Nothing,DataFrame}
end
# ----------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------- #
# Basic -
mutable struct VLBasicMinorityGameStrategy <: VLAbstractGameStrategy

    # data -
    strategy::Array{Int,1}
    score::Int64

    function VLBasicMinorityGameStrategy(strategy::Array{Int,1}, score::Int64)
        _ = new(strategy, score)
    end
end

mutable struct VLBasicMinorityGameAgent <: VLAbstractGameAgent

    # data -
    agentStrategyCollection::Array{VLBasicMinorityGameStrategy,1}
    bestAgentStrategy::Union{Nothing,VLBasicMinorityGameStrategy}
    score::Union{Int64,Float64}

    function VLBasicMinorityGameAgent(agentStrategyCollection::Array{VLBasicMinorityGameStrategy,1};
        bestStrategy::Union{Nothing,VLBasicMinorityGameStrategy}=nothing, score::Union{Int64,Float64}=1)
        _ = new(agentStrategyCollection, bestStrategy, score)
    end
end

struct VLBasicMinorityGameWorld <: VLAbstractGameWorld

    # data -
    numberOfAgents::Int64
    agentMemorySize::Int64
    gameAgentArray::Array{VLBasicMinorityGameAgent,1}

    function VLBasicMinorityGameWorld(numberOfAgents::Int64, agentMemorySize::Int64, agentArray::Array{VLBasicMinorityGameAgent,1})
        _ = new(numberOfAgents, agentMemorySize, agentArray)
    end
end
# ----------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------- #
# Thermal -
mutable struct VLThermalMinorityGameStrategy <: VLAbstractGameStrategy

    # data -
    strategy::Array{Int,1}
    score::Int64
    probability::Float64

    function VLThermalMinorityGameStrategy(strategy::Array{Int,1}, score::Int64, probability::Float64)
        _ = new(strategy, score, probability)
    end
end

mutable struct VLThermalMinorityGameAgent <: VLAbstractGameAgent

    # data -
    agentStrategyCollection::Array{VLThermalMinorityGameStrategy,1}
    score::Int64

    function VLThermalMinorityGameAgent(agents::Array{VLThermalMinorityGameStrategy,1}, score::Int64)
        _ = new(agents, score)
    end
end

struct VLThermalMinorityGameWorld <: VLAbstractGameWorld

    # data -
    gameAgentArray::Array{VLThermalMinorityGameAgent,1}
    numberOfAgents::Int64
    agentMemorySize::Int64
    temperature::Float64

    function VLThermalMinorityGameWorld(gameAgentArray::Array{VLThermalMinorityGameAgent,1}, numberOfAgents::Int64,
        agentMemorySize::Int64, temperature::Float64)
        _ = new(gameAgentArray, numberOfAgents, agentMemorySize, temperature)
    end
end
# ----------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------- #
# Grand cannonical  -
mutable struct VLGCMinorityGameStrategy <: VLAbstractGameStrategy

    # data -
    strategy::Array{Int,1}
    score::Int64

    function VLGCMinorityGameStrategy(strategy::Array{Int,1}, score::Int64)
        _ = new(strategy, score)
    end
end

mutable struct VLGCMinorityGameAgent <: VLAbstractGameAgent

    # data -
    agentStrategyCollection::Array{VLGCMinorityGameStrategy,1}
    bestAgentStrategy::Union{Nothing,VLGCMinorityGameStrategy}
    score::Union{Int64,Float64}

    function VLGCMinorityGameAgent(agentStrategyCollection::Array{VLGCMinorityGameStrategy,1};
        bestStrategy::Union{Nothing,VLGCMinorityGameStrategy}=nothing, score::Union{Int64,Float64}=0)
        _ = new(agentStrategyCollection, bestStrategy, score)
    end
end

struct VLGCMinorityGameWorld <: VLAbstractGameWorld

    # data -
    numberOfAgents::Int64
    agentMemorySize::Int64
    gameAgentArray::Array{VLGCMinorityGameAgent,1}

    function VLGCMinorityGameWorld(numberOfAgents::Int64, agentMemorySize::Int64, agentArray::Array{VLGCMinorityGameAgent,1})
        _ = new(numberOfAgents, agentMemorySize, agentArray)
    end
end
# ----------------------------------------------------------------------------- #


text:
import Base: +, -, *, /, ^, !=, <, >, ==
import Base: isequal, isless, isapprox, cmp, size, getindex, broadcastable

"""
The Point type holds two coordinates. It's immutable, you can't
change the values of the x and y values directly.
"""
struct Point
   x::Float64
   y::Float64
end

"""
O is a shortcut for the current origin, `0/0`
"""
const O = Point(0, 0)

# basics
+(z::Number, p1::Point)              = Point(p1.x + z,    p1.y + z)
+(p1::Point, z::Number)              = Point(p1.x + z,    p1.y + z)
+(p1::Point, p2::Point)              = Point(p1.x + p2.x, p1.y + p2.y)
-(p1::Point, p2::Point)              = Point(p1.x - p2.x, p1.y - p2.y)
-(p::Point)                          = Point(-p.x,        -p.y)
-(z::Number, p1::Point)              = Point(p1.x - z,    p1.y - z)
-(p1::Point, z::Number)              = Point(p1.x - z,    p1.y - z)
*(k::Number, p2::Point)              = Point(k * p2.x,    k * p2.y)
*(p2::Point, k::Number)              = Point(k * p2.x,    k * p2.y)
/(p2::Point, k::Number)              = Point(p2.x/k,      p2.y/k)
*(p1::Point, p2::Point)              = Point(p1.x * p2.x, p1.y * p2.y)
^(p::Point, e::Integer)              = Point(p.x^e,       p.y^e)
^(p::Point, e::Float64)              = Point(p.x^e,       p.y^e)

# some refinements
# modifying points with tuples
+(p1::Point, shift::NTuple{2, Real}) = Point(p1.x + shift[1], p1.y + shift[2])
-(p1::Point, shift::NTuple{2, Real}) = Point(p1.x - shift[1], p1.y - shift[2])
*(p1::Point, shift::NTuple{2, Real}) = Point(p1.x * shift[1], p1.y * shift[2])
/(p1::Point, shift::NTuple{2, Real}) = Point(p1.x / shift[1], p1.y / shift[2])

# for broadcasting
Base.size(::Point) = 2
Base.getindex(p::Point, i) = [p.x, p.y][i]

Base.broadcastable(x::Point) = Ref(x)

# for iteration
Base.eltype(::Point) = Float64
Base.iterate(p::Point, state = 1) = state > length(p) ? nothing : (p[state], state + 1)
Base.length(::Point) = 2

"""
    dotproduct(a::Point, b::Point)

Return the scalar dot product of the two points.
"""
function dotproduct(a::Point, b::Point)
    result = 0.0
    for i in [:x, :y]
        result += getfield(a, i) * getfield(b, i)
    end
    return result
end

# comparisons

isequal(p1::Point, p2::Point)         = isapprox(p1.x, p2.x, atol=0.00000001) && (isapprox(p1.y, p2.y, atol=0.00000001))
isapprox(p1::Point, p2::Point)        = isapprox(p1.x, p2.x, atol=0.00000001) && (isapprox(p1.y, p2.y, atol=0.00000001))
isless(p1::Point, p2::Point)          = (p1.x < p2.x || (isapprox(p1.x, p2.x) && p1.y < p2.y))
!=(p1::Point, p2::Point)              = !isequal(p1, p2)
<(p1::Point, p2::Point)               = isless(p1, p2)
>(p1::Point, p2::Point)               = p2 < p1
==(p1::Point, p2::Point)              = isequal(p1, p2)

cmp(p1::Point, p2::Point)             = (p1 < p2) ? -1 : (p2 < p1) ? 1 : 0

# a unique that works better on points?
# I think this uses ==
# TODO perhaps unique(x -> round(x, sigdigits=13), myarray) ?
# "any implementation of unique with a tolerance will have some odd behaviors" Steven GJohnson
function Base.unique(pts::Array{Point, 1})
    apts = Point[]
    for pt in pts
        if pt ∉ apts
            push!(apts, pt)
        end
    end
    return apts
end

"""
    distance(p1::Point, p2::Point)

Find the distance between two points (two argument form).
"""
function distance(p1::Point, p2::Point)
    dx = p2.x - p1.x
    dy = p2.y - p1.y
    return sqrt(dx*dx + dy*dy)
end

"""
    pointlinedistance(p::Point, a::Point, b::Point)

Find the distance between a point `p` and a line between two points `a` and `b`.
"""
function pointlinedistance(p::Point, a::Point, b::Point)
  dx = b.x - a.x
  dy = b.y - a.y
  return abs(p.x * dy - p.y * dx + b.x * a.y - b.y * a.x) / hypot(dx, dy);
end

"""
    getnearestpointonline(pt1::Point, pt2::Point, startpt::Point)

Given a line from `pt1` to `pt2`, and `startpt` is the start of a perpendicular heading
to meet the line, at what point does it hit the line?
"""
function getnearestpointonline(pt1::Point, pt2::Point, startpt::Point)
    # first convert line to normalized unit vector
    dx = pt2.x - pt1.x
    dy = pt2.y - pt1.y
    mag = hypot(dx, dy)
    dx /= mag
    dy /= mag

    # translate the point and get the dot product
    lambda = (dx * (startpt.x - pt1.x)) + (dy * (startpt.y - pt1.y))
    return Point((dx * lambda) + pt1.x, (dy * lambda) + pt1.y)
end

"""
    midpoint(p1, p2)

Find the midpoint between two points.
"""
midpoint(p1::Point, p2::Point) = Point((p1.x + p2.x) / 2.0, (p1.y + p2.y) / 2.0)

"""
    midpoint(a)

Find midpoint between the first two elements of an array of points.
"""
midpoint(pt::Array) = midpoint(pt[1], pt[2])

"""
    between(p1::Point, p2::Point, x)
    between((p1::Point, p2::Point), x)

Find the point between point `p1` and point `p2` for `x`, where `x` is typically between 0
and 1. `between(p1, p2, 0.5)` is equivalent to `midpoint(p1, p2)`.
"""
function between(p1::Point, p2::Point, x)
    return p1 + (x * (p2 - p1))
end

function between(couple::NTuple{2, Point}, x)
    p1, p2 = couple
    return p1 + (x * (p2 - p1))
end

"""
    perpendicular(p1, p2, k)

Return a point `p3` that is `k` units away from `p1`, such that a line `p1 p3`
is perpendicular to `p1 p2`.

Convention? to the right?
"""
function perpendicular(p1::Point, p2::Point, k)
    px = p2.x - p1.x
    py = p2.y - p1.y
    l = hypot(px, py)
    if l > 0.0
        ux = -py/l
        uy = px/l
        return Point(p1.x + (k * ux), p1.y + (k * uy))
    else
        error("these two points are the same")
    end
end

"""
    perpendicular(p::Point)

Returns point `Point(p.y, -p.x)`.
"""
function perpendicular(p::Point)
    return Point(p.y, -p.x)
end

"""
    crossproduct(p1::Point, p2::Point)

This is the *perp dot product*, really, not the crossproduct proper (which is 3D):
"""
function crossproduct(p1::Point, p2::Point)
    return dotproduct(p1, perpendicular(p2))
end

function randomordinate(low, high)
    low + rand() * abs(high - low)
end

"""
    randompoint(lowpt, highpt)

Return a random point somewhere inside the rectangle defined by the two points.
"""
function randompoint(lowpt, highpt)
  Point(randomordinate(lowpt.x, highpt.x), randomordinate(lowpt.y, highpt.y))
end

"""
    randompoint(lowx, lowy, highx, highy)

Return a random point somewhere inside a rectangle defined by the four values.
"""
function randompoint(lowx, lowy, highx, highy)
    Point(randomordinate(lowx, highx), randomordinate(lowy, highy))
end

"""
    randompointarray(lowpt, highpt, n)

Return an array of `n` random points somewhere inside the rectangle defined by two points.
"""
function randompointarray(lowpt, highpt, n)
  array = Point[]
  for i in 1:n
    push!(array, randompoint(lowpt, highpt))
  end
  array
end

"""
    randompointarray(lowx, lowy, highx, highy, n)

Return an array of `n` random points somewhere inside the rectangle defined by the four
coordinates.
"""
function randompointarray(lowx, lowy, highx, highy, n)
    array = Point[]
    for i in 1:n
        push!(array, randompoint(lowx, lowy, highx, highy))
    end
    array
end
"""
    ispointonline(pt::Point, pt1::Point, pt2::Point;
        extended = false,
        atol = 10E-5)

Return `true` if the point `pt` lies on a straight line between `pt1` and `pt2`.

If `extended` is false (the default) the point must lie on the line segment between
`pt1` and `pt2`. If `extended` is true, the point lies on the line if extended in
either direction.
"""
function ispointonline(pt::Point, pt1::Point, pt2::Point;
        atol=10E-5,
        extended=false)
    dxc = pt.x - pt1.x
    dyc = pt.y - pt1.y
    dxl = pt2.x - pt1.x
    dyl = pt2.y - pt1.y
    cpr = (dxc * dyl) - (dyc * dxl)

    # point not on line
    if !isapprox(cpr, 0.0, atol=atol)
        return false
    end

    # point somewhere on extended line
    if extended == true
        return true
    end

    # point on the line
    if (abs(dxl) >= abs(dyl))
        return dxl > 0 ?
            pt1.x <= pt.x && pt.x <= pt2.x :
            pt2.x <= pt.x && pt.x <= pt1.x
    else
        return dyl > 0 ?
            pt1.y <= pt.y && pt.y <= pt2.y :
            pt2.y <= pt.y && pt.y <= pt1.y
    end
end

# """
#     intersection(p1::Point, p2::Point, p3::Point, p4::Point;
#         commonendpoints = false,
#         crossingonly = false,
#         collinearintersect = false)
#
# Find intersection of two lines `p1`-`p2` and `p3`-`p4`
#
# Deprecated. Use `intersectionlines().`
# """
# function intersection(A::Point, B::Point, C::Point, D::Point;
#         commonendpoints = false,
#         crossingonly = false,
#         collinearintersect = false
#         )
#     # false if either line is undefined
#     if (A == B) || (C == D)
#         return (false, Point(0, 0))
#     end
#
#     # false if the lines share a common end point
#     if commonendpoints
#         if (A == C) || (B == C) || (A == D) || (B == D)
#             return (false, Point(0, 0))
#         end
#     end
#
#     # Cramer's ?
#     A1 = (A.y - B.y)
#     B1 = (B.x - A.x)
#     C1 = (A.x * B.y - B.x * A.y)
#
#     L1 = (A1, B1, -C1)
#
#     A2 = (C.y - D.y)
#     B2 = (D.x - C.x)
#     C2 = (C.x * D.y - D.x * C.y)
#
#     L2 = (A2, B2, -C2)
#
#     d  = L1[1] * L2[2] - L1[2] * L2[1]
#     dx = L1[3] * L2[2] - L1[2] * L2[3]
#     dy = L1[1] * L2[3] - L1[3] * L2[1]
#
#     # if you ask me collinear points don't really intersect
#     if (C1 == C2) && (collinearintersect == false)
#         return (false, Point(0, 0))
#     end
#
#     if d != 0
#         pt = Point(dx/d, dy/d)
#         if crossingonly == true
#             if ispointonline(pt, A, B) && ispointonline(pt, C, D)
#                return (true, pt)
#             else
#                return (false, pt)
#             end
#         else
#             if ispointonline(pt, A, B, extended=true) && ispointonline(pt, C, D, extended=true)
#                 return (true, pt)
#             else
#                return (false, pt)
#             end
#         end
#     else
#         return (false, Point(0, 0))
#     end
# end

"""
    slope(pointA::Point, pointB::Point)

Find angle of a line starting at `pointA` and ending at `pointB`.

Return a value between 0 and 2pi. Value will be relative to the current axes.

    slope(O, Point(0, 100)) |> rad2deg # y is positive down the page
    90.0

    slope(Point(0, 100), O) |> rad2deg
    270.0

"""
function slope(pointA, pointB)
    return mod2pi(atan(pointB.y - pointA.y, pointB.x - pointA.x))
end

function intersection_line_circle(p1::Point, p2::Point, cpoint::Point, r)
    a = (p2.x - p1.x)^2 + (p2.y - p1.y)^2
    b = 2.0 * ((p2.x - p1.x) * (p1.x - cpoint.x) +
               (p2.y - p1.y) * (p1.y - cpoint.y))
    c = ((cpoint.x)^2 + (cpoint.y)^2 + (p1.x)^2 + (p1.y)^2 - 2.0 *
        (cpoint.x * p1.x + cpoint.y * p1.y) - r^2)
    i = b^2 - 4.0 * a * c
    if i < 0.0
        number_of_intersections = 0
        intpoint1 = O
        intpoint2 = O
    elseif isapprox(i, 0.0)
        number_of_intersections = 1
        mu = -b / (2.0 * a)
        intpoint1 = Point(p1.x + mu * (p2.x - p1.x), p1.y + mu * (p2.y - p1.y))
        intpoint2 = O
    elseif i > 0.0
        number_of_intersections = 2
        # first intersection
        mu = (-b + sqrt(i)) / (2.0 * a)
        intpoint1 = Point(p1.x + mu * (p2.x - p1.x),
                     p1.y + mu * (p2.y - p1.y))
        # second intersection
        mu = (-b - sqrt(i)) / (2.0 * a)
        intpoint2 = Point(p1.x + mu * (p2.x - p1.x), p1.y + mu * (p2.y - p1.y))
    end
    return number_of_intersections, intpoint1, intpoint2
end

"""
    intersectionlinecircle(p1::Point, p2::Point, cpoint::Point, r)

Find the intersection points of a line (extended through points `p1` and `p2`) and a circle.

Return a tuple of `(n, pt1, pt2)`

where

- `n` is the number of intersections, `0`, `1`, or `2`
- `pt1` is first intersection point, or `Point(0, 0)` if none
- `pt2` is the second intersection point, or `Point(0, 0)` if none

The calculated intersection points won't necessarily lie on the line segment between `p1` and `p2`.
"""
intersectionlinecircle(p1::Point, p2::Point, cpoint::Point, r) =
    intersection_line_circle(p1, p2, cpoint, r)

"""
    @polar (p)

Convert a tuple of two numbers to a Point of x, y Cartesian coordinates.

    @polar (10, pi/4)
    @polar [10, pi/4]

produces

    Luxor.Point(7.0710678118654755, 7.071067811865475)
"""
macro polar(p)
    quote
      Point($(esc(p))[1] * cos($(esc(p))[2]), $(esc(p))[1] * sin($(esc(p))[2]))
    end
end

"""
    polar(r, theta)

Convert point in polar form (radius and angle) to a Point.

    polar(10, pi/4)

produces

    Luxor.Point(7.071067811865475, 7.0710678118654755)
"""
polar(r, theta) = Point(r * cos(theta), r * sin(theta))

"""
    intersectionlines(p0, p1, p2, p3,
        crossingonly=false)

Find point where two lines intersect.

If `crossingonly == true` the point of intersection must lie on both lines.

If `crossingonly == false` the point of intersection can be where the lines meet
if extended almost to 'infinity'.

Accordng to this function, collinear, overlapping, and parallel lines never
intersect. Ie, the line segments might be collinear but have no points in
common, or the lines segments might be collinear and have many points in
common, or the line segments might be collinear and one is entirely contained
within the other.

If the lines are collinear and share a point in common, that
is the intersection point.
"""
function intersectionlines(p0::Point, p1::Point, p2::Point, p3::Point;
        crossingonly = false)
    resultflag = false
    resultip = Point(0.0, 0.0)
    if p0 == p1 # no lines at all
    elseif p2 == p3
    elseif p0 == p2 && p1 == p3 # lines are the same
    elseif p0 == p3 && p1 == p2
    elseif p0 == p2 # common end points
        resultflag = true
        resultip = p0
    elseif p0 == p3
        resultflag = true
        resultip = p0
    elseif p1 == p2
        resultflag = true
        resultip = p1
    elseif p1 == p3
        resultflag = true
        resultip = p1
    else
        # Cramers rule
        a1 = p0.y - p1.y
        b1 = p1.x - p0.x
        c1 = p0.x * p1.y - p1.x * p0.y

        l1 = (a1, b1, -c1)

        a2 = p2.y - p3.y
        b2 = p3.x - p2.x
        c2 = p2.x * p3.y - p3.x * p2.y

        l2 = (a2, b2, -c2)

        d  = l1[1] * l2[2] - l1[2] * l2[1]
        dx = l1[3] * l2[2] - l1[2] * l2[3]
        dy = l1[1] * l2[3] - l1[3] * l2[1]

        if !iszero(d)
            resultip = pt = Point(dx/d, dy/d)
            if crossingonly == true
                if ispointonline(resultip, p0, p1) && ispointonline(resultip, p2, p3)
                    resultflag = true
                else
                    resultflag = false
                end
            else
                if ispointonline(resultip, p0, p1, extended=true) &&
                   ispointonline(resultip, p2, p3, extended=true)
                    resultflag = true
                else
                    resultflag = false
                end
            end
        else
            resultflag = false
            resultip =  Point(0, 0)
        end
    end
    return (resultflag, resultip)
end

"""
    pointinverse(A::Point, centerpoint::Point, rad)

Find `A′`, the inverse of a point A with respect to a circle `centerpoint`/`rad`, such that:

```
distance(centerpoint, A) * distance(centerpoint, A′) == rad^2
```

Return (true, A′) or (false, A).
"""
function pointinverse(A::Point, centerpoint, rad)
    A == centerpoint && throw(error("pointinverse(): point $A and centerpoint $centerpoint are the same"))
    result = (false, A)
    n, C, pt2 = intersectionlinecircle(centerpoint, A, centerpoint, rad)
    if n > 0
        B = polar(rad, 0.7)                # arbitrary point on circle
        h = getnearestpointonline(B, C, A) # perp
        d = between(A, h, 2)               # reflection
        flag, A′ = intersectionlines(B, d, centerpoint, A, crossingonly = false)
        if flag == true
            result = (true, A′)
        end
    end
    return result
end

"""
    currentpoint()

Return the current point.
"""
function currentpoint()
    x, y = Cairo.get_current_point(get_current_cr())
    return Point(x, y)
end

"""
    hascurrentpoint()

Return true if there is a current point. Obtain the current point
with `currentpoint()`.
"""
function hascurrentpoint()
    return Cairo.has_current_point(get_current_cr())
end


text:
using BringRadical
using Test

@testset "BringRadical.jl" begin
    f(x) = x^5 + x
    for T in (Float32, Float64, BigFloat)
        for _ in 2^8
            for x in (rand(T), 20*rand(T), 100*rand(T))
                @test abs(f(bringrad(x)) - x) < eps(x)
                @test abs(f(bringrad(-x)) + x) < eps(x)
            end
        end
    end
end


text:
"""
    britter_puff_factory(scenario)

Generates a Britter-McQuaid dispersion model on the given scenario and returns a
function giving the concentration of the form
c(x, y, z, t)

"""
function britter_puff_factory(scenario)
    # TODO
    error("Not Implemented")
end


text:
A = 0.01
B = 1.6
C = 0.005
D = 1.0

#Diabatic interactions
function diabatic_potential(x)
    V = zeros(2,2)
    V[1] = (-sign(x))*A*(1-exp(-B*x*sign(x)))
    V[2] = C*exp(-D*x^2)
    V[4] = -V[1]
    V[3] = V[2]

    a = 1; b = -(V[1]+V[4]); c = det(V)
    lambdas = ( (-b + sqrt(disc(a,b,c)))/(2a), (-b - sqrt(disc(a,b,c)))/(2a) )
    return lambdas
end

function disc(a,b,c)
    return b^2 - 4*a*c
end


text:
## Fluid statics

export Circle, Semicircle, Rectangle, Triangle

abstract type AbstractShape end
abstract type Circle <: AbstractShape end
abstract type Rectangle <: AbstractShape end
abstract type Semicircle <: AbstractShape end
abstract type Triangle <: AbstractShape end

#ThermofluidQuantities.SpecificWeight(sg::SpecificGravity) = SpecificWeight(sg*SpecificWeight(Water))

##### Some geometric formulas #####
"""
    Area(::Circle,radius)

Calculate the area of a circle, given the `radius`.
"""
ThermofluidQuantities.Area(::Type{Circle},radius) = Area(π*radius^2)

"""
    Area(::Semicircle,radius)

Calculate the area of a semicircle, given the `radius`.
"""
ThermofluidQuantities.Area(::Type{Semicircle},radius) = Area(π*radius^2/2)

"""
    Area(::Rectangle,base,height)

Calculate the area of a rectangle, given the `base` and `height` lengths.
"""
ThermofluidQuantities.Area(::Type{Rectangle},base,height) = Area(base*height)

"""
    Area(::Triangle,base,height,apex)

Calculate the area of a triangle, given the lengths `base`, `height`, and `apex`,
which describes the horizontal distance from the left vertex to the apex. The base is
assumed to lie along the ``x`` axis.
"""
ThermofluidQuantities.Area(::Type{Triangle},base,height) = Area(base*height/2)

"""
    CentroidX(::Circle,radius), CentroidY(::Circle,radius)

Calculate the ``x`` and ``y`` coordinates of the centroid of a circle, given the `radius`.
The circle is assumed to be centered at the origin, so these simply return zero.
"""
CentroidX(::Type{Circle},radius) = CentroidX(0)
CentroidY(::Type{Circle},radius) = CentroidY(0)

"""
    CentroidX(::Semicircle,radius), CentroidY(::Semicircle,radius)

Calculate the ``x`` and ``y`` coordinates of the centroid of a semicircle, given the `radius`.
The semicircle is assumed to lie in the region above the ``x`` axis.
"""
CentroidX(::Type{Semicircle},radius) = CentroidX(0)
CentroidY(::Type{Semicircle},radius) = CentroidY(4*radius/3/π)

"""
    CentroidX(::Rectangle,base,height), CentroidY(::Rectangle,,base,height)

Calculate the ``x`` and ``y`` coordinates of the centroid of a rectangle, given the `base`
and `height`. The rectangle is assumed to be centered at the origin, so these simply return zero.
"""
CentroidX(::Type{Rectangle},base,height) = CentroidX(0)
CentroidY(::Type{Rectangle},base,height) = CentroidY(0)

"""
    CentroidX(::Triangle,base,height,apex), CentroidY(::Triangle,base,height,apex)

Calculate the area of a triangle, given the lengths `base`, `height`, and `apex`,
which describes the horizontal distance from the left vertex to the apex. The base is
assumed to lie along the ``x`` axis.
"""
CentroidX(::Type{Triangle},base,height,apex) = CentroidX((base+apex)/3)
CentroidY(::Type{Triangle},base,height,apex) = CentroidX(height/3)

"""
    SecondAreaMomentX(::Circle,radius), SecondAreaMomentY(::Circle,radius)

Calculate the second area moments of a circle about the centroid, given the `radius`.
"""
SecondAreaMomentX(::Type{Circle},radius) = SecondAreaMomentX(π*radius^4/4)
SecondAreaMomentY(::Type{Circle},radius) = SecondAreaMomentY(π*radius^4/4)

"""
    SecondAreaMomentX(::Semicircle,radius), SecondAreaMomentY(::Semicircle,radius)

Calculate the second area moments of a semicircle about the centroid, given the `radius`.
The semicircle is assumed to lie in the region above the ``x`` axis.
"""
SecondAreaMomentX(::Type{Semicircle},radius) = SecondAreaMomentX((π/8-8/9/π)*radius^4)
SecondAreaMomentY(::Type{Semicircle},radius) = SecondAreaMomentY(π*radius^4/8)

"""
    SecondAreaMomentX(::Rectangle,base,height), SecondAreaMomentY(::Rectangle,base,height)

Calculate the second area moments of a rectangle about the centroid, given the `base`
and `height` lengths. Note that ``x`` is parallel to the base and ``y`` to the height.
"""
SecondAreaMomentX(::Type{Rectangle},base,height) = SecondAreaMomentX(base*height^3/12)
SecondAreaMomentY(::Type{Rectangle},base,height) = SecondAreaMomentY(base^3*height/12)

"""
    SecondAreaMomentX(::Triangle,base,height,apex), SecondAreaMomentY(::Triangle,base,height,apex)

Calculate the second area moments of a triangle about the centroid, given the lengths `base`, `height`, and `apex`,
which describes the horizontal distance from the left vertex to the apex. The base is
assumed to lie along the ``x`` axis.
"""
SecondAreaMomentX(::Type{Triangle},base,height,apex) = SecondAreaMomentX(base*height^3/36)
SecondAreaMomentY(::Type{Triangle},base,height,apex) = SecondAreaMomentY((base^3*height-base^2*height*apex+base*height*apex^2)/36)

## Parallel axis theorem
"""
    SecondAreaMomentX(S,a...;x,y), SecondAreaMomentY(S,a...;x,y)

Calculate the second area moments about ``x`` and ``y``, using the
parallel axis theorem. If these arguments are omitted, it is assumed that
the moments are to be calculated about the centroid.
"""
SecondAreaMomentX(::Type{S},a...;x=CentroidX(S,a...),y=CentroidY(S,a...)) where S <: AbstractShape =
      SecondAreaMomentX(SecondAreaMomentX(S,a...) + Area(S,a...)*(CentroidY(S,a...)-y)^2)
SecondAreaMomentY(::Type{S},a...;x=CentroidX(S,a...),y=CentroidY(S,a...)) where S <: AbstractShape =
      SecondAreaMomentY(SecondAreaMomentY(S,a...) + Area(S,a...)*(CentroidX(S,a...)-x)^2)