text:
"""
Given a flatbuffers metadata type definition (a Field instance from Schema.fbs),
translate to the appropriate Julia storage eltype
"""
function juliaeltype end

"""
Given a FlatBuffers.Builder and a Julia column eltype,
Write the field.type flatbuffer definition
"""
function arrowtype end

"""
There are a couple places when writing arrow buffers where
we need to write a "dummy" value; it doesn't really matter
what we write, but we need to write something of a specific
type. So each supported writing type needs to define `default`.
"""
function default end

default(T) = zero(T)

finaljuliatype(T) = T
finaljuliatype(::Type{Missing}) = Missing
finaljuliatype(::Type{Union{T, Missing}}) where {T} = Union{Missing, finaljuliatype(T)}

function juliaeltype(f::Meta.Field)
    T = juliaeltype(f, f.type)
    if f.custom_metadata !== nothing
        fm = Dict(kv.key => kv.value for kv in f.custom_metadata)
    else
        fm = nothing
    end
    return (f.nullable ? Union{T, Missing} : T), fm
end

juliaeltype(f::Meta.Field, ::Meta.Null) = Missing

function arrowtype(b, ::Type{Missing})
    Meta.nullStart(b)
    return Meta.Null, Meta.nullEnd(b), nothing
end

function juliaeltype(f::Meta.Field, int::Meta.Int)
    if int.is_signed
        if int.bitWidth == 8
            Int8
        elseif int.bitWidth == 16
            Int16
        elseif int.bitWidth == 32
            Int32
        elseif int.bitWidth == 64
            Int64
        elseif int.bitWidth == 128
            Int128
        else
            throw(InvalidMetadataError("$int is not valid arrow type metadata"))
        end
    else
        if int.bitWidth == 8
            UInt8
        elseif int.bitWidth == 16
            UInt16
        elseif int.bitWidth == 32
            UInt32
        elseif int.bitWidth == 64
            UInt64
        elseif int.bitWidth == 128
            UInt128
        else
            throw(InvalidMetadataError("$int is not valid arrow type metadata"))
        end
    end
end

function arrowtype(b, ::Type{T}) where {T <: Integer}
    Meta.intStart(b)
    Meta.intAddBitWidth(b, Int32(8 * sizeof(T)))
    Meta.intAddIsSigned(b, T <: Signed)
    return Meta.Int, Meta.intEnd(b), nothing
end

# primitive types
function juliaeltype(f::Meta.Field, fp::Meta.FloatingPoint)
    if fp.precision == Meta.Precision.HALF
        Float16
    elseif fp.precision == Meta.Precision.SINGLE
        Float32
    elseif fp.precision == Meta.Precision.DOUBLE
        Float64
    end
end

function arrowtype(b, ::Type{T}) where {T <: AbstractFloat}
    Meta.floatingPointStart(b)
    Meta.floatingPointAddPrecision(b, T === Float16 ? Meta.Precision.HALF : T === Float32 ? Meta.Precision.SINGLE : Meta.Precision.DOUBLE)
    return Meta.FloatingPoint, Meta.floatingPointEnd(b), nothing
end

juliaeltype(f::Meta.Field, b::Union{Meta.Utf8, Meta.LargeUtf8}) = String

function arrowtype(b, ::Type{String})
    # To support LargeUtf8, we'd need a way to flag/pass a max length from user/actual data
    Meta.utf8Start(b)
    return Meta.Utf8, Meta.utf8End(b), nothing
end

default(::Type{String}) = ""

datasizeof(x) = sizeof(x)
datasizeof(x::AbstractVector) = sum(datasizeof, x)

juliaeltype(f::Meta.Field, b::Union{Meta.Binary, Meta.LargeBinary}) = Vector{UInt8}

function arrowtype(b, ::Type{Vector{UInt8}})
    # To support LargeBinary, we'd need a way to flag/pass a max length from user/actual data
    Meta.binaryStart(b)
    return Meta.Binary, Meta.binaryEnd(b), nothing
end

function default(::Type{A}) where {A <: AbstractVector{T}} where {T}
    a = similar(A, 1)
    a[1] = default(T)
    return a
end

juliaeltype(f::Meta.Field, x::Meta.FixedSizeBinary) = NTuple{Int(x.byteWidth), UInt8}

function arrowtype(b, ::Type{NTuple{N, UInt8}}) where {N}
    Meta.fixedSizeBinaryStart(b)
    Meta.fixedSizeBinaryAddByteWidth(b, Int32(N))
    return Meta.FixedSizeBinary, Meta.fixedSizeBinaryEnd(b), nothing
end

default(::Type{NTuple{N, T}}) where {N, T} = ntuple(i -> default(T), N)
default(::Type{T}) where {T <: Tuple} = Tuple(default(fieldtype(T, i)) for i = 1:fieldcount(T))

juliaeltype(f::Meta.Field, x::Meta.Bool) = Bool

function arrowtype(b, ::Type{Bool})
    Meta.boolStart(b)
    return Meta.Bool, Meta.boolEnd(b), nothing
end

struct Decimal{P, S}
    bytes::NTuple{16, UInt8}
end

Base.zero(::Type{Decimal{P, S}}) where {P, S} = Decimal{P, S}(ntuple(i->0x00, 16))

function juliaeltype(f::Meta.Field, x::Meta.Decimal)
    return Decimal{x.precision, x.scale}
end

function arrowtype(b, ::Type{Decimal{P, S}}) where {P, S}
    Meta.decimalStart(b)
    Meta.decimalAddPrecision(b, Int32(P))
    Meta.decimalAddScale(b, Int32(S))
    return Meta.Decimal, Meta.decimalEnd(b), nothing
end

abstract type ArrowTimeType end
Base.write(io::IO, x::ArrowTimeType) = Base.write(io, x.x)

struct Date{U, T} <: ArrowTimeType
    x::T
end

Base.zero(::Type{Date{U, T}}) where {U, T} = Date{U, T}(T(0))

bitwidth(x::Meta.DateUnit) = x == Meta.DateUnit.DAY ? Int32 : Int64
Date{Meta.DateUnit.DAY}(days) = Date{Meta.DateUnit.DAY, Int32}(Int32(days))
Date{Meta.DateUnit.MILLISECOND}(ms) = Date{Meta.DateUnit.MILLISECOND, Int64}(Int64(ms))

juliaeltype(f::Meta.Field, x::Meta.Date) = Date{x.unit, bitwidth(x.unit)}
finaljuliatype(::Type{Date{Meta.DateUnit.DAY, Int32}}) = Dates.Date
Base.convert(::Type{Dates.Date}, x::Date{Meta.DateUnit.DAY, Int32}) = Dates.Date(Dates.UTD(Int64(x.x + UNIX_EPOCH_DATE)))
finaljuliatype(::Type{Date{Meta.DateUnit.MILLISECOND, Int64}}) = Dates.DateTime
Base.convert(::Type{Dates.DateTime}, x::Date{Meta.DateUnit.MILLISECOND, Int64}) = Dates.DateTime(Dates.UTM(Int64(x.x + UNIX_EPOCH_DATETIME)))

function arrowtype(b, ::Type{Date{U, T}}) where {U, T}
    Meta.dateStart(b)
    Meta.dateAddUnit(b, U)
    return Meta.Date, Meta.dateEnd(b), nothing
end

const UNIX_EPOCH_DATE = Dates.value(Dates.Date(1970))
Base.convert(::Type{Date{Meta.DateUnit.DAY, Int32}}, x::Dates.Date) = Date{Meta.DateUnit.DAY, Int32}(Int32(Dates.value(x) - UNIX_EPOCH_DATE))

const UNIX_EPOCH_DATETIME = Dates.value(Dates.DateTime(1970))
Base.convert(::Type{Date{Meta.DateUnit.MILLISECOND, Int64}}, x::Dates.DateTime) = Date{Meta.DateUnit.MILLISECOND, Int64}(Int64(Dates.value(x) - UNIX_EPOCH_DATETIME))

struct Time{U, T} <: ArrowTimeType
    x::T
end

Base.zero(::Type{Time{U, T}}) where {U, T} = Time{U, T}(T(0))

bitwidth(x::Meta.TimeUnit) = x == Meta.TimeUnit.SECOND || x == Meta.TimeUnit.MILLISECOND ? Int32 : Int64
Time{U}(x) where {U <: Meta.TimeUnit} = Time{U, bitwidth(U)}(bitwidth(U)(x))

juliaeltype(f::Meta.Field, x::Meta.Time) = Time{x.unit, bitwidth(x.unit)}
finaljuliatype(::Type{<:Time}) = Dates.Time
periodtype(U::Meta.TimeUnit) = U === Meta.TimeUnit.SECOND ? Dates.Second :
                               U === Meta.TimeUnit.MILLISECOND ? Dates.Millisecond :
                               U === Meta.TimeUnit.MICROSECOND ? Dates.Microsecond : Dates.Nanosecond
Base.convert(::Type{Dates.Time}, x::Time{U, T}) where {U, T} = Dates.Time(Dates.Nanosecond(Dates.tons(periodtype(U)(x.x))))

function arrowtype(b, ::Type{Time{U, T}}) where {U, T}
    Meta.timeStart(b)
    Meta.timeAddUnit(b, U)
    return Meta.Time, Meta.timeEnd(b), nothing
end

Base.convert(::Type{Time{Meta.TimeUnit.NANOSECOND, Int64}}, x::Dates.Time) = Time{Meta.TimeUnit.NANOSECOND, Int64}(Dates.value(x))

struct Timestamp{U, TZ} <: ArrowTimeType
    x::Int64
end

Base.zero(::Type{Timestamp{U, T}}) where {U, T} = Timestamp{U, T}(Int64(0))

function juliaeltype(f::Meta.Field, x::Meta.Timestamp)
    return Timestamp{x.unit, x.timezone === nothing ? nothing : Symbol(x.timezone)}
end

finaljuliatype(::Type{Timestamp{U, nothing}}) where {U} = Dates.DateTime
Base.convert(::Type{Dates.DateTime}, x::Timestamp{U, nothing}) where {U} =
    Dates.DateTime(Dates.UTM(Int64(Dates.toms(periodtype(U)(x.x)) + UNIX_EPOCH_DATETIME)))

function arrowtype(b, ::Type{Timestamp{U, TZ}}) where {U, TZ}
    tz = TZ !== nothing ? FlatBuffers.createstring!(b, String(TZ)) : FlatBuffers.UOffsetT(0)
    Meta.timestampStart(b)
    Meta.timestampAddUnit(b, U)
    Meta.timestampAddTimezone(b, tz)
    return Meta.Timestamp, Meta.timestampEnd(b), nothing
end

struct Interval{U, T} <: ArrowTimeType
    x::T
end

Base.zero(::Type{Interval{U, T}}) where {U, T} = Interval{U, T}(T(0))

bitwidth(x::Meta.IntervalUnit) = x == Meta.IntervalUnit.YEAR_MONTH ? Int32 : Int64
Interval{Meta.IntervalUnit.YEAR_MONTH}(x) = Interval{Meta.IntervalUnit.YEAR_MONTH, Int32}(Int32(x))
Interval{Meta.IntervalUnit.DAY_TIME}(x) = Interval{Meta.IntervalUnit.DAY_TIME, Int64}(Int64(x))

function juliaeltype(f::Meta.Field, x::Meta.Interval)
    return Interval{x.unit, bitwidth(x.unit)}
end

function arrowtype(b, ::Type{Interval{U, T}}) where {U, T}
    Meta.intervalStart(b)
    Meta.intervalAddUnit(b, U)
    return Meta.Interval, Meta.intervalEnd(b), nothing
end

struct Duration{U} <: ArrowTimeType
    x::Int64
end

Base.zero(::Type{Duration{U}}) where {U} = Duration{U}(Int64(0))

function juliaeltype(f::Meta.Field, x::Meta.Duration)
    return Duration{x.unit}
end

finaljuliatype(::Type{Duration{U}}) where {U} = periodtype(U)
Base.convert(::Type{P}, x::Duration{U}) where {P <: Dates.Period, U} = P(periodtype(U)(x.x))

function arrowtype(b, ::Type{Duration{U}}) where {U}
    Meta.durationStart(b)
    Meta.durationAddUnit(b, U)
    return Meta.Duration, Meta.durationEnd(b), nothing
end

arrowperiodtype(P) = Meta.TimeUnit.SECOND
arrowperiodtype(::Type{Dates.Millisecond}) = Meta.TimeUnit.MILLISECOND
arrowperiodtype(::Type{Dates.Microsecond}) = Meta.TimeUnit.MICROSECOND
arrowperiodtype(::Type{Dates.Nanosecond}) = Meta.TimeUnit.NANOSECOND

Base.convert(::Type{Duration{U}}, x::Dates.Period) where {U} = Duration{U}(Dates.value(periodtype(U)(x)))

# nested types; call juliaeltype recursively on nested children
function juliaeltype(f::Meta.Field, list::Union{Meta.List, Meta.LargeList})
    T, _ = juliaeltype(f.children[1])
    return Vector{T}
end

# arrowtype will call fieldoffset recursively for children
function arrowtype(b, ::Type{Vector{T}}) where {T}
    children = [fieldoffset(b, -1, "", T, nothing, nothing)]
    Meta.listStart(b)
    return Meta.List, Meta.listEnd(b), children
end

function juliaeltype(f::Meta.Field, list::Meta.FixedSizeList)
    type, _ = juliaeltype(f.children[1])
    return NTuple{Int(list.listSize), type}
end

function arrowtype(b, ::Type{NTuple{N, T}}) where {N, T}
    children = [fieldoffset(b, -1, "", T, nothing, nothing)]
    Meta.fixedSizeListStart(b)
    Meta.fixedSizeListAddListSize(b, Int32(N))
    return Meta.FixedSizeList, Meta.fixedSizeListEnd(b), children
end

function juliaeltype(f::Meta.Field, map::Meta.Map)
    K, _ = juliaeltype(f.children[1].children[1])
    V, _ = juliaeltype(f.children[1].children[2])
    return Pair{K, V}
end

function arrowtype(b, ::Type{Pair{K, V}}) where {K, V}
    children = [fieldoffset(b, -1, "entries", KeyValue{K, V}, nothing, nothing)]
    Meta.mapStart(b)
    return Meta.Map, Meta.mapEnd(b), children
end

default(::Type{Pair{K, V}}) where {K, V} = default(K) => default(V)

struct KeyValue{K, V}
    key::K
    value::V
end
Base.length(kv::KeyValue) = 1
Base.iterate(kv::KeyValue, st=1) = st === nothing ? nothing : (kv, nothing)
default(::Type{KeyValue{K, V}}) where {K, V} = KeyValue(default(K), default(V))

function arrowtype(b, ::Type{KeyValue{K, V}}) where {K, V}
    children = [fieldoffset(b, -1, "key", K, nothing, nothing), fieldoffset(b, -1, "value", V, nothing, nothing)]
    Meta.structStart(b)
    return Meta.Struct, Meta.structEnd(b), children
end

function juliaeltype(f::Meta.Field, list::Meta.Struct)
    names = Tuple(Symbol(x.name) for x in f.children)
    types = Tuple(juliaeltype(x)[1] for x in f.children)
    return NamedTuple{names, Tuple{types...}}
end

function arrowtype(b, ::Type{NamedTuple{names, types}}) where {names, types}
    children = [fieldoffset(b, -1, names[i], fieldtype(types, i), nothing, nothing) for i = 1:length(names)]
    Meta.structStart(b)
    return Meta.Struct, Meta.structEnd(b), children
end

default(::Type{NamedTuple{names, types}}) where {names, types} = NamedTuple{names}(Tuple(default(fieldtype(types, i)) for i = 1:length(names)))

# Unions
function juliaeltype(f::Meta.Field, u::Meta.Union)
    return UnionT{u.mode, u.typeIds !== nothing ? Tuple(u.typeIds) : u.typeIds, Tuple{(juliaeltype(x)[1] for x in f.children)...}}
end

# Note: nested Union types can't be represented using julia's builtin Union{...}
arrowtype(b, U::Union) = arrowtype(b, UnionT{Meta.UnionMode.Dense, nothing, Tuple{eachunion(U)...}})

function arrowtype(b, ::Type{UnionT{T, typeIds, U}}) where {T, typeIds, U}
    if typeIds !== nothing
        Meta.unionStartTypeIdsVector(b, length(typeIds))
        for id in Iterators.reverse(typeIds)
            FlatBuffers.prepend!(b, id)
        end
        TI = FlatBuffers.endvector!(b, length(typeIds))
    end
    children = [fieldoffset(b, -1, "", fieldtype(U, i), nothing, nothing) for i = 1:fieldcount(U)]
    Meta.unionStart(b)
    Meta.unionAddMode(b, T)
    if typeIds !== nothing
        Meta.unionAddTypeIds(b, TI)
    end
    return Meta.Union, Meta.unionEnd(b), children
end


text:
function update_weights(data, e, closest_same, closest_other, weights, weights_mult, m, k, max_f_vals, min_f_vals)

	for t = 1:size(data, 2)

		# Penalty term
		penalty = sum(abs.(e[t] .- closest_same[:, t])/((max_f_vals[t] .- min_f_vals[t]) .+ eps(Float64)))

		# Reward term
		reward = sum(weights_mult .* (abs.(e[t] .- closest_other[:, t])/((max_f_vals[t] .- min_f_vals[t] .+ eps(Float64)))))

		# Weights update
		weights[t] = weights[t] - penalty/(m*k) + reward/(m*k)
	end

	# Return updated weights.
	return weights
end

ext: jl

max_stars_repo_name: Architecton/diplomska-naloga-koda

max_issues_repo_name: Architecton/diplomska-naloga-koda

max_forks_repo_name: Architecton/diplomska-naloga-koda


text:
### Display

"""
    splitprint([io::IO], box)

Print a representation of all boxes below `box`, using parentheses prefaced by a number `n`
to denote a split along dimension `n`, and `l` to represent a leaf.
"""
function splitprint(io::IO, box::Box)
    if isleaf(box)
        print(io, 'l')
    else
        print(io, box.split.dims, '[')
        splitprint(io, box.split.self)
        for c in box.split.others.children
            print(io, ", ")
            splitprint(io, c)
        end
        print(io, ']')
    end
end
splitprint(box::Box) = splitprint(stdout, box)

"""
    splitprint_colored([io::IO], box, innerbox)

Like [`splitprint`](@ref), except that `innerbox` is highlighted in red, and the chain
of parents of `innerbox` are highlighted in cyan.
"""
function splitprint_colored(io::IO, box::Box, thisbox::Box, allparents=get_allparents(thisbox))
    if isleaf(box)
        box == thisbox ? printstyled(io, 'l', color=:light_red) : print(io, 'l')
    else
        if box == thisbox
            printstyled(io, box.split.dims, color=:light_red)
        elseif box ∈ allparents
            printstyled(io, box.split.dims, color=:cyan)
        else
            print(io, box.split.dims)
        end
        print(io, '[')
        splitprint_colored(io, box.split.self, thisbox, allparents)
        for c in box.split.others.children
            print(io, ", ")
            splitprint_colored(io, c, thisbox, allparents)
        end
        print(io, ']')
    end
end
splitprint_colored(box::Box, thisbox::Box) = splitprint_colored(stdout, box, thisbox)

function get_allparents(box)
    allparents = Set{typeof(box)}()
    p = box
    while !isroot(p)
        p = p.parent
        push!(allparents, p)
    end
    allparents
end

print_metabox(io::IO, box::Box) = print(io, box, " with metadata ", box.metabox)

### Geometry

"""
    x = position(box)

Return the n-dimensional position vector `x` at which this box was evaluated.
"""
Base.position(box::Box{p,T}) where {p,T} = position!(Vector{T}(undef, ndims(box)), box)

"""
    x = position(box, d::Integer)

Return the position of `box` in dimension `d`.
"""
function Base.position(box::Box, d::Integer)
    default = baseposition(box.world.position[d])
    while !isroot(box)
        p = box.parent
        childindex = box.childindex
        if !isself(box)
            split = p.split
            i = findfirst(isequal(d), split.dims)
            if i !== nothing && ((childindex >> (i-1))&0x01) != 0x00
                return split.xs[i]
            end
        end
        box = p
    end
    return default
end

"""
    position!(x::AbstractVector, box)
    position!(x::AbstractVector, filled::AbstractVector{Bool}, box)

Fill `x` with the n-dimensional position of `box`. Including `filled` avoids allocation of temporaries.
"""
position!(x, box::Box) = position!(x, Vector{Bool}(undef, ndims(box)), box)

function position!(x, filled, box::Box)
    N = ndims(box)
    @assert(LinearIndices(x) == LinearIndices(filled) == Base.OneTo(N))
    # Fill in default values (those corresponding to root)
    for i = 1:N
        x[i] = baseposition(box.world.position[i])
    end
    # Traverse the tree to see which entries have changed since root
    fill!(filled, false)
    nfilled = 0
    while !isroot(box) && nfilled < length(filled)
        p = box.parent
        split = p.split
        childindex = box.childindex
        for i = 1:degree(box)
            sd, xo = split.dims[i], split.xs[i]
            sd > ndims(box) && continue  # fictive dimension
            if ((childindex >> (i-1))&0x01) != 0x00 && !filled[sd]
                x[sd] = xo
                filled[sd] = true
                nfilled += 1
            end
        end
        box = p
    end
    return x
end

"""
    bbs = boxbounds(box)

Compute the bounds (edge positions) of `box`. `bbs[d] = (lower, upper)` for dimension `d`.
"""
function boxbounds(box::Box{p,T}) where {p,T}
    bbs = Vector{Tuple{T,T}}(undef, ndims(box))
    return boxbounds!(bbs, box)
end

"""
    lower, upper = boxbounds(box, d::Integer)

Compute the bounds of `box` along dimension `d`.
"""
function boxbounds(box::Box, d::Integer)
    lower, upper = box.world.lower[d], box.world.upper[d]
    box0 = box
    lfilled, ufilled = false, false
    while !isroot(box) && !(lfilled & ufilled)
        p = box.parent
        split = p.split
        i = findfirst(isequal(d), split.dims)
        if i !== nothing
            xs, xo = position(split.self, d), split.xs[i]
            x = position(box0, d)
            # @show p lower upper lfilled ufilled xs xo x
            xl, xu = xs < xo ? (xs, xo) : (xo, xs)
            # @show xl x xu
            if xl <= x <= xu
                if !lfilled && xl < x
                    lower, lfilled = max((xl+xu)/2, (xl+x)/2), true
                elseif !ufilled && xu > x
                    upper, ufilled = min((xl+xu)/2, (x+xu)/2), true
                end
            end
        end
        box = p
    end
    lower, upper
end

"""
    boxbounds!(bbs::AbstractVector, box)

Fill `bbs` with the bounds (edge positions) of `box`.
"""
function boxbounds!(bbs, box::Box)
    lfilled = Vector{Bool}(undef, ndims(box))
    ufilled = Vector{Bool}(undef, ndims(box))
    boxbounds!(bbs, lfilled, ufilled, box)
end

function boxbounds!(bbs, lfilled, ufilled, box::Box)
    fill!(lfilled, false)
    fill!(ufilled, false)
    n = ndims(box)
    for d = 1:n
        bbs[d] = (box.world.lower[d], box.world.upper[d])
    end
    box0 = box
    nlfilled = nufilled = 0
    while !isroot(box) && min(nlfilled, nufilled) < n
        p = box.parent
        split = p.split
        for (sd, xo) = zip(split.dims, split.xs)
            sd > n && continue   # fictive dimensions
            if !lfilled[sd] || !ufilled[sd]
                xs = position(split.self, sd)
                x = position(box0, sd)
                xl, xu = xs < xo ? (xs, xo) : (xo, xs)
                if xl <= x <= xu
                    if !lfilled[sd] && xl < x
                        bbs[sd], lfilled[sd] = (max((xl+xu)/2, (xl+x)/2), bbs[sd][2]), true
                        nlfilled += 1
                    elseif !ufilled[sd] && xu > x
                        bbs[sd], ufilled[sd] = (bbs[sd][1], min((xl+xu)/2, (x+xu)/2)), true
                        nufilled += 1
                    end
                end
            end
        end
        box = p
    end
    bbs
end

function Base.in(xs::AbstractVector, box::Box)
    for (bb, x) in zip(boxbounds(box), xs)
        bb[1] <= x <= bb[2] || return false
    end
    return true
end

"""
    nc = ncollinear(box, top=getroot(box))

`nc[i]` is the number of points contained inside `top` that are
collinear with `position(box)` along dimension `i`. Being collinear
along dimension `i` implies that `position(otherbox, j) == position(box, j)`
for any dimension `j != i`.
"""
ncollinear(box::Box, top::Box=getroot(box)) = ncollinear!(Vector{Int}(undef, ndims(box)), Vector{Bool}(undef, ndims(box)), box, top)

function ncollinear!(nc::AbstractVector{Int}, filled::AbstractVector{Bool}, box::Box, top::Box=getroot(box))
    fill!(nc, 0)
    fill!(filled, false)
    nfilled = 0
    while nfilled < 2 && box != top
        split = box.parent.split
        childindex = cindex = box.childindex
        wasbase = nfilled == 0
        for i = 1:degree(box)
            sd = split.dims[i]
            otherindex = childindex ⊻ (0x01 << (i-1))
            if cindex&0x01 == 0x00 # if this is the same as the parent point
                (wasbase || filled[sd]) && _ncollinear!(nc, split.others.children[otherindex], sd)
            else  # if we go up with shift relative to the basepoint
                otherchild = otherindex == 0x00 ? split.self : split.others.children[otherindex]
                (wasbase || filled[sd]) && _ncollinear!(nc, otherchild, sd)
                nfilled += !filled[sd]
                filled[sd] = true
            end
            cindex = cindex >> 1
        end
        box = box.parent
    end
    return nc
end

function _ncollinear!(nc, box, splitdim)
    if isleaf(box)
        nc[splitdim] += 1
        return nc
    end
    split = box.split
    i = findfirst(isequal(splitdim), split.dims)
    _ncollinear!(nc, split.self, splitdim)
    i === nothing && return nc
    childidx = 1 << (i-1)
    _ncollinear!(nc, split.others.children[childidx], splitdim)
end

function epswidth(bb::Tuple{T,T}) where T<:AbstractFloat
    w1 = isfinite(bb[1]) ? eps(bb[1]) : T(0)
    w2 = isfinite(bb[2]) ? eps(bb[2]) : T(0)
    return max(w1, w2)
end
epswidth(bb::Tuple{Real,Real}) = epswidth(Float64.(bb))

# 10*eps is pretty arbitrary, but not unreasonable
splittable(bb::Tuple{Real,Real}) = bb[2] - bb[1] > 10*epswidth(bb)

### Traversal

"""
    root = getroot(box)

Return the root node for `box`.
"""
function getroot(box::Box)
    while !isroot(box)
        box = box.parent
    end
    return box
end

"""
    leaf = getleaf(box)

Get the leaf node at the position of `box`.
"""
function getleaf(box::Box)
    while !isleaf(box)
        box = box.split.self
    end
    return box
end

"""
    m = meta(box)

Return the metadata associated with `box`.
"""
function meta(box::Box)
    while isself(box)
        box = box.parent
    end
    isroot(box) && return box.world.meta
    @assert(box.childindex > 0)
    return box.parent.split.others.metas[box.childindex]
end

getchild(split::Split, childidx) = childidx == 0 ? split.self :
                                                   split.others.children[childidx]

# """
#     boxp = find_parent_with_splitdim(box, splitdim::Integer)

# Return the first node at or above `box` who's parent box was split
# along dimension `splitdim`. If `box` has not yet been split along
# `splitdim`, returns the root box.
# """
# function find_parent_with_splitdim(box::Box, splitdim::Integer)
#     while !isroot(box)
#         p = box.parent
#         if p.split.dim == splitdim
#             return box
#         end
#         box = p
#     end
#     return box
# end

"""
    box = find_leaf_at(root, x)

Return the leaf-node `box` beneath `root` that contains `x`.
"""
function find_leaf_at(root::Box, x)
    @noinline throwbb(bb, x, dim) = error("$x not within $bb along dimension $dim")
    @noinline throwdmm(n, l) = throw(DimensionMismatch("tree has $n dimensions, got $l"))
    n = ndims(root)
    length(x) == n || throwdmm(n, length(x))
    # Check that x is in the interior of root
    for i = 1:n
        bb = boxbounds(root, i)
        ((bb[1] <= x[i]) && (x[i] < bb[2] || (x[i] == bb[2] && bb[2] == root.world.upper[i]))) || throwbb(bb, x[i], i)
    end
    isleaf(root) && return root
    while !isleaf(root)
        split = root.split
        childidx = 0  # assign using bitwise arithmetic
        for j = 1:degree(root)
            sd, xother = split.dims[j], split.xs[j]
            sd > n && continue
            xself = position(root, sd)
            xsd = x[sd]
            # We need to choose the child where, for bb = boxbounds(child, sd),
            # we have bb[1] <= xsd < bb[2]. This is a requirement if we plan
            # to split child using the positions in x.
            #
            # Consequently this bit should be set according to
            #
            #                  |   xother > xself   |    xother < xself
            #    --------------|--------------------|---------------------
            #    xsd >= xmid   |         1          |           0
            #    --------------|--------------------|---------------------
            #    xsd < xmid    |         0          |           1
            #
            # where xmid = (xself+xother)/2.
            childidx |= ((xsd >= (xself+xother)/2) ⊻ (xother < xself)) << (j-1)
        end
        root = getchild(split, childidx)
    end
    return root
end

### Iteration

function Base.iterate(root::Box)
    isleaf(root) && return root, VisitorState(root, maxchildren(root))
    return root, VisitorState(root, 0)
end

function Base.iterate(root::Box, state::VisitorState)
    box, childindex = state.box, state.childindex
    if childindex < maxchildren(box) && !isleaf(box)
        # descend
        item = getchild(box.split, childindex)
        return item, VisitorState(item, isleaf(item) ? maxchildren(item) : 0)
    end
    # ascend
    while box != root
        childindex = box.childindex + 1
        box = box.parent
        if childindex < maxchildren(box)
            item = getchild(box.split, childindex)
            return item, VisitorState(item, isleaf(item) ? maxchildren(item) : 0)
        end
    end
    return nothing
end

"""
    iter = leaves(box)

Return an iterator visiting the leaf-nodes below `box` in depth-first order.
"""
leaves(root::Box) = Iterators.filter(box->isleaf(box) && !isfake(box), root)

"""
    isfake(box)

True if `box` has been displaced from its parent along a fictive dimension.
"""
function isfake(box::Box)
    n = ndims(box)
    p = box.parent
    isleaf(p) && return false
    split = p.split
    tf = false
    cidx = box.childindex
    for d in split.dims
        tf |= ((cidx & 0x01) == 0x01) & (d > n)
        cidx = cidx >> 0x01
    end
    return tf
end

"""
    iter = splits(box)

Return an iterator for visiting the dimension-splits in the sub-tree
containing `box`. Branches below `box` are explored first; once
exhausted, a step is taken up the tree and branches below the other
children of `box`'s parent are explored. Search terminates when all
branches below the root node have been explored.

The iterator returns the `split` field of non-leaf boxes (see
[`CoordinateSplittingPTrees.Split`](@ref)), so a loop should be written as

    for split in splits(box)
        # do something with split
    end
"""
splits(box::Box) = SplitIterator(box)

function nobranch(box, skipchildindex, childindex)
    return childindex < maxchildren(box) &&
        (childindex == skipchildindex ||
         isleaf(getchild(box.split, childindex)))
end

function Base.iterate(splits::SplitIterator)
    box = splits.base
    if !isleaf(box)
        # Start at the current split if that's what was used to create the iterator
        branchhead = box.split.self
        branchiter = Iterators.filter(isnonleaf, branchhead)
        return box.split, ClimbingState(box, true, typemax(UInt8), 0, branchiter, VisitorState(branchhead, -1))
    end
    # The iterator was created from a leaf, which has no splits.
    # Go up and grab a sibling or parent.
    isroot(box) && return nothing
    skipchildindex, childindex = box.childindex, 0
    box = box.parent
    while nobranch(box, skipchildindex, childindex)
        childindex += 1
    end
    if childindex >= maxchildren(box)
        # This is a tree with only 1 split, so use the parent
        branchiter = Iterators.filter(isnonleaf, box)
        item, branchstate = iterate(branchiter)
        return item.split, ClimbingState(box, true, skipchildindex, childindex, branchiter, branchstate)
    end
    # Return a sibling
    branchhead = getchild(box.split, childindex)
    branchiter = Iterators.filter(isnonleaf, branchhead)
    item, branchstate = iterate(branchiter)  # this can't fail to return an item, because we already checked
    return item.split, ClimbingState(box, false, skipchildindex, childindex, branchiter, branchstate)
end

Base.iterate(splits::SplitIterator, ::Nothing) = iterate(splits)

function Base.iterate(splits::SplitIterator, state::ClimbingState)
    # First, try iterating on the existing branch and see if that returns a valid item
    branchiter, branchstate = state.branchiter, state.branchstate
    ret = branchstate.childindex < 0 ? iterate(branchiter) : iterate(branchiter, branchstate)
    if ret !== nothing
        box, branchstate = ret
        return box.split, ClimbingState(state, branchstate)
    end
    # The branch is done. Let's see if we can find another one.
    box = state.box
    skipchildindex, childindex = state.skipchildindex, state.childindex+1
    while nobranch(box, skipchildindex, childindex)
        childindex += 1
    end
    if childindex < maxchildren(box)
        # We still haven't exhausted the children of the current box
        branchhead = getchild(box.split, childindex)
        branchiter = Iterators.filter(isnonleaf, branchhead)
        item, branchstate = iterate(branchiter)  # this can't fail to return an item, because we already checked
        return item.split, ClimbingState(box, state.visited, skipchildindex, childindex, branchiter, branchstate)
    end
    if !state.visited
        # We haven't returned the current box's split, so do that now
        return box.split, ClimbingState(box, true, state.skipchildindex, childindex, branchiter, branchstate)
    end
    # Go up
    isroot(box) && return nothing
    skipchildindex = box.childindex
    box = box.parent
    childindex = 0
    while nobranch(box, skipchildindex, childindex)
        childindex += 1
    end
    if childindex >= maxchildren(box)
        return box.split, ClimbingState(box, true, state.skipchildindex, childindex, branchiter, branchstate)
    end
    branchhead = getchild(box.split, childindex)
    branchiter = Iterators.filter(isnonleaf, branchhead)
    item, branchstate = iterate(branchiter)  # this can't fail to return an item, because we already checked
    return item.split, ClimbingState(box, false, skipchildindex, childindex, branchiter, branchstate)
end

"""
    iter = chain(box)

Create an iterator that visits splits that minimally cover all dimensions
"""
chain(box::Box) = ChainIterator(box)

function Base.iterate(iter::ChainIterator)
    ndims(iter.base) == 0 && return nothing
    box, count = iter.base, 0
    item = box.split
    return item, (item.others.children[end], count+degree(box))
end

function Base.iterate(iter::ChainIterator, state)
    state[2] >= ndims(iter.base) && return nothing
    box, count = state
    item = box.split
    return item, (item.others.children[end], count+degree(box))
end

Base.length(iter::Union{Box,CSpTreeIterator}) = _length(iter)
Base.length(iter::Iterators.Filter{F,I}) where {F,I<:Union{Box,CSpTreeIterator}} =
    _length(iter)
Base.length(iter::Iterators.Filter{F,I}) where {F,I<:Iterators.Filter{F2,I2}} where {F2,I2<:Box} =
    _length(iter)

function _length(iter)
    len = 0
    ret = iterate(iter)
    while ret !== nothing
        len += 1
        _, state = ret
        ret = iterate(iter, state)
    end
    len
end

function skipsplits(iter, skip)
    skip == 0 && return nothing
    _, state = iterate(iter)
    for i = 1:skip-1
        _, state = iterate(iter, state)
    end
    return state
end

# useful for debugging
function Base.show(io::IO, state::ClimbingState)
    println(io, "box $(state.box), visited $(state.visited), splitdims $(state.box.split.dims), skipidx $(state.skipchildindex), childidx $(state.childindex)")
    print(io, "branchiter $(state.branchiter), branch-is-done $(iterate(state.branchiter, state.branchstate)===nothing)")
end


text:
using DifferentialEquations
using Roots
using PyCall
include("YfindBIG.jl")
pygui(:qt)
using PyPlot
pygui(true)

Energy=BigFloat("0.1249")
H=(2*Energy)^2
t_end=100

barrier=20
Hamil(XX,YY,QQ,PP)=( (QQ-XX)^2+(PP-YY)^2 )*( (QQ+XX)^2+(PP+YY)^2 )/((PP^4+2*PP^2*(XX^2-1)+(1+XX^2)^2 )*(QQ^4+2*QQ^2*(YY^2+1)+(YY^2-1)^2 ))
H_test(u)=( (u[3]-u[1])^2+(u[2]-u[4])^2 )*( (u[3]+u[1])^2+(u[2]+u[4])^2 )/ ((u[2]^4+2*u[2]^2*(u[1]^2-1)+(1+u[1]^2)^2 )*(u[3]^4+2*u[3]^2*(u[4]^2+1)+(u[4]^2-1)^2 ))
function g(resid,u,p,t)
  resid[1] = ( (u[3]-u[1])^2+(u[2]-u[4])^2 )*( (u[3]+u[1])^2+(u[2]+u[4])^2 )/ ((u[2]^4+2*u[2]^2*(u[1]^2-1)+(1+u[1]^2)^2 )*(u[3]^4+2*u[3]^2*(u[4]^2+1)+(u[4]^2-1)^2 ))-H
  resid[2] = 0
  resid[3] = 0
  resid[4] = 0
end
condition(u,t,integrator)=abs(u[1])>barrier || abs(u[2])>barrier || abs(u[3])>barrier||abs(u[4])>barrier
affect!(integrator) = terminate!(integrator)
cb1 = DiscreteCallback(condition,affect!)

# Restrict to energy surface
cb2 = ManifoldProjection(g)
cb=CallbackSet(cb1)

ODE1(z,w)=conj(  im * w.*( BigFloat("1") ./(z.^2-w.^2)+BigFloat("1") ./(1+w.^2) ))
ODE2(z,w)=conj(  im * z.*( BigFloat("1") ./(w.^2-z.^2)+BigFloat("1") ./(1+z.^2) ))
function Eq_of_M(du,u,p,t)
    du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    return
end
# refernce
# z=X+iP W=Q+iY . Plane Q=0, P=0
# z=u[1]+i u[2] ; w=u[3]+i u[4]

#Q=0, have P find Y

Q=BigFloat(".1")
P=BigFloat(".1")
Y=YfindBIG(Q,P,H)
u0=zeros(4)
u0[1]=0 #X
u0[2]=P #P
u0[3]=Q #Q
u0[4]=Y #Y
TEST1_0=H_test(u0)/H
prob = ODEProblem(Eq_of_M,u0,(0., t_end))
# t,A=solve(prob, RK4(),reltol=1e-15,abstol=1e-16,maxiters=1e15)

sol=solve(prob, Vern9(),maxiters=BigInt(1e10), reltol=BigFloat("1e-15"),abstol=BigFloat("1e-17"),callback=cb)
# figure()

A=sol[:,:]
A=convert(Array{Float64},A)
figure()
 plot3D(A[3,:],A[2,:],A[1,:],linewidth=.5)

N=length(A[1,:])
H_v=zeros(BigFloat,N)
umax=zeros(BigFloat ,length(A[1,:]))
dH=zeros(BigFloat ,length(A[1,:]))
for i=1:length(A[1,:])
    uf=zeros(BigFloat ,4)
    uf[1]=A[1,i] #X
    uf[2]=A[2,i] #P
    uf[3]=A[3,i]#Q
    uf[4]=A[4,i] #Y
    dH[i]=abs(H_test(uf)/H)
    umax[i]=maximum(uf)
end
MAX=maximum(umax)
worse=log(maximum(dH))
dH=convert(Array{Float64}, dH)

figure()
plot(1:N, log.(dH))


text:
apply_palette(p::AbstractVector, uv) = [cycle(p, idx) for idx in eachindex(uv)]
apply_palette(::Automatic, uv) = eachindex(uv)
apply_palette(p, uv) = map(p, eachindex(uv))

# TODO: add more customizations?
struct Wrap end

const wrap = Wrap()

function apply_palette(::Wrap, uv)
    ncols = ceil(Int, sqrt(length(uv)))
    return [fldmod1(idx, ncols) for idx in eachindex(uv)]
end

struct ContinuousScale{T, F}
    f::F
    extrema::Tuple{T, T}
end

rescale(values::AbstractArray{<:Number}, c::ContinuousScale) = values # Is this ideal?
function rescale(values::AbstractArray{<:Union{Date, DateTime}}, c::ContinuousScale)
    @assert c.f === identity
    min, max = c.extrema
    return @. convert(Millisecond, DateTime(values) - DateTime(min)) / Millisecond(1)
end

struct CategoricalScale{S, T}
    data::S
    plot::T
end

function rescale(values, c::CategoricalScale)
    idxs = indexin(values, c.data)
    return c.plot[idxs]
end

Base.length(c::CategoricalScale) = length(c.data)

rescale(values, ::Nothing) = values

default_scale(::Nothing, palette) = nothing

function default_scale(summary::Tuple, palette)
    f = palette isa Function ? palette : identity
    return ContinuousScale(f, summary)
end

function default_scale(summary::AbstractVector, palette)
    plot = apply_palette(palette, summary)
    return CategoricalScale(summary, plot)
end

# Logic to infer good scales
function default_scales(summaries, palettes)
    defaults = Dict{KeyType, Any}()
    for (key, val) in pairs(summaries)
        # ensure `palette === automatic` for integer keys
        palette = key in propertynames(palettes) ? palettes[key] : automatic
        defaults[key] = default_scale(val, palette)
    end
    return defaults
end

# Logic to create ticks from a scale
# Should take current tick to incorporate information
function ticks(scale::CategoricalScale)
    u = map(string, scale.data)
    return (axes(u, 1), u)
end

function ticks(scale::ContinuousScale)
    return continuousticks(scale.extrema...)
end

continuousticks(min, max) = automatic

function continuousticks(min::T, max::T) where T<:Union{Date, DateTime}
    min_ms::Millisecond, max_ms::Millisecond = DateTime(min), DateTime(max)
    min_pure, max_pure = min_ms/Millisecond(1), max_ms/Millisecond(1)
    dates, labels = optimize_datetime_ticks(min_pure, max_pure)
    return (dates .-  min_pure, labels)
end

## Scale helpers

const ArrayLike = Union{AbstractArray, Tuple}
const StringLike = Union{AbstractString, Symbol}

function cycle(v::AbstractVector, i::Int)
    ax = axes(v, 1)
    return v[first(ax) + mod(i - first(ax), length(ax))]
end

"""
    iscontinuous(v::AbstractArray)

Determine whether `v` should be treated as a continuous or categorical vector.
"""
iscontinuous(::AbstractArray) = false
iscontinuous(::AbstractArray{<:Number}) = true
iscontinuous(::AbstractArray{<:Union{Date, DateTime}}) = true

isgeometry(::AbstractArray{<:AbstractGeometry}) = true
isgeometry(::AbstractArray{T}) where {T} = eltype(T) <: AbstractGeometry

extend_extrema((l1, u1), (l2, u2)) = min(l1, l2), max(u1, u2)

push_different!(v, val) = !isempty(v) && isequal(last(v), val) || push!(v, val)

function mergesorted(v1, v2)
    issorted(v1) && issorted(v2) || throw(ArgumentError("arguments must be sorted"))
    T = promote_type(eltype(v1), eltype(v2))
    v = sizehint!(T[], length(v1) + length(v2))
    i1, i2 = 1, 1
    while i2 ≤ length(v2)
        while i1 ≤ length(v1) && isless(v1[i1], v2[i2])
            push_different!(v, v1[i1])
            i1 += 1
        end
        push_different!(v, v2[i2])
        i2 += 1
    end
    while i1 ≤ length(v1)
        push_different!(v, v1[i1])
        i1 += 1
    end
    return v
end


text:
using StokesDiffEq

## Change plot commands to GR

Δx = 1//2^3 #Super large for quick testing. Lower this for better results.
mesh = FDMMesh(Δx,mins=[-1;-1],maxs=[1;1])
prob = prob_stokes_homogenous

#DGS Convergence
@time sol = solve(prob,mesh,converrors=true,maxiters=200,alg=:DGS)

err1 = sol.converrors[:rul∞]
err2 = sol.converrors[:rvl∞]
err3 = sol.converrors[:rpl∞]

#TEST_PLOT && plot([1:length(err1) 1:length(err2) 1:length(err3)],[err1 err2 err3],yscale=:log10,xguide="Error",yguide="Iterations",title="DGS Convergence, Delta t = 2^{-3}",label=["u" "v" "p"])

#Multigrid Relative Convergence
@time sol = solve(prob,mesh,converrors=true,maxiters=20,alg=:Multigrid,levels=4)

err1 = sol.converrors[:rresul2]
err2 = sol.converrors[:rresvl2]
err3 = sol.converrors[:rrespl2]
#TEST_PLOT && plot([1:length(err1) 1:length(err2) 1:length(err3)],[err1 err2 err3],yscale=:log10,xguide="Error",yguide="Iterations",title="Multigrid Relative Residual Convergence, Delta t = 2^{-3}",label=["u" "v" "p"])

err1[end] < 1e-12


text:
export pole_info

"""
### function pole_info_header()

Print the header of the table that shows the information about the poles:

                      Pole |         Damping |   Freq. [rad/s] | Time Const. [s]
    -----------------------+-----------------+-----------------+----------------

"""

function pole_info_header()
    println(((@bold)) *
            lpad("Pole",            22, " ") * " |" *
            lpad("Damping",         16, " ") * " |" *
            lpad("Freq. [rad/s]",   16, " ") * " |" *
            lpad("Time Const. [s]", 16, " ") *
            ((@normal)))

    println(((@bold)) * lpad("", 22, "-") * "-+" *
            lpad("", 16, "-")           * "-+" *
            lpad("", 16, "-")           * "-+" *
            lpad("", 16, "-")           *
            ((@normal)))
end

"""
### function pole_info(pole::Complex128; header = true, highlight = false)

Print the pole information:

* Pole: The location of the pole in the complex plane.
* Damping: The damping of the pole.
* Frequency: The natural frequency of the pole [rad/s].
* Time constant: The time constant of the pole [s].

##### Args

* pole: The pole.

##### Parameters

* header: If `true`, then the header will be printed (see `pole_info_header()`).
* highlight: If `true`, then the text will be printed in yellow.

"""

function pole_info(pole::Complex128; header = true, highlight = false)
    # Compute the damping, frequency, and time constant of the pole.
    ω_n = norm(pole)
    ζ = -real(pole)/ω_n
    t = 1/(ω_n*ζ)

    # Convert the numbers into text.
    pole_str = lpad("$(round(pole,5))",      22, " ")
    ω_n_str  = lpad((@sprintf "%16.5g" ω_n), 16, " ")
    ζ_str    = lpad((@sprintf "%16.5g" ζ),   16, " ")
    t_str    = lpad((@sprintf "%16.5g" t),   16, " ")

    (header) && pole_info_header()

    println(( highlight ? (@yellow) : "") * pole_str * ((@boldn)) * " |" * ((@normal)) *
            ( highlight ? (@yellow) : "") * ω_n_str  * ((@boldn)) * " |" * ((@normal)) *
            ( highlight ? (@yellow) : "") * ζ_str    * ((@boldn)) * " |" * ((@normal)) *
            ( highlight ? (@yellow) : "") * t_str    * ((@treset)))

    nothing
end

"""
### function pole_info(poles::Array{Complex128,1}; highlight_ind = 0)

Print the information of an array of poles.

##### Args

* poles: The array of poles.

##### Parameters

* highlight_ind: The index of a pole to be highlighted.

"""

function pole_info(poles::Array{Complex128,1}; highlight_ind = 0)
    pole_info_header()

    for i = 1:length(poles)
        if (i == highlight_ind)
            pole_info(poles[i]; header = false, highlight = true)
        else
            pole_info(poles[i]; header = false)
        end
    end

    nothing
end


text:
module TestGBIFRecords

using GBIF
using Test

# Version using pairs
set1 = occurrences("scientificName" => "Mus musculus", "year" => 1999, "hasCoordinate" => true)
@test typeof(set1) == GBIFRecords
@test length(set1) == 20

# Version with no query parameters
set2 = occurrences()
@test typeof(set2) == GBIFRecords
@test length(set2) == 20

# Version using ranged pairs
set3 = occurrences("scientificName" => "Mus musculus", "year" => 1999, "hasCoordinate" => true, "decimalLatitude" => (0.0, 50.0))
@test typeof(set3) == GBIFRecords
@test length(set3) == 20

# Version with the full query - this one has about 250 records
serval = GBIF.taxon("Leptailurus serval", strict=true)
obs = occurrences(serval, "hasCoordinate" => "true", "continent" => "AFRICA", "decimalLongitude" => (-30, 40))
while length(obs) < size(obs)
    occurrences!(obs)
end
@test length(obs) == size(obs)

# Version with the full query AND a set page size - this one has about 250 records
obs = occurrences(serval, "hasCoordinate" => "true", "continent" => "AFRICA", "decimalLongitude" => (-30, 40), "limit" => 45)
while length(obs) < size(obs)
    occurrences!(obs)
end
@test length(obs) == size(obs)

end


text:
function plotTrajectory(u_hat::Vector{Float64}, o_hat::Vector{Float64})
    plot(1:N, o_hat, color="black", marker="o", linewidth=2, alpha=0.8) # Observed locations
    grid("on")
    xlabel("Time t")
    ylabel("Position o")
end

"""
To find the index of the slide message in the schedule:
schedule = sumProductSchedule(u[2])
slide_msg_idx = findMessageIndex(schedule, :equ_s_1_1, 3)
"""
function findMessageIndex(schedule::Schedule, node_id::Symbol, outbound_interface_id::Union{Int64, Symbol})
    condensed_schedule = ForneyLab.condense(ForneyLab.flatten(schedule))
    for (idx, entry) in enumerate(condensed_schedule)
        node = entry.interface.node
        if (node.id == node_id) && (entry.interface == node.i[outbound_interface_id])
            return idx
        end
    end

    error("Message $node_id[$outbound_interface_id] could not be found in schedule")
end

function inspectSnippet(algo::String)
    println(algo[8:620])
    println("\n...\n")
    println(algo[end-89:end-13])
end

max_stars_repo_name: Yikeru/BMLIP

max_issues_repo_name: Yikeru/BMLIP

max_forks_repo_name: Yikeru/BMLIP


text:
const SEED = 999
const NROWS = 128
const N_IDX_TESTS = 16

Random.seed!(SEED)

randdate() = Date(rand(0:4000), rand(1:12), rand(1:27))
randtime() = Dates.Time(rand(0:23), rand(0:59), rand(0:59))
randdatetime() = randdate() + randtime()

randstrings() = String[[randstring(rand(0:20)) for i ∈ 1:(NROWS-1)]; "a"]
function randstrings(::Missing)
    Union{String,Missing}[[rand(Bool) ? missing : randstring(rand(0:20)) for i ∈ 1:(NROWS-1)]; "a"]
end

convstring(str::AbstractString) = String(str)
convstring(::Missing) = missing

io = IOBuffer()

@testset "ArrowTests" begin
df = DataFrame(ints=rand(Int32,NROWS),
               floats=rand(Float64,NROWS),
               dates=Date[randdate() for i ∈ 1:NROWS],
               datetimes=DateTime[randdatetime() for i ∈ 1:NROWS],
               times=Dates.Time[randtime() for i ∈ 1:NROWS],
               missingints=Union{Int64,Missing}[rand(Bool) ? missing : rand(Int64) for i ∈ 1:NROWS],
               strings=randstrings(),
               missingstrings=randstrings(missing),
               catstrings=categorical(randstrings()),
               catstringsmissing=categorical(randstrings(missing))
              )

Feather.write(io, df)

ndf = Feather.read(io)

@test typeof(ndf[:ints]) == Primitive{Int32}
@test typeof(ndf[:floats]) == Primitive{Float64}
@test typeof(ndf[:dates]) == Primitive{Arrow.Datestamp}
@test typeof(ndf[:datetimes]) == Primitive{Arrow.Timestamp{Dates.Millisecond}}
@test typeof(ndf[:times]) == Primitive{Arrow.TimeOfDay{Dates.Nanosecond,Int64}}
@test typeof(ndf[:missingints]) == NullablePrimitive{Int64}
@test typeof(ndf[:strings]) == List{String,Arrow.DefaultOffset,Primitive{UInt8}}
@test typeof(ndf[:missingstrings]) == NullableList{String,Arrow.DefaultOffset,Primitive{UInt8}}
@test typeof(ndf[:catstrings]) == DictEncoding{String,Primitive{Int32},
                                               List{String,Arrow.DefaultOffset,Primitive{UInt8}}}
@test typeof(ndf[:catstringsmissing]) ==
        DictEncoding{Union{String,Missing},NullablePrimitive{Int32},List{String,Arrow.DefaultOffset,
                                                                         Primitive{UInt8}}}

for j ∈ 1:N_IDX_TESTS
    i = rand(1:NROWS)
    @test df[i, :ints] == ndf[i, :ints]
    @test df[i, :floats] == ndf[i, :floats]
    @test df[i, :dates] == convert(Date, ndf[i, :dates])
    @test df[i, :datetimes] == convert(DateTime, ndf[i, :datetimes])
    @test df[i, :times] == convert(Dates.Time, ndf[i, :times])
    @test df[i, :missingints] ≅ ndf[i, :missingints]
    @test df[i, :strings] == ndf[i, :strings]
    @test df[i, :missingstrings] ≅ ndf[i, :missingstrings]
    @test df[i, :catstrings] == String(ndf[i, :catstrings])
    @test df[i, :catstringsmissing] ≅ convstring(ndf[i, :catstringsmissing])
end
for j ∈ 1:N_IDX_TESTS
    a, b = extrema(rand(1:NROWS, 2))
    i = a:b
    @test df[i, :ints] == ndf[i, :ints]
    @test df[i, :floats] == ndf[i, :floats]
    @test df[i, :dates] == convert.(Date, ndf[i, :dates])
    @test df[i, :datetimes] == convert.(DateTime, ndf[i, :datetimes])
    @test df[i, :times] == convert.(Dates.Time, ndf[i, :times])
    @test df[i, :missingints] ≅ ndf[i, :missingints]
    @test df[i, :strings] == ndf[i, :strings]
    @test df[i, :missingstrings] ≅ ndf[i, :missingstrings]
    @test df[i, :catstrings] == String.(ndf[i, :catstrings])
    @test df[i, :catstringsmissing] ≅ convstring.(ndf[i, :catstringsmissing])
end
for j ∈ 1:N_IDX_TESTS
    i = rand(1:NROWS, rand(1:4))
    @test df[i, :ints] == ndf[i, :ints]
    @test df[i, :floats] == ndf[i, :floats]
    @test df[i, :dates] == convert.(Date, ndf[i, :dates])
    @test df[i, :datetimes] == convert.(DateTime, ndf[i, :datetimes])
    @test df[i, :times] == convert.(Dates.Time, ndf[i, :times])
    @test df[i, :missingints] ≅ ndf[i, :missingints]
    @test df[i, :strings] == ndf[i, :strings]
    @test df[i, :missingstrings] ≅ ndf[i, :missingstrings]
    @test df[i, :catstrings] == String.(ndf[i, :catstrings])
    @test df[i, :catstringsmissing] ≅ convstring.(ndf[i, :catstringsmissing])
end
end

ndf = nothing;
GC.gc()


text:
using BlackBoxOptim

# Another problem from Hans W. Borchers, constrained 15-points:
function points15(x)
    d = 2.0
    for i = 1:14, j = (i+1):15
        s = (x[i]-x[j])^2 + (x[i+15]-x[j+15])^2 + (x[i+30]-x[j+30])^2
        s = sqrt(s)
        if s < d
          d = s
        end
    end
    return -d
end

Dims = 15 * 3

# The best known minimum is around -0.62 but we do not clearly know if there are better
# ones so we aim somewhat lower...
p = BlackBoxOptim.fixeddim_problem(points15; range = (0.0, 1.0),
  dims = Dims, name = "constrained points 15", fmins = [-0.64])

# This is a hard problem so run for some time...
MaxMinutes = 0.50

# Run the AdaptiveDE once:
#bboptimize(points15, (0.0, 1.0); dimensions = 45, max_time = MaxMinutes * 60)

# Run several methods several times:
BlackBoxOptim.repeated_bboptimize(2, p, Dims, [
  :generating_set_search,
  :adaptive_de_rand_1_bin_radiuslimited,
  #:random_search,
  #:separable_nes # Can't apply separable NES since it does not yet respect box constraints...
  ],
  MaxMinutes * 60)

#compare_optimizers(points15; search_range = (0.0, 1.0), dimensions = 45,
#  max_time = MaxMinutes * 60,
#  methods = [
#    :generating_set_search,
#    :adaptive_de_rand_1_bin_radiuslimited,
#  ]
#)

# To run XNES on it we need to add a penalty for going outside the (0,1) box.
function penalty(x, range)
  lower = x .< range[1]
  higher = x .> range[2]
  100 * (norm(range[1] - x[lower]) + norm(x[higher] - range[2]))
end

penalized_points15(x) = points15(x) + penalty(x, (0.0, 1.0))

# Run XNES
#best, fitness = bboptimize(penalized_points15, (0.0, 1.0); dimensions = 45,
#  max_time = MaxMinutes * 60, method = :xnes)

#println("points15 fitness = ", points15(best))


text:
# imports
include("iron_line_profile_approx.jl")
using Measures
using DelimitedFiles
gr()

function iron_line_heatmap_approx(
                            tolerance=1e-12,
                            size_multiplier=1,
                            fov=12,
                            dtmax=10000,
                            nbins=100,
                            vars=LinRange(5.0, 85.0, 5),
                            ϵ_3=0,
                            α_13=0,
                            α_22=0,
                            α_52=0
                            )
    # nbins=600
    # vars = LinRange(5.0, 85.0, 400)
    # vars = 5.0:5.0:85.0

    # values to be plotted
    bins_vals = Array{Float64}[]
    x_vals_vals = []

    # generating data for different spins
    for var in vars
        fvar = Float64(var)
        x_vals, bins = iron_line_profile_approx(
                                        tolerance=tolerance,
                                        size_multiplier=size_multiplier,
                                        fov=fov,
                                        dtmax=dtmax,
                                        obs_angle=fvar,
                                        output="data",
                                        normalised=true,
                                        nbins=nbins,
                                        ϵ_3=ϵ_3,
                                        α_13=α_13,
                                        α_22=α_22,
                                        α_52=α_52
                                        )

        push!(bins_vals, bins)
        push!(x_vals_vals, x_vals)
    end

    new = mapreduce(permutedims, vcat, bins_vals)
    hmp = heatmap(
                new,
                xticks=(LinRange(0,nbins,5),string.(LinRange(0,10,5))),
                yticks=(LinRange(0,length(vars),5),string.(round.(vars))),
                xlabel="Energy (keV)",
                ylabel="Observation Angle",
                cb_title="Flux (Arbitrary Units)"
                )
end

hmp = iron_line_heatmap_approx()
writedlm("heatmap_data.txt", new)
display(hmp)
png(hmp, "iron_heatmap.png")


text:
using Test
using GeoParams
@testset "SeismicVelocity.jl" begin

# This tests the MaterialParameters structure
CharUnits_GEO   =   GEO_units(viscosity=1e19, length=1000km);

# Define constant velocities
x1      =   ConstantSeismicVelocity(Vp=8.05km/s, Vs=3.5km/s)
@test  Value(x1.Vp)==8.05km/s
@test  Value(x1.Vs)==3.5km/s

Nondimensionalize!(x1,CharUnits_GEO)
@test NumValue(x1.Vp) ≈ 8.050000000000001e9
@test NumValue(x1.Vs) ≈ 3.5e9

# Compute
@test ComputePwaveVelocity(1.0,1.0, x1) ≈ 8.050000000000001e9
@test ComputeSwaveVelocity(1.0,1.0, x1) ≈ 3.5e9

# Read Phase diagram interpolation object
fname   =   "./test_data/Peridotite.in"
PD_data =   PerpleX_LaMEM_Diagram(fname);
@test PD_data.Vp(1500,1e7) ≈ 6.5290725233303935
@test PD_data.Vs(1500,1e7) ≈ 2.4874400647487658

@test ComputePwaveVelocity(1e7, 1500, PD_data) ≈  6.5290725233303935
@test ComputeSwaveVelocity(1e7, 1500, PD_data) ≈  2.4874400647487658

# Do the same but non-dimensionalize the result
CharDim  =  GEO_units();
PD_data1 =  PerpleX_LaMEM_Diagram(fname, CharDim=CharDim);

rho_ND   =  PD_data1.Rho(Nondimensionalize(1500K,CharDim),Nondimensionalize(1e8*Pa,CharDim))
Vp_ND    =  PD_data1.Vp(Nondimensionalize(1500K,CharDim),Nondimensionalize(1e8*Pa,CharDim))
Vs_ND    =  PD_data1.Vs(Nondimensionalize(1500K,CharDim),Nondimensionalize(1e8*Pa,CharDim))

# redimensionalize and check with value from original structure that did not use non-dimensionalization
@test   ustrip(Dimensionalize(rho_ND,kg/m^3,CharDim)) ≈ PD_data.Rho(1500,1e8)
@test   ustrip(Dimensionalize(Vp_ND, km/s,  CharDim)) ≈ PD_data.Vp(1500,1e8)
@test   ustrip(Dimensionalize(Vs_ND, km/s,  CharDim)) ≈ PD_data.Vs(1500,1e8)

# Test computation of velocity for the whole computational domain, using arrays
MatParam    =   Array{MaterialParams, 1}(undef, 3);
MatParam[1] =   SetMaterialParams(Name="Mantle", Phase=0,
                        SeismicVelocity   = PerpleX_LaMEM_Diagram("test_data/Peridotite.in"));

MatParam[2] =   SetMaterialParams(Name="Crust", Phase=1,
                        SeismicVelocity   = ConstantSeismicVelocity());

MatParam[3] =   SetMaterialParams(Name="UpperCrust", Phase=2,
                        SeismicVelocity   = ConstantSeismicVelocity(Vp=10km/s, Vs=3km/s));

# test computing material properties
Phases              = ones(Int64,400,400)*0;
Phases[:,20:end] .= 1
Phases[:,300:end] .= 2

Vp      = zeros(size(Phases))
Vs      = zeros(size(Phases))
T       =  ones(size(Phases))
P       =  ones(size(Phases))*10

ComputePwaveVelocity!(Vp, Phases, P,T, MatParam)
@test sum(Vp)/400^2 ≈ 8.541562850000005

# test computing material properties when we have PhaseRatios, instead of Phase numbers
PhaseRatio  = zeros(size(Phases)...,length(MatParam));
for i in CartesianIndices(Phases)
    iz = Phases[i]
    I = CartesianIndex(i,iz+1)
    PhaseRatio[I] = 1.0
end

ComputeSwaveVelocity!(Vs, PhaseRatio, P,T, MatParam)
@test sum(Vs)/400^2 ≈ 4.0739837

end


text:
function kernelscore_test(n)
    X = rand(n)
    points = [-1.0, 2.0, 5.0]
    N = size(points,1)
    b = 0.5
    exps::Float64 = 0.0
    tic()
    for i in 1:length(X)
        d = -(X[i]-points).^2./(2*b^2)
        m = minimum(d)
        exps += m-log(b*N)+log(sum(exp(d-m)))
    end
    toc()
    return exps
end

n = 128
if length(ARGS)==1
    n = parse(Int, ARGS[1])
end

println(kernelscore_test(n))

max_stars_repo_name: ehsantn/data-analytics-benchmarks

max_issues_repo_name: ehsantn/data-analytics-benchmarks

max_forks_repo_name: ehsantn/data-analytics-benchmarks


text:
"""
Unicode Normalization, Category constants and documentation

Copyright 2017-2018 Gandalf Software, Inc. Scott P. Jones,
and other contributors to the UTF8proc library (for the documentation of the API)

Licensed under MIT License, see LICENSE.md
"""
module Uni

using ModuleInterfaceTools: m_eval

@enum(Category::Int8,Cn,Lu,Ll,Lt,Lm,Lo,Mn,Mc,Me,Nd,Nl,No,Pc,Pd,Ps,Pe,Pi,Pf,Po,Sm,Sc,Sk,So,
                     Zs,Zl,Zp,Cc,Cf,Cs,Co,Ci,Cm)

# strings corresponding to the category constants
const category_strings = [
    "Other, not assigned",
    "Letter, uppercase",
    "Letter, lowercase",
    "Letter, titlecase",
    "Letter, modifier",
    "Letter, other",
    "Mark, nonspacing",
    "Mark, spacing combining",
    "Mark, enclosing",
    "Number, decimal digit",
    "Number, letter",
    "Number, other",
    "Punctuation, connector",
    "Punctuation, dash",
    "Punctuation, open",
    "Punctuation, close",
    "Punctuation, initial quote",
    "Punctuation, final quote",
    "Punctuation, other",
    "Symbol, math",
    "Symbol, currency",
    "Symbol, modifier",
    "Symbol, other",
    "Separator, space",
    "Separator, line",
    "Separator, paragraph",
    "Other, control",
    "Other, format",
    "Other, surrogate",
    "Other, private use",
    "Invalid, too high",
    "Malformed, bad data",
]

Base.rem(c::Category, t::Type{T}) where {T<:Integer}  = T(Int(c))
Base.rem(c::Category, t::Type{T}) where {T<:Unsigned} = T(Int(c))
Base.:(==)(x::Category, t::Integer) = Int(x) == t
Base.:(==)(t::Integer, x::Category) = Int(x) == t
Base.:(+)(x::Category, t::Integer)  = Category(Int(x) + t)
Base.:(-)(x::Category, y::Integer)  = Category(Int(x) - y)
Base.:(-)(x::Category, y::Category) = Int(x) - Int(y)

const curmod = @static VERSION < v"0.7" ? current_module() : @__MODULE__

for i = Int(typemin(Category)):Int(typemax(Category))
    m_eval(curmod,
           Expr(:macrocall, Symbol("@doc"),
                string("Unicode Category: ", category_strings[i+1]),
                Symbol(Category(i))))
end

## Constants for calling the normalize function

"""The given UTF-8 input is NULL terminated"""
const NULLTERM  = (1<<0)

"""Unicode Versioning Stability has to be respected"""
const STABLE    = (1<<1)

"""Compatibility decomposition (i.e. formatting information is lost)"""
const COMPAT    = (1<<2)

"""Return a result with decomposed characters"""
const COMPOSE   = (1<<3)

"""Return a result with decomposed characters"""
const DECOMPOSE = (1<<4)

"""Strip "default ignorable characters" such as SOFT-HYPHEN or ZERO-WIDTH-SPACE"""
const IGNORE    = (1<<5)

"""Return an error, if the input contains unassigned codepoints"""
const REJECTNA  = (1<<6)

"""
Indicating that NLF-sequences (LF, CRLF, CR, NEL) are representing a line break,
and should be converted to the codepoint for line separation (LS)
"""
const NLF2LS    = (1<<7)

"""
Indicating that NLF-sequences are representing a paragraph break, and
should be converted to the codepoint for paragraph separation (PS)
"""
const NLF2PS    = (1<<8)

"""Indicating that the meaning of NLF-sequences is unknown"""
const NLF2LF    = (NLF2LS | NLF2PS)

"""
Strips and/or converts control characters.
NLF-sequences are transformed into space, except if one of the NLF2LS/PS/LF
options is given. HorizontalTab (HT) and FormFeed (FF) are treated as a
NLF-sequence in this case.  All other control characters are simply removed
"""
const STRIPCC   = (1<<9)

"""Performs unicode case folding, to be able to do a case-insensitive string comparison"""
const CASEFOLD  = (1<<10)

"""
Inserts 0xFF bytes at the beginning of each sequence which is representing
a single grapheme cluster (see UAX#29)
"""
const CHARBOUND = (1<<11)

"""
Lumps certain characters together.

E.g. HYPHEN U+2010 and MINUS U+2212 to ASCII "-". See lump.md for details.

If NLF2LF is set, this includes a transformation of paragraph and line separators
to ASCII line-feed (LF).
"""
const LUMP      = (1<<12)

"""
Strips all character markings.

This includes non-spacing, spacing and enclosing (i.e. accents).
@note This option works only with @ref COMPOSE or
@ref DECOMPOSE
"""
const STRIPMARK = (1<<13)

"""Strip unassigned codepoints"""
const STRIPNA    = (1<<14)

end # module Uni


text:
function solve()
    k = parse(Int, readline())
    ret = 0
    for a=1:k, b=1:k, c=1:k
        ret += gcd(a, b, c)
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using ImputeNaNs
using Base.Test

novalue = NaN
a = [novalue]
b = [novalue, novalue]
c = [1.0, 2.0, 3.0]
d = [1.0, novalue, 3.0]
e = [1.0, 2.0, novalue]
f = [novalue, 2.0, 3.0]

@test isnan(locf(a)[1])
@test isnan(locf(b)[1])
@test isnan(locf(b)[2])
@test locf(c) == c
@test locf(d) == [1.0, 1.0, 3.0]
@test locf(e) == [1.0, 2.0, 2.0]
@test locf(f) == [2.0, 2.0, 3.0]
@test locf(f, true) == [2.0, 2.0, 3.0]
@test isnan(locf(f, false)[1])

novalue = Nullable{Float64}()
a = [Nullable{Float64}(1.0), Nullable{Float64}(2.0), Nullable{Float64}(3.0)]
b = [Nullable{Float64}(1.0), novalue, Nullable{Float64}(3.0)]
c = [Nullable{Float64}(1.0), Nullable{Float64}(2.0), novalue]
d = [novalue, Nullable{Float64}(2.0), Nullable{Float64}(3.0)]

@test all(locf(a) .=== a)
@test all(locf(b) .=== [Nullable{Float64}(1.0), Nullable{Float64}(1.0), Nullable{Float64}(3.0)])
@test all(locf(c) .=== [Nullable{Float64}(1.0), Nullable{Float64}(2.0), Nullable{Float64}(2.0)])
@test all(locf(d) .=== [Nullable{Float64}(2.0), Nullable{Float64}(2.0), Nullable{Float64}(3.0)])
@test isnull(locf(d, false)[1])
@test !isnull(locf(d, true)[1])


text:
module kadi

import PyCall: PyObject, pyimport, PyNULL, pycall
import DataFrames: DataFrame

kadi_states = PyNULL()
kadi_commands = PyNULL()

function __init__()
    copy!(kadi_states, pyimport("kadi.commands.states"))
    copy!(kadi_commands, pyimport("kadi.commands"))
end

struct Commands
    cmds_table::PyObject
    function Commands(cmds_table:PyObject)
        new(cmds_table)
    end
end

struct CmdStates
    states_table::PyObject
    function CmdStates(states_table::PyObject)
        new(states_table)
    end
end

function CmdStates(start, stop; state_keys=nothing, continuity=nothing,
                   reduce=true, merge_identical=false)
    cmd_states = kadi_states.get_states(start=start, stop=stop,
                                        state_keys=state_keys, reduce=reduce,
                                        merge_identical=merge_identical)
    CmdStates(cmd_states)
end

function CmdStates(cmds; state_keys=nothing, continuity=nothing,
                   reduce=true, merge_identical=false)
    cmd_states = kadi_states.get_states(cmds=cmds, state_keys=state_keys,
                                        reduce=reduce,
                                        merge_identical=merge_identical)
    CmdStates(cmd_states)
end

function export_to_dataframe(cmd_states::CmdStates)
    colnames = cmd_states.states_table.colnames
    df = DataFrame()
    getitem = cmd_states.states_table."__getitem__"
    for colname in colnames
        if startswith(colname, "date")
            df[!, colname] = pycall(getitem, Array, colname)
        elseif colname == "trans_keys"
            df[!, colname] = [Set(tk) for tk in getitem(colname)]
        else
            df[!, colname] = getitem(colname)
        end

    end
    df
end

function export_to_dataframe(commands::Commands)
    colnames = commands.cmds_table.colnames
    df = DataFrame()
    getitem = commands.cmds_table."__getitem__"
    for colname in colnames
        if startswith(colname, "date")
            df[!, colname] = pycall(getitem, Array, colname)
        else
            df[!, colname] = getitem(colname)
        end

    end
    df
end


text:
function init_session(session::Session)
    for msg in session.message_queue
        for connection in session.connections
            serialize_websocket(connection, msg)
        end
    end
    empty!(session.message_queue)
    put!(session.js_fully_loaded, true)
end

function Session(connections::Vector{WebSocket}=WebSocket[])
    Session(
        Ref(false),
        connections,
        Dict{String, Tuple{Bool, Observable}}(),
        Dict{Symbol, Any}[],
        Set{Asset}(),
        JSCode[],
        string(uuid4()),
        Channel{Bool}(1),
        init_session
    )
end

session(session::Session) = session

function Base.close(session::Session)
    foreach(close, session.connections)
    empty!(session.connections)
    empty!(session.observables)
    empty!(session.on_document_load)
    empty!(session.message_queue)
    empty!(session.dependencies)
end

function Base.push!(session::Session, observable::Observable)
    if !haskey(session.observables, observable.id)
        session.observables[observable.id] = (true, observable)
        # Register on the JS side by sending the current value
        updater = JSUpdateObservable(session, observable.id)
        # Make sure we update the Javascript values!
        on(updater, observable)
        if isopen(session)
            # If websockets are already open, we need to also update the value
            # to register it with js
            updater(observable[])
        end
    end
end

function Base.push!(session::Session, dependency::Dependency)
    for asset in dependency.assets
        push!(session, asset)
    end
    return dependency
end

function Base.push!(session::Session, asset::Asset)
    push!(session.dependencies, asset)
    if asset.onload !== nothing
        on_document_load(session, asset.onload)
    end
    return asset
end

function Base.push!(session::Session, websocket::WebSocket)
    push!(session.connections, websocket)
    filter!(isopen, session.connections)
    return session
end

"""
    queued_as_script(session::Session)

Returns all queued messages as a script that can be included into html
"""
function queued_as_script(io::IO, session::Session)
    # send all queued messages
    # # first register observables
    observables = Dict{String, Any}()

    for (id, (registered, observable)) in session.observables
        observables[observable.id] = observable[]
    end

    data = Dict("observables" => observables, "messages" => session.message_queue)

    isdir(dependency_path("session_temp_data")) || mkdir(dependency_path("session_temp_data"))

    deps_path = dependency_path("session_temp_data", session.id * ".msgpack")
    open(io -> MsgPack.pack(io, serialize_js(data)), deps_path, "w")
    data_url = url(AssetRegistry.register(deps_path))
    println(io, js"""
    var url = $(data_url);

    var http_request = new XMLHttpRequest();
    http_request.open("GET", url, true);
    http_request.responseType = "arraybuffer";
    var t0 = performance.now();
    http_request.onload = function (event) {
        var t1 = performance.now();
        console.log("download done! " + (t1 - t0) + " milliseconds.");
        var arraybuffer = http_request.response; // Note: not oReq.responseText
        if (arraybuffer) {
            var bytes = new Uint8Array(arraybuffer);
            var data = msgpack.decode(bytes);
            for (let obs_id in data.observables) {
                registered_observables[obs_id] = data.observables[obs_id];
            }
            for (let message in data.messages) {
                process_message(data.messages[message]);
            }
            t1 = performance.now();
            console.log("msg process done! " + (t1 - t0) + " milliseconds.");
            websocket_send({msg_type: JSDoneLoading});
        }else{
            send_warning("Didn't receive any setup data from server.")
        }
    };
    http_request.send(null);
    """)
    empty!(session.message_queue)
end

"""
    send(session::Session; attributes...)

Send values to the frontend via JSON for now
"""
Sockets.send(session::Session; kw...) = send(session, Dict{Symbol, Any}(kw))

function Sockets.send(session::Session, message::Dict{Symbol, Any})
    if isopen(session) && !session.fusing[]
        for connection in session.connections
            serialize_websocket(connection, message)
        end
    else
        push!(session.message_queue, message)
    end
end

fuse(f, has_session) = fuse(f, session(has_session))
function fuse(f, session::Session)
    session.fusing[] = true
    result = f()
    session.fusing[] = false
    if !isempty(session.message_queue)
        send(session; msg_type=FusedMessage, payload=session.message_queue)
        empty!(session.message_queue)
    end
    return result
end

function Base.isopen(session::Session)
    return !isempty(session.connections) && any(isopen, session.connections)
end

"""
    onjs(session::Session, obs::Observable, func::JSCode)

Register a javascript function with `session`, that get's called when `obs` gets a new value.
If the observable gets updated from the JS side, the calling of `func` will be triggered
entirely in javascript, without any communication with the Julia `session`.
"""
function onjs(session::Session, obs::Observable, func::JSCode)
    # register the callback with the JS session
    register_resource!(session, (obs, func))

    send(
        session,
        msg_type = OnjsCallback,
        id = obs.id,
        # eval requires functions to be wrapped in ()
        payload = js"($func)"
    )
end

function onjs(has_session, obs::Observable, func::JSCode)
    onjs(session(has_session), obs, func)
end

"""
    onload(session::Session, node::Node, func::JSCode)

calls javascript `func` with node, once node has been displayed.
"""
function onload(session::Session, node::Node, func::JSCode)
    on_document_load(session, js"""
        ($(func))($node);
    """)
end

"""
    on_document_load(session::Session, js::JSCode)

executes javascript after document is loaded
"""
function on_document_load(session::Session, js::JSCode)
    register_resource!(session, js)
    push!(session.on_document_load, js)
end

"""
    linkjs(session::Session, a::Observable, b::Observable)

for an open session, link a and b on the javascript side. This will also
Link the observables in Julia, but only as long as the session is active.
"""
function linkjs(session::Session, a::Observable, b::Observable)
    # register the callback with the JS session
    onjs(
        session,
        a,
        js"""
        function (value){
            // update_obs will return false once b is gone,
            // so this will automatically deregister the link!
            return update_obs($b, value)
        }
        """
    )
end

function linkjs(has_session, a::Observable, b::Observable)
    linkjs(session(has_session), a, b)
end

"""
    evaljs(session::Session, jss::JSCode)

Evaluate a javascript script in `session`.
"""
function evaljs(session::Session, jss::JSCode)
    register_resource!(session, jss)
    send(session, msg_type = EvalJavascript, payload = jss)
end

function evaljs(has_session, jss::JSCode)
    evaljs(session(has_session), jss)
end

const JS_COMM_CHANNEL = Channel{Dict{String, Any}}(1)
const JS_COMM_OBSERVABLE = Observable(Dict{String, Any}())

"""
    get_js_comm()

Gets the JS communication channels and asserts, that they're in the right state!
"""
function get_js_comm()
    # on first run, we need to register our listener
    if isempty(JS_COMM_OBSERVABLE.listeners)
        # whenever we get a new value, put it in the channel
        on(JS_COMM_OBSERVABLE) do value
            put!(JS_COMM_CHANNEL, value)
        end
    end

    @assert length(JS_COMM_OBSERVABLE.listeners) >= 1
    js_listeners = @view JS_COMM_OBSERVABLE.listeners[2:end]
    # There should be only the listeners needed to forward values to the Frontend!
    @assert all(x-> x isa JSUpdateObservable, js_listeners) "Someone else registered to the comm channel. Nobody should do this!!"
    @assert isopen(JS_COMM_CHANNEL)
    @assert !isready(JS_COMM_CHANNEL) "Channel is still containing a value from another run!"
    return JS_COMM_OBSERVABLE, JS_COMM_CHANNEL
end

"""
    evaljs_value(session::Session, js::JSCode)

Evals `js` code and returns the jsonified value.
Blocks until value is returned. May block indefinitely, when called with a session
that doesn't have a connection to the browser.
"""
function evaljs_value(session::Session, js; error_on_closed=true, time_out=10.0)
    if error_on_closed && !isopen(session)
        error("Session is not open and would result in this function to indefinitely block.
        It may unblock, if the browser is still connecting and opening the session later on. If this is expected,
        you may try setting `error_on_closed=false`")
    end
    comm, comm_channel = get_js_comm()
    js_with_result = js"""
    try{
        var result = $(js);
        update_obs($(comm), {result: result});
    }catch(e){
        update_obs($(comm), {error: e.toString()});
    }
    """
    evaljs(session, js_with_result)
    # TODO, have an on error callback, that triggers when evaljs goes wrong
    # (e.g. because of syntax error that isn't caught by the above try catch!)
    wait_timeout(()-> isready(comm_channel), "Waited for $(time_out) seconds for JS to return, but it didn't!", time_out)
    value = take!(comm_channel)
    if haskey(value, "error")
        error(value["error"])
    else
        return value["result"]
    end
end

"""
    register_resource!(session::Session, domlike)

Walks dom like structures and registers all resources (Observables, Assets Depencies)
with the session.
"""
register_resource!(session::Session, @nospecialize(jss)) = nothing # do nothing for unknown type

function register_resource!(session::Session, list::Union{Tuple, AbstractVector, Pair})
    for elem in list
        register_resource!(session, elem)
    end
end

function register_resource!(session::Session, jss::JSCode)
    register_resource!(session, jss.source)
end

function register_resource!(session::Session, asset::Union{Asset, Dependency, Observable})
    push!(session, asset)
end

function register_resource!(session::Session, node::Node)
    walk_dom(session, node) do x
        register_resource!(session, x)
    end
end


text:
module MultilinearOpt

import JuMP, MathOptInterface, PiecewiseLinearOpt

using LinearAlgebra: Symmetric, eigvals

export HyperRectangle, MultilinearFunction, Discretization, outerapproximate, relaxbilinear!

const MOI = MathOptInterface
const MOIU = MOI.Utilities

struct HyperRectangle{D}
    l::NTuple{D,Float64}
    u::NTuple{D,Float64}
end

mutable struct MultilinearFunction{D}
    bounds::HyperRectangle{D}
    f::Function
end
(mlf::MultilinearFunction)(args...) = mlf.f(args...)

mutable struct Discretization{D}
    d::NTuple{D,Vector{Float64}}
end
Discretization(d...) = Discretization(map(t->convert(Vector{Float64},t), d))

mutable struct MultilinearData
    product_dict::Dict{JuMP.UnorderedPair{JuMP.VariableRef}, JuMP.VariableRef}

    MultilinearData() = new(Dict{JuMP.UnorderedPair{JuMP.VariableRef}, JuMP.VariableRef}())
end

function outerapproximate(m::JuMP.Model, x::NTuple{D,JuMP.VariableRef}, mlf::MultilinearFunction{D}, disc::Discretization, method) where {D}
    vform = (method in (:Logarithmic1D,:Logarithmic2D,:ZigZag1D,:ZigZag2D,:Unary))
    r = length(disc.d)
    @assert r == 2
    V = collect(Base.product(disc.d...))
    T = map(t -> mlf.f(t...), V)
    z = JuMP.@variable(m, base_name="z")
    if vform
        PiecewiseLinearOpt.initPWL!(m)
        λ = JuMP.@variable(m, [V], lower_bound=0, upper_bound=1, base_name="λ")
        JuMP.@constraint(m, sum(λ) == 1)
        for i in 1:r
            JuMP.@constraint(m, sum(λ[v]*v[i] for v in V) == x[i])
        end
        JuMP.@constraint(m, sum(λ[v]*mlf.f(v...) for v in V)  == z)
        for i in 1:r
            I = disc.d[i]
            n = length(I)-1
            n > 1 || continue
            γ = JuMP.@expression(m, [j=1:(n+1)], sum(λ[v] for v in V if v[i] == I[j]))
            if method == :Logarithmic1D || method == :Logarithmic2D
                PiecewiseLinearOpt.sos2_logarithmic_formulation!(m, γ)
            elseif method == :Unary
                PiecewiseLinearOpt.sos2_mc_formulation!(m, γ) # TODO: where's γ supposed to come from?
            elseif method == :ZigZag1D || method == :ZigZag2D
                PiecewiseLinearOpt.sos2_zigzag_general_integer_formulation!(m, γ)
            else
                throw(ArgumentError("Unrecognized method: $method"))
            end
        end
    else
        # Methods from Misener only work for bilinear terms, and discretized along only one dimension
        @assert minimum(map(length, disc.d)) == 2
        iˣ = length(disc.d[1]) == 2 ? 2 : 1 # direction we discretize along
        iʸ = iˣ == 1 ? 2 : 1
        I = disc.d[iˣ]
        NP = length(I)-1
        xˡ, xᵘ = minimum(I), maximum(I)
        yˡ, yᵘ = minimum(disc.d[iʸ]), maximum(disc.d[iʸ])
        a = I[2] - I[1] # should assert all lengths are the same
        x, y = x[iˣ], x[iʸ]
        if method == :MisenerLinear
            λ  = JuMP.@variable(m, [1:NP], Bin)
            Δy = JuMP.@variable(m, [1:NP], lower_bound=0, upper_bound=yᵘ-yˡ)
            JuMP.@constraint(m, sum(λ) == 1)
            JuMP.@constraints(m, begin
                xˡ + sum(a*(nP-1)*λ[nP] for nP in 1:NP) ≤ x
                x ≤ xˡ + sum(a*nP*λ[nP] for nP in 1:NP)
            end)
            JuMP.@constraints(m, begin
                y == yˡ + sum(Δy[nP] for nP in 1:NP)
                [nP=1:NP], Δy[nP] ≤ (yᵘ-yˡ)*λ[nP]
            end)
            JuMP.@constraints(m, begin
                z ≥ x*yˡ + sum((xˡ+a*(nP-1))* Δy[nP]                for nP in 1:NP)
                z ≥ x*yᵘ + sum((xˡ+a* nP   )*(Δy[nP]-(yᵘ-yˡ)*λ[nP]) for nP in 1:NP)
                z ≤ x*yˡ + sum((xˡ+a* nP   )* Δy[nP]                for nP in 1:NP)
                z ≤ x*yᵘ + sum((xˡ+a*(nP-1))*(Δy[nP]-(yᵘ-yˡ)*λ[nP]) for nP in 1:NP)
            end)
        elseif method == :MisenerLog1
            NL = ceil(Int, log2(NP))
            λ    = JuMP.@variable(m, [1:NL], Bin)
            Δy   = JuMP.@variable(m, [1:NP], lower_bound=0, upper_bound=yᵘ-yˡ)
            λhat = JuMP.@variable(m, [1:NP], lower_bound=0, upper_bound=1)
            JuMP.@constraints(m, begin
                xˡ + sum(2^(NL-nL)*a*λ[nL] for nL in 1:NL)     ≤ x
                xˡ + sum(2^(NL-nL)*a*λ[nL] for nL in 1:NL) + a ≥ x
            end)
            JuMP.@constraint(m, sum(λhat) == 1)
            for nL in 1:NL
                JuMP.@constraints(m, begin
                    sum(λhat[nP] for nP in 1:NP if mod(floor((nP-1)/2^(NL-nL)),2)==0) ≤ 1 - λ[nL]
                    sum(λhat[nP] for nP in 1:NP if mod(floor((nP-1)/2^(NL-nL)),2)==1) ≤     λ[nL]
                end)
            end
            for nP in 1:NP
                JuMP.@constraint(m, Δy[nP] ≤ (yᵘ-yˡ)*λhat[nP])
            end
            JuMP.@constraint(m, y == yˡ + sum(Δy))
            JuMP.@constraints(m, begin
                z ≥ x*yˡ + sum((xˡ+a*(nP-1))* Δy[nP]                   for nP in 1:NP)
                z ≥ x*yᵘ + sum((xˡ+a* nP   )*(Δy[nP]-(yᵘ-yˡ)*λhat[nP]) for nP in 1:NP)
                z ≤ x*yˡ + sum((xˡ+a* nP   )* Δy[nP]                   for nP in 1:NP)
                z ≤ x*yᵘ + sum((xˡ+a*(nP-1))*(Δy[nP]-(yᵘ-yˡ)*λhat[nP]) for nP in 1:NP)
            end)
        elseif method == :MisenerLog2
            NL = ceil(Int, log2(NP))
            λ  = JuMP.@variable(m, [1:NL], Bin)
            Δy = JuMP.@variable(m, [1:NL], lower_bound=0, upper_bound=yᵘ-yˡ)
            s  = JuMP.@variable(m, [1:NL], lower_bound=0, upper_bound=yᵘ-yˡ)
            JuMP.@constraints(m, begin
                xˡ +     sum(2^(nL-1)*a*λ[nL] for nL in 1:NL) ≤ x
                xˡ + a + sum(2^(nL-1)*a*λ[nL] for nL in 1:NL) ≥ x
                xˡ + a + sum(2^(nL-1)*a*λ[nL] for nL in 1:NL) ≤ xᵘ
            end)
            for nL in 1:NL
                JuMP.@constraints(m, begin
                    Δy[nL] ≤ (yᵘ-yˡ)*λ[nL]
                    Δy[nL] == (y-yˡ) - s[nL]
                    s[nL] ≤ (yᵘ-yˡ)*(1-λ[nL])
                end)
            end
            JuMP.@constraints(m, begin
                z ≥ x*yˡ +  xˡ   *(y-yˡ) + sum(a*2^(nL-1)* Δy[nL]             for nL in 1:NL)
                z ≥ x*yᵘ + (xˡ+a)*(y-yᵘ) + sum(a*2^(nL-1)*(Δy[nL]-(yᵘ-yˡ)*λ[nL]) for nL in 1:NL)
                z ≤ x*yˡ + (xˡ+a)*(y-yˡ) + sum(a*2^(nL-1)* Δy[nL]             for nL in 1:NL)
                z ≤ x*yᵘ +  xˡ   *(y-yᵘ) + sum(a*2^(nL-1)*(Δy[nL]-(yᵘ-yˡ)*λ[nL]) for nL in 1:NL)
            end)
        else
            throw(ArgumentError("Unrecognized method: $method"))
        end
    end
    z
end

function gramian(expr::JuMP.GenericQuadExpr{T}) where T
    vars = unique(Iterators.flatten((var1, var2) for (coeff, var1, var2) in JuMP.quad_terms(expr)))
    varindices = Dict(v => i for (i, v) in enumerate(vars))
    n = length(vars)
    gramian = zeros(T, n, n)
    for (coeff, var1, var2) in JuMP.quad_terms(expr)
        ind1, ind2 = varindices[var1], varindices[var2]
        row, col = extrema((ind1, ind2))
        gramian[row, col] += row == col ? coeff : coeff / 2
    end
    Symmetric(gramian), vars
end

ispossemidef(mat) = all(eigvals(mat) .>= -1e-10)
isconvex(x) = error("Could not determine convexity.")
isconvex(expr::JuMP.GenericAffExpr) = true
isconvex(expr::JuMP.GenericQuadExpr) = ispossemidef(first(gramian(expr)))
isconvex(set::MOI.LessThan) = true
isconvex(set::MOI.GreaterThan) = true
isconvex(set::MOI.EqualTo) = true
isconvex(set::MOI.Interval) = true

function isconvex(constr::JuMP.ScalarConstraint)
    f = JuMP.jump_function(constr)
    set = JuMP.moi_set(constr)
    if f isa JuMP.GenericAffExpr
        return isconvex(set)
    elseif f isa JuMP.GenericQuadExpr
        interval = MOI.Interval(set)
        l, u = interval.lower, interval.upper
        if l == -Inf && isfinite(u)
            return isconvex(f)
        elseif isfinite(l) && u == Inf
            return isconvex(-f)
        elseif isfinite(l) && isfinite(u)
            return all(iszero, coeff for (coeff, _, _) in JuMP.quad_terms(f))
        elseif l == -Inf && u == Inf
            return true
        else
            error("Should never get here.")
        end
    else
        error("Function type $(typeof(f)) not recognized.")
    end
end

function relaxbilinear!(m::JuMP.Model; method=:Logarithmic1D, disc_level::Int = 9, constraints_to_skip=())
    # replace each bilinear term in (nonconvex) quadratic constraints with outer approx
    product_dict = Dict{JuMP.UnorderedPair{JuMP.VariableRef}, Any}() #m.ext[:Multilinear].product_dict
    sense = JuMP.objective_sense(m)
    if sense !== MOI.FEASIBILITY_SENSE
        @assert sense ∈ (MOI.MIN_SENSE, MOI.MAX_SENSE)
        obj = JuMP.objective_function(m)
        min_obj = sense == MOI.MIN_SENSE ? obj : -obj
        if !isconvex(min_obj)
            min_aff = linearize_quadratic!(m, min_obj, product_dict, method, disc_level)
            aff = sense == MOI.MIN_SENSE ? min_aff : -min_aff
            JuMP.set_objective_function(m, aff)
        end
    end
    linearized_quad_constrs = JuMP.ConstraintRef[]
    for (F, S) in JuMP.list_of_constraint_types(m)
        if F <: JuMP.GenericQuadExpr
            for constr in JuMP.all_constraints(m, F, S)
                constr ∈ constraints_to_skip && continue
                q = JuMP.constraint_object(constr)
                if !isconvex(q)
                    f = JuMP.jump_function(q)
                    aff = linearize_quadratic!(m, f, product_dict, method, disc_level)
                    set = JuMP.moi_set(q)
                    interval = MOI.Interval(set)
                    lb = interval.lower - JuMP.constant(aff)
                    ub = interval.upper - JuMP.constant(aff)
                    aff.constant = 0
                    JuMP.@constraint(m, lb <= aff <= ub)
                    push!(linearized_quad_constrs, constr)
                end
            end
        end
    end
    for constr in linearized_quad_constrs
        JuMP.delete(m, constr)
    end
    nothing
end

function linearize_quadratic!(m::JuMP.Model, t::JuMP.QuadExpr, product_dict::Dict, method::Symbol, disc_level::Int)
    aff = copy(t.aff)
    for (coeff, x, y) in JuMP.quad_terms(t)
        z = get!(product_dict, JuMP.UnorderedPair(x, y)) do
            # TODO: better relaxation for the case x == y
            if JuMP.is_fixed(x) && JuMP.is_fixed(y)
                return JuMP.fix_value(x) * JuMP.fix_value(y)
            elseif JuMP.is_fixed(x)
                return JuMP.fix_value(x) * y
            elseif JuMP.is_fixed(y)
                return x * JuMP.fix_value(y)
            end
            lˣ, lʸ = JuMP.lower_bound(x), JuMP.lower_bound(y)
            uˣ, uʸ = JuMP.upper_bound(x), JuMP.upper_bound(y)
            @assert isfinite(lˣ) && isfinite(lʸ) && isfinite(uˣ) && isfinite(uʸ)
            hr = HyperRectangle((lˣ,lʸ), (uˣ,uʸ))
            mlf = MultilinearFunction(hr, (a,b) -> a*b)
            disc_levelˣ = lˣ == uˣ ? 1 : disc_level # if variable is fixed, no need to discretize
            disc_levelʸ = lʸ == uʸ ? 1 : (!(method in (:Logarithmic2D,:ZigZag2D)) ? 2 : disc_level)
            # disc = Discretization(range(lˣ, stop=uˣ, length=disc_levelˣ),  range(lʸ, stop=uʸ, length=disc_levelʸ))
            disc = Discretization(range(lˣ, stop=uˣ, length=disc_levelˣ),  range(lʸ, stop=uʸ, length=disc_levelʸ))
            outerapproximate(m, (x, y), mlf, disc, method)
        end
        JuMP.add_to_expression!(aff, coeff * z)
    end
    aff
end

end # module


text:
cd(@__DIR__)
using Pkg
Pkg.activate(".")
##
using StochasticPrograms
using GLPK
using Pipe
using Random
using Plots
using JuMP
using DataFrames
using CSV
using Cbc
##

include("samplers.jl")
include("plot_utils.jl")

##
# Define the time interval
t_s = 1
t_f = 168
offset = 6531
timesteps = t_s:t_f

F_range = (100., 200.)
f_samp = F_sampler(timesteps, F_range)
# Define weather and demand data
pv = CSV.read("../basic_example.csv", DataFrame)[timesteps.+offset, 3]
wind = CSV.read("../basic_example.csv", DataFrame)[timesteps.+offset, 4]
demand = CSV.read("../basic_example.csv", DataFrame)[timesteps.+offset, 2] ./80
##
# Note on signs: s_xi<0 means that energy is requested, s_xi>0 means an additional consumption
energy_model = @stochastic_model begin
    @stage 1 begin
        @parameters begin
            asset_lifetime = 10.

            c_pv = 100.
            c_wind = 1000.
            c_storage = 100.
            inv_budget = 500000000.
            flexible_demand = 0.
        end
        time_scale = asset_lifetime *365*24/length(timesteps)
        # Investments:
        @decision(model, u_pv >= 0)
        @decision(model, u_wind >= 0)
        @decision(model, u_storage >= 0)
        @constraint(model, c_pv*u_pv+c_wind*u_wind+c_storage*u_storage<=inv_budget)
        # Grid connection
        @decision(model, gci[t in timesteps] >= 0)
        @decision(model, gco[t in timesteps] >= 0)
        # Flexible demand
        @decision(model, fl_dem[t in timesteps] >= 0)
        @constraint(model, sum(fl_dem) == flexible_demand)
        # Storage model
        @decision(model, sto_in[t in timesteps] >= 0) #into the bus from storage
        @decision(model, sto_out[t in timesteps] >= 0)
        @constraint(model, [t in timesteps], -0.5*u_storage <= -sum(sto_in[t_s:t])+sum(sto_out[t_s:t]))
        @constraint(model, [t in timesteps], -sum(sto_in[t_s:t])+sum(sto_out[t_s:t]) <= 0.5*u_storage)
        @constraint(model, -sum(sto_in)+sum(sto_out) == 0)
        # Energy balance
        @constraint(model, [t in timesteps],
        gci[t]-gco[t]+u_pv*pv[t]+u_wind*wind[t]-demand[t]-fl_dem[t]+sto_in[t]-sto_out[t]==0)
        # Investment costs
        @objective(model, Min, (u_pv*c_pv+u_wind*c_wind+u_storage*c_storage)/time_scale)
    end
    @stage 2 begin
        @parameters begin
            c_i = .03
            c_o = .01
            c_flex = .5
            #F = 10.
            c_sto_op = 0.00001

        end
        @uncertain t_xi s_xi F_xi from F_scenario #t_xi the time of flexibility demand, s_xi - sign (±1 or 0)
        @known(model, u_pv)
        @known(model, u_wind)
        @known(model, u_storage)
        @known(model, gci)
        @known(model, gco)
        @known(model, sto_in)
        @known(model, sto_out)
        @known(model, fl_dem)
        # Post event components
        @recourse(model, gci2[t in t_s:t_f] >= 0)
        @recourse(model, gco2[t in t_s:t_f] >= 0)
        @recourse(model, sto_in2[t in timesteps] >= 0)
        @recourse(model, sto_out2[t in timesteps] >= 0)

        @constraint(model, [t in t_s:t_xi], gci[t] == gci2[t])
        @constraint(model, [t in t_s:t_xi], gco[t] == gco2[t])

        @constraint(model, [t in timesteps], -0.5*u_storage <= -sum(sto_in2[t_s:t])+sum(sto_out2[t_s:t]))
        @constraint(model, [t in timesteps], -sum(sto_in2[t_s:t])+sum(sto_out2[t_s:t]) <= 0.5*u_storage)

        @constraint(model, -sum(sto_in2)+sum(sto_out2) == 0)

        #Flexible demand
        @recourse(model, fl_dem2[t in timesteps]>=0)
        @constraint(model, [t in t_s:(t_xi-1)], fl_dem[t] == fl_dem2[t])
        @constraint(model, sum(fl_dem2) == sum(fl_dem))
        # Put storage at time of event into same state
        @constraint(model, [t in t_s:(t_xi-1)], sto_in[t] == sto_in2[t])
        @constraint(model, [t in t_s:(t_xi-1)], sto_out[t] == sto_out2[t])

        # Event energy balance
        # The storage and other fast acting components use the recourse variables here.
        # They provide the balance. Grid connection is not allowed, as we are suporting the grid here.
        @constraint(model, gci[t_xi]-gco[t_xi]+u_pv*pv[t_xi]+u_wind*wind[t_xi]-demand[t_xi]-fl_dem2[t_xi] +sto_in2[t_xi]-sto_out2[t_xi]+F_xi*s_xi==0)
        # Post event energy balance
        @constraint(model, [t in (t_xi+1):t_f],
        gci2[t]-gco2[t]+u_pv*pv[t]+u_wind*wind[t]-demand[t]-fl_dem2[t]+sto_in2[t]-sto_out2[t]==0)
        @objective(model, Min, c_i*sum(gci2)-c_o*sum(gco2)+c_sto_op*sum(sto_in2) + c_sto_op*sum(sto_out2) - c_flex*F_xi)
    end
end
##
#c_i = 0.5, c_o = 0.04, c_wind = 10., c_pv = 10., c_storage = 1. # very expensive grid connection, cheap storage
#xi = [simple_scenario(5,1), simple_scenario(5,-1)]

sp0 = instantiate(energy_model, [F_scenario(t_f, 0, 0.)], c_pv = 700., c_wind = 2000., c_i = 0.3, c_o = 0.05, c_storage = 600., flexible_demand = 0., optimizer = GLPK.Optimizer)

optimize!(sp0)

od0 = optimal_decision(sp0)
ov0 = objective_value(sp0)

plot_results(sp0, pv, wind, demand, s = 1, stage_1 = [:gci, :gco, :fl_dem], stage_2 = [:gci2, :gco2, :fl_dem2],)

##
sp = instantiate(energy_model, f_samp, 100,  c_pv = 700., c_wind = 2000., c_i = 0.3, c_o = 0.05, c_storage = 600., flexible_demand = 0., optimizer = GLPK.Optimizer)
optimize!(sp)

od = optimal_decision(sp)
ov = objective_value(sp)

##
#sp = instantiate(energy_model, [simple_scenario(1,0)], F=100., optimizer = Cbc.Optimizer)

##
for s in 1:5
    plot_results(sp, pv, wind, demand, s = s, stage_1 = [:gci, :gco, :fl_dem], stage_2 = [:gci2, :gco2, :fl_dem2],)
end
##
# Main result
println("Termination status: $(termination_status(sp))")
println("Objective value: $(objective_value(sp))")
println("Optimal decision: $(optimal_decision(sp))")

##
# Investment result
u_pv = sp[1, :u_pv]
u_wind = sp[1, :u_wind]
u_storage = sp[1, :u_storage]
# First stage
println("value(u_pv) = $(value(u_pv))")
println("value(u_wind) = $(value(u_wind))")
println("value(u_storage) = $(value(u_storage))")

# Second stage
for s in 1:3
    println("Objective value in scenario $s: $(objective_value(sp, s))")
    #println("Optimal recourse in scenario $s: $(optimal_recourse_decision(sp, s))")
    plot_results(sp, pv, wind, demand, s = s)
end

##

scen = scenarios(sp)
print(scen)
##
function test_decision(p, decision, timesteps)
    infeasible_count = 0
    costs = []
    #scenario_results = []
    for sign in [-1,1]
        for t in timesteps
            scenario = F_scenario(t,sign,150.)
            try evaluate_decision(p, decision, scenario);
            catch e;
                infeasible_count += 1;
                append!(costs, Inf)
                #append!(scenario_results, [Inf, scenario.t_xi, scenario.s_xi, scen.F_xi]);
            end
            append!(costs, evaluate_decision(p, decision, scenario))
            #append!(scenario_results, [evaluate_decision(p, decision, scenario), scenario.t_xi, scenario.s_xi, scenario.F_xi]);
        end
    end
    N = length(timesteps)
    print(infeasible_count/N/2.)

    scenario_results = hcat(costs, vcat(timesteps, timesteps), vcat(-1*ones(N), ones(N)))
    return scenario_results
end

function test_decision_variate_F(p, decision, timesteps; F_step = 50., F_max = 300.)
    costs = zeros(2*length(timesteps))
    F_potential = zeros(2*length(timesteps))
    #scenario_results = []
    for Sign in [-1,1]
        for t in timesteps
            i = t
            if Sign == 1
                i = t + length(timesteps)
            end
            for F in F_step:F_step:F_max
                scenario = F_scenario(t,Sign,F*1.)
                c = evaluate_decision(p, decision, scenario)
                if c!= Inf
                    F_potential[i] = F*Sign
                    costs[i] = c
                else
                    if costs[i] == 0
                        costs[i] = Inf
                        F_potential[i] = F*Sign
                        break
                    end
                end
            end
        end
    end
    print(F_potential)
    scenario_results = hcat(costs, vcat(timesteps, timesteps), F_potential)
    return scenario_results
end

##

scen_results = test_decision_variate_F(sp, od, timesteps)

scen_results0 = test_decision_variate_F(sp0, od0, timesteps)

positive_flexibility, negative_flexibility, positive_potential, negative_potential = plot_flexibility(scen_results0, ov0)

plot_flexibility_cost(scen_results0, ov0)

plot(positive_potential)
plot!(negative_potential)

plot_flexibility(scen_results, ov0)

plot_flexibility(scen_results0, ov0)

##
scen_results = test_decision_variate_F(sp, od, timesteps, F_step = 300., F_max = 3000.)

plt = plot()
plot!(plt, positive_flexibility./positive_potential, label = "price of positive flexibility")
plot!(plt, negative_flexibility./negative_potential, label = "price of negative flexibility", legend = :topleft)
plot!(twinx(), positive_potential, fillrange = 0, fillalpha = 0.35, label = "positive flexibility potential", ylimits = (-300., 300.))
plot!(twinx(), negative_potential, color = :red, fillrange = 0, fillalpha = 0.35, label = "negative flexibility potential", ylimits = (-300., 300.))
#display(plot(plt_cost, plt_pot, layout = (2,1)))
display(plt)


text:
# test page for checking if site is live
# testing post and get requests
function route_test_page_get_post!(progress_test,html_coast_introduction)

    # allow cors headers because chrome complains about these otherwise
    Genie.config.cors_headers["Access-Control-Allow-Methods"] ="GET, POST"
    Genie.config.cors_headers["Access-Control-Allow-Origin"] = "*"

    # simple message to say hi
    route("/") do
        print("Run get script COAST \n")
        html_coast_introduction
     end

     # use e.g. curl to test
     route("/", method = POST) do
      print(jsonpayload())
       aa = jsonpayload()["name"]
       print("testing progress script! \n")
       progress_test[1] += 1.0
       return Genie.Renderer.Json.json("Tested input: $(aa) ")
     end

     # use e.g. curl to test
     route("/testpost", method = GET) do

      return "<p>Registered GET. Current progress is $(progress_test)%<p>"
    end

  end


text:
using BenchmarkTools, NTupleKmers, BioSequences

SUITE = BenchmarkGroup()

SUITE["creation"] = BenchmarkGroup(["constructors"])

# Add some benchmarks to the "creation" group
dnaseq = randdnaseq(63)
SUITE["creation"]["63mer"] = @benchmarkable DNAKmer{63,2}($dnaseq)

onem   = NTupleKmers.DNAKmer{31, 1}(LongSequence{DNAAlphabet{2}}(randdnaseq(31)))
twom   = NTupleKmers.DNAKmer{63, 2}(LongSequence{DNAAlphabet{2}}(randdnaseq(63)))
threem = NTupleKmers.DNAKmer{95, 3}(LongSequence{DNAAlphabet{2}}(randdnaseq(95)))
fourm  = NTupleKmers.DNAKmer{127, 4}(LongSequence{DNAAlphabet{2}}(randdnaseq(127)))

m = NTupleKmers.DNAKmer{63, 2}(dnaseq)
oldm = BigDNAMer{63}(dnaseq)

@benchmark shiftright($m)
@benchmark shiftright($oldm)


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, Pkg

name = "SQLCipher"
version = v"4.4.3"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/sqlcipher/sqlcipher.git", "0663d8500204e14bd2bb0ca25162d91e4555528d")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/sqlcipher

# use same flags as SQLite_jll
# plus those provided at https://github.com/sqlcipher/sqlcipher
export CPPFLAGS="-DSQLITE_ENABLE_COLUMN_METADATA=1 \
                 -DSQLITE_ENABLE_UNLOCK_NOTIFY \
                 -DSQLITE_ENABLE_DBSTAT_VTAB=1 \
                 -DSQLITE_ENABLE_FTS3_TOKENIZER=1 \
                 -DSQLITE_SECURE_DELETE \
                 -DSQLITE_MAX_VARIABLE_NUMBER=250000 \
                 -DSQLITE_MAX_EXPR_DEPTH=10000"
./configure --prefix=${prefix} \
    --build=${MACHTYPE} \
    --host=${target} \
    --enable-tempstore=yes \
    --disable-static \
    --enable-fts3 \
    --enable-fts4 \
    --enable-fts5 \
    --enable-rtree \
    --enable-json1 \
    CFLAGS="-DSQLITE_HAS_CODEC" \
    LDFLAGS="-L${libdir}" \
    LDFLAGS="-lcrypto"

make -j${nproc}
make install

# SQLCipher and SQLite licenses
install_license LICENSE*
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)
filter!(p -> libc(p) != "musl", platforms)
filter!(p -> !Sys.isfreebsd(p) && !Sys.iswindows(p), platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libsqlcipher", :libsqlcipher),
    ExecutableProduct("sqlcipher", :sqlcipher)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # Required for amalgamation, could not build without it
    HostBuildDependency("Tcl_jll"),
    Dependency(PackageSpec(name="OpenSSL_jll", uuid="458c3c95-2e84-50aa-8efc-19380b2a3a95")),
]

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

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
__precompile__()

module MendelPlots

using Compose, Statistics, DataFrames, Distributions, Reexport, DataFrames, TextParse, Colors
@reexport using Gadfly
import Cairo

export
    qq,
    manhattan,
    ppoints

include("gwasplots.jl")

end


text:
module GetFinancialSeries

using Requests

export fred, truefx, google, quandl #, yahoo

include("fred/lookup.jl")
include("fred/request.jl")
include("fred/update.jl")

include("truefx/authorization.jl")
include("truefx/lookup.jl")
include("truefx/request.jl")
include("truefx/update.jl")

include("google/lookup.jl")
include("google/request.jl")
include("google/update.jl")

include("quandl/lookup.jl")
include("quandl/request.jl")
include("quandl/update.jl")

#=
include("yahoo/lookup.jl")
include("yahoo/request.jl")
include("yahoo/update.jl")
=#

end # module