text:
"""
    GCNConv([fg,] in => out, σ=identity; bias=true, init=glorot_uniform)

Graph convolutional layer.

# Arguments

- `fg`: Optionally pass a [`FeaturedGraph`](@ref).
- `in`: The dimension of input features.
- `out`: The dimension of output features.
- `σ`: Activation function.
- `bias`: Add learnable bias.
- `init`: Weights' initializer.

The input to the layer is a node feature array `X`
of size `(num_features, num_nodes)`.
"""
struct GCNConv{A<:AbstractMatrix, B, F, S<:AbstractFeaturedGraph} <: AbstractGraphLayer
    weight::A
    bias::B
    σ::F
    fg::S
end

function GCNConv(fg::AbstractFeaturedGraph, ch::Pair{Int,Int}, σ=identity;
                 init=glorot_uniform, bias::Bool=true)
    in, out = ch
    W = init(out, in)
    b = Flux.create_bias(W, bias, out)
    GCNConv(W, b, σ, fg)
end

GCNConv(ch::Pair{Int,Int}, σ = identity; kwargs...) =
    GCNConv(NullGraph(), ch, σ; kwargs...)

@functor GCNConv

function (l::GCNConv)(fg::FeaturedGraph, x::AbstractMatrix)
    L̃ = normalized_laplacian(fg, eltype(x); selfloop=true)
    l.σ.(l.weight * x * L̃ .+ l.bias)
end

(l::GCNConv)(fg::FeaturedGraph) = FeaturedGraph(fg, nf = l(fg, node_feature(fg)))

function Base.show(io::IO, l::GCNConv)
    out, in = size(l.weight)
    print(io, "GCNConv($in => $out")
    l.σ == identity || print(io, ", ", l.σ)
    print(io, ")")
end

"""
    ChebConv([fg,] in=>out, k; bias=true, init=glorot_uniform)

Chebyshev spectral graph convolutional layer.

# Arguments

- `fg`: Optionally pass a [`FeaturedGraph`](@ref).
- `in`: The dimension of input features.
- `out`: The dimension of output features.
- `k`: The order of Chebyshev polynomial.
- `bias`: Add learnable bias.
- `init`: Weights' initializer.
"""
struct ChebConv{A<:AbstractArray{<:Number,3}, B, S<:AbstractFeaturedGraph} <: AbstractGraphLayer
    weight::A
    bias::B
    fg::S
    k::Int
end

function ChebConv(fg::AbstractFeaturedGraph, ch::Pair{Int,Int}, k::Int;
                  init=glorot_uniform, bias::Bool=true)
    in, out = ch
    W = init(out, in, k)
    b = Flux.create_bias(W, bias, out)
    ChebConv(W, b, fg, k)
end

ChebConv(ch::Pair{Int,Int}, k::Int; kwargs...) =
    ChebConv(NullGraph(), ch, k; kwargs...)

@functor ChebConv

function (c::ChebConv)(fg::FeaturedGraph, X::AbstractMatrix{T}) where T
    check_num_nodes(fg, X)
    @assert size(X, 1) == size(c.weight, 2) "Input feature size must match input channel size."

    L̃ = scaled_laplacian(fg, eltype(X))

    Z_prev = X
    Z = X * L̃
    Y = view(c.weight,:,:,1) * Z_prev
    Y += view(c.weight,:,:,2) * Z
    for k = 3:c.k
        Z, Z_prev = 2*Z*L̃ - Z_prev, Z
        Y += view(c.weight,:,:,k) * Z
    end
    return Y .+ c.bias
end

(l::ChebConv)(fg::FeaturedGraph) = FeaturedGraph(fg, nf = l(fg, node_feature(fg)))

function Base.show(io::IO, l::ChebConv)
    out, in, k = size(l.weight)
    print(io, "ChebConv(", in, " => ", out)
    print(io, ", k=", k)
    print(io, ")")
end

"""
    GraphConv([fg,] in => out, σ=identity, aggr=+; bias=true, init=glorot_uniform)

Graph neural network layer.

# Arguments

- `fg`: Optionally pass a [`FeaturedGraph`](@ref).
- `in`: The dimension of input features.
- `out`: The dimension of output features.
- `σ`: Activation function.
- `aggr`: An aggregate function applied to the result of message function. `+`, `-`,
`*`, `/`, `max`, `min` and `mean` are available.
- `bias`: Add learnable bias.
- `init`: Weights' initializer.
"""
struct GraphConv{V<:AbstractFeaturedGraph, A<:AbstractMatrix, B} <: MessagePassing
    fg::V
    weight1::A
    weight2::A
    bias::B
    σ
    aggr
end

function GraphConv(fg::AbstractFeaturedGraph, ch::Pair{Int,Int}, σ=identity, aggr=+;
                   init=glorot_uniform, bias::Bool=true)
    in, out = ch
    W1 = init(out, in)
    W2 = init(out, in)
    b = Flux.create_bias(W1, bias, out)
    GraphConv(fg, W1, W2, b, σ, aggr)
end

GraphConv(ch::Pair{Int,Int}, σ=identity, aggr=+; kwargs...) =
    GraphConv(NullGraph(), ch, σ, aggr; kwargs...)

@functor GraphConv

message(gc::GraphConv, x_i, x_j::AbstractVector, e_ij) = gc.weight2 * x_j

update(gc::GraphConv, m::AbstractVector, x::AbstractVector) = gc.σ.(gc.weight1*x .+ m .+ gc.bias)

function (gc::GraphConv)(fg::FeaturedGraph, x::AbstractMatrix)
    check_num_nodes(fg, x)
    _, x = propagate(gc, adjacency_list(fg), Fill(0.f0, 0, ne(fg)), x, +)
    x
end

(l::GraphConv)(fg::FeaturedGraph) = FeaturedGraph(fg, nf = l(fg, node_feature(fg)))

function Base.show(io::IO, l::GraphConv)
    in_channel = size(l.weight1, ndims(l.weight1))
    out_channel = size(l.weight1, ndims(l.weight1)-1)
    print(io, "GraphConv(", in_channel, " => ", out_channel)
    l.σ == identity || print(io, ", ", l.σ)
    print(io, ", aggr=", l.aggr)
    print(io, ")")
end

"""
    GATConv([fg,] in => out;
            heads=1,
            concat=true,
            init=glorot_uniform
            bias=true,
            negative_slope=0.2)

Graph attentional layer.

# Arguments

- `fg`: Optionally pass a [`FeaturedGraph`](@ref).
- `in`: The dimension of input features.
- `out`: The dimension of output features.
- `bias::Bool`: Keyword argument, whether to learn the additive bias.
- `heads`: Number attention heads
- `concat`: Concatenate layer output or not. If not, layer output is averaged.
- `negative_slope::Real`: Keyword argument, the parameter of LeakyReLU.
"""
struct GATConv{V<:AbstractFeaturedGraph, T, A<:AbstractMatrix{T}, B} <: MessagePassing
    fg::V
    weight::A
    bias::B
    a::A
    negative_slope::T
    channel::Pair{Int, Int}
    heads::Int
    concat::Bool
end

function GATConv(fg::AbstractFeaturedGraph, ch::Pair{Int,Int};
                 heads::Int=1, concat::Bool=true, negative_slope=0.2f0,
                 init=glorot_uniform, bias::Bool=true)
    in, out = ch
    W = init(out*heads, in)
    b = Flux.create_bias(W, bias, out*heads)
    a = init(2*out, heads)
    GATConv(fg, W, b, a, negative_slope, ch, heads, concat)
end

GATConv(ch::Pair{Int,Int}; kwargs...) = GATConv(NullGraph(), ch; kwargs...)

@functor GATConv

# Here the α that has not been softmaxed is the first number of the output message
function message(gat::GATConv, x_i::AbstractVector, x_j::AbstractVector)
    x_i = reshape(gat.weight*x_i, :, gat.heads)
    x_j = reshape(gat.weight*x_j, :, gat.heads)
    x_ij = vcat(x_i, x_j+zero(x_j))
    e = sum(x_ij .* gat.a, dims=1)  # inner product for each head, output shape: (1, gat.heads)
    e_ij = leakyrelu.(e, gat.negative_slope)
    vcat(e_ij, x_j)  # shape: (n+1, gat.heads)
end

# After some reshaping due to the multihead, we get the α from each message,
# then get the softmax over every α, and eventually multiply the message by α
function apply_batch_message(gat::GATConv, i, js, X::AbstractMatrix)
    e_ij = mapreduce(j -> GeometricFlux.message(gat, _view(X, i), _view(X, j)), hcat, js)
    n = size(e_ij, 1)
    αs = Flux.softmax(reshape(view(e_ij, 1, :), gat.heads, :), dims=2)
    msgs = view(e_ij, 2:n, :) .* reshape(αs, 1, :)
    reshape(msgs, (n-1)*gat.heads, :)
end

update_batch_edge(gat::GATConv, adj, E::AbstractMatrix, X::AbstractMatrix, u) = update_batch_edge(gat, adj, X)

function update_batch_edge(gat::GATConv, adj, X::AbstractMatrix)
    n = size(adj, 1)
    # a vertex must always receive a message from itself
    Zygote.ignore() do
        GraphLaplacians.add_self_loop!(adj, n)
    end
    mapreduce(i -> apply_batch_message(gat, i, adj[i], X), hcat, 1:n)
end

# The same as update function in batch manner
update_batch_vertex(gat::GATConv, M::AbstractMatrix, X::AbstractMatrix, u) = update_batch_vertex(gat, M)

function update_batch_vertex(gat::GATConv, M::AbstractMatrix)
    M = M .+ gat.bias
    if !gat.concat
        N = size(M, 2)
        M = reshape(mean(reshape(M, :, gat.heads, N), dims=2), :, N)
    end
    return M
end

function (gat::GATConv)(fg::FeaturedGraph, X::AbstractMatrix)
    check_num_nodes(fg, X)
    _, X = propagate(gat, adjacency_list(fg), Fill(0.f0, 0, ne(fg)), X, +)
    X
end

(l::GATConv)(fg::FeaturedGraph) = FeaturedGraph(fg, nf = l(fg, node_feature(fg)))

function Base.show(io::IO, l::GATConv)
    in_channel = size(l.weight, ndims(l.weight))
    out_channel = size(l.weight, ndims(l.weight)-1)
    print(io, "GATConv(", in_channel, "=>", out_channel)
    print(io, ", LeakyReLU(λ=", l.negative_slope)
    print(io, "))")
end

"""
    GatedGraphConv([fg,] out, num_layers; aggr=+, init=glorot_uniform)

Gated graph convolution layer.

# Arguments

- `fg`: Optionally pass a [`FeaturedGraph`](@ref).
- `out`: The dimension of output features.
- `num_layers`: The number of gated recurrent unit.
- `aggr`: An aggregate function applied to the result of message function. `+`, `-`,
`*`, `/`, `max`, `min` and `mean` are available.
"""
struct GatedGraphConv{V<:AbstractFeaturedGraph, A<:AbstractArray{<:Number,3}, R} <: MessagePassing
    fg::V
    weight::A
    gru::R
    out_ch::Int
    num_layers::Int
    aggr
end

function GatedGraphConv(fg::AbstractFeaturedGraph, out_ch::Int, num_layers::Int;
                        aggr=+, init=glorot_uniform)
    w = init(out_ch, out_ch, num_layers)
    gru = GRUCell(out_ch, out_ch)
    GatedGraphConv(fg, w, gru, out_ch, num_layers, aggr)
end

GatedGraphConv(out_ch::Int, num_layers::Int; kwargs...) =
    GatedGraphConv(NullGraph(), out_ch, num_layers; kwargs...)

@functor GatedGraphConv

message(ggc::GatedGraphConv, x_i, x_j::AbstractVector, e_ij) = x_j

update(ggc::GatedGraphConv, m::AbstractVector, x) = m

function (ggc::GatedGraphConv)(fg::FeaturedGraph, H::AbstractMatrix{S}) where {T<:AbstractVector,S<:Real}
    check_num_nodes(fg, H)
    m, n = size(H)
    @assert (m <= ggc.out_ch) "number of input features must less or equals to output features."
    adj = adjacency_list(fg)
    if m < ggc.out_ch
        Hpad = similar(H, S, ggc.out_ch - m, n)
        H = vcat(H, fill!(Hpad, 0))
    end
    for i = 1:ggc.num_layers
        M = view(ggc.weight, :, :, i) * H
        _, M = propagate(ggc, adj, Fill(0.f0, 0, ne(fg)), M, +)
        H, _ = ggc.gru(H, M)  # BUG: FluxML/Flux.jl#1381
    end
    H
end

(l::GatedGraphConv)(fg::FeaturedGraph) = FeaturedGraph(fg, nf = l(fg, node_feature(fg)))

function Base.show(io::IO, l::GatedGraphConv)
    print(io, "GatedGraphConv(($(l.out_ch) => $(l.out_ch))^$(l.num_layers)")
    print(io, ", aggr=", l.aggr)
    print(io, ")")
end

"""
    EdgeConv([fg,] nn; aggr=max)

Edge convolutional layer.

# Arguments

- `fg`: Optionally pass a [`FeaturedGraph`](@ref).
- `nn`: A neural network (e.g. a Dense layer or a MLP).
- `aggr`: An aggregate function applied to the result of message function. `+`, `max` and `mean` are available.
"""
struct EdgeConv{V<:AbstractFeaturedGraph} <: MessagePassing
    fg::V
    nn
    aggr
end

EdgeConv(fg::AbstractFeaturedGraph, nn; aggr=max) = EdgeConv(fg, nn, aggr)
EdgeConv(nn; kwargs...) = EdgeConv(NullGraph(), nn; kwargs...)

@functor EdgeConv

message(ec::EdgeConv, x_i::AbstractVector, x_j::AbstractVector, e_ij) = ec.nn(vcat(x_i, x_j .- x_i))
update(ec::EdgeConv, m::AbstractVector, x) = m

function (ec::EdgeConv)(fg::FeaturedGraph, X::AbstractMatrix)
    check_num_nodes(fg, X)
    _, X = propagate(ec, adjacency_list(fg), Fill(0.f0, 0, ne(fg)), X, ec.aggr)
    X
end

(l::EdgeConv)(fg::FeaturedGraph) = FeaturedGraph(fg, nf = l(fg, node_feature(fg)))

function Base.show(io::IO, l::EdgeConv)
    print(io, "EdgeConv(", l.nn)
    print(io, ", aggr=", l.aggr)
    print(io, ")")
end

"""
    GINConv([fg,] nn, [eps=0])

    Graph Isomorphism Network.

# Arguments

- `fg`: Optionally pass in a FeaturedGraph as input.
- `nn`: A neural network/layer.
- `eps`: Weighting factor.

The definition of this is as defined in the original paper,
Xu et. al. (2018) https://arxiv.org/abs/1810.00826.
"""
struct GINConv{G,R} <: MessagePassing
    fg::G
    nn
    eps::R

    function GINConv(fg::G, nn, eps::R=0f0) where {G<:AbstractFeaturedGraph,R<:Real}
        new{G,R}(fg, nn, eps)
    end
end

function GINConv(nn, eps::Real=0f0)
    GINConv(NullGraph(), nn, eps)
end

Flux.trainable(g::GINConv) = (fg=g.fg, nn=g.nn)

message(g::GINConv, x_i::AbstractVector, x_j::AbstractVector) = x_j
update(g::GINConv, m::AbstractVector, x) = g.nn((1 + g.eps) * x + m)

@functor GINConv

function (g::GINConv)(fg::FeaturedGraph, X::AbstractMatrix)
    gf = graph(fg)
    GraphSignals.check_num_node(gf, X)
    _, X = propagate(g, adjacency_list(gf), Fill(0.f0, 0, ne(gf)), X, +)
    X
end

(l::GINConv)(fg::FeaturedGraph) = FeaturedGraph(fg.graph, nf = l(fg, node_feature(fg)))

"""
    CGConv([fg,] (node_dim, edge_dim), out, init, bias=true, as_edge=false)

Crystal Graph Convolutional network. Uses both node and edge features.

# Arguments

- `fg`: Optional [`FeaturedGraph`] argument(@ref)
- `node_dim`: Dimensionality of the input node features. Also is necessarily the output dimensionality.
- `edge_dim`: Dimensionality of the input edge features.
- `out`: Dimensionality of the output features.
- `init`: Initialization algorithm for each of the weight matrices
- `bias`: Whether or not to learn an additive bias parameter.
- `as_edge`: When call to layer `CGConv(M)`, accept input feature as node features or edge features.

# Usage

You can call `CGConv` in several different ways:

- Pass a FeaturedGraph: `CGConv(fg)`, returns `FeaturedGraph`
- Pass both node and edge features: `CGConv(X, E)`
- Pass one matrix, which is determined as node features or edge features by `as_edge` keyword argument.
"""
struct CGConv{E, V<:AbstractFeaturedGraph, A<:AbstractMatrix, B} <: MessagePassing
    fg::V
    Wf::A
    Ws::A
    bf::B
    bs::B
end

@functor CGConv

function CGConv(fg::G, dims::NTuple{2,Int};
                init=glorot_uniform, bias=true, as_edge=false) where {G<:AbstractFeaturedGraph}
    node_dim, edge_dim = dims
    Wf = init(node_dim, 2*node_dim + edge_dim)
    Ws = init(node_dim, 2*node_dim + edge_dim)
    bf = Flux.create_bias(Wf, bias, node_dim)
    bs = Flux.create_bias(Ws, bias, node_dim)
    T, S = typeof(Wf), typeof(bf)

    CGConv{as_edge,G,T,S}(fg, Wf, Ws, bf, bs)
end

function CGConv(dims::NTuple{2,Int}; init=glorot_uniform, bias=true, as_edge=false)
    CGConv(NullGraph(), dims; init=init, bias=bias, as_edge=as_edge)
end

message(c::CGConv,
        x_i::AbstractVector, x_j::AbstractVector, e::AbstractVector) = begin
    z = vcat(x_i, x_j, e)
    σ.(c.Wf * z + c.bf) .* softplus.(c.Ws * z + c.bs)
end
update(c::CGConv, m::AbstractVector, x) = x + m

function (c::CGConv)(fg::FeaturedGraph, X::AbstractMatrix, E::AbstractMatrix)
    check_num_nodes(fg, X)
    check_num_edges(fg, E)
    _, Y = propagate(c, adjacency_list(fg), E, X, +)
    Y
end

(l::CGConv)(fg::FeaturedGraph) = FeaturedGraph(fg, nf=l(fg, node_feature(fg),
                                                        edge_feature(fg)),
                                               ef=edge_feature(fg))

(l::CGConv)(X::AbstractMatrix, E::AbstractMatrix) = l(l.fg, X, E)

(l::CGConv{true})(M::AbstractMatrix) = l(l.fg, node_feature(l.fg), M)
(l::CGConv{false})(M::AbstractMatrix) = l(l.fg, M, edge_feature(l.fg))


text:
#include("C:/Users/10500508/Documents/julia/MDP_HSM/MDP_HSM_Model.jl")
using MDP_HSM
#global_logger(SimpleLogger(open(joinpath(path, "Logging.txt"),"w")))

if size(ARGS,1)==0
    myModel=MDP_HSM.MDP_HSM_Model("C:/tmp/20180801085842")
elseif size(ARGS,1)==1
    arg=ARGS[1]
    #error(arg)
    myModel=MDP_HSM.MDP_HSM_Model(arg)
else
    nArgs=size(ARGS,1)
    for i=1:nArgs
        println(ARGS[i])
    end
    error("0 or 1 arguments expected, received $nArgs")
end

MDP_HSM.RunModel(myModel) ==:Optimal #|| MDP_HSM.RunModel(myModel,false)
close(myModel.logFile)

ext: jl


text:
using GoogleCloud
using Base.Test
@test 1 == 1


text:
"""
    changed_traverse_strategy!(com::CS.CoM, old_traverse_strategy)

Changing from one traverse strategy to another. This means that the priority queue
needs to be rebuilt from scratch.
"""
function changed_traverse_strategy!(com::CS.CoM, old_traverse_strategy)
    old_backtrack_pq = deepcopy(com.backtrack_pq)
    com.backtrack_pq = PriorityQueue{Int,Priority}(Base.Order.Reverse)
    if old_traverse_strategy == Val(:DFS) && com.traverse_strategy == Val(:BFS)
        for elem in old_backtrack_pq
            idx = elem.first
            priority = elem.second
            com.backtrack_pq[idx] =
                PriorityBFS(priority.bound, priority.depth, priority.neg_idx)
        end
    elseif old_traverse_strategy == Val(:BFS) && com.traverse_strategy == Val(:DFS)
        for elem in old_backtrack_pq
            idx = elem.first
            priority = elem.second
            com.backtrack_pq[idx] =
                PriorityDFS(priority.depth, priority.bound, priority.neg_idx)
        end
    end
end

"""
    set_update_backtrack_pq!(com::CS.CoM, backtrack_obj::BacktrackObj; best_bound=backtrack_obj.best_bound)

Inserts or updates the value for `backtrack_obj` in the priority queue: `backtrack_pq`.
It uses the best bound of the object by default but can be overwritten with `; best_bound=other`
"""
function set_update_backtrack_pq!(
    com::CS.CoM,
    backtrack_obj::BacktrackObj;
    best_bound = backtrack_obj.best_bound,
)
    backtrack_obj.status == :Closed && return
    if com.traverse_strategy == Val(:DFS)
        if com.sense == MOI.MIN_SENSE
            com.backtrack_pq[backtrack_obj.idx] =
                PriorityDFS(backtrack_obj.depth, -best_bound, -backtrack_obj.idx)
        else
            com.backtrack_pq[backtrack_obj.idx] =
                PriorityDFS(backtrack_obj.depth, best_bound, -backtrack_obj.idx)
        end
    else
        if com.sense == MOI.MIN_SENSE
            com.backtrack_pq[backtrack_obj.idx] =
                PriorityBFS(-best_bound, backtrack_obj.depth, -backtrack_obj.idx)
        else
            com.backtrack_pq[backtrack_obj.idx] =
                PriorityBFS(best_bound, backtrack_obj.depth, -backtrack_obj.idx)
        end
    end
end

"""
    add2priorityqueue(com::CS.CoM, backtrack_obj::BacktrackObj)

Add the backtrack_obj to the priority queue `backtrack_pq`
"""
function add2priorityqueue(com::CS.CoM, backtrack_obj::BacktrackObj)
    set_update_backtrack_pq!(com, backtrack_obj)
end

"""
    close_node!(com::CS.CoM, node_idx::Int)

Close a backtrack object node if not already closed by setting the status to `:Closed`
and deleting it from the priority queue `backtrack_pq`
"""
function close_node!(com::CS.CoM, node_idx::Int)
    backtrack_vec = com.backtrack_vec
    backtrack_vec[node_idx].status == :Closed && return
    backtrack_vec[node_idx].status = :Closed

    delete!(com.backtrack_pq, node_idx)
end

function get_first_open(backtrack_vec)
    found = false
    l = length(backtrack_vec)
    backtrack_obj = backtrack_vec[l]
    while backtrack_obj.status != :Open
        l -= 1
        if l == 0
            break
        end
        backtrack_obj = backtrack_vec[l]
    end
    if l != 0
        found = true
    end
    return found, backtrack_obj
end

function get_next_node(
    com::CS.CoM,
    backtrack_vec::Vector{BacktrackObj{T}},
    sorting,
) where {T<:Real}
    found, backtrack_obj = get_next_node(com, com.traverse_strategy, backtrack_vec, sorting)
    com.c_step_nr += 1
    backtrack_obj.step_nr = com.c_step_nr

    # if we found the optimal solution or one feasible
    # => check whether all solutions are requested
    if !found && com.options.all_solutions
        l = length(backtrack_vec)
        backtrack_obj = backtrack_vec[l]
        while backtrack_obj.status == :Closed
            l -= 1
            if l == 0
                break
            end
            backtrack_obj = backtrack_vec[l]
        end
        if l != 0
            found = true
        end
    end

    if !found && com.options.all_optimal_solutions
        l = length(backtrack_vec)
        backtrack_obj = backtrack_vec[l]
        # get an obj which has the same bound as the optimal solution
        while backtrack_obj.best_bound != com.best_sol || backtrack_obj.status == :Closed
            l -= 1
            if l == 0
                break
            end
            backtrack_obj = backtrack_vec[l]
        end
        if l != 0
            found = true
        end
    end
    return found, backtrack_obj
end

"""
    get_next_node(com::CS.CoM, ::Val{:BFS}, backtrack_vec::Vector{BacktrackObj{T}}, sorting) where T <: Real

Get the next node we want to prune on if there is any. This uses best first search and if two
nodes have the same `best_bound` the deeper one is chosen.
Check other `get_next_node` functions for other possible traverse strategies.
Return whether a node was found and the corresponding backtrack_obj
"""
function get_next_node(
    com::CS.CoM,
    ::Val{:BFS},
    backtrack_vec::Vector{BacktrackObj{T}},
    sorting,
) where {T<:Real}
    found = false
    obj_factor = com.sense == MOI.MIN_SENSE ? 1 : -1
    backtrack_obj = backtrack_vec[1]

    # if there is no objective or sorting is set to false
    if com.sense == MOI.FEASIBILITY_SENSE || !sorting
        return get_first_open(backtrack_vec)
    end
    # sort for objective
    # don't actually sort => just get the best backtrack idx
    # the one with the best bound and if same best bound choose the one with higher depth
    l = 0
    best_fac_bound = typemax(Int)
    best_depth = 0
    for i in 1:length(backtrack_vec)
        bo = backtrack_vec[i]
        if bo.status == :Open
            if obj_factor * bo.best_bound < best_fac_bound ||
               (obj_factor * bo.best_bound == best_fac_bound && bo.depth > best_depth)
                l = i
                best_depth = bo.depth
                best_fac_bound = obj_factor * bo.best_bound
            end
        end
    end

    if l != 0
        backtrack_obj = backtrack_vec[l]
        found = true
    end
    return found, backtrack_obj
end

"""
    get_next_node(com::CS.CoM, :Val{:DFS}, backtrack_vec::Vector{BacktrackObj{T}}, sorting) where T <: Real

Get the next node we want to prune on if there is any. This uses depth first search and if two
nodes have the same depth the one with the better `best_bound` is chosen.
Check other `get_next_node` functions for other possible traverse strategies.
Return whether a node was found and the corresponding backtrack_obj
"""
function get_next_node(
    com::CS.CoM,
    ::Val{:DFS},
    backtrack_vec::Vector{BacktrackObj{T}},
    sorting,
) where {T<:Real}
    obj_factor = com.sense == MOI.MIN_SENSE ? 1 : -1
    backtrack_obj = backtrack_vec[1]

    # if sorting is set to false
    if !sorting
        return get_first_open(backtrack_vec)
    end # sort for depth
    # don't actually sort => just get the best backtrack idx
    # the one with the highest depth and if same choose better bound
    isempty(com.backtrack_pq) && return false, backtrack_obj

    backtrack_idx, _ = peek(com.backtrack_pq)

    backtrack_obj = backtrack_vec[backtrack_idx]

    return true, backtrack_obj
end


text:
# In England the currency is made up of pound, £, and pence, p, and there are
# eight coins in general circulation:
#
# 1p, 2p, 5p, 10p, 20p, 50p, £1 (100p) and £2 (200p).
# It is possible to make £2 in the following way:
#
# 1×£1 + 1×50p + 2×20p + 1×5p + 1×2p + 3×1p
# How many different ways can £2 be made using any number of coins?

using ProjectEulerSolutions

# Recursive solution count coin combinations.  In general, the idea is to
# consider each coin type and see how many ways there are to make a sum with
# that coin and lower denominations.  Different solutions can solve the
# recursion in slightly different ways.  Memoization is used to store lower
# level computatons.

# Uses a double recursion.
function p031solution_double(units::Array{Integer}=ones(Integer, 1), value::Integer=5)::Integer
    cache = zeros(Integer, value, length(units))
    return countcoins2(cache, units, length(units), value)
end

# Uses a single recursion, makes less total recursive calls.
function p031solution_single(units::Array{Integer}=ones(Integer, 1), value::Integer=5)::Integer
    cache = zeros(Integer, value, length(units))
    return countcoins1(cache, units, length(units), value)
end

# Uses no recursion, solution from https://projecteuler.net/thread=31.
function p031solution_dp(units::Array{Integer}=ones(Integer, 1), value::Integer=5)::Integer
    ways = zeros(Integer, value)
    ways[1] = 1
    for coin in units
        for j in (coin+1):value
            ways[j] += ways[j - coin]
        end
    end
    return ways[value]
end

# Most straightforward recursive solution.
function countcoins2(cache::Array{Integer,2}, units::Array{Integer}, coin::Integer, amount::Integer)::Integer
    if amount == 0 # Coin divides equally, add 1
        return 1
    end
    if coin < 1 || amount < 0 # Coins don't divide equally
        return 0
    end
    if cache[amount, coin] > 0
        return cache[amount, coin]
    end

    # Sum of solutions without coin and with coin
    val = countcoins2(cache, units, coin - 1, amount) +
          countcoins2(cache, units, coin, amount - units[coin])
    cache[amount, coin] = val
    return val
end

# Simple refactoring, less total recursive calls.
function countcoins1(cache::Array{Integer,2}, units::Array{Integer}, coin::Integer, amount::Integer)::Integer
    if amount == 0
        return 1
    end
    if cache[amount, coin] > 0
        return cache[amount, coin]
    end
    res = 0
    for i in coin:-1:2
        if amount - units[i] >= 0
            res += countcoins1(cache, units, i, amount - units[i])
        end
    end

    # Account for smallest coin evenly dividing into amount
    if amount % units[1] == 0
        res += 1
    end

    cache[amount, coin] = res
    return res
end

p031 = Problems.Problem(Dict("Double recusion" => p031solution_double,
                             "Single recusion" => p031solution_single,
                             "Dynamic programming" => p031solution_dp))

Problems.benchmark(p031, convert(Array{Integer}, [1, 2, 5, 10, 20, 50, 100, 200]), 200)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
using LinearAlgebra

"""
    lattice(s, N, ispbc=true)

Create a adjacency matrix of lattice `s` with size ``N``.
The periodic boundary condition is imposed by default.
Admissible values for `s` are:

| `s`                      | graph type                       |
| :------------------------|:---------------------------------|
| :chain                   | A one dimensional chain.         |
| :square                  | A two dimensional square lattice.|

# Examples
```julia-repl
julia> lattice(:chain, 3)
[0 1 1
 1 0 1
 1 1 0]

julia> lattice(:square, 4, ispbc=false)
[0 1 1 0
 1 0 0 1
 1 0 0 1
 0 1 1 0]
```
"""
function lattice(s::Symbol, N::Int; ispbc=true)
    graphmap = Dict(
        :chain => chain_graph,
        :square => square_graph,
    )

    if (s in keys(graphmap))
        return graphmap[s](N, ispbc)
    end
    throw(ArgumentError("$s is not a valid graph"))
end

"""
    chain_graph(N, ispbc)

An adjacency matrix of a chain graph with length ``N``.

# Examples
```julia-repl
julia> chain_graph(3)
[0 1 1
 1 0 1
 1 1 0]
```
"""
function chain_graph(N, ispbc)
    @assert N > 0 "N must be larger than 1."
    if N == 1
        return [0]
    end
    A = Matrix(SymTridiagonal(zeros(Int, N), ones(Int, N-1)))
    if ispbc
        A[1, N] += 1
        A[N, 1] += 1
    end
    A
end

"""
    issquare(N)

Check whether the number is perfect square or not.

# Examples
```julia-repl
julia> issquare(2)
false

julia> issquare(4)
true
```
"""
function issquare(N)
    isqrt(N)^2 == N
end

"""
    square_graph(N, ispbc)

An adjacency matrix of a square lattice whose number of vertices is ``N``.

# Examples
```julia-repl
julia> square_graph(4)
[0 2 2 0
 2 0 0 2
 2 0 0 2
 0 2 2 0]

julia> square_graph(4, ispbc=false)
[0 1 1 0
 1 0 0 1
 1 0 0 1
 0 1 1 0]
```
"""
function square_graph(N, ispbc)
    @assert N > 0 "N must be larger than 1."
    @assert issquare(N) "N must be a perfect square number."
    L = isqrt(N)
    if N == 1
        return [0]
    end

    A = transpose(reshape(1:N, L,L))
    h = [(A[CartesianIndex(i,j)], A[CartesianIndex(i,pbc(j+1,L))]) for i in 1:L for j in 1:L]
    v = [(A[CartesianIndex(i,j)], A[CartesianIndex(pbc(i+1,L),j)]) for i in 1:L for j in 1:L]
    bonds = vcat(v, h)

    h = [(A[CartesianIndex(i,1)], A[CartesianIndex(i,L)]) for i in 1:L]
    v = [(A[CartesianIndex(1,j)], A[CartesianIndex(L,j)]) for j in 1:L]
    pbcbonds = vcat(v, h)

    Res = zeros(Int, N, N)
    for b in bonds
        Res[b[1], b[2]] += 1
        Res[b[2], b[1]] += 1
    end
    if !ispbc
        for b in pbcbonds
            Res[b[1], b[2]] -= 1
            Res[b[2], b[1]] -= 1
        end
    end
    Res
end


text:
function outer_dGx_pos(loc_globaldata, loc_ghost_holder, dist_length, idx, configData, phi_i, phi_k)

    power::Float64 = configData["core"]["power"]
    limiter_flag::Float64 = configData["core"]["limiter_flag"]

    sum_delx_sqr = zero(Float64)
    sum_dely_sqr = zero(Float64)
    sum_delx_dely = zero(Float64)

    sum_delx_delf = zeros(Float64,4)
    sum_dely_delf = zeros(Float64,4)

    x_i = loc_globaldata[idx].x
    y_i = loc_globaldata[idx].y

    nx = loc_globaldata[idx].nx
    ny = loc_globaldata[idx].ny

    tx = ny
    ty = -nx

    G_i = zeros(Float64,4)
    G_k = zeros(Float64,4)
    result = zeros(Float64,4)

    for itm in loc_globaldata[idx].xpos_conn

        if itm <= dist_length
            globaldata_itm = loc_globaldata[itm]
        else
            globaldata_itm = loc_ghost_holder[1][itm]
        end
        x_k = globaldata_itm.x
        y_k = globaldata_itm.y

        delx = x_k - x_i
        dely = y_k - y_i

        dels = delx*tx + dely*ty
        deln = delx*nx + dely*ny

        dist = hypot(dels, deln)
        weights = dist ^ power

        dels_weights = dels*weights
        deln_weights = deln*weights

        sum_delx_sqr = sum_delx_sqr + dels*dels_weights
        sum_dely_sqr = sum_dely_sqr + deln*deln_weights

        sum_delx_dely = sum_delx_dely + dels*deln_weights

        qtilde_i = @. loc_globaldata[idx].q - 0.5*(delx*loc_globaldata[idx].dq[1] + dely*loc_globaldata[idx].dq[2])
        qtilde_k = @. globaldata_itm.q - 0.5*(delx*globaldata_itm.dq[1] + dely*globaldata_itm.dq[2])

        if limiter_flag == 1
            venkat_limiter(qtilde_i, loc_globaldata, idx, configData, phi_i)
            venkat_limiter(qtilde_k, loc_globaldata, idx, configData, phi_k)
            @. qtilde_i = loc_globaldata[idx].q - 0.5 * phi_i * (delx*loc_globaldata[idx].dq[1] + dely*loc_globaldata[idx].dq[2])
            @. qtilde_k = globaldata_itm.q - 0.5 * phi_k * (delx*globaldata_itm.dq[1] + dely*globaldata_itm.dq[2])
        end
                #if limiter_flag == 2
        #    maxi = max_q_values(globaldata, idx)
        #    mini = min_q_values(globaldata, idx)
#
        #    for i in 1:4
        #        if qtilde_i[i] > maxi[i]
        #            qtilde_i[i] = maxi[i]
        #        end
        #        if qtilde_i[i] < mini[i]
        #            qtilde_i[i] = mini[i]
        #        end
        #        if qtilde_k[i] > maxi[i]
        #            qtilde_k[i] = maxi[i]
        #        end
        #        if qtilde_k[i] < mini[i]
        #            qtilde_k[i] = mini[i]
        #        end
        #    end
        #end
        qtilde_to_primitive(result, qtilde_i, configData)
        flux_quad_GxIII(G_i, nx, ny, result[1], result[2], result[3], result[4])

        qtilde_to_primitive(result, qtilde_k, configData)
        flux_quad_GxIII(G_k, nx, ny, result[1], result[2], result[3], result[4])

        for i in 1:4
            sum_delx_delf[i] += (G_k[i] - G_i[i]) * dels_weights
            sum_dely_delf[i] += (G_k[i] - G_i[i]) * deln_weights
        end
    end

    det = @. sum_delx_sqr*sum_dely_sqr - sum_delx_dely*sum_delx_dely
    one_by_det = 1 / det

    return @. (sum_delx_delf*sum_dely_sqr - sum_dely_delf*sum_delx_dely)*one_by_det
end

function outer_dGx_neg(loc_globaldata, loc_ghost_holder, dist_length, idx, configData, phi_i, phi_k)

    power::Float64 = configData["core"]["power"]::Float64
    limiter_flag::Float64 = configData["core"]["limiter_flag"]::Float64

    sum_delx_sqr = zero(Float64)
    sum_dely_sqr = zero(Float64)
    sum_delx_dely = zero(Float64)

    sum_delx_delf = zeros(Float64, 4)
    sum_dely_delf = zeros(Float64, 4)

    x_i = loc_globaldata[idx].x
    y_i = loc_globaldata[idx].y

    nx = loc_globaldata[idx].nx
    ny = loc_globaldata[idx].ny

    tx = ny
    ty = -nx

    G_i = zeros(Float64,4)
    G_k = zeros(Float64,4)
    result = zeros(Float64,4)

    for itm in loc_globaldata[idx].xneg_conn

        if itm <= dist_length
            globaldata_itm = loc_globaldata[itm]
        else
            globaldata_itm = loc_ghost_holder[1][itm]
        end
        x_k = globaldata_itm.x
        y_k = globaldata_itm.y

        delx = x_k - x_i
        dely = y_k - y_i

        dels = delx*tx + dely*ty
        deln = delx*nx + dely*ny

        dist = hypot(dels, deln)
        weights = dist ^ power

        dels_weights = dels*weights
        deln_weights = deln*weights

        sum_delx_sqr = sum_delx_sqr + dels*dels_weights
        sum_dely_sqr = sum_dely_sqr + deln*deln_weights

        sum_delx_dely = sum_delx_dely + dels*deln_weights

        qtilde_i = @. loc_globaldata[idx].q - 0.5*(delx*loc_globaldata[idx].dq[1] + dely*loc_globaldata[idx].dq[2])
        qtilde_k = @. globaldata_itm.q - 0.5*(delx*globaldata_itm.dq[1] + dely*globaldata_itm.dq[2])

        if limiter_flag == 1
            venkat_limiter(qtilde_i, loc_globaldata, idx, configData, phi_i)
            venkat_limiter(qtilde_k, loc_globaldata, idx, configData, phi_k)
            @. qtilde_i = loc_globaldata[idx].q - 0.5 * phi_i * (delx*loc_globaldata[idx].dq[1] + dely*loc_globaldata[idx].dq[2])
            @. qtilde_k = globaldata_itm.q - 0.5 * phi_k * (delx*globaldata_itm.dq[1] + dely*globaldata_itm.dq[2])
        end

                #if limiter_flag == 2
        #    maxi = max_q_values(globaldata, idx)
        #    mini = min_q_values(globaldata, idx)
#
        #    for i in 1:4
        #        if qtilde_i[i] > maxi[i]
        #            qtilde_i[i] = maxi[i]
        #        end
        #        if qtilde_i[i] < mini[i]
        #            qtilde_i[i] = mini[i]
        #        end
        #        if qtilde_k[i] > maxi[i]
        #            qtilde_k[i] = maxi[i]
        #        end
        #        if qtilde_k[i] < mini[i]
        #            qtilde_k[i] = mini[i]
        #        end
        #    end
        #end
        qtilde_to_primitive(result, qtilde_i, configData)
        flux_quad_GxIV(G_i, nx, ny, result[1], result[2], result[3], result[4])

        qtilde_to_primitive(result, qtilde_k, configData)
        flux_quad_GxIV(G_k, nx, ny, result[1], result[2], result[3], result[4])

        for i in 1:4
            sum_delx_delf[i] += (G_k[i] - G_i[i]) * dels_weights
            sum_dely_delf[i] += (G_k[i] - G_i[i]) * deln_weights
        end
    end
    det = @. sum_delx_sqr*sum_dely_sqr - sum_delx_dely*sum_delx_dely
    one_by_det = 1 / det
    return @. (sum_delx_delf*sum_dely_sqr - sum_dely_delf*sum_delx_dely)*one_by_det
end

function outer_dGy_pos(loc_globaldata, loc_ghost_holder, dist_length, idx, configData, phi_i, phi_k)

    power::Float64 = configData["core"]["power"]::Float64
    limiter_flag::Float64 = configData["core"]["limiter_flag"]::Float64

    sum_delx_sqr = zero(Float64)
    sum_dely_sqr = zero(Float64)
    sum_delx_dely = zero(Float64)

    sum_delx_delf = zeros(Float64, 4)
    sum_dely_delf = zeros(Float64, 4)

    x_i = loc_globaldata[idx].x
    y_i = loc_globaldata[idx].y

    nx = loc_globaldata[idx].nx
    ny = loc_globaldata[idx].ny

    tx = ny
    ty = -nx

    G_i = zeros(Float64,4)
    G_k = zeros(Float64,4)
    result = zeros(Float64,4)

    for itm in loc_globaldata[idx].ypos_conn
        if itm <= dist_length
            globaldata_itm = loc_globaldata[itm]
        else
            globaldata_itm = loc_ghost_holder[1][itm]
        end
        x_k = globaldata_itm.x
        y_k = globaldata_itm.y

        delx = x_k - x_i
        dely = y_k - y_i

        dels = delx*tx + dely*ty
        deln = delx*nx + dely*ny

        dist = hypot(dels, deln)
        weights = dist ^ power

        dels_weights = dels*weights
        deln_weights = deln*weights

        sum_delx_sqr = sum_delx_sqr + dels*dels_weights
        sum_dely_sqr = sum_dely_sqr + deln*deln_weights

        sum_delx_dely = sum_delx_dely + dels*deln_weights

        qtilde_i = @. loc_globaldata[idx].q - 0.5*(delx*loc_globaldata[idx].dq[1] + dely*loc_globaldata[idx].dq[2])
        qtilde_k = @. globaldata_itm.q - 0.5*(delx*globaldata_itm.dq[1] + dely*globaldata_itm.dq[2])

        if limiter_flag == 1
            venkat_limiter(qtilde_i, loc_globaldata, idx, configData, phi_i)
            venkat_limiter(qtilde_k, loc_globaldata, idx, configData, phi_k)
            @. qtilde_i = loc_globaldata[idx].q - 0.5 * phi_i * (delx*loc_globaldata[idx].dq[1] + dely*loc_globaldata[idx].dq[2])
            @. qtilde_k = globaldata_itm.q - 0.5 * phi_k * (delx*globaldata_itm.dq[1] + dely*globaldata_itm.dq[2])
        end
                #if limiter_flag == 2
        #    maxi = max_q_values(globaldata, idx)
        #    mini = min_q_values(globaldata, idx)
#
        #    for i in 1:4
        #        if qtilde_i[i] > maxi[i]
        #            qtilde_i[i] = maxi[i]
        #        end
        #        if qtilde_i[i] < mini[i]
        #            qtilde_i[i] = mini[i]
        #        end
        #        if qtilde_k[i] > maxi[i]
        #            qtilde_k[i] = maxi[i]
        #        end
        #        if qtilde_k[i] < mini[i]
        #            qtilde_k[i] = mini[i]
        #        end
        #    end
        #end
        qtilde_to_primitive(result, qtilde_i, configData)
        flux_Gyp(G_i, nx, ny, result[1], result[2], result[3], result[4])

        qtilde_to_primitive(result, qtilde_k, configData)
        flux_Gyp(G_k, nx, ny, result[1], result[2], result[3], result[4])

        # print(sum_delx_delf)
        # print(" ")
        # print(G_k)
        # print(" ")
        # print(G_i * dels_weights)
        # print("\n")
        for i in 1:4
            sum_delx_delf[i] += (G_k[i] - G_i[i]) * dels_weights
            sum_dely_delf[i] += (G_k[i] - G_i[i]) * deln_weights
        end

    end
    det = @. sum_delx_sqr*sum_dely_sqr - sum_delx_dely*sum_delx_dely

    one_by_det = 1 / det

    return @. (sum_dely_delf*sum_delx_sqr - sum_delx_delf*sum_delx_dely)*one_by_det
end

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
export CDS_MVp_MT

"""
compute multi-threaded matrix vector product with vector x, output is vector y: y=A*x
MVP is in the compressed diagonal format.
R is a tall matrix N by ndiagonals, corresponding to a square matrix A
offsets indicate offset of diagonal compared to the main diagonal in A (which is 0)
"""
function CDS_MVp_MT(
                        N      ::Integer,
                        ndiags ::Integer,
                        R      ::Array{TF,2},
                        offset ::Vector{TI},
                        x      ::Vector{TF},
                        y      ::Vector{TF}) where {TF<:Real,TI<:Integer}

  for i = 1 : ndiags
      d  = offset[i]
      r0 = max(1, 1-d)
      r1 = min(N, N-d)
      c0 = max(1, 1+d)
      CDS_MVp_MT_subfunc(R,x,y,r0,c0,r1,i)
  end
  return y
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


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

mutable struct VirtualMachineStatusCodeCount <: SwaggerModel
    code::Any # spec type: Union{ Nothing, String } # spec name: code
    count::Any # spec type: Union{ Nothing, Int32 } # spec name: count

    function VirtualMachineStatusCodeCount(;code=nothing, count=nothing)
        o = new()
        validate_property(VirtualMachineStatusCodeCount, Symbol("code"), code)
        setfield!(o, Symbol("code"), code)
        validate_property(VirtualMachineStatusCodeCount, Symbol("count"), count)
        setfield!(o, Symbol("count"), count)
        o
    end
end # type VirtualMachineStatusCodeCount

const _property_map_VirtualMachineStatusCodeCount = Dict{Symbol,Symbol}(Symbol("code")=>Symbol("code"), Symbol("count")=>Symbol("count"))
const _property_types_VirtualMachineStatusCodeCount = Dict{Symbol,String}(Symbol("code")=>"String", Symbol("count")=>"Int32")
Base.propertynames(::Type{ VirtualMachineStatusCodeCount }) = collect(keys(_property_map_VirtualMachineStatusCodeCount))
Swagger.property_type(::Type{ VirtualMachineStatusCodeCount }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineStatusCodeCount[name]))}
Swagger.field_name(::Type{ VirtualMachineStatusCodeCount }, property_name::Symbol) =  _property_map_VirtualMachineStatusCodeCount[property_name]

function check_required(o::VirtualMachineStatusCodeCount)
    true
end

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


text:
module NDTensors

using Base.Threads
using Compat
using Dictionaries
using Random
using LinearAlgebra
using StaticArrays
using HDF5
using Requires
using Strided
using TimerOutputs
using TupleTools

using Base: @propagate_inbounds, ReshapedArray

using Base.Cartesian: @nexprs

using Base.Threads: @spawn

#####################################
# Imports and exports
#
include("exports.jl")
include("imports.jl")

#####################################
# DenseTensor and DiagTensor
#
include("aliasstyle.jl")
include("similar.jl")
include("tupletools.jl")
include("dims.jl")
include("tensorstorage.jl")
include("tensor.jl")
include("contraction_logic.jl")
include("dense.jl")
include("symmetric.jl")
include("linearalgebra.jl")
include("diag.jl")
include("combiner.jl")
include("truncate.jl")
include("svd.jl")

#####################################
# BlockSparseTensor
#
include("blocksparse/blockdims.jl")
include("blocksparse/block.jl")
include("blocksparse/blockoffsets.jl")
include("blocksparse/blocksparse.jl")
include("blocksparse/blocksparsetensor.jl")
include("blocksparse/diagblocksparse.jl")
include("blocksparse/combiner.jl")
include("blocksparse/linearalgebra.jl")

#####################################
# Empty
#
include("empty.jl")

#####################################
# Deprecations
#
include("deprecated.jl")

#####################################
# A global timer used with TimerOutputs.jl
#

const timer = TimerOutput()

#####################################
# Optional block sparse multithreading
#

include("blas_get_num_threads.jl")

const _using_threaded_blocksparse = Ref(false)

function enable_threaded_blocksparse_docstring(module_name)
  return """
             $(module_name).enable_threaded_blocksparse()
             $(module_name).disable_threaded_blocksparse()

         Enable or disable block sparse multithreading.

         Returns the current state of `$(module_name).using_threaded_blocksparse()`, i.e. `true` if threaded block sparse was previously enabled, and `false` if threaded block sparse was previously disabled. This is helpful for turning block sparse threading on or off temporarily. For example:
         ```julia
         using_threaded_blocksparse = $(module_name).enable_threaded_blocksparse()
         # Run code that you want to be threaded
         if !using_threaded_blocksparse
           $(module_name).disable_threaded_blocksparse()
         end
         ```

         Note that you need to start Julia with multiple threads. For example, to start Julia with 4 threads, you can use any of the following:
         ```
         \$ julia --threads=4

         \$ julia -t 4

         \$ JULIA_NUM_THREADS=4 julia
         ```

         In addition, we have found that it is best to disable `BLAS` and `Strided` multithreading when using block sparse multithreading. You can do that with the commands `using LinearAlgebra; BLAS.set_num_threads(1)` and `$(module_name).Strided.disable_threads()`.

         See also: `$(module_name).enable_threaded_blocksparse`, `$(module_name).disable_threaded_blocksparse`, `$(module_name).using_threaded_blocksparse`.
         """
end

function _enable_threaded_blocksparse()
  current_using_threaded_blocksparse = using_threaded_blocksparse()
  if !current_using_threaded_blocksparse
    if Threads.nthreads() == 1
      println(
        "WARNING: You are trying to enable block sparse multithreading, but you have started Julia with only a single thread. You can start Julia with `N` threads with `julia -t N`, and check the number of threads Julia can use with `Threads.nthreads()`. Your system has $(Sys.CPU_THREADS) threads available to use, which you can determine by running `Sys.CPU_THREADS`.\n",
      )
    end
    if blas_get_num_threads() > 1 && Threads.nthreads() > 1
      println(
        "WARNING: You are enabling block sparse multithreading, but BLAS $(BLAS.vendor()) is currently set to use $(blas_get_num_threads()) threads. When using block sparse multithreading, we recommend setting BLAS to use only a single thread, otherwise you may see suboptimal performance. You can set it with `using LinearAlgebra; BLAS.set_num_threads(1)`.\n",
      )
    end
    if Strided.get_num_threads() > 1
      println(
        "WARNING: You are enabling block sparse multithreading, but Strided.jl is currently set to use $(Strided.get_num_threads()) threads for performing dense tensor permutations. When using block sparse multithreading, we recommend setting Strided.jl to use only a single thread, otherwise you may see suboptimal performance. You can set it with `NDTensors.Strided.disable_threads()` and see the current number of threads it is using with `NDTensors.Strided.get_num_threads()`.\n",
      )
    end
    _using_threaded_blocksparse[] = true
  end
  return current_using_threaded_blocksparse
end

function _disable_threaded_blocksparse()
  current_using_threaded_blocksparse = using_threaded_blocksparse()
  if current_using_threaded_blocksparse
    _using_threaded_blocksparse[] = false
  end
  return current_using_threaded_blocksparse
end

"""
$(enable_threaded_blocksparse_docstring(@__MODULE__))
"""
using_threaded_blocksparse() = _using_threaded_blocksparse[]

"""
$(enable_threaded_blocksparse_docstring(@__MODULE__))
"""
enable_threaded_blocksparse() = _enable_threaded_blocksparse()

"""
$(enable_threaded_blocksparse_docstring(@__MODULE__))
"""
disable_threaded_blocksparse() = _disable_threaded_blocksparse()

#####################################
# Optional TBLIS contraction backend
#

const _using_tblis = Ref(false)

using_tblis() = _using_tblis[]

function enable_tblis()
  _using_tblis[] = true
  return nothing
end

function disable_tblis()
  _using_tblis[] = false
  return nothing
end

function __init__()
  @require TBLIS = "48530278-0828-4a49-9772-0f3830dfa1e9" begin
    enable_tblis()
    include("tblis.jl")
  end
  @require Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" begin
    include("octavian.jl")
  end
end

end # module NDTensors


text:
using TimedAutomata, GraphMakie
##

states = ["S0", "S1", "S2", "S3"]
initial_state = "S0"
clocks = [:x,:y,:z]
symbols = Set(Symbol[])
arcs = Vector([
    @arc("S0","S1",true,TAMessage(),[:z],Int64),
    @arc("S1","S2",y>2,TAMessage(),[:y],Int64),
    @arc("S2","S3",x-z<1 && z-y<1,TAMessage(),[],Int64)
])
invariants = Dict("S0"=>@constraint x-y==0 && y-z==0 && z-x==0 Int64)

automaton1 = TA(
    states,
    initial_state,
    clocks,
    symbols,
    arcs,
    invariants
)

##
states = ["start", "loop", "end"]
initial_state = "start"
clocks = [:x,:y]
symbols = Set(Symbol[])
arcs = Vector([
    @arc("start","loop",true,TAMessage(),[:x,:y],Int64),
    @arc("loop","loop",x==10,TAMessage(),[:x],Int64),
    @arc("loop","end",y≥20,TAMessage(),[:x,:y],Int64)
])
invariants = Dict("start"=>@constraint(x-y==0,Int64), "loop"=>@constraint(x≤10,Int64))

automaton2 = TA(
    states,
    initial_state,
    clocks,
    symbols,
    arcs,
    invariants
)

##
automaton4=TA(
    ["off","dim","bright"],
    "off",
    [:x],
    Set([:press]),
    Vector([
        @arc("off", "dim", true, TAMessage(MSG_IN, :press), [:x]),
        @arc("dim", "off", x>10//1, TAMessage(MSG_IN, :press), []),
        @arc("dim", "bright", x≤10//1, TAMessage(MSG_IN, :press), []),
        @arc("bright", "off", true, TAMessage(MSG_IN, :press), [])
    ]),
    Dict{String,TAConstraint{Rational{Int}}}()
)
zg4 = TTS(automaton4)

automaton5=TA(
    ["idle","t","study", "relax"],
    "idle",
    [:y],
    Set([:press]),
    Vector([
        @arc("t", "study", true, TAMessage(MSG_OUT, :press), []),
        @arc("study", "study", true, TAMessage(), []),
        @arc("study", "idle", true, TAMessage(MSG_OUT, :press), []),
        @arc("idle", "t", true, TAMessage(MSG_OUT, :press), [:y]),
        @arc("idle", "relax", true, TAMessage(MSG_OUT, :press), [:y]),
        @arc("relax", "idle", y>10//1, TAMessage(MSG_OUT, :press), [])
    ]),
    Dict(
        "t" => @constraint y<5//1
    )
)

automaton45 = automaton4 | automaton5
automaton44 = (automaton4 | automaton4)
automaton54 = automaton5 | automaton4

find_isomorphism(automaton45,automaton54)

automaton445 = |(automaton4, automaton4, automaton5)

##

automata = ("Example 1"=>automaton1, "Example 2"=>automaton2, "Lamp"=>automaton4, "Student"=>automaton5, "Lamp + Student"=>automaton45, "2 Lamps + Student"=>automaton445);
for (name,automaton) in automata
    for condition in ("automaton", "pruned automaton", "zone graph")
        obj=if condition=="automaton"
            automaton
        elseif condition=="pruned automaton"
            prune(automaton)
        elseif condition=="zone graph"
            TTS(automaton)
        end

        f,ax,_ = graphplot(obj)
        hidedecorations!(ax)
        ax.title[] = "$(name) ($(condition))"
        save(joinpath("examples","figures","$(name)_$(condition).png"), f)
    end
end


text:
srand(1)
M = 10_000_000; K = 100
svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:M÷K], M)
# memcmp("abc","def",3)
memcmp(a::String, b::String, sz) = ccall(:memcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}, Int32), pointer(a), pointer(b), sz)

@time sort!(svec1, lt=(x,y)->memcmp(x,y,8)<0, alg=QuickSort); # 2.7 seconds
issorted(svec1)


text:
function CLT(
    dist;
    boundingbox = BoundingBox(O + (-250, -120), O + (250, 120)),
    margin = 5,
    pathname = "central_limit_theorem",
)

    function ground(args...)
        background("black")
        sethue("white")
    end

    function fixed_gaussian(loc_hist; color)
        (values, i, lowpos, highpos, barwidth, scaledvalue) -> begin

            minvalue, maxvalue = extrema(loc_hist.weights)
            barchartheight = boxheight(boundingbox) - 2margin
            minbarrange = minvalue - abs(minvalue)
            maxbarrange = maxvalue + abs(maxvalue)
            @layer begin

                sethue(color)
                if i <= length(loc_hist.edges[1])
                    scaledgaussvalue =
                        rescale(
                            pdf(
                                gauss,
                                loc_hist.edges[1][i + 1] - loc_hist.edges[1].step.hi / 2,
                            ),
                            minbarrange,
                            maxbarrange,
                        ) * barchartheight
                    circle(lowpos - (0, scaledgaussvalue), 3, :fill)
                end
            end

            tickline(
                boxbottomleft(boundingbox) - Point(0, margin),
                boxbottomright(boundingbox) - Point(0, margin),
                startnumber = loc_hist.edges[1][1],
                finishnumber = loc_hist.edges[1][end],
            )
        end
    end

    function hist_bar(; color)
        (values, i, lowpos, highpos, barwidth, scaledvalue) -> begin
            @layer begin
                sethue(color)
                Luxor.setline(barwidth)
                line(lowpos, highpos, :stroke)
            end
        end
    end

    # Parameters for the barchart function
    boundingbox = BoundingBox(O + (-250, -120), O + (250, 120))
    margin = 5

    # number of frames
    n_frames = 700

    # Number of histograms shown. To make the animation
    # slower set this to a lower value, less histograms
    # will be shown but they will reach the same n for
    # convergence as set by n_frames
    n_hists = 700

    # number of samples at each n
    n_samples = 10000

    # Change to adjust bar colors and gauss colors
    barcolor = HSV(colorant"goldenrod1")
    gausscolor = HSV(240, barcolor.s, barcolor.v)

    # Sample the mean r.v. for increasing values of n (1 to n_hists)
    samples = map(1:n_hists) do n
        map(1:n_samples) do _
            (mean(rand(dist, n)) - mean(dist)) * sqrt(n)
        end
    end

    finalmin, finalmax = extrema(samples[end])

    # Turn the samples into histograms
    steps = map(samples) do sampling
        hist = fit(
            Histogram,
            sampling,
            weights(ones(length(sampling))),
            range(finalmin, finalmax, length = 100),
        )
        hist = StatsBase.normalize(hist, mode = :pdf)
        hist
    end

    # Define the gaussian distribution where they should converge
    gauss = Normal(0, StatsBase.std(dist))

    my_video = Video(800, 800)
    Background(1:n_frames, ground)

    step_size = n_frames ÷ n_hists
    frame_brakes = 1:step_size:n_frames

    # Fix some point where writings will be shown
    titlepoint = Point(0, -100)
    distpoint = Point(-200, 0)
    counterpoint = Point(200, 0)

    # The last and thus hopefully most closely converged histogram in our sequence
    final_hist = steps[end]

    for (frame_n, hist) in zip(frame_brakes, steps)

        Object(
            frame_n:(frame_n + step_size - 1),
            @JShape begin
                barchart(
                    hist.weights,
                    boundingbox = boundingbox,
                    margin = margin,
                    labels = true,

                    # Provide the hist_bar we defined as the barfunction
                    barfunction = hist_bar(color = barcolor),

                    # Provide the fixed_gaussian we defined as the labelfunction
                    # it will plot the gaussian dots and the ticks on the bottom
                    labelfunction = fixed_gaussian(final_hist, color = gausscolor),
                )
            end
        )

        # The counter digits
        Object(
            frame_n:(frame_n + step_size - 1),
            @JShape begin
                sethue(barcolor)
                fontsize(15)
                text(string(frame_n ÷ step_size), counterpoint, halign = :center)
            end
        )
    end

    # All the writings except the changing digit.
    Object(
        1:n_frames,
        @JShape begin
            fontsize(40)
            text("Central Limit Theorem", titlepoint, halign = :center)

            fontsize(20)
            label("N", :N, counterpoint, offset = 20)
            label("Distribution", :N, distpoint, offset = 20)

            sethue(barcolor)
            fontsize(15)
            text(
                # Gather the distribution name and parameters
                # only works if dist is from Distributions.jl
                join([string(typeof(dist).name.name); string(params(dist))], ""),
                distpoint,
                halign = :center,
            )
        end
    )

    render(my_video, pathname = pathname)
end

function pendulum(;
    shift = O,
    basepoint = shift + O,
    radius = 300,
    margin = 10,
    k = 0.7π,
    first_frame = 1,
    last_frame = 500,
    speed_up = 3,
)

    ## Constants
    n_frames = last_frame - first_frame

    ## Constant shapes
    Object(first_frame:last_frame, (args...) -> arc(O, radius, 0, π, :stroke))
    Object(
        first_frame:last_frame,
        (args...) -> Javis.latex(
            L"\Huge{\theta(t) = \pi e^{-\frac{t}{\tau}} \cos(\omega t)}",
            Point(150, -300),
            :middle,
            :center,
        ),
    )
    Object(
        first_frame:last_frame,
        JLine(basepoint - (radius, 0), basepoint + (radius, 0), color = "white"),
    )
    Object(
        first_frame:last_frame,
        JLine(basepoint, basepoint + (0, radius), color = "white"),
    )
    Object(
        first_frame:last_frame,
        JLine(
            Point(-radius, -radius / 2 - margin),
            Point(radius, -radius / 2 - margin),
            linewidth = 1,
            color = "white",
        ),
    )

    Object(
        first_frame:last_frame,
        @JShape begin
            setline(2)
            tickline(
                Point(-radius, -radius - margin),
                Point(-radius, -margin),
                startnumber = k,
                finishnumber = -k,
            )
        end
    )

    motion(t) = k * ℯ^(-3 * (t) / 2) * cos(10t)

    points = map(first_frame:last_frame) do idx
        # Frame constants
        t = speed_up * (idx - first_frame + 1) / n_frames
        p = shift + radius * Point(sin(motion(t)), cos(motion(t)))
        endpoint = shift + p

        # Lines
        Object(idx:idx, @JShape begin
            setline(3)
            sethue("white")
            line(basepoint, shift + p, :stroke)
            setline(1)
            sethue("orange")
            line(endpoint, Point(endpoint.x, 0), :stroke)
            line(endpoint, Point(0, endpoint.y), :stroke)
        end)
        Object(idx:idx, JCircle(endpoint, 20, action = :fill, color = "red"))

        # Labels
        Object(
            idx:idx,
            (args...) -> begin
                fontsize(25)
                cos_x_val = midpoint(basepoint, endpoint).x
                Javis.latex(
                    L"\sin(\theta)",
                    Point(cos_x_val, endpoint.y) + (0, -10),
                    :middle,
                    :center,
                )
                sin_y_val = midpoint(basepoint, endpoint).y
                Javis.latex(
                    L"\cos(\theta)",
                    Point(endpoint.x, sin_y_val) + (-4, 0),
                    :middle,
                    :center,
                )
            end,
        )

        # Arcs
        Object(
            idx:idx,
            (args...) -> begin
                mypoint = Point(sin(motion(t)), cos(motion(t)))
                if motion(t) < 0
                    move(Point(0, 40))
                    arc(basepoint, 40, π / 2, π / 2 - motion(t), :stroke)
                else
                    move(40 * mypoint)
                    arc(basepoint, 40, π / 2 - motion(t), π / 2, :stroke)
                end
                fontsize(25)
                Javis.latex(L"\theta", basepoint + (-2margin, 2margin), :middle, :center)
            end,
        )

        Object(
            idx:idx,
            (args...) -> begin
                # Tangent component
                Luxor.arrow(
                    endpoint,
                    endpoint +
                    radius / 3 *
                    cos(motion(t) - π / 2) *
                    Point(sin(motion(t) - π / 2), cos(motion(t) - π / 2)),
                    linewidth = 2,
                )

                # Perpendicular component
                Luxor.arrow(
                    endpoint,
                    endpoint +
                    radius / 3 *
                    sin(motion(t) - π / 2) *
                    Point(-cos(motion(t) - π / 2), sin(motion(t) - π / 2)),
                    linewidth = 2,
                )

                # Full force
                sethue("grey")
                Luxor.arrow(endpoint, endpoint + radius / 3 * Point(0, 1))
            end,
        )

        # Plot
        Object(
            idx:last_frame,
            JCircle(
                Point(
                    -radius +
                    margin +
                    (idx - first_frame + 1) * (2 * radius - margin) / n_frames,
                    -radius / 2 - margin - radius * motion(t) / 2k,
                ),
                2,
                action = :fill,
                color = "red",
            ),
        )
    end

    Object(first_frame:last_frame, JCircle(basepoint, 10, action = :fill, color = "blue"))

end


text:
using Distributions

# In the literature Cauchy distributions have been used for sampling the
# f and cr constants used in DE.
function sample_bimodal_cauchy_once(cauchyDistrs, cutoffProb = 0.5)
  index = (rand() < cutoffProb) ? 1 : 2
  rand(cauchyDistrs[index])
end

# When sampling it is common to truncate in either or both ends.
function sample_bimodal_cauchy(cauchyDistrs; cutoffProb = 0.5,
  truncateAbove1 = true, truncateBelow0 = true)
  value = sample_bimodal_cauchy_once(cauchyDistrs, cutoffProb)
  if value > 1.0
    if truncateAbove1
      return 1.0
    else
      return sample_bimodal_cauchy(cauchyDistrs; cutoffProb = cutoffProb, truncateAbove1 = truncateAbove1, truncateBelow0 = truncateBelow0)
    end
  elseif value < 0.0
    if truncateBelow0
      return 0.0
    else
      return sample_bimodal_cauchy(cauchyDistrs; cutoffProb = cutoffProb, truncateAbove1 = truncateAbove1, truncateBelow0 = truncateBelow0)
    end
  else
    return value
  end
end

# For sampling f in DE, bimodal_cauchy(0.65, 0.1, 1.0, 0.1) has been proposed.
# For sampling cr in DE, bimodal_cauchy(0.1, 0.1, 0.95, 0.1) has been proposed.
function bimodal_cauchy(location1, scale1, location2, scale2)
  (Cauchy(location1, scale1), Cauchy(location2, scale2))
end


text:
#Exponential Loss

using LaTeXStrings, LossFunctions, Measures, Latexify

Plots.plot(ExpLoss(), c=:olive, leg=:topright, linewidth=2, annotation=(1, 1.5, L"L\left( a \right) = e^{-a}", 18))

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
export  parsexml, readxml,
        # Debug, remove later
        xmlparser, parse_node, parse_element

tagstrip(tag::AbstractString) = strip(tag, ['<', '>', '/'])

function match_closing_tag(n::ElementNode, t::Token)
    close_tagname = tagstrip(t.lexeme)
    if n.name != close_tagname
        error("opening tag '$(n.name)' does not match closing tag '$close_tagname'")
    end
end

function add_parsed_nodes!(parser::Parser, parent::ElementNode)
    t = peek_token(parser)
    while t.kind == BEGIN_TAG || t.kind == TEXT
        addchild!(parent, parse_node(parser))
        t = peek_token(parser)
    end
end

function add_parsed_attribute_nodes!(parser::Parser, parent::ElementNode)
    t = peek_token(parser)
    while t.kind == IDENT
        t = expect(parser, IDENT)
        name = t.lexeme
        expect(parser, EQUAL)
        t = expect(parser, STRING)
        value = t.lexeme
        parent[name] = value
        t = peek_token(parser)
    end
end

function parse_text(parser::Parser)
    t = expect(parser, TEXT)
    TextNode(t.lexeme)
end

function parse_element(parser::Parser)
    t = peek_token(parser)
    expect(parser, BEGIN_TAG)
    tagname = tagstrip(t.lexeme)
    n = ElementNode(tagname)

    t = peek_token(parser)
    while t.kind in [END_TAG, IDENT]
        if t.kind == END_TAG
            expect(parser, END_TAG)
            add_parsed_nodes!(parser, n)
        elseif t.kind == IDENT
            add_parsed_attribute_nodes!(parser, n)
        end
        t = peek_token(parser)
    end

    t = next_token(parser)
    if t.kind == CLOSE_TAG
        match_closing_tag(n, t)
    end

    if t.kind ∉ [CLOSE_TAG, END_AND_CLOSE_TAG]
        error("Element node needs to end with /> or </$(n.name)> not '$t'")
    end

    return n
end

function parse_node(parser::Parser)
    t = peek_token(parser)
    if t.kind == BEGIN_TAG
        parse_element(parser)::Node
    elseif t.kind == TEXT
        parse_text(parser)::Node
    else
        error("Had not expected token '$t' while looking for start of new XML node")
    end
end

xmlparser(s::AbstractString) = Parser(lex_xml(s))

function strip_xml_header(xmlstring::AbstractString, ignore_declaration::Bool)
    # Get the XML declaration. It is not part of the XML DOM, so we want
    # to exclude it.
    r = findfirst(r"<\?xml.*\?>", xmlstring)
    s = xmlstring # assume there is no XML declaration until proven otherwise
    if r == nothing
        ignore_declaration || @warn "Did not find any XML declaration such as <?xml version=\"1.0\" encoding=\"UTF-8\"?>"
    else
        s = xmlstring[last(r)+1:end]
        # Check encoding used
        decl = xmlstring[r]
        m = match(r"encoding=\"([\w-]+)\"", decl)
        if isempty(m.captures)
            @warn "Could not determine encoding, will assume UTF-8"
        else
            encoding = m.captures[1]
            if uppercase(encoding) != "UTF-8"
                @warn "XML parser is not made to handle $encoding encoding. XML file should have UTF-8 encoding"
            end
        end
    end

    # Skip DOC type as we don't handle it
    r = findfirst(r"<!DOCTYPE[^>]+>", s)
    if r != nothing
        s = s[last(r)+1:end]
    end
    return s
end

"""
    parsexml(xmlstring::AbstractString; ignore_declaration=false)

Parse a text string containing XML, and return an XML document object
"""
function parsexml(xmlstring::AbstractString; ignore_declaration=false)
    s = strip_xml_header(xmlstring, ignore_declaration)
    l = lex_xml(s)
    p = Parser(l)
    Document(parse_element(p))
end

"""
    readxml(stream::IO=stdin)
    readxml(filename::AbstractString)

Read and parse XML from an I/O stream or file. Result is returned as an XML document object.
"""
function readxml(stream::IO=stdin)
    text = read(stream, String)
    parsexml(text)
end

function readxml(filename::AbstractString)
    open(readxml, filename)
end


text:
module NamedVectors

using Printf: @printf
using ..Factories: Inference,TypeFactory,FunctionFactory,addargs!,extendbody!
using ..Factories: MixEscaped,Escaped,UnEscaped
using ..TypeTraits: efficientoperations

export NamedVector,@namedvector
export HomoNamedVector,@homonamedvector

"""
    NamedVector

Abstract type for all named vectors.
"""
abstract type NamedVector end

"""
    getindex(nv::NamedVector,index::Int)

Get the value by the `[]` syntax.
"""
Base.getindex(nv::NamedVector,index::Int)=getfield(nv,index)

"""
    setindex!(nv::NamedVector,value,index::Int)

Set the value by the `[]` syntax if mutable.
"""
Base.setindex!(nv::NamedVector,value,index::Int)=setfield!(nv,index,value)

"""
    ==(nv1::NamedVector,nv2::NamedVector) -> Bool
    isequal(nv1::NamedVector,nv2::NamedVector) -> Bool

Overloaded equivalent operator. Two named vector are equal to each other if and only if their keys as well as their values are equal to each other.
!!! note
    It is not necessary for two named vectors to be of the same concrete type to be equal to each other.
"""
Base.:(==)(nv1::NamedVector,nv2::NamedVector)=nv1|>keys==nv2|>keys && nv1|>values==nv2|>values
Base.isequal(nv1::NamedVector,nv2::NamedVector)=isequal(nv1|>keys,nv2|>keys) && isequal(nv1|>values,nv2|>values)

"""
    <(nv1::NamedVector,nv2::NamedVector) -> Bool
    isless(nv1::NamedVector,nv2::NamedVector) -> Bool

Compare two named vectors and judge whether the first is less than the second.
"""
Base.:<(nv1::NamedVector,nv2::NamedVector) = <(efficientoperations,nv1,nv2)
Base.isless(nv1::NamedVector,nv2::NamedVector)=isless(efficientoperations,nv1,nv2)

"""
    show(io::IO,nv::NamedVector)

Show a concrete `NamedVector`.
"""
Base.show(io::IO,nv::NamedVector)=@printf io "%s(%s)" nv|>typeof|>nameof join(repr.(nv|>values),',')

"""
    hash(nv::NamedVector,h::UInt)

Hash a concrete `NamedVector`.
"""
Base.hash(nv::NamedVector,h::UInt)=hash(nv|>values,h)

"""
    convert(::Type{Tuple},nv::NamedVector) -> Tuple
    convert(::Type{NV},nv::Tuple) where NV<:NamedVector -> NV

Convert a named vector to tuple and vice versa.
"""
@generated Base.convert(::Type{Tuple},nv::NamedVector)=Expr(:tuple,(:(getfield(nv,$i)) for i=1:(nv|>fieldcount))...)
function Base.convert(::Type{NV},nv::Tuple) where NV<:NamedVector
    @assert NV|>fieldcount==nv|>length "convert error: dismatched length between $NV($(NV|>fieldcount)) and input tuple($(nv|>length))."
    return NV(nv...)
end

"""
    length(::Type{NV}) where NV<:NamedVector -> Int
    length(nv::NamedVector) -> Int

Get the length of a concrete `NamedVector`.
"""
Base.length(::Type{NV}) where NV<:NamedVector=NV|>fieldcount
Base.length(nv::NamedVector)=nv|>typeof|>length

"""
    zero(::Type{NV}) where NV<:NamedVector -> NV
    zero(nv::NamedVector) -> typeof(nv)

Get a concrete `NamedVector` with all values being zero.
"""
@generated Base.zero(::Type{NV}) where NV<:NamedVector=(zeros=(zero(fieldtype(NV,i)) for i=1:(NV|>fieldcount));:(NV($(zeros...))))
Base.zero(nv::NamedVector)=nv|>typeof|>zero

"""
    iterate(nv::NamedVector,state=1)
    iterate(rv::Iterators.Reverse{<:NamedVector},state=length(rv.itr))

Iterate or reversely iterate over the values of a concrete `NamedVector`.
"""
Base.iterate(nv::NamedVector,state=1)=state>length(nv) ? nothing : (getfield(nv,state),state+1)
Base.iterate(rv::Iterators.Reverse{<:NamedVector},state=length(rv.itr))=state<1 ? nothing : (getfield(rv.itr,state),state-1)

"""
    keys(nv::NamedVector) -> NTuple(nv|>fieldcount,Symbol)
    values(nv::NamedVector) -> Tuple
    pairs(nv::NamedVector) -> Base.Generator

Iterate over the names.
"""
Base.keys(nv::NamedVector)=nv|>typeof|>fieldnames
Base.values(nv::NamedVector)=convert(Tuple,nv)
Base.pairs(nv::NamedVector)=Base.Generator(=>,keys(nv),values(nv))

"""
    replace(nv::NamedVector;kwargs...) -> typeof(nv)

Return a copy of a concrete `NamedVector` with some of the field values replaced by the keyword arguments.
"""
Base.replace(nv::NamedVector;kwargs...)=replace(efficientoperations,nv;kwargs...)

"""
    map(f,nvs::NV...) where NV<:NamedVector -> NV

Apply function `f` elementwise on the input named vectors.
"""
@generated function Base.map(f,nvs::NV...) where NV<:NamedVector
    exprs=Vector{Expr}(undef,NV|>fieldcount)
    for i=1:length(exprs)
        tmp=[:(nvs[$j][$i]) for j=1:length(nvs)]
        exprs[i]=:(f($(tmp...)))
    end
    return :(($NV)($(exprs...)))
end

"""
    @namedvector structdef::Expr

Decorate a "raw" struct to be a subtype of `NamedVector`. Here, "raw" means that the input struct has no explicit supertype and no inner constructors.
"""
macro namedvector(structdef::Expr)
    tf=TypeFactory(structdef)
    @assert tf.supertype==Inference(:Any) "@namedvector error: no explicit supertype except `Any` is allowed."
    @assert length(tf.constructors)==0 "@namedvector error: no inner constructor is allowed."
    tf.supertype=Inference(:NamedVector)
    paramnames=tuple((param.name for param in tf.params)...)
    fieldnames=tuple((field.name for field in tf.fields)...)
    fldnm=FunctionFactory(name=:(Base.fieldnames))
    addargs!(fldnm,:(::Type{<:$(tf.name)}))
    extendbody!(fldnm,Expr(:tuple,QuoteNode.(fieldnames)...))
    structdef=tf(MixEscaped(Escaped(tf.name),UnEscaped(paramnames...,:NamedVector)))
    fldnmdef=fldnm(MixEscaped(Escaped(:tuple)))
    return Expr(:block,:(Base.@__doc__($structdef)),fldnmdef)
end

"""
    HomoNamedVector{T}

Abstract type for all homogeneous named vectors.
"""
abstract type HomoNamedVector{T} <: NamedVector end

"""
    eltype(::Type{<:HomoNamedVector{T}}) where T
    eltype(nv::HomoNamedVector)

Get the type parameter of a concrete `HomoNamedVector`.
"""
Base.eltype(::Type{<:HomoNamedVector{T}}) where T=T
Base.eltype(nv::HomoNamedVector)=nv|>typeof|>eltype

"""
    @homonamedvector typename fieldnames dtype::Union{Expr,Symbol}=:nothing mutable::Union{Expr,Bool}=false

Construct a concrete homogeneous named vector with the type name being `typename` and the fieldnames specified by `fieldnames`, and optionally, the type parameters specified by `dtype`.`mutable` can be used as a keyword argument to determine whether the concrete type is mutable.
"""
macro homonamedvector(typename,fieldnames,dtype::Union{Expr,Symbol}=:nothing,mutable::Union{Expr,Bool}=false)
    typename=Symbol(typename)
    fieldnames=tuple(eval(fieldnames)...)
    @assert all(isa(name,Symbol) for name in fieldnames) "homonamedvector error: every field name should be a `Symbol`."
    isa(dtype,Expr) && dtype.head==:(<:) && (@assert dtype.args|>length==1 "homonamedvector error: wrong `dtype`.")
    dname,dscope=dtype==:nothing ? (:T,:Any) : isa(dtype,Expr) && dtype.head==:(<:) ? (:T,dtype.args[1]) : (dtype,:concrete)
    if isa(mutable,Expr)
        @assert mutable.head==:(=) && mutable.args[1]==:mutable && isa(mutable.args[2],Bool) "homonamedvector error: wrong `mutable`."
        mutable=mutable.args[2]
    end
    if dscope==:concrete
        new=:($(esc(typename)))
        super=:(HomoNamedVector{$(esc(dname))})
        body=(:($field::$(esc(dname))) for field in fieldnames)
    else
        new=:($(esc(typename)){$dname<:$(esc(dscope))})
        super=:(HomoNamedVector{$dname})
        body=(:($field::$dname) for field in fieldnames)
    end
    structdef=Expr(:struct,mutable,Expr(:<:,new,super),Expr(:block,body...))
    functions=:(Base.fieldnames(::Type{<:$(esc(typename))})=$fieldnames)
    return Expr(:block,:(Base.@__doc__($structdef)),functions)
end

end #module


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

const oldversion = VERSION < v"0.6.0-dev.2390"

# Test whether r[i] has full precision
function test_linspace{T<:AbstractFloat}(r::Ranges.LinSpace{T})
    isempty(r) && return nothing
    n = length(r)
    f, l = first(r), last(r)
    a, b, d = BigFloat(f), BigFloat(l), max(n-1,1)
    Δ = 2.1*max(eps(f), eps(l))  # FIXME: tighten these bounds
    for i = 1:n
        c = b*(BigFloat(i-1)/d) + a*(BigFloat(d-i+1)/d)
        @test abs(r[i]-T(c)) <= Δ
    end
    nothing
end

L32 = Ranges.linspace(Int32(1), Int32(4), 4)
L64 = Ranges.linspace(Int64(1), Int64(4), 4)
@test L32[1] == 1 && L64[1] == 1
@test L32[2] == 2 && L64[2] == 2
@test L32[3] == 3 && L64[3] == 3
@test L32[4] == 4 && L64[4] == 4

@test Ranges.range(Int64(1), Int64(1), 4) === L64

test_linspace(Ranges.linspace(0.1,0.3,3))
test_linspace(Ranges.linspace(0.0,0.3,4))
test_linspace(Ranges.linspace(0.3,-0.1,5))
test_linspace(Ranges.linspace(0.1,-0.3,5))
test_linspace(Ranges.linspace(0.0,1.0,11))
test_linspace(Ranges.linspace(0.0,1.0,0))
test_linspace(Ranges.linspace(0.0,-1.0,0))
test_linspace(Ranges.linspace(0.0,-1.0,11))
test_linspace(Ranges.linspace(1.0,27.0,1275))
test_linspace(Ranges.linspace(0.0,2.1,4))
test_linspace(Ranges.linspace(0.0,3.3,4))
test_linspace(Ranges.linspace(0.1,3.4,4))
test_linspace(Ranges.linspace(0.0,3.9,4))
test_linspace(Ranges.linspace(0.1,4.0,4))
test_linspace(Ranges.linspace(1.1,3.3,3))
test_linspace(Ranges.linspace(0.3,1.1,9))
test_linspace(Ranges.linspace(0.0,0.0,1))
test_linspace(Ranges.linspace(0.0,0.0,1))

for T = (Float32, Float64,),# BigFloat),
    a = -5:25, s = [-5:-1;1:25;], d = 1:25, n = -1:15
    den   = convert(T,d)
    start = convert(T,a)/den
    stop  = convert(T,(a+(n-1)*s))/den
    test_linspace(Ranges.linspace(start, stop, max(n,0)))
end

# linspace & ranges with very small endpoints
for T = (Float32, Float64)
    z = zero(T)
    u = eps(z)
    @test first(Ranges.linspace(u,u,0)) == u
    @test last(Ranges.linspace(u,u,0)) == u
    @test first(Ranges.linspace(-u,u,0)) == -u
    @test last(Ranges.linspace(-u,u,0)) == u
    @test [Ranges.linspace(-u,u,0);] == []
    @test [Ranges.linspace(-u,-u,1);] == [-u]
    @test [Ranges.linspace(-u,u,2);] == [-u,u]
    @test [Ranges.linspace(-u,u,3);] == [-u,0,u]
    @test first(Ranges.linspace(-u,-u,0)) == -u
    @test last(Ranges.linspace(-u,-u,0)) == -u
    @test first(Ranges.linspace(u,-u,0)) == u
    @test last(Ranges.linspace(u,-u,0)) == -u
    @test [Ranges.linspace(u,-u,0);] == []
    @test [Ranges.linspace(u,u,1);] == [u]
    @test [Ranges.linspace(u,-u,2);] == [u,-u]
    @test [Ranges.linspace(u,-u,3);] == [u,0,-u]
    v = [Ranges.linspace(-u,u,12);]
    @test length(v) == 12
end

# linspace with very large endpoints
for T = (Float32, Float64)
    a = realmax()
    for i = 1:5
        @test [Ranges.linspace(a,a,1);] == [a]
        @test [Ranges.linspace(-a,-a,1);] == [-a]
        b = realmax()
        for j = 1:5
            @test [Ranges.linspace(-a,b,0);] == []
            @test [Ranges.linspace(-a,b,2);] == [-a,b]
            @test [Ranges.linspace(-a,b,3);] == [-a,(b-a)/2,b]
            @test [Ranges.linspace(a,-b,0);] == []
            @test [Ranges.linspace(a,-b,2);] == [a,-b]
            @test [Ranges.linspace(a,-b,3);] == [a,(a-b)/2,-b]
            for c = maxintfloat(T)-3:maxintfloat(T)
                s = Ranges.linspace(-a,b,c)
                @test first(s) == -a
                @test last(s) == b
                c <= typemax(Int) && @test length(s) == c
                @test s.len == c
                s = Ranges.linspace(a,-b,c)
                @test first(s) == a
                @test last(s) == -b
                c <= typemax(Int) && @test length(s) == c
                @test s.len == c
            end
            b = prevfloat(b)
        end
        a = prevfloat(a)
    end
end

# linspace with 1 or 0 elements (whose step length is NaN)
@test issorted(Ranges.linspace(1,1,0))
@test issorted(Ranges.linspace(1,1,1))

# comparing and hashing ranges
let
    Rs = Range[1:2, map(Int32,1:3:17), map(Int64,1:3:17), 1:0, 17:-3:0,
               0.0:0.1:1.0, map(Float32,0.0:0.1:1.0),
               Ranges.linspace(0, 1, 20), map(Float32, Ranges.linspace(0, 1, 20))]
    for r in Rs
        ar = collect(r)
        @test r != ar
        @test !isequal(r,ar)
        for s in Rs
            as = collect(s)
            @test !isequal(r,s) || hash(r)==hash(s)
            @test (r==s) == (ar==as)
        end
    end
end

r = Ranges.linspace(1/3,5/7,6)
@test length(r) == 6
@test r[1] == 1/3
@test abs(r[end] - 5/7) <= eps(5/7)
r = Ranges.linspace(0.25,0.25,1)
@test length(r) == 1
@test_throws ErrorException Ranges.linspace(0.25,0.5,1)

# Issue #11245
let io = IOBuffer()
    show(io, Ranges.linspace(1, 2, 3))
    @static if VERSION < v"0.6-"
        str = takebuf_string(io)
    else
        str = String(take!(io))
    end
    @test str == "Ranges.linspace(1.0,2.0,3)"
end

# stringmime/show should display the range or linspace nicely
# to test print_range in range.jl
replstrmime(x) = sprint((io,x) -> show(IOContext(io, limit=true), MIME("text/plain"), x), x)
@test stringmime("text/plain", Ranges.linspace(1,5,7)) == "7-element Ranges.LinSpace{Float64}:\n 1.0,1.66667,2.33333,3.0,3.66667,4.33333,5.0"
@test repr(Ranges.linspace(1,5,7)) == "Ranges.linspace(1.0,5.0,7)"
# next is to test a very large range, which should be fast because print_range
# only examines spacing of the left and right edges of the range, sufficient
# to cover the designated screen size.
@test replstrmime(Ranges.linspace(0,100, 10000)) == "10000-element Ranges.LinSpace{Float64}:\n 0.0,0.010001,0.020002,0.030003,0.040004,…,99.95,99.96,99.97,99.98,99.99,100.0"

@test sprint(io -> show(io,UnitRange(1,2))) == "1:2"
@test sprint(io -> show(io,StepRange(1,2,5))) == "1:2:5"

# Issue 11049 and related
@test promote(Ranges.linspace(0f0, 1f0, 3), Ranges.linspace(0., 5., 2)) ===
    (Ranges.linspace(0., 1., 3), Ranges.linspace(0., 5., 2))
@test convert(Ranges.LinSpace{Float64}, Ranges.linspace(0., 1., 3)) === Ranges.linspace(0., 1., 3)
@test convert(Ranges.LinSpace{Float64}, Ranges.linspace(0f0, 1f0, 3)) === Ranges.linspace(0., 1., 3)

@test promote(Ranges.linspace(0., 1., 3), 0:5) === (Ranges.linspace(0., 1., 3),
                                             Ranges.linspace(0., 5., 6))
@test convert(Ranges.LinSpace{Float64}, 0:5) === Ranges.linspace(0., 5., 6)
@test convert(Ranges.LinSpace{Float64}, 0:1:5) === Ranges.linspace(0., 5., 6)
@test convert(Ranges.LinSpace, 0:5) === Ranges.linspace(0., 5., 6)
@test convert(Ranges.LinSpace, 0:1:5) === Ranges.linspace(0., 5., 6)

function test_range_index(r, s)
    @test typeof(r[s]) == typeof(r)
    @test [r;][s] == [r[s];]
end
test_range_index(Ranges.linspace(0.1, 0.3, 3), 1:2)
test_range_index(Ranges.linspace(0.1, 0.3, 3), 1:0)
test_range_index(Ranges.linspace(1.0, 1.0, 1), 1:1)
test_range_index(Ranges.linspace(1.0, 1.0, 1), 1:0)
test_range_index(Ranges.linspace(1.0, 2.0, 0), 1:0)

function test_linspace_identity{T}(r::Ranges.LinSpace{T}, mr::Ranges.LinSpace{T})
    @test -r == mr
    @test -collect(r) == collect(mr)
    @test isa(-r, Ranges.LinSpace)

    oldversion || return

    @test 1 + r + (-1) == r
    @test 1 .+ r .+ (-1) == r
    @test isa(1 + r + (-1), Ranges.LinSpace)
    @test 1 - r - 1 == mr
    @test 1 .- r .- 1 == mr
    @test isa(1 - r - 1, Ranges.LinSpace)

    @test 1 * r * 1 == r
    @test 2 * r * T(0.5) == r
    @test isa(1 * r * 1, Ranges.LinSpace)
    @test r / 1 == r
    @test r / 2 * 2 == r
    @test r / T(0.5) * T(0.5) == r
    @test isa(r / 1, Ranges.LinSpace)

    @test (2 * collect(r) == collect(r * 2) == collect(2 * r) ==
           collect(r * T(2.0)) == collect(T(2.0) * r) ==
           collect(r / T(0.5)) == -collect(mr * T(2.0)))
end

test_linspace_identity(Ranges.linspace(1.0, 27.0, 10), Ranges.linspace(-1.0, -27.0, 10))
test_linspace_identity(Ranges.linspace(1f0, 27f0, 10), Ranges.linspace(-1f0, -27f0, 10))

test_linspace_identity(Ranges.linspace(1.0, 27.0, 0), Ranges.linspace(-1.0, -27.0, 0))
test_linspace_identity(Ranges.linspace(1f0, 27f0, 0), Ranges.linspace(-1f0, -27f0, 0))

test_linspace_identity(Ranges.linspace(1.0, 1.0, 1), Ranges.linspace(-1.0, -1.0, 1))
test_linspace_identity(Ranges.linspace(1f0, 1f0, 1), Ranges.linspace(-1f0, -1f0, 1))

@test reverse(Ranges.linspace(1.0, 27.0, 1275)) == Ranges.linspace(27.0, 1.0, 1275)
@test [reverse(Ranges.linspace(1.0, 27.0, 1275));] ==
    reverse([Ranges.linspace(1.0, 27.0, 1275);])

# PR 12200 and related
for _r in (1:2:100, 1:100, 1f0:2f0:100f0, 1.0:2.0:100.0,
           Ranges.linspace(1, 100, 10), Ranges.linspace(1f0, 100f0, 10))
    float_r = float(_r)
    big_r = @static oldversion ? big(_r) : big.(_r)
    @static if oldversion
        @test typeof(big_r).name === typeof(_r).name
    end
    if eltype(_r) <: AbstractFloat
        @test isa(float_r, typeof(_r))
        @test eltype(big_r) === BigFloat
    else
        @test isa(float_r, Range)
        @test eltype(float_r) <: AbstractFloat
        @test eltype(big_r) === BigInt
    end
end

@test_throws DimensionMismatch Ranges.linspace(1.,5.,5) + Ranges.linspace(1.,5.,6)
@test_throws DimensionMismatch Ranges.linspace(1.,5.,5) - Ranges.linspace(1.,5.,6)
@test_throws DimensionMismatch Ranges.linspace(1.,5.,5) .* Ranges.linspace(1.,5.,6)
@test_throws DimensionMismatch Ranges.linspace(1.,5.,5) ./ Ranges.linspace(1.,5.,6)

# linspace of other types
r = Ranges.linspace(0, 3//10, 4)
@test eltype(r) == Rational{Int}
@test r[2] === 1//10

a, b = 1.0, nextfloat(1.0)
ba, bb = BigFloat(a), BigFloat(b)
r = Ranges.linspace(ba, bb, 3)
@test eltype(r) == BigFloat
@test r[1] == a && r[3] == b
@test r[2] == (ba+bb)/2

a, b = rand(10), rand(10)
acopy = copy(a)
@static if oldversion
    ba, bb = big(a), big(b)
else
    ba, bb = big.(a), big.(b)
end
r = Ranges.linspace(a, b, 5)
@test r[1] == a && r[5] == b
for i = 2:4
    x = ((5-i)//4)*ba + ((i-1)//4)*bb
    @test r[i] ≈ Float64.(x)
end
@test first(r .+ 2) == a+2
@test first(2 .+ r) == a+2
@test first(r) === a && a == acopy

r = Ranges.linspace(3+im, 1+3im, 3)
@test r[2] === 2.0+2.0im
@test first(r .+ (5-3im)) === 8.0-2.0im
@test first((5-3im) .+ r) == 8-2im
@test first(r .- (5-3im)) == -2+4im
@test first((5-3im) .- r) == 2-4im

@test Ranges.logspace(0, 3, 4) == [1,10,100,1000]

if oldversion
    list = [Ranges.linspace(1,3,3), 0.0:0.1:0.3]
    @test eltype(list) == Ranges.LinSpace{Float64}
    @test isa(convert(Ranges.LinSpace, 0.0:0.1:0.3), Ranges.LinSpace{Float64})
end


text:
@testset "planetoid" begin
    graph = graphdata(Planetoid(), :cora)
    @test typeof(graph) == SimpleGraph{UInt32}
    @test nv(graph) == 2708
    @test ne(graph) == 5275

    train_X, train_y = traindata(Planetoid(), :cora)
    @test typeof(train_X) == SparseMatrixCSC{Float32,Int64}
    @test size(train_X) == (1433, 140)
    @test typeof(train_y) == SparseMatrixCSC{Int32,Int64}
    @test size(train_y) == (7, 140)

    test_X, test_y = testdata(Planetoid(), :cora)
    @test typeof(test_X) == SparseMatrixCSC{Float32,Int64}
    @test size(test_X) == (1433, 1000)
    @test typeof(test_y) == SparseMatrixCSC{Int32,Int64}
    @test size(test_y) == (7, 1000)

    all_X, all_y = alldata(Planetoid(), :cora)
    @test typeof(all_X) == SparseMatrixCSC{Float32,Int64}
    @test size(all_X) == (1433, 1708)
    @test typeof(all_y) == SparseMatrixCSC{Int32,Int64}
    @test size(all_y) == (7, 1708)

    raw = rawdata(Planetoid(), :cora)
    @test typeof(raw[:graph]) == Dict{Any,Any}
    @test typeof(raw[:train_X]) == SparseMatrixCSC{Float32,Int64}
    @test typeof(raw[:train_y]) == SparseMatrixCSC{Int32,Int64}
    @test typeof(raw[:test_X]) == SparseMatrixCSC{Float32,Int64}
    @test typeof(raw[:test_y]) == SparseMatrixCSC{Int32,Int64}
    @test typeof(raw[:all_X]) == SparseMatrixCSC{Float32,Int64}
    @test typeof(raw[:all_y]) == SparseMatrixCSC{Int32,Int64}

    meta = metadata(Planetoid(), :cora)
    @test meta.graph.num_V == nv(graph)
    @test meta.graph.num_E == ne(graph)
end


text:
module GCode

@enum POSITION absolute relative

export absolute, relative
export linear_moves, feed_rate, positioning, set_position

function print_code(io::IO, code::String, args::Dict; range=nothing, rfrom=1)
	function w(i)
		print(io, code)
		for k in filter(k->k!=';', keys(args))
			print(io, ' ', uppercase(k), args[k][i])
		end
		if haskey(args, ';')
			print(io, " ; ", args[';'][i])
		end
		print(io, "\n")
	end

	if range == nothing
		range=rfrom:length(dict_entry(args))
	end

	w.(range)
	return range
end

dict_entry(d; i=1) = d[collect(keys(d))[i]]
entry_isa(d, t; i=1) = isa(dict_entry(d;i=i), t)
arrayize_dict(d) = Dict(k=>[d[k]] for k in keys(d))
dict_of_arrays(d) = entry_isa(d, Array)

function linear_moves(io::IO, points::Dict)
	print_code(io, "G00", points, range=1:1)
	print_code(io, "G01", points, rfrom=2)
end

function feed_rate(io::IO, r)
	print_code(io, "G01", Dict('F'=>[r], ';'=>["Feed rate"]))
end

function positioning(io::IO, p::POSITION)
	if p == absolute
		print_code(io, "G90", Dict(';'=>["Absolute Positioning"]))
	else
		print_code(io, "G91", Dict(';'=>["Relative Positioning"]))
	end
end

function set_position(io::IO, points::Dict)
	if dict_of_arrays(points)
		print_code(io, "G92", points)
	else
		print_code(io, "G92", arrayize_dict(points))
	end
end

###
end


text:
#' ---
#' title: Facility Location
#' ---

#' **Originally Contributed by**: Mathieu Tanneau and Alexis Montoison

#+
using Random
using LinearAlgebra

using JuMP
import GLPK
using Plots

#' ## Uncapacitated facility location

#' ### Problem description
#'
#' We are given
#' * $M=\{1, \dots, m\}$ clients
#' * $N=\{ 1, \dots, n\}$ sites where a facility can be built
#'
#' **Decision variables**
#' Decision variables are split into two categories:
#' * Binary variable $y_{j}$ indicates whether facility $j$ is built or not
#' * Binary variable $x_{i, j}$ indicates whether client $i$ is assigned to facility $j$
#'
#' **Objective**
#' The objective is to minimize the total cost of serving all clients.
#' This costs breaks down into two components:
#' * Fixed cost of building a facility.
#' In this example, this cost is $f_{j} = 1, \ \forall j$.
#' * Cost of serving clients from the assigned facility.
#' In this example, the cost $c_{i, j}$
#' of serving client $i$ from facility $j$
#' is the Euclidean distance between the two.
#'
#' **Constraints**
#' * Each customer must be served by exactly one facility
#' * A facility cannot serve any client unless it is open

#' ### MILP formulation
#'
#' The problem can be formulated as the following MILP:
#'
#' \begin{align}
#' \min_{x, y} \ \ \ &
#' \sum_{i, j} c_{i, j} x_{i, j} +
#' \sum_{j} f_{j} y_{j}\\
#' s.t. &
#' \sum_{j} x_{i, j} = 1, && \forall i \in M\\
#' & x_{i, j} \leq y_{j}, && \forall i \in M, j \in N\\
#' & x_{i, j}, y_{j} \in \{0, 1\}, && \forall i \in M, j \in N
#' \end{align}
#'
#' where the first set of constraints ensures
#' that each client is served exactly once,
#' and the second set of constraints ensures
#' that no client is served from an unopened facility.

#' ### Problem data

#+
Random.seed!(314)

m = 12  # number of clients
n = 5  # number of facility locations

# Clients' locations
Xc = rand(m)
Yc = rand(m)

# Facilities' potential locations
Xf = rand(n)
Yf = rand(n)

# Fixed costs
f = ones(n);

# Distance
c = zeros(m, n)
for i in 1:m
    for j in 1:n
        c[i, j] = norm([Xc[i] - Xf[j], Yc[i] - Yf[j]], 2)
    end
end

#+
# Display the data
scatter(Xc, Yc, label = "Clients", markershape=:circle, markercolor=:blue)
scatter!(Xf, Yf, label="Facility",
    markershape=:square, markercolor=:white, markersize=6,
    markerstrokecolor=:red, markerstrokewidth=2
)

#' ### JuMP implementation

#+
# Create a JuMP model
ufl = Model(GLPK.Optimizer)

# Variables
@variable(ufl, y[1:n], Bin);
@variable(ufl, x[1:m, 1:n], Bin);

# Each client is served exactly once
@constraint(ufl, client_service[i in 1:m],
    sum(x[i, j] for j in 1:n) == 1
);

# A facility must be open to serve a client
@constraint(ufl, open_facility[i in 1:m, j in 1:n],
    x[i, j] <= y[j]
)

# Objective
@objective(ufl, Min, f'y + sum(c .* x));

#+
# Solve the uncapacitated facility location problem with GLPK
optimize!(ufl)
println("Optimal value: ", objective_value(ufl))

#' ### Visualizing the solution

#+
# The threshold 1e-5 ensure that edges between clients and facilities are drawn when x[i, j] ≈ 1.
x_ = value.(x) .> 1 - 1e-5
y_ = value.(y) .> 1 - 1e-5

# Display clients
p = scatter(Xc, Yc, markershape=:circle, markercolor=:blue, label=nothing)

# Show open facility
mc = [(y_[j] ? :red : :white) for j in 1:n]
scatter!(Xf, Yf,
    markershape=:square, markercolor=mc, markersize=6,
    markerstrokecolor=:red, markerstrokewidth=2,
    label=nothing
)

# Show client-facility assignment
for i in 1:m
    for j in 1:n
        if x_[i, j] == 1
           plot!([Xc[i], Xf[j]], [Yc[i], Yf[j]], color=:black, label=nothing)
        end
    end
end

display(p)

#' ## Capacitated Facility location

#' ### Problem formulation
#'
#' The capacitated variant introduces a capacity constraint on each facility, i.e., clients have a certain level of demand to be served, while each facility only has finite capacity which cannot be exceeded.
#'
#' Specifically, let
#' * $a_{i} \geq 0$ denote the demand of client $i$
#' * $q_{j} \geq 0$ denote the capacity of facility $j$
#'
#' The capacity constraints then write
#' \begin{align}
#' \sum_{i} a_{i} x_{i, j} &\leq q_{j} y_{j} && \forall j \in N
#' \end{align}
#'
#' Note that, if $y_{j}$ is set to $0$, the capacity constraint above automatically forces $x_{i, j}$ to $0$.

#' Thus, the capacitated facility location can be formulated as follows
#'
#' \begin{align}
#' \min_{x, y} \ \ \ &
#' \sum_{i, j} c_{i, j} x_{i, j} +
#' \sum_{j} f_{j} y_{j}\\
#' s.t. &
#' \sum_{j} x_{i, j} = 1, && \forall i \in M\\
#' & \sum_{i} a_{i} x_{i, j} \leq q_{j} y_{j}, && \forall j \in N\\
#' & x_{i, j}, y_{j} \in \{0, 1\}, && \forall i \in M, j \in N
#' \end{align}
#'
#' For simplicity, we will assume that there is enough capacity to serve the demand,
#'  i.e., there exists at least one feasible solution.

#+
# Demands
a = rand(1:3, m);

# Capacities
q = rand(5:10, n);

#+
# Display the data
scatter(Xc, Yc, label=nothing,
    markershape=:circle, markercolor=:blue, markersize= 2 .*(2 .+ a)
)

scatter!(Xf, Yf, label=nothing,
    markershape=:rect, markercolor=:white, markersize= q,
    markerstrokecolor=:red, markerstrokewidth=2
)

#' ### JuMP implementation

#+
# Create a JuMP model
cfl = Model(GLPK.Optimizer)

# Variables
@variable(cfl, y[1:n], Bin);
@variable(cfl, x[1:m, 1:n], Bin);

# Each client is served exactly once
@constraint(cfl, client_service[i in 1:m], sum(x[i, :]) == 1)

# Capacity constraint
@constraint(cfl, capacity, x'a .<= (q .* y))

# Objective
@objective(cfl, Min, f'y + sum(c .* x));

#+
# Solve the problem
optimize!(cfl)
println("Optimal value: ", objective_value(cfl))

#' ### Visualizing the solution

#+
# The threshold 1e-5 ensure that edges between clients and facilities are drawn when x[i, j] ≈ 1.
x_ = value.(x) .> 1 - 1e-5
y_ = value.(y) .> 1 - 1e-5

# Display the solution
p = scatter(Xc, Yc, label=nothing,
    markershape=:circle, markercolor=:blue, markersize= 2 .*(2 .+ a)
)

mc = [(y_[j] ? :red : :white) for j in 1:n]
scatter!(Xf, Yf, label=nothing,
    markershape=:rect, markercolor=mc, markersize=q,
    markerstrokecolor=:red, markerstrokewidth=2
)

# Show client-facility assignment
for i in 1:m
    for j in 1:n
        if x_[i, j] == 1
            plot!([Xc[i], Xf[j]], [Yc[i], Yf[j]], color=:black, label=nothing)
            break
        end
    end
end

display(p)

#' ## Further
#' * [Benders decomposition](https://github.com/JuliaOpt/JuMPTutorials.jl/blob/master/script/optimization_concepts/benders_decomposition.jl)
#' is a method of choice for solving facility location problems.
#' * Benchmark instances can be found
#' [here](https://resources.mpi-inf.mpg.de/departments/d1/projects/benchmarks/UflLib/).


text:
"""
    transform(eq::Equation) -> DifferentialEquations.ODEFunction

Transform `eq` into an ordinary differential equation (ODE) defined in terms of the Boltzmann variable ϕ.

Returns an ODE with independent variable ϕ and two components, where the first is the solution itself and
the second component is the ϕ-derivative of the solution. The ODE is optimized for components stored in
`StaticArrays.SVector`s.

See also: [`DifferentialEquations`](https://diffeq.sciml.ai/stable/), [`StaticArrays.SVector`](https://juliaarrays.github.io/StaticArrays.jl/stable/pages/api/#SVector-1)
"""
function transform(eq::DiffusionEquation{1})
    let D = eq.D
        function f((θ, dθ_dϕ), ::NullParameters, ϕ)
            try
                D_, dD_dθ = value_and_derivative(D, typeof(θ), θ)

                d²θ_dϕ² = -((ϕ/2 + dD_dθ*dθ_dϕ)/D_)*dθ_dϕ

                return @SVector [dθ_dϕ, d²θ_dϕ²]
            catch e
                isa(e, DomainError) || rethrow()
                return @SVector [dθ_dϕ, typeof(dθ_dϕ)(NaN)]
            end
        end
        return ODEFunction{false}(f, syms=(eq.symbol, :d_dϕ))
    end
end

function transform(eq::DiffusionEquation{m}) where m
    @assert 2 ≤ m ≤ 3
    let D = eq.D, k = m-1
        function f((θ, dθ_dϕ), ::NullParameters, ϕ)
            try
                D_, dD_dθ = value_and_derivative(D, typeof(θ), θ)

                d²θ_dϕ² = -((ϕ/2 + dD_dθ*dθ_dϕ)/D_ + k/ϕ)*dθ_dϕ

                return @SVector [dθ_dϕ, d²θ_dϕ²]
            catch e
                isa(e, DomainError) || rethrow()
                return @SVector [dθ_dϕ, typeof(dθ_dϕ)(NaN)]
            end
        end
        return ODEFunction{false}(f, syms=(eq.symbol, :d_dϕ))
    end
end

d_dϕ(eq::DiffusionEquation{2}, θ, ϕ; flux_mul_r) = flux_mul_r/(-eq.D(θ)*ϕ)

function transform(eq::RichardsEquation{1})
    let C = eq.C, K = eq.K
        function f((h, dh_dϕ), ::NullParameters, ϕ)
            try
                K_, dK_dh = value_and_derivative(K, typeof(h), h)

                d²h_dϕ² = -((C(h)*ϕ/2 + dK_dh*dh_dϕ)/K_)*dh_dϕ

                return @SVector [dh_dϕ, d²h_dϕ²]
            catch e
                isa(e, DomainError) || rethrow()
                return @SVector [dh_dϕ, typeof(dh_dϕ)(NaN)]
            end
        end
        return ODEFunction{false}(f, syms=(eq.symbol, :d_dϕ))
    end
end

function transform(eq::RichardsEquation{m}) where m
    @assert 2 ≤ m ≤ 3
    let C = eq.C, K = eq.K, k=m-1
        function f((h, dh_dϕ), ::NullParameters, ϕ)
            try
                K_, dK_dh = value_and_derivative(K, typeof(h), h)

                d²h_dϕ² = -((C(h)*ϕ/2 + dK_dh*dh_dϕ)/K_ + k/ϕ)*dh_dϕ

                return @SVector [dh_dϕ, d²h_dϕ²]
            catch e
                isa(e, DomainError) || rethrow()
                return @SVector [dh_dϕ, typeof(dh_dϕ)(NaN)]
            end
        end
        return ODEFunction{false}(f, syms=(eq.symbol, :d_dϕ))
    end
end

d_dϕ(eq::RichardsEquation{2}, h, ϕ; flux_mul_r) = flux_mul_r/(-eq.K(h)*ϕ)


text:
module InterfaceTest

using JET.JETInterface, JET, Test
import JET: get_reports, BasicPass, GlobalUndefVarErrorReport

# customized report pass
# ======================

struct IgnoreAllPass <: ReportPass end
(::IgnoreAllPass)(::Type{<:InferenceErrorReport}, @nospecialize(_...)) = return
let
    result = @report_call report_pass=IgnoreAllPass() sum("julia")
    @test isempty(get_reports(result))
end

struct IgnoreAllExceptGlobalUndefVarPass <: ReportPass end
(::IgnoreAllExceptGlobalUndefVarPass)(::Type{<:InferenceErrorReport}, @nospecialize(_...)) = return
(::IgnoreAllExceptGlobalUndefVarPass)(::Type{GlobalUndefVarErrorReport}, @nospecialize(args...)) = BasicPass()(GlobalUndefVarErrorReport, args...)
let
    result = report_call(; report_pass=IgnoreAllExceptGlobalUndefVarPass()) do
        sum("julia") # should be ignored
        undefvar
    end
    @test length(get_reports(result)) == 1
    let r = first(get_reports(result))
        @test isa(r, GlobalUndefVarErrorReport)
        @test r.name === :undefvar
    end
end

# AbstractAnalyzer API
# ====================

# NOTE actual functionalities are tested by documentation

struct APIValidator <: AbstractAnalyzer
    state::AnalyzerState
end

function compute_sins(i)
    n = 1000000
    out = Vector{Float64}(undef, n)
    Threads.@threads for itr in 1:n
        i = itr
        out[i] = sin(i)
    end
    return out
end

@test_throws ErrorException @report_call analyzer=APIValidator compute_sins(10)

# interface 1: `APIValidator(; jetconfigs...) -> APIValidator`
APIValidator(; jetconfigs...) = APIValidator(AnalyzerState(; jetconfigs...))

@test_throws ErrorException @report_call analyzer=APIValidator compute_sins(10)

# interface 2: `AnalyzerState(analyzer::APIValidator) -> AnalyzerState`
JETInterface.AnalyzerState(analyzer::APIValidator) = analyzer.state

@test_throws ErrorException @report_call analyzer=APIValidator compute_sins(10)

# interface 3: `AbstractAnalyzer(analyzer::APIValidator, state::AnalyzerState) -> APIValidator`
JETInterface.AbstractAnalyzer(analyzer::APIValidator, state::AnalyzerState) = APIValidator(state)

@test_throws ErrorException @report_call analyzer=APIValidator compute_sins(10)

# interface 4: `ReportPass(analyzer::APIValidator) -> ReportPass`
JETInterface.ReportPass(analyzer::APIValidator) = IgnoreAllPass()

@test_throws ErrorException @report_call analyzer=APIValidator compute_sins(10)

# interface 5: `get_cache_key(analyzer::APIValidator) -> UInt`
JETInterface.get_cache_key(analyzer::APIValidator) = AnalyzerState(analyzer).param_key

# because `APIValidator` uses `IgnoreAllPass`, we won't get any reports
let
    result = @report_call analyzer=APIValidator compute_sins(10)
    @test isempty(get_reports(result))
end

end # module InterfaceTest