text:
module ConstraintSolver

using MatrixNetworks
using JSON
using MathOptInterface
using JuMP: @variable, @constraint, @objective, Model, with_optimizer, VariableRef, backend
import JuMP.sense_to_set

const MOI = MathOptInterface
const MOIU = MOI.Utilities

include("options.jl")

CS = ConstraintSolver

mutable struct Variable
    idx                 :: Int
    lower_bound         :: Int # inital lower and
    upper_bound         :: Int # upper bound of the variable see min, max otherwise
    first_ptr           :: Int
    last_ptr            :: Int
    values              :: Vector{Int}
    indices             :: Vector{Int}
    offset              :: Int
    min                 :: Int # the minimum value during the solving process
    max                 :: Int # for initial see lower/upper_bound
    changes             :: Vector{Vector{Tuple{Symbol,Int64,Int64,Int64}}}
    has_upper_bound     :: Bool # must be true to work
    has_lower_bound     :: Bool # must be true to work
    is_fixed            :: Bool
    is_integer          :: Bool # must be true to work

end

Variable(idx) = Variable(idx,0,0,0,0,[],[],0,0,0,Vector{Vector{Tuple{Symbol,Int64,Int64,Int64}}}(), false, false, false, false)

mutable struct CSInfo
    pre_backtrack_calls :: Int
    backtracked         :: Bool
    backtrack_fixes     :: Int
    in_backtrack_calls  :: Int
    backtrack_reverses  :: Int
end

abstract type Constraint
end

abstract type ObjectiveFunction
end

mutable struct SingleVariableObjective <: ObjectiveFunction
    fct     :: Function
    index   :: Int # index of the variable
    indices :: Vector{Int}
end

mutable struct BasicConstraint <: Constraint
    idx                 :: Int
    fct                 :: Function
    indices             :: Vector{Int}
    pvals               :: Vector{Int}
    hash                :: UInt64
    BasicConstraint() = new()
end

# support for a <= b constraint
mutable struct SingleVariableConstraint <: Constraint
    idx                 :: Int
    fct                 :: Function
    indices             :: Vector{Int}
    pvals               :: Vector{Int}
    lhs                 :: Int
    rhs                 :: Int
    hash                :: UInt64
    SingleVariableConstraint() = new()
end

struct AllDifferentSet <: MOI.AbstractVectorSet
    dimension :: Int64
end

struct NotEqualSet{T} <: MOI.AbstractScalarSet
    value :: T
end

mutable struct LinearVariables
    indices             :: Vector{Int}
    coeffs              :: Vector{Int}
end

mutable struct LinearConstraint <: Constraint
    idx                 :: Int
    fct                 :: Function
    indices             :: Vector{Int}
    pvals               :: Vector{Int}
    coeffs              :: Vector{Int}
    operator            :: Symbol
    rhs                 :: Int
    in_all_different    :: Bool
    mins                :: Vector{Int}
    maxs                :: Vector{Int}
    pre_mins            :: Vector{Int}
    pre_maxs            :: Vector{Int}
    hash                :: UInt64
    LinearConstraint() = new()
end

mutable struct BacktrackObj
    idx                 :: Int
    parent_idx          :: Int
    depth               :: Int
    status              :: Symbol
    variable_idx        :: Int
    pval                :: Int
    best_bound          :: Int

    BacktrackObj() = new()
end

mutable struct TreeLogNode
    id              :: Int
    status          :: Symbol
    best_bound      :: Int
    step_nr         :: Int
    var_idx         :: Int
    set_val         :: Int
    var_states      :: Dict{Int64,Vector{Int64}}
    var_changes     :: Dict{Int64,Vector{Tuple{Symbol, Int64, Int64, Int64}}}
    children        :: Vector{TreeLogNode}
    TreeLogNode() = new()
    TreeLogNode(id, status, best_bound, step_nr, var_idx, set_val, var_states, var_changes, children) = new(id, status, best_bound, step_nr, var_idx, set_val, var_states, var_changes, children)
end

mutable struct ConstraintSolverModel
    init_search_space   :: Vector{Variable}
    search_space        :: Vector{Variable}
    subscription        :: Vector{Vector{Int}}
    constraints         :: Vector{Constraint}
    bt_infeasible       :: Vector{Int}
    c_backtrack_idx     :: Int
    backtrack_vec       :: Vector{BacktrackObj}
    sense               :: MOI.OptimizationSense
    objective           :: ObjectiveFunction
    best_sol            :: Int # Objective of the best solution
    best_bound          :: Int # Overall best bound
    solutions           :: Vector{Int} # saves only the id to the BacktrackObj
    info                :: CSInfo
    input               :: Dict{Symbol,Any}
    logs                :: Vector{TreeLogNode}

    ConstraintSolverModel() = new()
end

const CoM = ConstraintSolverModel

include("MOI_wrapper/MOI_wrapper.jl")
include("printing.jl")
include("logs.jl")
include("hashes.jl")
include("Variable.jl")
include("objective.jl")
include("linearcombination.jl")
include("all_different.jl")
include("eq_sum.jl")
include("svc.jl")
include("equal.jl")
include("not_equal.jl")

"""
    init()

Initialize the constraint model object
"""
function init()
    com = CoM()
    com.constraints         = Vector{Constraint}()
    com.subscription        = Vector{Vector{Int}}()
    com.init_search_space   = Vector{Variable}()
    com.search_space        = Vector{Variable}()
    com.bt_infeasible       = Vector{Int}()
    com.c_backtrack_idx     = 1
    com.sense               = MOI.FEASIBILITY_SENSE
    com.backtrack_vec       = Vector{BacktrackObj}()
    com.solutions           = Vector{Int}()
    com.info                = CSInfo(0, false, 0, 0, 0)
    com.input               = Dict{Symbol, Bool}()
    com.logs                = Vector{TreeLogNode}()
    return com
end

"""
    add_var!(com::CS.CoM, from::Int, to::Int; fix=nothing)

Adding a variable to the constraint model `com`. The variable is discrete and has the possible values from,..., to.
If the variable should be fixed to something one can use the `fix` keyword i.e `add_var!(com, 1, 9; fix=5)`
"""
function add_var!(com::CS.CoM, from::Int, to::Int; fix=nothing)
    ind = length(com.search_space)+1
    changes = Vector{Vector{Tuple{Symbol,Int64,Int64,Int64}}}()
    push!(changes, Vector{Tuple{Symbol,Int64,Int64,Int64}}())
    var = Variable(ind, from, to, 1, to-from+1, from:to, 1:to-from+1, 1-from, from, to, changes, true, true, fix !== nothing, true)
    if fix !== nothing
        fix!(com, var, fix)
    end
    push!(com.search_space, var)
    push!(com.subscription, Int[])
    push!(com.bt_infeasible, 0)
    return var
end

"""
    fulfills_constraints(com::CS.CoM, index, value)

Return whether the model is still feasible after setting the variable at position `index` to `value`.
"""
function fulfills_constraints(com::CS.CoM, index, value)
    # variable doesn't have any constraint
    if index > length(com.subscription)
        return true
    end
    feasible = true
    for ci in com.subscription[index]
        constraint = com.constraints[ci]
        feasible = constraint.fct(com, constraint, value, index)
        if !feasible
            break
        end
    end
    return feasible
end

"""
    fixed_vs_unfixed(search_space, indices)

Return the fixed_vals as well as the unfixed_indices
"""
function fixed_vs_unfixed(search_space, indices)
    # get all values which are fixed
    fixed_vals = Int[]
    unfixed_indices = Int[]
    for (i,ind) in enumerate(indices)
        if isfixed(search_space[ind])
            push!(fixed_vals, CS.value(search_space[ind]))
        else
            push!(unfixed_indices, i)
        end
    end
    return fixed_vals, unfixed_indices
end

"""
    set_pvals!(com::CS.CoM, constraint::Constraint)

Compute the possible values inside this constraint and set it as constraint.pvals
"""
function set_pvals!(com::CS.CoM, constraint::Constraint)
    indices = constraint.indices
    variables = Variable[v for v in com.search_space[indices]]
    pvals_intervals = Vector{NamedTuple}()
    push!(pvals_intervals, (from = variables[1].lower_bound, to = variables[1].upper_bound))
    for (i,ind) in enumerate(indices)
        extra_from = variables[i].min
        extra_to   = variables[i].max
        comp_inside = false
        for cpvals in pvals_intervals
            if extra_from >= cpvals.from && extra_to <= cpvals.to
                # completely inside the interval already
                comp_inside = true
                break
            elseif extra_from >= cpvals.from && extra_from <= cpvals.to
                extra_from = cpvals.to+1
            elseif extra_to <= cpvals.to && extra_to >= cpvals.from
                extra_to = cpvals.from-1
            end
        end
        if !comp_inside && extra_to >= extra_from
            push!(pvals_intervals, (from = extra_from, to = extra_to))
        end
    end
    pvals = collect(pvals_intervals[1].from:pvals_intervals[1].to)
    for interval in pvals_intervals[2:end]
        pvals = vcat(pvals, collect(interval.from:interval.to))
    end
    constraint.pvals = pvals
end

"""
    add_constraint!(com::CS.CoM, constraint::Constraint)

Add a constraint to the model i.e `add_constraint!(com, a != b)`
"""
function add_constraint!(com::CS.CoM, constraint::Constraint)
    constraint.idx = length(com.constraints)+1
    push!(com.constraints, constraint)
    set_pvals!(com, constraint)
    for (i,ind) in enumerate(constraint.indices)
        push!(com.subscription[ind], constraint.idx)
    end
end

"""
    get_weak_ind(com::CS.CoM)

Get the next weak index for backtracking. This will be the next branching variable.
Return whether there is an unfixed variable and a best index
"""
function get_weak_ind(com::CS.CoM)
    lowest_num_pvals = typemax(Int)
    biggest_inf = -1
    best_ind = -1
    biggest_dependent = typemax(Int)
    found = false

    for ind in 1:length(com.search_space)
        if !isfixed(com.search_space[ind])
            num_pvals = nvalues(com.search_space[ind])
            inf = com.bt_infeasible[ind]
            if inf >= biggest_inf
                if inf > biggest_inf || num_pvals < lowest_num_pvals
                    lowest_num_pvals = num_pvals
                    biggest_inf = inf
                    best_ind = ind
                    found = true
                end
            end
        end
    end
    return found, best_ind
end

"""
    prune!(com::CS.CoM, prune_steps::Vector{Int})

Prune the search space based on a list of backtracking indices `prune_steps`.
"""
function prune!(com::CS.CoM, prune_steps::Vector{Int})
    search_space = com.search_space
    for backtrack_idx in prune_steps
        for var in search_space
            for change in var.changes[backtrack_idx]
                fct_symbol = change[1]
                val = change[2]
                if fct_symbol == :fix
                    fix!(com, var, val; changes=false)
                elseif fct_symbol == :rm
                    rm!(com, var, val; changes=false)
                elseif fct_symbol == :remove_above
                    remove_above!(com, var, val; changes=false)
                elseif fct_symbol == :remove_below
                    remove_below!(com, var, val; changes=false)
                else
                    throw(ErrorException("There is no pruning function for $fct_symbol"))
                end
            end
        end
    end
end

function open_possibilities(search_space, indices)
    open = 0
    for vi in indices
        if !isfixed(search_space[vi])
            open += nvalues(search_space[vi])
        end
    end
    return open
end

function find_best_constraint(com::CS.CoM, constraint_idxs_vec)
    best_ci = 0
    best_open = typemax(Int64)
    best_hash = typemax(UInt64)
    for ci = 1:length(constraint_idxs_vec)
        if constraint_idxs_vec[ci] <= best_open
            if constraint_idxs_vec[ci] < best_open || com.constraints[ci].hash < best_hash
                best_ci = ci
                best_open = constraint_idxs_vec[ci]
                best_hash = com.constraints[ci].hash
            end
        end
    end
    return best_open, best_ci
end

"""
    prune!(com::CS.CoM; pre_backtrack=false, all=false, only_once=false, initial_check=false)

Prune based on changes by initial solve or backtracking. The information for it is stored in each variable.
There are several parameters:
`pre_backtrack` when set to true `com.info.in_backtrack_calls` is incremented
`all` instead of only looking at changes each constraint is check at least once (if there are open values)
`only_once` Just run on the changed constraints or `all` once instead of repeatedly until nothing can be pruned
`initial_check` Checks on `all` constraints and also checks if variables are set for the whole constraint
    whether the constraint is fulfilled or the problem is infeasible.
Return whether it's still feasible
"""
function prune!(com::CS.CoM; pre_backtrack=false, all=false, only_once=false, initial_check=false)
    feasible = true
    N = typemax(Int64)
    search_space = com.search_space
    prev_var_length = zeros(Int, length(search_space))
    constraint_idxs_vec = fill(N, length(com.constraints))
    # get all constraints which need to be called (only once)
    current_backtrack_id = com.c_backtrack_idx
    for var in search_space
        new_var_length = length(var.changes[current_backtrack_id])
        if new_var_length > 0 || all || initial_check
            prev_var_length[var.idx] = new_var_length
            inner_constraints = com.constraints[com.subscription[var.idx]]
            for ci in com.subscription[var.idx]
                inner_constraint = com.constraints[ci]
                constraint_idxs_vec[inner_constraint.idx] = open_possibilities(search_space, inner_constraint.indices)
            end
        end
    end

    # while we haven't called every constraint
    while true
        b_open_constraint = false
        # will be changed or b_open_constraint => false
        open_pos, ci = find_best_constraint(com, constraint_idxs_vec)
        # no open values => don't need to call again
        if open_pos == 0 && !initial_check
            constraint_idxs_vec[ci] = N
            continue
        end
        # checked all
        if open_pos == N
            break
        end
        constraint_idxs_vec[ci] = N
        constraint = com.constraints[ci]

        feasible = constraint.fct(com, constraint; logs = false)
        if !pre_backtrack
            com.info.in_backtrack_calls += 1
        else
            com.info.pre_backtrack_calls += 1
        end

        if !feasible
            break
        end

        # if we changed another variable increase the level of the constraints to call them later
        for var_idx in constraint.indices
            var = search_space[var_idx]
            new_var_length = length(var.changes[current_backtrack_id])
            if new_var_length > prev_var_length[var.idx]
                prev_var_length[var.idx] = new_var_length
                inner_constraints = com.constraints[com.subscription[var.idx]]
                for ci in com.subscription[var.idx]
                    if ci != constraint.idx
                        inner_constraint = com.constraints[ci]
                        # if initial check or don't add constraints => update only those which already have open possibilities
                        if (only_once || initial_check) && constraint_idxs_vec[inner_constraint.idx] == N
                            continue
                        end
                        constraint_idxs_vec[inner_constraint.idx] = open_possibilities(search_space, inner_constraint.indices)
                    end
                end
            end
        end
    end
    return feasible
end

"""
    single_reverse_pruning!(search_space, index::Int, prune_int::Int, prune_fix::Int)

Reverse a single variable using `prune_int` (number of value removals) and `prune_fix` (new last_ptr if not 0).
"""
function single_reverse_pruning!(search_space, index::Int, prune_int::Int, prune_fix::Int)
    if prune_int > 0
        var = search_space[index]
        l_ptr = max(1,var.last_ptr)

        new_l_ptr = var.last_ptr + prune_int
        min_val, max_val = extrema(var.values[l_ptr:new_l_ptr])
        if min_val < var.min
            var.min = min_val
        end
        if max_val > var.max
            var.max = max_val
        end
        var.last_ptr = new_l_ptr
    end
    if prune_fix > 0
        var = search_space[index]
        var.last_ptr = prune_fix

        min_val, max_val = extrema(var.values[1:prune_fix])
        if min_val < var.min
            var.min = min_val
        end
        if max_val > var.max
            var.max = max_val
        end
        var.first_ptr = 1
    end
end

"""
    reverse_pruning!(com, backtrack_idx)

Reverse the changes made by a specific backtrack object
"""
function reverse_pruning!(com::CS.CoM, backtrack_idx::Int)
    search_space = com.search_space
    for var in search_space
        v_idx = var.idx

        for change in Iterators.reverse(var.changes[backtrack_idx])
            single_reverse_pruning!(search_space, v_idx, change[4], change[3])
        end
    end
end

"""
    get_best_bound(com::CS.CoM; var_idx=0, val=0)

Return the best bound if setting the variable with idx: `var_idx` to `val` if `var_idx != 0`.
Without an objective function return 0.
"""
function get_best_bound(com::CS.CoM; var_idx=0, val=0)
    if com.sense == MOI.FEASIBILITY_SENSE
        return 0
    end
    return com.objective.fct(com, var_idx, val)
end

function checkout_from_to!(com::CS.CoM, from_idx::Int, to_idx::Int)
    backtrack_vec = com.backtrack_vec
    from = backtrack_vec[from_idx]
    to = backtrack_vec[to_idx]
    if to.parent_idx == from.idx
        return
    end
    reverse_pruning!(com, from.idx)

    prune_steps = Vector{Int}()
    # first go to same level if new is higher in the tree
    if to.depth < from.depth
        depth = from.depth
        parent_idx = from.parent_idx
        parent = backtrack_vec[parent_idx]
        while to.depth < depth
            reverse_pruning!(com, parent_idx)
            parent = backtrack_vec[parent_idx]
            parent_idx = parent.parent_idx
            depth -= 1
        end
        if parent_idx == to.parent_idx
            return
        else
            from = parent
        end
    elseif from.depth < to.depth
        depth = to.depth
        parent_idx = to.parent_idx
        parent = backtrack_vec[parent_idx]
        while from.depth < depth
            pushfirst!(prune_steps, parent_idx)
            parent = backtrack_vec[parent_idx]
            parent_idx = parent.parent_idx
            depth -= 1
        end

        to = parent
        if backtrack_vec[prune_steps[1]].parent_idx == from.parent_idx
            prune!(com, prune_steps)
            return
        end
    end
    @assert from.depth == to.depth
    # same diff but different parent
    # => level up until same parent
     while from.parent_idx != to.parent_idx
        reverse_pruning!(com, from.parent_idx)
        from = backtrack_vec[from.parent_idx]

        pushfirst!(prune_steps, to.parent_idx)
        to = backtrack_vec[to.parent_idx]
    end

    prune!(com, prune_steps)
end

"""
    update_best_bound!(backtrack_obj::BacktrackObj, com::CS.CoM, constraints)

Check all constraints which change the objective and update the best bound of the backtrack_obj accordingly.
Pruning should not be continued if the new best bound has changed.
Return feasible and if pruning should be continued.
"""
function update_best_bound!(backtrack_obj::BacktrackObj, com::CS.CoM, constraints)
    further_pruning = true
    feasible = true
    for constraint in constraints
        relevant = false
        for obj_index in com.objective.indices
            if obj_index in constraint.indices
                relevant = true
                break
            end
        end
        if relevant
            feasible = constraint.fct(com, constraint; logs = false)
            if !feasible
                return false, false
            end
        end
    end
    if !feasible
        com.info.backtrack_reverses += 1
        return false, false
    end

    # check best_bound again
    # if best bound unchanged => continue pruning
    # otherwise try another path but don't close the current
    # -> means open new paths from here even if not pruned til the end
    new_bb = get_best_bound(com)
    if backtrack_obj.best_bound != new_bb
        further_pruning = false
    end
    return true, further_pruning
end

"""
    backtrack!(com::CS.CoM, max_bt_steps; sorting=true)

Start backtracking and stop after `max_bt_steps`.
If `sorting` is set to `false` the same ordering is used as when used without objective this has only an effect when an objective is used.
Return :Solved or :Infeasible if proven or `:NotSolved` if interrupted by `max_bt_steps`.
"""
function backtrack!(com::CS.CoM, max_bt_steps; sorting=true)
    found, ind = get_weak_ind(com)
    com.info.backtrack_fixes   = 1

    pvals = reverse!(values(com.search_space[ind]))
    dummy_backtrack_obj = BacktrackObj()
    dummy_backtrack_obj.status = :Close
    dummy_backtrack_obj.idx = 1
    dummy_backtrack_obj.variable_idx = 1

    backtrack_vec = com.backtrack_vec
    push!(backtrack_vec, dummy_backtrack_obj)

    # the first solve (before backtrack) has idx 1
    num_backtrack_objs = 1
    step_nr = 1
    for pval in pvals
        backtrack_obj = BacktrackObj()
        num_backtrack_objs += 1
        backtrack_obj.idx = num_backtrack_objs
        backtrack_obj.parent_idx = 1
        backtrack_obj.depth = 1
        backtrack_obj.status = :Open
        backtrack_obj.best_bound = get_best_bound(com; var_idx=ind, val=pval)
        backtrack_obj.variable_idx = ind
        backtrack_obj.pval = pval
        push!(backtrack_vec, backtrack_obj)
        for v in com.search_space
            push!(v.changes, Vector{Tuple{Symbol,Int64,Int64,Int64}}())
        end
        if com.input[:logs]
            push!(com.logs, log_one_node(com, length(com.search_space), num_backtrack_objs, -1))
        end
    end
    last_backtrack_id = 0

    started = true
    obj_factor = com.sense == MOI.MIN_SENSE ? 1 : -1

    while length(backtrack_vec) > 0
        step_nr += 1
        # get next open backtrack object
        l = 1

        # if there is no objective or sorting is set to false
        if com.sense == MOI.FEASIBILITY_SENSE || !sorting
            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
        else # 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(Int64)
            best_depth = 0
            found_sol = length(com.solutions) > 0
            nopen_nodes = 0
            for i=1:length(backtrack_vec)
                bo = backtrack_vec[i]
                if bo.status == :Open
                    nopen_nodes += 1
                    if found_sol
                        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
                    else
                        if bo.depth > best_depth || (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
            end

            if l != 0
                backtrack_obj = backtrack_vec[l]
            end
        end

        # no open node => Infeasible
        if l <= 0 || l > length(backtrack_vec)
            break
        end
        if com.sense == MOI.MIN_SENSE
            max_val = typemax(Int64)
            com.best_bound = minimum([bo.status == :Open ? bo.best_bound : max_val for bo in backtrack_vec])
        elseif com.sense == MOI.MAX_SENSE
            min_val = typemin(Int64)
            com.best_bound = maximum([bo.status == :Open ? bo.best_bound : min_val for bo in backtrack_vec])
        else
            com.best_bound = 0
        end

        ind = backtrack_obj.variable_idx

        com.c_backtrack_idx = backtrack_obj.idx

        if !started
            com.c_backtrack_idx = 0
            checkout_from_to!(com, last_backtrack_id, backtrack_obj.idx)
            com.c_backtrack_idx = backtrack_obj.idx
        end

        if com.input[:logs]
            com.logs[backtrack_obj.idx].step_nr = step_nr
        end

        started = false
        last_backtrack_id = backtrack_obj.idx

        backtrack_obj.status = :Closed

        pval = backtrack_obj.pval

        # check if this value is still possible
        constraints = com.constraints[com.subscription[ind]]
        feasible = true
        for constraint in constraints
            feasible = constraint.fct(com, constraint, pval, ind)
            if !feasible
                break
            end
        end
        if !feasible
            continue
        end
        # value is still possible => set it
        fix!(com, com.search_space[ind], pval)
        com.info.backtrack_fixes   += 1

        further_pruning = true
        # first update the best bound (only constraints which have an index in the objective function)
        if com.sense != MOI.FEASIBILITY_SENSE
            feasible, further_pruning = update_best_bound!(backtrack_obj, com, constraints)

            if !feasible
                com.info.backtrack_reverses += 1
                continue
            end
        end

        if further_pruning
            # prune completely start with all that changed by the fix or by updating best bound
            feasible = prune!(com)

            if !feasible
                com.info.backtrack_reverses += 1
                continue
            end
        end

        found, ind = get_weak_ind(com)
        # no index found => solution found
        if !found
            new_sol = get_best_bound(com)
            if length(com.solutions) == 0
                com.best_sol = new_sol
            elseif obj_factor*new_sol < obj_factor*com.best_sol
                com.best_sol = new_sol
            end
            push!(com.solutions, backtrack_obj.idx)
            if com.best_sol == com.best_bound
                return :Solved
            else
                # set all nodes to :Worse if they can't achieve a better solution
                for bo in backtrack_vec
                    if bo.status == :Open && obj_factor*bo.best_bound >= com.best_sol
                        bo.status = :Worse
                    end
                end
                continue
            end
        end

        if com.info.backtrack_fixes > max_bt_steps
            return :NotSolved
        end

        if com.input[:logs]
            com.logs[backtrack_obj.idx] = log_one_node(com, length(com.search_space), backtrack_obj.idx, step_nr)
        end

        pvals = reverse!(values(com.search_space[ind]))
        last_backtrack_obj = backtrack_vec[last_backtrack_id]
        for pval in pvals
            backtrack_obj = BacktrackObj()
            num_backtrack_objs += 1
            backtrack_obj.parent_idx = last_backtrack_obj.idx
            backtrack_obj.depth = last_backtrack_obj.depth + 1
            backtrack_obj.idx = num_backtrack_objs
            backtrack_obj.status = :Open
            backtrack_obj.best_bound = get_best_bound(com; var_idx=ind, val=pval)
            backtrack_obj.variable_idx = ind
            backtrack_obj.pval = pval

            # only include nodes which have a better objective than the current best solution if one was found already
            if backtrack_obj.best_bound*obj_factor < com.best_sol || length(com.solutions) == 0
                push!(backtrack_vec, backtrack_obj)
                for v in com.search_space
                    push!(v.changes, Vector{Tuple{Symbol,Int64,Int64,Int64}}())
                end
                if com.input[:logs]
                    push!(com.logs, log_one_node(com, length(com.search_space), num_backtrack_objs, -1))
                end
            else
                num_backtrack_objs -= 1
            end
        end
    end
    if length(com.solutions) > 0
        # find one of the best solutions
        sol, sol_id = findmin([backtrack_vec[sol_id].best_bound*obj_factor for sol_id in com.solutions])
        backtrack_id = com.solutions[sol_id]
        checkout_from_to!(com, last_backtrack_id, backtrack_id)
        # prune the last step as checkout_from_to! excludes the to part
        prune!(com, [backtrack_id])
        return :Solved
    end
    return :Infeasible
end

function arr2dict(arr)
    d = Dict{Int,Bool}()
    for v in arr
        d[v] = true
    end
    return d
end

"""
    simplify!(com)

Simplify constraints i.e by adding new constraints which uses an implicit connection between two constraints.
i.e an `all_different` does sometimes include information about the sum.
"""
function simplify!(com)
    # check if we have all_different and sum constraints
    # (all different where every value is used)
    b_all_different = false
    b_all_different_sum = false
    b_sum = false
    for constraint in com.constraints
        if nameof(constraint.fct) == :all_different
            b_all_different = true
            if length(constraint.indices) == length(constraint.pvals)
                b_all_different_sum = true
            end
        elseif nameof(constraint.fct) == :eq_sum
            b_sum = true
        end
    end
    if b_all_different_sum && b_sum
        # for each all_different constraint
        # which can be formulated as a sum constraint
        # check which sum constraints are completely inside all different
        # which are partially inside
        # compute inside sum and total sum
        n_constraints_before = length(com.constraints)
        for constraint_idx in 1:length(com.constraints)
            constraint = com.constraints[constraint_idx]

            if nameof(constraint.fct) == :all_different
                add_sum_constraint = true
                if length(constraint.indices) == length(constraint.pvals)
                    all_diff_sum = sum(constraint.pvals)
                    in_sum = 0
                    total_sum = 0
                    outside_indices = Int[]
                    cons_indices_dict = arr2dict(constraint.indices)
                    for variable_idx in keys(cons_indices_dict)
                        found_sum_constraint = false
                        for sub_constraint_idx in com.subscription[variable_idx]
                            # don't mess with constraints added later on
                            if sub_constraint_idx > n_constraints_before
                                continue
                            end
                            sub_constraint = com.constraints[sub_constraint_idx]
                            # it must be an equal constraint and all coefficients must be 1 otherwise we can't add a constraint
                            if nameof(sub_constraint.fct) == :eq_sum && all(c->c==1, sub_constraint.coeffs)
                                found_sum_constraint = true
                                total_sum += sub_constraint.rhs
                                all_inside = true
                                for sub_variable_idx in sub_constraint.indices
                                    if !haskey(cons_indices_dict, sub_variable_idx)
                                        all_inside = false
                                        push!(outside_indices, sub_variable_idx)
                                    else
                                        delete!(cons_indices_dict, sub_variable_idx)
                                    end
                                end
                                if all_inside
                                    in_sum += sub_constraint.rhs
                                end
                                break
                            end
                        end
                        if !found_sum_constraint
                            add_sum_constraint = false
                            break
                        end
                    end

                    # make sure that there are not too many outside indices
                    if add_sum_constraint && length(outside_indices) < 3
                        add_constraint!(com, sum(com.search_space[outside_indices]) == total_sum - all_diff_sum)
                    end
                end
            end
        end
    end
end

"""
    set_in_all_different!(com::CS.CoM)

Set `constraint.in_all_different` if all variables in the constraint are part of the same `all_different` constraint.
"""
function set_in_all_different!(com::CS.CoM)
    for constraint in com.constraints
        if :in_all_different in fieldnames(typeof(constraint))
            if !constraint.in_all_different
                subscriptions_idxs = [[i for i in com.subscription[v]] for v in constraint.indices]
                intersects = intersect(subscriptions_idxs...)

                for i in intersects
                    if nameof(com.constraints[i].fct) == :all_different
                        constraint.in_all_different = true
                        break
                    end
                end
            end
        end
    end
end

"""
    solve!(com::CS.CoM, options::SolverOptions)

Solve the constraint model based on the given settings.
"""
function solve!(com::CS.CoM, options::SolverOptions)
    # no constraints => Solved
    if length(com.constraints) == 0
        return :Solved
    end

    backtrack = options.backtrack
    max_bt_steps = options.max_bt_steps
    backtrack_sorting = options.backtrack_sorting
    keep_logs = options.keep_logs

    com.input[:logs] = keep_logs
    if keep_logs
        com.init_search_space = deepcopy(com.search_space)
    end

    set_in_all_different!(com)

    # check for better constraints
    simplify!(com)

    # check if all feasible even if for example everything is fixed
    feasible = prune!(com; pre_backtrack=true, initial_check=true)

    if !feasible
        return :Infeasible
    end

    if all(v->isfixed(v), com.search_space)
        com.best_bound = get_best_bound(com)
        com.best_sol = com.best_bound
        return :Solved
    end
    feasible = prune!(com; pre_backtrack=true)

    com.best_bound = get_best_bound(com)
    if keep_logs
        push!(com.logs, log_one_node(com, length(com.search_space), 1, 1))
    end

    if !feasible
        return :Infeasible
    end

    if all(v->isfixed(v), com.search_space)
        com.best_sol = com.best_bound
        return :Solved
    end
    if backtrack
        com.info.backtracked = true
        return backtrack!(com, max_bt_steps; sorting=backtrack_sorting)
    else
        @info "Backtracking is turned off."
        return :NotSolved
    end
end

"""
    solve!(com::Optimizer)

Solve the constraint model based on the given settings.
"""
function solve!(model::Optimizer)
    com = model.inner
    return solve!(com, model.options)
end

end # module


text:
using Mux
using Test
using Lazy
using HTTP
import HTTP.ExceptionRequest: StatusError

@test Mux.notfound()(d())[:status] == 404
d1 = Dict("one"=> "1", "two"=> "2")
d2 = Dict("one"=> "1", "two"=> "")
# Test basic server
@app test = (
  Mux.defaults,
  page("/",respond("
Hello World!
")),
  page("/about", respond("
Boo!
")),
  page("/user/:user", req -> "
Hello, $(req[:params][:user])!
"),
  query(d1, respond("
query1
")),
  query(d2, respond("
query2
")),
  Mux.notfound())
serve(test)
@test String(HTTP.get("http://localhost:8000").body) ==
            "
Hello World!
"
@test String(HTTP.get("http://localhost:8000/about").body) ==
            "
Boo!
"
@test String(HTTP.get("http://localhost:8000/user/julia").body) ==
            "
Hello, julia!
"

# Issue #68
@test Mux.fileheaders("foo.css")["Content-Type"] == "text/css"
@test Mux.fileheaders("foo.html")["Content-Type"] == "text/html"
@test Mux.fileheaders("foo.js")["Content-Type"] == "application/javascript"

function f()
  @app foo = (Mux.defaults)
end

@test f() == nothing

# Query based routing
@test String(HTTP.get("http://localhost:8000/dum?one=1&two=2").body) ==
            "
query1
"
@test_throws StatusError String(HTTP.get("http://localhost:8000/dum?one=1").body)
@test_throws StatusError String(HTTP.get("http://localhost:8000/dum?one=1&two=2&sarv=boo").body)
@test_throws StatusError String(HTTP.get("http://localhost:8000/dum?one=1").body)
@test String(HTTP.get("http://localhost:8000/dum?one=1&two=56").body) ==
            "
query2
"
@test String(HTTP.get("http://localhost:8000/dum?one=1&two=hfjd").body) ==
            "
query2
"
@test_throws StatusError String(HTTP.get("http://localhost:8000/dum?one=1&two=2&sarv=boo").body)

throwapp() = (_...) -> error("An error!")

# Test production server
@app test = (
  Mux.prod_defaults,
  page("/",respond("
Hello World!
")),
  page("/about", respond("
Boo!
")),
  page("/user/:user", req -> "
Hello, $(req[:params][:user])!
"),
  throwapp(),
  Mux.notfound())
serve(test, 8001)
@test String(HTTP.get("http://localhost:8001").body) ==
            "
Hello World!
"
@test String(HTTP.get("http://localhost:8001/about").body) ==
            "
Boo!
"
@test String(HTTP.get("http://localhost:8001/user/julia").body) ==
            "
Hello, julia!
"
@test String(HTTP.get("http://localhost:8001/badurl";
                      status_exception=false).body) ==
             "Internal server error"

ext: jl


text:
using DataStructures

### List

# TODO. Find where rules for List are put this doc there.
@sjdoc List """
    [e1, e2, ...],  List(e1, e2, ...)

creates a list. Lists are printed delimited by square brackets. The head
of the internal form is `List`.
"""

### Args

@sjdoc Args """
    Args(ex)

replace the `Head` of expression `ex` with `List`.

`Args(ex)` is equivalent to `Apply(List,ex)`.

    Args(pyobj::PyCall.PyObject)

return the arguments of the sympy expression `pyobj`.
"""
@mkapprule Args :nargs => 1
@doap Args(ex::Mxpr) = MList(copy(margs(ex)))

### First

@mkapprule First :nargs => 1
@doap First(x::Mxpr) = isempty(x) ? mx : x[1]
@doap First(x) = mx

### Join

@mkapprule Join

@sjdoc Join """
    Join(expr1, expr2, ...)

concatenate arguments of expressions with the same `Head`, returning an expression with the same `Head`.
"""

@doap function Join(args::Mxpr{T}...) where T
    nargs = newargs()
    foreach( x -> append!(nargs,margs(x)), args)
    mxpra(mhead(args[1]),nargs)
end

### Rest

@mkapprule Rest :nargs => 1

@doap function Rest(x::Mxpr)
    isempty(x) && return mx
    nargs = newargs(length(x)-1)
    for i in 1:length(nargs)
        nargs[i] = deepcopy(x[i+1])
    end
    mxpra(mhead(x),nargs)
end

@doap Rest(x) = mx

### Most

@mkapprule Most :nargs => 1

@doap function Most(x::Mxpr)
    isempty(x) && return mx
    nargs = newargs(length(x)-1)
    for i in 1:length(nargs)
        nargs[i] = deepcopy(x[i])
    end
    mxpra(mhead(x),nargs)
end

@doap Most(x) = mx

### Last

@mkapprule Last :nargs => 1:2

@doap function Last(x::Mxpr)
    isempty(x) && return mx  # FIXME add warning
    x[end]
end

@doap function Last(x::Mxpr,default)
    isempty(x) && return doeval(default)  # Last has Attribute HoldRest
    x[end]
end

### Fold

@mkapprule Fold :nargs => 2:3

@sjdoc Fold """
    Fold(f,x,[a,b,c,...])

return  `f(f(f(x,a),b),c)...`

    Fold(f,lst)

return  `f(First(lst),Rest(lst))`.

`f` may be a `Symbol`, or a function, or Julia function. Pure functions are not yet implemented.
"""

@mkapprule FoldList :nargs => 2:3

for head in (:Fold, :FoldList)
    fl = (head == :FoldList)
@eval begin
    @doap function ($head)(f, x, lst::Mxpr)
                local res
                n = length(lst)
                n == 0 && return x
                $(fl ? :(nargs = newargs(n+1)) : nothing)
                $(fl ? :(nargs[1] = x) : nothing)
                if isa(f,Function)
                    res = doeval(Base.invokelatest(f,x,lst[1]))
                    $(fl ? :(nargs[2] = res) : nothing)
                    for i in 2:n
                        res = invokelatest(f,res,lst[i])
                        $(fl ? :(nargs[i+1] = res) : nothing)
                    end
                else
                    res = doeval(mxpr(f,x,lst[1]))
                    $(fl ? :(nargs[2] = res) : nothing)
                    for i in 2:n
                        res = doeval(mxpr(f,res,lst[i]))
                        $(fl ? :(nargs[i+1] = res) : nothing)
                    end
                end
                $(fl ? :(mxpra(mhead(lst), nargs)) : :( res ))
         end

   @doap function ($head)(f,lst::Mxpr)
           isempty(lst) && return mx
           mxpr($(QuoteNode(head)), f, lst[1], mxpr(mhead(lst),lst[2:end]))
         end
end
end

@mkapprule Nest :nargs => 3
@mkapprule NestList :nargs => 3

## FIXME: Using invokelatest makes a typical call 10x slower.
## But, it seems we have to do this in v0.6

for head in (:Nest, :NestList)
    nl = (head == :NestList)
@eval begin
    @doap function ($head)(f, x, n::Integer)
                local res
                n == 0 && return $(nl ? :(mxpr(:List,x)) : :(x))
                $(nl ? :(nargs = newargs(n+1)) : nothing)
                $(nl ? :(nargs[1] = x) : nothing)
                if isa(f,Function)
                    res = doeval(Base.invokelatest(f,x))
                    if is_throw() return res end
                    $(nl ? :(nargs[2] = res) : nothing)
                    for i in 2:n
                        res = invokelatest(f,res)
                        if is_throw() return res end
                        $(nl ? :(nargs[i+1] = res) : nothing)
                    end
                else
                    res = doeval(mxpr(f,x))
                    if is_throw() return res end
                    $(nl ? :(nargs[2] = res) : nothing)
                    for i in 2:n
                        res = doeval(mxpr(f,res))
                        if is_throw()
                            return res
                        end
                        $(nl ? :(nargs[i+1] = res) : nothing)
                    end
                end
                $(nl ? :(mxpra(:List, nargs)) : :( res ))
         end
end
end

### Range

@sjdoc Range """
    Range(n)

return the `List` of integers from `1` through `n`.

    Range(n1,n2)

return the `List` of numbers from `n1` through `n2`.

    Range(n1,n2,di)

return the `List` of numbers from `n1` through `n2` in steps of `di`.
`di` may be negative. Floats and some symbolic arguments are supported.

Similar `Lists` can be created by importing Julia `Array`s with `Unpack(:([1.0:10^5]))`.
This uses embedded Julia to create a typed `Array` and then unpacks it to a List.
"""

# Need to check for uprules for free symbols
function apprules(mx::Mxpr{:Range})
    iter = make_sjitera(margs(mx))
    args = do_Range(iter)
    r = mxpra(:List,args)
    setfixed(r)
    setcanon(r)
    mergesyms(r,:nothing)  # not correct if we have symbols.
    return r
end

function do_Range(iter::SJIterA1)  # iter is parameterized, so we hope type of n is inferred.
    n = iter.num_iters
    args = newargs(n);
    _do_Range_fill(args,n,typeof(iter.imax))
    return args
end

function _do_Range_fill(args, n, ::Type{Int})
    copyto!(args,1:n)
  # @inbounds @simd for i in 1:n
  #                    args[i] = i
  #                 end
    return args
end

function _do_Range_fill(args, n, ::Type{T}) where T<:Real
#    j = one(T)
    n1 = one(T)
    n2 = convert(T,n)
    copyto!(args,n1:n2)
    # @inbounds for i in 1:n
    #    args[i] = j
    #    j += 1
    # end
    return args
end

# Fails for rationals. and counting is wrong
function do_Range(iter::SJIterA2{T,V}) where {T<:Real,V<:Real}
    nd = round(Int,mplus(iter.imax,-iter.imin) + 1)  # Try bug fix!
    if nd > 1
        args = newargs(iter.num_iters)
        @inbounds for i in 0:nd-1
            args[i+1] = mplus(i,iter.imin)  # Bug here, if iter.imin is not of type Int.
        end
    else  # Mma does not allow this second branch: eg Range(5,1) implies di = -1
        nd = -nd + 2
        args = newargs(iter.num_iters)
        @inbounds for i in zero(iter.imin):(nd - 1)
            args[i+1] = mplus(iter.imin, -i)
        end
    end
    return args
end

# Symbolic values
# FIXME. We don't record free symbols and check for upvalues.
# This is about as fast as Mma 3 (running on a somewhat slower cpu)
function do_Range(iter::SJIterA2)
    args = newargs(iter.num_iters)
    imin = iter.imin
    args[1] = imin
    s = imin
    if is_Mxpr(imin,:Plus)  # imin is a sum
        if is_Number(imin[1])  # number is always first in canon order.
            b = imin[1]  # extract number
            r = imin[2:end]  # the rest of the sum
            _do_Range_A2(args,b,r,iter.num_iters)
        else  # imin is a sum with no numbers, so we put a number in front
            sargs = margs(s)
            for i in 2:iter.num_iters
                args[i] = mxpr(:Plus,i-1,sargs...)
                setfixed(args[i])
            end
        end
    else  # imin is not a sum
        @inbounds for i in 2:iter.num_iters
            args[i] = mxpr(:Plus,i-1,s)
            setfixed(args[i])
        end
    end
    #  we don't handle counting down case.
    return args
end

function _do_Range_A2(args,b,r,n)
    @inbounds for i in 2:n
        args[i] = mxpr(:Plus,b+i-1,r...) # only a little slower than Mma if disable gc
        setfixed(args[i])
    end
    return args
end

# seems to be little penalty for mplus instead of +
function do_Range(iter::SJIterA3{T,V,W}) where {T<:Real,V<:Real,W<:Real}
    n = iter.num_iters
    args = newargs(n)
    j = iter.imin
    @inbounds for i in 1:n
        args[i] = j
        j = mplus(j,iter.di)
    end
    return args
end

# Symbolic again
function do_Range(iter::SJIterA3)
    args = newargs(iter.num_iters)
    imin = iter.imin
    args[1] = imin
    s = imin
    if is_Number(iter.di)
        if true
            if is_Mxpr(imin,:Plus)
                if is_Number(imin[1])  # number is always first in canon order.
                    b = imin[1]        # extract number
                    r = imin[2:end]    # the rest of the sum
                    for i in 2:iter.num_iters
                        b = b + iter.di
                        if b == 0  # more efficient to move this branch out
                            if length(r) == 1
                                args[i] = r[1]
                            else
                                args[i] = mxpr(:Plus,r...)
                            end
                        else
                            args[i] = mxpr(:Plus,b,r...)
                        end
                        setfixed(args[i])
                    end
                else  # imin is a sum with no numbers, so we put a number in front
                    sargs = margs(s)
                    j = zero(iter.di)
                    for i in 2:iter.num_iters
                        j += iter.di
                        args[i] = mxpr(:Plus,j,sargs...)
                        setfixed(args[i])
                    end
                end
            else # imin is not a sum, so just create one
                j = zero(iter.di)
                @inbounds for i in 2:iter.num_iters
                    j += iter.di
                    args[i] = mxpr(:Plus,j,s)
                    setfixed(args[i])
                end
            end
        else  #  iter.di < 0
            error("unimplemented")
        end
    else # di is not a number
        error("unimplemented")
    end
    return args
end

# Some is implemented here that is not in the new Range yet
function apprules(mx::Mxpr{:OldRange})
    if length(mx) == 1
        n = mx[1]
        args = range_args1(n) # use function for optimization on type
    elseif length(mx) == 2
        n0 = mx[1] - 1
        n = mx[2]
        args = range_args2(n0,n)
    elseif length(mx) == 3
        n0 = mx[1]
        n = mx[2]
        di = mx[3]
        off = n > n0 ? 1 : -1
        args = range_args3(n0,n,di,off)
    else
        return mx
    end
    r = mxpra(:List,args)
    setfixed(r)
    setcanon(r)
    mergesyms(r,:nothing)
    return r
end

# separate functions are *essential* for type stability and efficiency.
function range_args1(n::T) where T<:Integer
    args = newargs(n);
    @inbounds for i in one(n):n
        args[i] = i
    end
    return args
end

function range_args1(n::T) where T<:AbstractFloat
    ni = floor(Int,n)
    args = newargs(ni);
    @inbounds for i in 1:ni
        args[i] = convert(T,i)
    end
    return args
end

function range_args2(n0,n)
    nd = n - n0
    args = newargs(nd);
    @inbounds for i in one(n0):nd
        args[i] = i+n0
    end
    return args
end

function range_args3(n0,n,di,off)
    args = newargs(div(n-n0+off,di));
    len = length(args) # cheap
    s = n0
    @inbounds for i in one(n0):len
        args[i] = s
        s += di
    end
    args
end

### ConstantArray

@sjdoc ConstantArray """
    ConstantArray(expr,n)

create a `List` of `n` copies of `expr`.
"""

# We take only attribute to be Protected. So expr is evaled already

@mkapprule ConstantArray

# function apprules(mx::Mxpr{:ConstantArray})
#     do_ConstantArray(mx,margs(mx)...)
# end
#do_ConstantArray(mx,args...) = mx

# The annotation for Number is needed, because deepcopy tries
# to do something very slow with numbers.
# Copying a small Mxpr is extremely slow
# 'c^2' is 200 times slower than Symbol 'c'.

@doap function ConstantArray(expr,n::Integer)
    _constantarray(expr,n)
end

function _constantarray(expr,n::Integer)
    nargs = newargs(n)
    @inbounds for i in 1:n
        nargs[i] = recursive_copy(expr)
    end
    mxpra(:List,nargs)
end

function _constantarray(expr::Mxpr,n::Integer)
    nargs = newargs(n)
     @inbounds for i in 1:n
        nargs[i] = setfixed(recursive_copy(expr))
     end
    mxpra(:List,nargs)
end

@doap function ConstantArray(expr,ns::Mxpr{:List})
    all(x->isa(x,Integer), margs(ns)) || return mx
    _constantarray(expr,reverse(margs(ns)))
end

function _constantarray(expr,ns::Array)
    isempty(ns) && return Null
    a = _constantarray(expr,ns[1])
    length(ns) == 1 && return a
    _constantarray(a,ns[2:end])
end

function _constantarray(expr::T,n::Integer) where T<:Union{Number,SJSym,String}
    nargs = newargs(n)
    fill!(nargs,expr)
    a = mxpra(:List,nargs)
    expr == :Nothing && return a
    setfixed(a)
end

### Nothing

# Nothing is removed where it appears in a List, but not from expressions with other Heads.
# This is too expensive to implement as a rule. Better to write special code in evaluation.jl

@mkapprule Nothing :nodefault => true

@sjdoc Nothing """
     Nothing(args...)

is replaced by `Nothing`.

Instances of `Nothing` as arguments to `List` are removed.

```
[a,b,Nothing,c] == [a,b,c].
```

`Nothing` is only removed from the arguments of expressions with `Head` `List`.
"""

@doap Nothing(args...) = :Nothing

### Keys

@sjdoc Keys """
    Keys(d)

returns a list of the keys in `Dict` `d`.
"""

apprules(mx::Mxpr{:Keys}) = do_keys(mx,mx[1])
do_keys(mx, d::AbstractDict) = mxpr(:List,collect(Any,keys(d))...)
do_keys(mx, x) = (@symwarn("Can't return keys of $x"); mx)

### Values

@sjdoc Values """
    Values(d)

return a list of the values in `Dict` `d`.
"""

apprules(mx::Mxpr{:Values}) = do_values(mx,mx[1])
do_values(mx, d::AbstractDict) = mxpr(:List,collect(Any,values(d))...)
do_values(mx, x) = (@symwarn("Can't return values of $mx"); mx)

### Splat

@mkapprule Splat :nargs => 1

@sjdoc Splat """
    Splat(expr)

equivalent to `Apply(Sequence, expr)`.

In a list of arguments `Splat(expr)` is replaced by the arguments of `expr`.

```
symata> f(a,b, Splat([c,d]))
 f(a,b, c, d)
```
"""

@doap Splat(x::Mxpr) = mxprcf(:Sequence, margs(x))
@doap Splat(x) = x

### Sort

@mkapprule Sort

@sjdoc Sort """
    Sort(expr)

sort the elements of `expr`.
"""

@doap Sort(expr::Mxpr{:List}) = mxpr(:List,sort(margs(expr)))
@doap Sort(a::AbstractArray) = sort(a)

### Split

@sjdoc Split """
    Split(lst)

splits `lst` into runs of identical elements.

    Split(lst,test)

apply `test` to determine identical elements.
"""
@mkapprule Split

@doap function Split(lst::ListT)
    isempty(lst) && return MList()
    a = margs(lst)
    a0 = newargs()
    a1 = newargs()
    for i in 1:length(a)
        if length(a1) != 0 && a[i-1] != a[i]
            push!(a0, MList(a1))
            a1 = newargs()
        end
        push!(a1, a[i])
    end
    push!(a0,MList(a1))
    MList(a0)
end

@doap function Split(lst::ListT, test)
    isempty(lst) && return MList()
    a = margs(lst)
    a0 = newargs()
    a1 = newargs()
    for i in 1:length(a)
        if length(a1) != 0
            t = isa(doeval(mxpr(test, a[i-1], a[i])),Bool)
            if ! t
                push!(a0, MList(a1))
                a1 = newargs()
            end
        end
        push!(a1, a[i])
    end
    push!(a0,MList(a1))
    MList(a0)
end

### Partition

@sjdoc Partition """
    Partition(list,n)

partitions `list` into lists of `n` elements.

    Partition(list,n,m)

partitions `list` into lists of `n` elements with offset `m`.

`Partition(list,n)` is equivalent to `Partition(list,n,n)`
"""
@mkapprule Partition :nargs => 2:3

@doap Partition(x::Mxpr, n::Integer) = partition(x,n,n)
@doap Partition(x::Mxpr, n::Integer, m::Integer) = partition(x,n,m)

function partition(x::Mxpr, n, m)
    m < 1 && symerror("Partition, m must be greater than zero")
    h = mhead(x)
    a = margs(x)
    na = newargs()
    i = 1
    while true
        na1 = newargs(n)
        if i + n - 1 > length(a) break end
        for j in 1:n
            na1[j] = a[i]
            i += 1
        end
        i -= (n - m)
        push!(na, mxpr(h,na1))
    end
    mxpra(h,na)
end

### Riffle

@mkapprule Riffle

## TODO: implement all features. document
@doap function Riffle(x::Mxpr{:List},ex)
    n = 2 * length(x)
    nargs = newargs(n)
    j = 0
    for i in 1:2:n
        j += 1
        nargs[i] = x[j]
        nargs[i+1] = ex
    end
    MListA(nargs)
end

@doap function Riffle(x::Mxpr{:List},ex::Mxpr{:List})
    nargs = newargs()
    j = 1
    n = length(ex)
    for i in 1:length(x)
        push!(nargs,x[i])
        push!(nargs,ex[j])
        j = j >= n ? 1 : j + 1
    end
    MListA(nargs)
end

### Union

@mkapprule Union
@doap function Union(exprs::Mxpr...)
    isempty(exprs) && return mxpr(:List)
    return mxpra(mhead(exprs[1]), union(map(margs, exprs)...))
end

@mkapprule Intersection
@doap function Intersection(exprs::Mxpr...)
    isempty(exprs) && return mxpr(:List)
    return mxpra(mhead(exprs[1]), intersect(map(margs, exprs)...))
end

@mkapprule Complement
@doap function Complement(exprs::Mxpr...)
    isempty(exprs) && return mxpr(:List)
    return mxpra(mhead(exprs[1]), setdiff(map(margs, exprs)...))
end

### Thread

@sjdoc Thread """
    Thread(f(args))

threads `f` over any lists appearing in `args`.

    Thread(f(args),h)

threads `f` over any expressions with head `h`

    Thread(f(args),h,spec)

threads only over elments specified by the standard seuence specification `spec.
"""
@mkapprule Thread :nargs => 1:3

@doap Thread(x::Mxpr) = threadlistable(x)
@doap Thread(x::Mxpr,head) = threadlistable(x,head)

## TODO: generate this with a macro or something, there is a copy in evaluation.jl with :List hardcoded
function threadlistable(mx::Mxpr, head)
    pos = Array{Int}(undef, 0)      # should avoid this
    lenmx = length(mx)
    lenlist::Int = -1
    h = mhead(mx)
    @inbounds for i in 1:lenmx
        if is_Mxpr(mx[i],head)
            nlen = length(mx[i])
            if lenlist >= 0 && nlen != lenlist
                error("Can't thread over lists of different lengths.")
            end
            lenlist = nlen
            push!(pos,i)
        end
    end
    lenp = length(pos)
    lenp == 0 && return mx      # Nothing to do. return input array
    largs = newargs(lenlist)
    @inbounds for i in 1:lenlist
        nargs = newargs(lenmx)
        p = 1
        @inbounds for j in 1:lenmx
            if p <= lenp && pos[p] == j
                nargs[j] = mx[j][i]
                p += 1
            else
                nargs[j] = mx[j]
            end
        end
        largs[i] = mxpr(h,nargs)
    end
    nmx = mxpr(head,largs)
    return nmx
end

@doap function Thread(x::Mxpr,head,inspec)
    seqspec = sequencespec(inspec,length(x))
    threadlistable(x,head,seqspec)
end

threadlistable(x,head,s::SequenceAll) = threadlistable(x,head)
threadlistable(x,head,s::SequenceNone) = x  # copy ?

function threadlistable(mx, head, seqspec::SequenceSpec)
    pos = Array{Int}(undef, 0)      # should avoid this
    lenmx = length(mx)
    lenlist::Int = -1  # length of lists that we will thread over
    h = mhead(mx)
    (n1,n2,di) = seqiter(seqspec,lenmx)
    @inbounds for i in n1:di:n2
        if is_Mxpr(mx[i],head)
            nlen = length(mx[i])
            if lenlist >= 0 && nlen != lenlist
                error("Can't thread over lists of different lengths.")
            end
            lenlist = nlen
            push!(pos,i)
        end
    end
    lenp = length(pos)
    lenp == 0 && return mx      # Nothing to do. return input array
    largs = newargs(lenlist)
    @inbounds for i in 1:lenlist
        nargs = newargs(lenmx)
        p = 1
        @inbounds for j in 1:lenmx
            if p <= lenp && pos[p] == j
                nargs[j] = mx[j][i]
                p += 1
            else
                nargs[j] = mx[j]
            end
        end
        largs[i] = mxpr(h,nargs)
    end
    nmx = mxpr(head,largs)
    return nmx
end

### MapThread

@mkapprule MapThread nargs => 2  """
    MapThread(f,list)

threads `f` over lists in `list` and evaluates the result.
"""

@doap MapThread(f,x) = threadlistable(mxpr(f,margs(x)))

### Counts

@mkapprule Counts nargs => 1 """
    Counts(list)

return a dictionary of the number of times each distinct element of list
occurs.
"""

@doap function Counts(list::ListT)
    d = OrderedDict{Any,Int}()
    for el in list
        val = get!(d,el,0)
        d[el] = val + 1
    end
    return d
end

### Select

@mkapprule Select nargs => 1:3

@doap function Select(x::Mxpr, crit)
    nargs = newargs()
    for y in x
        if trueq(doeval(mxpr(crit,y))) push!(nargs, y) end
    end
    mxpra(mhead(x),nargs)
end

@doap function Select(x::Mxpr, crit, n::Integer)
    nargs = newargs()
    c = 1
    for y in x
        if trueq(doeval(mxpr(crit,y)))
            push!(nargs, y)
            c += 1
            c > n && break
        end
    end
    mxpra(mhead(x),nargs)
end

@curry_second Select


text:
include("rediscommon.jl")

basedir = "/home/kfischer/cs262project/redis/utils/create-cluster/traces"

reads = Dict{Int,Vector{Any}}()
writes = Dict{Int,Vector{Any}}()
accepts = Dict{Int,Vector{Any}}()
closes = Dict{Int,Vector{Any}}()

timeline, modules = replay(joinpath(basedir,"redis-cli-0"));
client_id = 0
client_reads, client_writes, client_accepts, client_closes, client_nticks = trace_rw(timeline, modules);
reads[0] = client_reads
writes[0] = client_writes
client_ends = [start_time(timeline); end_time(timeline)]
client_layers = timeline_layers([
client_ends;
map(x->TimelineEvent(x[2],:read), client_reads);
map(x->TimelineEvent(x[2],:write), client_writes);
]; timeline_id = client_id);

server_layers = Any[]
nticks = Dict{Any,Int}()
for i = 0:5
  server_id = i+1
  timeline, modules = replay(joinpath(basedir,"redis-server-$i"))
  server_reads, server_writes, server_accepts, server_closes, server_nticks = trace_rw(timeline, modules);
  reads[server_id] = server_reads
  writes[server_id] = server_writes
  accepts[server_id] = server_accepts
  closes[server_id] = server_closes
  nticks[server_id] = server_nticks
  server_ends = [start_time(timeline); end_time(timeline);]
  append!(server_layers,timeline_layers([
  server_ends;
  map(x->TimelineEvent(x[2],:read), server_reads);
  map(x->TimelineEvent(x[2],:write), server_writes);
  ]; timeline_id = server_id));
end

messages = Dict{Any,Any}()
for i=0:6, j=0:6
    @show (i,j)
    messages[i=>j] = correlate_messages(writes[i], reads[j], i, j)
end

ranges = Dict{Int,Any}()
for i = 1:6
    ranges[i] = compute_fd_ranges(accepts[i], closes[i]);
end

# Filter messages by range
for i = 1:6
    if isempty(messages[0=>i])
        messages[i=>0] = Message[]
        continue
    end
    m = messages[0=>i][1]
    server_fd = m.recv_id
    range = compute_range(ranges[i], m.send_t, server_fd)
    messages[i=>0] = collect(filter(x->((x.send_fd == server_fd) && x.send_t in range), messages[i=>0]))
end

all_messages = reduce(vcat,map(collect,values(messages)))

m_colors = [colorant"red",colorant"blue",colorant"green",colorant"orange",colorant"purple",colorant"black",colorant"magenta"]

m_offset = messages[0=>6][1].recv_t

client_layers = timeline_layers([
client_ends;
map(x->TimelineEvent(x[2]+m_offset,:read), client_reads);
map(x->TimelineEvent(x[2]+m_offset,:write), client_writes);
]; timeline_id = client_id);

layers = [layer(x=[m.send_t + ((m.send_id == 0) ? m_offset : 0), m.recv_t + ((m.recv_id == 0) ? m_offset : 0)],
    y=[m.send_id, m.recv_id], Geom.line, Theme(default_color=m_colors[m.send_id+1],line_width=0.5mm)) for m in all_messages];

preplication(theme) = plot(client_layers..., server_layers..., layers...,
    theme,Guide.xlabel("Time"),
    Guide.ylabel("Timeline"),
    Coord.Cartesian(xmin = m_offset),
    Scale.color_discrete_manual(colors...), Scale.y_continuous(minvalue=-1, maxvalue=7))

display(preplication(repl_theme))
draw(PDF("preplication.pdf", 17Compose.cm, 2inch), preplication(paper_theme))


text:
# Draws histogram for a given one dimesional data
function hist2D(data,binsizex)
    xmax = maximum(data[:,1])
    xmin = minimum(data[:,1])
    xmax = xmax+0.15*(xmax-xmin)
    xmin = xmin-0.15*(xmax-xmin)

    # X grid point locations
    xpoints = linspace(xmin, xmax, binsizex+1)
    # Bins
    bins = zeros(binsizex,1)
    # Bin allocation
    for i = 1:size(data)[1]
        tmp1 = data[i,1]*ones(binsizex,1) - xpoints[1:(end-1),1]
        [tmp1[i]=(tmp1[i] >= 0 ? 1 : 0) for i in 1:binsizex]
        xbinindx = Int(sum(tmp1))
        # print(tmp1)
        # if xbinindx==0
        #     error("Zero in column 1 not possible")
        # end
        bins[xbinindx,1] = bins[xbinindx,1] + 1
    end
    return bins , xpoints
end

max_stars_repo_name: niladridas/UKF

max_issues_repo_name: niladridas/UKF

max_forks_repo_name: niladridas/UKF


text:
Base.@deprecate MockModel() mock_model()

function mock_model()
    optimizer = MOIU.MockOptimizer(Parametron.ParametronMOIModel{Float64}())
    Parametron.Model(optimizer)
end


text:
function simulate!(dataset::Dataset, p)
    simulate!.(dataset.ponctual,Ref(p))
    simulate!.(dataset.timeseries, Ref(p))
end

##################################
### SIMULATE PONCTUAL DATASETS ###
##################################

function simulate!(pds::PonctualDataset, p)
    data = pds.data
    for i in 1:nrow(data)
        data_params = data[i,Not(pds.target)] # Dataframerow without target variable
        df_sol = p.f.ponctual(p, data_params, pds.target)
        (i == 1) && insert_sim_cols!(pds, names(df_sol))
        data[i,sol_colnames] .= df_sol[1,!]
    end
    return nothing
end

####################################
### SIMULATE TIMESERIES DATASETS ###
####################################

function simulate!(tsds::TimeseriesDataset, p)
    data = tsds.data
    gdata = groupby(data,tsds.var)
    for i in 1:length(gdata)
        data_params = gdata[i][1,Not([tsds.target ; :t])] # Dataframerow without target variable
        df_sol = p.f.timeseries(p, data_params, tsds.target, gdata[i].t) #
        sol_colnames = names(df_sol)
        (i == 1) && insert_sim_cols!(tsds, sol_colnames)
        gdata[i][1:size(df_sol,1),sol_colnames] .= df_sol[:,sol_colnames]
    end
    return nothing
end

### HELPER ###
function insert_sim_cols!(ds::DatasetType, sim_colnames)
    colnames = names(ds.data)
    for i in 1:length(sim_colnames)
        if string(sim_colnames[i]) ∉ colnames
            insertcols!(ds.data, last(names(ds.data)),
                sim_colnames[i] => Vector{Union{Missing,Float64}}(missing,size(ds.data,1)),
                after=true)
        end
    end
end

export simulate!

max_stars_repo_name: Leooop/DataFormatter

max_issues_repo_name: Leooop/DataFormatter

max_forks_repo_name: Leooop/DataFormatter


text:
## Transformations
probit(x) = quantile(Normal(), x)

function transform!(trans_samples::AbstractMatrix, samples::AbstractMatrix, lb::AbstractVector, ub::AbstractVector)
    nd, ns = size(samples)
    for j = 1:ns, i = 1:nd
        trans_samples[i,j] = transform(samples[i,j], lb[i], ub[i])
    end
    return trans_samples
end
transform(samples::AbstractMatrix, lb::AbstractVector, ub::AbstractVector) = transform!(similar(samples), samples, lb, ub)
transform(x, lb::Real, ub::Real) = probit((x - lb)/(ub - lb))
transform(x, ::Missing, ::Missing) = x
transform(x, lb::Real, ::Missing) = log(x - lb)
transform(x, ::Missing, ub::Real) = log(ub - x)

# With bijectors
transform(x, dist::Distribution, trans::Bool) = trans ? Bijectors.link(dist,x) : x

# Inverse transformations
function invtransform!(samples::AbstractMatrix, trans_samples::AbstractMatrix, lb::AbstractVector, ub::AbstractVector)
    nd, ns = size(samples)
    for j = 1:ns, i = 1:nd
        samples[i,j] = invtransform(trans_samples[i,j], lb[i], ub[i])
    end
    return samples
end
invtransform(trans_samples::AbstractMatrix, lb::AbstractVector, ub::AbstractVector) = invtransform!(similar(trans_samples), trans_samples, lb, ub)
invtransform(x, lb::Real, ub::Real) = lb + (ub - lb) * cdf(Normal(), x)
invtransform(x, ::Missing, ::Missing) = x
invtransform(x, lb::Real, ::Missing) = exp(x) + lb
invtransform(x, ::Missing, ub::Real) = ub - exp(x)

# With bijectors
invtransform(x, dist::Distribution, trans::Bool) = trans ? Bijectors.invlink(dist,x) : x

# Log jacobian contribution : x is the transformed variable
logjacobian_contribution(x, lb, ub) = log(ub-lb) + logpdf(Normal(), x)
logjacobian_contribution(x, ::Missing, ::Missing) = 0.0
logjacobian_contribution(x, lb, ::Missing) = x
logjacobian_contribution(x, ::Missing, ub) = x

# With bijectors
function logjacobian_contribution(x, dist::Distribution, trans::Bool)
    b =  Bijectors.bijector(dist)
    return trans ? Bijectors.logabsdetjac(b, Bijectors.invlink(dist, x)) : 0.0
end


text:
using Serialization

function julia_func(z)
    c = complex(-0.5, 0.61)
    z = z^2 + c
    return z
end

function gen_julia()
    im_width, im_height = 5000, 5000
    max_it = 300
    max_z2 = 100
    min_x, max_x = -1, 1
    min_y, max_y = -1, 1

    julia_set = zeros(im_width,im_height)
    for ix in 1:im_width
        for iy in 1:im_height
            real_part = ((max_x-min_x)/im_width)*ix + min_x
            im_part = ((max_y-min_y)/im_height)*iy + min_y
            z = complex(real_part, im_part)
            it = 0
            while(abs2(z) < max_z2 && it < max_it)
                z = julia_func(z)
                it += 1
            end
            ratio = it/max_it
            julia_set[ix,iy] = ratio
        end
    end
    return julia_set
end

function main()
    @time julia_set = gen_julia()
    serialize("julia_data",julia_set)
end

main()

max_stars_repo_name: ofbennett/julia-set-fractal

max_issues_repo_name: ofbennett/julia-set-fractal

max_forks_repo_name: ofbennett/julia-set-fractal


text:
export  Protocol, ProtocolParams, name, description, scanner, params, runProtocol,
        init, execute, cleanup, timeEstimate, ProtocolEvent, InfoQueryEvent,
        InfoEvent, DecisionEvent, AnswerEvent, StopEvent, ResumeEvent, CancelEvent, RestartEvent, ProgressQueryEvent,
        ProgressEvent, UndefinedEvent, DataQueryEvent, DataAnswerEvent, FinishedNotificationEvent, FinishedAckEvent,
        ExceptionEvent, IllegaleStateEvent, DatasetStoreStorageRequestEvent, StorageSuccessEvent, StorageRequestEvent,
        OperationSuccessfulEvent, OperationUnsuccessfulEvent, OperationNotSupportedEvent, MultipleChoiceEvent, ChoiceAnswerEvent

abstract type ProtocolParams end

export ProtocolState, PS_UNDEFINED, PS_INIT, PS_RUNNING, PS_PAUSED, PS_FINISHED, PS_FAILED
@enum ProtocolState PS_UNDEFINED PS_INIT PS_RUNNING PS_PAUSED PS_FINISHED PS_FAILED

name(protocol::Protocol)::AbstractString = protocol.name
description(protocol::Protocol)::AbstractString = protocol.description
scanner(protocol::Protocol)::MPIScanner = protocol.scanner
params(protocol::Protocol)::ProtocolParams = protocol.params
biChannel(protocol::Protocol) = protocol.biChannel

"General constructor for all concrete subtypes of Protocol."
function Protocol(protocolDict::Dict{String, Any}, scanner::MPIScanner)
  if haskey(protocolDict, "name")
    name = pop!(protocolDict, "name")
  else
    throw(ProtocolConfigurationError("There is no protocol name given in the configuration."))
  end

  if haskey(protocolDict, "description")
    description = pop!(protocolDict, "description")
  else
    throw(ProtocolConfigurationError("There is no protocol description given in the configuration."))
  end

  if haskey(protocolDict, "targetScanner")
    targetScanner = pop!(protocolDict, "targetScanner")
    if targetScanner != scannerName(scanner)
      throw(ProtocolConfigurationError("The target scanner (`$targetScanner`) for the protocol does not match the given scanner (`$(scannerName(scanner))`)."))
    end
  else
    throw(ProtocolConfigurationError("There is no target scanner for the protocol given in the configuration."))
  end

  if haskey(protocolDict, "type")
    protocolType = pop!(protocolDict, "type")
  else
    throw(ProtocolConfigurationError("There is no protocol type given in the configuration."))
  end

  paramsType = getConcreteType(ProtocolParams, protocolType*"Params")
  params = paramsType(protocolDict, scanner)
  ProtocolImpl = getConcreteType(Protocol, protocolType)

  return ProtocolImpl(name=name, description=description, scanner=scanner, params=params)
end

function Protocol(protocolName::AbstractString, scanner::MPIScanner)
  configDir_ = configDir(scanner)
  filename = joinpath(configDir_, "Protocols", "$protocolName.toml")

  if isfile(filename)
    protocolDict = TOML.parsefile(filename)
  else
    throw(ProtocolConfigurationError("Could not find a valid configuration for protocol with name `$protocolName` and the derived path `$filename`."))
  end

  protocolDict["name"] = protocolName

  return Protocol(protocolDict, scanner)
end

Protocol(protocolName::AbstractString, scannerName::AbstractString) = Protocol(protocolName, MPIScanner(scannerName))
Protocol(protocolDict::Dict{String, Any}, scannerName::AbstractString) = Protocol(protocolDict, MPIScanner(scannerName))

function runProtocol(protocol::Protocol)
  # TODO: Error handling
  # TODO command line "handler"
  channel = init(protocol)
  @async begin
    execute(protocol)
    cleanup(protocol)
  end
  return channel
end

abstract type ProtocolEvent end

@mustimplement _init(protocol::Protocol) # Prepare necessary data structures, ex. buffer for samples + background meas
@mustimplement _execute(protocol::Protocol)
@mustimplement enterExecute(protocol::Protocol) # Prepare execution tracking flags, ex. execution done or cancelled flags
@mustimplement cleanup(protocol::Protocol)
@mustimplement stop(protocol::Protocol)
@mustimplement resume(protocol::Protocol)
@mustimplement cancel(protocol::Protocol)

function init(protocol::Protocol)
  checkRequiredDevices(protocol)
  _init(protocol)
end

function checkRequiredDevices(protocol::Protocol)
  missingDevices = []
  scanner = protocol.scanner
  for device in requiredDevices(protocol)
    if length(getDevices(scanner, device)) == 0
      push!(missingDevices, device)
    end
  end
  if length(missingDevices) > 0
    message = "Protocol $(typeof(protocol)) is missing the following required devices: " * join(string.(missingDevices), ", ", " and ")
    throw(IllegalStateException(message))
  end
end

requiredDevices(protocol::Protocol) = []

timeEstimate(protocol::Protocol) = "Unknown"

function execute(protocol::Protocol, threadID::Integer = 1)
  if isnothing(protocol.executeTask) || istaskdone(protocol.executeTask)
    protocol.biChannel = BidirectionalChannel{ProtocolEvent}(32)
    @tspawnat threadID executionTask(protocol)
    return BidirectionalChannel{ProtocolEvent}(protocol.biChannel)
  else
    @warn "Protocol cannot be executed again as it is still running"
    return nothing
  end
end

function executionTask(protocol::Protocol)
  protocol.executeTask = current_task()
  try
    enterExecute(protocol)
    _execute(protocol)
  catch ex
    if ex isa CancelException
      put!(protocol.biChannel, OperationSuccessfulEvent(CancelEvent()))
      close(protocol.biChannel)
    elseif ex isa IllegalStateException
      put!(protocol.biChannel, IllegaleStateEvent(ex.message))
      close(protocol.biChannel)
    else
      # Let task fail
      put!(protocol.biChannel, ExceptionEvent(ex))
      close(protocol.biChannel)
      throw(ex)
    end
  end
  close(protocol.biChannel)
end

struct UndefinedEvent <: ProtocolEvent
  event::ProtocolEvent
end
# Interaction Events, only necessary for interactive protocolts
struct DecisionEvent <: ProtocolEvent
  message::AbstractString
end
struct AnswerEvent <: ProtocolEvent
  answer::Bool
  question::DecisionEvent
end
struct MultipleChoiceEvent <: ProtocolEvent
  message::AbstractString
  choices::Vector{AbstractString}
end
struct ChoiceAnswerEvent <: ProtocolEvent
  answer::Int64
  question::MultipleChoiceEvent
end

# (Mandatory) Control flow events for all protocols
struct StopEvent <: ProtocolEvent end
struct ResumeEvent <: ProtocolEvent end
struct CancelEvent <: ProtocolEvent end
struct RestartEvent <: ProtocolEvent end
struct OperationNotSupportedEvent <: ProtocolEvent
  operation::ProtocolEvent
end
struct OperationSuccessfulEvent <: ProtocolEvent
  operation::ProtocolEvent
end
struct OperationUnsuccessfulEvent <: ProtocolEvent
  operation::ProtocolEvent
end
struct IllegaleStateEvent <: ProtocolEvent
  message::AbstractString
end
struct ExceptionEvent <: ProtocolEvent
  exception::Exception
end
struct FinishedNotificationEvent <: ProtocolEvent end
struct FinishedAckEvent <: ProtocolEvent end
#Maybe a status (+ query) event and all Protocols have the states: UNKNOWN, INIT, EXECUTING, PAUSED, FINISHED

# Display/Information Events
struct ProgressQueryEvent <: ProtocolEvent end
struct ProgressEvent <: ProtocolEvent
  done::Int
  total::Int
  unit::AbstractString
  query::ProgressQueryEvent
end
struct DataQueryEvent <: ProtocolEvent
  message::AbstractString
end
struct DataAnswerEvent <: ProtocolEvent
  data::Any
  query::DataQueryEvent
end
struct UnknownDataQueryEvent <: ProtocolEvent
  query::DataQueryEvent
end

abstract type StorageRequestEvent <: ProtocolEvent end
struct DatasetStoreStorageRequestEvent <: StorageRequestEvent
  datastore::DatasetStore
  params::Dict
end
struct StorageSuccessEvent <: ProtocolEvent
  filename::AbstractString
end

function askConfirmation(protocol::Protocol, message::AbstractString)
  channel = biChannel(protocol)
  question = DecisionEvent(message)
  put!(channel, question)
  while isopen(channel) || isready(channel)
    event = take!(channel)
    # Note that for immutable objects the '==' does not guarantee that the reply is to the actual current question
    if event isa AnswerEvent && event.question == question
      return event.answer
    else
      handleEvent(protocol, event)
    end
    sleep(0.001)
  end
end

function askChoices(protocol::Protocol, message::AbstractString, choices::Vector{<:AbstractString})
  channel = biChannel(protocol)
  question = MultipleChoiceEvent(message, choices)
  put!(channel, question)
  while isopen(channel) || isready(channel)
    event = take!(channel)
    # Note that for immutable objects the '==' does not guarantee that the reply is to the actual current question
    if event isa ChoiceAnswerEvent && event.question == question
      return event.answer
    else
      handleEvent(protocol, event)
    end
    sleep(0.001)
  end
end

handleEvent(protocol::Protocol, event::StopEvent) = stop(protocol)
handleEvent(protocol::Protocol, event::ResumeEvent) = resume(protocol)
handleEvent(protocol::Protocol, event::CancelEvent) = cancel(protocol)
handleEvent(protocol::Protocol, event::ProtocolEvent) = put!(biChannel(protocol), UndefinedEvent(event))
#handleEvent(protocol::Protocol, event::InfoQueryEvent) =

function handleEvents(protocol::Protocol)
  while isready(protocol.biChannel)
    event = take!(protocol.biChannel)
    handleEvent(protocol, event)
  end
end

# Traits
abstract type ProtocolInteractivity end
struct NonInteractive <: ProtocolInteractivity end
struct Interactive <: ProtocolInteractivity end
@mustimplement protocolInteractivity(protocol::Protocol)

include("DAQMeasurementProtocol.jl")
include("MPIMeasurementProtocol.jl")
include("RobotBasedProtocol.jl")
include("RobotBasedSystemMatrixProtocol.jl")
include("ContinousMeasurementProtocol.jl")
#include("TransferFunctionProtocol.jl")


text:
using WAV

FILENAME = "british_speaking.wav"

require("k_svd")

wav, FS = wavread(FILENAME)

midpoint = int(size(wav)[1]/2)

# We want to build a dictionary off of the first half of the file
# to see how well we can emulate the sound in the second half of the file

function training_data(num_samples, sample_dimension)
    # Make sure to only get samples from the first half.
    # Second half is what we will test on.
    max_index = midpoint - sample_dimension
    starting_values = rand(1:max_index, num_samples)
    data = [wav[i:(i+sample_dimension-1)] for i in starting_values]

    hcat(data...)
end

function encode_half_2(D, sparsity, seconds=8)
    sample_length = size(D)[1]
    output = zeros(seconds * FS)
    num_frames = int(seconds * FS / sample_length)
    for i = 0:(num_frames-1)
        start_offset = i*sample_length
        end_offset = start_offset + sample_length
        signal_piece = wav[midpoint + start_offset:midpoint + end_offset-1]
        sparse_encoding = orthogonal_matching_persuit(signal_piece, D, sparsity)
        sparsely_encoded = D*sparse_encoding
        output[start_offset+1:end_offset] = sparsely_encoded
    end
    output /= maximum(abs(output[:,1]))
    wavwrite(wav[midpoint:end], "original_second_half.wav", Fs=FS)
    wavwrite(output, "encoded_second_half.wav", Fs=FS)
end

function demo()
    S = training_data(1000, 200);
    sparsity = 40
    # Slow step...
    D = k_svd(S, 30, 500, sparsity);
    encode_half_2(D, sparsity)
end

max_stars_repo_name: justinvf-zz/algorithmshop

max_issues_repo_name: justinvf-zz/algorithmshop

max_forks_repo_name: justinvf-zz/algorithmshop


text: #Websocket tester: http://jsfiddle.net/EAVvQ/24/ or chrome extension: https://chrome.google.com/webstore/detail/simple-websocket-client/

max_stars_repo_name: albapa/ulti

max_issues_repo_name: albapa/ulti

max_forks_repo_name: albapa/ulti


text:
# ---
# title: 79. Word Search
# id: problem79
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Backtracking
# link: <https://leetcode.com/problems/word-search/description/>
# hidden: true
# ---
#
# Given an `m x n` `board` and a `word`, find if the word exists in the grid.
#
# The word can be constructed from letters of sequentially adjacent cells, where
# "adjacent" cells are horizontally or vertically neighboring. The same letter
# cell may not be used more than once.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/04/word2.jpg)
#
#
#
#     Input: board = [["A","B","C","E"],["S","F","C","S"],["A","D","E","E"]], word = "ABCCED"
#     Output: true
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/04/word-1.jpg)
#
#
#
#     Input: board = [["A","B","C","E"],["S","F","C","S"],["A","D","E","E"]], word = "SEE"
#     Output: true
#
#
# **Example 3:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/15/word3.jpg)
#
#
#
#     Input: board = [["A","B","C","E"],["S","F","C","S"],["A","D","E","E"]], word = "ABCB"
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `m == board.length`
#   * `n = board[i].length`
#   * `1 <= m, n <= 200`
#   * `1 <= word.length <= 103`
#   * `board` and `word` consists only of lowercase and uppercase English letters.
#
#
## @lc code=start
using LeetCode

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


text:
module TestKVStore
using MXNet
using Base.Test

using ..Main: rand_dims

SHAPE = rand_dims()
KEYS  = [5,7,11]

function init_kv()
  kv = mx.KVStore()
  mx.init!(kv, 3, mx.zeros(SHAPE))

  vals = [mx.zeros(SHAPE) for k in KEYS]
  mx.init!(kv, KEYS, vals)
  return kv
end

function test_kv_basic()
  info("KVStore::basic")

  kv = init_kv()
  @test mx.get_type(kv) == :local
  @test mx.get_rank(kv) == 0
  @test mx.get_num_workers(kv) == 1
end

function test_single_kv_pair()
  info("KVStore::single")

  kv = init_kv()
  mx.push!(kv, 3, mx.ones(SHAPE))
  val = mx.empty(SHAPE)
  mx.pull!(kv, 3, val)
  @test maximum(abs(copy(val) - 1)) == 0
end

function test_aggregator()
  info("KVStore::aggregator")

  kv = init_kv()

  num_devs = 4
  devs = [mx.Context(mx.CPU, i) for i=0:num_devs-1]
  vals = [mx.ones(SHAPE, dev) for dev in devs]

  mx.push!(kv, 3, vals)
  mx.pull!(kv, 3, vals)
  for v in vals
    @test maximum(abs(copy(v)) - num_devs) == 0
  end

  # list
  vals = [mx.NDArray[mx.ones(SHAPE, dev)*2 for dev in devs] for k in KEYS]
  mx.push!(kv, KEYS, vals)
  mx.pull!(kv, KEYS, vals)

  for vv in vals
    for v in vv
      @test maximum(abs(copy(v)) - 2num_devs) == 0
    end
  end
end

test_kv_basic()
test_single_kv_pair()
test_aggregator()

end


text:
using TropicalNumbers
using Test, Documenter

@testset "tropical" begin
    include("tropical.jl")
end

@testset "counting_tropical" begin
    include("counting_tropical.jl")
end

doctest(TropicalNumbers)


text:
struct VerticallyStretchedRectilinearGrid{FT, TX, TY, TZ, R, A, Arch} <: AbstractRectilinearGrid{FT, TX, TY, TZ}

    architecture :: Arch

    # Number of grid points in (x,y,z).
    Nx :: Int
    Ny :: Int
    Nz :: Int

    # Halo size in (x,y,z).
    Hx :: Int
    Hy :: Int
    Hz :: Int

    # Domain size [m].
    Lx :: FT
    Ly :: FT
    Lz :: FT

    # Grid spacing [m].
       Δx :: FT
       Δy :: FT
    Δzᵃᵃᶜ :: A
    Δzᵃᵃᶠ :: A

    # Range of coordinates at the centers of the cells.
    xᶜᵃᵃ :: R
    yᵃᶜᵃ :: R
    zᵃᵃᶜ :: A

    # Range of grid coordinates at the faces of the cells.
    # Note: there are Nx+1 faces in the x-dimension, Ny+1 in the y, and Nz+1 in the z.
    xᶠᵃᵃ :: R
    yᵃᶠᵃ :: R
    zᵃᵃᶠ :: A
end

"""
    VerticallyStretchedRectilinearGrid([FT=Float64]; architecture=CPU(), size, z_faces,
                                        x = nothing, y = nothing,
                                        topology = (Periodic, Periodic, Bounded), halo = nothing)

Creates a `VerticallyStretchedRectilinearGrid` with `size = (Nx, Ny, Nz)` grid points and
a vertical grid specified by `z_faces`.

Keyword arguments
=================

- `size` (required): A tuple prescribing the number of grid points in non-`Flat` directions.
                     `size` is a 3-tuple for 3D models, a 2-tuple for 2D models, and either a
                     scalar or 1-tuple for 1D models.

- `topology`: A 3-tuple `(Tx, Ty, Tz)` specifying the topology of the domain.
              `Tz` must be `Bounded` for `VerticallyStretchedRectilinearGrid`.
              `Tx` and `Ty` specify whether the `x`- and `y`-` directions are
              `Periodic`, `Bounded`, or `Flat`. The topology `Flat` indicates that a model does
              not vary in that directions so that derivatives and interpolation are zero.
              The default is `topology=(Periodic, Periodic, Bounded)`.

- `architecture`: Specifies whether the array of vertical coordinates, interfaces, and spacings
                  are stored on the CPU or GPU. Default: `architecture = CPU()`.

- `z_faces`: An array or function of vertical index `k` that specifies the location of cell faces
        in the `z-`direction for indices `k=1` through `k=Nz+1`, where `Nz` is the third element
        of `size`.

- `x`, `y`: Each of `x, y` are 2-tuples that specify the end points of the domain
            in their respect directions. Scalar values may be used in `Flat` directions.

- `halo`: A tuple of integers that specifies the size of the halo region of cells surrounding
          the physical interior for each non-`Flat` direction.

The physical extent of the domain can be specified via `x` and `y` keyword arguments
indicating the left and right endpoints of each dimensions, e.g. `x=(-π, π)`.

A grid topology may be specified via a tuple assigning one of `Periodic`, `Bounded, and `Flat`
to each dimension. By default, a horizontally periodic grid topology `(Periodic, Periodic, Bounded)`
is assumed.

Constants are stored using floating point values of type `FT`. By default this is `Float64`.
Make sure to specify the desired `FT` if not using `Float64`.

Grid properties
===============

- `(Nx, Ny, Nz)::Int`: Number of physical points in the (x, y, z)-direction

- `(Hx, Hy, Hz)::Int`: Number of halo points in the (x, y, z)-direction

- `(Lx, Ly, Lz)::FT`: Physical extent of the grid in the (x, y, z)-direction

- `(Δx, Δy)::FT`: Grid spacing (distance between grid nodes) in the (x, y)-direction

- `Δzᵃᵃᶜ`: Grid spacing in the z-direction between cell faces.
           Defined at cell centers in `z` and independent of cell location in (x, y).

- `Δzᵃᵃᶠ`: Grid spacing in the z-direction between cell centers, and defined at cell faces in z.
           Defined at cell faces in `z` and independent of cell location in (x, y).

- `(xᶜᵃᵃ, yᵃᶜᵃ, zᵃᵃᶜ)`: (x, y, z) coordinates of cell centers.

- `(xᶠᵃᵃ, yᵃᶠᵃ, zᵃᵃᶠ)`: (x, y, z) coordinates of cell faces.
"""
function VerticallyStretchedRectilinearGrid(FT = Float64;
                                            architecture = CPU(),
                                            size,
                                            z_faces,
                                            x = nothing,
                                            y = nothing,
                                            halo = nothing,
                                            topology = (Periodic, Periodic, Bounded))

    TX, TY, TZ = validate_topology(topology)
    size = validate_size(TX, TY, TZ, size)
    halo = validate_halo(TX, TY, TZ, halo)
    x = validate_dimension_specification(TX, x, :x)
    y = validate_dimension_specification(TY, y, :y)
    Lx, Ly, x, y = validate_vertically_stretched_grid_xy(TX, TY, FT, x, y)

    Nx, Ny, Nz = size
    Hx, Hy, Hz = halo

    # Initialize vertically-stretched arrays on CPU
    Lz, zᵃᵃᶠ, zᵃᵃᶜ, Δzᵃᵃᶜ, Δzᵃᵃᶠ = generate_stretched_vertical_grid(FT, topology[3], Nz, Hz, z_faces)

    # Construct uniform horizontal grid
    Lh, Nh, Hh, X₁ = (Lx, Ly), size[1:2], halo[1:2], (x[1], y[1])
    Δx, Δy = Δh = Lh ./ Nh

    # Face-node limits in x, y, z
    xF₋, yF₋ = XF₋ = @. X₁ - Hh * Δh
    xF₊, yF₊ = XF₊ = @. XF₋ + total_extent(topology[1:2], Hh, Δh, Lh)

    # Center-node limits in x, y, z
    xC₋, yC₋ = XC₋ = @. XF₋ + Δh / 2
    xC₊, yC₊ = XC₊ = @. XC₋ + Lh + Δh * (2Hh - 1)

    # Total length of Center and Face quantities
    TFx, TFy, TFz = total_length.(Face, topology, size, halo)
    TCx, TCy, TCz = total_length.(Center, topology, size, halo)

    # Include halo points in coordinate arrays
    xᶠᵃᵃ = range(xF₋, xF₊; length = TFx)
    yᵃᶠᵃ = range(yF₋, yF₊; length = TFy)

    xᶜᵃᵃ = range(xC₋, xC₊; length = TCx)
    yᵃᶜᵃ = range(yC₋, yC₊; length = TCy)

    xᶜᵃᵃ = OffsetArray(xᶜᵃᵃ,  -Hx)
    yᵃᶜᵃ = OffsetArray(yᵃᶜᵃ,  -Hy)
    zᵃᵃᶜ = OffsetArray(zᵃᵃᶜ,  -Hz)

    xᶠᵃᵃ = OffsetArray(xᶠᵃᵃ,  -Hx)
    yᵃᶠᵃ = OffsetArray(yᵃᶠᵃ,  -Hy)
    zᵃᵃᶠ = OffsetArray(zᵃᵃᶠ,  -Hz)

    Δzᵃᵃᶠ = OffsetArray(Δzᵃᵃᶠ, -Hz)
    Δzᵃᵃᶜ = OffsetArray(Δzᵃᵃᶜ, -Hz)

    # Seems needed to avoid out-of-bounds error in viscous dissipation
    # operators wanting to access Δzᵃᵃᶠ[Nz+2].
    Δzᵃᵃᶠ = OffsetArray(cat(Δzᵃᵃᶠ[0], Δzᵃᵃᶠ..., Δzᵃᵃᶠ[Nz], dims=1), -Hz-1)

    # Convert to appropriate array type for arch
    zᵃᵃᶠ  = OffsetArray(arch_array(architecture,  zᵃᵃᶠ.parent),  zᵃᵃᶠ.offsets...)
    zᵃᵃᶜ  = OffsetArray(arch_array(architecture,  zᵃᵃᶜ.parent),  zᵃᵃᶜ.offsets...)
    Δzᵃᵃᶜ = OffsetArray(arch_array(architecture, Δzᵃᵃᶜ.parent), Δzᵃᵃᶜ.offsets...)
    Δzᵃᵃᶠ = OffsetArray(arch_array(architecture, Δzᵃᵃᶠ.parent), Δzᵃᵃᶠ.offsets...)

    R = typeof(xᶠᵃᵃ)
    A = typeof(zᵃᵃᶠ)
    Arch = typeof(architecture)

    return VerticallyStretchedRectilinearGrid{FT, TX, TY, TZ, R, A, Arch}(architecture,
        Nx, Ny, Nz, Hx, Hy, Hz, Lx, Ly, Lz, Δx, Δy, Δzᵃᵃᶜ, Δzᵃᵃᶠ, xᶜᵃᵃ, yᵃᶜᵃ, zᵃᵃᶜ, xᶠᵃᵃ, yᵃᶠᵃ, zᵃᵃᶠ)
end

#####
##### Vertically stretched grid utilities
#####

get_z_face(z::Function, k) = z(k)
get_z_face(z::AbstractVector, k) = CUDA.@allowscalar z[k]

lower_exterior_Δzᵃᵃᶜ(z_topo,          zFi, Hz) = [zFi[end - Hz + k] - zFi[end - Hz + k - 1] for k = 1:Hz]
lower_exterior_Δzᵃᵃᶜ(::Type{Bounded}, zFi, Hz) = [zFi[2]  - zFi[1] for k = 1:Hz]

upper_exterior_Δzᵃᵃᶜ(z_topo,          zFi, Hz) = [zFi[k + 1] - zFi[k] for k = 1:Hz]
upper_exterior_Δzᵃᵃᶜ(::Type{Bounded}, zFi, Hz) = [zFi[end]   - zFi[end - 1] for k = 1:Hz]

function generate_stretched_vertical_grid(FT, z_topo, Nz, Hz, z_faces)

    # Ensure correct type for zF and derived quantities
    interior_zF = zeros(FT, Nz+1)

    for k = 1:Nz+1
        interior_zF[k] = get_z_face(z_faces, k)
    end

    Lz = interior_zF[Nz+1] - interior_zF[1]

    # Build halo regions
    ΔzF₋ = lower_exterior_Δzᵃᵃᶜ(z_topo, interior_zF, Hz)
    ΔzF₊ = upper_exterior_Δzᵃᵃᶜ(z_topo, interior_zF, Hz)

    z¹, zᴺ⁺¹ = interior_zF[1], interior_zF[Nz+1]

    zF₋ = [z¹   - sum(ΔzF₋[k:Hz]) for k = 1:Hz] # locations of faces in lower halo
    zF₊ = reverse([zᴺ⁺¹ + sum(ΔzF₊[k:Hz]) for k = 1:Hz]) # locations of faces in width of top halo region

    zF = vcat(zF₋, interior_zF, zF₊)

    # Build cell centers, cell center spacings, and cell interface spacings
    TCz = total_length(Center, z_topo, Nz, Hz)
     zC = [ (zF[k + 1] + zF[k]) / 2 for k = 1:TCz ]
    ΔzC = [  zC[k] - zC[k - 1]      for k = 2:TCz ]

    # Trim face locations for periodic domains
    TFz = total_length(Face, z_topo, Nz, Hz)
    zF = zF[1:TFz]

    ΔzF = [zF[k + 1] - zF[k] for k = 1:TFz-1]

    return Lz, zF, zC, ΔzF, ΔzC
end

"""
    with_halo(new_halo, old_grid::VerticallyStretchedRectilinearGrid)

Returns a new `VerticallyStretchedRectilinearGrid` with the same properties as
`old_grid` but with halos set to `new_halo`.

Note that in contrast to the constructor for `VerticallyStretchedRectilinearGrid`,
`new_halo` is expected to be a 3-`Tuple` by `with_halo`. The elements
of `new_halo` corresponding to `Flat` directions are removed (and are
therefore ignored) prior to constructing the new `VerticallyStretchedRectilinearGrid`.
"""
function with_halo(new_halo, old_grid::VerticallyStretchedRectilinearGrid)

    Nx, Ny, Nz = size = (old_grid.Nx, old_grid.Ny, old_grid.Nz)
    topo = topology(old_grid)

    x = x_domain(old_grid)
    y = y_domain(old_grid)
    z = z_domain(old_grid)

    # Remove elements of size and new_halo in Flat directions as expected by grid
    # constructor
    size = pop_flat_elements(size, topo)
    new_halo = pop_flat_elements(new_halo, topo)

    new_grid = VerticallyStretchedRectilinearGrid(eltype(old_grid);
                                                  architecture = old_grid.architecture,
                                                  size = size,
                                                  x = x, y = y,
                                                  z_faces = old_grid.zᵃᵃᶠ,
                                                  topology = topo,
                                                  halo = new_halo)

    return new_grid
end

@inline x_domain(grid::VerticallyStretchedRectilinearGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} = domain(TX, grid.Nx, grid.xᶠᵃᵃ)
@inline y_domain(grid::VerticallyStretchedRectilinearGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} = domain(TY, grid.Ny, grid.yᵃᶠᵃ)
@inline z_domain(grid::VerticallyStretchedRectilinearGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} = domain(TZ, grid.Nz, grid.zᵃᵃᶠ)

short_show(grid::VerticallyStretchedRectilinearGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} =
    "VerticallyStretchedRectilinearGrid{$FT, $TX, $TY, $TZ}(Nx=$(grid.Nx), Ny=$(grid.Ny), Nz=$(grid.Nz))"

function show(io::IO, g::VerticallyStretchedRectilinearGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ}
    Δz_min = minimum(view(parent(g.Δzᵃᵃᶜ), g.Hz+1:g.Nz+g.Hz))
    Δz_max = maximum(view(parent(g.Δzᵃᵃᶜ), g.Hz+1:g.Nz+g.Hz))
    print(io, "VerticallyStretchedRectilinearGrid{$FT, $TX, $TY, $TZ}\n",
              "                   domain: $(domain_string(g))\n",
              "                 topology: ", (TX, TY, TZ), '\n',
              "  resolution (Nx, Ny, Nz): ", (g.Nx, g.Ny, g.Nz), '\n',
              "   halo size (Hx, Hy, Hz): ", (g.Hx, g.Hy, g.Hz), '\n',
              "grid spacing (Δx, Δy, Δz): (", g.Δx, ", ", g.Δy, ", [min=", Δz_min, ", max=", Δz_max,"])",)
end

Adapt.adapt_structure(to, grid::VerticallyStretchedRectilinearGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} =
    VerticallyStretchedRectilinearGrid{FT, TX, TY, TZ,
                                       typeof(grid.xᶠᵃᵃ),
                                       typeof(Adapt.adapt(to, grid.zᵃᵃᶠ)),
                                       Nothing}(
        nothing,
        grid.Nx, grid.Ny, grid.Nz,
        grid.Hx, grid.Hy, grid.Hz,
        grid.Lx, grid.Ly, grid.Lz,
        grid.Δx, grid.Δy,
        Adapt.adapt(to, grid.Δzᵃᵃᶜ),
        Adapt.adapt(to, grid.Δzᵃᵃᶠ),
        grid.xᶜᵃᵃ, grid.yᵃᶜᵃ,
        Adapt.adapt(to, grid.zᵃᵃᶜ),
        grid.xᶠᵃᵃ, grid.yᵃᶠᵃ,
        Adapt.adapt(to, grid.zᵃᵃᶠ))

#####
##### Should merge with grid_utils.jl at some point
#####

@inline xnode(::Center, i, grid::VerticallyStretchedRectilinearGrid) = @inbounds grid.xᶜᵃᵃ[i]
@inline xnode(::Face, i, grid::VerticallyStretchedRectilinearGrid) = @inbounds grid.xᶠᵃᵃ[i]

@inline ynode(::Center, j, grid::VerticallyStretchedRectilinearGrid) = @inbounds grid.yᵃᶜᵃ[j]
@inline ynode(::Face, j, grid::VerticallyStretchedRectilinearGrid) = @inbounds grid.yᵃᶠᵃ[j]

@inline znode(::Center, k, grid::VerticallyStretchedRectilinearGrid) = @inbounds grid.zᵃᵃᶜ[k]
@inline znode(::Face, k, grid::VerticallyStretchedRectilinearGrid) = @inbounds grid.zᵃᵃᶠ[k]

all_x_nodes(::Type{Center}, grid::VerticallyStretchedRectilinearGrid) = grid.xᶜᵃᵃ
all_x_nodes(::Type{Face}, grid::VerticallyStretchedRectilinearGrid) = grid.xᶠᵃᵃ
all_y_nodes(::Type{Center}, grid::VerticallyStretchedRectilinearGrid) = grid.yᵃᶜᵃ
all_y_nodes(::Type{Face}, grid::VerticallyStretchedRectilinearGrid) = grid.yᵃᶠᵃ
all_z_nodes(::Type{Center}, grid::VerticallyStretchedRectilinearGrid) = grid.zᵃᵃᶜ
all_z_nodes(::Type{Face}, grid::VerticallyStretchedRectilinearGrid) = grid.zᵃᵃᶠ

#
# Get minima of grid
#

function min_Δx(grid::VerticallyStretchedRectilinearGrid)
    topo = topology(grid)
    if topo[1] == Flat
        return Inf
    else
        return grid.Δx
    end
end

function min_Δy(grid::VerticallyStretchedRectilinearGrid)
    topo = topology(grid)
    if topo[2] == Flat
        return Inf
    else
        return grid.Δy
    end
end

function min_Δz(grid::VerticallyStretchedRectilinearGrid)
    topo = topology(grid)
    if topo[3] == Flat
        return Inf
    else
        return minimum(parent(grid.Δzᵃᵃᶜ))
    end
end