text:
using Tokenize: Tokens

function is_template_expr(expr)
    isexpr(expr, CSTParser.UnarySyntaxOpCall) || return (false, nothing, false)
    if isexpr(children(expr)[2], OPERATOR, Tokens.DDDOT)
        isexpr(children(expr)[1], CSTParser.UnarySyntaxOpCall) || return (false, nothing, false)
        isexpr(children(children(expr)[1])[1], OPERATOR, Tokens.EX_OR) || return (false, nothing, false)
        return (true, Symbol( children(children(expr)[1])[2].val), true)
    end
    isexpr(children(expr)[1], OPERATOR, Tokens.EX_OR) || return (false, nothing, false)
    return (true, Symbol(children(expr)[2].val), false)
end

function matches_template(x, y)
    if typeof(x) != typeof(y)
        return false
    end
    if typeof(x) == CSTParser.IDENTIFIER
        return strip(x.val) == strip(y.val)
    end
    if typeof(x) in (CSTParser.KEYWORD, CSTParser.PUNCTUATION, CSTParser.OPERATOR)
        return x.kind == y.kind
    end
    true
end

matches_template(x::OverlayNode, y::OverlayNode) = matches_template(x.expr, y.expr)

function matches_template2(x, y)
    if (typeof(x) == BinarySyntaxOpCall && typeof(y) == CSTParser.WhereOpCall) ||
            (typeof(x) == CSTParser.WhereOpCall && typeof(y) == BinarySyntaxOpCall)
        return true
    end
    if typeof(x) != typeof(y)
        return false
    end
    if typeof(x) in (CSTParser.KEYWORD, CSTParser.PUNCTUATION, CSTParser.OPERATOR)
        if x.kind != y.kind
            return false
        end
    end
    if typeof(x) == CSTParser.IDENTIFIER
        if x.val != y.val
            return false
        end
    end
    true
end
matches_template2(x::OverlayNode, y::OverlayNode) = matches_template2(x.expr, y.expr)

struct EmptyMatch
    parent
end

function match_parameters(template, match, result)
    !matches_template2(template, match) && error("Shouldn't have gotten here $template, $match")
    j = 1
    for (i,x) in enumerate(children(template))
        if j > length(children(match)) && i == length(children(template))
            ret, sym, slurp = is_template_expr(x)
            if ret
                result[sym] = (("", "", ""), (EmptyMatch(match),))
                break
            end
        end
        if j > length(children(match))
            return false
        end
        y = children(match)[j]
        ret, sym, slurp = is_template_expr(x)
        if !ret && matches_template(x, y)
            ok = match_parameters(x, y, result)
            ok || return ok
            j += 1
        else
            if ret
                ws = is_last_leaf(x), string(prev_node_ws(x), leading_trivia(x)), string(trailing_trivia(x), next_node_ws(x))
                if !slurp
                    result[sym] =  (ws, (y,))
                    j += 1
                else
                    matched_exprs = Any[]
                    if i == length(children(template))
                        result[sym] = (ws, children(match)[j:end])
                        j = length(children(match))
                    else
                        nextx = children(template)[i + 1]
                        startj = j
                        while j <= length(children(match)) && !matches_template2(nextx, children(match)[j])
                            push!(matched_exprs, children(match)[j])
                            j += 1
                        end
                        result[sym] = (ws, isempty(matched_exprs) ? (EmptyMatch(children(match)[j]),) : matched_exprs)
                    end
                end
            else
                return false
            end
        end
    end
    if j < length(children(match))
        return false
    end
    return true
end

function leaf_is_template_expr(x)
    isempty(children(x)) && return (false, nothing, false)
    ret, sym, slurp = is_template_expr(x)
    ret || return leaf_is_template_expr(children(x)[1])
    ret, sym, slurp
end
is_template_expr(x::OverlayNode) = is_template_expr(x.expr)
leaf_is_template_expr(x::OverlayNode) = leaf_is_template_expr(x.expr)

using AbstractTrees: prevsibling, nextsibling

function prev_node_ws(node)
    ws = ""
    while true
        sib = prevsibling(node)
        if sib != nothing
            ws = string(trailing_trivia(sib), ws)
            !isempty(sib.span) && return ws
            node = sib
            continue
        end
        node.parent == nothing && return ws
        node = node.parent
    end
end

function next_node_ws(node)
    while true
        sib = nextsibling(node)
        sib != nothing && return leading_trivia(sib)
        node.parent == nothing && return ""
        node = node.parent
    end
end

function next_is_template(node)
    while true
        sib = nextsibling(node)
        sib != nothing && return leaf_is_template_expr(sib)
        node.parent == nothing && return (false, nothing, false)
        node = node.parent
    end
end

function is_last_leaf(node)
    while true
        sib = nextsibling(node)
        sib != nothing && return false
        node.parent == nothing && return true
        node = node.parent
    end
end

function reassemble(out::IO, replacement, matches)
    if isempty(children(replacement))
        ret, sym, _ = next_is_template(replacement)
        skip_trailing_trivia = ret && endswith(String(sym), "!")
        rt = text(replacement, true, !skip_trailing_trivia)
        write(out, rt)
    end
    i = 1
    while i <= length(children(replacement))
        x = children(replacement)[i]
        ret, sym, slurp = is_template_expr(x)
        if !ret
            reassemble(out, x, matches)
        else
            endswith(String(sym), "!") && (sym = Symbol(String(sym)[1:end-1]))
            exprs = matches[sym]
            write(out, prev_node_ws(first(exprs)))
            for expr in exprs
                write(out, fullspan_text(expr))
            end
        end
        i += 1
    end
    replacement
end

function skip_next_ws(c)
    ret, sym, _ = next_is_template(c)
    ret
end

"""
    For whitespace before or after a template parameter, given the whitespace:
        1. In the template
        2. In the corresponding replacement
        3. In the match

    Compute the whitespace that should end up in the final replacement.
"""
function process_whitespace(template, replacement, match)
    # The idea here to apply the "diff" between template and replacement
    # to match. For now, we require one to be a substring of the other.
    if template == replacement
        return match
    elseif length(replacement) < length(template)
        @assert startswith(template, replacement)
        # Any characters at the start of match that match `replacement` are fine
        i = 0
        while i < min(length(replacement), length(match))
            i = nextind(match, i)
            if replacement[i] != match[j]
                return match
            end
        end
        i > length(match) && return match
        j = nextind(match, i)
        # Then we cut any characters of the match that match the template
        while j <= min(length(template), length(match))
            if template[j] != match[j]
                break
            end
            j = nextind(match, j)
        end
        return string(match[1:i],match[j:end])
    elseif length(template) < length(replacement)
        @assert endswith(replacement, template)
        return string(replacement[end-sizeof(template):end], match)
    end
end

function reassemble_tree(replacement, matches, parent = nothing)
    # Whitespace hanlding: For any match, we need to process both the leading
    # and the trailing whitespace. Sometimes one template exprs trailing
    # whitespace is the next's leading whitespace. In that case, we break the
    # ties in favor of always processing leading ws (i.e. we skip trailing whitespace)
    if isempty(children(replacement))
        return skip_next_ws(replacement) ?
            TriviaReplacementNode(parent, replacement, leading_trivia(replacement), "") :
            replacement
    end
    ret = ChildReplacementNode(parent, Any[], replacement)
    for (i, x) in enumerate(children(replacement))
        istemp, sym, slurp = is_template_expr(x)
        if !istemp
            push!(ret.children, reassemble_tree(x, matches, ret))
        else
            use_template_trailing_trivia = true
            if endswith(String(sym), "!")
                use_template_trailing_trivia = false
                sym = Symbol(String(sym)[1:end-1])
            end
            (no_trailing_trivia, template_leading_trivia, template_trailing_trivia), exprs = matches[sym]
            expr = first(exprs)
            if typeof(expr) == EmptyMatch
                push!(ret.children, TriviaInsertionNode(prev_node_ws(expr.parent)))
                continue
            end
            lws = process_whitespace(template_leading_trivia,
                string(prev_node_ws(x), leading_trivia(x)),
                string(prev_node_ws(expr), leading_trivia(expr)),
            )
            tws = process_whitespace(template_trailing_trivia,
                string(next_node_ws(x), trailing_trivia(x)),
                string(next_node_ws(last(exprs)), trailing_trivia(last(exprs))),
            )
            if length(exprs) == 1
                push!(ret.children, TriviaReplacementNode(ret, expr,
                    lws, skip_next_ws(x) ? "" : tws))
            else
                push!(ret.children, TriviaReplacementNode(ret, expr,
                    lws, trailing_trivia(expr)))
                append!(ret.children, exprs[2:end-1])
                lexpr = last(exprs)
                push!(ret.children, TriviaReplacementNode(ret, lexpr,
                    leading_trivia(lexpr), skip_next_ws(x) ? "" : tws))
            end
        end
    end
    ret
end

function inspect_matches(result, text)
    for (sym, ((pre_ws, offset), expr)) in result
        print(stdout, sym, ": ")
        span = isa(expr, Array) ? (first(first(expr).fullspan):last(last(expr).fullspan)) : expr.fullspan
        show(stdout, text[1 + span])
        println(stdout)
    end
end


text:
Base.cmp(x::T, y::T) where {T<:Kmer} = cmp(packed_data(x), packed_data(y))
Base.:(==)(x::T, y::T) where {T<:Kmer} = packed_data(x) == packed_data(y)
Base.isless(x::T, y::T) where {T<:Kmer} = isless(packed_data(x), packed_data(y))

function Base.hash(x::Kmer{A,K,N}, h::UInt) where {A<:NucleicAcidAlphabet{2},K,N}
    return Base.hash(packed_data(x) ⊻ K ⊻ h)
end

ext: jl


text:
function breakpoint_hook(::typeof(JuliaInterpreter.breakpoint), bp)
    @debug "Breakpoint created" bp
    make_interpretted_mode_overdub(bp)
end

function breakpoint_hook(::typeof(JuliaInterpreter.remove), bp)
    @debug "Breakpoint Removed" bp
    remove_interpretted_mode_overdub(bp)
end

function breakpoint_hook(::typeof(JuliaInterpreter.update_states!), bp)
    is_enabled = bp.enabled[]
    @debug "Breakpoint update" bp is_enabled
    if is_enabled
        make_interpretted_mode_overdub(bp)
    else
        remove_interpretted_mode_overdub(bp)
    end
end


text:
# ---
# title: 208. Implement Trie (Prefix Tree)
# id: problem208
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Design, Trie
# link: <https://leetcode.com/problems/implement-trie-prefix-tree/description/>
# hidden: true
# ---
#
# Implement a trie with `insert`, `search`, and `startsWith` methods.
#
# **Example:**
#
#
#
#     Trie trie = new Trie();
#
#     trie.insert("apple");
#     trie.search("apple");   // returns true
#     trie.search("app");     // returns false
#     trie.startsWith("app"); // returns true
#     trie.insert("app");
#     trie.search("app");     // returns true
#
#
# **Note:**
#
#   * You may assume that all inputs are consist of lowercase letters `a-z`.
#   * All inputs are guaranteed to be non-empty strings.
#
#
## @lc code=start
using LeetCode

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


text:
# Interface function: nnmf

function nnmf(X::AbstractMatrix{T}, k::Integer;
              init::Symbol=:nndsvdar,
              alg::Symbol=:alspgrad,
              maxiter::Integer=100,
              tol::Real=cbrt(eps(T)/100),
              verbose::Bool=false) where T

    p, n = size(X)
    k <= min(p, n) || throw(ArgumentError("The value of k should not exceed min(size(X))."))

    # determine whether H needs to be initialized
    initH = alg == :projals

    # perform initialization
    if init == :random
        W, H = randinit(X, k; zeroh=!initH, normalize=true)
    elseif init == :nndsvd
        W, H = nndsvd(X, k; zeroh=!initH)
    elseif init == :nndsvda
        W, H = nndsvd(X, k; variant=:a, zeroh=!initH)
    elseif init == :nndsvdar
        W, H = nndsvd(X, k; variant=:ar, zeroh=!initH)
    else
        throw(ArgumentError("Invalid value for init."))
    end

    # choose algorithm
    if alg == :projals
        alginst = ProjectedALS{T}(maxiter=maxiter, tol=tol, verbose=verbose)
    elseif alg == :alspgrad
        alginst = ALSPGrad{T}(maxiter=maxiter, tol=tol, verbose=verbose)
    elseif alg == :multmse
        alginst = MultUpdate{T}(obj=:mse, maxiter=maxiter, tol=tol, verbose=verbose)
    elseif alg == :multdiv
        alginst = MultUpdate{T}(obj=:div, maxiter=maxiter, tol=tol, verbose=verbose)
    elseif alg == :cd
        alginst = CoordinateDescent{T}(maxiter=maxiter, tol=tol, verbose=verbose)
    else
        throw(ArgumentError("Invalid algorithm."))
    end

    # run optimization
    solve!(alginst, X, W, H)
end


text:
#=
    The File contains tests for some project utils and etc..
=#

#------------------------------------------------------------------------------

#using Suppressor

#------------------------------------------------------------------------------

@testset "Gizmos -- construct_jacobians" begin

    # simple tests #

    R, (x₁, x₂, x₃, x₄) = Nemo.QQ["x₁", "x₂", "x₃", "x₄"]

    system = [
        x₁^2 + 2*x₁*x₂,
        x₂^2 + x₃ + x₄,
        x₂ + x₄,
        x₁ + x₃
    ]
    true_js = [
        from_dense([
            2 0 0 0;
            2 0 0 0;
            0 0 0 0;
            0 0 0 0;
        ], Nemo.QQ),
        from_dense([
            2 0 0 0;
            0 2 0 0;
            0 0 0 0;
            0 0 0 0;
        ], Nemo.QQ),
        from_dense([
            0 0 0 1;
            0 0 1 0;
            0 1 0 1;
            0 1 1 0;
        ], Nemo.QQ)
    ]

    js = construct_jacobians(system)

    @test all(J in true_js for J in js) && length(true_js) == length(js)

    # test #2
    R, (x₁, x₂, x₃) = Nemo.QQ["x₁", "x₂", "x₃"]

    system = [
        x₂^2 + 4x₂*x₃ + 4x₃^2,  # derivative of x1
        4x₃ - 2x₁,              # derivative of x2
        x₁ + x₂                 # derivative of x3
    ]
    true_js = [
        from_dense([
            0 0 0;
            2 0 0;
            4 0 0;
        ], Nemo.QQ),
        from_dense([
            0 0 0;
            4 0 0;
            8 0 0;
        ], Nemo.QQ),
        from_dense([
            0 -2 1;
            0 0 1;
            0 4 0;
        ], Nemo.QQ)
    ]

    js = construct_jacobians(system)

    @test all(J in true_js for J in js) && length(true_js) == length(js)

end

@testset "General -- randomness" begin

    ε = 0.01
    for n in (100, )
        for d in (0.1, 0.5, 1)
            for 𝔽 in (QQ, GF(2^31-1))
                vector = random_sparsik(n, 𝔽, density=d)
                @test abs(density(vector) - d) < ε

                matrix = random_sparsik((n, n), 𝔽, density=d)
                @test abs(density(matrix) - d) < 2ε
            end
        end
    end

end

@testset "Subspaces -- check invariance" begin

    A = @sparse [1 0 0; 0 2 0; 0 0 3;]

    @test check_invariance!([], [A])

    v1 = @sparse [1,0,0]
    v2 = @sparse [0,1,0]
    v3 = @sparse [0,0,1]

    @test check_invariance!([A], [v1, v2, v3])
    @test check_invariance!([A], [v1])
    @test check_invariance!([A], [v2])
    @test check_invariance!([A], [v3])

    A = @sparse [1 2 0; 0 1 0; 0 0 3;]
    @test check_invariance!([A], [v1, v3])
    @test check_invariance!([A], [v1 + v2, v2])
    @test !check_invariance!([A], [v2])

    A = @sparse [0 4 0 1; 0 0 1 1; 0 0 0 2; 0 0 0 4]
    v1 = @sparse [1, 0, 0, 0]
    v2 = @sparse [0, 3//8, 1//2, 1]
    @test check_invariance!([A], [v1, v2])
    @test !check_invariance!([A], [v1 + v2])

end

@info "OK"


text:
import Base.rand

export BallInf,
       volume

"""
    BallInf{N, VN<:AbstractVector{N}} <: AbstractHyperrectangle{N}

Type that represents a ball in the infinity norm.

### Fields

- `center` -- center of the ball as a real vector
- `radius` -- radius of the ball as a real scalar (``≥ 0``)

### Notes

Mathematically, a ball in the infinity norm is defined as the set

```math
\\mathcal{B}_∞^n(c, r) = \\{ x ∈ \\mathbb{R}^n : ‖ x - c ‖_∞ ≤ r \\},
```
where ``c ∈ \\mathbb{R}^n`` is its center and ``r ∈ \\mathbb{R}_+`` its radius.
Here ``‖ ⋅ ‖_∞`` denotes the infinity norm, defined as
``‖ x ‖_∞ = \\max\\limits_{i=1,…,n} \\vert x_i \\vert`` for any
``x ∈ \\mathbb{R}^n``.

### Examples

Create the two-dimensional unit ball and compute its support function along the
positive ``x=y`` direction:

```jldoctest
julia> B = BallInf(zeros(2), 1.0)
BallInf{Float64, Vector{Float64}}([0.0, 0.0], 1.0)

julia> dim(B)
2

julia> ρ([1., 1.], B)
2.0
```
"""
struct BallInf{N, VN<:AbstractVector{N}} <: AbstractHyperrectangle{N}
    center::VN
    radius::N

    # default constructor with domain constraint for radius
    function BallInf(center::VN, radius::N) where {N, VN<:AbstractVector{N}}
        @assert radius >= zero(N) "radius must not be negative"
        return new{N, VN}(center, radius)
    end
end

isoperationtype(::Type{<:BallInf}) = false
isconvextype(::Type{<:BallInf}) = true

# --- AbstractHyperrectangle interface functions ---

"""
    radius_hyperrectangle(B::BallInf, i::Int)

Return the box radius of a ball in the infinity norm in a given dimension.

### Input

- `B` -- ball in the infinity norm
- `i` -- dimension of interest

### Output

The box radius of the ball in the infinity norm in the given dimension.
"""
function radius_hyperrectangle(B::BallInf, i::Int)
    return B.radius
end

"""
    radius_hyperrectangle(B::BallInf)

Return the box radius of a ball in the infinity norm, which is the same in every
dimension.

### Input

- `B` -- ball in the infinity norm

### Output

The box radius of the ball in the infinity norm.
"""
function radius_hyperrectangle(B::BallInf)
    return fill(B.radius, dim(B))
end

"""
    isflat(B::BallInf)

Determine whether a ball in the infinity norm is flat, i.e. whether its radius
is zero.

### Input

- `B` -- ball in the infinity norm

### Output

`true` iff the ball is flat.

### Notes

For robustness with respect to floating-point inputs, this function relies on
the result of `isapproxzero` when applied to the radius of the ball.
Hence, this function depends on the absolute zero tolerance `ABSZTOL`.
"""
function isflat(B::BallInf)
    return isapproxzero(B.radius)
end

function load_genmat_ballinf_static()
return quote
    function genmat(B::BallInf{N, SVector{L, N}}) where {L, N}
        r = B.radius
        if isapproxzero(B.radius)
            return SMatrix{L, 0, N, 0}()
        else
            gens = zeros(MMatrix{L, L})
            @inbounds for i in 1:L
                gens[i, i] = r
            end
            return SMatrix(gens)
        end
    end
end
end

# --- AbstractCentrallySymmetric interface functions ---

"""
    center(B::BallInf)

Return the center of a ball in the infinity norm.

### Input

- `B` -- ball in the infinity norm

### Output

The center of the ball in the infinity norm.
"""
function center(B::BallInf)
    return B.center
end

# --- LazySet interface functions ---

"""
    σ(d::AbstractVector, B::BallInf)

Return the support vector of a ball in the infinity norm in a given direction.

### Input

- `d` -- direction
- `B` -- ball in the infinity norm

### Output

The support vector in the given direction.
If the direction has norm zero, the vertex with biggest values is returned.
"""
function σ(d::AbstractVector, B::BallInf)
    @assert length(d) == dim(B) "a $(length(d))-dimensional vector is " *
                                "incompatible with a $(dim(B))-dimensional set"
    return center(B) .+ sign.(d) .* B.radius
end

# Particular dispatch for SingleEntryVector
function σ(d::SingleEntryVector, B::BallInf)
    return _σ_sev_hyperrectangle(d, B)
end

"""
    ρ(d::AbstractVector, B::BallInf)

Evaluate the support function of a ball in the infinity norm in a given
direction.

### Input

- `d` -- direction
- `B` -- ball in the infinity norm

### Output

Evaluation of the support function in the given direction.

### Algorithm

Let ``B`` be a ball in the infinity norm with center ``c`` and radius ``r`` and
let `d` be the direction of interest.
For balls with dimensions less than 30 we use the implementation for
`AbstractHyperrectangle`, taylored to a `BallInf`, which computes

```math
    ∑_{i=1}^n d_i * (c_i + \\textrm{sgn}(d_i) * r)
```

where ``\\textrm{sgn}(α) = 1`` if ``α ≥ 0`` and ``\\textrm{sgn}(α) = -1`` if ``α < 0``.

For balls of higher dimension, we instead exploit that for a support vector
``v = σ(d, B) = c + \\textrm{sgn}(d) * (r, …, r)ᵀ`` we have

```math
    ρ(d, B) = ⟨d, v⟩ = ⟨d, c⟩ + ⟨d, \\textrm{sgn}(d) * (r, …, r)ᵀ⟩ = ⟨d, c⟩ + r · ∑_{i=1}^n |d_i|
```

where ``⟨·, ·⟩`` denotes the dot product.
"""
function ρ(d::AbstractVector, B::BallInf)
    @assert length(d) == dim(B) "a $(length(d))-dimensional vector is " *
                                "incompatible with a $(dim(B))-dimensional set"
    c = center(B)
    r = B.radius
    if length(d) > 30
        # more efficient for higher dimensions
        return dot(d, c) + r * sum(abs, d)
    end
    N = promote_type(eltype(d), eltype(B))
    res = zero(N)
    @inbounds for (i, di) in enumerate(d)
        if di < zero(N)
            res += di * (c[i] - r)
        elseif di > zero(N)
            res += di * (c[i] + r)
        end
    end
    return res
end

# Particular dispatch for SingleEntryVector
function ρ(d::SingleEntryVector, B::BallInf)
    return _ρ_sev_hyperrectangle(d, B)
end

"""
    radius(B::BallInf, [p]::Real=Inf)

Return the radius of a ball in the infinity norm.

### Input

- `B` -- ball in the infinity norm
- `p` -- (optional, default: `Inf`) norm

### Output

A real number representing the radius.

### Notes

The radius is defined as the radius of the enclosing ball of the given
``p``-norm of minimal volume with the same center.
"""
function radius(B::BallInf, p::Real=Inf)
    return (p == Inf) ? B.radius : norm(fill(B.radius, dim(B)), p)
end

"""
    rand(::Type{BallInf}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing)

Create a random ball in the infinity norm.

### Input

- `BallInf` -- type for dispatch
- `N`       -- (optional, default: `Float64`) numeric type
- `dim`     -- (optional, default: 2) dimension
- `rng`     -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`    -- (optional, default: `nothing`) seed for reseeding

### Output

A random ball in the infinity norm.

### Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1.
Additionally, the radius is nonnegative.
"""
function rand(::Type{BallInf};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing)
    rng = reseed(rng, seed)
    center = randn(rng, N, dim)
    radius = abs(randn(rng, N))
    return BallInf(center, radius)
end

"""
    translate(B::BallInf, v::AbstractVector)

Translate (i.e., shift) a ball in the infinity norm by a given vector.

### Input

- `B` -- ball in the infinity norm
- `v` -- translation vector

### Output

A translated ball in the infinity norm.

### Algorithm

We add the vector to the center of the ball.

### Notes

See also [`translate!(::BallInf, ::AbstractVector)`](@ref) for the in-place version.
"""
function translate(B::BallInf, v::AbstractVector)
    return translate!(copy(B), v)
end

"""
    translate!(B::BallInf, v::AbstractVector)

Translate (i.e., shift) a ball in the infinity norm by a given vector, in-place.

### Input

- `B` -- ball in the infinity norm
- `v` -- translation vector

### Output

The ball `B` translated by `v`.

### Notes

See also [`translate(::BallInf, ::AbstractVector)`](@ref) for the out-of-place version.
"""
function translate!(B::BallInf, v::AbstractVector)
    @assert length(v) == dim(B) "cannot translate a $(dim(B))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    c = center(B)
    c .+= v
    return B
end

@inline function _vol_prod(B::BallInf{N}, n) where {N}
    vol = one(N)
    diam = 2 * B.radius
    for i in 1:n
        vol *= diam
    end
    return vol
end

# computation for floating-point results
function volume(B::BallInf{N}) where {N<:AbstractFloat}
    n = dim(B)
    if n < 50
        vol = _vol_prod(B, n)
    else
        r = B.radius
        vol = exp(n*log(2r))
    end
    return vol
end

function volume(B::BallInf)
    n = dim(B)
    vol = _vol_prod(B, n)
    return vol
end

function project(B::BallInf, block::AbstractVector{Int}; kwargs...)
    return BallInf(B.center[block], B.radius)
end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

fname = "sector-test.pdf"
Drawing(1500, 1500, fname) # points/pixels, 2000 points is 70.56cm × 70.56cm

setline(1)
setopacity(0.8)
origin()
background("black")

# Test endangle < angle
randomhue()
sector(O, 0, 20, pi/12, 0, :fill)
sethue("black")
sector(0, 20, pi/12, 0, :stroke)

# normal
for inner in 0:20:600
    for a in 0:pi/12:2pi
        randomhue()
        sector(O, inner, inner + 20, a, a + pi/12, :fill)
        sethue("black")
        sector(inner, inner + 20, a, a + pi/12, :stroke)
    end
end

# rounded
setopacity(0.75)
for i in 1:10
    gsave()
    translate(rand(-500:500), rand(-500:500))
    for inner in 90:30:210
        for a in 0:pi/12:2pi
            randomhue()
            sector(O, inner, inner + 20, a, a + pi/12, 15, :fill)
            sethue("black")
            sector(inner, inner + 20, a, a + pi/12, 10, :stroke)
        end
    end
    grestore()
end

@test finish() == true
println("...finished test: output in $(fname)")


text:
__precompile__()
module BayesLM

import Base.show
import Base.summary

using Distributions

export lm, summary, dic, show, cpo, Hyper, glm, bclogit

include("MCMC.jl")
include("common.jl")
include("lm.jl")
include("anova.jl")
include("glm.jl")
include("manova.jl")
include("survreg.jl")
include("bclogit.jl")

end # module BayesLM


text:
using eFEMpart
include("MeshTransform.jl")

point = [rand(),rand()]

p1 = 2*[-1.0,-1.0]
p2 = 2*[1.0,-1.0]
p3 = 2*[1.0,1.0]
p4 = 2*[-1.0,1.0]
quad = [p1,p2,p3,p4]

s,t = reverseQuadMap(quad,point)


text:
using ColorTypes, FixedPointNumbers
using Base.Test

if VERSION >= v"0.5.0-dev+2396"
    macro inferred5(ex)
        Expr(:macrocall, Symbol("@inferred"), esc(ex))
    end
else
    macro inferred5(ex)
        esc(ex)
    end
end

@test ColorTypes.to_top(AGray32(.8)) == ColorTypes.Colorant{FixedPointNumbers.UFixed{UInt8,8},2}
@test @inferred(eltype(Color{U8})) == U8
@test @inferred(eltype(RGB{Float32})) == Float32
@test @inferred(eltype(RGBA{Float64})) == Float64
# @test eltype(RGB) == TypeVar(:T, Fractional)
@inferred(eltype(RGB))      # just test that it doesn't error

@test length(RGB)    == 3
@test length(RGB1)   == 3
@test length(Gray)   == 1
@test length(ARGB)   == 4
@test length(RGB24)  == 3
@test length(ARGB32) == 4
@test length(AGray{Float32}) == 2

@test @inferred(color_type(RGB{U8})) == RGB{U8}
@test @inferred(color_type(RGB)) == RGB
@test @inferred(color_type(RGBA{Float32})) == RGB{Float32}
@test @inferred(color_type(GrayA{U8})) == Gray{U8}
@test @inferred(color_type(RGBA)) == RGB
@test @inferred(color_type(RGB24) ) == RGB24
@test @inferred(color_type(ARGB32)) == RGB24
@test @inferred(color_type(TransparentColor{RGB})) == RGB
@test @inferred(color_type(TransparentColor{RGB,Float64})) == RGB
@test @inferred(color_type(TransparentColor{RGB{Float64},Float64})) == RGB{Float64}
@test color_type(TransparentColor) <: Color
@test Color <: color_type(TransparentColor)
@test_throws MethodError color_type(Colorant{U8})

@test @inferred(base_color_type(RGBA{Float32})) == RGB
@test @inferred(base_color_type(ARGB{Float32})) == RGB
@test @inferred(base_color_type(BGR{U8})      ) == BGR
@test @inferred(base_color_type(HSV) ) == HSV
@test @inferred(base_color_type(HSVA)) == HSV
@test @inferred(base_color_type(TransparentColor{RGB{Float64},Float64})) == RGB

@test @inferred(base_colorant_type(RGBA{Float32})) == RGBA
@test @inferred(base_colorant_type(ARGB{Float32})) == ARGB
@test @inferred(base_colorant_type(BGR{U8})      ) == BGR
@test @inferred(base_colorant_type(HSV) ) == HSV
@test @inferred(base_colorant_type(HSVA)) == HSVA

@test @inferred(ccolor(Colorant{U8,3}, BGR{U8})) == BGR{U8}

@test @inferred(ccolor(RGB{Float32}, HSV{Float32})) == RGB{Float32}
@test @inferred(ccolor(RGB{U8},      HSV{Float32})) == RGB{U8}
@test @inferred(ccolor(RGB,          HSV{Float32})) == RGB{Float32}
@test @inferred(ccolor(ARGB{Float32}, HSV{Float32})) == ARGB{Float32}
@test @inferred(ccolor(ARGB{U8},      HSV{Float32})) == ARGB{U8}
@test @inferred(ccolor(ARGB,          HSV{Float32})) == ARGB{Float32}

@test @inferred(ccolor(Gray{U8}, Bool)) === Gray{U8}
@test @inferred(ccolor(Gray,     Bool)) === Gray{Bool}

@test @inferred(ccolor(RGB,  RGB))  === RGB
@test @inferred(ccolor(Gray, Gray)) === Gray

# Traits for instances (and their constructors)
@test @inferred(eltype(RGB{U8}(1,0,0))) == U8
@test @inferred(eltype(RGB(1.0,0,0))) == Float64
@test @inferred(eltype(ARGB(1.0,0.8,0.6,0.4))) == Float64
@test @inferred(eltype(RGBA{Float32}(1.0,0.8,0.6,0.4))) == Float32
@test @inferred(eltype(RGB(0x01,0x00,0x00))) == U8

@test length(ARGB(1.0,0.8,0.6,0.4)) == 4

@test @inferred(color_type(RGB{U8}(1,0,0))) == RGB{U8}
@test @inferred(color_type(ARGB(1.0,0.8,0.6,0.4))) == RGB{Float64}
@test @inferred(color_type(RGBA{Float32}(1.0,0.8,0.6,0.4))) == RGB{Float32}

@test @inferred(base_color_type(RGB{U8}(1,0,0))) == RGB
@test @inferred(base_color_type(ARGB(1.0,0.8,0.6,0.4))) == RGB
@test @inferred(base_color_type(RGBA{Float32}(1.0,0.8,0.6,0.4))) == RGB

@test @inferred(base_colorant_type(RGB{U8}(1,0,0))) == RGB
@test @inferred(base_colorant_type(ARGB(1.0,0.8,0.6,0.4))) == ARGB
@test @inferred(base_colorant_type(RGBA{Float32}(1.0,0.8,0.6,0.4))) == RGBA

# Constructors
@test eltype(Gray()) == U8
@test Gray(Gray()) == Gray()  # no StackOverflowError
for C in ColorTypes.parametric3
    @test eltype(C{Float32}) == Float32
    et = (C <: AbstractRGB) ? U8 : Float32
    @test eltype(C(1,0,0)) == et
    @test color_type(C(1,0,0)) == C{et}
    @test color_type(C) == C
    @test color_type(C{Float32}) == C{Float32}
    @test eltype(C{Float32}(1,0,0)) == Float32
    @test C(C()) == C()  # no StackOverflowError
end

# Specifically test the AbstractRGB types
# This checks that the constructor order is the same, even if the
# storage order is not
for C in subtypes(AbstractRGB)
    c = C(1, 0.5, 0)
    C == RGB24 && continue
    @test red(c)   == c.r == 1
    @test green(c) == c.g == 0.5
    @test blue(c)  == c.b == 0
end

c = Gray(0.8)
@test gray(c) == 0.8
@test gray(0.8) == 0.8
c = convert(Gray, 0.8)
@test c === Gray{Float64}(0.8)

c = AGray(0.8)
@test gray(c) == 0.8
@test color(c) == Gray(0.8)

c = convert(Gray, true)
@test c === Gray{Bool}(true)
@test gray(c) === true
@test gray(false) === false

# Transparency
@test alphacolor(Gray24(.2), .8) == AGray32(.2,.8)
@test alphacolor(RGB24(1,0,0), .8) == ARGB32(1,0,0,.8)
@test alphacolor(RGB(1,0,0), .8) == ARGB{U8}(1,0,0,.8)
@test coloralpha(RGB(1,0,0), .8) == RGBA{U8}(1,0,0,.8)
@test alphacolor(RGBA(1,0,0,.8)) == ARGB{Float64}(1,0,0,.8)
@test coloralpha(ARGB(1,0,0,.8)) == RGBA{Float64}(1,0,0,.8)
@test alphacolor(RGBA(1,0,0,.8)) == ARGB{Float64}(1,0,0,.8)
@test coloralpha(ARGB(1,0,0,.8)) == RGBA{Float64}(1,0,0,.8)
for C in setdiff(ColorTypes.parametric3, [RGB1,RGB4])
    for A in (alphacolor(C), coloralpha(C))
        @test eltype(A{Float32}) == Float32
        @test color_type(A) == C
        @test color_type(A{Float32}) == C{Float32}
        @test eltype(A(1,0.8,0.6,0.4)) == Float64
        c = A{Float64}(1,0.8,0.6,0.4)
        @test color_type(c) == C{Float64}
        cc = color(c)
        @test cc == C{Float64}(1,0.8,0.6)
        if VERSION >= v"0.4.0-dev"
            @test A(cc) == A{Float64}(1,0.8,0.6,1)
            @test A(cc, 0.4)  == c
            @test A(cc, 0x01) == A{Float64}(1,0.8,0.6,1)
            @test A{Float32}(cc, 0x01) == A{Float32}(1,0.8,0.6,1)
            @test C(c         ) == C{Float64}(1,0.8,0.6)
            @test C{Float32}(c) == C{Float32}(1,0.8,0.6)
        end
        @test convert(A, cc) == A{Float64}(1,0.8,0.6,1)
        @test convert(A, cc, 0.4)  == c
        @test convert(A, cc, 0x01) == A{Float64}(1,0.8,0.6,1)
        @test convert(A{Float32}, cc, 0x01) == A{Float32}(1,0.8,0.6,1)
        @test convert(C,          c) == C{Float64}(1,0.8,0.6)
        @test convert(C{Float32}, c) == C{Float32}(1,0.8,0.6)
    end
end
ac = ARGB32(rand(UInt32))
@test convert(ARGB32, ac) == ac
c = convert(RGB24, ac)
@test convert(RGB24, c) == c

crgb   = convert(RGB, c)
acargb = convert(ARGB, ac)
@test convert(Colorant,       crgb) === crgb
@test convert(Colorant{U8},   crgb) === crgb
@test convert(Colorant{Float32},       crgb) === convert(RGB{Float32}, crgb)
@test convert(Color,          crgb) === crgb
@test convert(Color{U8},      crgb) === crgb
@test convert(Color{Float32}, crgb) === convert(RGB{Float32}, crgb)
@test_throws ErrorException convert(TransparentColor, crgb)
@test_throws ErrorException convert(TransparentColor{RGB{U8}}, crgb)
@test_throws ErrorException convert(TransparentColor{RGB{U8},U8}, crgb)
@test_throws ErrorException convert(TransparentColor{RGB{U8},U8,2}, crgb)
@test convert(AlphaColor, crgb) === alphacolor(crgb)
@test convert(AlphaColor{RGB{U8}}, crgb) === alphacolor(crgb)
@test convert(AlphaColor{RGB{U8},U8}, crgb) === alphacolor(crgb)
@test convert(AlphaColor{RGB{U8},U8,4}, crgb) === alphacolor(crgb)
@test convert(ColorAlpha, crgb) === coloralpha(crgb)
@test convert(ColorAlpha{RGB{U8}}, crgb) === coloralpha(crgb)
@test convert(ColorAlpha{RGB{U8},U8}, crgb) === coloralpha(crgb)
@test convert(ColorAlpha{RGB{U8},U8,4}, crgb) === coloralpha(crgb)

@test convert(Colorant, acargb) === acargb
@test convert(Colorant{U8},   acargb) === acargb
@test convert(Colorant{U8,3}, acargb) === crgb
@test convert(TransparentColor, acargb) == acargb
@test convert(Color, acargb) == crgb

@test red(c)   == red(ac)
@test green(c) == green(ac)
@test blue(c)  == blue(ac)
ac2 = convert(ARGB32, c)
@test ac2.color == (c.color | 0xff000000)
@test color(c) == c
@test color(ac) == c
@test alpha(c) == U8(1)
@test alpha(ac) == UFixed8(ac.color>>24, 0)
@test alpha(ac2) == U8(1)
@test convert(RGB24,  0xff020304).color == 0xff020304
@test convert(ARGB32, 0x01020304).color == 0x01020304
ac3 = convert(RGBA, ac)
@test convert(RGB24, ac3) == c
ac4 = AGray32(.2,.8)
@test alpha(ac4) == .8
@test gray(ac4) == .2

for C in subtypes(AbstractRGB)
    rgb = convert(C, c)
    C == RGB24 && continue
    @test ccolor(Gray24, C) == Gray24
    @test ccolor(AGray32, C) == AGray32
    @test convert(AbstractRGB, c) == c
    @test convert(AbstractRGB{Float64}, rgb) === convert(C{Float64}, c)
    argb = convert(alphacolor(C), ac)
    rgba = convert(coloralpha(C), ac)
    @test rgb.r == red(c)
    @test rgb.g == green(c)
    @test rgb.b == blue(c)
    @test argb.alpha == alpha(ac)
    @test argb.r == red(ac)
    @test argb.g == green(ac)
    @test argb.b == blue(ac)
    @test rgba.alpha == alpha(ac)
    @test rgba.r == red(ac)
    @test rgba.g == green(ac)
    @test rgba.b == blue(ac)
end

@test convert(Float64, Gray(.3)) === .3
@test convert(GrayA{U8}, .2) == GrayA{U8}(.2)
@test convert(AGray{U8}, .2) == AGray{U8}(.2)
@test Gray{U8}(0.37).val           == U8(0.37)
@test convert(Gray{U8}, 0.37).val  == U8(0.37)
@test Gray24(0x0duf8).color           == 0x000d0d0d
@test convert(Gray24, 0x0duf8).color  == 0x000d0d0d
@test AGray32(0x0duf8).color          == 0xff0d0d0d
@test convert(AGray32, 0x0duf8).color == 0xff0d0d0d
@test AGray32(0x0duf8, 0x80uf8).color    == 0x800d0d0d
@test convert(Gray{UFixed16}, Gray24(0x0duf8)) == Gray{UFixed16}(0.05098)

iob = IOBuffer()
cf = RGB{Float32}(0.32218,0.14983,0.87819)
c  = convert(RGB{U8}, cf)
show(iob, c)
@test takebuf_string(iob) == "RGB{U8}(0.322,0.149,0.878)"
c = RGB{UFixed16}(0.32218,0.14983,0.87819)
show(iob, c)
@test takebuf_string(iob) == "RGB{UFixed16}(0.32218,0.14983,0.87819)"
c = RGBA{UFixed8}(0.32218,0.14983,0.87819,0.99241)
show(iob, c)
@test takebuf_string(iob) == "RGBA{U8}(0.322,0.149,0.878,0.992)"
showcompact(iob, c)
@test takebuf_string(iob) == "RGBA{U8}(0.322,0.149,0.878,0.992)"
show(iob, cf)
@test takebuf_string(iob) == "RGB{Float32}(0.32218f0,0.14983f0,0.87819f0)"
showcompact(iob, cf)
@test takebuf_string(iob) == "RGB{Float32}(0.32218,0.14983,0.87819)"

@test one(Gray{U8}) == Gray{U8}(1)
@test zero(Gray{U8}) == Gray{U8}(0)

c = Gray(0.8)
@test c == 0.8
show(iob, c)
@test takebuf_string(iob) == "Gray{Float64}(0.8)"
show(iob, AGray(0.8))
@test takebuf_string(iob) == "AGray{Float64}(0.8,1.0)"
show(iob, AGray32(0.8))
@test takebuf_string(iob) == "AGray32{U8}(0.8,1.0)"

# if the test below fails, please extend the list of types at the call to
# make_alpha in types.jl (this is the price of making that list explicit)
@test Set(ColorTypes.ctypes) ==
Set([DIN99d, DIN99o, DIN99, HSI, HSL, HSV, LCHab, LCHuv,
     LMS, Lab, Luv, XYZ, YCbCr, YIQ, xyY, BGR, RGB, Gray])

## operations
for T in (Gray{U8}, AGray{Float32}, GrayA{Float64},
          RGB{U8}, ARGB{U16}, RGBA{Float32},
          BGR{Float16}, RGB1{U8}, RGB4{Float64}, ABGR{U8})
    a = rand(T)
    @test isa(a, T)
    a = rand(T, (3, 5))
    @test isa(a, Array{T,2})
    @test size(a) == (3,5)
end

# colorfields
@test ColorTypes.colorfields(AGray32(.2)) == (:color,:alpha)
@test ColorTypes.colorfields(Gray) == (:val,)
@test ColorTypes.colorfields(RGB1) == (:r, :g, :b)
@test ColorTypes.colorfields(RGB4) == (:r, :g, :b)
@test ColorTypes.colorfields(BGR) == (:r, :g, :b)

# UInt32 comparison
@test Gray24() == 0x00000000
@test Gray24(.2) == 0x00333333
@test Gray24(0x23uf8) == 0x00232323
@test convert(UInt32, Gray24(0x23uf8)) === 0x00232323
@test AGray32() == 0xff000000
@test AGray32(.2) == 0xff333333
@test convert(AGray32, .2, 0.) == 0x00333333
@test AGray32(0x23uf8) == 0xff232323
@test convert(UInt32, AGray32(0x23uf8)) === 0xff232323
@test RGB24() == 0x00000000
@test RGB24(0x00232323) == 0x00232323
@test convert(UInt32, RGB24(0x00232323)) === 0x00232323
@test ARGB32() == 0xff000000
@test ARGB32(0xff232323) == 0xff232323
@test ARGB32(1,.2,.3) == 0xffff334c
@test convert(UInt32, ARGB32(1,.2,.3)) === 0xffff334c

@test @inferred5(mapc(sqrt, Gray{U8}(0.04))) == Gray(sqrt(U8(0.04)))
@test @inferred5(mapc(sqrt, AGray{U8}(0.04, 0.4))) == AGray(sqrt(U8(0.04)), sqrt(U8(0.4)))
@test @inferred5(mapc(sqrt, GrayA{U8}(0.04, 0.4))) == GrayA(sqrt(U8(0.04)), sqrt(U8(0.4)))
@test @inferred5(mapc(x->2x, RGB{U8}(0.04,0.2,0.3))) == RGB(map(x->2*U8(x), (0.04,0.2,0.3))...)
@test @inferred5(mapc(sqrt, RGBA{U8}(0.04,0.2,0.3,0.7))) == RGBA(map(x->sqrt(U8(x)), (0.04,0.2,0.3,0.7))...)
@test @inferred5(mapc(x->1.5f0x, RGBA{U8}(0.04,0.2,0.3,0.4))) == RGBA(map(x->1.5f0*U8(x), (0.04,0.2,0.3,0.4))...)

@test @inferred5(mapc(max, Gray{U8}(0.2), Gray{U8}(0.3))) == Gray{U8}(0.3)
@test @inferred5(mapc(-, AGray{Float32}(0.3), AGray{Float32}(0.2))) == AGray{Float32}(0.3f0-0.2f0,0.0)
@test @inferred5(mapc(min, RGB{U8}(0.2,0.8,0.7), RGB{U8}(0.5,0.2,0.99))) == RGB{U8}(0.2,0.2,0.7)
@test @inferred5(mapc(+, RGBA{U8}(0.2,0.8,0.7,0.3), RGBA{Float32}(0.5,0.2,0.99,0.5))) == RGBA(0.5f0+U8(0.2),0.2f0+U8(0.8),0.99f0+U8(0.7),0.5f0+U8(0.3))
@test_throws ArgumentError mapc(min, RGB{U8}(0.2,0.8,0.7), BGR{U8}(0.5,0.2,0.99))


text:
using Sugar: LazyMethod
using Transpiler: CLMethod
using Iterators
using GPUArrays
ctx = CLBackend.init()

f = Transpiler.CLFunction(
    Base._ind2sub, ((2,2,2), 2), ctx.queue
)

f, gltypes = Base._ind2sub, map(typeof, ((2,2,2), 2))
decl = CLMethod((f, gltypes))
funcsource = Sugar.getast!(decl)

for elem in funcsource.args[end].args[1].expressions
    println(elem)
end

tuple((int)((ind - l * indnext) + f)[1], _ind2sub_4(tail_2(inds), indnext)[1], _ind2sub_4(tail_2(inds), indnext)[2])

(tuple)(
    tuple((_3 - (_7 * _5)) + _6)[1],
    _ind2sub(tail(_2), _5)[1],
    _ind2sub(tail(_2), _5)[2]
)

tuple((+)((-)(_3, (*)(_7, _5)), _6))

tuple(
(int)((ind - l * indnext) + f)[1],
_ind2sub_4(tail_2(inds), indnext)[1],
_ind2sub_4(tail_2(inds), indnext)[2]
)


text:
using .GZip

# if GZip is loaded in, overwrite openvcf to this method
function openvcf(infile::String)
    if endswith(infile, ".gz")
        return GZip.open(infile, "r")
    else
        return open(infile, "r")
    end
end


text:
using Test, ArgParse, MAT

#=
The point of this script is to run a component of the test suite while not running the whole thing as the start up overhead is unreasonably long.

Environment requirements: Test, MAT

Example:
julia tinytest.jl --f test/test_utils.jl
=#

###########
## INPUT ##
###########

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--f"
        arg_type = String
        help = "input data file"
        required = true
    end

    return parse_args(s)
end

##########
## MAIN ##
##########

args = parse_commandline()

base = split(args["f"], "_")[end]

@testset "Tiny Test" begin
    include(string(pwd(), "/src/", base)) # source file
    include(string(pwd(), "/", args["f"])) # test file

end


text:
using Plots, MAT

nprocs = (2, 2) # nprocs (x, y) dim

@views function vizme2D_mpi(nprocs)
    H  = []
    ip = 1
    for ipx = 1:nprocs[1]
        for ipy = 1:nprocs[2]
            file = matopen("H_$(ip-1).mat"); H_loc = read(file, "H"); close(file)
            nx_i, ny_i = size(H_loc,1)-2, size(H_loc,2)-2
            ix1, iy1   = 1+(ipx-1)*nx_i, 1+(ipy-1)*ny_i
            if (ip==1)  H = zeros(nprocs[1]*nx_i, nprocs[2]*ny_i)  end
            H[ix1:ix1+nx_i-1,iy1:iy1+ny_i-1] .= H_loc[2:end-1,2:end-1]
            ip += 1
        end
    end
    fontsize = 12
    opts = (aspect_ratio=1, yaxis=font(fontsize, "Courier"), xaxis=font(fontsize, "Courier"),
        ticks=nothing, framestyle=:box, titlefontsize=fontsize, titlefont="Courier",
        xlabel="Lx", ylabel="Ly", xlims=(1, size(H,1)), ylims=(1, size(H,2)) )
    display(heatmap(H'; c=:davos, title="diffusion 2D MPI", opts...))
    return
end

vizme2D_mpi(nprocs)


text:
__precompile__()
module MDP_HSM
using DataFrames
using JuMP
using Cbc
using Query
using DataValues
using CSV
using Dates

export MDP_HSM_Model,RunModel

struct MDP_HSM_Model
    workFolder::String
    logFile::IOStream
    dateFrmt::DateFormat
    dfOrders::DataFrame
    dfRounds::DataFrame# for each round min max volume + min max occurence
    dfFlows::DataFrame# for each flow min max volume
    params::Dict{String,String}
end

CSV2DF(path::AbstractString)= DataFrame(CSV.File(path))

#would be great to precomplie this one ...
function df2ParamDict(dfParams::DataFrame)
    params=@from rds in dfParams begin
        @select rds.Key=>rds.Value
        @collect Dict
    end
    return params::Dict{String,String}
end
toList(x)=split(x,"#")

#function MDP_HSM_Model(path::AbstractString; orderFile::AbstractString="HSMOrders.csv", roundFile::AbstractString="HSMRounds.csv",flowFile::AbstractString="HSMFlows.csv",paramFile::AbstractString="HSMParams.csv",dateFrmt::DateFormat = DateFormat("mm/dd/yyyy"))
function MDP_HSM_Model(path::String; orderFile::String="HSMOrders.csv", roundFile::String="HSMRounds.csv",flowFile::String="HSMFlows.csv",paramFile::String="HSMParams.csv",dateFrmt::DateFormat = DateFormat("m/d/y"))
    logFile=open(joinpath(path, "Logging.txt"),"w")

    write(logFile, "Reading $path $roundFile\r\n")
    dfRounds=CSV2DF(joinpath(path,roundFile))

    write(logFile,"Reading $path $flowFile\r\n")
    dfFlows=CSV2DF(joinpath(path,flowFile))

    if(false && :Active in names(dfFlows))
        dfFlows=@from fl in dfFlows begin
        @where fl.Active==1
        @select fl
        @collect DataFrame
        end
    end

    write(logFile, "Reading $path $paramFile\r\n")
    dfParams=CSV2DF(joinpath(path,paramFile))

    params=df2ParamDict(dfParams)

    write(logFile, "Reading $path $orderFile\r\n")
    dfOrders=DataFrame(CSV.File(joinpath(path,orderFile);delim=";",types=Dict("Works_Order_No"=>Union{String,Missing},"Expedite_Level"=>Union{String,Missing},"Furnace_Group"=>Union{String,Missing},"Galv_Options"=>Union{String,Missing},"CULPST"=>Union{Date,Missing},"HSM_LPST"=>Union{Date,Missing},"HSM_PSD"=>Union{Date,Missing},"SO_Due_Date"=>Union{Date,Missing},"Slab_Grade"=>Union{String,Missing}) ,dateformat=dateFrmt))
    #showall(dfOrders)

    #remove trials , invalid ROUND_ID
    dfOrders=@from ord in dfOrders begin
        @where (ord.Customer_Name != "AM/NS Calvert Quality Internal Trials") && (!(ord.Round_ID in("NOK_Width" ,"NOK_MES Unavailable","NOK_No Location","NOK_River Terminal")))
        @select ord
        @collect DataFrame
    end

    if (!(:Flows in propertynames(dfOrders)) && (:Flow in propertynames(dfOrders)))
        write(logFile,"Using <Flow> for <Flows>\r\n")
        rename!(dfOrders,(:Flow=>:Flows))
    end

    dfOrders.FlowList = map( (x) -> toList(x),dfOrders.Flows)
    dfOrders.RoundList= map( (x) -> toList(ismissing(x) ? "???" : replace(x,"IF_BH" => "IF")),dfOrders.Round_Type)
    dfOrders.RoundList= map( (x) -> ("IF" in(x)) ? push!(x,"IF_Exposed") : x ,dfOrders.RoundList)

    if (!(:Volume in propertynames(dfOrders)) && (:Slab_Weight in propertynames(dfOrders)))
        dfOrders.Volume = map( (x) -> x/1000,dfOrders.Slab_Weight)
        dfOrders.WidthGroup =map( (x) -> floor(Int,x/100),dfOrders.Aim_Width)
    end
    #show(aMDPModel.dfOrders, allrows=true,allcols=true)
   #show(aMDPModel.dfOrders, allcols=true)
    #println(eltypes(dfOrders))
    #showall(dfOrders)

    res= MDP_HSM_Model(path,logFile,dateFrmt,dfOrders,dfRounds,dfFlows,params)
    write(logFile,"Model created\r\n")
    return res
end

#precomp hint

CSV2DF(joinpath(@__DIR__,"../data/HSMParams.csv"))
#no gain
#toList("GI")
# not working... [1] query(::DataFrames.DataFrame) at C:\Users\10500508\.julia\v0.6\QueryOperators\src\source_iterable.jl:
#df2ParamDict(CSV2DF(joinpath(@__DIR__,"../data/HSMParams.csv")))

function orderInFlow(orderFlow::Array, flow::AbstractString)
#println(typeof(orderFlow))
if in(flow, orderFlow)
    return 1
end
flowSum=split(flow,"+")# we have a case like GI_3+GI_34
if length(flowSum)==1
    return 0
end
for f in flowSum
    if in(f, orderFlow)
        #println(orderFlow, " in ", f)
        return 1
   end
end
return 0
end
#orderInFlow("GI3","GI3+GI34")#precomp hint

function RunModel(aMDPModel::MDP_HSM_Model;RoundLimitsAsConstraint::Bool=true)
#write(aMDPModel.logFile,"Building Model\r\n")
m=Model(optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1))
nOrders=size(aMDPModel.dfOrders,1)

@variable(m,VolInRd[i=1:nOrders,r in aMDPModel.dfOrders[i,:RoundList]]>=0)
@variable(m,VolOuterBayInRd[i=1:nOrders,r in aMDPModel.dfOrders[i,:RoundList]]>=0)

@variable(m,RdPerWidthGroupLength[r in aMDPModel.dfRounds.RoundName,w=1:20]>=0)

@variable(m,Flow[f in aMDPModel.dfFlows.FlowName]>=0)

@variable(m,FlowShortage[f in aMDPModel.dfFlows.FlowName]>=0)
@variable(m,FlowExcess[f in aMDPModel.dfFlows.FlowName]>=0)

minOcc=@from rds in aMDPModel.dfRounds begin
    @select rds.RoundName=>  rds.MinOccurence
    @collect Dict
end

maxOcc=@from rds in aMDPModel.dfRounds begin
    @select rds.RoundName=> rds.MaxOccurence
    @collect Dict
end

MinVolumePerRound=@from rds in aMDPModel.dfRounds begin
    @select rds.RoundName=>rds.MinVolPerRound
    @collect Dict
end

MaxVolumePerRound=@from rds in aMDPModel.dfRounds begin
    @select rds.RoundName=>rds.MaxVolPerRound
    @collect Dict
end

OuterBayVolumePerRound=@from rds in aMDPModel.dfRounds begin
    @select rds.RoundName=>rds.OuterbayVol
    @collect Dict
end

#show(aMDPModel.dfFlows, allrows=true,allcols=true)
minFlow=@from fl in aMDPModel.dfFlows begin
    @select fl.FlowName=>get(fl.FlowMin, 0)
    @collect Dict
end
#show(minFlow)

maxFlow=@from fl in aMDPModel.dfFlows begin
    @select fl.FlowName=>get(fl.FlowMax,0)
    @collect Dict
end

@variable(m, Rd[r in aMDPModel.dfRounds.RoundName] <=maxOcc[r], lower_bound=minOcc[r], Int)

@variable(m,RdShortage[r in aMDPModel.dfRounds.RoundName]>=0)
@variable(m,RdExcess[r in aMDPModel.dfRounds.RoundName]>=0)

for i=1:nOrders
    @constraint(m,sum(VolInRd[i,r] for r in aMDPModel.dfOrders[i,:RoundList]) <= aMDPModel.dfOrders[i,:Volume])
#unknown flow/round
    for r in aMDPModel.dfOrders[i,:RoundList]
        if !haskey( minOcc,r)
            println("Order $i has invalid round $r")
            @constraint(m,VolInRd[i,r] ==0 )
        end
        @constraint(m,VolOuterBayInRd[i,r]==(aMDPModel.dfOrders[i,:Yard_Location]  =="OuterBay" ? VolInRd[i,r] : 0))
    end
    #for f in aMDPModel.dfOrders[i,:FlowList]
    #    if !haskey(minFlow,f)
    #        @constraint(m,sum(VolInRd[i,r] for r in aMDPModel.dfOrders[i,:RoundList]) ==0 )
    #    end
    #end
end

if haskey(aMDPModel.params,"OuterbayVol")
    @constraint(m,sum(VolOuterBayInRd[i,r] for i=1:nOrders, r in (aMDPModel.dfOrders[i,:RoundList]) ) <= parse(Float64, aMDPModel.params["OuterbayVol"]))
end

@constraint(m, sum(Rd[r] for r in aMDPModel.dfRounds.RoundName) >= parse(Int32, aMDPModel.params["Min_Rounds"]))
@constraint(m, sum(Rd[r] for r in aMDPModel.dfRounds.RoundName) <= parse(Int32, aMDPModel.params["Max_Rounds"]))

@constraint(m, conRdPerWidthGroupLength[r in (aMDPModel.dfRounds.RoundName) ,w=1:20] ,RdPerWidthGroupLength[r,w] == sum(VolInRd[i,r]/aMDPModel.dfOrders[i,:Volume]*aMDPModel.dfOrders[i,:Coil_Length] for i=1:nOrders if aMDPModel.dfOrders[i,:WidthGroup]==w && r in aMDPModel.dfOrders[i,:RoundList]))

if haskey(aMDPModel.params,"WidthGroupLimit") && aMDPModel.params["WidthGroupLimit"]!=0
    write(aMDPModel.logFile,"Adding Width Group limit constraint.\n")
    @constraint(m, conRdPerWidthGroupLengthLimit[r in (aMDPModel.dfRounds.RoundName) ,w=1:20] ,RdPerWidthGroupLength[r,w] <=55*Rd[r])
end

@variable(m,RdVol[r in aMDPModel.dfRounds.RoundName])
@variable(m,totalVol)

for r in aMDPModel.dfRounds.RoundName
    @constraint(m,RdVol[r]==sum(VolInRd[i,r] for i=1:nOrders if r in aMDPModel.dfOrders[i,:RoundList]))
    @constraint(m,RdShortage[r] >= MinVolumePerRound[r]*Rd[r] - RdVol[r])
    @constraint(m,RdExcess[r]>=RdVol[r] - Rd[r]*MaxVolumePerRound[r] )
    if RoundLimitsAsConstraint
        @constraint(m,RdShortage[r]==0)
        @constraint(m,RdExcess[r]==0)
    end
    @constraint(m,sum(VolOuterBayInRd[i,r] for i=1:nOrders if r in aMDPModel.dfOrders[i,:RoundList])<=OuterBayVolumePerRound[r])
end

@constraint(m,totalVol==sum(RdVol[r] for r in (aMDPModel.dfRounds.RoundName) ))
if haskey(aMDPModel.params,"MinVolume")
    @constraint(m,totalVol >= parse(Float64, aMDPModel.params["MinVolume"]))
end
if haskey(aMDPModel.params,"MaxVolume")
    @constraint(m,totalVol <= parse(Float64, aMDPModel.params["MaxVolume"]))
end

for f in aMDPModel.dfFlows.FlowName
    #@constraint(m,Flow[f]==sum(VolInRd[i,r]*in(f,aMDPModel.dfOrders[i,:FlowList]) for i=1:nOrders, r in aMDPModel.dfOrders[i,:RoundList]))
    @constraint(m,Flow[f]==sum(VolInRd[i,r]* orderInFlow(aMDPModel.dfOrders[i,:FlowList],f) for i=1:nOrders, r in aMDPModel.dfOrders[i,:RoundList]))
    @constraint(m,FlowShortage[f]>=minFlow[f]-Flow[f])
    @constraint(m,FlowExcess[f]>=Flow[f]-maxFlow[f])
end

@variable(m,flowShortagecost[f in aMDPModel.dfFlows.FlowName])
@variable(m,flowExcesscost[f in aMDPModel.dfFlows.FlowName])
@variable(m,totalflowShortagecost)
@variable(m,totalflowExcesscost)

FlowShortageCostFactor=@from fl in aMDPModel.dfFlows begin
    @select fl.FlowName=>fl.FlowShortageCostFactor*fl.Active
    @collect Dict
end

FlowExcessCostFactor=@from fl in aMDPModel.dfFlows begin
    @select fl.FlowName=>fl.FlowExcessCostFactor*fl.Active
    @collect Dict
end

for f in aMDPModel.dfFlows.FlowName
    @constraint(m,flowShortagecost[f]==sum(FlowShortage[f]*FlowShortageCostFactor[f] ))
    @constraint(m,flowExcesscost[f]==sum(FlowExcess[f]*FlowExcessCostFactor[f] ))
end

@constraint(m,totalflowShortagecost==sum(flowShortagecost[f] for f in aMDPModel.dfFlows.FlowName))
@constraint(m,totalflowExcesscost==sum(flowExcesscost[f] for f in aMDPModel.dfFlows.FlowName))

@variable(m,totalSelectionCost)
@constraint(m,totalSelectionCost==sum(VolInRd[i,r]*aMDPModel.dfOrders[i,:SelectionCost] for i=1:nOrders,r in aMDPModel.dfOrders[i,:RoundList]))
@variable(m,totalOuterbayCost)
if ! haskey(aMDPModel.params,"OuterbayCost")
    @constraint(m,totalOuterbayCost==0)
else
    @constraint(m,totalOuterbayCost==sum(VolOuterBayInRd[i,r])*aMDPModel.params["OuterbayCost"])
end

@objective(m,Min,totalflowExcesscost + totalflowShortagecost + totalSelectionCost + totalOuterbayCost)
write(aMDPModel.logFile,"Solving Model\r\n")

modelPrint=open(joinpath(aMDPModel.workFolder,"model.txt"),"w")

println(modelPrint,m)
optimize!(m)
stat = termination_status(m)
write(aMDPModel.logFile,"Solved Model: $stat\r\n")
result=open(joinpath(aMDPModel.workFolder,"Result.csv"),"w")
    println(result,"Status,$stat\r")

if stat==MOI.OPTIMAL
    println(result,"Cost,",objective_value(m),"\r")
    println(result,"FlowExcessCost,",value(totalflowExcesscost),"\r")
    println(result,"FlowShortageCost,",value(totalflowShortagecost),"\r")
    println(result,"SelectionCost,",value(totalSelectionCost),"\r")
    #for r in aMDPModel.dfRounds.RoundName
    #    if value(Rd[r])>0
    #        write(result,"$r\r\n")
            #for w = 1:20
            #    val=value(RdPerWidthGroupLength[r,w])
            #    if val>0
            #        write(result," $w,$val\r\n")
            #    end
     #       end
     #   end
    #end

    close(result)

    flowResult=open(joinpath(aMDPModel.workFolder,"HSMFlows_Result.csv"),"w")
    #write Flows
    for f in aMDPModel.dfFlows.FlowName
        val =value(Flow[f])
        shortage=value(FlowShortage[f])
        excess=value(FlowExcess[f])
        write(flowResult,"$f,$val,$shortage,$excess\r\n")
    end
    close(flowResult)
    #write rounds
    roundResult=open(joinpath(aMDPModel.workFolder,"HSMRounds_Result.csv"),"w")
    for r in aMDPModel.dfRounds.RoundName
        val=value(Rd[r])
        vol=value(RdVol[r])
        write(roundResult,"$r,$val,$vol\r\n")
        #println(r,",",value(Rd[r]))
    end
    close(roundResult)
    #write Orders
    orderResult=open(joinpath(aMDPModel.workFolder,"HSMOrders_Result.csv"),"w")
    for i=1:nOrders
        vol=0
        for r in aMDPModel.dfOrders[i,:RoundList]
            if (value(VolInRd[i,r])>0)
               #println("\t$i,$r," , aMDPModel.dfOrders[i,:Works_Order_No] , ",$vol")
                vol += value(VolInRd[i,r])
            end
        end
        if vol>0
            println(orderResult,aMDPModel.dfOrders[i,:Works_Order_No],",$vol\r")
        end
    end
    close(orderResult)
end
println(result,"Cost,",999999999999999999999,"\r")
close(result)
close(aMDPModel.logFile)
return stat
end
end


text:
"""
Subset quality metrics
"""
function coverage(fileCount::Int, subset::Array{Int}, totalN::Int)
  println("Coverage metric:")
  reachableNeighbours = @parallel (vcat) for i in subset
    tempNeighbours = Array(Array{Int},0)
    for j in 1:fileCount
      index = findfirst(fileData[j][1], i)
      if index != 0
        neighbours = fileData[j][3][:, index]
        found = findfirst(neighbours, -1)
        if found != 0
          deleteat!(neighbours, found:length(neighbours))
        end
        push!(tempNeighbours, neighbours)
      end
    end
    tempNeighbours = vcat(tempNeighbours...)
  end
  C = length(unique(reachableNeighbours))/totalN
  return C
end

function coverageApprox(fileCount::Int, subset::Array{Int}, minN::Int, maxN::Int)
  println("Coverage metric:")
  totalNodes = maxN-minN+1
  reachableNeighbours = 0
  # BitArray to keep track of unique set of neighbours
  foundNeighbours = falses(totalNodes)
  @inbounds for j in 1:fileCount
    println("File $(j)")
    subsetIndices = find(x -> x in subset, fileData[j][1])'
    @inbounds for nodeIndex in subsetIndices
      nodeNeighbours = fileData[j][3][:, nodeIndex]
      @inbounds for node in nodeNeighbours
        if node == -1
          break
        end
        if !foundNeighbours[node-minN+1]
          reachableNeighbours += 1
          foundNeighbours[node-minN+1] = true
        end
      end
    end
  end
  C = reachableNeighbours / totalNodes
  return C, totalNodes
end

"""
Clustering quality metrics
"""
function modularity(fileCount::Int, k::Int, weighted::Int)
  println("Modularity metric:")
  clusterResult = fill(0.0, k, 2)
  @inbounds for i in 1:k
    println("Cluster $(i)")
    tic()
    # nodes belonging to cluster i
    clusterNodes = qData[(qData[:,2] .== i), 1]
    println("Cluster size $(length(clusterNodes))")
    clusterValues = @parallel (.+) for node in clusterNodes
      clusterAssoc = 0.0
      clusterDeg = 0.0
      foundNeighbours = Int[]
      foundClusterNeighbours = Int[]
      @inbounds for j in 1:fileCount
        indexNode = findfirst(fileData[j][1], node)
        if indexNode != 0
          # degree of new node neighbours
          nodeNeighbours = find(x -> x ∉ foundNeighbours, fileData[j][3][:, indexNode])
          if length(nodeNeighbours) != 0
            push!(foundNeighbours, fileData[j][3][nodeNeighbours, indexNode]...)
            # if unweighted
            if weighted == 0
              clusterDeg += float(length(nodeNeighbours))
            else
              clusterDeg += sum(fileData[j][4][nodeNeighbours, indexNode])
            end
          end
          # degree of new node cluster neighbours
          nodeNeighbours = find(x -> x ∉ foundClusterNeighbours, fileData[j][3][:, indexNode])
          clusterNeighbours = find(x -> x in clusterNodes, fileData[j][3][nodeNeighbours, indexNode])
          if length(clusterNeighbours) != 0
            push!(foundClusterNeighbours, fileData[j][3][nodeNeighbours, indexNode][clusterNeighbours]...)
            # if unweighted
            if weighted == 0
              clusterAssoc += float(length(clusterNeighbours))
            else
              clusterAssoc += sum(fileData[j][4][nodeNeighbours, indexNode][clusterNeighbours])
            end
          end
        end
      end
      # end parallel loop
      [clusterAssoc, clusterDeg]
    end
    clusterResult[i, :] = clusterValues
    toc()
  end
  # total sum of weights of all edges in network
  edgesSum = sum(clusterResult[:, 2])
  Q = 0.0
  @inbounds for i in 1:k
    Q += clusterResult[i, 1]/edgesSum - (clusterResult[i, 2]/edgesSum)^2
  end
  return Q
end

function modularityApprox(fileCount::Int, k::Int, weighted::Int)
  println("Modularity metric:")
  clusterResult = fill(0.0, k, 2)
  @inbounds for i in 1:k
    println("Cluster $(i)")
    # nodes belonging to cluster i
    clusterNodes = qData[(qData[:,2] .== i), 1]
    println("Cluster size $(length(clusterNodes))")
    clusterDeg, clusterAssoc = @parallel (.+) for j in 1:fileCount
      tic()
      fileDeg = 0.0
      fileAssoc = 0.0
      fileNodes = find(x -> x in clusterNodes, fileData[j][1])
      fileNeighbours = fileData[j][3][:, fileNodes]
      # sum of degree of cluster nodes' neighbours
      if weighted == 0
        fileDeg = float(length(fileNeighbours[fileNeighbours .>= 0]))
      else
        fileDeg = sum(fileData[j][4][:, fileNodes])
      end
      # sum of degree of cluster nodes' in cluster neighbours
      clusterNeighbours = find(x -> x in clusterNodes, fileNeighbours)
      if length(clusterNeighbours) != 0
        if weighted == 0
          fileAssoc = float(length(clusterNeighbours))
        else
          fileAssoc = sum(fileData[j][4][:, fileNodes][clusterNeighbours])
        end
      end
      println("File $(j) done")
      toc()
      [fileDeg, fileAssoc]
    end
    clusterResult[i, :] = [clusterDeg, clusterAssoc]
  end
  # total sum of weights of all edges in network
  edgesSum = sum(clusterResult[:, 2])
  Q = 0.0
  @inbounds for i in 1:k
    Q += clusterResult[i, 1]/edgesSum - (clusterResult[i, 2]/edgesSum)^2
  end
  return Q
end

#=
function modularityApprox(fileCount::Int, k::Int, weighted::Int)
  println("Modularity metric:")
  clusterResult = fill(0.0, k, 2)
  @inbounds for i in 1:k
    println("Cluster $(i)")
    # nodes belonging to cluster i
    clusterNodes = qData[(qData[:,2] .== i), 1]
    println("Cluster size $(length(clusterNodes))")
    @inbounds for j in 1:fileCount
      println("File $(j)")
      tic()
      fileNodes = find(x -> x in clusterNodes, fileData[j][1])
      fileNeighbours = fileData[j][3][:, fileNodes]
      # sum of degree of cluster nodes' neighbours
      if weighted == 0
        clusterDeg = float(length(fileNeighbours[fileNeighbours .>= 0]))
      else
        clusterDeg = sum(fileData[j][4][:, fileNodes])
      end
      # sum of degree of cluster nodes' in cluster neighbours
      clusterNeighbours = find(x -> x in clusterNodes, fileNeighbours)
      if length(clusterNeighbours) != 0
        if weighted == 0
          clusterAssoc = float(length(clusterNeighbours))
        else
          clusterAssoc = sum(fileData[j][4][:, fileNodes][clusterNeighbours])
        end
      end
      clusterResult[i, :] .+= [clusterDeg clusterAssoc]
      toc()
    end
  end
  # total sum of weights of all edges in network
  edgesSum = sum(clusterResult[:, 2])
  Q = 0.0
  @inbounds for i in 1:k
    Q += clusterResult[i, 1]/edgesSum - (clusterResult[i, 2]/edgesSum)^2
  end
  return Q
end
=#

function cut_conductance(fileCount::Int)
end


text:
using StatisticalRethinking, CmdStan
#gr(size=(600,600));

ProjDir = rel_path("..", "scripts", "04")
cd(ProjDir)

howell1 = CSV.read(rel_path("..", "data", "Howell1.csv"), delim=';')
df = convert(DataFrame, howell1);

df2 = filter(row -> row[:age] >= 18, df)
first(df2, 5)

heightsmodel = "
// Inferring a Rate
data {
  int N;
  real<lower=0> h[N];
}
parameters {
  real<lower=0> sigma;
  real<lower=0,upper=250> mu;
}
model {
  // Priors for mu and sigma
  mu ~ normal(178, 20);
  sigma ~ uniform( 0 , 50 );

  // Observed heights
  h ~ normal(mu, sigma);
}
";

stanmodel = Stanmodel(name="heights", monitors = ["mu", "sigma"],model=heightsmodel,
  output_format=:mcmcchains);

heightsdata = Dict("N" => length(df2[:height]), "h" => df2[:height]);

rc, chn, cnames = stan(stanmodel, heightsdata, ProjDir, diagnostics=false,
  CmdStanDir=CMDSTAN_HOME);

MCMCChains.describe(chn)

serialize("m4.1s.jls", chn)
chn2 = deserialize("m4.1s.jls")

MCMCChains.describe(chn2)

# end of m4.1s


text:
# # Robust approximate fitting
# Section 6.4.2
# Boyd & Vandenberghe "Convex Optimization"
# Original by Lieven Vandenberghe
# Adapted for Convex by Joelle Skaf - 10/03/05
#
# Adapted for Convex.jl by Karanveer Mohan and David Zeng - 26/05/14
# Original cvx code and plots here:
# <http://web.cvxr.com/cvx/examples/cvxbook/Ch06_approx_fitting/html/fig6_15.html>
#
# Consider the least-squares problem:
#       minimize $\|(A + tB)x - b\|_2$
# where $t$ is an uncertain parameter in [-1,1]
# Three approximate solutions are found:
#
#   1. nominal optimal (i.e. letting t=0)
#   2. stochastic robust approximation:
#           minimize $\mathbb{E}\|(A+tB)x - b\|_2$
#      assuming $u$ is uniformly distributed on [-1,1].
#      (reduces to minimizing $\mathbb{E} \|(A+tB)x-b\|^2 = \|A*x-b\|^2  + x^TPx$
#        where $P = \mathbb{E}(t^2) B^TB = (1/3) B^TB$ )
#   3. worst-case robust approximation:
#           minimize $\mathrm{sup}_{-1\leq u\leq 1} \|(A+tB)x - b\|_2$
#      (reduces to minimizing $\max\{\|(A-B)x - b\|_2, \|(A+B)x - b\|_2\}$ ).
#
using Convex, LinearAlgebra, SCS

# Input Data
m = 20;
n = 10;
A = randn(m, n);
(U, S, V) = svd(A);
S = diagm(exp10.(range(-1, stop = 1, length = n)));
A = U[:, 1:n] * S * V';

B = randn(m, n);
B = B / norm(B);

b = randn(m, 1);
x = Variable(n)

# Case 1: Nominal optimal solution
p = minimize(norm(A * x - b, 2))
solve!(p, SCS.Optimizer; silent_solver = true)
x_nom = evaluate(x)

# Case 2: Stochastic robust approximation
P = 1 / 3 * B' * B;
p = minimize(square(pos(norm(A * x - b))) + quadform(x, Symmetric(P)))
solve!(p, SCS.Optimizer; silent_solver = true)
x_stoch = evaluate(x)

# Case 3: Worst-case robust approximation
p = minimize(max(norm((A - B) * x - b), norm((A + B) * x - b)))
solve!(p, SCS.Optimizer; silent_solver = true)
x_wc = evaluate(x)

# Plot residuals:
parvals = range(-2, stop = 2, length = 100);

errvals(x) = [norm((A + parvals[k] * B) * x - b) for k in eachindex(parvals)]
errvals_ls = errvals(x_nom)
errvals_stoch = errvals(x_stoch)
errvals_wc = errvals(x_wc)

using Plots
plot(parvals, errvals_ls, label = "Nominal problem")
plot!(parvals, errvals_stoch, label = "Stochastic Robust Approximation")
plot!(parvals, errvals_wc, label = "Worst-Case Robust Approximation")
plot!(
    title = "Residual r(u) vs a parameter u for three approximate solutions",
    xlabel = "u",
    ylabel = "r(u) = ||A(u)x-b||_2",
)


text:
#!/cvmfs/soft.computecanada.ca/easybuild/software/2017/avx2/Compiler/gcc7.3/julia/1.3.0/bin/julia
#SBATCH --mail-user=mkschleg@ualberta.ca
#SBATCH --mail-type=ALL
#SBATCH -o reproduce_test.out # Standard output
#SBATCH -e reproduce_test.err # Standard error
#SBATCH --mem-per-cpu=512M # Memory request of 2 GB
#SBATCH --time=0:10:00 # Running time of 12 hours
#SBATCH --ntasks=8
#SBATCH --account=rrg-whitem

using Pkg
Pkg.activate(".")

using Reproduce

const save_loc = "test_exp"
const exp_file = "examples/experiment.jl"
const exp_module_name = :Main
const exp_func_name = :main_experiment

# What arguments get swept over.
arg_dict = Dict(
    "opt1"=>1:50,
    "opt2"=>[5,6,7,8]
)
# The order of the arguments in the sweep (optional)
arg_order = ["opt2", "opt1"]

static_args = Dict("steps" => 102902)

args_iterator = ArgIterator(arg_dict,
                            static_args;
                            arg_order=arg_order)

# Don't name anything in global scope `experiment`!!!
exp = Experiment(save_loc,
                 exp_file,
                 exp_module_name,
                 exp_func_name,
                 args_iterator)

pre_experiment(exp; tldr="Hello tldr")
ret = job(exp; num_workers=6, extra_args=[1])
post_experiment(exp, ret)


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# common move-up functions

function setMoveUpModule!(sim::Simulation, moveUpModule::MoveUpModule)
	sim.moveUpData.moveUpModule = moveUpModule
	sim.moveUpData.useMoveUp = moveUpModule != nullMoveUpModule
	if isdefined(sim, :backup)
		setMoveUpModule!(sim.backup, moveUpModule)
	end
end

function isAmbMovable(ambulance::Ambulance)
	return isFree(ambulance.status)
end

# for an ambulance, gives travel time for it to travel to each station
function ambMoveUpTravelTimes!(sim::Simulation, ambulance::Ambulance;
	stations::Vector{Station} = sim.stations)

	# shorthand:
	@unpack net, travel = sim
	priority = lowPriority # default travel priority for this function

	travelMode = getTravelMode!(travel, sim.responseTravelPriorities[priority], sim.time)
	(node1, time1) = getRouteNextNode!(sim, ambulance.route, travelMode.index, sim.time) # next/nearest node in ambulance route

	# get travel times to each station
	numStations = length(stations)
	ambToStationTimes = Vector{Float}(undef, numStations)
	for (i,station) in enumerate(stations)
		if ambulance.stationIndex == station.index && ambulance.status == ambIdleAtStation
			# ambulance already idle at this station, travel time is 0
			ambToStationTimes[i] = 0.0
		else
			# lookup travel time from node1 to station
			pathTravelTime = shortestPathTravelTime(net, travelMode.index, node1, station.nearestNodeIndex)
			time2 = offRoadTravelTime(travelMode, station.nearestNodeDist)
			ambToStationTimes[i] = time1 + pathTravelTime + time2 # on and off road travel times
		end
	end

	return ambToStationTimes
end

function createStationPairs(sim::Simulation, travelMode::TravelMode;
	maxPairsPerStation::Int = Inf, maxPairSeparation::Float = Inf)

	# shorthand:
	@unpack net, travel, map, stations, numStations = sim

	maxPairsPerStation = min(maxPairsPerStation, numStations)

	# calculate travel times between each pair of stations
	stationToStationTimes = zeros(Float, numStations, numStations) #Array{Float,2}(undef,numStations,numStations)
	for i = 1:numStations, j = 1:numStations
		if i != j
			# check travel time from station i to j
			travelTime = shortestPathTravelTime(net, travelMode.index, stations[i].nearestNodeIndex, stations[j].nearestNodeIndex)
			travelTime += offRoadTravelTime(travelMode, stations[i].nearestNodeDist)
			travelTime += offRoadTravelTime(travelMode, stations[j].nearestNodeDist)
			stationToStationTimes[i,j] = travelTime
		end
	end

	areStationsPaired = Array{Bool,2}(undef,numStations,numStations) # areStationsPaired[i,j] = false if stations i, j, should not be paired
	areStationsPaired[:] .= true
	sortedTimes = sort(stationToStationTimes, dims = 2)
	areStationsPaired .*= (stationToStationTimes .<= sortedTimes[:, maxPairsPerStation])
	areStationsPaired .*= (stationToStationTimes .<= maxPairSeparation)

	# make sure that areStationsPaired is symmetric
	areStationsPaired .*= areStationsPaired'

	# convert areStationsPaired to stationPairs
	stationPairs = Vector{Vector{Int}}() # stationPairs[i][k] gives kth station in pair i (k = 1,2)
	for i = 1:numStations, j = i+1:numStations
		if areStationsPaired[i,j]
			push!(stationPairs, [i,j])
		end
	end

	return stationPairs
end

function moveUpNull()
	return Ambulance[], Station[]
end


text:
# Script to test the use of cross-correlation beamforming

using Plots

using Beamforming
using Seis, Seis.Plot
using SeisTau

"Return a random selection of `n` items from `set`"
function random_selection(set, n)
    n > length(set) &&
        throw(ArgumentError("cannot select more than all the items in set"))
    out = Vector{eltype(set)}(undef, n)
    indices = collect(eachindex(set))
    for i in 1:n
        index = rand(1:length(indices))
        out[i] = set[indices[index]]
        deleteat!(indices, index)
    end
    out
end

# Plotting defaults
default(fontfamily="Helvetica")

# Choose case to look at
case = :synth
case = :real_array
case = :real_local

# Synthetic case
if case == :synth
    # Synthetic slowness of arrival
    sx, sy = 1, 2
    t = Beamforming.Synth.synthetic_arrival(sx, sy, nsta=8, lon=2rand(8), lat=2rand(8),
        distance=90)

    # Beamforming parameters
    t1, t2 = 40, 60
    window_time = 1
    s_max = 3
    ds = 0.05

    # Plotting parameters
    phases = []
# Real data example
elseif case == :real_array
    # Find a random selection of 6 traces from the UK network
    tt = sample_data(:array)
    t = random_selection(tt, 6)
    remove_trend!.(t)
    taper!.(t)
    bandpass!.(t, 0.05, 2, twopass=true)

    # Beamforming parameters
    t1, t2 = 1110, 1150
    window_time = 0.5
    s_max = 4
    ds = 0.05
    sx = sy = NaN

    phases = ["PKiKP", "PKIKP"]
elseif case == :real_local
    # Vertical stations only
    t = filter(x -> endswith(x.sta.sta, r"z"i), sample_data(:local))
    for tt in t
        remove_trend!(tt)
        taper!(tt)
        bandpass!(tt, 0.05, 1, twopass=true)
    end
    t1, t2 = 1, 10
    window_time = 2
    s_max = 20
    ds = 0.2
    sx = sy = NaN

    phases = ["p", "s"]
else
    error("choose a case to look at")
end

# Generic BF options for both CCBF and BF
method = :linear
n = 3

# Cross-correlation beamforming
bf = Beamforming.crosscorrelation_beamform(t, t1, t2,
    -s_max, s_max, -s_max, s_max, ds, t=window_time, method=method, n=n)

# Conventional beamforming for comparison
bf′ = Beamforming.beamform(t, t1, t2, s_max, ds, method=method, n=n)

# Plot parameters
colors = cgrad([:white, :blue, :green, :yellow, :orange, :red])
bf_plot_kwargs = (powscale=:dB, c=colors, clim=(-10,0),
    xlabel="X-slowness (s/°)", ylabel="Y-slowness (s/°)",
    phases=phases
    )

# Assemble plots
begin
    # Input trace plot
    p_in = section(bf.traces, title="Input traces")
    vspan!(p_in, [t1, t2], label="", fillcolor=:blue, fillalpha=0.1)

    # Cross-correlations plots
    p_cc_d = section(bf.corrs, distance_km, ylabel="Inter-station distance / km",
        title="Cross-correlations")
    vspan!(p_cc_d, [-window_time, window_time], label="", fillcolor=:blue, fillalpha=0.1)

    p_cc_az = section(bf.corrs, azimuth, ylabel="Inter-station azimuth / °",
        title="Cross-correlations")
    vspan!(p_cc_az, [-window_time, window_time], label="", fillcolor=:blue, fillalpha=0.1)

    # Map of stations with pair lines plotted
    p_map = plot([[x.evt.lon, x.sta.lon] for x in bf.corrs],
        [[x.evt.lat, x.sta.lat] for x in bf.corrs], linecolor=:black, primary=false,
        framestyle=:box, grid=false)
    scatter!(p_map, bf.traces.sta.lon, bf.traces.sta.lat, markershape=:dtriangle,
        markercolor=:blue, aspect_ratio=:equal, xlabel="Longitude", ylabel="Latitude",
        title="Array map", label="")

    # CC beamforming plot
    p_bf = plot(bf; title="CC beam power (dB)", bf_plot_kwargs...)
    if case == :synth
        scatter!(p_bf, [bf.sx_max[1] sx], [bf.sy_max[1] sy],
            label=["Recovered" "True"], markercolor=[:black :white], legend=:bottomleft)
    end

    # BF plot
    p_bf′ = plot(bf′; title="Beam power (dB)", bf_plot_kwargs...)
    if case == :synth
        scatter!(p_bf′, [sx], [sy], label="", markercolor=:black)
    end

    # Put all plots together (and display)
    plot(p_in, p_cc_d, p_cc_az, p_map, p_bf, p_bf′, layout=(2, 3), size=(1200,800))
end

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming


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

using BAT
using Test

using ArraysOfArrays, ElasticArrays

struct _SampleInfo
    x::Int
end

struct _SampleAux
    x::Float32
end

_SampleInfo() = _SampleInfo(0)
_SampleAux() = _SampleInfo(0)

@testset "density_sample" begin
    param1 = [Float64(1.),-0.1, 0.5]
    param2 = [Float64(-2.5), 0.2, 2.7]
    param4 = ones(Float64, 3)
    ds1 = @inferred DensitySample(param1, Float32(-3.3868156), 1, _SampleInfo(7), _SampleAux(0.378f0))
    ds2 = @inferred DensitySample(param2, Float32(-2.8723492), 2, _SampleInfo(8), _SampleAux(0.435f0))
    ds4 = @inferred DensitySample(param4, Float32(-4.2568156), 4, _SampleInfo(9), _SampleAux(0.612f0))

    @testset "DensitySample" begin
        @test typeof(ds1)  <: DensitySample{Float64,Float32,Int,_SampleInfo}

        @test nparams(ds1) == 3
        @test nparams.(ds1) == 3

        @test typeof(ds2)  <: DensitySample{Float64,Float32,Int,_SampleInfo,_SampleAux}

        @test ds1 != ds4
        @test ds1 != ds2
        ds3 = @inferred DensitySample(param1, Float32(-3.3868156), 1, _SampleInfo(7), _SampleAux(0.378f0))
        @test ds1 == ds3
    end

    @testset "DensitySampleVector" begin
        dsv1 = @inferred DensitySampleVector{Float64,Float32,Int,_SampleInfo,_SampleAux}(undef, 0, 3)
        @test typeof(dsv1) <: DensitySampleVector{Float64,Float32,Int,_SampleInfo,_SampleAux}

        @test size(dsv1) == (0,)
        push!(dsv1, ds1)
        @test size(dsv1) == (1,)
        @test dsv1[1] == ds1
        @test IndexStyle(dsv1) == IndexLinear()

        push!(dsv1, ds2)
        dsv2 = @inferred DensitySampleVector{Float64,Float64,Float32,_SampleInfo,_SampleAux}(undef, 0, 3)
        push!(dsv2, ds2)
        push!(dsv2, ds4)
        append!(dsv1, dsv2)
        @test dsv1[4] == ds4
        @test dsv1[2] == ds2
    end
end


text:
if VERSION < v"1.5"
    # https://github.com/JuliaLang/julia/blob/70cc57cb36d839afe6ce56ea48ff6ed01bc262c4/base/int.jl#L524-L551
    function bitrotate(x::T, k::Integer) where {T <: Integer}
        (x << ((sizeof(T) << 3 - 1) & k)) | (x >>> ((sizeof(T) << 3 - 1) & -k))
    end
end

if VERSION < v"1.1"
    isnothing(x) = x === nothing
end


text:
"""
    SimpleUnary{f,TimeAgnostic,T}

Represents a stateless unary operator that will always emit a value.

The value of the `TimeAgnostic` type parmater is coupled to [`time_agnostic`](@ref).
"""
struct SimpleUnary{f,TimeAgnostic,T} <: UnaryNodeOp{T} end

always_ticks(::SimpleUnary) = true
stateless_operator(::SimpleUnary) = true
time_agnostic(::SimpleUnary{f,false}) where {f} = false
time_agnostic(::SimpleUnary{f,true}) where {f} = true
operator!(::SimpleUnary{f,false}, t, x) where {f} = f(t, x)
operator!(::SimpleUnary{f,true}, x) where {f} = f(x)

"""
    SimpleBinary{f,TimeAgnostic,T,A}

Represents a stateless binary operator that will always emit a value.

The value of the `TimeAgnostic` type parmater is coupled to [`time_agnostic`](@ref).
"""
struct SimpleBinary{f,TimeAgnostic,T,A} <: BinaryNodeOp{T,A} end

always_ticks(::SimpleBinary) = true
stateless_operator(::SimpleBinary) = true
time_agnostic(::SimpleBinary{f,false}) where {f} = false
time_agnostic(::SimpleBinary{f,true}) where {f} = true
operator!(::SimpleBinary{f,false}, t, x, y) where {f} = f(t, x, y)
operator!(::SimpleBinary{f,true}, x, y) where {f} = f(x, y)

"""
    SimpleBinaryUnionInitial{f,TimeAgnostic,T,L,R}

Represents a stateless binary operator that will always emit a value.

The value of the `TimeAgnostic` type parmater is coupled to [`time_agnostic`](@ref).

Unlike [`SimpleBinary`](@ref), this also contains initial values for its parent nodes.
See [Initial values](@ref) for more details.
"""
struct SimpleBinaryUnionInitial{f,TimeAgnostic,T,L,R} <: BinaryNodeOp{T,UnionAlignment}
    initial_l::L
    initial_r::R
end

always_ticks(::SimpleBinaryUnionInitial) = true
stateless_operator(::SimpleBinaryUnionInitial) = true
time_agnostic(::SimpleBinaryUnionInitial{f,false}) where {f} = false
time_agnostic(::SimpleBinaryUnionInitial{f,true}) where {f} = true
operator!(::SimpleBinaryUnionInitial{f,false}, t, x, y) where {f} = f(t, x, y)
operator!(::SimpleBinaryUnionInitial{f,true}, x, y) where {f} = f(x, y)
has_initial_values(::SimpleBinaryUnionInitial) = true
initial_left(op::SimpleBinaryUnionInitial) = op.initial_l
initial_right(op::SimpleBinaryUnionInitial) = op.initial_r

"""
    SimpleBinaryLeftInitial{f,T,R}

Represents a stateless binary operator that will always emit a value.

The value of the `TimeAgnostic` type parmater is coupled to [`time_agnostic`](@ref).

Unlike [`SimpleBinary`](@ref), this also contains initial values for its right parent.
See [Initial values](@ref) for more details.
"""
struct SimpleBinaryLeftInitial{f,TimeAgnostic,T,R} <: BinaryNodeOp{T,LeftAlignment}
    initial_r::R
end

always_ticks(::SimpleBinaryLeftInitial) = true
stateless_operator(::SimpleBinaryLeftInitial) = true
time_agnostic(::SimpleBinaryLeftInitial{f,false}) where {f} = false
time_agnostic(::SimpleBinaryLeftInitial{f,true}) where {f} = true
operator!(::SimpleBinaryLeftInitial{f,false}, t, x, y) where {f} = f(t, x, y)
operator!(::SimpleBinaryLeftInitial{f,true}, x, y) where {f} = f(x, y)
has_initial_values(::SimpleBinaryLeftInitial) = true
initial_right(op::SimpleBinaryLeftInitial) = op.initial_r

"""
    SimpleNary{f,TimeAgnostic,N,T,A}

Represents a stateless `N`ary operator that will always emit a value.

The value of the `TimeAgnostic` type parmater is coupled to [`time_agnostic`](@ref).
"""
struct SimpleNary{f,TimeAgnostic,N,T,A} <: NaryNodeOp{N,T,A} end

always_ticks(::SimpleNary) = true
stateless_operator(::SimpleNary) = true
time_agnostic(::SimpleNary{f,false}) where {f} = false
time_agnostic(::SimpleNary{f,true}) where {f} = true
operator!(::SimpleNary{f,false}, t, values...) where {f} = f(t, values...)
operator!(::SimpleNary{f,true}, values...) where {f} = f(values...)

"""
    SimpleNaryInitial{f,TimeAgnostic,N,T,A,Types}

Represents a stateless binary operator that will always emit a value.

The value of the `TimeAgnostic` type parmater is coupled to [`time_agnostic`](@ref).

Unlike [`SimpleNary`](@ref), this also contains initial values.
See [Initial values](@ref) for more details.
"""
struct SimpleNaryInitial{f,TimeAgnostic,N,T,A,Types} <: NaryNodeOp{N,T,A}
    initial::Types
end

always_ticks(::SimpleNaryInitial) = true
stateless_operator(::SimpleNaryInitial) = true
time_agnostic(::SimpleNaryInitial{f,false}) where {f} = false
time_agnostic(::SimpleNaryInitial{f,true}) where {f} = true
operator!(::SimpleNaryInitial{f,false}, t, values...) where {f} = f(t, values...)
operator!(::SimpleNaryInitial{f,true}, values...) where {f} = f(values...)
has_initial_values(::SimpleNaryInitial) = true
initial_values(op::SimpleNaryInitial) = op.initial

"""
    apply(f::Function, x; out_type=nothing, time_agnostic=true)
    apply(
        f::Function, x, y[, z, ..., alignment=DEFAULT_ALIGNMENT];
        out_type=nothing, time_agnostic=true, initial_values=nothing
    )

Obtain a node with values constructed by applying the pure function `f` to the input values.

Iff `time_agnostic` is `false`, `f` will be passed the time of the current knot as the first
argument, in addition to any values.

With more than one nodal argument, alignment will be performed. In this case, the
`alignment` argument can be specified as one of [`INTERSECT`](@ref), [`LEFT`](@ref) or
[`UNION`](@ref). If unspecified, `DEFAULT_ALIGNMENT` will be used.

Internally this will infer the output type of `f` applied to the arguments, and will also
ensure that subgraph elimination occurs when possible.

If `out_type` is not specified, we attempt to infer the value type of the resulting node
automatically, using [`output_type`](@ref). Alternatively, if `out_type` is given as
anything other than `nothing`, it will be used instead.
"""
function apply(
    f::Function, x; out_type::Union{Nothing,Type}=nothing, time_agnostic::Bool=true
)
    x = _ensure_node(x)
    T = if isnothing(out_type)
        arg_types = (value_type(x),)
        arg_types = time_agnostic ? arg_types : (DateTime, arg_types...)
        output_type(f, arg_types...)
    else
        out_type
    end
    return obtain_node((x,), SimpleUnary{f,time_agnostic,T}())
end

function _get_op(f, time_agnostic, T, ::UnionAlignment, l::L, r::R) where {L,R}
    return SimpleBinaryUnionInitial{f,time_agnostic,T,L,R}(l, r)
end
function _get_op(f, time_agnostic, T, ::LeftAlignment, l, r::R) where {R}
    return SimpleBinaryLeftInitial{f,time_agnostic,T,R}(r)
end
function _get_op(f, time_agnostic, T, ::IntersectAlignment, l, r)
    return SimpleBinary{f,time_agnostic,T,IntersectAlignment}()
end

function apply(
    f::Function,
    x,
    y,
    alignment::Alignment;
    out_type::Union{Nothing,Type}=nothing,
    initial_values::Union{Nothing,Tuple{<:Any,<:Any}}=nothing,
    time_agnostic::Bool=true,
)
    x = _ensure_node(x)
    y = _ensure_node(y)
    A = typeof(alignment)
    T = if isnothing(out_type)
        arg_types = (value_type(x), value_type(y))
        arg_types = time_agnostic ? arg_types : (DateTime, arg_types...)
        output_type(f, arg_types...)
    else
        out_type
    end
    op = if isnothing(initial_values)
        SimpleBinary{f,time_agnostic,T,A}()
    else
        initial_l, initial_r = initial_values
        L = value_type(x)
        R = value_type(y)
        isa(initial_l, L) || throw(ArgumentError("$initial_l should be of type $L"))
        isa(initial_r, R) || throw(ArgumentError("$initial_r should be of type $R"))
        _get_op(f, time_agnostic, T, alignment, initial_l, initial_r)
    end
    return obtain_node((x, y), op)
end

function apply(f::Function, x::Node, y::Node; kwargs...)
    return apply(f, x, y, DEFAULT_ALIGNMENT; kwargs...)
end
apply(f::Function, x::Node, y; kwargs...) = apply(f, x, y, DEFAULT_ALIGNMENT; kwargs...)
apply(f::Function, x, y::Node; kwargs...) = apply(f, x, y, DEFAULT_ALIGNMENT; kwargs...)

function apply(
    f::Function,
    x,
    y,
    z,
    rest...;
    out_type::Union{Nothing,Type}=nothing,
    initial_values::Union{Nothing,Tuple{Vararg{Any}}}=nothing,
    time_agnostic::Bool=true,
)
    # The last argument *might* be an alignment. If it isn't, we should use the default
    # alignment.
    args = (x, y, z, rest...)
    head = args[1:(end - 1)]
    tail = last(args)
    args, alignment = if isa(tail, Alignment)
        head, tail
    else
        args, DEFAULT_ALIGNMENT
    end

    inputs = map(_ensure_node, args)
    A = typeof(alignment)
    N = length(inputs)
    input_types = map(value_type, inputs)
    T = if isnothing(out_type)
        arg_types = time_agnostic ? input_types : (DateTime, input_types...)
        output_type(f, arg_types...)
    else
        out_type
    end
    # Note that initial values should always be ignored for intersect alignment, since by
    # definition they will never be used.
    op = if A <: IntersectAlignment || isnothing(initial_values)
        SimpleNary{f,time_agnostic,N,T,A}()
    else
        # Sanity check the initial values.
        Types = Tuple{input_types...}
        if !isa(initial_values, Types)
            throw(ArgumentError("$initial_values should be of type $Types"))
        end
        SimpleNaryInitial{f,time_agnostic,N,T,A,Types}(initial_values)
    end
    return obtain_node(inputs, op)
end

"""
    BCast(f::Function)

Represent a function which should be broadcasted when called.

That is, `BCast{f}(x, y...)` is identical to `f.(x, y...)`.
"""
struct BCast{f} <: Function
    BCast(f::Function) = new{f}()
end
@inline (::BCast{f})(args...; kwargs...) where {f} = f.(args...; kwargs...)

"""
    Wrapped{f}

Represent a function that, when called, will expect its arguments to be nodes and try to
convert them as such.
"""
struct Wrapped{f} <: Function
    time_agnostic::Bool
end
@inline function (wrapped::Wrapped{f})(x, rest...; kwargs...) where {f}
    return apply(f, x, rest...; time_agnostic=wrapped.time_agnostic, kwargs...)
end

# TODO Would we like an identity map? Not important right now, but could be in the future
#   if we allow stateful things?

"""
    wrap(f::Function; time_agnostic=true)

Return a callable object that acts on nodes, and returns a node.

It is assumed that `f` is stateless (this will therefore not work with closures). We also
assume that we will always emit a knot when the alignment semantics say we should — thus `f`
must always return a valid output value.

Iff `time_agnostic` is `false`, `f` will be passed the time of the current knot as the first
argument, followed by the value of every input.

If the object is called with more than one node, alignment will be performed.
If an alignment other than the default should be used, provide it as the final argument.

Internally this will call `TimeDag.apply(f, args...; kwargs...)`; see there for further
details.
"""
wrap(f::Function; time_agnostic::Bool=true) = Wrapped{f}(time_agnostic)

"""
    wrapb(f::Function; time_agnostic=true)

`wrapb` is like [`wrap`](@ref), however `f` will be broadcasted over all input values.
"""
wrapb(f::Function; time_agnostic::Bool=true) = Wrapped{BCast(f)}(time_agnostic)

"""
    @simple_unary(f)

Define a method `f(::Node)` that will obtain the correct instance of `SimpleUnary{f}`.

This will internally infer the output value, and perform subgraph elimination.
"""
macro simple_unary(f, self_inverse::Bool=false)
    f = esc(f)
    return quote
        """
            $($f)(x::Node)

        Obtain a node with values constructed by applying `$($f)` to each input value.
        """
        $f(x::Node) = apply($f, x)
    end
end

"""
    @simple_unary_self_inverse(f)

Define a method `f(::Node)` that will obtain the correct instance of `SimpleUnary{f}`.

This must ONLY be used if `f(f(x)) == x` for all nodes `x`.

This will internally infer the output value, and perform subgraph elimination.
"""
macro simple_unary_self_inverse(f)
    f = esc(f)
    return quote
        """
            $($f)(x::Node)

        Obtain a node with values constructed by applying `$($f)` to each input value.
        """
        function $f(x::Node)
            # Optimisation: self-inverse
            isa(x.op, SimpleUnary{$f}) && return only(parents(x))
            return apply($f, x)
        end
    end
end

"""
    @simple_binary(f)

Define methods `f(x, y)` that will obtain the correct instance of `SimpleBinary{f}`.

These will internally infer the output value, and perform subgraph elimination.
"""
macro simple_binary(f)
    f = esc(f)
    return quote
        """
            $($f)(x, y[, alignment=DEFAULT_ALIGNMENT; kwargs...])

        Obtain a node with values constructed by applying `$($f)` to the input values.

        An `alignment` can optionally be specified. `x` and `y` should be nodes, or
        constants that can be converted to nodes.

        Other keyword arguments are passed to [`apply`](@ref).
        """
        $f(x, y, alignment::Alignment; kwargs...) = apply($f, x, y, alignment; kwargs...)
        $f(x::Node, y::Node; kwargs...) = $f(x, y, DEFAULT_ALIGNMENT; kwargs...)
        $f(x::Node, y; kwargs...) = $f(x, y, DEFAULT_ALIGNMENT; kwargs...)
        $f(x, y::Node; kwargs...) = $f(x, y, DEFAULT_ALIGNMENT; kwargs...)
    end
end


text:
using DelimitedFiles
using Sundials
using SteadyStateDiffEq
import BioMASS: isinstalled

@testset "Bifurcation analysis" begin
    for model in ["restriction_point", "g1s_transition", "mitotic_control"]
        @testset "$model" begin
            MODEL_PATH = "../examples/bifurcation/" * model
            create_diffeq(MODEL_PATH)
            include(MODEL_PATH * "/diagram.jl")
            fp, br = calc_fixed_point_vec(MODEL_PATH)
            for file in ["/data/fp.dat", "/data/ev.dat"]
                @test isfile(MODEL_PATH * file)
            end
            if isinstalled("matplotlib")
                using PyPlot
                bifurcation_diagram(MODEL_PATH, fp, br)
                @test isfile(MODEL_PATH * "/bifurcation_diagram.pdf")
            end
            rm(MODEL_PATH * "/forwarddiff.jl")
            rm(MODEL_PATH * "/data", recursive=true, force=true)
        end
    end
end


text:
"""
Bridges `CP.AllDifferentExceptConstants` to reifications.
"""
struct AllDifferentExceptConstants2ReificationBridge{T} <: MOIBC.AbstractBridge
    vars_compare::Matrix{MOI.VariableIndex} # First index: value in the array; second index: excepted variable
    vars_compare_bin::Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}
    cons_compare_reif::Matrix{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}}}
    # An upper-triangular matrix (i.e. nothing if i < j, constraint if i >= j).
    # Standard sparse matrices cannot store anything that is not a Number
    # (more specifically, anything that does not implement `zero(T)`).
    # https://github.com/JuliaLang/julia/issues/30573
    vars_different::Dict{Tuple{Int, Int}, MOI.VariableIndex}
    vars_different_bin::Dict{
        Tuple{Int, Int},
        MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne},
    }
    cons_different_reif::Dict{
        Tuple{Int, Int},
        MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Reification{CP.DifferentFrom{T}}},
    }
    cons::Dict{
        Tuple{Int, Int},
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}},
    }
end

function MOIBC.bridge_constraint(
    ::Type{AllDifferentExceptConstants2ReificationBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.AllDifferentExceptConstants{T},
) where {T}
    return MOIBC.bridge_constraint(
        AllDifferentExceptConstants2ReificationBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{AllDifferentExceptConstants2ReificationBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.AllDifferentExceptConstants{T},
) where {T}
    f_scalars = MOIU.scalarize(f)
    # dim = MOI.output_dimension(f) # TODO: remove everywhere

    values = collect(s.k)

    vars_compare = Matrix{MOI.VariableIndex}(undef, s.dimension, length(values))
    vars_compare_bin = Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}(undef, s.dimension, length(values))
    for i in 1:s.dimension
        vars_compare[i, :], vars_compare_bin[i, :] = MOI.add_constrained_variables(model, [MOI.ZeroOne() for _ in 1:length(values)])
    end

    cons_compare_reif = Matrix{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}}}(undef, s.dimension, length(values))
    for i in 1:s.dimension
        for j in 1:length(values)
            cons_compare_reif[i, j] = MOI.add_constraint(
                model,
                MOIU.vectorize(
                    [
                        one(T) * vars_compare[i, j],
                        f_scalars[i] - T(values[j]),
                    ]
                ),
                CP.Reification(MOI.EqualTo(zero(T))),
            )
        end
    end

    # Upper-triangular matrix of constraints: i >= j, i.e. d(d-1)/2 elements:
    #     \sum_{i=2}^{d} (n - i + 1) = d (d - 1) / 2
    vars_different = Dict{Tuple{Int, Int}, MOI.VariableIndex}()
    vars_different_bin = Dict{Tuple{Int, Int}, MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}()
    cons_different_reif = Dict{Tuple{Int, Int}, MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Reification{CP.DifferentFrom{T}}}}()
    cons = Dict{Tuple{Int, Int}, MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}}()

    sizehint!(vars_different, s.dimension * (s.dimension - 1) / 2)
    sizehint!(vars_different_bin, s.dimension * (s.dimension - 1) / 2)
    sizehint!(cons_different_reif, s.dimension * (s.dimension - 1) / 2)
    sizehint!(cons, s.dimension * (s.dimension - 1) / 2)

    for i in 1:s.dimension
        for j in (i+1):s.dimension
            vars_different[i, j], vars_different_bin[i, j] = MOI.add_constrained_variable(model, MOI.ZeroOne())

            cons_different_reif[i, j] = MOI.add_constraint(
                model,
                MOIU.vectorize(
                    [
                        one(T) * vars_different[i, j],
                        f_scalars[i] - f_scalars[j],
                    ]
                ),
                CP.Reification(CP.DifferentFrom(zero(T))),
            )

            cons[i, j] = MOI.add_constraint(
                model,
                sum(one(T) .* vars_compare[i, :]) + sum(one(T) .* vars_compare[j, :]) + vars_different[i, j],
                MOI.GreaterThan(one(T)),
            )
        end
    end

    return AllDifferentExceptConstants2ReificationBridge{T}(vars_compare, vars_compare_bin, cons_compare_reif, vars_different, vars_different_bin, cons_different_reif, cons)
end

function MOI.supports_constraint(
    ::Type{AllDifferentExceptConstants2ReificationBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.AllDifferentExceptConstants{T}},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{AllDifferentExceptConstants2ReificationBridge{T}}) where {T}
    return [(MOI.ZeroOne,)]
end

function MOIB.added_constraint_types(::Type{AllDifferentExceptConstants2ReificationBridge{T}}) where {T}
    return [
        (MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}),
        (MOI.ScalarAffineFunction{T}, CP.Reification{CP.DifferentFrom{T}}),
        (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}),
    ]
end

function MOI.get(b::AllDifferentExceptConstants2ReificationBridge, ::MOI.NumberOfVariables)
    return length(b.vars_compare) + length(b.vars_different)
end

function MOI.get(
    b::AllDifferentExceptConstants2ReificationBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return length(b.vars_compare_bin) + length(b.vars_different_bin)
end

function MOI.get(
    b::AllDifferentExceptConstants2ReificationBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}},
    },
) where {T}
    return length(b.cons_compare_reif)
end

function MOI.get(
    b::AllDifferentExceptConstants2ReificationBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, CP.Reification{CP.DifferentFrom{T}},
    },
) where {T}
    return length(b.cons_different_reif)
end

function MOI.get(
    b::AllDifferentExceptConstants2ReificationBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T},
    },
) where {T}
    return length(b.cons)
end

function MOI.get(
    b::AllDifferentExceptConstants2ReificationBridge{T},
    ::MOI.ListOfVariableIndices,
) where {T}
    return vcat(vec(collect(values(b.vars_compare))), vec(collect(values(b.vars_different))))
end

function MOI.get(
    b::AllDifferentExceptConstants2ReificationBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return vcat(vec(collect(values(b.vars_compare_bin))), vec(collect(values(b.vars_different_bin))))
end

function MOI.get(
    b::AllDifferentExceptConstants2ReificationBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}},
    },
) where {T}
    return vec(b.cons_compare_reif)
end

function MOI.get(
    b::AllDifferentExceptConstants2ReificationBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, CP.Reification{CP.DifferentFrom{T}},
    },
) where {T}
    return collect(values(b.cons_different_reif))
end

function MOI.get(
    b::AllDifferentExceptConstants2ReificationBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T},
    },
) where {T}
    return collect(values(b.cons))
end