text:
using ..Cos
import Base: get
using ..Cos: CosDocImpl

mutable struct PDDocImpl <: PDDoc
    cosDoc::CosDocImpl
    catalog::CosIndirectObject{CosDict}
    pages::CosObject
    structTreeRoot::CosObject
    isTagged::Symbol #Valid values :tagged, :none and :suspect
    fonts::Dict{CosObject, PDFont}
    xobjs::Dict{CosObject, PDXObject}
    function PDDocImpl(fp::AbstractString)
        cosDoc = cosDocOpen(fp)
        catalog = cosDocGetRoot(cosDoc)
        new(cosDoc,catalog,CosNull,CosNull,:none,
            Dict{CosObject, PDFont}(), Dict{CosObject, PDXObject}())
    end
end

function pdDocGetPage(doc::PDDocImpl, num::Int)
    cosobj = find_page_from_treenode(doc.pages, num)
    return create_pdpage(doc, cosobj)
end

"""
```
    show(io::IO, doc::PDDoc)
```
Prints the PDDoc. The intent is to print lesser information from the structure.
"""
function Base.show(io::IO, doc::PDDoc)
    print(io, "\nPDDoc ==>\n")
    print(io, doc.cosDoc)
    print(io, '\n')
    print(io, "Catalog:")
    print(io, doc.catalog)
    print(io, "isTagged: $(doc.isTagged)\n")
end

"""
Recursively reads the page object and populates the indirect objects
Ensures indirect objects are read and updated in the xref Dictionary.
"""
populate_doc_pages(doc::PDDocImpl, dict::CosIndirectObject{CosDict}) =
    populate_doc_pages(doc, dict.obj)

function populate_doc_pages(doc::PDDocImpl, dict::CosDict)
    if (cn"Pages" == get(dict, cn"Type"))
        kids = cosDocGetObject(doc.cosDoc, dict, cn"Kids")
        arr = get(kids)
        len = length(arr)
        for i=1:len
            ref = splice!(arr, 1)
            obj = cosDocGetObject(doc.cosDoc, ref)
            if obj !== CosNull
                push!(arr, obj)
                populate_doc_pages(doc, obj)
            end
        end
    end
    parent = get(dict, cn"Parent")
    if (parent === CosNull)
        obj = cosDocGetObject(doc.cosDoc, parent)
        set!(dict,CosName("Parent"),obj)
    end
    return nothing
end

populate_doc_pages(doc::PDDocImpl, dict::CosObject) = nothing

function update_page_tree(doc::PDDocImpl)
    pagesref = get(doc.catalog, cn"Pages")::CosIndirectObjectRef
    pages = cosDocGetObject(doc.cosDoc, pagesref)::CosIndirectObject{CosDict}
    populate_doc_pages(doc, pages)
    doc.pages = pages
    return nothing
end

#=
This implementation may seem non-intuitive to some due to recursion and also
seemingly non-standard way of computing page count. However, PDF spec does not
discount the possibility of an intermediate node having page and pages nodes
in the kids array. Hence, this implementation.
=#
function find_page_from_treenode(node::CosObject, pageno::Int)::CosObject
    mytype = get(node, cn"Type")
    #If this is a page object the pageno has to be 1
    if mytype == cn"Page"
        pageno == 1 && return node
        throw(ErrorException(E_INVALID_PAGE_NUMBER))
    end
    kids = get(node, cn"Kids")
    kidsarr = get(kids)
    sum = 0
    for kid in kidsarr
        cnt = Cos.get_internal_pagecount(kid)
        (sum + cnt) >= pageno &&
            return find_page_from_treenode(kid, pageno-sum)::CosObject
        sum += cnt
    end
    throw(ErrorException(E_INVALID_PAGE_NUMBER))
end

# The structure tree is not fully loaded but the object linkages are established for future
# correlations during text extraction.

function update_structure_tree!(doc::PDDocImpl)
    catalog = pdDocGetCatalog(doc)
    marking = get(catalog, cn"MarkInfo")

    if (marking !== CosNull)
        tagged = get(marking, cn"Marked")
        suspect = get(marking, cn"Suspect")
        doc.isTagged = (suspect === CosTrue) ? (:suspect) :
                       (tagged  === CosTrue) ? (:tagged)  : (:none)
    end

    structTreeRef = get(catalog, cn"StructTreeRoot")
    doc.structTreeRoot = cosDocGetObject(doc.cosDoc, structTreeRef)
    return nothing
end

get_pd_font!(doc::PDDocImpl, cosfont::CosObject) =
    get!(doc.fonts, cosfont, PDFont(doc, cosfont))

get_pd_xobject!(doc::PDDocImpl, cosxobj::CosObject) =
    get!(doc.xobjs, cosxobj, createPDXObject(doc, cosxobj))


text:
@agent SchellingAgent GridAgent{2} begin
    mood::Bool # whether the agent is happy in its position. (true = happy)
    group::Int # The group of the agent,  determines mood as it interacts with neighbors
end

"""
``` julia
schelling(;
    numagents = 320,
    griddims = (20, 20),
    min_to_be_happy = 3,
)
```
Same as in [Schelling's segregation model](@ref).
"""
function schelling(; numagents = 320, griddims = (20, 20), min_to_be_happy = 3)
    @assert numagents < prod(griddims)
    space = GridSpace(griddims, periodic = false)
    properties = Dict(:min_to_be_happy => min_to_be_happy)
    model = ABM(SchellingAgent, space; properties, scheduler = Schedulers.randomly)
    for n in 1:numagents
        agent = SchellingAgent(n, (1, 1), false, n < numagents / 2 ? 1 : 2)
        add_agent_single!(agent, model)
    end
    return model, schelling_agent_step!, dummystep
end

function schelling_agent_step!(agent, model)
    agent.mood == true && return # do nothing if already happy
    count_neighbors_same_group = 0
    for neighbor in nearby_agents(agent, model)
        if agent.group == neighbor.group
            count_neighbors_same_group += 1
        end
    end
    if count_neighbors_same_group ≥ model.min_to_be_happy
        agent.mood = true
    else
        move_agent_single!(agent, model)
    end
end

ext: jl


text:
#This is an additional wrapper to generate julia function definitions
# from the netcdf.h header file. This is necessary because it looks like
# Clang.jl had problems translating size_t and ptrdiff_t to the correct
# Julia types

#A small translator for c types:
tdict = Dict(
    "int"=>"Cint",
    "size_t"=>"Csize_t",
    "char"=>"Cchar",
    "NC_Dispatch"=>"Cvoid",
    "nc_type"=>"Cint",
    "nc_vlen_t"=>"nc_vlen_t",
    "void"=>"Cvoid",
    "long long"=>"Clonglong",
    "ptrdiff_t"=>"Cptrdiff_t",
    "unsigned int"=>"Cuint",
    "float"=>"Cfloat",
    "unsigned char"=>"Cuchar",
    "signed char"=>"Cchar",
    "short"=>"Cshort",
    "long"=>"Clong",
    "double"=>"Cdouble",
    "unsigned short"=>"Cushort",
    "unsigned long"=>"Culong",
    "unsigned long long"=>"Culonglong",
)

"Parses the signature from the c function argument and returns the Julia type as well as the argument name"
function parsesig(s,tdict)
    s = strip(replace(s,"const"=>""))
    nstar = count(isequal('*'),s)
    s = replace(s,"*"=>"")
    if endswith(s,"[]")
        nstar += 1
        s = replace(s,"[]"=>"")
    end
    spl = split(strip(s)," ", keepempty=false)
    if length(spl) == 1
        return tdict[spl[1]],"x"
    end
    if length(spl)>2
        spl = (join(spl[1:end-1]," "), spl[end])
    end
    t, argname = spl
    if nstar == 2
        return "Ptr{Ptr{$(tdict[t])}}", argname
    elseif nstar == 1
        return "Ptr{$(tdict[t])}", argname[2:end]
    else
        return tdict[t], argname
    end
end

"Write the function with name `fname` and signature `sig` to the IO stream `outp`"
function writefun(outp,fname,sig)
    print(outp,"function ")
    print(outp,fname)
    print(outp,"(")
    for (t,name) in sig
        print(outp,name)
#         if parse(t) <: Integer
#             print(outp,"::Integer")
#         end
        print(outp,", ")
    end
    print(outp, ")\n")
    println(outp, "    check(ccall(")
    println(outp, "        (:$(fname), libnetcdf),")
    println(outp, "        Cint,")
    print(outp, "        (")
    foreach(i->print(outp,string(i[1],", ")),sig)
    println(outp,"),")
    foreach(i->println(outp,"        $(i[2]),"),sig)
    println(outp, "    ))")
    println(outp, "end")
end

include(joinpath(@__DIR__,"..","src","netcdf_constants.jl"))

function parseheader(input, output)
    cheader = eachline(input) |> collect;
    open(output,"w") do outp
        ii = 1
        while true
        i1 = findnext(startswith("EXTERNL int"), cheader,ii)
        i1 === nothing && break
        i2 = findnext(contains(";"), cheader,i1)
        cdef = join(cheader[i1:i2])
        sig = match(r"int\s*(\w+)\((.*)\)",cdef)
        if sig !== nothing
            fname = sig.captures[1]
            sig = parsesig.(split(sig.captures[2],","), Ref(tdict))
            writefun(outp,fname,sig)
        end
        ii = i2+1
        end
    end
end

input = eachline(`locate netcdf.h`) |> first
output = joinpath(@__DIR__, "..", "src","netcdf_c.jl")
parseheader(input, output)


text:
module TestWiringDiagramCore

using Test
using Catlab.Doctrines, Catlab.WiringDiagrams

# For testing purposes, check equality of port symbols.
function WiringDiagramCore.validate_ports(source_port::Symbol, target_port::Symbol)
  if source_port != target_port
    throw(PortValueError(source_port, target_port))
  end
end

# Low-level graph interface
###########################

A, B, C, D = Ob(FreeSymmetricMonoidalCategory, :A, :B, :C, :D)
f = Hom(:f, A, B)
g = Hom(:g, B, C)
h = Hom(:h, C, D)

# Operations on boxes
d = WiringDiagram(A, C)
@test nboxes(d) == 0
@test_throws KeyError box(d,input_id(d))
@test_throws KeyError box(d,output_id(d))

fv = add_box!(d, f)
@test nboxes(d) == 1
@test box(d, fv).value == :f
@test box(d, fv) == Box(f)
rem_box!(d, fv)
@test nboxes(d) == 0

fv = add_box!(d, f)
gv = add_box!(d, g)
hv = add_box!(d, h)
@test nboxes(d) == 3
@test [b.value for b in boxes(d)] == [:f,:g,:h]
@test boxes(d) == [Box(f),Box(g),Box(h)]
rem_boxes!(d, [fv,hv])
@test nboxes(d) == 1
@test boxes(d) == [Box(g)]

# Operations on ports
d = WiringDiagram(A, C)
fv = add_box!(d, f)
gv = add_box!(d, g)

@test input_ports(d, fv) == [:A]
@test output_ports(d, fv) == [:B]
@test input_ports(d, output_id(d)) == [:C]
@test output_ports(d, input_id(d)) == [:A]
@test_throws ErrorException input_ports(d, input_id(d))
@test_throws ErrorException output_ports(d, output_id(d))

@test port_value(d, Port(input_id(d),OutputPort,1)) == :A
@test port_value(d, Port(output_id(d),InputPort,1)) == :C
@test port_value(d, Port(fv,InputPort,1)) == :A
@test port_value(d, Port(fv,OutputPort,1)) == :B

# Operations on wires
@test nwires(d) == 0
@test !has_wire(d, fv, gv)
@test !has_wire(d, (fv,1) => (gv,1))
add_wire!(d, (input_id(d),1) => (fv,1))
add_wire!(d, (fv,1) => (gv,1))
add_wire!(d, (gv,1) => (output_id(d),1))
@test nwires(d) == 3
@test has_wire(d, fv, gv)
@test has_wire(d, (fv,1) => (gv,1))
@test_throws PortValueError add_wire!(d, (gv,1) => (fv,1))
@test wires(d) == map(Wire, [
  (input_id(d),1) => (fv,1),
  (fv,1) => (gv,1),
  (gv,1) => (output_id(d),1),
])

# Shallow copies.
d_copy = copy(d)
rem_boxes!(d_copy, [fv,gv])
@test nboxes(d) == 2
@test nwires(d) == 3

# Graph properties.
@test Set(all_neighbors(d, fv)) == Set([input_id(d),gv])
@test Set(all_neighbors(d, gv)) == Set([fv,output_id(d)])
@test neighbors(d, fv) == [gv]
@test outneighbors(d, fv) == [gv]
@test inneighbors(d, gv) == [fv]
@test wires(d, input_id(d)) == [ Wire((input_id(d),1) => (fv,1)) ]
@test wires(d, fv) == map(Wire, [
  ((input_id(d),1) => (fv,1)),
  ((fv,1) => (gv,1))
])
@test out_wires(d, fv) == [ Wire((fv,1) => (gv,1)) ]
@test out_wires(d, Port(fv,OutputPort,1)) == [ Wire((fv,1) => (gv,1)) ]
@test in_wires(d, gv) == [ Wire((fv,1) => (gv,1)) ]
@test in_wires(d, Port(gv,InputPort,1)) == [ Wire((fv,1) => (gv,1)) ]

rem_wires!(d, fv, gv)
@test nwires(d) == 2
@test !has_wire(d, fv, gv)
rem_wire!(d, (input_id(d),1) => (fv,1))
@test wires(d) == [ Wire((gv,1) => (output_id(d),1)) ]

# Substitution
f, g, h = Hom(:f,A,B), Hom(:g,B,C), Hom(:h,C,D)
sub = WiringDiagram(B,D)
gv = add_box!(sub, g)
hv = add_box!(sub, h)
add_wires!(sub, Pair[
  (input_id(sub),1) => (gv,1),
  (gv,1) => (hv,1),
  (hv,1) => (output_id(sub),1),
])
d = WiringDiagram(A,D)
fv = add_box!(d, f)
subv = add_box!(d, sub)
add_wires!(d, Pair[
  (input_id(d),1) => (fv,1),
  (fv,1) => (subv,1),
  (subv,1) => (output_id(d),1),
])
@test boxes(d) == [ Box(f), sub ]
@test boxes(sub) == [ Box(g), Box(h) ]
substitute!(d, subv)
@test nboxes(d) == 3
@test Set(boxes(d)) == Set([ Box(f), Box(g), Box(h) ])
box_map = Dict(box(d,v).value => v for v in box_ids(d))
@test nwires(d) == 4
@test Set(wires(d)) == Set(map(Wire, [
  (input_id(d),1) => (box_map[:f],1),
  (box_map[:f],1) => (box_map[:g],1),
  (box_map[:g],1) => (box_map[:h],1),
  (box_map[:h],1) => (output_id(d),1),
]))

# Encapsulation
d0 = WiringDiagram(A,D)
fv = add_box!(d0, f)
gv = add_box!(d0, g)
hv = add_box!(d0, h)
add_wires!(d0, Pair[
  (input_id(d),1) => (fv,1),
  (fv,1) => (gv,1),
  (gv,1) => (hv,1),
  (hv,1) => (output_id(d),1)
])

d = deepcopy(d0)
encapsulate!(d, [fv,gv])
@test nboxes(d) == 2
@test sum(isa(b, WiringDiagram) for b in boxes(d)) == 1
@test nwires(d) == 3
sub = first(b for b in boxes(d) if isa(b, WiringDiagram))
@test nboxes(sub) == 2
@test Set(boxes(sub)) == Set([ Box(f), Box(g) ])
box_map = Dict(box(sub,v).value => v for v in box_ids(sub))
@test Set(wires(sub)) == Set(map(Wire, [
  (input_id(sub),1) => (box_map[:f],1),
  (box_map[:f],1) => (box_map[:g],1),
  (box_map[:g],1) => (output_id(sub),1),
]))

d = deepcopy(d0)
encapsulate!(d, [fv,gv], discard_boxes=true, value=:e)
@test Set(boxes(d)) == Set([ Box(:e, [:A], [:C]), Box(h) ])
box_map = Dict(box(d,v).value => v for v in box_ids(d))
@test Set(wires(d)) == Set(map(Wire, [
  (input_id(d),1) => (box_map[:e],1),
  (box_map[:e],1) => (box_map[:h],1),
  (box_map[:h],1) => (output_id(d),1),
]))

d0 = WiringDiagram(A,B)
v1 = add_box!(d0, f)
v2 = add_box!(d0, f)
add_wires!(d0, Pair[
  (input_id(d),1) => (v1,1),
  (input_id(d),1) => (v2,1),
  (v1,1) => (output_id(d),1),
  (v2,1) => (output_id(d),1),
])
d = deepcopy(d0)
encapsulate!(d, [v1,v2])
@test nboxes(d) == 1
@test nwires(d) == 2
sub = first(boxes(d))
@test sub == d0

# High-level categorical interface
##################################

# Category
#---------

# Generators
f = singleton_diagram(Box(Hom(:f,A,B)))
g = singleton_diagram(Box(Hom(:g,B,A)))
@test nboxes(f) == 1
@test boxes(f) == [ Box(Hom(:f,A,B)) ]
@test nwires(f) == 2

# Composition
@test nboxes(compose(f,g)) == 2
@test boxes(compose(f,g)) == [ Box(Hom(:f,A,B)), Box(Hom(:g,B,A)) ]
@test nwires(compose(f,g)) == 3

# Domains and codomains
@test dom(f) == Ports([:A])
@test codom(f) == Ports([:B])
@test dom(compose(f,g)) == Ports([:A])
@test codom(compose(f,g)) == Ports([:A])
@test_throws Exception compose(f,f)

# Associativity
@test compose(compose(f,g),f) == compose(f,compose(g,f))

# Identity
@test compose(id(dom(f)), f) == f
@test compose(f, id(codom(f))) == f

# Symmetric monoidal category
#----------------------------

# Domains and codomains
@test dom(otimes(f,g)) == otimes(dom(f),dom(g))
@test codom(otimes(f,g)) == otimes(codom(f),codom(g))

# Associativity and unit
X, Y = Ports([:A,:B]), Ports([:C,:D])
I = munit(Ports)
@test otimes(X,I) == X
@test otimes(I,X) == X
@test otimes(otimes(X,Y),X) == otimes(X,otimes(Y,X))
@test otimes(otimes(f,g),f) == otimes(f,otimes(g,f))

# Braiding
@test compose(braid(X,Y),braid(Y,X)) == id(otimes(X,Y))

# Permutations
W = otimes(X,Y)
@test permute(W, [1,2,3,4]) == id(W)
@test permute(W, [1,2,3,4], inverse=true) == id(W)
@test permute(W, [3,4,1,2]) == braid(X,Y)
@test permute(W, [3,4,1,2], inverse=true) == braid(Y,X)
@test_throws AssertionError permute(W, [1,2])

# Diagonals
#----------

# Basic composition
d = WiringDiagram(dom(f), otimes(codom(f),codom(f)))
fv1 = add_box!(d, first(boxes(f)))
fv2 = add_box!(d, first(boxes(f)))
add_wires!(d, [
  (input_id(d),1) => (fv1,1),
  (input_id(d),1) => (fv2,1),
  (fv1,1) => (output_id(d),1),
  (fv2,1) => (output_id(d),2),
])
@test compose(mcopy(dom(f)), otimes(f,f)) == d

# Domains and codomains
@test dom(mcopy(Ports([:A]))) == Ports([:A])
@test codom(mcopy(Ports([:A]))) == Ports([:A,:A])
@test dom(mcopy(Ports([:A,:B]),3)) == Ports([:A,:B])
@test codom(mcopy(Ports([:A,:B]),3)) == Ports([:A,:B,:A,:B,:A,:B])

# Associativity
X = Ports([:A])
@test compose(mcopy(X), otimes(id(X),mcopy(X))) == mcopy(X,3)
@test compose(mcopy(X), otimes(mcopy(X),id(X))) == mcopy(X,3)

# Commutativity
@test compose(mcopy(X), braid(X,X)) == mcopy(X)

# Unit
@test compose(mcopy(X), otimes(id(X),delete(X))) == id(X)

# Codiagonals
#------------

# Domains and codomains
@test dom(mmerge(Ports([:A]))) == Ports([:A,:A])
@test codom(mmerge(Ports([:A]))) == Ports([:A])
@test dom(mmerge(Ports([:A,:B]),3)) == Ports([:A,:B,:A,:B,:A,:B])
@test codom(mmerge(Ports([:A,:B]),3)) == Ports([:A,:B])

# Associativity
X = Ports([:A])
@test compose(otimes(id(X),mmerge(X)), mmerge(X)) == mmerge(X,3)
@test compose(otimes(mmerge(X),id(X)), mmerge(X)) == mmerge(X,3)

# Commutativity
@test compose(braid(X,X), mmerge(X)) == mmerge(X)

# Unit
@test compose(otimes(id(X),create(X)), mmerge(X)) == id(X)

end


text:
function read_request(channel::gRPCChannel, controller::gRPCController, services)
    connection = channel.session
    evt = Session.take_evt!(connection)

    if isa(evt, EvtGoaway)
        @debug("received EvtGoaway")
        close(channel)
        return (nothing, nothing, nothing)
    end

    if !isa(evt, EvtRecvHeaders)
        @debug("unexpected event while reading request (closing channel)", evt)
        close(channel)
        return (nothing, nothing, nothing)
    end

    channel.stream_id = evt.stream_identifier
    headers = Headers()
    while isa(evt, EvtRecvHeaders)
        merge!(headers, evt.headers)
        evt.is_end_stream && break
        evt = Session.take_evt!(connection)
    end

    method = headers[":method"]
    path = headers[":path"]
    pathcomps = split(path, "/"; keepempty=false)
    if length(pathcomps) != 2
        @debug("unexpected path while reading request (closing channel)", path)
        close(channel)
        return (nothing, nothing, nothing)
    end

    servicename, methodname = split(path, "/"; keepempty=false)

    data = UInt8[]
    while isa(evt, EvtRecvData)
        append!(data, evt.data)
        evt.is_end_stream && break
        evt = Session.take_evt!(connection)
    end

    @debug("received request", method, path, stream_id=channel.stream_id, servicename, methodname, nbytes=length(data))

    service = services[servicename]
    method = find_method(service, methodname)
    request_type = get_request_type(service, method)
    request = request_type()
    from_delimited_message_bytes(data, request)

    service, method, request
end

# TODO: use http connection settings to delimit frame size
function write_response(channel::gRPCChannel, controller::gRPCController, response)
    sending_headers = Dict(":status" => "200", "content-type" => "application/grpc")
    data_buff = to_delimited_message_bytes(response)

    Session.put_act!(channel.session, Session.ActSendHeaders(channel.stream_id, sending_headers, false))
    Session.put_act!(channel.session, Session.ActSendData(channel.stream_id, data_buff, false))
    Session.put_act!(channel.session, Session.ActSendHeaders(channel.stream_id, Dict("grpc-status" => "0"), true))
    @info "done with request"
    nothing
end

function write_error_response(channel::gRPCChannel, controller::gRPCController, err::Exception)
    sending_headers = Dict(":status" => "200", "content-type" => "application/grpc")
    ex = wrap(err)

    Session.put_act!(channel.session, Session.ActSendHeaders(channel.stream_id, sending_headers, false))
    Session.put_act!(channel.session, Session.ActSendData(channel.stream_id, to_delimited_message_bytes(""), false))
    Session.put_act!(channel.session, Session.ActSendHeaders(channel.stream_id, Dict("grpc-status" => string(ex.code), "grpc-message" => ex.message), true))
    @info "done with request"
    nothing
end

function call_method(channel::gRPCChannel, service::ServiceDescriptor, method::MethodDescriptor, controller::gRPCController, request)
    write_request(channel, controller, service, method, request)
    response_type = get_response_type(method)
    response = response_type()
    read_response(channel, controller, response)
end

# gRPC server implementation
#
mutable struct gRPCServer
    sock::TCPServer
    services::Dict{String, ProtoService}
    run::Bool

    gRPCServer(services::Tuple{ProtoService}, ip::IPv4, port::Integer) = gRPCServer(services, listen(ip, port))
    gRPCServer(services::Tuple{ProtoService}, port::Integer) = gRPCServer(services, listen(port))
    function gRPCServer(services::Tuple{ProtoService}, sock::TCPServer)
        svcdict = Dict{String,ProtoService}()
        for svc in services
            svcdict[svc.desc.name] = svc
        end
        new(sock, svcdict, true)
    end
end

# TODO: close all channels, wait for/interrupt processors
close(srvr::gRPCServer) = close(srvr.sock)

function process(controller::gRPCController, srvr::gRPCServer, channel::gRPCChannel)
    @info("start processing channel")
    try
        while(!channel.session.closed)
            wait(channel.session.channel_evt)
            service, method, request = read_request(channel, controller, srvr.services)
            (service === nothing) && continue

            try
                response = call_method(service, method, controller, request)
                write_response(channel, controller, response)
            catch e
                write_error_response(channel, controller, e)
            end
        end
    catch ex
        @warn("channel stopped with exception", ex)
    end
    # TODO: close channel if not closed, remove reference from srvr
    @info("stopped processing channel")
end

function run(srvr::gRPCServer)
    controller = gRPCController()
    try
        while(srvr.run)
            buffer = accept(srvr.sock)
            connection = Session.new_connection(buffer; isclient=false)

            # TODO: keep channel reference in server
            channel = gRPCChannel(connection)
            @async process(controller, srvr, channel)
        end
    catch ex
        @warn("server stopped with exception", ex)
    end
    @info("stopped server")
end


text:
using JsonGrinder
using Test

include("extractors.jl")
include("multirepresentation.jl")
include("pipelinetest.jl")
include("schema.jl")
include("show_html.jl")


text:
# setup filepaths
path_source = string(@__DIR__,"\\..\\source");
path_dev = string(@__DIR__,"\\..\\dev");

# LOAD IN PACKAGES AND DATA
using DataFrames, CSV, Plots, StatsPlots, StatsBase, Missings
df = DataFrame(CSV.read(joinpath(path_dev,"19-21VAERSCOMB.csv"), DataFrame));

# Filter dataset to reports with non-missing Age > 16
df = df[(.!ismissing.(df.AGE_YRS)) .& (df.AGE_YRS .>= 16), :] #filter age >= 16
first(df, 10)

#labeling of events as serious or non-serious
a = (df.DIED .== "Y") .| (df.L_THREAT .== "Y") .| (df.HOSPITAL .== "Y") .| (df.X_STAY .== "Y") .| (df.DISABLE .== "Y") .| (df.BIRTH_DEFECT .== "Y")
df.SERIOUS_EVENT = replace(a, missing => 0)
show(names(df))

colnames = ["VAERS_ID", "RECVDATE", "AGE_YRS", "VAX_NAME", "VAX_TYPE", "VAX_MANU", "SYMPTOMS", "SERIOUS_EVENT"]
df_cleaned = select(df, colnames)

# write to dev
CSV.write(joinpath(path_dev,"19-21VAERSCOMB_clean.csv"), df_cleaned);


text:
module C_
    using MultiInterface

    abstract type Subby{T} end
    abstract type Subby2{N,T} <: Subby{T} end
    struct SubS{N,T} <: Subby2{N,T} end
    struct SubS2{N,T} <: Subby2{2,SubS{N,T}} end

    Base.print_without_params(x::Type{Subby}) = false
    Base.print_without_params(x::Type{Subby2}) = false
    Base.print_without_params(x::Type{SubS}) = false
    Base.print_without_params(x::Type{SubS2}) = false

    @interface sub(s::Subby{<:T}, t::T)::T where T
    @interface sub2(s::Subby2{2,<:SubS{N,<:T}}, t::T) where {N,T}

    @impl begin
        struct Subby2Sub end
        function sub(s::Subby2{N,T}, t::T)::T where {N,T}
            t
        end
    end

    @impl begin
        struct Subby2Sub1 end
        # This should not break
        function sub(s::Subby2{N,Float64}, t::AbstractFloat)::AbstractFloat where {N}
            t
        end

    end

    error = nothing
    try
        @impl begin
            # this should break
            struct Subby2Sub2 end
            function sub(s::Subby2{N,AbstractFloat}, t::Float64)::Float64 where {N}
                t
            end
        end
    catch e
        global error = e
    end

    @impl begin
        struct Subby2Sub3 end
        function sub2(s::SubS2{N,T}, t::T)::T where {N,T}
            t
        end
    end
end

@testset "subtype testing" begin
    @test C_.sub(C_.SubS{1,Int64}(), 10) == 10
    @test C_.sub(C_.SubS{1,Symbol}(), :a) == :a
    @test C_.sub(C_.SubS{1, Float64}(), 1.1) == 1.1
    @test C_.sub2(C_.SubS2{10,Int64}(), 2) == 2
    @test C_.error isa ErrorException
    @test occursin("Invalid typing--interface typevars misbound", C_.error.msg)
    @test occursin(r"rule        :  var\"#s\d+\" <: T", C_.error.msg)
    @test occursin("bound rule  :  AbstractFloat <: Float64", C_.error.msg)
end

module D_
    using MultiInterface

    abstract type SupD{S,T} end
    struct SubD{T} <: SupD{Nothing, T} end

    Base.print_without_params(x::Type{SupD}) = false
    Base.print_without_params(x::Type{SubD}) = false

    @interface retval1(x::SubD{T})::T where T

    @impl begin
        struct SubDRetVal1 end
        function retval1(t::SubD{T})::T where T
            T isa Type{<:Integer} && return 5
            T isa Type{Symbol} && return :b
            T
        end
    end

    @interface retval2(x::SubD{<:T})::T where T

    @impl begin
        struct SubDSubRetVal2 end
        function retval2(t::SubD{T})::T where T
            T isa Type{<:Integer} && return 10
            T isa Type{Symbol} && return :a
            T
        end
    end
end

@testset "return value testing" begin
    @test D_.retval1(D_.SubD{Int64}()) == 5
    @test D_.retval1(D_.SubD{Symbol}()) == :b
    @test D_.retval2(D_.SubD{Int64}()) == 10
    @test D_.retval2(D_.SubD{Symbol}()) == :a
end

module E_
    using MultiInterface

    @interface ten(t::T)::T where T
    @interface eleven(x::Vector{T}, t::T) where T
    @interface twelve(x::(Vector{T} where T), t::(T where T))

    tenerror = nothing
    try
        @impl begin
            struct Ten_ end
            function ten(t::Int64)::Float64
                1.0
            end
        end
    catch e
        global tenerror = e
    end

    elevenerror1 = nothing
    try
        @impl begin
            struct Eleven1_ end
            function eleven(x::Vector{T}, t::Float64)::T where T
                t
            end
        end
    catch e
        global elevenerror1 = e
    end

    elevenerror2 = nothing
    try
        @impl begin
            struct Eleven2_ end
            function eleven(x::Vector{Int64}, t::Float64)::Number
                t
            end
        end
    catch e
        global elevenerror2 = e
    end

    @impl begin
        struct Twelve_ end
        function twelve(x::Vector, t::Int64)::Int64
            t
        end
        function twelve(x::Vector{Float64}, t::Int64)
            x
        end
    end
end

@testset "matching and non-matching typevars" begin
    @test E_.tenerror isa ErrorException
    @test contains(E_.tenerror.msg, "Invalid typing--interface typevar bound to multiple types")
    @test contains(E_.tenerror.msg, "typevar:  T")
    @test contains(E_.tenerror.msg, "types  :  Any[:Int64, :Float64]")

    @test E_.elevenerror1 isa ErrorException
    @test contains(E_.elevenerror1.msg, "Invalid typing--interface typevar bound to multiple types")
    @test contains(E_.elevenerror1.msg, "typevar:  T")
    @test contains(E_.elevenerror1.msg, "types  :  Any[:T, :Float64]")

    @test E_.elevenerror2 isa ErrorException
    @test contains(E_.elevenerror2.msg, "Invalid typing--interface typevar bound to multiple types")
    @test contains(E_.elevenerror2.msg, "typevar:  T")
    @test contains(E_.elevenerror2.msg, "types  :  Any[:Int64, :Float64]")

    @test E_.twelve([:a,:b,"a",1], 10) == 10
    @test E_.twelve([1.0,2.0], 10) == [1.0,2.0]
end

module F_
    using MultiInterface

    @interface thirteen(t::Tuple{AbstractFloat, Signed})
    @interface fourteen(x::Tuple{T}, t::T)::T where T
    @interface fifteen(x::Tuple{S,T}, s::S, t::T) where {S,T}
    @interface sixteen(x::Vector{<:Tuple{S,T}}, s::S, t::T) where {S,T}

    @impl begin
        struct Thirteen_ end
        function thirteen(t::Tuple{Float64, Int64})
            t
        end
    end

    @impl begin
        struct Fourteen_ end
        function fourteen(x::Tuple{Symbol}, t::Symbol)::Symbol
            Symbol(string(x[1], t))
        end

        function fourteen(x::Tuple{Tuple{T}}, t::Tuple{T})::Tuple{T} where T
            x[1]
        end
    end

    fourteenerror1 = nothing
    try
        @impl begin
            struct Fourteen_1 end
            function fourteen(x::Tuple{Tuple{T}}, t::T)::T where T
                Symbol(string(x[1], t))
            end
        end
    catch e
        global fourteenerror1 = e
    end

    @impl begin
        struct Fifteen_ end
        function fifteen(x::Tuple{AbstractFloat, Signed}, s::AbstractFloat, t::Signed)
            x
        end
        function fifteen(x::Tuple{Symbol, Symbol}, s::Symbol, t::Symbol)
            x
        end
    end

    fifteenerror1 = nothing
    try
        @impl begin
            struct Fifteen_1 end
            function fifteen(x::Tuple{Float64, Signed}, s::AbstractFloat, t::Signed)
                x
            end
        end
    catch e
        global fifteenerror1 = e
    end

    fifteenerror2 = nothing
    try
        @impl begin
            struct Fifteen_1 end
            function fifteen(x::Tuple, s::Float64, t::Symbol)
                x
            end
        end
    catch e
        global fifteenerror2 = e
    end

    @impl begin
        struct Sixteen_ end
        function sixteen(x::Vector{<:Tuple{AbstractFloat, Signed}}, s::AbstractFloat, t::Signed)
            (s,t)
        end

        function sixteen(x::Vector{Tuple{Float64, Int64}}, s::Float64, t::Int64)
            s
        end

        function sixteen(x::Vector{Tuple{Float32, Int32}}, s::Float32, t::Int32)
            t
        end

        function sixteen(x::Vector{<:Tuple{Float64, B}}, s::Float64, t::B) where B
            B
        end
    end

    sixteenerror1 = nothing
    try
        @impl begin
            struct Sixteen_1 end
            function sixteen(x::Vector{Tuple{Symbol,Symbol}}, s::Symbol, t::Integer)
                x
            end
        end
    catch e
        global sixteenerror1 = e
    end

    sixteenerror2 = nothing
    try
        @impl begin
            struct Sixteen_2 end
            function sixteen(x::Vector{Tuple{A,B}}, s::A, t::Integer) where {A,B}
                x
            end
        end
    catch e
        global sixteenerror2 = e
    end

end

@testset "tuples" begin
    @test F_.thirteen((1.0, 1)) == (1.0, 1)
    @test F_.fourteen((:a,), :b) == :ab
    @test F_.fourteen(((1,),), (2,)) == (1,)

    @test F_.fourteenerror1 isa ErrorException
    @test contains(F_.fourteenerror1.msg, "Invalid typing--interface typevar bound to multiple types")
    @test contains(F_.fourteenerror1.msg, "typevar:  T")
    @test contains(F_.fourteenerror1.msg, "types  :  Any[:(Tuple{T}), :T]")

    @test F_.fifteen((1.0, 1), 2.0, 2) == (1.0, 1)
    @test F_.fifteen((:a, :b), :c, :d) == (:a, :b)

    @test F_.fifteenerror1 isa ErrorException
    @test contains(F_.fifteenerror1.msg, "Invalid typing--interface typevar bound to multiple types")
    @test contains(F_.fifteenerror1.msg, "typevar:  S")
    @test contains(F_.fifteenerror1.msg, "types  :  Any[:Float64, :AbstractFloat]")

    @test F_.fifteenerror2 isa ErrorException
    @test contains(F_.fifteenerror2.msg, "Invalid definition of fifteen")
    @test contains(F_.fifteenerror2.msg, "defined :  DataType[Tuple, Float64, Symbol]")
    @test contains(F_.fifteenerror2.msg, "expected:  Type[Tuple{S, T} where {S, T}, Any, Any]")
    @test contains(F_.fifteenerror2.msg, "at      :  Bool[0, 1, 1]")

    @test F_.sixteen(Tuple{Float16, Int16}[(1.0,1)], convert(Float16, 2.0), convert(Int16, 2)) == (2.0,2)
    @test F_.sixteen([(1.0,1)], 2.0, 2) == 2.0
    @test F_.sixteen(Tuple{Float32, Int32}[(1.0,1)], 2.0f0, convert(Int32, 2)) == 2
    @test F_.sixteen([(1.0,"bar")], 2.0, "foo") == String
    @test typeof(F_.sixteenerror1) == ErrorException
    @test contains(F_.sixteenerror1.msg, "Invalid typing--interface typevars misbound")
    @test contains(F_.sixteenerror1.msg, r"var\"#s\d+\" <: T")
    @test contains(F_.sixteenerror1.msg, "bound rule  :  Tuple{Symbol, Symbol} <: Tuple{Symbol, Integer}")

    @test typeof(F_.sixteenerror2) == ErrorException
    @test contains(F_.sixteenerror2.msg, "Invalid typing--interface typevars misbound")
    @test contains(F_.sixteenerror2.msg, r"var\"#s\d+\" <: Tuple{S, T}")
    @test contains(F_.sixteenerror2.msg, "bound rule  :  Tuple{A, B} <: Tuple{A, Integer}")
end

module G_
    using MultiInterface

    abstract type A{I,O,P} end
    abstract type B{I,O,P<:AbstractFloat} <: A{I,O,P} end
    abstract type C{I,O} <: B{I,O,Float64} end

    struct D <: C{Symbol, Int64} end
    struct E <: C{Tuple{Symbol, Int64}, String} end
    struct A1{I,O,P} <: A{I,O,P} end

    abstract type AA{I,O,Q,P<:A{I,O,Q}} <: A{I,O,Q} end
    abstract type BB{I,O,Q,P<:A{I,O,Q}} <: A{I,O,P} end
    struct AAA{I,O,Q,P} <: AA{I,O,Q,P} end
    struct CC{D} <: AA{Symbol, Int64, Float64, D} end
    struct DD{D} <: BB{Symbol, Int64, Float64, D} end

    Base.print_without_params(x::Type{A}) = false
    Base.print_without_params(x::Type{B}) = false
    Base.print_without_params(x::Type{C}) = false
    Base.print_without_params(x::Type{D}) = false
    Base.print_without_params(x::Type{E}) = false
    Base.print_without_params(x::Type{A1}) = false
    Base.print_without_params(x::Type{AA}) = false
    Base.print_without_params(x::Type{BB}) = false
    Base.print_without_params(x::Type{AAA}) = false
    Base.print_without_params(x::Type{CC}) = false
    Base.print_without_params(x::Type{DD}) = false

    @interface seventeen(a::A{I,O,P})::O where {I,O,P}
    @interface eighteen(a::B{I,O,P})::P where {I,O,P}
    @interface nineteen(a::AA{I,O,Q,P}, b::I, c::O, d::Q) where {I,O,Q,P<:A{I,O,Q}}

    @impl begin
        struct Seventeen_1 end
        function seventeen(a::B{I,Int64})::Int64 where I
            16
        end
        function seventeen(a::D)::Int64
            17
        end
    end

    @impl begin
        struct Eighteen_1 end
        function eighteen(a::D)::Float64
            18.1
        end
        function eighteen(a::E)::Float64
            18.2
        end
    end

    @impl begin
        struct Nineteen_1 end
        function nineteen(a::AAA{I,O,Q,P}, b::I, c::O, d::Q) where {I,O,Q,P<:A{I,O,Q}}
            (b,c,d)
        end

        function nineteen(a::CC{D}, b::Symbol, c::Int64, d::Float64)::Tuple where {D<:A{Symbol,Int64,Float64}}
            (b,c,d)
        end
    end

    seventeenerror1 = nothing
    try
        @impl begin
            struct Seventeen_2 end
            function seventeen(a::D)::Float64
                16.0
            end
        end
    catch e
        global seventeenerror1 = e
    end

    eightteenerror1 = nothing
    try
        @impl begin
            struct Eightteen_1 end
            function eighteen(a::D)::Int64
                18
            end
        end
    catch e
        global eightteenerror1 = e
    end

    nineteenerror1 = nothing
    try
        @impl begin
            function nineteen(a::AAA, b, c, d)
                (b,c,d)
            end
        end
    catch e
        global nineteenerror1 = e
    end
end

@testset "Complex subtying" begin
    struct BStruct <: G_.B{Symbol, Int64, Float64} end
    struct BStruct2 <: G_.B{Symbol, Symbol, Float64} end

    Base.print_without_params(x::Type{BStruct}) = false
    Base.print_without_params(x::Type{BStruct2}) = false

    @test G_.seventeen(BStruct()) == 16
    @test G_.seventeen(G_.D()) == 17
    @test_throws MethodError G_.seventeen(BStruct2())
    @test typeof(G_.seventeenerror1) == ErrorException

    @test contains(G_.seventeenerror1.msg, "Invalid typing--interface typevar bound to multiple types")
    @test contains(G_.seventeenerror1.msg, "typevar:  O")
    @test contains(G_.seventeenerror1.msg, "types  :  Any[:Int64, :Float64]")
    @test G_.eighteen(G_.D()) == 18.1
    @test G_.eighteen(G_.E()) == 18.2

    @test typeof(G_.eightteenerror1) == ErrorException
    @test contains(G_.eightteenerror1.msg, "Invalid typing--interface typevar bound to multiple types")
    @test contains(G_.eightteenerror1.msg, "typevar:  P")
    @test contains(G_.eightteenerror1.msg, "types  :  Any[:Float64, :Int64]")

    @test G_.nineteen(G_.AAA{Int64, Int64, Int64, G_.A1{Int64,Int64,Int64}}(), 1, 2, 3) == (1,2,3)
    @test G_.nineteen(G_.CC{G_.D}(), :a, 1, 1.0) == (:a, 1, 1.0)
    @test_throws MethodError G_.nineteen(G_.CC{G_.D}(), 1, 1, 1.0)

    @test typeof(G_.nineteenerror1) == ErrorException
    @test contains(G_.nineteenerror1.msg, "Invalid definition of nineteen")
    @test contains(G_.nineteenerror1.msg, "defined :  Type[Main.G_.AAA{I, O, Q, P} where {I, O, Q, P}, Any, Any, Any]")
    @test contains(G_.nineteenerror1.msg, "expected:  Type[Main.G_.AA{I, O, Q, P} where {I, O, Q, P<:Main.G_.A{I, O, Q}}, Any, Any, Any]")
    @test contains(G_.nineteenerror1.msg, "defined :  Union{Expr, Symbol}[:(AAA{I, O, Q, P} where {I, O, Q, P}), :Any, :Any, :Any]")
    @test contains(G_.nineteenerror1.msg, "expected:  Union{Expr, Symbol}[:(AA{I, O, Q, P} where P <: Main.G_.A{I, O, Q}), :(I where {I, O, Q, P <: A{I, O, Q}}), :(O where {I, O, Q, P <: A{I, O, Q}}), :(Q where {I, O, Q, P <: A{I, O, Q}})]")
    @test contains(G_.nineteenerror1.msg, "at      :  Bool[0, 1, 1, 1]")

end

module H_
    using MultiInterface

    abstract type A{I,O,P} end
    abstract type B{I,O,P} end
    struct AA{I,O,P} <: A{I,O,P} end
    struct BB{I,O,P} <: B{I,O,P} end

    Base.print_without_params(x::Type{A}) = false
    Base.print_without_params(x::Type{B}) = false
    Base.print_without_params(x::Type{AA}) = false
    Base.print_without_params(x::Type{BB}) = false

    @interface do_a(a::Union{A{Float64}, AA{Float64}})
    @interface do_b(a::Union{B{Float64}, BB{Float64}})
    @interface do_ab(ab::Union{A{Float64}, B{Float64}})
    @interface do_aabb(ab::Union{AA{Float64}, BB{Float64}})

    @interface do_t(t::Union{A{T}, B{T}}) where {T <: Real}

    @impl begin
        struct DoA1 end

        function do_a(a::A{Float64, Int64, Nothing})
            1.0
        end

        function do_a(a::Union{A{Float64, Float64}, A{Float64, Symbol}})
            2.0
        end

        function do_a(a::A{Float64})
            3.0
        end
    end

    doaerror1 = nothing
    try
        @impl begin
            struct DoA2 end

            function do_a(a::Union{A{Float64}, B{Int64}})
                3.0
            end
        end
    catch e
        global doaerror1 = e
    end

    doaerror2 = nothing
    try
        @impl begin
            struct DoA3 end

            function do_a(a::Union{A{Int64}, A{Float64}})
                4.0
            end
        end
    catch e
        global doaerror2 = e
    end

    @impl begin
        struct DoT1 end

        function do_t(t::AA{Int64})
            1.0
        end

        function do_t(t::B{Float64})
            2.0
        end

        function do_t(t::A{T}) where {T<:Real}
            3.0
        end

        function do_t(t::B{T}) where {T<:Real}
            4.0
        end
    end

    @impl begin
        struct DoT2 end

        function do_t(t::Union{A{Float64}, B{Int64}})
            5.0
        end
    end

    doterror1 = nothing
    try
        @impl begin
            struct DoT3 end

            function do_t(t::Union{A{Symbol}, B{Int64}})
                5.0
            end
        end
    catch e
        global doterror1 = e
    end
end

@testset "Unions" begin
    @test H_.do_a(H_.AA{Float64, Int64, Nothing}()) == 1.0
    @test H_.do_a(H_.AA{Float64, Float64, Nothing}()) == 2.0
    @test H_.do_a(H_.AA{Float64, Symbol, Nothing}()) == 2.0
    @test H_.do_a(H_.AA{Float64, Nothing, Nothing}()) == 3.0

    @test typeof(H_.doaerror1) == ErrorException
    @test contains(H_.doaerror1.msg, "Invalid definition of do_a")
    @test contains(H_.doaerror1.msg, "defined :  Union[Union{Main.H_.A{Float64, O, P} where {O, P}, Main.H_.B{Int64, O, P} where {O, P}}]")
    @test contains(H_.doaerror1.msg, "expected:  UnionAll[Main.H_.A{Float64, O, P} where {O, P}]")
    @test contains(H_.doaerror1.msg, "defined :  Union{Expr, Symbol}[:(Union{Main.H_.A{Float64, O, P} where {O, P}, Main.H_.B{Int64, O, P} where {O, P}})]")
    @test contains(H_.doaerror1.msg, "expected:  Union{Expr, Symbol}[:(A{Float64, O, P} where {O, P})]")
    @test contains(H_.doaerror1.msg, "at      :  Bool[0]")

    @test typeof(H_.doaerror2) == ErrorException
    @test contains(H_.doaerror2.msg, "Invalid definition of do_a")
    @test contains(H_.doaerror2.msg, "defined :  Union[Union{Main.H_.A{Float64, O, P} where {O, P}, Main.H_.A{Int64, O, P} where {O, P}}]")
    @test contains(H_.doaerror2.msg, "expected:  UnionAll[Main.H_.A{Float64, O, P} where {O, P}]")
    @test contains(H_.doaerror2.msg, "defined :  Union{Expr, Symbol}[:(Union{Main.H_.A{Float64, O, P} where {O, P}, Main.H_.A{Int64, O, P} where {O, P}})]")
    @test contains(H_.doaerror2.msg, "expected:  Union{Expr, Symbol}[:(A{Float64, O, P} where {O, P})]")
    @test contains(H_.doaerror2.msg, "at      :  Bool[0]")

    @test H_.do_t(H_.AA{Int64, Symbol, Nothing}()) == 1.0
    @test H_.do_t(H_.BB{Float64, Nothing, Nothing}()) == 2.0
    @test H_.do_t(H_.AA{Integer, Nothing, Nothing}()) == 3.0
    @test H_.do_t(H_.BB{Integer, Nothing, Nothing}()) == 4.0
    @test H_.do_t(H_.BB{Int64, Nothing, Nothing}()) == 5.0

    @test typeof(H_.doterror1) == ErrorException
    @test contains(H_.doterror1.msg, "Invalid definition of do_t")
    @test contains(H_.doterror1.msg, "defined :  Union[Union{Main.H_.A{Symbol, O, P} where {O, P}, Main.H_.B{Int64, O, P} where {O, P}}]")
    @test contains(H_.doterror1.msg, "expected:  UnionAll[Union{Main.H_.A{T, O, P} where {O, P}, Main.H_.B{T, O, P} where {O, P}} where T<:Real]")
    @test contains(H_.doterror1.msg, "defined :  Union{Expr, Symbol}[:(Union{Main.H_.A{Symbol, O, P} where {O, P}, Main.H_.B{Int64, O, P} where {O, P}})]")
    @test contains(H_.doterror1.msg, "expected:  Union{Expr, Symbol}[:(Union{Main.H_.A{T, O, P} where {O, P}, Main.H_.B{T, O, P} where {O, P}} where T <: Real)]")
    @test contains(H_.doterror1.msg, "at      :  Bool[0]")

end

module I_
    __OptVec{T} = Union{Vector{Union{}}, Vector{T}}
    __Opt{T} = Union{Nothing, T}

    using MultiInterface

    abstract type A{I,O,P} end
    abstract type AA{T,P} <: A{Tuple{}, T, P} end
    mutable struct AAA <: AA{Float64, Tuple{Float64, Float64}}
        params::Tuple{Float64, Float64}
    end

    Base.print_without_params(x::Type{A}) = false
    Base.print_without_params(x::Type{AA}) = false
    Base.print_without_params(x::Type{AAA}) = false

    @interface bounded_probs(sf::A{I,O},
                            range::__OptVec{<:O},
                            parranges::NTuple{N,Vector})::Tuple{Vector{<:AbstractFloat}, Vector{<:AbstractFloat}} where {I,O,P,N}

    @impl begin
        struct BoundedProbs_1 end

        function bounded_probs(n::AAA,
            range::Union{Vector{Union{}}, Vector{Float64}},
            parranges::NTuple{N,Vector})::Tuple{Vector{<:AbstractFloat}, Vector{<:AbstractFloat}} where {N}

            ([n.params[1]],[n.params[2]])
        end
    end
end

@testset "Unions with Aliasing" begin
    @test I_.bounded_probs(I_.AAA((1.0,2.0)), Union{}[], ([1.0],)) == ([1.0],[2.0])
end

module K_
    using MultiInterface
    struct X end

    @interface myx(x::K_.X)::Vector{<:K_.X}

    @impl begin
        struct MyX_1 end
        function myx(x::X)::Vector{<:K_.X}
            [x]
        end
    end

    abstract type SFunc{I,O} end
    Base.print_without_params(x::Type{<:SFunc}) = false

    abstract type Dist{T} <: SFunc{Tuple{}, T} end
    abstract type Score{I} <: SFunc{I,Nothing} end
    struct SC <: Score{Float64} end
    struct D <: Dist{Float64} end

    @interface send_lambda(sf::SFunc{I,O},
                        lambda::Score{<:O},
                        range::Union{Vector{Union{}}, Vector{<:O}},
                        parranges::NTuple{N,Vector},
                        incoming_pis::Union{Vector{Union{}}, Vector{<:Dist}},
                        parent_idx::Integer)::Score where {N,I,O}

    @impl begin
        function send_lambda(sf::Dist{T},
                            lambda::Score{<:T},
                            range::Union{Vector{Union{}}, Vector{<:T}},
                            parranges::NTuple{N,Vector},
                            incoming_pis::Union{Vector{Union{}}, Vector{<:Dist}},
                            parent_idx::Integer)::Score where {N,T}
            lambda
        end
    end

@testset "Arguments with module specification" begin
    @test K_.myx(K_.X()) == [K_.X()]
    @test K_.send_lambda(K_.D(), K_.SC(), Float64[], (), Union{}[], 1) == K_.SC()
end


text:
@testset "RSparseBitSet" begin
    @testset "RSparseBitSet(size, trailer)" begin
        trailer = SeaPearl.Trailer()
        n = 540

        sets = [
            SeaPearl.RSparseBitSet{UInt8}(n, trailer)
            SeaPearl.RSparseBitSet{UInt16}(n, trailer)
            SeaPearl.RSparseBitSet{UInt32}(n, trailer)
            SeaPearl.RSparseBitSet{UInt64}(n, trailer)
            SeaPearl.RSparseBitSet{UInt128}(n, trailer)
        ]

        strings = [join([bitstring(word.value) for word in set.words]) for set in sets]
        bools = [occursin(Regex("^1{$(n)}0*\$"), str) for str in strings]
        @test all(bools)
        @test all([set.limit.value == length(set.words) for set in sets])
    end

    @testset "Base.getindex" begin
        trailer = SeaPearl.Trailer()
        n = 540

        set = SeaPearl.RSparseBitSet(n, trailer)
        for word in set.words
            SeaPearl.setValue!(word, UInt64(0))
        end
        SeaPearl.setValue!(set.words[2], UInt64(1) << 63 | UInt64(1))

        @test !any([set[i] for i in 1:64])
        @test set[65]
        @test !any([set[i] for i in 66:127])
        @test set[128]
    end

    @testset "clearMask!" begin
        trailer = SeaPearl.Trailer()
        n = 540

        set = SeaPearl.RSparseBitSet(n, trailer)
        set.mask .= typemax(UInt64)
        SeaPearl.setValue!(set.limit, 6)

        SeaPearl.clearMask!(set)
        @test all(set.mask[1:6] .== UInt64(0))
        @test all(set.mask[7:end] .== typemax(UInt64))
    end

    @testset "reverseMask!" begin
        trailer = SeaPearl.Trailer()
        n = 540

        set = SeaPearl.RSparseBitSet(n, trailer)
        SeaPearl.clearMask!(set)
        set.mask[1:3] .= typemax(UInt64)
        SeaPearl.setValue!(set.limit, 6)

        SeaPearl.reverseMask!(set)
        @test all(set.mask[1:3] .== UInt64(0))
        @test all(set.mask[4:6] .== typemax(UInt64))
        @test all(set.mask[7:end] .== UInt64(0))
    end

    @testset "addToMask!" begin
        trailer = SeaPearl.Trailer()
        n = 20

        set = SeaPearl.RSparseBitSet{UInt8}(n, trailer)
        SeaPearl.clearMask!(set)
        x = [
            UInt8(1) << 7 | UInt8(1) << 2,
            UInt8(1) << 4,
            UInt8(1)
        ]
        y = [
            UInt8(1) << 3,
            UInt8(1) << 4,
            UInt8(1) << 2
        ]

        SeaPearl.addToMask!(set, x)
        @test all(set.mask .== x)

        SeaPearl.addToMask!(set, y)
        @test all(set.mask .== x .| y)
    end

    @testset "intersectWithMask!" begin
        trailer = SeaPearl.Trailer()
        n = 20

        set = SeaPearl.RSparseBitSet{UInt8}(n, trailer)
        SeaPearl.clearMask!(set)
        set.mask[1:end] = [
            UInt8(1) << 7 | UInt8(1) << 2,
            UInt8(1) << 4,
            UInt8(1)
        ]

        SeaPearl.intersectWithMask!(set)
        @test set[1]
        @test set[6]
        @test set[12]
        @test !any([set[i] for i in setdiff(1:24, [1,6,12])])
    end

    @testset "intersectIndex" begin
        trailer = SeaPearl.Trailer()
        n = 20

        set = SeaPearl.RSparseBitSet{UInt8}(n, trailer)
        SeaPearl.clearMask!(set)
        set.mask[1:end] = [
            UInt8(1) << 7 | UInt8(1) << 2,
            UInt8(1) << 4,
            UInt8(1) << 7
        ]
        SeaPearl.intersectWithMask!(set)

        x = [UInt8(1) << 2, typemax(UInt8), typemax(UInt8)]
        @test SeaPearl.intersectIndex(set, x) == 1

        y = [UInt8(1) << 3, UInt8(0), typemax(UInt8)]
        @test SeaPearl.intersectIndex(set, y) == 3

        z = [UInt8(1) << 3, UInt8(0), UInt8(0)]
        @test SeaPearl.intersectIndex(set, z) == -1

    end
end


text:
using DiffEqBayes, OrdinaryDiffEq, ParameterizedFunctions, RecursiveArrayTools
using Test, Distributions, SteadyStateDiffEq
println("Four parameter case")
f2 = @ode_def begin
    dx = a*x - b*x*y
    dy = -c*y + d*x*y
end a b c d
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob2 = ODEProblem(f2,u0,tspan,p)
sol = solve(prob2,Tsit5())
t = collect(range(1,stop=10,length=10))
randomized = VectorOfArray([(sol(t[i]) + .5randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)

priors = [truncated(Normal(1.5,0.5),0.5,2.5),truncated(Normal(1.2,0.5),0,2.0),
          truncated(Normal(1.0,0.5),1.0,4),truncated(Normal(1.0,0.5),0,2)]

@time bayesian_result = turing_inference(prob2,Tsit5(),t,data,priors;num_samples=1000,
                                   syms = [:a,:b,:c,:d])

@show bayesian_result

@test mean(get(bayesian_result,:a)[1]) ≈ 1.5 atol=3e-1
@test mean(get(bayesian_result,:b)[1]) ≈ 1.0 atol=3e-1
@test mean(get(bayesian_result,:c)[1]) ≈ 3.0 atol=3e-1
@test mean(get(bayesian_result,:d)[1]) ≈ 1.0 atol=3e-1

#=

# Ben's settings

parameters {
  real<lower = 0.5, upper = 2.5> alpha;
  real<lower = 0.0, upper = 2.0> beta;
  real<lower = 1.0, upper = 4.0> gamma;
  real<lower = 0.0, upper = 2.0> delta;
  real<lower = 0> sigma;
}
  alpha ~ normal(1.5, 0.5);
  beta ~ normal(1.2, 0.5);
  gamma ~ normal(3.0, 0.5);
  delta ~ normal(1.0, 0.5);

=#


text:
#= Grid types for OceanTurb.jl.

See fields.jl for the geometry of an OceanTurb.Field.
=#

import Base: eltype, length, size

# Staggered grid lengths and sizes for fields
cell_length(N) = N+2
face_length(N) = N+3
  cell_size(N) = (cell_length(N),)
  face_size(N) = (face_length(N),)

     height(g::Grid) = g.H
     length(g::Grid) = g.N
       size(g::Grid) = (g.N,)
cell_length(g::Grid) = cell_length(g.N)
face_length(g::Grid) = face_length(g.N)
  cell_size(g::Grid) = cell_size(g.N)
  face_size(g::Grid) = face_size(g.N)

height(m::AbstractModel) = height(m.grid)
length(m::AbstractModel) = length(m.grid)
  size(m::AbstractModel) = size(m.grid)

eltype(::Grid{T}) where T = T

"""
    arraytype(grid::Grid)

Return the array type corresponding to data
that lives on `grid`. Defaults to `Array`.
New data types (for example, grids that exist on GPUs) must
implement new array types.
"""
arraytype(grid::Grid{T}) where T = Array{T, 1} # default array type

"""
    UniformGrid([T], H, N)

Construct a 1D finite-volume grid with type `T` and array type `A`,
with `N` cells on the domain `z = [0, -H]`.
A `Grid` has two type parameters: an element type `T`,
and and a range type `R`. `T` defaults to Float64.
"""
struct UniformGrid{T, R} <: Grid{T, R}
    N  :: Int
    H  :: T
    Δc :: T
    Δf :: T
    zc :: R
    zf :: R
end

function UniformGrid(T, N::Int, H::Number)
    Δ = convert(T, H/N)
    H = convert(T, H)
    half_Δ = convert(T, 0.5Δ)

    zc = range(-H+half_Δ; length=N, stop=-half_Δ)
    zf = range(-H; length=N+1, stop=zero(T))

    UniformGrid(N, H, Δ, Δ, zc, zf)
end

# Defaults
UniformGrid(N, H) = UniformGrid(Float64, N, H)
UniformGrid(T=Float64; N=3, H=1) = UniformGrid(T, N, H)

"Return the cell spacing at index i."
@inline Δc(grid::UniformGrid, i) = grid.Δc

"Return the face spacing at index i."
@inline Δf(grid::UniformGrid, i) = grid.Δf


text:
using RecipesBase

function circle(x, y, r, θ)
    u = x + r * cos(θ) + 0.5
    v = y + r * sin(θ) + 0.5
    return u, v
end

circle(ap::CircularAperture, θ) = circle(ap.x, ap.y, ap.r, θ)

function ellipse(x, y, a, b, θ, t=0)
    u = x + a * cos(θ) * cosd(t) - b * sin(θ) * sind(t) + 0.5
    v = y + a * cos(θ) * sind(t) + b * sin(θ) * cosd(t) + 0.5
    return u, v
end

ellipse(ap::EllipticalAperture, θ) = ellipse(ap.x, ap.y, ap.a, ap.b, θ, ap.theta)

@recipe function f(c::CircularAperture, npoints = 100)
    seriestype := :shape
    aspect_ratio --> :equal
    seriescolor --> :red
    linecolor --> :match
    fillcolor --> nothing
    label --> ""

    t = range(0, 2pi, length = npoints)
    circle.((c,), t)
end

@recipe function f(c::CircularAnnulus, npoints = 100)
    seriestype := :shape
    aspect_ratio --> :equal
    seriescolor --> :red
    linecolor --> :match
    fillcolor --> nothing
    label --> ""

    t = range(0, 2pi, length = npoints)
    # outer ring
    @series begin
        circle.(c.x, c.y, c.r_out, t)
    end

    # inner ring
    @series begin
        circle.(c.x, c.y, c.r_in, t)
    end
end

@recipe function f(e::EllipticalAperture, npoints = 100)
    seriestype := :shape
    aspect_ratio --> :equal
    seriescolor --> :red
    linecolor --> :match
    fillcolor --> nothing
    label --> ""

    t = range(0, 2pi, length = npoints)
    ellipse.((e,), t)
end

@recipe function f(e::EllipticalAnnulus, npoints = 100)
    seriestype := :shape
    aspect_ratio --> :equal
    seriescolor --> :red
    linecolor --> :match
    fillcolor --> nothing
    label --> ""

    t = range(0, 2pi, length = npoints)

    # outer ring
    @series begin
        ellipse.(e.x, e.y, e.a_out, e.b_out, t, e.theta)
    end

    # inner ring
    @series begin
        ellipse.(e.x, e.y, e.a_in, e.b_in, t, e.theta)
    end
end

function rectangle(x, y, w, h, θ, t=0)
    half_w, half_h = w / 2, h / 2
    sinth, costh = sincos(θ)
    u = half_w * (abs(costh) * costh + abs(sinth) * sinth)
    v = half_h * (abs(costh) * costh - abs(sinth) * sinth)
    sint, cost = sincos(deg2rad(t))
    j = x + u * cost - v * sint + 0.5
    k = y + u * sint + v * cost + 0.5
    return j, k
end

rectangle(ap::RectangularAperture, θ) = rectangle(ap.x, ap.y, ap.w, ap.h, θ, ap.theta)

@recipe function f(ap::RectangularAperture, npoints = 100)
    seriestype := :shape
    aspect_ratio --> :equal
    seriescolor --> :red
    linecolor --> :match
    fillcolor --> nothing
    label --> ""

    t = range(0, 2π, length = npoints)
    return rectangle.((ap,), t)
end

@recipe function f(ap::RectangularAnnulus, npoints = 100)
    seriestype := :shape
    aspect_ratio --> :equal
    seriescolor --> :red
    linecolor --> :match
    fillcolor --> nothing
    label --> ""

    t = range(0, 2π, length = npoints)

    # outer box
    @series begin
        rectangle.(ap.x, ap.y, ap.w_out, ap.h_out, t, ap.theta)
    end

    # inner box
    @series begin
        rectangle.(ap.x, ap.y, ap.w_in, ap.h_in, t, ap.theta)
    end
end

@recipe function f(aps::AbstractVector{<:AbstractAperture})
    seriestype := :shape
    aspect_ratio --> :equal
    seriescolor --> :red
    linecolor --> :match
    fillcolor --> nothing
    label --> ""

    for ap in aps
        @series begin ap end
    end
end


text:
struct Mbplsr
    fm
    T::Matrix{Float64}
    scal
    weights::Vector{Float64}
end

"""
    mbplsr(X, Y, weights = ones(size(X[1], 1)); nlv, scal = nothing)
Multiblock PLSR (scaling and concatenation of blocks).
* `X` : List (vector) of blocks (matrices) of X-data.
    Each component of the list is a block.
* `Y` : Y-data.
* `weights` : Weights of the observations (rows).
* `nlv` : Nb. latent variables (LVs) to compute.
* `scal` : If `nothing`, each block is "autoscaled" (i.e. divided) by the
    the square root of the sum of the variances of each column of the block.
    Else, `scal` must be a vector defining the scaling values dividing the blocks.

Vector `weights` (row-weighting) is internally normalized to sum to 1.
See the help of `plskern` for details.
"""
function mbplsr(X, Y, weights = ones(size(X[1], 1)); nlv, scal = nothing)
    Y = ensure_mat(Y)
    weights = mweights(weights)
    res = blockscal(X, weights; scal = scal)
    zX = reduce(hcat, res.X)
    fm = plskern(zX, Y, weights; nlv = nlv)
    Mbplsr(fm, fm.T, res.scal, weights)
end

"""
    transform(object::Mbplsr, X; nlv = nothing)
Compute LVs ("scores" T) from a fitted model.
* `object` : The maximal fitted model.
* `X` : A list (vector) of blocks (matrices) of X-data for which LVs are computed.
* `nlv` : Nb. LVs to consider. If nothing, it is the maximum nb. LVs.
"""
function transform(object::Mbplsr, X; nlv = nothing)
    res = blockscal(X; scal = object.scal)
    zX = reduce(hcat, res.X)
    transform(object.fm, zX; nlv = nlv)
end

"""
    predict(object::Mbplsr, X; nlv = nothing)
Compute Y-predictions from a fitted model.
* `object` : The fitted model.
* `X` : A list (vector) of X-data for which predictions are computed.
* `nlv` : Nb. LVs, or collection of nb. LVs, to consider.
    If nothing, it is the maximum nb. LVs.
"""
function predict(object::Mbplsr, X; nlv = nothing)
    res = blockscal(X; scal = object.scal)
    zX = reduce(hcat, res.X)
    pred = predict(object.fm, zX; nlv = nlv).pred
    (pred = pred,)
end


text:
###################################################
## geo/timings/geotimings.jl
## abstract geotimings definition
###################################################

"""
    `GeoTimings`, abstract type that represents an estimator of timings in continuous setting
    must contain:
    - `trips`, reference to trip data that may be used
    must implement:
    - `estimateTime`, estimate time between two points, or for a trip
"""
abstract type GeoTimings end

function Base.show(io::IO, p::GeoTimings)
    typeName = split(string(typeof(p)),".")[end]
    println(io,"Continuous timings: $(typeName)")
end

"""
    `estimateTime`, estimate time of continuous trip
"""
estimateTime(gt::GeoTimings, o::Tuple{Float32, Float32}, d::Tuple{Float32, Float32})=
estimateTime(gt,o[1],o[2],d[1],d[2])
estimateTime(gt::GeoTimings, t::GeoTrip) =
estimateTime(gt, t.pLon, t.pLat, t.dLon, t.dLat)
estimateTime(gt::GeoTimings, tId::Int) =
estimateTime(gt, gt.trips[tId])
estimateTime(gt::GeoTimings, list::AbstractArray{Int,1}) =
Float64[estimateTime(gt, tId) for tId in list]


text:
function test_jaccard_distance()
  v = 6

  u = [1; 1; 1; 1; 1; 1]
  nu = [6; 0; 0; 0; 0; 0]
  ku = 1

  a = [1; 1; 1; 1; 2; 2]
  na = [4; 2; 0; 0; 0; 0]
  ka = 2

  b = [1; 2; 3; 1; 2; 3]
  nb = [2; 2; 2; 0; 0; 0]
  kb = 3

  c = [1; 1; 2; 3; 4; 3]
  nc = [2; 1; 2; 1; 0; 0]
  kc = 4

  d = [1; 1; 2; 3; 4; 5]
  nd = [2; 1; 1; 1; 1; 0]
  kd = 5

  t = [1; 2; 3; 4; 5; 6]
  nt = [1; 1; 1; 1; 1; 1]
  kt = 6

  @test Kpax3.jaccard(u, ku, u, ku, v) == 0.0
  @test Kpax3.jaccard_index_basic(u, u, v) == 1.0
  @test Kpax3.jaccard_index_table(u, ku, u, ku, v) == 1.0

  @test isapprox(Kpax3.jaccard(u, ku, a, ka, v), 0.5333333333333333)
  @test isapprox(Kpax3.jaccard(u, ku, b, kb, v), 0.8)
  @test isapprox(Kpax3.jaccard(u, ku, c, kc, v), 0.8666666666666667)
  @test isapprox(Kpax3.jaccard(u, ku, d, kd, v), 0.9333333333333333)
  @test isapprox(Kpax3.jaccard(u, ku, t, kt, v), 1.0)

  @test isapprox(Kpax3.jaccard_index_basic(u, a, v), 0.4666666666666667)
  @test isapprox(Kpax3.jaccard_index_basic(u, b, v), 0.2)
  @test isapprox(Kpax3.jaccard_index_basic(u, c, v), 0.1333333333333333)
  @test isapprox(Kpax3.jaccard_index_basic(u, d, v), 0.0666666666666667)
  @test isapprox(Kpax3.jaccard_index_basic(u, t, v), 0.0)

  @test isapprox(Kpax3.jaccard_index_table(u, ku, a, ka, v), 0.4666666666666667)
  @test isapprox(Kpax3.jaccard_index_table(u, ku, b, kb, v), 0.2)
  @test isapprox(Kpax3.jaccard_index_table(u, ku, c, kc, v), 0.1333333333333333)
  @test isapprox(Kpax3.jaccard_index_table(u, ku, d, kd, v), 0.0666666666666667)
  @test isapprox(Kpax3.jaccard_index_table(u, ku, t, kt, v), 0.0)

  @test isapprox(Kpax3.jaccardlower(na, ka, v), 0.5333333333333333)
  @test isapprox(Kpax3.jaccardlower(nb, kb, v), 0.8)
  @test isapprox(Kpax3.jaccardlower(nc, kc, v), 0.8666666666666667)
  @test isapprox(Kpax3.jaccardlower(nd, kd, v), 0.9333333333333333)
  @test isapprox(Kpax3.jaccardlower(nt, kt, v), 1.0)

  @test Kpax3.jaccard(a, ka, a, ka, v) == 0.0
  @test Kpax3.jaccard_index_basic(a, a, v) == 1.0
  @test Kpax3.jaccard_index_table(a, ka, a, ka, v) == 1.0

  @test isapprox(Kpax3.jaccard(a, ka, u, ku, v), 0.5333333333333333)
  @test isapprox(Kpax3.jaccard(a, ka, b, kb, v), 0.8888888888888889)
  @test isapprox(Kpax3.jaccard(a, ka, c, kc, v), 0.875)
  @test isapprox(Kpax3.jaccard(a, ka, d, kd, v), 0.8571428571428572)
  @test isapprox(Kpax3.jaccard(a, ka, t, kt, v), 1.0)

  @test isapprox(Kpax3.jaccard_index_basic(a, u, v), 0.4666666666666667)
  @test isapprox(Kpax3.jaccard_index_basic(a, b, v), 0.1111111111111111)
  @test isapprox(Kpax3.jaccard_index_basic(a, c, v), 0.125)
  @test isapprox(Kpax3.jaccard_index_basic(a, d, v), 0.1428571428571428)
  @test isapprox(Kpax3.jaccard_index_basic(a, t, v), 0.0)

  @test isapprox(Kpax3.jaccard_index_table(a, ka, u, ku, v), 0.4666666666666667)
  @test isapprox(Kpax3.jaccard_index_table(a, ka, b, kb, v), 0.1111111111111111)
  @test isapprox(Kpax3.jaccard_index_table(a, ka, c, kc, v), 0.125)
  @test isapprox(Kpax3.jaccard_index_table(a, ka, d, kd, v), 0.1428571428571428)
  @test isapprox(Kpax3.jaccard_index_table(a, ka, t, kt, v), 0.0)

  @test Kpax3.jaccard(b, kb, b, kb, v) == 0.0
  @test Kpax3.jaccard_index_basic(b, b, v) == 1.0
  @test Kpax3.jaccard_index_table(b, kb, b, kb, v) == 1.0

  @test isapprox(Kpax3.jaccard(b, kb, u, ku, v), 0.8)
  @test isapprox(Kpax3.jaccard(b, kb, a, ka, v), 0.8888888888888889)
  @test isapprox(Kpax3.jaccard(b, kb, c, kc, v), 1.0)
  @test isapprox(Kpax3.jaccard(b, kb, d, kd, v), 1.0)
  @test isapprox(Kpax3.jaccard(b, kb, t, kt, v), 1.0)

  @test isapprox(Kpax3.jaccard_index_basic(b, u, v), 0.2)
  @test isapprox(Kpax3.jaccard_index_basic(b, a, v), 0.1111111111111111)
  @test isapprox(Kpax3.jaccard_index_basic(b, c, v), 0.0)
  @test isapprox(Kpax3.jaccard_index_basic(b, d, v), 0.0)
  @test isapprox(Kpax3.jaccard_index_basic(b, t, v), 0.0)

  @test isapprox(Kpax3.jaccard_index_table(b, kb, u, ku, v), 0.2)
  @test isapprox(Kpax3.jaccard_index_table(b, kb, a, ka, v), 0.1111111111111111)
  @test isapprox(Kpax3.jaccard_index_table(b, kb, c, kc, v), 0.0)
  @test isapprox(Kpax3.jaccard_index_table(b, kb, d, kd, v), 0.0)
  @test isapprox(Kpax3.jaccard_index_table(b, kb, t, kt, v), 0.0)

  @test Kpax3.jaccard(c, kc, c, kc, v) == 0.0
  @test Kpax3.jaccard_index_basic(c, c, v) == 1.0
  @test Kpax3.jaccard_index_table(c, kc, c, kc, v) == 1.0

  @test isapprox(Kpax3.jaccard(c, kc, u, ku, v), 0.8666666666666667)
  @test isapprox(Kpax3.jaccard(c, kc, a, ka, v), 0.875)
  @test isapprox(Kpax3.jaccard(c, kc, b, kb, v), 1.0)
  @test isapprox(Kpax3.jaccard(c, kc, d, kd, v), 0.5)
  @test isapprox(Kpax3.jaccard(c, kc, t, kt, v), 1.0)

  @test isapprox(Kpax3.jaccard_index_basic(c, u, v), 0.1333333333333333)
  @test isapprox(Kpax3.jaccard_index_basic(c, a, v), 0.125)
  @test isapprox(Kpax3.jaccard_index_basic(c, b, v), 0.0)
  @test isapprox(Kpax3.jaccard_index_basic(c, d, v), 0.5)
  @test isapprox(Kpax3.jaccard_index_basic(c, t, v), 0.0)

  @test isapprox(Kpax3.jaccard_index_table(c, kc, u, ku, v), 0.1333333333333333)
  @test isapprox(Kpax3.jaccard_index_table(c, kc, a, ka, v), 0.125)
  @test isapprox(Kpax3.jaccard_index_table(c, kc, b, kb, v), 0.0)
  @test isapprox(Kpax3.jaccard_index_table(c, kc, d, kd, v), 0.5)
  @test isapprox(Kpax3.jaccard_index_table(c, kc, t, kt, v), 0.0)

  @test Kpax3.jaccard(d, kd, d, kd, v) == 0.0
  @test Kpax3.jaccard_index_basic(d, d, v) == 1.0
  @test Kpax3.jaccard_index_table(d, kd, d, kd, v) == 1.0

  @test isapprox(Kpax3.jaccard(d, kd, u, ku, v), 0.9333333333333333)
  @test isapprox(Kpax3.jaccard(d, kd, a, ka, v), 0.8571428571428572)
  @test isapprox(Kpax3.jaccard(d, kd, b, kb, v), 1.0)
  @test isapprox(Kpax3.jaccard(d, kd, c, kc, v), 0.5)
  @test isapprox(Kpax3.jaccard(d, kd, t, kt, v), 1.0)

  @test isapprox(Kpax3.jaccard_index_basic(d, u, v), 0.0666666666666667)
  @test isapprox(Kpax3.jaccard_index_basic(d, a, v), 0.1428571428571428)
  @test isapprox(Kpax3.jaccard_index_basic(d, b, v), 0.0)
  @test isapprox(Kpax3.jaccard_index_basic(d, c, v), 0.5)
  @test isapprox(Kpax3.jaccard_index_basic(d, t, v), 0.0)

  @test isapprox(Kpax3.jaccard_index_table(d, kd, u, ku, v), 0.0666666666666667)
  @test isapprox(Kpax3.jaccard_index_table(d, kd, a, ka, v), 0.1428571428571428)
  @test isapprox(Kpax3.jaccard_index_table(d, kd, b, kb, v), 0.0)
  @test isapprox(Kpax3.jaccard_index_table(d, kd, c, kc, v), 0.5)
  @test isapprox(Kpax3.jaccard_index_table(d, kd, t, kt, v), 0.0)

  @test Kpax3.jaccard(t, kt, t, kt, v) == 0.0
  @test Kpax3.jaccard_index_basic(t, t, v) == 1.0
  @test Kpax3.jaccard_index_table(t, kt, t, kt, v) == 1.0

  @test isapprox(Kpax3.jaccard(t, kt, u, ku, v), 1.0)
  @test isapprox(Kpax3.jaccard(t, kt, a, ka, v), 1.0)
  @test isapprox(Kpax3.jaccard(t, kt, b, kb, v), 1.0)
  @test isapprox(Kpax3.jaccard(t, kt, c, kc, v), 1.0)
  @test isapprox(Kpax3.jaccard(t, kt, d, kd, v), 1.0)

  @test isapprox(Kpax3.jaccard_index_basic(t, u, v), 0.0)
  @test isapprox(Kpax3.jaccard_index_basic(t, a, v), 0.0)
  @test isapprox(Kpax3.jaccard_index_basic(t, b, v), 0.0)
  @test isapprox(Kpax3.jaccard_index_basic(t, c, v), 0.0)
  @test isapprox(Kpax3.jaccard_index_basic(t, d, v), 0.0)

  @test isapprox(Kpax3.jaccard_index_table(t, kt, u, ku, v), 0.0)
  @test isapprox(Kpax3.jaccard_index_table(t, kt, a, ka, v), 0.0)
  @test isapprox(Kpax3.jaccard_index_table(t, kt, b, kb, v), 0.0)
  @test isapprox(Kpax3.jaccard_index_table(t, kt, c, kc, v), 0.0)
  @test isapprox(Kpax3.jaccard_index_table(t, kt, d, kd, v), 0.0)

  @test Kpax3.jaccard(u, ku, t, kt, v) <= (Kpax3.jaccard(u, ku, b, kb, v) + Kpax3.jaccard(b, kb, t, kt, v))
  @test Kpax3.jaccard(a, ka, b, kb, v) <= (Kpax3.jaccard(a, ka, c, kc, v) + Kpax3.jaccard(c, kc, b, kb, v))
  @test Kpax3.jaccard(b, kb, c, kc, v) <= (Kpax3.jaccard(b, kb, d, kd, v) + Kpax3.jaccard(d, kd, c, kc, v))
  @test Kpax3.jaccard(c, kc, d, kd, v) <= (Kpax3.jaccard(c, kc, a, ka, v) + Kpax3.jaccard(a, ka, d, kd, v))
  @test Kpax3.jaccard(a, ka, b, kb, v) <= (Kpax3.jaccard(a, ka, u, ku, v) + Kpax3.jaccard(u, ku, b, kb, v))
  @test Kpax3.jaccard(a, ka, b, kb, v) <= (Kpax3.jaccard(a, ka, t, kt, v) + Kpax3.jaccard(t, kt, b, kb, v))
  @test Kpax3.jaccard(a, ka, d, kd, v) <= (Kpax3.jaccard(a, ka, u, ku, v) + Kpax3.jaccard(u, ku, d, kd, v))
  @test Kpax3.jaccard(a, ka, d, kd, v) <= (Kpax3.jaccard(a, ka, t, kt, v) + Kpax3.jaccard(t, kt, d, kd, v))

  nothing
end

test_jaccard_distance()


text:
using JLD2, JSON, CSV, DataFrames, CodecZlib

edir = "data/SpanishLang/embeddings"
k = 33
common = CSV.read("data/SpanishLang/common-tokens-per-region.tsv.gz", DataFrame, delim='\t')
subset!(common, :n_regions => n -> n .> 10)
cclist = split(first(subset(common, :n_regions => n -> n .== 26).country_codes), ',')
@info sort!(cclist)

function encode_model(vocab, umapmodel)
    #  ivocab = Dict(token => i for (i, token) in enumerate(vocab))
    obj = Dict{String,Vector}()

    for (i, (idcol, distcol)) in enumerate(zip(eachcol(umapmodel.knns), eachcol(umapmodel.dists)))
        obj[vocab[i]] = [vocab[j] for j in idcol]
    end

    obj
end

#db = Dict()
for modelname in sort(glob(joinpath(edir, "umap+index-common-tokens-*.k=$k.jld2")))
    model_, vocab_ = load(modelname, "U2", "vocab")

    cc = replace(modelname, r".+tokens-" => "", r".k=.+" => "")
    v = encode_model(vocab_, model_)
    #v = db[cc]
    open("data/SpanishLang/common-voc-semantic-knn-$cc.json.gz", "w") do f
        gz = GzipCompressorStream(f)
        println(gz, JSON.json(v))
        close(gz)
    end
end

max_stars_repo_name: sadit/RegionalSpanish

max_issues_repo_name: sadit/RegionalSpanish

max_forks_repo_name: sadit/RegionalSpanish


text:
"""
    module TestLinearIntegerEquations - Test ModiaBase/src/LinearIntegerEquations.jl
"""
module TestLinearIntegerEquations

using Test
using ModiaBase

println("... Test LinearIntegerEquations.jl")

@testset "Test LinearIntegerEquations.jl" begin

    @testset "... Test Voltage source and resistor without ground" begin
        println("\n    --- Test Voltage source and resistor without ground")
        #=
            # Resistor
            1: R.v = R.p.v - R.n.v
            2: 0 = R.p.i + R.n.i
            3: R.v = R.R*R.p.i

            # Voltage source
            4: V.v = V.p.v - V.n.v
            5: 0   = V.p.i + V.n.i
            6: V.v = 10

            # Connect equations
            connect(V.p, R.p)
            connect(R.n, V.n)
            ->  7: V.p.v = R.p.v
                8: 0 = V.p.i + R.p.i
                9: R.n.v = V.n.v
            10: 0 = R.n.i + V.n.i

            # Variables:
            1: R.v
            2: R.p.v
            3: R.p.i
            4: R.n.v
            5: R.n.i
            6: V.v
            7: V.p.v
            8: V.p.i
            9: V.n.v
            10: V.n.i
        =#
        G = Vector{Int}[[1, 2, 4],
                        [3, 5],
                        [1, 3],
                        [6, 7, 9],
                        [8, 10],
                        [6],
                        [7, 2],
                        [8, 3],
                        [4, 9],
                        [5, 10]]

        eInt = Int[1,2,4,5,7,8,9,10]

        GcInt = Vector{Int}[[-1, 1, -1],
                            [1,1],
                            [-1,1,-1],
                            [1,1],
                            [-1,1],
                            [1,1],
                            [-1,1],
                            [1,1]]

        Avar = fill(0,10)

        vNames = ["R.v",
                "R.p.v",
                "R.p.i",
                "R.n.v",
                "R.n.i",
                "V.v",
                "V.p.v",
                "V.p.i",
                "V.n.v",
                "V.n.i"]

        name(v) = vNames[v]

        (vEliminated, vProperty, nvArbitrary, redundantEquations) = simplifyLinearIntegerEquations!(G, eInt, GcInt, Avar)

        printSimplifiedLinearIntegerEquations(G, eInt, GcInt, vEliminated, vProperty, nvArbitrary,
                                              redundantEquations, name, printTest=false)

        @test nvArbitrary == 1
        @test vEliminated == [9, 6, 2, 4, 10, 7, 8, 5]
        @test vProperty[vEliminated] == [0, 1, 1, 0, 3, 1, -3, -3]
        @test redundantEquations == [10]
        @test eInt  == [2, 5, 7, 8, 9, 1, 4, 10]
        @test G[eInt] == [Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[]]
        @test GcInt == [Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[]]

    end

    @testset "... Test Voltage source and resistor with ground" begin
        println("\n    --- Test Voltage source and resistor with ground")
        #=
            # Resistor
            1: R.v = R.p.v - R.n.v
            2: 0 = R.p.i + R.n.i
            3: R.v = R.R*R.p.i

            # Voltage source
            4: V.v = V.p.v - V.n.v
            5: 0   = V.p.i + V.n.i
            6: V.v = 10

            # Ground
            7: ground.p.v = 0

            # Connect equations
            connect(V.p, R.p)
            connect(R.n, V.n)
            connect(R.n, ground.p)
            ->   8: V.p.v = R.p.v
                 9: 0 = V.p.i + R.p.i
                10: R.n.v = V.n.v
                11: R.n.v = ground.p.v
                12: 0 = R.n.i + V.n.i + ground.p.i

            # Variables:
            1: R.v
            2: R.p.v
            3: R.p.i
            4: R.n.v
            5: R.n.i
            6: V.v
            7: V.p.v
            8: V.p.i
            9: V.n.v
            10: V.n.i
            11: ground.p.v
            12: ground.p.i
        =#
        G = Vector{Int}[[1, 2, 4],
                        [3, 5],
                        [1, 3],
                        [6, 7, 9],
                        [8, 10],
                        [6],
                        [11],
                        [7, 2],
                        [8, 3],
                        [4, 9],
                        [4,11],
                        [5, 10, 12]]

        eInt = Int[1,2,4,5,7,8,9,10,11,12]

        GcInt = Vector{Int}[[-1, 1, -1],
                            [1,1],
                            [-1,1,-1],
                            [1,1],
                            [1],
                            [-1,1],
                            [1,1],
                            [-1,1],
                            [-1,1],
                            [1,1,1]]

        Avar = fill(0,12)

        vNames = ["R.v",
                "R.p.v",
                "R.p.i",
                "R.n.v",
                "R.n.i",
                "V.v",
                "V.p.v",
                "V.p.i",
                "V.n.v",
                "V.n.i",
                "ground.p.v",
                "ground.p.i"]

        name(v) = vNames[v]

        (vEliminated, vProperty, nvArbitrary, redundantEquations) = simplifyLinearIntegerEquations!(G, eInt, GcInt, Avar)

        printSimplifiedLinearIntegerEquations(G, eInt, GcInt, vEliminated, vProperty, nvArbitrary,
                                              redundantEquations, name, printTest=false)

        @test nvArbitrary == 0
        @test vEliminated == [6, 12, 5, 10, 7, 2, 8, 9, 4, 11]
        @test vProperty[vEliminated] == [1, 0, -3, 3, 1, 1, -3, 0, 0, 0]
        @test redundantEquations == Int64[]
        @test eInt  == [7, 11, 10, 5, 1, 8, 9, 2, 12, 4]
        @test G[eInt] == [Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[]]
        @test GcInt == [Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[]]

    end

    @testset "... Test Voltage source and capacitor without ground" begin
        println("\n    --- Test Voltage source and capacitor without ground")
        #=
            # Capacitor
            1: C.v = C.p.v - C.n.v
            2: 0 = C.p.i + C.n.i
            3: C.C*der(C.v) = C.p.i

            # Voltage source
            4: V.v = V.p.v - V.n.v
            5: 0   = V.p.i + V.n.i
            6: V.v = 10

            # Connect equations
            connect(V.p, C.p)
            connect(C.n, V.n)
            ->  7: V.p.v = C.p.v
                8: 0 = V.p.i + C.p.i
                9: C.n.v = V.n.v
               10: 0 = C.n.i + V.n.i

            # Variables:
             1: C.v
             2: C.p.v
             3: C.p.i
             4: C.n.v
             5: C.n.i
             6: V.v
             7: V.p.v
             8: V.p.i
             9: V.n.v
            10: V.n.i
            11: der(C.v)
        =#
        G = Any[[1, 2, 4],
                [3, 5],
                [11, 3],
                [6, 7, 9],
                [8, 10],
                [6],
                [7, 2],
                [8, 3],
                [4, 9],
                [5, 10]]

        eInt = Int[1,2,4,5,7,8,9,10]

        GcInt = Any[[-1, 1, -1],
                    [1,1],
                    [-1,1,-1],
                    [1,1],
                    [-1,1],
                    [1,1],
                    [-1,1],
                    [1,1]]

        Avar = fill(0,10)
        pushfirst!(Avar,11)

        vNames = ["C.v",
                  "C.p.v",
                  "C.p.i",
                  "C.n.v",
                  "C.n.i",
                  "V.v",
                  "V.p.v",
                  "V.p.i",
                  "V.n.v",
                  "V.n.i",
                  "der(C.v)"]

        name(v) = vNames[v]

        (vEliminated, vProperty, nvArbitrary, redundantEquations) = simplifyLinearIntegerEquations!(G, eInt, GcInt, Avar)

        printSimplifiedLinearIntegerEquations(G, eInt, GcInt, vEliminated, vProperty, nvArbitrary,
                                              redundantEquations, name, printTest=false)

        @test nvArbitrary == 1
        @test vEliminated == [9, 6, 2, 4, 10, 7, 8, 5]
        @test vProperty[vEliminated] == [0, 1, 1, 0, 3, 1, -3, -3]
        @test redundantEquations == [10]
        @test eInt  == [2, 5, 7, 8, 9, 1, 4, 10]
        @test G[eInt] == Any[Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[]]
        @test GcInt == Any[Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[]]
    end

    @testset "... Test two inductances in series with parallel resistors without ground" begin
        println("\n    --- Test two inductances in series with parallel resistors without ground")
        #=
            # Inductor 1
            1: L1.v = L1.p.v - L1.n.v
            2: 0 = L1.p.i + L1.n.i
            3: L1.L*der(L1.p.i) = L1.v

            # Inductor 2
            4: L2.v = L2.p.v - L2.n.v
            5: 0 = L2.p.i + L2.n.i
            6: L2.L*der(L2.p.i) = L2.v

            # Resistor 1
            7: R1.v = R1.p.v - R1.n.v
            8:    0 = R1.p.i + R1.n.i
            9: R1.v = R1.R*R1.p.i

            # Resistor 2
            10: R2.v = R2.p.v - R2.n.v
            11:    0 = R2.p.i + R2.n.i
            12: R2.v = R2.R*R2.p.i

            # Voltage source
            13: V.v = V.p.v - V.n.v
            14: 0   = V.p.i + V.n.i
            15: V.v = 10

            # Connect equations
            connect(V.p, L1.p)
            connect(L1.n, R1.p)
            connect(L1.n, R2.p)
            connect(R1.n, L2.p)
            connect(R2.n, L2.p)
            connect(L2.n, V.n)
            ->  16: V.p.v = L1.p.v
                17: 0 = V.p.i + L1.p.i
                18: L1.n.v = R1.p.v
                19: L1.n.v = R2.p.v
                20: 0 = L1.n.i + R1.p.i + R2.p.i
                21: R1.n.v = L2.p.v
                22: R2.n.v = L2.p.v
                23: 0 = R1.n.i + R2.n.i + L2.p.i
                24: L2.n.v = V.n.v
                25: 0 = L2.n.i + V.n.i

            # Variables:
             1: L1.v
             2: L1.p.v
             3: L1.p.i
             4: L1.n.v
             5: L1.n.i
             6: L2.v
             7: L2.p.v
             8: L2.p.i
             9: L2.n.v
            10: L2.n.i
            11: R1.v
            12: R1.p.v
            13: R1.p.i
            14: R1.n.v
            15: R1.n.i
            16: R2.v
            17: R2.p.v
            18: R2.p.i
            19: R2.n.v
            20: R2.n.i
            21: V.v
            22: V.p.v
            23: V.p.i
            24: V.n.v
            25: V.n.i
            26: der(L1.p.i)
            27: der(L2.p.i)
        =#
        G = [[1, 2, 4],
             [3, 5],
             [26, 1],

             [6, 7, 9],
             [8,10],
             [27, 6],

             [11, 12, 14],
             [13, 15],
             [11, 13],

             [16, 17, 19],
             [18, 20],
             [16, 18],

             [21, 22, 24],
             [23, 25],
             [21],

             [22, 2],
             [23, 3],

             [4, 12],
             [4, 17],
             [5, 13, 18],

             [14, 7],
             [19, 7],
             [15, 20, 8],

             [9, 24],
             [10, 25]]

        eInt = Int[1,2,4,5,7,8,10,11,13,14,16,17,18,19,20,21,22,23,24,25]

        GcInt = [[-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]]

        Avar = fill(0,27)
        Avar[3] = 26
        Avar[8] = 27

        vNames = ["L1.v",
                  "L1.p.v",
                  "L1.p.i",
                  "L1.n.v",
                  "L1.n.i",
                  "L2.v",
                  "L2.p.v",
                  "L2.p.i",
                  "L2.n.v",
                  "L2.n.i",
                  "R1.v",
                  "R1.p.v",
                  "R1.p.i",
                  "R1.n.v",
                  "R1.n.i",
                  "R2.v",
                  "R2.p.v",
                  "R2.p.i",
                  "R2.n.v",
                  "R2.n.i",
                  "V.v",
                  "V.p.v",
                  "V.p.i",
                  "V.n.v",
                  "V.n.i",
                  "der(L1.p.i)",
                  "der(L2.p.i)"]

        var_name(v) = vNames[v]
        (vEliminated, vProperty, nvArbitrary, redundantEquations) = simplifyLinearIntegerEquations!(G, eInt, GcInt, Avar; log=true, var_name = var_name)

        printSimplifiedLinearIntegerEquations(G, eInt, GcInt, vEliminated, vProperty, nvArbitrary, redundantEquations, var_name, printTest=false)

        @test nvArbitrary == 1
        @test vEliminated == [7, 16, 12, 9, 19, 14, 17, 4, 25, 22, 23, 20, 15, 10, 5, 6]
        @test vProperty[vEliminated] == [0, 11, 11, 24, 0, 0, 11, 11, 3, 2, -3, -18, -13, -8, -3, -24]
        @test redundantEquations == [25]
        @test eInt  == [2, 5, 8, 11, 14, 16, 17, 18, 19, 21, 22, 24, 1, 13, 7, 10, 4, 23, 20, 25]
        @test G[eInt] == [Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], [2, 1, 11], [21, 2, 24], Int64[], Int64[], Int64[], [13, 8, 18], [3, 8], Int64[]]
        @test GcInt == [Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], [1, -1, -1], [-1, 1, -1], Int64[], Int64[], Int64[], [-1, 1, -1], [1, -1], Int64[]]

    end
end

end


text:
module DA_NlayerkernelTests

using Test
using Parameters, SpecialFunctions

abstract type DiffusionParameters end

include(joinpath(dirname(@__FILE__), "..", "..","src/forwardmodels/Diffusion Approximation/DAcylinder_layered.jl"))

## Test the derivations with original

@inline function α_coeff!(α, μa, D, sn)
    @inbounds for ind in 1:length(μa)
        α[ind] = sqrt(μa[ind] / D[ind] + sn^2)
    end
    return α
end

# greens function for 2 layer cylinder as presented by Liemert
function green_2layer(α, sn, μa, D, z, z0, zb, l, n)
    α = α_coeff!(α, μa, D, sn)

    β3 = sinh(α[2]*(l[2] + zb[2]))
    γ3 = cosh(α[2]*(l[2] + zb[2]))

    g1 = sinh(α[1]*(z0 + zb[1]))*sinh(α[1]*(z + zb[1])) / (D[1]*α[1]*exp(α[1]*(l[1] + zb[1])))
    g1 *= (D[1]*α[1]*n[1]^2*β3 - D[2]*α[2]*n[2]^2*γ3) / (D[1]*α[1]*n[1]^2*β3*cosh(α[1]*(l[1] + zb[1])) + (D[2]*α[2]*n[2]^2*γ3*sinh(α[1]*(l[1] + zb[1]))))

    g = (exp(-α[1]*abs(z - z0)) - exp(-α[1]*(z + z0 + 2zb[1]))) / (2*D[1]*α[1])

    return g + g1
end
# greens function for 3 layer cylinder as presented by Liemert
function green_3layer(α, sn, μa, D, z, z0, zb, l, n)
    α = α_coeff!(α, μa, D, sn)

    β3 = D[2]*α[2]*n[2]^2*cosh(α[2]*l[2])*sinh(α[3]*(l[3] + zb[2])) + D[3]*α[3]*n[3]^2*sinh(α[2]*l[2])*cosh(α[3]*(l[3] + zb[2]))
    γ3 = D[2]*α[2]*n[2]^2*sinh(α[2]*l[2])*sinh(α[3]*(l[3] + zb[2])) + D[3]*α[3]*n[3]^2*cosh(α[2]*l[2])*cosh(α[3]*(l[3] + zb[2]))

    g1 = sinh(α[1]*(z0 + zb[1]))*sinh(α[1]*(z + zb[1])) / (D[1]*α[1]*exp(α[1]*(l[1] + zb[1])))
    g1 *= (D[1]*α[1]*n[1]^2*β3 - D[2]*α[2]*n[2]^2*γ3) / (D[1]*α[1]*n[1]^2*β3*cosh(α[1]*(l[1] + zb[1])) + (D[2]*α[2]*n[2]^2*γ3*sinh(α[1]*(l[1] + zb[1]))))

    g = (exp(-α[1]*abs(z - z0)) - exp(-α[1]*(z + z0 + 2zb[1]))) / (2*D[1]*α[1])

    return g + g1
end

### test the 2 layer case
α = [1.2, 1.5]
sn = 2.2
D = [0.2, 0.5]
μa = [0.1, 0.14]
z0 = 0.2
zb = [0.2, 0.3]
l = [1.0, 2.0]
n = [1.0, 1.0]
n2 = @. D * n^2
z = 0.0

@test _green_Nlaycylin_top(sn, μa, D, z, z0, zb, l, n2, 2) ≈ green_2layer(α, sn, μa, D, z, z0, zb, l, n)

### test the 3 layer case
α = [1.2, 1.5, 1.3]
sn = 5.8
D = [0.2, 0.5, 0.7]
μa = [0.1, 0.14, 0.2]
z0 = 0.2
zb = [0.2, 0.3, 0.3]
l = [1.0, 2.0, 1.2]
n = [1.0, 1.0, 1.1]
n2 = @. D * n^2
z = 0.0

@test _green_Nlaycylin_top(sn, μa, D, z, z0, zb, l, n2, 3) ≈ green_3layer(α, sn, μa, D, z, z0, zb, l, n)

end # module


text:
using CausalForest
using Test

@testset "oob classification random forest test" begin

    features, labels = load_data("iris")
    features = float.(features)
    labels   = string.(labels)
    model = build_forest_oob(labels, features, -1, 100)
    predictions = apply_forest_oob(model)
    comparison = predictions .== labels
    len = length(labels)
    @test (sum(comparison)/len) > 0.9

end


text:
      # Tests the routine dh17 jacobian2

#include("../src/ttv.jl")

# Next, try computing three-body Keplerian Jacobian:

@testset "dh17" begin

#n = 8
n = 3
H = [3,1,1]
#n = 2
t0 = 7257.93115525
h  = 0.05
hbig = big(h)
tmax = 600.0
dlnq = big(1e-20)

#nstep = 8000
nstep = 500
#nstep = 1

elements = readdlm("elements.txt",',')
# Increase mass of inner planets:
elements[2,1] *= 100.
elements[3,1] *= 100.

m =zeros(n)
x0=zeros(3,n)
v0=zeros(3,n)

# Define which pairs will have impulse rather than -drift+Kepler:
pair = zeros(Bool,n,n)
# We want Keplerian between star & planets, and impulses between
# planets.  Impulse is indicated with 'true', -drift+Kepler with 'false':
#for i=2:n
#  pair[1,i] = false
#  # We don't need to define this, but let's anyways:
#  pair[i,1] = false
#end

# Initialize with identity matrix:
jac_step = Matrix{Float64}(I,7*n,7*n)
jac_error = zeros(7*n,7*n)

for k=1:n
  m[k] = elements[k,1]
end
m0 = copy(m)

init = ElementsIC(elements,H,t0)
x0,v0,_ = init_nbody(init)

# Tilt the orbits a bit:
x0[2,1] = 5e-1*sqrt(x0[1,1]^2+x0[3,1]^2)
x0[2,2] = -5e-1*sqrt(x0[1,2]^2+x0[3,2]^2)
x0[2,3] = -5e-1*sqrt(x0[1,2]^2+x0[3,2]^2)
v0[2,1] = 5e-1*sqrt(v0[1,1]^2+v0[3,1]^2)
v0[2,2] = -5e-1*sqrt(v0[1,2]^2+v0[3,2]^2)
v0[2,3] = -5e-1*sqrt(v0[1,2]^2+v0[3,2]^2)
xtest = copy(x0); vtest=copy(v0)
xbig = big.(x0); vbig = big.(v0)
# Take a single step (so that we aren't at initial coordinates):
@time for i=1:nstep; dh17!(x0,v0,h,m,n,pair); end
# Take a step with big precision:
ah18!(xbig,vbig,big(h),big.(m),n,pair)
# Take a single AH18 step:
@time for i=1:nstep; ah18!(xtest,vtest,h,m,n,pair);end
println("AH18 vs. DH17 x/v difference: ",x0-xtest,v0-vtest)

# Now, copy these to compute Jacobian (so that I don't step
# x0 & v0 forward in time):
x = copy(x0)
v = copy(v0)
m = copy(m0)
xerror = zeros(3,n); verror = zeros(3,n)
xbig = big.(x0)
vbig = big.(v0)
mbig = big.(m0)
xerr_big = zeros(BigFloat,3,n); verr_big=zeros(BigFloat,3,n);

# Compute jacobian exactly over nstep steps:
for istep=1:nstep
  dh17!(x,v,xerror,verror,h,m,n,jac_step,pair)
end
#println(typeof(h)," ",jac_step)
#read(STDIN,Char)

# The following lines have a Julia bug - jac_big gets
# misaligned or shifted when returned.  If I modify dh17! to output
# jac_step, then it works.
# Initialize with identity matrix:
## Compute jacobian exactly over nstep steps:
#for istep=1:nstep
#  jac_copy = dh17!(xbig,vbig,hbig,mbig,n,jac_big,pair)
#end
#println(typeof(hbig)," ",convert(Array{Float64,2},jac_copy))
#println("Comparing x & xbig: ",maximum(abs,x-xbig))
#println("Comparing v & vbig: ",maximum(abs,v-vbig))
#println("Comparing jac_step and jac_big: ",maxabs(jac_step-jac_copy))
#println("Comparing jac_step and jac_big: ",jac_step./convert(Array{Float64,2},jac_big))

# Test that both versions of dh17 give the same answer:
#xtest = copy(x0)
#vtest = copy(v0)
#m = copy(m0)
#for istep=1:nstep
#  dh17!(xtest,vtest,h,m,n,pair)
#end
#println("x/v difference: ",x-xtest,v-vtest)

# Now compute numerical derivatives:
jac_step_num = zeros(BigFloat,7*n,7*n)
# Vary the initial parameters of planet j:
for j=1:n
  # Vary the initial phase-space elements:
  for jj=1:3
  # Initial positions, velocities & masses:
    xm = big.(x0)
    vm = big.(v0)
    mm = big.(m0)
    dq = dlnq * xm[jj,j]
    if xm[jj,j] != 0.0
      xm[jj,j] -=  dq
    else
      dq = dlnq
      xm[jj,j] = -dq
    end
    for istep=1:nstep
      dh17!(xm,vm,hbig,mm,n,pair)
    end
    xp = big.(x0)
    vp = big.(v0)
    mp = big.(m0)
    dq = dlnq * xp[jj,j]
    if xp[jj,j] != 0.0
      xp[jj,j] +=  dq
    else
      dq = dlnq
      xp[jj,j] = dq
    end
    for istep=1:nstep
      dh17!(xp,vp,hbig,mp,n,pair)
    end
  # Now x & v are final positions & velocities after time step
    for i=1:n
      for k=1:3
        jac_step_num[(i-1)*7+  k,(j-1)*7+ jj] = .5*(xp[k,i]-xm[k,i])/dq
        jac_step_num[(i-1)*7+3+k,(j-1)*7+ jj] = .5*(vp[k,i]-vm[k,i])/dq
      end
    end
  # Next velocity derivatives:
    xm= big.(x0)
    vm= big.(v0)
    mm= big.(m0)
    dq = dlnq * vm[jj,j]
    if vm[jj,j] != 0.0
      vm[jj,j] -=  dq
    else
      dq = dlnq
      vm[jj,j] = -dq
    end
    for istep=1:nstep
      dh17!(xm,vm,hbig,mm,n,pair)
    end
    xp= big.(x0)
    vp= big.(v0)
    mp= big.(m0)
    dq = dlnq * vp[jj,j]
    if vp[jj,j] != 0.0
      vp[jj,j] +=  dq
    else
      dq = dlnq
      vp[jj,j] = dq
    end
    for istep=1:nstep
      dh17!(xp,vp,hbig,mp,n,pair)
    end
    for i=1:n
      for k=1:3
        jac_step_num[(i-1)*7+  k,(j-1)*7+3+jj] = .5*(xp[k,i]-xm[k,i])/dq
        jac_step_num[(i-1)*7+3+k,(j-1)*7+3+jj] = .5*(vp[k,i]-vm[k,i])/dq
      end
    end
  end
# Now vary mass of planet:
  xm= big.(x0)
  vm= big.(v0)
  mm= big.(m0)
  dq = mm[j]*dlnq
  mm[j] -= dq
  for istep=1:nstep
    dh17!(xm,vm,hbig,mm,n,pair)
  end
  xp= big.(x0)
  vp= big.(v0)
  mp= big.(m0)
  dq = mp[j]*dlnq
  mp[j] += dq
  for istep=1:nstep
    dh17!(xp,vp,hbig,mp,n,pair)
  end
  for i=1:n
    for k=1:3
      jac_step_num[(i-1)*7+  k,j*7] = .5*(xp[k,i]-xm[k,i])/dq
      jac_step_num[(i-1)*7+3+k,j*7] = .5*(vp[k,i]-vm[k,i])/dq
    end
  end
  # Mass unchanged -> identity
  jac_step_num[j*7,j*7] = 1.0
end

# Now, compare the results:
#println(jac_step)
#println(jac_step_num)

#for j=1:n
#  for i=1:7
#    for k=1:n
#      println(i," ",j," ",k," ",jac_step[(j-1)*7+i,(k-1)*7+1:k*7]," ",jac_step_num[(j-1)*7+i,(k-1)*7+1:k*7]," ",jac_step[(j-1)*7+i,(k-1)*7+1:7*k]./jac_step_num[(j-1)*7+i,(k-1)*7+1:7*k]-1.)
#    end
#  end
#end

jacmax = 0.0; jac_diff = 0.0
imax = 0; jmax = 0; kmax = 0; lmax = 0
for i=1:7, j=1:3, k=1:7, l=1:3
  if jac_step[(j-1)*7+i,(l-1)*7+k] != 0
    diff = abs(jac_step_num[(j-1)*7+i,(l-1)*7+k]/jac_step[(j-1)*7+i,(l-1)*7+k]-1.0)
    if diff > jacmax
      jac_diff = diff; imax = i; jmax = j; kmax = k; lmax = l; jacmax = jac_step[(j-1)*7+i,(l-1)*7+k]
    end
  end
end

println("Maximum fractional error: ",jac_diff," ",imax," ",jmax," ",kmax," ",lmax," ",jacmax)
#println(jac_step./jac_step_num)
println("Maximum error jac_step:   ",maximum(abs.(jac_step-jac_step_num)))
println("Maximum diff asinh(jac_step):   ",maximum(abs.(asinh.(jac_step)-asinh.(jac_step_num))))

# Compute dqdt:
# Set a short timestep so that particles will drift:
h = 0.001; hbig = big(h)
dqdt = zeros(7*n)
dqdt_num = zeros(BigFloat,7*n)
x = copy(x0)
v = copy(v0)
m = copy(m0)
xerror = zeros(3,n); verror = zeros(3,n)
dh17!(x,v,xerror,verror,h,m,n,dqdt,pair)
xm= big.(x0)
vm= big.(v0)
mm= big.(m0)
dq = hbig*dlnq
hbig -= dq
dh17!(xm,vm,hbig,mm,n,pair)
xp= big.(x0)
vp= big.(v0)
mp= big.(m0)
hbig += 2dq
dh17!(xp,vp,hbig,mp,n,pair)
for i=1:n, k=1:3
  dqdt_num[(i-1)*7+  k] = (xp[k,i]-xm[k,i])/dq/2
  dqdt_num[(i-1)*7+3+k] = (vp[k,i]-vm[k,i])/dq/2
end
dqdt_num = convert(Array{Float64,1},dqdt_num)
println("dqdt: ",dqdt," dqdt_num: ",dqdt_num," diff: ",dqdt-dqdt_num)
println("dqdt-dqdt_num: ",maxabs(dqdt-convert(Array{Float64,1},dqdt_num)))

#@test isapprox(jac_step,jac_step_num)
#@test isapprox(jac_step,jac_step_num;norm=maxabs)
@test isapprox(asinh.(jac_step),asinh.(jac_step_num);norm=maxabs)
@test isapprox(dqdt,dqdt_num;norm=maxabs)
end


text:
status,len,point,buffer = Grisu.bignumdtoa(1.0, Grisu.SHORTEST, 0, buffer)
@test "1" == trimrep(buffer)
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(1.0, Grisu.FIXED, 3, buffer)
@test 3 >= len - 1 - point
@test "1" == trimrep(buffer)
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(1.0, Grisu.PRECISION, 3, buffer)
@test 3 >= len - 1
@test "1" == trimrep(buffer)
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(1.5, Grisu.SHORTEST, 0, buffer)
@test "15" == trimrep(buffer)
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(1.5, Grisu.FIXED, 10, buffer)
@test 10 >= len - 1 - point
@test "15" == trimrep(buffer)
@test 1 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(1.5, Grisu.PRECISION, 10, buffer)
@test 10 >= len - 1
@test "15" == trimrep(buffer)
@test 1 == point
fill!(buffer,0);

min_double = 5e-324
status,len,point,buffer = Grisu.bignumdtoa(min_double, Grisu.SHORTEST, 0, buffer)
@test "5" == trimrep(buffer)
@test -323 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(min_double, Grisu.FIXED, 5, buffer)
@test 5 >= len - 1 - point
@test "" == trimrep(buffer)

status,len,point,buffer = Grisu.bignumdtoa(min_double, Grisu.PRECISION, 5, buffer)
@test 5 >= len - 1
@test "49407" == trimrep(buffer)
@test -323 == point
fill!(buffer,0);

max_double = 1.7976931348623157e308
status,len,point,buffer = Grisu.bignumdtoa(max_double, Grisu.SHORTEST, 0, buffer)
@test "17976931348623157" == trimrep(buffer)
@test 309 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(max_double, Grisu.PRECISION, 7, buffer)
@test 7 >= len - 1
@test "1797693" == trimrep(buffer)
@test 309 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(4294967272.0, Grisu.SHORTEST, 0, buffer)
@test "4294967272" == bytestring(pointer(buffer))
@test 10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(4294967272.0, Grisu.FIXED, 5, buffer)
@test "429496727200000" == bytestring(pointer(buffer))
@test 10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(4294967272.0, Grisu.PRECISION, 14, buffer)
@test 14 >= len - 1
@test "4294967272" == trimrep(buffer)
@test 10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(4.1855804968213567e298, Grisu.SHORTEST, 0,buffer)
@test "4185580496821357" == trimrep(buffer)
@test 299 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(4.1855804968213567e298, Grisu.PRECISION, 20,buffer)
@test 20 >= len - 1
@test "41855804968213567225" == trimrep(buffer)
@test 299 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(5.5626846462680035e-309, Grisu.SHORTEST, 0, buffer)
@test "5562684646268003" == trimrep(buffer)
@test -308 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(5.5626846462680035e-309, Grisu.PRECISION, 1, buffer)
@test 1 >= len - 1
@test "6" == trimrep(buffer)
@test -308 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(2147483648.0, Grisu.SHORTEST, 0, buffer)
@test "2147483648" == trimrep(buffer)
@test 10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(2147483648.0, Grisu.FIXED, 2, buffer)
@test 2 >= len - 1 - point
@test "2147483648" == trimrep(buffer)
@test 10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(2147483648.0, Grisu.PRECISION, 5, buffer)
@test 5 >= len - 1
@test "21475" == trimrep(buffer)
@test 10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(3.5844466002796428e+298, Grisu.SHORTEST, 0, buffer)
@test "35844466002796428" == trimrep(buffer)
@test 299 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(3.5844466002796428e+298, Grisu.PRECISION, 10, buffer)
@test 10 >= len - 1
@test "35844466" == trimrep(buffer)
@test 299 == point
fill!(buffer,0);

v = reinterpret(Float64,0x0010000000000000)
status,len,point,buffer = Grisu.bignumdtoa(v, Grisu.SHORTEST, 0, buffer)
@test "22250738585072014" == trimrep(buffer)
@test -307 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(v, Grisu.PRECISION, 20, buffer)
@test 20 >= len - 1
@test "22250738585072013831" == trimrep(buffer)
@test -307 == point
fill!(buffer,0);

v = reinterpret(Float64,0x000FFFFFFFFFFFFF)
status,len,point,buffer = Grisu.bignumdtoa(v, Grisu.SHORTEST, 0, buffer)
@test "2225073858507201" == trimrep(buffer)
@test -307 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(v, Grisu.PRECISION, 20, buffer)
@test 20 >= len - 1
@test "2225073858507200889" == trimrep(buffer)
@test -307 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(4128420500802942e-24, Grisu.SHORTEST, 0, buffer)
@test "4128420500802942" == trimrep(buffer)
@test -8 == point
fill!(buffer,0);

v = 3.9292015898194142585311918e-10;
status,len,point,buffer = Grisu.bignumdtoa(v, Grisu.SHORTEST, 0, buffer)
@test "39292015898194143" == trimrep(buffer)

v = 4194304.0;
status,len,point,buffer = Grisu.bignumdtoa(v, Grisu.FIXED, 5, buffer)
@test 5 >= len - 1 - point
@test "4194304" == trimrep(buffer)

v = 3.3161339052167390562200598e-237;
status,len,point,buffer = Grisu.bignumdtoa(v, Grisu.PRECISION, 19, buffer)
@test 19 >= len - 1
@test "3316133905216739056" == trimrep(buffer)
@test -236 == point
fill!(buffer,0);

v = 7.9885183916008099497815232e+191;
status,len,point,buffer = Grisu.bignumdtoa(v, Grisu.PRECISION, 4, buffer)
@test 4 >= len - 1
@test "7989" == trimrep(buffer)
@test 192 == point
fill!(buffer,0);

v = 1.0000000000000012800000000e+17;
status,len,point,buffer = Grisu.bignumdtoa(v, Grisu.FIXED, 1, buffer)
@test 1 >= len - 1 - point
@test "100000000000000128" == trimrep(buffer)
@test 18 == point
fill!(buffer,0);

min_float = float32(1e-45)
status,len,point,buffer = Grisu.bignumdtoa(min_float, Grisu.SHORTEST, 0, buffer)
@test "1" == trimrep(buffer)
@test -44 == point
fill!(buffer,0);

max_float = float32(3.4028234e38)
status,len,point,buffer = Grisu.bignumdtoa(max_float, Grisu.SHORTEST, 0, buffer)
@test "34028235" == trimrep(buffer)
@test 39 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(float32(4294967272.0), Grisu.SHORTEST, 0, buffer)
@test "42949673" == trimrep(buffer)
@test 10 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(float32(3.32306998946228968226e+35), Grisu.SHORTEST, 0, buffer)
@test "332307" == trimrep(buffer)
@test 36 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(float32(1.2341e-41), Grisu.SHORTEST, 0, buffer)
@test "12341" == trimrep(buffer)
@test -40 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(float32(3.3554432e7), Grisu.SHORTEST, 0, buffer)
@test "33554432" == trimrep(buffer)
@test 8 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(float32(3.26494756798464e14), Grisu.SHORTEST, 0, buffer)
@test "32649476" == trimrep(buffer)
@test 15 == point
fill!(buffer,0);

status,len,point,buffer = Grisu.bignumdtoa(float32(3.91132223637771935344e37), Grisu.SHORTEST, 0, buffer)
@test "39113222" == trimrep(buffer)
@test 38 == point
fill!(buffer,0);

v = reinterpret(Float32,0x00800000)
status,len,point,buffer = Grisu.bignumdtoa(v, Grisu.SHORTEST, 0, buffer)
@test "11754944" == trimrep(buffer)
@test -37 == point
fill!(buffer,0);

v = reinterpret(Float32,0x007FFFFF)
status,len,point,buffer = Grisu.bignumdtoa(v, Grisu.SHORTEST, 0, buffer)
@test "11754942" == trimrep(buffer)
@test -37 == point
fill!(buffer,0);


text:
@testset "Pipeline" begin
  @test true
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
# include("interpolate.jl") # not work
include("resample.jl")
include("aggregate.jl")

include("flipud.jl")

export aggregate, resample