text:
#  fdp.jl
#  fused dot product and related methods.
#  this code was funded by Etaphase, inc, as part of DARPA TRADES program
#  award number BAA-16-39

doc"""
  SigmoidNumbers::Quire

  represents the 'quire' data structure.  This is a scratchpad of unsigned
  integers which represents a mutable fixed point value, useful as an
  accumulator for fused floating point operations.

  Currently, no shortcut operations are supported.
"""

type Quire
  fixed_point_value::Vector{UInt64}
  infinity::Bool
end

#a few properties:
maxpos(::Type{Quire}) = 2047
minpos(::Type{Quire}) = -2048
Base.isinf(q::Quire) = q.infinity
iszero(q::Quire) = mapreduce((n) -> (n == zero(UInt64)), &, true, q.fixed_point_value )

#a few convenience functions.
maximum_exponent{N,ES}(::Type{Posit{N,ES}}) = (2^(ES) * N - 2)
minimum_exponent{N,ES}(::Type{Posit{N,ES}}) = -(2^(ES) * N - 2)

#create a zeroed quire.
function (::Type{Quire}){N,ES}(::Type{Posit{N,ES}})
  # IN THE FUTURE:
  # calculate the  number of digits we'll need.  The top exponent is going to be
  # 2^(ES) * (N - 2).  Double that and you get 2^(ES + 1) * (N - 2)..  Then add padding.
  # bottom exponent is going to be 2^(ES) * (N - 2).

  #for now, default to accomodating an accumulator for {64, 4}.  Layout is as follows:
  # 32 words . 32 words
  Quire(zeros(UInt64, 64), false)
end

type fquire{F}; q::F; end
(::Type{Quire})(::Type{Float64}) = fquire{Float64}(zero(Float64))
(::Type{Quire})(::Type{Float32}) = fquire{Float32}(zero(Float32))

zero!{T}(q::fquire{T}) = q.q = zero(T)
#zeroes out the fused dot product accumulator.  Just throw out the old array.
function zero!(q::Quire)
  q.infinity = false
  for idx = 1:64
    q.fixed_point_value[idx] = zero(UInt64)
  end
end

doc"""
  inf!(q::Quire) forces the quire to carry an infinite value.
"""
inf!(q::Quire) = (q.infinity = true; q)
inf!{T}(q::fquire{T}) = (q.q = T(Inf); q)

doc"""
  isnegative(q::Quire) checks if the quire contains a negative value.
"""
isnegative(q::Quire) = (last(q.fixed_point_value) & 0x8000_0000_0000_0000) != 0

#this function builds a posit value based on some passed sign/exponent/fraction
#parameters.

function (::Type{Posit{N,0}}){N}(sign::Bool, exp::Int64, frac::UInt64)
  frac = frac >> 3
  #In both positive and negative cases, there's a straightforward translation of
  #exponent to the requisite shift.
  exp -= sign

  shift = (exp < 0) ? (-exp - 1) : exp
  frac |= (sign ⊻ (exp < 0)) ? 0x2000_0000_0000_0000 : 0xC000_0000_0000_0000

  sfrac = reinterpret(Int64, frac)
  frac = reinterpret(UInt64, sfrac >> shift)

  frac = (sign) ? (frac | 0x8000_0000_0000_0000) : (frac & 0x7FFF_FFFF_FFFF_FFFF)

  __round(reinterpret(Posit{N,0}, frac))
end

function (::Type{Posit{N,ES}}){N,ES}(sign::Bool, exp::Int64, frac::UInt64)

  (exp < -2048) && return zero(Posit{N,ES})

  frac = frac >> (3 + ES)
  #In both positive and negative cases, there's a straightforward translation of
  #exponent to the requisite shift.

  exp -= sign
  reg = exp >> ES

  exp = ((1 << ES) - 1) & exp
  exp = (sign ? ((1 << ES) - 1) - exp : exp)

  shift = (reg < 0) ? (-reg - 1) : reg
  frac |= (sign ⊻ (reg < 0)) ? 0x2000_0000_0000_0000 : 0xC000_0000_0000_0000
  exp_shifted = exp << (61 - ES)
  frac |= exp_shifted

  sfrac = reinterpret(Int64, frac)
  frac = reinterpret(UInt64, sfrac >> shift)

  frac = (sign) ? (frac | 0x8000_0000_0000_0000) : (frac & 0x7FFF_FFFF_FFFF_FFFF)

  __round(reinterpret(Posit{N,ES}, frac))
end

# an iterator that iterates over the quire indices and returns the fixed point
# cell contents.

immutable __p2quire_iter
  sgn::Bool
  exp::Int64
  frc::UInt64
  zero::Bool
end

__p2quire_iter{N,ES}(p::Posit{N,ES}) = __p2quire_iter(posit_components(p)..., p == zero(Posit{N,ES}))

#return the starting state, which is just index 1.
Base.start(pi::__p2quire_iter) = 1
#the quire iterator is finished when the index exceeds the array length.
Base.done(pi::__p2quire_iter, index::Int64) = (index > 64)
Base.length(pi::__p2quire_iter) = 64

#takes the index, and the iterator object, and returns the corresponding cell.
function Base.next(pi::__p2quire_iter, index::Int64)

  pi.zero && return (zero(UInt64), index + 1)

  #return the current item and the state.
  if pi.exp < (64 * (index - 33))

    #           34              33                32              31
    #            1               0                -1              -2
    #   |________________|________________|________________|________________|
    #    127           64 63             0 -1           -64 -65         -128
    #                     1XXXXXXXXXXXXXXX X
    #
    #   63 -> 1           64 -> 2

    (pi.sgn ? -one(UInt64) : zero(UInt64), index + 1)

  elseif pi.exp > (64 * (index - 31) - 1)

    #           34              33                32              31
    #            1               0                -1              -2
    #   |________________|________________|________________|________________|
    #    127           64 63             0 -1           -64 -65         -128
    #                   1 XXXXXXXXXXXXXXXX
    #
    #   -1 -> 64          -2 -> 63

    (zero(UInt64), index + 1)

  elseif pi.exp < (64 * (index - 32))  # top half case.
    #           34              33                32              31
    #            1               0                -1              -2
    #   |________________|________________|________________|________________|
    #    127           64 63             0 -1           -64 -65         -128
    #                     1XXXXXXXXXXXXXXX X
    #

    # first determine the shift.  Use this value to append the front value

    shift = (reinterpret(UInt64, pi.exp) & 0x3F)

    prefix           = (pi.sgn ? -one(UInt64) : one(UInt64)) << shift
    shifted_fraction = pi.frc >> (64 - shift + pi.sgn)

    (prefix | shifted_fraction, index + 1)

  elseif pi.exp == (64 * (index - 32)) #the strange case

    (pi.sgn ? pi.frc >> 1 : pi.frc, index + 1)

  else

    shift = reinterpret(UInt64, pi.exp) & 0x3F - pi.sgn

    (pi.frc << shift, index + 1)
  end
end

function positvals_from(q::Quire)
  #look at the top fixed point to get important info.
  msword = last(q.fixed_point_value)

  #determine the sign.
  sign = (msword & 0x8000_0000_0000_0000) != 0
  exponent = 2047 + sign
  fractionbits = 0
  fraction = zero(UInt64)
  summary = false

  for idx = 64:-1:1
    current_cell = q.fixed_point_value[idx]

    if fractionbits == 0
      #we haven't seen an exponent yet.
      exp_delta = (sign ? leading_ones(current_cell) : leading_zeros(current_cell))
      exponent -= exp_delta
      fractionbits = 64 - exp_delta
      fraction = current_cell << (exp_delta + 1)
    elseif fractionbits == -1
      #special case where we can transfer the entire fraction over.
      fraction = current_cell
      fractionbits = 64
    elseif fractionbits == 64
      #case where we've seen everything
      (current_cell != zero(UInt64)) && (return (sign, exponent, fraction, true))
    else  #we have an incomplete fraction.  This cycle WILL complete it.
      fraction |= current_cell >> (fractionbits - 1)
      current_cell = current_cell & ((one(UInt64) << (fractionbits - 1)) - one(UInt64))
      fractionbits = 64
      (current_cell != zero(UInt64)) && (return (sign, exponent, fraction, true))
    end
  end
  (sign, exponent, fraction, false)
end

function (::Type{Posit{N,ES}}){N,ES}(q::Quire)
  #exceptional value handling:
  isinf(q) && return Posit{N,ES}(Inf)

  (sgn,exp,frc,summ) = positvals_from(q)

  #check for zero
  (exp < -2048) && return zero(Posit{N,ES})

  Posit{N,ES}(sgn,exp,frc)  #for now.  We will have to manage inexacts later.

end

doc"""
  posit_components breaks up a posit into components (sign, exp, frac)

  frac is in 2's complement.  This representation corresponds to: (-1^sign) * 2^exp + (2^exp * (1.frac))
"""
function posit_components{N,ES}(x::Posit{N,ES})
  sign = (@u(x) & @signbit) != 0
  inverted = (@u(x) & @invertbit) == 0

  shift = max(leading_zeros(@u(x) & (~@signbit)),leading_ones(@u(x) | @signbit))
  regime = shift - 2 #(sign ? 2 : 3)
  regime = (sign ⊻ inverted) ? -(regime + 1) : regime

  fraction = @u(x) << (shift + 1)

  #next extract the esize from the remainder of the fraction
  exp = ((-((@signbit) >> ES)) & fraction) >> (64 - ES)

  fraction = fraction << ES

  exponent = (regime + sign) << ES + (sign ? -@s(exp) : @s(exp))

  (sign, exponent, fraction)
end

################################################################################
## QUIRE FUNCTIONS

set!{T}(q::fquire{T}, x::T) = (q.q = x; x)

function set!{N,ES}(acc::Quire, x::Posit{N,ES})
  #shortcut evaluation of infinity.
  isinf(x) && return inf!(acc)
  acc.infinity = false

  for (idx, cell) in zip(1:64, __p2quire_iter(x))
    acc.fixed_point_value[idx] = cell
  end
  x
end

add!{T}(q::fquire{T}, x::T) = (q.q += x; q.q)
function add!{N,ES}(acc::Quire, x::Posit{N,ES})

  if isinf(x)
    isinf(acc) && throw(NaNError)
    return inf!(acc)
  end

  carry = false
  for (idx, cell) in zip(1:64, __p2quire_iter(x))
    old = acc.fixed_point_value[idx]
    acc.fixed_point_value[idx] += cell + carry
    carry = (acc.fixed_point_value[idx] <= old) & ((cell != zero(UInt64)) | (carry))
  end

  #return the cumulative sum.
  Posit{N,ES}(acc)
end

fdp!{T}(q::fquire{T}, a::T, b::T) = (q.q = fma(a, b, q.q); q.q)
function fdp!{N,ES}(acc::Quire, a::Posit{N,ES}, b::Posit{N,ES})

  #intercept exceptional cases.
  if isinf(a)
    (reinterpret(UInt64, a) == zero(UInt64)) && throw(NaNError(fdp!, a, b))
    isinf(acc) && throw(NaNError(fdp!, acc, [a, b]))
    inf!(acc)
    return Posit{N,ES}(Inf)
  end

  if isinf(b)
    (reinterpret(UInt64, a) == zero(UInt64)) && throw(NaNError(fdp!, [a, b]))
    isinf(acc) && throw(NaNError(fdp!, acc, [a, b]))
    inf!(acc)
    return Posit{N,ES}(Inf)
  end

  #intercept zero times anything, without changing the quire value.
  (reinterpret(UInt64, a) == zero(UInt64)) && return Posit{N,ES}(acc)
  (reinterpret(UInt64, b) == zero(UInt64)) && return Posit{N,ES}(acc)

  #the fused dot product will be calculated manually here.  For now, we'll only
  #support 32 bit x 32 bit fdp.

  (a_sgn, a_exp, a_frc) = posit_components(a)
  (b_sgn, b_exp, b_frc) = posit_components(b)

  #first, take the a_frc and b_frc and shift them right by 32 bits.
  a_frc_shft = a_frc >> 32
  b_frc_shft = b_frc >> 32

  #next actually do the multiplication.
  old_fraction = fraction = (a_frc_shft * b_frc_shft)#as a 64-bit unsigned integer.

  #the carry can be up to negative 2.
  mulcarry = 1 * (a_sgn ? -2 : 1) * (b_sgn ? -2 : 1)

  #next, add in the a_frc value
  if b_sgn
    #in the case that b is negative, invert and right shift the fraction addend.
    deltacarry = ((a_frc & 0x8000_0000_0000_0000) != 0) ? 1 : 0
    fraction = fraction - (a_frc << 1)
    mulcarry -= (fraction > old_fraction)
    mulcarry -= ((a_frc & 0x8000_0000_0000_0000) != 0)
  else
    fraction = fraction + a_frc
    mulcarry += (fraction < old_fraction)
  end

  old_fraction = fraction

  #then, add the b_frc value
  if a_sgn
    deltacarry = ((b_frc & 0x8000_0000_0000_0000) != 0) ? 1 : 0
    fraction = fraction - (b_frc << 1)
    mulcarry -= (fraction > old_fraction)
    mulcarry -= ((b_frc & 0x8000_0000_0000_0000) != 0)
  else
    fraction = fraction + b_frc
    mulcarry += (fraction < old_fraction)
  end

  #exponents in multiplication are additive.
  exponent = a_exp + b_exp
  #finally, look at the carry and shift the product fraction accordingly.

  sign = a_sgn ⊻ b_sgn

  (fraction, exponent) =
  if     (mulcarry == -4)
    (fraction >> 1, exponent + 1)
  elseif (mulcarry == -3)
    (fraction >> 1 | 0x8000_0000_0000_0000, exponent + 1)
  elseif (mulcarry == -2)
    (fraction, exponent)
  elseif (mulcarry == -1)
    (fraction << 1, exponent - 1)
  elseif (mulcarry == 0)
    (fraction << 1, exponent - 1)

  #one, two and three can be accessible both from the product of two positives and
  #the product of two negatives.
  elseif (mulcarry == 1)
    (fraction, exponent - (2 * a_sgn))
  elseif (mulcarry == 2)
    (fraction >> 1, exponent + (a_sgn ? -1 : 1))
  elseif (mulcarry == 3)
    ((fraction >> 1) | 0x8000_0000_0000_0000 , exponent + (a_sgn ? -1 : 1))

  #four can only be attained by two negatives.
  elseif (mulcarry == 4)
    (fraction >> 2, exponent)
  end

  #set up a quire iterator.
  qi = __p2quire_iter(sign, exponent, fraction, false)

  carry = false
  for (idx, cell) in zip(1:64, qi)
    old = acc.fixed_point_value[idx]
    acc.fixed_point_value[idx] += cell + carry
    carry = (acc.fixed_point_value[idx] <= old) & ((cell != zero(UInt64)) | (carry))
  end

  Posit{N,ES}(acc)
end

################################################################################
## some utility functions.

function find_lsb(q::Quire)
  digits_so_far = -2048
  for idx in 1:64
    lbits = trailing_zeros(q.fixed_point_value[idx])
    digits_so_far += lbits
    (lbits != 64) && break
  end
  digits_so_far
end

function fdp{N,ES}(v1::Vector{Posit{N,ES}}, v2::Vector{Posit{N,ES}}, q = Quire(Posit{N,ES}))
  zero!(q)
  res = zero(Posit{N,ES})
  for (x,y) in zip(v1, v2)
    res = fdp!(q, x, y)
  end
  res
end

function exact_sum{N,ES}(v1::Vector{Posit{N,ES}}, q = Quire(Posit{N,ES}))
  res = zero(Posit{N,ES})
  for x in v1
    res = add!(q, x)
  end
  res
end

function exact_mmult{N,ES}(M::Matrix{Posit{N,ES}}, v::Vector{Posit{N,ES}}, q = Quire(Posit{N,ES}))
  [fdp(M[idx, :], v, q) for idx in 1:size(M,1)]
end

function exact_mmult{N,ES}(M1::Matrix{Posit{N,ES}}, M2::Matrix{Posit{N,ES}}, q = Quire(Posit{N,ES}))
  [fdp(M1[idx, :], M2[:, jdx], q) for idx in 1:size(M1,1), jdx in 1:size(M2,2)]
end

export Quire, set!, add!, fdp!, fdp, exact_sum, exact_mmult

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAppsV1StatefulSetSpec <: SwaggerModel
    podManagementPolicy::Any # spec type: Union{ Nothing, String } # spec name: podManagementPolicy
    replicas::Any # spec type: Union{ Nothing, Int32 } # spec name: replicas
    revisionHistoryLimit::Any # spec type: Union{ Nothing, Int32 } # spec name: revisionHistoryLimit
    selector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: selector
    serviceName::Any # spec type: Union{ Nothing, String } # spec name: serviceName
    template::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PodTemplateSpec } # spec name: template
    updateStrategy::Any # spec type: Union{ Nothing, IoK8sApiAppsV1StatefulSetUpdateStrategy } # spec name: updateStrategy
    volumeClaimTemplates::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1PersistentVolumeClaim} } # spec name: volumeClaimTemplates

    function IoK8sApiAppsV1StatefulSetSpec(;podManagementPolicy=nothing, replicas=nothing, revisionHistoryLimit=nothing, selector=nothing, serviceName=nothing, template=nothing, updateStrategy=nothing, volumeClaimTemplates=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1StatefulSetSpec, Symbol("podManagementPolicy"), podManagementPolicy)
        setfield!(o, Symbol("podManagementPolicy"), podManagementPolicy)
        validate_property(IoK8sApiAppsV1StatefulSetSpec, Symbol("replicas"), replicas)
        setfield!(o, Symbol("replicas"), replicas)
        validate_property(IoK8sApiAppsV1StatefulSetSpec, Symbol("revisionHistoryLimit"), revisionHistoryLimit)
        setfield!(o, Symbol("revisionHistoryLimit"), revisionHistoryLimit)
        validate_property(IoK8sApiAppsV1StatefulSetSpec, Symbol("selector"), selector)
        setfield!(o, Symbol("selector"), selector)
        validate_property(IoK8sApiAppsV1StatefulSetSpec, Symbol("serviceName"), serviceName)
        setfield!(o, Symbol("serviceName"), serviceName)
        validate_property(IoK8sApiAppsV1StatefulSetSpec, Symbol("template"), template)
        setfield!(o, Symbol("template"), template)
        validate_property(IoK8sApiAppsV1StatefulSetSpec, Symbol("updateStrategy"), updateStrategy)
        setfield!(o, Symbol("updateStrategy"), updateStrategy)
        validate_property(IoK8sApiAppsV1StatefulSetSpec, Symbol("volumeClaimTemplates"), volumeClaimTemplates)
        setfield!(o, Symbol("volumeClaimTemplates"), volumeClaimTemplates)
        o
    end
end # type IoK8sApiAppsV1StatefulSetSpec

const _property_map_IoK8sApiAppsV1StatefulSetSpec = Dict{Symbol,Symbol}(Symbol("podManagementPolicy")=>Symbol("podManagementPolicy"), Symbol("replicas")=>Symbol("replicas"), Symbol("revisionHistoryLimit")=>Symbol("revisionHistoryLimit"), Symbol("selector")=>Symbol("selector"), Symbol("serviceName")=>Symbol("serviceName"), Symbol("template")=>Symbol("template"), Symbol("updateStrategy")=>Symbol("updateStrategy"), Symbol("volumeClaimTemplates")=>Symbol("volumeClaimTemplates"))
const _property_types_IoK8sApiAppsV1StatefulSetSpec = Dict{Symbol,String}(Symbol("podManagementPolicy")=>"String", Symbol("replicas")=>"Int32", Symbol("revisionHistoryLimit")=>"Int32", Symbol("selector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("serviceName")=>"String", Symbol("template")=>"IoK8sApiCoreV1PodTemplateSpec", Symbol("updateStrategy")=>"IoK8sApiAppsV1StatefulSetUpdateStrategy", Symbol("volumeClaimTemplates")=>"Vector{IoK8sApiCoreV1PersistentVolumeClaim}")
Base.propertynames(::Type{ IoK8sApiAppsV1StatefulSetSpec }) = collect(keys(_property_map_IoK8sApiAppsV1StatefulSetSpec))
Swagger.property_type(::Type{ IoK8sApiAppsV1StatefulSetSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1StatefulSetSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1StatefulSetSpec }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1StatefulSetSpec[property_name]

function check_required(o::IoK8sApiAppsV1StatefulSetSpec)
    (getproperty(o, Symbol("selector")) === nothing) && (return false)
    (getproperty(o, Symbol("serviceName")) === nothing) && (return false)
    (getproperty(o, Symbol("template")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAppsV1StatefulSetSpec }, name::Symbol, val)
end

ext: jl


text:
macro log(exs...)
    blk = Expr(:block)
    loc = string("\t", __source__.file, "#", __source__.line)
    for ex in exs
        push!(blk.args, :(println(getsocket(),
                                  $(sprint(show_unquoted,ex)*" = "),
                                  repr(begin value=$(esc(ex)) end),
                                  "\n",
                                  $loc,
                                  "\n",
                                  )))
    end
    isempty(exs) || push!(blk.args, :value)
    return blk
end

# macro log(str)
#     blk = Expr(:block)
#     loc = string("\t", __source__.file, "#", __source__.line)
#     push!(blk.args, :(println(getsocket(),
#                               # repr(begin value=$(esc(str)) end),
#                               $str,
#                               "\n",
#                               $loc,
#                               "\n",
#                               )))
#     push!(blk.args, :value)
#     return blk
# end

function log(s...)
    println(getsocket(), s...)
end

function test_bind()

end

# 'j' => (s::LE.MIState, o...)->LE.edit_move_down(s),

macro bindkey(c)
    # char = ($(esc(c)))
    # char = esc(c)
    # @show sym = Symbol(char)
    return :(
        ()->eval(Expr(:call, Symbol($(esc(c))))))
         # pair item 2
         # (s::LE.MIState, o...)->begin
         # @show Symbol($(esc(c)))
         # Expr(:call, Symbol($(esc(c))),(VB.mode, esc(s)))
end

# buffer(s :: LE.MIState) = LE.buffer(s)

# alpha_keymap = AnyDict()


text:
export continuity

Space(d::IntervalOrSegment) = Chebyshev(d)
Space(d::FullSpace{<:Real}) = Chebyshev(Line())

Fun(::typeof(identity), d::IntervalOrSegment{T}) where {T<:Number} =
    Fun(Chebyshev(d), [mean(d), complexlength(d)/2])

## Calculus

# the default domain space is higher to avoid negative ultraspherical spaces
Integral(d::IntervalOrSegment,n::Integer) = Integral(Ultraspherical(1,d),n)


text:
#=
test_parser:
- Julia version: 1.0.3
- Author: rhdekker
- Date: 2019-02-17
=#
using Test

@testset "parser" begin
    # using HyperGraphTools.PartiallyOrderedParser
    using HyperGraphTools.SimpleTAGMLTokenizer
    using HyperGraphTools.PartiallyOrderedParser

    r = my_parse()
    @test r == "Hello World"

    function has_no_errors(x)
        isempty(filter(x -> x.type == SimpleTAGMLTokenizer.ERROR, x))
    end

    # use case 1: 2 layers, full overlap
    tagml1 = "[root>[tag|a>[other|b>one two three<other]<tag]<root]"
    # grammar1 = """
    # S -> ROOT
    # ROOT -> [root> TAG <root]
    # TAG -> [tag|a> OTHER <tag]
    # OTHER -> [other|b> TEXT <other]
    # TEXT -> "one two three"
    # """
    t1 = tokenize(tagml1)
    @test has_no_errors(t1)
end

max_stars_repo_name: rhdekker/HyperGraphTools

max_issues_repo_name: rhdekker/HyperGraphTools

max_forks_repo_name: rhdekker/HyperGraphTools


text:
using URIs

export uri_add_path

"""
    uri_add_path(::URI, ::String...)
Return a new `URI` with the original path extended by the
specified additional components.
Also see URIs.resolvereference, which only accepts a single
additional path component.
"""
function uri_add_path(uri::URI, subs::String...)::URI
    path = ["", URIs.splitpath(uri)..., subs...]
    URI(uri; path=join(path, "/"))
end

function uri_add_path(uri::AbstractString, subs::String...)::URI
    uri_add_path(URI(uri), subs...)
end


text:
import Base.length
import Base.start
import Base.next
import Base.done

type StatefulIterator{T}
    iter::T
    state
end
"Captures the state of the iterator so the user does not need to."
statefuliterator(iter) = StatefulIterator(iter, start(iter))
function Base.next(si::StatefulIterator)
    item,state = next(si.iter, si.state)
    si.state = state
    item
end
Base.done(si::StatefulIterator) = done(si.iter, si.state)
Base.length(si::StatefulIterator) = length(si.iter)

"Distributes the given number of samples proportionally."
function allocate_samples(proportions, nsamples, maxCounts)
    counts = round(Int, nsamples*proportions/sum(proportions))
    counts[counts .> maxCounts] = maxCounts[counts .> maxCounts]
    total = sum(counts)
    changed = true
    while changed
        changed = false
        for ind in shuffle(find(counts .!= maxCounts))
            total != nsamples || break

            if total < nsamples && counts[ind] < maxCounts[ind]
                counts[ind] += 1
                total += 1
                changed = true
            elseif (total > nsamples && counts[ind] > 0) || (counts[ind] > maxCounts[ind])
                counts[ind] -= 1
                total -= 1
                changed = true
            end
        end
    end
    counts
end

# http://www.nowozin.net/sebastian/blog/streaming-mean-and-variance-computation.html
type MeanVarianceAccumulator
    sumw::Float64
    wmean::Float64
    t::Float64
    n::Int

    function MeanVarianceAccumulator()
        new(0.0, 0.0, 0.0, 0)
    end
end
function observe!(mvar::MeanVarianceAccumulator, value, weight)
    @assert weight >= 0.0
    q = value - mvar.wmean
    temp_sumw = mvar.sumw + weight
    r = q*weight / temp_sumw

    mvar.wmean += r
    mvar.t += q*r*mvar.sumw
    mvar.sumw = temp_sumw
    mvar.n += 1

    nothing
end
Base.length(mvar::MeanVarianceAccumulator) = mvar.n
Base.mean(mvar::MeanVarianceAccumulator) = mvar.wmean
Base.var(mvar::MeanVarianceAccumulator) = (mvar.t*mvar.n)/(mvar.sumw*(mvar.n-1))
Base.std(mvar::MeanVarianceAccumulator) = sqrt(var(mvar))

"Identify which feature groups vary."
function varying_groups(x, X, featureGroups)
    varying = zeros(length(featureGroups))
    for (i,inds) in enumerate(featureGroups)
        varying[i] = sum(vec(sum(x[inds] .== X[inds,:],1) .!= length(inds)))
    end
    find(varying)
end


text:
### Julia OpenStreetMap Package ###
### MIT License                 ###
### Copyright 2014              ###

type Highway
    @compat class::AbstractString       # Type of highway
    lanes::Int          # Number of lanes (1 if unspecified)
    oneway::Bool        # True if road is one-way
    @compat sidewalk::AbstractString    # Sidewalk classifier, if available
    @compat cycleway::AbstractString    # Cycleway classifier, if available
    @compat bicycle::AbstractString     # Bicycle classifier, if available
    @compat name::AbstractString        # Name, if available
    maxspeed::Float64       # Speed limit in kph
    nodes::Vector{Int}  # List of nodes
end

type Segment
    node0::Int          # Source node ID
    node1::Int          # Target node ID
    nodes::Vector{Int}  # List of nodes falling within node0 and node1
    dist::Real          # Length of the segment
    class::Int          # Class of the segment
    parent::Int         # ID of parent highway
    oneway::Bool        # True if road is one-way
end

type Feature
    @compat class::AbstractString       # Shop, amenity, crossing, etc.
    @compat detail::AbstractString      # Class qualifier
    @compat name::AbstractString        # Name
end

type Building
    @compat class::AbstractString       # Building type (usually "yes")
    @compat name::AbstractString        # Building name (usually unavailable)
    nodes::Vector{Int}  # List of nodes
end

type Intersection
    highways::Set{Int}  # Set of highway IDs
end
Intersection() = Intersection(Set{Int}())

type HighwaySet # Multiple highways representing a single "street" with a common name
    highways::Set{Int}
end

# Transporation network graph data and helpers to increase routing speed
type Network
    g                                   # Graph object
    v::Dict{Int,Graphs.KeyVertex{Int}}  # (node id) => (graph vertex)
    w::Vector{Float64}                  # Edge weights, indexed by edge id
    class::Vector{Int}                 # Road class of each edge
end

### Rendering style data
type Style
    @compat color::UInt32
    width::Real
    @compat spec::AbstractString
end
Style(x, y) = Style(x, y, "-")


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using RadMod2D

include("./cases_epseff.jl")

# epsilon_effective_round_groove()
# epsilon_effective_v_groove()
# epsilon_effective_square_groove()
# epsilon_effective_elbow_piece()
# epsilon_effective_trapez()
epsilon_effective_cosinus()

# plot_epseff_round_groove()
# plot_epseff_v_groove()
# plot_epseff_square_groove()
# plot_epseff_elbow_piece()
# plot_epseff_trapez()
plot_epseff_cosinus()

# include("./test/run_epseff.jl")


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

##########################################################################################################################################
# The reserves.jl module contains functions to creates decision variables related to frequency regulation and reserves provision
# and constraints setting overall system requirements for regulation and operating reserves.
##########################################################################################################################################

@doc raw"""
	reserves(EP::Model, inputs::Dict, UCommit::Int)

This function sets up reserve decisions and constraints, using the reserves_core()` and reserves_contingency()` functions.
"""
function reserves(EP::Model, inputs::Dict, UCommit::Int)

	if inputs["pStatic_Contingency"] > 0 || (UCommit >= 1 && inputs["pDynamic_Contingency"] >= 1)
		EP = reserves_contingency(EP, inputs, UCommit)
	end

	EP = reserves_core(EP, inputs, UCommit)

	return EP
end

@doc raw"""
	reserves_contingency(EP::Model, inputs::Dict, UCommit::Int)

This function establishes several different versions of contingency reserve requirement expression, $CONTINGENCY$ used in the reserves_core() function below.

Contingency reserves represent requirements for upward ramping capability within a specified time frame to compensated for forced outages or unplanned failures of generators or transmission lines (e.g. N-1 contingencies).

There are three options for the $Contingency$ expression, depending on user settings:
	1. a static contingency, in which the contingency requirement is set based on a fixed value (in MW) specified in the '''reserves.csv''' input file;
	2. a dynamic contingency based on installed capacity decisions, in which the largest 'installed' generator is used to determine the contingency requirement for all time periods; and
	3. dynamic unit commitment based contingency, in which the largest 'committed' generator in any time period is used to determine the contingency requirement in that time period.

Note that the two dynamic contigencies are only available if unit commitment is being modeled.

**Static contingency**
Option 1 (static contingency) is expressed by the following constraint:
```math
\begin{aligned}
	Contingency = \epsilon^{contingency}
\end{aligned}
```
where $\epsilon^{contingency}$ is static contingency requirement in MWs.

**Dynamic capacity-based contingency**
Option 2 (dynamic capacity-based contingency) is expressed by the following constraints:
```math
\begin{aligned}
	&Contingency \geq \Omega^{size}_{y,z} \times \alpha^{Contingency,Aux}_{y,z} & \forall y \in \mathcal{UC}, z \in \mathcal{Z}\\
	&\alpha^{Contingency,Aux}_{y,z} \leq \Delta^{\text{total}}_{y,z} & \forall y \in \mathcal{UC}, z \in \mathcal{Z}\\
	&\alpha^{Contingency,Aux}_{y,z} \geq M_y \times \Delta^{\text{total}}_{y,z} & \forall y \in \mathcal{UC}, z \in \mathcal{Z}\\
\end{aligned}
```

where $M_y$ is a `big M' constant equal to the largest possible capacity that can be installed for generation cluster $y$, and $\alpha^{Contingency,Aux}_{y,z} \in [0,1]$ is a binary auxiliary variable that is forced by the second and third equations above to be 1 if the total installed capacity $\Delta^{\text{total}}_{y,z} > 0$ for any generator $y \in \mathcal{UC}$ and zone $z$, and can be 0 otherwise. Note that if the user specifies contingency option 2, and is also using the linear relaxation of unit commitment constraints, the capacity size parameter for units in the set $\mathcal{UC}$ must still be set to a discrete unit size for this contingency to work as intended.

**Dynamic commitment-based contingency**
Option 3 (dynamic commitment-based contingency) is expressed by the following set of constraints:
```math
\begin{aligned}
	& Contingency \geq \Omega^{size}_{y,z} \times Contingency\_Aux_{y,z,t} & \forall y \in \mathcal{UC}, z \in \mathcal{Z}\\
	& Contingency\_Aux_{y,z,t} \leq \nu_{y,z,t} & \forall y \in \mathcal{UC}, z \in \mathcal{Z}\\
	& Contingency\_Aux_{y,z,t} \geq M_y \times \nu_{y,z,t} & \forall y \in \mathcal{UC}, z \in \mathcal{Z}\\
\end{aligned}
```

where $M_y$ is a `big M' constant equal to the largest possible capacity that can be installed for generation cluster $y$, and $Contingency\_Aux_{y,z,t} \in [0,1]$ is a binary auxiliary variable that is forced by the second and third equations above to be 1 if the commitment state for that generation cluster $\nu_{y,z,t} > 0$ for any generator $y \in \mathcal{UC}$ and zone $z$ and time period $t$, and can be 0 otherwise. Note that this dynamic commitment-based contingency can only be specified if discrete unit commitment decisions are used (e.g. it will not work if relaxed unit commitment is used).
"""
function reserves_contingency(EP::Model, inputs::Dict, UCommit::Int)

	println("Reserves Contingency Module")

	dfGen = inputs["dfGen"]

	T = inputs["T"]     # Number of time steps (hours)
	COMMIT = inputs["COMMIT"]

	if UCommit >= 1
		pDynamic_Contingency = inputs["pDynamic_Contingency"]
	end

	### Variables ###

	# NOTE: If Dynamic_Contingency == 0, then contingency is a fixed parameter equal the value specified in reserves.csv via pStatic_Contingency.
	if UCommit == 1 && pDynamic_Contingency == 1
		# Contingency = largest installed thermal unit
		@variable(EP, vLARGEST_CONTINGENCY >= 0)
		# Auxiliary variable that is 0 if vCAP = 0, 1 otherwise
		@variable(EP, vCONTINGENCY_AUX[y in COMMIT], Bin)
	elseif UCommit == 1 && pDynamic_Contingency == 2
		# Contingency = largest committed thermal unit in each time period
		@variable(EP, vLARGEST_CONTINGENCY[t=1:T] >= 0)
		# Auxiliary variable that is 0 if vCOMMIT = 0, 1 otherwise
		@variable(EP, vCONTINGENCY_AUX[y in COMMIT, t=1:T], Bin)
	end

	### Expressions ###
	if UCommit == 1 && pDynamic_Contingency == 1
		# Largest contingency defined as largest installed generator
		println("Dynamic Contingency Type 1: Modeling the largest contingency as the largest installed generator")
		@expression(EP, eContingencyReq[t=1:T], vLARGEST_CONTINGENCY)
	elseif UCommit == 1 && pDynamic_Contingency == 2
		# Largest contingency defined for each hour as largest committed generator
		println("Dynamic Contingency Type 2: Modeling the largest contingency as the largest largest committed generator")
		@expression(EP, eContingencyReq[t=1:T], vLARGEST_CONTINGENCY[t])
	else
		# Largest contingency defined fixed as user-specifed static contingency in MW
		println("Static Contingency: Modeling the largest contingency as user-specifed static contingency")
		@expression(EP, eContingencyReq[t=1:T], inputs["pStatic_Contingency"])
	end

	### Constraints ###

	# Dynamic contingency related constraints
		# option 1: ensures vLARGEST_CONTINGENCY is greater than the capacity of the largest installed generator
	if UCommit == 1 && pDynamic_Contingency == 1
		@constraint(EP, cContingency[y in COMMIT], vLARGEST_CONTINGENCY >=
			inputs["dfGen"][!,:Cap_Size][y]*vCONTINGENCY_AUX[y] )
		# Ensure vCONTINGENCY_AUX = 0 if total capacity = 0
		@constraint(EP, cContAux1[y in COMMIT], vCONTINGENCY_AUX[y] <= EP[:eTotalCap][y])
		# Ensure vCONTINGENCY_AUX = 1 if total capacity > 0
		@constraint(EP, cContAux2[y in COMMIT], EP[:eTotalCap][y] <= inputs["pContingency_BigM"][y]*vCONTINGENCY_AUX[y])

		# option 2: ensures vLARGEST_CONTINGENCY is greater than the capacity of the largest commited generator in each hour
	elseif UCommit == 1 && pDynamic_Contingency == 2
		@constraint(EP, cContingency[y in COMMIT, t=1:T], vLARGEST_CONTINGENCY[t] >=
			inputs["dfGen"][!,:Cap_Size][y]*vCONTINGENCY_AUX[y,t] )
		# Ensure vCONTINGENCY_AUX = 0 if vCOMMIT = 0
		@constraint(EP, cContAux[y in COMMIT, t=1:T], vCONTINGENCY_AUX[y,t] <= EP[:vCOMMIT][y,t])
		# Ensure vCONTINGENCY_AUX = 1 if vCOMMIT > 0
		@constraint(EP, cContAux2[y in COMMIT, t=1:T], EP[:vCOMMIT][y, t] <= inputs["pContingency_BigM"][y]*vCONTINGENCY_AUX[y,t])
	end

	return EP
end

@doc raw"""
	reserves_core(EP::Model, inputs::Dict, UCommit::Int)

This function creates decision variables related to frequency regulation and reserves provision and constraints setting overall system requirements for regulation and operating reserves.

**Regulation and reserves decisions**
$f_{y,t,z} \geq 0$ is the contribution of generation or storage resource $y \in Y$ in time $t \in T$ and zone $z \in Z$ to frequency regulation

$r_{y,t,z} \geq 0$ is the contribution of generation or storage resource $y \in Y$ in time $t \in T$ and zone $z \in Z$ to operating reserves up

We assume frequency regulation is symmetric (provided in equal quantity towards both upwards and downwards regulation). To reduce computational complexity, operating reserves are only modeled in the upwards direction, as downwards reserves requirements are rarely binding in practice.

Storage resources ($y \in \mathcal{O}$) have two pairs of auxilary variables to reflect contributions to regulation and reserves when charging and discharging, where the primary variables ($f_{y,z,t}$ \& $r_{y,z,t}$) becomes equal to sum of these auxilary variables.

**Unmet operating reserves**

$unmet\_rsv_{t} \geq 0$ denotes any shortfall in provision of operating reserves during each time period $t \in T$

There is a penalty $C^{rsv}$ added to the objective function to penalize reserve shortfalls, equal to:

```math
\begin{aligned}
	C^{rvs} = \sum_{t \in T} \omega_t \times unmet\_rsv_{t}
\end{aligned}
```

**Frequency regulation requirements**

Total requirements for frequency regulation (aka primary reserves) in each time step $t$ are specified as fractions of hourly demand (to reflect demand forecast errors) and variable renewable avaialblity in the time step (to reflect wind and solar forecast errors).

```math
\begin{aligned}
	\sum_{y \in Y, z \in Z} f_{y,t,z} \geq \epsilon^{load}_{reg} \times \sum_{z \in Z} \mathcal{D}_{z,t} + \epsilon^{vre}_{reg} \times \sum_{z \in Z} \rho^{max}_{y,z,t} \times \Delta^{\text{total}}_{y,z} \quad \forall t \in T
\end{aligned}
```
where $\mathcal{D}_{z,t}$ is the forecasted electricity demand in zone $z$ at time $t$ (before any demand flexibility); $\rho^{max}_{y,z,t}$ is the forecasted capacity factor for variable renewable resource $y \in VRE$ and zone $z$ in time step $t$; $\Delta^{\text{total}}_{y,z}$ is the total installed capacity of variable renewable resources $y \in VRE$ and zone $z$; and $\epsilon^{load}_{reg}$ and $\epsilon^{vre}_{reg}$ are parameters specifying the required frequency regulation as a fraction of forecasted demand and variable renewable generation.

**Operating reserve requirements**

Total requirements for operating reserves in the upward direction (aka spinning reserves or contingency reserces or secondary reserves) in each time step $t$ are specified as fractions of time step's demand (to reflect demand forecast errors) and variable renewable avaialblity in the time step (to reflect wind and solar forecast errors) plus the largest planning contingency (e.g. potential forced generation outage).

```math
\begin{aligned}
	\sum_{y \in Y, z \in Z} r_{y,z,t} + r^{unmet}_{t} \geq \epsilon^{load}_{rsv} \times \sum_{z \in Z} \mathcal{D}_{z,t} + \epsilon^{vre}_{rsv} \times \sum_{z \in Z} \rho^{max}_{y,z,t} \times \Delta^{\text{total}}_{y,z} + Contingency \quad \forall t \in T
\end{aligned}
```

where $\mathcal{D}_{z,t}$ is the forecasted electricity demand in zone $z$ at time $t$ (before any demand flexibility); $\rho^{max}_{y,z,t}$ is the forecasted capacity factor for variable renewable resource $y \in VRE$ and zone $z$ in time step $t$; $\Delta^{\text{total}}_{y,z}$ is the total installed capacity of variable renewable resources $y \in VRE$ and zone $z$; and $\epsilon^{load}_{rsv}$ and $\epsilon^{vre}_{rsv}$ are parameters specifying the required contingency reserves as a fraction of forecasted demand and variable renewable generation. $Contingency$ is an expression defined in the reserves\_contingency() function meant to represent the largest `N-1` contingency (unplanned generator outage) that the system operator must carry operating reserves to cover and depends on how the user wishes to specify contingency requirements.
"""
function reserves_core(EP::Model, inputs::Dict, UCommit::Int)

	# DEV NOTE: After simplifying reserve changes are integrated/confirmed, should we revise such that reserves can be modeled without UC constraints on?
	# Is there a use case for economic dispatch constraints with reserves?

	println("Reserves Core Module")

	dfGen = inputs["dfGen"]

	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones

	REG = inputs["REG"]
	RSV = inputs["RSV"]

	### Variables ###

	## Integer Unit Commitment configuration for variables

	## Decision variables for reserves
	@variable(EP, vREG[y in REG, t=1:T] >= 0) # Contribution to regulation (primary reserves), assumed to be symmetric (up & down directions equal)
	@variable(EP, vRSV[y in RSV, t=1:T] >= 0) # Contribution to operating reserves (secondary reserves or contingency reserves); only model upward reserve requirements

	# Storage techs have two pairs of auxilary variables to reflect contributions to regulation and reserves
	# when charging and discharging (primary variable becomes equal to sum of these auxilary variables)
	@variable(EP, vREG_discharge[y in intersect(inputs["STOR_ALL"], REG), t=1:T] >= 0) # Contribution to regulation (primary reserves) (mirrored variable used for storage devices)
	@variable(EP, vRSV_discharge[y in intersect(inputs["STOR_ALL"], RSV), t=1:T] >= 0) # Contribution to operating reserves (secondary reserves) (mirrored variable used for storage devices)
	@variable(EP, vREG_charge[y in intersect(inputs["STOR_ALL"], REG), t=1:T] >= 0) # Contribution to regulation (primary reserves) (mirrored variable used for storage devices)
	@variable(EP, vRSV_charge[y in intersect(inputs["STOR_ALL"], RSV), t=1:T] >= 0) # Contribution to operating reserves (secondary reserves) (mirrored variable used for storage devices)

	@variable(EP, vUNMET_RSV[t=1:T] >= 0) # Unmet operating reserves penalty/cost

	### Expressions ###
	## Total system reserve expressions
	# Regulation requirements as a percentage of load and scheduled variable renewable energy production in each hour
	# Reg up and down requirements are symmetric
	@expression(EP, eRegReq[t=1:T], inputs["pReg_Req_Load"]*sum(inputs["pD"][t,z] for z=1:Z) +
		inputs["pReg_Req_VRE"]*sum(inputs["pP_Max"][y,t]*EP[:eTotalCap][y] for y in intersect(inputs["VRE"], inputs["MUST_RUN"])))
	# Operating reserve up / contingency reserve requirements as ˚a percentage of load and scheduled variable renewable energy production in each hour
	# and the largest single contingency (generator or transmission line outage)
	@expression(EP, eRsvReq[t=1:T], inputs["pRsv_Req_Load"]*sum(inputs["pD"][t,z] for z=1:Z) +
				inputs["pRsv_Req_VRE"]*sum(inputs["pP_Max"][y,t]*EP[:eTotalCap][y] for y in intersect(inputs["VRE"], inputs["MUST_RUN"])))

	# N-1 contingency requirement is considered only if Unit Commitment is being modeled
	if UCommit >= 1 && (inputs["pDynamic_Contingency"] >= 1 || inputs["pStatic_Contingency"] > 0)
		EP[:eRsvReq] = EP[:eRsvReq] + EP[:eContingencyReq]
	end

	## Objective Function Expressions ##

	# Penalty for unmet operating reserves
	@expression(EP, eCRsvPen[t=1:T], inputs["omega"][t]*inputs["pC_Rsv_Penalty"]*vUNMET_RSV[t])
	@expression(EP, eTotalCRsvPen, sum(eCRsvPen[t] for t=1:T) +
		sum(dfGen[!,:Reg_Cost][y]*vRSV[y,t] for y in RSV, t=1:T) +
		sum(dfGen[!,:Rsv_Cost][y]*vREG[y,t] for y in REG, t=1:T) )
	EP[:eObj] += eTotalCRsvPen

	### Constraints ###

	## Total system reserve constraints
	# Regulation requirements as a percentage of load and scheduled variable renewable energy production in each hour
	# Note: frequencty regulation up and down requirements are symmetric and all resources contributing to regulation are assumed to contribute equal capacity to both up and down directions
	@constraint(EP, cReg[t=1:T], sum(vREG[y,t] for y in REG) >= EP[:eRegReq][t])

	@constraint(EP, cRsvReq[t=1:T], sum(vRSV[y,t] for y in RSV) + vUNMET_RSV[t] >= EP[:eRsvReq][t])

	return EP
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
include("forecast_types.jl")
include("helper_functions.jl")

function calculate_bounds(local_cases, τ, state)

    # observation period
    T = length(local_cases)

    idxs_limits =  Vector{UnitRange{Int64}}()

    idxs_start = 1
    while idxs_start - 1 != length(local_cases)
        idxs_end = idxs_start + τ - 1
        if length(local_cases) - idxs_end < τ
            idxs_end = length(local_cases)
        end
        push!(idxs_limits, idxs_start:idxs_end)
        idxs_start = idxs_end + 1
    end

    Cₜ = [sum(local_cases[idx]) for idx in idxs_limits]

    # multipliers on the n-day average
    # (ℓ, u) = (0.25, 4.0)
    (ℓ, u) = (0.1, 4.0)
    Lₜ = ceil.(Int, ℓ * Cₜ)
    Uₜ = ceil.(Int, u * Cₜ)

    # remove restrictions over last τ * 2 days
    (ℓ, u) = (0.1, 1.5)
    Lₜ[(end - 4):end] = ceil.(Int, ℓ * Cₜ[(end - 4):end])
    Uₜ[(end - 4):end] = ceil.(Int, u * Cₜ[(end - 4):end])

    # remove restrictions over last τ * 2 days
    (ℓ, u) = (0.25, 10.0)
    Lₜ[end] = ceil.(Int, ℓ * Cₜ[end])
    Uₜ[end] = ceil.(Int, u * Cₜ[end])

    # (ℓ, u) = (0.7, 2.0)
    # Lₜ[end-1:end] = ceil.(Int, ℓ * Cₜ[end-1:end])
    # Uₜ[end-1:end] = ceil.(Int, u * Cₜ[end-1:end])

    min_limit = τ * 50

    Lₜ[Lₜ .< min_limit] .= 0
    Uₜ[Uₜ .< min_limit] .= min_limit

    return (Lₜ, Uₜ, idxs_limits)

end

function count_cases!(
    case_counts,
    forecast::Forecast,
)

    idxs_limits = forecast.sim_features.idxs_limits
    D = forecast.sim_realisation.D
    case_counts .= 0

    for (i, idx) in enumerate(idxs_limits)
        case_counts[i] = sum(@view D[idx, 1:2])
    end

    return nothing

end

function get_simulation_limits(
    local_cases,
    forecast_start_date,
    TP_indices,
    N,
    state;
    τ = 5,
    data_truncation = 7,
    nowcast_truncation = 30,
    fitting_truncation = 14,
)
    """
    Using the observed cases, determine the limits of cases over the backcast and
    nowcast. This assumes consistency over windows of fixed length and a nowcast period
    of 14 days.
    """

    # length of observed time series
    T_observed = length(local_cases)
    # duration of forecast simulation including initial day
    T_end = sum(ind >= 0 for ind in TP_indices)

    (min_cases, max_cases, idxs_limits) = calculate_bounds(local_cases, τ, state)

    # the maximum allowable cases over the forecast period is the population size
    max_forecast_cases = N ÷ 2

    omicron_start_date = "2021-11-15"
    omicron_only_date = "2022-02-01"
    omicron_start_day = Dates.value(
        Dates.Date(omicron_start_date) - Dates.Date(forecast_start_date)
    )
    omicron_only_day = Dates.value(
        Dates.Date(omicron_only_date) - Dates.Date(forecast_start_date)
    )

    sim_features = Features(
        max_forecast_cases,
        N,
        T_observed,
        T_end,
        τ,
        min_cases,
        max_cases,
        idxs_limits,
        state,
        omicron_start_day,
        omicron_only_day,
    )

    return sim_features

end

function check_sim!(
    forecast::Forecast,
    case_counts,
    local_cases,
    reinitialise_allowed;
    day = 0,
    sim = 0,
)
    """
    Check for consistency of the simulations against the data. This will also
    check for instances of superspreading events and will insert cases if conditions
    are met.
    """

    # guard clause for if we shouldn't check the simulation
    day < 0 && return (false, false)

    print_status = false
    print_inject_status = false

    Z = forecast.sim_realisation.Z
    D = forecast.sim_realisation.D
    U = forecast.sim_realisation.U

    # days forecast observed for
    T_observed = forecast.sim_features.T_observed
    T_end = forecast.sim_features.T_end
    τ = forecast.sim_features.τ
    max_forecast_cases = forecast.sim_features.max_forecast_cases
    consistency_multiplier = forecast.sim_constants.consistency_multiplier

    min_cases = forecast.sim_features.min_cases
    max_cases = forecast.sim_features.max_cases
    idxs_limits = forecast.sim_features.idxs_limits

    # initialise as bad sim
    bad_sim = false
    injected_cases = false

    # if day == 0
    #     # count how many cases each day 0
    #     count_cases!(case_counts, forecast)
    #     if any(case_counts .< min_cases)
    #         bad_sim = true
    #         return (bad_sim, injected_cases)
    #     end
    # end

    # only check for consistency if we don't inject cases.
    ζ = ceil(map_day_to_index_cases(day) ÷ τ)

    if 1 <= ζ <= length(case_counts)
        # count how many cases each day 0
        count_cases!(case_counts, forecast)
        if any(case_counts[1:ζ] .> max_cases[1:ζ])
            if print_status
                bad_id = case_counts[1:ζ] .> max_cases[1:ζ]
                println(
                    "Sim: ", sim, "\n",
                    "Too many cases: ", case_counts[1:ζ][bad_id], "\n",
                    "Maximum: ", max_cases[1:ζ][bad_id], "\n",
                    "Window: ", findfirst(bad_id), "\n",
                )
            end
            bad_sim = true
        end

        if !bad_sim && reinitialise_allowed
            # calculate the number of detections missed over the τ day period
            X_day_range = min(T_observed - 1, day - τ + 1):min(T_observed - 1, day)
            actual_X_day_cases = sum(
                @view local_cases[map_day_to_index_cases.(X_day_range)]
            )
            # calculate total number of cases over week
            sim_X_day_cases = sum(@view D[map_day_to_index_UD.(X_day_range), 1:2])
            missing_detections = 0
            threshold = ceil(Int, consistency_multiplier * max(1, sim_X_day_cases))

            if (actual_X_day_cases > threshold)
                if print_status && print_inject_status
                    println(
                        "======\n",
                        "Sim: ", sim, "\n",
                        "Injected cases...\n",
                        "Actual cases: ", actual_X_day_cases, "\n",
                        "Sim cases: ", sim_X_day_cases, "\n",
                        "Day added: ", day, "\n",
                        "======",
                    )
                end

                # uniformly sample a number of missing detections to add in
                missing_detections = ceil(
                    Int,
                    actual_X_day_cases - sim_X_day_cases,
                )

                injected_cases = true

                inject_cases!(
                    forecast::Forecast,
                    missing_detections,
                    X_day_range,
                )
            end
        end

        # count updated cases
        count_cases!(case_counts, forecast)

        if any(case_counts[1:ζ] .< min_cases[1:ζ])
            if print_status
                bad_id = case_counts[1:ζ] .< min_cases[1:ζ]
                println(
                    "Sim: ", sim, "\n",
                    "Too few cases: ", case_counts[1:ζ][bad_id], "\n",
                    "Minimum:  ", min_cases[1:ζ][bad_id], "\n",
                    "Window: ", findfirst(bad_id), "\n",
                )
            end
            bad_sim = true
        end

    elseif ζ > length(case_counts)
        # this is just the total cases over the forecast horizon
        D_forecast = sum(@view D[T_observed + 1:end, 1:2])
        # if we've exceeded the max cases over a given period
        if D_forecast > max_forecast_cases
            bad_sim = true
        end
    end

    return (bad_sim, injected_cases)

end

function inject_cases!(
    forecast::Forecast,
    missing_detections,
    X_day_range,
)
    """
    A model for injecting cases into the simulation following long periods of low
    incidence.
    """

    Z = forecast.sim_realisation.Z
    D = forecast.sim_realisation.D
    U = forecast.sim_realisation.U
    individual_type_map = forecast.individual_type_map

    total_injected = 0
    # assume that by mid December the adjusted GI and IP are noticeable for the case insertion.
    # small approximation but should be Ok overall.
    omicron_dominant_day = forecast.sim_features.omicron_start_day + 30
    # initialise the parameters
    p_symp = 0.0
    p_detect_given_symp = 0.0
    p_detect_given_asymp = 0.0
    p_symp_given_detect = 0.0

    ind = map_day_to_index_p(X_day_range[end])

    p_symp = forecast.sim_constants.p_symp[ind]
    p_detect_given_symp = forecast.sim_constants.p_detect_given_symp[ind]
    p_detect_given_asymp = forecast.sim_constants.p_detect_given_asymp[ind]
    p_symp_given_detect = forecast.sim_constants.p_symp_given_detect[ind]

    # sample number detected
    num_symptomatic_detected = sample_binomial_limit(missing_detections, p_symp_given_detect)
    num_asymptomatic_detected = missing_detections - num_symptomatic_detected

    # infer some undetected symptomatic
    num_symptomatic_undetected = sample_negative_binomial_limit(
        missing_detections, p_detect_given_symp
    )

    # infer some undetected asumptomatic
    num_symptomatic_total = num_symptomatic_detected + num_symptomatic_undetected

    num_asymptomatic_undetected = 0
    if num_symptomatic_total == 0
        num_asymptomatic_undetected = sample_negative_binomial_limit(1, p_symp)
    else
        num_asymptomatic_undetected = sample_negative_binomial_limit(num_symptomatic_total, p_symp)
    end

    # infection time is then 1 generation interval earlier
    for _ in 1:num_symptomatic_detected
        total_injected += 1
        onset_time = rand(X_day_range)
        inf_time = onset_time - sample_onset_time(onset_time >= omicron_dominant_day)
        Z[map_day_to_index_Z(inf_time), individual_type_map.S] += 1
        D[map_day_to_index_UD(onset_time), individual_type_map.S] += 1
    end

    # infection time is then 1 generation interval earlier
    for _ in 1:num_asymptomatic_detected
        total_injected += 1
        onset_time = rand(X_day_range)
        # -1 as the sampled onset time is ceiling and so technically the infection day
        # would be rounded down
        inf_time = onset_time - sample_onset_time(onset_time >= omicron_dominant_day)
        Z[map_day_to_index_Z(inf_time), individual_type_map.A] += 1
        D[map_day_to_index_UD(onset_time), individual_type_map.A] += 1
    end

    # now add in the inferred undetected cases in the same way. Note that the
    # onset times are kinda irrelevant here, but just help with consistent
    for _ in 1:num_symptomatic_undetected
        total_injected += 1
        onset_time = rand(X_day_range)
        inf_time = onset_time - sample_onset_time(onset_time >= omicron_dominant_day)
        Z[map_day_to_index_Z(inf_time), individual_type_map.S] += 1
        U[map_day_to_index_UD(onset_time), individual_type_map.S] += 1
    end

    # now add in the inferred undetected cases
    for _ in 1:num_asymptomatic_undetected
        total_injected += 1
        onset_time = rand(X_day_range)
        inf_time = onset_time - sample_onset_time(onset_time >= omicron_dominant_day)
        Z[map_day_to_index_Z(inf_time), individual_type_map.A] += 1
        U[map_day_to_index_UD(onset_time), individual_type_map.A] += 1
    end

    return nothing

end


text:
function leading_boundary(state::InfiniteMPS, H::InfiniteMPO, alg::GradientGrassmann, envs = environments(state,H))
    (multi,envs,err) = leading_boundary(convert(MPSMultiline,state),convert(MPOMultiline,H),alg,envs)
    state = convert(InfiniteMPS,multi)
    return (state,envs,err)
end

function leading_boundary(state::MPSMultiline, H, alg::GradientGrassmann,envs=environments(state, H))
    res = optimize(GrassmannMPS.fg, (state, envs), alg.method;
                   transport! = GrassmannMPS.transport!,
                   retract = GrassmannMPS.retract,
                   inner = GrassmannMPS.inner,
                   scale! = GrassmannMPS.scale!,
                   add! = GrassmannMPS.add!,
                   finalize! = alg.finalize!,
                   precondition = GrassmannMPS.precondition,
                   isometrictransport = true)
    (x, fx, gx, numfg, normgradhistory) = res
    (state, envs) = x
    return state, envs, normgradhistory[end]
end


text:
abstract type GraphEntity end
abstract type IPFSEntity <: GraphEntity end

struct GQLQuery
    args::Dict{AbstractString,Any}
    fields::AbstractVector{AbstractString}
end

struct Allocation <: IPFSEntity
    id::AbstractString
    ipfshash::AbstractString
    amount::Real
    created_at_epoch::Integer

    function Allocation(
        ipfshash::AbstractString, amount::AbstractString, created_at_epoch::Integer
    )
        return new("", ipfshash, togrt(amount), created_at_epoch)
    end
    function Allocation(ipfshash::AbstractString, amount::Real, created_at_epoch::Integer)
        return new("", ipfshash, amount, created_at_epoch)
    end
    function Allocation(id::AbstractString, ipfshash::AbstractString)
        return new(id, ipfshash, 0.0, 0)
    end
end

struct Indexer <: GraphEntity
    id::AbstractString
    stake::Real
    allocations::AbstractVector{Allocation}

    function Indexer(id, delegation::AbstractString, stake::AbstractString, allocation)
        return new(
            id,
            togrt(stake) + togrt(delegation),
            map(
                x -> Allocation(
                    x["subgraphDeployment"]["ipfsHash"],
                    x["allocatedTokens"],
                    x["createdAtEpoch"],
                ),
                allocation,
            ),
        )
    end
    function Indexer(id, stake::Real, allocation)
        return new(id, stake, allocation)
    end
    function Indexer(allocation)
        return new(
            "",
            0.0,
            map(x -> Allocation(x["id"], x["subgraphDeployment"]["ipfsHash"]), allocation),
        )
    end
end

struct SubgraphDeployment <: IPFSEntity
    id::AbstractString
    ipfshash::AbstractString
    signal::Real

    function SubgraphDeployment(id, ipfshash, signal::AbstractString)
        return new(id, ipfshash, togrt(signal))
    end
    SubgraphDeployment(id, ipfshash, signal::Real) = new(id, ipfshash, signal)
end

struct Repository
    indexers::AbstractVector{Indexer}
    subgraphs::AbstractVector{SubgraphDeployment}
end

verify_ipfshash(x::AbstractString) = startswith(x, "Qm") && length(x) == 46

function togrt(x)::Float64
    return parse(Float64, x) / 1e18
end

signal(s::SubgraphDeployment) = s.signal

ipfshash(x::IPFSEntity) = x.ipfshash

id(x::GraphEntity) = x.id

allocation(i::Indexer) = i.allocations

allocated_stake(a::Allocation) = a.amount


text:
# This test file is mostly about performance
# and probably will not be FactTested

using ProfileView
using Distributions
include(Pkg.dir("MethylUtils", "src","mlml-bs-ox-typed.jl"))
include(Pkg.dir("MethylUtils", "src","Features.jl"))

function rand_CytosineCounts(num_cgs, max_cov, min_cov=0)
   coverage_distribution = [ x for x= 1:max_cov]
   cgs = Array(CytosineCount,num_cgs)
   k=1
   while k <= num_cgs
       (Cs, Ts) = rand_c_counts(0.8, coverage_distribution)
       if ( Cs + Ts ) >= min_cov
           cgs[k] = CytosineCount(GenomePos("chr1",1,2),true,int16(Cs) ,int16(Ts))
           k +=1
       end
   end
   return cgs
end

function rand_bs_ox_counts(num_cgs, max_cov, min_cov=0)
   coverage_distribution = [ x for x= 1:max_cov]
   cgs = Array(Array{Int16},num_cgs)
   k=1
   while k <= num_cgs
       (M_B, u_B) = rand_c_counts(0.8, coverage_distribution)
       (m_O, U_B) = rand_c_counts(0.77,coverage_distribution)

       if( ( M_B + u_B ) >= min_cov && (m_O + U_B ) >= min_cov )
           meth_counts = map( int16, [M_B, u_B, m_O, U_B ] )
           cgs[k] = meth_counts
           k += 1
       end
   end
   return cgs
end

function rand_c_counts(prob, coverage_distribution)
    num_reads = sample(coverage_distribution)
    binom = Binomial(num_reads,prob)
    Cs = rand(binom)
    Ts = num_reads - Cs
    return (Cs,Ts)
end

function start_point(M_B,u_B,m_O,U_O)
        starting_estimate_p_m = M_B / (M_B + u_B)
        starting_estimate_p_h = 1 - starting_estimate_p_m
        return ( starting_estimate_p_m, starting_estimate_p_h )
end

function mlml_performance(meth_counts)
    num_cgs = length(meth_counts)
    new_array = Array(Float64, num_cgs)
    for k = 1:num_cgs
        (M_B,u_B,m_O,U_O) = meth_counts[k]
        (starting_estimate_p_m, starting_estimate_p_h) = start_point(M_B,u_B,m_O,U_O)
        # if we have overshoot then calculate
        #p_m_O = m_O / ( m_O + U_O )
        #if (starting_estimate_p_m - p_m_O) < 0
        (em_p_m, em_p_h, iter) = EM(M_B, u_B, m_O, U_O, starting_estimate_p_m, starting_estimate_p_h)
        new_array[k] = em_p_m
        #end
    end
    return new_array
end

# THIS WILL HAVE TO BE IN ANOTHER FILE TO REQUIRE
function mlml_parallel_performance(d::DArray, meth_counts)
   DArray(size(d),procs(d)) do I
         yrange = I[1]
         ymin=yrange[1]
         #xrange = I[2]
         #xmin=xrange[1]

         # HERE
         #new_array = Array(Float64, size(yrange,1) )
         new_array = Array(Array{Float64}, size(yrange,1))
         # here should be 2-d for both em_p_m and em_p_h and column for conflict or iteration
         # or rather than matrix perhaps just push on array of arrays?

         num_cgs = length(meth_counts)
         for y in yrange
            (M_B,u_B,m_O,U_O) = meth_counts[y]
            (starting_estimate_p_m, starting_estimate_p_h) = start_point(M_B,u_B,m_O,U_O)
            (em_p_m, em_p_h, iter) = EM(M_B, u_B, m_O, U_O, starting_estimate_p_m, starting_estimate_p_h)
            #put onto new array
            new_array[ y - ymin + 1 ] = [ em_p_m, em_p_h]
         end
          return new_array
    end
end

function mlml_parallel_performance(d::DArray, bs_counts , ox_counts )
   DArray(size(d),procs(d)) do I
         yrange = I[1]
         ymin=yrange[1]
         #xrange = I[2]
         #xmin=xrange[1]

         # HERE
         #new_array = Array(Float64, size(yrange,1) )
         len = size(yrange,1)
         println("size of array is $len")
         new_array = Array(CytosineMethylation, len)
         # here should be 2-d for both em_p_m and em_p_h and column for conflict or iteration
         # or rather than matrix perhaps just push on array of arrays?

         #num_cgs = length(bs_counts)
         for y in yrange
            (M_B,u_B) = ( bs_counts[y].C_count, bs_counts[y].T_count)
            (m_O,U_O) = ( ox_counts[y].C_count, ox_counts[y].T_count)
            (starting_estimate_p_m, starting_estimate_p_h) = start_point(M_B,u_B,m_O,U_O)
            (em_p_m, em_p_h, iter) = EM(M_B, u_B, m_O, U_O, starting_estimate_p_m, starting_estimate_p_h)
            #put onto new array
            new_array[ y - ymin + 1 ] = CytosineMethylation( bs_counts[y].pos, bs_counts[y].strand, em_p_m, em_p_h )
         end
         return new_array
    end
end


text:
using Slabs
using Test

# For testing `datamerge`
struct Merger
    val::Int
end
Slabs.datamerge(a::Merger, b::Merger) = Merger(a.val + b.val)

@testset "NamedTuple Utilities" begin
    # NT Sorting
    nt = (b = 1, c = 2, a = 3)
    @test Slabs.ntsort(nt) == (a = 3, b = 1, c = 2)

    #####
    ##### Matching
    #####

    # Expect a match
    a = (a = 1, b = nothing)
    @test Slabs.ismatch(a, (a = 1,))
    @test Slabs.ismatch(a, (a = 1, b = nothing))
    @test Slabs.ismatch(a, (b = nothing, a = 1))
    @test Slabs.ismatch(a, (b = nothing,))

    # Don't expect a match
    @test !Slabs.ismatch(a, (c = 1,))
    @test !Slabs.ismatch(a, (a = 1, b = 1))

    # Make sure the `data` field is not checked.
    a = (a = 1, b = nothing, data = -1234)
    @test Slabs.ismatch(a, (a = 1, data = "hello wold"))

    # Empty NamedTuple should always match
    @test Slabs.ismatch(a, NamedTuple())

    #####
    ##### Expanding
    #####

    # Make sure the result is sorted
    a = (c = 1, b = 2)
    b = Slabs.expand(a, (:a, :d))
    @test b == (a = nothing, b = 2, c = 1, d = nothing)

    # Try again with a different callback.
    b = Slabs.expand(a, (:a, :d); f = () -> 0)
    @test b == (a = 0, b = 2, c = 1, d = 0)

    #####
    ##### pushconvert
    #####

    x = [1,2,3]
    @test Slabs.pushconvert(x, nothing) == [1, 2, 3, nothing]
    @test Slabs.pushconvert(x, "hi") == [1, 2, 3, "hi"]

    #####
    ##### Slab
    #####

    S = Slab()
    println(S)
    @test length(S) == 0
    @test eltype(S) == typeof(NamedTuple())

    data1 = Dict(:merger => Merger(10), :meta => 1)
    S[(a = 1, b = 2)] = data1

    @test length(S) == 1

    # Test some of the indexing methods
    @test S[1] == (a = 1, b = 2, data = data1)
    @test S[nothing] == S
    @test S[:a] == [1]
    @test S[(:a, :b)] == [(1, 2)]

    # Add another entry
    data2 = Dict(:merger => Merger(100), :meta => 2)
    S[(c = 10,)] = data2
    @test propertynames(S.slab) == (:a, :b, :c, :data)
    @test S[1] == (a = 1, b = 2, c = nothing, data = data1)

    # Try filtering
    filt = S[(c = nothing,)]
    @test isa(filt, Slab)
    @test length(filt) == 1
    @test filt[1] == S[1]

    # Do the same thing with the keyword path.
    filt = S[c = nothing]
    @test isa(filt, Slab)
    @test length(filt) == 1
    @test filt[1] == S[1]

    # Some more getindex testing
    @test S[:a] == [1, nothing]
    @test S[(:a, :b)] == [(1, 2), (nothing, nothing)]
    println(S)

    for i in S
        @test isa(i, NamedTuple)
    end

    # Finally, try datamerging.
    data3 = Dict(
        :merger => Merger(100),
        :meta => 2,
        :misc => 0,
    )

    S[(a = 1, b = 2)] = data3
    row = S[(a = 1, b = 2)][1]

    # The mergers should have added their values together
    @test row.data[:merger] == Merger(110)

    # Normal collisions will have been overwritten.
    @test row.data[:meta] == 2
    @test row.data[:misc] == 0

    # Finally, test that we get an exception if we have multiple matches.
    S[(a = 1, b = 3)] = data3
    @test_throws ErrorException S[a = 1] = data3
end


text:
module DynamicGridsInteract
# Use the README as the module docs
@doc read(joinpath(dirname(@__DIR__), "README.md"), String) DynamicGridsInteract

using Blink,
      DynamicGrids,
      ImageShow,
      Interact,
      InteractBase,
      InteractModels,
      Mux,
      WebIO,
      WebSockets

const DG = DynamicGrids

using DynamicGrids: AbstractSimData, SimData, rules, fps, init, tspan, isrunning, setrunning!

import Base: length, size, firstindex, lastindex, getindex, setindex!, push!, append!, parent, show

export AbstractWebOutput, InteractOutput, ElectronOutput, ServerOutput

include("interact.jl")
include("electron.jl")
include("server.jl")

end


text:
module DCAlign
using FastaIO
using Printf, Random, Logging, Distributed, DelimitedFiles
using ExtractMacro, OffsetArrays
export Seq, palign

import Base.show

include("types.jl")     # general types (all marginals and BP/MS msg)
include("iterate_bplc.jl")   # BPlc algoritm at chosen temperature
include("utils.jl")     # utilities
include("palign_bplc.jl")    # the main alignment function (BPlc at chosen temperature)
include("dataset.jl")   # read full length seqs or generate them
include("iterate_bplc0.jl") # BPlc update at zero temperature
include("decimation.jl") # all decimation routines
include("align_all.jl") # align a set of sequences

end #end module


text:
using Test
using Plots

opts = TrajectoryOptimization.SolverOptions()
opts.verbose = false
opts.cost_tolerance = 1e-6
opts.cost_tolerance_intermediate = 1e-6
opts.constraint_tolerance = 1e-8
opts.square_root = true
opts.outer_loop_update_type = :default
opts.constraint_tolerance_second_order_dual_update = sqrt(opts.constraint_tolerance)
opts.use_second_order_dual_update = true
opts.penalty_max = 1e8
opts.iterations_outerloop = 25
opts.live_plotting = false
model, obj = TrajectoryOptimization.Dynamics.double_integrator
u_min = -.15
u_max = 0.2
obj_con = TrajectoryOptimization.ConstrainedObjective(obj, tf=5.0, u_min=u_min, u_max=u_max)#, x_min=x_min, x_max=x_max)

integrator = :rk4
dt = 0.05
solver = TrajectoryOptimization.Solver(model,obj_con,integration=integrator,dt=dt,opts=opts)
U = rand(solver.model.m, solver.N)

results, stats = TrajectoryOptimization.solve(solver,U)

@show stats["iterations"]
@show stats["outer loop iterations"]
@show stats["c_max"][end]
@show max_violation(results) < opts.constraint_tolerance

plot(stats["c_max"],yscale=:log10,ylabel="log(c_max)",label="1st")

plot(to_array(results.U)',label="")
plot(to_array(results.X)',label="")


text:
scatter(sin, 0, 2π; msize=0.1)
xlim(0,2π)


text:
function Init_cold_4D(NX,NY,NZ,NT,Nwing,NC)
    Dim = 4

    U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")

    temp1 = similar(U[1])
    temp2 = similar(U[1])

    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)

    #factor = 2/(U[1].NV*4*3*U[1].NC)
    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    println("plaq_t = $plaq_t")
    poly = calculate_Polyakov_loop(U,temp1,temp2)
    println("polyakov loop = $(real(poly)) $(imag(poly))")
    return plaq_t

end

function Init_hot_4D(NX,NY,NZ,NT,Nwing,NC)
    Random.seed!(123)
    Dim = 4

    U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "hot") #for debug
    U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "hot",randomnumber="Reproducible")

    temp1 = similar(U[1])
    temp2 = similar(U[1])

    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)

    #factor = 2/(U[1].NV*4*3*U[1].NC)
    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    println("plaq_t = $plaq_t")
    poly = calculate_Polyakov_loop(U,temp1,temp2)
    println("polyakov loop = $(real(poly)) $(imag(poly))")
    return plaq_t

end

function Init_ildg_4D(NX,NY,NZ,NT,Nwing,NC,filename)
    Random.seed!(123)
    Dim = 4

    U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "hot")
    save_binarydata(U,filename)

    U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")
    ildg = ILDG(filename)
    i = 1
    L = [NX,NY,NZ,NT]
    load_gaugefield!(U,i,ildg,L,NC)

    temp1 = similar(U[1])
    temp2 = similar(U[1])

    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)
    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    poly = calculate_Polyakov_loop(U,temp1,temp2)
    println("polyakov loop = $(real(poly)) $(imag(poly))")
    return plaq_t

end

function Init_cold_2D(NX,NT,Nwing,NC)
    Dim = 2

    U = Initialize_Gaugefields(NC,Nwing,NX,NT,condition = "cold")

    temp1 = similar(U[1])
    temp2 = similar(U[1])

    if Dim == 4
        comb = 6 #4*3/2
    elseif Dim == 3
        comb = 3
    elseif Dim == 2
        comb = 1
    else
        error("dimension $Dim is not supported")
    end
    factor = 1/(comb*U[1].NV*U[1].NC)

    #factor = 2/(U[1].NV*4*3*U[1].NC)
    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    println("plaq_t = $plaq_t")
    poly = calculate_Polyakov_loop(U,temp1,temp2)
    println("polyakov loop = $(real(poly)) $(imag(poly))")
    return plaq_t

end

function Init_hot_2D(NX,NT,Nwing,NC)
    Random.seed!(123)
    Dim = 2

    #U = Initialize_Gaugefields(NC,Nwing,NX,NT,condition = "hot")
    U = Initialize_Gaugefields(NC,Nwing,NX,NT,condition = "hot",randomnumber="Reproducible")

    temp1 = similar(U[1])
    temp2 = similar(U[1])

    if Dim == 4
        comb = 6 #4*3/2
    elseif Dim == 3
        comb = 3
    elseif Dim == 2
        comb = 1
    else
        error("dimension $Dim is not supported")
    end
    factor = 1/(comb*U[1].NV*U[1].NC)

    #factor = 2/(U[1].NV*4*3*U[1].NC)
    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    println("plaq_t = $plaq_t")
    poly = calculate_Polyakov_loop(U,temp1,temp2)
    println("polyakov loop = $(real(poly)) $(imag(poly))")
    return plaq_t

end

@testset "cold start" begin
    println("cold start")
    println("4D system")
    @testset "4D" begin
        NX = 4
        NY = 4
        NZ = 4
        NT = 4
        Nwing = 1

        @testset "NC=2" begin
            NC = 2
            println("NC = $NC")
            plaq_t = Init_cold_4D(NX,NY,NZ,NT,Nwing,NC)

            @test plaq_t == one(plaq_t)
        end

        @testset "NC=3" begin
            NC = 3
            println("NC = $NC")
            plaq_t = Init_cold_4D(NX,NY,NZ,NT,Nwing,NC)
            @test plaq_t == one(plaq_t)
        end

        @testset "NC=4" begin
            NC = 4
            println("NC = $NC")
            plaq_t = Init_cold_4D(NX,NY,NZ,NT,Nwing,NC)
            @test plaq_t == one(plaq_t)
        end

        @testset "NC=5" begin
            NC = 5
            println("NC = $NC")
            plaq_t = Init_cold_4D(NX,NY,NZ,NT,Nwing,NC)
            @test plaq_t == one(plaq_t)
        end
    end

    println("2D system")
    @testset "2D" begin
        NX = 4
        #NY = 4
        #NZ = 4
        NT = 4
        Nwing = 1

        @testset "NC=2" begin
            NC = 2
            println("NC = $NC")
            plaq_t = Init_cold_2D(NX,NT,Nwing,NC)

            @test plaq_t == one(plaq_t)
        end

        @testset "NC=3" begin
            NC = 3
            println("NC = $NC")
            plaq_t = Init_cold_2D(NX,NT,Nwing,NC)
            @test plaq_t == one(plaq_t)
        end

        @testset "NC=4" begin
            NC = 4
            println("NC = $NC")
            plaq_t = Init_cold_2D(NX,NT,Nwing,NC)
            @test plaq_t == one(plaq_t)
        end

        @testset "NC=5" begin
            NC = 5
            println("NC = $NC")
            plaq_t = Init_cold_2D(NX,NT,Nwing,NC)
            @test plaq_t == one(plaq_t)
        end
    end
end

eps = 1e-8

@testset "hot start" begin
    println("hot start")
    println("4D system")
    @testset "4D" begin
        NX = 4
        NY = 4
        NZ = 4
        NT = 4
        Nwing = 1

        @testset "NC=2" begin
            NC = 2
            println("NC = $NC")
            plaq_t = Init_hot_4D(NX,NY,NZ,NT,Nwing,NC)
            #val = -0.007853743153861802
            val = 0.0002565633289190507
            @test abs(plaq_t-val)/abs(val) < eps
        end

        @testset "NC=3" begin
            NC = 3
            println("NC = $NC")
            plaq_t = Init_hot_4D(NX,NY,NZ,NT,Nwing,NC)
            #val = 0.0015014233929744197
            val = 0.008449494077606137
            @test abs(plaq_t-val)/abs(val) < eps
        end

        @testset "NC=4" begin
            NC = 4
            println("NC = $NC")
            plaq_t = Init_hot_4D(NX,NY,NZ,NT,Nwing,NC)
            #val = -0.004597227507817238
            val =  0.0029731967994215515
            @test abs(plaq_t-val)/abs(val) < eps
        end

        @testset "NC=5" begin
            NC = 5
            println("NC = $NC")
            plaq_t = Init_hot_4D(NX,NY,NZ,NT,Nwing,NC)
            #val = 0.0037580826460029506
            val = -0.0005506443879866896
            @test abs(plaq_t-val)/abs(val) < eps
        end
    end

    println("2D system")
    @testset "2D" begin
        NX = 4
        #NY = 4
        #NZ = 4
        NT = 4
        Nwing = 1

        @testset "NC=2" begin
            NC = 2
            println("NC = $NC")
            plaq_t = Init_hot_2D(NX,NT,Nwing,NC)
            #val = 0.022601163616639157
            val  =-0.09978762099361757
            @test abs(plaq_t-val)/abs(val) < eps
        end

        @testset "NC=3" begin
            NC = 3
            println("NC = $NC")
            plaq_t = Init_hot_2D(NX,NT,Nwing,NC)
            #val = -0.04647124293538649
            val = 0.046263787757668186
            @test abs(plaq_t-val)/abs(val) < eps
        end

        @testset "NC=4" begin
            NC = 4
            println("NC = $NC")
            plaq_t = Init_hot_2D(NX,NT,Nwing,NC)
            #val = 0.07457370324173362
            val = 0.03150732915218983
            @test abs(plaq_t-val)/abs(val) < eps
        end

        @testset "NC=5" begin
            NC = 5
            println("NC = $NC")
            plaq_t = Init_hot_2D(NX,NT,Nwing,NC)
            #val = -0.013511504030861661
            val = -0.03564420387393341
            @test abs(plaq_t-val)/abs(val) < eps
        end
    end

end

@testset "one instanton" begin
    @testset "4D" begin
        NX = 4
        NY = 4
        NZ = 4
        NT = 4
        NC = 2
        Nwing = 1
        U = Oneinstanton(NC,NX,NY,NZ,NT,Nwing)

        comb = 6
        factor = 1/(comb*U[1].NV*U[1].NC)

        temp1 = similar(U[1])
        temp2 = similar(U[1])

        @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        println("plaq_t = $plaq_t")
        val =  0.9796864531099871
        @test abs(plaq_t-val)/abs(val) < eps
    end

    @testset "2D" begin
        Dim = 2
        NX = 4
        #NY = 4
        #NZ = 4
        NT = 4
        NC = 2
        Nwing = 1
        U = Oneinstanton(NC,NX,NT,Nwing)

        if Dim == 4
            comb = 6 #4*3/2
        elseif Dim == 3
            comb = 3
        elseif Dim == 2
            comb = 1
        else
            error("dimension $Dim is not supported")
        end
        factor = 1/(comb*U[1].NV*U[1].NC)

        temp1 = similar(U[1])
        temp2 = similar(U[1])

        @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        println("plaq_t = $plaq_t")
        val = 0.9300052284270868
        @test abs(plaq_t-val)/abs(val) < eps
    end
end

#=
@testset "File start" begin
    @testset "4444 SU(2)" begin
        Dim = 4
        NX = 4
        NY = 4
        NZ = 4
        NT = 4
        NC = 2
        Nwing = 1
        #U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "hot")

        filename = "./data/conf_00000100_4444_test.ildg"
        plaq_t = Init_ildg_4D(NX,NY,NZ,NT,Nwing,NC,filename)

        println("plaq_t = $plaq_t")
        val = -0.007853743153861798
        @test abs(plaq_t-val)/abs(val) < eps
    end
end
=#


text:
# Clear all objects in the scene and delete all meshes
function clear_objects!()
    length(data.meshes) == 0 && return
    object_ops.mode_set(mode="OBJECT")
    object_ops.select_by_type(; :type => "MESH")
    object_ops.delete(use_global = false)

    for m in data.meshes
        data.meshes.remove(m, do_unlink = true)
    end
end

# https://blender.stackexchange.com/a/91687
function smooth!(ob)
    for f in ob.data.polygons
        f.use_smooth = true
    end
end

function add_object!(name::String, vertices, normals, faces; smooth=false)
    mesh = data.meshes.new("$(name)_mesh")
    object = data.objects.new("$(name)_obj", mesh)

    scene = context.scene
    scene.collection.objects.link(object)
    context.view_layer.objects.active = object
    # object.select = true

    facesm1 = map(f -> f .- 1, faces)

    mesh.from_pydata(vertices,
                     [],
                     facesm1)
    mesh.update()

    smooth && smooth!(object)
    object
end

function add_object!(name::String, m::Mesh{3,T};
                     trf::Transformation=IdentityTransformation(),
                     kwargs...) where T
    verts = map(trf.(decompose(Point3{T}, m))) do v
        v[1],v[2],v[3]
    end

    # Need to transform normals properly
    norms = map(decompose(Normal(), m)) do n
        n[1],n[2],n[3]
    end

    faces = map(decompose(TriangleFace{UInt32}, m)) do f
        f[1],f[2],f[3]
    end

    add_object!(name, verts, norms, faces; kwargs...)
end

function add_object!(name::String, p::GeometryBasics.GeometryPrimitive;
                     nvertices=24, kwargs...)
    m = Mesh(collect(coordinates(p, nvertices)),
             collect(faces(p, nvertices)))
    add_object!(name, m; kwargs...)
end

setcolor!(o, c::Colorant) =
    o.color = (red(c), green(c), blue(c), alpha(c))

function vertex_colors!(fun::Function, ob)
    # https://blender.stackexchange.com/a/911
    mesh = ob.data
    isempty(mesh.vertex_colors) && mesh.vertex_colors.new()
    cl = mesh.vertex_colors.active
    for (i,p) in enumerate(mesh.polygons)
        for (j,vi) in zip(p.loop_indices,p.vertices)
            setcolor!(cl.data[j+1], fun(vi+1))
        end
    end
end

function append_material!(ob, mat; activate=true)
    mesh = ob.data
    mesh.materials.append(mat)
    if activate
        ob.active_material_index = length(mesh.materials)-1
    end
end

export clear_objects!, smooth!, add_object!,
    vertex_colors!, append_material!


text:
module Lexical

# REFERENCES
#
#    [1] ISO 8879:1986 Information processing -- Text and office systems -- Standard Generalized Markup Language (SGML)
#        https://www.iso.org/standard/16387.html
#    [2] Extensible Markup Language (XML) 1.0 (Fifth Edition),
#        https://www.w3.org/TR/xml/

# ----------------------------------------
# EXPORTED INTERFACE
# ----------------------------------------

export tokens
export location_of

# ----------------------------------------
# TOKEN TYPES
# ----------------------------------------

@enum TokenType begin
    # These token names and values come from the SGML specification. See [1].
    #
    mdo     # markup delimiter open ... <!
    mdc     # markup delimiter close ... > ... this one will never show up ... we'll use tagc instead.
    dso     # declaration subset open ... [
    dsc     # declaration subset close ... ]
    msc     # marked section close ... ]]
    com     # comment ... --
    rni     # reserved name indicator ... #
    lit     # literal ... "
    lita    # alternative literal ... '
    grpo    # group open ... (
    grpc    # group close ... )
    and     # and connector ... & ... this is not actually used in XML, and will never show up because ero will be found instead.
    or      # or connector ... |
    seq     # sequence connector ... ,
    opt     # optional occurrence indicator ... ?
    rep     # zero-or-more occurrence indicator ... *
    plus    # one-or-more occurrence indicator ... +
    minus   # exclusion/omission flag ... - ... this is not actually used in XML.
    cro     # character reference open ... &#
    ero     # entity reference open ... &
    pero    # parameter entity reference open ... %
    refc    # reference close ... ;
    pio     # processing instruction open ... <?
    pic     # processing instruction close ... > ... this is actually ?> in XML ... but see below in TwoCharacterTokens.
    stago   # start tag open ... <
    etago   # end tag open ... </
    tagc    # tag close ... >
    net     # null end tag ... /
    vi      # value indicator ... =

    # The rest of these don't represent traditional SGML tokens.
    #
    ws      # XML white-space ... see [2], § 2.3
    eoi     # end of input
    text    # anything else, except ...

    illegal # this is not actually used in XML.
end

# ----------------------------------------
# TYPE DECLARATIONS
# ----------------------------------------

struct Location
   identification ::String
   line_number    ::Int64
end

mutable struct State
   io          ::Union{IOStream, IOBuffer}
   line_number ::Int64
   last_match  ::Union{Nothing, Array{Char, 1}}
   length_of   ::Int64
end

struct Token
   token_type ::TokenType
   value      ::String
   location   ::Location
end

# ----------------------------------------
# CONSTANTS
# ----------------------------------------

# In XML, pic should be ?>, as compared to > in SGML. However, we can't actually tokenize ?> as a single token, because
# then element declarations that happen to be optional (e.g., <!ELEMENT a (b)?>) would be tokenized incorrectly (unless
# there happens to be space between the two characters) and then we'd have to bend over backwards in the element
# declaration parser. What kills me is that according to [1], § 1.1, "It shall be easy to write programs which process
# XML documents.". If they'd left pic as >, like in SGML, it would've been easier to parse: you'd end a PI by looking
# for tagc. Otherwise as it stands you'd need a context-sensitive lexical layer to distinguish the token ?> from the two
# tokens ? and >. And in the same section, "XML shall be compatible with SGML." ... well it isn't: PIs in SGML by
# default end with a tagc (i.e., >), not with the sequence ?>. Anyway, that's my tirade for today.
#
const TwoCharacterTokens = Dict([ '<', '!' ] => mdo,
                                [ ']', ']' ] => msc,
                                [ '-', '-' ] => com,
                                [ '&', '#' ] => cro,
                                [ '<', '?' ] => pio,
                                [ '<', '/' ] => etago)
const OneCharacterTokens = Dict('>'  => mdc, # Again, this will never show up ... we'll always emit tagc instead.
                                '['  => dso,
                                ']'  => dsc,
                                # '#'  => rni, ... Including this as a token in XML makes parsing reserved names a pain.
                                '\"' => lit,
                                '\'' => lita,
                                '('  => grpo,
                                ')'  => grpc,
                                '&'  => and, # Again, this will never show up ... we'll always emit ero instead.
                                '|'  => or,
                                ','  => seq,
                                '?'  => opt,
                                '*'  => rep,
                                '+'  => plus,
                                # '-'  => minus, ... This is not a token in XML, and including it causes too many headaches.
                                '&'  => ero,
                                '%'  => pero,
                                ';'  => refc,
                                '<'  => stago,
                                '>'  => tagc,
                                '/'  => net,
                                '='  => vi)
const IllegalCharacters = Set(vcat('\u00':'\u08', '\u0e':'\u1f', '\ud800':'\udfff', [ '\u0b', '\u0c', '\ufffe', '\uffff' ]))
const TokenStarts       = Set([ '\"', '%', '&', '\'', '(', ')', '*', '+', ',', '-', '/', ';', '<', '=', '>', '?', '[', ']', '|' ])
const WhiteSpaces       = Set([ '\u20', '\u09', '\u0a', '\u0d' ])

# ----------------------------------------
# FUNCTIONS
# ----------------------------------------

State(io) = State(io, -1, nothing, -1)
Token(token_type::TokenType, token_value::String, state::State) = Token(token_type, token_value, location_of(state))
location_of(token::Token) = token.location

function consume_last_match(state::State)
    skip(state.io, state.length_of)

    return state.last_match
end

function consume_until(state::State, end_markers::Set{Char})
    consumed   = Array{Char, 1}()
    start      = position(state.io)
    bytes_read = 0

    while !eof(state.io)
        c = read(state.io, Char)
        if c ∈ end_markers
            if c == '-'
                # Be careful here: '-' is a token in SGML (minus), but in XML we need to look ahead to see if it is
                # followed by another '-'. If so, then "--" is a token (com), but only if it, in turn, is followed by
                # '>' (mdc/tagc). Otherwise it is part of the current token (text).
                #
                # Ultimately, it might be easier to handle this in an other layer, by consuming a sequence of tokens and
                # grouping then. But this will do for now.
                #
                if eof(state.io)
                    bytes_read = position(state.io) - start
                    push!(consumed, c)
                    break
                end

                current_position = position(state.io)

                if read(state.io, Char) == '-' && !eof(state.io) && read(state.io, Char) ∈ [ '-', '>' ]
                    break

                else
                    # Keep going.
                    #
                    seek(state.io, current_position)
                end

            else
                break
            end
        end
        bytes_read = position(state.io) - start
        push!(consumed, c)
    end

    return ( bytes_read, consumed )
end

function consume_while(state::State, set::Set{Char})
    consumed   = Array{Char, 1}()
    start      = position(state.io)
    bytes_read = 0

    while !eof(state.io)
        c = read(state.io, Char)
        if c ∉ set
            break
        end
        bytes_read = position(state.io) - start
        push!(consumed, c)
    end

    return ( bytes_read, consumed )
end

function location_of(state::State)
    name_of(io::IOStream) = io.name
    name_of(io) = "a buffer"

    return Location(name_of(state.io), -1)
end

function is_delimiter_one(state::State)::Bool
    start = mark(state.io)
    token_value = read(state.io, Char)
    current = position(state.io)
    reset(state.io)

    if haskey(OneCharacterTokens, token_value)
        state.last_match = [ token_value ]
        state.length_of  = current - start

        return true

    else
        return false
    end
end

function is_delimiter_two(state::State)::Bool
    start = mark(state.io)
    first = read(state.io, Char)

    if eof(state.io)
        reset(state.io)

        return false

    else
        second = read(state.io, Char)
        current = position(state.io)
        reset(state.io)

        if haskey(TwoCharacterTokens, [ first, second ])
            state.last_match = [ first, second ]
            state.length_of  = current - start

            return true

        else
            return false
        end
    end
end

function is_text(state::State)::Bool
    start = mark(state.io)
    ( bytes_read, token_value ) = consume_until(state, ∪(TokenStarts, WhiteSpaces, IllegalCharacters))
    reset(state.io)

    if length(token_value) > 0
        state.last_match = token_value
        state.length_of  = bytes_read

        return true

    else
        return false
    end
end

function is_white_space(state::State)::Bool
    start = mark(state.io)
    ( bytes_read, token_value ) = consume_while(state, WhiteSpaces)
    reset(state.io)

    if length(token_value) > 0
        state.last_match = token_value
        state.length_of  = bytes_read

        return true

    else
        return false
    end
end

function tokens(state::State)
    function next(state::State)::Union{Token, Nothing}
        if is_delimiter_two(state)
            token_value = consume_last_match(state)

            return Token(TwoCharacterTokens[token_value], String(token_value), state)

        elseif is_delimiter_one(state)
            token_value = consume_last_match(state)

            return Token(OneCharacterTokens[token_value[1]], String(token_value), state)

        elseif is_text(state)
            return Token(text, String(consume_last_match(state)), state)

        elseif is_white_space(state)
            return Token(ws, String(consume_last_match(state)), state)

        else
            return Token(illegal, string(read(state.io, Char)), state)
        end
    end

    function tokenized(channel::Channel)
        while !eof(state.io)
            put!(channel, next(state))
        end
    end

    return Channel(tokenized; ctype = Token, csize = 1)
end

end

max_stars_repo_name: jlegare/NativeXML

max_issues_repo_name: jlegare/NativeXML

max_forks_repo_name: jlegare/NativeXML


text:
module DLLBFGS

using JSOSolvers, NLPModels, SolverTools, LinearOperators

using LinearAlgebra, Logging

include("solver.jl")
end # module


text:
abstract type AbstractScene{T} end
abstract type AbstractLayer{T} end
abstract type AbstractWindow{T} end

struct Event{TYPE}
    fields::Dict{Symbol, Any}
end
Event{TYPE}() where {TYPE} = Event{TYPE}(Dict{Symbol, Any}())
Event{TYPE}(pairs::Pair...) where {TYPE}  = Event{TYPE}(Dict{Symbol, Any}(pairs...))

Base.getproperty(e::Event, name::Symbol) = getfield(e, :fields)[name]
Base.setproperty!(e::Event, name::Symbol, x) = setindex!(getfield(e, :fields), x, name)

mutable struct Loop
    name::String
    target_speed::Float64
    paused::Bool
    task::Task
    Loop(name::String, target_speed=50.0) = new(name, target_speed, false)
end

Base.schedule(loop::Loop) = schedule(loop.task)
Base.wait(loop::Loop) = wait(loop.task)

struct Timeout
    t1::Float64
    fn::Function
end
Base.isless(a::Timeout, b::Timeout) = isless(b.t1, a.t1)

"""

A Window.

"""
mutable struct Window{T} <: AbstractWindow{T}
    window_ptr::Ptr{SDL.Window}
    render_ptr::Ptr{SDL.Renderer}
    scene_stack::T
    loops::Dict{String, Loop}
    event_queue::Vector{Event}
    timer_queue::Vector{Timeout}

    function Window{T}(title, width, height, scene_stack::T; fullscreen=false) where {T}
        window_ptr = SDL.CreateWindow(title, Int32(SDL.WINDOWPOS_CENTERED_MASK), Int32(SDL.WINDOWPOS_CENTERED_MASK), Int32(width), Int32(height), fullscreen ? SDL.WINDOW_FULLSCREEN : UInt32(0))
        render_ptr = SDL.CreateRenderer(window_ptr, Int32(-1), UInt32(0))
        if render_ptr == C_NULL
            SDL.DestroyWindow(window_ptr)
            SDL.DestroyRenderer(render_ptr)
            error(unsafe_string(SDL.GetError()))
        end
        SDL.SetRenderDrawBlendMode(render_ptr, SDL.BLENDMODE_BLEND)
        loops = Dict{String,Loop}()
        event_queue = Event[]
        timer_queue = Timeout[]
        self = new(window_ptr, render_ptr, scene_stack, loops, event_queue, timer_queue)
        finalizer(destroy!, self)
        return self
    end
end

Window(title, width, height; kwargs...) = Window(title, width, height, AbstractScene[]; kwargs...)
function Window(title, width, height, scene_stack::T; fullscreen=false) where {T}
    return Window{T}(title, width, height, scene_stack, fullscreen=fullscreen)
end

function Window(title, width, height, scenes::AbstractScene...; kwargs...)
    return Window(title, width, height, [scenes...]; kwargs...)
end

Base.show(io::IO, window::Window) = print(io, "Window(\"", unsafe_title(window), "\", ", join(size(window), ", "), ")")

function destroy!(window::Window)
    empty!(window.scene_stack)
    SDL.DestroyWindow(window.window_ptr)
    SDL.DestroyRenderer(window.render_ptr)
    return nothing
end

function Base.getproperty(window::Window, name::Symbol)
    name == :width && return size(window)[1]
    name == :height && return size(window)[2]
    getfield(window, name)
end

activescene(window::Window) = window.scene_stack[end]

"""

`Scene` objects run inside a `Window`.  Only one scene can be active
inside a given `Window` at any time.

"""
mutable struct Scene{T <: AbstractVector{<:AbstractLayer}} <: AbstractScene{T}
    layers::T
    color::Colors.RGB
end

Scene(layer::AbstractLayer, layers::AbstractLayer...; color=Colors.colorant"darkgray") = Scene([layer, layers...], color)
Scene(layers::T) where {T<:AbstractVector{<:AbstractLayer}} = Scene(layers, Colors.colorant"darkgray")

struct RenderTask{T}
    source::T
    x::Float64
    y::Float64
    θ::Float64
    scale::Float64
    offset_x::Int
    offset_y::Int
    flip::Symbol
    color::Colors.RGB
end

RenderTask(source::T, args...) where {T} = RenderTask{T}(source, args...)

"""

A Layer object runs inside a Scene.

"""
mutable struct Layer{T} <: AbstractLayer{T}
    objects::OrderedSet{T}
    x::Float64
    y::Float64
    θ::Float64
    scale::Float64
    show::Bool
    render_tasks::Vector{RenderTask}
    new_objects::OrderedSet{T}
    dead_objects::Set{T}

    function Layer{T}(objects, x = 0.0, y = 0.0, θ = 0.0; show = true, scale = 1.0) where {T}
        return new(objects, x, y, θ, scale, show, RenderTask[], OrderedSet{T}(), Set{T}())
    end
end
Layer(objects::OrderedSet{T}, args...; kwargs...) where {T} = Layer{T}(objects, args...; kwargs...)
Layer(objects::Vector, args...; kwargs...) = Layer(OrderedSet(objects), args...; kwargs...)

visible(layer::Layer) = layer.show

Base.show(io::IO, layer::AbstractLayer) = print(io, "Layer(", join([layer.objects, layer.x, layer.y, layer.θ], ", "), ")")

Base.iterate(layer::AbstractLayer) = iterate(layer.objects)
Base.iterate(layer::AbstractLayer, state) = iterate(layer.objects, state)
Base.filter(f, layer::AbstractLayer) = filter(f, layer.objects)

"""
    wait(window::Window)

Block the current task until `window` has finished executing.

"""
function Base.wait(window::Window)
    for (name, loop) in window.loops
        wait(loop)
    end
    return nothing
end

"""
    size(window::Window)

Return a tuple containing the dimensions of `window`.

"""
function Base.size(window::Window)
    width, height = Int32[0], Int32[0]
    SDL.GetWindowSize(window.window_ptr, pointer(width), pointer(height))
    return width[], height[]
end

unsafe_title(window::Window) = unsafe_string(SDL.GetWindowTitle(window.window_ptr))

Base.push!(window::Window, scenes::AbstractScene...) = push!(window.scene_stack, scenes...)
Base.pop!(window::Window) = pop!(window.scene_stack)
Base.append!(window::Window, scenes::AbstractVector{<:AbstractScene}) = append!(window.scene_stack, scenes)
Base.push!(scene::AbstractScene, layers::AbstractLayer...) = push!(scene.layers, layers...)
Base.pop!(scene::AbstractScene) = pop!(scene.layers)
Base.append!(scene::AbstractScene, layers::AbstractVector{<:AbstractLayer}) = append!(scene.layers, layers)
Base.push!(layer::AbstractLayer, objs...) = (push!(layer.objects, objs...); layer)
Base.delete!(layer::AbstractLayer, obj) = (delete!(layer.objects, obj); layer)

Base.getindex(layer::Layer, index) = getindex(layer.objects, index)

isalive(layer::Layer, obj) = (obj in layer.objects && !(obj in layer.dead_objects))

"""

"""
add!(layer::AbstractLayer, objs...) = (push!(layer.new_objects, objs...); layer)

"""

"""
function populate!(layer::AbstractLayer)
    union!(layer.objects, layer.new_objects)
    empty!(layer.new_objects)
    return layer
end

"""

"""
kill!(layer::AbstractLayer, objs...) = (push!(layer.dead_objects, objs...); layer)
killall!(layer::AbstractLayer) = (union!(layer.dead_objects, layer.objects); layer)

"""

"""
function cleanup!(layer::AbstractLayer)
    setdiff!(layer.objects, layer.dead_objects)
    empty!(layer.dead_objects)
    return layer
end

tocoordinates(layer::Layer, x, y) = x - layer.x, y - layer.y