text:
##################################################################################
#
#  ORTHOFILTER TRANSFORMS
#  Periodic boundaries, Orthogonal, dyadic length (powers of 2)
#
##################################################################################

function makeqmf(h::AbstractVector, fw::Bool, T::Type=eltype(h))
	scfilter, dcfilter = makereverseqmf(h, fw, T)
    return reverse(scfilter), reverse(dcfilter)
end
function makereverseqmf(h::AbstractVector, fw::Bool, T::Type=eltype(h))
	h = convert(Vector{T}, h)
    if fw
    	scfilter = reverse(h)
        dcfilter = mirror(h)
    else
        scfilter = h
        dcfilter = reverse(mirror(h))
    end
    return scfilter, dcfilter
end

for (Xwt) in (:dwt!, :wpt!)
@eval begin
    # pseudo "inplace" by copying
    function $Xwt{T<:FloatingPoint}(x::AbstractArray{T}, filter::OrthoFilter, L::Integer, fw::Bool)
        y = Array(T, size(x))
        $Xwt(y, x, filter, L, fw)
        copy!(x,y)
        return x
    end
end
end

# DWT
# 1-D
# writes to y
function dwt!{T<:FloatingPoint}(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, L::Integer, fw::Bool)
    si = Array(T, length(filter)-1)       # tmp filter vector
    scfilter, dcfilter = makereverseqmf(filter.qmf, fw, T)

    dwt!(y, x, filter, L, fw, dcfilter, scfilter, si)
    return y
end
function dwt!{T<:FloatingPoint}(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, L::Integer, fw::Bool, dcfilter::Vector{T}, scfilter::Vector{T}, si::Vector{T}, snew::Vector{T} = Array(T, ifelse(L>1, length(x)>>1, 0)))
    n = length(x)
    J = nscales(n)
    @assert size(x) == size(y)
    @assert isdyadic(y)
    @assert 0 <= L <= J
    is(y,x) && error("input vector is output vector")

    L == 0 && return copy!(y,x)     # do nothing
    L > 1 && length(snew) != n>>1 && error("length of snew incorrect")
    s = x                           # s is currect scaling coefs location
    filtlen = length(filter)

    if fw
        jrange = (J-1):-1:(J-L)
    else
        jrange = (J-L):(J-1)
    end
    for j in jrange
        if fw
            # detail coefficients
            filtdown!(dcfilter, si, y, detailindex(j,1), detailn(j), s, 1,-filtlen+1, true)
            # scaling coefficients
            filtdown!(scfilter, si, y,                1, detailn(j), s, 1, 0, false)
        else
            # scaling coefficients
            filtup!(false, scfilter, si, y, 1, detailn(j+1), s, 1, -filtlen+1, false)
            # detail coefficients
            filtup!(true,  dcfilter, si, y, 1, detailn(j+1), x, detailindex(j,1), 0, true)
        end
        # if not final iteration: copy to tmp location
        fw  && j != jrange[end] && copy!(snew,1,y,1,detailn(j))
        !fw && j != jrange[end] && copy!(snew,1,y,1,detailn(j+1))
        L > 1 && (s = snew)
    end
    return y
end
function unsafe_dwt1level!{T<:FloatingPoint}(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, fw::Bool, dcfilter::Vector{T}, scfilter::Vector{T}, si::Vector{T})
    n = length(x)
    j = nscales(n) - 1
    filtlen = length(filter)

    if fw
        # detail coefficients
        filtdown!(dcfilter, si, y, detailindex(j,1), detailn(j), x, 1,-filtlen+1, true)
        # scaling coefficients
        filtdown!(scfilter, si, y,                1, detailn(j), x, 1, 0, false)
    else
        # scaling coefficients
        filtup!(false, scfilter, si, y, 1, detailn(j+1), x, 1, -filtlen+1, false)
        # detail coefficients
        filtup!(true,  dcfilter, si, y, 1, detailn(j+1), x, detailindex(j,1), 0, true)
    end
    return y
end

# 2-D
# writes to y
function dwt!{T<:FloatingPoint}(y::Matrix{T}, x::AbstractMatrix{T}, filter::OrthoFilter, L::Integer, fw::Bool)
    n = size(x,1)
    si = Array(T, length(filter)-1)       # tmp filter vector
    tmpvec = Array(T,n<<1)             # tmp storage vector
    scfilter, dcfilter = makereverseqmf(filter.qmf, fw, T)

    dwt!(y, x, filter, L, fw, dcfilter, scfilter, si, tmpvec)
    return y
end
function dwt!{T<:FloatingPoint}(y::Matrix{T}, x::AbstractMatrix{T}, filter::OrthoFilter, L::Integer, fw::Bool, dcfilter::Vector{T}, scfilter::Vector{T}, si::Vector{T}, tmpvec::Vector{T})

    n = size(x,1)
    J = nscales(n)
    @assert size(x) == size(y)
    @assert iscube(y)
    @assert isdyadic(y)
    @assert 0 <= L <= J
    @assert length(tmpvec) >= n<<1
    is(y,x) && error("input matrix is output matrix")

    L == 0 && return copy!(y,x)        # do nothing
    xs = n
    #s = x

    if fw
        jrange = (J-1):-1:(J-L)
        nsub = n
    else
        jrange = (J-L):(J-1)
        nsub = int(2^(J-L+1))
        copy!(y,x)
    end

    for j in jrange
        tmpsub = unsafe_vectorslice(tmpvec, 1, nsub)
        tmpsub2 = unsafe_vectorslice(tmpvec, nsub+1, nsub)
        if fw
            # rows
            for i=1:nsub
                xi = i
                if j != jrange[1]
                    stridedcopy!(tmpsub, y, xi, xs, nsub)
                else  # use x in first iteration
                    stridedcopy!(tmpsub, x, xi, xs, nsub)
                end
                unsafe_dwt1level!(tmpsub2, tmpsub, filter, fw, dcfilter, scfilter, si)
                stridedcopy!(y, xi, xs, tmpsub2, nsub)
            end
            # columns
            for i=1:nsub
                xi = 1+(i-1)*n
                ya = unsafe_vectorslice(y, xi, nsub)
                copy!(tmpsub,1,ya,1,nsub)
                unsafe_dwt1level!(ya, tmpsub, filter, fw, dcfilter, scfilter, si)
            end
        else
            # columns
            for i=1:nsub
                xi = 1+(i-1)*n
                ya = unsafe_vectorslice(y, xi, nsub)
                if j != jrange[1]
                    copy!(tmpsub,1,ya,1,nsub)
                    unsafe_dwt1level!(ya, tmpsub, filter, fw, dcfilter, scfilter, si)
                else  # use x in first iteration
                    #xsub = unsafe_vectorslice(x, xi, nsub)  # not safe for AbstractArray
                    copy!(tmpsub,1,x,xi,nsub)
                    unsafe_dwt1level!(ya, tmpsub, filter, fw, dcfilter, scfilter, si)
                end
            end
            # rows
            for i=1:nsub
                xi = i
                stridedcopy!(tmpsub, y, xi, xs, nsub)
                unsafe_dwt1level!(tmpsub2, tmpsub, filter, fw, dcfilter, scfilter, si)
                stridedcopy!(y, xi, xs, tmpsub2, nsub)
            end

        end

        fw  && (nsub = nsub>>1)
        !fw && (nsub = nsub<<1)
    end
    return y
end

# WPT
# 1-D
# writes to y
function wpt!{T<:FloatingPoint}(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, L::Integer, fw::Bool)
    wpt!(y, x, filter, maketree(length(y), L, :full), fw)
end
function wpt!{T<:FloatingPoint}(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, tree::BitVector, fw::Bool)
    si = Array(T, length(filter)-1)
    ns = ifelse(fw, length(x)>>1, length(x))
    snew = Array(T, ns)
    scfilter, dcfilter = makereverseqmf(filter.qmf, fw, T)

    wpt!(y, x, filter, tree, fw, dcfilter, scfilter, si, snew)
    return y
end
function wpt!{T<:FloatingPoint}(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, tree::BitVector, fw::Bool, dcfilter::Vector{T}, scfilter::Vector{T}, si::Vector{T}, snew::Vector{T})
    n = length(x)
    J = nscales(n)
    @assert size(x) == size(y)
    @assert isdyadic(y)
    @assert isvalidtree(y, tree)
    is(y, x) && error("input vector is output vector")
    tree[1] || return copy!(y,x)

    @assert length(snew) >= ifelse(fw, length(x)>>1, length(x))
    first = true
    L = J
    while L > 0
        ix = 1
        k = 1
        fw  && (Lfw = J-L)
        !fw && (Lfw = L-1)
        nj = detailn(tl2level(n, Lfw))
        treeind = 2^(Lfw)-1
        dx = unsafe_vectorslice(snew, 1, nj)

        while ix <= n
            if tree[treeind+k]
                dy = unsafe_vectorslice(y, ix, nj)
                if first
                    dx = unsafe_vectorslice(x, ix, nj)
                else
                    copy!(dx, dy)
                end
                unsafe_dwt1level!(dy, dx, filter, fw, dcfilter, scfilter, si)
            elseif first
                dy = unsafe_vectorslice(y, ix, nj)
                dx = unsafe_vectorslice(x, ix, nj)
                copy!(dy, dx)
            end
            ix += nj
            k += 1
        end
        L -= 1
        first = false
    end

    return y
end

macro filtermainloop(si, silen, b, val)
    quote
        @inbounds for j=2:$(esc(silen))
            $(esc(si))[j-1] = $(esc(si))[j] + $(esc(b))[j]*$(esc(val))
        end
        @inbounds si[$(esc(silen))] = $(esc(b))[$(esc(silen))+1]*$(esc(val))
    end
end
macro filtermainloopzero(si, silen)
    quote
        @inbounds for j=2:$(esc(silen))
            $(esc(si))[j-1] = $(esc(si))[j]
        end
        @inbounds si[$(esc(silen))] = 0.0
    end
end

# periodic filter and downsampling (by 2)
# apply f to x[ix:ix+nx-1] and write to out[iout:iout+nout-1]
# f : filter
# si: tmp array of length 1 less than f
# out : result gets written to out[iout:iout+nout-1]
# x : filter convolved with x[ix:ix+nx-1], where nx=nout*2 (shifted by shift)
# ss : shift downsampling
# based on Base.filt
function filtdown!{T<:FloatingPoint}(f::Vector{T}, si::Vector{T},
                              out::AbstractVector{T},  iout::Integer, nout::Integer,
                              x::AbstractVector{T}, ix::Integer, shift::Integer=0, ss::Bool=false)
    nx = nout<<1
    silen = length(si)
    flen = length(f)
    @assert length(x) >= ix+nx-1
    @assert length(out) >= iout+nout-1
    @assert (nx-1)>>1 + iout <= length(out)
    @assert silen == flen-1
    @assert shift <= 0

    fill!(si,0.0)
    istart = flen + int(ss)
    dsshift = (flen%2 + int(ss))%2  # is flen odd, and shift downsampling

    rout1, rin, rout2 = splitdownrangeper(istart, ix, nx, shift)
    @inbounds begin
        for i in rout1  # rout1 assumed to be in 1:istart-1
            # periodic in the range [ix:ix+nx-1]
            xatind = x[mod(i-1+shift, nx) + ix]
            @filtermainloop(si, silen, f, xatind)
        end
        # rin is inbounds in [ix:ix+nx-1]
        ixsh = -1 + shift + ix
        for i in rin
            xatind = x[i + ixsh]
            # take every other value after istart
            if (i+dsshift)%2 == 0 && i >= istart
                out[(i-istart)>>1 + iout] = si[1] + f[1]*xatind
            end
            @filtermainloop(si, silen, f, xatind)
        end
        for i in rout2
            # periodic in the range [ix:ix+nx-1]
            xatind = x[mod(i-1+shift, nx) + ix]
            # take every other value after istart
            if (i+dsshift)%2 == 0 && i >= istart
                out[(i-istart)>>1 + iout] = si[1] + f[1]*xatind
            end
            @filtermainloop(si, silen, f, xatind)
        end
    end

    return nothing
end
# find part of range which is inbounds for [ix:ix+nx-1] (ixsh = -1 + shift + ix)
# where mod(i-1+shift, nx) + ix == i + ixsh
function splitdownrangeper(istart, ix, nx, shift)
    inxi = 0
    ixsh = -1 + shift + ix
    if mod(shift, nx) + ix == 1 + ixsh   # shift likely 0
        inxi = 1
        iend = nx-1
        while mod(iend - 1 + shift, nx) + ix != iend + ixsh
            iend -= 1
        end
        return (0:-1, inxi:iend, (iend+1):(nx-1+istart))
    elseif mod(istart - 1 + shift, nx) + ix == istart + ixsh
        inxi = istart
        iend = nx-1+istart
        while mod(iend - 1 + shift, nx) + ix != iend + ixsh
            iend -= 1
        end
        return (1:inxi-1, inxi:iend, (iend+1):(nx-1+istart))
    else
        return (0:-1, 0:-1, 1:(nx-1+istart))
    end
end

# periodic filter and upsampling (by 2)
# apply f to x[ix:ix+nx-1] upsampled and write to out[iout:iout+nout-1]
# f : filter
# si: tmp array of length 1 less than f
# out : result gets written to out[iout:iout+nout-1]
# x : filter convolved with x[ix:ix+nx-1] upsampled, where nout==nx*2 (then shifted by shift)
# ss : shift upsampling
# based on Base.filt
function filtup!{T<:FloatingPoint}(add2out::Bool, f::Vector{T}, si::Vector{T},
                              out::AbstractVector{T},  iout::Integer, nout::Integer,
                              x::AbstractVector{T}, ix::Integer, shift::Integer=0, ss::Bool=false)
    nx = nout>>1
    silen = length(si)
    flen = length(f)
    @assert length(x) >= ix+nx-1                # check array size
    @assert length(out) >= iout+nout-1          # check array size
    @assert (nx<<1-1)>>1 + iout <= length(out)  # max for out index
    @assert silen == flen-1
    @assert shift <= 0

    fill!(si,0.0)
    istart = flen - shift%2
    dsshift = int(ss)%2  # shift upsampling

    rout1, rin, rout2 = splituprangeper(istart, ix, nx, nout, shift)
    @inbounds begin
        xatind = 0.0
        for i in rout1  # rout1 assumed to be in 1:istart-1
            if (i+dsshift)%2 == 0
                @filtermainloopzero(si, silen)
            else
                # periodic in the range [ix:ix+nx-1]
                xindex = mod((i-1)>>1+shift>>1,nx) + ix    #(i-1)>>1 increm. every other
                xatind = x[xindex]
                @filtermainloop(si, silen, f, xatind)
            end
        end
        ixsh = shift>>1 + ix
        for i in rin
            if (i+dsshift)%2 == 0
                xatind = 0.0
            else
                xindex = (i-1)>>1 + ixsh
                xatind = x[xindex]
            end
            if i >= istart
                if add2out
                    out[(i-istart) + iout] += si[1] + f[1]*xatind
                else
                    out[(i-istart) + iout] = si[1] + f[1]*xatind
                end
            end
            if (i+dsshift)%2 == 0
                @filtermainloopzero(si, silen)
            else
                @filtermainloop(si, silen, f, xatind)
            end
        end
        for i in rout2
            if (i+dsshift)%2 == 0
                xatind = 0.0
            else
                xindex = mod((i-1)>>1+shift>>1,nx) + ix
                xatind = x[xindex]
            end
            if i >= istart
                if add2out
                    out[(i-istart) + iout] += si[1] + f[1]*xatind
                else
                    out[(i-istart) + iout] = si[1] + f[1]*xatind
                end
            end
            if (i+dsshift)%2 == 0
                @filtermainloopzero(si, silen)
            else
                @filtermainloop(si, silen, f, xatind)
            end
        end
    end

    return nothing
end
# find part of range which is inbounds for [ix:ix+nx-1] (ixsh = shift>>1 + ix)
# where mod((i-1)>>1+shift>>1, nx) + ix == (i-1)>>1 + ixsh
function splituprangeper(istart, ix, nx, nout, shift)
    inxi = 0
    ixsh = shift>>1 + ix
    if mod(shift>>1, nx) + ix == ixsh   # shift likely 0
        inxi = 1
        iend = nout-1
        while mod((iend-1)>>1+shift>>1, nx) + ix != (iend-1)>>1 + ixsh
            iend -= 1
        end
        return (0:-1, inxi:iend, (iend+1):(nout-1+istart))
    elseif mod((istart-1)>>1+shift>>1, nx) + ix == (istart-1)>>1 + ixsh
        inxi = istart
        iend = nout-1+istart
        while mod((iend-1)>>1+shift>>1, nx) + ix != (iend-1)>>1 + ixsh
            iend -= 1
        end
        return (1:inxi-1, inxi:iend, (iend+1):(nout-1+istart))
    else
        return (0:-1, 0:-1, 1:(nout-1+istart))
    end
end


text:
"""
Handles Ajax communication logic.
"""
module WebThreads

import HTTP, Distributed, Logging, Dates
import Genie, Genie.Renderer

const MESSAGE_QUEUE = Dict{UInt,Vector{String}}()

const ClientId = UInt # session id
const ChannelName = String

mutable struct ChannelClient
  client::UInt
  channels::Vector{ChannelName}
  last_active::Dates.DateTime
end

const ChannelClientsCollection = Dict{ClientId,ChannelClient} # { uint => { :client => ws, :channels => ["foo", "bar", "baz"] } }
const ChannelSubscriptionsCollection = Dict{ChannelName,Vector{ClientId}}  # { "foo" => ["4", "12"] }
const MessagePayload = Union{Nothing,Dict}

mutable struct ChannelMessage
  channel::ChannelName
  client::ClientId
  message::String
  payload::MessagePayload
end

const CLIENTS = ChannelClientsCollection()
const SUBSCRIPTIONS = ChannelSubscriptionsCollection()

clients() = collect(values(CLIENTS))
subscriptions() = SUBSCRIPTIONS
webthreads() = map(c -> c.client, clients())
channels() = collect(keys(SUBSCRIPTIONS))

function connected_clients(channel::ChannelName) :: Vector{ChannelClient}
  clients = ChannelClient[]
  for client_id in SUBSCRIPTIONS[channel]
    ((Dates.now() - CLIENTS[client_id].last_active) <= Genie.config.webthreads_connection_threshold) && push!(clients, CLIENTS[client_id])
  end

  clients
end
function connected_clients() :: Vector{ChannelClient}
  clients = ChannelClient[]
  for ch in channels()
    clients = vcat(clients, connected_clients(ch))
  end

  clients
end

function disconnected_clients(channel::ChannelName) :: Vector{ChannelClient}
  clients = ChannelClient[]
  for client_id in SUBSCRIPTIONS[channel]
    ((Dates.now() - CLIENTS[client_id].last_active) > Genie.config.webthreads_connection_threshold) && push!(clients, CLIENTS[client_id])
  end

  clients
end
function disconnected_clients() :: Vector{ChannelClient}
  clients = ChannelClient[]
  for ch in channels()
    clients = vcat(clients, disconnected_clients(ch))
  end

  clients
end

"""
Subscribes a web thread client `wt` to `channel`.
"""
function subscribe(wt::UInt, channel::ChannelName) :: ChannelClientsCollection
  if haskey(CLIENTS, wt)
    in(channel, CLIENTS[wt].channels) || push!(CLIENTS[wt].channels, channel)
  else
    CLIENTS[wt] = ChannelClient(wt, ChannelName[channel], Dates.now())
  end

  push_subscription(wt, channel)

  haskey(MESSAGE_QUEUE, wt) || (MESSAGE_QUEUE[wt] = String[])

  CLIENTS
end

"""
Unsubscribes a web socket client `wt` from `channel`.
"""
function unsubscribe(wt::UInt, channel::ChannelName) :: ChannelClientsCollection
  haskey(CLIENTS, wt) && deleteat!(CLIENTS[wt].channels, CLIENTS[wt].channels .== channel)
  pop_subscription(wt, channel)

  CLIENTS
end
function unsubscribe(channel_client::ChannelClient, channel::ChannelName) :: ChannelClientsCollection
  unsubscribe(channel_client.client, channel)
end

"""
Unsubscribes a web socket client `wt` from all the channels.
"""
function unsubscribe_client(wt::UInt) :: ChannelClientsCollection
  if haskey(CLIENTS, wt)
    for channel_id in CLIENTS[wt].channels
      pop_subscription(wt, channel_id)
    end

    delete!(CLIENTS, wt)
  end

  CLIENTS
end
function unsubscribe_client(channel_client::ChannelClient) :: ChannelClientsCollection
  unsubscribe_client(channel_client.client)

  CLIENTS
end

"""
unsubscribe_disconnected_clients() :: ChannelClientsCollection

Unsubscribes clients which are no longer connected.
"""
function unsubscribe_disconnected_clients() :: ChannelClientsCollection
  for channel_client in disconnected_clients()
    unsubscribe_client(channel_client)
  end

  CLIENTS
end
function unsubscribe_disconnected_clients(channel::ChannelName) :: ChannelClientsCollection
  for channel_client in disconnected_clients(channel)
    unsubscribe(channel_client, channel)
  end

  CLIENTS
end

function unsubscribe_clients()
  empty!(CLIENTS)
  empty!(SUBSCRIPTIONS)
end

function timestamp_client(client_id::ClientId) :: Nothing
  haskey(CLIENTS, client_id) && (CLIENTS[client_id].last_active = Dates.now())

  nothing
end

"""
Adds a new subscription for `client` to `channel`.
"""
function push_subscription(client_id::ClientId, channel::ChannelName) :: ChannelSubscriptionsCollection
  if haskey(SUBSCRIPTIONS, channel)
    ! in(client_id, SUBSCRIPTIONS[channel]) && push!(SUBSCRIPTIONS[channel], client_id)
  else
    SUBSCRIPTIONS[channel] = ClientId[client_id]
  end

  timestamp_client(client_id)

  SUBSCRIPTIONS
end
function push_subscription(channel_client::ChannelClient, channel::ChannelName) :: ChannelSubscriptionsCollection
  push_subscription(channel_client.client, channel)
end

"""
Removes the subscription of `client` to `channel`.
"""
function pop_subscription(client::ClientId, channel::ChannelName) :: ChannelSubscriptionsCollection
  if haskey(SUBSCRIPTIONS, channel)
    filter!(SUBSCRIPTIONS[channel]) do (client_id)
      client_id != client
    end
    isempty(SUBSCRIPTIONS[channel]) && delete!(SUBSCRIPTIONS, channel)
  end

  timestamp_client(client)

  SUBSCRIPTIONS
end
function pop_subscription(channel_client::ChannelClient, channel::ChannelName) :: ChannelSubscriptionsCollection
  pop_subscription(channel_client.client, channel)
end

"""
Removes all subscriptions of `client`.
"""
function pop_subscription(channel::ChannelName) :: ChannelSubscriptionsCollection
  if haskey(SUBSCRIPTIONS, channel)
    delete!(SUBSCRIPTIONS, channel)
  end

  SUBSCRIPTIONS
end

"""
Pushes `msg` (and `payload`) to all the clients subscribed to the channels in `channels`.
"""
function broadcast(channels::Union{ChannelName,Vector{ChannelName}}, msg::String;
                    except::Union{HTTP.WebSockets.WebSocket,Nothing,UInt} = nothing) :: Bool
  isa(channels, Array) || (channels = ChannelName[channels])

  for channel in channels
    haskey(SUBSCRIPTIONS, channel) || throw(Genie.WebChannels.ChannelNotFoundException(channel))

    for client in SUBSCRIPTIONS[channel]
      except !== nothing && client == except && continue

      try
        message(client, msg)
      catch ex
        @error ex
      end
    end
  end

  true
end
function broadcast(channels::Union{ChannelName,Vector{ChannelName}}, msg::String, payload::Dict) :: Bool
  isa(channels, Array) || (channels = [channels])

  for channel in channels
    haskey(SUBSCRIPTIONS, channel) || throw(Genie.WebChannels.ChannelNotFoundException(channel))

    for client in SUBSCRIPTIONS[channel]
      try
        message(client, ChannelMessage(channel, client, msg, payload) |> Renderer.Json.JSONParser.json)
      catch ex
        @error ex
      end
    end
  end

  true
end

"""
Pushes `msg` (and `payload`) to all the clients subscribed to all the channels.
"""
function broadcast(msg::String, payload::Union{Dict,Nothing} = nothing) :: Bool
  payload === nothing ?
    broadcast(collect(keys(SUBSCRIPTIONS)), msg) :
    broadcast(collect(keys(SUBSCRIPTIONS)), msg, payload)
end

"""
Pushes `msg` (and `payload`) to `channel`.
"""
function message(channel::ChannelName, msg::String, payload::Union{Dict,Nothing} = nothing) :: Bool
  payload === nothing ?
    broadcast(channel, msg) :
    broadcast(channel, msg, payload)
end

"""
Writes `msg` to message queue for `client`.
"""
function message(wt::UInt, msg::String)
  push!(MESSAGE_QUEUE[wt], msg)
end
function message(client::ChannelClient, msg::String)
  message(client.client, msg)
end

function pull(wt::UInt, channel::ChannelName)
  output = ""
  if haskey(MESSAGE_QUEUE, wt) && ! isempty(MESSAGE_QUEUE[wt])
    output = MESSAGE_QUEUE[wt] |> Renderer.Json.JSONParser.json
    empty!(MESSAGE_QUEUE[wt])
  end

  timestamp_client(wt)

  output
end

function push(wt::UInt, channel::ChannelName, message::String)
  timestamp_client(wt)

  Genie.Router.route_ws_request(Genie.Router.@params(Genie.PARAMS_REQUEST_KEY), message, wt)
end

end

ext: jl


text:
module MPI

using Libdl, Serialization
using Requires
using DocStringExtensions

export mpiexec, UBuffer, VBuffer

function serialize(x)
    s = IOBuffer()
    Serialization.serialize(s, x)
    take!(s)
end

function deserialize(x)
    s = IOBuffer(x)
    Serialization.deserialize(s)
end

primitive type SentinelPtr Sys.WORD_SIZE
end

primitive type MPIPtr Sys.WORD_SIZE
end
@assert sizeof(MPIPtr) == sizeof(Ptr{Cvoid})
Base.cconvert(::Type{MPIPtr}, x::SentinelPtr) = x
Base.unsafe_convert(::Type{MPIPtr}, x::SentinelPtr) = reinterpret(MPIPtr, x)

function _doc_external(fname)
"""
- `$fname` man page: [OpenMPI](https://www.open-mpi.org/doc/current/man3/$fname.3.php), [MPICH](https://www.mpich.org/static/docs/latest/www3/$fname.html)
"""
end

import MPIPreferences

if MPIPreferences.binary == "MPICH_jll"
    import MPICH_jll: libmpi, mpiexec
    const libmpiconstants = nothing
elseif MPIPreferences.binary == "OpenMPI_jll"
    import OpenMPI_jll: libmpi, mpiexec
    const libmpiconstants = nothing
elseif MPIPreferences.binary == "MicrosoftMPI_jll"
    import MicrosoftMPI_jll: libmpi, mpiexec
    const libmpiconstants = nothing
elseif MPIPreferences.binary == "MPItrampoline_jll"
    import MPItrampoline_jll: MPItrampoline_jll, libmpi, mpiexec
    const libmpiconstants = MPItrampoline_jll.libload_time_mpi_constants_path
elseif MPIPreferences.binary == "system"
    import MPIPreferences.System: libmpi, mpiexec
    const libmpiconstants = nothing
else
    error("Unknown MPI binary: $(MPIPreferences.binary)")
end

include("consts/consts.jl")
using .Consts

# These functions are run after reading the values of the constants above)
const _mpi_load_time_hooks = Any[]
const _finished_loading = Ref(false)
function add_load_time_hook!(f)
    @assert !_finished_loading[]
    push!(_mpi_load_time_hooks, f)
end
function run_load_time_hooks()
    @assert !_finished_loading[]
    _finished_loading[] = true
    for hook in _mpi_load_time_hooks
        hook()
    end
    empty!(_mpi_load_time_hooks)
    nothing
end

include("implementations.jl")
include("error.jl")
include("info.jl")
include("group.jl")
include("comm.jl")
include("environment.jl")
include("datatypes.jl")
include("buffers.jl")
include("operators.jl")
include("nonblocking.jl")
include("pointtopoint.jl")
include("collective.jl")
include("topology.jl")
include("onesided.jl")
include("io.jl")
include("errhandler.jl")
include("mpiexec_wrapper.jl")

include("deprecated.jl")

function __init__()
    @static if Sys.isunix()
        # dlopen the MPI library before any ccall:
        # - RTLD_GLOBAL is required for Open MPI
        #   https://www.open-mpi.org/community/lists/users/2010/04/12803.php
        # - also allows us to ccall global symbols, which enables profilers
        #   which use LD_PRELOAD
        # - don't use RTLD_DEEPBIND; this leads to issues with multiple MPI
        #   libraries:
        #   https://github.com/JuliaParallel/MPI.jl/pull/109
        #   https://github.com/JuliaParallel/MPI.jl/issues/587
        Libdl.dlopen(libmpi, Libdl.RTLD_LAZY | Libdl.RTLD_GLOBAL)
    end

    # Needs to be called after `dlopen`. Use `invokelatest` so that `cglobal`
    # calls don't trigger early `dlopen`-ing of the library.
    Base.invokelatest(Consts.init_consts)

    # disable UCX memory cache, since it doesn't work correctly
    # https://github.com/openucx/ucx/issues/5061
    if !haskey(ENV, "UCX_MEMTYPE_CACHE")
        ENV["UCX_MEMTYPE_CACHE"] = "no"
    end

    # Julia multithreading uses SIGSEGV to sync threads
    # https://docs.julialang.org/en/v1/devdocs/debuggingtips/#Dealing-with-signals-1
    # By default, UCX will error if this occurs (issue #337)
    if !haskey(ENV, "UCX_ERROR_SIGNALS")
        # default is "SIGILL,SIGSEGV,SIGBUS,SIGFPE"
        ENV["UCX_ERROR_SIGNALS"] = "SIGILL,SIGBUS,SIGFPE"
    end

    if MPIPreferences.binary == "MPItrampoline_jll" && !haskey(ENV, "MPITRAMPOLINE_MPIEXEC")
        ENV["MPITRAMPOLINE_MPIEXEC"] = MPItrampoline_jll.mpich_mpiexec_path
    end

    run_load_time_hooks()

    @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" include("cuda.jl")
end

end


text:
#=
 = This module returns a dictionary of symbols and alphanumeric values, as given by a single
 = ASCII line from an International Maritime Meteorological Archive file (IMMA Release 2.5;
 = following http://rda.ucar.edu/datasets/ds540.0/docs/R2.5-imma_short.pdf).  This format
 = consists of a core entry and up to six attachments, some of which are still in the process
 = of being standardized.  Note that the dictionary returns only defined values (that is, the
 = original string if a scale parameter is zero, a scaled float if the scale is positive, and
 = a float converted from a different base if the scale is negative) and for the time being,
 = attachments 5 and 6 are not parsed (at least until 5 is formalized) - RD March 2016, May 2017
 =#

module ICOADSDict                                       # first define the values contained in the core
export imma                                             # and attachments by symbol, location, and scale

const symb0 = [:YR, :MO, :DY,  :HR, :LAT, :LON, :IM, :ATTC, :TI, :LI, :DS, :VS, :NID, :II, :ID, :C1, :DI, :D, :WI,  :W, :VI, :VV, :WW, :W1, :SLP, :A, :PPP, :IT, :AT, :WBTI, :WBT, :DPTI, :DPT, :SI, :SST, :N, :NH, :CL, :HI, :H, :CM, :CH, :WD, :WP, :WH, :SD, :SP, :SH]
const leng0 = [  4,   2,   2,    4,    5,    6,   2,     1,   1,   1,   1,   1,    2,   2,   9,   2,   1,  3,   1,   3,   1,   2,   2,   1,    5,  1,    3,   1,   4,     1,    4,     1,    4,   2,    4,  1,   1,   1,   1,  1,   1,   1,   2,   2,   2,   2,   2,   2]
const scal0 = [  1,   1,   1, 0.01, 0.01, 0.01,   1,     1,   1,   1,   1,   1,    1,   1,   0,   0,   1,  1,   1, 0.1,   1,   1,   1,   1,  0.1,  1,  0.1,   1, 0.1,     1,  0.1,     1,  0.1,   1,  0.1,  1,   1,   0,   1,  0,   0,   0,   1,   1,   1,   1,   1,   1]

const symb1 = [:ATTI1, :ATTL1, :BSI, :B10, :B1, :DCK, :SID, :PT, :DUPS, :DUPC, :TC, :PB, :WX, :SX, :C2, :SQZ, :SQA, :AQZ, :AQA, :UQZ, :UQA, :VQZ, :VQA, :PQZ, :PQA, :DQZ, :DQA, :ND, :SF, :AF, :UF, :VF, :PF, :RF, :ZNC, :WNC, :BNC, :XNC, :YNC, :PNC, :ANC, :GNC, :DNC, :SNC, :CNC, :ENC, :FNC, :TNC, :QCE, :LZ, :QCZ]
const leng1 = [     2,      2,    1,    3,   2,    3,    3,   2,     2,     1,   1,   1,   1,   1,   2,    1,    1,    1,    1,    1,    1,    1,    1,    1,    1,    1,    1,   1,   1,   1,   1,   1,   1,   1,    1,    1,    1,    1,    1,    1,    1,    1,    1,    1,    1,    1,    1,    1,    2,   1,    2]
const scal1 = [     1,      1,    0,    1,   1,    1,    1,   1,     1,     1,   1,   1,   1,   1,   1,  -36,  -36,  -36,  -36,  -36,  -36,  -36,  -36,  -36,  -36,  -36,  -36,   1, -36, -36, -36, -36, -36, -36,  -36,  -36,  -36,  -36,  -36,  -36,  -36,  -36,  -36,  -36,  -36,  -36,  -36,  -36,    1,   1,    1]

const symb2 = [:ATTI2, :ATTL2, :OS, :OP, :FM, :IX, :W2, :SGN, :SGT, :SGH, :WMI, :SD2, :SP2, :SH2, :IS, :ES, :RS, :IC1, :IC2, :IC3, :IC4, :IC5, :IR, :RRR, :TR, :QCI, :QI1, :QI2, :QI3, :QI4, :QI5, :QI6, :QI7, :QI8, :QI9, :QI10, :QI11, :QI12, :QI13, :QI14, :QI15, :QI16, :QI17, :QI18, :QI19, :QI20, :QI21, :HDG, :COG, :SOG, :SLL, :SLHH, :RWD, :RWS]
const leng2 = [     2,      2,   1,   1,   2,   1,   1,    1,    1,    2,    1,    2,    2,    2,   1,   2,   1,    1,    1,    1,    1,    1,   1,    3,   1,    1,    1,    1,    1,    1,    1,    1,    1,    1,    1,     1,     1,     1,     1,     1,     1,     1,     1,     1,     1,     1,     1,    3,    3,    2,    2,     3,    3,    3]
const scal2 = [     1,      1,   1,   1,   1,   1,   1,    1,    0,    1,    1,    1,    1,    1,   1,   1,   1,    0,    0,    0,    0,    0,   1,    1,   1,    1,    1,    1,    1,    1,    1,    1,    1,    1,    1,     1,     1,     1,     1,     1,     1,     1,     1,     1,     1,     1,     1,    1,    1,    1,    1,     1,    1,    1]

const symb3 = [:ATTI3, :ATTL3, :CCCC, :BUID, :BMP, :BSWU, :SWU, :BSWV, :SWV, :BSAT, :BSRH, :SRH, :SIX, :BSST, :MST, :MSH, :BY, :BM, :BD, :BH, :BFL]
const leng3 = [     2,      2,     4,     6,    5,     4,    4,     4,    4,     4,     3,    3,    1,     4,    1,    3,   4,   2,   2,   2,    2]
const scal3 = [     1,      1,     0,     0,  0.1,   0.1,  0.1,   0.1,  0.1,   0.1,     1,    1,    1,   0.1,    1,    1,   1,   1,   1,   1,    1]

const symb4 = [:ATTI4, :ATTL4, :C1M, :OPM, :KOV, :COR, :TOB, :TOT, :EOT, :LOT, :TOH, :EOH, :SIM, :LOV, :DOS, :HOP, :HOT, :HOB, :HOA, :SMF, :SME, :SMV]
const leng4 = [     2,      2,    2,    2,    2,    2,    3,    3,    2,    2,    1,    2,    3,    3,    2,    3,    3,    3,    3,    5,    5,    2]
const scal4 = [     1,      1,    0,    1,    0,    0,    0,    0,    0,    0,    0,    0,    0,    1,    1,    1,    1,    1,    1,    1,    1,    1]

const symb5 = [:ATTI5, :ATTL5, :WFI, :WF, :XWI, :XW, :XDI, :XD, :SLPI, :TAI, :TA, :XNI, :XN, :OTHERSTUFF5]
const leng5 = [     2,      2,    1,   2,    1,   3,    1,   2,     1,    1,   4,    1,   2,           71]
const scal5 = [     1,      1,    0,   1,    0,   1,    0,   0,     0,    0,   1,    0,   0,            0]

const symb6 = [:ATTI6, :ATTL6, :ATTE, :OTHERSTUFF6]
const leng6 = [     2,      2,     1,          999]
const scal6 = [     1,      1,     0,            0]

function imma(line::AbstractString)
  vals = Dict{Symbol, Any}()                                                  # initialize a "vals" dictionary
  linlen = length(line)

  subsec = 0                                                                  # then loop through the core and
  substa = 1                                                                  # all available attachments until
  while substa < linlen                                                       # we're past the end of the line
    if substa > 1                                                             # (sta/end are section start/end)
      subsec = parse(Int, line[substa:substa+1])
    end

    symb = eval(Symbol("symb", subsec))                                       # use the relevant section arrays
    leng = eval(Symbol("leng", subsec))
    scal = eval(Symbol("scal", subsec))
    sumlen = sum(leng)
    subend = substa + sumlen - 1

    if subsec < 5                                                             # ignore historical/supplemental
      if subend <= linlen                                                     # isolate this section and pad if
        sublin =      line[substa:subend]                                     # at the end of a truncated line
      else
        sublin = rpad(line[substa:   end], sumlen, ' ')
      end
#     println(subsec, " :", sublin, ": ")

      index = [0; cumsum(leng)] .+ 1                                          # add any valid entries from each
      for (a, sym) in enumerate(symb)                                         # section and return the dictionary
        tmpval = sublin[index[a]:index[a]+leng[a]-1]
        if tmpval != " "^leng[a]
          if scal[a] > 0
            vals[sym] = parse(Float64, tmpval) * scal[a]
          elseif scal[a] == 0
            vals[sym] = strip(tmpval)
          else
            vals[sym] = float(parse(Int, tmpval, base = 36))
          end
#         println(sym, " ", vals[sym])
        end
      end
    end
    substa = subend + 1
  end
  return vals
end

end                                                                           # end ICOADSDict module


text:
function Base.length(x::NASM)
  return length(fieldnames(typeof(x)))
end

function Base.getindex(x::NASM, i::Int64)
  return getproperty(x, fieldnames(typeof(x))[i])
end

@recipe function f(iter,
                   x::Vector{RiskParameters{T}}) where {T <: DiseaseStateSequence}
  xguide --> "Iteration"
  yguide --> "Value"
  legend --> :none
  seriestype := :line
  # label --> permutedims(_parameters(x[1]))
  iter, convert(Array{Float64, 2}, x[iter])
end

@recipe function f(iter, x::Vector{T}) where {T <: NASM}
  xguide --> "Iteration"
  yguide --> "Value"
  legend --> :none
  seriestype := :line
  # label --> permutedims(_parameters(x[1]))
  iter, convert(Array{Float64, 2}, x[iter])
end

@recipe function f(iter, p::Int64, x::MCMC{S, M}) where {S <: DiseaseStateSequence, M <: ILM}
  np = length(x.markov_chains[1].risk_parameters[1])
  sp = 0
  if M == PhyloILM
    sp += length(x.markov_chains[1].substitution_model[1])
  end
  xguide --> "Iteration"
  yguide --> "Value"
  guidefontsize --> 8
  seriestype --> :line
  if 0 <= p <= np
    for i in eachindex(x.markov_chains)
      @series begin
        label --> "Markov chain $i"
        iter, getindex.(x.markov_chains[i].risk_parameters[iter], p)
      end
    end
  elseif np < p <= (np + sp)
    for i in eachindex(x.markov_chains)
      @series begin
        label --> "Markov chain $i"
        iter, getindex.(x.markov_chains[i].substitution_model[iter], p-np)
      end
    end
  else
    error("Invalid parameter index")
  end
end

@recipe function f(iter, x::MCMC{S, M}) where {S <: DiseaseStateSequence, M <: ILM}
  np = length(x.markov_chains[1].risk_parameters[1])
  if M == PhyloILM
    np += length(x.markov_chains[1].substitution_model[1])
  end
  layout --> (np, 1)
  legend --> :none
  link --> :y
  for i = 1:np
    @series begin
      subplot := i
      iter, i, x
    end
  end
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, HammerScalpel

makedocs(modules = [HammerScalpel], sitename = "HammerScalpel.jl")

deploydocs(repo = "github.com/JulianStreyczek/HammerScalpel.jl.git", devbranch = "main")


text:
# Autogenerated wrapper script for Dierckx_jll for x86_64-w64-mingw32-libgfortran4
export libddierckx

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Dierckx")
JLLWrappers.@declare_library_product(libddierckx, "libddierckx.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libddierckx,
        "bin\\libddierckx.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function addSavingCSV_!(p::Program_,abm::Agent,Space::SimulationSpace,platform::String)

    addSavingRAM_!(p,abm,Space,platform)

    push!(p.declareVar.args,:(saveCounter=1))

    p.execInit = postwalk(x->@capture(x,push!(commRAM_,ob)) ? :(saveCSV(ob,string(saveFile,"_",saveCounter)); saveCounter += 1) : x, p.execInit)
    p.execInloop = postwalk(x->@capture(x,push!(commRAM_,ob)) ? :(saveCSV(ob,string(saveFile,"_",saveCounter)); saveCounter += 1) : x, p.execInloop)

    p.returning = quote end

    return
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
# Removal of LinearKernel
Base.@deprecate LinearKernel(a::Real=1, c::Real=1) PolynomialKernel(a, c, 1)

# Renaming of GammaRational to GammaRationalQuadratic
Base.@deprecate GammaRationalKernel GammaRationalQuadraticKernel

# Removal of PeriodicKernel
struct SineSquared <: PreMetric end
@inline base_aggregate(::SineSquared, s::T, x::T, y::T) where {T} = s + sin(x-y)^2
@inline isstationary(::SineSquared) = true

struct PeriodicKernel{T<:AbstractFloat} <: MercerKernel{T}
    α::T
    function PeriodicKernel{T}(α::Real) where {T<:AbstractFloat}
        Base.depwarn("PeriodicKernel will be removed in the next major release", :PeriodicKernel)
        @check_args(PeriodicKernel, α, α > zero(α), "α > 0")
        new{T}(α)
    end
end
PeriodicKernel(α::T₁ = 1.0) where {T₁<:Real} = PeriodicKernel{promote_float(T₁)}(α)

@inline basefunction(::PeriodicKernel) = SineSquared()

@inline function kappa(κ::PeriodicKernel{T}, z::T) where {T}
    return exp(-κ.α*z)
end


text:
struct Frame
    image::Image
    view::ImageView
    framebuffer::Framebuffer
end

struct FrameSynchronization
    image_acquired::Semaphore
    image_rendered::Semaphore
    has_rendered::Fence
end

FrameSynchronization(device) = FrameSynchronization(Semaphore(device), Semaphore(device), Fence(device; flags = FENCE_CREATE_SIGNALED_BIT))

struct FrameState
    device::Device
    swapchain::Ref{Created{SwapchainKHR,SwapchainCreateInfoKHR}}
    render_pass::RenderPass
    frames::Vector{Frame}
    current_frame::Ref{Frame}
    frame_count::Ref{Int}
    syncs::Dictionary{Frame,FrameSynchronization}
end

Vulkan.device(fs::FrameState) = fs.device

function FrameState(swapchain::Created{SwapchainKHR,SwapchainCreateInfoKHR}, render_pass)
    max_in_flight = info(swapchain).min_image_count
    fs = FrameState(device(render_pass), Ref(swapchain), render_pass, [], Ref{Frame}(), Ref(0), Dictionary{Frame,FrameSynchronization}())
    update!(fs)
end

function Vulkan.SurfaceCapabilitiesKHR(fs::FrameState)
    unwrap(get_physical_device_surface_capabilities_khr(device(fs).physical_device, info(fs.swapchain[]).surface))
end

function recreate_swapchain!(fs::FrameState, new_extent::Extent2D)
    swapchain = fs.swapchain[]
    swapchain_info = setproperties(info(swapchain), old_swapchain = handle(swapchain), image_extent = new_extent)
    swapchain_handle = unwrap(create_swapchain_khr(device(fs), swapchain_info))
    fs.swapchain[] = Created(swapchain_handle, swapchain_info)
    fs
end

function get_frames(fs::FrameState)
    swapchain = fs.swapchain[]
    extent = info(swapchain).image_extent
    map(unwrap(get_swapchain_images_khr(device(fs), swapchain))) do img
        view = ImageView(
            device(fs),
            img,
            IMAGE_VIEW_TYPE_2D,
            info(swapchain).image_format,
            ComponentMapping(fill(COMPONENT_SWIZZLE_IDENTITY, 4)...),
            ImageSubresourceRange(IMAGE_ASPECT_COLOR_BIT, 0, 1, 0, 1),
        )

        fb = Framebuffer(device(fs), fs.render_pass, [view], extent.width, extent.height, 1)
        Frame(img, view, fb)
    end
end

function update!(fs::FrameState)
    empty!(fs.frames)
    empty!(fs.syncs)

    # TODO: make fields of returned-only structs high-level in Vulkan.jl
    _extent = SurfaceCapabilitiesKHR(fs).current_extent.vks
    new_extent = Extent2D(_extent.width, _extent.height)
    if new_extent ≠ info(fs.swapchain[]).image_extent
        recreate_swapchain!(fs, new_extent)
    end

    append!(fs.frames, get_frames(fs))
    fs.current_frame[] = first(fs.frames)

    foreach(fs.frames) do frame
        insert!(fs.syncs, frame, FrameSynchronization(device(fs.render_pass)))
    end
    fs
end

"""
    command_buffers(frame_state, app)

Command buffers generated by the application.
"""
function command_buffers end

function next_frame!(fs::FrameState, dispatch::QueueDispatch, app)
    swapchain = fs.swapchain[]

    # acquire next image
    image_acquired = fs.syncs[fs.current_frame[]].image_acquired
    status = @timeit to "Acquire next image" acquire_next_image_khr(device(fs), swapchain, typemax(UInt64); semaphore = image_acquired)
    if !iserror(status)
        idx, result = unwrap(status)
        @assert result in (SUCCESS, SUBOPTIMAL_KHR) "$result: Could not retrieve next swapchain image"
        frame = fs.frames[idx + 1]
        sync = fs.syncs[frame]

        # submit rendering commands
        @timeit to "Submit rendering commands" begin
            cbuffs = @timeit to "Create command buffers" command_buffers(fs, app)
            img_acquired_info = SemaphoreSubmitInfoKHR(image_acquired, 0, 0; stage_mask = PIPELINE_STAGE_2_COLOR_ATTACHMENT_OUTPUT_BIT_KHR)
            img_rendered_info = SemaphoreSubmitInfoKHR(sync.has_rendered, 0, 0; stage_mask = PIPELINE_STAGE_2_COLOR_ATTACHMENT_OUTPUT_BIT_KHR)
            render_info = SubmitInfo2KHR([img_acquired_info], CommandBufferSubmitInfoKHR.(cbuffs, 0), [img_rendered_info])
            @timeit to "Wait for frame being rendered" wait_for_fences(device(fs), [sync.has_rendered], false, typemax(UInt64))
            reset_fences(device(fs), [sync.has_rendered])
            submit(dispatch, QUEUE_GRAPHICS_BIT, [render_info]; fence = sync.has_rendered)
        end

        # submit presentation commands
        @timeit to "Submit presentation commands" begin
            present_info = PresentInfoKHR([sync.image_rendered], [swapchain], [idx])
            unwrap(present(dispatch, present_info))
        end

        fs.frame_count[] += 1
        fs.current_frame[] = frame
    else
        err = unwrap_error(status)
        if err.code == ERROR_OUT_OF_DATE_KHR
            # recreate swapchain and start over
            @timeit to "Recreate swapchain" update!(fs)
            next_frame!(fs, rdr, app)
        end
    end
end

function wait_hasrendered(fs::FrameState)
    wait_for_fences(device(fs), getproperty.(fs.syncs, :has_rendered), true, typemax(UInt64))
end

"""
Encode a lifetime context that allows the deletion of staged resources once expired.
"""
abstract type LifeTimeContext end

"""
Work between frames is assumed to be independent,
and so is any resource attached to a particular piece of work.
"""
struct FrameContext <: LifeTimeContext
    state::FrameState
    pending::Vector{Any}
end

function next_frame!(ctx::FrameContext)
    next_frame!(ctx.state)
    foreach(ctx.pending) do resource
        finalize(resource)
    end
    empty!(ctx.pending)
end

function stage!(ctx::FrameContext, resource)
    push!(ctx.pending, resource)
end


text:
module Forecast

using CSV, Distributions, DataFrames, GZip, LinearAlgebra, Plots,
      RecipesBase, TimeSeries, Statistics

# types
export CCF, STL
# methods
export acf, ccf, d, loess, sma, stl, hma, hmaSymmetricWeights
# datasets
export co2

# source files

## ccf
include("ccf.jl")
include("acf.jl")

include("d.jl")

## stl
include("loess.jl")
include("sma.jl")
include("stl.jl")

## hma
include("hma.jl")

include("utils.jl")

include("datasets.jl")

## recipes
include("plotrecipes.jl")

"""
Collection of methods for Time Series analysis

Methods implemented:

    acf:        Auto-correlation or auto-covariance of univariate series.
    ccf:        Cross-correlation or cros-covariance of two univariate series.
    d:          Lagged differences of a given order for Vector, Array and TimeSeries.
    hma:        Henderson moving average filters.
    loess:      Locally estimated scatterplot smoothing.
    sma:        Simple moving average.
    stl:        Seasonal and Trend decomposition using loess.
"""
Forecast

end


text:
# This is test of the code for computing a transit model integrated over a time step.
# Computes derivatives over the timestep.
using PyPlot

if VERSION >= v"0.7"
  using Statistics
end

#include("../src/integrate_transit_cubature.jl")
#include("../src/integrate_transit_simpson_vec.jl")
include("../src/integrate_lightcurve.jl")
#include("integrate_transit_simpson_vec.jl")

# Test it out:

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

t1 = -1.5; t2 = 1.5; nt = 10000; dt = 0.3
t = zeros(nt)
t .= linearspace(t1,t2,nt)
#t = [t[1235]] ; nt =1
# The following compares different tolerances and maxdepths:
r = 0.1; b0 = 0.5; u_n = [0.3,0.3]; nu = length(u_n)
favg0 = zeros(5+nu,nt)
neval0= ones(Int64,nt)
favg1 = zeros(5+nu,nt)
neval1= zeros(Int64,nt)
depthmax1= zeros(Int64,nt)
favg2 = zeros(5+nu,nt)
neval2= zeros(Int64,nt)
depthmax2= zeros(Int64,nt)
favg3 = zeros(5+nu,nt)
neval3= zeros(Int64,nt)
depthmax3= zeros(Int64,nt)
trans = transit_init(r,b0,u_n,true)
param = [0.0,1.0,b0]   # [t_0,v,b_0]

# Now, carry out speed test:

function compute_lightcurve!(trans::Transit_Struct{T},param::Array{T,1},t::Array{T,1},favg0::Array{T,2},nt::Int64,neval::Array{Int64,1}) where {T <: Real}
@inbounds for i=1:nt
  if neval[i] == 0
    neval[i] = 1
  end
  for j=1:neval[i]
    trans.b = sqrt(param[3]^2+(param[2]*(t[i]-param[1]))^2)
    favg0[1,i] =  transit_poly!(trans)
    binv = inv(trans.b)
    favg0[2,i] =  trans.dfdrb[1]
    favg0[3,i] =  trans.dfdrb[2]*binv*param[2]^2*(param[1]-t[i])
    favg0[4,i] =  trans.dfdrb[2]*param[2]*binv*(t[i]-param[1])^2
    favg0[5,i] =  trans.dfdrb[2]*param[3]*binv
    favg0[6:5+trans.n,i] = trans.dfdu
  end
end
return
end

compute_lightcurve!(trans,param,t,favg0,nt,neval0)
@time compute_lightcurve!(trans,param,t,favg0,nt,neval0)

#integrate_lightcurve!(trans,param,t,dt,favg1,nt,1e-4,3,neval1)
#@time integrate_lightcurve!(trans,param,t,dt,favg1,nt,1e-4,8,neval1)

ntol = 9
tol = logarithmspace(-10.0,-2.0,ntol)
K = 1
neval_mean = zeros(K,ntol)
precision = zeros(K,ntol)
timing = zeros(K,ntol)
for k=1:K
  texp = 0.01*10.0^(rand())
  q1,q2 = rand(2)
  u_n=[2sqrt(q1)*q2,q1*(1-2*q2)]
  r=0.01*10.0^rand()
  b=rand()*(1+r)
  trans = transit_init(r,b,u_n,true)
  tau = 1.0 # Duration of transit
  t .= linearspace(-0.5*(tau+texp),0.5*(tau+texp),nt)
  v = 2/tau*sqrt((1+r)^2-b^2)
#  t  .= x/v
  param = [0.0,v,b]   # [t_0,v,b_0]
#  t_nobin1 = time_ns()
#  compute_lightcurve!(trans,param,t,favg0,nt,neval0)
#  t_nobin = time_ns()-t_nobin1
  integrate_lightcurve!(trans,param,t,texp,favg2,nt,1e-12,30,neval2,depthmax2)
  println("r: ",r," b: ",b," texp: ",texp)
  for j=1:ntol
    t_bin1 = time_ns()
    integrate_lightcurve!(trans,param,t,texp,favg1,nt,tol[j],30,neval1,depthmax1)
    t_bin = time_ns()-t_bin1
    t_nobin1 = time_ns()
    compute_lightcurve!(trans,param,t,favg0,nt,neval1)
    t_nobin = time_ns()-t_nobin1
    # Now time an "unbinned" version:
    precision[k,j] += maximum(abs.(favg1-favg2))
    neval_mean[k,j] += mean(neval1)
    println("k: ",k," j: ",j," tol: ",tol[j]," neval: ",mean(neval1)," relative time: ",t_bin/t_nobin," max depth: ",maximum(depthmax1))
  end
  clf()
  plot(t,favg2[1,:])
  plot(reverse(t),favg2[1,:])
  if VERSION >= v"0.7"
    read(stdin,Char)
  else
    read(STDIN,Char)
  end
end
precision_tol = zeros(ntol)
neval_tol = zeros(ntol)
for j=1:ntol
  precision_tol[j]=median(precision[:,j])
  neval_tol[j] = median(neval_mean[:,j])
end
clf()
loglog(precision_tol,neval_tol)


text:
immutable BrownianMotion <: ContinuousUnivariateStochasticProcess
  t::Vector{Float64}
  n::Int64

  function BrownianMotion(t::Vector{Float64}, n::Int64)
    t[1] == 0. || error("Starting time point must be equal to 0.0.")
    issorted(t, lt=<=) || error("The time points must be strictly sorted.")
    length(t) == n || error("Number of time points must be equal to the vector holding the time points.")
    new(t, n)
  end
end

BrownianMotion(t::Vector{Float64}) = BrownianMotion(t, Int64(length(t)))
BrownianMotion(t::Range) = BrownianMotion(collect(t), Int64(length(t)))
BrownianMotion(t::Float64, n::Int64) = BrownianMotion(0.0:t/n:t-t/n)
BrownianMotion(t::Float64) = BrownianMotion([t], 1)

BrownianMotion(t::Matrix{Float64}) = BrownianMotion[BrownianMotion(t[:, i]) for i = 1:size(t, 2)]
BrownianMotion(t::Range, np::Int) = BrownianMotion[BrownianMotion(t) for i = 1:np]
BrownianMotion(t::Float64, n::Int64, np::Int) = BrownianMotion[BrownianMotion(t, n) for i = 1:np]

function rand!(p::BrownianMotion, x::Vector{Float64})
  x[1] = 0.
  for i = 2:p.n
    x[i] = rand(Normal(0.0, sqrt(p.t[i]-p.t[i-1])))+x[i-1]
  end
  x
end

function rand!(p::Vector{BrownianMotion}, x::Matrix{Float64})
  for j = 1:length(p)
    x[1, j] = 0.
    for i = 2:p[j].n
      x[i, j] = rand(Normal(0.0, sqrt(p[j].t[i]-p[j].t[i-1])))+x[i-1, j]
    end
  end
  x
end

rand(p::BrownianMotion) = rand!(p, Array(Float64, p.n))

rand(p::Vector{BrownianMotion}) = rand!(p, Array(Float64, p[1].n, length(p)))


text:
module TestString
using Test
using Unicode
using CategoricalArrays

@testset "AbstractString operations on values of CategoricalPool{String}" begin
    pool = CategoricalPool(["", "café"])

    v1 = CategoricalValue(1, pool)
    v2 = CategoricalValue(2, pool)

    @test typeof(promote("a", v1)) == Tuple{String,String}

    @test String(v1)::String == ""
    @test String(v2)::String == "café"

    @test Symbol(v1) === Symbol("")
    @test Symbol(v2) === :café

    @test v1 == ""
    @test v2 == "café"
    @test v2 != ""

    @test isequal(v1, "")
    @test isequal(v2, "café")
    @test !isequal(v2, "")

    @test isempty(v1)
    @test !isempty(v2)

    @test eltype(v1) === Char
    @test eltype(v2) === Char

    @test length(v1) === 0
    @test length(v2) === 4

    @test sizeof(v1) === 0
    @test sizeof(v2) === 5

    @test_throws BoundsError nextind(v1, 1)
    @test nextind(v2, 4) === 6

    @test prevind(v1, 1) === 0
    @test prevind(v2, 6) === 4

    @test firstindex(v1) === 1
    @test firstindex(v2) === 1

    @test lastindex(v1) === 0
    @test lastindex(v2) === 4

    @test collect(v1) == []
    @test collect(v2) == collect("café")

    @test v2[2] === 'a'
    @test v2[4] === 'é'
    @test_throws BoundsError v1[1]
    @test_throws StringIndexError v2[5]

    @test codeunit(v2, 2) === 0x61
    @test codeunit(v2, 5) === 0xa9
    @test_throws BoundsError codeunit(v1, 1)
    @test_throws BoundsError codeunit(v2, 6)

    @test ascii(v1)::String == ""
    @test_throws ArgumentError ascii(v2)

    @test Unicode.normalize(v1) == ""
    @test Unicode.normalize(v2) == "café"
    @test Unicode.normalize(v2, :NFKD) == "café"

    @test isempty(collect(graphemes(v1)))
    @test collect(graphemes(v2)) == collect(graphemes("café"))

    @test isvalid(v1)
    @test isvalid(v2)
    @test !isvalid(v1, 1)
    @test isvalid(v2, 4)
    @test !isvalid(v2, 5)

    @test_throws BoundsError length(v1, 0, 0)
    @test length(v2, 1, 4) === 4

    @test_throws BoundsError nextind(v1, 1, 1)
    @test nextind(v2, 1, 2) === 3

    @test string(v1) == ""
    @test string(v2) == "café"
    @test string(v1, "a") == "a"
    @test string(v1, "a", v2) == "acafé"
    @test string(v2, 1) == "café1"

    @test sprint(print, v1) == ""
    @test sprint(print, v2) == "café"

    @test repr(v1) == "\"\""
    @test repr(v2) == "\"café\""

    @test "a" * v1 == "a"
    @test "a" * v1 * "b" == "ab"
    @test "a" * v1 * v2 == "acafé"

    @test v1^1 == ""
    @test v2^1 == "café"
    @test v1^2 == ""
    @test v2^2 == "cafécafé"

    @test repeat(v1, 10) == ""
    @test repeat(v2, 2) == "cafécafé"

    @test isempty(collect(eachmatch(r"af", v1)))
    @test first(eachmatch(r"af", v2)).offset == 2

    @test match(r"af", v1) === nothing
    @test match(r"af", v2).offset === 2
    @test match(r"af", v2, 2).offset === 2
    @test match(r"af", v2, 2, UInt32(0)).offset === 2

    @test collect((m.match for m ∈ eachmatch(r"af", v1))) == []
    @test collect((m.match for m ∈ eachmatch(r"af", v2))) == ["af"]
    @test collect((m.match for m ∈ eachmatch(r"af", v2, overlap=true))) == ["af"]

    @test lpad(v1, 1) == " "
    @test lpad(v2, 1) == "café"
    @test lpad(v2, 5) == " café"

    @test rpad(v1, 1) == " "
    @test rpad(v2, 1) == "café"
    @test rpad(v2, 5) == "café "

    @test findfirst("", v1) === 1:0
    @test findfirst("a", v2) === 2:2
    @test findfirst(==('a'), v2) === 2
    @test findnext(==('a'), v2, 3) === nothing

    @test findlast("a", v1) === nothing
    @test findlast("a", v2) === 2:2
    @test findlast(==('a'), v2) === 2
    @test findprev(==('a'), v2, 1) === nothing

    @test !occursin("a", v1)
    @test occursin("", v1)
    @test occursin("fé", v2)

    @test !occursin(r"af", v1)
    @test occursin(r"af", v2)

    @test startswith(v1, "")
    @test !startswith(v1, "a")
    @test startswith(v2, "caf")

    @test endswith(v1, "")
    @test !endswith(v1, "a")
    @test endswith(v2, "fé")

    @test reverse(v1) == ""
    @test reverse(v2) == "éfac"

    @test replace(v1, "a"=>"b") == ""
    @test replace(v2, 'a'=>'b') == "cbfé"
    @test replace(v2, "ca"=>"b", count=1) == "bfé"

    @test isempty(split(v1))
    @test split(v1, "a") == [""]
    @test split(v2) == ["café"]
    @test split(v2, "f") == ["ca", "é"]

    @test rsplit(v1, "a") == [""]
    @test rsplit(v2, "f") == ["ca", "é"]

    @test strip(v1) == ""
    @test strip(v2) == v2
    @test strip(v2, 'é') == "caf"

    @test lstrip(v1) == ""
    @test lstrip(v2) == v2
    @test lstrip(v2, 'é') == "café"

    @test rstrip(v1) == ""
    @test rstrip(v2) == v2
    @test rstrip(v2, 'é') == "caf"

    @test uppercase(v1) == ""
    @test uppercase(v2) == "CAFÉ"

    @test lowercase(v1) == ""
    @test lowercase(v2) == "café"

    @test titlecase(v1) == ""
    @test titlecase(v2) == "Café"

    @test uppercasefirst(v1) == ""
    @test uppercasefirst(v2) == "Café"

    @test lowercasefirst(v1) == ""
    @test lowercasefirst(v2) == "café"

    @test join([v1, "a"]) == "a"
    @test join([v1, "a"], v2) == "caféa"

    if VERSION >= v"1.1.1" # Previous versions throw a BoundsError
        @test chop(v1) == ""
    end
    @test chop(v2) == "caf"

    @test chomp(v1) == ""
    @test chomp(v2) == "café"

    @test textwidth(v1) === 0
    @test textwidth(v2) === 4

    @test isascii(v1)
    @test !isascii(v2)

    @test escape_string(v1) == ""
    @test escape_string(v2) == "café"

    @test collect(v1) == Char[]
    @test collect(v2) == Char['c', 'a', 'f', 'é']
end

end


text:
using Test

@testset "Tonti Diagrams" begin
  include("TontiDiagrams.jl")
end


text:
"""
    convolve(inp, kernel)

Convolves an inpL x inpL array with the FFT of a centered kernel
of size kerL x kerL to produce an output of size (kerL - inpL) x (kerL - inpL).

TODO: maybe this should return a real array
"""
function convolve(inp, kernel)
    inpL = size(inp)[1]
    kerL = size(kernel)[1]
    outL = kerL - inpL

    arr_pad = [inp zeros(inpL, outL); zeros(outL, inpL) zeros(outL, outL)]
    #out_pad = P \ ((P * arr_pad) .* kernel)
    out_pad = planned_ifft(planned_fft(arr_pad) .* kernel)
    out = out_pad[inpL+1:kerL, inpL+1:kerL]

    out
end

function convolveT(out, kernel)
    outL = size(out)[1]
    kerL = size(kernel)[1]
    inpL = kerL - outL

    out_pad = [zeros(inpL, inpL) zeros(inpL, outL); zeros(outL, inpL) out]
    #arr_pad = P * ((P \ out_pad) .* kernel)
    arr_pad = planned_fft(planned_ifft(out_pad) .* kernel)
    arr = arr_pad[1:inpL, 1:inpL]

    arr
end

# padded stores the intermediate padded result
function convolve!(inp, kernel, padded)
    inpL = size(inp)[1]
    kerL = size(kernel)[1]
    outL = kerL - inpL

    padded .= 0
    padded[1:inpL, 1:inpL] .= inp
    planned_fft!(padded)
    padded .*= kernel
    planned_ifft!(padded)
    out = @view padded[inpL+1:kerL, inpL+1:kerL]

    out
end

function convolveT!(out, kernel, padded)
    outL = size(out)[1]
    kerL = size(kernel)[1]
    inpL = kerL - outL

    padded .= 0
    padded[inpL+1:inpL+outL, inpL+1:inpL+outL] .= out
    planned_ifft!(padded)
    padded .*= kernel
    planned_fft!(padded)
    inp = @view padded[1:inpL, 1:inpL]

    inp
end


text:
"""
    cascademodel(S::Int64, Co::Float64)

Return matrix of the type `UnipartiteNetwork` randomly assembled according to
the cascade model for a given nuber of `S` and connectivity `Co`.

See also: `nichemodel`, `mpnmodel`, `nestedhierarchymodel`

#### References

Cohen, J.E., Newman, C.M., 1985. A stochastic theory of community food webs I.
Models and aggregated data. Proceedings of the Royal Society of London. Series
B. Biological Sciences 224, 421–448. https://doi.org/10.1098/rspb.1985.0042
"""
function cascademodel(S::Int64, Co::Float64)

    # Evaluate input
    maxconnectance = ((S^2-S)/2)/(S*S)
    Co >= maxconnectance && throw(ArgumentError("Connectance for $(S) species cannot be larger than $(maxconnectance) "))
    Co <= 0 && throw(ArgumentError("Connectance C must be positive"))
    S <= 0 && throw(ArgumentError("Number of species L must be positive"))

    # Initiate matrix
    A = UnipartiteNetwork(zeros(Bool, (S, S)))

    # For each species randomly asscribe rank e
    e = sort(rand(S), rev=false)

    # Probability for linking two species
    p = 2*Co*S/(S - 1)

    for consumer in  1:S

        # Rank for a consumer
        rank = e[consumer]

        # Get a set of resources smaller than consumer
        potentialresources = findall(x -> x > rank, e)  # indices
        #resourcesranks = η[potentialresources]         # ranks
        #R = length(potentialresources)                 # length

        # Each gets a link with a probability p
        for resource in potentialresources

            # Take the resources and for each of them check a random number
            # if it is smaller than probability of linking two speceis create
            # a link in the A matrix

            # Random number for a potential resource
            rand() < p && (A[consumer, resource] = true)

        end

    end

    return A

end

"""
    cascademodel(N::T) where {T <: UnipartiteNetwork}

Applied to a `UnipartiteNetwork` return its randomized version.

#### References

Cohen, J.E., Newman, C.M., 1985. A stochastic theory of community food webs I.
Models and aggregated data. Proceedings of the Royal Society of London. Series
B. Biological Sciences 224, 421–448. https://doi.org/10.1098/rspb.1985.0042
"""
function cascademodel(N::T) where {T <: UnipartiteNetwork}

    return cascademodel(richness(N), connectance(N))

end

"""
    cascademodel(S::Int64, L::Int64)

Number of links can be specified instead of connectance.

#### References

Cohen, J.E., Newman, C.M., 1985. A stochastic theory of community food webs I.
Models and aggregated data. Proceedings of the Royal Society of London. Series
B. Biological Sciences 224, 421–448. https://doi.org/10.1098/rspb.1985.0042
"""
function cascademodel(S::Int64, L::Int64)
    Co = (L/(S*S))
    return cascademodel(S, Co)

end

"""

    cascademodel(parameters::Tuple)

Parameters tuple can also be provided in the form (S::Int64, Co::Float64)
or (S::Int64, L::Int64).

#### References

Cohen, J.E., Newman, C.M., 1985. A stochastic theory of community food webs I.
Models and aggregated data. Proceedings of the Royal Society of London. Series
B. Biological Sciences 224, 421–448. https://doi.org/10.1098/rspb.1985.0042
"""
function cascademodel(parameters::Tuple)
    return cascademodel(parameters[1], parameters[2])
end


text:
module SymbolServer
module LoadingBay
end
export SymbolServerProcess, change_env, load_manifest_packages, load_project_packages, get_context, getstore

using Serialization, Pkg, SHA
using Base: UUID
@static if VERSION < v"1.1"
    const PackageEntry = Vector{Dict{String,Any}}
else
    using Pkg.Types: PackageEntry
end
include("symbols.jl")

mutable struct SymbolServerProcess
    process::Base.Process
    context::Union{Nothing,Pkg.Types.Context}
    depot::Dict{String,ModuleStore}
    process_stderr::Union{IOBuffer,Nothing}

    caching_packages::Set{UUID}
    newly_cached_packages::Vector{UUID}

    function SymbolServerProcess(;environment = nothing, depot = nothing)
        jl_cmd = joinpath(Sys.BINDIR, Base.julia_exename())
        server_script = joinpath(@__DIR__, "server.jl")

        env_to_use = copy(ENV)
        if depot !== nothing
            if depot == ""
                delete!(env_to_use, "JULIA_DEPOT_PATH")
            else
                env_to_use["JULIA_DEPOT_PATH"] = depot
            end
        end

        stderr_for_client_process = VERSION < v"1.1.0" ? nothing : IOBuffer()

        p = if environment === nothing
            open(pipeline(Cmd(`$jl_cmd $server_script`, env = env_to_use), stderr = stderr_for_client_process), read = true, write = true)
        else
            open(pipeline(Cmd(`$jl_cmd --project=$environment $server_script`, dir = environment, env = env_to_use), stderr = stderr_for_client_process), read = true, write = true)
        end
        ssp = new(p, nothing, deepcopy(stdlibs), stderr_for_client_process, Set{UUID}(), UUID[])
        get_context(ssp)
        return ssp
    end
end

function Base.show(io::IO, ssp::SymbolServerProcess)
    println(io, "SymbolServerProcess with $(length(ssp.depot)) ($(sum(!isempty(v.vals) for (k,v) in ssp.depot))) packages")

    print(join(sort!([string(isempty(v.vals) ? " ** " : "    ", k) for (k, v) in ssp.depot], lt = (a,b) ->a[5:end]<b[5:end]), "\n"))
end

function request(ssp::SymbolServerProcess, message::Symbol, payload)
    serialize(ssp.process, (message, payload))
    ret_val = try
        deserialize(ssp.process)
    catch err
        # Only Julia 1.1 and newer support capturing stderr into an IOBuffer
        if ssp.process_stderr !== nothing
            stderr_from_client_process = String(take!(ssp.process_stderr))

            complete_error_message = string(sprint(showerror, err), "\n\nstderr from client process:\n\n", stderr_from_client_process)

            error(complete_error_message)
        else
            complete_error_message = string(sprint(showerror, err), "\n\nCouldn't capture stderr from client process on julia 1.0.\n\n")

            error(complete_error_message)
        end
    end

    !(ret_val isa Tuple{Symbol,<:Any}) && error("Invalid response:\n", ret_val)
    return ret_val
end

"""
    load_manifest_packages(ssp)
Load all packages from the current active environments manifest into the client
side depot.
"""
function load_manifest_packages(ssp::SymbolServerProcess)
    # for uuid in keys(manifest(ssp.context))
    for pkg in manifest(ssp.context)
        load_package_cache(ssp, packageuuid(pkg))
    end
    update(ssp)
end

function load_project_packages(ssp::SymbolServerProcess)
    for uuid in values(deps(project(ssp.context)))
        load_package_cache(ssp, uuid)
    end
    update(ssp)
end

function getstore(ssp::SymbolServerProcess)
    load_project_packages(ssp)
    return ssp.depot
end

"""
    get_context(ssp)
Retrieves the active context (environment) from the server.
"""
function get_context(ssp::SymbolServerProcess)
    status, payload = request(ssp, :get_context, nothing)
    if status == :success
        ssp.context = payload
        return
    else
        error(payload)
    end
end
"""
    change_env(ssp, env_path)
Activates a new environment on the server. The new active context must then be retrieved separately.
"""
function change_env(ssp::SymbolServerProcess, env_path::String)
    status, payload = request(ssp, :change_env, env_path)
    if status == :success
        return payload
    else
        error(payload)
    end
end

"""
    load_package_cache(ssp::SymbolServerProcess, uuid::UUID)
Tries to load the on-disc stored cache for a package (uuid). Attempts to generate (and save to disc) a new cache if the file does not exist or is unopenable.
"""
function load_package_cache(ssp::SymbolServerProcess, uuid::UUID)
    storedir = abspath(joinpath(@__DIR__, "..", "store"))
    cache_path = joinpath(storedir, string(uuid, ".jstore"))

    if !isinmanifest(ssp.context, uuid)
        @info "Tried to load $uuid but failed to find it in the manifest."
        return
    end

    pe = frommanifest(ssp.context, uuid)
    pe_name = packagename(ssp.context, uuid)
    if isfile(cache_path)
        try
            store = open(cache_path) do io
                deserialize(io)
            end
            if version(pe) != store.ver || (store.ver isa String && endswith(store.ver, "+") && sha_pkg(pe) != store.sha)
                @info "$pe_name changed, updating cache."
                cache_package(ssp, uuid)
                store = open(cache_path) do io
                    deserialize(io)
                end
            end
            ssp.depot[pe_name] = store.val
            for dep in deps(pe)
                load_dependency_cache(ssp, packageuuid(dep))
            end
        catch err
            @info "Tried to load $pe_name but failed to load from disc, re-caching."
            if err isa UndefVarError && err.var in (:structStore, :abstractStore)
                @info "Package cache pre-v3.0"
            end
            rm(cache_path)
            cache_package(ssp, uuid)
        end
    else
        @info "$(pe_name) not stored on disc"
        cache_package(ssp, uuid)
    end
end

function load_dependency_cache(ssp::SymbolServerProcess, uuid::UUID)
    storedir = abspath(joinpath(@__DIR__, "..", "store"))
    cache_path = joinpath(storedir, string(uuid, ".jstore"))
    if !isinmanifest(ssp.context, uuid)
        @info "Tried to load $uuid cache as a dependency but failed to find it in the manifest."
        return
    end
    pe = frommanifest(ssp.context, uuid)
    pe_name = packagename(ssp.context, uuid)
    haskey(ssp.depot, pe_name) && return
    @info "loading dependency $pe_name"
    if isfile(cache_path)
        try
            store = open(cache_path) do io
                deserialize(io)
            end
            if version(pe) != store.ver || (store.ver isa String && endswith(store.ver, "+") && sha_pkg(pe) != store.sha)
                @info "Tried to load $pe_name cache as a dependency but failed."
            end
            ssp.depot[pe_name] = store.val
            for dep in deps(pe)
                load_dependency_cache(ssp, packageuuid(dep))
            end
        catch err
            if err isa UndefVarError && err.var in (:structStore, :abstractStore)
                @info "Package cache pre-v3.0"
            end
            # rm(cache_path)
            # cache_package(ssp, uuid)
        end
    else
        @info "Tried to load $pe_name cache as a dependency but no file found."
        # cache_package(ssp, uuid)
    end

end

load_package_cache(ssp::SymbolServerProcess, uuid::String) = load_package_cache(ssp, UUID(uuid))

function Base.kill(server::SymbolServerProcess)
    serialize(server.process, (:close, nothing))
end
"""
    cache_package(ssp, uuid::Union{UUID, Vector{UUID}})
Sends a request to the server to cache a package or collection of packages.
Requested packages are added to the list `ssp.caching_packages` to prevent the
sending multiple requests for the same package.

The server returns a list of packages that it has loaded and adds it to
`ssp.newly_cached_packages`. This list should be used to `update` the client side
depot `ssp.depot`.
"""
cache_package(ssp::SymbolServerProcess, uuid::UUID) = cache_package(ssp, [uuid])
function cache_package(ssp::SymbolServerProcess, uuid::Vector{UUID})
    if uuid in ssp.caching_packages
        return
    end
    union!(ssp.caching_packages, uuid)
    status, payload = request(ssp, :cache_package, string.(uuid))
    if status == :success
        for u in uuid
            delete!(ssp.caching_packages, u)
        end
        append!(ssp.newly_cached_packages, UUID.(first(p) for p in payload))
        return payload
    else
        error(payload)
    end
end

"""
    update(ssp)
Update the client side depot with newly cached packages.
"""
function update(ssp::SymbolServerProcess)
    for uuid in ssp.newly_cached_packages
        load_package_cache(ssp, uuid)
    end
    empty!(ssp.newly_cached_packages)
end

function clear_disc_store()
    storedir = abspath(joinpath(@__DIR__, "..", "store"))
    for f in readdir(storedir)
        if endswith(f, ".jstore")
            rm(joinpath(storedir, f))
        end
    end
end
include("utils.jl")
const stdlibs = load_core()
end # module


text:
module Gherkin
export tokenize

using JSON

"""Takes a .feature file and outputs a .feature.tokens file.

For example, given the `.feature` file:

```
Feature: Background

  Background: a simple background
    Given the minimalism inside a background

  Scenario: minimalistic
    Given the minimalism
```
output the `.feature.tokens` file:
```
(1:1)FeatureLine:Feature/Background/
(2:1)Empty://
(3:3)BackgroundLine:Background/a simple background/
(4:5)StepLine:Given /the minimalism inside a background/
(5:1)Empty://
(6:1)Empty://
(7:3)ScenarioLine:Scenario/minimalistic/
(8:5)StepLine:Given /the minimalism/
EOF
```

and with a language:
```
#language:no
Egenskap: i18n support

  Scenario: Parsing many languages
    Gitt Gherkin supports many languages
    Når Norwegian keywords are parsed
    Så they should be recognized
```
tokenizes to
```
(1:1)Language:/no/
(2:1)FeatureLine:Egenskap/i18n support/
(3:1)Empty://
(4:3)ScenarioLine:Scenario/Parsing many languages/
(5:5)StepLine:Gitt /Gherkin supports many languages/
(6:5)StepLine:Når /Norwegian keywords are parsed/
(7:5)StepLine:Så /they should be recognized/
EOF
```

"""

const DEFAULT_LANG = "en"
const I18N = JSON.parsefile(joinpath(dirname(@__FILE__), "glangs.json"))
const LANGUAGE_RE = r"^\s*#\s*language\s*:(.*)$"
const TITLE_KEYWORDS = ["feature", "scenario", "scenarioOutline",
                        "background", "examples"]
const STEPS = ["but", "and", "when", "given", "then"]

type Parser
    lines::Vector{UTF8String}
    lang::UTF8String
    keywords::Dict{AbstractString, Any} # TODO remove this
    tokens::Vector{UTF8String}
end
Parser(lines) = Parser(lines, DEFAULT_LANG, I18N[DEFAULT_LANG], UTF8String[])

function tokenize(filename::AbstractString)
    # TODO stream this rather than read in at once
    tokenize(readlines(open(filename)))
end

function tokenize(lines::Vector)
    p = Parser(lines)
    nextline = get_lang!(p)
    s = no_docstring()
    for i in nextline:length(p.lines)
        s = parseline!(i, p, s)
    end
    push!(p.tokens, "EOF")
    p.tokens
end

"""The language may be specified in the first line of the `.features` file.
If the language is not present use `DEFAULT_LANG` (`"en"`).

Return line to continue parsing, i.e. 2 if there is a language line otherwise 1.
"""
function get_lang!(p::Parser)
    length(p.lines) === 0 && return 1
    m = match(LANGUAGE_RE, p.lines[1])
    m === nothing && return 1

    p.lang = strip(m.captures[1])
    # TODO pushing to tokens here feels like a complete hack
    # This should really be done as part of the next section (more generically)
    # BUT maybe it needs to be here, as can only be at the start...
    @assert p.lang in keys(I18N) "The language \"$lang\" is not supported."
    push!(p.tokens, "(1:1)Language:/$(p.lang)/")
    p.keywords = I18N[p.lang]
    return 2
end

type DocstringState
    in_docstring::Bool
    indent::Int
    sep::Char
end
no_docstring() = DocstringState(false, 0, ' ')
start_docstring(i::Int, c::Char) = DocstringState(true, i, c)
in_docstring(s::DocstringState) = s.in_docstring

function parseline!(i::Int, p::Parser, s::DocstringState)
    in_docstring(s) ? parsedocstring!(i, p, s) : parseline!(i, p)
end

function parseline!(i::Int, p::Parser)
    line = p.lines[i]
    j = _get_indent(line)
    line = strip(line)
    # match against things from p.lang....
    # is there a hierarchy for this? feature -> scenario -> given ...?
    # TODO specifically how to distinguish between empty Other and Empty.
    if line == ""
        push!(p.tokens, "($i:1)Empty://")
        return no_docstring()
    end
    for t in TITLE_KEYWORDS
        b, k = startswith_titlekeyword(line, t, p)
        if b
            rest = split(line, k * ":"; limit=2)[2] # TODO remove this hack
            if startswith(rest, ' ')
                rest = rest[2:end]
            end
            push!(p.tokens, "($i:$j)$(_titlecase(t))Line:$k/$rest/")
            return no_docstring()
        end
    end
    for s in STEPS
        b, k= startswith_step(line, s, p)
        if b
            rest = split(line, k; limit=2)[2] # TODO remove this hack
            push!(p.tokens, "($i:$j)StepLine:$k/$rest/")
            return no_docstring()
        end
    end
    if startswith(line, "#")
        push!(p.tokens, "($i:1)Comment:/" * " " ^ (j - 1) * line * "/")
        return no_docstring()
    end
    if startswith(line, "|")
        # TODO test escaping better
        entries = split(line, r"(?!<\\)\|")
        columns = (j-1) .+ cumsum(map(length, entries) .+ 1)[1:end-2]
        entries = entries[2:end-1]
        row = join(["$(c + _get_indent(v)):$(strip(v))"
                    for (c, v) in zip(columns, entries)], ",")

        push!(p.tokens, "($i:$j)TableRow://$row")
        return no_docstring()
    end
    if startswith(line, "@")
        inds = findin(line, "@")
        tags = join(["$(j-1+x[1]):$(line[x[1]:x[2]-1])" for x in zip(inds, inds[2:end])], ",")
        tags = tags * (tags == "" ? "" : ",") * "$(j-1+inds[end]):$(line[inds[end]:end])"
        push!(p.tokens, "($i:$j)TagLine://$tags")
        return no_docstring()
    end
    if startswith(line, "\"\"\"") || startswith(line, "```")
        content_type = strip(line[4:end])
        push!(p.tokens, "($i:$j)DocStringSeparator:/$content_type/")
        return start_docstring(j, line[1])
    end

    push!(p.tokens, "($i:1)Other:/" * " "^ (j - 1) * line * "/")
    no_docstring()
end

function parsedocstring!(i::Int, p::Parser, s::DocstringState)
    @assert in_docstring(s)
    line = rstrip(p.lines[i], '\n')
    if startswith(strip(line), "$(s.sep)" ^ 3)
        push!(p.tokens, "($i:$(s.indent))DocStringSeparator://")
        return no_docstring()
    end
    if !startswith(line, " " ^ (s.indent - 1))
        # This is a bad docstring indent
        push!(p.tokens, "($i:1)Other:/$(strip(line))/")
    else
        push!(p.tokens, "($i:1)Other:/$(line[s.indent:end])/")
    end
    s
end

function startswith_titlekeyword(line, title_keyword, p)
    for k in p.keywords[title_keyword]
        startswith(line, k * ":") && return true, k
    end
    return false, ""
end

function startswith_step(line, keyword, p)
    for k in p.keywords[keyword]
        startswith(line, k) && return true, k
    end
    return false, ""
end

function _titlecase(word)
    # Note this isn't portable
    string(uppercase(word[1])) * word[2:end]
end

function _get_indent(line)
    i = 1
    for c in line
        if c != ' '
            return i
        end
        i += 1
    end
    # potentially we could return 0
    return i
end

end # module


text:
using Flux
using WindMixing
using JLD2
using FileIO
using OceanParameterizations
using OrdinaryDiffEq
using Random
using GalacticOptim
using LinearAlgebra

train_files = [
        "-1e-3",
]

𝒟train = WindMixing.data(train_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)
#
PATH = pwd()

OUTPUT_PATH = joinpath(PATH, "training_output")
# OUTPUT_PATH = "D:\\University Matters\\MIT\\CLiMA Project\\OceanParameterizations.jl\\training_output"

EXTRACTED_OUTPUT_PATH = joinpath(PATH, "extracted_training_output")

FILE_NAME = "test"
FILE_PATH = joinpath(OUTPUT_PATH, "$(FILE_NAME).jld2")
EXTRACTED_FILE_PATH = joinpath(EXTRACTED_OUTPUT_PATH, "$(FILE_NAME)_extracted.jld2")
@assert !isfile(FILE_PATH)

N_inputs = 96
hidden_units = 400
N_outputs = 31

weights, re = Flux.destructure(Chain(Dense(N_inputs, hidden_units, relu), Dense(hidden_units, N_outputs)))

uw_NN = re(weights ./ 1f5)
vw_NN = re(weights ./ 1f5)
wT_NN = re(weights ./ 1f5)

gradient_scaling = 5f-3
train_parameters = Dict("ν₀" => 1f-4, "ν₋" => 0.1f0, "Riᶜ" => 0.25f0, "ΔRi" => 1f-1, "Pr" => 1f0, "κ" => 10f0,
                        "modified_pacanowski_philander" => true, "convective_adjustment" => false,
                        "smooth_profile" => false, "smooth_NN" => false, "smooth_Ri" => false, "train_gradient" => true,
                        "zero_weights" => true, "unscaled" => false, "gradient_scaling" => gradient_scaling)

train_epochs = [1]
train_tranges = [1:20:200]
train_iterations = [3]
train_optimizers = [[ADAM(2e-4)]]

timestepper = Rosenbrock23()

function train(FILE_PATH, train_files, train_epochs, train_tranges, train_parameters, train_optimizers, train_iterations, uw_NN, vw_NN, wT_NN, 𝒟train, timestepper, unscaled)
    write_metadata_NDE_training(FILE_PATH, train_files, train_epochs, train_tranges, train_parameters, train_optimizers, uw_NN, vw_NN, wT_NN)
    if unscaled
        for i in 1:length(train_epochs)
            @info "iteration $i/$(length(train_epochs)), time range $(train_tranges[i])"
            # uw_NN, vw_NN, wT_NN = train_NDE_convective_adjustment(uw_NN, vw_NN, wT_NN, 𝒟train, train_tranges[i], timestepper, train_optimizers[i], train_epochs[i], FILE_PATH, 1, 1, 10f0, 5)
            if train_parameters["modified_pacanowski_philander"]
                uw_NN, vw_NN, wT_NN = train_NDE_unscaled(uw_NN, vw_NN, wT_NN, 𝒟train, train_tranges[i], timestepper, train_optimizers[i], train_epochs[i], FILE_PATH, i, n_simulations=length(train_files), maxiters=train_iterations[i],
                modified_pacanowski_philander=train_parameters["modified_pacanowski_philander"], convective_adjustment=train_parameters["convective_adjustment"],
                ν₀=train_parameters["ν₀"], ν₋=train_parameters["ν₋"], ΔRi=train_parameters["ΔRi"], Riᶜ=train_parameters["Riᶜ"],
                κ=train_parameters["κ"],
                smooth_profile=train_parameters["smooth_profile"], smooth_NN=train_parameters["smooth_NN"], smooth_Ri=train_parameters["smooth_Ri"], train_gradient=train_parameters["train_gradient"],
                gradient_scaling=train_parameters["gradient_scaling"])
            else
                uw_NN, vw_NN, wT_NN = train_NDE_unscaled(uw_NN, vw_NN, wT_NN, 𝒟train, train_tranges[i], timestepper, train_optimizers[i], train_epochs[i], FILE_PATH, i, n_simulations=length(train_files), maxiters=train_iterations[i],
                modified_pacanowski_philander=train_parameters["modified_pacanowski_philander"], convective_adjustment=train_parameters["convective_adjustment"],
                κ=train_parameters["κ"],
                smooth_profile=train_parameters["smooth_profile"], smooth_NN=train_parameters["smooth_NN"], smooth_Ri=train_parameters["smooth_Ri"], train_gradient=train_parameters["train_gradient"],
                gradient_scaling=train_parameters["gradient_scaling"])
            end
        end
    else
        for i in 1:length(train_epochs)
            @info "iteration $i/$(length(train_epochs)), time range $(train_tranges[i])"
            # uw_NN, vw_NN, wT_NN = train_NDE_convective_adjustment(uw_NN, vw_NN, wT_NN, 𝒟train, train_tranges[i], timestepper, train_optimizers[i], train_epochs[i], FILE_PATH, 1, 1, 10f0, 5)
            if train_parameters["modified_pacanowski_philander"]
                uw_NN, vw_NN, wT_NN = train_NDE(uw_NN, vw_NN, wT_NN, 𝒟train, train_tranges[i], timestepper, train_optimizers[i], train_epochs[i], FILE_PATH, i, n_simulations=length(train_files), maxiters=train_iterations[i],
                modified_pacanowski_philander=train_parameters["modified_pacanowski_philander"], convective_adjustment=train_parameters["convective_adjustment"],
                ν₀=train_parameters["ν₀"], ν₋=train_parameters["ν₋"], ΔRi=train_parameters["ΔRi"], Riᶜ=train_parameters["Riᶜ"],
                κ=train_parameters["κ"],
                smooth_profile=train_parameters["smooth_profile"], smooth_NN=train_parameters["smooth_NN"], smooth_Ri=train_parameters["smooth_Ri"], train_gradient=train_parameters["train_gradient"],
                zero_weights = train_parameters["zero_weights"],
                gradient_scaling=train_parameters["gradient_scaling"])
            else
                uw_NN, vw_NN, wT_NN = train_NDE(uw_NN, vw_NN, wT_NN, 𝒟train, train_tranges[i], timestepper, train_optimizers[i], train_epochs[i], FILE_PATH, i, n_simulations=length(train_files), maxiters=train_iterations[i],
                modified_pacanowski_philander=train_parameters["modified_pacanowski_philander"], convective_adjustment=train_parameters["convective_adjustment"],
                κ=train_parameters["κ"],
                smooth_profile=train_parameters["smooth_profile"], smooth_NN=train_parameters["smooth_NN"], smooth_Ri=train_parameters["smooth_Ri"], train_gradient=train_parameters["train_gradient"],
                zero_weights = train_parameters["zero_weights"],
                gradient_scaling=train_parameters["gradient_scaling"])
            end
        end
    end
    return uw_NN, vw_NN, wT_NN
end

uw_NN_res, vw_NN_res, wT_NN_res = train(FILE_PATH, train_files, train_epochs, train_tranges, train_parameters, train_optimizers, train_iterations, uw_NN, vw_NN, wT_NN, 𝒟train, timestepper, train_parameters["unscaled"])


text:
"Parameters for right-hand side function of variational partitioned Runge-Kutta methods."
mutable struct ParametersVPRKpStandard{DT, TT, ET <: IODE{DT,TT}, D, S} <: AbstractParametersVPRK{DT,TT,ET,D,S}
    equ::ET
    tab::TableauVPRK{TT}
    Δt::TT

    RU::Vector{TT}
    RG::Vector{TT}

    RU1::Vector{TT}
    RU2::Vector{TT}
    RG1::Vector{TT}
    RG2::Vector{TT}

    t::TT
    q::Vector{DT}
    p::Vector{DT}

    function ParametersVPRKpStandard{DT,TT,ET,D,S}(equ, tab, Δt, RU, RG, R∞) where {DT,TT,ET,D,S}
        RU1 = [RU[1], zero(TT)]
        RU2 = [zero(TT), R∞ * RU[2]]
        RG1 = [RG[1], zero(TT)]
        RG2 = [zero(TT), R∞ * RG[2]]
        new(equ, tab, Δt, RU, RG, RU1, RU2, RG1, RG2, zero(TT), zeros(DT,D), zeros(DT,D))
    end
end

function ParametersVPRKpStandard(equ::ET, tab::TableauVPRK{TT}, Δt::TT, R::Vector, R∞=tab.R∞) where {DT, TT, ET <: IODE{DT,TT}}
    R = Vector{TT}(R)
    ParametersVPRKpStandard{DT, TT, ET, equ.d, tab.s}(equ, tab, Δt, R, R, R∞)
end

function ParametersVPRKpStandard(equ::ET, tab::TableauVPRK{TT}, Δt::TT, RU::Vector, RG::Vector, R∞=tab.R∞) where {DT, TT, ET <: IODE{DT,TT}}
    RU = Vector{TT}(RU)
    RG = Vector{TT}(RG)
    ParametersVPRKpStandard{DT, TT, ET, equ.d, tab.s}(equ, tab, Δt, RU, RG, R∞)
end

function update_params!(params::ParametersVPRKpStandard, cache::IntegratorCacheVPRK)
    # set time for nonlinear solver and copy previous solution
    params.t  = cache.t
    params.q .= cache.q
    params.p .= cache.p
end

function compute_projection!(
                x::Vector{ST}, q::SolutionVector{ST}, p::SolutionVector{ST}, λ::SolutionVector{ST},
                U::Vector{Vector{ST}}, G::Vector{Vector{ST}},
                params::ParametersVPRKpStandard{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    @assert D == length(q) == length(p) == length(λ)
    @assert D == length(U[1]) == length(U[2])
    @assert D == length(G[1]) == length(G[2])

    # copy x to q, λ
    for k in 1:D
        q[k] = x[0*D+k]
        λ[k] = x[1*D+k]
    end

    # compute u=λ and g=∇ϑ(q)⋅λ
    U[1] .= λ
    U[2] .= λ

    params.equ.g(params.t, q, λ, G[1])
    G[2] .= G[1]

    # compute p̅=ϑ(q)
    params.equ.ϑ(params.t, q, λ, p)
end

"Compute stages of projected variational partitioned Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST},
                params::ParametersVPRKpStandard{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    cache = NonlinearFunctionCacheVPRKprojection{ST}(D,S)

    quote
        @assert length(x) == length(b)

        compute_projection!(x, $cache.q, $cache.p, $cache.λ, $cache.U, $cache.G, params)

        # compute b = - [q̅-q-U]
        for k in 1:D
            b[0*D+k] = - ($cache.q[k] - params.q[k]) + params.Δt * params.RU[2] * $cache.U[2][k]
        end

        # compute b = - [p̅-p-G]
        for k in 1:D
            b[1*D+k] = - ($cache.p[k] - params.p[k]) + params.Δt * params.RG[2] * $cache.G[2][k]
        end
    end
end

"Variational partitioned Runge-Kutta integrator."
mutable struct IntegratorVPRKpStandard{DT, TT,
                SPT <: ParametersVPRK{DT,TT},
                PPT <: ParametersVPRKpStandard{DT,TT},
                SST <: NonlinearSolver{DT},
                STP <: NonlinearSolver{DT},
                IT  <: InitialGuessPODE{DT,TT}} <: AbstractIntegratorVPRKwProjection{DT,TT}

    sparams::SPT
    pparams::PPT

    solver::SST
    projector::STP
    iguess::IT
end

function IntegratorVPRKpStandard(args...)
    IntegratorVPRKpStandardConstructor(args..., [0,1]; R∞=1)
end

function IntegratorVPRKpSymplectic(args...)
    IntegratorVPRKpStandardConstructor(args..., [1,1])
end

function IntegratorVPRKpVariationalQ(args...)
    IntegratorVPRKpStandardConstructor(args..., [1,0], [0,1]; R∞=1)
end

function IntegratorVPRKpVariationalP(args...)
    IntegratorVPRKpStandardConstructor(args..., [0,1], [1,0]; R∞=1)
end

function IntegratorVPRKpStandardConstructor(equation, tableau, Δt, R; R∞=tableau.R∞)
    IntegratorVPRKpStandardConstructor(equation, tableau, Δt, R, R; R∞=R∞)
end

function IntegratorVPRKpStandardConstructor(equation::ET, tableau::TableauVPRK{TT}, Δt::TT,
                                    RU::Vector, RG::Vector; R∞::Int=tableau.R∞) where {DT, TT, ET <: IODE{DT,TT}}
    D = equation.d
    M = equation.n
    S = tableau.s

    # create solver params
    sparams = ParametersVPRK(equation, tableau, Δt)

    # create projector params
    pparams = ParametersVPRKpStandard(equation, tableau, Δt, RU, RG, R∞)

    # create nonlinear solver
    solver = create_nonlinear_solver(DT, D*S, sparams)

    # create projector
    projector = create_nonlinear_solver(DT, 2*D, pparams)

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

    # create integrator
    IntegratorVPRKpStandard{DT, TT, typeof(sparams), typeof(pparams), typeof(solver), typeof(projector), typeof(iguess)}(
                sparams, pparams, solver, projector, iguess)
end

equation(integrator::IntegratorVPRKpStandard) = integrator.sparams.equ
timestep(integrator::IntegratorVPRKpStandard) = integrator.sparams.Δt
tableau(integrator::IntegratorVPRKpStandard) = integrator.sparams.tab
nstages(integrator::IntegratorVPRKpStandard) = integrator.sparams.tab.s

function initialize!(int::IntegratorVPRKpStandard{DT,TT}, cache::IntegratorCacheVPRK) where {DT,TT}
    cache.t̅ = cache.t - timestep(int)

    equation(int).v(cache.t, cache.q, cache.p, cache.v)
    equation(int).f(cache.t, cache.q, cache.p, cache.f)

    initialize!(int.iguess, cache.t, cache.q, cache.p, cache.v, cache.f,
                            cache.t̅, cache.q̅, cache.p̅, cache.v̅, cache.f̅)

    # initialise projector
    cache.U[1] .= cache.λ
    equation(int).g(cache.t, cache.q, cache.λ, cache.G[1])
end

function initial_guess!(int::IntegratorVPRKpStandard, cache::IntegratorCacheVPRK)
    for i in 1:nstages(int)
        evaluate!(int.iguess, cache.q, cache.p, cache.v, cache.f,
                              cache.q̅, cache.p̅, cache.v̅, cache.f̅,
                              cache.q̃, cache.ṽ,
                              tableau(int).q.c[i])

        for k in 1:ndims(int)
            int.solver.x[ndims(int)*(i-1)+k] = cache.ṽ[k]
        end
    end
end

function initial_guess_projection!(int::IntegratorVPRKpStandard, cache::IntegratorCacheVPRK)
    offset_q = 0
    offset_λ = ndims(int)

    for k in 1:ndims(int)
        int.projector.x[offset_q+k] = cache.q[k]
        int.projector.x[offset_λ+k] = 0
    end
end

"Integrate ODE with variational partitioned Runge-Kutta integrator."
function integrate_step!(int::IntegratorVPRKpStandard{DT,TT}, cache::IntegratorCacheVPRK{DT,TT}) where {DT,TT}
    # add perturbation for next time step to solution
    # (same vector field as previous time step)
    project_solution!(int, cache, int.pparams.RU1, int.pparams.RG1)

    # update nonlinear solver parameters from cache
    update_params!(int.sparams, cache)

    # compute initial guess
    initial_guess!(int, cache)

    # reset cache
    reset!(cache, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params, cache.n)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params, cache.n)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, cache.Q, cache.V, cache.P, cache.F, int.sparams)

    # compute unprojected solution
    update_solution!(int, cache)

    # set time and solution for projection solver
    update_params!(int.pparams, cache)

    # set initial guess for projection
    initial_guess_projection!(int, cache)

    # call projection solver
    solve!(int.projector)

    # print solver status
    print_solver_status(int.projector.status, int.projector.params, cache.n)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.projector.status, int.projector.params, cache.n)

    # compute projection vector fields
    compute_projection!(int.projector.x, cache.q̃, cache.p̃, cache.λ, cache.U, cache.G, int.pparams)

    # add projection to solution
    project_solution!(int, cache, int.pparams.RU2, int.pparams.RG2)

    # copy solution to initial guess
    update!(int.iguess, cache.t, cache.q, cache.p, cache.v, cache.f)

    # take care of periodic solutions
    cut_periodic_solution!(cache, equation(int).periodicity)
end


text:
using Test
using JLD2
using Statistics
using LinearAlgebra

using Oceananigans
using Oceananigans.Units
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize
using Oceananigans: fields

using ParameterEstimocean
using ParameterEstimocean.EnsembleKalmanInversions: iterate!
using ParameterEstimocean.EnsembleKalmanInversions: FullEnsembleDistribution, Resampler
using ParameterEstimocean.EnsembleKalmanInversions: resample!
using ParameterEstimocean.InverseProblems: inverting_forward_map

data_path = "convective_adjustment_test.jld2"
Nensemble = 3
architecture = CPU()

@testset "EnsembleKalmanInversions tests" begin
    @info "  Testing EnsembleKalmanInversion..."

    #####
    ##### Build two InverseProblem
    #####

    observation = SyntheticObservations(data_path, field_names=(:u, :v, :b))

    file = jldopen(observation.path)
    closure = file["serialized/closure"]
    close(file)

    simulation = ensemble_column_model_simulation(observation;
                                                  Nensemble,
                                                  architecture,
                                                  closure,
                                                  tracers = :b)

    batched_observations = BatchedSyntheticObservations([observation, observation], weights=[0.5, 1.0])
    batched_simulation = ensemble_column_model_simulation(batched_observations;
                                                          Nensemble,
                                                          architecture,
                                                          closure,
                                                          tracers = :b)

    for sim in [simulation, batched_simulation]
        sim.model.velocities.u.boundary_conditions.top.condition .= observation.metadata.parameters.Qᵘ
        sim.model.tracers.b.boundary_conditions.top.condition .= observation.metadata.parameters.Qᵇ
        sim.model.tracers.b.boundary_conditions.bottom.condition .= observation.metadata.parameters.N²
        sim.Δt = observation.metadata.parameters.Δt
    end

    lower_bound, upper_bound = bounds = [0.9, 1.1]
    priors = (convective_κz = ScaledLogitNormal(; bounds),
              background_κz = ScaledLogitNormal(bounds=(5e-5, 2e-4)))

    Nparams = length(priors)

    free_parameters = FreeParameters(priors)
    calibration = InverseProblem(observation, simulation, free_parameters)
    batched_calibration = InverseProblem(batched_observations, batched_simulation, free_parameters)

    #####
    ##### Test EKI
    #####

    for eki in [EnsembleKalmanInversion(calibration; pseudo_stepping=ConstantConvergence(0.9)),
                EnsembleKalmanInversion(batched_calibration; pseudo_stepping=ConstantConvergence(0.9))]

        batch_str = string("(Nbatch = ", size(eki.inverse_problem.simulation.model.grid, 2), ")")
        @testset "EnsembleKalmanInversions construction and iteration tests $batch_str" begin
            @info "  Testing EnsembleKalmanInversion construcation and basic iteration $batch_str..."

            iterations = 5
            iterate!(eki; iterations = iterations, show_progress = false)

            @test length(eki.iteration_summaries) == iterations + 1
            @test eki.iteration == iterations

            θθ_std_arr = sqrt.(hcat(diag.(getproperty.(eki.iteration_summaries, :ensemble_cov))...))
            @test size(θθ_std_arr) == (Nparams, iterations + 1)

            # Test that parameters stay within bounds
            parameters = eki.iteration_summaries[0].parameters
            convective_κzs = getproperty.(parameters, :convective_κz)
            @test all(convective_κzs .> lower_bound)
            @test all(convective_κzs .< upper_bound)

            # Test that parameters change
            @test convective_κzs[1] != convective_κzs[2]

            iterate!(eki; iterations = 1, show_progress = false)

            @test length(eki.iteration_summaries) == iterations + 2
            @test eki.iteration == iterations + 1
        end

        #####
        ##### Test Resampler
        #####

        @testset "Resampler tests $batch_str" begin
            @info "  Testing resampling and NaN handling $batch_str"

            # Test resample!
            resampler = Resampler(acceptable_failure_fraction = 1.0,
                                  distribution = FullEnsembleDistribution())

            θ = rand(Nparams, Nensemble)
            θ1 = deepcopy(θ[:, 1])
            θ2 = deepcopy(θ[:, 2])
            θ3 = deepcopy(θ[:, 3])

            # Fake a forward map output with NaNs
            norm_exceeds_median = NormExceedsMedian(Inf)
            G = inverting_forward_map(eki.inverse_problem, θ)
            view(G, :, 2) .= NaN
            @test any(isnan.(G)) == true

            @test sum(norm_exceeds_median(G)) == 1
            @test norm_exceeds_median(G)[1] == false
            @test norm_exceeds_median(G)[2] == true
            @test norm_exceeds_median(G)[3] == false

            resample!(resampler, θ, G, eki)

            @test sum(norm_exceeds_median(G)) == 0

            @test any(isnan.(G)) == false
            @test θ[:, 1] == θ1
            @test θ[:, 2] != θ2
            @test θ[:, 3] == θ3

            # Test that model fields get overwritten without NaN
            G = inverting_forward_map(eki.inverse_problem, θ)

            # Particle 2
            view(G, :, 2) .= NaN
            model = eki.inverse_problem.simulation.model
            time_series_collector = eki.inverse_problem.time_series_collector

            # Fill one batch...
            for field_name in keys(fields(model))
                field = fields(model)[field_name]
                collector = time_series_collector.field_time_serieses[field_name]

                field2 = view(parent(field), 2, 1, :)
                collector2 = view(parent(collector), 2, 1, :, :)
                fill!(field2, NaN)
                fill!(collector2, NaN)
            end

            @test any(isnan.(model.tracers.b))
            @test any(isnan.(time_series_collector.field_time_serieses.b))
            @test any(isnan.(G))

            resample!(resampler, θ, G, eki)

            @test !any(isnan.(model.tracers.b))
            @test !any(isnan.(time_series_collector.field_time_serieses.b))

            # Resample all particles, not just failed ones

            resampler = Resampler(acceptable_failure_fraction = 1.0,
                                  only_failed_particles = false,
                                  distribution = FullEnsembleDistribution())

            θ = rand(Nparams, Nensemble)
            θ1 = deepcopy(θ[:, 1])
            θ2 = deepcopy(θ[:, 2])
            θ3 = deepcopy(θ[:, 3])

            # Fake a forward map output with NaNs
            G = inverting_forward_map(eki.inverse_problem, θ)
            Gcopy = deepcopy(G)

            resample!(resampler, θ, G, eki)
            @test G != Gcopy
            @test θ[:, 1] != θ1
            @test θ[:, 2] != θ2
            @test θ[:, 3] != θ3

            # Resample particles with SuccessfulEnsembleDistribution.
            # NaN out 2 or 3 columns so that all particles end up identical
            # after resampling.

            resampler = Resampler(acceptable_failure_fraction = 1.0,
                                  only_failed_particles = false,
                                  distribution = SuccessfulEnsembleDistribution())

            θ = rand(Nparams, Nensemble)
            θ1 = deepcopy(θ[:, 1])
            θ2 = deepcopy(θ[:, 2])
            θ3 = deepcopy(θ[:, 3])

            # Fake a forward map output with NaNs
            G = inverting_forward_map(eki.inverse_problem, θ)
            view(G, :, 1) .= NaN
            view(G, :, 2) .= NaN

            @test sum(norm_exceeds_median(G)) == 2
            @test norm_exceeds_median(G)[1]
            @test norm_exceeds_median(G)[2]
            @test !(norm_exceeds_median(G)[3])

            resample!(resampler, θ, G, eki)

            @test !any(isnan.(G))
            @test θ[:, 1] != θ3
            @test θ[:, 2] != θ3
            @test θ[:, 3] != θ3
        end
    end
end