text:
using OrderedCollections

"""
    last_point(abaco, ne::String, var::String)::Union{Nothing, Missing, Tuple{Int,Float64}}

Returns the most recent in time value for the `ne` node metric `var` as a tuple (time, value).

If the node `ne` is unknown returns `nothing` and if there are not values
for metric `var` returns `missing`.
"""
function last_point(abaco, ne::String, var::String)::Union{Nothing, Missing, Tuple{Int,Float64}}
    if haskey(abaco.node, ne)
        cursor = abaco.node[ne].currsnap
        ts = abaco.node[ne].snap[cursor].ts
        vals = abaco.node[ne].snap[cursor].vals
        if haskey(vals, var)
            return (ts, vals[var].value)
        end
    else
        # unknown ne node
        return nothing
    end
    # missing var for ne node
    missing
end

"""
    last_value(abaco, ne::String, var::String)::Union{Nothing, Missing, Float64}

Returns the most recent in time value for the `ne` node metric `var`

If the node `ne` is unknown returns `nothing` and if there are not values
for metric `var` returns `missing`.
"""
function last_value(abaco, ne::String, var::String)::Union{Nothing, Missing, Float64}
    if haskey(abaco.node, ne)
        cursor = abaco.node[ne].currsnap
        ts = abaco.node[ne].snap[cursor].ts
        vals = abaco.node[ne].snap[cursor].vals
        if haskey(vals, var)
            return vals[var].value
        end
    else
        # unknown ne node
        return nothing
    end
    # missing var for ne node
    missing
end

"""
    get_values(abaco, ne::String, var::String)::Dict{Int,Float64}

Returns the ordered by time sequence of `var` values for `ne` node.

The returned values dictionary is ordered by descending time, most
recent value first. The number of entries are at most equal to the
value of ages.
"""
function get_values(abaco, ne::String, var::String)::Dict{Int,Float64}
    result = OrderedDict{Int, Float64}()

    if haskey(abaco.node, ne)
        cursor = abaco.node[ne].currsnap
        for i = 1:abaco.ages
            vals = abaco.node[ne].snap[cursor].vals
            ts = abaco.node[ne].snap[cursor].ts
            if haskey(vals, var)
                result[ts] = vals[var].value
            end
            cursor = cursor == 1 ? abaco.ages : cursor - 1
        end
    end
    result
end

"""
    ingest!(abaco, payload)

Adds the input variables included in the `payload` dictionary.

The Dict `msg` must contains the keys `ts` and `ne` and a numbers
of other keys managed as input variables.

This function modifies the `payload` dictionary: `ne` and `ts` keys are popped out.

```julia
    payload = Dict(
        "ts" => nowts(),
        "ne" => "trento.castello",
        "x" => 23.2,
        "y" => 100
    )
    ingest!(abaco, ts, ne, payload)
```
"""
function ingest!(abaco, payload)
    ts = pop!(payload, "ts")
    ne = pop!(payload, "ne")
    ingest(abaco, ts, ne, payload)
end

"""
    ingest(abaco, ts, ne, values)

Adds the input variables include in the dictionary `values`.

```julia
    # now timestamp
    ts = nowts()

    # short name of network node
    ne = "trento.castello"

    values = Dict(
        "x" => 23.2,
        "y" => 100
    )
    ingest(abaco, ts, ne, values)
```
"""
function ingest(abaco, ts, ne, values)
    snap = getsnap(abaco, ts, ne)
    if snap !== nothing
        for (var, val) in values
            snap_add(snap, ne, var, val)
        end
        trigger_formulas(abaco, snap, ne, keys(values))
    end
end

function ingest(abaco, payload::Dict{String, Any})
    ts = payload["ts"]
    ne = payload["ne"]
    snap = getsnap(abaco, ts, ne)
    vars = [var for var in keys(payload) if !(var in ["ne", "ts"])]
    if snap !== nothing
        for var in vars
            snap_add(snap, ne, var, payload[var])
        end
        trigger_formulas(abaco, snap, ne, vars)
    end
end

"""
    ingest(abaco, ts::Int, ne::String, var::String, val::Real)

Adds the input variable `var` with value `val`.

* `ts`: timestamp with seconds resolution
* `ne`: scope name
* `var`: variable name
* `val`: variable value

"""
function ingest(abaco, ts::Int, ne::String, var::String, val::Real)
    snap = getsnap(abaco, ts, ne)
    snap_add(snap, ne, var, val)
    trigger_formulas(abaco, snap, ne, var)
end

function ingest(abaco, ts::Int, ne::String, var::String, val::Vector{<:Real})
    snap = getsnap(abaco, ts, ne)
    snap_add(snap, ne, var, val)
    trigger_formulas(abaco, snap, ne, var)
end

function sum_collected(abaco::Context, ne, variable, ts)
    result = get_collected(abaco, ne, variable, ts)
    if result === nothing
        return nothing
    else
        return PValue(length(result[1]), result[2], sum(result[1]), result[3])
    end
end

function get_collected(abaco::Context, ne, variable, ts)
    result = nothing

    (origin_type, var) = split(variable, ".")

    index = abaco.interval == -1 ? 1 : mvindex(ts, abaco.interval, abaco.ages)

    if haskey(abaco.origins, ne)
        result = LValue(length(abaco.origins[ne]), [])
        for origin_elem in abaco.origins[ne]
            snap = origin_elem.snap[index]
            if haskey(snap.vals, var) && origin_type == origin_elem.tag
                ropts = span(ts, abaco.interval)
                # @debug "[get_collected] interval: $(abaco.interval) - snap ts ts: $(snap.ts), ts: $ts, ropts: $ropts"
                if (ropts === snap.ts)
                    append!(result.value, snap.vals[var].value)
                    # updated is the last updated timestamp of origins
                    if result.updated < snap.vals[var].updated
                        result.updated = snap.vals[var].updated
                    end
                end
            end
        end
    end
    if result === nothing
        return nothing
    else
        return (result.value, result.contribs, result.updated)
    end
end

function sum_collected(abaco, ne, variable)
    ts = nowts()
    sum_collected(abaco, ne, variable, ts)
end

function get_collected(abaco, ne, variable)
    ts = nowts()
    get_collected(abaco, ne, variable, ts)
end

etype(abaco, ne) = haskey(abaco.node, ne) ? abaco.node[ne].tag : DEFAULT_TYPE

snapsetting(abaco, tag) = haskey(abaco.cfg, tag) ? abaco.cfg[tag] : abaco.cfg[DEFAULT_TYPE]

"""
    getsnap(abaco::Context, ts, ne)

Returns the `ne` node snapshot relative to timestamp `ts`.
"""
function getsnap(abaco::Context, ts, ne)
    tag = etype(abaco, ne)
    ## cfg = snapsetting(abaco, tag)

    if !haskey(abaco.node, ne)
        node(abaco, ne, tag)
    end

    if abaco.interval == -1
        snap =  abaco.node[ne].snap[1]
        snap.ts = ts
        return snap
    end

    # get the span
    ropts = span(ts, abaco.interval)
    index = mvindex(ts, abaco.interval, abaco.ages)

    elem = abaco.node[ne]
    snap = elem.snap[index]

    if snap.ts < ropts

        last_ts = elem.snap[elem.currsnap].ts
        if last_ts < ropts
            elem.currsnap = index
        end
        snap.ts  = ropts

        # reset the values of snap
        map!(v->begin
            isa(v, SValue) ? v.value=NaN : v.value = []
            v
        end, values(snap.vals))

        for fstate in values(snap.outputs)
            fstate.done = false
        end
    end
    if ropts < snap.ts
        return nothing
    end

    snap
end

"""
    snap_add(snap::Snap, ne, var::String, val::Real)

Adds the variable value of `ne` node to the `snap` snapshot.
"""
function snap_add(snap::Snap, ne, var::String, val::Real)
    if !haskey(snap.vals, var)
        # first arrival of variable var
        snap.vals[var] = SValue(val)
    else
        snap.vals[var].value = val
        snap.vals[var].updated = nowts()
    end
end

function snap_add(snap::Snap, ne, var::String, val::PValue)
    if val.contribs === val.expected
        # Add only completed qvalues to the snapshot
        if !haskey(snap.vals, var)
            # first arrival of variable var
            snap.vals[var] = SValue(val.value)
        else
            snap.vals[var].value = val.value
            snap.vals[var].updated = nowts()
        end
    end
end

function snap_add(snap::Snap, ne, var::String, val::Vector{<:Real})
    if !haskey(snap.vals, var) || isa(snap.vals[var], SValue)
        # first arrival of variable var
        snap.vals[var] = LValue(ne, val)
    else
        snap.vals[var].value = val
        snap.vals[var].updated = nowts()
    end
end

function snap_add(::Nothing, ne, var::String, values::Any)
    # the target snap is older of any current snap, do nothing
end

function snap_add(::Snap, ne, var::String, value::Any)
    # it is not a real value, do nothing
end

function setup_settings(abaco::Context,
                        tag;
                        handle=missing,
                        oncomplete=missing,
                        emitone=missing)
    if !haskey(abaco.cfg, tag)
        abaco.cfg[tag] = SnapsSetting(handle === missing ? nothing : handle,
                                       emitone === missing ? false : emitone,
                                       oncomplete === missing ? nothing : oncomplete)
    else
        if handle !== missing
            abaco.cfg[tag].handle = handle
        end
        if oncomplete !== missing
            abaco.cfg[tag].oncomplete = oncomplete
        end
        if emitone !== missing
            abaco.cfg[tag].emitone = emitone
        end
    end
end


text:
import Base: unique, unique!, ==, deleteat!

"""
$(TYPEDEF)

A basis over the states with parameters, independent variable, and possible exogenous controls.
It extends an `AbstractSystem` as defined in `ModelingToolkit.jl`. `f` can either be a Julia function which is able to use ModelingToolkit variables or
a vector of `eqs`.
It can be called with the typical SciML signature, meaning out of place with `f(u,p,t)`
or in place with `f(du, u, p, t)`. If control inputs are present, it is assumed that no control corresponds to
zero for all inputs. The corresponding function calls are `f(u,p,t,inputs)` and `f(du,u,p,t,inputs)` and need to
be specified fully.

The optional `implicits` declare implicit variables in the `Basis`, meaning variables representing the (measured) target of the system.
Right now only supported with the use of `ImplicitOptimizer`s.

If `linear_independent` is set to `true`, a linear independent basis is created from all atom functions in `f`.

If `simplify_eqs` is set to `true`, `simplify` is called on `f`.

Additional keyworded arguments include `name`, which can be used to name the basis, and
`observed` for defining observables.

# Fields
$(FIELDS)

# Example

```julia
using ModelingToolkit
using DataDrivenDiffEq

@parameters w[1:2] t
@variables u[1:2](t)

Ψ = Basis([u; sin.(w.*u)], u, parameters = p, iv = t)
```

## Note

The keyword argument `eval_expression` controls the function creation
behavior. `eval_expression=true` means that `eval` is used, so normal
world-age behavior applies (i.e. the functions cannot be called from
the function that generates them). If `eval_expression=false`,
then construction via GeneralizedGenerated.jl is utilized to allow for
same world-age evaluation. However, this can cause Julia to segfault
on sufficiently large basis functions. By default eval_expression=false.

"""
mutable struct Basis <: AbstractBasis
    """The equations of the basis"""
    eqs::Vector{Equation}
    """Dependent (state) variables"""
    states::Vector
    """Control variables"""
    ctrls::Vector
    """Parameters"""
    ps::Vector
    """Observed"""
    observed::Vector
    """Independent variable"""
    iv::Num
    """Implicit variables of the basis"""
    implicit::Vector
    """Internal function representation of the basis"""
    f::Function
    """Name of the basis"""
    name::Symbol
    """Internal systems"""
    systems::Vector{Basis}
end

## Constructors

function Basis(eqs::AbstractVector, states::AbstractVector;
    parameters::AbstractVector = [], iv = nothing,
    controls::AbstractVector = [], implicits = [],
    observed::AbstractVector = [],
    name = gensym(:Basis),
    simplify = false, linear_independent = false,
    eval_expression = false,
    kwargs...)

    iv === nothing && (iv = Symbolics.variable(:t))
    iv = value(iv)
    eqs = scalarize(eqs)

    states, controls, parameters, implicits, observed = value.(scalarize(states)), value.(scalarize(controls)), value.(scalarize(parameters)), value.(scalarize(implicits)), value.(scalarize(observed))

    eqs = [eq for eq in eqs if ~isequal(Num(eq),zero(Num))]

    if linear_independent
        eqs_ = create_linear_independent_eqs(eqs, simplify)
    else
        eqs_ = simplify ? ModelingToolkit.simplify.(eqs) : eqs
    end

    unique!(eqs, !simplify)

    f = _build_ddd_function(eqs, [states;implicits], parameters, iv, controls, eval_expression)

    eqs = [Symbolics.variable(:φ,i) ~ eq for (i,eq) ∈ enumerate(eqs_)]

    return Basis(collect(eqs), states, controls, parameters, observed, iv, implicits, f, name, Basis[])
end

function Basis(eqs::AbstractVector{Equation}, states::AbstractVector;
    parameters::AbstractVector = [], iv = nothing,
    controls::AbstractVector = [], implicits = [],
    observed::AbstractVector = [],
    name = gensym(:Basis),
    simplify = false, linear_independent = false,
    eval_expression = false,
    kwargs...)

    iv === nothing && (iv = Symbolics.variable(:t))
    iv = value(iv)
    eqs = scalarize(eqs)
    states, controls, parameters, implicits, observed = value.(scalarize(states)), value.(scalarize(controls)), value.(scalarize(parameters)), value.(scalarize(implicits)), value.(scalarize(observed))

    lhs = [x.lhs for x in eqs]
    # We filter out 0s
    eqs_ = [Num(x.rhs) for x in eqs if ~isequal(Num(x),zero(Num))]
    if linear_independent
        eqs_ = create_linear_independent_eqs(eqs_, simplify)
    else
        eqs_ = simplify ? ModelingToolkit.simplify.(eqs_) : eqs_
    end

    isnothing(iv) && (iv = Num(variable(:t)))
    unique!(eqs_, !simplify)

    f = _build_ddd_function(eqs_, [states;implicits], parameters, iv, controls, eval_expression)

    eqs = [lhs[i] ~ eq for (i,eq) ∈ enumerate(eqs_)]

    return Basis(collect(eqs), states, controls, parameters, observed, iv, implicits, f, name, Basis[])

end

function Basis(f::Function, states::AbstractVector; parameters::AbstractVector = [], controls::AbstractVector = [], implicits::AbstractVector = [],
     iv = nothing, kwargs...)

    isnothing(iv) && (iv = Num(Variable(:t)))

    try
        eqs = isempty(controls) ? f([states;implicits], parameters, iv) : f([states;implicits], parameters, iv, controls)
        return Basis(eqs, states, parameters = parameters, iv = iv, controls = controls, implicits = implicits; kwargs...)
    catch e
        rethrow(e)
    end
end

function Basis(n::Int)
    @parameters t
    @variables u[1:n](t)
    return Basis(u, u, iv = t)
end

## Printing

@inline function Base.print(io::IO, x::AbstractBasis)
    state = states(x)
    ps = parameters(x)
    Base.printstyled(io, "Model $(nameof(x)) with $(length(x)) equations\n"; bold=true)
    print(io, "States :")
    if length(state) < 5
        for xi in state
            print(io, " ", xi)
        end
    else
        print(io, " ", length(state))
    end

    if !isempty(ps)
         print(io, "\nParameters :")
         if length(ps) < 5
             for p_ in ps
                 print(io, " ", p_)
             end
         else
             print(io, " ", length(ps))
         end
     end

    println(io, "\nIndependent variable: $(get_iv(x))")
    println(io, "Equations")
    for (i,eq) ∈ enumerate(equations(x))
        if i < 5 || i == length(x)
        println(io, "$(eq.lhs) = $(eq.rhs)")
        elseif i == 5
            println(io, "...")
        else
            continue
        end
    end
end

@inline function Base.println(io::IO, x::AbstractBasis, fullview::DataType = Val{false})
    fullview == Val{false} && return print(io, x)

    state = states(x)
    ps = parameters(x)
    Base.printstyled(io, "Model $(nameof(x)) with $(length(x)) equations\n"; bold=true)
    print(io, "States :")
    if length(state) < 5
        for xi in state
            print(io, " ", xi)
        end
    else
        print(io, " ", length(state))
    end

    if !isempty(ps)
         print(io, "\nParameters :")
         if length(ps) < 5
             for p_ in ps
                 print(io, " ", p_)
             end
         else
             print(io, " ", length(ps))
         end
     end

    println(io, "\nIndependent variable: $(get_iv(x))")
    println(io, "Equations")
    for (i,eq) ∈ enumerate(equations(x))
        println(io, "$(eq.lhs) = $(eq.rhs)")
    end
end

## Getters

"""
    $(SIGNATURES)

    Returns the internal function representing the dynamics of the `Basis`. This can be called either inplace or out-of-place
    with the typical SciML signature `f(u,p,t)` or `f(du,u,p,t)`. If control variables are defined, the function can also be called
    by `f(u,p,t,control)` or `f(du,u,p,t,control)` and assumes `control .= 0` if no control is given.
"""
function dynamics(b::AbstractBasis)
    return get_f(b)
end

"""
$(SIGNATURES)

Return the implicit variables of the basis.
"""
function implicit_variables(b::AbstractBasis)
    return getfield(b, :implicit)
end

# For internal use
is_implicit(b::AbstractBasis) = !isempty(implicit_variables(b))

## Callable
get_f(b::AbstractBasis) = getfield(b, :f)

# Fallback
(b::AbstractBasis)(args...) = get_f(b)(args...)

# OOP
function (b::AbstractBasis)(x::AbstractVector{T} where T, p::AbstractVector{T} where T = parameters(b),
    t::T where T <: Number = get_iv(b))
    return get_f(b)(x,p,t)
end

function (b::AbstractBasis)(x::AbstractVector{T} where T, p::AbstractVector{T} where T,
    t::T where T <: Number , u::AbstractVector{T} where T)
    return get_f(b)(x,p,t, u)
end

function (b::AbstractBasis)(x::AbstractMatrix{T} where T)
    t = get_iv(b)
    return get_f(b)(x,parameters(b),[t for i in 1:size(x,2)])
end

function (b::AbstractBasis)(x::AbstractMatrix{T} where T, p::AbstractVector{T} where T,)
    t = get_iv(b)
    return get_f(b)(x,p,[t for i in 1:size(x,2)])
end

function (b::AbstractBasis)(x::AbstractMatrix{T} where T, p::AbstractVector{T} where T,
    t::AbstractVector{T} where T <: Number)
    return get_f(b)(x,p,t)
end

function (b::AbstractBasis)(x::AbstractMatrix{T} where T, p::AbstractVector{T} where T,
    t::AbstractVector{T} where T <: Number, u::AbstractMatrix{T} where T)
    return get_f(b)(x,p,t, u)
end

# IIP
function (b::AbstractBasis)(y::AbstractMatrix{T} where T, x::AbstractMatrix{T} where T)
    t = get_iv(b)
    return get_f(b)(y,x,parameters(b),[t for i in 1:size(x,2)])
end

function (b::AbstractBasis)(y::AbstractMatrix{T} where T, x::AbstractMatrix{T} where T, p::AbstractVector{T} where T,)
    t = get_iv(b)
    return get_f(b)(y,x,p,[t for i in 1:size(x,2)])
end

function (b::AbstractBasis)(y::AbstractMatrix{T} where T, x::AbstractMatrix{T} where T, p::AbstractVector{T} where T,
    t::AbstractVector{T} where T <: Number, )
    return get_f(b)(y,x,p,t)
end

function (b::AbstractBasis)(y::AbstractMatrix{T} where T, x::AbstractMatrix{T} where T, p::AbstractVector{T} where T,
    t::AbstractVector{T} where T <: Number, u::AbstractMatrix{T} where T)
    return get_f(b)(y,x,p,t,u)
end

## Information and Iteration

Base.length(x::AbstractBasis) = length(equations(x))
Base.size(x::AbstractBasis) = size(equations(x))

Base.getindex(x::AbstractBasis, idx) = getindex(equations(x), idx)
Base.firstindex(x::AbstractBasis) = firstindex(equations(x))
Base.lastindex(x::AbstractBasis) = lastindex(equations(x))
Base.iterate(x::AbstractBasis) = iterate(equations(x))
Base.iterate(x::AbstractBasis, id) = iterate(equations(x), id)

## Internal update

function update!(b::AbstractBasis, eval_expression = false)

    ff = _build_ddd_function([bi.rhs for bi in collect(equations(b))],
        states(b), parameters(b), [get_iv(b)],
        controls(b), eval_expression)
    Core.setfield!(b, :f, ff)

    return
end

function Base.setindex!(x::AbstractBasis, idx, val, eval_expression = false)
    setindex!(equations(x), idx, val)
    update!(x, eval_expression)
    return
end

## Derivatives

"""
    $(SIGNATURES)

    Returns a function representing the jacobian matrix / gradient of the `Basis` with respect to the
    states as a function with the common signature `f(u,p,t)` for out of place and `f(du, u, p, t)` for in place computation.
    If control variables are defined, the function can also be called by `f(u,p,t,control)` or `f(du,u,p,t,control)` and assumes `control .= 0` if no control is given.

    If the jacobian with respect to other variables is needed, it can be passed via a second argument.
"""
jacobian(x::Basis, eval_expression::Bool = false) = jacobian(x, states(x), eval_expression)

function jacobian(x::Basis, s, eval_expression::Bool = false)

    j = Symbolics.jacobian([xi.rhs for xi in equations(x)], s)

    return _build_ddd_function(j,
        states(x), parameters(x), get_iv(x),
        controls(x), eval_expression)
end

## Utilities
function Base.deleteat!(b::Symbolics.Arr{T,N}, idxs) where {T,N}
    deleteat!(Symbolics.unwrap(b), idxs)
end

function Base.unique(b::AbstractArray{SymbolicUtils.Sym}, simplify_eqs::Bool)
    b = simplify_eqs ? simplify.(b) : b
    returns = ones(Bool, size(b)...)
    N = maximum(eachindex(b))
    for i ∈ eachindex(b)
        returns[i] = !any([isequal(b[i], b[j]) for j in i+1:N])
    end
    return Num.(b[returns])
end

function Base.unique(b::AbstractArray{Num}, y::Bool)
    unique(collect(map(x->x.val, b)), y)
end

function Base.unique!(b::AbstractArray{T}, simplify_eqs = false) where {T <: Union{Num, SymbolicUtils.Sym}}
    bs = simplify_eqs ? simplify.(b) : b
    removes = zeros(Bool, size(bs)...)
    N = maximum(eachindex(bs))
    for i ∈ eachindex(bs)
        removes[i] = any([isequal(bs[i], bs[j]) for j in i+1:N])
    end
    deleteat!(b, removes)
end

function Base.unique(b::AbstractArray{Equation}, simplify_eqs::Bool)
    b = simplify_eqs ? simplify.(b) : b
    returns = ones(Bool, size(b)...)
    N = maximum(eachindex(b))
    for i ∈ eachindex(b)
        returns[i] = !any([isequal(b[i].rhs, b[j].rhs) for j in i+1:N])
    end
    return b[returns]
end

function Base.unique!(b::AbstractArray{Equation}, simplify_eqs::Bool)
    bs = [bi.rhs for bi in b]
    bs = simplify_eqs ? simplify.(bs) : bs
    removes = zeros(Bool, size(bs)...)
    N = maximum(eachindex(bs))
    for i ∈ eachindex(bs)
        removes[i] = any([isequal(bs[i], bs[j]) for j in i+1:N])
    end
    deleteat!(b, removes)
end

## Interfacing && merging

function Base.unique!(b::Basis, simplify_eqs = false; eval_expression = false)
    unique!(equations(b), simplify_eqs)
    update!(b, eval_expression)
end

function Base.unique(b::Basis; kwargs...)
    eqs = unique(equations(b))
    return Basis(eqs, states(b), parameters = parameters(b), iv = get_iv(b), kwargs...)
end

"""
    $(SIGNATURES)

    Delete the entries specified by `inds` and update the `Basis` accordingly.
"""
function Base.deleteat!(b::Basis, inds; eval_expression = false)
    deleteat!(equations(b), inds)
    update!(b, eval_expression)
    return
end

"""
    $(SIGNATURES)

    Push the equations(s) in `eq` into the basis and update all internal fields accordingly.
    `eq` can either be a single equation or an array. If `simplify_eq` is true, the equation will be simplified.
"""
function Base.push!(b::Basis, eqs::AbstractArray, simplify_eqs = true; eval_expression = false)
    @inbounds for eq ∈ eqs
        push!(b, eq, false)
    end
    unique!(b, simplify_eqs, eval_expression = eval_expression)
    return
end

function Base.push!(b::Basis, eq, simplify_eqs = true; eval_expression = false)
    push!(equations(b), variable(:φ, length(b)+1)~eq)
    unique!(b, simplify_eqs, eval_expression = eval_expression)
    return
end

function Base.push!(b::Basis, eq::Equation, simplify_eqs = true; eval_expression = false)
    push!(equations(b), eq)
    unique!(b, simplify_eqs, eval_expression = eval_expression)
    return
end

"""
    $(SIGNATURES)

    Return a new `Basis`, which is defined via the union of `x` and `y` .
"""
function Base.merge(x::Basis, y::Basis; eval_expression = false)
    b =  unique(vcat([xi.rhs  for xi ∈ equations(x)], [xi.rhs  for xi ∈ equations(y)]))
    vs = unique(vcat(states(x), states(y)))
    ps = unique(vcat(parameters(x), parameters(y)))

    c = unique(vcat(controls(x), controls(y)))
    observed = unique(vcat(get_observed(x), get_observed(y)))

    return Basis(Num.(b), vs, parameters = ps, controls = c, observed = observed, eval_expression = eval_expression)
end

"""
    $(SIGNATURES)

    Updates `x` to include the union of both `x` and `y`.
"""
function Base.merge!(x::Basis, y::Basis; eval_expression = false)
    push!(x, map(x->x.rhs, equations(y)))
    Core.setfield!(x, :states, unique(vcat(states(x),states(y))))
    Core.setfield!(x, :ps, unique(vcat(parameters(x), parameters(y))))
    Core.setfield!(x, :ctrls, unique(vcat(controls(x), controls(y))))
    Core.setfield!(x, :observed, unique(vcat(get_observed(x), get_observed(y))))
    update!(x, eval_expression)
    return
end

## Additional functionalities

function (==)(x::Basis, y::Basis)
    length(x) == length(y) || return false
    n = zeros(Bool, length(x))
    yrhs = [yi.rhs for yi in equations(y)]
    xrhs = [xi.rhs for xi in equations(x)]
    @inbounds for (i, xi) in enumerate(xrhs)
        n[i] = any(isequal.([xi], yrhs))
        !n[i] && break
    end
    return all(n)
end

free_parameters(b::AbstractBasis; operations = [+]) = count_operation([xi.rhs for xi in equations(b)], operations) + length(b)

ext: jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "GSW"
version = v"0.0.1"

# Collection of sources required to build GSW
sources = [
    "http://www.teos-10.org/software/GSW-C-3.05-4.zip" =>
    "ea72c3311f35de42c5e838d05709c3dfd36570c7780629586353ee670ad5f145",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd GSW-C-3.05.0-4/
sed -i 's#/usr#../../destdir/#' GNUmakefile
make
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libgswteos", :libgswteos10),
    ExecutableProduct(prefix, "gsw_check", :libgswcheck)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
module Aquarius

using VegaLite

export histogram

include("histogram.jl")
end # module


text:
using Test

import PlutoTool: plutotool, commands, lookup_command, get_notebook, test_context, interactive_context,
                  CommandException, CommandNotFound, BadOption, BadOption, CellNotEmpty,
                  new_notebook, new_cell, find_empty, find, delete, set_contents, set_workspace

HELP_DOC_FILE = "../PlutoTool_Commands.md"

@testset "Verify help documentation" begin
  ctx = test_context()
  plutotool(ctx, "help")
  out = String(take!(ctx.stdout))
  err = String(take!(ctx.stderr))
  @test length(err) == 0
  # At least two lines of output per command
  @test length(split(out, "\n")) >= 2 * length(commands)
  # Every command is documented:
  for cmd in commands
    @test occursin(string(cmd), out)
  end
  open(f -> write(f, out),
       HELP_DOC_FILE; truncate=true)
end

@testset "Build a sample notebook from scratch and test with it" begin
  ctx = test_context()
  # Missing command detection:
  @test_throws CommandNotFound lookup_command("no-such-command")
  let notebook_file = "SampleNotebook"
    rm(notebook_file, force=true)
    # Look for notebook that doesn't exist:
    @test_throws SystemError find_empty(ctx, notebook_file)
    # Create notebook:
    new_notebook(ctx, notebook_file)
    @test isfile(notebook_file)
    @test length(get_notebook(notebook_file).cells) == 1
    # Find empty cell
    empty = find_empty(ctx, notebook_file)
    @test length(empty) == 1
    # Set its content:
    original_cell = empty[1]
    set_contents(ctx, notebook_file, original_cell, "# This is a cell.")
    found = find(ctx, notebook_file, "This", "cell")
    @test length(found) == 1
    @test found[1] == original_cell
    # Add a cell before it:
    before = new_cell(ctx, "before", notebook_file, original_cell)
    # Add a cell after it:
    after = new_cell(ctx, "after", notebook_file, original_cell)
    @test length(get_notebook(notebook_file).cells) == 3
    # Did they get added?
    set_contents(ctx, notebook_file, before, "# Before.")
    set_contents(ctx, notebook_file, after, "# After.")
    # negative find cell test
    found = find(ctx, notebook_file, "notfound")
    @test length(found) == 0
    # negative enpty test
    @test length(find_empty(ctx, notebook_file)) == 0
    # Deletion:
    @test length(find_empty(ctx, notebook_file)) == 0
    delete_me = new_cell(ctx, "after", notebook_file, original_cell)
    @test length(find_empty(ctx, notebook_file)) == 1
    delete(ctx, notebook_file, delete_me)
    @test length(find_empty(ctx, notebook_file)) == 0
  end
end

@testset "test set_workspace" begin
    ctx = test_context()
    let notebook_file = "NotebookWithWorkspace"
        rm(notebook_file, force=true)
        # Create notebook:
        new_notebook(ctx, notebook_file)
        set_workspace(ctx, notebook_file, "..")
    end
end


text:
# Description:
# Dimensions: d

# The Trid function has no local minimum except the global one. It is shown here in its two-dimensional form.

# Input Domain:
# The function is usually evaluated on the hypercube xi ∈ [-d2, d2], for all i = 1, …, d.

# Global Minimum:
# https://www.sfu.ca/~ssurjano/trid3.png

function trid(xx)
    d = length(xx);
    sum1 = (xx[1] - 1)^2;
    sum2 = 0;

    for ii = 2:d
        xi = xx[ii];
        xold = xx[ii - 1];
        sum1 = sum1 + (xi - 1)^2;
        sum2 = sum2 + xi * xold;
    end

    y = sum1 - sum2;

    return [y]
end


text:
module Twiddle

include("nibbles.jl")
include("bitpairs.jl")

end


text:
using ThreeBodyProblem
using LinearAlgebra
using Test

sys = sun_earth()

Az = 0.001
Lpt = 2 # Which libration point will I center about
NS = 1
npts = 100
t, rvs, T, Ax = rich3(sys, Az, Lpt, NS, npts)

rv₀, T = differential_corrector(sys, rvs[1,:], myconst=3, tf=T)

Wsp, Wsn, Wup, Wun, Φₜ, ww = invariant_manifolds(sys, rv₀, T, tf=5*T, nPts=20)


text:
"""
Runs the semi-distributed bucket model for a given area with one precipitation and temperature input.

$(SIGNATURES)

The function returns the modeled discharge for each timestep, the modelled snow extend for each elevation and timestep and the Waterbalance.
"""
function run_model(Area, Evaporation_Mean::Array{Float64,1}, Precipitation::Array{Float64,2}, Temp::Array{Float64,2},
                bare_input::HRU_Input, forest_input::HRU_Input, grass_input::HRU_Input, rip_input::HRU_Input,
                bare_storage::Storages, forest_storage::Storages, grass_storage::Storages, rip_storage::Storages, Slowstorage::Float64,
                bare_parameters::Parameters, forest_parameters::Parameters, grass_parameters::Parameters, rip_parameters::Parameters, slow_parameters::Slow_Paramters, Total_Elevationbands, Elevation_Percentage)
    # the function takes as input the parameters of each HRU, the inital storage values of each HRU, the inital value of the slow storage
    # KS, ratio riparian, all inputs

    # define the maximum time
    tmax::Int128 = length(Precipitation[:,1])

    Discharge::Array{Float64,1} = zeros(tmax)
    WBtotal::Array{Float64,1} = zeros(tmax)
    Snow_Extend::Array{Float64,2} = zeros(tmax, Total_Elevationbands)

    # store the initial storage values, Does not take into account GW storage!!
    # Assumption_ GW storage is 0 at start
    #TO DO: average values over area
    # as form now only works if storage input is 0 at start
    Initial_Storage_bare::Float64 = bare_storage.Fast + sum(bare_storage.Interception) + sum(bare_storage.Snow) + bare_storage.Soil
    Initial_Storage_forest::Float64 = forest_storage.Fast + sum(forest_storage.Interception) + sum(forest_storage.Snow) + forest_storage.Soil
    Initial_Storage_grass::Float64 = grass_storage.Fast + sum(grass_storage.Interception) + sum(grass_storage.Snow) + grass_storage.Soil
    Initial_Storage_rip::Float64 = rip_storage.Fast + sum(rip_storage.Interception) + sum(rip_storage.Snow) + rip_storage.Soil
    Initial_Storage::Float64 = Initial_Storage_bare + Initial_Storage_forest + Initial_Storage_grass + Initial_Storage_rip + Slowstorage

    for t in 1:tmax
        # at each timestep new temp, precipitation and Epot values have to be delivered
        # areas don't change
        # riparian discharge from former timestep has to be used
        # gives the current precipitation, evaporation and temperature
        Evaporation_Mean_Current = Evaporation_Mean[t]
        Precipitation_Current = Precipitation[t, :]
        Temperature_Current = Temp[t, :]

        bare_input::HRU_Input = input_timestep(bare_input, Evaporation_Mean_Current, Precipitation_Current, Temperature_Current)
        forest_input::HRU_Input = input_timestep(forest_input, Evaporation_Mean_Current, Precipitation_Current, Temperature_Current)
        grass_input::HRU_Input = input_timestep(grass_input, Evaporation_Mean_Current, Precipitation_Current, Temperature_Current)
        rip_input::HRU_Input = input_timestep(rip_input, Evaporation_Mean_Current, Precipitation_Current, Temperature_Current)

        Riparian_Discharge::Float64, Total_Discharge::Float64, Total_Interception_Evaporation::Float64, Total_Soil_Evaporation::Float64, bare_storage::Storages, forest_storage::Storages, grass_storage::Storages, rip_storage::Storages, Slowstorage::Float64, WB::Float64, Total_Prec::Float64 = allHRU(bare_input, forest_input, grass_input, rip_input,
                                                                                bare_storage, forest_storage, grass_storage, rip_storage,
                                                                                bare_parameters, forest_parameters, grass_parameters, rip_parameters,
                                                                                Slowstorage, slow_parameters)
        # give new riparian discharge as input for next timestep
        rip_input.Riparian_Discharge = Riparian_Discharge
        # save the fluxes of the current timestep
        Discharge[t]::Float64 = Total_Discharge/1000 * Area / (3600 * 24) # mm convert it to meter and than * area / seconds in one day
        # OPTIONAL: store the total amount of snow in each elevation (considering whole catchment)
        Snow_Elevations, Snow_Extend[t,:]::Array{Float64,1} = snowperelevation(bare_input, forest_input, grass_input, rip_input, bare_storage, forest_storage, grass_storage, rip_storage, Total_Elevationbands)
        #Bare_Snow[t,:]::Array{Float64,1} = bare_storage.Snow
        Snow_Extend[t,:] = Snow_Extend[t,:] ./ Elevation_Percentage

        WBtotal[t]::Float64 = WB
    end

    # Check Water Balance
    Waterbalance = sum(WBtotal)::Float64
    @assert Waterbalance <= 10^(-10)
    return Discharge::Array{Float64,1}, Snow_Extend::Array{Float64,2}, Waterbalance::Float64
end

"""
Runs the semi-distributed bucket model for a given area with one precipitation and temperature input.

$(SIGNATURES)

The function returns the modeled discharge for each timestep, the modelled snow extend for each elevation and timestep and the Waterbalance.
"""
function run_model_glacier(Area, Evaporation_Mean::Array{Float64,1},Glacier::Array{Float64,2}, Precipitation::Array{Float64,2}, Temp::Array{Float64,2},
                bare_input::HRU_Input, forest_input::HRU_Input, grass_input::HRU_Input, rip_input::HRU_Input,
                bare_storage::Storages, forest_storage::Storages, grass_storage::Storages, rip_storage::Storages, Slowstorage::Float64,
                bare_parameters::Parameters, forest_parameters::Parameters, grass_parameters::Parameters, rip_parameters::Parameters, slow_parameters::Slow_Paramters, Total_Elevationbands, Elevation_Percentage)
    # the function takes as input the parameters of each HRU, the inital storage values of each HRU, the inital value of the slow storage
    # KS, ratio riparian, all inputs

    # define the maximum time
    tmax::Int128 = length(Precipitation[:,1])

    Discharge::Array{Float64,1} = zeros(tmax)
    WBtotal::Array{Float64,1} = zeros(tmax)
    Snow_Extend::Array{Float64,2} = zeros(tmax, Total_Elevationbands)

    # store the initial storage values, Does not take into account GW storage!!
    # Assumption_ GW storage is 0 at start
    #TO DO: average values over area
    # as form now only works if storage input is 0 at start
    Initial_Storage_bare::Float64 = bare_storage.Fast + sum(bare_storage.Interception) + sum(bare_storage.Snow) + bare_storage.Soil
    Initial_Storage_forest::Float64 = forest_storage.Fast + sum(forest_storage.Interception) + sum(forest_storage.Snow) + forest_storage.Soil
    Initial_Storage_grass::Float64 = grass_storage.Fast + sum(grass_storage.Interception) + sum(grass_storage.Snow) + grass_storage.Soil
    Initial_Storage_rip::Float64 = rip_storage.Fast + sum(rip_storage.Interception) + sum(rip_storage.Snow) + rip_storage.Soil
    Initial_Storage::Float64 = Initial_Storage_bare + Initial_Storage_forest + Initial_Storage_grass + Initial_Storage_rip + Slowstorage

    for t in 1:tmax
        # at each timestep new temp, precipitation and Epot values have to be delivered
        # areas don't change
        # riparian discharge from former timestep has to be used
        # gives the current precipitation, evaporation and temperature
        Evaporation_Mean_Current = Evaporation_Mean[t]
        Precipitation_Current = Precipitation[t, :]
        Temperature_Current = Temp[t, :]
        Glacier_Current = Glacier[:,t]

        bare_input::HRU_Input = input_timestep_glacier(bare_input, Evaporation_Mean_Current, Glacier_Current, Precipitation_Current, Temperature_Current)
        forest_input::HRU_Input = input_timestep(forest_input, Evaporation_Mean_Current, Precipitation_Current, Temperature_Current)
        grass_input::HRU_Input = input_timestep(grass_input, Evaporation_Mean_Current, Precipitation_Current, Temperature_Current)
        rip_input::HRU_Input = input_timestep(rip_input, Evaporation_Mean_Current, Precipitation_Current, Temperature_Current)

        Riparian_Discharge::Float64, Total_Discharge::Float64, Total_Interception_Evaporation::Float64, Total_Soil_Evaporation::Float64, bare_storage::Storages, forest_storage::Storages, grass_storage::Storages, rip_storage::Storages, Slowstorage::Float64, WB::Float64, Total_Prec::Float64 = allHRU(bare_input, forest_input, grass_input, rip_input,
                                                                                bare_storage, forest_storage, grass_storage, rip_storage,
                                                                                bare_parameters, forest_parameters, grass_parameters, rip_parameters,
                                                                                Slowstorage, slow_parameters)
        # give new riparian discharge as input for next timestep
        rip_input.Riparian_Discharge = Riparian_Discharge
        # save the fluxes of the current timestep
        Discharge[t]::Float64 = Total_Discharge/1000 * Area / (3600 * 24) # mm convert it to meter and than * area / seconds in one day
        # OPTIONAL: store the total amount of snow in each elevation (considering whole catchment)
        Snow_Elevations, Snow_Extend[t,:]::Array{Float64,1} = snowperelevation(bare_input, forest_input, grass_input, rip_input, bare_storage, forest_storage, grass_storage, rip_storage, Total_Elevationbands)
        #Bare_Snow[t,:]::Array{Float64,1} = bare_storage.Snow
        Snow_Extend[t,:] = Snow_Extend[t,:] ./ Elevation_Percentage

        WBtotal[t]::Float64 = WB
    end

    # Check Water Balance
    Waterbalance = sum(WBtotal)::Float64
    @assert Waterbalance <= 10^(-10)
    return Discharge::Array{Float64,1}, Snow_Extend::Array{Float64,2}, Waterbalance::Float64
end

"""
Runs the semi-distributed bucket model for a given area with one precipitation and temperature input.

$(SIGNATURES)

The function returns the modeled discharge for each timestep, the modelled snow extend for each elevation and timestep and the Waterbalance and the snowstorage.
"""
function run_model_future(Area, Evaporation_Mean::Array{Float64,1}, Precipitation::Array{Float64,2}, Temp::Array{Float64,2},
                bare_input::HRU_Input, forest_input::HRU_Input, grass_input::HRU_Input, rip_input::HRU_Input,
                bare_storage::Storages, forest_storage::Storages, grass_storage::Storages, rip_storage::Storages, Slowstorage::Float64,
                bare_parameters::Parameters, forest_parameters::Parameters, grass_parameters::Parameters, rip_parameters::Parameters, slow_parameters::Slow_Paramters, Total_Elevationbands, Elevation_Percentage)
    # the function takes as input the parameters of each HRU, the inital storage values of each HRU, the inital value of the slow storage
    # KS, ratio riparian, all inputs

    # define the maximum time
    tmax::Int128 = length(Precipitation[:,1])

    Discharge::Array{Float64,1} = zeros(tmax)
    Total_Melt::Array{Float64,1} = zeros(tmax)
    WBtotal::Array{Float64,1} = zeros(tmax)
    Snow_Extend::Array{Float64,2} = zeros(tmax, Total_Elevationbands)
    Snowstorage::Array{Float64,1} = zeros(tmax)

    # store the initial storage values, Does not take into account GW storage!!
    # Assumption_ GW storage is 0 at start
    #TO DO: average values over area
    # as form now only works if storage input is 0 at start
    Initial_Storage_bare::Float64 = bare_storage.Fast + sum(bare_storage.Interception) + sum(bare_storage.Snow) + bare_storage.Soil
    Initial_Storage_forest::Float64 = forest_storage.Fast + sum(forest_storage.Interception) + sum(forest_storage.Snow) + forest_storage.Soil
    Initial_Storage_grass::Float64 = grass_storage.Fast + sum(grass_storage.Interception) + sum(grass_storage.Snow) + grass_storage.Soil
    Initial_Storage_rip::Float64 = rip_storage.Fast + sum(rip_storage.Interception) + sum(rip_storage.Snow) + rip_storage.Soil
    Initial_Storage::Float64 = Initial_Storage_bare + Initial_Storage_forest + Initial_Storage_grass + Initial_Storage_rip + Slowstorage

    for t in 1:tmax
        # at each timestep new temp, precipitation and Epot values have to be delivered
        # areas don't change
        # riparian discharge from former timestep has to be used
        # gives the current precipitation, evaporation and temperature
        Evaporation_Mean_Current = Evaporation_Mean[t]
        Precipitation_Current = Precipitation[t, :]
        Temperature_Current = Temp[t, :]

        bare_input::HRU_Input = input_timestep(bare_input, Evaporation_Mean_Current, Precipitation_Current, Temperature_Current)
        forest_input::HRU_Input = input_timestep(forest_input, Evaporation_Mean_Current, Precipitation_Current, Temperature_Current)
        grass_input::HRU_Input = input_timestep(grass_input, Evaporation_Mean_Current, Precipitation_Current, Temperature_Current)
        rip_input::HRU_Input = input_timestep(rip_input, Evaporation_Mean_Current, Precipitation_Current, Temperature_Current)

        Riparian_Discharge::Float64, Total_Discharge::Float64, Total_Interception_Evaporation::Float64, Total_Soil_Evaporation::Float64, bare_storage::Storages, forest_storage::Storages, grass_storage::Storages, rip_storage::Storages, Slowstorage::Float64, WB::Float64, Total_Prec::Float64, Melt::Float64 = allHRU_future(bare_input, forest_input, grass_input, rip_input,
                                                                                bare_storage, forest_storage, grass_storage, rip_storage,
                                                                                bare_parameters, forest_parameters, grass_parameters, rip_parameters,
                                                                                Slowstorage, slow_parameters)
        # give new riparian discharge as input for next timestep
        rip_input.Riparian_Discharge = Riparian_Discharge
        # save the fluxes of the current timestep
        Discharge[t]::Float64 = Total_Discharge/1000 * Area / (3600 * 24) # mm convert it to meter and than * area / seconds in one day
        # OPTIONAL: store the total amount of snow in each elevation (considering whole catchment)
        Snow_Elevations, Snow_Extend[t,:]::Array{Float64,1} = snowperelevation(bare_input, forest_input, grass_input, rip_input, bare_storage, forest_storage, grass_storage, rip_storage, Total_Elevationbands)
        #Bare_Snow[t,:]::Array{Float64,1} = bare_storage.Snow
        Snow_Extend[t,:] = Snow_Extend[t,:] ./ Elevation_Percentage
        Total_Melt[t] = Melt

        Bare_Interceptionstorage::Float64, Bare_Snowstorage::Float64 = Storage_Total(bare_storage, bare_input)
        Forest_Interceptionstorage::Float64, Forest_Snowstorage::Float64 = Storage_Total(forest_storage, forest_input)
        Grass_Interceptionstorage::Float64, Grass_Snowstorage::Float64 = Storage_Total(grass_storage, grass_input)
        Rip_Interceptionstorage::Float64, Rip_Snowstorage::Float64 = Storage_Total(rip_storage, rip_input)

        Snowstorage[t] = Bare_Snowstorage * bare_input.Area_HRU + Forest_Snowstorage * forest_input.Area_HRU + Grass_Snowstorage * grass_input.Area_HRU + Rip_Snowstorage * rip_input.Area_HRU

        WBtotal[t]::Float64 = WB
    end

    # Check Water Balance
    Waterbalance = sum(WBtotal)::Float64
    @assert Waterbalance <= 10^(-10)
    return Discharge::Array{Float64,1}, Snow_Extend::Array{Float64,2}, Snowstorage::Array{Float64,1}, Waterbalance::Float64, Total_Melt::Array{Float64,1}
end

"""
Runs the semi-distributed bucket model for a given area with one precipitation and temperature input.

$(SIGNATURES)

The function returns the modeled discharge for each timestep, the modelled snow extend for each elevation and timestep and the Waterbalance and the snowstorage.
"""
function run_model_glacier_future(Area, Evaporation_Mean::Array{Float64,1}, Glacier::Array{Float64,2}, Precipitation::Array{Float64,2}, Temp::Array{Float64,2},
                bare_input::HRU_Input, forest_input::HRU_Input, grass_input::HRU_Input, rip_input::HRU_Input,
                bare_storage::Storages, forest_storage::Storages, grass_storage::Storages, rip_storage::Storages, Slowstorage::Float64,
                bare_parameters::Parameters, forest_parameters::Parameters, grass_parameters::Parameters, rip_parameters::Parameters, slow_parameters::Slow_Paramters, Total_Elevationbands, Elevation_Percentage)
    # the function takes as input the parameters of each HRU, the inital storage values of each HRU, the inital value of the slow storage
    # KS, ratio riparian, all inputs

    # define the maximum time
    tmax::Int128 = length(Precipitation[:,1])

    Discharge::Array{Float64,1} = zeros(tmax)
    Total_Melt::Array{Float64,1} = zeros(tmax)
    WBtotal::Array{Float64,1} = zeros(tmax)
    Snow_Extend::Array{Float64,2} = zeros(tmax, Total_Elevationbands)
    Snowstorage::Array{Float64,1} = zeros(tmax)

    # store the initial storage values, Does not take into account GW storage!!
    # Assumption_ GW storage is 0 at start
    #TO DO: average values over area
    # as form now only works if storage input is 0 at start
    Initial_Storage_bare::Float64 = bare_storage.Fast + sum(bare_storage.Interception) + sum(bare_storage.Snow) + bare_storage.Soil
    Initial_Storage_forest::Float64 = forest_storage.Fast + sum(forest_storage.Interception) + sum(forest_storage.Snow) + forest_storage.Soil
    Initial_Storage_grass::Float64 = grass_storage.Fast + sum(grass_storage.Interception) + sum(grass_storage.Snow) + grass_storage.Soil
    Initial_Storage_rip::Float64 = rip_storage.Fast + sum(rip_storage.Interception) + sum(rip_storage.Snow) + rip_storage.Soil
    Initial_Storage::Float64 = Initial_Storage_bare + Initial_Storage_forest + Initial_Storage_grass + Initial_Storage_rip + Slowstorage

    for t in 1:tmax
        # at each timestep new temp, precipitation and Epot values have to be delivered
        # areas don't change
        # riparian discharge from former timestep has to be used
        # gives the current precipitation, evaporation and temperature
        Evaporation_Mean_Current = Evaporation_Mean[t]
        Precipitation_Current = Precipitation[t, :]
        Temperature_Current = Temp[t, :]
        Glacier_Current = Glacier[:,t]

        bare_input::HRU_Input = input_timestep_glacier(bare_input, Evaporation_Mean_Current, Glacier_Current, Precipitation_Current, Temperature_Current)
        forest_input::HRU_Input = input_timestep(forest_input, Evaporation_Mean_Current, Precipitation_Current, Temperature_Current)
        grass_input::HRU_Input = input_timestep(grass_input, Evaporation_Mean_Current, Precipitation_Current, Temperature_Current)
        rip_input::HRU_Input = input_timestep(rip_input, Evaporation_Mean_Current, Precipitation_Current, Temperature_Current)
        Riparian_Discharge::Float64, Total_Discharge::Float64, Total_Interception_Evaporation::Float64, Total_Soil_Evaporation::Float64, bare_storage::Storages, forest_storage::Storages, grass_storage::Storages, rip_storage::Storages, Slowstorage::Float64, WB::Float64, Total_Prec::Float64, Melt::Float64 = allHRU_future(bare_input, forest_input, grass_input, rip_input,
                                                                                bare_storage, forest_storage, grass_storage, rip_storage,
                                                                                bare_parameters, forest_parameters, grass_parameters, rip_parameters,
                                                                                Slowstorage, slow_parameters)
        # give new riparian discharge as input for next timestep
        rip_input.Riparian_Discharge = Riparian_Discharge
        # save the fluxes of the current timestep
        Discharge[t]::Float64 = Total_Discharge/1000 * Area / (3600 * 24) # mm convert it to meter and than * area / seconds in one day
        # OPTIONAL: store the total amount of snow in each elevation (considering whole catchment)
        Snow_Elevations, Snow_Extend[t,:]::Array{Float64,1} = snowperelevation(bare_input, forest_input, grass_input, rip_input, bare_storage, forest_storage, grass_storage, rip_storage, Total_Elevationbands)
        #Bare_Snow[t,:]::Array{Float64,1} = bare_storage.Snow
        Snow_Extend[t,:] = Snow_Extend[t,:] ./ Elevation_Percentage
        Total_Melt[t] = Melt

        Bare_Interceptionstorage::Float64, Bare_Snowstorage::Float64 = Storage_Total(bare_storage, bare_input)
        Forest_Interceptionstorage::Float64, Forest_Snowstorage::Float64 = Storage_Total(forest_storage, forest_input)
        Grass_Interceptionstorage::Float64, Grass_Snowstorage::Float64 = Storage_Total(grass_storage, grass_input)
        Rip_Interceptionstorage::Float64, Rip_Snowstorage::Float64 = Storage_Total(rip_storage, rip_input)

        Snowstorage[t] = Bare_Snowstorage * bare_input.Area_HRU + Forest_Snowstorage * forest_input.Area_HRU + Grass_Snowstorage * grass_input.Area_HRU + Rip_Snowstorage * rip_input.Area_HRU

        WBtotal[t]::Float64 = WB
    end

    # Check Water Balance
    Waterbalance = sum(WBtotal)::Float64
    #@assert Waterbalance <= 10^(-10)
    return Discharge::Array{Float64,1}, Snow_Extend::Array{Float64,2}, Snowstorage::Array{Float64,1}, Waterbalance::Float64, Total_Melt::Array{Float64,1}
end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
using IncrementalInference

@testset "test fourdoor early example" begin

global N=100
global fg = initfg()
global doors = reshape(Float64[-100.0;0.0;100.0;300.0],1,4)
global pd = kde!(doors,[3.0])
global pd = resample(pd,N);
global bws = getBW(pd)[:,1]
global doors2 = getPoints(pd);

addVariable!(fg,:x0,ContinuousScalar,N=N)
addFactor!(fg,[:x0], Prior( pd ) )

# tem = 2.0*randn(1,N)+getVal(v1)+50.0
addVariable!(fg,:x2, ContinuousScalar, N=N)
addFactor!(fg, [:x0; :x2], LinearConditional(Normal(50.0,2.0)) )
# addFactor!(fg, [v1;v2], Odo(50.0*ones(1,1),[2.0]',[1.0]))

# monocular sighting would look something like
#addFactor!(fg, Mono, [:x3,:l1], [14.0], [1.0], [1.0])
#addFactor!(fg, Mono, [:x4,:l1], [11.0], [1.0], [1.0])

addVariable!(fg,:x3,ContinuousScalar, N=N)
addFactor!(fg,[:x2;:x3], LinearConditional( Normal(50.0,4.0)) )
addFactor!(fg,[:x3], Prior( pd ))

addVariable!(fg,:x4,ContinuousScalar, N=N)
addFactor!(fg,[:x3;:x4], LinearConditional( Normal(50.0,2.0)) )

addVariable!(fg, :l1, ContinuousScalar, N=N)
addFactor!(fg, [:x3,:l1], Ranged([64.0],[0.5],[1.0]))
addFactor!(fg, [:x4,:l1], Ranged([16.0],[0.5],[1.0]))

addVariable!(fg,:x5,ContinuousScalar, N=N)
addFactor!(fg,[:x4;:x5], LinearConditional( Normal(50.0,2.0)) )

addVariable!(fg,:x6,ContinuousScalar, N=N)
addFactor!(fg,[:x5;:x6], LinearConditional( Normal(40.0,1.20)) )

addVariable!(fg,:x7,ContinuousScalar, N=N)
addFactor!(fg,[:x6;:x7], LinearConditional( Normal(60.0,2.0)) )

# ensureAllInitialized!(fg)

global mlc = MixturePrior(Normal.(doors[1,:], bws[1]), 0.25*ones(4))

# getSample(mlc)

addFactor!(fg,[:x7], mlc )

# HMM computed ground truth, extended for 7 poses with landmark
global gt = Dict{Symbol, Array{Float64,2}}()
gt[:x0]=reshape(Float64[0.0;1.97304 ],2,1) # -0.0342366
gt[:x2]=reshape(Float64[50.0; 2.83153 ],2,1) # 49.8797
gt[:x3]=reshape(Float64[100.0; 1.65557 ],2,1) # 99.8351
gt[:x4]=reshape(Float64[150.0; 1.64945 ],2,1) # 148.637
gt[:x5]=reshape(Float64[200.0; 1.77992 ],2,1) # 198.62
gt[:x6]=reshape(Float64[240.0; 2.20466 ],2,1) # 238.492
gt[:x7]=reshape(Float64[300.0; 2.14353 ],2,1) # 298.467
gt[:l1]=reshape(Float64[165.0; 1.17284 ],2,1) # 164.102

tree, smt, hist = solveTree!(fg)

# global tree = prepBatchTree!(fg, drawpdf=false);
#
# # list vertices in fg
# @show xx,ll = IncrementalInference.ls(fg)
#
# # do belief propagation inference over tree once
# # using recursive single core approach (better stack trace for development)
# # inferOverTreeR!(fg, tree)
# inferOverTree!(fg, tree, N=N, dbg=true, treeinit=true)
#
#  #
# # test multi-processor solve (operational fast solving)
#  inferOverTree!(fg, tree)
# # inferOverTree!(fg, tree, dbg=true)
#

end

#


text:
using MacroTools: @q

function lower_loop(is, out, box, body)
  @q begin
    $out = $box
    $(Expr(:vertex, Loop(), :(($(is...),) -> $body), out))
  end
end

function lower_stencil(is, body)
  @gensym out
  lower_loop(is, out, :(Array{Any,$(length(is))}()), :($out[$(is...)] = $body))
end

function desugar(ex)
  MacroTools.postwalk(ex) do x
    @capture(x, [is__] -> (out_ = box_) => body_) ? lower_loop(is, out, box, body) :
    @capture(x, [is__] -> body_) ? lower_stencil(is, body) :
    @capture(x, c_[is__] = body_) ? :($c = $(lower_stencil(is, body))) :
      x
  end
end

graphm(args, body) = DataFlow.graphm(MacroTools.flatten(body), args = args)

struct Func
  name::Symbol
  graph::IVertex{Any}
end

Base.show(io::IO, f::Func) = print(io, f.name, "::Func")

macro tk(ex)
  @capture(shortdef(desugar(ex)), f_(args__) = body_)
  ex = esc(DataFlow.constructor(graphm(args, body)))
  :($(esc(f)) = Func($(Expr(:quote, f)), $ex); nothing)
end

# Output cleanup

function inline_blocks(ex)
  MacroTools.postwalk(ex) do ex
    @capture(ex, x_ = (body__; y_)) || return ex
    unblock(:($(body...); $x = $y))
  end
end

function rm_aliases(ex)
  aliases = Dict()
  ex = MacroTools.postwalk(ex) do ex
    @capture(ex, (body__; x_Symbol = y_Symbol)) || return ex
    # Avoid for loop assignments
    isexpr(ex, :block) && !(all(isexpr.(ex.args, :(=)))) || return ex
    aliases[x] = y
    :($(body...);)
  end
  MacroTools.postwalk(x -> get(aliases, x, x), ex)
end

code(v::IVertex) =
  v |> DataFlow.syntax |> inline_blocks |>
    rm_aliases |> MacroTools.prettify

max_stars_repo_name: JuliaTagBot/Tokamak

max_issues_repo_name: JuliaTagBot/Tokamak

max_forks_repo_name: FluxML/Tokamak


text:
@testset "667.beautiful-arrangement-ii.jl" begin
    @test beautiful_arrangement(3, 1) == [1, 2, 3]
    @test beautiful_arrangement(3, 2) == [1, 3, 2]
    @test beautiful_arrangement(6, 3) == [1, 2, 3, 6, 4, 5]
    res1 = beautiful_arrangement(60, 32)
    @test res1 ==
          [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22,
           23, 24, 25, 26, 27, 28, 60, 29, 59, 30, 58, 31, 57, 32, 56, 33, 55, 34, 54, 35,
           53, 36, 52, 37, 51, 38, 50, 39, 49, 40, 48, 41, 47, 42, 46, 43, 45, 44]
    @test length(Set(abs(res1[i] - res1[i + 1]) for i in 1:(length(res1) - 1))) == 32
end


text:
import Base: print, show, showcompact
import Base.Dates: value, DateFormat, Slot, slotparse, slotformat

print(io::IO, tz::TimeZone) = print(io, tz.name)
function print(io::IO, tz::FixedTimeZone)
    name = string(tz.name)
    isempty(name) ? print(io, "UTC", tz.offset) : print(io, name)
end
print(io::IO, zdt::ZonedDateTime) = print(io, localtime(zdt), zdt.zone.offset)

showcompact(io::IO, tz::TimeZone) = print(io, string(tz))

function show(io::IO, t::Transition)
    name_str = string(t.zone.name)
    print(io, t.utc_datetime, " ")
    show(io, t.zone.offset)
    !isempty(name_str) && print(io, " (", name_str, ")")
end

function show(io::IO, tz::FixedTimeZone)
    offset_str = "UTC" * offset_string(tz.offset, true)  # Use ISO 8601 for comparision
    name_str = string(tz.name)
    if isempty(name_str)
        print(io, offset_str)
    elseif name_str != offset_str && !(value(tz.offset) == 0 && name_str in ("UTC", "GMT"))
        print(io, name_str, " (UTC", offset_string(tz.offset), ")")
    else
        print(io, name_str)
    end
end

function show(io::IO,tz::VariableTimeZone)
    trans = tz.transitions

    # Retrieve the "modern" time zone transitions. We'll treat the latest transitions as
    # the same as the transitions for `now()` since these future transitions should be
    # based upon the same rules.
    if isnull(tz.cutoff) || length(trans) == 1
        trans = trans[end:end]
    else
        trans = trans[end-1:end]

        # Attempt to show a standard time offset before daylight saving time offset. Sorting
        # should work as long as the DST adjustment is always positive. Fixes differences
        # between the north and south hemispheres.
        sort!(trans, by=el -> el.zone.offset)
    end

    # Show standard time offset before daylight saving time offset.
    print(
        io,
        string(tz.name),
        " (", join(["UTC" * offset_string(t.zone.offset) for t in trans], "/"), ")",
    )
end

show(io::IO,dt::ZonedDateTime) = print(io, string(dt))

# DateTime Parsing
function slotparse(slot::Slot{TimeZone},x,locale)
    if slot.letter == 'z'
        # TODO: Should 'z' only parse numeric UTC offsets? e.g. disallow "UTC-7"
        if ismatch(r"\d", x)
            return FixedTimeZone(x)
        else
            throw(ArgumentError("Time zone offset contains no digits"))
        end
    elseif slot.letter == 'Z'
        # First attempt to create a timezone from the string. An error will be thrown if the
        # time zone is unrecognized.
        tz = TimeZone(x)

        # If the time zone is recognized make sure that it is well-defined. For our purposes
        # we'll treat all abbreviations except for UTC and GMT as ambiguous.
        # e.g. "MST": "Mountain Standard Time" (UTC-7) or "Moscow Summer Time" (UTC+3:31).
        if contains(x, "/") || x in ("UTC", "GMT")
            return tz
        else
            throw(ArgumentError("Time zone is ambiguous"))
        end
    end
end

function slotformat(slot::Slot{TimeZone},zdt::ZonedDateTime,locale)
    if slot.letter == 'z'
        return string(zdt.zone.offset)
    elseif slot.letter == 'Z'
        return string(zdt.zone)  # In most cases will be an abbreviation.
    end
end

# Note: ISOZonedDateTimeFormat is defined in the module __init__ which means that this
# function can not be called from within this module. TODO: Ignore linting for this line
ZonedDateTime(dt::AbstractString,df::DateFormat=ISOZonedDateTimeFormat) = ZonedDateTime(Base.Dates.parse(dt,df)...)
ZonedDateTime(dt::AbstractString,format::AbstractString;locale::AbstractString="english") = ZonedDateTime(dt,DateFormat(format,locale))


text: sec_per_year = 60. * 60. * 24. * 365.25


text:
#=
    Contents: Getters for ConnectionPool.
=#

"""
Gets a connection from the pool if one is available, else returns 0.

Notes:
- Set the target number of connections in the pool under typical usage.
- Set a (finite) peak number of connections that the pool can scale to. When these extra connections are released back to the pool they are deleted until the pool size reaches the target.
- The constraint `target_pool_size <= peak_size` is automatically enforced.
- Reset the target and peak numbers as desired.
- When requesting a connection from the pool:
    - If the target has not been reached, get a connection from the `unoccupied` set if there is one, otherwise create a new one.
    - If the target upper bound has been reached:
        - If the peak number has not been reached, create a new connection (and manually delete it when finished).
        - If the peak number has been reached, wait `ms_wait` milliseconds and try again to acquire a connection.
        - Try a maximum of `n_tries` times to acquire a connection from the pool.
        - If all attempts to acquire a connection fail, return the connection pool's connection prototype, which is a disconnected instance of the connection used to instantiate the connection pool.
"""
function get_connection!(cp::ConnectionPool)
    result = cp.connection_prototype
    if get_n_connections(cp) < get_peak(cp)
	if get_n_unoccupied(cp) == 0                          # If no connection available, create a new one
	    result = new_connection(cp.connection_prototype)
	    push!(cp.occupied, result)
	else                                                  # If connection available, use it
	    result = pop!(cp.unoccupied)
	    push!(cp.occupied, result)
	end
    else    # get_n_connections(cp) == get_peak(cp)
	s = 0.001 * get_n_tries(cp)
	for i = 1:cp.n_tries
	    sleep(s)
	    if get_n_unoccupied(cp) > 0
		result = pop!(cp.unoccupied)
		push!(cp.occupied, result)
		break
	    end
	end
    end
    result
end

function get_n_unoccupied(cp::ConnectionPool)
    length(cp.unoccupied)
end

function get_n_occupied(cp::ConnectionPool)
    length(cp.occupied)
end

"Number of connections, occupied or not."
function get_n_connections(cp::ConnectionPool)
    get_n_unoccupied(cp) + get_n_occupied(cp)
end

function get_target(cp::ConnectionPool)
    cp.target_pool_size
end

function get_peak(cp::ConnectionPool)
    cp.peak_size
end

function get_wait(cp::ConnectionPool)
    cp.ms_wait
end

function get_n_tries(cp::ConnectionPool)
    cp.n_tries
end

# EOF


text:
export update_IminusTVp!

# dense version
function update_IminusTVp!(IminusTVp::Matrix{T}, Π::Matrix{T}, β::T, q0::AbstractVector{T}) where {T<:AbstractFloat}
  n, m = size(IminusTVp)
  n == m || throw(error(DimensionMismatch()))
  (n, n,) == size(Π) || throw(error(DimensionMismatch()))

  # make_eye!(IminusTVp)
  IminusTVp .= zero(T)
  @inbounds for j = 1:n
    for i = 1:n
      IminusTVp[i,j] =       i == j     ?     one(T) - Π[i, j] * β * q0[j]  :       - Π[i, j] * β * q0[j]
    end
  end
end

# sparse version
function update_IminusTVp!(IminusTVp::SparseMatrixCSC{T}, Π::SparseMatrixCSC{T}, β::T, q0::AbstractVector{T}) where {T<:AbstractFloat}
  n, m = size(IminusTVp)
  n == m             || throw(error(DimensionMismatch()))
  (n, n,) == size(Π) || throw(error(DimensionMismatch()))

  Π_rows = rowvals(Π)
  Π_vals = nonzeros(Π)
  IminusTVp_rows = rowvals( IminusTVp)
  IminusTVp_vals = nonzeros(IminusTVp)

  # consider eliminmating this check?
  length(Π_rows) == length(IminusTVp_rows) == length(Π_vals) == length(IminusTVp_vals) || throw(error("Π and IminusTVp not same pattern"))

  @inbounds for j = 1:n
    for nzi in nzrange(IminusTVp, j)
      IminusTVp_vals[nzi] =  j == IminusTVp_rows[nzi]  ?  one(T) - Π_vals[nzi] * β * q0[j]   :   - Π_vals[nzi] * β * q0[j]
    end
  end

end


text:
module MSC
export removeMSC, gapFillMSC, getMSC, getMedSC
using ..Cubes
using ..DAT
using ..Proc
import Statistics: quantile!

function removeMSC(aout,ain,NpY::Integer,tmsc,tnmsc)
    #Start loop through all other variables
    fillmsc(1,tmsc,tnmsc,ain,NpY)
    subtractMSC(tmsc,ain,aout,NpY)
    nothing
end

"""
    removeMSC(c::AbstractCubeData)

Removes the mean annual cycle from each time series of a data cube.

**Input Axis** `Time`axis

**Output Axis** `Time`axis
"""
function removeMSC(c::AbstractCubeData;kwargs...)
    NpY = getNpY(c)
    mapCube(
        removeMSC,
        c,
        NpY,
        zeros(Union{Missing,Float64},NpY),
        zeros(Int,NpY);
        indims  = InDims("Time"),
        outdims = OutDims("Time"),
        kwargs...
    )
end

"""
    gapFillMSC(c::AbstractCubeData)

Fills missing values of each time series in a cube with the mean annual cycle.

**Input Axis** `Time`axis

**Output Axis** `Time`axis
"""
function gapFillMSC(c::AbstractCubeData;kwargs...)
  NpY=getNpY(c)
  mapCube(gapFillMSC,c,NpY,zeros(Union{Float64,Missing},NpY),zeros(Int,NpY);indims=InDims("Time"),outdims=OutDims("Time"),kwargs...)
end

function gapFillMSC(aout::AbstractVector,ain::AbstractVector,NpY::Integer,tmsc,tnmsc)
  fillmsc(1,tmsc,tnmsc,ain,NpY)
  replaceMisswithMSC(tmsc,ain,aout,NpY)
end

"""
    getMSC(c::AbstractCubeData)

Returns the mean annual cycle from each time series.

**Input Axis** `Time`axis

**Output Axis** `MSC`axis

"""
function getMSC(c::AbstractCubeData;kwargs...)
  outdims = OutDims(MSCAxis(getNpY(c)))
  indims = InDims(TimeAxis)
  mapCube(getMSC,c,zeros(Int,getNpY(c)),zeros(Union{Float64,Missing},getNpY(c));indims=indims,outdims=outdims,kwargs...)
end

function getMSC(aout::AbstractVector,ain::AbstractVector,nmsc,fmsc;imscstart::Int=1,NpY=length(aout))
    NpY=length(aout)
    fillmsc(imscstart,fmsc,nmsc,ain,NpY)
    copyto!(aout,fmsc)
end

"Subtracts given msc from input vector"
function subtractMSC(msc::AbstractVector,xin2::AbstractVector,xout2,NpY)
    imsc=1
    ltime=length(xin2)
    for i in 1:ltime
        xout2[i] = xin2[i]-msc[imsc]
        imsc =imsc==NpY ? 1 : imsc+1 # Increase msc time step counter
    end
end

"Replaces missing values with mean seasonal cycle"
function replaceMisswithMSC(msc::AbstractVector,xin::AbstractArray,xout::AbstractArray,NpY::Integer)
  imsc=1
  for i in eachindex(xin)
    if ismissing(xin[i]) && !ismissing(msc[imsc])
      xout[i]=msc[imsc]
    else
      xout[i]=xin[i]
    end
    imsc= imsc==NpY ? 1 : imsc+1 # Increase msc time step counter
  end
end

"""
    getMedMSC(c::AbstractCubeData)

Returns the median annual cycle from each time series.

**Input Axis** `Time`axis

**Output Axis** `MSC`axis
"""
function getMedSC(c::AbstractCubeData;kwargs...)
  outdims = OutDims(MSCAxis(getNpY(c)))
  indims = InDims(TimeAxis)
  mapCube(getMedSC,c;indims=indims,outdims=outdims,kwargs...)
end

function getMedSC(aout::AbstractVector{Union{T,Missing}},ain::AbstractVector) where T
    #Reshape the cube to squeeze unimportant variables
    NpY=length(aout)
    yvec=T[]
    q=[convert(eltype(yvec),0.5)]
    for doy=1:length(aout)
        empty!(yvec)
        for i=doy:NpY:length(ain)
            ismissing(ain[i]) || push!(yvec,ain[i])
        end
        aout[doy] = isempty(yvec) ? missing : quantile!(yvec,q)[1]
    end
    aout
end

"Calculates the mean seasonal cycle of a vector"
function fillmsc(imscstart::Integer,msc::AbstractVector{T1},nmsc::AbstractVector{Int},xin::AbstractVector,NpY) where T1
  imsc=imscstart
  fill!(msc,zero(T1))
  fill!(nmsc,0)
  for v in xin
    if !ismissing(v)
      msc[imsc]  += v
      nmsc[imsc] += 1
    end
    imsc=imsc==NpY ? 1 : imsc+1 # Increase msc time step counter
  end
  for i in 1:NpY msc[i] = nmsc[i] > 0 ? msc[i]/nmsc[i] : missing end # Get MSC by dividing by number of points
end

end


text:
using Documenter, CfGRIB

makedocs(
    modules = [CfGRIB],
    format = Documenter.HTML(
        # Use clean URLs, unless built as a "local" build
        prettyurls = !("local" in ARGS),
        canonical = "https://ecmwf.github.io/cfgrib.jl/stable/",
    ),
    clean = false,
    sitename="CfGRIB.jl",
    doctest = false,
    pages = [
        "Home"    => "index.md",
        "Manual"  => "man/guide.md",
        "Library" => Any[
            # "Public" => "lib/public.md",
            "Internals" => map(
                s -> "lib/internals/$(s)",
                sort(readdir(joinpath(@__DIR__, "src/lib/internals")))
            ),
        ],
    ]
)

deploydocs(
    repo = "github.com/ecmwf/cfgrib.jl",
    devbranch="dev",
    push_preview=true,
)


text:
using NumericalMethodsforEngineers

a = [2. -3. 1.; 1. -1. -2.; 3. 1. -1.]
b = [7., -2., 0.]

(lower, upper) = lufac(a)
println("NumericalMethodsforEngineers lufac(a) - lower:")
println(lower)
println()
println("NumericalMethodsforEngineers lufac(a) - upper:")
println(upper)
println()
@show lower * upper
println()

subfor!(lower, b)
println("Updated RHS:")
println(b)
println()

subbac!(upper, b)
println("Solution vector:")
println(b)
println()

c = a\b

@show round(c, 14)
@show round(a * c, 14) == b
println()

println("Julia lufact(a):")
f = lufact(a);
println(f[:L])
println()
println(f[:U])
println()


text:
include("htmlescape.jl")

include("HTMLCell.jl")
include("HTMLStyle.jl")
include("HTMLStyleElementConverter.jl")

include("render.jl")


text:
using Colors
using FixedPointNumbers
using Images
using FileIO
using DataStructures
using Lazy: @>>
using LRUCache

import Base: display, close, +, convert, promote_rule, convert, push!,
  filter!, length, collect, copy

 # importing solely to allow their use in user code
import Colors: @colorant_str, RGB

export visual, window, font, display, close, @colorant_str, RGB,
  clear_image_cache

@static if is_windows()
  const font_dirs = [".",joinpath(ENV["WINDIR"],"fonts")]
elseif is_apple()
  const font_dirs = [".","/Library/Fonts"]
elseif is_linux()
  const font_dirs = [".","/usr/share/fonts","/usr/local/share/fonts"]
end

function find_font(name,dirs)
  reg = Regex(lowercase(name)*"\\.ttf\$")
  for d in dirs
    for file in readdir(d)
      if ismatch(reg,lowercase(file))
        return joinpath(d,file)
      end
    end
  end
  error("No font matching the name \"$name\" in the directories "*
        join(dirs,", "," and ")*".")
end

abstract type SDLRendered end
abstract type SDLSimpleRendered <: SDLRendered end
timed(r) = 0.0 < display_duration(r) < Inf
visual(x::SDLRendered;kwds...) = update_arguments(x;kwds...)

struct RenderItem
  r::SDLRendered
  delete_at::Float64
end
display_priority(x::RenderItem) = display_priority(x.r)
display_duration(x::RenderItem) = display_duration(x.r)
timed(x::RenderItem) = timed(x.r)

mutable struct DisplayStack
  data::OrderedSet{RenderItem}
  next_change::Float64
end
DisplayStack() = DisplayStack(OrderedSet{RenderItem}(),Inf)

function push!(x::DisplayStack,r::SDLRendered)
  item = RenderItem(r,(timed(r) ? Weber.tick() + display_duration(r) : Inf))
  push!(x.data,item)
  if timed(r)
    x.next_change = min(x.next_change,item.delete_at)
  end
  x
end

function delete_untimed!(x::DisplayStack)
  filter!(r -> display_duration(r) > 0.0,x.data)
end

function delete_timed!(x::DisplayStack)
  filter!(r -> display_duration(r) <= 0.0 || isinf(display_duration(r)),x.data)
  x
end

const change_resolution = 0.001
function delete_expired!(stack::DisplayStack,tick=Weber.tick())
  next_change = Inf
  stack.data = filter!(stack.data) do item
    if item.delete_at + change_resolution <= tick
      next_change = min(next_change,item.delete_at)
      false
    else
      true
    end
  end
  stack.next_change = next_change

  stack
end
length(x::DisplayStack) = length(x.data)
collect(x::DisplayStack) = collect(x.data)
copy(x::DisplayStack) = DisplayStack(copy(x.data),x.next_change)
ischanging(x::DisplayStack,tick) = x.next_change + change_resolution <= tick

abstract type ExperimentWindow end
mutable struct NullWindow <: ExperimentWindow
  w::Cint
  h::Cint
  closed::Bool
end
visual(win::NullWindow,args...;kwds...) = nothing
display(win::NullWindow,r;kwds...) = nothing

mutable struct SDLWindow <: ExperimentWindow
  data::Ptr{Void}
  renderer::Ptr{Void}
  w::Cint
  h::Cint
  closed::Bool
  stack::DisplayStack
end

const SDL_WINDOWPOS_CENTERED = 0x2fff0000
const SDL_WINDOW_FULLSCREEN_DESKTOP = 0x00001001
const SDL_WINDOW_INPUT_GRABBED = 0x00000100
const SDL_WINDOW_INPUT_FOCUS = 0x00000200
const SDL_WINDOW_MOUSE_FOCUS = 0x00000400
const SDL_RENDERER_SOFTWARE = 0x00000001
const SDL_RENDERER_ACCELERATED = 0x00000002
const SDL_RENDERER_PRESENTVSYNC = 0x00000004
const SDL_RENDERER_TARGETTEXTURE = 0x00000008

const display_is_setup = Array{Bool}()
display_is_setup[] = false
"""
    window([width=1024],[height=768];[fullscreen=true],[title="Experiment"],
           [accel=true])

Create a window to which various objects can be rendered. See the `visual`
method.
"""
function window(width=1024,height=768;fullscreen=true,
                title="Experiment",accel=true,null=false)

  if !display_is_setup[] && !null
    setup_display()
    display_is_setup[] = true
  elseif null
    return NullWindow(width,height,false)
  end

  if !ccall((:SDL_SetHint,weber_SDL2),Bool,(Cstring,Cstring),
            "SDL_RENDER_SCALE_QUALITY","1")
    warn("Linear texture filtering not enabled.")
  end

  if ccall((:TTF_Init,weber_SDL2_ttf),Cint,()) == -1
    error("Failed to initialize SDL_ttf: "*TTF_GetError())
  end

  x = y = SDL_WINDOWPOS_CENTERED
  flags = (fullscreen ? SDL_WINDOW_FULLSCREEN_DESKTOP : 0x0) |
    SDL_WINDOW_INPUT_GRABBED | SDL_WINDOW_INPUT_FOCUS | SDL_WINDOW_MOUSE_FOCUS

  win = ccall((:SDL_CreateWindow,weber_SDL2),Ptr{Void},
              (Cstring,Cint,Cint,Cint,Cint,UInt32),
              pointer(title),x,y,width,height,flags)
  if win == C_NULL
    error("Failed to create a window: "*SDL_GetError())
  end

  flags = SDL_RENDERER_ACCELERATED |
    SDL_RENDERER_PRESENTVSYNC |
    SDL_RENDERER_TARGETTEXTURE

  fallback_flags = SDL_RENDERER_SOFTWARE

  rend = ccall((:SDL_CreateRenderer,weber_SDL2),Ptr{Void},
               (Ptr{Void},Cint,UInt32),win,-1,(accel ? flags : fallback_flags))
  if rend == C_NULL
    accel_error = SDL_GetError()
    if accel
      rend = ccall((:SDL_CreateRenderer,weber_SDL2),Ptr{Void},
                   (Ptr{Void},Cint,UInt32),win,-1,fallback_flags)
      if rend == C_NULL
        error("Failed to create a renderer: "*SDL_GetError())
      end
      warn("Failed to create accelerated graphics renderer: "*accel_error)
    else
      error("Failed to create a renderer: "*accel_error())
    end
  end

  wh = Array{Cint}(2)
  ccall((:SDL_GetWindowSize,weber_SDL2),Void,
        (Ptr{Void},Ptr{Cint},Ptr{Cint}),win,pointer(wh,1),pointer(wh,2))
  ccall((:SDL_ShowCursor,weber_SDL2),Void,(Cint,),0)

  x = SDLWindow(win,rend,wh[1],wh[2],false,DisplayStack())
  finalizer(x,x -> (x.closed ? nothing : close(x)))

  x
end

"""
    close(win::SDLWindow)

Closes a visible SDLWindow window.
"""
function close(win::SDLWindow)
  ccall((:SDL_DestroyRenderer,weber_SDL2),Void,(Ptr{Void},),win.renderer)
  ccall((:SDL_DestroyWindow,weber_SDL2),Void,(Ptr{Void},),win.data)
  ccall((:SDL_ShowCursor,weber_SDL2),Void,(Cint,),1)

  win.closed = true
end
close(win::NullWindow) = win.closed = true

function clear(window,color::Color)
  clear(window,convert(RGB{N0f8},color))
end

function clear(window::SDLWindow,color::RGB{N0f8}=colorant"gray")
  ccall((:SDL_SetRenderDrawColor,weber_SDL2),Void,
        (Ptr{Void},UInt8,UInt8,UInt8),window.renderer,
        reinterpret(UInt8,red(color)),
        reinterpret(UInt8,green(color)),
        reinterpret(UInt8,blue(color)))
  ccall((:SDL_RenderClear,weber_SDL2),Void,(Ptr{Void},),window.renderer)
  nothing
end
clear(win::NullWindow,color) = nothing

mutable struct SDLFont
  data::Ptr{Void}
  color::RGBA{N0f8}
end

"""
    font(name,size,[dirs=os_default],[color=colorant"white"])

Creates an `SDLFont` object to be used for for rendering text as an image.

By default this function looks in the current directory and then an os specific
default font directory for a font with the given name (case insensitive). You
can specify a different list of directories using the `dirs` parameter.

"""
function font(name::String,size;dirs=font_dirs,color=colorant"white")
  file = find_font(name,dirs)
  font = ccall((:TTF_OpenFont,weber_SDL2_ttf),Ptr{Void},(Cstring,Cint),
               pointer(file),size)
  if font == C_NULL
    error("Failed to load the font $file: "*TTF_GetError())
  end

  x = SDLFont(font,color)
  finalizer(x,x -> ccall((:TTF_CloseFont,weber_SDL2_ttf),Void,
                         (Ptr{Void},),x.data))
  x
end

const forever = -1

"""
    visual(obj,[duration=0s],[priority=0],keys...)

Render an object, allowing `display` to show the object in current experiment's
window.

# Arguments
* duration: A positive duration means the object is
  displayed for the given duration, otherwise the object displays until a new
  object is displayed.
* priority: Higher priority objects are always visible above lower priority ones.
  Newer objects display over same-priority older objects.

If coordinates are used they are in units of half screen widths (for x)
and heights (for y), with (0,0) at the center of the screen.

!!! note

    By using using the `+` operator, multiple visual objects can be composed
    into one object, so that they are displayed together
"""
visual(x,args...;keys...) = visual(win(get_experiment()),x,args...;keys...)
visual(win::SDLWindow,x,args...;keys...) =
  throw(MethodError(visual,typeof((win,x))))

function draw(r::SDLRendered)
  draw(win(get_experiment()),r)
end

struct SDLClear <: SDLSimpleRendered
  color::Color
  duration::Float64
  priority::Float64
end
display_duration(clear::SDLClear) = clear.duration
display_priority(clear::SDLClear) = clear.priority
draw(window::SDLWindow,cl::SDLClear) = clear(window,cl.color)
function update_arguments(cl::SDLClear;duration=cl.duration*s,
                          priority=cl.priority,kwds...)
  SDLClear(cl.color,ustrip(inseconds(duration)),priority)
end

"""
    visual(color,[duration=0s],[priority=0])

Render a color, across the entire screen.
"""
function visual(window::SDLWindow,color::Color;duration=0s,priority=0)
  SDLClear(color,ustrip(inseconds(duration)),priority)
end

immutable SDLRect
  x::Cint
  y::Cint
  w::Cint
  h::Cint
end

function as_screen_coordinates(window,x,y,w,h)
  max(0,min(window.w,round(Cint,window.w/2 + x*window.w/4 - w / 2))),
  max(0,min(window.h,round(Cint,window.h/2 - y*window.h/4 - h / 2)))
end

function as_screen_coordinate_x(window,x,w)
  max(0,min(window.w,round(Cint,window.w/2 + x*window.w/4 - w / 2)))
end
function as_screen_coordinate_y(window,y,w)
  max(0,min(window.h,round(Cint,window.h/2 - y*window.h/4 - h / 2)))
end

function update_arguments(rect::SDLRect;w=NaN,h=NaN,x=NaN,y=NaN,kwds...)
  if !isnan(x) || !isnan(y) || !isnan(w) || !isnan(h)
    if isnan(w)
      w = rect.w
    end
    if isnan(h)
      h = rect.h
    end

    if !isnan(x)
      newx = as_screen_coordinate_x(window,x,w)
    else
      newx = rect.x
    end

    if !isnan(y)
      newy = as_screen_coordinate_y(window,y,h)
    else
      newy = rect.y
    end

    SDLRect(x,y,w,h)
  else
    rect
  end
end

abstract type SDLTextured <: SDLSimpleRendered end

function draw(window::SDLWindow,texture::SDLTextured)
  ccall((:SDL_RenderCopy,weber_SDL2),Void,
        (Ptr{Void},Ptr{Void},Ptr{Void},Ptr{SDLRect}),
        window.renderer,data(texture),C_NULL,Ref(rect(texture)))
  nothing
end

mutable struct SDLText <: SDLTextured
  str::String
  data::Ptr{Void}
  rect::SDLRect
  duration::Float64
  priority::Float64
  color::Color
end
display_duration(text::SDLText) = text.duration
display_priority(text::SDLText) = text.priority
data(text::SDLText) = text.data
rect(text::SDLText) = text.rect
function update_arguments(text::SDLText;w=NaN,h=NaN,color=nothing,
                          duration=text.duration*s,priority=text.priority,
                          kwds...)
  if !isnan(w) || !isnan(h)
    error("Cannot update the height or width of text on display. These",
          "properties are automatically determiend by the original string",
          " used to render the text")
  end

  if color != nothing
    error("You cannot change the color of text on display. Create"*
          " a new text object using the desired color.")
  end

  rect = update_arguments(text.rect;kwds...)
  SDLText(text.str,text.data,rect,ustrip(inseconds(duration)),priority,text.color)
end

fonts = Dict{Tuple{String,Int},SDLFont}()

image_formats = [
  "BMP",
  "AVI",
  "CRW",
  "CUR",
  "DCX",
  "DOT",
  "EPS",
  "GIF",
  "HDR",
  "ICO",
  "INFO",
  "JP2",
  "JPEG",
  "PCX",
  "PDB",
  "PDF",
  "PGM",
  "PNG",
  "PSD",
  "RGB",
  "TIFF",
  "WMF",
  "WPG",
  "TGA"
]
function isimage(str::String)
  ismatch(r".*\.(\w{3,4})$",str) &&
    uppercase(match(r".*\.(\w{3,4})$",str)[1]) in image_formats
end

"""
    visual(str::String, [font=nothing], [font_name="arial"], [size=32],
           [color=colorant"white"],
           [wrap_width=0.8],[clean_whitespace=true],[x=0],[y=0],[duration=0s],
           [priority=0])

Render the given string as an image that can be displayed. An optional second
argument can specify a font, loaded using the `font` function.

!!! note "Strings treated as files..."

    If the string passed refers to an image file--becasue the string ends in a
    file type, like .bmp or .png---it will be treated as an image to be loaded
    and displayed, rather than as a string to be printed to the screen.
    Refer to the documentation of `visual` for image objects.

# Arguments
* wrap_width: the proporition of the screen that the text can utilize
  before wrapping.
* clean_whitespace: if true, replace all consecutive white space with a single
  space.
"""
function visual(window::SDLWindow,str::String,cache=true;keys...)
  if isimage(str)
    image_cache(cache,str) do
      visual(window,load(str),false;keys...)
    end
  else
    visual_text(window,str;keys...)
  end
end

function visual_text(window::SDLWindow,str::String;
                     font=nothing,font_name="arial",size=32,info...)
  if font == nothing
    f = get!(fonts,(font_name,size)) do
      Weber.font(font_name,size)
    end
    visual(window,str,f;info...)
  else
    if font_name != "arial" || size != 32
      error("Cannot specify both a font object and font_name or size. Either ",
            "use a font object, or specify the font using a name and size.")
    end
    visual(window,str,font;info...)
  end
end

function visual(window::SDLWindow,str::String,font::SDLFont;
                color::RGB{N0f8}=colorant"white",
                wrap_width=0.8,clean_whitespace=true,x=0,y=0,
                duration=0s,priority=0)
  if clean_whitespace
    str = replace(str,r"^\s+","")
    str = replace(str,r"\s+"," ")
  end
  visual(window,x,y,font,color,round(UInt32,window.w*wrap_width),str,
         inseconds(duration),priority)
end

const w_ptr = 0x0000000000000010 # icxx"offsetof(SDL_Surface,w);"
const h_ptr = 0x0000000000000014 # icxx"offsetof(SDL_Surface,h);"

const text_cache = LRU{Tuple{String,UInt32,SDLFont},SDLText}(256)

function visual(window::SDLWindow,x::Real,y::Real,font::SDLFont,color::RGB{N0f8},
                wrap_width::UInt32,str::String,duration=0s,priority=0)
  get!(text_cache,(str,wrap_width,font)) do
    surface = ccall((:TTF_RenderUTF8_Blended_Wrapped,weber_SDL2_ttf),Ptr{Void},
                    (Ptr{Void},Cstring,RGBA{N0f8},UInt32),
                    font.data,pointer(str),color,wrap_width)
    if surface == C_NULL
      error("Failed to render text: "*TTF_GetError())
    end

    texture = ccall((:SDL_CreateTextureFromSurface,weber_SDL2),Ptr{Void},
                    (Ptr{Void},Ptr{Void}),window.renderer,surface)
    if texture == C_NULL
      error("Failed to create texture for render text: "*SDL_GetError())
    end

    w = at(surface,Cint,w_ptr)
    h = at(surface,Cint,h_ptr)

    xint,yint = as_screen_coordinates(window,x,y,w,h)

    result = SDLText(str,texture,SDLRect(xint,yint,w,h),
                     ustrip(inseconds(duration)),priority,color)
    finalizer(result,x ->
              ccall((:SDL_DestroyTexture,weber_SDL2),Void,(Ptr{Void},),x.data))
    ccall((:SDL_FreeSurface,weber_SDL2),Void,(Ptr{Void},),surface)

    result
  end
end

mutable struct SDLImage <: SDLTextured
  data::Ptr{Void}
  img::Array{RGBA{N0f8}}
  rect::SDLRect
  duration::Float64
  priority::Float64
end
display_duration(img::SDLImage) = img.duration
display_priority(img::SDLImage) = img.priority
data(img::SDLImage) = img.data
rect(img::SDLImage) = img.rect
function update_arguments(img::SDLImage;w=NaN,h=NaN,duration=img.duration*s,
                          priority=img.priority,kwds...)
  if !isnan(w) || !isnan(h)
    error("Cannot update the height or width of an image on display. These",
          "properties are automatically determiend by the original image.")
  end

  rect = update_arguments(img.rect;kwds...)
  SDLImage(img.data,img.img,rect,ustrip(inseconds(duration)),priority)
end

const convert_cache = LRU{UInt,Array{RGBA{N0f8}}}(256)
const _image_cache = LRU{UInt,SDLImage}(256)

function clear_image_cache()
  empty!(convert_cache)
  empty!(_imge_cache)
end

function image_cache(fn,usecache,x)
  if usecache
    get!(fn,_image_cache,object_id(x))
  else
    fn()
  end
end

function image_cache(fn,usecache::Union{File,String},x)
  if usecache
    get!(fn,_image_cache,x)
  else
    fn()
  end
end

"""
    visual(img, [x=0],[y=0],[duration=0s],[priority=0])

Prepare the color or gray scale image to be displayed to the screen.

For a string or file reference, this loads and prepares for display the given
image file. For an array this utilizes all the conventions in the `Images`
package for representing images. Internally, real-number 2d arrays are
interpreted as gray scale images, and real-number 3d arrays as an RGB image or
RGBA image, depending on whether size(img,1) is of size 3 or 4. A 3d array with
a size(img,1) ∉ [3,4] results in an error.
"""
function visual(window::SDLWindow,img::Array{<:AbstractFloat},cache=true;keys...)
  image_cache(cache,img) do
    converted = if length(size(img)) == 3
      if size(img,1) == 3
        n0f8.(colorview(RGB,img))
      elseif size(img,1) == 4
        n0f8.(colorview(RGBA,img))
      else
        error("Could not interpret array of size $(size(imge)) as a color image.")
      end
    elseif length(size(img)) == 2
      n0f8.(colorview(Gray,img))
    end
    visual(window,converted,false;keys...)
  end
end

function visual(window::SDLWindow,img::Array,cache=true;keys...)
  image_cache(cache,img) do
    visual(window,convert(RGBA,n0f8.(img)),false;keys...)
  end
end

function visual(window::SDLWindow,img::Array{RGBA{N0f8}},cache=true;
                x=0,y=0,duration=0s,priority=0)
  image_cache(cache,img) do
    surface = ccall((:SDL_CreateRGBSurfaceFrom,weber_SDL2),Ptr{Void},
                    (Ptr{Void},Cint,Cint,Cint,Cint,UInt32,UInt32,UInt32,UInt32),
                    pointer(copy(img')),size(img,2),size(img,1),32,
                    4size(img,2),0x000000ff,0x0000ff00,0x00ff0000,0xff000000)
    if surface == C_NULL
      error("Failed to create image surface: "*SDL_GetError())
    end

    texture = ccall((:SDL_CreateTextureFromSurface,weber_SDL2),Ptr{Void},
                    (Ptr{Void},Ptr{Void}),window.renderer,surface)
    if texture == C_NULL
      error("Failed to create texture from image surface: "*SDL_GetError())
    end

    h,w = size(img)
    xint,yint = as_screen_coordinates(window,x,y,w,h)

    result = SDLImage(texture,img,SDLRect(xint,yint,w,h),
                      ustrip(inseconds(duration)),priority)
    finalizer(result,x ->
              ccall((:SDL_DestroyTexture,weber_SDL2),Void,(Ptr{Void},),x.data))
    ccall((:SDL_FreeSurface,weber_SDL2),Void,(Ptr{Void},),surface)

    result
  end
end

# Below is a scratch implementation that uses Compose rendering. My
# implementation currently throws a mysterious MacOS GUI error on my machine. I
# suspect there is some interaction between Cairo and SDL that is causing the
# problem, but I haven't narrowed it down yet.

# Also note that this currently wont' work well because Pkg.installed isn't a
# good way to handle conditional code. Julia doesn't yet have a better
# solution for this that I know of.

#=
if Pkg.installed("Compose") != nothing && Pkg.installed("Cairo") != nothing
  import Compose
  import Cairo

  const sdl_compose_implemented = true

  type SDLComposed <: SDLTextured
    data::Ptr{Void}
    surface::Cairo.CairoSurface
    rect::SDLRect
    duration::Float64
    priority::Float64
  end
  display_duration(img::SDLComposed) = img.duration
  display_priority(img::SDLComposed) = img.priority
  data(img::SDLComposed) = img.data
  rect(img::SDLComposed) = img.rect

"""
  visual(comp::Compose.Context,[x=0],[y=0],[w=1],[h=1],
                           [dpi=72],[priority=0],[duration=0s])

Renders the given Compose.jl object so it can be displayed on screen.
Width and height are specified as a proportion of the full width and height.
"""
  function visual(comp::Compose.Context;keys...)
    visual(get_experiment().win,comp;keys...)
  end
  function visual(win::SDLWindow,comp::Compose.Context;
                  x=0,y=0,w=1,h=1,dpi=72)
    w,h = round(Cint,w*win.w),round(Cint,h*win.h)
    png = Compose.PNG(w,h,dpi)
    Compose.draw(png,comp)
    cairo_surface = png.surface

    surface = ccall((:SDL_CreateRGBSurfaceFrom,weber_SDL2),Ptr{Void},
                    (Ptr{Void},Cint,Cint,Cint,Cint,UInt32,UInt32,UInt32,UInt32),
                    cairo_surface.ptr,w,h,8,4w,0,0,0,0)
    if surface == C_NULL
      error("Failed to create surface for Compose.Context: "*SDL_GetError())
    end

    texture = ccall((:SDL_CreateTextureFromSurface,weber_SDL2),Ptr{Void},
                    (Ptr{Void},Ptr{Void}),window.renderer,surface)
    if texture == C_NULL
      error("Failed to create texture from Compose.Context surface: "*SDL_GetError())
    end

    xint,yint = as_screen_coordinates(window,x,y,w,h)

    result = SDLComposed(texture,cairo_surface,SDLRect(xint,yint,w,h),
                         duraiton,priority)
    finalizer(result,x ->
              ccall((:SDL_DestroyTexture,weber_SDL2),Void,(Ptr{Void},),x.data))

    result
  end
elseif Pkg.installed("Compose") != nothing
  import Compose

  const sdl_compose_implemented = false
  function visual(win::SDLWindow,comp::Compose.Context;keys...)
    error("To render Compose.jl objects you must have Cairo.jl installed!\n"*
          "Please call Pkg.add(\"Cairo\") from the command line")
  end
else
  const sdl_compose_implemented = false
end
=#

function show_drawn(window::SDLWindow)
  ccall((:SDL_RenderPresent,weber_SDL2),Void,(Ptr{Void},),window.renderer)
  nothing
end
show_drawn(win::NullWindow) = nothing

const SDL_INIT_VIDEO = 0x00000020
function setup_display()
  init = ccall((:SDL_Init,weber_SDL2),Cint,(UInt32,),SDL_INIT_VIDEO)
  if init < 0
    error("Failed to initialize SDL: "*SDL_GetError())
  end
  if !sdl_is_setup[]
    sdl_is_setup[] = true
    atexit(() -> ccall((:SDL_Quit,weber_SDL2),Void,()))
  end

  #=
  if !sdl_compose_implemented && Pkg.installed("Compose") != nothing
    warn("The Compose package was installed after Weber. "*
         "To render Compose objects using Weber you will need to run\n "*
         "rm(Pkg.dir(\"Weber\"),recursive=true,force=true); "*
         "Pkg.add(\"Weber\")\n before you call `using Weber`.")
  end
  =#
end

"""
    display(r::SDLRendered;kwds...)

Displays anything rendered by `visual` onto the current experiment window.

Any keyword arguments, available from [`visual`](@ref) are also available here.
They overload the arguments as specified during visual (but do not change them).

    display(x;kwds...)

Short-hand for `display(visual(x);kwds...)`. This is the most common way to use
display. For example:

    moment(0.5s,display,"Hello, World!")

This code will show the text "Hello, World!" on the screen 0.5 seconds after
the start of the previous moment.

!!! warning

    Assuming your hardware and video drivers permit it, `display` sycnrhonizes
    to the screen refresh rate so long as the experiment window uses accelerated
    graphics (true by default). The display of a visual can be no more accurate
    than that permitted by this refresh rate. In particular, display can block
    for up to the length of an entire refresh cycle. If you want accurate timing
    in your experiment, make sure that there is nothing you want to occur
    immediately after calling display. If you want to display multiple visuals
    at once remember that you can compose visuals using the `+` operator, do not
    call display multiple times and expect these visual to all display at the
    same time (also note that the default behavior of visuals is to disappear
    when the next visual is shown).

"""
function display(r::SDLRendered;kwds...)
  if in_experiment() && !experiment_running()
    error("You cannot call `display` during experiment `setup`. During `setup`",
          " you should add `dispaly` to a trial (e.g. ",
          "`addtrial(moment(display,my_visual))`).")
  end
  display(win(get_experiment()),r;kwds...)
end

"""
    display(fn::Function;kwds...)

Display the visual returned by calling `fn`.
"""
display(fn::Function;kwds...) = display(fn();kwds...)

update_arguments(r) = r
function display(window::SDLWindow,r::SDLRendered;kwds...)
  r = update_arguments(r;kwds...)
  update_stack!(window,r)

  # if there's no experiment, handle display duration
  if timed(r) && !in_experiment()
    eventually_remove!(window,stack,r)
    warn("calling display outside of an experiment")
  end
end

function update_stack!(window,r)
  update_stack_helper!(window,r)
  draw_stack(window)
end

function eventually_remove!(window::SDLWindow,stack::DisplayStack,r::SDLRendered)
  Timer(duration(r)) do t
    window.stack = delete_expired!(window.stack,precise_time())
    draw_stack(window)
  end
end

function display(w::SDLWindow,r)
  if experiment_running()
    warn("Visual was not precomputed! To minimize latency* call "*
         "`x = visual(obj)` before running an experiment* then call"*
         " `display(x)` during the experiment.",moment_trace_string())
  end
  display(w,visual(w,r))
end

function update_stack_helper!(window,r::SDLRendered)
  delete_untimed!(window.stack)
  push!(window.stack,r)
end

function draw_stack(window::SDLWindow)
  clear(window)
  for item in sort(collect(window.stack),by=display_priority,alg=MergeSort)
    draw(window,item.r)
  end
  show_drawn(window)
end

function refresh_display(window::SDLWindow,tick=Weber.tick())
  if ischanging(window.stack,tick)
    window.stack = delete_expired!(window.stack,tick)
    draw_stack(window)
  end
end

refresh_display(window::NullWindow) = nothing

struct SDLCompound <: SDLRendered
  data::Array{SDLRendered}
end
function +(a::SDLSimpleRendered,b::SDLSimpleRendered)
  SDLCompound([a,b])
end
function +(a::SDLCompound,b::SDLCompound)
  SDLCompound(vcat(a.data,b.data))
end
+(a::SDLRendered,b::SDLRendered) = +(promote(a,b)...)
promote_rule{T <: SDLSimpleRendered}(::Type{SDLCompound},::Type{T}) = SDLCompound
convert(::Type{SDLCompound},x::SDLSimpleRendered) = SDLCompound([x])
function update_stack_helper!(window,rs::SDLCompound)
  delete_untimed!(window.stack)
  for r in rs.data
    push!(window.stack,r)
  end
end
function eventually_remove!(window::SDLWindow,stack::DisplayStack,
                            rs::SDLCompound)
  for r in rs
    eventually_remove!(window,stack,r)
  end
end
timed(r::SDLCompound) = any(timed,r.data)

function update_arguments(rs::SDLCompound;kwds...)
  SDLCompound(map(r -> update_arguments(r;kwds...),rs.data))
end

struct RestoreDisplay <: SDLRendered
  x::DisplayStack
end
function update_stack_helper!(window,restore::RestoreDisplay)
  # only restore objects that aren't timed. Otherwise the timed objects will
  # never get deleted. There is no well defined way to set their duration again
  # and remove them, since some unknown amount of time has passed since they
  # started being displayed. TODO: this is no longer true,
  # their duration could be calculated
  window.stack = delete_timed!(restore.x)
end

const saved_display = Array{DisplayStack}()
saved_display[] = DisplayStack()
function save_display(window::SDLWindow)
  saved_display[] = copy(window.stack)
end
save_display(win::NullWindow) = nothing

function restore_display(window::SDLWindow)
  update_stack!(window,RestoreDisplay(saved_display[]))
end
restore_display(win::NullWindow) = nothing

function focus(window::SDLWindow)
  ccall((:SDL_RaiseWindow,weber_SDL2),Void,(Ptr{Void},),window.data)
end
focus(win::NullWindow) = nothing