text:
const EMPTYSTRING = ""

# Implementation of MOI for AbstractModel
abstract type AbstractModelLike{T} <: MOI.ModelLike end
abstract type AbstractOptimizer{T} <: MOI.AbstractOptimizer end
const AbstractModel{T} = Union{AbstractModelLike{T},AbstractOptimizer{T}}

# Variables
function MOI.get(model::AbstractModel, ::MOI.NumberOfVariables)::Int64
    if model.variable_indices === nothing
        return model.num_variables_created
    else
        return length(model.variable_indices)
    end
end

function MOI.add_variable(model::AbstractModel{T}) where {T}
    vi = VI(model.num_variables_created += 1)
    push!(model.single_variable_mask, 0x0)
    push!(model.lower_bound, zero(T))
    push!(model.upper_bound, zero(T))
    if model.variable_indices !== nothing
        push!(model.variable_indices, vi)
    end
    return vi
end

function MOI.add_variables(model::AbstractModel, n::Integer)
    return [MOI.add_variable(model) for i in 1:n]
end

"""
    remove_variable(f::MOI.AbstractFunction, s::MOI.AbstractSet, vi::MOI.VariableIndex)

Return a tuple `(g, t)` representing the constraint `f`-in-`s` with the
variable `vi` removed. That is, the terms containing the variable `vi` in the
function `f` are removed and the dimension of the set `s` is updated if
needed (e.g. when `f` is a `VectorOfVariables` with `vi` being one of the
variables).
"""
remove_variable(f, s, vi::VI) = remove_variable(f, vi), s
function remove_variable(f::MOI.VectorOfVariables, s, vi::VI)
    g = remove_variable(f, vi)
    if length(g.variables) != length(f.variables)
        t = MOI.update_dimension(s, length(g.variables))
    else
        t = s
    end
    return g, t
end

filter_variables(keep::F, f, s) where {F<:Function} = filter_variables(keep, f), s
function filter_variables(keep::F, f::MOI.VectorOfVariables, s) where {F<:Function}
    g = filter_variables(keep, f)
    if length(g.variables) != length(f.variables)
        t = MOI.update_dimension(s, length(g.variables))
    else
        t = s
    end
    return g, t
end
function _delete_variable(
    model::AbstractModel{T},
    vi::MOI.VariableIndex,
) where {T}
    MOI.throw_if_not_valid(model, vi)
    model.single_variable_mask[vi.value] = 0x0
    if model.variable_indices === nothing
        model.variable_indices =
            Set(MOI.get(model, MOI.ListOfVariableIndices()))
    end
    delete!(model.variable_indices, vi)
    model.name_to_var = nothing
    delete!(model.var_to_name, vi)
    model.name_to_con = nothing
    delete!(
        model.con_to_name,
        MOI.ConstraintIndex{MOI.SingleVariable,MOI.EqualTo{T}}(vi.value),
    )
    delete!(
        model.con_to_name,
        MOI.ConstraintIndex{MOI.SingleVariable,MOI.GreaterThan{T}}(vi.value),
    )
    delete!(
        model.con_to_name,
        MOI.ConstraintIndex{MOI.SingleVariable,MOI.LessThan{T}}(vi.value),
    )
    delete!(
        model.con_to_name,
        MOI.ConstraintIndex{MOI.SingleVariable,MOI.Interval{T}}(vi.value),
    )
    delete!(
        model.con_to_name,
        MOI.ConstraintIndex{MOI.SingleVariable,MOI.Integer}(vi.value),
    )
    delete!(
        model.con_to_name,
        MOI.ConstraintIndex{MOI.SingleVariable,MOI.ZeroOne}(vi.value),
    )
    delete!(
        model.con_to_name,
        MOI.ConstraintIndex{MOI.SingleVariable,MOI.Semicontinuous{T}}(vi.value),
    )
    return delete!(
        model.con_to_name,
        MOI.ConstraintIndex{MOI.SingleVariable,MOI.Semiinteger{T}}(vi.value),
    )
end
function MOI.delete(model::AbstractModel, vi::MOI.VariableIndex)
    vis = [vi]
    broadcastcall(model) do constrs
        _throw_if_cannot_delete(constrs, vis, vis)
    end
    _delete_variable(model, vi)
    broadcastcall(model) do constrs
        _deleted_constraints(constrs, vi) do ci
            delete!(model.con_to_name, ci)
        end
    end
    model.objective = remove_variable(model.objective, vi)
    model.name_to_con = nothing
    return
end

function MOI.delete(model::AbstractModel, vis::Vector{MOI.VariableIndex})
    if isempty(vis)
        # In `keep`, we assume that `model.variable_indices !== nothing` so
        # at least one variable need to be deleted.
        return
    end
    fast_in_vis = Set(vis)
    broadcastcall(model) do constrs
        _throw_if_cannot_delete(constrs, vis, fast_in_vis)
    end
    broadcastcall(model) do constrs
        _deleted_constraints(constrs, vis) do ci
            delete!(model.con_to_name, ci)
        end
    end
    for vi in vis
        _delete_variable(model, vi)
    end
    keep(vi::MOI.VariableIndex) = vi in model.variable_indices
    model.objective = filter_variables(keep, model.objective)
    model.name_to_con = nothing
    return
end

function MOI.is_valid(
    model::AbstractModel,
    ci::CI{MOI.SingleVariable,S},
) where {S}
    return 1 ≤ ci.value ≤ length(model.single_variable_mask) &&
           !iszero(
        model.single_variable_mask[ci.value] & single_variable_flag(S),
    )
end

function MOI.is_valid(model::AbstractModel, ci::CI{F,S}) where {F,S}
    if MOI.supports_constraint(model, F, S)
        return MOI.is_valid(constraints(model, ci), ci)
    else
        return false
    end
end

function MOI.is_valid(model::AbstractModel, vi::VI)
    if model.variable_indices === nothing
        return 1 ≤ vi.value ≤ model.num_variables_created
    else
        return in(vi, model.variable_indices)
    end
end

function MOI.get(model::AbstractModel, ::MOI.ListOfVariableIndices)
    if model.variable_indices === nothing
        return VI.(1:model.num_variables_created)
    else
        vis = collect(model.variable_indices)
        sort!(vis, by = vi -> vi.value) # It needs to be sorted by order of creation
        return vis
    end
end

# Names
MOI.supports(::AbstractModel, ::MOI.Name) = true
function MOI.set(model::AbstractModel, ::MOI.Name, name::String)
    return model.name = name
end
MOI.get(model::AbstractModel, ::MOI.Name) = model.name

MOI.supports(::AbstractModel, ::MOI.VariableName, vi::Type{VI}) = true
function MOI.set(model::AbstractModel, ::MOI.VariableName, vi::VI, name::String)
    model.var_to_name[vi] = name
    model.name_to_var = nothing # Invalidate the name map.
    return
end

function MOI.get(model::AbstractModel, ::MOI.VariableName, vi::VI)
    return get(model.var_to_name, vi, EMPTYSTRING)
end

"""
    build_name_to_var_map(con_to_name::Dict{MOI.VariableIndex, String})

Create and return a reverse map from name to variable index, given a map from
variable index to name. The special value `MOI.VariableIndex(0)` is used to
indicate that multiple variables have the same name.
"""
function build_name_to_var_map(var_to_name::Dict{VI,String})
    name_to_var = Dict{String,VI}()
    for (var, var_name) in var_to_name
        if haskey(name_to_var, var_name)
            # 0 is a special value that means this string does not map to
            # a unique variable name.
            name_to_var[var_name] = VI(0)
        else
            name_to_var[var_name] = var
        end
    end
    return name_to_var
end

function throw_multiple_name_error(::Type{MOI.VariableIndex}, name::String)
    return error("Multiple variables have the name $name.")
end
function throw_multiple_name_error(::Type{<:MOI.ConstraintIndex}, name::String)
    return error("Multiple constraints have the name $name.")
end
function throw_if_multiple_with_name(::Nothing, ::String) end
function throw_if_multiple_with_name(index::MOI.Index, name::String)
    if iszero(index.value)
        throw_multiple_name_error(typeof(index), name)
    end
end

function MOI.get(model::AbstractModel, ::Type{VI}, name::String)
    if model.name_to_var === nothing
        # Rebuild the map.
        model.name_to_var = build_name_to_var_map(model.var_to_name)
    end
    result = get(model.name_to_var, name, nothing)
    throw_if_multiple_with_name(result, name)
    return result
end

function MOI.get(
    model::AbstractModel,
    ::MOI.ListOfVariableAttributesSet,
)::Vector{MOI.AbstractVariableAttribute}
    return isempty(model.var_to_name) ? [] : [MOI.VariableName()]
end

MOI.supports(model::AbstractModel, ::MOI.ConstraintName, ::Type{<:CI}) = true
function MOI.set(
    model::AbstractModel,
    ::MOI.ConstraintName,
    ci::CI,
    name::String,
)
    model.con_to_name[ci] = name
    model.name_to_con = nothing # Invalidate the name map.
    return
end

function MOI.get(model::AbstractModel, ::MOI.ConstraintName, ci::CI)
    return get(model.con_to_name, ci, EMPTYSTRING)
end

"""
    build_name_to_con_map(con_to_name::Dict{MOI.ConstraintIndex, String})

Create and return a reverse map from name to constraint index, given a map from
constraint index to name. The special value
`MOI.ConstraintIndex{Nothing, Nothing}(0)` is used to indicate that multiple
constraints have the same name.
"""
function build_name_to_con_map(con_to_name::Dict{CI,String})
    name_to_con = Dict{String,CI}()
    for (con, con_name) in con_to_name
        if haskey(name_to_con, con_name)
            name_to_con[con_name] = CI{Nothing,Nothing}(0)
        else
            name_to_con[con_name] = con
        end
    end
    return name_to_con
end

function MOI.get(model::AbstractModel, ConType::Type{<:CI}, name::String)
    if model.name_to_con === nothing
        # Rebuild the map.
        model.name_to_con = build_name_to_con_map(model.con_to_name)
    end
    ci = get(model.name_to_con, name, nothing)
    throw_if_multiple_with_name(ci, name)
    return ci isa ConType ? ci : nothing
end

function MOI.get(
    model::AbstractModel,
    ::MOI.ListOfConstraintAttributesSet,
)::Vector{MOI.AbstractConstraintAttribute}
    return isempty(model.con_to_name) ? [] : [MOI.ConstraintName()]
end

# Objective
MOI.get(model::AbstractModel, ::MOI.ObjectiveSense) = model.sense
MOI.supports(model::AbstractModel, ::MOI.ObjectiveSense) = true
function MOI.set(
    model::AbstractModel{T},
    ::MOI.ObjectiveSense,
    sense::MOI.OptimizationSense,
) where {T}
    if sense == MOI.FEASIBILITY_SENSE
        model.objectiveset = false
        model.objective = zero(MOI.ScalarAffineFunction{T})
    end
    model.senseset = true
    model.sense = sense
    return
end

function MOI.get(model::AbstractModel, ::MOI.ObjectiveFunctionType)
    return MOI.typeof(model.objective)
end
function MOI.get(model::AbstractModel, ::MOI.ObjectiveFunction{T})::T where {T}
    return model.objective
end
function MOI.supports(
    model::AbstractModel{T},
    ::MOI.ObjectiveFunction{<:Union{
        MOI.SingleVariable,
        MOI.ScalarAffineFunction{T},
        MOI.ScalarQuadraticFunction{T},
    }},
) where {T}
    return true
end
function MOI.set(
    model::AbstractModel,
    attr::MOI.ObjectiveFunction{F},
    f::F,
) where {F<:MOI.AbstractFunction}
    if !MOI.supports(model, attr)
        throw(MOI.UnsupportedAttribute(attr))
    end
    model.objectiveset = true
    # f needs to be copied, see #2
    model.objective = copy(f)
    return
end

function MOI.modify(
    model::AbstractModel,
    ::MOI.ObjectiveFunction,
    change::MOI.AbstractFunctionModification,
)
    model.objective = modify_function(model.objective, change)
    model.objectiveset = true
    return
end

function MOI.get(::AbstractModel, ::MOI.ListOfOptimizerAttributesSet)
    return MOI.AbstractOptimizerAttribute[]
end
function MOI.get(
    model::AbstractModel,
    ::MOI.ListOfModelAttributesSet,
)::Vector{MOI.AbstractModelAttribute}
    listattr = MOI.AbstractModelAttribute[]
    if model.senseset
        push!(listattr, MOI.ObjectiveSense())
    end
    if model.objectiveset
        push!(listattr, MOI.ObjectiveFunction{typeof(model.objective)}())
    end
    if !isempty(model.name)
        push!(listattr, MOI.Name())
    end
    return listattr
end

# Constraints
single_variable_flag(::Type{<:MOI.EqualTo}) = 0x1
single_variable_flag(::Type{<:MOI.GreaterThan}) = 0x2
single_variable_flag(::Type{<:MOI.LessThan}) = 0x4
single_variable_flag(::Type{<:MOI.Interval}) = 0x8
single_variable_flag(::Type{MOI.Integer}) = 0x10
single_variable_flag(::Type{MOI.ZeroOne}) = 0x20
single_variable_flag(::Type{<:MOI.Semicontinuous}) = 0x40
single_variable_flag(::Type{<:MOI.Semiinteger}) = 0x80
# If a set is added here, a line should be added in
# `MOI.delete(::AbstractModel, ::MOI.VariableIndex)`

function flag_to_set_type(flag::UInt8, ::Type{T}) where {T}
    if flag == 0x1
        return MOI.EqualTo{T}
    elseif flag == 0x2
        return MOI.GreaterThan{T}
    elseif flag == 0x4
        return MOI.LessThan{T}
    elseif flag == 0x8
        return MOI.Interval{T}
    elseif flag == 0x10
        return MOI.Integer
    elseif flag == 0x20
        return MOI.ZeroOne
    elseif flag == 0x40
        return MOI.Semicontinuous{T}
    else
        @assert flag == 0x80
        return MOI.Semiinteger{T}
    end
end

# Julia doesn't infer `S1` correctly, so we use a function barrier to improve
# inference.
function _throw_if_lower_bound_set(variable, S2, mask, T)
    S1 = flag_to_set_type(mask, T)
    throw(MOI.LowerBoundAlreadySet{S1,S2}(variable))
    return
end

function throw_if_lower_bound_set(variable, S2, mask, T)
    lower_mask = mask & LOWER_BOUND_MASK
    if iszero(lower_mask)
        return  # No lower bound set.
    elseif iszero(single_variable_flag(S2) & LOWER_BOUND_MASK)
        return  # S2 isn't related to the lower bound.
    end
    return _throw_if_lower_bound_set(variable, S2, lower_mask, T)
end

# Julia doesn't infer `S1` correctly, so we use a function barrier to improve
# inference.
function _throw_if_upper_bound_set(variable, S2, mask, T)
    S1 = flag_to_set_type(mask, T)
    throw(MOI.UpperBoundAlreadySet{S1,S2}(variable))
    return
end

function throw_if_upper_bound_set(variable, S2, mask, T)
    upper_mask = mask & UPPER_BOUND_MASK
    if iszero(upper_mask)
        return  # No upper bound set.
    elseif iszero(single_variable_flag(S2) & UPPER_BOUND_MASK)
        return  # S2 isn't related to the upper bound.
    end
    return _throw_if_upper_bound_set(variable, S2, upper_mask, T)
end

# Sets setting lower bound:
extract_lower_bound(set::MOI.EqualTo) = set.value
function extract_lower_bound(
    set::Union{MOI.GreaterThan,MOI.Interval,MOI.Semicontinuous,MOI.Semiinteger},
)
    return set.lower
end
# 0xcb = 0x80 | 0x40 | 0x8 | 0x2 | 0x1
const LOWER_BOUND_MASK = 0xcb

# Sets setting upper bound:
extract_upper_bound(set::MOI.EqualTo) = set.value
function extract_upper_bound(
    set::Union{MOI.LessThan,MOI.Interval,MOI.Semicontinuous,MOI.Semiinteger},
)
    return set.upper
end
# 0xcd = 0x80 | 0x40 | 0x8 | 0x4 | 0x1
const UPPER_BOUND_MASK = 0xcd

const SUPPORTED_VARIABLE_SCALAR_SETS{T} = Union{
    MOI.EqualTo{T},
    MOI.GreaterThan{T},
    MOI.LessThan{T},
    MOI.Interval{T},
    MOI.Integer,
    MOI.ZeroOne,
    MOI.Semicontinuous{T},
    MOI.Semiinteger{T},
}
function MOI.supports_constraint(
    ::AbstractModel{T},
    ::Type{MOI.SingleVariable},
    ::Type{<:SUPPORTED_VARIABLE_SCALAR_SETS{T}},
) where {T}
    return true
end

function MOI.add_constraint(
    model::AbstractModel{T},
    f::MOI.SingleVariable,
    s::SUPPORTED_VARIABLE_SCALAR_SETS{T},
) where {T}
    flag = single_variable_flag(typeof(s))
    index = f.variable.value
    mask = model.single_variable_mask[index]
    throw_if_lower_bound_set(f.variable, typeof(s), mask, T)
    throw_if_upper_bound_set(f.variable, typeof(s), mask, T)
    # No error should be thrown now, we can modify `model`.
    if !iszero(flag & LOWER_BOUND_MASK)
        model.lower_bound[index] = extract_lower_bound(s)
    end
    if !iszero(flag & UPPER_BOUND_MASK)
        model.upper_bound[index] = extract_upper_bound(s)
    end
    model.single_variable_mask[index] = mask | flag
    return CI{MOI.SingleVariable,typeof(s)}(index)
end

function MOI.add_constraint(model::AbstractModel, func::MOI.AbstractFunction, set::MOI.AbstractSet)
    if MOI.supports_constraint(model, typeof(func), typeof(set))
        return MOI.add_constraint(constraints(model, typeof(func), typeof(set)), func, set)
    else
        throw(MOI.UnsupportedConstraint{typeof(func),typeof(set)}())
    end
end
function constraints(
    model::AbstractModel,
    ci::MOI.ConstraintIndex{F,S}
) where {F,S}
    if !MOI.supports_constraint(model, F, S)
        throw(MOI.InvalidIndex(ci))
    end
    return constraints(model, F, S)
end
function MOI.get(model::AbstractModel, attr::Union{MOI.AbstractFunction, MOI.AbstractSet}, ci::MOI.ConstraintIndex)
    return MOI.get(constraints(model, ci), attr, ci)
end
function MOI.modify(model::AbstractModel, ci::MOI.ConstraintIndex, change)
    return MOI.modify(constraints(model, ci), ci, change)
end

function _delete_constraint(
    model::AbstractModel,
    ci::MOI.ConstraintIndex{MOI.SingleVariable,S},
) where {S}
    MOI.throw_if_not_valid(model, ci)
    model.single_variable_mask[ci.value] &= ~single_variable_flag(S)
    return
end

function _delete_constraint(model::AbstractModel, ci::MOI.ConstraintIndex)
    MOI.delete(constraints(model, ci), ci)
    return
end

function MOI.delete(model::AbstractModel, ci::MOI.ConstraintIndex)
    _delete_constraint(model, ci)
    model.name_to_con = nothing
    delete!(model.con_to_name, ci)
    return
end

function MOI.modify(
    model::AbstractModel,
    ci::MOI.ConstraintIndex,
    change::MOI.AbstractFunctionModification,
)
    MOI.modify(constraints(model, ci), ci, change)
    return
end

function MOI.set(
    ::AbstractModel,
    ::MOI.ConstraintFunction,
    ::MOI.ConstraintIndex{MOI.SingleVariable,<:MOI.AbstractScalarSet},
    ::MOI.SingleVariable,
)
    return throw(MOI.SettingSingleVariableFunctionNotAllowed())
end
function MOI.set(
    model::AbstractModel{T},
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{MOI.SingleVariable,S},
    set::S,
) where {T,S<:SUPPORTED_VARIABLE_SCALAR_SETS{T}}
    MOI.throw_if_not_valid(model, ci)
    flag = single_variable_flag(typeof(set))
    if !iszero(flag & LOWER_BOUND_MASK)
        model.lower_bound[ci.value] = extract_lower_bound(set)
    end
    if !iszero(flag & UPPER_BOUND_MASK)
        model.upper_bound[ci.value] = extract_upper_bound(set)
    end
    return
end

function MOI.set(
    model::AbstractModel,
    attr::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{<:MOI.AbstractFunction,S},
    set::S,
) where {S<:MOI.AbstractSet}
    MOI.set(constraints(model, ci), attr, ci, set)
    return
end

function MOI.set(
    model::AbstractModel,
    attr::MOI.ConstraintFunction,
    ci::MOI.ConstraintIndex{F,<:MOI.AbstractSet},
    func::F,
) where {F<:MOI.AbstractFunction}
    MOI.set(constraints(model, ci), attr, ci, func)
    return
end

function MOI.get(
    model::AbstractModel,
    ::MOI.NumberOfConstraints{MOI.SingleVariable,S},
) where {S}
    flag = single_variable_flag(S)
    return count(mask -> !iszero(flag & mask), model.single_variable_mask)
end
function MOI.get(model::AbstractModel, noc::MOI.NumberOfConstraints{F,S}) where {F,S}
    if MOI.supports_constraint(model, F, S)
        return MOI.get(constraints(model, F, S), noc)
    else
        return 0
    end
end

function _add_contraint_type(
    list,
    model::AbstractModel,
    S::Type{<:MOI.AbstractScalarSet},
)
    flag = single_variable_flag(S)
    if any(mask -> !iszero(flag & mask), model.single_variable_mask)
        push!(list, (MOI.SingleVariable, S))
    end
    return
end
function MOI.get(model::AbstractModel{T}, loc::MOI.ListOfConstraints) where {T}
    list = broadcastvcat(model) do v
        MOI.get(v, loc)
    end
    for S in (
        MOI.EqualTo{T},
        MOI.GreaterThan{T},
        MOI.LessThan{T},
        MOI.Interval{T},
        MOI.Semicontinuous{T},
        MOI.Semiinteger{T},
        MOI.Integer,
        MOI.ZeroOne,
    )
        _add_contraint_type(list, model, S)
    end
    return list
end

function MOI.get(
    model::AbstractModel,
    ::MOI.ListOfConstraintIndices{MOI.SingleVariable,S},
) where {S}
    list = CI{MOI.SingleVariable,S}[]
    flag = single_variable_flag(S)
    for (index, mask) in enumerate(model.single_variable_mask)
        if !iszero(mask & flag)
            push!(list, CI{MOI.SingleVariable,S}(index))
        end
    end
    return list
end

function MOI.get(model::AbstractModel, loc::MOI.ListOfConstraintIndices{F,S}) where {F,S}
    if MOI.supports_constraint(model, F, S)
        return MOI.get(constraints(model, F, S), loc)
    else
        return MOI.ConstraintIndex{F,S}[]
    end
end

function MOI.get(
    model::AbstractModel,
    ::MOI.ConstraintFunction,
    ci::CI{MOI.SingleVariable},
)
    MOI.throw_if_not_valid(model, ci)
    return MOI.SingleVariable(MOI.VariableIndex(ci.value))
end
function MOI.get(
    model::AbstractModel,
    attr::Union{MOI.ConstraintFunction, MOI.ConstraintSet},
    ci::MOI.ConstraintIndex
)
    return MOI.get(constraints(model, ci), attr, ci)
end

function _get_single_variable_set(
    model::AbstractModel,
    S::Type{<:MOI.EqualTo},
    index,
)
    return MOI.EqualTo(model.lower_bound[index])
end
function _get_single_variable_set(
    model::AbstractModel,
    S::Type{<:Union{MOI.GreaterThan,MOI.EqualTo}},
    index,
)
    # Lower and upper bounds are equal for `EqualTo`, we can take either of them.
    return S(model.lower_bound[index])
end
function _get_single_variable_set(
    model::AbstractModel,
    S::Type{<:MOI.LessThan},
    index,
)
    return S(model.upper_bound[index])
end
function _get_single_variable_set(
    model::AbstractModel,
    S::Type{<:Union{MOI.Interval,MOI.Semicontinuous,MOI.Semiinteger}},
    index,
)
    return S(model.lower_bound[index], model.upper_bound[index])
end
function _get_single_variable_set(
    ::AbstractModel,
    S::Type{<:Union{MOI.Integer,MOI.ZeroOne}},
    index,
)
    return S()
end
function MOI.get(
    model::AbstractModel,
    ::MOI.ConstraintSet,
    ci::CI{MOI.SingleVariable,S},
) where {S}
    MOI.throw_if_not_valid(model, ci)
    return _get_single_variable_set(model, S, ci.value)
end

function MOI.is_empty(model::AbstractModel)
    return isempty(model.name) &&
               !model.senseset &&
               !model.objectiveset &&
               isempty(model.objective.terms) &&
               iszero(model.objective.constant) &&
               iszero(model.num_variables_created) &&
               mapreduce_constraints(MOI.is_empty, &, model, true)
end
function MOI.empty!(model::AbstractModel{T}) where {T}
    model.name = ""
    model.senseset = false
    model.sense = MOI.FEASIBILITY_SENSE
    model.objectiveset = false
    model.objective = zero(MOI.ScalarAffineFunction{T})
    model.num_variables_created = 0
    model.variable_indices = nothing
    model.single_variable_mask = UInt8[]
    model.lower_bound = T[]
    model.upper_bound = T[]
    empty!(model.var_to_name)
    model.name_to_var = nothing
    empty!(model.con_to_name)
    model.name_to_con = nothing
    broadcastcall(MOI.empty!, model)
    return
end

function MOI.copy_to(dest::AbstractModel, src::MOI.ModelLike; kws...)
    return automatic_copy_to(dest, src; kws...)
end
supports_default_copy_to(model::AbstractModel, copy_names::Bool) = true

# Allocate-Load Interface
# Even if the model does not need it and use default_copy_to, it could be used
# by a layer that needs it
supports_allocate_load(model::AbstractModel, copy_names::Bool) = true

function allocate_variables(model::AbstractModel, nvars)
    return MOI.add_variables(model, nvars)
end
allocate(model::AbstractModel, attr...) = MOI.set(model, attr...)
function allocate_constraint(
    model::AbstractModel,
    f::MOI.AbstractFunction,
    s::MOI.AbstractSet,
)
    return MOI.add_constraint(model, f, s)
end

function load_variables(::AbstractModel, nvars) end
function load(::AbstractModel, attr...) end
function load_constraint(
    ::AbstractModel,
    ::CI,
    ::MOI.AbstractFunction,
    ::MOI.AbstractSet,
) end

# Can be used to access constraints of a model
"""
broadcastcall(f::Function, model::AbstractModel)

Calls `f(contrs)` for every vector `constrs::Vector{ConstraintIndex{F, S}, F, S}` of the model.

# Examples

To add all constraints of the model to a solver `solver`, one can do
```julia
_addcon(solver, ci, f, s) = MOI.add_constraint(solver, f, s)
function _addcon(solver, constrs::Vector)
    for constr in constrs
        _addcon(solver, constr...)
    end
end
MOIU.broadcastcall(constrs -> _addcon(solver, constrs), model)
```
"""
function broadcastcall end

"""
broadcastvcat(f::Function, model::AbstractModel)

Calls `f(contrs)` for every vector `constrs::Vector{ConstraintIndex{F, S}, F, S}` of the model and concatenate the results with `vcat` (this is used internally for `ListOfConstraints`).

# Examples

To get the list of all functions:
```julia
_getfun(ci, f, s) = f
_getfun(cindices::Tuple) = _getfun(cindices...)
_getfuns(constrs::Vector) = _getfun.(constrs)
MOIU.broadcastvcat(_getfuns, model)
```
"""
function broadcastvcat end

function mapreduce_constraints end

# Macro to generate Model
abstract type Constraints{F} end

abstract type SymbolFS end
struct SymbolFun <: SymbolFS
    s::Union{Symbol,Expr}
    typed::Bool
    cname::Expr # `esc(scname)` or `esc(vcname)`
end
struct SymbolSet <: SymbolFS
    s::Union{Symbol,Expr}
    typed::Bool
end

# QuoteNode prevents s from being interpolated and keeps it as a symbol
# Expr(:., MOI, s) would be MOI.s
# Expr(:., MOI, $s) would be Expr(:., MOI, EqualTo)
# Expr(:., MOI, :($s)) would be Expr(:., MOI, :EqualTo)
# Expr(:., MOI, :($(QuoteNode(s)))) is Expr(:., MOI, :(:EqualTo)) <- what we want

# (MOI, :Zeros) -> :(MOI.Zeros)
# (:Zeros) -> :(MOI.Zeros)
_set(s::SymbolSet) = esc(s.s)
_fun(s::SymbolFun) = esc(s.s)
function _typedset(s::SymbolSet)
    if s.typed
        :($(_set(s)){T})
    else
        _set(s)
    end
end
function _typedfun(s::SymbolFun)
    if s.typed
        :($(_fun(s)){T})
    else
        _fun(s)
    end
end

# Base.lowercase is moved to Unicode.lowercase in Julia v0.7
using Unicode

_field(s::SymbolFS) = Symbol(replace(lowercase(string(s.s)), "." => "_"))

_getC(s::SymbolSet) = :(VectorOfConstraints{F,$(_typedset(s))})
_getC(s::SymbolFun) = _typedfun(s)

_getCV(s::SymbolSet) = :($(_getC(s))())
_getCV(s::SymbolFun) = :($(s.cname){T,$(_getC(s))}())

_callfield(f, s::SymbolFS) = :($f(model.$(_field(s))))
_broadcastfield(b, s::SymbolFS) = :($b(f, model.$(_field(s))))
_mapreduce_field(s::SymbolFS) = :(cur = $MOIU.mapreduce_constraints(f, op, model.$(_field(s)), cur))
_mapreduce_constraints(s::SymbolFS) = :(cur = op(cur, f(model.$(_field(s)))))

# This macro is for expert/internal use only. Prefer the concrete Model type
# instantiated below.
"""
    macro model(
        model_name,
        scalar_sets,
        typed_scalar_sets,
        vector_sets,
        typed_vector_sets,
        scalar_functions,
        typed_scalar_functions,
        vector_functions,
        typed_vector_functions,
        is_optimizer = false
    )

Creates a type `model_name` implementing the MOI model interface and containing
`scalar_sets` scalar sets `typed_scalar_sets` typed scalar sets, `vector_sets`
vector sets, `typed_vector_sets` typed vector sets, `scalar_functions` scalar
functions, `typed_scalar_functions` typed scalar functions, `vector_functions`
vector functions and `typed_vector_functions` typed vector functions.
To give no set/function, write `()`, to give one set `S`, write `(S,)`.

The function [`MathOptInterface.SingleVariable`](@ref) should not be given in
`scalar_functions`. The model supports [`MathOptInterface.SingleVariable`](@ref)-in-`F`
constraints where `F` is [`MathOptInterface.EqualTo`](@ref),
[`MathOptInterface.GreaterThan`](@ref), [`MathOptInterface.LessThan`](@ref),
[`MathOptInterface.Interval`](@ref), [`MathOptInterface.Integer`](@ref),
[`MathOptInterface.ZeroOne`](@ref), [`MathOptInterface.Semicontinuous`](@ref)
or [`MathOptInterface.Semiinteger`](@ref). The sets supported
with the [`MathOptInterface.SingleVariable`](@ref) cannot be controlled from the
macro, use the [`UniversalFallback`](@ref) to support more sets.

This macro creates a model specialized for specific types of constraint,
by defining specialized structures and methods. To create a model that,
in addition to be optimized for specific constraints, also support arbitrary
constraints and attributes, use [`UniversalFallback`](@ref).

This implementation of the MOI model certifies that the constraint indices, in
addition to being different between constraints `F`-in-`S` for the same types
`F` and `S`, are also different between constraints for different types `F` and
`S`. This means that for constraint indices `ci1`, `ci2` of this model,
`ci1 == ci2` if and only if `ci1.value == ci2.value`. This fact can be used to
use the the value of the index directly in a dictionary representing a mapping
between constraint indices and something else.

If `is_optimizer = true`, the resulting struct is a subtype of
of `MOIU.AbstractOptimizer`, which is a subtype of
[`MathOptInterface.AbstractOptimizer`](@ref), otherwise, it is a subtype of
`MOIU.AbstractModelLike`, which is a subtype of
[`MathOptInterface.ModelLike`](@ref).

### Examples

The model describing an linear program would be:
```julia
@model(LPModel,                                                   # Name of model
      (),                                                         # untyped scalar sets
      (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval), #   typed scalar sets
      (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives),            # untyped vector sets
      (),                                                         #   typed vector sets
      (),                                                         # untyped scalar functions
      (MOI.ScalarAffineFunction,),                                #   typed scalar functions
      (MOI.VectorOfVariables,),                                   # untyped vector functions
      (MOI.VectorAffineFunction,),                                #   typed vector functions
      false
    )
```

Let `MOI` denote `MathOptInterface`, `MOIU` denote `MOI.Utilities` and
`MOIU.ConstraintEntry{F, S}` be defined as `MOI.Tuple{MOI.ConstraintIndex{F, S}, F, S}`.
The macro would create the types:
```julia
struct LPModelScalarConstraints{T, F <: MOI.AbstractScalarFunction} <: MOIU.Constraints{F}
    equalto::Vector{MOIU.ConstraintEntry{F, MOI.EqualTo{T}}}
    greaterthan::Vector{MOIU.ConstraintEntry{F, MOI.GreaterThan{T}}}
    lessthan::Vector{MOIU.ConstraintEntry{F, MOI.LessThan{T}}}
    interval::Vector{MOIU.ConstraintEntry{F, MOI.Interval{T}}}
end
struct LPModelVectorConstraints{T, F <: MOI.AbstractVectorFunction} <: MOIU.Constraints{F}
    zeros::Vector{MOIU.ConstraintEntry{F, MOI.Zeros}}
    nonnegatives::Vector{MOIU.ConstraintEntry{F, MOI.Nonnegatives}}
    nonpositives::Vector{MOIU.ConstraintEntry{F, MOI.Nonpositives}}
end
mutable struct LPModel{T} <: MOIU.AbstractModel{T}
    name::String
    sense::MOI.OptimizationSense
    objective::Union{MOI.SingleVariable, MOI.ScalarAffineFunction{T}, MOI.ScalarQuadraticFunction{T}}
    num_variables_created::Int64
    # If nothing, no variable has been deleted so the indices of the
    # variables are VI.(1:num_variables_created)
    variable_indices::Union{Nothing, Set{MOI.VariableIndex}}
    # Union of flags of `S` such that a `SingleVariable`-in-`S`
    # constraint was added to the model and not deleted yet.
    single_variable_mask::Vector{UInt8}
    # Lower bound set by `SingleVariable`-in-`S` where `S`is
    # `GreaterThan{T}`, `EqualTo{T}` or `Interval{T}`.
    lower_bound::Vector{T}
    # Lower bound set by `SingleVariable`-in-`S` where `S`is
    # `LessThan{T}`, `EqualTo{T}` or `Interval{T}`.
    upper_bound::Vector{T}
    var_to_name::Dict{MOI.VariableIndex, String}
    # If `nothing`, the dictionary hasn't been constructed yet.
    name_to_var::Union{Dict{String, MOI.VariableIndex}, Nothing}
    con_to_name::Dict{MOI.ConstraintIndex, String}
    name_to_con::Union{Dict{String, MOI.ConstraintIndex}, Nothing}
    scalaraffinefunction::LPModelScalarConstraints{T, MOI.ScalarAffineFunction{T}}
    vectorofvariables::LPModelVectorConstraints{T, MOI.VectorOfVariables}
    vectoraffinefunction::LPModelVectorConstraints{T, MOI.VectorAffineFunction{T}}
end
```
The type `LPModel` implements the MathOptInterface API except methods specific
to solver models like `optimize!` or `getattribute` with `VariablePrimal`.
"""
macro model(
    model_name,
    ss,
    sst,
    vs,
    vst,
    sf,
    sft,
    vf,
    vft,
    is_optimizer = false,
)
    scalar_sets = [SymbolSet.(ss.args, false); SymbolSet.(sst.args, true)]
    vector_sets = [SymbolSet.(vs.args, false); SymbolSet.(vst.args, true)]

    scname = esc(Symbol(string(model_name) * "ScalarConstraints"))
    vcname = esc(Symbol(string(model_name) * "VectorConstraints"))

    esc_model_name = esc(model_name)
    header = if is_optimizer
        :($(esc(model_name)){T} <: AbstractOptimizer{T})
    else
        :($(esc(model_name)){T} <: AbstractModelLike{T})
    end

    scalar_funs = [
        SymbolFun.(sf.args, false, Ref(scname))
        SymbolFun.(sft.args, true, Ref(scname))
    ]
    vector_funs = [
        SymbolFun.(vf.args, false, Ref(vcname))
        SymbolFun.(vft.args, true, Ref(vcname))
    ]
    funs = [scalar_funs; vector_funs]

    scalarconstraints = :(
        struct $scname{T,F<:$MOI.AbstractScalarFunction} <: Constraints{F} end
    )
    vectorconstraints = :(
        struct $vcname{T,F<:$MOI.AbstractVectorFunction} <: Constraints{F} end
    )
    for (c, sets) in
        ((scalarconstraints, scalar_sets), (vectorconstraints, vector_sets))
        for s in sets
            field = _field(s)
            push!(c.args[3].args, :($field::$(_getC(s))))
        end
    end

    modeldef = quote
        mutable struct $header
            name::String
            senseset::Bool
            sense::$MOI.OptimizationSense
            objectiveset::Bool
            objective::Union{
                $MOI.SingleVariable,
                $MOI.ScalarAffineFunction{T},
                $MOI.ScalarQuadraticFunction{T},
            }
            num_variables_created::Int64
            # If nothing, no variable has been deleted so the indices of the
            # variables are VI.(1:num_variables_created)
            variable_indices::Union{Nothing,Set{$VI}}
            # Union of flags of `S` such that a `SingleVariable`-in-`S`
            # constraint was added to the model and not deleted yet.
            single_variable_mask::Vector{UInt8}
            # Lower bound set by `SingleVariable`-in-`S` where `S`is
            # `GreaterThan{T}`, `EqualTo{T}` or `Interval{T}`.
            lower_bound::Vector{T}
            # Lower bound set by `SingleVariable`-in-`S` where `S`is
            # `LessThan{T}`, `EqualTo{T}` or `Interval{T}`.
            upper_bound::Vector{T}
            var_to_name::Dict{$VI,String}
            # If `nothing`, the dictionary hasn't been constructed yet.
            name_to_var::Union{Dict{String,$VI},Nothing}
            con_to_name::Dict{$CI,String}
            name_to_con::Union{Dict{String,$CI},Nothing}
            # A useful dictionary for extensions to store things. These are
            # _not_ copied between models!
            ext::Dict{Symbol,Any}
        end
    end
    for f in funs
        cname = f.cname
        field = _field(f)
        push!(modeldef.args[2].args[3].args, :($field::$cname{T,$(_getC(f))}))
    end

    code = quote
        function $MOIU.broadcastcall(f::F, model::$esc_model_name) where {F<:Function}
            return $(Expr(:block, _broadcastfield.(Ref(:(broadcastcall)), funs)...))
        end
        function $MOIU.broadcastvcat(f::F, model::$esc_model_name) where {F<:Function}
            return vcat($(_broadcastfield.(Ref(:(broadcastvcat)), funs)...))
        end
        function $MOIU.mapreduce_constraints(f::Function, op::Function, model::$esc_model_name, cur)
            return $(Expr(:block, _mapreduce_field.(funs)...))
        end
    end
    for (cname, sets) in ((scname, scalar_sets), (vcname, vector_sets))
        code = quote
            $code
            function $MOIU.broadcastcall(f::F, model::$cname) where {F<:Function}
                return $(Expr(:block, _callfield.(:f, sets)...))
            end
            function $MOIU.broadcastvcat(f::F, model::$cname) where {F<:Function}
                return vcat($(_callfield.(:f, sets)...))
            end
            function $MOIU.mapreduce_constraints(f::Function, op::Function, model::$cname, cur)
                return $(Expr(:block, _mapreduce_constraints.(sets)...))
            end
        end
    end

    for (c, sets) in ((scname, scalar_sets), (vcname, vector_sets))
        for s in sets
            set = _set(s)
            field = _field(s)
            code = quote
                $code
                function $MOIU.constraints(
                    model::$c,
                    ::Type{<:$set},
                ) where {F}
                    return model.$field
                end
            end
        end
    end

    for f in funs
        fun = _fun(f)
        field = _field(f)
        code = quote
            $code
            function $MOIU.constraints(
                model::$esc_model_name,
                ::Type{<:$fun},
                ::Type{S}
            ) where S
                return $MOIU.constraints(model.$field, S)
            end
        end
    end

    code = quote
        $scalarconstraints
        function $scname{T,F}() where {T,F}
            return $scname{T,F}($(_getCV.(scalar_sets)...))
        end

        $vectorconstraints
        function $vcname{T,F}() where {T,F}
            return $vcname{T,F}($(_getCV.(vector_sets)...))
        end

        $modeldef
        function $esc_model_name{T}() where {T}
            return $esc_model_name{T}(
                "",
                false,
                $MOI.FEASIBILITY_SENSE,
                false,
                $SAF{T}($MOI.ScalarAffineTerm{T}[], zero(T)),
                0,
                nothing,
                UInt8[],
                T[],
                T[],
                Dict{$VI,String}(),
                nothing,
                Dict{$CI,String}(),
                nothing,
                Dict{Symbol,Any}(),
                $(_getCV.(funs)...),
            )
        end

        function $MOI.supports_constraint(
            model::$esc_model_name{T},
            ::Type{<:Union{$(_typedfun.(scalar_funs)...)}},
            ::Type{<:Union{$(_typedset.(scalar_sets)...)}},
        ) where {T}
            return true
        end
        function $MOI.supports_constraint(
            model::$esc_model_name{T},
            ::Type{<:Union{$(_typedfun.(vector_funs)...)}},
            ::Type{<:Union{$(_typedset.(vector_sets)...)}},
        ) where {T}
            return true
        end

        $code
    end
    return code
end

const LessThanIndicatorSetOne{T} =
    MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE,MOI.LessThan{T}}
const LessThanIndicatorSetZero{T} =
    MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO,MOI.LessThan{T}}

@model(
    Model,
    (MOI.ZeroOne, MOI.Integer),
    (
        MOI.EqualTo,
        MOI.GreaterThan,
        MOI.LessThan,
        MOI.Interval,
        MOI.Semicontinuous,
        MOI.Semiinteger,
    ),
    (
        MOI.Reals,
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.Complements,
        MOI.NormInfinityCone,
        MOI.NormOneCone,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.GeometricMeanCone,
        MOI.ExponentialCone,
        MOI.DualExponentialCone,
        MOI.RelativeEntropyCone,
        MOI.NormSpectralCone,
        MOI.NormNuclearCone,
        MOI.PositiveSemidefiniteConeTriangle,
        MOI.PositiveSemidefiniteConeSquare,
        MOI.RootDetConeTriangle,
        MOI.RootDetConeSquare,
        MOI.LogDetConeTriangle,
        MOI.LogDetConeSquare,
    ),
    (
        MOI.PowerCone,
        MOI.DualPowerCone,
        MOI.SOS1,
        MOI.SOS2,
        LessThanIndicatorSetOne,
        LessThanIndicatorSetZero,
    ),
    (),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction)
)

@doc raw"""

An implementation of `ModelLike` that supports all functions and sets defined
in MOI. It is parameterized by the coefficient type.

# Examples
```jl
model = Model{Float64}()
x = add_variable(model)
```
"""
Model


text:
using YaoBlocks
using LuxurySparse: IMatrix
using BitBasis: readbit,bcat,@bit_str

Locs = Union{Vector{Int}, Int}
"""
Note: according to Figure 2, page 4 of the paper:

    U(nbits::Int, locs::Union{Vector{Int}, Int}, θ::Real, ϕ::Real, λ::Real)
"""
function U(nbits::Int, locs::Locs, θ::Real, ϕ::Real, λ::Real)
    put(nbits, locs=>Rz(ϕ)*Ry(θ)*Rz(λ))
end

"""
Note: according to the section before Figure 2, page 4 of the paper:

    CX(nbits::Int, locs_c::Union{Vector{Int}, Int}, locs_x::Union{Vector{Int}, Int})

CNOT gate, `locs_c` are control qubits locations and `locs_x` are target qubits locations.
"""
function CX(nbits::Int, locs_c::Locs, locs_x::Locs)
    control(nbits, (locs_c...,), (locs_x...,)=>X)
end

"""
    RESET(nbits::Int, locs::Locs)

The reset operation in Open QASM, note here, instead of keeping the remaining bits in mixed ensemble,
we measure the target qubit directly. See `measure_collapseto!` in `Yao.jl` for detail.
"""
function RESET(nbits::Int, locs::Locs)
    Measure(nbits; locs=locs, resetto=bcat(bit"0" for i in 1:nbits))
end

"""
    measure(nbits::Int, locs::Locs, cmem)

Measure qubits in location `locs`, and store results to `cmem`.
Here `cmem` is a classical memory that can be specified like `view(bitarray, [2,3])`.
"""
function MEASURE(nbits::Int, locs::Locs, cmem)
    FlushResults(Measure(nbits; locs=locs), cmem)
end

"""
    Barrier{N} <: PrimitiveBlock{N}

The barrier instruction prevents optimizations from reordering gates across its source
line.
"""
struct Barrier{N} <: PrimitiveBlock{N} end
YaoBlocks.mat(b::Barrier{N}) where N = IMatrix{1<<N}()
YaoBlocks.getiparams(b::Barrier) = ()

include("Bag.jl")
"""
    FlushResults{N, C, MT<:Measure{N,C}} <: AbstractBag{MT,N}

The measure wrapper that flushs the output in a `Measure` block to classical memory.
"""
struct FlushResults{N,C,TT, MT<:Measure{N,C}} <: AbstractBag{MT,N}
    content::MT
    target_space::TT
    function FlushResults(mblock::MT, target_space::TT) where {N,C,TT,MT<:Measure{N,C}}
        length(target_space) != C && throw(QubitMismatchError("Expect $C, got target space size $(length(target_space))"))
        new{N,C,TT,MT}(mblock, target_space)
    end
end

function YaoBlocks.apply!(reg::AbstractRegister, fl::FlushResults{N,C}) where {N, C}
    throw(NotImplementedError())
end
function YaoBlocks.apply!(reg::AbstractRegister{1}, fl::FlushResults{N,C}) where {N, C}
    out = apply!(reg, fl.content)
    res = fl.content.results[]
    for i=1:C
        fl.target_space[i] = readbit(res, i)
    end
    return out
end

function YaoBlocks.print_annotation(io::IO, fl::FlushResults)
    printstyled(io, "$(fl.target_space) ← "; bold=true, color=:cyan)
end

function Base.show(io::IO, ::MIME"plain/text", blk::FlushResults)
    return print_tree(io, blk; title=false, compact=false)
end

ext: jl


text:
###
# Block
###

Base.in(K::Block, B::BroadcastVector{<:Block,Type{Block}}) = Int(K) in B.args[1]

###
# BlockBanded
###

bandedblockbandedbroadcaststyle(::LazyArrayStyle{2}) = LazyArrayStyle{2}()
bandedblockbandedcolumns(::LazyLayout) = BandedBlockBandedColumns{LazyLayout}()
bandedblockbandedcolumns(::ApplyLayout) = BandedBlockBandedColumns{LazyLayout}()
bandedblockbandedcolumns(::BroadcastLayout) = BandedBlockBandedColumns{LazyLayout}()

isblockbanded(K::BlockKron) = isbanded(first(K.args))
isbandedblockbanded(K::BlockKron) = all(isbanded, K.args)
blockbandwidths(K::BlockKron) = bandwidths(first(K.args))
subblockbandwidths(K::BlockKron) = bandwidths(last(K.args))

"""
    DiagTrav(A::AbstractMatrix)

converts a matrix to a block vector by traversing the anti-diagonals.
"""
struct DiagTrav{T, N, AA<:AbstractArray{T,N}} <: AbstractBlockVector{T}
    array::AA
end

function axes(A::DiagTrav{<:Any,2})
    m,n = size(A.array)
    mn = min(m,n)
    (blockedrange(Vcat(OneTo(mn), Fill(mn,max(m,n)-mn))),)
end

function axes(A::DiagTrav{<:Any,3})
    m,n,p = size(A.array)
    @assert m == n == p
    (blockedrange(cumsum(OneTo(m))),)
end

function getindex(A::DiagTrav{<:Any,2}, K::Block{1})
    k = Int(K)
    m,n = size(A.array)
    mn = min(m,n)
    st = stride(A.array,2)
    if k ≤ m
        A.array[range(k; step=st-1, length=min(k,mn))]
    else
        A.array[range(m+(k-m)*st; step=st-1, length=min(k,mn))]
    end
end

function getindex(A::DiagTrav{T,3}, K::Block{1}) where T
    k = Int(K)
    m,n,p = size(A.array)
    @assert m == n == p
    st = stride(A.array,2)
    st3 = stride(A.array,3)
    ret = A.array[range(k; step=st-1, length=k)]
    for j = 1:k-1
        append!(ret, view(A.array, range(j*st3 + (k-j); step=st-1, length=k-j)))
    end
    ret
end

getindex(A::DiagTrav, k::Int) = A[findblockindex(axes(A,1), k)]

function resize!(A::DiagTrav{<:Any,2}, K::Block{1})
    k = Int(K)
    DiagTrav(A.array[1:k, 1:k])
end

struct InvDiagTrav{T, AA<:AbstractVector{T}} <: AbstractMatrix{T}
    vector::AA
end

size(A::InvDiagTrav) = (blocksize(A.vector,1),blocksize(A.vector,1))

function getindex(A::InvDiagTrav{T}, k::Int, j::Int)  where T
    if k+j-1 ≤ blocksize(A.vector,1)
        A.vector[Block(k+j-1)][j]
    else
        zero(T)
    end
end

struct KronTrav{T, N, AA<:AbstractArray{T,N}, BB<:AbstractArray{T,N}, AXES} <: AbstractBlockArray{T, N}
    A::AA
    B::BB
    axes::AXES
end

KronTrav(A::AbstractArray{T,N}, B::AbstractArray{V,N}, axes) where {T,V,N} =
    KronTrav{promote_type(T,V),N,typeof(A),typeof(B),typeof(axes)}(A, B, axes)

KronTrav(A::AbstractArray{T,N}, B::AbstractArray{V,N}) where {T,V,N} =
    KronTrav(A, B, _krontrav_axes(axes(A), axes(B)))

function _krontrav_axes(A::NTuple{N,OneTo{Int}}, B::NTuple{N,OneTo{Int}}) where N
    m,n = length.(A), length.(B)
    mn = min.(m,n)
    @. blockedrange(Vcat(OneTo(mn), Fill(mn,max(m,n)-mn)))
end

axes(A::KronTrav) = A.axes

function getindex(A::KronTrav{<:Any,1}, K::Block{1})
    m,n = length(A.A), length(A.B)
    mn = min(m,n)
    k = Int(K)
    if k ≤ mn
        A.A[1:k] .* A.B[k:-1:1]
    elseif m < n
        A.A .* A.B[k:-1:(k-m+1)]
    else # n < m
        A.A[(k-n+1):k] .* A.B[end:-1:1]
    end
end

function getindex(A::KronTrav{<:Any,2}, K::Block{2})
    m,n = size(A.A), size(A.B)
    @assert m == n
    k,j = K.n
    # Following is equivalent to A.A[1:k,1:j] .* A.B[k:-1:1,j:-1:1]
    # layout_getindex to avoid SparseArrays from Diagonal
    # rot180 to preserve bandedness
    layout_getindex(A.A,1:k,1:j) .* rot180(layout_getindex(A.B,1:k,1:j))
end
getindex(A::KronTrav{<:Any,N}, kj::Vararg{Int,N}) where N =
    A[findblockindex.(axes(A), kj)...]

# A.A[1:k,1:j] has A_l,A_u
# A.B[k:-1:1,j:-1:1] has bandwidths (B_u + k-j, B_l + j-k)
subblockbandwidths(A::KronTrav) = bandwidths(A.A)
blockbandwidths(A::KronTrav) = bandwidths(A.A) .+ bandwidths(A.B)
isblockbanded(A::KronTrav) = isbanded(A.A) && isbanded(A.B)
isbandedblockbanded(A::KronTrav) = isbanded(A.A) && isbanded(A.B)

struct KronTravBandedBlockBandedLayout <: AbstractBandedBlockBandedLayout end

krontravlayout(_, _) = UnknownLayout()
krontravlayout(::AbstractBandedLayout, ::AbstractBandedLayout) = KronTravBandedBlockBandedLayout()
MemoryLayout(::Type{KronTrav{T,N,AA,BB,AXIS}}) where {T,N,AA,BB,AXIS} = krontravlayout(MemoryLayout(AA), MemoryLayout(BB))

sublayout(::KronTravBandedBlockBandedLayout, ::Type{<:NTuple{2,BlockSlice1}}) = BroadcastBandedLayout{typeof(*)}()

call(b::BroadcastLayout{typeof(*)}, a::KronTrav) = *
call(b::BroadcastBandedLayout{typeof(*)}, a::SubArray) = *

function _broadcast_sub_arguments(::KronTravBandedBlockBandedLayout, A, V)
    K,J = parentindices(V)
    k,j = Int(K.block),Int(J.block)
    view(A.A,1:k,1:j), ApplyMatrix(rot180,view(A.B,1:k,1:j))
end

krontavbroadcaststyle(::BandedStyle, ::BandedStyle) = BandedBlockBandedStyle()
krontavbroadcaststyle(::BandedStyle, ::StructuredMatrixStyle{<:Diagonal}) = BandedBlockBandedStyle()
krontavbroadcaststyle(::StructuredMatrixStyle{<:Diagonal}, ::BandedStyle) = BandedBlockBandedStyle()
krontavbroadcaststyle(::LazyArrayStyle{2}, ::BandedStyle) = LazyArrayStyle{2}()
krontavbroadcaststyle(::BandedStyle, ::LazyArrayStyle{2}) = LazyArrayStyle{2}()
krontavbroadcaststyle(::LazyArrayStyle{2}, ::LazyArrayStyle{2}) = LazyArrayStyle{2}()
BroadcastStyle(::Type{KronTrav{T,N,AA,BB,AXIS}}) where {T,N,AA,BB,AXIS} =
    krontavbroadcaststyle(BroadcastStyle(AA), BroadcastStyle(BB))


text:
export Analysis

abstract type Options end

"""
    Analysis
Set-up analysis parameters.
### Fields:
* `filename`    - unique filename of the model.
* `analyse`     - R, C, I, G, B stand for ray, coheren, incoherence, geometri and gaussian
* `num_ray`     - number of rays.
* `alpha`       - ray lauch angles between -90 and 90 deg.
* `box`         - analysis box x- range, km and y-depth, m.
### Available Constructors:
* `Source(;frequency::Real = 1000.0f0, height::Real = 80.0f0))`
* `Source(frequency, height)`
"""
struct Analysis{T<:AbstractArray} <: Options
    filename::String
    analyse::String
    option1::String
    option2::String
    num_ray::Int
    alpha::Vec2{T}
    box::Vec2{T}
    step::Real
end

Analysis(;filename = "Model_f1000Hz",
            analyse = "CB",
            option1 = "SVW",
            option2 = "A",
            num_ray = 80,
            alpha = Vec2(-80.0f0,80.0f0),
            box = Vec2(10f0,1000.0f0),
            step= 0) =
    Analysis(filename, analyse, option1, option2, num_ray, alpha, box, step)


text:
# import TopicModelHealthWebsites
include("/publichealth-websites/TopicModelHealthWebsites/src/TopicModelHealthWebsites.jl") # update path
TopicModelHealthWebsites.HALT_LDA(K = 10, datapath = "/data/dataanddict.jld", savepath = "/data/test.jld") # update path

max_stars_repo_name: jwanghb/publichealth-websites

max_issues_repo_name: jwanghb/publichealth-websites

max_forks_repo_name: jwanghb/publichealth-websites


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/PkgTestSuite.jl/blob/master/LICENSE

using PkgTestSuite: checktabs
using Base.Test

new_file_path = Pkg.dir("PkgTestSuite", "src", "new_file_with_tabs.jl")
isfile(new_file_path) && rm(new_file_path)

@test checktabs("PkgTestSuite")

test_file_contents = """
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/PkgTestSuite.jl/blob/master/LICENSE

function foo()
\treturn 1
end
"""
open(new_file_path, "w") do fid
    write(fid, test_file_contents)
end
@test_throws Exception checktabs("PkgTestSuite")
isfile(new_file_path) && rm(new_file_path)


text:
module ConwayPolyhedra

include("operators.jl")

end


text:
module WebApp

using HTTP, Sockets
using ..Gameplay, ..GameSession, ..Wikipedia, ..Wikipedia.Articles

# Configuration
const HOST = ip"0.0.0.0"
const PORT = 8888
const ROUTER = HTTP.Router()
const SERVER = HTTP.Server(ROUTER)

# Functions
function wikiarticle(game, article)
  html = """
  <!DOCTYPE html>

  $(head())

    $(objective(game))

    $(
      if losinggame(game)
        "
You Lost :(
"
      else
        puzzlesolved(game, article) ? "
You Won!
" : ""
      end
    )

$(article.title)

      $(replace(article.content, "/wiki/"=>"/$(game.id)/wiki/"))

  """
end

function history(game)
  html = """
"""
  iter = 0
  for a in game.history
    html *= """

        $(a.title)

    """
    iter += 1
  end

  html * "
"
end

function objective(game)
  """

Go from
      $(game.articles[1].title)
      to
      $(game.articles[end].title)

      Progress:
      $(size(game.history, 1) - 1)
      out of maximum
      $(size(game.articles, 1) - 1)
      links in
      $(game.steps_taken)
      steps

    $(history(game))

      Solution? |
      New game

  """
end

function head()
  """

    <meta charset="utf-8" />
    <link rel="stylesheet" href="https://stackpath.bootstrapcdn.com/bootstrap/4.1.3/css/bootstrap.min.css" integrity="sha384-MCw98/SFnGE8fJT3GXwEOngsV7Zt27NXFoaoApmYm81iuXoPkFOJwJ8ERdknLPMO" crossorigin="anonymous">
    <title>6 Degrees of Wikipedia</title>

  """
end

function puzzlesolved(game, article)
  article.url == game.articles[end].url
end

function losinggame(game)
  game.steps_taken >= Gameplay.MAX_NUMBER_OF_STEPS
end

function parseuri(uri)
  map(x -> String(x), split(uri, "/", keepempty = false))
end

# Routes handlers
const landingpage = HTTP.HandlerFunction() do req
  html = """
  <!DOCTYPE html>

  $(head())

Six degrees of Wikipedia

        The goal of the game is to find the shortest path between two random Wikipedia articles.

        Depending on the difficulty level you choose, the Wiki pages will be further apart and less related.

        If you can’t find the solution, you can always go back up the articles chain, but you need to find the solution within the maximum number of steps, otherwise you lose.

        If you get stuck, you can always check the solution, but you'll lose.

        Good luck and enjoy!

New game

          Easy ($(Gameplay.DIFFICULTY_EASY) links away) |
          Medium ($(Gameplay.DIFFICULTY_MEDIUM) links away) |
          Hard ($(Gameplay.DIFFICULTY_HARD) links away)

  """

  HTTP.Messages.Response(200, html)
end

const newgamepage = HTTP.HandlerFunction() do req
  game = parse(UInt8, (replace(req.target, "/new/"=>""))) |> newgamesession
  article = game.articles[1]
  push!(game.history, article)

  HTTP.Messages.Response(200, wikiarticle(game, article))
end

const articlepage = HTTP.HandlerFunction() do req
  uri_parts = parseuri(req.target)
  game = gamesession(uri_parts[1])
  article_uri = "/wiki/$(uri_parts[end])"

  existing_articles = Articles.find(article_uri)
  article = isempty(existing_articles) ? persistedarticle(fetchpage(article_uri)...) : existing_articles[1]

  push!(game.history, article)
  game.steps_taken += 1

  puzzlesolved(game, article) && destroygamesession(game.id)

  HTTP.Messages.Response(200, wikiarticle(game, article))
end

const backpage = HTTP.HandlerFunction() do req
  uri_parts = parseuri(req.target)
  game = gamesession(uri_parts[1])
  history_index = parse(UInt8, uri_parts[end])

  article = game.history[history_index]
  game.history = game.history[1:history_index]

  HTTP.Messages.Response(200, wikiarticle(game, article))
end

const solutionpage = HTTP.HandlerFunction() do req
  uri_parts = parseuri(req.target)
  game = gamesession(uri_parts[1])
  game.history = game.articles
  game.steps_taken = Gameplay.MAX_NUMBER_OF_STEPS
  article = game.articles[end]

  HTTP.Messages.Response(200, wikiarticle(game, article))
end

const notfoundpage = HTTP.HandlerFunction() do req
  HTTP.Messages.Response(404, "Sorry, this can't be found")
end

# Routes definitions
HTTP.register!(ROUTER, "/", landingpage) # root page
HTTP.register!(ROUTER, "/new/*", newgamepage) # /new/$difficulty_level -- new game
HTTP.register!(ROUTER, "/*/wiki/*", articlepage) # /$session_id/wiki/$wikipedia_article_url -- article page
HTTP.register!(ROUTER, "/*/back/*", backpage) # /$session_id/back/$number_of_steps -- go back the navigation history
HTTP.register!(ROUTER, "/*/solution", solutionpage) # /$session_id/solution -- display the solution
HTTP.register!(ROUTER, "*", notfoundpage) # everything else -- not found

# Start server
HTTP.serve(SERVER, HOST, PORT)

end

max_stars_repo_name: PacktPublishing/Julia-Projects

max_issues_repo_name: PacktPublishing/Julia-Projects

max_forks_repo_name: PacktPublishing/Julia-Projects


text:
### Fetch required packages and reaction networks ###
using Catalyst, Test
using ModelingToolkit: get_ps, get_states, get_eqs, get_systems, get_iv
include("test_networks.jl")

using StableRNGs
rng = StableRNG(12345)

function unpacksys(sys)
    get_eqs(sys),get_iv(sys),get_ps(sys),nameof(sys),get_systems(sys)
end

### Tests construction of empty reaction networks ###
empty_network_1 = @reaction_network
eqs,iv,ps,name,systems = unpacksys(empty_network_1)
@test length(eqs) == 0
@test nameof(iv) == :t
@test length(get_states(empty_network_1)) == 0
@test length(ps) == 0

empty_network_2 = @reaction_network
@variables p1 p2 p3 p4 p5
addparam!(empty_network_2, p1)
addparam!(empty_network_2, p2)
addparam!(empty_network_2, p3)
addparam!(empty_network_2, p4)
addparam!(empty_network_2, p5)
eqs,iv,ps,name,systems = unpacksys(empty_network_2)
@test length(eqs) == 0
@test nameof(iv) == :t
@test length(get_states(empty_network_2)) == 0
@test length(ps) == 5
@test all(getproperty.(ps,:name) .== [:p1,:p2,:p3,:p4,:p5])

### Tests accessing parameters and species added with network API ###
empty_network_3 = @reaction_network
@variables x p
addspecies!(empty_network_3, x)
addparam!(empty_network_3, p)
@test isequal(empty_network_3.x, states(empty_network_3, x))
@test isequal(empty_network_3.p, parameters(empty_network_3, p))

### Tests creating a network and adding reactions ###
unfinished_network = @reaction_network begin
    (k1,k2), X1 ↔ X2
    (k3,k4), X3 ↔ X4
end k0 k1 k2 k3 k4
@add_reactions unfinished_network begin
    (k3,k4), X3 ↔ X4
    (k5,k6), X5 ↔ X6
end k0 k3 k4 k5 k6
@add_reactions unfinished_network begin
    (k5,k6), X5 ↔ X6
    (k7,k8), X7 ↔ X8
end k0 k5 k6 k7 k8
@test length(get_states(unfinished_network)) == 8
@test length(get_ps(unfinished_network)) == 9

### Compares test network to identical network constructed via @add_reactions ###
identical_networks = Vector{Pair}()

step_by_step_network_1 = @reaction_network rns5 begin
    p, ∅ → X1
    (k1,k2), X1 ⟷ X2
    (k3,k4), X2 ⟷ X3
end p k1 k2 k3 k4
@add_reactions step_by_step_network_1 begin
    (k5,k6), X3 ⟷ X4
end k5 k6
@add_reactions step_by_step_network_1 begin
    d, X4 → ∅
end d
push!(identical_networks, reaction_networks_standard[5] => step_by_step_network_1)

step_by_step_network_2 = @reaction_network rns7 begin
    (p1,p2,p3), ∅ → (X1,X2,X3)
end p1 p2 p3
@add_reactions step_by_step_network_2 begin
    (k1,k2), X1 + 2X2 ⟷ X4
    (mm(X3,v1,K1),k3), X4 ⟷ X5
    (d1,d2,d3,d4,d5), (X1,X2,X3,X4,X5) → ∅
end k1 k2 k3 v1 K1 d1 d2 d3 d4 d5
push!(identical_networks, reaction_networks_standard[7] => step_by_step_network_2)

step_by_step_network_3 = @reaction_network rns10 begin
    p, ∅ ⟶ X1
    (k1,k2), X1 → X2
end p k1 k2 k3 k4
@add_reactions step_by_step_network_3 begin
    (k3,k4), X2 → X3
    (k5,k6), X3 → X4
end p k1 k2 k3 k4 k5 k6
@add_reactions step_by_step_network_3 begin
    (k7,k8), X4 → X5
    d, X5  ⟶ ∅
end p k1 k2 k3 k4 k5 k6 k7 k8 d
push!(identical_networks, reaction_networks_standard[10] => step_by_step_network_3)

step_by_step_network_4 = @reaction_network rnh7 begin
    v/10 + hill(X1,v,K,n), ∅ → X1 + X2
end v K n k1 k2 k3 d
@add_reactions step_by_step_network_4 begin
    k1, X1 + X2 → X3
    k2, X3 → X1 + X2
end k1 k2
@add_reactions step_by_step_network_4 begin
    k3, X3 → X1
    d, (X1,X2,X3) → ∅
end k2 k3 d
push!(identical_networks, reaction_networks_hill[7] => step_by_step_network_4)

step_by_step_network_5 = @reaction_network rnc1 begin
    (k1,k2), X1 ↔ X2
    (k3,k4), X2 ↔ X3
end k1 k2 k3 k4 k5 k6
@add_reactions step_by_step_network_5 begin
    k5, X3 → X1
end k5
@add_reactions step_by_step_network_5 begin
    k6, X1 → X3
end k6
push!(identical_networks, reaction_networks_constraint[1] => step_by_step_network_5)

step_by_step_network_6 = @reaction_network rnc5 begin
    (k1,k2), X1 ↔ 2X2
end k1 k2 k3 k4 k5 k6
@add_reactions step_by_step_network_6 begin
    (k3,k4), 2X2 ↔ 3X3
end k3 k4
@add_reactions step_by_step_network_6 begin
    (k5,k6), 3X3 ↔ 4X4
end k5 k6
push!(identical_networks, reaction_networks_constraint[5] => step_by_step_network_6)

step_by_step_network_7 = @reaction_network rnr3 begin
    k2p, Y → 0
end k1 k2p k2pp k3p k3pp A J3 k4 m J4
@add_reactions step_by_step_network_7 begin
    (k3p+k3pp*A)/(J3+Po), Po → P
    k2pp*P, Y → 0
end k3p k3pp A J3 k2pp
@add_reactions step_by_step_network_7 begin
    k1, 0 → Y
end k1
@add_reactions step_by_step_network_7 begin
    (k4*m)/(J4+P), Y + P → Y + Po
end k4 m J4
push!(identical_networks, reaction_networks_real[3] => step_by_step_network_7)

step_by_step_network_8 = @reaction_network rnw7
@parameters k1
addparam!(step_by_step_network_8,k1)
@add_reactions step_by_step_network_8 begin
    k1, X1 → X2
    0, X2 → X3
end k1
@add_reactions step_by_step_network_8 begin
    k2, X3 → X4
    k3, X4 → X5
end k2 k3
push!(identical_networks, reaction_networks_weird[7] => step_by_step_network_8)

step_by_step_network_9 = @reaction_network rnw10
@add_reactions step_by_step_network_9 begin
    d, 5X1 → 4X1
end d
push!(identical_networks, reaction_networks_weird[10] => step_by_step_network_9)

function permute_ps(pvals, rn1, rn2)
    ps1    = parameters(rn1)
    ps2    = parameters(rn2)
    pvals2 = similar(pvals)
    for (i,p) in enumerate(ps2)
        pidx      = findfirst(isequal(p), ps1)
        pvals2[i] = pvals[pidx]
    end
    pvals2
end

for networks in identical_networks
    f1 = ODEFunction(convert(ODESystem,networks[1]),jac=true)
    f2 = ODEFunction(convert(ODESystem,networks[2]),jac=true)
    g1 = SDEFunction(convert(SDESystem,networks[1]))
    g2 = SDEFunction(convert(SDESystem,networks[2]))
    @test networks[1] == networks[2]
    for factor in [1e-2, 1e-1, 1e0, 1e1]
        u0 = factor*rand(rng,length(get_states(networks[1])))
        pp = factor*rand(rng,length(get_ps(networks[1])))

        # needed because this code assumes an ordering of the parameters and species...
        pp2 = permute_ps(pp, networks[1], networks[2])
        t = rand(rng)
        @test all(abs.(f1(u0,pp,t) .- f2(u0,pp2,t)) .< 1000*eps())
        @test all(abs.(f1.jac(u0,pp,t) .- f2.jac(u0,pp2,t)) .< 1000*eps())
        @test all(abs.(g1(u0,pp,t) .- g2(u0,pp2,t)) .< 1000*eps())
    end
end


text:
function julia_info()
    return (
        version = string(VERSION),
        git = (
            commit = Base.GIT_VERSION_INFO.commit,
            branch = Base.GIT_VERSION_INFO.branch,
        ),
        is_debugbuild = ccall(:jl_is_debugbuild, Cint, ()) != 0,
        libllvm_version = string(Base.libllvm_version),
        Sys = (
            WORD_SIZE = Sys.WORD_SIZE,
            JIT = Sys.JIT,
            # CPU_NAME = Sys.CPU_NAME,
            # CPU_THREADS = Sys.CPU_THREADS,
        ),
    )
end

const JuliaInfo = typeof(julia_info())
const JuliaInfoWithCommand = typeof((; command = String[], julia_info()...))

function julia_info(julia::JuliaConfig)
    # TODO: Don't assume project-compatibility
    code = """
    $(Base.load_path_setup_code())

    using BenchmarkConfigSweeps
    using BenchmarkConfigSweeps.Internal: JSON3, julia_info
    JSON3.write(stdout, julia_info())
    """

    cmd = Cmd(julia.cmd)
    cmd = `$cmd --startup-file=no -e 'include_string(Main, read(stdin, String))'`
    outio = communicate(cmd; input = code).stdout
    output = String(take!(outio))
    return JSON3.read(output, JuliaInfo)
end

function asinfo(configs::Vector{Config})
    julia_configs =
        unique!(JuliaConfig[cfg.julia for cfg in configs if cfg.julia !== nothing])
    if isempty(julia_configs)
        julia = [(; command = String[], julia_info()...)]
    else
        julia = [(; command = jl.cmd, julia_info(jl)...) for jl in julia_configs]
    end
    return (
        BenchmarkConfigSweeps = (
            version = "0.0.1",  # TODO: load it from Project.toml
        ),
        julia = julia::Vector{JuliaInfoWithCommand},
        configs = configs,
        # TODO: add more metadata? begin/end time?
    )
end

const Info = typeof(asinfo(Config[]))

function dumpinfo(path, configs::Vector{Config})
    info = asinfo(configs)
    JSON3.write(path, info)
end

loadinfo(path) = JSON3.read(read(path, String), Info)


text:
using DataFrames, TimeZones, Statistics

const DESIGN_COLS = Dict("speed" => 2,
                         "climb" => 3)

struct PreDesign
    hr::Vector{Float64}
    time::Vector{Float64}
    date::Float64
    X::Matrix{Float64} # [speed climb]

    function PreDesign(rs::RunSummary)
        return new(rs.hr, rs.time, decimaldate(rs.start_time),
                   [rs.dist ./ rs.time, rs.alt ./ rs.time])
    end
end

function covariate_dimension()
    return 2
end

# function speed_col()
#     return
function exponential_decay(rs::RunSummary, raw::AbstractArray{Float64, 1},
                           α::Float64;
                           initial_value::Float64 = 0.0,
                           x::AbstractArray{Float64, 1} = zeros(length(rs)))
    n = length(rs)
    @assert length(x) == n

    val = initial_value

    for i = 1:n
        val = exp(-α * rs.time[i]) * val + raw[i]
        x[i] = val
    end

    return x
end

function sliding_scale(rs::RunSummary, raw::AbstractArray{Float64, 1},
                        window::Float64;
                        initial_value::Float64 = 0.0,
                        x::AbstractArray{Float64, 1} = zeros(length(rs)))

    # @assert rs.unit_time

    # t_unit = rs.time[1]
    # m = Int(round(window / t_unit))
    # new_window = m * t_unit

    n = length(rs)
    @assert length(x) == n

    val = initial_value

    previous_raw = 0.0
    to_last = window
    last_value = initial_value
    last_ind = 0

    for i = 1:n
        if rs.time[i] > window
            val = window * raw[i]
            x[i] = val
        else
            t_pre = rs.time[i]
            t_post = rs.time[i]
            while to_last < t_pre
                val -= to_last * last_value

                last_ind += 1
                last_value = raw[last_ind]
                t_pre -= to_last
                to_last = rs.time[last_ind]
            end

            to_last -= t_pre
            val -= t_pre * last_value
            val += rs.time[i] * raw[i]
        end
        x[i] = val
    end

    x ./= window

    return x
end

function design_mat(run_sums::Array{RunSummary}, α::Float64;
                    standardize::Bool = false,
                    decay_function::Function = exponential_decay)
    n_runs = length(run_sums)
    ns = [length(rs) for rs in run_sums]
    n = sum(ns)

    n_covars = covariate_dimension()

    X = zeros(n, n_covars + n_runs + 1)
    y = zeros(n)

    offset = 0
    for i = 1:n_runs
        rs = run_sums[i]
        speed = rs.dist ./ rs.time # TODO: buffer
        climb = rs.alt ./ rs.time # TODO: buffer
        x = findall(isnan, speed)
        # @show rs.dist[x]
        # @show rs.time[x]
        speed[x] .= 0
        climb[x] .= 0
        # @show findall(isnan, speed)
        # @show findall(isnan, climb)

        run_rng = (offset + 1):(offset + ns[i])

        speed_decay = @view X[run_rng, 2]
        climb_decay = @view X[run_rng, 3]

        decay_function(rs, speed, α, x = speed_decay)
        decay_function(rs, climb, α, x = climb_decay)

        date_ind = i + n_covars + 1
        X[run_rng, date_ind] .= 1.0
        y[run_rng] .= rs.hr

        t = 0.0

        for j = 1:ns[i]
            t += rs.time[j]
            X[offset + j, 1] = t
        end

        offset += ns[i]
    end

    if standardize
        for j = 2:3
            μ = 0.0
            σ2 = 0.0
            for i = 1:n
                μ += X[i, j]
                σ2 += X[i, j] * X[i, j]
            end
            μ = μ / n
            σ2 = σ2 / n - μ * μ

            X[:, j] .-= μ
            X[:, j] ./= sqrt(σ2)
        end
    end
    return y, X
end


text:
export
    d_to_dms,
    diff_angle,
    dms_to_d,
    fday_to_hms,
    radianstime_to_fday,
    hangle_to_dec_deg_alt,
    hangle_to_dec_deg,
    hms_to_fday,
    interpolate_angle3,
    interpolate3,
    @polynomial_horner,
    polynomial,
    quadratic_interpolation,
    quadratic_roots,
    set_latitude,
    set_longitude,
    set_obliquity

function set_latitude(lat)
     global latitude = lat
end

function set_longitude(long)
     global longitude = long
end

function set_obliquity(ob)
     global currentobliquity = ob
end

# Meeus chapter 4, quadratic curve fitting:

function quadratic_interpolation(x_list, y_list)
    (P, Q, R, S, T, U, V) = 0, 0, 0, 0, 0, 0, 0
    for (i, v) in enumerate(x_list)
    	v2 = v * v
    	v3 = v2 * v
    	v4 = v3 * v
    	P = P + v
    	Q = Q + v2
    	R = R + v3
    	S = S + v4
    	T = T + y_list[i]
    	U = U + v * y_list[i]
    	V = V + v2 * y_list[i]
    end
    N = length(x_list)
    D =  N * Q * S + 2 * P * Q * R - Q * Q * Q - P * P * S - N * R * R
    a = (N * Q * V + P * R * T + P * Q * U - Q * Q * T - P * P * V - N * R * U) / D
    b = (N * S * U + P * Q * V + Q * R * T - Q * Q * U - P * S * T - N * R * V) / D
    c = (Q * S * T + Q * R * U + P * R * V - Q * Q * V - P * S * U - R * R * T) / D
    return (a, b, c)
end

function quadratic_roots(a, b, c)
    delta = b * b - 4 * a * c
    if delta < 0
        return (false, false)
    end
    if delta == 0
        return (-b/2/a, -b/2/a)
    end
    sqd = sqrt(delta)
    return ((-b + sqd)/2/a, (-b - sqd)/2/a)
end

"""
    d_to_dms(x)

Convert an angle in decimal degrees to degree components.

Return a tuple (degrees, minutes, seconds). Degrees and minutes
will be integers, seconds may be floating.

If the argument is negative:
    The return value of degrees will be negative.
    If degrees is 0, minutes will be negative.
    If minutes is 0, seconds will be negative.

Parameters:
    x : degrees

Returns:
    degrees
    minutes
    seconds
"""
function d_to_dms(x)
    if x < 0
        negative = true
    else
        negative = false
    end
    x = abs(x)
    deg =  floor(Integer, x)
    x = x-deg
    mn =  floor(Integer, x * 60)
    x = x - mn / 60.0
    second = x * 3600
    if negative
        if deg > 0
            deg = -deg
        elseif mn > 0
            mn = -mn
        else
            second = -second
        end
    end
    return (deg, mn, second)
end

"""
    diff_angle(a, b)

Return angle b - a, accounting for circular values.

Parameters a and b should be in the range 0..pi*2. The
result will be in the range -pi..pi.

This allows us to directly compare angles which cross through 0:

    359 degress... 0 degrees... 1 degree... etc

Parameters:
    a : first angle, in radians
    b : second angle, in radians

Returns:
    b - a, in radians
"""
function diff_angle(a, b)
    if b < a
      result = (b + (2 * pi)) - a
    else
      result = b - a
    end
    return mod2pi(result)
end

"""
    dms_to_d(deg, minute, second)

Convert an angle in degree components to decimal degrees.

If any of the components are negative the result will also be negative.

Parameters:
    deg : degrees
    min : minutes
    sec : seconds

Returns:
    decimal degrees
"""
function dms_to_d(deg, minute, second)
    result = abs(deg) + (abs(minute) / 60) + (abs(second) / 3600)
    if deg < 0 || minute < 0 || second < 0
      result = -result
    end
    return result
end

"""
    hangle_to_dec_deg(hour, minute, second)

Convert an hour angle in hour, minute, seconds to decimal degrees.

?if any of the components are negative the result will also be negative.?

Parameters:
    hour : hours
    min : minutes
    sec : seconds

Returns:
    decimal degrees

"""
function hangle_to_dec_deg(hour, minute, second)
    t = dms_to_d(hour, minute, second)
    return mod((t * 15), 360)
end

# function hangle_to_dec_deg(hour, minute, second)
#     global seconds_per_day
#     s = abs(hour) * 3600 + abs(minute) * 60 + abs(second)
#     s = mod(s, seconds_per_day)
#     if hour < 0 || minute < 0 || second < 0
#     	s = -s
#     end
#     return mod(s, seconds_per_day)
# end

"""
    fday_to_hms(day)

Convert fractional day (0.0..1.0) to integral hours, minutes, seconds.

Parameters:
    day : a fractional day in the range 0.0..1.0

Returns:
    hour : 0..23
    minute : 0..59
    second : 0..59
"""
function fday_to_hms(day)
    global seconds_per_day
    tsec = day * seconds_per_day
    tmin = tsec / 60
    thour = tmin / 60
    hour = thour % 24
    minutes = tmin % 60
    seconds = tsec % 60
    return (floor(Integer, hour),  floor(Integer, minutes), seconds)
end

"""
    hms_to_fday(hr, mn, seconds)

Convert hours-minutes-seconds into a fractional day 0.0..1.0.

Parameters:
    hr : hours, 0..23
    mn : minutes, 0..59
    sec : seconds, 0..59

Returns:
    fractional day, 0.0..1.0
"""
function hms_to_fday(hr, mn, seconds)
    global minutes_per_day, seconds_per_day
    return ((hr / 24.0) + (mn / minutes_per_day) + (seconds / seconds_per_day))
end

"""
    radianstime_to_fday(tr)

Convert a time of day in radians to fractional day 0 through 1.
"""
function radianstime_to_fday(tr)
    return (tr / (2 * pi))
end

"""
    interpolate3(n, y)

Interpolate from three equally spaced tabular values.

[Meeus-1998 equation 3.3]

Parameters:
    n : the interpolating factor, must be between -1 and 1
    y : a sequence of three values

Results:
    the interpolated value of y
"""
function interpolate3(n, y)
    if (n < -1) || (n > 1)
    	error("Interpolating factor out of range: $n")
    end
    a = y[2] - y[1]
    b = y[3] - y[2]
    c = b - a
    return y[2] + n/2 * (a + b + n*c)
end

"""
    interpolate_angle3(n, y)

Interpolate from three equally spaced tabular angular values.

[Meeus-1998 equation 3.3]

This version is suitable for interpolating from a table of
angular values which may cross the origin of the circle,
for example: 359 degrees...0 degrees...1 degree.

Parameters:
    n : the interpolating factor, must be between -1 and 1
    y : a sequence of three values

Results:
    the interpolated value of y
"""
function interpolate_angle3(n, y)
    if (n < -1) || (n > 1)
    	error("Interpolating factor $n out of range")
    end
    a = diff_angle(y[1], y[2])
    b = diff_angle(y[2], y[3])
    c = diff_angle(a, b)
    return y[2] + n/2 * (a + b + n * c)
end

"""
    polynomial(coefficients, x)

Evaluate a simple polynomial.

Where: terms[1] is constant, terms[2] is for x, terms[3] is for x^2, etc.

Example:

```
y = polynomial((1.1, 2.2, 3.3, 4.4), t)
```

returns the value of:

```
1.1 + 2.2 * t + 3.3 * t^2 + 4.4 * t^3
```

Parameters:
    terms : sequence of coefficients
    x : variable value

Results:
    value of the polynomial
"""
function polynomial(coefficients, x)
    result = 0
    for (power, coefficient) in enumerate(coefficients)
        result += (x ^ (power - 1)) * coefficient
    end
    return result
end

"""

Alternatively use this from Julia:

Example

```
y = @polynomial_horner(x, p1, p2, p3)
```

evaluates

```
p[1] + x *  (p[2] + x * (....)),
```

Notice the arguments are differently ordered from `polynomial()` above!

Parameters:
    x : variable value
    term1 : coefficient 1
    ...

Results:
    value of the polynomial
"""
macro polynomial_horner(x, p...)
    ex = esc(p[end])
    for i = length(p)-1:-1:1
        ex = :($(esc(p[i])) + t * $ex)
    end
    Expr(:block, :(t = $(esc(x))), ex)
end


text:
# push!(LOAD_PATH,joinpath("..","src"))
using Multilane
using Test
using POMDPs
# using Gallium
using Debugger
using Random
using LinearAlgebra

# import Iterators.product
# import Base.assert

##############
##UNIT TESTS##
##############

function Multilane.MLState(pos::Real, vel::Real, cars::Array{CarState,1}, x::Real=50.)
    #x = mdp.phys_param.lane_length/2.
    insert!(cars,1,CarState(x, pos, vel, 0, Multilane.NORMAL, 0))
    return MLState(0.0, 0.0, cars, nothing)
end
# function MLState(crashed::Bool,pos::Real,vel::Real,cars::Array{CarState,1},x::Real=50.)
#     insert!(cars,1,CarState(x, pos, vel, 0, NORMAL, 0))
#     return MLState(false, crashed, 0.0, 0.0, cars)
# end
function Multilane.MLState(crashed::Bool,x::Float64,t::Float64,pos::Real,vel::Real,cars::Array{CarState,1},x_ego::Real=50.)
    insert!(cars,1,CarState(x_ego, pos, vel, 0, Multilane.NORMAL, 0))
    if crashed
        return MLState(x, t, cars, :crashed)
    else
        return MLState(x, t, cars)
    end
end

function test_hashing(s::AbstractString, ps)
	##TODO: make a more comprehensive test that uses more of the created objects
	println("\t\tTesting $s Hashing")
	if length(ps) < 2
		error("Make More objects for testing $s hashing")
	end
	d1 = Dict{typeof(ps[1]),Int}()
	d1[ps[1]] = 1
	d1[ps[1]] = 1
	@assert(length(d1)==1)
	d2 = Dict{typeof(ps[1]),Int}(p=>3 for p in [ps[2], ps[2]])
	@assert(length(d2)==1)
	@assert(get(d1,ps[1],0)==1)
end

function test_equality(s::AbstractString,ps)
	println("\t\tTesting $s Equality")
	if length(ps) < 3
		error("Make More objects for testing $s equality")
	end
	for i = 1:length(ps)
		for j = 1:length(ps)
			if i == j
				@assert(ps[i]==deepcopy(ps[j]))
			else
				@assert(ps[i] != ps[j])
				@assert(ps[j] != ps[i])
			end
		end
	end
	#@assert(ps[1] == ps[1])
	#@assert(ps[2] != ps[3])
	#@assert(ps[3] != ps[2])
end

include("test_idm.jl")
include("test_mobil.jl")
include("test_types.jl")
include("test_crash.jl")
include("test_model.jl")
include("test_actions.jl")

println("Running Unit Tests...")
test_idm()
test_mobil()
test_mdp_types()
test_crash()
test_model()
test_actions()

# println("Running an OriginalMDP example...")
# include("run_example.jl")

# include("change_into.jl") # broken in 0.6 :(

println("Running a NoCrashMDP example...")
include("run_nocrash.jl")
# write_tmp_gif(mdp, sim)
# println("Running a NoCrashMDP example with the robust solver...")
# include("run_with_robust_solver.jl")


text:
function get_total_profit(u::String="*")
    tmp = glob(joinpath(Pinguu.base_dir, "users/" * u * "/*/closed_orders/*/*.json"));
    total_profit = 0
    #
    ## Pinguus
    for file in tmp
        month = JSON.parsefile(file)
        for (day, info) in month
            if !(day in ["profit" "trade_num"])
                for (id_str, trade) in info["trades"]
                    total_profit += haskey(trade, "profit_EUR") ? trade["profit_EUR"] : trade["profit"]
                end
            end
        end
    end
    #
    ## Grid bots
    tmp = glob(joinpath(Pinguu.base_dir, "grid_bots/user/" * u * "/*/config.json"));
    for grid_bot_path in tmp
        tmp_config = JSON.parsefile(grid_bot_path)
        if tmp_config["testing"] == "0"
            tmp_history = glob(joinpath(dirname(grid_bot_path), "history/*/*json"))
            for tmp_history_file in tmp_history
                month = JSON.parsefile(tmp_history_file)
                for (day, trades) in month
                    for trade in trades
                        total_profit += trade["profit"]
                    end
                end
            end

        end
    end
    return total_profit
end
export get_total_profit


text:
mapinteraction!(multidata::MultiSimData, interaction::Interaction) = begin
    nrows, ncols = framesize(multidata)

    # Only pass in the data that the interaction wants, in that order
    interactiondata = map(key -> multidata[key], keys(interaction))

    for j in 1:ncols, i in 1:nrows
        ismasked(multidata, i, j) && continue
        state = map(d -> source(d)[i, j], interactiondata)
        newstate = applyinteraction(interaction, interactiondata, state, (i, j))
        map(interactiondata, newstate) do d, s
            @inbounds dest(d)[i, j] = s
        end
    end
end

mapinteraction!(multidata::MultiSimData, interaction::PartialInteraction) = begin
    nrows, ncols = framesize(multidata)

    # Only pass in the data that the interaction wants, in that order
    keys_ = keys(interaction)
    interactiondata = @set multidata.data =
        NamedTuple{keys_}(map(key -> WritableSimData(multidata[key]), keys_))

    # Copy all the source and dest arrays
    map(data(interactiondata)) do d
        @inbounds parent(dest(d)) .= parent(source(d))
    end

    # Loop manually. TODO: use blockrun?
    for j in 1:ncols, i in 1:nrows
        ismasked(multidata, i, j) && continue
        state = map(d -> source(d)[i, j], data(interactiondata))
        applyinteraction!(interaction, interactiondata, state, (i, j))
    end

    map(d -> updatestatus!(sourcestatus(d), deststatus(d)), data(interactiondata))

end


text:
export CalArray

function make_test_cal_array(name:: Symbol, N :: Integer;
    allFixed = false, offset = 0.0)

    rng = MersenneTwister(23);
    sizeV = 2 : (1 + N);
    sz = (sizeV..., );

    T1 = Float64;
    defaultValue = randn(rng, T1, sz) .+ offset;
    lb = defaultValue .- one(T1);
    ub = defaultValue .+ one(T1);
    if allFixed
        isCalM = trues(size(lb));
    else
        isCalM = Matrix{Bool}(defaultValue .> 0.0);
    end

    ca = CalArray(name, string(name), string(name), defaultValue .+ 0.1, defaultValue, lb, ub, isCalM);
    return ca
end

## -------------  Basics

Base.size(ca :: CalArray{T1, N}) where {T1, N} =
    size(ca.defaultValue);
Base.show(io :: IO,  p :: CalArray{T1, N}) where {T1, N} =
    print(io, "CalArray:  " * show_string(p));

is_calibrated(ca :: CalArray{T1, N}) where {T1, N} = any(ca.isCalM);

# Returns a Vector with the calibrated parameter values (or empty).
calibrated_value(ca :: CalArray{T1, N}) where {T1, N} =
    calibrated_elements(ca, ca.value);

calibrated_lb(ca :: CalArray{T1, N}) where {T1, N} =
    calibrated_elements(ca, ca.lb);

calibrated_ub(ca :: CalArray{T1, N}) where {T1, N} =
    calibrated_elements(ca, ca.ub);

function calibrated_elements(ca :: CalArray{T1, N}, m :: AbstractArray{T1, N}) where {T1, N}
    if is_calibrated(ca)
        v = m[map_indices(ca)];
    else
        v = Array{T1,N}();
    end
    return v
end

## --------  Validate

function validate(ca :: CalArray{T1, N}; silent = true) where {T1, N}
    isValid = true;
    isValid = isValid  && (size(ca.defaultValue) == size(ca.lb) == size(ca.ub) == size(ca.isCalM));
    # isValid = isValid  &&  check_values(ca);
    return isValid
end

# function check_values(ca :: CalArray{T1, N}) where {T1, N}
#     isValid = check_calibrated_values(ca)  &&  check_fixed_values(ca);
#     return isValid
# end

# function check_calibrated_values(ca :: CalArray{T1, N}) where {T1, N}
#     idxV = map_indices(ca);
#     isempty(idxV)  &&  return true;

#     calValueV = ca.valueM[idxV];
#     isValid = isapprox(calValueV, ca.calValueV);

#     pvec = get_pvector(ca);
#     isValid = isValid  &&  isapprox(calValueV, param_value(pvec, :calValueV));
#     return isValid
# end

# function check_fixed_values(ca :: CalArray{T1, N}) where {T1, N}
#     valueM = values(ca);
#     isValid = true;
#     for (j, isCal) in enumerate(ca.isCalM)
#         if !isCal
#             isValid = isValid && (valueM[j] ≈ ca.defaultValue[j]);
#         end
#     end
#     return isValid
# end

## ------------  Retrieve

function value(ca :: CalArray{T1, N}) where {T1, N}
    # update_values!(ca);
    return ca.value;
end

# function update_values!(ca :: CalArray{T1, N}) where {T1, N}
#     ca.value .= ca.defaultValue;
#     idxV = map_indices(ca);
#     for (j, idx) in enumerate(idxV)
#         ca.value[idx] = ca.value[];
#     end
# end

# Linear index into arrays for each calibrated value.
function map_indices(ca :: CalArray{T1, N}) where {T1, N}
    nCal = sum(ca.isCalM);
    if nCal > 0
        idxV = zeros(Int, nCal);
        j2 = 0;
        for (j, isCal) in enumerate(ca.isCalM)
            if isCal
                j2 += 1;
                idxV[j2] = j;
            end
        end
        @assert j2 == nCal;
    else
        idxV = Vector{Int}();
    end
    return idxV
end

## ----------  Change / update

"""
    $(SIGNATURES)

Change calibration status to `true`
"""
function calibrate!(p :: CalArray{T1, N}) where {T1, N}
    p.isCalM .= true
    return nothing
end

"""
    $(SIGNATURES)

Change calibration status to `false`
"""
function fix!(p :: CalArray{T1, N}; pValue = nothing) where {T1, N}
    p.isCalM .= false;
    if !isnothing(pValue)
        # set_value!(p, pValue);
        set_default_value!(p, pValue);
    end
    return nothing
end

# Used in calibration. Input should match dimension of calibrated params.
function set_calibrated_value!(ca :: CalArray{T1, N}, vIn :: AbstractVector{T1};
    skipInvalidSize = false) where {T1, N}

    idxV = map_indices(ca);
    if size(idxV) == size(vIn)
        ca.value[idxV] .= vIn;
    else
        @warn """
            Invalid size of new calibrated values for $ca
            Given: $(size(vIn))
            Expected: $(size(idxV))
            """
        skipInvalidSize  ||  error("Stopped");
    end
end

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


text:
#==
Example taken from
Messac et. al., Amir Ismail-Yahaya, and Christopher A Mattson.
The normalized normal constraint method for generating the pareto frontier.
Structural and multidisciplinary optimization, 25(2):86–98, 2003.

For the NC method, it is necessary to scale the objectives to get
an evenly spaced Pareto front. This example shows that this is
not necessary with the NBI method.
==#

using MultiJuMP, JuMP
using Ipopt

m = multi_model(Ipopt.Optimizer)
@variable(m, x[i=1:2])
@NLexpression(m, f1, x[1])
@NLexpression(m, f2, x[2])
@NLconstraint(m, ((x[1]-20)/20)^8+(x[2]-1)^8 <= 1)

obj1 = SingleObjective(f1, sense = MOI.MIN_SENSE)
obj2 = SingleObjective(f2, sense = MOI.MIN_SENSE)

multim = get_multidata(m)
multim.objectives = [obj1, obj2]
multim.pointsperdim = 20
optimize!(m, method = NBI())

using Plots
plot(multim)


text:
# ProductSpace
#--------------
# Tensor product of several ElementarySpace objects
immutable ProductSpace{S<:ElementarySpace,N} <: TensorSpace{S,N}
    spaces::NTuple{N, S}
end

# Additional constructors
ProductSpace{S<:ElementarySpace}(V::S,Vlist::S...) = ProductSpace(tuple(V,Vlist...))
ProductSpace(P::ProductSpace) = P

# This one is needed for unambiguous treatment of zero-dimensional ProductSpaces
ProductSpace{S<:ElementarySpace,N}(::Type{S}, spaces::NTuple{N,S}) = ProductSpace{S,N}(spaces)

# Corresponding methods
dim(P::ProductSpace) = (d::Int=1;for V in P;d*=dim(V);end;return d)
iscnumber(P::ProductSpace) = all(iscnumber,P)

dim{S<:UnitaryRepresentationSpace,G<:Sector,N}(P::ProductSpace{S,N},s::NTuple{N,G})=_dim(P.spaces,s)
sectors{S<:UnitaryRepresentationSpace}(P::ProductSpace{S})=_sectors(P.spaces)
sectortype{S<:UnitaryRepresentationSpace}(P::ProductSpace{S})=sectortype(S)
sectortype{S<:UnitaryRepresentationSpace,N}(::Type{ProductSpace{S,N}})=sectortype(S)

# Convention on dual, conj, transpose and ctranspose of tensor product spaces
dual{S,N}(P::ProductSpace{S,N}) = ProductSpace{S,N}(ntuple(n->dual(P[n]),N))
Base.conj{S,N}(P::ProductSpace{S,N}) = ProductSpace{S,N}(ntuple(n->conj(P[n]),N))

Base.transpose(P::ProductSpace) = reverse(P)
Base.ctranspose(P::ProductSpace) = reverse(conj(P))

# Default construction from product of spaces:
⊗{S<:ElementarySpace}(V1::S, V2::S) = ProductSpace(S, (V1, V2))
⊗{S<:ElementarySpace}(P1::ProductSpace{S}, V2::S) = ProductSpace(S, tuple(P1.spaces..., V2))
⊗{S<:ElementarySpace}(V1::S, P2::ProductSpace{S}) = ProductSpace(S, tuple(V1, P2.spaces...))
⊗{S<:ElementarySpace}(P1::ProductSpace{S}, P2::ProductSpace{S}) = ProductSpace(S, tuple(P1.spaces..., P2.spaces...))

# Functionality for extracting and iterating over spaces
Base.length{S,N}(P::ProductSpace{S,N}) = N
Base.endof(P::ProductSpace) = length(P)
Base.getindex(P::ProductSpace, n::Integer) = P.spaces[n]
Base.getindex{S<:ElementarySpace}(P::ProductSpace{S}, r) = ProductSpace(S, P.spaces[r])

Base.reverse(P::ProductSpace)=ProductSpace(reverse(P.spaces))
Base.map(f::Base.Callable,P::ProductSpace) = map(f,P.spaces) # required to make map(dim,P) efficient

Base.start(P::ProductSpace) = start(P.spaces)
Base.next(P::ProductSpace, state) = next(P.spaces, state)
Base.done(P::ProductSpace, state) = done(P.spaces, state)

# Comparison
==(P1::ProductSpace,P2::ProductSpace) = P1.spaces==P2.spaces
==(P::ProductSpace,V::ElementarySpace) = length(P)==1 && P[1] == V
==(V::ElementarySpace,P::ProductSpace) = length(P)==1 && P[1] == V

# Show method
function Base.show(io::IO, P::ProductSpace)
    print(io,"ProductSpace(")
    for i in 1:length(P)
        i==1 || print(io," ⊗ ")
        show(io, P[i])
    end
    print(io, ")")
end

# # Promotion and conversion
# Base.convert{S<:ElementarySpace}(::Type{ProductSpace{S,1}}, V::S) = ProductSpace(V)
# Base.convert{S<:ElementarySpace}(::Type{ProductSpace{S}}, V::S) = ProductSpace(V)
# Base.convert(::Type{ProductSpace}, V::ElementarySpace) = ProductSpace(V)
#
# Base.promote_rule{S<:ElementarySpace,N}(::Type{ProductSpace{S,N}},::Type{S}) = ProductSpace{S}
# Base.promote_rule{S<:ElementarySpace}(::Type{ProductSpace{S}},::Type{S}) = ProductSpace{S}
#

# # basis and basisvector
# typealias ProductBasisVector{S,N} BasisVector{ProductSpace{S,N},Int} # use integer from 1 to dim as identifier
# typealias ProductBasis{S,N} Basis{ProductSpace{S,N}}
#
# Base.length{S,N}(B::ProductBasis{S,N}) = dim(space(B))
# Base.start{S,N}(B::ProductBasis{S,N}) = 1
# Base.next{S,N}(B::ProductBasis{S,N}, state::Int) = (ProductBasisVector{S,N}(space(B),state),state+1)
# Base.done{S,N}(B::ProductBasis{S,N}, state::Int) = state>length(B)
#
# Base.to_index{S,N}(b::ProductBasisVector{S,N}) = b.identifier # use linear indexing as long as we cannot efficiently generate a cartesian iterator
# Base.show{S,N}(io::IO,b::ProductBasisVector{S,N}) = print(io, "BasisVector($(b.space),$(ind2sub(map(dim,b.space),b.identifier)))")


text:
module TestDataFrame
    using Base.Test
    using DataFrames

    #test_group("Operations on DataFrames that have column groupings")

    x = DataFrame(a = [1, 2, 3], b = [4, 5, 6])
    y = DataFrame(c = [1, 2, 3], d = [4, 5, 6])

    #
    # Equality
    #
    @test isequal(DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])), DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])))
    @test !isequal(DataFrame(a=@data([1, 2]), b=@data([4, 5])), DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])))
    @test !isequal(DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])), DataFrame(a=@data([1, 2, 3])))
    @test !isequal(DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])), DataFrame(a=@data([1, 2, 3]), c=@data([4, 5, 6])))
    @test !isequal(DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])), DataFrame(b=@data([4, 5, 6]), a=@data([1, 2, 3])))
    @test !isequal(DataFrame(a=@data([1, 2, 2]), b=@data([4, 5, 6])), DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])))
    @test isequal(DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6])), DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6])))

    @test DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])) == DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6]))
    @test DataFrame(a=@data([1, 2]), b=@data([4, 5])) != DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6]))
    @test DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])) != DataFrame(a=@data([1, 2, 3]))
    @test DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])) != DataFrame(a=@data([1, 2, 3]), c=@data([4, 5, 6]))
    @test DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])) != DataFrame(b=@data([4, 5, 6]), a=@data([1, 2, 3]))
    @test DataFrame(a=@data([1, 2, 2]), b=@data([4, 5, 6])) != DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6]))
    @test DataFrame(a=@data([1, 3, NA]), b=@data([4, 5, 6])) != DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6]))
    @test isna(DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6])) == DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6])))
    @test isna(DataFrame(a=@data([1, 2, NA]), b=@data([4, 5, 6])) == DataFrame(a=@data([1, 2, 3]), b=@data([4, 5, 6])))

    z = deepcopy(x)

    z = hcat(x, y)

    v = DataFrame(
        a = [5,6,7],
        b = [8,9,10]
    )
    z = vcat(DataFrame[v, x])

    z = vcat(v, x)

    ## del calls ref, which properly deals with groupings
    z2 = z[:,[1,1,2]]
    @test names(z2) == [:a, :a_1, :b]

    #test_group("DataFrame assignment")
    # Insert single column
    x0 = x[[], :]
    @test_throws ErrorException x0[:d] = [1]
    @test_throws ErrorException x0[:d] = 1:3

    # Insert single value
    x[:d] = 3
    @test x[:d] == [3, 3, 3]

    x0[:d] = 3
    @test x0[:d] == Int[]

    # Associative methods

    df = DataFrame(a=[1, 2], b=[3., 4.])
    @test haskey(df, :a)
    @test !haskey(df, :c)
    @test get(df, :a, -1) === df.columns[1]
    @test get(df, :c, -1) == -1
    @test keys(df) == [:a, :b]
    @test values(df) == {df[:a], df[:b]}
    @test !isempty(df)

    @test empty!(df) === df
    @test isempty(df.columns)
    @test isempty(df)

    df = DataFrame(a=[1, 2], b=[3., 4.])
    @test_throws BoundsError insert!(df, 5, ["a", "b"], :newcol)
    @test_throws ErrorException insert!(df, 1, ["a"], :newcol)
    @test insert!(df, 1, ["a", "b"], :newcol) == df
    @test isequal(df, DataFrame(newcol=["a", "b"], a=[1, 2], b=[3., 4.]))
    df = DataFrame(a=[1, 2], b=[3., 4.])
    @test insert!(df, 3, ["a", "b"], :newcol) == df
    @test isequal(df, DataFrame(a=[1, 2], b=[3., 4.], newcol=["a", "b"]))

    df = DataFrame(a=[1, 2], b=[3., 4.])
    df2 = DataFrame(c=["a", "b"], d=[:c, :d])
    @test insert!(df, df2) == df
    @test isequal(df, DataFrame(a=[1, 2], b=[3., 4.], c=["a", "b"], d=[:c, :d]))

    #test_group("Empty DataFrame constructors")
    df = DataFrame(Int, 10, 3)
    @test size(df, 1) == 10
    @test size(df, 2) == 3
    @test typeof(df[:, 1]) == DataVector{Int}
    @test typeof(df[:, 2]) == DataVector{Int}
    @test typeof(df[:, 3]) == DataVector{Int}

    df = DataFrame({Int, Float64, ASCIIString}, 100)
    @test size(df, 1) == 100
    @test size(df, 2) == 3
    @test typeof(df[:, 1]) == DataVector{Int}
    @test typeof(df[:, 2]) == DataVector{Float64}
    @test typeof(df[:, 3]) == DataVector{ASCIIString}

    df = DataFrame({Int, Float64, ASCIIString}, [:A, :B, :C], 100)
    @test size(df, 1) == 100
    @test size(df, 2) == 3
    @test typeof(df[:, 1]) == DataVector{Int}
    @test typeof(df[:, 2]) == DataVector{Float64}
    @test typeof(df[:, 3]) == DataVector{ASCIIString}

    df = convert(DataFrame, zeros(10, 5))
    @test size(df, 1) == 10
    @test size(df, 2) == 5
    @test typeof(df[:, 1]) == Vector{Float64}

    df = convert(DataFrame, ones(10, 5))
    @test size(df, 1) == 10
    @test size(df, 2) == 5
    @test typeof(df[:, 1]) == Vector{Float64}

    df = convert(DataFrame, eye(10, 5))
    @test size(df, 1) == 10
    @test size(df, 2) == 5
    @test typeof(df[:, 1]) == Vector{Float64}

    #test_group("Other DataFrame constructors")
    df = DataFrame([{:a=>1, :b=>'c'}, {:a=>3, :b=>'d'}, {:a=>5}])
    @test size(df, 1) == 3
    @test size(df, 2) == 2
    @test typeof(df[:,:a]) == DataVector{Int}
    @test typeof(df[:,:b]) == DataVector{Char}

    df = DataFrame([{:a=>1, :b=>'c'}, {:a=>3, :b=>'d'}, {:a=>5}], [:a, :b])
    @test size(df, 1) == 3
    @test size(df, 2) == 2
    @test typeof(df[:,:a]) == DataVector{Int}
    @test typeof(df[:,:b]) == DataVector{Char}

    data = {:A => [1, 2], :C => [:1, :2], :B => [3, 4]}
    df = DataFrame(data)
    # Specify column_names
    df = DataFrame(data, [:C, :A, :B])

    # This assignment was missing before
    df = DataFrame(Column = [:A])
    df[1, :Column] = "Testing"

    # zero-row dataframe and subdataframe test
    df = DataFrame(x=[], y=[])
    @test nrow(df) == 0
    df = DataFrame(x=[1:3], y=[3:5])
    sdf = sub(df, df[:x] .== 4)
    @test size(sdf, 1) == 0

    @test hash(convert(DataFrame, [1 2; 3 4])) == hash(convert(DataFrame, [1 2; 3 4]))
    @test hash(convert(DataFrame, [1 2; 3 4])) != hash(convert(DataFrame, [1 3; 2 4]))

    # push!(df, row)
    df=DataFrame( first=[1,2,3], second=["apple","orange","pear"] )

    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    push!(dfb, {3,"pear"})
    @test df==dfb

    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    push!(dfb, (3,"pear"))
    @test df==dfb

    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    @test_throws ArgumentError push!(dfb, (33.33,"pear"))

    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    @test_throws ArgumentError push!(dfb, ("coconut",22))

    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    push!(dfb, [ :first=>3,:second=>"pear" ])
    @test df==dfb

    df=DataFrame( first=[1,2,3], second=["apple","orange","banana"] )
    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    push!(dfb, [ "first"=>3,"second"=>"banana" ])
    @test df==dfb

    df0= DataFrame( first=[1,2], second=["apple","orange"] )
    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    @test_throws ArgumentError push!(dfb, [ :first=>true,:second=>false ])
    @test df0==dfb

    df0= DataFrame( first=[1,2], second=["apple","orange"] )
    dfb= DataFrame( first=[1,2], second=["apple","orange"] )
    @test_throws ArgumentError push!(dfb, ["first"=>"chicken", "second"=>"stuff"] )
    @test df0==dfb

    # delete!
    df = DataFrame(a=1, b=2, c=3, d=4, e=5)
    @test_throws ArgumentError delete!(df, 0)
    @test_throws ArgumentError delete!(df, 6)
    @test_throws KeyError delete!(df, :f)

    d = copy(df)
    delete!(d, [:a, :e, :c])
    @test names(d) == [:b, :d]
    delete!(d, :b)
    @test isequal(d, df[[:d]])

    d = copy(df)
    delete!(d, [2, 5, 3])
    @test names(d) == [:a, :d]
    delete!(d, 2)
    @test isequal(d, df[[:a]])

    # deleterows!
    df = DataFrame(a=[1, 2], b=[3., 4.])
    @test deleterows!(df, 1) === df
    @test isequal(df, DataFrame(a=[2], b=[4.]))

    df = DataFrame(a=[1, 2], b=[3., 4.])
    @test deleterows!(df, 2) === df
    @test isequal(df, DataFrame(a=[1], b=[3.]))

    df = DataFrame(a=[1, 2, 3], b=[3., 4., 5.])
    @test deleterows!(df, 2:3) === df
    @test isequal(df, DataFrame(a=[1], b=[3.]))

    df = DataFrame(a=[1, 2, 3], b=[3., 4., 5.])
    @test deleterows!(df, [2, 3]) === df
    @test isequal(df, DataFrame(a=[1], b=[3.]))

    df = DataFrame(a=@data([1, 2]), b=@data([3., 4.]))
    @test deleterows!(df, 1) === df
    @test isequal(df, DataFrame(a=@data([2]), b=@data([4.])))

    df = DataFrame(a=@data([1, 2]), b=@data([3., 4.]))
    @test deleterows!(df, 2) === df
    @test isequal(df, DataFrame(a=@data([1]), b=@data([3.])))

    df = DataFrame(a=@data([1, 2, 3]), b=@data([3., 4., 5.]))
    @test deleterows!(df, 2:3) === df
    @test isequal(df, DataFrame(a=@data([1]), b=@data([3.])))

    df = DataFrame(a=@data([1, 2, 3]), b=@data([3., 4., 5.]))
    @test deleterows!(df, [2, 3]) === df
    @test isequal(df, DataFrame(a=@data([1]), b=@data([3.])))
end


text:
HODIndex(client::HODClient, name::AbstractString) = HODIndex(client, utf8(name))

Base.length(index::HODIndex) = length(index.documents)

function Base.push!(index::HODIndex, document::UTF8String)
    push!(index.documents, document)
end

function Base.append!(index::HODIndex, documents::Vector{UTF8String})
    append!(index.documents, documents)
end

function commit(index::HODIndex; async::Bool = false)
    response = add_documents(index, async = async)
    empty!(index.documents)
    return response
end

add_document(index::HODIndex; async::Bool = false) = add_documents([doc], async)

function add_documents(index::HODIndex; async::Bool = false)
    documents = Dict("documents" => index.documents)
    data = Dict("json" => json(documents), "index" => index.name)
    response = post(
        index.client,
        "addtotextindex",
        data = data,
        files = Dict("fake" => ""),
        async = async
    )
end

Base.delete!(index::HODIndex) = delete_index(index.client, index.name)


text:
using MyInterpolations
using Base.Test

include("test_lin_interp.jl")