text:
log(message) = @info "$(now()) | $message"

const ALPHANUMERIC_SNAKE_CASE_CHARACTERS = Char['_',
                                                '0':'9'...,
                                                'a':'z'...]

function is_lower_snake_case_alphanumeric(x::AbstractString, also_allow=())
    return !isempty(x) && !startswith(x, '_') && !endswith(x, '_') &&
           all(i -> i in ALPHANUMERIC_SNAKE_CASE_CHARACTERS || i in also_allow, x)
end

# TODO port a generic version of this + notion of primary key to Legolas.jl
function _fully_validate_legolas_table(table, schema::Legolas.Schema, primary_key)
    Legolas.validate(table, schema)
    primary_counts = Dict{Any,Int}()
    for (i, row) in enumerate(Tables.rows(table))
        local validated_row
        try
            validated_row = Legolas.Row(schema, row)
        catch err
            log("Encountered invalid row $i when validating table's compliance with $schema:")
            rethrow(err)
        end
        primary = Tables.getcolumn(validated_row, primary_key)
        primary_counts[primary] = get(primary_counts, primary, 0) + 1
    end
    filter!(>(1) ∘ last, primary_counts)
    if !isempty(primary_counts)
        throw(ArgumentError("duplicate $primary_key values found in given $schema table: $primary_counts"))
    end
    return table
end

#####
##### arrrrr i'm a pirate
#####
# The Onda Format defines `span` elements to correspond to the Arrow-equivalent of `(start=Nanosecond(...), stop=Nanosecond(...))`.
# Here we define the generic `TimeSpans` interface on this type in order to ensure that this structure can be treated like a
# `TimeSpan` anywhere. This way, callers don't need to do any fiddling if e.g. they're working with an Onda file written from
# a source that wasn't using `TimeSpans` (e.g. if it was written out by a non-Julia process).

const NamedTupleTimeSpan = NamedTuple{(:start, :stop),Tuple{Nanosecond,Nanosecond}}

TimeSpans.istimespan(::NamedTupleTimeSpan) = true
TimeSpans.start(x::NamedTupleTimeSpan) = x.start
TimeSpans.stop(x::NamedTupleTimeSpan) = x.stop

const TIME_SPAN_ARROW_NAME = Symbol("JuliaLang.TimeSpan")

Arrow.ArrowTypes.arrowname(::Type{TimeSpan}) = TIME_SPAN_ARROW_NAME
ArrowTypes.JuliaType(::Val{TIME_SPAN_ARROW_NAME}) = TimeSpan

#####
##### zstd_compress/zstd_decompress
#####

function zstd_compress(bytes::Vector{UInt8}, level=3)
    compressor = ZstdCompressor(; level=level)
    TranscodingStreams.initialize(compressor)
    compressed_bytes = transcode(compressor, bytes)
    TranscodingStreams.finalize(compressor)
    return compressed_bytes
end

zstd_decompress(bytes::Vector{UInt8}) = transcode(ZstdDecompressor, bytes)

#####
##### read/write/bytes/streams
#####

jump(io::IO, n) = (read(io, n); nothing)
jump(io::IOStream, n) = (skip(io, n); nothing)
jump(io::IOBuffer, n) = ((io.seekable ? skip(io, n) : read(io, n)); nothing)

unsafe_vec_uint8(x::AbstractVector{UInt8}) = convert(Vector{UInt8}, x)
unsafe_vec_uint8(x::Base.ReinterpretArray{UInt8,1}) = unsafe_wrap(Vector{UInt8}, pointer(x), length(x))

"""
    read_byte_range(path, byte_offset, byte_count)

Return the equivalent `read(path)[(byte_offset + 1):(byte_offset + byte_count)]`,
but try to avoid reading unreturned intermediate bytes. Note that the
effectiveness of this method depends on the type of `path`.
"""
function read_byte_range(path, byte_offset, byte_count)
    return open(path, "r") do io
        jump(io, byte_offset)
        return read(io, byte_count)
    end
end

read_byte_range(path, ::Missing, ::Missing) = read(path)

write_full_path(path::AbstractString, bytes) = (mkpath(dirname(path)); write(path, bytes))
write_full_path(path, bytes) = write(path, bytes)


text:
"""Calculate the number of grains on square `square`."""
function on_square(square)
    check_square_input(square)
    2 ^ (square - 1)
end

"""Calculate the total number of grains after square `square`."""
function total_after(square)
    check_square_input(square)
    sum(map(on_square, 1:square))
end

"""Validate an arbitrary square."""
function check_square_input(square)
    square == 0 && throw(DomainError(square, "Square input of zero is invalid."))
    square <  0 && throw(DomainError(square, "Negative square input is invalid."))
    square > 64 && throw(DomainError(square, "Square input greater than 64 is invalid."))
end

if VERSION < v"0.7"  # backwards compatibility
    Base.DomainError(val) = DomainError()
    Base.DomainError(val, msg) = DomainError()
end

ext: jl


text:
# Compressor Codecs
# ==================

abstract type CompressorCodec <: TranscodingStreams.Codec end

function Base.show(io::IO, codec::CompressorCodec)
    print(io, summary(codec), "(level=$(codec.level), windowbits=$(codec.windowbits))")
end

# Gzip
# ----

struct GzipCompressor <: CompressorCodec
    zstream::ZStream
    level::Int
    windowbits::Int
end

"""
    GzipCompressor(;level=$(Z_DEFAULT_COMPRESSION), windowbits=$(Z_DEFAULT_WINDOWBITS))

Create a gzip compression codec.

Arguments
---------
- `level`: compression level (-1..9)
- `windowbits`: size of history buffer (8..15)
"""
function GzipCompressor(;level::Integer=Z_DEFAULT_COMPRESSION,
                         windowbits::Integer=Z_DEFAULT_WINDOWBITS)
    if !(-1 ≤ level ≤ 9)
        throw(ArgumentError("compression level must be within -1..9"))
    elseif !(8 ≤ windowbits ≤ 15)
        throw(ArgumentError("windowbits must be within 8..15"))
    end
    return GzipCompressor(ZStream(), level, windowbits+16)
end

const GzipCompressorStream{S} = TranscodingStream{GzipCompressor,S} where S<:IO

"""
    GzipCompressorStream(stream::IO; kwargs...)

Create a gzip compression stream (see `GzipCompressor` for `kwargs`).
"""
function GzipCompressorStream(stream::IO; kwargs...)
    x, y = splitkwargs(kwargs, (:level, :windowbits))
    return TranscodingStream(GzipCompressor(;x...), stream; y...)
end

# Zlib
# ----

struct ZlibCompressor <: CompressorCodec
    zstream::ZStream
    level::Int
    windowbits::Int
end

"""
    ZlibCompressor(;level=$(Z_DEFAULT_COMPRESSION), windowbits=$(Z_DEFAULT_WINDOWBITS))

Create a zlib compression codec.

Arguments
---------
- `level`: compression level (-1..9)
- `windowbits`: size of history buffer (8..15)
"""
function ZlibCompressor(;level::Integer=Z_DEFAULT_COMPRESSION,
                         windowbits::Integer=Z_DEFAULT_WINDOWBITS)
    if !(-1 ≤ level ≤ 9)
        throw(ArgumentError("compression level must be within -1..9"))
    elseif !(8 ≤ windowbits ≤ 15)
        throw(ArgumentError("windowbits must be within 8..15"))
    end
    return ZlibCompressor(ZStream(), level, windowbits)
end

const ZlibCompressorStream{S} = TranscodingStream{ZlibCompressor,S} where S<:IO

"""
    ZlibCompressorStream(stream::IO)

Create a zlib compression stream (see `ZlibCompressor` for `kwargs`).
"""
function ZlibCompressorStream(stream::IO; kwargs...)
    x, y = splitkwargs(kwargs, (:level, :windowbits))
    return TranscodingStream(ZlibCompressor(;x...), stream; y...)
end

# Deflate
# -------

struct DeflateCompressor <: CompressorCodec
    zstream::ZStream
    level::Int
    windowbits::Int
end

"""
    DeflateCompressor(;level=$(Z_DEFAULT_COMPRESSION), windowbits=$(Z_DEFAULT_COMPRESSION))

Create a deflate compression codec.

Arguments
---------
- `level`: compression level (-1..9)
- `windowbits`: size of history buffer (8..15)
"""
function DeflateCompressor(;level::Integer=Z_DEFAULT_COMPRESSION,
                        windowbits::Integer=Z_DEFAULT_WINDOWBITS)
    if !(-1 ≤ level ≤ 9)
        throw(ArgumentError("compression level must be within -1..9"))
    elseif !(8 ≤ windowbits ≤ 15)
        throw(ArgumentError("windowbits must be within 8..15"))
    end
    return DeflateCompressor(ZStream(), level, -Int(windowbits))
end

const DeflateCompressorStream{S} = TranscodingStream{DeflateCompressor,S} where S<:IO

"""
    DeflateCompressorStream(stream::IO; kwargs...)

Create a deflate compression stream (see `DeflateCompressor` for `kwargs`).
"""
function DeflateCompressorStream(stream::IO; kwargs...)
    x, y = splitkwargs(kwargs, (:level, :windowbits))
    return TranscodingStream(DeflateCompressor(;x...), stream; y...)
end

# Methods
# -------

function TranscodingStreams.initialize(codec::CompressorCodec)
    code = deflate_init!(codec.zstream, codec.level, codec.windowbits)
    if code != Z_OK
        zerror(codec.zstream, code)
    end
    return
end

function TranscodingStreams.finalize(codec::CompressorCodec)
    zstream = codec.zstream
    if zstream.state != C_NULL
        code = deflate_end!(zstream)
        if code != Z_OK
            zerror(zstream, code)
        end
    end
    return
end

function TranscodingStreams.startproc(codec::CompressorCodec, state::Symbol, error::Error)
    code = deflate_reset!(codec.zstream)
    if code == Z_OK
        return :ok
    else
        error[] = ErrorException(zlib_error_message(codec.zstream, code))
        return :error
    end
end

function TranscodingStreams.process(codec::CompressorCodec, input::Memory, output::Memory, error::Error)
    zstream = codec.zstream
    zstream.next_in = input.ptr
    zstream.avail_in = input.size
    zstream.next_out = output.ptr
    zstream.avail_out = output.size
    code = deflate!(zstream, input.size > 0 ? Z_NO_FLUSH : Z_FINISH)
    Δin = Int(input.size - zstream.avail_in)
    Δout = Int(output.size - zstream.avail_out)
    if code == Z_OK
        return Δin, Δout, :ok
    elseif code == Z_STREAM_END
        return Δin, Δout, :end
    else
        error[] = ErrorException(zlib_error_message(zstream, code))
        return Δin, Δout, :error
    end
end


text:
include("../src/fdfd_local_field.jl")

using Plots
pyplot()

ps = transpose(0.1:0.01:0.9)

E_local_field = [get_local_field(ps[:,i])[3] for i=1:length(ps)]

Plots.plot(real.(E_local_field))
Plots.plot!(imag.(E_local_field))
gui()

Plots.plot(angle.(E_local_field))
gui()


text:
Qx, x = PolynomialRing(FlintQQ, "x")
QasNumberField, _ = NumberField(x - 1)
Kt, t = PolynomialRing(QasNumberField, "t")
K1, a1 = NumberField(x^3 - 2)
K2, (a2, ) = NumberField([x^3 - 2])
K3, a3 = NumberField(t^3 - 2)
K4, (a4, ) = NumberField([t^3 - 2])

@testset for (K, a) in [(K1, a1), (K2, a2), (K3, a3), (K4, a4)]
  b = one(K)
  fl = @inferred isintegral(b)
  @test fl
  fl = @inferred isintegral(a)

  B = @inferred basis(K)
  c = @inferred basis_matrix([one(K), a^4])
  @assert nrows(c) == 2
  @assert ncols(c) == 3
  @assert sum(B[i] * c[1, i] for i in 1:3) == one(K)
  @assert sum(B[i] * c[2, i] for i in 1:3) == a^4

  f = @inferred charpoly(one(K()))
  @test f == (gen(parent(f)) - 1)^3
  f = @inferred minpoly(one(K))
  @test f == (gen(parent(f)) - 1)

  f = @inferred charpoly(a)
  @test f == gen(parent(f))^3 - 2
  f = @inferred minpoly(a)
  @test f == gen(parent(f))^3 - 2

end


text:
"""
CovMoment(f; m = 0)

Calculate the spectural densitiy matrix - p.190 equation (145) in Cochrane (2000)

INPUT
`f`: T x q array of q moment conditions
`cov_method': method for calculating the variance-covariance matrix of the error terms:
            - time iid: error terms are iid across time but cross-sectionally correlated
            - iid: error terms are iiid both across time and assets
            - NW: Newey-West coveriance matrix with m lags
            - HH: Hansen-Hodrick covariance matrix with m lag
`m': number of lags for the Newey-West or Hansen-Hodrick covariance matrix
OUTPUT
- `Ef2`: q x q covariance matrix
"""

function CovMoment(f, cov_method = "time iid", m = 0)
    Ef2 = (1 / (T-1)) * (f' * f);
    if cov_method == "time iid"
        Ef2 = Ef2; # Spectural densitiy matrix under time iid errors
    elseif cov_method == "iid"
        Ef2 = diagm(diag(Ef2)); # Spectural densitiy matrix under iid errors
    elseif cov_method == "NW"
        for i = 1:m
            w = 1 - i / (m + 1);
            Ef2 = Ef2 + 2 * w * (1 / (T - 1 - i)) * f[1:end - i, :]' * f[1 + i:end, :]; # Spectural densitiy matrix under Newey-West errors
        end
    elseif cov_method == "HH"
        for i = 1:m
            Ef2 = Ef2 + (2 / (T - 1 - i)) * f[1:end - i, :]' * f[1 + i:end, :]; # Spectural densitiy matrix under Hansen-Hodrick errors
        end
    end

    return Ef2

end


text:
@testset "HTTP.Client" begin

@testset "HTTP.Connection" begin
    conn = HTTP.Connection(IOBuffer())
    @test conn.state == HTTP.Busy
    HTTP.idle!(conn)
    @test conn.state == HTTP.Idle
    HTTP.busy!(conn)
    @test conn.state == HTTP.Busy
    HTTP.dead!(conn)
    @test conn.state == HTTP.Dead
    HTTP.idle!(conn)
    @test conn.state == HTTP.Dead
    HTTP.busy!(conn)
    @test conn.state == HTTP.Dead
end

for sch in ("http", "https")
    println("runnin $sch client tests...")

    println("simple GET, HEAD, POST, DELETE, etc.")
    @test HTTP.status(HTTP.get("$sch://httpbin.org/ip")) == 200
    @test HTTP.status(HTTP.head("$sch://httpbin.org/ip")) == 200
    @test HTTP.status(HTTP.options("$sch://httpbin.org/ip")) == 200
    @test HTTP.status(HTTP.post("$sch://httpbin.org/ip"; statusraise=false)) == 405
    @test HTTP.status(HTTP.post("$sch://httpbin.org/post")) == 200
    @test HTTP.status(HTTP.put("$sch://httpbin.org/put")) == 200
    @test HTTP.status(HTTP.delete("$sch://httpbin.org/delete")) == 200
    @test HTTP.status(HTTP.patch("$sch://httpbin.org/patch")) == 200

    # Testing within tasks, see https://github.com/JuliaWeb/HTTP.jl/issues/18
    println("async client request")
    @test HTTP.status(wait(@schedule HTTP.get("$sch://httpbin.org/ip"))) == 200

    @test HTTP.status(HTTP.get("$sch://httpbin.org/encoding/utf8")) == 200

    println("pass query to uri")
    r = HTTP.get("$sch://httpbin.org/response-headers"; query=Dict("hey"=>"dude"))
    h = HTTP.headers(r)
    @test (haskey(h, "Hey") ? h["Hey"] == "dude" : h["hey"] == "dude")

    println("cookie requests")
    r = HTTP.get("$sch://httpbin.org/cookies")
    body = String(take!(r))
    @test (body == "{\n  \"cookies\": {}\n}\n" || body == "{\n  \"cookies\": {\n    \"hey\": \"\"\n  }\n}\n" || body == "{\n  \"cookies\": {\n    \"hey\": \"sailor\"\n  }\n}\n")
    r = HTTP.get("$sch://httpbin.org/cookies/set?hey=sailor")
    @test HTTP.status(r) == 200
    body = String(take!(r))
    @test (body == "{\n  \"cookies\": {\n    \"hey\": \"sailor\"\n  }\n}\n" || body == "{\n  \"cookies\": {\n    \"hey\": \"\"\n  }\n}\n")

    # r = HTTP.get("$sch://httpbin.org/cookies/delete?hey")
    # @test String(take!(r)) == "{\n  \"cookies\": {\n    \"hey\": \"\"\n  }\n}\n"

    # stream
    println("client streaming tests")
    r = HTTP.post("$sch://httpbin.org/post"; body="hey")
    @test HTTP.status(r) == 200
    # stream, but body is too small to actually stream
    r = HTTP.post("$sch://httpbin.org/post"; body="hey", stream=true)
    @test HTTP.status(r) == 200
    r = HTTP.get("$sch://httpbin.org/stream/100")
    @test HTTP.status(r) == 200
    totallen = length(HTTP.body(r)) # number of bytes to expect
    bytes = take!(r)
    begin
        r = HTTP.get("$sch://httpbin.org/stream/100"; stream=true)
        @test HTTP.status(r) == 200
        len = length(HTTP.body(r))
        HTTP.@timeout 15.0 begin
            while !eof(HTTP.body(r))
                b = take!(r)
            end
        end throw(HTTP.TimeoutException(15.0))
    end

    # body posting: Vector{UInt8}, String, IOStream, IOBuffer, FIFOBuffer
    println("client body posting of various types")
    @test HTTP.status(HTTP.post("$sch://httpbin.org/post"; body="hey")) == 200
    @test HTTP.status(HTTP.post("$sch://httpbin.org/post"; body=UInt8['h','e','y'])) == 200
    io = IOBuffer("hey"); seekstart(io)
    @test HTTP.status(HTTP.post("$sch://httpbin.org/post"; body=io)) == 200
    tmp = tempname()
    open(f->write(f, "hey"), tmp, "w")
    io = open(tmp)
    @test HTTP.status(HTTP.post("$sch://httpbin.org/post"; body=io)) == 200
    close(io); rm(tmp)
    f = HTTP.FIFOBuffer("hey")
    @test HTTP.status(HTTP.post("$sch://httpbin.org/post"; body=f)) == 200

    # chunksize
    println("client transfer-encoding chunked")
    @test_broken HTTP.status(HTTP.post("$sch://httpbin.org/post"; body="hey", chunksize=2)) == 200
    @test_broken HTTP.status(HTTP.post("$sch://httpbin.org/post"; body=UInt8['h','e','y'], chunksize=2)) == 200
    io = IOBuffer("hey"); seekstart(io)
    @test_broken HTTP.status(HTTP.post("$sch://httpbin.org/post"; body=io, chunksize=2)) == 200
    tmp = tempname()
    open(f->write(f, "hey"), tmp, "w")
    io = open(tmp)
    @test_broken HTTP.status(HTTP.post("$sch://httpbin.org/post"; body=io, chunksize=2)) == 200
    close(io); rm(tmp)
    f = HTTP.FIFOBuffer("hey")
    @test_broken HTTP.status(HTTP.post("$sch://httpbin.org/post"; body=f, chunksize=2)) == 200

    # multipart
    println("client multipart body")
    r = HTTP.post("$sch://httpbin.org/post"; body=Dict("hey"=>"there"))
    @test HTTP.status(r) == 200
    @test startswith(String(take!(r)), "{\n  \"args\": {}, \n  \"data\": \"\", \n  \"files\": {}, \n  \"form\": {\n    \"hey\": \"there\"\n  }")

    r = HTTP.post("$sch://httpbin.org/post"; body=Dict("hey"=>"there"), chunksize=1000)
    @test HTTP.status(r) == 200
    @test startswith(String(take!(r)), "{\n  \"args\": {}, \n  \"data\": \"\", \n  \"files\": {}, \n  \"form\": {\n    \"hey\": \"there\"\n  }")

    tmp = tempname()
    open(f->write(f, "hey"), tmp, "w")
    io = open(tmp)
    r = HTTP.post("$sch://httpbin.org/post"; body=Dict("hey"=>"there", "iostream"=>io))
    close(io); rm(tmp)
    @test HTTP.status(r) == 200
    str = String(take!(r))
    @test startswith(str, "{\n  \"args\": {}, \n  \"data\": \"\", \n  \"files\": {\n    \"iostream\": \"hey\"\n  }, \n  \"form\": {\n    \"hey\": \"there\"\n  }")

    tmp = tempname()
    open(f->write(f, "hey"), tmp, "w")
    io = open(tmp)
    r = HTTP.post("$sch://httpbin.org/post"; body=Dict("hey"=>"there", "iostream"=>io), chunksize=1000)
    close(io); rm(tmp)
    @test HTTP.status(r) == 200
    @test startswith(String(take!(r)), "{\n  \"args\": {}, \n  \"data\": \"\", \n  \"files\": {\n    \"iostream\": \"hey\"\n  }, \n  \"form\": {\n    \"hey\": \"there\"\n  }")

    tmp = tempname()
    open(f->write(f, "hey"), tmp, "w")
    io = open(tmp)
    m = HTTP.Multipart("mycoolfile.txt", io)
    r = HTTP.post("$sch://httpbin.org/post"; body=Dict("hey"=>"there", "multi"=>m))
    close(io); rm(tmp)
    @test HTTP.status(r) == 200
    @test startswith(String(take!(r)), "{\n  \"args\": {}, \n  \"data\": \"\", \n  \"files\": {\n    \"multi\": \"hey\"\n  }, \n  \"form\": {\n    \"hey\": \"there\"\n  }")

    tmp = tempname()
    open(f->write(f, "hey"), tmp, "w")
    io = open(tmp)
    m = HTTP.Multipart("mycoolfile", io, "application/octet-stream")
    r = HTTP.post("$sch://httpbin.org/post"; body=Dict("hey"=>"there", "multi"=>m), chunksize=1000)
    close(io); rm(tmp)
    @test HTTP.status(r) == 200
    @test startswith(String(take!(r)), "{\n  \"args\": {}, \n  \"data\": \"\", \n  \"files\": {\n    \"multi\": \"hey\"\n  }, \n  \"form\": {\n    \"hey\": \"there\"\n  }")

    # asynchronous
    println("asynchronous client request body")
    begin
        f = HTTP.FIFOBuffer()
        write(f, "hey")
        t = @async HTTP.post("$sch://httpbin.org/post"; body=f)
        wait(f) # wait for the async call to write it's first data
        write(f, " there ") # as we write to f, it triggers another chunk to be sent in our async request
        write(f, "sailor")
        close(f) # setting eof on f causes the async request to send a final chunk and return the response
        @test HTTP.status(wait(t)) == 200
    end

    # redirects
    println("client redirect following")
    r = HTTP.get("$sch://httpbin.org/redirect/1")
    @test HTTP.status(r) == 200
    @test length(HTTP.history(r)) == 1
    @test_throws HTTP.RedirectException HTTP.get("$sch://httpbin.org/redirect/6")
    @test HTTP.status(HTTP.get("$sch://httpbin.org/relative-redirect/1")) == 200
    @test HTTP.status(HTTP.get("$sch://httpbin.org/absolute-redirect/1")) == 200
    @test HTTP.status(HTTP.get("$sch://httpbin.org/redirect-to?url=http%3A%2F%2Fexample.com")) == 200

    @test HTTP.status(HTTP.post("$sch://httpbin.org/post"; body="√")) == 200
    println("client basic auth")
    @test HTTP.status(HTTP.get("$sch://user:pwd@httpbin.org/basic-auth/user/pwd")) == 200
    @test HTTP.status(HTTP.get("$sch://user:pwd@httpbin.org/hidden-basic-auth/user/pwd")) == 200

    # readtimeout
    # @test_throws HTTP.TimeoutException HTTP.get("$sch://httpbin.org/delay/3"; readtimeout=1.0)

    # custom client & other high-level entries
    println("high-level client request methods")
    buf = IOBuffer()
    cli = HTTP.Client(buf)
    HTTP.get(cli, "$sch://httpbin.org/ip")
    seekstart(buf)
    @test length(String(take!(buf))) > 0

    r = HTTP.request("$sch://httpbin.org/ip")
    @test HTTP.status(r) == 200

    uri = HTTP.URI("$sch://httpbin.org/ip")
    r = HTTP.request(uri)
    @test HTTP.status(r) == 200
    r = HTTP.get(uri)
    @test HTTP.status(r) == 200
    r = HTTP.get(cli, uri)
    @test HTTP.status(r) == 200

    r = HTTP.request(HTTP.GET, "$sch://httpbin.org/ip")
    @test HTTP.status(r) == 200

    uri = HTTP.URI("$sch://httpbin.org/ip")
    r = HTTP.request("GET", uri)
    @test HTTP.status(r) == 200

    req = HTTP.Request(HTTP.GET, uri, HTTP.Headers(), HTTP.FIFOBuffer())
    r = HTTP.request(req)
    @test HTTP.status(r) == 200
    @test !HTTP.isnull(HTTP.request(r))
    @test length(take!(r)) > 0

    for c in HTTP.DEFAULT_CLIENT.httppool["httpbin.org"]
        HTTP.dead!(c)
    end

    r = HTTP.get(cli, "$sch://httpbin.org/ip")
    @test isempty(HTTP.cookies(r))
    @test isempty(HTTP.history(r))

    r = HTTP.get("$sch://httpbin.org/image/png")
    @test HTTP.status(r) == 200

    # r = HTTP.connect("http://47.89.41.164:80")
    # gzip body = "hey"
    # body = UInt8[0x1f,0x8b,0x08,0x00,0x00,0x00,0x00,0x00,0x00,0x03,0xcb,0x48,0xad,0x04,0x00,0xf0,0x15,0xd6,0x88,0x03,0x00,0x00,0x00]
    # r = HTTP.post("$sch://httpbin.org/post"; body=body, chunksize=1)
end

end # @testset "HTTP.Client"


text:
function make_grid(Ny, Nz, yLength, zLength)
    return RegularCartesianGrid(topology = (Bounded, Bounded, Bounded),
                                size=(1, Ny, Nz),
                                x=(-π, π),
                                y=(-yLength/2, yLength/2),
                                z=(-zLength, 0))
end


text:
using KdotP
using Test
using Crystalline
using LinearAlgebra

# ---------------------------------------------------------------------------------------- #
# Plane group 17, K-point (Dirac cone)
lgirsd = lgirreps(17, Val(2))
lgirs = lgirsd["K"][end] # K₃ (a 2D irrep)

H = kdotp(lgirs)
H((.1,.2))

using PlotlyJS
kxys = range(-.2, .2, 50)
Es = [eigvals(H((kx,ky))) for kx in kxys, ky in kxys]
plot(map(i->PlotlyJS.surface(z=getindex.(Es, i)), 1:irdim(H)))

# ---------------------------------------------------------------------------------------- #
# Space group 222, R-point (6-fold degeneracy)
lgirsd = lgirreps(222, Val(3))
lgirs = lgirsd["R"][end] # R₄ (a 6D irrep)
H = kdotp(lgirs)

using PlotlyJS
kxys = range(-.2, .2, 50)
Es = [eigvals(H((kx,ky,ky), 1)) for kx in kxys, ky in kxys]
plot(map(i->PlotlyJS.surface(z=getindex.(Es, i), surfacecolor=1000*i/N*ones(50^2)'), 1:irdim(H)))

using Brillouin
pts = [-[1,1,1].*.2, [0,0,0], [-1,0,-1]*.2] # from Γ, to R, to X (relative to R)
ks = Brillouin.interpolate(pts, 50)
kdists = Brillouin.cumdists(ks)

Es = [eigvals(H(k, [1,0])) for k in ks]
plot(map(i->PlotlyJS.scatter(x=kdists, y=getindex.(Es, i)), 1:irdim(H)))

@testset begin "functor"
    # 3D
    lgir = lgirreps(222,3)["R"][end]
    H = kdotp(lgir)
    k = [.1,.2,-.3]
    @test (.9H(k, 1) - .7H(k,2)) ≈ H(k, [.9,-.7])
end

# ---------------------------------------------------------------------------------------- #
# Space group 139, P-point (2D irrep)
# TODO: BROKEN FOR BOTH P and Γ
lgirsd = lgirreps(139, Val(3))
lgirs = lgirsd["P"][end] # P₅ (a 2D irrep)
H = kdotp(lgirs) # BROKEN: only depends on k₃

lgirs = lgirsd["Γ"][end] # Γ₅⁻ (a 2D irrep)
H = kdotp(lgirs) # BROKEN: no elements

# ---------------------------------------------------------------------------------------- #
# Space group 225, W-point (2D irrep)
# TODO: BROKEN FOR BOTH W and Γ
lgirsd = lgirreps(225, Val(3))
lgirs = lgirsd["W"][end] # P₅ (a 2D irrep)
H = kdotp(lgirs) # BROKEN: only depends on k₁

lgirs = lgirsd["Γ"][end] # Γ₅⁻ (a 2D irrep)
H = kdotp(lgirs) # BROKEN: no elements

# ---------------------------------------------------------------------------------------- #
lgirs = lgirreps(92, Val(3))["A"]
lgir = lgirs[1] # P₅ (a 2D irrep)
H = kdotp(lgir) # BROKEN: only depends on k₁


text:
"""
    gpw(ζ, φ; k=1)

Closure that defines a generalized plane wave in polar coordinates.

The closure captures the values of the evanescence parameter ``ζ``, the
angle ``φ`` and the wavenumber ``k``.
It does not include any normalization.
The function can then be evaluated at any ``(r,θ)`` point.

The generalized plane wave in polar coordinates is defined by
```math
ϕ := (r,θ) ↦ e^{ı k \\mathbf{d} ⋅ \\mathbf{x}},
```
where
```math
\\mathbf{d} = (\\cos[φ+ıζ], \\sin[φ+ıζ]),
\\qquad\\text{and}\\qquad
\\mathbf{x} = (r \\cos θ, r\\sin θ).
```
"""
function gpw(ζ,φ;k=1)
    d = [cos(φ+im*ζ), sin(φ+im*ζ)] # Direction
    return (r,θ) -> exp(im * k * [r*cos(θ), r*sin(θ)] ⋅ d)
end

"""
    approximation_set(Y, W; k=1)

Generate a set of generalized plane waves.

The parameters `Y` and `W` are respectively the sets of complex angles
``(ζ_j,φ_j)`` and associated weights ``ω_j`` corresponding to the generalized
plane waves, see [`gpw`](@ref)
```math
ϕ_j := (r,θ) ↦ e^{ı k \\mathbf{d}_j ⋅ \\mathbf{x}},
```
where
```math
\\mathbf{d}_j = (\\cos[φ_j+ıζ_j], \\sin[φ_j+ıζ_j]),
\\qquad\\text{and}\\qquad
\\mathbf{x} = (r \\cos θ, r\\sin θ).
```
Approximations are constructed in the form of
```math
(r, θ) ↦ \\sum_j ξ_j ω_j ϕ_j(r, θ),
```
where ``(ξ_j)_j`` is the set of unknown coefficients.

The value of the wavenumber `k` defaults to ``1`` and can be provided as an
optional parameter.
"""
function approximation_set(Y, W; k=1)
    @assert length(Y) == length(W)
    return [(r, θ) -> w * gpw(y...;k=k)(r, θ) for (y, w) in zip(Y, W)]
end

"""
    approximation_set(N; k=1)

Generate a set of ``N`` propagative plane waves with equispaced angles.

The equispaced angles are defined by
```math
θ_n := 2π s / R, \\qquad n = 0,…,N-1.
```

The value of the wavenumber `k` defaults to ``1`` and can be provided as an
optional parameter.
"""
function approximation_set(N; k=1)
    Y = [(0, 2π * n / N) for n in 0:N-1]
    W = ones(Int64, N) ./ sqrt(N)
    return approximation_set(Y, W; k=k)
end

"""
    approximation_set(N, qs, smpl_type; k=1)
    approximation_set(N, Q::Integer, smpl_type; k=1)

Generate a set of ``N`` evanescent plane waves according to some sampling type.

The second parameter can be an integer ``Q`` which controls the truncation
parameter of the [`TruncKernel`](@ref), or a set of integers `qs` to specify
only some specific terms in the expansion.
The parameter `smpl_type` should be a function name: see
[`uniform_sampling`](@ref), [`sobol_sampling`](@ref) and
[`random_sampling`](@ref).

The value of the wavenumber `k` defaults to ``1`` and can be provided as an
optional parameter.
"""
function approximation_set(N, qs, smpl_type; k=1)
    K = TruncKernel(Tuple(qs), logwz(;k=k));
    smpl = Sampler(K)
    Y, W = smpl_type(smpl, N)
    W ./= sqrt(N)
    return approximation_set(Y, W; k=k)
end
function approximation_set(N, Q::Int64, smpl_type; k=1)
    return approximation_set(N, -Q:Q, smpl_type; k=k)
end

max_stars_repo_name: EmileParolin/STrAW

max_issues_repo_name: EmileParolin/STrAW

max_forks_repo_name: EmileParolin/STrAW


text:
using Test

include("../src/common.jl")
include("../src/uninformed_search.jl")
include("../src/a_star.jl")

@testset "AStar" begin
  include("reachnumber.jl")
  include("maze.jl")
  include("15puzzle.jl")
  include("collatz.jl")
  include("no_path.jl")
end


text:
"""
    writeMetric( mesh::Mesh2D, mesher::MesherFEFLOA,
                      metric::Vector{Float64}, flname::String )

Writes 2D isotropic metric to feflo.a format in `flname`.
"""
function writeMetric( mesh::Mesh2D, mesher::MesherFEFLOA,
                      metric::Vector{Float64}, flname::String )

    np = size(mesh.p,1)
    if np != length(metric)
        error(" Isotropic metric error: size of metric does not correspond with size of mesh")
    end

    fid = open( flname, "w" )

    @printf( fid, "MeshVersionFormatted 2\n\n" )

    @printf( fid, "Dimension 2\n\n" )

    @printf( fid, "SolAtVertices\n" )

    @printf( fid, "%i\n", np )
    @printf( fid, "1 1\n" )

    # Vertices
    for jj = 1:np
        @printf(fid, "%16.15e\n", metric[jj])
    end

    close(fid)

    # TODO: Corners

end

"""
    writeMetric( mesh::Mesh2D, mesher::MesherFEFLOA,
                      metric::Matrix{Float64}, flname::String )

Writes 2D anisotropic metric to feflo.a format in `flname`.
"""
function writeMetric( mesh::Mesh2D, mesher::MesherFEFLOA,
                      metric::Matrix{Float64}, flname::String )

    np = size(mesh.p,1)
    if np != size(metric,1)
        error(" Ansotropic metric error: size of metric does not correspond with size of mesh")
    end
    if size(metric,2) != 3
        error(" Ansotropic metric error: wrong number of entries for metric")
    end

    fid = open( flname, "w" )

    @printf( fid, "MeshVersionFormatted 2\n\n" )

    @printf( fid, "Dimension 2\n\n" )

    @printf( fid, "SolAtVertices\n" )

    @printf( fid, "%i\n", np )
    @printf( fid, "1 3\n" )

    # Vertices
    for jj = 1:np
        @printf(fid, "%16.15e\t%16.15e\t%16.15e\n", metric[jj,1], metric[jj,2], metric[jj,3])
    end

    close(fid)

    # TODO: Corners

end

"""
    writeMetric( mesh::Mesh3D, mesher::MesherFEFLOA,
                      metric::Vector{Float64}, flname::String )

Writes 3D isotropic metric to feflo.a format in `flname`.
"""
function writeMetric( mesh::Mesh3D, mesher::MesherFEFLOA,
                      metric::Vector{Float64}, flname::String )

    np = size(mesh.p,1)
    if np != length(metric)
        error(" Isotropic metric error: size of metric does not correspond with size of mesh")
    end

    fid = open( flname, "w" )

    @printf( fid, "MeshVersionFormatted 2\n\n" )

    @printf( fid, "Dimension 3\n\n" )

    @printf( fid, "SolAtVertices\n" )

    @printf( fid, "%i\n", np )
    @printf( fid, "1 1\n" )

    # Vertices
    for jj = 1:np
        @printf(fid, "%16.15e\n", metric[jj])
    end

    close(fid)

    # TODO: Corners

end

"""
    writeMetric( mesh::Mesh3D, mesher::MesherFEFLOA,
                      metric::Matrix{Float64}, flname::String )

Writes 3D anisotropic metric to feflo.a format in `flname`.
"""
function writeMetric( mesh::Mesh3D, mesher::MesherFEFLOA,
                      metric::Matrix{Float64}, flname::String )

    np = size(mesh.p,1)
    if np != size(metric,1)
        error(" Anisotropic metric error: size of metric does not correspond with size of mesh")
    end
    if size(metric,2) != 6
        error(" Anisotropic metric error: wrong number of entries for metric")
    end

    fid = open( flname, "w" )

    @printf( fid, "MeshVersionFormatted 2\n\n" )

    @printf( fid, "Dimension 3\n\n" )

    @printf( fid, "SolAtVertices\n" )

    @printf( fid, "%i\n", np )
    @printf( fid, "1 3\n" )

    # Vertices
    for jj = 1:np
        @printf(fid, "%16.15e\t%16.15e\t%16.15e\t%16.15e\t%16.15e\t%16.15e\n",
         metric[jj,1], metric[jj,2], metric[jj,3],
         metric[jj,4], metric[jj,5], metric[jj,6])
         # Ordered (1,1), (2,1), (2,2), (3,1), (3,2), (3,3)
    end

    close(fid)

    # TODO: Corners

end


text:
####################
### Split output ###
####################

type OutputSplit{F} <: Output
    cellhandles::Dict{UInt,FASTQ.Writer}
    umihandles::Dict{UInt,IOStream}
    cellpaths::Dict{UInt,String}
    namegen::F
    maxopenfiles::Int
    outputdir::String
end

function OutputSplit(protocol::Protocol;
                     outputdir::String = ".",
                     maxopenfiles::Int = 1000,
                     kwargs...)

    mkpath(outputdir)

    namegen(cellid, unmatched) =
        joinpath(outputdir,
                 unmatched ? "unmatched" : String(cellid))

    return OutputSplit(Dict{UInt,FASTQ.Writer}(),
                       Dict{UInt,IOStream}(),
                       Dict{UInt,String}(),
                       namegen,
                       maxopenfiles,
                       outputdir)
end

"""

    Writes a read (`InterpretedRecord`) to corresponding output file(s).

    """
function Base.write(oh::OutputSplit, ir::InterpretedRecord)
    if !ir.unmatched
        closetoomanyfiles!(oh)
        writeinsert(oh,ir)
        writeumi(oh,ir)
        addpath(oh,ir)
    end
end

"""

    Counts the number of open files and closes them if necessary.  Should
    be made more intelligent in the future.

    """
function closetoomanyfiles!{N}(oh::OutputSplit{N})
    if length(oh.cellhandles) + length(oh.umihandles) >= oh.maxopenfiles
        # close only the cell-like handles, leave the raw handles open
        map(close, values(oh.cellhandles))
        map(close, values(oh.umihandles))

        oh.cellhandles=Dict()
        oh.umihandles=Dict()
    end
    return nothing
end

function writeinsert(oh::OutputSplit,ir::InterpretedRecord)

    writeid = ir.unmatched ? UInt(0) : ir.cellid.val

    celldesc = get!(oh.cellhandles, writeid) do
        filename = oh.namegen(ir.cell, ir.unmatched)*".fastq"
        FASTQ.Writer(open(filename, "w+"))
    end
    write(celldesc,ir.output)

    return nothing
end

function writeumi(oh::OutputSplit,ir::InterpretedRecord)

    writeid = ir.unmatched ? UInt(0) : ir.cellid.val
    umidesc = get!(oh.umihandles, writeid) do
        open(oh.namegen(ir.cell, ir.unmatched)*".umi", "w+")
    end
    unsafe_write(umidesc,pointer(ir.umi),length(ir.umi))
    write(umidesc,'\n')

    return nothing
end

function addpath(oh::OutputSplit,ir::InterpretedRecord)
    writeid = ir.unmatched ? UInt(0) : ir.cellid.val
    get!(oh.cellpaths, writeid) do
        oh.namegen(ir.cell, ir.unmatched)
    end
    return nothing
end

function mergeoutput{N}(outputs::Vector{OutputSplit{N}};
                        outputdir::String = ".",
                        kwargs...)

    # 1) merge the cell barcodes
    paths = Dict{UInt,Vector{String}}()
    for o in outputs
        for (id,path) in o.cellpaths
            pbase = get!(paths, id) do
                String[]
            end
            push!(pbase,path)
        end
    end

    # 2) for each cellid concatenate the contents of all files into
    # one and remove the directories
    pmap(values(paths)) do ps
        mergecellid(outputdir,ps,".umi")
        mergecellid(outputdir,ps,".fastq")
    end

    # 3) clean up the directories
    pmap(outputs) do oh
        rm(oh.outputdir,force=true,recursive=true)
    end
end

function mergecellid(outputdir,path,ext)
    outfile = joinpath(outputdir,basename(path[1])*ext)
    infiles = map(f->f*ext,path)
    catfiles(outfile,infiles)
    map(rm, infiles)
end

function Base.close(oh::OutputSplit)
    map(close, values(oh.cellhandles))
    map(close, values(oh.umihandles))
end

max_stars_repo_name: pwl/FASTQDemultiplexer

max_issues_repo_name: pwl/FASTQDemultiplexer

max_forks_repo_name: pwl/FASTQDemultiplexer


text:
# From data (which is an array with planet number,
# epoch, time of observation, error), return
# the epochs, times, errors, and number of observations
# for the ith planet:
function extract_planet(data,i)
eobs = zeros(0)
tobs = zeros(0)
sobs = zeros(0)
nobs = 0
for j=1:size(data,1)
  if data[j,1] == i
# Extract data:
    if data[j,4] != 0.0
      push!(eobs,data[j,2])
      push!(tobs,data[j,3])
      push!(sobs,data[j,4])
      nobs += 1
    end
  end
end
return eobs,tobs,sobs,nobs
end


text:
using Test, Omniscape, Circuitscape

@testset "Internals" begin
### Unit tests for components
## source target matching
source_subset = [1.0 1 1 1 1 1 1 1 1 1; # 6x10 array
                   1 1 1 1 1 1 1 1 1 1;
                   1 1 1 1 1 1 1 1 1 1;
                   1 1 1 1 1 1 1 1 1 1;
                   1 1 1 1 1 1 1 1 1 1;
                   1 1 1 1 1 1 1 1 1 1]

con1pres =   [1.0 4 3 2 7 6 9 2 3 0;
                4 2 6 5 3 9 4 5 3 3;
                1 1 5 2 7 7 8 8 9 1;
                8 7 6 5 5 5 0 0 1 1;
                3 5 4 4 4 8 5 0 1 8;
                1 2 8 2 1 6 5 5 8 7]

con1fut =   [2.0 1 3 2 7 6 7 1 3 0;
               4 2 4 5 3 9 4 5 5 3;
               3 8 3 2 6 3 8 5 9 1;
               8 7 6 5 1 5 0 0 1 1;
               3 1 4 5 4 8 5 7 1 8;
               5 2 3 2 1 6 4 5 3 7]
x = 5
y = 4
con1_lower = -1.0
con1_upper = 1.0
source_subset = convert(Array{Union{Float64, Missing}, 2}, source_subset)

condition_layers = Omniscape.ConditionLayers(
    convert(Array{Union{Float64, Missing}, 2}, con1pres),
    convert(Array{Union{Float64, Missing}, 2}, con1fut),
    Array{Union{Float64, Missing}, 2}(undef, 1, 1),
    Array{Union{Float64, Missing}, 2}(undef, 1, 1)
)

Omniscape.source_target_match!(source_subset,
                     1,
                     condition_layers,
                     Omniscape.Conditions("within", "within", con1_lower, con1_upper, 0.0, 0.0),
                     y,
                     y,
                     x,
                     x,
                     1,
                     6,
                     1,
                     10
                     )
target_val = con1fut[y, x]
# Make sure no present day vals in con1 are outside of range
@test sum((con1pres[source_subset .== 1] .< (target_val + con1_lower)) .|
    (con1pres[source_subset .== 1] .> (target_val + con1_upper))) == 0
@info "conditional connectivity tests passed"

## Check that targets are IDed properly
source_strength = Omniscape.read_raster("input/source.asc", Float64)[1]

int_arguments = Dict{String, Int64}()
int_arguments["block_size"] = 7
int_arguments["block_radius"] = 3 # must be (size - 1) / 2
int_arguments["nrows"] = size(source_strength)[1]
int_arguments["ncols"] = size(source_strength)[2]
targets = Omniscape.get_targets(source_strength,
                                int_arguments,
                                Float64)

n_targets = floor(int_arguments["nrows"] / int_arguments["block_size"]) *
    floor(int_arguments["ncols"] / int_arguments["block_size"])

# Correct number of targets
@test size(targets)[1] == n_targets # would be less than if any of the sources
                                    # had 0 strength, but sources.asc does
                                    # not have 0's

# correct source strength for a target
block_sources = source_strength[Int(targets[1,2] - int_arguments["block_radius"]):Int(targets[1,2] + int_arguments["block_radius"]),
                            Int(targets[1,1] - int_arguments["block_radius"]):Int(targets[1,1] + int_arguments["block_radius"])]
@test targets[1,3] ≈ sum(block_sources)
@info "target tests passed"

# Test error throws
@info "Testing error throws"
@test run_omniscape("input/config7.ini") === nothing
@test run_omniscape("input/bad_config.ini") === nothing

end

@testset "run_omnsicape()" begin
### Tests for run_omniscape()
l, f, p = run_omniscape("input/config4.ini")
l_verify = Omniscape.read_raster("output_verify/test4/cum_currmap.tif", Float64)[1]
f_verify = Omniscape.read_raster("output_verify/test4/flow_potential.tif", Float64)[1]
p_verify = Omniscape.read_raster("output_verify/test4/normalized_cum_currmap.tif", Float64)[1]
@test Omniscape.arrays_equal(l, l_verify)
@test Omniscape.arrays_equal(f, f_verify)
@test Omniscape.arrays_equal(p, p_verify)

g = run_omniscape("input/config5.ini")
g_verify = Omniscape.read_raster("output_verify/test5/cum_currmap.tif", Float64)[1]
@test Omniscape.arrays_equal(g, g_verify)

h = run_omniscape("input/config6.ini")
h_verify = Omniscape.read_raster("output_verify/test6/cum_currmap.tif", Float64)[1]
@test Omniscape.arrays_equal(h, h_verify)

a, b, c = run_omniscape("input/config.ini")
a_verify = Omniscape.read_raster("output_verify/test1/cum_currmap.tif", Float64)[1]
b_verify = Omniscape.read_raster("output_verify/test1/flow_potential.tif", Float64)[1]
c_verify = Omniscape.read_raster("output_verify/test1/normalized_cum_currmap.tif", Float64)[1]
@test Omniscape.arrays_equal(a, a_verify)
@test Omniscape.arrays_equal(b, b_verify)
@test Omniscape.arrays_equal(c, c_verify)

a1, b1, c1 = run_omniscape("input/config_32bit.ini")
a1_verify = Omniscape.read_raster("output_verify/test1_32bit/cum_currmap.tif", Float32)[1]
b1_verify = Omniscape.read_raster("output_verify/test1_32bit/flow_potential.tif", Float32)[1]
c1_verify = Omniscape.read_raster("output_verify/test1_32bit/normalized_cum_currmap.tif", Float32)[1]
@test Omniscape.arrays_equal(a1, a1_verify)
@test Omniscape.arrays_equal(b1, b1_verify)
@test Omniscape.arrays_equal(c1, c1_verify)

q, e = run_omniscape("input/config3.ini")
q_verify = Omniscape.read_raster("output_verify/test3/cum_currmap.asc", Float64)[1]
e_verify = Omniscape.read_raster("output_verify/test3/flow_potential.asc", Float64)[1]
@test Omniscape.arrays_equal(q, q_verify)
@test Omniscape.arrays_equal(e, e_verify)

d = run_omniscape("input/config2.ini")
d_verify = Omniscape.read_raster("output_verify/test2/cum_currmap.asc", Float64)[1]
@test Omniscape.arrays_equal(d, d_verify)
d_1 = run_omniscape("input/config2.ini")
d_2 = run_omniscape("input/config2.ini")

reclassed = run_omniscape("input/config_reclass.ini")
reclassed_verify = Omniscape.read_raster("output_verify/test_reclass/cum_currmap.asc", Float64)[1]
@test Omniscape.arrays_equal(reclassed, reclassed_verify)

@test typeof(f) == Array{Union{Float64, Missing},2}
@test typeof(g) == Array{Union{Float64, Missing},2}
@test typeof(h) == Array{Union{Float64, Missing},2}
@test typeof(a) == Array{Union{Float64, Missing},2}
@test typeof(b) == Array{Union{Float64, Missing},2}
@test typeof(c) == Array{Union{Float64, Missing},2}
@test typeof(a1) == Array{Union{Float32, Missing},2}
@test typeof(b1) == Array{Union{Float32, Missing},2}
@test typeof(c1) == Array{Union{Float32, Missing},2}
@test typeof(d) == Array{Union{Float64, Missing},2}
@test typeof(e) == Array{Union{Float64, Missing},2}
@test typeof(reclassed) == Array{Union{Float64, Missing},2}
@test Omniscape.arrays_equal(a, d) #parallel and serial produce same result

# Single and double produce similar results
@test Omniscape.arrays_equal(a, a1)
@test Omniscape.arrays_equal(b, b1)
@test Omniscape.arrays_equal(c, c1)

GC.gc()

rm("test1", recursive = true)
rm("test1_32bit", recursive = true)
rm("test2", recursive = true)
rm("test2_1", recursive = true)
rm("test2_2", recursive = true)
rm("test3", recursive = true)
rm("test4", recursive = true)
rm("test5", recursive = true)
rm("test6", recursive = true)
rm("test_reclass", recursive = true)

end


text:
module Generate

using ..Models: ClusterParameters, DataParameters, GaussianParameters, GeometricArrivals, NtlParameters
using ..Models: FiniteTopicModelParameters
using ..Models: GaussianWishartParameters, PoissonArrivals
using ..Models: ArrivalDistribution, MultinomialParameters
using ..Models: Model, Mixture, Changepoint, HiddenMarkovModel
using Distributions
using LinearAlgebra
using DataFrames

struct LocationScale
    location::Vector{Float64}
    scale::Matrix{Float64}
end

function generate_log_pmfs(psi_variates::Vector{Float64})
    num_clusters = length(psi_variates)
    complement_psi_variates = 1 .- psi_variates
    log_psi_variates = log.(psi_variates)
    log_comple_psi_variates = log.(complement_psi_variates)

    log_pmfs = -ones(num_clusters, num_clusters)./0
    log_pmfs_diag = diagind(log_pmfs)
    log_pmfs[log_pmfs_diag] = log_psi_variates

    for j = 1:num_clusters
        log_pmfs[(j+1):num_clusters, j] .= log_psi_variates[j]
    end

    for i = 1:num_clusters
        log_pmfs[i:num_clusters, 1:(i-1)] .+= log_comple_psi_variates[i]
    end

    return log_pmfs
end

function generate_log_pmfs(psi_variates::Matrix{Float64})
    num_states = size(psi_variates)[1]
    log_transition_matrices = Array{Float64}(undef, num_states, num_states, num_states)
    for state = 1:num_states
        log_transition_matrices[state, :, :] = generate_log_pmfs(vec(psi_variates[state, :]))
    end
    return log_transition_matrices
end

function sample_clusters(arrivals::Array{Int64}, model::Mixture)
    cluster_parameters = model.cluster_parameters
    n = length(arrivals)
    cluster_assignments = Array{Int64}(undef, n)
    num_clusters = sum(arrivals)
    new_cluster_indices = arrivals .=== 1
    arrival_times = (1:n)[new_cluster_indices]
    psi_prior = cluster_parameters.prior
    psi_variates  = rand(Beta(psi_prior[1], psi_prior[2]), num_clusters)
    psi_variates[1] = 1

    log_pmfs = generate_log_pmfs(psi_variates)
    pmfs = exp.(log_pmfs)

    cluster_assignments[new_cluster_indices] = 1:num_clusters
    for i = 1:num_clusters
        arrival_time = arrival_times[i]
        if i < num_clusters
            next_arrival_time = arrival_times[i+1]
        else
            next_arrival_time = n + 1
        end
        observation_range = (arrival_time + 1):(next_arrival_time - 1)
        num_obs = length(observation_range)
        pmf = vec(pmfs[i, :])
        observation_assignments = rand(Categorical(pmf), num_obs)
        cluster_assignments[observation_range] = observation_assignments
    end

    return cluster_assignments
end

function sample_clusters(arrivals::Array{Int64}, ::Changepoint)
    n = length(arrivals)
    changepoint_assignments = Array{Int64}(undef, n)
    num_changepoints = sum(arrivals)
    changepoint_indices = arrivals .=== 1
    arrival_times = (1:n)[changepoint_indices]

    for index = 1:num_changepoints
        arrival_time = arrival_times[index]
        if index === num_changepoints
            next_arrival_time = n + 1
        else
            next_arrival_time = arrival_times[index+1]
        end
        changepoint_assignments[arrival_time:(next_arrival_time-1)] .= index
    end

    return changepoint_assignments
end

function sample_clusters(arrivals::Array{Int64}, model::HiddenMarkovModel)
    cluster_parameters = model.cluster_parameters
    n = length(arrivals)
    cluster_assignments = Array{Int64}(undef, n)
    num_clusters = sum(arrivals)
    new_cluster_indices = arrivals .=== 1
    arrival_times = (1:n)[new_cluster_indices]
    psi_prior = cluster_parameters.prior

    psi_variates = rand(Beta(psi_prior[1], psi_prior[2]), num_clusters*num_clusters)
    psi_variates_matrix = reshape(psi_variates, num_clusters, num_clusters)
    psi_variates_matrix[:, 1] .= 1

    log_pmfs = generate_log_pmfs(psi_variates_matrix)
    pmfs = exp.(log_pmfs)

    cluster_assignments[new_cluster_indices] = 1:num_clusters
    for observation = 2:n
        if !(observation in arrival_times)
            previous_assignment = cluster_assignments[observation-1]
            newest_state = maximum(cluster_assignments[1:observation-1])
            pmf = vec(pmfs[previous_assignment, newest_state, :])
            assignment = rand(Categorical(pmf), 1)[1]
            cluster_assignments[observation] = assignment
        end
    end
    return cluster_assignments
end

function sample_cluster_parameters(num_clusters, data_parameters::GaussianParameters)
    dist = MvNormal(data_parameters.prior_mean, data_parameters.prior_covariance)
    return rand(dist, num_clusters)
end

function sample_cluster_parameters(num_clusters, data_parameters::GaussianWishartParameters)
    location_scale_parameters = Vector{LocationScale}(undef, num_clusters)
    dim = data_parameters.dim
    for i = 1:num_clusters
        covariance_matrix = rand(InverseWishart(data_parameters.dof, data_parameters.scale_matrix), 1)[1]
        covariance_matrix = reshape(covariance_matrix, dim, dim)
        scale = data_parameters.scale
        data_mean = vec(rand(MvNormal(data_parameters.prior_mean, (1/scale)covariance_matrix), 1))
        location_scale = LocationScale(data_mean, covariance_matrix)
        location_scale_parameters[i] = location_scale
    end
    return location_scale_parameters
end

function sample_cluster_parameters(num_clusters, data_parameters::MultinomialParameters)
    dist = Dirichlet(data_parameters.prior_dirichlet_scale)
    return rand(dist, num_clusters)
end

function sample_cluster_parameters(num_clusters, data_parameters::FiniteTopicModelParameters)
    # Generates the word distribution for each topic
    num_topics = data_parameters.num_topics
    num_words = data_parameters.num_words
    word_parameter = data_parameters.word_parameter
    topic_word_distributions = Matrix{Float64}(undef, num_words, num_topics)
    for topic = 1:num_topics
        word_prior = Dirichlet(num_words, word_parameter)
        word_distribution_parameter = rand(word_prior, 1)
        topic_word_distributions[:, topic] = word_distribution_parameter
    end
    return topic_word_distributions
end

function sample_data(cluster, num_assigned, topic_word_distributions, data_parameters::FiniteTopicModelParameters)
    # For a particular cluster, samples topics from the topic prior for each word, and then samples word token from
    # topic distribution
    num_topics = data_parameters.num_topics
    topic_parameter = data_parameters.topic_parameter
    topic_prior = Dirichlet(num_topics, topic_parameter)
    topic_distribution_parameter = vec(rand(topic_prior, 1))
    num_words = data_parameters.num_words
    topic_distribution = Multinomial(data_parameters.length, topic_distribution_parameter)
    topic_frequencies = rand(topic_distribution, num_assigned)

    words = zeros(Int64, num_words, num_assigned)
    for observation = 1:num_assigned
        for topic = 1:num_topics
            num_words_assigned_to_topic = topic_frequencies[topic, observation]
            word_distribution_parameter = vec(topic_word_distributions[:, topic])
            word_distribution = Multinomial(num_words_assigned_to_topic, word_distribution_parameter)
            words_from_topic = vec(rand(word_distribution, 1))
            words[:, observation] += words_from_topic
        end
    end
    return words
end

function sample_data(cluster, num_assigned, location_scale_parameters::Vector{LocationScale},
                     ::GaussianWishartParameters)
    location_scale = location_scale_parameters[cluster]
    location = location_scale.location
    scale = location_scale.scale
    normal_distribution = MvNormal(location, scale)
    cluster_data = rand(normal_distribution, num_assigned)
    return cluster_data
end

function sample_data(cluster::Int64, num_assigned::Int64, cluster_means::Matrix{Float64},
                     data_parameters::GaussianParameters)
    cluster_mean = vec(cluster_means[:, cluster])
    normal_distribution = MvNormal(cluster_mean, data_parameters.data_covariance)
    cluster_data = rand(normal_distribution, num_assigned)
    return cluster_data
end

function sample_data(cluster::Int64, num_assigned::Int64, cluster_probabilities::Matrix{Float64},
                     data_parameters::MultinomialParameters)
    cluster_probability = vec(cluster_probabilities[:, cluster])
    multinomial_distribution = Multinomial(data_parameters.n, cluster_probability)
    cluster_data = rand(multinomial_distribution, num_assigned)
    return cluster_data
end

function generate_arrival_times(n::Int64, arrival_distribution::GeometricArrivals)
    arrival_prior = arrival_distribution.prior
    phi = reshape(rand(Beta(arrival_prior[1], arrival_prior[2]), 1), 1)[1]
    arrivals = rand(Binomial(1, phi), n-1)
    prepend!(arrivals, [1])
    return arrivals
end

function generate_arrival_times(n::Int64, arrival_distribution::PoissonArrivals)
    alpha = arrival_distribution.alpha
    beta = arrival_distribution.beta
    rate = rand(Gamma(alpha, 1/beta), 1)[1]
    arrivals = zeros(Int64, n)
    arrival_time = 1
    while arrival_time <= n
        arrivals[arrival_time] = 1
        arrival_time += rand(Poisson(rate), 1)[1] + 1
    end
    return arrivals
end

function prepare_data_matrix(data_parameters::MultinomialParameters, n)
    dim = data_parameters.dim
    return Array{Int64}(undef, dim, n)
end

function prepare_data_matrix(data_parameters::GaussianParameters, n)
    dim = data_parameters.dim
    return Array{Float64}(undef, dim, n)
end

function prepare_data_matrix(data_parameters::GaussianWishartParameters, n)
    dim = data_parameters.dim
    return Array{Float64}(undef, dim, n)
end

function prepare_data_matrix(data_parameters::FiniteTopicModelParameters, n)
    dim = data_parameters.num_words
    return Array{Int64}(undef, dim, n)
end

get_cluster_parameters(model::Union{Mixture, HiddenMarkovModel}) = model.cluster_parameters
get_cluster_parameters(model::Changepoint) = model.cluster_parameters

function generate(model::Model; n::Int64=100)
    data_parameters = model.data_parameters

    cluster_prior_parameters = get_cluster_parameters(model)
    arrivals = generate_arrival_times(n, cluster_prior_parameters.arrival_distribution)
    num_clusters = sum(arrivals)
    assignments = sample_clusters(arrivals, model)
    data = prepare_data_matrix(data_parameters, n)

    cluster_parameters = sample_cluster_parameters(num_clusters, data_parameters)

    for cluster = 1:num_clusters
        assigned_to_cluster = (assignments .=== cluster)
        num_assigned = count(assigned_to_cluster)
        cluster_data = sample_data(cluster, num_assigned, cluster_parameters, data_parameters)
        data[:, assigned_to_cluster] = cluster_data
    end

    assignments = DataFrame(reshape(assignments, length(assignments), 1), :auto)
    assignments = select(assignments, "x1" => "cluster")
    data = DataFrame(transpose(data), :auto)

    data = hcat(assignments, data)
    return data
end

end

max_stars_repo_name: realseanla/ntl-stick-breaking-julia

max_issues_repo_name: realseanla/ntl-stick-breaking-julia

max_forks_repo_name: realseanla/ntl-stick-breaking-julia


text:
# Uncomment the following
# line to run this example
# with more julia workers:
addprocs(4)

@everywhere begin
    using StochasticSearch

    function tour_cost(x::Configuration, parameters::Dict{Symbol, Any})
        result = float(readall(`./tour_cost $(x["Tour"].value)`))
        result
    end

end

tour = ["1"]
for i = 2:48
    push!(tour, string(i))
end
shuffle!(tour)

configuration = Configuration([PermutationParameter(tour ,"Tour")],
                               "TSP Solution")

tuning_run = Run(cost                = tour_cost,
                 starting_point      = configuration,
                 methods             = [[iterated_local_search 4];],
                 stopping_criterion  = elapsed_time_criterion,
                 duration            = 600,
                 reporting_criterion = elapsed_time_reporting_criterion,
                 report_after        = 30)

search_task = @task optimize(tuning_run)

result = consume(search_task)
print("$(result.current_time) ")
println(result.cost_minimum)
while result.is_final == false
    result = consume(search_task)
    print("$(result.current_time) ")
    println(result.cost_minimum)
end


text:
abstract type Constraints end

"""
    constraints template
"""
function constraints!(c, Z, con::Constraints, model, idx, h, T)
    return nothing
end

function constraints_jacobian!(∇c, Z, con::Constraints, model, idx, h, T)
   return nothing
end

function constraints_sparsity(con::Constraints, model, idx, T;
    shift_row = 0, shift_col = 0)
    row = []
    col = []
    return collect(zip(row, col))
end

function constraints_inequality(con::Constraints)
    return (1:0)
end

"""
    empty constraints
"""

struct EmptyConstraints <: Constraints
    n
    ineq
end

function EmptyConstraints()
    EmptyConstraints(0, (1:0))
end

"""
    multiple constraints
"""
struct MultiConstraints <: Constraints
    con::Vector{Constraints}
    n
    ineq
end

function multiple_constraints(con)
    n = sum([_con.n for _con in con])
    ineq = vcat([i == 1 ? _con.ineq : sum([con[k].n for k = 1:i-1]) .+ _con.ineq for (i,_con) in enumerate(con)]...)

    return MultiConstraints(con, n, ineq)
end

function add_constraints(c1::Constraints, c2::Constraints)
    multiple_constraints([c1, c2])
end

function add_constraints(c1::MultiConstraints, c2::Constraints)
    multiple_constraints([c1.con..., c2])
end

function constraints!(c, Z, con::MultiConstraints, model, idx, h, T)
    shift = 0
    for _con in con.con
        c_idx = shift .+ (1:_con.n)
        constraints!(view(c, c_idx), Z, _con, model, idx, h, T)
        shift += _con.n
    end
    return nothing
end

function constraints_jacobian!(∇c, Z, con::MultiConstraints, model, idx, h, T)
    con_shift = 0
    jac_shift = 0
    for _con in con.con
        jac_len = length(constraints_sparsity(_con, model, idx, T, shift_row = con_shift))
        jac_idx = jac_shift .+ (1:jac_len)
        constraints_jacobian!(view(∇c, jac_idx), Z, _con, model, idx, h, T)
        con_shift += _con.n
        jac_shift += jac_len
    end
    return nothing
end

function constraints_sparsity(con::MultiConstraints, model, idx, T;
    shift_row = 0, shift_col = 0)

    sparsity = []
    for _con in con.con
        spar = constraints_sparsity(_con, model, idx, T,
            shift_row = shift_row, shift_col = shift_col)
        sparsity = collect([sparsity..., spar...])
        shift_row += _con.n
    end
    return sparsity
end

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
using Oceananigans.Grids: AbstractGrid
using Oceananigans.Operators: ∂xᶠᵃᵃ, ∂yᵃᶠᵃ
using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions
using Oceananigans.Fields

using Adapt

struct ImplicitFreeSurface{E, G, B, VF, I}
    η :: E
    gravitational_acceleration :: G
    barotropic_volume_flux :: B
    vertically_integrated_lateral_face_areas :: VF
    implicit_step_solver :: I
end

# User interface to ImplicitFreeSurface
ImplicitFreeSurface(; gravitational_acceleration=g_Earth) =
    ImplicitFreeSurface(nothing, gravitational_acceleration, nothing, nothing, nothing)

### ExplicitFreeSurface(; gravitational_acceleration=g_Earth) =
###     ExplicitFreeSurface(nothing, gravitational_acceleration)

### Adapt.adapt_structure(to, free_surface::ExplicitFreeSurface) =
###     ExplicitFreeSurface(Adapt.adapt(to, free_surface.η), free_surface.gravitational_acceleration)

# Internal function for HydrostaticFreeSurfaceModel
function FreeSurface(free_surface::ImplicitFreeSurface{Nothing}, velocities, arch, grid)
    η = CenterField(arch, grid, TracerBoundaryConditions(grid))
    g = convert(eltype(grid), free_surface.gravitational_acceleration)

    barotropic_x_volume_flux = ReducedField(Face, Center, Nothing, arch, grid; dims=(3), boundary_conditions=nothing)
    barotropic_y_volume_flux = ReducedField(Center, Face, Nothing, arch, grid; dims=(3), boundary_conditions=nothing)
    barotropic_volume_flux = (u=barotropic_x_volume_flux, v=barotropic_y_volume_flux)

    Ax_zintegral = ReducedField(Face, Center, Nothing, arch, grid; dims=(3), boundary_conditions=nothing)
    Ay_zintegral = ReducedField(Center, Face, Nothing, arch, grid; dims=(3), boundary_conditions=nothing)
    vertically_integrated_lateral_face_areas = (Ax = Ax_zintegral, Ay=Ay_zintegral)
    compute_vertically_integrated_lateral_face_areas!(vertically_integrated_lateral_face_areas, grid, arch)

    ## implicit_step_solver = ImplicitFreeSurfaceSolver(arch, η, vertically_integrated_lateral_face_areas; maxit=100)
    implicit_step_solver = ImplicitFreeSurfaceSolver(arch, η, vertically_integrated_lateral_face_areas)

    return ImplicitFreeSurface(η, g,
                               barotropic_volume_flux,
                               vertically_integrated_lateral_face_areas,
                               implicit_step_solver)
end

### In final form the two functions below will return 0 (these are invoked in hydrostatic_free_surface_tendency_kernel_functions.jl )
explicit_barotropic_pressure_x_gradient(i, j, k, grid, free_surface::ImplicitFreeSurface) = 0

explicit_barotropic_pressure_y_gradient(i, j, k, grid, free_surface::ImplicitFreeSurface) = 0

@kernel function _compute_vertically_integrated_lateral_face_areas!(grid, A )
    i, j, k = @index(Global, NTuple)
    # U.w[i, j, 1] = 0 is enforced via halo regions.
    A.Ax[i, j, 1] = 0
    A.Ay[i, j, 1] = 0
    @unroll for k in 1:grid.Nz
        @inbounds A.Ax[i, j, 1] += Δyᶠᶜᵃ(i, j, k, grid)*Δzᵃᵃᶜ(i, j, k, grid)
        @inbounds A.Ay[i, j, 1] += Δxᶜᶠᵃ(i, j, k, grid)*Δzᵃᵃᶜ(i, j, k, grid)
    end
end

function compute_vertically_integrated_lateral_face_areas!(vertically_integrated_lateral_face_areas, grid, arch)

    event = launch!(arch,
                    grid,
                    :xyz,
                    _compute_vertically_integrated_lateral_face_areas!,
                    grid,
                    vertically_integrated_lateral_face_areas,
                    dependencies=Event(device(arch)))

    wait(device(arch), event)

    return
end


text:
ENV["GKSwstype"]="100"

using  Literate
using  Plots
import Remark
using Remark, FileWatching

# files =  filter( f -> startswith(f, "0"), readdir("src")) |> collect

files = [ "01.Introduction.jl",
          "02.RungeKuttaMethods.jl",
          "03.PoissonEquation.jl",
          "04.HOODESolver.jl",
          "05.Conclusion.jl"]

run(pipeline(`cat src/$files`; stdout="slides/src/index.jl" ))

slideshowdir = Remark.slideshow("slides",
                                options = Dict("ratio" => "16:9"),
                                title = "Differential equations with Julia")

# Open presentation in default browser.
# Remark.open(slideshowdir)

Literate.notebook("slides/src/index.jl", execute=false)

cp("slides/build", "docs", force=true)


text:
using ChemicalReactionNetworks

###
chems = ["CH₄", "O₂", "CO₂", "H₂O"]
rn = ReactionNetwork([Reaction([1,2],[1,2],[3,4],[1,2],1.0,0.0)], chems)
print(rn)

###

reactions = [Reaction([1,2],[1,1],[3],[1],1.0,1.0),
            Reaction([1,1],[1,1],[2],[2],1.0,1.0)]

cl = conservation_laws(ReactionNetwork(reactions))

@assert convert(Array,cl) == [1 1 2]

reactions = [Reaction([1,2],[1,1],[3],[1],1.0,1.0),
            Reaction([1,1],[1,1],[2],[1],1.0,1.0)]

cl = conservation_laws(ReactionNetwork(reactions))

@assert convert(Array,cl) == [1 2 3]

# Figure 1 in Rao paper

reactions = [Reaction([1],[1],[2,3],[2,1],1.0,1.0),
            Reaction([3,4],[1,1],[5],[1],1.0,1.0)]

cmm=conservation_laws(ReactionNetwork(reactions))

@assert convert(Array,cmm) ==
[ 2  1  0  0  0
 0  0  0  1  1
 1  0  1  0  1]

# Example 5 in Rao paper, without chemostatting
# species: [Y1, Xa, Xb, Xc, Y2]
reactions = [Reaction([1,2],[1,1],[3],[1],1.0,1.0),
            Reaction([3],[1],[2,5],[1,1],1.0,1.0),
            Reaction([2,5],[1,1],[4],[1],1.0,1.0),
            Reaction([4],[1],[1,2],[1,1],1.0,1.0)]

clm = cycles(ReactionNetwork(reactions))

@assert vec(convert(Array,clm)) == [1,1,1,1]

# Example 5 in Rao paper, with chemostatting.
# A chemostatted species can be represented as a 'reaction',
# with no reactants, and the species as product.
# The equilibrium concentration is k+/k-,
# And the speed at which it reaches equilibrium is proportional to k+
reactions = [Reaction([1,2],[1,1],[3],[1],1.0,1.0),
            Reaction([3],[1],[2,5],[1,1],1.0,1.0),
            Reaction([2,5],[1,1],[4],[1],1.0,1.0),
            Reaction([4],[1],[1,2],[1,1],1.0,1.0),
            Reaction([],[],[1],[1],100.0,100.0),
            Reaction([],[],[5],[1],100.0,100.0)]

clm = cycles(ReactionNetwork(reactions))

# Note that Rao's conclusion of 'one linearly independent emergent cycle (20) arises'
# is actually wrong; there are two emergent cycles, but they can be combined into one.
@assert convert(Array,clm[1:4,:]) ==
[ 0   1
  0   1
  1   0
  1   0]

# There is also only one conservation law in the chemostatted system
clm = conservation_laws(ReactionNetwork(reactions))
@assert vec(convert(Array,clm)) == [0,1,1,1,0]

#-----
# Test edge cases

# Case where one species does not react
reactions = [Reaction([1],[1],[2],[1],1.0,1.0),
             Reaction([2],[1],[4],[1],1.0,1.0)]
clm = conservation_laws(ReactionNetwork(reactions))
@assert convert(Array,clm) ==
[ 1  1  0  1
  0  0  1  0]

#-----
# Test complex decomposition
cm,imm=complex_decomposition(ReactionNetwork(reactions))
dp,dr=stoichiometric_matrix(ReactionNetwork(reactions))
@assert cm*imm == dr-dp

# Test (29) in paper
# [Xa,Xb,Xc]
reactions = [Reaction([1],[1],[2],[1],1.0,1.0),
            Reaction([1,2],[1,1],[2],[2],1.0,1.0),
            Reaction([2],[2],[3],[1],1.0,1.0)]
cm,imm=complex_decomposition(ReactionNetwork(reactions))
dp,dr=stoichiometric_matrix(ReactionNetwork(reactions))
@assert cm ==
[ 1  0  1  0  0
  0  1  1  2  0
  0  0  0  0  1
]
@assert imm ==
[ -1   0   0
   1   0   0
   0  -1   0
   0   1  -1
   0   0   1
]
@assert cm*imm == dr-dp

# Example 9 in paper.
# Here we 'chemostat' species by simply removing them from the equations;
# incorporating the chemostatted concentrations into the rate equations.
# This procedure is described:
# "This regrouping corresponds to the equivalent CRN only made of
#  internal species with effective rate constant ruling each reaction"

# Here the chemostatted species are 1 and 3.
# The only species left is 2, which we re-label as 1.
reactions = [Reaction([ ],[ ],[1],[1],1.0,1.0),
            Reaction([1],[1],[1],[2],1.0,1.0),
            Reaction([1],[2],[ ],[ ], 1.0,1.0)]
# we can also get the above by doing:
# chemostatted!(reactions, [1,3], [1.0, 1.0])
cm,imm=complex_decomposition(ReactionNetwork(reactions))
@assert cm ==  [0  1  2]
@assert imm ==
[ -1   0   1
   1  -1   0
   0   1  -1]

# ---------------
# Test chemostatting, with previous example.
####
reactions = [Reaction([1],[1],[2],[1],1.0,1.0),
            Reaction([1,2],[1,1],[2],[2],1.0,1.0),
            Reaction([2],[2],[3],[1],1.0,1.0)]
chemostatted!(reactions, [1,3],[1.0,1.0])
relabel!(reactions, [0,1,0])
cm,imm=complex_decomposition(ReactionNetwork(reactions))
@assert cm ==  [0  1  2]
@assert imm ==
[ -1   0   1
   1  -1   0
   0   1  -1]

#----------------
# [Ya,Xb,Xc,Xd,Ye]
reactions = [Reaction([1],[1],[2],[1],1.0,1.0),
            Reaction([2],[1],[3,4],[1,1],1.0,1.0),
            Reaction([3,4],[1,1],[5],[1],1.0,1.0)]

@assert deficiency(ReactionNetwork(reactions)) == 0

chemostatted!(reactions, [1,5], [1.0, 1.0])

relabel!(reactions, [0,1,2,3,0])

@assert deficiency(ReactionNetwork(reactions)) == 0

#------------------
# Futile cycle enzyme example.
# The deficiency is 1.
n_species=4
reactions = [Reaction([1,2],[1,1],[3,2],[1,1],1.0,1.0),
            Reaction([3,4],[1,1],[1,4],[1,1],1.0,1.0)]

@assert deficiency(ReactionNetwork(reactions)) == 1

#-------------------
# Examples given on:
# https://reaction-networks.net/wiki/Reaction_graph#Deficiency

reactions = [Reaction([],[],[1],[1],1.0,0.0),
            Reaction([1],[1],[2],[1],1.0,0.0),
            Reaction([2],[1],[],[],1.0,0.0)]

@assert deficiency(ReactionNetwork(reactions)) == 0

reactions = [Reaction([1],[2],[2],[2],1.0,0.0),
            Reaction([2],[2],[1,2],[1,1],1.0,0.0),
            Reaction([1,2],[1,1],[1],[2],1.0,0.0)]
@assert deficiency(ReactionNetwork(reactions)) == 1

reactions = [Reaction([1,2],[2,1],[1],[3],1.0,0.0),
            Reaction([1],[3],[1,2],[1,2],1.0,0.0),
            Reaction([1,2],[1,2],[2],[3],1.0,0.0),
            Reaction([2],[3],[1,2],[2,1],1.0,0.0)]
@assert deficiency(ReactionNetwork(reactions)) == 2

# ----------

# example taken from Sarkar paper
reactions = [Reaction([1,1],[1,1],[2],[1],1.0,1.0),
            Reaction([1,2],[1,1],[3],[1],1.0,1.0),
            Reaction([1,3],[1,1],[4],[1],1.0,1.0),
            Reaction([1,3],[1,1],[2],[2],2.0,1.0),
            Reaction([1,4],[1,1],[2,3],[1,1],1.0,1.0),
            Reaction([2,4],[1,1],[3],[2],1.0,1.0)]
chem_potentials=ones(4)
fe=free_energies(chem_potentials, reactions)
kf,kr = reaction_rates(fe, ones(6), 1.0)
for i=1:6
    reactions[i].kf = kf[i]
    reactions[i].kr = kr[i]
end
z=equilibrium_state(ReactionNetwork(reactions))
dz=mass_action(reactions, z)
@assert all(dz .< 1e-9)   # the network should admit an equilibrium state

# ------------------------------------
# Testing mass action: the Oregonator

chems = ["HBrO₂", "Br⁻", "Ce(IV)", "BrO₃⁻", "CH₂(COOH)₂", "HOBr"]
#chems = ["X", "Y", "Z", "A", "B", "P"]
k1 = 1.0
k2 = 1.0
k3 = 1.0
k4 = 1.0
k5 = 1.0

reactions = [Reaction([4,2],[1,1],[1,6],[1,1],k1,0),
            Reaction([1,2],[1,1],[6],[2],k2,0),
            Reaction([4,1],[1,1],[1,3],[2,2],k3,0),
            Reaction([1],[2],[4,6],[1,1],k4,0),
            Reaction([5,3],[2,2],[2],[1],k5,0)]

@show reactions

# -----------------------------------
# Autocatalysis reaction
chems = ["A₂","B","AB","Z","AZ","BZ"]
n_species=length(chems)
reactions = [Reaction([1,4],[1,2],[5],[2],1.0,1.0),
            Reaction([2,4],[1,1],[6],[1],1.0,1.0),
            Reaction([5,6],[1,1],[3,4],[1,2],1.0,1e-4)]
# Chemostat A₂, B, and AB
chemostatted!(reactions, [1,2,3], [1.0,1.0,1.0])
f(z,p,t) = mass_action(reactions,z)
#=function eom!(znew,z,p,t)
    znew[:] = mass_action(reactions,z)
end
ds=ContinuousDynamicalSystem(eom!, ones(n_species), nothing)
=#
function eom(z,p,t)
    dz = mass_action(reactions,z)
    for i=1:length(z)
        if z[i] < 0
            dz[i] = exp(-z[i])
        end
    end
    return dz
end
# ssz=solve(SteadyStateProblem(f, equilibrium_state(6, reactions)))

# ---------------------------------
# Test specificity

reactions = [Reaction([1],[1],[2],[1],1.0,1.0),
             Reaction([1],[1],[4],[1],1.0,1.0)]
spmat = [0.5 1.0
         0.5 1.0]
isapprox(specificity(ReactionNetwork(reactions), rand(n_species)), spmat)

reactions = [Reaction([1],[1],[2],[1],1.0,1.0),
             Reaction([1],[1],[4],[1],0.0,1.0)]
spmat = [1.0 1.0
         0.0 1.0]
isapprox(specificity(ReactionNetwork(reactions), rand(n_species)), spmat)

# ----------------------------------
# Test normal form
reactions = [Reaction([1,1],[1,2],[3],[1],1.0,1.0),
             Reaction([1],[2],[4],[1],1.0,1.0)]

rnf = normal_form(ReactionNetwork(reactions))

@assert Reaction([1],[3],[3],[1],1.0,1.0) == rnf.reactions[1]


text:
"""
"""
function fit!(
        transformer::ImmutablePredictProbaSingleLabelInt2StringTransformer,
        varargs...;
        kwargs...
        )
    if length(varargs) == 1
        return varargs[1]
    else
        return varargs
    end
end

"""
"""
function predict(
        transformer::ImmutablePredictProbaSingleLabelInt2StringTransformer,
        varargs...;
        kwargs...
        )
    if length(varargs) == 1
        return varargs[1]
    else
        return varargs
    end
end

"""
"""
function predict_proba(
        transformer::ImmutablePredictProbaSingleLabelInt2StringTransformer,
        single_labelprobabilities::AbstractDict;
        kwargs...
        )
    labelint2stringmap = getlabelint2stringmap(
        transformer.levels,
        transformer.index,
        )
    result = Dict()
    for key in keys(single_labelprobabilities)
        result[labelint2stringmap[key]] = single_labelprobabilities[key]
    end
    result = fix_type(result)
    return result
end


text:
"Used to compile to a function to machine code."
struct Native <: AbstractMachine
end

compile(::Type{Native}, program::Program) = func(program)

# Based on Mike's IRTools.jl
function func(ir::Program)
    @eval @generated function $(gensym())($([v.symbol for v in input_variables(ir)]...))
        return build_function_body_expr($ir)
    end
end

tmpsym(i::Integer) = Symbol(:t, i)

const BB = Union{BasicBlock, CompiledBasicBlock}

function build_function_body_expr(ir::BB)
    block = Expr(:block)

    itr = eachinstruction(ir)
    for (i, instruction) in enumerate(itr)

        expr = julia_expr(ir, instruction)

        if expr !== nothing
            push!(block.args, Expr(:(=), tmpsym(i), expr))
        end
    end

    push!(block.args, return_expr(ir))

    return block
end

function return_expr(bb::BB)
    out_vars = output_variables(bb)

    if isempty(out_vars)

        # function does not return values
        return :(return nothing)

    elseif length(out_vars) == 1

        # function returns a single value
        return :(return $(julia_expr(bb, out_vars[1])))

    else

        # function returns a tuple
        ret_tuple = Expr(:tuple)

        for v in out_vars
            push!(ret_tuple.args, julia_expr(bb, v))
        end

        return Expr(:return, ret_tuple)
    end
end

julia_expr(bb::BB, c::CallUnary) = Expr(:call, c.op, julia_expr(bb, c.arg))

function julia_expr(bb::BB, c::CallBinary)
    Expr(:call, c.op, julia_expr(bb, c.arg1), julia_expr(bb, c.arg2))
end

function julia_expr(bb::BB, c::Call3Args)
    Expr(:call, c.op, julia_expr(bb, c.arg1), julia_expr(bb, c.arg2), julia_expr(bb, c.arg3))
end

function julia_expr(bb::BB, c::Call4Args)
    Expr(:call, c.op, julia_expr(bb, c.arg1), julia_expr(bb, c.arg2), julia_expr(bb, c.arg3), julia_expr(bb, c.arg4))
end

function julia_expr(bb::BB, c::Call5Args)
    Expr(:call, c.op, julia_expr(bb, c.arg1), julia_expr(bb, c.arg2), julia_expr(bb, c.arg3), julia_expr(bb, c.arg4), julia_expr(bb, c.arg5))
end

function julia_expr(bb::BB, c::Call6Args)
    Expr(:call, c.op, julia_expr(bb, c.arg1), julia_expr(bb, c.arg2), julia_expr(bb, c.arg3), julia_expr(bb, c.arg4), julia_expr(bb, c.arg5), julia_expr(bb, c.arg6))
end

julia_expr(bb::BB, c::ImpureInstruction) = julia_expr(bb, c.call)
julia_expr(bb::BB, c::PureInstruction) = julia_expr(bb, c.call)

julia_expr(bb::BB, constant::Const) = constant.val
julia_expr(bb::BB, ssa::SSAValue) = tmpsym(ssa.address)

function julia_expr(bb::BB, ::Assignment{V}) where {V<:ImmutableVariable}
    # no-op
    return nothing
end

function julia_expr(bb::BB, op::Assignment{V}) where {V<:MutableVariable}
    Expr(:(=), op.lhs.symbol, julia_expr(bb, op.rhs))
end

julia_expr(bb::BB, variable::MutableVariable) = variable.symbol

function julia_expr(bb::BB, variable::ImmutableVariable)
    if is_input(bb, variable)
        return variable.symbol
    else
        return julia_expr(bb, bb.immutable_locals[variable])
    end
end


text:
#
# Level Sets
#
#    Functions for computing level curves and surfaces in simplicial complexes
#
#    Author: Carlos Valero Valdés
#    e-mail: valerocar@gmail.com
#
#

"""
    levelcurve(ps::Array{<:Real,2}, fs::Array{<:Int,2}, vals::::Array{<:Real,1}, level::Real)

Computes the "level curve" for a function defined as a collection of values on the vertices of a triangular mesh.
# Example

```julia-repl
julia> levelcurve(points,faces,values,level)
2×4 Array{Float64,2}:
 0.0  0.5  0.5  0.0
 0.5  1.0  1.0  0.5

julia> points = [[0 0];[1 0];[0 1];[1 1]]; faces = [[1 2 3];[2 4 3]]; values = [0,1,1,0]; level =1/2;

julia> levelcurve(points,faces,values,level)
2×4 Array{Float64,2}:
 0.0  0.5  0.5  0.0
 1.0  0.5  0.5  1.0
```

In the example above, each row of the matrix correspond to a level edge.
Columns 1 and 2  are the coordinates of the first vertex, and
columns 3 and 4 are the coordinates of the second vertex.
The edges are oriented so that the gradient of the function points to their left side.
"""
function levelcurve(ps::Array{<:Real,2}, fs::Array{<:Int,2}, vals::Array{<:Real,1}, level::Real)
    @assert size(fs,2) == 3
    dim = size(ps,2)
    law =[[[3 1],[3 2]], [[2 3],[2 1]], [[1 3],[1 2]],[[1 2],[1 3]],[[2 1],[2 3]],[[3 2],[3 1]]]
    cnt = 0
    es = zeros(size(fs,1),2*dim) # make enough space to hold edges (edges as pairs of points in 3-space)
    for i in 1:size(fs,1)
        fvs = fs[i,:] # Face vertices
        vals_fvs= vals[fvs] .> level
        ss = 1*vals_fvs[3] + 2*vals_fvs[2] + 4*vals_fvs[1]
        if 0 < ss < 7 # When intersection occurs
            ea, eb = law[ss]
            a1, a2 = fvs[ea]
            b1, b2 = fvs[eb]
            val_a1 = vals[a1]; val_a2 = vals[a2]
            val_b1 = vals[b1]; val_b2 = vals[b2]
            λa = (level-val_a1)/(val_a2-val_a1)
            λb = (level-val_b1)/(val_b2-val_b1)
            pa = (1-λa) *ps[a1,:] + λa*ps[a2,:]
            pb = (1-λb) *ps[b1,:] + λb*ps[b2,:]
            es[cnt+1,1:dim] = pa
            es[cnt+1,dim+1:2*dim] = pb
            cnt = cnt + 1
        end
    end
    es[1:cnt,:]
end

"""
    levelsurface(ps::Array{<:Real,2}, ts::Array{<:Int,2}, vals::Array{<:Real,1}, level::Real)

Computes the "level surface" for a function defined as a collection of values on the vertices of a tetrahedral mesh.
# Example

```julia-repl
julia> points = [[0 0 0];[1 0 0];[0 1 0];[0 0 1]]; tetrahedrons= [[1 2 3 4];]; values = [0,1,1,1]; level =1/2;

julia> levelsurface(points, tetrahedrons, val
valtype values
julia> levelsurface(points, tetrahedrons, values, level)
1×9 Array{Float64,2}:
 0.5  0.0  0.0  0.0  0.5  0.0  0.0  0.0  0.5
```

In the example above, each row of the matrix correspond to a level triangle.
Columns 1, 2 and 3 are the coordinates of the first vertex,
columns 4, 5 and 6 are the coordinates of the second vertex, and
columns 7, 8 and 9 are the coordinates of the third vertex.
The triangles are oriented so that the gradient of the function points to their left side.
"""
function levelsurface(ps::Array{<:Real,2}, ts::Array{<:Int,2}, vals::Array{<:Real,1}, level::Real)
    @assert size(ts,2) == 4
    dim = size(ps,2)
    # vertex1+
    l08 = [[[1 2],[1 4],[1 3]]]
    l07 = [[[1 2],[1 3],[1 4]]]
    # vertex2+
    l04 = [[[2 1],[2 3],[2 4]]]
    l11 = [[[2 1],[2 4],[2 3]]]
    # vertex3+
    l02 = [[[3 1],[3 4],[3 2]]]
    l13 = [[[3 1],[3 2],[3 4]]]
    # vertex4+
    l01 = [[[4 1],[4 2],[4 3]]]
    l14 = [[[4 1],[4 3],[4 2]]]

    #vertices??+
    l12 = [[[1 3],[2 3],[2 4]], [[1 3],[2 4],[1 4]]]
    l03 = [[[1 3],[2 4],[2 3]], [[1 3],[1 4],[2 4]]]
    #vertices??+
    l10 = [[[1 2],[1 4],[4 3]], [[1 2],[4 3],[2 3]]]
    l05 = [[[1 2],[4 3],[1 4]], [[1 2],[2 3],[4 3]]]
    #vertices??+
    l09 = [[[1 2],[2 4],[1 3]], [[1 3],[2 4],[4 3]]]
    l06 = [[[1 2],[1 3],[2 4]], [[1 3],[4 3],[2 4]]]

    ll = [l01,l02,l03,l04,l05,l06,l07,l08,l09,l10,l11,l12,l13,l14]
    cnt = 0
    fs = zeros(size(ts,1),3*dim) # make enough space to hold edges (edges as triplets of points in 3-space)

    for i in 1:size(ts,1)
        tvs = ts[i,:] # Tetrahedron vertices
        vals_tvs= vals[tvs] .> level
        c0, c1, c2, c3 = vals[tvs] .> level
        ss = 1*vals_tvs[4] + 2*vals_tvs[3] + 4*vals_tvs[2] + 8*vals_tvs[1]
        if 0 < ss < 15 # When intersection occurs
            for tr in ll[ss]
                ea, eb, ec = tr # edges having intersection points forming a level face
                a1, a2 = tvs[ea]
                b1, b2 = tvs[eb]
                c1, c2 = tvs[ec]
                val_a1 = vals[a1]; val_a2 = vals[a2]
                val_b1 = vals[b1]; val_b2 = vals[b2]
                val_c1 = vals[c1]; val_c2 = vals[c2]
                λa = (level-val_a1)/(val_a2-val_a1)
                λb = (level-val_b1)/(val_b2-val_b1)
                λc = (level-val_c1)/(val_c2-val_c1)
                pa = (1-λa) *ps[a1,:] + λa*ps[a2,:]
                pb = (1-λb) *ps[b1,:] + λb*ps[b2,:]
                pc = (1-λc) *ps[c1,:] + λc*ps[c2,:]
                fs[cnt+1,1:dim] = pa
                fs[cnt+1,dim+1:2*dim] = pb
                fs[cnt+1,2*dim+1:3*dim] = pc
                cnt = cnt + 1
            end
        end
    end
    fs[1:cnt,:]
end


text:
__precompile__()

module BedgraphFiles

using FileIO
using Requires

using Bedgraph

using IteratorInterfaceExtensions, TableTraits, TableTraitsUtils
using TableShowUtils

import IterableTables

function __init__()
    @require DataFrames="a93c6f00-e57d-5684-b7b6-d8193f3e46c0" include(joinpath(@__DIR__, "integrations","DataFrames.jl"))
end

const BedgraphFileFormat = File{format"bedGraph"}

struct BedgraphFile
    filename::String
    keywords
end

function Base.show(io::IO, source::BedgraphFile)
    TableShowUtils.printtable(io, getiterator(source), "bedGraph file")
end

function Base.read(file::BedgraphFile) :: Vector{Bedgraph.Record}
    # Read file using Bedgraph package.
    return open(file.filename, "r") do io
        Bedgraph.read(io, Vector{Bedgraph.Record})
    end

end

function load(f::BedgraphFileFormat; args...)
    return BedgraphFile(f.filename, args)
end

IteratorInterfaceExtensions.isiterable(x::BedgraphFile) = true
TableTraits.isiterabletable(x::BedgraphFile) = true
IteratorInterfaceExtensions.isiterable(x::Vector{Bedgraph.Record}) = true #Note: Vector{Bedgraph.Record} is iterable by default.
TableTraits.isiterabletable(x::Vector{Bedgraph.Record}) = true

function IteratorInterfaceExtensions.getiterator(records::Vector{Bedgraph.Record})

    columns = [
        Bedgraph.chrom.(records),
        Bedgraph.first.(records),
        Bedgraph.last.(records),
        Bedgraph.value.(records)
    ]

    names = Symbol[:chrom, :first, :last, :value]

    it = TableTraitsUtils.create_tableiterator(columns, names)

    return it
end

function IteratorInterfaceExtensions.getiterator(file::BedgraphFile)

    records = read(file) #TODO: Generate iterator from first record?

    it = getiterator(records)

    return it
end

function Base.collect(x::BedgraphFile)
    return collect(getiterator(x))
end

function _Records(x) :: Vector{Bedgraph.Record} #TODO: consider formalising Records function in bedgraph (e.g. Bedgraph.Records, Bedgraph.Bedgraph.Records) that returns Vector{Bedgraph.Record}.
    cols, names = create_columns_from_iterabletable(x, na_representation=:missing)

    return Bedgraph.Record.(cols[1], cols[2], cols[3], cols[4])
end

function Vector{Bedgraph.Record}(x::AbstractVector{T}) :: Vector{Bedgraph.Record} where {T<:NamedTuple}
    @debug "Vector{Bedgraph.Record}(x::AbstractVector{T})"
    return  _Records(x)
end

function Vector{Bedgraph.Record}(file::B) :: Vector{Bedgraph.Record} where {B<:BedgraphFile}
    @debug "Vector{Bedgraph.Record}(file::BedgraphFile)"
    return read(file)
end

function Vector{Bedgraph.Record}(x::T) :: Vector{Bedgraph.Record} where {T} #TODO: consider formalising Records function in bedgraph (e.g. Bedgraph.Records, Bedgraph.Bedgraph.Records) that returns Vector{Bedgraph.Record}.

    if TableTraits.isiterabletable(x)
        @debug "Vector{Bedgraph.Record}(x) - isiterabletable"
        return _Records(x)
    end

    @debug "Vector{Bedgraph.Record}(x) - converting"
    return convert(Vector{Bedgraph.Record}, x)

end

function save(file::BedgraphFileFormat, header::Bedgraph.BedgraphHeader, records::Vector{Bedgraph.Record}) :: Vector{Bedgraph.Record}

    write(file.filename, header, records)

    return records #Note: this return is useful when piping (e.g., records = some_operation | save(file)).
end

function save(file::BedgraphFileFormat, records::Vector{Bedgraph.Record}; bump_forward = true) :: Vector{Bedgraph.Record}

    sort!(records)

    header = Bedgraph.generateBasicHeader(records, bump_forward = bump_forward) #TODO: consolidate header generation and determine whether there is a need for bump_forward.

    return save(file, header, records)
end

function save(file::BedgraphFileFormat, data; bump_forward = true)

    it = getiterator(data)

    records = Vector{Bedgraph.Record}(it)

    save(file, records, bump_forward = bump_forward)

    return data #Note: this return is usful when piping.
end

end # module


text:
include("DM1_1.jl")
include("DM2_1.jl")
include("DM1_2.jl")
include("DM2_2.jl")
include("loadSPP.jl")
include("setSPP.jl")
include("getfname.jl")

target = "B:/Cours/Nantes/Metaheuristique/DM1_Metaheuristique/solveSPP-master/Data"
fnames = getfname(target)

rep = 10
alpha = 0.8
ite = 15
alphaset = 0
temps = 2
nbAlpha = 5

function UpgradeGRASP(cost, matrix, n, m, alpha)
    (SOL, z, desactive_condition) = GRASP(cost, matrix, n, m, alpha)
    (SOL, z) = exchange1_2(SOL,n,m,cost,desactive_condition,matrix)
    (SOL, z) = exchange1_1(SOL,n,m,cost,desactive_condition,matrix)
    return(SOL, z)
end

function UpgradeGlouton(cost, matrix, n, m, alpha)
    (SOL, z, desactive_condition) = Glouton(cost, matrix, n, m)
    (SOL, z) = exchange1_2(SOL,n,m,cost,desactive_condition,matrix)
    (SOL, z) = exchange1_1(SOL,n,m,cost,desactive_condition,matrix)
    return(SOL, z)
end

function expGRASP(fnames)
    for f in fnames
        println("====================================================================")
        println(f)
        cost, matrix, n, m = loadSPP(f)
        (SOLBest, zBest, desactive_conditionBestBest) = GRASP(cost, matrix, n, m, alpha)
        @time for i in 1:rep
#            (SOL, z, desactive_condition) = GRASP(cost, matrix, n, m, alpha)
#            (SOL, z) = exchange1_2(SOL,n,m,cost,desactive_condition,matrix)
#            (SOL, z) = exchange1_1(SOL,n,m,cost,desactive_condition,matrix)
            (SOL, z) = UpgradeGRASP(cost, matrix, n, m, alpha)
            if z > zBest
                zBest = z
                SOLBest = SOL
            end
        end
        println("Meilleur soluttion pour ", rep , " répétitions, et alpha = ", alpha, " :  ", zBest)
    end
    cd("../")
end

function expGlouton(fnames)
    for f in fnames
        println("====================================================================")
        println(f)
        cost, matrix, n, m = loadSPP(f)

        @time (SOL, z) = UpgradeGlouton(cost, matrix, n, m, alpha)

        println("Solution glouton après amélioration: ", z)
    end
    cd("../")
end

function expReact(fnames)
    for f in fnames
        println("====================================================================")
        println(f)
        cost, matrix, n, m = loadSPP(f)
        @time for i in 1:rep
            ReactiveGRASP(matrix,cost, n, m, ite, nbAlpha, alphaset, temps)
        end
    end
    cd("../")
end

#expReact(fnames)
#expGRASP(fnames)
#expGlouton(fnames)


text:
function render_response(response)
    if Tables.istable(response)
        render_df_response(response)
    else
        error("Don't know how to render result of type `$(typeof(response))`")
    end
end

render_response(response::AbstractString) = response
render_response(response::Bool) = response
render_response(response::Number) = response
render_response(response::Nothing) = nothing
render_response(response::DataFrames.AbstractDataFrame) = render_df_response(response)

function render_response(response::Plots.Plot)
    io = IOBuffer()
    show(io, "image/png", response)
    png = read(seekstart(io), String)

    "data:image/png;base64," * base64encode(png)
end

function jsonify_table_row(row)
    Dict(columnname => getproperty(row, columnname) for columnname in propertynames(row))
end

function render_df_response(response)
    data = [jsonify_table_row(row) for row in Tables.rows(response)]
    if length(data) == 0
         data = String[]
     end
    colnames = string.(Tables.schema(response).names)
    header = [Dict(:text => colname, :value => colname) for colname in colnames]

    [header, data...]
end

function render_response(response::AbstractArray)
    response
end

function handle_request(id, input_dict, server_module, session)
    fn = get_handler(id, server_module)
    names = argument_names(fn)
    args = ()
    for argname in names
        if argname == :session
            args = (args..., session)
        else
            try
                args = (args..., input_dict[string(argname)])
            catch e
                msg = sprint(showerror, e)
                if typeof(e) <: KeyError
                    msg = "Bad Server module configuration: Server-side function `$id` takes `$argname` as an input, but there is no UI element with that id"
                    @error msg
                end
                return JSON.json(Dict("matte_error_msg" => msg))
            end
        end
    end
    try
        response = fn(args...)
        response = render_response(response)
        JSON.json(Dict("id" => id, "value" => response))
    catch e
        msg = sprint(showerror, e)
        @error string(msg, "\nError occurred calling function `$fn` with arguments: `$(args...)`")

        JSON.json(Dict("matte_error_msg" => msg))
    end
end


text:
using Documenter, MCVI

makedocs(modules=MCVI)

deploydocs(
           repo   = "github.com/JuliaPOMDP/MCVI.jl.git",
           julia  = "0.5",
           osname = "linux"
           )


text:
#=

Benchmark simple models

=#

using DiscreteChoiceModels, CSV, DataFrames, BenchmarkTools, Logging, Dagger, Distributed, CodecZlib

const SAMPLES = 10

get_test_file_path(file) = joinpath(Base.source_dir(), "../test/data/", file)

function benchmark_mnl()
    data = CSV.read(get_test_file_path("biogeme_swissmetro.dat"), DataFrame, delim='\t')
    data = data[in.(data.PURPOSE, [Set([1, 3])]) .& (data.CHOICE .!= 0), :]

    data.avtr = (data.TRAIN_AV .== 1) .& (data.SP .!= 0)
    data.avsm = data.SM_AV .== 1
    data.avcar = (data.CAR_AV .== 1) .& (data.SP .!= 0)

    # TODO ensure compilation time for the macro is included!!!! for a fair comparison.
    @benchmarkable(
        multinomial_logit(
            @utility(begin
                1 ~ αtrain + βtravel_time * TRAIN_TT / 100 + βcost * (TRAIN_CO * (GA == 0)) / 100
                2 ~ αswissmetro + βtravel_time * SM_TT / 100 + βcost * SM_CO * (GA == 0) / 100
                3 ~ αcar + βtravel_time * CAR_TT / 100 + βcost * CAR_CO / 100

                αswissmetro = 0, fixed  # fix swissmetro ASC to zero
            end),
            :CHOICE,
            $data,
            availability=[
                1 => :avtr,
                2 => :avsm,
                3 => :avcar,
            ]
        ),
        samples=SAMPLES,
        seconds=1e6,
        evals=1
     )
end

function benchmark_mnl_dtable()
    data = DTable(Tables.rowtable(CSV.File(get_test_file_path("biogeme_swissmetro.dat"), delim='\t')), 2000)

    # filter to wanted cases
    data = filter(data) do row
        ((row.PURPOSE == 1) | (row.PURPOSE == 3)) && row.CHOICE != 0
    end

    data = map(data) do r
        merge((
            avtr=(r.TRAIN_AV == 1) && (r.SP != 0),
            avsm=r.SM_AV == 1,
            avcar=(r.CAR_AV == 1) && (r.SP != 0)
        ), pairs(r))
    end

    @benchmarkable(
        multinomial_logit(
            @utility(begin
                1 ~ αtrain + βtravel_time * TRAIN_TT / 100 + βcost * (TRAIN_CO * (GA == 0)) / 100
                2 ~ αswissmetro + βtravel_time * SM_TT / 100 + βcost * SM_CO * (GA == 0) / 100
                3 ~ αcar + βtravel_time * CAR_TT / 100 + βcost * CAR_CO / 100

                αswissmetro = 0, fixed  # fix swissmetro ASC to zero
            end),
            :CHOICE,
            $data,
            availability=[
                1 => :avtr,
                2 => :avsm,
                3 => :avcar,
            ]
        ),
        samples=SAMPLES,
        seconds=1e6,
        evals=1
    )
end

# benchmark the actual macroexpand, since it would be only compiled once during the other benchmarking
# https://discourse.julialang.org/t/benchmark-macro-with-benchmarktools/72723/2
function benchmark_mnl_macroexpand()
    @benchmarkable(@macroexpand(@utility(begin
                1 ~ αtrain + βtravel_time * TRAIN_TT / 100 + βcost * (TRAIN_CO * (GA == 0)) / 100
                2 ~ αswissmetro + βtravel_time * SM_TT / 100 + βcost * SM_CO * (GA == 0) / 100
                3 ~ αcar + βtravel_time * CAR_TT / 100 + βcost * CAR_CO / 100

                αswissmetro = 0, fixed  # fix swissmetro ASC to zero
            end)),
            samples=SAMPLES,
            seconds=1e6,
            evals=1)
end

#########################################################
# NHTS Model                                            #
#########################################################

# convenience functions
# to save space in the repo, data are gzipped. helper function to read a
# gzipped csv "straight" into JuliaDB or DataFrames (via a temporary file)
function read_gzipped_csv(file, T; kwargs...)
    open(GzipDecompressorStream, file) do is
        mktemp() do path, os
            write(os, read(is))
            if T <: DTable
                DTable(CSV.File(path), 17500, kwargs...)
            elseif T <: DataFrame
                CSV.read(path, DataFrame, kwargs...)
            else
                error("Don't know how to read table into type $T")
            end
        end
    end
end

function benchmark_mnl_nhts_dtable()
    data = read_gzipped_csv(get_test_file_path("nhts/hhpub.csv.gz"), DTable)

    # topcode hhvehcnt
    data = map(r -> merge((hhveh_topcode=min(r.HHVEHCNT, 4),), pairs(r)), data)

    @benchmarkable(
        multinomial_logit(
            @utility(begin
                0 ~ α0
                1 ~ α1 + β1homeown * (HOMEOWN == 2) + @dummy_code(β1, HH_RACE, [2, 3, 4, 5, 6, 97]) + β1_hhsize * HHSIZE
                2 ~ α2 + β2homeown * (HOMEOWN == 2) + @dummy_code(β2, HH_RACE, [2, 3, 4, 5, 6, 97]) + β2_hhsize * HHSIZE
                3 ~ α3 + β3homeown * (HOMEOWN == 2) + @dummy_code(β3, HH_RACE, [2, 3, 4, 5, 6, 97]) + β3_hhsize * HHSIZE
                4 ~ α4plus + β4plushomeown * (HOMEOWN == 2) + @dummy_code(β4plus, HH_RACE, [2, 3, 4, 5, 6, 97]) + β4plus_hhsize * HHSIZE

                α0 = 0, fixed
            end),
            :hhveh_topcode,
            $data
        ),
        samples=SAMPLES,
        evals=1,
        seconds=1e6
    )
end

function benchmark_mnl_nhts_dataframe()
    data = read_gzipped_csv(get_test_file_path("nhts/hhpub.csv.gz"), DataFrame)

    data.hhveh_topcode = min.(data.HHVEHCNT, 4)

    @benchmarkable(
        multinomial_logit(
            @utility(begin
                0 ~ α0
                1 ~ α1 + β1homeown * (HOMEOWN == 2) + @dummy_code(β1, HH_RACE, [2, 3, 4, 5, 6, 97]) + β1_hhsize * HHSIZE
                2 ~ α2 + β2homeown * (HOMEOWN == 2) + @dummy_code(β2, HH_RACE, [2, 3, 4, 5, 6, 97]) + β2_hhsize * HHSIZE
                3 ~ α3 + β3homeown * (HOMEOWN == 2) + @dummy_code(β3, HH_RACE, [2, 3, 4, 5, 6, 97]) + β3_hhsize * HHSIZE
                4 ~ α4plus + β4plushomeown * (HOMEOWN == 2) + @dummy_code(β4plus, HH_RACE, [2, 3, 4, 5, 6, 97]) + β4plus_hhsize * HHSIZE

                α0 = 0, fixed
            end),
            :hhveh_topcode,
            $data
        ),
        samples=SAMPLES,
        evals=1,
        seconds=1e6
    )
end

function benchmark_mnl_nhts_macroexpand()
    @benchmarkable(@macroexpand(@utility(begin
                0 ~ α0
                1 ~ α1 + β1homeown * (HOMEOWN == 2) + @dummy_code(β1, HH_RACE, [2, 3, 4, 5, 6, 97]) + β1_hhsize * HHSIZE
                2 ~ α2 + β2homeown * (HOMEOWN == 2) + @dummy_code(β2, HH_RACE, [2, 3, 4, 5, 6, 97]) + β2_hhsize * HHSIZE
                3 ~ α3 + β3homeown * (HOMEOWN == 2) + @dummy_code(β3, HH_RACE, [2, 3, 4, 5, 6, 97]) + β3_hhsize * HHSIZE
                4 ~ α4plus + β4plushomeown * (HOMEOWN == 2) + @dummy_code(β4plus, HH_RACE, [2, 3, 4, 5, 6, 97]) + β4plus_hhsize * HHSIZE

                α0 = 0, fixed
            end)),
            samples=SAMPLES,
            seconds=1e6,
            evals=1)
end

function should_run(benchmark_name)
    run = length(ARGS) == 0 || benchmark_name in ARGS
    # INFO level debugging is disabled. If I were fancy I'd filter log messages by module, but... I'm not
    run || @warn "skipping benchmark $benchmark_name"
    return run
end

function main()
    "--verbose" in ARGS || Logging.disable_logging(Logging.Info)
    suite = BenchmarkGroup()

    should_run("MNL_Swissmetro_DTable") && (suite["MNL_Swissmetro_DTable"] = benchmark_mnl_dtable())
    should_run("MNL_Swissmetro_DataFrame") && (suite["MNL_Swissmetro_DataFrame"] = benchmark_mnl())
    should_run("MNL_Swissmetro_macroexpand") && (suite["MNL_Swissmetro_macroexpand"] = benchmark_mnl_macroexpand())
    should_run("MNL_NHTS_DTable") && (suite["MNL_NHTS_DTable"] = benchmark_mnl_nhts_dtable())
    should_run("MNL_NHTS_DataFrame") && (suite["MNL_NHTS_DataFrame"] = benchmark_mnl_nhts_dataframe())
    should_run("MNL_NHTS_macroexpand") && (suite["MNL_NHTS_macroexpand"] = benchmark_mnl_nhts_macroexpand())

    #tune!(suite)
    results = run(suite, verbose=true)

    show(IOContext(stdout, :compact => false), results)
end

# ugly that this has to be done at top level
addprocs()

@everywhere begin
    using Pkg
    Pkg.activate(Base.source_dir())
    using DiscreteChoiceModels
end

main()


text:
@testset "mdalternative interface" begin
    mpat = "Harwell-Boeing/psadmit/1138_bus"
    spat = "HB/1138_bus"
    @test mdlist(mm(spat)) == [mpat]
    @test mdlist(sp(mpat)) == [spat]
    @test mdlist(mm(mpat)) == [mpat]
    @test mdlist(sp(spat)) == [spat]
    @test mdlist(sp("vand")) == String[]
    @test mdlist(mm("Sandia/adder_dcop_17")) == String[]

    s1 = mdlist(sp(mm(:)))
    @test length(s1) == 455
    m1 = mdlist(mm(s1))
    @test length(m1) == 455
    @test m1 ⊆ mdlist(mm(:))

    m2 = mdlist(mm(sp(:)))
    @test length(m2) == 455
    s2 = mdlist(sp(m2))
    @test length(s2) == 455
    @test s2 ⊆ mdlist(sp(:))

    @test m1 == m2
end


text:
include("../src/HMat2d.jl");

function laplacian2d(nx,ny)
    M = speye(nx*ny);
    for x=1:nx
        for y=1:ny
            s = x+(y-1)*nx;
            M[s,s] = 4;
            if x > 1
                M[s,s-1] = -1;
            end
            if x < nx
                M[s,s+1] = -1;
            end
            if y > 1
                M[s,s-nx] = -1;
            end
            if y < ny
                M[s,s+nx] = -1;
            end
        end
    end
    return M;
end

EPS = 1e-6;
MaxRank = 2;
minn = 8;

for n = [16 32 64]
    A = full(laplacian2d(n,n));
    Z2C = Z2Cmapper([n,n],minn,reshape(1:n^2,n,n));
    A = full(A[Z2C,Z2C]);
    AH = HMat2dd2h(A, [n,n], [n,n], WEAK, [0,0], [0,0], 1, EPS, MaxRank, minn);
    BH = hcopy(AH);
    v = rand(n^2,5);
    tic();
    uhmat = hmatvec(AH,v);
    toc();
    tic();
    CH = hmul(AH,BH);
    toc();
    CD = hh2d(CH);
    @printf("N^2 = %4d^2, hmul error: %.3e\n",n,norm(CD-A*A)/norm(A*A));
end

max_stars_repo_name: YingzhouLi/HMat

max_issues_repo_name: YingzhouLi/HMat


text:
__precompile__()

module DistLearn

export try_fetch
include("utils.jl")

export Dataset
export load_datasets, get_num_samples, get_num_features
include("dataset.jl")

export Metric, MetricAccumulator, Accuracy
export evaluate, merge
include("metric.jl")

export Worker
export spawn_workers, invoke_on_workers
export accumulate_metric
include("worker.jl")

export make_arg_parser
include("arg_parse.jl")

end  # module DistLearn


text:
const INCLUDE_BRANCHES = Regex[
    r"^default$"i,
    r"^dev$"i,
    r"^develop$"i,
    r"^development$"i,
    r"^latest$"i,
    r"^main$"i,
    r"^master$"i,
    r"^production$"i,
    r"^stable$"i,
    r"^trunk$"i,
    r"^unstable$"i,
    ]

const EXCLUDE_BRANCHES = Regex[
    r"^gh-pages$"i,
    r"^gh-pages1$"i,
    ]