text:
puzzleinput = readlines("inputs/day22.txt")

Deck = Vector{Int}

function readinput(lines)

    player₁ = map(lines[2:26]) do s
        parse(Int, s)
    end
    player₂ = map(lines[29:53]) do s
        parse(Int, s)
    end

    return player₁, player₂
end

# First Half

function combat!(player₁::Deck, player₂::Deck)

    round = 1
    while !(isempty(player₁) || isempty(player₂))
        # Games end wheh the deck of one of the players is empty

        println("--- Round $round --")
        println("Player 1's deck: $player₁")
        println("Player 2's deck: $player₂")
        draw₁ = popfirst!(player₁)
        draw₂ = popfirst!(player₂)
        println("Plater 1 plays: $draw₁")
        println("Plater 2 plays: $draw2")
        if draw₁ > draw₂
            player₁ = [player₁; draw₁; draw₂]
            println("Player 1 wins the round!")
        else
            player₂ = [player₂; draw₂; draw₁]
            println("Player 2 wins the round!")
        end
        round += 1
        println()
    end

    println("== Post-game results ==")
    println("Player 1's deck: $player₁")
    println("Player 2's deck: $player₂")
    winner = isempty(player₁) ? player₂ : player₁
    score = sum(reverse(collect(1:length(winner))) .* winner)
    println("Score = $score")
end

player₁, player₂ = readinput(puzzleinput)
combat!(player1, player2)

# Second Half

function recursivecombat!(player₁::Deck, player₂::Deck)

    function playgame(calling_game, p₁, p₂)

        local game = calling_game + 1

        round = 1
        local history = []
        while !(isempty(p₁) || isempty(p₂))
            if (p₁, p₂) ∈ history
                # Check if there was a previous round in this game that had
                # exactly the same cards in the same order in the same players' decks

                return 1 # If the case, player 1 wins
            end
            push!(history, (deepcopy(p₁), deepcopy(p₂)))

            d₁ = popfirst!(p₁)
            d₂ = popfirst!(p₂)

            winner = if length(p1) >= d₁ && length(p₂) >= d₂
                # If both players have at least as many cards remaining in their
                # deck as the value of the card they just drew, the winner of
                # the round is determined by playing a new game of Recursive Combat

                local p₁copy = deepcopy(p1[1:d₁])
                local p₂copy = deepcopy(p₂[1:d₂])
                # To play a sub-game of Recursive Combat, each player creates
                # a new deck by making a copy of the next cards in their deck
                # (the quantity of cards copied is equal to the number on the card
                #  they drew to trigger the sub-game).

                playgame(game, p₁copy, p₂copy)
                # The winner is the winner of the sub-game

            else
                # The winner is the player with the higher-value card

                d₁ > d₂ ? 1 : 2
            end

            if winner == 1
                p₁ = [p₁; d₁; d₂]
            else
                p₂ = [p₂; d₂; d₁]
            end
            round += 1
        end

        winner = isempty(p₁) ? 2 : 1
        if game > 1
            return winner
        else
            println()
            println()
            println("== Post-game results ==")
            println("Player 1's deck: $p₁")
            println("Player 2's deck: $p₂")
            wdeck = winner == 1 ? p₁ : p₂
            score = sum(reverse(collect(1:length(wdeck))) .* wdeck)
            println("Score = $score")
        end
    end

    game = 0
    playgame(game, player₁, player₂)

end

player₁, player₂ = readinput(puzzleinput)
@time recursivecombat!(player₁, player₂)


text:
using EconPDEs, Distributions

Base.@kwdef struct WachterModel{T<: Distribution}
    # consumption process parameters
    μ::Float64 = 0.025
    σ::Float64 = 0.02

    # idiosyncratic
    λbar::Float64 = 0.0355
    κλ::Float64 = 0.08
    νλ::Float64 = 0.067
    ZDistribution::T = Normal(-0.4, 0.25) #from Ian Martin higher order cumulants paper

    # utility parameters
    ρ::Float64 = 0.012
    γ::Float64 = 3.0
    ψ::Float64 = 1.1
    ϕ::Float64 = 2.6
end

function (m::WachterModel)(state::NamedTuple, y::NamedTuple)
    (; μ, σ, λbar, κλ, νλ, ZDistribution, ρ, γ, ψ, ϕ) = m
    (; λ) = state
    (; p, pλ_up, pλ_down, pλλ) = y

    # Drift and volatility of λ, p
    μλ = κλ * (λbar - λ)
    σλ = νλ * sqrt(λ)
    pλ = (μλ >= 0) ? pλ_up : pλ_down
    μp = pλ / p * μλ + 0.5 * pλλ / p * σλ^2
    σp_Zλ = pλ / p * σλ

    # Market price of risk
    κ_Zc = γ * σ
    κ_Zλ = (γ * ψ - 1) / (ψ - 1) * σp_Zλ
    η =  λ * (mgf(ZDistribution, 1) - 1 + mgf(ZDistribution, -γ) - 1 - (mgf(ZDistribution, 1 - γ) - 1))

    # Interest rate r
    r = ρ + μ / ψ - (1 + 1 / ψ) / 2 * γ * σ^2 - λ * (mgf(ZDistribution, -γ) - 1) + (1 / ψ - γ) / (1 - γ) * λ * (mgf(ZDistribution, 1 - γ) - 1) - (γ * ψ - 1) / (2 * (ψ - 1)) * σp_Zλ^2

    # Market Pricing
    pt = - p * (1 / p  + μ + μp + λ * (mgf(ZDistribution, 1) - 1) - r - κ_Zc * σ - κ_Zλ * σp_Zλ - η)

    return (; pt)
end

m = WachterModel()
λn = 30
stategrid = OrderedDict(:λ => range(0.0, 0.1, length = λn))
yend =  OrderedDict(:p => ones(λn))
result = pdesolve(m, stategrid, yend)
@assert result.residual_norm <= 1e-5

#========================================================================================

Solve for levered equity claim

========================================================================================#

function pde_levered(m, state, y, r, κ_Zc, κ_Zλ)
    μ = m.μ ; σ = m.σ ; λbar = m.λbar ; κλ = m.κλ ; νλ = m.νλ ; ZDistribution = m.ZDistribution ; ρ = m.ρ ; γ = m.γ ; ψ = m.ψ ; ϕ = m.ϕ
    λ = state.λ
    pe, peλ_up, peλ_down, peλλ = y.pe, y.peλ, y.peλλ
    μλ = κλ * (λbar - λ)
    σλ = νλ * sqrt(λ)
    peλ = (μλ >= 0) ? peλ_up : peλ_down
    μpe = peλ / pe * μλ + 0.5 * peλλ / pe * σλ^2
    σpe_Zλ = peλ / pe * σλ
    η =  λ * (mgf(ZDistribution, ϕ) - 1 + mgf(ZDistribution, -γ) - 1 - (mgf(ZDistribution, ϕ - γ) - 1))
    pet = pe * (1 / pe  + μ + μpe + λ * (mgf(ZDistribution, ϕ) - 1) - r(λ) - κ_Zc(λ) * ϕ * σ - κ_Zλ(λ) * σpe_Zλ - η)
    return (; pet)
end
#using Interpolations
#r = interpolate((state[:λ],), result[:r], Gridded(Linear()))
#κ_Zc = interpolate((state[:λ],), result[:κ_Zc], Gridded(Linear()))
#κ_Zλ = interpolate((state[:λ],), result[:κ_Zλ], Gridded(Linear()))
#y2, _, distance =  pdesolve((state, y) -> pde_levered(m, state, y, r, κ_Zc, κ_Zλ), stategrid, OrderedDict(:pe => y[:p]))

ext: jl


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
Support for sockets. Provides [`IPAddr`](@ref) and subtypes, [`TCPSocket`](@ref), and [`UDPSocket`](@ref).
"""
module Sockets

export
    accept,
    bind,
    connect,
    getaddrinfo,
    getalladdrinfo,
    getnameinfo,
    getipaddr,
    getipaddrs,
    getpeername,
    getsockname,
    listen,
    listenany,
    recv,
    recvfrom,
    send,
    TCPSocket,
    UDPSocket,
    @ip_str,
    IPAddr,
    IPv4,
    IPv6

import Base: isless, show, print, parse, bind, convert, isreadable, iswritable, alloc_buf_hook, _uv_hook_close

using Base: LibuvStream, LibuvServer, PipeEndpoint, @handle_as, uv_error, associate_julia_struct, uvfinalize,
    notify_error, uv_req_data, uv_req_set_data, preserve_handle, unpreserve_handle, _UVError, IOError,
    eventloop, StatusUninit, StatusInit, StatusConnecting, StatusOpen, StatusClosing, StatusClosed, StatusActive,
    preserve_handle, unpreserve_handle, iolock_begin, iolock_end,
    uv_status_string, check_open, OS_HANDLE, RawFD,
    UV_EINVAL, UV_ENOMEM, UV_ENOBUFS, UV_EAGAIN, UV_ECONNABORTED, UV_EADDRINUSE, UV_EACCES, UV_EADDRNOTAVAIL,
    UV_EAI_ADDRFAMILY, UV_EAI_AGAIN, UV_EAI_BADFLAGS,
    UV_EAI_BADHINTS, UV_EAI_CANCELED, UV_EAI_FAIL,
    UV_EAI_FAMILY, UV_EAI_NODATA, UV_EAI_NONAME,
    UV_EAI_OVERFLOW, UV_EAI_PROTOCOL, UV_EAI_SERVICE,
    UV_EAI_SOCKTYPE, UV_EAI_MEMORY

include("IPAddr.jl")
include("addrinfo.jl")

"""
    TCPSocket(; delay=true)

Open a TCP socket using libuv. If `delay` is true, libuv delays creation of the
socket's file descriptor till the first [`bind`](@ref) call. `TCPSocket` has various
fields to denote the state of the socket as well as its send/receive buffers.
"""
mutable struct TCPSocket <: LibuvStream
    handle::Ptr{Cvoid}
    status::Int
    buffer::IOBuffer
    cond::Base.ThreadSynchronizer
    readerror::Any
    sendbuf::Union{IOBuffer, Nothing}
    lock::ReentrantLock # advisory lock
    throttle::Int

    function TCPSocket(handle::Ptr{Cvoid}, status)
        tcp = new(
                handle,
                status,
                PipeBuffer(),
                Base.ThreadSynchronizer(),
                nothing,
                nothing,
                ReentrantLock(),
                Base.DEFAULT_READ_BUFFER_SZ)
        associate_julia_struct(tcp.handle, tcp)
        finalizer(uvfinalize, tcp)
        return tcp
    end
end

# kw arg "delay": if true, libuv delays creation of the socket fd till the first bind call
function TCPSocket(; delay=true)
    tcp = TCPSocket(Libc.malloc(Base._sizeof_uv_tcp), StatusUninit)
    af_spec = delay ? 0 : 2   # AF_UNSPEC is 0, AF_INET is 2
    iolock_begin()
    err = ccall(:uv_tcp_init_ex, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Cuint),
                eventloop(), tcp.handle, af_spec)
    uv_error("failed to create tcp socket", err)
    tcp.status = StatusInit
    iolock_end()
    return tcp
end

function TCPSocket(fd::OS_HANDLE)
    tcp = TCPSocket()
    iolock_begin()
    err = ccall(:uv_tcp_open, Int32, (Ptr{Cvoid}, OS_HANDLE), pipe.handle, fd)
    uv_error("tcp_open", err)
    tcp.status = StatusOpen
    iolock_end()
    return tcp
end
if OS_HANDLE != RawFD
    TCPSocket(fd::RawFD) = TCPSocket(Libc._get_osfhandle(fd))
end

mutable struct TCPServer <: LibuvServer
    handle::Ptr{Cvoid}
    status::Int
    cond::Base.ThreadSynchronizer

    function TCPServer(handle::Ptr{Cvoid}, status)
        tcp = new(
            handle,
            status,
            Base.ThreadSynchronizer())
        associate_julia_struct(tcp.handle, tcp)
        finalizer(uvfinalize, tcp)
        return tcp
    end
end

# Keyword arg "delay": if true, libuv delays creation of socket fd till bind.
# It can be set to false if there is a need to set socket options before
# further calls to `bind` and `listen`, e.g. `SO_REUSEPORT`.
function TCPServer(; delay=true)
    tcp = TCPServer(Libc.malloc(Base._sizeof_uv_tcp), StatusUninit)
    af_spec = delay ? 0 : 2   # AF_UNSPEC is 0, AF_INET is 2
    iolock_begin()
    err = ccall(:uv_tcp_init_ex, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Cuint),
                eventloop(), tcp.handle, af_spec)
    uv_error("failed to create tcp server", err)
    tcp.status = StatusInit
    iolock_end()
    return tcp
end

isreadable(io::TCPSocket) = isopen(io) || bytesavailable(io) > 0
iswritable(io::TCPSocket) = isopen(io) && io.status != StatusClosing

"""
    accept(server[, client])

Accepts a connection on the given server and returns a connection to the client. An
uninitialized client stream may be provided, in which case it will be used instead of
creating a new stream.
"""
accept(server::TCPServer) = accept(server, TCPSocket())

function accept(callback, server::LibuvServer)
    task = @async try
            while true
                client = accept(server)
                callback(client)
            end
        catch ex
            # accept below may explicitly throw UV_ECONNABORTED:
            # filter that out since we expect that error
            if !(ex isa IOError && ex.code == UV_ECONNABORTED) || isopen(server)
                rethrow()
            end
        end
    return task # caller is responsible for checking for errors
end

# UDP
"""
    UDPSocket()

Open a UDP socket using libuv. `UDPSocket` has various
fields to denote the state of the socket.
"""
mutable struct UDPSocket <: LibuvStream
    handle::Ptr{Cvoid}
    status::Int
    recvnotify::Base.ThreadSynchronizer
    cond::Base.ThreadSynchronizer

    function UDPSocket(handle::Ptr{Cvoid}, status)
        cond = Base.ThreadSynchronizer()
        udp = new(handle, status, Base.ThreadSynchronizer(cond.lock), cond)
        associate_julia_struct(udp.handle, udp)
        finalizer(uvfinalize, udp)
        return udp
    end
end
function UDPSocket()
    this = UDPSocket(Libc.malloc(Base._sizeof_uv_udp), StatusUninit)
    iolock_begin()
    err = ccall(:uv_udp_init, Cint, (Ptr{Cvoid}, Ptr{Cvoid}),
                eventloop(), this.handle)
    uv_error("failed to create udp socket", err)
    this.status = StatusInit
    iolock_end()
    return this
end

show(io::IO, stream::UDPSocket) = print(io, typeof(stream), "(", uv_status_string(stream), ")")

function _uv_hook_close(sock::UDPSocket)
    sock.handle = C_NULL
    lock(sock.cond)
    try
        sock.status = StatusClosed
        notify(sock.cond)
        notify_error(sock.recvnotify, EOFError())
    finally
        unlock(sock.cond)
    end
    nothing
end

# Disables dual stack mode.
const UV_TCP_IPV6ONLY = 1

# Disables dual stack mode. Only available when using ipv6 binf
const UV_UDP_IPV6ONLY = 1

# Indicates message was truncated because read buffer was too small. The
# remainder was discarded by the OS.
const UV_UDP_PARTIAL = 2

# Indicates if SO_REUSEADDR will be set when binding the handle in uv_udp_bind. This sets
# the SO_REUSEPORT socket flag on the BSDs and OS X. On other Unix platforms, it sets the
# SO_REUSEADDR flag. What that means is that multiple threads or processes can bind to the
# same address without error (provided they all set the flag) but only the last one to bind
# will receive any traffic, in effect "stealing" the port from the previous listener.
const UV_UDP_REUSEADDR = 4

##

function _bind(sock::TCPServer, host::Union{IPv4, IPv6}, port::UInt16, flags::UInt32=UInt32(0))
    host_in = Ref(hton(host.host))
    return ccall(:jl_tcp_bind, Int32, (Ptr{Cvoid}, UInt16, Ptr{Cvoid}, Cuint, Cint),
            sock, hton(port), host_in, flags, host isa IPv6)
end

function _bind(sock::UDPSocket, host::Union{IPv4, IPv6}, port::UInt16, flags::UInt32=UInt32(0))
    host_in = Ref(hton(host.host))
    return ccall(:jl_udp_bind, Int32, (Ptr{Cvoid}, UInt16, Ptr{Cvoid}, Cuint, Cint),
            sock, hton(port), host_in, flags, host isa IPv6)
end

"""
    bind(socket::Union{UDPSocket, TCPSocket}, host::IPAddr, port::Integer; ipv6only=false, reuseaddr=false, kws...)

Bind `socket` to the given `host:port`. Note that `0.0.0.0` will listen on all devices.

* The `ipv6only` parameter disables dual stack mode. If `ipv6only=true`, only an IPv6 stack is created.
* If `reuseaddr=true`, multiple threads or processes can bind to the same address without error
  if they all set `reuseaddr=true`, but only the last to bind will receive any traffic.
"""
function bind(sock::Union{TCPServer, UDPSocket}, host::IPAddr, port::Integer; ipv6only = false, reuseaddr = false, kws...)
    if sock.status != StatusInit
        error("$(typeof(sock)) is not in initialization state")
    end
    flags = 0
    if isa(host, IPv6) && ipv6only
        flags |= isa(sock, UDPSocket) ? UV_UDP_IPV6ONLY : UV_TCP_IPV6ONLY
    end
    if isa(sock, UDPSocket) && reuseaddr
        flags |= UV_UDP_REUSEADDR
    end
    iolock_begin()
    err = _bind(sock, host, UInt16(port), UInt32(flags))
    if err < 0
        iolock_end()
        if err != UV_EADDRINUSE && err != UV_EACCES && err != UV_EADDRNOTAVAIL
            #TODO: this codepath is not currently tested
            throw(_UVError("bind", err))
        else
            return false
        end
    end
    sock.status = StatusOpen
    isa(sock, UDPSocket) && setopt(sock; kws...)
    iolock_end()
    return true
end

bind(sock::TCPServer, addr::InetAddr) = bind(sock, addr.host, addr.port)

"""
    setopt(sock::UDPSocket; multicast_loop=nothing, multicast_ttl=nothing, enable_broadcast=nothing, ttl=nothing)

Set UDP socket options.

* `multicast_loop`: loopback for multicast packets (default: `true`).
* `multicast_ttl`: TTL for multicast packets (default: `nothing`).
* `enable_broadcast`: flag must be set to `true` if socket will be used for broadcast
  messages, or else the UDP system will return an access error (default: `false`).
* `ttl`: Time-to-live of packets sent on the socket (default: `nothing`).
"""
function setopt(sock::UDPSocket; multicast_loop=nothing, multicast_ttl=nothing, enable_broadcast=nothing, ttl=nothing)
    iolock_begin()
    if sock.status == StatusUninit
        error("Cannot set options on uninitialized socket")
    end
    if multicast_loop !== nothing
        uv_error("multicast_loop", ccall(:uv_udp_set_multicast_loop, Cint, (Ptr{Cvoid}, Cint), sock.handle, multicast_loop) < 0)
    end
    if multicast_ttl !== nothing
        uv_error("multicast_ttl", ccall(:uv_udp_set_multicast_ttl, Cint, (Ptr{Cvoid}, Cint), sock.handle, multicast_ttl))
    end
    if enable_broadcast !== nothing
        uv_error("enable_broadcast", ccall(:uv_udp_set_broadcast, Cint, (Ptr{Cvoid}, Cint), sock.handle, enable_broadcast))
    end
    if ttl !== nothing
        uv_error("ttl", ccall(:uv_udp_set_ttl, Cint, (Ptr{Cvoid}, Cint), sock.handle, ttl))
    end
    iolock_end()
    nothing
end

"""
    recv(socket::UDPSocket)

Read a UDP packet from the specified socket, and return the bytes received. This call blocks.
"""
function recv(sock::UDPSocket)
    addr, data = recvfrom(sock)
    return data
end

"""
    recvfrom(socket::UDPSocket) -> (host_port, data)

Read a UDP packet from the specified socket, returning a tuple of `(host_port, data)`, where
`host_port` will be an InetAddr{IPv4} or InetAddr{IPv6}, as appropriate.

!!! compat "Julia 1.3"
    Prior to Julia version 1.3, the first returned value was an address (`IPAddr`).
    In version 1.3 it was changed to an `InetAddr`.
"""
function recvfrom(sock::UDPSocket)
    iolock_begin()
    # If the socket has not been bound, it will be bound implicitly to ::0 and a random port
    if sock.status != StatusInit && sock.status != StatusOpen && sock.status != StatusActive
        error("UDPSocket is not initialized and open")
    end
    if ccall(:uv_is_active, Cint, (Ptr{Cvoid},), sock.handle) == 0
        err = ccall(:uv_udp_recv_start, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
                    sock, Base.uv_jl_alloc_buf::Ptr{Cvoid}, uv_jl_recvcb::Ptr{Cvoid})
        uv_error("recv_start", err)
    end
    sock.status = StatusActive
    lock(sock.recvnotify)
    iolock_end()
    try
        From = Union{InetAddr{IPv4}, InetAddr{IPv6}}
        Data = Vector{UInt8}
        from, data = wait(sock.recvnotify)::Tuple{From, Data}
        return (from, data)
    finally
        unlock(sock.recvnotify)
    end
end

alloc_buf_hook(sock::UDPSocket, size::UInt) = (Libc.malloc(size), size) # size is always 64k from libuv

function uv_recvcb(handle::Ptr{Cvoid}, nread::Cssize_t, buf::Ptr{Cvoid}, addr::Ptr{Cvoid}, flags::Cuint)
    sock = @handle_as handle UDPSocket
    lock(sock.recvnotify)
    try
        buf_addr = ccall(:jl_uv_buf_base, Ptr{UInt8}, (Ptr{Cvoid},), buf)
        if nread == 0 && addr == C_NULL
            Libc.free(buf_addr)
        elseif nread < 0
            Libc.free(buf_addr)
            notify_error(sock.recvnotify, _UVError("recv", nread))
        elseif flags & UV_UDP_PARTIAL > 0
            Libc.free(buf_addr)
            notify_error(sock.recvnotify, "Partial message received")
        else
            buf_size = Int(ccall(:jl_uv_buf_len, Csize_t, (Ptr{Cvoid},), buf))
            if buf_size - nread < 16384 # waste at most 16k (note: buf_size is currently always 64k)
                buf = unsafe_wrap(Array, buf_addr, nread, own=true)
            else
                buf = Vector{UInt8}(undef, nread)
                GC.@preserve buf unsafe_copyto!(pointer(buf), buf_addr, nread)
                Libc.free(buf_addr)
            end
            # need to check the address type in order to convert to a Julia IPAddr
            host = IPv4(0)
            port = UInt16(0)
            if ccall(:jl_sockaddr_is_ip4, Cint, (Ptr{Cvoid},), addr) == 1
                host = IPv4(ntoh(ccall(:jl_sockaddr_host4, UInt32, (Ptr{Cvoid},), addr)))
                port = ntoh(ccall(:jl_sockaddr_port4, UInt16, (Ptr{Cvoid},), addr))
            elseif ccall(:jl_sockaddr_is_ip6, Cint, (Ptr{Cvoid},), addr) == 1
                tmp = Ref{UInt128}(0)
                scope_id = ccall(:jl_sockaddr_host6, UInt32, (Ptr{Cvoid}, Ptr{UInt128}), addr, tmp)
                host = IPv6(ntoh(tmp[]))
                port = ntoh(ccall(:jl_sockaddr_port6, UInt16, (Ptr{Cvoid},), addr))
            end
            from = InetAddr(host, port)
            notify(sock.recvnotify, (from, buf), all=false)
        end
        if sock.status == StatusActive && isempty(sock.recvnotify)
            sock.status = StatusOpen
            ccall(:uv_udp_recv_stop, Cint, (Ptr{Cvoid},), sock)
        end
    finally
        unlock(sock.recvnotify)
    end
    nothing
end

function _send_async(sock::UDPSocket, ipaddr::Union{IPv4, IPv6}, port::UInt16, buf)
    req = Libc.malloc(Base._sizeof_uv_udp_send)
    uv_req_set_data(req, C_NULL) # in case we get interrupted before arriving at the wait call
    host_in = Ref(hton(ipaddr.host))
    err = ccall(:jl_udp_send, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, UInt16, Ptr{Cvoid}, Ptr{UInt8}, Csize_t, Ptr{Cvoid}, Cint),
            req, sock, hton(port), host_in, buf, sizeof(buf), Base.uv_jl_writecb_task::Ptr{Cvoid}, ipaddr isa IPv6)
    if err < 0
        Libc.free(req)
        uv_error("send", err)
    end
    return req
end

"""
    send(socket::UDPSocket, host::IPAddr, port::Integer, msg)

Send `msg` over `socket` to `host:port`.
"""
function send(sock::UDPSocket, ipaddr::IPAddr, port::Integer, msg)
    # If the socket has not been bound, it will be bound implicitly to ::0 and a random port
    iolock_begin()
    if sock.status != StatusInit && sock.status != StatusOpen && sock.status != StatusActive
        error("UDPSocket is not initialized and open")
    end
    uvw = _send_async(sock, ipaddr, UInt16(port), msg)
    ct = current_task()
    preserve_handle(ct)
    Base.sigatomic_begin()
    uv_req_set_data(uvw, ct)
    iolock_end()
    status = try
        Base.sigatomic_end()
        wait()::Cint
    finally
        Base.sigatomic_end()
        iolock_begin()
        ct.queue === nothing || list_deletefirst!(ct.queue, ct)
        if uv_req_data(uvw) != C_NULL
            # uvw is still alive,
            # so make sure we won't get spurious notifications later
            uv_req_set_data(uvw, C_NULL)
        else
            # done with uvw
            Libc.free(uvw)
        end
        iolock_end()
        unpreserve_handle(ct)
    end
    uv_error("send", status)
    nothing
end

#from `connect`
function uv_connectcb(conn::Ptr{Cvoid}, status::Cint)
    hand = ccall(:jl_uv_connect_handle, Ptr{Cvoid}, (Ptr{Cvoid},), conn)
    sock = @handle_as hand LibuvStream
    lock(sock.cond)
    try
        if status >= 0 # success
            if !(sock.status == StatusClosed || sock.status == StatusClosing)
                sock.status = StatusOpen
            end
        else
            sock.readerror = _UVError("connect", status) # TODO: perhaps we should not reuse readerror for this
            if !(sock.status == StatusClosed || sock.status == StatusClosing)
                ccall(:jl_forceclose_uv, Cvoid, (Ptr{Cvoid},), hand)
                sock.status = StatusClosing
            end
        end
        notify(sock.cond)
    finally
        unlock(sock.cond)
    end
    Libc.free(conn)
    nothing
end

function connect!(sock::TCPSocket, host::Union{IPv4, IPv6}, port::Integer)
    iolock_begin()
    if sock.status != StatusInit
        error("TCPSocket is not in initialization state")
    end
    if !(0 <= port <= typemax(UInt16))
        throw(ArgumentError("port out of range, must be 0 ≤ port ≤ 65535, got $port"))
    end
    host_in = Ref(hton(host.host))
    uv_error("connect", ccall(:jl_tcp_connect, Int32, (Ptr{Cvoid}, Ptr{Cvoid}, UInt16, Ptr{Cvoid}, Cint),
                              sock, host_in, hton(UInt16(port)), uv_jl_connectcb::Ptr{Cvoid}, host isa IPv6))
    sock.status = StatusConnecting
    iolock_end()
    nothing
end

connect!(sock::TCPSocket, addr::InetAddr) = connect!(sock, addr.host, addr.port)

function wait_connected(x::LibuvStream)
    iolock_begin()
    check_open(x)
    isopen(x) || x.readerror === nothing || throw(x.readerror)
    preserve_handle(x)
    lock(x.cond)
    try
        while x.status == StatusConnecting
            iolock_end()
            wait(x.cond)
            unlock(x.cond)
            iolock_begin()
            lock(x.cond)
        end
        isopen(x) || x.readerror === nothing || throw(x.readerror)
    finally
        unlock(x.cond)
        unpreserve_handle(x)
    end
    iolock_end()
    nothing
end

# Default Host to localhost

"""
    connect([host], port::Integer) -> TCPSocket

Connect to the host `host` on port `port`.
"""
connect(sock::TCPSocket, port::Integer) = connect(sock, localhost, port)
connect(port::Integer) = connect(localhost, port)

# Valid connect signatures for TCP
connect(host::AbstractString, port::Integer) = connect(TCPSocket(), host, port)
connect(addr::IPAddr, port::Integer) = connect(TCPSocket(), addr, port)
connect(addr::InetAddr) = connect(TCPSocket(), addr)

function connect!(sock::TCPSocket, host::AbstractString, port::Integer)
    if sock.status != StatusInit
        error("TCPSocket is not in initialization state")
    end
    ipaddr = getaddrinfo(host)
    connect!(sock, ipaddr, port)
    return sock
end

function connect(sock::LibuvStream, args...)
    connect!(sock, args...)
    wait_connected(sock)
    return sock
end

"""
    nagle(socket::Union{TCPServer, TCPSocket}, enable::Bool)

Enables or disables Nagle's algorithm on a given TCP server or socket.
"""
function nagle(sock::Union{TCPServer, TCPSocket}, enable::Bool)
    # disable or enable Nagle's algorithm on all OSes
    Sockets.iolock_begin()
    Sockets.check_open(sock)
    err = ccall(:uv_tcp_nodelay, Cint, (Ptr{Cvoid}, Cint), sock.handle, Cint(!enable))
    # TODO: check err
    Sockets.iolock_end()
    return err
end

"""
    quickack(socket::Union{TCPServer, TCPSocket}, enable::Bool)

On Linux systems, the TCP_QUICKACK is disabled or enabled on `socket`.
"""
function quickack(sock::Union{TCPServer, TCPSocket}, enable::Bool)
    Sockets.iolock_begin()
    Sockets.check_open(sock)
    @static if Sys.islinux()
        # tcp_quickack is a linux only option
        if ccall(:jl_tcp_quickack, Cint, (Ptr{Cvoid}, Cint), sock.handle, Cint(enable)) < 0
            @warn "Networking unoptimized ( Error enabling TCP_QUICKACK : $(Libc.strerror(Libc.errno())) )" maxlog=1
        end
    end
    Sockets.iolock_end()
    nothing
end

##

const BACKLOG_DEFAULT = 511

"""
    listen([addr, ]port::Integer; backlog::Integer=BACKLOG_DEFAULT) -> TCPServer

Listen on port on the address specified by `addr`.
By default this listens on `localhost` only.
To listen on all interfaces pass `IPv4(0)` or `IPv6(0)` as appropriate.
`backlog` determines how many connections can be pending (not having
called [`accept`](@ref)) before the server will begin to
reject them. The default value of `backlog` is 511.
"""
function listen(addr; backlog::Integer=BACKLOG_DEFAULT)
    sock = TCPServer()
    bind(sock, addr) || error("cannot bind to port; may already be in use or access denied")
    listen(sock; backlog=backlog)
    return sock
end
listen(port::Integer; backlog::Integer=BACKLOG_DEFAULT) = listen(localhost, port; backlog=backlog)
listen(host::IPAddr, port::Integer; backlog::Integer=BACKLOG_DEFAULT) = listen(InetAddr(host, port); backlog=backlog)

function listen(sock::LibuvServer; backlog::Integer=BACKLOG_DEFAULT)
    uv_error("listen", trylisten(sock))
    return sock
end

# from `listen`
function uv_connectioncb(stream::Ptr{Cvoid}, status::Cint)
    sock = @handle_as stream LibuvServer
    lock(sock.cond)
    try
        if status >= 0
            notify(sock.cond)
        else
            notify_error(sock.cond, _UVError("connection", status))
        end
    finally
        unlock(sock.cond)
    end
    nothing
end

function trylisten(sock::LibuvServer; backlog::Integer=BACKLOG_DEFAULT)
    iolock_begin()
    check_open(sock)
    err = ccall(:uv_listen, Cint, (Ptr{Cvoid}, Cint, Ptr{Cvoid}),
                sock, backlog, uv_jl_connectioncb::Ptr{Cvoid})
    sock.status = StatusActive
    iolock_end()
    return err
end

##

function accept_nonblock(server::TCPServer, client::TCPSocket)
    iolock_begin()
    if client.status != StatusInit
        error("client TCPSocket is not in initialization state")
    end
    err = ccall(:uv_accept, Int32, (Ptr{Cvoid}, Ptr{Cvoid}), server.handle, client.handle)
    if err == 0
        client.status = StatusOpen
    end
    iolock_end()
    return err
end

function accept_nonblock(server::TCPServer)
    client = TCPSocket()
    uv_error("accept", accept_nonblock(server, client))
    return client
end

function accept(server::LibuvServer, client::LibuvStream)
    iolock_begin()
    if server.status != StatusActive && server.status != StatusClosing && server.status != StatusClosed
        throw(ArgumentError("server not connected, make sure \"listen\" has been called"))
    end
    while isopen(server)
        err = accept_nonblock(server, client)
        if err == 0
            iolock_end()
            return client
        elseif err != UV_EAGAIN
            uv_error("accept", err)
        end
        preserve_handle(server)
        lock(server.cond)
        iolock_end()
        try
            wait(server.cond)
        finally
            unlock(server.cond)
            unpreserve_handle(server)
        end
        iolock_begin()
    end
    uv_error("accept", UV_ECONNABORTED)
    nothing
end

## Utility functions

const localhost = ip"127.0.0.1"

"""
    listenany([host::IPAddr,] port_hint) -> (UInt16, TCPServer)

Create a `TCPServer` on any port, using hint as a starting point. Returns a tuple of the
actual port that the server was created on and the server itself.
"""
function listenany(host::IPAddr, default_port)
    addr = InetAddr(host, default_port)
    while true
        sock = TCPServer()
        if bind(sock, addr) && trylisten(sock) == 0
            if default_port == 0
                _addr, port = getsockname(sock)
                return (port, sock)
            end
            return (addr.port, sock)
        end
        close(sock)
        addr = InetAddr(addr.host, addr.port + 1)
        if addr.port == default_port
            error("no ports available")
        end
    end
end

listenany(default_port) = listenany(localhost, default_port)

"""
    getsockname(sock::Union{TCPServer, TCPSocket}) -> (IPAddr, UInt16)

Get the IP address and port that the given socket is bound to.
"""
getsockname(sock::Union{TCPSocket, TCPServer}) = _sockname(sock, true)

"""
    getpeername(sock::TCPSocket) -> (IPAddr, UInt16)

Get the IP address and port of the remote endpoint that the given
socket is connected to. Valid only for connected TCP sockets.
"""
getpeername(sock::TCPSocket) = _sockname(sock, false)

function _sockname(sock, self=true)
    rport = Ref{Cushort}(0)
    raddress = zeros(UInt8, 16)
    rfamily = Ref{Cuint}(0)

    iolock_begin()
    if self
        r = ccall(:jl_tcp_getsockname, Int32,
                (Ptr{Cvoid}, Ref{Cushort}, Ptr{Cvoid}, Ref{Cuint}),
                sock.handle, rport, raddress, rfamily)
    else
        r = ccall(:jl_tcp_getpeername, Int32,
                (Ptr{Cvoid}, Ref{Cushort}, Ptr{Cvoid}, Ref{Cuint}),
                sock.handle, rport, raddress, rfamily)
    end
    iolock_end()
    uv_error("cannot obtain socket name", r)
    port = ntoh(rport[])
    af_inet6 = @static if Sys.iswindows() # AF_INET6 in <sys/socket.h>
        23
    elseif Sys.isapple()
        30
    elseif Sys.KERNEL ∈ (:FreeBSD, :DragonFly)
        28
    elseif Sys.KERNEL ∈ (:NetBSD, :OpenBSD)
        24
    else
        10
    end

    if rfamily[] == 2 # AF_INET
        addrv4 = raddress[1:4]
        naddr = ntoh(unsafe_load(Ptr{Cuint}(pointer(addrv4)), 1))
        addr = IPv4(naddr)
    elseif rfamily[] == af_inet6
        naddr = ntoh(unsafe_load(Ptr{UInt128}(pointer(raddress)), 1))
        addr = IPv6(naddr)
    else
        error(string("unsupported address family: ", rfamily[]))
    end
    return addr, port
end

# domain sockets

include("PipeServer.jl")

# libuv callback handles

function __init__()
    global uv_jl_getaddrinfocb = @cfunction(uv_getaddrinfocb, Cvoid, (Ptr{Cvoid}, Cint, Ptr{Cvoid}))
    global uv_jl_getnameinfocb = @cfunction(uv_getnameinfocb, Cvoid, (Ptr{Cvoid}, Cint, Cstring, Cstring))
    global uv_jl_recvcb        = @cfunction(uv_recvcb, Cvoid, (Ptr{Cvoid}, Cssize_t, Ptr{Cvoid}, Ptr{Cvoid}, Cuint))
    global uv_jl_connectioncb  = @cfunction(uv_connectioncb, Cvoid, (Ptr{Cvoid}, Cint))
    global uv_jl_connectcb     = @cfunction(uv_connectcb, Cvoid, (Ptr{Cvoid}, Cint))
end

end

max_stars_repo_name: Lun-Liu/SC-by-Default-Julia

max_issues_repo_name: Lun-Liu/SC-by-Default-Julia

max_forks_repo_name: Lun-Liu/SC-by-Default-Julia


text:
using Colors
using FileIO
using GLMakie
using GeoJSON
using GeoMakie.GeoInterface
using XLSX

inzidenz_colors = [
    RGB(244/255, 225/255, 96/255),
    RGB(224/255, 164/255, 74/255),
    RGB(220/255, 117/255, 56/255),
    RGB(195/255, 75/255, 64/255),
    RGB(117/255, 47/255, 46/255)
]

read_kreise() = GeoJSON.read(read("landkreise.geojson"))

function color_for_inzidenz(x)
    if x <= 25
        inzidenz_colors[1]
    elseif x <= 50
        inzidenz_colors[2]
    elseif x <= 100
        inzidenz_colors[3]
    elseif x <= 250
        inzidenz_colors[4]
    elseif x <= 500
        inzidenz_colors[5]
    else
        RGB(0, 0, 0)
    end
end

function read_inzidenzen(rki_excel_file)
    xf = XLSX.readxlsx(rki_excel_file)
    rows = xf["LK_7-Tage-Inzidenz (fixiert)"][:][6:end, [3, end]]

    rows_with_content = filter(collect(eachrow(rows))) do row
        !ismissing(row[1])
    end
    Dict(map(rows_with_content) do row
         @sprintf("%05i", row[1]) => row[2]
    end)
end

function colors_from_rki(kreise, rki_excel_file)
    inz_dict = read_inzidenzen(rki_excel_file)
    map(kreise.features) do kreis
        ags = kreis.properties["AGS"]

        try
            color_for_inzidenz(inz_dict[ags])
        catch
            RGB(1, 1, 1) # white for not known data
        end
    end
end

function inzidenz_legend(scene::Scene, date)
    u25 = PolyElement(color = inzidenz_colors[1], strokecolor=:black)
    u50 = PolyElement(color = inzidenz_colors[2], strokecolor=:black)
    u100 = PolyElement(color = inzidenz_colors[3], strokecolor=:black)
    u250 = PolyElement(color = inzidenz_colors[4], strokecolor=:black)
    u500 = PolyElement(color = inzidenz_colors[5], strokecolor=:black)
    above = PolyElement(color = RGB(0, 0, 0), strokecolor=:black)

    leg = Legend(scene,
                 [u25, u50, u100, u250, u500, above],
                 ["<=25", "25<x<=50", "50<x<=100", "100<x<=250", "250<x<=500", "x>500"],
                 "7T-Inzidenz, $date",
                 patchsize=(35,35),
                 orientation=:vertical,
                 tellwidth=false,
                 tellheight=false,
                 halign=:center,
                 margin=(350, 0, 350, 0),
                )
end

function create_plot(kreise, colors, date)
    fig = Figure(resolution=(1500, 1920), title="11. November 2021")
    ax = fig[1,1] = Axis(fig)

    poly!(ax, kreise, color=colors, strokewidth=1, strokecolor=:white, )
    inzidenz_legend(fig.scene, date)

    display(fig)
    return fig
end

function rki_plot(rki_excel_file, date)
    kreise = read_kreise()
    colors = colors_from_rki(kreise, rki_excel_file)

    create_plot(kreise, colors, date)
end

max_stars_repo_name: EightSQ/inzidenzkarte

max_issues_repo_name: EightSQ/inzidenzkarte

max_forks_repo_name: EightSQ/inzidenzkarte


text:
using ClusterTrees

"""
    Sort objects to lie on a space filling curve
"""
function sort_sfc(points)

    ct, sz = CompScienceMeshes.boundingbox(points)
    tree = ClusterTrees.LevelledTrees.LevelledTree(ct, sz, Int[])

    smb = sz / 2^log(length(points)+1)
    for (i,pt) = enumerate(points)
        dest = (smallest_box_size=smb, target_point=pt)
        state = ClusterTrees.LevelledTrees.rootstate(tree, dest)
        ClusterTrees.update!(tree, state, i, dest) do tree, node, i
            push!(data(tree, node).values, i)
        end
    end

    sorted = Vector{Int}()
    for node in ClusterTrees.DepthFirstIterator(tree, root(tree))
        append!(sorted, data(tree,node).values)
    end

    return sorted
end


text:
"""
Residual model to pick up the slack between approximation of -∂z(wT) and the truth.
This example uses the ResidualData struct and GP.

    # temperature eq.
    # ∂t(T) + ∂x(uT) + ∂y(vT) + ∂z(wT) = 𝓀∇²T
    # --> horizontal average -->
    #
    #           diffusive  advective
    # ∂t(T) = - ∂z(wT) +  𝓀∂z(∂z(T))
    #
    # GOAL: express wT in terms of only large-scale terms wT = F(T,h,𝓀,μ)
    #  - 𝓀: diffusivity
    #  - μ: coefficient of viscosity

    # T(n+1) = T(n) + Δt * (- ∂z(wT) +  𝓀∂z(∂z(T)))
    # Gp = -∂z(wT)
    # y = T(n+1) - T(n) - Δt * ∂z(∂z(T))

    use GPR to capture the difference between the truth
        Δt(-∂z(wT))
    and the approximation
        T(n+1) - T(n) - Δt * ∂z(∂z(T))

    # target
        approx - truth
"""

# using Statistics, LinearAlgebra, Plots
# include("GP.jl")

struct ResidualData
    truth
    approx
end

# get_truth(𝒟::ProfileData, 𝒢::GP)

δ(ϕ, z) = diff(ϕ) ./ diff(z)

function construct_residual_data(filename, D; N=4)

    # return ProfileData(V, vavg, x, y, x_train, y_train, verification_set, z, zavg, t, Nt, n_train)
end

# filename = "general_strat_16_profiles.jld2"
filename = "dns_profiles.jld2"
D = 16
N = 4

𝒟_wT = construct_profile_data(filename, "wT", D; N=N)
𝒟_T  = construct_profile_data(filename, "T",  D; N=N)
Nt = 𝒟_wT.Nt
t = 𝒟_wT.t
Δt = t[2]-t[1]
κₑ = 𝒟_wT.κₑ

avgd = false
smooth = false

if avgd
    wT = 𝒟_wT.vavg
    T = 𝒟_T.vavg
    z = 𝒟_wT.zavg
    ∂wT∂z = [δ(q, z) for q in wT] # Δt(∂z(wT))
    # Δt * ∂z(∂z(T)) + T(n) - T(n+1)
    Δt∂²T∂z² = [δ( δ(q,z), z[1:end-1] )*Δt for q in T]
    approx = [κₑ*Δt∂²T∂z²[i] .+ T[i][1:end-2] .- T[i+1][1:end-2] for i in 1:(Nt-1)]
else
    wT = 𝒟_wT.v
    if smooth
        wT = smooth_window(wT)
    end

    T  = 𝒟_T.v
    z  = 𝒟_wT.z
    ∂wT∂z = [δ(wT[:,i+1], z)*Δt for i in 1:Nt-1] # Δt(∂z(wT))
    # Δt * ∂z(∂z(T)) + T(n) - T(n+1)
    Δt∂²T∂z² = [δ( δ(T[:,i],z), z[1:end-1] )*Δt for i in 1:Nt]
    approx = [κₑ.*2*Δt∂²T∂z²[i] .+ T[:,i][1:end-2] .- T[:,i+2][1:end-2] for i in 1:(Nt-2)]
end

# Δt(-∂z(wT)) = T(n+1) - T(n) - Δt * ∂z(∂z(T))
# Δt(∂z(wT)) = Δt * ∂z(∂z(T)) + T(n) - T(n+1)

# animation_set = 1:20:(Nt)
animation_set = 1:50
anim = @animate for i in animation_set
    day_string = string(floor(Int, t[i]/86400))
    p1 = scatter(∂wT∂z[i], z, label = "Δt*∂z(wT)", xlims=(-2e-6,2e-6))
    # scatter!(approx[i], z[1:end-2], legend = :topleft, label = "∂z(∂z(T)) + T(n) - T(n+1)", ylabel = "depth", title = "day " * day_string)
    scatter!(approx[i], z[1:end-2], legend = :topright, label = "Δt*κₑ*∂z(∂z(T)) + T(n) - T(n+1)", ylabel = "depth", title = "t = $(𝒟_wT.t[i])")
end
gif(anim, pwd() * "/try2_dns.gif", fps = 2)

function smooth_window(wT)
    # takes in an Nz x Nt array
    Nz, Nt = size(wT)
    smooth_wT = similar(wT)
    for i in 2:(Nz-1)
        for j in 1:Nt
            smooth_wT[i,j] = (wT[i-1,j] + wT[i+1,j])/2
        end
    end
    return smooth_wT
end

smooth_window(wT)

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
Sphere(x...) = HyperSphere(x...)

Circle(x...) = HyperSphere(x...)

widths(c::HyperSphere{N, T}) where {N, T} = Vec{N, T}(radius(c)*2)
radius(c::HyperSphere) = c.r
origin(c::HyperSphere) = c.center
minimum(c::HyperSphere{N, T}) where {N, T} = Vec{N, T}(origin(c)) - Vec{N, T}(radius(c))
maximum(c::HyperSphere{N, T}) where {N, T} = Vec{N, T}(origin(c)) + Vec{N, T}(radius(c))
function isinside(c::Circle, x::Real, y::Real)
    @inbounds ox, oy = origin(c)
    xD = abs(ox - x)
    yD = abs(oy - y)
    xD <= c.r && yD <= c.r
end

centered(S::Type{HyperSphere{N, T}}) where {N, T} = S(Vec{N,T}(0), T(0.5))
centered(::Type{T}) where {T <: HyperSphere} = centered(HyperSphere{ndims_or(T, 3), eltype_or(T, Float32)})


text:
# Do pairwise tests for effects of one model's surprisal estimates over and above the other,
# comparing SRN, GRU and LSTM after averaging surprisal over training repetitions (fully trained only)

using Distributions, CategoricalArrays

function add_mean_surp(new_col_name, surp_col_names)
    surp = surprisal[surp_col_names]
    mean_surp = Array{Union{Missing, Float64}}(undef, size(surp,1), 1)
    for i = collect(1:1:size(surp,1))
        mean_surp[i,:] .= mean(surp[i,:])
    end
    surprisal[new_col_name] = vec(mean_surp)
    surprisal
end

function benjamini(pvals,fdr=.05)
    #=Benjamini & Hochberg's (J. Royal Statistical Society B,1995) method for
    controlling the false discovery rate for multiple comparisons.

    pvals  : a list of (uncorrected) p-values for the individual comparisons
    fdr    : the desired false discovery rate. Default: .05

    Returns the indices of "pvals" corresponding to the null hypotheses that
    should be rejected (i.e., significant effects).
    =#
    nrcomp = length(pvals)
    ind    = sortperm(pvals)    # Sort indices of p-values from lowest to highest p-value
    pvals  = pvals[ind]         # Sort p-values
    sort(ind[1:findall(pvals .<= fdr*collect(1:nrcomp)/nrcomp)[end]])
end

function create_data_surp(data)
    # Create data frame with human data and average surprisal of each RNN type

    # Add average surprisals
    surprisal = add_mean_surp(:srn, [ :SRN_6470000_0_surprisal, :SRN_6470000_1_surprisal, :SRN_6470000_2_surprisal, :SRN_6470000_3_surprisal, :SRN_6470000_4_surprisal, :SRN_6470000_5_surprisal])
    surprisal = add_mean_surp(:gru, [ :GRU_6470000_0_surprisal, :GRU_6470000_1_surprisal, :GRU_6470000_2_surprisal, :GRU_6470000_3_surprisal, :GRU_6470000_4_surprisal, :GRU_6470000_5_surprisal])
    surprisal = add_mean_surp(:lstm,[:LSTM_6470000_0_surprisal,:LSTM_6470000_1_surprisal,:LSTM_6470000_2_surprisal,:LSTM_6470000_3_surprisal,:LSTM_6470000_4_surprisal,:LSTM_6470000_5_surprisal])
    surprisal = add_mean_surp(:srn_prev, [ :SRN_6470000_0_surprisal_prev, :SRN_6470000_1_surprisal_prev, :SRN_6470000_2_surprisal_prev, :SRN_6470000_3_surprisal_prev, :SRN_6470000_4_surprisal_prev, :SRN_6470000_5_surprisal_prev])
    surprisal = add_mean_surp(:gru_prev, [ :GRU_6470000_0_surprisal_prev, :GRU_6470000_1_surprisal_prev, :GRU_6470000_2_surprisal_prev, :GRU_6470000_3_surprisal_prev, :GRU_6470000_4_surprisal_prev, :GRU_6470000_5_surprisal_prev])
    surprisal = add_mean_surp(:lstm_prev,[:LSTM_6470000_0_surprisal_prev,:LSTM_6470000_1_surprisal_prev,:LSTM_6470000_2_surprisal_prev,:LSTM_6470000_3_surprisal_prev,:LSTM_6470000_4_surprisal_prev,:LSTM_6470000_5_surprisal_prev])

    # Join human data with surprisal values and normalize suprisals
    surp_df = surprisal[[:srn, :gru, :lstm, :srn_prev, :gru_prev, :lstm_prev, :item]]
    data_surp = join(data, surp_df, kind=:left, on=:item, makeunique=true)
    data_surp[[:srn, :gru, :lstm, :srn_prev, :gru_prev, :lstm_prev]] = zscore(data_surp[[:srn, :gru, :lstm, :srn_prev, :gru_prev, :lstm_prev]])
    data_surp[!,1] = CategoricalArray(data_surp[!,1]) # subject
    data_surp[!,4] = CategoricalArray(data_surp[!,1]) # item
    data_surp
end

function do_model_comparisons(lmm_srn,lmm_gru,lmm_lstm,lmm_srn_gru,lmm_srn_lstm,lmm_gru_lstm,df)
    #=
    Compare deviance of models with one versus two sets of surprisal estimates. Returns
    χ²-statistics and uncorrected p-values
    =#
    chi2_srn_gru   = deviance(lmm_gru)  - deviance(lmm_srn_gru)
    chi2_srn_lstm  = deviance(lmm_lstm) - deviance(lmm_srn_lstm)
    chi2_gru_srn   = deviance(lmm_srn)  - deviance(lmm_srn_gru)
    chi2_gru_lstm  = deviance(lmm_lstm) - deviance(lmm_gru_lstm)
    chi2_lstm_srn  = deviance(lmm_srn)  - deviance(lmm_srn_lstm)
    chi2_lstm_gru  = deviance(lmm_gru)  - deviance(lmm_gru_lstm)

    p_srn_gru  = 1-cdf(Chisq(df), chi2_srn_gru)
    p_srn_lstm = 1-cdf(Chisq(df), chi2_srn_lstm)
    p_gru_srn  = 1-cdf(Chisq(df), chi2_gru_srn)
    p_gru_lstm = 1-cdf(Chisq(df), chi2_gru_lstm)
    p_lstm_srn = 1-cdf(Chisq(df), chi2_lstm_srn)
    p_lstm_gru = 1-cdf(Chisq(df), chi2_lstm_gru)

    (Array([missing      chi2_srn_gru   chi2_srn_lstm;
           chi2_gru_srn  missing        chi2_gru_lstm;
           chi2_lstm_srn chi2_lstm_gru missing]),
     Array([missing    p_srn_gru  p_srn_lstm;
            p_gru_srn  missing    p_gru_lstm;
            p_lstm_srn p_lstm_gru missing]))
end

function print_model_comparisons(chi2,df,sign)
    # Print results table of pairwise comparisons (χ²-statistics and asterisks indicating corrected significance levels)

    function asterisks(sign,idx)
        if idx in sign[1]
            stars = "***"
        elseif idx in sign[2]
            stars = "**"
        elseif idx in sign[3]
            stars = "*"
        else
            stars = ""
        end
    end

    chi2text = " χ²(" * repr(df) * ") = "
    println("SRN over GRU :" * chi2text * repr(round(chi2[1,2],digits=3)) * asterisks(sign,3))
    println("SRN over LSTM:" * chi2text * repr(round(chi2[1,3],digits=3)) * asterisks(sign,5))
    println("GRU over SRN :" * chi2text * repr(round(chi2[2,1],digits=3)) * asterisks(sign,1))
    println("GRU over LSTM:" * chi2text * repr(round(chi2[2,3],digits=3)) * asterisks(sign,6))
    println("LSTM over SRN:" * chi2text * repr(round(chi2[3,1],digits=3)) * asterisks(sign,2))
    println("LSTM over GRU:" * chi2text * repr(round(chi2[3,2],digits=3)) * asterisks(sign,4))
end

function do_fit(data_type, surp_from)
    # Fit LMM with all factors named in 'surp_from' and other factors depending on data_type
    # ("SPR", "ET" or "EEG)

    function fixef(surp_from)
        fixef_str = ""
        for surp = surp_from
            fixef_str *= surp * "+"
        end
        fixef_str
    end

    function ranef(surp_from)
        ranef_str = ""
        for surp = surp_from
            ranef_str *= "(0+" * surp * "|subj_nr) +"
        end
        ranef_str
    end

    if data_type == "SPR"
        formula_str = "(logRT ~" * fixef(surp_from) *
                      "logRT_prev * log_freq * log_freq_prev * nr_char * nr_char_prev * word_pos + (1|subj_nr) +" * ranef(surp_from) *
                      "(0+logRT_prev|subj_nr) + (0+log_freq|subj_nr) + (0+log_freq_prev|subj_nr) + (0+nr_char|subj_nr) + (0+nr_char_prev|subj_nr) + (0+word_pos|subj_nr) + (1|item))"
    elseif data_type == "ET"
        formula_str = "(logRTfirstpass ~" * fixef(surp_from) *
                      "log_freq * log_freq_prev * nr_char * nr_char_prev * word_pos * prevfix + (1|subj_nr) +" * ranef(surp_from) *
                      "(0+log_freq|subj_nr) + (0+log_freq_prev|subj_nr) + (0+nr_char|subj_nr) + (0+nr_char_prev|subj_nr) + (0+word_pos|subj_nr) +(0+prevfix|subj_nr) + (1|item) + (0+prevfix|item))"
    elseif data_type == "EEG"
        formula_str = "(N400 ~" * fixef(surp_from) *
                      "baseline + log_freq * nr_char * word_pos+ (1|subj_nr) +" * ranef(surp_from) *
                      "(0+baseline|subj_nr) + (0+log_freq|subj_nr) + (0+nr_char|subj_nr) + (0+word_pos|subj_nr) + (1|item))"
    else error("Unknown data type"*data_type)
    end

    fit!(LinearMixedModel(eval(Meta.parse("@formula" * formula_str)), data_surp))
end
#-------------------------------------------------------------------------------------------------------------------
println("Doing pairwise model comparisons")

#=
data_SPR = read_data("SPR")
data_ET  = read_data("ET")
data_EEG = read_data("EEG")
=#

# Fit regression models with average surprisals from each RNN and from each pair
println("Running SPR analyses")
data_surp = create_data_surp(data_SPR)

lmm_srn_SPR  = do_fit("SPR", ["srn", "srn_prev"])
lmm_gru_SPR  = do_fit("SPR", ["gru", "gru_prev"])
lmm_lstm_SPR = do_fit("SPR", ["lstm","lstm_prev"])

lmm_srn_gru_SPR  = do_fit("SPR", ["srn", "srn_prev", "gru", "gru_prev"])
lmm_srn_lstm_SPR = do_fit("SPR", ["srn", "srn_prev", "lstm", "lstm_prev"])
lmm_gru_lstm_SPR = do_fit("SPR", ["gru", "gru_prev", "lstm", "lstm_prev"])

println("Running ET analyses")
data_surp = create_data_surp(data_ET)

lmm_srn_ET  = do_fit("ET", ["srn", "srn_prev"])
lmm_gru_ET  = do_fit("ET", ["gru", "gru_prev"])
lmm_lstm_ET = do_fit("ET", ["lstm","lstm_prev"])

lmm_srn_gru_ET  = do_fit("ET", ["srn", "srn_prev", "gru", "gru_prev"])
lmm_srn_lstm_ET = do_fit("ET", ["srn", "srn_prev", "lstm", "lstm_prev"])
lmm_gru_lstm_ET = do_fit("ET", ["gru", "gru_prev", "lstm", "lstm_prev"])

println("Running EEG analyses")
data_surp = create_data_surp(data_EEG)

lmm_srn_EEG = do_fit("EEG", ["srn"])
lmm_gru_EEG = do_fit("EEG", ["gru"])
lmm_lstm_EEG = do_fit("EEG", ["lstm"])

lmm_srn_gru_EEG = do_fit("EEG", ["srn",  "gru"])
lmm_srn_lstm_EEG = do_fit("EEG", ["srn", "lstm"])
lmm_gru_lstm_EEG = do_fit("EEG", ["gru", "lstm"])

# Do the model comparisons
chi2_SPR, p_SPR = do_model_comparisons(lmm_srn_SPR, lmm_gru_SPR, lmm_lstm_SPR, lmm_srn_gru_SPR, lmm_srn_lstm_SPR, lmm_gru_lstm_SPR, 4)
chi2_ET, p_ET   = do_model_comparisons(lmm_srn_ET,  lmm_gru_ET,  lmm_lstm_ET,  lmm_srn_gru_ET,  lmm_srn_lstm_ET,  lmm_gru_lstm_ET,  4)
chi2_EEG, p_EEG = do_model_comparisons(lmm_srn_EEG, lmm_gru_EEG, lmm_lstm_EEG, lmm_srn_gru_EEG, lmm_srn_lstm_EEG, lmm_gru_lstm_EEG, 2)

#=
lmm_srn_SPR  = fit(LinearMixedModel, @formula(logRT ~ srn + srn_prev + logRT_prev * log_freq * log_freq_prev * nr_char * nr_char_prev * word_pos
                                         + (1|subj_nr) + (0+srn|subj_nr) + (0+srn_prev|subj_nr) + (0+logRT_prev|subj_nr) + (0+log_freq|subj_nr) + (0+log_freq_prev|subj_nr)
                                         + (0+nr_char|subj_nr) + (0+nr_char_prev|subj_nr) + (0+word_pos|subj_nr) + (1|item)), data_surp)
lmm_gru_SPR  = fit(LinearMixedModel, @formula(logRT ~ gru + gru_prev + logRT_prev * log_freq * log_freq_prev * nr_char * nr_char_prev * word_pos
                                         + (1|subj_nr) + (0+gru|subj_nr) + (0+gru_prev|subj_nr) + (0+logRT_prev|subj_nr) + (0+log_freq|subj_nr) + (0+log_freq_prev|subj_nr)
                                         + (0+nr_char|subj_nr) + (0+nr_char_prev|subj_nr) + (0+word_pos|subj_nr) + (1|item)), data_surp)
lmm_lstm_SPR = fit(LinearMixedModel, @formula(logRT ~ lstm + lstm_prev + logRT_prev * log_freq * log_freq_prev * nr_char * nr_char_prev * word_pos
                                         + (1|subj_nr) + (0+lstm|subj_nr) + (0+lstm_prev|subj_nr) + (0+logRT_prev|subj_nr) + (0+log_freq|subj_nr) + (0+log_freq_prev|subj_nr)
                                         + (0+nr_char|subj_nr) + (0+nr_char_prev|subj_nr) + (0+word_pos|subj_nr) + (1|item)), data_surp)

lmm_srn_gru_SPR  = fit(LinearMixedModel, @formula(logRT ~ srn + srn_prev + gru  + gru_prev  + logRT_prev * log_freq * log_freq_prev * nr_char * nr_char_prev * word_pos
                                             + (1|subj_nr) + (0+srn|subj_nr) + (0+srn_prev|subj_nr) + (0+gru|subj_nr)  + (0+gru_prev|subj_nr)  + (0+logRT_prev|subj_nr) + (0+log_freq|subj_nr) + (0+log_freq_prev|subj_nr)
                                             + (0+nr_char|subj_nr) + (0+nr_char_prev|subj_nr) + (0+word_pos|subj_nr) + (1|item)), data_surp)
lmm_srn_lstm_SPR = fit(LinearMixedModel, @formula(logRT ~ srn + srn_prev + lstm + lstm_prev + logRT_prev * log_freq * log_freq_prev * nr_char * nr_char_prev * word_pos
                                             + (1|subj_nr) + (0+srn|subj_nr) + (0+srn_prev|subj_nr) + (0+lstm|subj_nr) + (0+lstm_prev|subj_nr) + (0+logRT_prev|subj_nr) + (0+log_freq|subj_nr) + (0+log_freq_prev|subj_nr)
                                             + (0+nr_char|subj_nr) + (0+nr_char_prev|subj_nr) + (0+word_pos|subj_nr) + (1|item)), data_surp)
lmm_gru_lstm_SPR = fit(LinearMixedModel, @formula(logRT ~ gru + gru_prev + lstm + lstm_prev + logRT_prev * log_freq * log_freq_prev * nr_char * nr_char_prev * word_pos
                                             + (1|subj_nr) + (0+gru|subj_nr) + (0+gru_prev|subj_nr) + (0+lstm|subj_nr) + (0+lstm_prev|subj_nr) + (0+logRT_prev|subj_nr) + (0+log_freq|subj_nr) + (0+log_freq_prev|subj_nr)
                                             + (0+nr_char|subj_nr) + (0+nr_char_prev|subj_nr) + (0+word_pos|subj_nr) + (1|item)), data_surp)

#-------------------------------------------------------------------------------------------------------------------
println("Running ET analyses")
data_surp = create_data_surp(data_ET)

lmm_srn_ET  = fit(LinearMixedModel, @formula(logRTfirstpass ~ srn + srn_prev + log_freq * log_freq_prev * nr_char * nr_char_prev * word_pos * prevfix
                                          + (1|subj_nr) + (0+srn|subj_nr) + (0+srn_prev|subj_nr) + (0+log_freq|subj_nr) + (0+log_freq_prev|subj_nr)
                                          + (0+nr_char|subj_nr) + (0+nr_char_prev|subj_nr) + (0+word_pos|subj_nr) +(0+prevfix|subj_nr) + (1|item) + (0+prevfix|item)), data_surp)
lmm_gru_ET  = fit(LinearMixedModel, @formula(logRTfirstpass ~ gru + gru_prev + log_freq * log_freq_prev * nr_char * nr_char_prev * word_pos * prevfix
                                          + (1|subj_nr) + (0+gru|subj_nr) + (0+gru_prev|subj_nr) + (0+log_freq|subj_nr) + (0+log_freq_prev|subj_nr)
                                          + (0+nr_char|subj_nr) + (0+nr_char_prev|subj_nr) + (0+word_pos|subj_nr) +(0+prevfix|subj_nr) + (1|item) + (0+prevfix|item)), data_surp)
lmm_lstm_ET = fit(LinearMixedModel, @formula(logRTfirstpass ~ lstm + lstm_prev + log_freq * log_freq_prev * nr_char * nr_char_prev * word_pos * prevfix
                                          + (1|subj_nr) + (0+lstm|subj_nr) + (0+lstm_prev|subj_nr) + (0+log_freq|subj_nr) + (0+log_freq_prev|subj_nr)
                                          + (0+nr_char|subj_nr) + (0+nr_char_prev|subj_nr) + (0+word_pos|subj_nr) +(0+prevfix|subj_nr) + (1|item) + (0+prevfix|item)), data_surp)

lmm_srn_gru_ET  = fit(LinearMixedModel, @formula(logRTfirstpass ~ srn + srn_prev + gru + gru_prev + log_freq * log_freq_prev * nr_char * nr_char_prev * word_pos * prevfix
                                              + (1|subj_nr) + (0+srn|subj_nr) + (0+srn_prev|subj_nr) + (0+gru|subj_nr) + (0+gru_prev|subj_nr) + (0+log_freq|subj_nr) + (0+log_freq_prev|subj_nr)
                                              + (0+nr_char|subj_nr) + (0+nr_char_prev|subj_nr) + (0+word_pos|subj_nr) +(0+prevfix|subj_nr) + (1|item) + (0+prevfix|item)), data_surp)
lmm_srn_lstm_ET = fit(LinearMixedModel, @formula(logRTfirstpass ~ srn + srn_prev + lstm + lstm_prev + log_freq * log_freq_prev * nr_char * nr_char_prev * word_pos * prevfix
                                              + (1|subj_nr) + (0+srn|subj_nr) + (0+srn_prev|subj_nr) + (0+lstm|subj_nr) + (0+lstm_prev|subj_nr) + (0+log_freq|subj_nr) + (0+log_freq_prev|subj_nr)
                                              + (0+nr_char|subj_nr) + (0+nr_char_prev|subj_nr) + (0+word_pos|subj_nr) +(0+prevfix|subj_nr) + (1|item) + (0+prevfix|item)), data_surp)
lmm_gru_lstm_ET = fit(LinearMixedModel, @formula(logRTfirstpass ~ gru + gru_prev + lstm + lstm_prev + log_freq * log_freq_prev * nr_char * nr_char_prev * word_pos * prevfix
                                              + (1|subj_nr) + (0+gru|subj_nr) + (0+gru_prev|subj_nr) + (0+lstm|subj_nr) + (0+lstm_prev|subj_nr) + (0+log_freq|subj_nr) + (0+log_freq_prev|subj_nr)
                                              + (0+nr_char|subj_nr) + (0+nr_char_prev|subj_nr) + (0+word_pos|subj_nr) +(0+prevfix|subj_nr) + (1|item) + (0+prevfix|item)), data_surp)

#-------------------------------------------------------------------------------------------------------------------
println("Running EEG analyses")
data_surp = create_data_surp(data_EEG)

lmm_srn_EEG  = fit(LinearMixedModel, @formula(N400 ~ srn + baseline + log_freq * nr_char * word_pos
                                       + (1|subj_nr) + (0+srn|subj_nr) + (0+baseline|subj_nr) + (0+log_freq|subj_nr) + (0+nr_char|subj_nr) + (0+word_pos|subj_nr) + (1|item)), data_surp)
lmm_gru_EEG  = fit(LinearMixedModel, @formula(N400 ~ gru + baseline + log_freq * nr_char * word_pos
                                       + (1|subj_nr) + (0+gru|subj_nr) + (0+baseline|subj_nr) + (0+log_freq|subj_nr) + (0+nr_char|subj_nr) + (0+word_pos|subj_nr) + (1|item)), data_surp)
lmm_lstm_EEG = fit(LinearMixedModel, @formula(N400 ~ lstm + baseline + log_freq * nr_char * word_pos
                                       + (1|subj_nr) + (0+lstm|subj_nr) + (0+baseline|subj_nr) + (0+log_freq|subj_nr) + (0+nr_char|subj_nr) + (0+word_pos|subj_nr) + (1|item)), data_surp)

lmm_srn_gru_EEG  = fit(LinearMixedModel, @formula(N400 ~ srn + gru + baseline + log_freq * nr_char * word_pos
                                       + (1|subj_nr) + (0+srn|subj_nr) + (0+gru|subj_nr) + (0+baseline|subj_nr) + (0+log_freq|subj_nr) + (0+nr_char|subj_nr) + (0+word_pos|subj_nr)
                                       + (1|item)), data_surp)
lmm_srn_lstm_EEG = fit(LinearMixedModel, @formula(N400 ~ srn + lstm + baseline + log_freq * nr_char * word_pos
                                       + (1|subj_nr) + (0+srn|subj_nr) + (0+lstm|subj_nr) + (0+baseline|subj_nr) + (0+log_freq|subj_nr) + (0+nr_char|subj_nr) + (0+word_pos|subj_nr)
                                       + (1|item)), data_surp)
lmm_gru_lstm_EEG = fit(LinearMixedModel, @formula(N400 ~ gru + lstm + baseline + log_freq * nr_char * word_pos
                                       + (1|subj_nr) + (0+gru|subj_nr) + (0+lstm|subj_nr) + (0+baseline|subj_nr) + (0+log_freq|subj_nr) + (0+nr_char|subj_nr) + (0+word_pos|subj_nr)
                                       + (1|item)), data_surp)
=#
#-----------------------------------------

#---------------------------------------------------
# Perform FDR correction at three different α levels
#--------------------------------------------------

p_vals = [p_SPR[:]; p_ET[:]; p_EEG[:]]   # concatenate all p-values column wise
p_vals = p_vals[.~ismissing.(p_vals)]    # remove missing values

sign001 = benjamini(p_vals,.001)                             # indices of p_vals with corrected p<.001
sign01  = setdiff(benjamini(p_vals,.01), sign001)            # indices of p_vals with corrected .001<p<.01
sign05  = setdiff(benjamini(p_vals,.05), [sign01;sign001])   # indices of p_vals with corrected .01<p<.05

# Split by data type
sign_SPR = (sign001[sign001.<=6],                      sign01[ sign01 .<=6],                      sign05[ sign05 .<=6])
sign_ET  = (sign001[(sign001.>6) .& (sign001.<=12)].-6, sign01[ (sign01 .>6) .& (sign01 .<=12)].-6, sign05[ (sign05 .>6) .& (sign05 .<=12)].-6)
sign_EEG = (sign001[sign001.>12].-12,                   sign01[ sign01 .>12].-12,                   sign05[ sign05 .>12].-12)

println("\nSelf-paced reading")
println("------------------")
print_model_comparisons(chi2_SPR,4,sign_SPR)

println("\nEye tracking")
println("------------------")
print_model_comparisons(chi2_ET, 4, sign_ET)

println("\nEEG")
println("------------------")
print_model_comparisons(chi2_EEG,4,sign_EEG)

max_stars_repo_name: caurnhammer/gated-rnns

max_issues_repo_name: caurnhammer/gated-rnns

max_forks_repo_name: caurnhammer/gated-rnns


text:
correct = sqrt(π/sqrt(ℯ)) / 2
let
    x, w = hermite(10)
    Q = sum(w .* ( x .- 1.0) .* sin.(x))
    @test Q ≈ correct
end


text:
export rule

@rule Categorical(:p, Marginalisation) (q_out::Any, ) = Dirichlet(probvec(q_out) .+ one(eltype(probvec(q_out))))


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{Type{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Axes, F, Args} where Args<:Tuple where F where Axes}, typeof(Base._views), Tuple{Array{Any, 1}}})
    precompile(Tuple{Type{Base.Pair{A, B} where B where A}, Pkg.BinaryPlatforms.FreeBSD, Base.Dict{String, Any}})
    precompile(Tuple{Type{NamedTuple{(:libgfortran_version, :libstdcxx_version, :cxxstring_abi), T} where T<:Tuple}, Tuple{Nothing, Nothing, Symbol}})
end


text:
include("utils.jl")

const tutorialtests = [
  "line"
]

const samplertests = [
  "amm",
  "amwg",
  "bhmc",
  "bia",
  "bmc3",
  "bmg",
  "hmc",
  "mala",
  "nuts",
  "rwm",
  "slice",
  "slicesimplex"
]

const mcmctests = [
  "readcoda"
]

const extensiontests = [
  "newunivardist",
  "newmultivardist"
]

println("Running tests:")

for t in tutorialtests
  @everywhere srand(123)
  @runtest "../doc/tutorial/" t
end

for t in samplertests
  @everywhere srand(123)
  @runtest "../doc/samplers/" t
end

for t in mcmctests
  @runtest "../doc/mcmc/" t
end

for t in extensiontests
  @everywhere srand(123)
  @runtest "../doc/mcmc/" t
end


text:
using Documenter, ReTest

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

deploydocs(
    repo = "github.com/JuliaTesting/ReTest.jl.git",
)


text:
module Chatter

include("message.jl")
include("source.jl")
include("client.jl")
include(joinpath("sources", "IRC.jl"))
include(joinpath("sources", "subreddit.jl"))

export ChatterClient, start, get!, flush!, send, IRC, Subreddit

end


text:
include("../src/read/read_datasets.jl")
include("../src/utils.jl")
include("../src/opq/OPQ.jl");
include("../src/linscan/Linscan.jl");

function demo_opq(
  dataset_name="SIFT1M",
  nread::Integer=Int(1e4)) # Increase this to 1e5 to use the full dataset

  # === Hyperparams ===
  m       = 8
  h       = 256
  verbose = true
  nquery  = Int(1e4)
  knn     = Int(1e3) # Compute recall up to
  b       = Int( log2(h) * m )
  niter   = 10

  # === Train ===
  x_train              = read_dataset(dataset_name, nread )
  C, B, R, train_error = train_opq(x_train, m, h, niter, "natural", verbose)
  @printf("Error in training is %e\n", train_error[end])

  # === Encode the base set ===
  nread_base = Int(1e6)
  x_base     = read_dataset(dataset_name * "_base", nread_base )
  B_base     = quantize_opq( x_base, R, C, verbose )
  base_error = qerror_opq( x_base, B_base, C, R )
  @printf("Error in base is %e\n", base_error)

  # === Compute recall ===
  x_query = read_dataset( dataset_name * "_query", nquery, verbose )
  gt      = read_dataset( dataset_name * "_groundtruth", nquery, verbose )
  if dataset_name == "SIFT1M" || dataset_name == "GIST1M"
    gt = gt + 1;
  end
  gt = convert( Vector{UInt32}, gt[1,1:nquery] )
  B_base  = convert( Matrix{UInt8}, B_base-1 )

  print("Querying m=$m ... ");
  @time dists, idx = linscan_opq( B_base, x_query[:,1:nquery], C, b, R, knn )
  println("done");

  rec = eval_recall( gt, idx, knn );

end # function train_process

# train
demo_opq()

max_issues_repo_name: xiaosu-zhu/local-search-quantization

max_forks_repo_name: xiaosu-zhu/local-search-quantization


text:
struct SinkhornEpsilonScaling{A<:Sinkhorn,T<:Real} <: Sinkhorn
    alg::A
    factor::T
    steps::Int
end

function Base.show(io::IO, alg::SinkhornEpsilonScaling)
    return print(
        io, alg.alg, " with ε-scaling (steps: ", alg.steps, ", factor: ", alg.factor, ")"
    )
end

"""
    SinkhornEpsilonScaling(algorithm::Sinkhorn; factor=1//2, steps=5)

Construct an ε-scaling Sinkhorn algorithm for solving an entropically regularized optimal
transport problem.

The function uses the specified Sinkhorn `algorithm` with `steps` ε-scaling steps
with scaling factor `factor`. It sequentially solves the entropically regularized
optimal transport with regularization parameters
```math
\\varepsilon_i := \\varepsilon \\lambda^{i-k}, \\qquad (i = 1,\\ldots,k),
```
where ``\\lambda`` is the scaling factor and ``k`` the number of scaling steps.
"""
function SinkhornEpsilonScaling(alg::Sinkhorn; factor::Real=1//2, steps::Int=5)
    zero(factor) < factor < one(factor) ||
        throw(ArgumentError("the scaling `factor` must be in (0, 1)"))
    steps > 1 || throw(ArgumentError("number of scaling `steps` must be greater than 1"))
    return SinkhornEpsilonScaling(alg, factor, steps)
end

function sinkhorn(μ, ν, C, ε, alg::SinkhornEpsilonScaling; kwargs...)
    # unpack
    scaling_factor = alg.factor
    scaling_steps = alg.steps

    # initial regularization parameter
    εstep = ε * scaling_factor^(1 - scaling_steps)
    @debug "ε-scaling Sinkhorn algorithm: ε = $εstep"

    # initialize solver and perform Sinkhorn algorithm
    solver = build_solver(μ, ν, C, εstep, alg.alg; kwargs...)
    solve!(solver)

    for _ in 2:(scaling_steps - 1)
        εstep *= scaling_factor
        @debug "ε-scaling Sinkhorn algorithm: ε = $εstep"

        # re-run Sinkhorn algorithm with smaller regularization parameter
        solver = update_epsilon!(solver, εstep)
        solve!(solver)
    end

    # re-run Sinkhorn algorithm with desired regularization parameter
    @debug "ε-scaling Sinkhorn algorithm: ε = $ε"
    solver = update_epsilon!(solver, oftype(εstep, ε))
    solve!(solver)

    # compute final plan
    γ = sinkhorn_plan(solver)

    return γ
end

# support for `SinkhornGibbs` and `SinkhornStabilized`
function update_epsilon!(solver::SinkhornSolver{SinkhornGibbs}, ε::Real)
    cache = solver.cache
    @. cache.K = exp(-solver.C / ε)
    newsolver = SinkhornSolver(
        solver.source,
        solver.target,
        solver.C,
        ε,
        solver.alg,
        solver.atol,
        solver.rtol,
        solver.maxiter,
        solver.check_convergence,
        cache,
        solver.convergence_cache,
    )
    return newsolver
end

function update_epsilon!(solver::SinkhornSolver{<:SinkhornStabilized}, ε::Real)
    newsolver = SinkhornSolver(
        solver.source,
        solver.target,
        solver.C,
        ε,
        solver.alg,
        solver.atol,
        solver.rtol,
        solver.maxiter,
        solver.check_convergence,
        solver.cache,
        solver.convergence_cache,
    )
    absorb!(newsolver)
    return newsolver
end

# deprecations

"""
    sinkhorn_stabilized_epsscaling(
        μ, ν, C, ε;
        scaling_factor=1//2, scaling_steps=5, absorb_tol=1_000, kwargs...
    )

This method is deprecated, please use
```julia
sinkhorn(
    μ,
    ν,
    C,
    ε,
    SinkhornEpsilonScaling(
        SinkhornStabilized(; absorb_tol=absorb_tol);
        factor=scaling_factor,
        steps=scaling_steps,
    );
    kwargs...,
)
```
instead.

See also: [`sinkhorn`](@ref), [`SinkhornEpsilonScaling`](@ref)
"""
function sinkhorn_stabilized_epsscaling(
    μ,
    ν,
    C,
    ε;
    lambda=nothing,
    scaling_factor=lambda,
    k=nothing,
    scaling_steps=k,
    absorb_tol=1_000,
    kwargs...,
)
    if lambda !== nothing
        Base.depwarn(
            "keyword argument `lambda` is deprecated, please use `scaling_factor`",
            :sinkhorn_stabilized_epsscaling,
        )
    end
    if k !== nothing
        Base.depwarn(
            "keyword argument `k` is deprecated, please use `scaling_steps`",
            :sinkhorn_stabilized_epsscaling,
        )
    end

    Base.depwarn(
        "`sinkhorn_stabilized_epsscaling` is deprecated, " *
        "please use `sinkhorn` with `SinkhornEpsilonScaling` instead",
        :sinkhorn_stabilized_epsscaling,
    )

    # construct ε-scaling algorithm
    factor = scaling_factor === nothing ? 1//2 : scaling_factor
    steps = scaling_steps === nothing ? 5 : scaling_steps
    algorithm = SinkhornEpsilonScaling(
        SinkhornStabilized(; absorb_tol=absorb_tol); factor=factor, steps=steps
    )

    return sinkhorn(μ, ν, C, ε, algorithm; kwargs...)
end


text:
function cornerhome!(r)
    x=0
    z=0
    while isborder(r,West)==0
        move!(r,West)
        x=x+1
    end
    while isborder(r,Nord)==0
        move!(r,Nord)
        z=z+1
    end
    Large!(r)
    movetoside!(r,West)
    for t in 1:z
        move!(r,Sud)
    end
    for tt in 1:x
        move!(r,Ost)
    end
end

function movetoside!(r,side)
    while isborder(r,side)==0
        move!(r,side)
    end
end

function Large!(r)
    s=0
    side=Ost
    movetoside!(r,side)
    side=Sud
    movetoside!(r,side)
    while isborder(r,West)==0
        putmarker!(r)
        move!(r,West)
        putmarker!(r)
        s=s+1
    end
    move!(r,Nord)
    s=s-1
    while isborder(r,Nord)==0
        for i in 1:s
            putmarker!(r)
            move!(r,Ost)
            putmarker!(r)
        end
        s=s-1
        side=West
        movetoside!(r,side)
        move!(r,Nord)
    end
    for i in 1:s
        putmarker!(r)
        move!(r,Ost)
        putmarker!(r)
    end
end

cornerhome!(r)


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct AttrValue_ListValue <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function AttrValue_ListValue(; kwargs...)
        obj = new(meta(AttrValue_ListValue), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct AttrValue_ListValue (has cyclic type dependency)
const __meta_AttrValue_ListValue = Ref{ProtoMeta}()
function meta(::Type{AttrValue_ListValue})
    ProtoBuf.metalock() do
        if !isassigned(__meta_AttrValue_ListValue)
            __meta_AttrValue_ListValue[] = target = ProtoMeta(AttrValue_ListValue)
            fnum = Int[2,3,4,5,6,7,8,9]
            pack = Symbol[:i,:f,:b,:_type]
            allflds = Pair{Symbol,Union{Type,String}}[:s => Base.Vector{Vector{UInt8}}, :i => Base.Vector{Int64}, :f => Base.Vector{Float32}, :b => Base.Vector{Bool}, :_type => Base.Vector{Int32}, :shape => Base.Vector{TensorShapeProto}, :tensor => Base.Vector{TensorProto}, :func => "Base.Vector{NameAttrList}"]
            meta(target, AttrValue_ListValue, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_AttrValue_ListValue[]
    end
end
function Base.getproperty(obj::AttrValue_ListValue, name::Symbol)
    if name === :s
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Vector{UInt8}}
    elseif name === :i
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int64}
    elseif name === :f
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float32}
    elseif name === :b
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Bool}
    elseif name === :_type
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int32}
    elseif name === :shape
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{TensorShapeProto}
    elseif name === :tensor
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{TensorProto}
    elseif name === :func
        return (obj.__protobuf_jl_internal_values[name])::Any
    else
        getfield(obj, name)
    end
end

mutable struct AttrValue <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function AttrValue(; kwargs...)
        obj = new(meta(AttrValue), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct AttrValue (has cyclic type dependency)
const __meta_AttrValue = Ref{ProtoMeta}()
function meta(::Type{AttrValue})
    ProtoBuf.metalock() do
        if !isassigned(__meta_AttrValue)
            __meta_AttrValue[] = target = ProtoMeta(AttrValue)
            fnum = Int[2,3,4,5,6,7,8,1,10,9]
            allflds = Pair{Symbol,Union{Type,String}}[:s => Vector{UInt8}, :i => Int64, :f => Float32, :b => Bool, :_type => Int32, :shape => TensorShapeProto, :tensor => TensorProto, :list => AttrValue_ListValue, :func => "NameAttrList", :placeholder => AbstractString]
            oneofs = Int[1,1,1,1,1,1,1,1,1,1]
            oneof_names = Symbol[Symbol("value")]
            meta(target, AttrValue, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_AttrValue[]
    end
end
function Base.getproperty(obj::AttrValue, name::Symbol)
    if name === :s
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :i
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :f
        return (obj.__protobuf_jl_internal_values[name])::Float32
    elseif name === :b
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :_type
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :shape
        return (obj.__protobuf_jl_internal_values[name])::TensorShapeProto
    elseif name === :tensor
        return (obj.__protobuf_jl_internal_values[name])::TensorProto
    elseif name === :list
        return (obj.__protobuf_jl_internal_values[name])::AttrValue_ListValue
    elseif name === :func
        return (obj.__protobuf_jl_internal_values[name])::Any
    elseif name === :placeholder
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct NameAttrList_AttrEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function NameAttrList_AttrEntry(; kwargs...)
        obj = new(meta(NameAttrList_AttrEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct NameAttrList_AttrEntry (mapentry) (has cyclic type dependency)
const __meta_NameAttrList_AttrEntry = Ref{ProtoMeta}()
function meta(::Type{NameAttrList_AttrEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_NameAttrList_AttrEntry)
            __meta_NameAttrList_AttrEntry[] = target = ProtoMeta(NameAttrList_AttrEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => AttrValue]
            meta(target, NameAttrList_AttrEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_NameAttrList_AttrEntry[]
    end
end
function Base.getproperty(obj::NameAttrList_AttrEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::AttrValue
    else
        getfield(obj, name)
    end
end

mutable struct NameAttrList <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function NameAttrList(; kwargs...)
        obj = new(meta(NameAttrList), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct NameAttrList (has cyclic type dependency)
const __meta_NameAttrList = Ref{ProtoMeta}()
function meta(::Type{NameAttrList})
    ProtoBuf.metalock() do
        if !isassigned(__meta_NameAttrList)
            __meta_NameAttrList[] = target = ProtoMeta(NameAttrList)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :attr => "Base.Dict{AbstractString,AttrValue}"]
            meta(target, NameAttrList, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_NameAttrList[]
    end
end
function Base.getproperty(obj::NameAttrList, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :attr
        return (obj.__protobuf_jl_internal_values[name])::Any
    else
        getfield(obj, name)
    end
end

export AttrValue_ListValue, AttrValue, NameAttrList_AttrEntry, NameAttrList, AttrValue_ListValue, AttrValue, NameAttrList_AttrEntry, NameAttrList
# mapentries: "NameAttrList_AttrEntry" => ("AbstractString", "AttrValue")


text:
# The primes 3, 7, 109, and 673, are quite remarkable. By taking any two primes
# and concatenating them in any order the result will always be prime. For
# example, taking 7 and 109, both 7109 and 1097 are prime. The sum of these four
# primes, 792, represents the lowest sum for a set of four primes with this
# property.
#
# Find the lowest sum for a set of five primes for which any two primes
# concatenate to produce another prime.

using ProjectEulerSolutions

# Fairly brute-force approach, start with Seive of Eratosthenes to get an
# initial set of primes, then start building up a combination of 5 that each
# satisfy the concatenated condition given all of the previous numbers.
# Short-cutting can be aided by upper bounding some prime number ranges based
# on previous best results.  Finally, the order you pass through the pairs
# matters, looping with the largest value on the outside loop finds the minimum
# set faster.  A recursive solution would look cleaner, but it wouldn't
# necessarily be faster.
#
# This seems like the kind of thing you can do by building a graph with edges
# representing legal concatenated pairs and then solving for the fully connected
# cliques of size 5.
function p060solution(n::Integer=1_000)::Integer
    primes = sieve_eratosthenes(n, 3)

    min_sum = 10_000_000
    n = length(primes)

    # Use a bitarray, but might need to like a Dict of Sets if the number of
    # primes is very large
    pairs = falses(n, n)
    ashift = typeof(primes[1])(10)

    # Find concatenated prime pairs
    for a in 1:n
        # Break early if the sum is too large
        if 5 * primes[a] > min_sum
            break
        end

        # Fill in the pairs bitarray
        prime_a = primes[a]
        jshift = typeof(prime_a)(10)
        if prime_a > ashift
            ashift *= 10
        end
        for j in a+1:n
            prime_j = primes[j]
            if prime_j > jshift
                jshift *= 10
            end
            if isprime(ashift * prime_j + prime_a) && isprime(jshift * prime_a + prime_j)
                pairs[a,j] = true
            end
        end

        # Start looking for legal combinations of 2, 3, 4, and 5 primes
        for b in 4:a-1
            if !pairs[b,a]
                continue
            end

            for c in 3:b-1
                if !pairs[c,a] || !pairs[c,b]
                    continue
                end

                for d in 2:c-1
                    if !pairs[d,a] || !pairs[d,b] || !pairs[d,c]
                        continue
                    end

                    for e in 1:d-1
                        if pairs[e,a] && pairs[e,b] && pairs[e,c] && pairs[e,d]
                            min_sum = primes[a] + primes[b] + primes[c] + primes[d] + primes[e]
                        end
                    end
                end
            end
        end
    end

    return min_sum
end

p060 = Problems.Problem(p060solution)

Problems.benchmark(p060, 10_000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
"""
Handles Http server related functionality, manages requests and responses and their logging.
"""
module AppServer

using Revise, HTTP, HTTP.IOExtras, HTTP.Sockets, Millboard, MbedTLS, WebSockets, URIParser, Sockets, Distributed
using Genie, Genie.Router, Genie.Loggers, Genie.Sessions, Genie.Configuration, Genie.WebChannels

"""
    startup(port::Int = 8000)

Starts the web server on the configurated port.
Automatically invoked when Genie is started with the `s` or the `server:start` command line params.

# Examples
```julia
julia> AppServer.startup()
Listening on 0.0.0.0:8000...
```
"""
function startup(port::Int = 8000, host = "10.52.247.40"; ws_port = port + 1, async = !Genie.config.run_as_server)
  web_server = HTTP.Servers.Server((req) -> begin
    setup_http_handler(req, req.response)
  end, devnull)

  if async
    @async HTTP.Servers.serve(web_server, host, port)
  else
    HTTP.Servers.serve(web_server, host, port)
  end

  log("Web server running at $host:$port")

  if Genie.config.websocket_server
    @async HTTP.listen(host, ws_port) do req
      if HTTP.WebSockets.is_upgrade(req.message)
        HTTP.WebSockets.upgrade(req) do ws
          setup_ws_handler(req.message, ws)
        end
      end
      log("WebSockets server running at $host:$ws_port")
    end
  end
end

"""
"""
function setup_http_handler(req, res)
  try
    response = @fetch handle_request(req, res)
    if response.status >= 500 && response.status < 600
      error(response)
    end

    response
  catch ex
    log(string(ex), :critical)
    log(sprint(io->Base.show_backtrace(io, catch_backtrace() )), :critical)
    log("$(@__FILE__):$(@__LINE__)", :critical)

    message = Genie.Configuration.is_prod() ?
                "The error has been logged and we'll look into it ASAP." :
                string(ex, " in $(@__FILE__):$(@__LINE__)", "\n\n", sprint(io->Base.show_backtrace(io, catch_backtrace())))

    Genie.Router.serve_error_file(500, message, Genie.Router.@params)
  end
end

"""
"""
function setup_ws_handler(req, ws_client)
  while ! eof(ws_client)
    write(ws_client, String(@fetch handle_ws_request(req, String(readavailable(ws_client)), ws_client)))
  end
end

"""
    handle_request(req::Request, res::Response, ip::IPv4 = ip"0.0.0.0") :: Response

Http server handler function - invoked when the server gets a request.
"""
function handle_request(req::HTTP.Request, res::HTTP.Response, ip::IPv4 = ip"0.0.0.0") :: HTTP.Response
  Genie.config.server_signature != "" && sign_response!(res)
  set_headers!(req, res, Genie.Router.route_request(req, res, ip))
end

function set_headers!(req::HTTP.Request, res::HTTP.Response, app_response::HTTP.Response) :: HTTP.Response
  if req.method == Genie.Router.OPTIONS
    Genie.config.cors_headers["Access-Control-Allow-Origin"] = strip(Genie.config.cors_headers["Access-Control-Allow-Origin"])

    ! isempty(Genie.config.cors_allowed_origins) &&
      in(req.headers["Origin"], Genie.config.cors_allowed_origins) &&
      (Genie.config.cors_headers["Access-Control-Allow-Origin"] == "" ||
        Genie.config.cors_headers["Access-Control-Allow-Origin"] == "*") &&
      (Genie.config.cors_headers["Access-Control-Allow-Origin"] = req.headers["Origin"])

    app_response.headers = merge(res.headers, Genie.config.cors_headers)
  end

  app_response.headers = [d for d in merge(Dict(res.headers), Dict(app_response.headers))]

  # app_response.cookies = merge(res.cookies, app_response.cookies)

  app_response
end

"""
    handle_ws_request(req::Request, client::Client, ip::IPv4 = ip"0.0.0.0") :: String

Http server handler function - invoked when the server gets a request.
"""
function handle_ws_request(req, msg::String, ws_client, ip::IPv4 = ip"0.0.0.0") :: String
  msg == "" && return "" # keep alive
  Genie.Router.route_ws_request(req, msg, ws_client, ip)
end

"""
    sign_response!(res::Response) :: Response

Adds a signature header to the response using the value in `Genie.config.server_signature`.
If `Genie.config.server_signature` is empty, the header is not added.
"""
function sign_response!(res::HTTP.Response) :: HTTP.Response
  headers = Dict(res.headers)
  isempty(Genie.config.server_signature) || (headers["Server"] = Genie.config.server_signature)

  res.headers = [k for k in headers]
  res
end

end

max_stars_repo_name: kiranchowdary/Genie

max_issues_repo_name: kiranchowdary/Genie

max_forks_repo_name: kiranchowdary/Genie


text:
reciprocal_distretization(x, L) = reciprocal_distretization(lattice(x), L)
function reciprocal_distretization(lattice::AbstractLattice, L)
    r1, r2 = reciprocal_vectors(lattice, L)
    reciprocal_distretization(r1, r2, L)
end
function reciprocal_distretization(v1, v2, L)
    pgram = map(((i, j) for i in 0:L-1 for j in 0:L-1)) do (i, j)
        out = [0., 0.]
        temp = [0., 0.]
        min = Inf
        for n in -1:1, m in -1:1
            temp .= v1 * (n + i/L) .+ v2 * (m + j/L)
            if norm(temp) < min
                min = norm(temp)
                out .= temp
            end
        end
        out
    end
    return pgram
end

structure_factor(args...) = fourier_transform(args...)
function fourier_transform(mc::DQMC, key::Symbol)
    fourier_transform(mc, mc[key])
end
function fourier_transform(mc::DQMC, m::DQMCMeasurement)
    qs = reciprocal_distretization(mc)
    fourier_transform(mc, qs, m)
end
function fourier_transform(mc::DQMC, qs::Vector, m::DQMCMeasurement)
    dirs = directions(mc)
    values = mean(m)
    fourier_transform(qs, dirs, values)
end

function fourier_transform(qs::Vector, dirs::Vector, value::Vector)
    map(qs) do q
        sum(cis(dot(q, v)) * o for (v, o) in zip(dirs, values))
    end
end

uniform_fourier(mc::DQMC, key::Symbol) = uniform_fourier(mc[key])
uniform_fourier(m::DQMCMeasurement) = sum(mean(m))

apply_symmetry(mc::DQMC, key::Symbol, weights=(1)) = apply_symmetry(mc[key], weights)
apply_symmetry(m::DQMCMeasurement, weights=(1)) = apply_symmetry(mean(m), weights)
function apply_symmetry(data::AbstractArray{T, 3}, weights=(1)) where {T}
    out = zeros(T, size(data, 1))
    @avx for i in eachindex(out)
        for j in eachindex(weights), k in eachindex(weights)
            out[i] += weights[j] * weights[k] * data[i,j,k]
        end
    end
    out
end
function apply_symmetry(data::AbstractArray{T, 1}, weights=(1)) where {T}
    out = zero(T)
    @avx for i in eachindex(weights)
        out += weights[i] * data[i]
    end
    out
end
function apply_symmetry(data::AbstractArray{T, 2}, weights=(1)) where {T}
    out = zeros(T, size(data, 1))
    @avx for i in eachindex(out)
        for j in eachindex(weights)
            out[i] += weights[j] * data[i,j]
        end
    end
    out
end

# cc2superfluid density
superfluid_density(mc, key::Symbol) = superfluid_density(mc, mc[key])
function superfluid_density(mc, m::DQMCMeasurement)
    superfluid_density(mean(m), lattice(mc))
end
function superfluid_density(data::Array{T, 2}, lattice)
    qx, qy = reciprocal_vectors(lattice)
    dirs = directions(mc)
    output = ComplexF64(0)
    for i in axes(data, 1)
        for j in axes(data, 2)
            output += (cis(dot(qy, dirs[j])) - cis(dot(qx, dirs[j]))) * data[i, j]
        end
    end
    output
end


text:
"""
Example: a univariate time series (local linear trend) with an increasing variance, going from 1 to 5.
"""

# Packages
using Plots

T = 1_000
y = zeros(T, 1)
F = [1, 0]
G = Matrix(I, 2, 2)
m = zeros(2, 1)
P = Matrix(I, 2, 2) * 1000
S = Matrix(I, 1, 1)

#s = StateSpaceModel(y, F, G, m, P, S, 0.95, 0.95)
s = LocalLevelTrend(y, m, P, S, 0.98, 0.98)

# Simulate
simulate!(s)

# Estimate
m = estimate(s)

plot([s.y m.predicted.μ m.filtered.μ], labels = ["Simulated" "Predicted" "Filtered"], legend = :topleft)

#plot([y m.predicted.μ m.filtered.μ], labels = ["True" "Predicted" "Filtered"], legend = :topleft)
#plot([fill(Σ, T) m.predicted.Σ[:, 1, 1] m.filtered.Σ[:, 1, 1]][2:end, :], labels = ["True" "Predicted" "Filtered"], legend = :topleft)


text:
using PyPlot
using DelimitedFiles

if !isdir("figures_summary")
  mkdir("figures_summary")
end

m = 15
n = 30
h = 30.0
dz = 3.0 # meters
dx = 3.0
nz = Int64(round((m * h) / dz)) + 1
nx = Int64(round((n * h) / dx)) + 1
z_src = (collect(5:10:nz-5) .- 1 ) .* dz .+ dz/2.0
x_src = (5-1)ones(Int64, size(z_src)) .* dx .+ dx/2.0
z_rec = (collect(5:1:nz-5) .- 1) .* dz .+ dz/2.0
x_rec = (nx-5-1) .* ones(Int64, size(z_rec)) .*dx .+ dx/2.0

z_inj = (9-1)*h + h/2.0
x_inj = (3-1)*h + h/2.0
z_prod = (9-1)*h + h/2.0
x_prod = (28-1)*h + h/2.0

iter = 100
Prj_names = ["CO2", "CO2_1src", "CO2_2surveys", "Brie_3_nocoefupdate", "Brie_tune_coef_true3_start2", "Brie_true3_set2_noupdate", "Brie_tune_coef_true3_start2_scale30"]
K_name = "/K$iter.txt"

rc("axes", titlesize=20)
rc("axes", labelsize=18)
rc("xtick", labelsize=18)
rc("ytick", labelsize=18)
rc("legend", fontsize=20)
# true model
figure()
K = 20.0 .* ones(m,n)
K[8:10,:] .= 120.0
imshow(K, extent=[0,n*h,m*h,0]);
xlabel("Distance (m)")
ylabel("Depth (m)")
cb = colorbar()
clim([20, 120])
cb.set_label("Permeability (md)")
shot_inds = collect(1:length(z_src))
scatter(x_src[shot_inds], z_src[shot_inds], c="w", marker="*")
scatter(x_rec, z_rec, s=16.0, c="r", marker="v")
scatter(x_inj, z_inj, c="r", marker=">", s=64)
scatter(x_prod, z_prod, c="r", marker="<", s=64)
savefig("figures_summary/K_true.pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);

# init model
figure()
K = 20.0 .* ones(m,n)
imshow(K, extent=[0,n*h,m*h,0]);
xlabel("Distance (m)")
ylabel("Depth (m)")
cb = colorbar()
clim([20, 120])
cb.set_label("Permeability (md)")
shot_inds = collect(1:length(z_src))
scatter(x_src[shot_inds], z_src[shot_inds], c="w", marker="*")
scatter(x_rec, z_rec, s=16.0, c="r", marker="v")
scatter(x_inj, z_inj, c="r", marker=">", s=64)
scatter(x_prod, z_prod, c="r", marker="<", s=64)
savefig("figures_summary/K_init.pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);

figure()
iPrj = 1
K = readdlm(Prj_names[iPrj] * K_name)
imshow(K, extent=[0,n*h,m*h,0]);
xlabel("Distance (m)")
ylabel("Depth (m)")
cb = colorbar()
clim([20, 120])
cb.set_label("Permeability (md)")
shot_inds = collect(1:length(z_src))
scatter(x_src[shot_inds], z_src[shot_inds], c="w", marker="*")
scatter(x_rec, z_rec, s=16.0, c="r", marker="v")
scatter(x_inj, z_inj, c="r", marker=">", s=64)
scatter(x_prod, z_prod, c="r", marker="<", s=64)
savefig("figures_summary/K_$(Prj_names[iPrj]).pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);

figure()
iPrj = 3
K = readdlm(Prj_names[iPrj] * K_name)
imshow(K, extent=[0,n*h,m*h,0]);
xlabel("Distance (m)")
ylabel("Depth (m)")
cb = colorbar()
clim([20, 120])
cb.set_label("Permeability (md)")
shot_inds = collect(1:length(z_src))
scatter(x_src[shot_inds], z_src[shot_inds], c="w", marker="*")
scatter(x_rec, z_rec, s=16.0, c="r", marker="v")
scatter(x_inj, z_inj, c="r", marker=">", s=64)
scatter(x_prod, z_prod, c="r", marker="<", s=64)
savefig("figures_summary/K_$(Prj_names[iPrj]).pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);

figure()
iPrj = 4
K = readdlm(Prj_names[iPrj] * K_name)
imshow(K, extent=[0,n*h,m*h,0]);
xlabel("Distance (m)")
ylabel("Depth (m)")
cb = colorbar()
clim([20, 120])
cb.set_label("Permeability (md)")
shot_inds = collect(1:length(z_src))
scatter(x_src[shot_inds], z_src[shot_inds], c="w", marker="*")
scatter(x_rec, z_rec, s=16.0, c="r", marker="v")
scatter(x_inj, z_inj, c="r", marker=">", s=64)
scatter(x_prod, z_prod, c="r", marker="<", s=64)
savefig("figures_summary/K_$(Prj_names[iPrj]).pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);

figure()
iPrj = 5
K = readdlm(Prj_names[iPrj] * K_name)
imshow(K, extent=[0,n*h,m*h,0]);
xlabel("Distance (m)")
ylabel("Depth (m)")
cb = colorbar()
clim([20, 120])
cb.set_label("Permeability (md)")
shot_inds = collect(1:length(z_src))
scatter(x_src[shot_inds], z_src[shot_inds], c="w", marker="*")
scatter(x_rec, z_rec, s=16.0, c="r", marker="v")
scatter(x_inj, z_inj, c="r", marker=">", s=64)
scatter(x_prod, z_prod, c="r", marker="<", s=64)
savefig("figures_summary/K_$(Prj_names[iPrj]).pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);

figure()
iPrj = 6
K = readdlm(Prj_names[iPrj] * K_name)
imshow(K, extent=[0,n*h,m*h,0]);
xlabel("Distance (m)")
ylabel("Depth (m)")
cb = colorbar()
clim([20, 120])
cb.set_label("Permeability (md)")
shot_inds = collect(1:length(z_src))
scatter(x_src[shot_inds], z_src[shot_inds], c="w", marker="*")
scatter(x_rec, z_rec, s=16.0, c="r", marker="v")
scatter(x_inj, z_inj, c="r", marker=">", s=64)
scatter(x_prod, z_prod, c="r", marker="<", s=64)
savefig("figures_summary/K_$(Prj_names[iPrj]).pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);

figure()
iPrj = 7
K = readdlm(Prj_names[iPrj] * K_name)
imshow(K, extent=[0,n*h,m*h,0]);
xlabel("Distance (m)")
ylabel("Depth (m)")
cb = colorbar()
clim([20, 120])
cb.set_label("Permeability (md)")
shot_inds = collect(1:length(z_src))
scatter(x_src[shot_inds], z_src[shot_inds], c="w", marker="*")
scatter(x_rec, z_rec, s=16.0, c="r", marker="v")
scatter(x_inj, z_inj, c="r", marker=">", s=64)
scatter(x_prod, z_prod, c="r", marker="<", s=64)
savefig("figures_summary/K_$(Prj_names[iPrj]).pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);

figure()
z_src = [75].* dz .+ dz/2.0 # single source
x_src = (5-1)ones(Int64, size(z_src)) .* dx .+ dx/2.0
iPrj = 2
K = readdlm(Prj_names[iPrj] * K_name)
imshow(K, extent=[0,n*h,m*h,0]);
xlabel("Distance (m)")
ylabel("Depth (m)")
cb = colorbar()
clim([20, 120])
cb.set_label("Permeability (md)")
shot_inds = collect(1:length(z_src))
scatter(x_src[shot_inds], z_src[shot_inds], c="w", marker="*")
scatter(x_rec, z_rec, s=16.0, c="r", marker="v")
scatter(x_inj, z_inj, c="r", marker=">", s=64)
scatter(x_prod, z_prod, c="r", marker="<", s=64)
savefig("figures_summary/K_$(Prj_names[iPrj]).pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);

rc("axes", titlesize=14)
rc("axes", labelsize=14)
rc("xtick", labelsize=14)
rc("ytick", labelsize=14)
rc("legend", fontsize=14)
figure()
iPrj = 7
brie_coef = readdlm(Prj_names[iPrj] * "/brie_coef.txt")[:,2]./30.0
plot(0:100,[2;brie_coef], "k");grid(ls="--")
# plot(1:100, 3ones(100))
xlabel("Iterations")
ylabel("Brie model coefficient")
savefig("figures_summary/brie_coef_curve.pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);


text:
include("utils.jl")
@testset "Symmetric Positive Definite Matrices" begin
    M1 = Manifolds.SymmetricPositiveDefinite(3)
    M2 = MetricManifold(Manifolds.SymmetricPositiveDefinite(3), Manifolds.LinearAffineMetric())
    M3 = MetricManifold(Manifolds.SymmetricPositiveDefinite(3), Manifolds.LogCholeskyMetric())
    M4 = MetricManifold(Manifolds.SymmetricPositiveDefinite(3), Manifolds.LogEuclideanMetric())

    @test is_default_metric(M2) == Val{true}()
    @test is_default_metric(M1, Manifolds.LinearAffineMetric()) == Val{true}()
    @test is_default_metric(M1, Manifolds.LogCholeskyMetric()) == Val{false}()
    @test is_default_metric(M3) == Val{false}()

    metrics = [M1, M2, M3]
    types = [ Matrix{Float32},
            Matrix{Float64},
            MMatrix{3,3,Float32},
        ]
    for M in metrics
        @testset "$(typeof(M))" begin
            @test representation_size(M) == (3,3)
            for T in types
                if M == M3 && T <: MMatrix
                    # Cholesky or something does not work in vector_transport yet for MMatrix
                    continue
                end
                A(α) = [1. 0. 0.; 0. cos(α) sin(α); 0. -sin(α) cos(α)]
                ptsF = [#
                    [1. 0. 0.; 0. 1. 0.; 0. 0. 1],
                    [2. 0. 0.; 0. 2. 0.; 0. 0. 1],
                    A(π/6) * [1. 0. 0.; 0. 2. 0.; 0. 0. 1] * transpose(A(π/6)),
                ]
                pts = [convert(T, a) for a in ptsF]
                test_manifold(M, pts;
                    test_vector_transport = true,
                    test_forward_diff = false,
                    test_reverse_diff = false,
                    exp_log_atol_multiplier = 8
                )
            end
            @testset "Test Error cases in is_manifold_point and is_tangent_vector" begin
                pt1f = zeros(2,3); # wrong size
                pt2f = [1. 0. 0.; 0. 0. 0.; 0. 0. 1.]; # not positive Definite
                pt3f = [2. 0. 1.; 0. 1. 0.; 0. 0. 4.]; # not symmetric
                pt4 = [2. 1. 0.; 1. 2. 0.; 0. 0. 4.]
                @test !is_manifold_point(M,pt1f)
                @test !is_manifold_point(M,pt2f)
                @test !is_manifold_point(M,pt3f)
                @test is_manifold_point(M, pt4)
                @test !is_tangent_vector(M,pt4, pt1f)
                @test is_tangent_vector(M,pt4, pt2f)
                @test !is_tangent_vector(M,pt4, pt3f)
            end
        end
    end
    x = [1. 0. 0.; 0. 1. 0.; 0. 0. 1]
    y = [2. 0. 0.; 0. 2. 0.; 0. 0. 1]
    @testset "Convert SPD to Cholesky" begin
        v = log(M1,x,y)
        (l,w) = Manifolds.spd_to_cholesky(x,v)
        (xs,vs) = Manifolds.cholesky_to_spd(l,w)
        @test isapprox(xs,x)
        @test isapprox(vs,v)
    end
    @testset "Preliminary tests for LogEuclidean" begin
        @test representation_size(M4) == (3,3)
        @test isapprox( distance(M4,x,y), sqrt(2)log(2))
        @test manifold_dimension(M4) == manifold_dimension(M1)
    end
    @testset "Test for tangent ONB on LinearAffineMetric" begin
        v = log(M2,x,y)
        (X,k) = tangent_orthonormal_basis(base_manifold(M2),x,v)
        @test isapprox(0.0,first(k))
        for i = 1:length(X)
            @test isapprox(1.0, norm(M2,x,X[1]))
            for j=i+1:length(X)
                @test isapprox(0.0, inner(M2,x,X[1],X[2]))
            end
        end
    end
end