text:
module Gates

export Gate, ID, X, Y, Z, H, PHASE, S, T, CZ, RX, RY, RZ, CNOT, CCNOT, SWAP, CSWAP, CCSWAP, controlled, noisify, damp_amplitude
export damping_residual_kraus, damping_kraus, damping_kraus_map

const empty_matrix = Matrix{Complex{Float64}}(undef, 0, 0)
const empty_kraus = Vector{Matrix{Complex{Float64}}}(undef, 0)

struct Gate
	# Gate has either matrix or SWAP bit indexes set. Control bits can be set in either case.
	matrix::Matrix{Complex{Float64}}
	matrix_func::Function
	qubit_index::Integer
	control_bit_indexes::AbstractVector{Int64}
	swap_bit_index_a::Integer
	swap_bit_index_b::Integer
	kraus_operators::AbstractVector{Matrix{Complex{Float64}}}

	# Constructor for ideal gates except for SWAP.
	function Gate(matrix::Matrix{Complex{Float64}},
		qubit_index::Integer,
		control_bit_indexes::AbstractArray{Int64,1} = Array{Int64,1}([]))

		return new(matrix, ()->(), qubit_index, control_bit_indexes, -1, -1, empty_kraus)
	end

	# Constructor for parametric gates.
	function Gate(matrix_func::Function,
		qubit_index::Integer,
		control_bit_indexes::AbstractArray{Int64,1} = Array{Int64,1}([]))

		return new(empty_matrix, matrix_func, qubit_index, control_bit_indexes, -1, -1, empty_kraus)
	end

	# Constructor for SWAP gates, including controlled once.
	function Gate(swap_bit_index_a::Integer,
		swap_bit_index_b::Integer,
		control_bit_indexes::AbstractArray{Int64,1} = Array{Int64,1}([]))

		return new(empty_matrix, ()->(), -1, control_bit_indexes, swap_bit_index_a, swap_bit_index_b, empty_kraus)
	end

	# Constructor for noisy gates.
	function Gate(matrix::Matrix{Complex{Float64}},
		qubit_index::Integer,
		kraus_operators::AbstractVector{Matrix{Complex{Float64}}},
		control_bit_indexes::AbstractArray{Int64,1} = Array{Int64,1}([]))

		return new(matrix, ()->(), qubit_index, control_bit_indexes, -1, -1, kraus_operators)
	end

	# Constructor that populates all properties; used to create a new gate based on existing one.
	function Gate(matrix::Matrix{Complex{Float64}},
		matrix_func::Function,
		qubit_index::Integer,
		control_bit_indexes::AbstractArray{Int64,1},
		swap_bit_index_a::Integer,
		swap_bit_index_b::Integer,
		kraus_operators::AbstractVector{Matrix{Complex{Float64}}})

		return new(matrix, matrix_func, qubit_index, control_bit_indexes, swap_bit_index_a, swap_bit_index_b, kraus_operators)
	end
end

# Identity gate.
# Named ID instead of I to avoid name conflict with LinearAlgebra.I.
function ID(qubit_index::Integer)::Gate
	return Gate(Matrix{Complex{Float64}}([1 0; 0 1]), qubit_index)
end

function X(qubit_index::Integer)::Gate
	return Gate(Matrix{Complex{Float64}}([0 1; 1 0]), qubit_index)
end

function Y(qubit_index::Integer)::Gate
	return Gate(Matrix{Complex{Float64}}([0 complex(0, -1); complex(0, 1) 0]), qubit_index)
end

function Z(qubit_index::Integer)::Gate
	return Gate(Matrix{Complex{Float64}}([1 0; 0 -1]), qubit_index)
end

function H(qubit_index::Integer)::Gate
	return Gate(Matrix{Complex{Float64}}(1 / √2 * [1 1; 1 -1]), qubit_index)
end

function PHASE(angle::Real, qubit_index::Integer)::Gate
	return Gate(Matrix{Complex{Float64}}([1 0; 0 exp(complex(0, 1) * angle)]), qubit_index)
end

function S(qubit_index::Integer)::Gate
	return Gate(Matrix{Complex{Float64}}([1 0; 0 complex(0, 1)]), qubit_index)
end

function T(qubit_index::Integer)::Gate
	return Gate(Matrix{Complex{Float64}}([1 0; 0 1 / √2 + 1im / √2]), qubit_index)
end

function CZ(control_index::Integer, target_index::Integer)::Gate
	return Gate(Matrix{Complex{Float64}}([1 0; 0 -1]), target_index, [control_index])
end

function RX(angle::Real, qubit_index::Integer)::Gate
	return Gate(Matrix{Complex{Float64}}([cos(angle / 2) complex(0, -1) * sin(angle / 2); complex(0, -1) * sin(angle / 2) cos(angle / 2)]), qubit_index)
end

function RX(param_name::String, qubit_index::Integer)::Gate
	matrix_func = function (params)
		angle = params[param_name]
		return Matrix{Complex{Float64}}([cos(angle / 2) complex(0, -1) * sin(angle / 2); complex(0, -1) * sin(angle / 2) cos(angle / 2)])
	end
	return Gate(matrix_func, qubit_index)
end

function RY(angle::Real, qubit_index::Integer)::Gate
	return Gate(Matrix{Complex{Float64}}([cos(angle / 2) -sin(angle / 2); sin(angle / 2) cos(angle / 2)]), qubit_index)
end

function RY(param_name::String, qubit_index::Integer)::Gate
	matrix_func = function (params)
		angle = params[param_name]
		return Matrix{Complex{Float64}}([cos(angle / 2) -sin(angle / 2); sin(angle / 2) cos(angle / 2)])
	end
	return Gate(matrix_func, qubit_index)
end

function RZ(angle::Real, qubit_index::Integer)::Gate
	return Gate(Matrix{Complex{Float64}}([exp(complex(0, -1) * angle / 2) 0; 0 exp(complex(0, 1) * angle / 2)]), qubit_index)
end

function RZ(param_name::String, qubit_index::Integer)::Gate
	matrix_func = function (params)
		angle = params[param_name]
		return Matrix{Complex{Float64}}([exp(complex(0, -1) * angle / 2) 0; 0 exp(complex(0, 1) * angle / 2)])
	end
	return Gate(matrix_func, qubit_index)
end

function CNOT(control_index::Integer, target_index::Integer)::Gate
	return Gate(Matrix{Complex{Float64}}([0 1; 1 0]), target_index, [control_index])
end

function CCNOT(control_index_a::Integer, control_index_b::Integer, target_index::Integer)::Gate
	return Gate(Matrix{Complex{Float64}}([0 1; 1 0]), target_index, [control_index_a, control_index_b])
end

function SWAP(swap_bit_index_a::Integer, swap_bit_index_b::Integer)::Gate
	return Gate(swap_bit_index_a, swap_bit_index_b)
end

function CSWAP(control_index::Integer, swap_bit_index_a::Integer, swap_bit_index_b::Integer)::Gate
	return Gate(swap_bit_index_a, swap_bit_index_b, [control_index])
end

function CCSWAP(control_index_a::Integer, control_index_b::Integer, swap_bit_index_a::Integer, swap_bit_index_b::Integer)::Gate
	return Gate(swap_bit_index_a, swap_bit_index_b, [control_index_a, control_index_b])
end

function controlled(control_bit_index::Integer, g::Gate)::Gate
	return controlled([control_bit_index], g)
end

function controlled(control_bit_indexes::AbstractArray{Int64,1}, g::Gate)::Gate
	return Gate(g.matrix,
		g.matrix_func,
		g.qubit_index,
		vcat(g.control_bit_indexes, control_bit_indexes),
		g.swap_bit_index_a,
		g.swap_bit_index_b,
		g.kraus_operators)
end

function noisify(g::Gate, kraus_operators::AbstractVector{Matrix{Complex{Float64}}})::Gate
	if size(g.kraus_operators)[1] > 0
		error("cannot apply Kraus operators to a gate that already has Kraus operators applied to")
	end
	return Gate(g.matrix,
		g.matrix_func,
		g.qubit_index,
		g.control_bit_indexes,
		g.swap_bit_index_a,
		g.swap_bit_index_b,
		kraus_operators)
end

damping_residual_kraus(decay_1_to_0_probability = .1) = Matrix{Complex{Float64}}([1 0; 0 √(1 - decay_1_to_0_probability)])
damping_kraus(decay_1_to_0_probability = .1) = Matrix{Complex{Float64}}([0 √decay_1_to_0_probability; 0 0])
damping_kraus_map(decay_1_to_0_probability = .1) = [damping_residual_kraus(decay_1_to_0_probability), damping_kraus(decay_1_to_0_probability)]

# damp_amplitude adds noise via amplitude damping kraus operators with given probability of decay |1> to |0>.
function damp_amplitude(g::Gate, ket1_to_ket0_decay_probability::Real)
	return noisify(g, damping_kraus_map(ket1_to_ket0_decay_probability))
end

end # module


text:
function chebyshev_extrema(n::S,domain = [1.0,-1.0]) where S <: Integer

  if n <= 0
    error("The number of nodes must be positive.")
  end

  nodes = zeros(n)

  if isodd(n)
    nodes[Int((n-1)/2)+1] = (domain[1]+domain[2])/2.0
  end

  for i = 1:div(n,2)
    x = -cos((i-1)*pi/(n-1))*(domain[1]-domain[2])/2.0
    nodes[i]       = (domain[1]+domain[2])/2.0  + x
    nodes[end-i+1] = (domain[1]+domain[2])/2.0  - x
  end

  return nodes

end

ext: jl

max_stars_repo_name: RJDennis/ChebyshevApprox

max_issues_repo_name: RJDennis/ChebyshevApprox

max_forks_repo_name: RJDennis/ChebyshevApprox


text:
"""
Instance that is internally stored as an `MathOptInterface` model.

    abstract type AbstractMOIBackedInstance <: AbstractInstance end

To interact with the underlying `MathOptInterface` model, use:

    get_moi_model(model::AbstractMOIBackedInstance)::MOI.ModelLike
"""
abstract type AbstractMOIBackedInstance <: AbstractInstance end

function get_moi_model end

"""
A "simple" `AbstractMOIBackedInstance` where the underlying `MathOptInterface`
model is stored in a field.

    struct SimpleMOIBackedInstance{T <: MOI.ModelLike} <: AbstractInstance
        model::T
    end
"""
struct SimpleMOIBackedInstance{T<:MOI.ModelLike} <: AbstractInstance
    model::T
end

get_moi_model(inst::SimpleMOIBackedInstance) = inst.model


text:
@inline front(t::Tuple, m::Integer) = ntuple(i -> t[i], m)
@inline front(t::Tuple, ::StaticOrVal{M}) where {M} = ntuple(i -> t[i], Val(M))
@inline front(t::NTuple{L,Any}) where {L} = front(t, Val(L - 1))

@inline function tail(t::Tuple, n::Integer)
    m = length(t) - n
    ntuple(i -> t[i+m], n)
end
@inline function tail(t::NTuple{L,Any}, ::StaticOrVal{N}) where {L,N}
    M = L - N
    ntuple(i -> t[i+M], Val(N))
end
@inline tail(t::NTuple{L,Any}) where {L} = tail(t, Val(L - 1))

@inline function tuple_split(t::NTuple{L,Any}, ::StaticOrVal{M}) where {L,M}
    front(t, Val(M)), tail(t, Val(L - M))
end
@inline function tuple_split(t::NTuple{L,Any}, m::Integer) where {L}
    front(t, Val(m)), tail(t, L - m)
end

@pure tuple_length(T::Type{<:Tuple}) = length(T.parameters)
@pure tuple_length(t::Tuple) = length(t)

@pure tuple_prod(T::Type{<:Tuple}) = length(T.parameters) == 0 ? 1 : *(T.parameters...)
@pure tuple_prod(t::Tuple) = prod(t)

@pure function tuple_minimum(T::Type{<:Tuple})
    length(T.parameters) == 0 ? 0 : minimum(tuple(T.parameters...))
end
@pure tuple_minimum(t::Tuple) = minimum(t)


text:
const LiteralAMQP = UInt8[65, 77, 81, 80] # "AMQP"
const ProtocolId = UInt8(0)
const ProtocolVersion = UInt8[0, 9, 1]

const ProtocolHeader = vcat(LiteralAMQP, ProtocolId, ProtocolVersion)

const ContentWeight = 0x0000
const FrameEnd = 0xCE
const HeartBeat = UInt8[8, 0, 0, FrameEnd]

abstract type TAMQPLengthPrefixed end

#const TAMQPBit                  = UInt8
const TAMQPBool                 = UInt8 # 0 = FALSE, else TRUE
const TAMQPScale                = UInt8 # number of decimal digits
const TAMQPOctet                = UInt8
const TAMQPShortShortInt        = UInt8
const TAMQPShortShortUInt       = UInt8
const TAMQPShortInt             = Int16
const TAMQPShortUInt            = UInt16
const TAMQPLongInt              = Int32
const TAMQPLongUInt             = UInt32
const TAMQPLongLongInt          = Int64
const TAMQPLongLongUInt         = UInt64
const TAMQPFloat                = Float32
const TAMQPDouble               = Float64
const TAMQPTimeStamp            = TAMQPLongLongUInt

struct TAMQPBit
    val::UInt8
end

function TAMQPBit(b::TAMQPBit, pos::Int)
    TAMQPBit((b.val >> (pos-1)) & 0x1)
end

function TAMQPBit(b::TAMQPBit, setbit::TAMQPBit, pos::Int)
    TAMQPBit(b.val | (setbit.val << (pos-1)))
end
TAMQPBit(b::Bool) = TAMQPBit(UInt8(b))
TAMQPBit(b::T) where {T<:Integer} = TAMQPBit(Bool(b))

struct TAMQPDecimalValue
    scale::TAMQPScale
    val::TAMQPLongUInt
end

struct TAMQPShortStr <: TAMQPLengthPrefixed
    len::TAMQPOctet
    data::Vector{UInt8}
end
TAMQPShortStr(d::Vector{UInt8}) = TAMQPShortStr(length(d), d)
TAMQPShortStr(s::AbstractString) = TAMQPShortStr(Vector{UInt8}(codeunits(String(s))))

struct TAMQPLongStr <: TAMQPLengthPrefixed
    len::TAMQPLongUInt
    data::Vector{UInt8}
end
TAMQPLongStr(d::Vector{UInt8}) = TAMQPLongStr(length(d), d)
TAMQPLongStr(s::AbstractString) = TAMQPLongStr(Vector{UInt8}(codeunits(String(s))))

struct TAMQPByteArray <: TAMQPLengthPrefixed
    len::TAMQPLongUInt
    data::Vector{UInt8}
end
TAMQPByteArray(d::Vector{UInt8}) = TAMQPByteArray(length(d), d)
TAMQPByteArray(s::AbstractString) = TAMQPByteArray(Vector{UInt8}(codeunits(String(s))))

const TAMQPFieldName = TAMQPShortStr
const TAMQPFV = Union{Real, TAMQPDecimalValue, TAMQPLengthPrefixed, Nothing}

struct TAMQPFieldValue{T <: TAMQPFV}
    typ::Char  # as in FieldValueIndicatorMap
    fld::T
end
TAMQPFieldValue(v::T) where {T} = TAMQPFieldValue{T}(FieldIndicatorMap[T], v)
TAMQPFieldValue(v::Dict) = TAMQPFieldValue(TAMQPFieldTable(v))
TAMQPFieldValue(v::String) = TAMQPFieldValue(TAMQPLongStr(v))
TAMQPFieldValue(v::Bool) = TAMQPFieldValue('b', TAMQPBool(v))

struct TAMQPFieldValuePair{T <: TAMQPFV}
    name::TAMQPFieldName
    val::TAMQPFieldValue{T}
end

struct TAMQPFieldArray <: TAMQPLengthPrefixed
    len::TAMQPLongInt
    data::Vector{TAMQPFieldValue}
end
TAMQPFieldArray(data::Vector{TAMQPFieldValue}) = TAMQPFieldArray(length(data), data)

struct TAMQPFieldTable <: TAMQPLengthPrefixed
    len::TAMQPLongUInt
    data::Vector{TAMQPFieldValuePair}
end
TAMQPFieldTable(data::Vector{TAMQPFieldValuePair}) = TAMQPFieldTable(length(data), data)
TAMQPFieldTable(dict::Dict) = TAMQPFieldTable(TAMQPFieldValuePair[TAMQPFieldValuePair(TAMQPShortStr(String(n)), TAMQPFieldValue(v)) for (n,v) in dict])

const TAMQPField = Union{TAMQPBit, Integer, TAMQPShortStr, TAMQPLongStr, TAMQPFieldTable}

const FieldValueIndicatorMap = Dict{Char,DataType}(
    't' => TAMQPBool,
    'b' => TAMQPShortShortInt,
    'B' => TAMQPShortShortUInt,
    'U' => TAMQPShortInt,
    'u' => TAMQPShortUInt,
    'I' => TAMQPLongInt,
    'i' => TAMQPLongUInt,
    'L' => TAMQPLongLongInt,
    'l' => TAMQPLongLongUInt,
    'f' => TAMQPFloat,
    'd' => TAMQPDouble,
    'D' => TAMQPDecimalValue,
    's' => TAMQPShortStr,
    'S' => TAMQPLongStr,
    'x' => TAMQPByteArray,
    'A' => TAMQPFieldArray,
    'T' => TAMQPTimeStamp,
    'F' => TAMQPFieldTable,
    'V' => Nothing
)

const FieldIndicatorMap = Dict{DataType,Char}(v=>n for (n,v) in FieldValueIndicatorMap)

const TAMQPChannel          = TAMQPShortUInt
const TAMQPPayloadSize      = TAMQPLongUInt
const TAMQPContentBodySize  = TAMQPLongLongUInt
const TAMQPClassId          = UInt16
const TAMQPMethodId         = UInt16
const TAMQPContentClass     = TAMQPClassId

struct TAMQPFrameProperties
    channel::TAMQPChannel
    payloadsize::TAMQPPayloadSize
end

struct TAMQPPropertyFlags
    flags::UInt16
    nextval::Union{TAMQPPropertyFlags, Nothing}
end
TAMQPPropertyFlags(flags::UInt16) = TAMQPPropertyFlags(flags, nothing)

struct TAMQPBodyPayload
    # TODO: may be better to allow sub arrays, for efficient writing of large messages
    data::Vector{TAMQPOctet}
end

struct TAMQPMethodPayload
    class::TAMQPClassId
    method::TAMQPMethodId
    fields::Vector{Pair{Symbol,TAMQPField}}

    TAMQPMethodPayload(p::TAMQPBodyPayload) = TAMQPMethodPayload(p.data)
    TAMQPMethodPayload(b::Vector{TAMQPOctet}) = TAMQPMethodPayload(IOBuffer(b))
    function TAMQPMethodPayload(io)
        class = ntoh(read(io, TAMQPClassId))
        method = ntoh(read(io, TAMQPMethodId))
        args = methodargs(class, method)
        fields = Vector{Pair{Symbol,TAMQPField}}(undef, length(args))
        @debug("reading method payload", class, method, nargs=length(args))
        bitpos = 0
        bitval = TAMQPBit(0)
        for idx in 1:length(fields)
            fld = args[idx]
            @debug("reading", field=fld.first, type=fld.second)
            if fld.second === TAMQPBit
                bitpos += 1
                (bitpos == 1) && (bitval = read(io, fld.second))
                v = TAMQPBit(bitval, bitpos)
                (bitpos == 8) && (bitpos == 0)
            else
                bitpos = 0
                v = read(io, fld.second)
            end
            (fld.second <: Integer) && (v = ntoh(v))
            fields[idx] = Pair{Symbol,TAMQPField}(fld.first, v)
        end
        new(class, method, fields)
    end
    function TAMQPMethodPayload(class_name::Symbol, method_name::Symbol, fldvals)
        class = CLASSNAME_MAP[class_name]
        method = CLASSMETHODNAME_MAP[(class_name,method_name)]
        fields = Pair{Symbol,TAMQPField}[]
        for idx in 1:length(method.args)
            (argname,argtype) = method.args[idx]
            argval = fldvals[idx]
            push!(fields, Pair{Symbol,TAMQPField}(argname, isa(argval, argtype) ? argval : argtype(argval)))
        end
        new(class.id, method.id, fields)
    end
end

struct TAMQPHeaderPayload
    class::TAMQPContentClass
    weight::UInt16  # must be ContentWeight
    bodysize::TAMQPContentBodySize
    propflags::TAMQPPropertyFlags
    proplist::Dict{Symbol,TAMQPField}

    TAMQPHeaderPayload(p::TAMQPBodyPayload) = TAMQPHeaderPayload(p.data)
    TAMQPHeaderPayload(b::Vector{TAMQPOctet}) = TAMQPHeaderPayload(IOBuffer(b))
    function TAMQPHeaderPayload(io)
        class = ntoh(read(io, TAMQPClassId))
        wt = ntoh(read(io, UInt16))
        @assert wt === ContentWeight
        bodysize = ntoh(read(io, TAMQPContentBodySize))
        propflags = TAMQPPropertyFlags(ntoh(read(io, UInt16)))
        proplist = Dict{Symbol,TAMQPField}()

        flags = propflags.flags
        for prop in SORTED_PROPERTIES
            if (flags & prop.mask) > 0x0000
                proplist[prop.name] = read(io, prop.typ)
            end
        end
        new(class, ContentWeight, bodysize, propflags, proplist)
    end
    function TAMQPHeaderPayload(class::TAMQPContentClass, message)
        bodysize = length(message.data)
        flags = 0x0000
        for name in keys(message.properties)
            flags = flags | PROPERTIES[name].mask
        end
        new(class, ContentWeight, bodysize, TAMQPPropertyFlags(flags), message.properties)
    end
end

# Generic frame, used to read any frame
struct TAMQPGenericFrame
    hdr::UInt8
    props::TAMQPFrameProperties
    payload::TAMQPBodyPayload
    fend::UInt8 # must be FrameEnd
end

# Type = 1, "METHOD": method frame
struct TAMQPMethodFrame
    props::TAMQPFrameProperties
    payload::TAMQPMethodPayload
end

function TAMQPMethodFrame(f::TAMQPGenericFrame)
    @debug("Frame Conversion: generic => method")
    @assert f.hdr == FrameMethod
    TAMQPMethodFrame(f.props, TAMQPMethodPayload(f.payload))
end

function TAMQPGenericFrame(f::TAMQPMethodFrame)
    @debug("Frame Conversion method => generic")
    iob = IOBuffer()
    methpayload = f.payload
    write(iob, hton(methpayload.class))
    write(iob, hton(methpayload.method))
    bitpos = 0
    bitval = TAMQPBit(0)
    for (n,v) in methpayload.fields
        if isa(v, TAMQPBit)
            bitpos += 1
            bitval = TAMQPBit(bitval, v, bitpos)
            if bitpos == 8
                write(iob, bitval)
                bitpos = 0
                bitval = TAMQPBit(0)
            end
        else
            if bitpos > 0
                write(iob, bitval)
                bitpos = 0
                bitval = TAMQPBit(0)
            end
            (typeof(v) <: Integer) && (v = hton(v))
            write(iob, v)
        end
    end
    if bitpos > 0
        write(iob, bitval)
    end
    bodypayload = TAMQPBodyPayload(take!(iob))
    TAMQPGenericFrame(FrameMethod, TAMQPFrameProperties(f.props.channel, length(bodypayload.data)), bodypayload, FrameEnd)
end

# Type = 2, "HEADER": content header frame.
struct TAMQPContentHeaderFrame
    props::TAMQPFrameProperties
    hdrpayload::TAMQPHeaderPayload
end

function TAMQPContentHeaderFrame(f::TAMQPGenericFrame)
    @debug("Frame Conversion: generic => contentheader")
    @assert f.hdr == FrameHeader
    TAMQPContentHeaderFrame(f.props, TAMQPHeaderPayload(f.payload))
end

function TAMQPGenericFrame(f::TAMQPContentHeaderFrame)
    @debug("Frame Conversion contentheader => generic")
    iob = IOBuffer()
    hdrpayload = f.hdrpayload
    propflags = hdrpayload.propflags
    proplist = hdrpayload.proplist

    write(iob, hton(hdrpayload.class))
    write(iob, hton(hdrpayload.weight))
    write(iob, hton(hdrpayload.bodysize))
    write(iob, hton(propflags.flags))

    flags = propflags.flags
    for prop in SORTED_PROPERTIES
        if (flags & prop.mask) > 0x0000
            write(iob, proplist[prop.name])
        end
    end
    bodypayload = TAMQPBodyPayload(take!(iob))
    TAMQPGenericFrame(FrameHeader, TAMQPFrameProperties(f.props.channel, length(bodypayload.data)), bodypayload, FrameEnd)
end

# Type = 3, "BODY": content body frame.
struct TAMQPContentBodyFrame
    props::TAMQPFrameProperties
    payload::TAMQPBodyPayload
end

function TAMQPContentBodyFrame(f::TAMQPGenericFrame)
    @debug("Frame Conversion: generic => contentbody")
    @assert f.hdr == FrameBody
    TAMQPContentBodyFrame(f.props, f.payload)
end

function TAMQPGenericFrame(f::TAMQPContentBodyFrame)
    @debug("Frame Conversion contentbody => generic")
    TAMQPGenericFrame(FrameBody, TAMQPFrameProperties(f.props.channel, length(f.payload.data)), f.payload, FrameEnd)
end

# Type = 4, "HEARTBEAT": heartbeat frame.
struct TAMQPHeartBeatFrame
end

function TAMQPHeartBeatFrame(f::TAMQPGenericFrame)
    @assert f.hdr == FrameHeartbeat
    TAMQPHeartBeatFrame()
end

function TAMQPGenericFrame(f::TAMQPHeartBeatFrame)
    @debug("Frame Conversion heartbeat => generic")
    TAMQPGenericFrame(FrameHeartbeat, TAMQPFrameProperties(DEFAULT_CHANNEL, 0), TAMQPBodyPayload(TAMQPOctet[]), FrameEnd)
end

struct TAMQPContent
    hdr::TAMQPContentHeaderFrame
    body::Vector{TAMQPContentBodyFrame}
end

struct TAMQPMethod
    frame::TAMQPMethodFrame
    content::Union{TAMQPContent, Nothing}
end

# Exceptions
mutable struct AMQPProtocolException <: Exception
    msg::String
end

mutable struct AMQPClientException <: Exception
    msg::String
end

# Spec code gen types
struct MethodSpec
    id::Int
    name::Symbol
    respname::Symbol
    args::Vector{Pair{Symbol,DataType}}
end

struct ClassSpec
    id::Int
    name::Symbol
    method_map::Dict{Int,MethodSpec}
end

struct CloseReason
    code::Int16
    msg::TAMQPShortStr
    classid::TAMQPClassId
    methodid::TAMQPMethodId
end

# Utility Methods for Types
method(classid::TAMQPClassId, methodid::TAMQPMethodId) = CLASS_MAP[classid].method_map[methodid]
methodargs(classid::TAMQPClassId, methodid::TAMQPMethodId) = method(classid, methodid).args
function displayname(classid::TAMQPClassId, methodid::TAMQPMethodId)
    c = CLASS_MAP[classid]
    m = c.method_map[methodid]
    "$(c.name).$(m.name)"
end


text:
module Spinner

# Imports
using IterTools

# Constants
# Could have supplied Inf , but Inf in Julia is a float!
MAX_CYCLES = 1000000000

# List of avaliable styles for spinner
circles = "○◔◑◕●"
corners = "⌜⌝⌟⌞"
dots = "⠋⠙⠹⠸⠼⠴⠦⠧⠇⠏"
moon = "🌑🌒🌓🌔🌕🌖🌗🌘"
quadrants = "▖▘▝▗'"
vertical_bars = "▁▂▃▄▅▆▇█"

"""
Spinner function

Sample usage:-
julia> using Spinner
julia> Spinner.spinner(Spinner.corners)

"""
function spinner(style)
	try
		for i in ncycle(style, MAX_CYCLES)
	           write(STDOUT, i);
	           flush
	           sleep(0.07)
	           write(STDOUT, "\b")
	    end
	catch
		nothing
	end
end

# Exports
export spinner
end


text:
#=
Created on Wednesday 9 Janurary 2020
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Trait matching model based on data from Olito et al. 2015.
=#

using CSV, Plots, LaTeXStrings
using OTSIN

interactions = CSV.read("data/Olito2015_SM_A2_rawdata_full/Aggregate network-Table 1.csv", header=2)

plant_morph = CSV.read("data/Olito2015_SM_A2_rawdata_full/Plant morphology-Table 1.csv")
pol_morph = CSV.read("data/Olito2015_SM_A2_rawdata_full/Pollinator morphology-Table 1.csv")

plant_densities = CSV.read("data/Olito2015_SM_A2_rawdata_full/Plant densities-Table 1.csv")

Y = Matrix(interactions[:,2:end])
P = Y / sum(Y)

# match of class
M = Matrix(plant_morph[:,2:5]) * Matrix(pol_morph[:,2:5])'

a, b = marginals(P)

histplants = histogram(sum(Y, dims=2)[:], color=:green, title="interactions plants",
    xlabel="number of interactions", ylabel=:counts, label="");
histpols = histogram(sum(Y, dims=1)[:], color=:yellow,  title="interactions pollinators",
    xlabel="number of interactions", ylabel=:counts, label="");
pmarg = plot(histplants, histpols, layout=(2,1))

λs = 2 .^(-10:0.1:0)

klvals = [KL(P, optimaltransport(M, a, b, λ=λ)) for λ in λs]

λopt = λs[argmin(klvals)]

pλ = plot(λs, klvals, xscale=:log10, lw=2, label="optimal transport")
ylabel!("Kullback-Leibler\ndivergence (nats)")
xlabel!("\$ \\lambda \$")

Q = optimaltransport(M, a, b, λ=λopt)

kl_optimal = minimum(klvals)
kl_neutral = KL(P, optimaltransport(a,b))

hline!([kl_neutral], lw=2, label="neutral", ls=:dash)

pM = heatmap(M, xlabel="plant sp.", ylabel="pol. sp.", title="M", color=:grays) #, colorbar=:none
pP = heatmap(P, xlabel="plant sp.", ylabel="pol. sp.", title="P", color=:heat)
pQ = heatmap(Q, xlabel="plant sp.", ylabel="pol. sp.", title="Q", color=:heat)

l = @layout [  a{0.5w} b{0.5w};
               c{0.3w} d{0.7w}]
plot(pP, pQ, pM, pλ, figsize=(300, 600))

savefig("figures/traitmatch.pdf")

max_stars_repo_name: PoisotLab/OTSIN

max_issues_repo_name: PoisotLab/OTSIN

max_forks_repo_name: PoisotLab/OTSIN


text:
include("examples/examples.jl")
include("3rd_party/3rd_party.jl")
include("wrappers/wrappers.jl")


text:
using Plots
using Dash, DashHtmlComponents, DashCoreComponents

plotly()

app = dash()

app.layout = html_div() do
    dcc_slider(
        id="angle_selector",
        min=0,
        max=2pi,
        step=0.01,
        value=pi/2,
    ),
    dcc_graph(id="plot")
end

callback!(
        app,
        Output("plot", "figure"),
        Input("angle_selector", "value")) do angle
    angles = range(0, angle, length=100)
    x = cos.(angles)
    y = sin.(angles)
    fig = plot(
        x, y,
        size=(500, 500),
        xlims=(-1.5, 1.5), ylims=(-1.5, 1.5),
        legend=:none,
        framestyle=:none)
    Plots.plotlybase_syncplot(fig)
    Plots.backend_object(fig)
end

run_server(app)


text:
# create genomic relationship matrix based on the genotype
using JWAS,DataFrames,CSV,Statistics,JWAS.Datasets,DelimitedFiles,Random, LinearAlgebra, StatsBase
Random.seed!(1)

for rep in 1:20
    genofile="/Users/tianjing/Box/Omicsdata/andrew1055/data$rep/geno.csv"
    genotypes=get_genotypes(genofile;separator=',',method="RR-BLUP");

    sum2pq=genotypes.sum2pq

    geno_centered=genotypes.genotypes
    round.(mean(geno_centered,dims=1),digits=2)

    GRM=(geno_centered*geno_centered')/sum2pq
    writedlm("/Users/tianjing/Box/Omicsdata/andrew1055/data$rep/GRM_geno.txt",GRM)
end

max_stars_repo_name: zhaotianjing/NN-LMM

max_issues_repo_name: zhaotianjing/NN-LMM

max_forks_repo_name: zhaotianjing/NN-LMM


text:
using ArgParse, DelimitedFiles

function nextParaRun()
    s = ArgParseSettings(description = "arg parse here for parallel experimentation")
    @add_arg_table! s begin
        "--location", "-l"  # location
    end
    parsed_args = parse_args(s) # the result is a Dict{String,Any}
    dir = parsed_args["location"]

    if isfile("$(dir)/counter.txt") == false
        mkdir("$(dir)")
        writedlm("$(dir)/counter.txt", "1")
    end
    c = Int(readdlm("$(dir)/counter.txt")[1])
    writedlm("$(dir)/counter.txt", c + 1)
    newdir = "$(dir)/$(c)"
    mkdir(newdir)
    println(newdir)
end

nextParaRun()

max_stars_repo_name: hipersys-team/arrow

max_issues_repo_name: arrow-bow/arrow

max_forks_repo_name: arrow-bow/arrow


text:
export find

import Base: show

"""
    find(conditionFn::F) where { F <: Function }

Creates a find operator, which emits only the first value emitted by the source Observable that meets some condition.

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Arguments
- `conditionFn::F`: condition function with `(data::T) -> Bool` signature

# Examples
```jldoctest
using Rocket

source = from([ 1, 2, 3, 4, 5, 6, 7, 8 ])
subscribe!(source |> find((d) -> d % 2 == 0), logger())
;

# output

[LogActor] Data: 2
[LogActor] Completed
```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`logger`](@ref)
"""
find(conditionFn::F) where { F <: Function } = filter(conditionFn) + take(1)


text:
using Bolt
using BenchmarkTools

par = CosmoParams()
bg = Background(par)
ih = IonizationHistory(Peebles(), par, bg)

k_grid = quadratic_k(0.1bg.H₀, 1000bg.H₀, 100)
k_i = 10
hierarchy = Hierarchy(BasicNewtonian(), par, bg, ih, k_grid[k_i])

xᵢ = log(1e-8)
u₀ = Bolt.initial_conditions(xᵢ, hierarchy)
udummy = deepcopy(u₀)
@btime Bolt.hierarchy!(udummy, u₀, hierarchy, xᵢ)

##
@time perturb = boltsolve(hierarchy);

##
@time sf = source_grid(par, bg, ih, k_grid, BasicNewtonian());

##


text:
using Test
using BSplines
using BSplines: KnotVector, StandardBasisVector, check_intervalindex
using LinearAlgebra: ldiv!, lmul!, rdiv!, rmul!
using OffsetArrays: OffsetArray

include("knotvector.jl")
include("standardbasisvector.jl")

@testset "intervalindex" begin
    function test_intervalindex(basis::BSplineBasis, x, leftknot)
        @test intervalindex(basis, x) === leftknot
        for start = eachindex(knots(basis))
            @test intervalindex(basis, x, start) === leftknot
        end
        @test_throws BoundsError intervalindex(basis, x, firstindex(knots(basis))-1)
        @test_throws BoundsError intervalindex(basis, x, lastindex(knots(basis))+1)
    end
    function test_intervalindex(vec, x, leftknot)
        @test intervalindex(vec, x) === leftknot
        for start = eachindex(vec)
            @test intervalindex(vec, x, start) === leftknot
        end
        @test_throws BoundsError intervalindex(vec, x, firstindex(vec)-1)
        @test_throws BoundsError intervalindex(vec, x, lastindex(vec)+1)
    end
    test_intervalindex(-10:10, -Inf32,         nothing)
    test_intervalindex(-10:10, -11,            nothing)
    test_intervalindex(-10:10, -10.0,          1)
    test_intervalindex(-10:10, -9.5f0,         1)
    test_intervalindex(-10:10, big(π),         14)
    test_intervalindex(-10:10, prevfloat(9.0), 19)
    test_intervalindex(-10:10, 9,              20)
    test_intervalindex(-10:10, 10.0,           20)
    test_intervalindex(-10:10, 10.1,           nothing)
    test_intervalindex(-10:10, NaN16,          nothing)
    test_intervalindex(-10:10, big(Inf),       nothing)
    test_intervalindex(-5//1:1//2:5//1, -Inf16,         nothing)
    test_intervalindex(-5//1:1//2:5//1, -5.1,           nothing)
    test_intervalindex(-5//1:1//2:5//1, -5.0,           1)
    test_intervalindex(-5//1:1//2:5//1, -4.5f0,         2)
    test_intervalindex(-5//1:1//2:5//1, 1.6,            14)
    test_intervalindex(-5//1:1//2:5//1, prevfloat(4.5), 19)
    test_intervalindex(-5//1:1//2:5//1, 9//2,           20)
    test_intervalindex(-5//1:1//2:5//1, 5.0,            20)
    test_intervalindex(-5//1:1//2:5//1, 5.1,            nothing)
    test_intervalindex(-5//1:1//2:5//1, NaN32,          nothing)
    test_intervalindex(-5//1:1//2:5//1, Inf64,          nothing)
    test_intervalindex(10_000_000.1:0.001:10_000_000.2, big(-Inf),                 nothing)
    test_intervalindex(10_000_000.1:0.001:10_000_000.2, 10_000_000,                nothing)
    test_intervalindex(10_000_000.1:0.001:10_000_000.2, 10_000_000.1,              1)
    test_intervalindex(10_000_000.1:0.001:10_000_000.2, prevfloat(10_000_000.123), 23)
    test_intervalindex(10_000_000.1:0.001:10_000_000.2, 10_000_000.123,            24)
    test_intervalindex(10_000_000.1:0.001:10_000_000.2, 10_000_000.15,             51)
    test_intervalindex(10_000_000.1:0.001:10_000_000.2, 10_000_000.1987,           99)
    test_intervalindex(10_000_000.1:0.001:10_000_000.2, 10_000_000.199,            100)
    test_intervalindex(10_000_000.1:0.001:10_000_000.2, 10_000_000.19999,          100)
    test_intervalindex(10_000_000.1:0.001:10_000_000.2, 10_000_000.2,              100)
    test_intervalindex(10_000_000.1:0.001:10_000_000.2, nextfloat(10_000_000.2),   nothing)
    test_intervalindex(10_000_000.1:0.001:10_000_000.2, big(NaN),                  nothing)
    test_intervalindex(10_000_000.1:0.001:10_000_000.2, Inf16,                     nothing)
    test_intervalindex([0, 1, 2, 3, 4, 4, 5, 6, 6], -Inf16, nothing)
    test_intervalindex([0, 1, 2, 3, 4, 4, 5, 6, 6], -1,     nothing)
    test_intervalindex([0, 1, 2, 3, 4, 4, 5, 6, 6], 0.0,    1)
    test_intervalindex([0, 1, 2, 3, 4, 4, 5, 6, 6], 0.5,    1)
    test_intervalindex([0, 1, 2, 3, 4, 4, 5, 6, 6], 1.0,    2)
    test_intervalindex([0, 1, 2, 3, 4, 4, 5, 6, 6], 4.0,    6)
    test_intervalindex([0, 1, 2, 3, 4, 4, 5, 6, 6], 4.5,    6)
    test_intervalindex([0, 1, 2, 3, 4, 4, 5, 6, 6], 6.0,    7)
    test_intervalindex([0, 1, 2, 3, 4, 4, 5, 6, 6], 6.1,    nothing)
    test_intervalindex([0, 1, 2, 3, 4, 4, 5, 6, 6], Inf64,  nothing)
    test_intervalindex([0, 1, 2, 3, 4, 4, 5, 6, 6], NaN,    nothing)
    test_intervalindex([5, 5, 5, 5, 5], 5, nothing)
    test_intervalindex(KnotVector(-5:5, 2), -Inf32,         nothing)
    test_intervalindex(KnotVector(-5:5, 2), -10,            nothing)
    test_intervalindex(KnotVector(-5:5, 2), -5,             3)
    test_intervalindex(KnotVector(-5:5, 2), -4.5,           3)
    test_intervalindex(KnotVector(-5:5, 2), 0,              8)
    test_intervalindex(KnotVector(-5:5, 2), 4.5,            12)
    test_intervalindex(KnotVector(-5:5, 2), 5,              12)
    test_intervalindex(KnotVector(-5:5, 2), nextfloat(5.0), nothing)
    test_intervalindex(KnotVector(-5:5, 2), NaN32,          nothing)
    test_intervalindex(KnotVector([0, 1, 2, 3, 4, 4, 5, 6], 4, 6), -1,  nothing)
    test_intervalindex(KnotVector([0, 1, 2, 3, 4, 4, 5, 6], 4, 6), 0.0, 5)
    test_intervalindex(KnotVector([0, 1, 2, 3, 4, 4, 5, 6], 4, 6), 0.5, 5)
    test_intervalindex(KnotVector([0, 1, 2, 3, 4, 4, 5, 6], 4, 6), 1.0, 6)
    test_intervalindex(KnotVector([0, 1, 2, 3, 4, 4, 5, 6], 4, 6), 4.0, 10)
    test_intervalindex(KnotVector([0, 1, 2, 3, 4, 4, 5, 6], 4, 6), 4.5, 10)
    test_intervalindex(KnotVector([0, 1, 2, 3, 4, 4, 5, 6], 4, 6), 6.0, 11)
    test_intervalindex(KnotVector([0, 1, 2, 3, 4, 4, 5, 6], 4, 6), 6.1, nothing)
    test_intervalindex(BSplineBasis(3, 0:5), -Inf, nothing)
    test_intervalindex(BSplineBasis(3, 0:5), -0.5, nothing)
    test_intervalindex(BSplineBasis(3, 0:5), 0//1, 3)
    test_intervalindex(BSplineBasis(3, 0:5), 0.5,  3)
    test_intervalindex(BSplineBasis(3, 0:5), 1,    4)
    test_intervalindex(BSplineBasis(3, 0:5), 4.0,  7)
    test_intervalindex(BSplineBasis(3, 0:5), 4.5,  7)
    test_intervalindex(BSplineBasis(3, 0:5), 5.0,  7)
    test_intervalindex(BSplineBasis(3, 0:5), 5.1,  nothing)
    test_intervalindex(BSplineBasis(3, 0:5), Inf,  nothing)
    test_intervalindex(BSplineBasis(3, 0:5), NaN,  nothing)
    test_intervalindex(BSplineBasis(4, [1,2,3,3,4,5,6]), 0.5,  nothing)
    test_intervalindex(BSplineBasis(4, [1,2,3,3,4,5,6]), 1,    4)
    test_intervalindex(BSplineBasis(4, [1,2,3,3,4,5,6]), 3,    7)
    test_intervalindex(BSplineBasis(4, [1,2,3,3,4,5,6]), 3.5,  7)
    test_intervalindex(BSplineBasis(4, [1,2,3,3,4,5,6]), 4.5,  8)
    test_intervalindex(BSplineBasis(4, [1,2,3,3,4,5,6]), 5.0,  9)
    test_intervalindex(BSplineBasis(4, [1,2,3,3,4,5,6]), 6.0,  9)
    test_intervalindex(BSplineBasis(4, [1,2,3,3,4,5,6]), 6.1,  nothing)
    test_intervalindex(BSplineBasis(4, [1,2,3,3,4,5,6]), NaN,  nothing)
    test_intervalindex(BSplineBasis(4, [1,2,3,3,4,5,6]), Inf,  nothing)
    test_intervalindex(BSplineBasis(4, [1,2,3,3,4,5,6]), -Inf, nothing)
    # check_intervalindex
    for vec = [[0, 1, 2, 3, 4, 4, 5, 6, 6], BSplineBasis(3, 0:5), BSplineBasis(4, [1,2,3,3,4,5,6])]
        for x = Real[-1, -0.5, 0, 0.5, 1, 3, 3.5, 4, 9//2, 5.0, 5.1, 6.0, 6.1, -Inf, Inf, NaN]
            for index = [-1:20; nothing]
                if intervalindex(vec, x) == index
                    @test check_intervalindex(Bool, vec, x, index)
                    @test check_intervalindex(vec, x, index) === nothing
                else
                    @test !check_intervalindex(Bool, vec, x, index)
                    @test_throws ArgumentError check_intervalindex(vec, x, index)
                end
            end
        end
    end
end

@testset "intervalindices" begin
    support = [(1, 3:3), (2, 3:4), (3, 3:5), (4, 4:6), (5, 5:7), (6, 6:7), (7, 7:7)]
    for basis = [BSplineBasis(3, 0:5), BSplineBasis(3, 1//2:1//2:3//1), BSplineBasis(3, 0.0:0.2:1.0)]
        for itr = [intervalindices(basis), intervalindices(basis, :), intervalindices(basis, 1:7),
                   intervalindices(basis, 1), intervalindices(basis, 1, 2)]
            @test itr isa UnitRange{Int}
        end
        @test intervalindices(basis)      == 3:7
        @test intervalindices(basis, :)   == 3:7
        @test intervalindices(basis, 1:7) == 3:7
        @test intervalindices(basis, 1:6) == 3:7
        @test intervalindices(basis, 3:7) == 3:7
        @test intervalindices(basis, 1:5) == 3:7
        @test intervalindices(basis, 2:4) == 3:6
        @test isempty(intervalindices(basis, 1:0))
        @test isempty(intervalindices(basis, 4:3))
        @test isempty(intervalindices(basis, 100:99))
        for (i, irange) = support
            @test intervalindices(basis, i:i) == irange
            @test intervalindices(basis, i) == irange
            for (j, jrange) = support
                @test intervalindices(basis, i, j) == irange ∩ jrange
                for (k, krange) = support
                    @test intervalindices(basis, i, j, k) == irange ∩ jrange ∩ krange
                end
            end
        end
        @test_throws BoundsError intervalindices(basis, 0:1)
        @test_throws BoundsError intervalindices(basis, 1:8)
        @test_throws BoundsError intervalindices(basis, 0)
        @test_throws BoundsError intervalindices(basis, 8)
        @test_throws BoundsError intervalindices(basis, 7, 0)
        @test_throws BoundsError intervalindices(basis, 1, 2, 8)
    end
    support = [(1, [4]), (2, [4,5]), (3, [4,5]), (4, [4,5,7]), (5, [5,7,8]),
               (6, [7,8,9]), (7, [7,8,9]), (8, [8,9]), (9, [9])]
    basis = BSplineBasis(4, [1,2,3,3,4,5,6])
    for itr = [intervalindices(basis), intervalindices(basis, :), intervalindices(basis, 1:7),
               intervalindices(basis, 1), intervalindices(basis, 1, 2)]
        @test eltype(typeof(itr)) === Int
        @test Base.IteratorEltype(typeof(itr)) === Base.HasEltype()
        @test Base.IteratorSize(typeof(itr)) === Base.SizeUnknown()
    end
    @test collect(intervalindices(basis)) == [4,5,7,8,9]
    @test collect(intervalindices(basis, :)) == [4,5,7,8,9]
    @test collect(intervalindices(basis, 1:9)) == [4,5,7,8,9]
    @test collect(intervalindices(basis, 1:8)) == [4,5,7,8,9]
    @test collect(intervalindices(basis, 4:8)) == [4,5,7,8,9]
    @test collect(intervalindices(basis, 1:6)) == [4,5,7,8,9]
    @test collect(intervalindices(basis, 5:6)) == [5,7,8,9]
    @test isempty(intervalindices(basis, 1:0))
    @test isempty(intervalindices(basis, 7:6))
    @test isempty(intervalindices(basis, -100:-101))
    for (i, irange) = support
        @test collect(intervalindices(basis, i:i)) == irange
        @test collect(intervalindices(basis, i)) == irange
        for (j, jrange) = support
            @test collect(intervalindices(basis, i, j)) == irange ∩ jrange
            for (k, krange) = support
                @test collect(intervalindices(basis, i, j, k)) == irange ∩ jrange ∩ krange
            end
        end
    end

    @testset "Iterators.reverse" begin
        support = [(1, 3:3), (2, 3:4), (3, 3:5), (4, 4:6), (5, 5:7), (6, 6:7), (7, 7:7)]
        for basis = [BSplineBasis(3, 0:5), BSplineBasis(3, 1//2:1//2:3//1), BSplineBasis(3, 0.0:0.2:1.0)]
            for itr = Iterators.reverse.([intervalindices(basis), intervalindices(basis, :),
                                          intervalindices(basis, 1:7), intervalindices(basis, 1),
                                          intervalindices(basis, 1, 2)])
                @test itr isa StepRange{Int}
            end
            @test Iterators.reverse(intervalindices(basis))      == 7:-1:3
            @test Iterators.reverse(intervalindices(basis, :))   == 7:-1:3
            @test Iterators.reverse(intervalindices(basis, 1:7)) == 7:-1:3
            @test Iterators.reverse(intervalindices(basis, 1:6)) == 7:-1:3
            @test Iterators.reverse(intervalindices(basis, 3:7)) == 7:-1:3
            @test Iterators.reverse(intervalindices(basis, 1:5)) == 7:-1:3
            @test Iterators.reverse(intervalindices(basis, 2:4)) == 6:-1:3
            @test isempty(Iterators.reverse(intervalindices(basis, 1:0)))
            @test isempty(Iterators.reverse(intervalindices(basis, 4:3)))
            @test isempty(Iterators.reverse(intervalindices(basis, 100:99)))
            for (i, irange) = support
                @test Iterators.reverse(intervalindices(basis, i:i)) == reverse(irange)
                @test Iterators.reverse(intervalindices(basis, i)) == reverse(irange)
                for (j, jrange) = support
                    @test Iterators.reverse(intervalindices(basis, i, j)) == reverse(irange ∩ jrange)
                    for (k, krange) = support
                        @test Iterators.reverse(intervalindices(basis, i, j, k)) == reverse(irange ∩ jrange ∩ krange)
                    end
                end
            end
        end
        support = [(1, [4]), (2, [4,5]), (3, [4,5]), (4, [4,5,7]), (5, [5,7,8]),
                   (6, [7,8,9]), (7, [7,8,9]), (8, [8,9]), (9, [9])]
        basis = BSplineBasis(4, [1,2,3,3,4,5,6])
        for itr = Iterators.reverse.([intervalindices(basis), intervalindices(basis, :),
                                      intervalindices(basis, 1:7), intervalindices(basis, 1),
                                      intervalindices(basis, 1, 2)])
            @test eltype(typeof(itr)) === Int
            @test Base.IteratorEltype(typeof(itr)) === Base.HasEltype()
            @test Base.IteratorSize(typeof(itr)) === Base.SizeUnknown()
        end
        @test collect(Iterators.reverse(intervalindices(basis))) == [9,8,7,5,4]
        @test collect(Iterators.reverse(intervalindices(basis, :))) == [9,8,7,5,4]
        @test collect(Iterators.reverse(intervalindices(basis, 1:9))) == [9,8,7,5,4]
        @test collect(Iterators.reverse(intervalindices(basis, 1:8))) == [9,8,7,5,4]
        @test collect(Iterators.reverse(intervalindices(basis, 4:8))) == [9,8,7,5,4]
        @test collect(Iterators.reverse(intervalindices(basis, 1:6))) == [9,8,7,5,4]
        @test collect(Iterators.reverse(intervalindices(basis, 5:6))) == [9,8,7,5]
        @test isempty(Iterators.reverse(intervalindices(basis, 1:0)))
        @test isempty(Iterators.reverse(intervalindices(basis, 7:6)))
        @test isempty(Iterators.reverse(intervalindices(basis, -100:-101)))
        for (i, irange) = support
            @test collect(Iterators.reverse(intervalindices(basis, i:i))) == reverse(irange)
            @test collect(Iterators.reverse(intervalindices(basis, i))) == reverse(irange)
            for (j, jrange) = support
                @test collect(Iterators.reverse(intervalindices(basis, i, j))) == reverse(irange ∩ jrange)
                for (k, krange) = support
                    @test collect(Iterators.reverse(intervalindices(basis, i, j, k))) == reverse(irange ∩ jrange ∩ krange)
                end
            end
        end
    end
end

@testset "support" begin
    # Basis
    @test support(BSplineBasis(3, 0:5)) === (0, 5)
    @test support(BSplineBasis(4, 0//1:1//10:5//1)) === (0//1, 5//1)
    @test support(BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])) === (-5.0, 3.6)
    # BSpline
    @test support(BSplineBasis(3, 0:5)[1]) === (0, 1)
    @test support(BSplineBasis(3, 0:5)[2]) === (0, 2)
    @test support(BSplineBasis(3, 0:5)[3]) === (0, 3)
    @test support(BSplineBasis(3, 0:5)[4]) === (1, 4)
    @test support(BSplineBasis(3, 0:5)[5]) === (2, 5)
    @test support(BSplineBasis(3, 0:5)[6]) === (3, 5)
    @test support(BSplineBasis(3, 0:5)[7]) === (4, 5)
    @test support(BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])[1]) === (-5.0, -3.5)
    @test support(BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])[2]) === (-5.0, -2.2)
    @test support(BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])[3]) === (-5.0, 1.0)
    @test support(BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])[4]) === (-5.0, 2.0)
    @test support(BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])[5]) === (-5.0, 3.6)
    @test support(BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])[6]) === (-3.5, 3.6)
    @test support(BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])[7]) === (-2.2, 3.6)
    @test support(BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])[8]) === (1.0, 3.6)
    @test support(BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])[9]) === (2.0, 3.6)
    # Spline
    @test support(Spline(BSplineBasis(3, 0:5), ones(7))) === (0, 5)
    @test support(Spline(BSplineBasis(4, 0//1:1//10:5//1), zeros(53))) === (0//1, 5//1)
    @test support(Spline(BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6]), 1:9)) === (-5.0, 3.6)
end

@testset "Equality and hashing" begin
    bases = [BSplineBasis(3, 0:5),
             BSplineBasis(3, [0.0,1.0,2.0,3.0,4.0,5.0]),
             BSplineBasis(3, 0//1:1//1:5//1),
             BSplineBasis(3, 0:0.5:5),
             BSplineBasis(4, 0:5),
             BSplineBasis(3, [1,2,3,3,4,5])]
    for b1 = bases, b2 = bases
        @test (b1 == b2) == (order(b1) == order(b2) && breakpoints(b1) == breakpoints(b2))
        @test (b1 == b2) == (hash(b1) == hash(b2))
    end
end

@testset "BSplineBasis" begin
    @testset "Construction" begin
        @test_throws DomainError BSplineBasis(0, 0:5) # order must be positive
        @test_throws InexactError BSplineBasis(5.5, 0:5) # order must be integer
        @test_throws TypeError BSplineBasis(5.5, [1+2im, 3+4im, 5+6im]) # breakpoints must be real
        @test_throws ArgumentError BSplineBasis(5, 1:0) # breakpoints must have length ≥ 2
        @test_throws ArgumentError BSplineBasis(5, [1]) # breakpoints must have length ≥ 2
        @test_throws ArgumentError BSplineBasis(5, OffsetArray(1:5, 1)) # breakpoints must not have offset axes
    end

    @testset "Dimensions" begin
        @test length(BSplineBasis(3, 0:5)) === 7
        @test length(BSplineBasis(3, big(0):big(5))) === 7
        @test length(BSplineBasis(5, [1,2,3.5,6,10])) === 8
    end

    @testset "Indices" begin
        @test firstindex(BSplineBasis(3, 0:5)) == 1
        @test firstindex(BSplineBasis(5, [1,2,3.5,6,10])) == 1
        @test lastindex(BSplineBasis(3, 0:5)) == 7
        @test lastindex(BSplineBasis(5, [1,2,3.5,6,10])) == 8
        @test eachindex(BSplineBasis(3, 0:5)) == 1:7
        @test eachindex(BSplineBasis(5, [1,2,3.5,6,10])) == 1:8
        @test keys(BSplineBasis(3, 0:5)) == 1:7
        @test keys(BSplineBasis(5, [1,2,3.5,6,10])) == 1:8
    end

    @testset "checkbounds" begin
        @test checkbounds(BSplineBasis(3, 0:5), :) === nothing
        @test checkbounds(BSplineBasis(3, 0:5), 2) === nothing
        @test checkbounds(BSplineBasis(3, 0:5), 1:7) === nothing
        @test checkbounds(BSplineBasis(3, 0:5), 2:7) === nothing
        @test_throws BoundsError checkbounds(BSplineBasis(3, 0:5), 0)
        @test_throws BoundsError checkbounds(BSplineBasis(3, 0:5), 8)
        @test_throws BoundsError checkbounds(BSplineBasis(3, 0:5), 2:8)
        @test checkbounds(BSplineBasis(5, [1,2,3.5,6,10]), :) === nothing
        @test checkbounds(BSplineBasis(5, [1,2,3.5,6,10]), 2) === nothing
        @test checkbounds(BSplineBasis(5, [1,2,3.5,6,10]), 1:8) === nothing
        @test checkbounds(BSplineBasis(5, [1,2,3.5,6,10]), 2:8) === nothing
        @test_throws BoundsError checkbounds(BSplineBasis(5, [1,2,3.5,6,10]), 0)
        @test_throws BoundsError checkbounds(BSplineBasis(5, [1,2,3.5,6,10]), 9)
        @test_throws BoundsError checkbounds(BSplineBasis(5, [1,2,3.5,6,10]), 2:9)
        @test checkbounds(Bool, BSplineBasis(3, 0:5), :)
        @test checkbounds(Bool, BSplineBasis(3, 0:5), 2)
        @test checkbounds(Bool, BSplineBasis(3, 0:5), 1:7)
        @test checkbounds(Bool, BSplineBasis(3, 0:5), 2:7)
        @test checkbounds(Bool, BSplineBasis(3, 0:5), 0) === false
        @test checkbounds(Bool, BSplineBasis(3, 0:5), 8) === false
        @test checkbounds(Bool, BSplineBasis(3, 0:5), 2:8) === false
        @test checkbounds(Bool, BSplineBasis(5, [1,2,3.5,6,10]), :)
        @test checkbounds(Bool, BSplineBasis(5, [1,2,3.5,6,10]), 2)
        @test checkbounds(Bool, BSplineBasis(5, [1,2,3.5,6,10]), 1:8)
        @test checkbounds(Bool, BSplineBasis(5, [1,2,3.5,6,10]), 2:8)
        @test checkbounds(Bool, BSplineBasis(5, [1,2,3.5,6,10]), 0) === false
        @test checkbounds(Bool, BSplineBasis(5, [1,2,3.5,6,10]), 9) === false
        @test checkbounds(Bool, BSplineBasis(5, [1,2,3.5,6,10]), 2:9) === false
    end

    @testset "Indexing and iteration" begin
        b1 = BSplineBasis(5, 0:5)
        @test eltype(typeof(b1)) === BSpline{BSplineBasis{UnitRange{Int}}}
        @test Base.IteratorEltype(typeof(b1)) === Base.HasEltype()
        @test Base.IteratorSize(typeof(b1)) === Base.HasLength()
        @test collect(b1) == [b1[i] for i=1:9]
        @test @inferred(first(b1)) isa BSpline{BSplineBasis{UnitRange{Int}}}
        @test @inferred(last(b1)) isa BSpline{BSplineBasis{UnitRange{Int}}}
        @test @inferred(b1[1]) isa BSpline{BSplineBasis{UnitRange{Int}}}
        @test first(b1) == Spline(b1, [1,0,0,0,0,0,0,0,0])
        @test last(b1) == Spline(b1, [0,0,0,0,0,0,0,0,1])
        @test b1[5] == Spline(b1, [0,0,0,0,1,0,0,0,0])
        @test b1[end] == Spline(b1, [0,0,0,0,0,0,0,0,1])
        @test_throws BoundsError b1[0]
        @test_throws BoundsError b1[10]
        b2 = BSplineBasis(3, [1,2,3.5,6,10])
        @test eltype(typeof(b2)) === BSpline{BSplineBasis{Vector{Float64}}}
        @test Base.IteratorEltype(typeof(b2)) === Base.HasEltype()
        @test Base.IteratorSize(typeof(b2)) === Base.HasLength()
        @test collect(b2) == [b2[i] for i=1:6]
        @test @inferred(first(b2)) isa BSpline{BSplineBasis{Vector{Float64}}}
        @test @inferred(last(b2)) isa BSpline{BSplineBasis{Vector{Float64}}}
        @test @inferred(b2[1]) isa BSpline{BSplineBasis{Vector{Float64}}}
        @test first(b2) == Spline(b2, [1,0,0,0,0,0])
        @test last(b2) == Spline(b2, [0,0,0,0,0,1])
        @test b2[5] == Spline(b2, [0,0,0,0,1,0])
        @test b2[end] == Spline(b2, [0,0,0,0,0,1])
        @test_throws BoundsError b2[0]
        @test_throws BoundsError b2[7]

        @testset "Iterators.reverse" begin
            b1 = BSplineBasis(5, 0:5)
            rb1 = Iterators.reverse(b1)
            @test eltype(typeof(rb1)) === BSpline{BSplineBasis{UnitRange{Int}}}
            @test Base.IteratorEltype(typeof(rb1)) === Base.HasEltype()
            @test Base.IteratorSize(typeof(rb1)) === Base.HasLength()
            @test collect(rb1) == [b1[i] for i=9:-1:1]
            @test @inferred(first(rb1)) isa BSpline{BSplineBasis{UnitRange{Int}}}
            @test @inferred(last(rb1)) isa BSpline{BSplineBasis{UnitRange{Int}}}
            @test first(rb1) == Spline(b1, [0,0,0,0,0,0,0,0,1])
            @test last(rb1) == Spline(b1, [1,0,0,0,0,0,0,0,0])
            b2 = BSplineBasis(3, [1,2,3.5,6,10])
            rb2 = Iterators.reverse(b2)
            @test eltype(typeof(rb2)) === BSpline{BSplineBasis{Vector{Float64}}}
            @test Base.IteratorEltype(typeof(rb2)) === Base.HasEltype()
            @test Base.IteratorSize(typeof(rb2)) === Base.HasLength()
            @test collect(rb2) == [b2[i] for i=6:-1:1]
            @test @inferred(first(rb2)) isa BSpline{BSplineBasis{Vector{Float64}}}
            @test @inferred(last(rb2)) isa BSpline{BSplineBasis{Vector{Float64}}}
            @test first(rb2) == Spline(b2, [0,0,0,0,0,1])
            @test last(rb2) == Spline(b2, [1,0,0,0,0,0])
        end
    end

    @testset "breakpoints" begin
        @test breakpoints(BSplineBasis(5, 0:5)) == 0:5
        @test breakpoints(BSplineBasis(3, [1,2,3,4,4,5])) == [1,2,3,4,4,5]
    end

    @testset "knots" begin
        @test knots(BSplineBasis(3, 0:5)) == [0, 0, 0, 1, 2, 3, 4, 5, 5, 5]
        @test knots(BSplineBasis(2, [1.5, 4.0])) == [1.5, 1.5, 4.0, 4.0]
    end

    @testset "order" begin
        @test order(BSplineBasis(5, 0:5)) == 5
        @test order(BSplineBasis(3, 0:1//2:10)) == 3
        @test order(BSplineBasis(5, 0:5)[1]) == 5
        @test order(BSplineBasis(3, 0:1//2:10)[2]) == 3
    end

    @testset "Printing" begin
        @test summary(BSplineBasis(5, 0:5)) == "9-element $(BSplineBasis{UnitRange{Int}})"
        @test summary(BSplineBasis(8, [1.0:0.1:3.0;])) == "27-element $(BSplineBasis{Vector{Float64}})"
        @test summary(Spline(BSplineBasis(5, 0:5), zeros(9))) == string(Spline{BSplineBasis{UnitRange{Int}},Vector{Float64}})
        @test summary(Spline(BSplineBasis(8, [1.0:0.1:3.0;]), 1:27)) == string(Spline{BSplineBasis{Vector{Float64}},UnitRange{Int}})
        @test summary(BSpline(BSplineBasis(5, 0:5), 1)) == "BSpline{$(BSplineBasis{UnitRange{Int}})}"
        @test summary(BSpline(BSplineBasis(8, [1.0:0.1:3.0;]), 1)) == "BSpline{$(BSplineBasis{Vector{Float64}})}"
        @test repr(MIME"text/plain"(), BSplineBasis(5, [1,π])) ==
            """
            5-element $(BSplineBasis{Vector{Float64}}):
             order: 5
             breakpoints: [1.0, 3.14159]"""
        @test repr(MIME"text/plain"(), Spline(BSplineBasis(5, [1,π]), sqrt.(0:4))) ==
            """
            $(Spline{BSplineBasis{Vector{Float64}},Vector{Float64}}):
             basis: 5-element $(BSplineBasis{Vector{Float64}}):
              order: 5
              breakpoints: [1.0, 3.14159]
             coeffs: [0.0, 1.0, 1.41421, 1.73205, 2.0]"""
        @test repr(MIME"text/plain"(), BSpline(BSplineBasis(5, [1,π]), 2)) ==
            """
            BSpline{$(BSplineBasis{Vector{Float64}})}:
             basis: 5-element $(BSplineBasis{Vector{Float64}}):
              order: 5
              breakpoints: [1.0, 3.14159]
             index: 2 (knots: [1.0, 1.0, 1.0, 1.0, 3.14159, 3.14159])"""
    end
end

@testset "Spline" begin
    @testset "Construction" begin
        @test_throws DimensionMismatch Spline(BSplineBasis(3, 0:5), zeros(8))
        @test_throws DimensionMismatch Spline(BSplineBasis(3, 0:5), OffsetArray{Int}(undef, 2:8))
        @test_throws DomainError BSpline(BSplineBasis(3, 0:5), 0)
        @test_throws DomainError BSpline(BSplineBasis(3, 0:5), 8)
        @test BSplineBasis(3, 0:5)[4] === BSpline(BSplineBasis(3, 0:5), 4)
    end

    @testset "basis" begin
        @test (b = BSplineBasis(4, [-1, 0, 1.5, 5, 10]); basis(Spline(b, ones(7))) === b)
        @test (b = BSplineBasis(3, 0:5); basis(Spline(b, 1:7)) === b)
        @test (b = BSplineBasis(5, 0:5); basis(Spline(b, zeros(9))) === b)
        @test (b = BSplineBasis(5, 0:5); basis(BSpline(b, 5)) === b)
    end

    @testset "coeffs" begin
        @test (c = ones(7); coeffs(Spline(BSplineBasis(4, [-1, 0, 1.5, 5, 10]), c)) === c)
        @test (c = 1:7; coeffs(Spline(BSplineBasis(3, 0:5), c)) === c)
        @test coeffs(BSpline(BSplineBasis(5, 0:10), 5)) == [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    end

    @testset "order" begin
        @test order(Spline(BSplineBasis(5, 0:5), 1:9)) == 5
        @test order(Spline(BSplineBasis(3, 0:1//2:10), zeros(22))) == 3
        @test order(BSpline(BSplineBasis(4, [1:10;]), 3)) == 4
    end

    @testset "Equality and hashing" begin
        splines = [Spline(BSplineBasis(3, 0:5), Float64[1:7;]),
                   Spline(BSplineBasis(3, Float64[0:5;]), 1:7),
                   BSplineBasis(4, 0:4)[1],
                   BSpline(BSplineBasis(4, 0:4), 1),
                   Spline(BSplineBasis(4, [0,1,2,3,4]), [1,0,0,0,0,0,0]),
                   BSplineBasis(4, 0:5)[1],
                   BSplineBasis(4, Float64[0:5;])[2]]
        for s1 = splines, s2 = splines
            @test (s1 == s2) == (basis(s1) == basis(s2) && coeffs(s1) == coeffs(s2))
            @test (s1 == s2) == (hash(s1) == hash(s2))
        end
    end

    @testset "Broadcasting" begin
        # Splines act as scalars for broadcasting
        @test size(Spline(BSplineBasis(3, 0:5), 1:7) .== BSplineBasis(3, 0:5)[1]) == ()
    end

    @testset "Arithmetic" begin
        @testset "Unary +/-" begin
            @test +(BSplineBasis(3, 0:5)[4]) == Spline(BSplineBasis(3, 0:5), [0,0,0,1,0,0,0])
            @test +(Spline(BSplineBasis(4, [1:10;]), 1:12)) == Spline(BSplineBasis(4, 1:10), 1:12)
            @test +(Spline(BSplineBasis(4, [1:10;]), [1:12;])) == Spline(BSplineBasis(4, 1:10), 1:12)
            @test -(BSplineBasis(3, 0:5)[4]) == Spline(BSplineBasis(3, 0:5), [0,0,0,-1,0,0,0])
            @test -(Spline(BSplineBasis(4, [1:10;]), 1:12)) == Spline(BSplineBasis(4, 1:10), -1:-1:-12)
            @test -(Spline(BSplineBasis(4, [1:10;]), [1:12;])) == Spline(BSplineBasis(4, 1:10), -1:-1:-12)
        end

        @testset "Binary +/-" begin
            # Addition
            @test (BSplineBasis(3, 0:5)[1]) + (BSplineBasis(3, 0:5)[4]) == Spline(BSplineBasis(3, 0:5), [1,0,0,1,0,0,0])
            @test (BSplineBasis(3, 0:5)[2]) + (BSplineBasis(3, [0:5;])[2]) == Spline(BSplineBasis(3, 0:5), [0,2,0,0,0,0,0])
            @test (BSplineBasis(4, [1:10;])[1]) + Spline(BSplineBasis(4, 1:10), ones(12)) ==
                Spline(BSplineBasis(4, 1:10), [2,1,1,1,1,1,1,1,1,1,1,1])
            @test Spline(BSplineBasis(4, [1:10;]), -11:0) + (BSplineBasis(4, [1:10;])[12]) ==
                Spline(BSplineBasis(4, 1:10), [-11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,1])
            @test Spline(BSplineBasis(5, -6:2:6), sqrt.(1:10)) + Spline(BSplineBasis(5, -6:2:6), log.(1:10)) ==
                Spline(BSplineBasis(5, -6:2:6), [sqrt(i)+log(i) for i=1:10])
            @test Spline(BSplineBasis(5, -3:1:3), 1.0:10.0) + Spline(BSplineBasis(5, -3:3), -10.0:-1.0) ==
                Spline(BSplineBasis(5, -3:3), -9.0:2.0:9.0)
            # Subtraction
            @test (BSplineBasis(3, 0:5)[1]) - (BSplineBasis(3, 0:5)[4]) == Spline(BSplineBasis(3, 0:5), [1,0,0,-1,0,0,0])
            @test (BSplineBasis(3, 0:5)[2]) - (BSplineBasis(3, [0:5;])[2]) == Spline(BSplineBasis(3, 0:5), zeros(7))
            @test (BSplineBasis(4, [1:10;])[1]) - Spline(BSplineBasis(4, 1:10), ones(12)) ==
                Spline(BSplineBasis(4, 1:10), [0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1])
            @test Spline(BSplineBasis(4, [1:10;]), -11:0) - (BSplineBasis(4, [1:10;])[12]) ==
                Spline(BSplineBasis(4, 1:10), [-11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,-1])
            @test Spline(BSplineBasis(5, -6:2:6), sqrt.(1:10)) - Spline(BSplineBasis(5, -6:2:6), log.(1:10)) ==
                Spline(BSplineBasis(5, -6:2:6), [sqrt(i)-log(i) for i=1:10])
            @test Spline(BSplineBasis(5, -3:1:3), 1.0:10.0) - Spline(BSplineBasis(5, -3:3), -10.0:-1.0) ==
                Spline(BSplineBasis(5, -3:3), fill(11.0, 10))
            # Different bases
            @test_throws ArgumentError (BSplineBasis(3, 0:5)[2]) + (BSplineBasis(4, 0:4)[2])
            @test_throws ArgumentError (BSplineBasis(3, 0:5)[2]) + Spline(BSplineBasis(4, 0:5), ones(8))
            @test_throws ArgumentError Spline(BSplineBasis(3, 0:5), zeros(7)) + Spline(BSplineBasis(4, 1:10), ones(12))
            @test_throws ArgumentError (BSplineBasis(3, 0:5)[2]) - (BSplineBasis(4, 0:4)[2])
            @test_throws ArgumentError (BSplineBasis(3, 0:5)[2]) - Spline(BSplineBasis(4, 0:5), ones(8))
            @test_throws ArgumentError Spline(BSplineBasis(3, 0:5), zeros(7)) - Spline(BSplineBasis(4, 1:10), ones(12))
        end

        @testset "Multiplication/division" begin
            # Base.*
            @test (BSplineBasis(3, 0:5)[5]) * 2 == Spline(BSplineBasis(3, 0:5), [0,0,0,0,2,0,0])
            @test Spline(BSplineBasis(4, [1:10;]), [sqrt(i) for i=1:12]) * 0.1 ==
                Spline(BSplineBasis(4, [1:10;]), [sqrt(i)*0.1 for i=1:12])
            @test 3.0 * (BSplineBasis(3, 0:5)[5])  == Spline(BSplineBasis(3, 0:5), [0,0,0,0,3,0,0])
            @test (5//2) * Spline(BSplineBasis(4, [1:10;]), [sqrt(i) for i=1:12]) ==
                Spline(BSplineBasis(4, [1:10;]), [(5//2)*sqrt(i) for i=1:12])
            # Base./
            @test (BSplineBasis(3, 0:5)[5]) / 2 == Spline(BSplineBasis(3, 0:5), [0,0,0,0,.5,0,0])
            @test Spline(BSplineBasis(4, [1:10;]), [sqrt(i) for i=1:12]) / 3.0 ==
                Spline(BSplineBasis(4, [1:10;]), [sqrt(i)/3.0 for i=1:12])
            # Base.\
            @test 10 \ (BSplineBasis(3, 0:5)[5])  == Spline(BSplineBasis(3, 0:5), [0,0,0,0,1/10,0,0])
            @test (5//2) \ Spline(BSplineBasis(4, [1:10;]), [sqrt(i) for i=1:12]) ==
                Spline(BSplineBasis(4, [1:10;]), [(5//2)\sqrt(i) for i=1:12])
            # LinearAlgebra.lmul!
            @test begin
                spl = Spline(BSplineBasis(4, [1:10;]), [sqrt(i) for i=1:12])
                lmul!(5//2, spl) === spl && spl == Spline(BSplineBasis(4, [1:10;]), [(5//2)*sqrt(i) for i=1:12])
            end
            # LinearAlgebra.rmul!
            @test begin
                spl = Spline(BSplineBasis(3, 0:5), [0,0,0,0,1,0,0])
                rmul!(spl, 2) === spl && spl == Spline(BSplineBasis(3, 0:5), [0,0,0,0,2,0,0])
            end
            @static if VERSION ≥ v"1.2"
                # LinearAlgebra.ldiv!
                @test begin
                    spl = Spline(BSplineBasis(3, 0:5), [0,0,0,0,1//1,0,0])
                    ldiv!(10, spl) === spl && spl == Spline(BSplineBasis(3, 0:5), [0,0,0,0,1//10,0,0])
                end
                # LinearAlgebra.rdiv!
                @test begin
                    spl = Spline(BSplineBasis(4, [1:10;]), [sqrt(i) for i=1:12])
                    rdiv!(spl, 3.0) === spl && spl ==  Spline(BSplineBasis(4, [1:10;]), [sqrt(i)/3.0 for i=1:12])
                end
            end
        end
    end
end

@testset "Derivatives" begin
    # Derivative
    @test Derivative{1}() isa Derivative
    @test_throws DomainError Derivative{-1}()
    @test_throws ArgumentError Derivative{UInt(1)}()
    @test Derivative(2) === Derivative{2}()
    @test Derivative(Int8(3)) === Derivative{3}()
    @test Derivative(UInt(4)) === Derivative{4}()
    @test_throws DomainError Derivative(-1)
    @test_throws DomainError Derivative(Int8(-2))
    @test Derivative(2) == Derivative(2)
    @test Derivative(2) != Derivative(3)
    # AllDerivatives
    @test AllDerivatives{2}() isa AllDerivatives
    @test_throws DomainError AllDerivatives{0}()
    @test_throws DomainError AllDerivatives{-1}()
    @test_throws ArgumentError AllDerivatives{UInt(1)}()
    @test AllDerivatives(2) === AllDerivatives{2}()
    @test AllDerivatives(Int8(3)) === AllDerivatives{3}()
    @test AllDerivatives(UInt(4)) === AllDerivatives{4}()
    @test_throws DomainError AllDerivatives(0)
    @test_throws DomainError AllDerivatives(-1)
    @test_throws DomainError AllDerivatives(Int8(-2))
    @test_throws DomainError AllDerivatives(UInt(0))
    @test AllDerivatives(2) == AllDerivatives(2)
    @test AllDerivatives(2) != AllDerivatives(3)
    # Derivative/AllDerivatives act as scalars for broadcasting
    @test size(Derivative(1) .== Derivative(1)) == ()
    @test size(AllDerivatives(4) .== AllDerivatives(4)) == ()
end

include("bases.jl")
include("bsplines.jl")
include("splinevalue.jl")

@testset "SplineFunctionWrapper" begin
    function test_splinefunctionwrapper(spline)
        a,b = support(spline)
        f  = Function(spline)
        f⁰ = Function(spline, Derivative(0))
        f¹ = Function(spline, Derivative(1))
        f² = Function(spline, Derivative(2))
        g  = Function(spline, false)
        g⁰ = Function(spline, Derivative(0), false)
        g¹ = Function(spline, Derivative(1), false)
        g² = Function(spline, Derivative(2), false)
        xs = range(a, stop=b, length=10)
        @test all(f.(xs) .== spline.(xs))
        @test all(f⁰.(xs) .== spline.(xs))
        @test all(f¹.(xs) .== spline.(xs, Derivative(1)))
        @test all(f².(xs) .== spline.(xs, Derivative(2)))
        @test all(g.(xs) .== spline.(xs))
        @test all(g⁰.(xs) .== spline.(xs))
        @test all(g¹.(xs) .== spline.(xs, Derivative(1)))
        @test all(g².(xs) .== spline.(xs, Derivative(2)))
        for x = (a-1, b+1)
            for fun = (f, f⁰, f¹, f²)
                @test isnan(fun(x))
            end
            for fun = (g, g⁰, g¹, g²)
                @test iszero(fun(x))
            end
        end
    end
    for spline = [BSpline(BSplineBasis(4, 0:2:20), 8),
                  Spline(BSplineBasis(3, 0//1:5//1), 1:7),
                  Spline(BSplineBasis(5, [0:0.1:1;]), ones(14))]
        test_splinefunctionwrapper(spline)
    end
end

@testset "basismatrix" begin
    basis = BSplineBasis(3, 0:5)
    xvalues = [0.3, 1.5, 3.2]
    @test basismatrix(basis, xvalues) == [basis[i](x) for x=xvalues, i=eachindex(basis)]
    @test basismatrix(basis, xvalues, indices=2:5) == [basis[i](x) for x=xvalues, i=2:5]
    @test_throws ArgumentError basismatrix(basis, xvalues, indices=0:5)
end

@testset "basismatrix!" begin
    basis = BSplineBasis(3, 0:5)
    xvalues = [0.3, 1.5, 3.2]
    @test begin
        dest = Matrix{Float64}(undef, 3, length(basis))
        basismatrix!(dest, basis, xvalues)
        dest == [basis[i](x) for x=xvalues, i=eachindex(basis)]
    end
    @test begin
        dest = Matrix{Float64}(undef, 3, 4)
        basismatrix!(dest, basis, xvalues, indices=2:5)
        dest == [basis[i](x) for x=xvalues, i=2:5]
    end
    @test_throws ArgumentError basismatrix!(Array{Float64}(undef, 3, 6), basis, xvalues, indices=0:5)
    @test_throws ArgumentError basismatrix!(OffsetArray(Array{Float64}(undef, 3, length(basis)),1,1), basis, xvalues)
    @test_throws ArgumentError basismatrix!(OffsetArray(Array{Float64}(undef, 3, 4),1,1), basis, xvalues, indices=2:5)
    @test_throws DimensionMismatch basismatrix!(Array{Float64}(undef, 3, 4), basis, xvalues)
    @test_throws DimensionMismatch basismatrix!(Array{Float64}(undef, 2, 4), basis, xvalues, indices=2:5)
end

@testset "approximate" begin
    basis = BSplineBasis(3, 0:5)
    knotavgs = [0.0, 0.5, 1.5, 2.5, 3.5, 4.5, 5.0]
    @test begin
        spline = @inferred(approximate(sqrt, basis))
        all(spline(x) ≈ sqrt(x) for x=knotavgs)
    end
    @test @inferred(approximate(sqrt, basis, indices=:)) == approximate(sqrt, basis)
    @test @inferred(approximate(sqrt, basis, indices=1:7)) == approximate(sqrt, basis)
    @test_throws ArgumentError approximate(sqrt, basis, indices=0:7)

    @test begin
        spline = approximate(sqrt, basis, indices=2:7)
        all(spline(x) ≈ sqrt(x) for x=knotavgs) && iszero(coeffs(spline)[1])
    end

    basis = BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])
    knotavgs = [-5.0, -4.625, -3.925, -2.425, -0.675, 1.1, 2.55, 3.2, 3.6]
    # basis = BSplineBasis(5, [-5, -7//2, -11//5, 1, 2, 18//5])
    # knotavgs = [-5, -37//8, -157//40, -97//40, -27//40, 11//10, 51//20, 16//5, 18//5]
    @test begin
        spline = @inferred(approximate(sin, basis))
        all(spline(x) ≈ sin(x) for x = knotavgs)
    end
    @test @inferred(approximate(sin, basis, indices=:)) == approximate(sin, basis)
    @test @inferred(approximate(sin, basis, indices=1:9)) == approximate(sin, basis)
    @test_throws ArgumentError approximate(sin, basis, indices=2:10)

    @test begin
        spline = approximate(one, BSplineBasis(3, 0:5))
        all(spline(x) ≈ 1 for x=range(0, stop=5, length=100))
    end
    @test begin
        basis = BSplineBasis(9, 0:0.1:1)
        spline = approximate(sinpi, basis, indices=2:length(basis)-1)
        all(spline(x) ≈ sinpi(x) for x=range(0, stop=1, length=100))
    end
end

@testset "interpolate" begin
    basis = BSplineBasis(3, 0:5)
    xvalues = [0.0, 0.5, 1.5, 2.5, 3.5, 4.5, 5.0]
    yvalues = sqrt.(xvalues)
    @test begin
        spline = @inferred(interpolate(basis, xvalues, yvalues))
        all(spline(xvalues[i]) ≈ yvalues[i] for i=eachindex(xvalues))
    end
    @test @inferred(interpolate(basis, xvalues, yvalues, indices=:)) == interpolate(basis, xvalues, yvalues)
    @test @inferred(interpolate(basis, xvalues, yvalues, indices=1:7)) == interpolate(basis, xvalues, yvalues)
    @test_throws ArgumentError interpolate(basis, xvalues, yvalues, indices=0:7)
    @test_throws DimensionMismatch interpolate(basis, xvalues, [yvalues; 0])
    @test_throws DimensionMismatch interpolate(basis, [xvalues; 0], yvalues, indices=:)
    @test_throws DimensionMismatch interpolate(basis, [xvalues; 0], yvalues, indices=1:7)

    @test begin
        spline = interpolate(basis, xvalues, yvalues, indices=2:7)
        all(spline(xvalues[i]) ≈ yvalues[i] for i=eachindex(xvalues)) && iszero(coeffs(spline)[1])
    end

    xvalues = xvalues[2:end]
    yvalues = yvalues[2:end]
    @test begin
        spline = interpolate(basis, xvalues, yvalues, indices=2:7)
        all(spline(xvalues[i]) ≈ yvalues[i] for i=eachindex(xvalues)) && iszero(coeffs(spline)[1])
    end

    basis = BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])
    xvalues = [-5.0, -4.625, -3.925, -2.425, -0.675, 1.1, 2.55, 3.2, 3.6]
    # basis = BSplineBasis(5, [-5, -7//2, -11//5, 1, 2, 18//5])
    # xvalues = [-5, -37//8, -157//40, -97//40, -27//40, 11//10, 51//20, 16//5, 18//5]
    yvalues = sin.(xvalues)
    @test begin
        spline = @inferred(interpolate(basis, xvalues, yvalues))
        all(spline(xvalues[i]) ≈ yvalues[i] for i=eachindex(xvalues))
    end
    @test @inferred(interpolate(basis, xvalues, yvalues, indices=:)) == interpolate(basis, xvalues, yvalues)
    @test @inferred(interpolate(basis, xvalues, yvalues, indices=1:9)) == interpolate(basis, xvalues, yvalues)
    @test_throws ArgumentError interpolate(basis, xvalues, yvalues, indices=2:10)
    @test_throws DimensionMismatch interpolate(basis, xvalues, [yvalues; 0])
    @test_throws DimensionMismatch interpolate(basis, [xvalues; 0], yvalues, indices=:)
    @test_throws DimensionMismatch interpolate(basis, [xvalues; 0], yvalues, indices=1:9)

    yvalues[end] = 0
    @test begin
        spline = interpolate(basis, xvalues, yvalues, indices=1:8)
        all(spline(xvalues[i]) ≈ yvalues[i] for i=eachindex(xvalues)) && iszero(coeffs(spline)[end])
    end

    xvalues = xvalues[1:end-1]
    yvalues = yvalues[1:end-1]
    @test begin
        spline = interpolate(basis, xvalues, yvalues, indices=1:8)
        all(spline(xvalues[i]) ≈ yvalues[i] for i=eachindex(xvalues)) && iszero(coeffs(spline)[end])
    end
end

@testset "averagebasis" begin
    ≈ₜₑₛₜ(x::BSplineBasis, y::BSplineBasis; kwargs...) =
        order(x) == order(y) && isapprox(breakpoints(x), breakpoints(y); kwargs...)

    @test averagebasis(3, [0.0, 0.5, 1.5, 2.5, 3.5, 4.5, 5.0]) ≈ₜₑₛₜ BSplineBasis(3, 0:5)
    @test averagebasis(3, [0//1, 1//2, 3//2, 5//2, 7//2, 9//2, 5//1]) == BSplineBasis(3, 0:5)
    @test averagebasis(4, [0.0, 1/3, 1.0, 2.0, 3.0, 4.0, 14/3, 5.0]) ≈ₜₑₛₜ
        BSplineBasis(4, [0.0, 10/9, 2.0, 3.0, 35/9, 5])
    @test averagebasis(4, [0//1, 1//3, 1//1, 2//1, 3//1, 4//1, 14//3, 5//1]) ==
        BSplineBasis(4, [0//1, 10//9, 2//1, 3//1, 35//9, 5//1])
    @test averagebasis(5, 1:10) ≈ₜₑₛₜ BSplineBasis(5, [1.0, 3.5, 4.5, 5.5, 6.5, 7.5, 10.0])
    @test averagebasis(5, 1//1:10//1) == BSplineBasis(5, [1//1, 7//2, 9//2, 11//2, 13//2, 15//2, 10//1])
    @test averagebasis(6, 1:1:10) ≈ₜₑₛₜ BSplineBasis(6, [1, 4, 5, 6, 7, 10])
    @test averagebasis(6, 1//1:1//1:10//1) == BSplineBasis(6, [1, 4, 5, 6, 7, 10])
    @test averagebasis(7, 1:7) ≈ₜₑₛₜ BSplineBasis(7, [1, 7])
    @test averagebasis(7, 1//1:7//1) == BSplineBasis(7, [1, 7])
    @test_throws ArgumentError averagebasis(8, 1:7) # too few datapoints
end

@testset "knotaverages" begin
    @test knotaverages(BSplineBasis(3, 0//1:5//1))              == [0, 1//2, 3//2, 5//2, 7//2, 9//2, 5]
    @test knotaverages(BSplineBasis(3, 0//1:5//1), indices=:)   == [0, 1//2, 3//2, 5//2, 7//2, 9//2, 5]
    @test knotaverages(BSplineBasis(3, 0//1:5//1), indices=2:6) == [1//2, 3//2, 5//2, 7//2, 9//2]
    @test_throws ArgumentError knotaverages(BSplineBasis(3, 0//1:5//1), indices=0:3)
    @test_throws ArgumentError knotaverages(BSplineBasis(3, 0//1:5//1), indices=1:8)
    @test knotaverages(BSplineBasis(4, 0:1:5))              ≈ [0, 1/3, 1, 2, 3, 4, 14/3, 5]
    @test knotaverages(BSplineBasis(4, 0:1:5), indices=:)   ≈ [0, 1/3, 1, 2, 3, 4, 14/3, 5]
    @test knotaverages(BSplineBasis(4, 0:1:5), indices=2:8) ≈ [1/3, 1, 2, 3, 4, 14/3, 5]
    @test knotaverages(BSplineBasis(4, 0.0:5.0))              ≈ [0, 1/3, 1, 2, 3, 4, 14/3, 5]
    @test knotaverages(BSplineBasis(4, 0.0:5.0), indices=:)   ≈ [0, 1/3, 1, 2, 3, 4, 14/3, 5]
    @test knotaverages(BSplineBasis(4, 0.0:5.0), indices=2:8) ≈ [1/3, 1, 2, 3, 4, 14/3, 5]
    @test knotaverages(BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6])) ≈
        [-5.0, -4.625, -3.925, -2.425, -0.675, 1.1, 2.55, 3.2, 3.6]
    @test knotaverages(BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6]), indices=:) ≈
        [-5.0, -4.625, -3.925, -2.425, -0.675, 1.1, 2.55, 3.2, 3.6]
    @test knotaverages(BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6]), indices=1:8) ≈
        [-5.0, -4.625, -3.925, -2.425, -0.675, 1.1, 2.55, 3.2]
    @test knotaverages(BSplineBasis(5, [-5//1, -7//2, -11//5, 1//1, 2//1, 18//5])) ==
        [-5//1, -37//8, -157//40, -97//40, -27//40, 11//10, 51//20, 16//5, 18//5]
    @test knotaverages(BSplineBasis(5, [-5//1, -7//2, -11//5, 1//1, 2//1, 18//5]), indices=:) ==
        [-5//1, -37//8, -157//40, -97//40, -27//40, 11//10, 51//20, 16//5, 18//5]
    @test knotaverages(BSplineBasis(5, [-5//1, -7//2, -11//5, 1//1, 2//1, 18//5]), indices=1:8) ==
        [-5//1, -37//8, -157//40, -97//40, -27//40, 11//10, 51//20, 16//5]
end

@testset "knotaverages!" begin
    @test begin
        dest = Vector{Float64}(undef, 7)
        knotaverages!(dest, BSplineBasis(3, 0//1:5//1))
        dest ≈ [0, 1/2, 3/2, 5/2, 7/2, 9/2, 5]
    end
    @test begin
        dest = Vector{Float64}(undef, 7)
        knotaverages!(dest, BSplineBasis(3, 0//1:5//1))
        dest ≈ [0, 1/2, 3/2, 5/2, 7/2, 9/2, 5]
    end
    @test begin
        dest = Vector{Rational{Int}}(undef, 5)
        knotaverages!(dest, BSplineBasis(3, 0//1:5//1), indices=2:6)
        dest == [1//2, 3//2, 5//2, 7//2, 9//2]
    end
    @test_throws ArgumentError knotaverages!(zeros(4), BSplineBasis(3, 0//1:5//1), indices=0:3)
    @test_throws ArgumentError knotaverages!(zeros(8), BSplineBasis(3, 0//1:5//1), indices=1:8)
    @test_throws ArgumentError knotaverages!(OffsetArray{Float64}(undef, 2:8), BSplineBasis(3, 0//1:5//1))
    @test_throws ArgumentError knotaverages!(OffsetArray{Float64}(undef, 2:8), BSplineBasis(3, 0//1:5//1), indices=:)
    @test_throws ArgumentError knotaverages!(OffsetArray{Float64}(undef, 0:5), BSplineBasis(3, 0//1:5//1), indices=2:7)
    @test_throws DimensionMismatch knotaverages!(zeros(6), BSplineBasis(3, 0//1:5//1))
    @test_throws DimensionMismatch knotaverages!(zeros(6), BSplineBasis(3, 0//1:5//1), indices=:)
    @test_throws DimensionMismatch knotaverages!(zeros(7), BSplineBasis(3, 0//1:5//1), indices=2:7)
    @test begin
        dest = Vector{Float32}(undef, 8)
        knotaverages!(dest, BSplineBasis(4, 0:1:5))
        dest ≈ [0, 1/3, 1, 2, 3, 4, 14/3, 5]
    end
    @test begin
        dest = Vector{Float64}(undef, 8)
        knotaverages!(dest, BSplineBasis(4, 0:1:5), indices=:)
        dest ≈ [0, 1/3, 1, 2, 3, 4, 14/3, 5]
    end
    @test begin
        dest = Vector{Rational{Int}}(undef, 5)
        knotaverages!(dest, BSplineBasis(4, 0:1:5), indices=2:6)
        dest == [1//3, 1, 2, 3, 4]
    end
    @test begin
        dest = Vector{Float32}(undef, 9)
        knotaverages!(dest, BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6]))
        dest ≈ [-5.0, -4.625, -3.925, -2.425, -0.675, 1.1, 2.55, 3.2, 3.6]
    end
    @test begin
        dest = Vector{Float64}(undef, 9)
        knotaverages!(dest, BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6]), indices=:)
        dest ≈ [-5.0, -4.625, -3.925, -2.425, -0.675, 1.1, 2.55, 3.2, 3.6]
    end
    @test begin
        dest = Vector{BigFloat}(undef, 8)
        knotaverages!(dest, BSplineBasis(5, [-5.0, -3.5, -2.2, 1.0, 2.0, 3.6]), indices=1:8)
        dest ≈ [-5.0, -4.625, -3.925, -2.425, -0.675, 1.1, 2.55, 3.2]
    end
end


text:
for N in [Float64, Rational{Int}, Float32]

    # ==================================
    # Constructor and interface methods
    # ==================================

    B = BallInf(zeros(N, 3), N(1))
    v = N[1, 0, 0] # translation along dimension 1
    M = Diagonal(N[1, 2, 3])
    am = AffineMap(M, B, v)

    # dimension check
    @test dim(am) == 3

    # dimension assertion
    @test_throws AssertionError AffineMap(M, B, N[0, 0])
    @test_throws AssertionError AffineMap(M, B × B, v)

    # linear map of an affine map is automatically simplified to an affine map
    Mam = M * am
    @test Mam isa AffineMap && Mam.M == M * am.M && Mam.v == M * am.v

    # scaling of an affine map is an affine map
    α = N(2)
    αam = α * am
    @test αam isa AffineMap && αam.M == α * am.M && αam.v == α * am.v

    # support vector
    @test σ(N[1, 0, 0], am) == N[2, 2, 3]

    # support function
    @test ρ(N[1, 0, 0], am) == N(2)

    # boundedness
    @test isbounded(am)
    @test !isbounded(AffineMap(M, Universe{N}(3), v))

    # function to get an element
    @test (an_element(am) - am.v) ∈ (am.M * am.X)
    @test an_element(am) ∈ am.M * am.X ⊕ am.v

    # emptiness check
    @test !isempty(am)
    # @test isempty(AffineMap(M, EmptySet{N}(), v)) # requires #1201

    # ==================================
    # Type-specific methods
    # ==================================

    # an affine map of the form I*X + b where I is the identity matrix is a pure translation
    #v = N[1, 0, 2]
    #am_tr = AffineMap(I, B, v) # crashes, see #1544
    #@test am_tr isa Translation && am_tr.v == v

    # two-dimensional case
    B2 = BallInf(zeros(N, 2), N(1))
    M = N[1 0; 0 2]
    v = N[-1, 0]
    am = AffineMap(M, B2, v)

    # list of vertices check
    vlist = vertices_list(am)
    @test ispermutation(vlist, [N[0, 2], N[-2, 2], N[0, -2], N[-2, -2]])

    # inclusion check
    h = Hyperrectangle(N[-1, 0], N[1, 2])
    @test h ⊆ am && am ⊆ h
end

# tests that only work with Float64 and Float32
for N in [Float64, Float32]
    B = BallInf(zeros(N, 3), N(1))

    # the translation is the origin and the linear map is the identity => constraints remain unchanged
    Id3 = Matrix(one(N) * I, 3, 3)
    @test ispermutation(constraints_list(AffineMap(Id3, B, zeros(N, 3))),
                        constraints_list(B))
end


text:
export Technology,
Timestep,
Variable,
Carrier,
Region,
Scenario,
Mask,
PivotResult,
pivotresult,
StackedResult,
stackedresult

const colnames = (
    Region = "region_dispatch_",
    Technology = "technology_",
    Timestep = "timestep_superordinate_dispatch_",
    Carrier = "carrier_",
    Variable = "variable",
    Scenario = "scenario"
)

abstract type ResultDimension{T} end
abstract type ResultWithDimension{T} <: ResultDimension{T} end
abstract type ResultWithoutDimension{T} <: ResultDimension{T} end
const SuperType = Union{AbstractVector, AbstractString, Nothing, Int}

for dimension = (:Technology, :Timestep, :Carrier, :Region)
    eval(quote
        struct $dimension{T<:SuperType} <: ResultWithDimension{T}
            vals::T
            mode::Symbol
            dim::Int
            colname::AbstractString

            function $dimension()
                colname = colnames.$dimension * "1"
                return new{Nothing}(nothing, :equal, 1, colname)
            end

            function $dimension(p::T; mode = :in) where T<:Pair{<:Int,<:AbstractVector}
                dim, val = p
                colname = colnames.$dimension * "$dim"
                return new{typeof(val)}(val, mode, dim, colname)
            end

            function $dimension(p::T; mode = :equal) where T<:Pair{<:Int,<:AbstractString}
                dim, val = p
                colname = colnames.$dimension * "$dim"
                int_val = tryparse(Int, val)
                if isnothing(int_val)
                    return new{typeof(val)}(val, mode, dim, colname)
                else
                    return new{typeof(int_val)}(int_val, mode, dim, colname)
                end
            end

            function $dimension(p::T; mode = :equal) where T<:Pair{<:Int,<:Int}
                dim, val = p
                colname = colnames.$dimension * "$dim"
                return new{typeof(val)}(val, mode, dim, colname)
            end

            function $dimension(x::T; mode = :equal, dim=1) where T<:AbstractString
                colname = colnames.$dimension * "$dim"
                int_val = tryparse(Int, x)
                if isnothing(int_val)
                    return new{T}(x, mode, dim, colname)
                else
                    return new{typeof(int_val)}(int_val, mode, dim, colname)
                end
            end

            function $dimension(i::T; mode = :equal) where T<:Integer
                colname = colnames.$dimension * "$i"
                return new{Nothing}(nothing, mode, i, colname)
            end

            function $dimension(vec::T; mode = :in, dim=1) where T<:AbstractVector
                colname = colnames.$dimension * "$dim"
                return new{T}(vec, mode, dim, colname)
            end

            function $dimension(args::T...; mode = :in, dim=1) where T<:AbstractString
                colname = colnames.$dimension * "$dim"
                args_int = tryparse.(Int, args)
                @show args_int
                if eltype(args_int) <: Int
                    vec = [args_int...]
                else
                    vec = [args...]
                end
                return new{typeof(vec)}(vec, mode, dim, colname)
            end
        end
    end)
end

for dimension = (:Variable, :Scenario)
    eval(quote
        struct $dimension{T<:SuperType} <: ResultWithoutDimension{T}
            vals::T
            mode::Symbol
            colname::AbstractString

            $dimension() = new{Nothing}(nothing, :equal, colnames.$dimension)

            function $dimension(x::T; mode = :equal) where T<:AbstractString
                colname = colnames.$dimension
                return new{T}(x, mode, colname)
            end

            function $dimension(vec::T; mode = :in) where T<:AbstractVector
                colname = colnames.$dimension
                return new{T}(vec, mode, colname)
            end

            function $dimension(args::T...; mode = :in) where T<:AbstractString
                colname = colnames.$dimension
                vec = [args...]
                return new{typeof(vec)}(vec, mode, colname)
            end
        end
    end)
end

struct Mask
    row
    col
    filters::Vector{ResultDimension}

    function Mask(row, col::Nothing, filters::ResultDimension...)
        new(row, col, [filters...])
    end

    function Mask(row, col::ResultDimension, filters::ResultDimension...)
        new(row, col, [filters...])
    end
end

function filter_df!(df::AbstractDataFrame, af::ResultDimension{Nothing})
    return df
end

function filter_df!(df::AbstractDataFrame, af::ResultDimension{T}) where T<:AbstractString

    if af.colname in names(df)
        if af.mode == :equal
            filter!(x-> isequal(af.vals, x[af.colname]), df)
        elseif af.mode == :occursin
            filter!(x-> occursin(af.vals, x[af.colname]), df)
        else
            @warn "Mode $(af.mode) is not supported for filtering based on a string" af.vals
        end
    else
        @warn "Column $(af.colname) does not exists in DataFrame" names(df)
    end

    return df
end

function filter_df!(df::AbstractDataFrame, af::ResultDimension{T}) where T<:Int

    if af.colname in names(df)
        if af.mode == :equal
            coltype = eltype(df[:,af.colname])
            if coltype <: Int
                filter!(x-> isequal(af.vals, x[af.colname]), df)
            elseif coltype <: AbstractString
                filter!(x-> string(af.vals) == x[af.colname], df)
            else
                @warn "Type $coltype in column $(af.colname) could not be filtered!"
            end
        else
            @warn "Mode $(af.mode) is not supported for filtering based on a string" af.vals
        end
    else
        @warn "Column $(af.colname) does not exists in DataFrame" names(df)
    end

    return df
end

function filter_df!(df::AbstractDataFrame, af::ResultDimension{T}) where T<:AbstractVector

    if af.colname in names(df)
        if af.mode == :in
            filter!(x-> x[af.colname] in af.vals, df)
        else
            @warn "Mode $(af.mode) is not supported for filtering based on a vector" af.vals
        end
    else
        @warn "Column $(af.colname) does not exists in DataFrame" names(df)
    end

    return df
end

function filter_df!(df::AbstractDataFrame, vec::T) where T<:AbstractVector
    for fil in vec
        filter_df!(df, fil)
    end
    return df
end

function filter_df!(df::AbstractDataFrame, vec::Nothing)
    return df
end

function filter_results(ar::AnymodResult, m::Mask)
    df = copy(ar.summarytable)
    # filter summarytable
    filter_df!(df, m.filters)
    filter_df!(df, m.row)
    filter_df!(df, m.col)

    return df
end

function get_colnames(row::T, col) where T<:AbstractArray
    colnames = [rt.colname for rt in row]
    groupkeys = [colnames..., col.colname]
    return colnames, groupkeys
end

function get_colnames(row::T, col::Nothing) where T<:AbstractArray
    colnames = groupkeys = [rt.colname for rt in row]
    return colnames, groupkeys
end

function get_colnames(row, col::Nothing)
    return row.colname, row.colname
end

function get_colnames(row, col)
    colnames = row.colname
    groupkeys = [row.colname, col.colname]
    return colnames, groupkeys
end

function prepare_df_result(ar::AnymodResult, m::Mask)
    df = filter_results(ar, m)
    # to do: handle case when df is empty
    colnames, groupkeys = get_colnames(m.row, m.col)
    df = combine(groupby(df, groupkeys), "value" => sum => "value")
    return colnames, df
end

struct PivotResult
    result::AnymodResult
    mask::Mask
    df::AbstractDataFrame

    function PivotResult(ar::AnymodResult, m::Mask)
        colnames, df = prepare_df_result(ar, m)
        if isnothing(m.col)
            return new(ar, m, df)
        else
            df = unstack(df, colnames, m.col.colname, "value")
        end
        return new(ar, m, df)
    end
end

function pivotresult(ar::AnymodResult, m::Mask)
    return PivotResult(ar, m).df
end

struct StackedResult
    result::AnymodResult
    mask::Mask
    df::AbstractDataFrame

    function StackedResult(ar::AnymodResult, m::Mask)
        _, df = prepare_df_result(ar, m)
        return new(ar, m, df)
    end
end

function stackedresult(ar::AnymodResult, m::Mask)
    return StackedResult(ar, m).df
end


text:
struct Pos
    x::Float64
    y::Float64
end

function jarvis_cross(point1::Pos, point2::Pos, point3::Pos)
    vec1 = Pos(point2.x - point1.x, point2.y - point1.y)
    vec2 = Pos(point3.x - point2.x, point3.y - point2.y)
    ret_cross = vec1.x*vec2.y - vec1.y*vec2.x
    return ret_cross*ret_cross
end

function jarvis_march(points::Vector{Pos})
    hull = Vector{Pos}()

    # sorting array based on leftmost point
    sort!(points, by = item -> item.x)
    push!(hull, points[1])

    i = 1
    curr_point = points[2]

    # Find cross product between points
    curr_product = jarvis_cross(Pos(0,0), hull[1], curr_point)
    while (curr_point != hull[1])
        for point in points
                product = 0.0
            if (i == 1)
                if (hull[i] != point)
                    product = jarvis_cross(Pos(0,0), hull[i], point)
                end
            else
                if (hull[i] != point && hull[i-1] != point)
                    product = jarvis_cross(hull[i-1], hull[i], point)
                end
            end
            if (product > curr_product)
                curr_point = point
                curr_product = product
            end
        end
        push!(hull, curr_point)
        curr_product = 0
        i += 1
    end

    return hull
end

function main()

    points = [Pos(2,1.5), Pos(1, 1), Pos(2, 4), Pos(3, 1)]
    hull = jarvis_march(points)
    println(hull)
end

main()

max_stars_repo_name: atocil/algorithm-archive

max_issues_repo_name: atocil/algorithm-archive

max_forks_repo_name: atocil/algorithm-archive


text:
############################
### Base Issue Functions ###
############################

listIssue(c::Client, page::Int64, state::String = "") = getParsedResponse(Vector{Issue},c,Requests.get,"/issues?page=$(page)")

listUserIssue(c::Client, page::Int64, state::String = "") = getParsedResponse(Vector{Issue},c,Requests.get,"/user/issues?page=$(page)")

listRepoIssues(c::Client,owner::String, repo::String, page::Int64, state::String = "") = getParsedResponse(Vector{Issue},c,Requests.get,"/repos/$(owner)/$(repo)/issues?page=$(page)")

getIssue(c::Client,owner::String, repo::String, index::Int64) = getParsedResponse(Issue,c,Requests.get,"/repos/$(owner)/$(repo)/issues/$(index)")

FieldTags.@tag immutable CreateIssueOption
	title::String
	body::Nullable{String}
	assignee::Nullable{String}
	milestone::Nullable{Int64}
	labels::Nullable{Vector{Int64}}
	closed::Nullable{Bool}
end

function CreateIssueOption(title; kwargs...)
	kwargsDict = Dict{Symbol,Any}(k=>v for (k,v) in kwargs)

	args = Any[title]

	for field in fieldnames(CreateIssueOption)[2:end]
		push!(args,get(kwargsDict,field, fieldtype(CreateIssueOption,field)()))
	end

	CreateIssueOption(args...)
end

createIssue(c::Client,owner::String,repo::String,opt::CreateIssueOption) = getParsedResponse(Issue,c,Requests.post,"/repos/$(owner)/$(repo)/issues"; json = marshalJSON(opt))

createIssue(c::Client,repo::Repository,opt::CreateIssueOption) = createIssue(c,repo.owner.userName,repo.name,opt)

FieldTags.@tag immutable EditIssueOption
	title::Nullable{String}
	body::Nullable{String}
	assignee::Nullable{String}
	milestone::Nullable{Int64}
	state::Nullable{String}
end

function EditIssueOption(;kwargs...)
	kwargsDict = Dict{Symbol,Any}(k=>v for (k,v) in kwargs)
	args = Any[]

	for field in fieldnames(EditIssueOption)
		push!(args,get(kwargsDict,field, fieldtype(EditIssueOption,field)()))
	end
	EditIssueOption(args...)
end

editIssue(c::Client,owner::String,repo::String, index::Int64, opt::EditIssueOption) = getParsedResponse(Issue,c,Requests.patch,"/repos/$(owner)/$(repo)/issues/$(index)"; json = marshalJSON(opt))

######################
### Issue Comments ###
######################

listIssueComments(c::Client,owner::String, repo::String, index::Int64) = getParsedResponse(Vector{Comment},c,Requests.get,"/repos/$(owner)/$(repo)/issues/$(index)/comments")

listRepoComments(c::Client,owner::String, repo::String) = getParsedResponse(Vector{Comment},c,Requests.get,"/repos/$(owner)/$(repo)/issues/comments")

createIssueComment(c::Client,owner::String, repo::String, index::Int64,body::String) = getParsedResponse(Comment,c,Requests.post,"/repos/$(owner)/$(repo)/issues/$(index)/comments"; json = Dict("body" => body))

createIssueComment(c::Client,repo::Repository,issue::Issue,body) = createIssueComment(c,repo.owner.userName,repo.name,issue.index,body)

editIssueComment(c::Client,owner::String, repo::String, index::Int64, commentID::Int64, body::String) = getParsedResponse(Comment,c,Requests.patch,"/repos/$(owner)/$(repo)/issues/$(index)/comments/$(commentID)"; json = Dict("body" => body))

deleteIssueComment(c::Client,owner::String, repo::String, index::Int64, commentID::Int64) = ( getResponse(c,Requests.delete,"/repos/$(owner)/$(repo)/issues/$(index)/comments/$(commentID)"); nothing)

####################
### Issue Labels ###
####################

listRepoLabels(c::Client, owner::String, repo::String) = getParsedResponse(Vector{Label},c,Requests.get,"/repos/$(owner)/$(repo)/labels")

listRepoLabels(c::Client, repo::Repository) = listRepoLabels(c,repo.owner.userName,repo.name)

getRepoLabel(c::Client, owner::String, repo::String, index::Int64) = getParsedResponse(Label,c,Requests.get,"/repos/$(owner)/$(repo)/labels/$(index)")

createLabel(c::Client, owner::String, repo::String, name::String, color::String) = getParsedResponse(Label,c,Requests.post,"/repos/$(owner)/$(repo)/labels"; json = Dict("name" => name, "color" => color))

createLabel(c::Client, repo::Repository, name::String, color::String) = createLabel(c,repo.owner.userName,repo.name, name, color)

FieldTags.@tag immutable EditLabelOption
	name::Nullable{String}
	color::Nullable{String}
end

function EditLabelOption(;kwargs...)
	kwargsDict = Dict{Symbol,Any}(k=>v for (k,v) in kwargs)
	args = Any[]

	for field in fieldnames(EditLabelOption)
		push!(args,get(kwargsDict,field, fieldtype(EditLabelOption,field)()))
	end
	EditLabelOption(args...)
end

editLabel(c::Client, owner::String, repo::String, index::Int64, opt::EditLabelOption) = getParsedResponse(Issue,c,Requests.patch,"/repos/$(owner)/$(repo)/labels/$(index)"; json = marshalJSON(opt))

deleteLabel(c::Client,owner::String, repo::String, index::Int64, commentID::Int64) = ( getResponse(c,Requests.delete,"/repos/$(owner)/$(repo)/labels/$(index)"); nothing)

getIssueLabels(c::Client, owner::String, repo::String, index::Int64) = getParsedResponse(Vector{Label},c,Requests.get,"/repos/$(owner)/$(repo)/issues/$(index)/labels")

addIssueLabels(c::Client, owner::String, repo::String, index::Int64, labels::Vector{Int64}) = getParsedResponse(Vector{Label},c,Requests.post,"/repos/$(owner)/$(repo)/issues/$(index)/labels"; json = Dict("labels" => labels))

replaceIssueLabels(c::Client, owner::String, repo::String, index::Int64, labels::Vector{Int64}) = getParsedResponse(Vector{Label},c,Requests.put,"/repos/$(owner)/$(repo)/issues/$(index)/labels"; json = Dict("labels" => labels))

deleteIssueLabel(c::Client,owner::String, repo::String, index::Int64, label::Int64) = ( getResponse(c,Requests.delete,"/repos/$(owner)/$(repo)/issues/$(index)/labels/$(issues)"); nothing)

deleteIssueLabels(c::Client,owner::String, repo::String, index::Int64) = ( getResponse(c,Requests.delete,"/repos/$(owner)/$(repo)/issues/$(index)/labels"); nothing)

########################
### Issue Milestones ###
########################

listRepoMilestones(c::Client, owner::String, repo::String) = getParsedResponse(Vector{Milestone},c,Requests.get,"/repos/$(owner)/$(repo)/milestones")

listRepoMilestones(c::Client, repo::Repository) = listRepoMilestones(c,repo.owner.userName,repo.name)

listRepoMilestone(c::Client, owner::String, repo::String, index::Int64) = getParsedResponse(Milestone,c,Requests.get,"/repos/$(owner)/$(repo)/milestones/$(index)")

FieldTags.@tag immutable CreateMilestoneOption
	title::String
	description::Nullable{String}
	deadline::Nullable{DateTime} => json:"due_on,format:y-m-dTH:M:SZ"
end

function CreateMilestoneOption(title;kwargs...)
	kwargsDict = Dict{Symbol,Any}(k=>v for (k,v) in kwargs)
	args = Any[title]

	for field in fieldnames(CreateMilestoneOption)[2:end]
		push!(args,get(kwargsDict,field, fieldtype(CreateMilestoneOption,field)()))
	end
	CreateMilestoneOption(args...)
end

createMilestone(c::Client, owner::String, repo::String, opt::CreateMilestoneOption) = getParsedResponse(Milestone,c,Requests.post,"/repos/$(owner)/$(repo)/milestones"; json = marshalJSON(opt))

createMilestone(c::Client, repo::Repository, opt) = createMilestone(c,repo.owner.userName,repo.name,opt)

FieldTags.@tag immutable EditMilestoneOption
	title::Nullable{String}
	description::Nullable{String}
	state::String
	deadline::Nullable{DateTime} => json:"due_on,format:y-m-dTH:M:SZ"
end

function EditMilestoneOption(;kwargs...)
	kwargsDict = Dict{Symbol,Any}(k=>v for (k,v) in kwargs)
	args = Any[]

	for field in fieldnames(EditMilestoneOption)
		push!(args,get(kwargsDict,field, fieldtype(EditMilestoneOption,field)()))
	end
	EditMilestoneOption(args...)
end

editMilestone(c::Client, owner::String, repo::String, index::Int64, opt::EditMilestoneOption) = getParsedResponse(Issue,c,Requests.patch,"/repos/$(owner)/$(repo)/milestones/$(index)"; json = marshalJSON(opt))

deleteMilestone(c::Client,owner::String, repo::String, index::Int64, commentID::Int64) = ( getResponse(c,Requests.delete,"/repos/$(owner)/$(repo)/milestones/$(index)"); nothing)


text:
#
# Support for automatic ordering of components
#

"""
    dependencies(md::ModelDef, comp_path::ComponentPath)

Return the set of component names that `comp_path` in `md` depends one, i.e.,
sources for which `comp_name` is the destination of an internal connection.
"""
function dependencies(md::ModelDef, comp_path::ComponentPath)
    conns = internal_param_conns(md)
    # For the purposes of the DAG, we don't treat dependencies on [t-1] as an ordering constraint
    deps = Set(c.src_comp_path for c in conns if (c.dst_comp_path == comp_path && c.offset == 0))
    return deps
end

"""
    comp_graph(md::ModelDef)

Return a MetaGraph containing a directed (LightGraph) graph of the components of
ModelDef `md`. Each vertex has a :name property with its component name.
"""
function comp_graph(md::ModelDef)
    comp_paths = [c.comp_path for c in compdefs(md)]
    graph = MetaDiGraph()

    for comp_path in comp_paths
        add_vertex!(graph, :path, comp_path)
    end

    set_indexing_prop!(graph, :path)

    for comp_path in comp_paths
        for dep_path in dependencies(md, comp_path)
            src = graph[dep_path,  :path]
            dst = graph[comp_path, :path]
            add_edge!(graph, src, dst)
        end
    end

    #TODO:  for now we can allow cycles since we aren't using the offset
    # if is_cyclic(graph)
    #     error("Component graph contains a cycle")
    # end

    return graph
end

"""
    _topological_sort(md::ModelDef)

Build a directed acyclic graph referencing the positions of the components in
the OrderedDict of model `md`, tracing dependencies to create the DAG.
Perform a topological sort on the graph for the given model and return a vector
of component paths in the order that will ensure dependencies are processed
prior to dependent components.
"""
function _topological_sort(md::ModelDef)
    graph = comp_graph(md)
    ordered = topological_sort_by_dfs(graph)
    paths = map(i -> graph[i, :path], ordered)
    return paths
end