text:
"""
    filter_matrix(f, domain, M, N)

Assemble discrete filtering matrix from a continuous filter `f`.
"""
function filter_matrix end

function filter_matrix(::IdentityFilter, domain::AbstractIntervalDomain, M, N; kwargs...)
    x = discretize(domain, M)
    ξ = discretize(domain, N)
    interpolation_matrix(x, ξ)
end

function filter_matrix(
    f::TopHatFilter,
    domain::ClosedIntervalDomain,
    M,
    N;
    degmax = 10,
    λ = 0,
)
    x = discretize(domain, M)
    ξ = discretize(domain, N)
    Δx = x[2] - x[1]
    ival_domain = (domain.left - Δx / 1000)..(domain.right + Δx / 1000)
    h = f.width
    W = spzeros(M + 1, N + 1)
    for m = 1:M+1
        # Filtering point
        xₘ = x[m]

        # Filter width at point
        hₘ = h(xₘ)

        # Indices of integration points inside domain
        ival = (xₘ ± hₘ) ∩ ival_domain
        ival_double = (xₘ ± 2hₘ) ∩ ival_domain
        n = ξ .∈ [ival_double]
        Nₘ = sum(n)

        # Polynomial basis
        P = min(degmax + 1, max(1, Nₘ == 2 ? 2 : 3 ≤ Nₘ < 9 ? 3 : floor(Int, √Nₘ)))
        ϕ = chebyshevt(0:P-1, ival_double)
        ϕ_int = integrate.(ϕ)

        # Assemble vandermonde matrix and moment vector
        Vₘ = spzeros(P, Nₘ)
        μₘ = zeros(P)
        for p = 1:P
            # Polynomials evaluated at integration points
            Vₘ[p, :] = ϕ[p].(ξ[n])

            # Polynomial moments around point
            μₘ[p] = (ϕ_int[p](ival.right) - ϕ_int[p](ival.left)) / (ival.right - ival.left)
        end

        # Fit weights
        # wₘ = ridge(Vₘ, μₘ, λ)
        # wₘ = Vₘ \ μₘ
        wₘ = nonneg_lsq(Vₘ, μₘ)

        # Store row of weights
        W[m, n] .= wₘ[:] ./ sum(wₘ)
    end
    dropzeros!(W)

    W
end

function filter_matrix(
    f::TopHatFilter,
    domain::PeriodicIntervalDomain,
    M,
    N;
    degmax = 10,
    λ = 0,
)
    x = discretize(domain, M)
    ξ = discretize(domain, N)
    Δx = x[2] - x[1]
    L = (domain.right - domain.left)
    h = f.width
    W = spzeros(M, N)
    for m = 1:M
        # Filtering point
        xₘ = x[m]

        # Filter width at point
        hₘ = h(xₘ)

        # Indices of integration points inside domain
        ival = xₘ ± hₘ
        n_left = ξ .+ L .∈ [ival]
        n_mid = ξ .∈ [ival]
        n_right = ξ .- L .∈ [ival]
        n = n_left .| n_mid .| n_right
        n = mapreduce(s -> circshift(n, s), .|, [-1, 0, 1])
        ξₙ = (ξ.+(xₘ.-ξ.>L/2)L.-(ξ.-xₘ.>L/2)L)[n]
        Nₘ = length(ξₙ)

        # Polynomial basis
        P = min(degmax + 1, Nₘ)
        # P = min(degmax + 1, max(1, Nₘ == 2 ? 2 : 3 ≤ Nₘ < 9 ? 3 : floor(Int, √Nₘ)))
        # ϕ = chebyshevt(0:P-1, ival)
        ϕ = chebyshevt(0:P-1, xₘ ± 3hₘ)
        ϕ_int = integrate.(ϕ)

        # Assemble vandermonde matrix and moment vector
        Vₘ = spzeros(P, Nₘ)
        zₘ = zeros(P)
        for p = 1:P
            # Polynomials evaluated at integration points
            Vₘ[p, :] = ϕ[p].(ξₙ)

            # Polynomial moments around point
            zₘ[p] = (ϕ_int[p](ival.right) - ϕ_int[p](ival.left)) / 2hₘ
        end

        # Fit weights
        wₘ = ridge(Vₘ, zₘ, λ)
        # wₘ = Vₘ \ zₘ
        # wₘ = nonneg_lsq(Vₘ, zₘ)

        # Store row of weights
        W[m, n] .= wₘ[:] ./ sum(wₘ)
    end
    dropzeros!(W)

    W
end

function filter_matrix(
    f::ConvolutionalFilter,
    domain::ClosedIntervalDomain,
    M,
    N;
    degmax = 10,
    λ = 0,
)
    x = discretize(domain, M)
    ξ = discretize(domain, N)
    Δx = x[2] - x[1]
    ival_domain = (domain.left - Δx / 1000)..(domain.right + Δx / 1000)
    h = f.width
    G = f.kernel
    W = spzeros(M + 1, N + 1)
    for m = 1:(M+1)
        # Filtering point
        xₘ = x[m]

        # Filter width at point
        hₘ = h(xₘ)

        # Indices of integration points inside domaain
        ival = (xₘ ± hₘ) ∩ ival_domain
        n = ξ .∈ [ival]
        ξₙ = ξ[n]
        Nₘ = length(ξₙ)

        # Polynomial basis
        P = min(degmax + 1, Nₘ)
        ϕ = chebyshevt(0:P-1, ival)

        # Normalized kernel Fun
        kern = Fun(ξ -> G(ξ - xₘ), ival)
        kern /= sum(kern)

        # Polynomial moments around point
        zₘ = sum.(kern .* ϕ)

        # Polynomials evaluated at integration points
        Vₘ = spzeros(P, Nₘ)
        for p = 1:P
            Vₘ[p, :] = ϕ[p].(ξₙ)
        end

        # Fit weights
        wₘ = ridge(Vₘ, zₘ, λ)
        # wₘ = Vₘ \ zₘ
        # wₘ = nonneg_lsq(Vₘ, zₘ)

        # Store row of weights
        W[m, n] .= wₘ[:] ./ sum(wₘ)
    end
    dropzeros!(W)

    W
end

function filter_matrix(
    f::ConvolutionalFilter,
    domain::PeriodicIntervalDomain,
    M,
    N;
    degmax = 10,
    λ = 0,
)
    x = discretize(domain, M)
    ξ = discretize(domain, N)
    L = (domain.right - domain.left)
    h = f.width
    G = f.kernel
    W = spzeros(M, N)
    for m = 1:M
        # Filtering point
        xₘ = x[m]

        # Filter width at point
        hₘ = h(xₘ)

        # Indices of integration points inside domain
        ival = (xₘ ± hₘ)
        n_left = ξ .+ L .∈ [ival]
        n_mid = ξ .∈ [ival]
        n_right = ξ .- L .∈ [ival]
        n = n_left .| n_mid .| n_right
        ξₙ = (ξ+L*n_left-L*n_right)[n]
        Nₘ = length(ξₙ)

        # Polynomial basis
        P = min(degmax + 1, Nₘ)
        # P = min(degmax + 1, max(1, floor(Int, √Nₘ)))
        ϕ = chebyshevt(0:P-1, ival)

        # Normalized kernel Fun
        kern = Fun(ξ -> G(ξ - xₘ), ival)
        kern /= sum(kern)

        # Assemble vandermonde matrix and moment vector
        Vₘ = spzeros(P, Nₘ)
        zₘ = zeros(P)
        # Fill in
        for p = 1:P
            # Polynomials evaluated at integration points
            Vₘ[p, :] = ϕ[p].(ξₙ)

            # Polynomial moments around point
            zₘ[p] = sum(kern * ϕ[p])
        end

        # Fit weights
        wₘ = ridge(Vₘ, zₘ, λ)
        # wₘ = Vₘ \ zₘ
        # wₘ = nonneg_lsq(Vₘ, zₘ)

        # Store row of weights
        W[m, n] .= wₘ[:] ./ sum(wₘ)
    end
    dropzeros!(W)

    W
end


text:
# This file is a part of EasyPkg, licensed under the MIT License (MIT).

import EasyPkg

try
    EasyPkg.run_all_tests("should_fail")
    error("Some tests were expected to fail, but didn't")
catch
    info("Ok, some tests were expected to fail")
end

include(joinpath("should_pass", "runtests.jl"))

ext: jl


text:
# Exception handler for RTDM interface

type RTDMError <: Exception
  errorcode :: Int
  tries     :: Int
end

function Base.showerror(io::IO, e::RTDMError)
  errortext = ""
  if e.errorcode > 0
    errortext = "targer replied with error code $(e.errorcode)"
  elseif e.errorcode == -1
    errortext = "malformed response"
  elseif e.errorcode == -2
    errortext = "error response crc fail"
  elseif e.errorcode == -3
    errortext = "communication link sanity check fail"
  elseif e.errorcode == -4
    errortext = "read fail"
  elseif e.errorcode == -5
    errortext = "write fail"
  end
  print(io, "RTDMError: ",errortext," (",e.tries," retries)")
end

function checkerrorcode(errorcode::Int, tries::Int)
  if errorcode != 0
    throw(RTDMError(errorcode,tries))
  end
  return nothing
end


text:
include("init.jl")

@testset "Bad Input message federate message" begin
    broker = createBroker(1)
    mFed1, fedinfo = createMessageFederate(1, "test")

    ept1 = h.helicsFederateRegisterGlobalEndpoint(mFed1, "ept1", "")
    @test_throws h.HELICSErrorRegistrationFailure h.helicsFederateRegisterGlobalEndpoint(mFed1, "ept1", "")

    h.helicsFederateEnterExecutingMode(mFed1)
    h.helicsEndpointSetDefaultDestination(ept1, "ept1")

    mess0 = h.helicsEndpointGetMessage(ept1)
    @test mess0.length == 0

    mess0 = h.helicsFederateGetMessage(mFed1)
    @test mess0.length == 0

    h.helicsEndpointSendMessage(ept1, mess0)

    h.helicsFederateRequestNextStep(mFed1)
    cnt = h.helicsEndpointPendingMessages(ept1)
    @test cnt == 1

    h.helicsFederateFinalize(mFed1)
    @test_throws h.HELICSErrorInvalidFunctionCall h.helicsEndpointSendMessage(ept1, mess0)

    destroyFederate(mFed1, fedinfo)
    destroyBroker(broker)

end

@testset "Bad Input filter test4" begin

    broker = createBroker(1)
    mFed1, fedinfo = createMessageFederate(1, "test")

    filt1 = h.helicsFederateRegisterCloningFilter(mFed1, "filt1")
    @test_throws h.HELICSErrorRegistrationFailure filt2 = h.helicsFederateRegisterCloningFilter(mFed1, "filt1")

    @test_throws h.HELICSErrorInvalidArgument h.helicsFilterSetString(filt1, "unknown", "string")

    h.helicsFederateRegisterGlobalEndpoint(mFed1, "ept1", "")

    h.helicsFilterAddDeliveryEndpoint(filt1, "ept1")
    h.helicsFilterAddSourceTarget(filt1, "ept1")
    h.helicsFilterAddDestinationTarget(filt1, "ept1")
    h.helicsFilterRemoveTarget(filt1, "ept1")

    @test_throws h.HELICSErrorInvalidArgument h.helicsFilterSet(filt1, "unknown", 10.0)
    h.helicsFederateFinalize(mFed1)

    destroyFederate(mFed1, fedinfo)
    destroyBroker(broker)

end

@testset "Bad Input filter core tests" begin

    broker = createBroker(1)
    mFed1, fedinfo = createMessageFederate(1, "test")

    cr = h.helicsFederateGetCoreObject(mFed1)

    filt1 = h.helicsCoreRegisterFilter(cr, h.HELICS_FILTER_TYPE_DELAY, "filt1")
    @test_throws h.HELICSErrorRegistrationFailure filt2 = h.helicsCoreRegisterFilter(cr, h.HELICS_FILTER_TYPE_DELAY, "filt1")
    h.helicsFilterSetOption(filt1, h.HELICS_HANDLE_OPTION_CONNECTION_OPTIONAL, true)
    @test h.helicsFilterGetOption(filt1, h.HELICS_HANDLE_OPTION_CONNECTION_OPTIONAL)
    h.helicsFederateFinalize(mFed1)
    h.helicsCoreDestroy(cr)

    destroyFederate(mFed1, fedinfo)
    destroyBroker(broker)
end

@testset "Bad Input type publication 2 tests" begin
    broker = createBroker(1)
    vFed1, fedinfo = createValueFederate(1, "test")

    pubid = h.helicsFederateRegisterGlobalTypePublication(vFed1, "pub1", "string", "")
    @test_throws h.HELICSErrorRegistrationFailure pubid2 = h.helicsFederateRegisterGlobalTypePublication(vFed1, "pub1", "string", "")

    @test_throws h.HELICSErrorInvalidArgument h.helicsFederateRegisterFromPublicationJSON(vFed1, "unknownfile.json")

    @test_throws h.HELICSErrorExternalType h.helicsFederateRegisterInterfaces(vFed1, "unknownfile.json")

    subid = h.helicsFederateRegisterTypeInput(vFed1, "inp1", "string", "")
    @test_throws h.HELICSErrorRegistrationFailure subid2 = h.helicsFederateRegisterTypeInput(vFed1, "inp1", "string", "")

    h.helicsInputAddTarget(subid, "pub1")

    h.helicsFederateSetTimeProperty(vFed1, h.HELICS_PROPERTY_TIME_PERIOD, 1.0)

    h.helicsFederateEnterExecutingModeIterativeAsync(vFed1, h.HELICS_ITERATION_REQUEST_NO_ITERATION)
    res = h.helicsFederateEnterExecutingModeIterativeComplete(vFed1)
    @test res == h.HELICS_ITERATION_RESULT_NEXT_STEP

    h.helicsPublicationPublishTime(pubid, 27.0)

    @test_throws h.HELICSErrorInvalidArgument h.helicsFederatePublishJSON(vFed1, "unknownfile.json")

    h.helicsFederateRequestNextStep(vFed1)
    str = h.helicsInputGetString(subid)
    @test str[1] == '2'
    @test str[2] == '7'
    h.helicsFederateClearUpdates(vFed1)

    h.helicsFederateFinalize(vFed1)

    @test_throws h.HELICSErrorInvalidFunctionCall h.helicsPublicationPublishRaw(pubid, str)
    @test_throws h.HELICSErrorInvalidFunctionCall h.helicsPublicationPublishString(pubid, str)
    @test_throws h.HELICSErrorInvalidFunctionCall h.helicsPublicationPublishInteger(pubid, 5)
    @test_throws h.HELICSErrorInvalidFunctionCall h.helicsPublicationPublishBoolean(pubid, true)
    @test_throws h.HELICSErrorInvalidFunctionCall h.helicsPublicationPublishDouble(pubid, 39.2)
    @test_throws h.HELICSErrorInvalidFunctionCall h.helicsPublicationPublishTime(pubid, 19.2)
    @test_throws h.HELICSErrorInvalidFunctionCall h.helicsPublicationPublishChar(pubid, 'a')

    @test_throws h.HELICSErrorInvalidFunctionCall h.helicsPublicationPublishComplex(pubid, 2.5 + -9.8im)

    @test_throws h.HELICSErrorInvalidFunctionCall h.helicsPublicationPublishVector(pubid, [1.3, 2.9])

    @test_throws h.HELICSErrorInvalidFunctionCall h.helicsPublicationPublishNamedPoint(pubid, "hello world", 2.0)

    destroyFederate(vFed1, fedinfo)
    destroyBroker(broker)
end

@testset "Bad Input tests raw tests" begin

    broker = createBroker(1)
    vFed1, fedinfo = createValueFederate(1, "test")

    pubid = h.helicsFederateRegisterPublication(vFed1, "pub1", h.HELICS_DATA_TYPE_RAW, "")

    subid = h.helicsFederateRegisterGlobalInput(vFed1, "inp1", h.HELICS_DATA_TYPE_RAW, "")

    h.helicsPublicationAddTarget(pubid, "inp1")

    h.helicsFederateSetTimeProperty(vFed1, h.HELICS_PROPERTY_TIME_PERIOD, 1.0)

    h.helicsFederateEnterExecutingMode(vFed1)

    h.helicsPublicationPublishDouble(pubid, 27.0)
    h.helicsFederateRequestNextStep(vFed1)
    h.helicsInputGetRawValue(subid)

    s = h.helicsInputGetString(subid)
    @test s == "0.000000"
    val = h.helicsInputGetComplexObject(subid)
    @test val == 0.0 + 0.0im

    h.helicsFederateFinalize(vFed1)

    t, res = h.helicsFederateRequestTimeIterative(vFed1, 1.0, h.HELICS_ITERATION_REQUEST_NO_ITERATION)
    @test res == h.HELICS_ITERATION_RESULT_HALTED

    destroyFederate(vFed1, fedinfo)
    destroyBroker(broker)
end

@testset "Bad Input duplicate publication and input pathways" begin

    broker = createBroker(1)
    vFed1, fedinfo = createValueFederate(1, "fed0")

    pubid = h.helicsFederateRegisterTypePublication(vFed1, "pub1", "string", "")

    @test_throws h.HELICSErrorRegistrationFailure pubid2 = h.helicsFederateRegisterTypePublication(vFed1, "pub1", "string", "")

    subid = h.helicsFederateRegisterGlobalTypeInput(vFed1, "inp1", "string", "")
    @test_throws h.HELICSErrorRegistrationFailure subid2 = h.helicsFederateRegisterGlobalTypeInput(vFed1, "inp1", "string", "")

    @test_throws h.HELICSErrorInvalidObject ept = h.helicsFederateRegisterEndpoint(vFed1, "ept1", "")

    h.helicsInputAddTarget(subid, "Testfed0/pub1")

    h.helicsFederateSetTimeProperty(vFed1, h.HELICS_PROPERTY_TIME_PERIOD, 1.0)

    h.helicsFederateEnterExecutingMode(vFed1)

    h.helicsPublicationPublishDouble(pubid, 27.0)
    h.helicsFederateRequestNextStep(vFed1)
    str = h.helicsInputGetString(subid)
    @test str[1] == '2'
    @test str[2] == '7'

    messages = h.helicsFederatePendingMessages(vFed1)
    @test messages == 0

    h.helicsFederateFinalize(vFed1)

    destroyFederate(vFed1, fedinfo)
    destroyBroker(broker)

end

@testset "Bad input init error" begin
    broker = createBroker(1)
    vFed1, fedinfo = createValueFederate(1, "fed0")

    # register the publications

    # the types here don't match which causes an error when initializing the federation
    h.helicsFederateRegisterGlobalTypePublication(vFed1, "pub1", "custom1", "")

    subid = h.helicsFederateRegisterTypeInput(vFed1, "inp1", "custom2", "")
    k1 = h.helicsInputGetKey(subid)

    # check some other calls
    inp2 = h.helicsFederateGetInput(vFed1, "inp1")
    k2 = h.helicsInputGetKey(inp2)
    @test k1 == k2

    inp3 = h.helicsFederateGetInputByIndex(vFed1, 0)
    k3 = h.helicsInputGetKey(inp3)
    @test k1 == k3

    h.helicsInputAddTarget(subid, "pub1")

    h.helicsFederateSetTimeProperty(vFed1, h.HELICS_PROPERTY_TIME_PERIOD, 1.0)

    # unknown publication
    @test_throws h.HELICSErrorInvalidArgument pub3 = h.helicsFederateGetPublication(vFed1, "unknown")

    # error in this call from the mismatch
    @test_throws h.HELICSErrorConnectionFailure h.helicsFederateEnterInitializingMode(vFed1)

    @test_throws h.HELICSErrorInvalidFunctionCall h.helicsFederateRequestTimeAdvance(vFed1, 0.1)

    # unknown input
    @test_throws h.HELICSErrorInvalidArgument inp4 = h.helicsFederateGetInput(vFed1, "unknown")

    # invalid input index
    # TODO: does this test segfault some times?
    @test_throws h.HELICSErrorInvalidArgument inp5 = h.helicsFederateGetInputByIndex(vFed1, 4)

    h.helicsFederateFinalize(vFed1)

    destroyFederate(vFed1, fedinfo)
    destroyBroker(broker)

end

@testset "Bad Inputs input tests" begin

    broker = createBroker(1)
    vFed1, fedinfo = createValueFederate(1, "fed0")

    # register the publications

    pubid = h.helicsFederateRegisterGlobalPublication(vFed1, "pub1", h.HELICS_DATA_TYPE_DOUBLE, "")

    subid = h.helicsFederateRegisterInput(vFed1, "inp1", h.HELICS_DATA_TYPE_DOUBLE, "")
    @test_throws h.HELICSErrorRegistrationFailure subid2 = h.helicsFederateRegisterInput(vFed1, "inp1", h.HELICS_DATA_TYPE_DOUBLE, "")

    h.helicsInputAddTarget(subid, "pub1")

    vf2 = h.helicsFederateClone(vFed1)
    @test h.helicsFederateGetName(vFed1) == h.helicsFederateGetName(vf2)

    h.helicsFederateSetTimeProperty(vFed1, h.HELICS_PROPERTY_TIME_PERIOD, 1.0)

    @test_throws h.HELICSErrorInvalidObject ept = h.helicsFederateRegisterEndpoint(vFed1, "ept1", "")

    h.helicsFederateEnterInitializingMode(vFed1)

    h.helicsPublicationPublishDouble(pubid, 27.0)

    comp = h.helicsFederateEnterExecutingModeIterative(vFed1, h.HELICS_ITERATION_REQUEST_FORCE_ITERATION)
    @test comp == h.HELICS_ITERATION_RESULT_ITERATING
    val = h.helicsInputGetDouble(subid)
    @test val == 27.0
    valt = h.helicsInputGetTime(subid)
    @test valt == 27.0

    comp = h.helicsFederateEnterExecutingModeIterative(vFed1, h.HELICS_ITERATION_REQUEST_ITERATE_IF_NEEDED)

    @test comp == h.HELICS_ITERATION_RESULT_NEXT_STEP

    val2 = h.helicsInputGetDouble(subid)
    @test val2 == val
    #expect error entering initializing Mode again
    @test_throws h.HELICSErrorInvalidFunctionCall h.helicsFederateEnterInitializingMode(vFed1)

    #expect error entering initializing Mode again
    @test_throws h.HELICSErrorInvalidFunctionCall h.helicsFederateEnterInitializingModeAsync(vFed1)

    #expect error entering initializing Mode again
    @test_throws h.HELICSErrorInvalidFunctionCall h.helicsFederateEnterInitializingModeComplete(vFed1)

    h.helicsFederateFinalize(vFed1)

    destroyFederate(vFed1, fedinfo)
    destroyBroker(broker)

end

@testset "Bad Inputs Core Link" begin
    broker = createBroker(1)
    vFed1, fedinfo = createValueFederate(1, "fed0")

    # register the publications

    h.helicsFederateRegisterGlobalTypePublication(vFed1, "pub1", "custom1", "")

    h.helicsFederateRegisterTypeInput(vFed1, "inp1", "custom2", "")

    fed2 = h.helicsGetFederateByName(h.helicsFederateGetName(vFed1))
    @test h.helicsFederateGetName(fed2) == h.helicsFederateGetName(vFed1)

    @test_throws h.HELICSErrorInvalidArgument fed3 = h.helicsGetFederateByName("fed_unknown")

    cr = h.helicsFederateGetCoreObject(vFed1)
    h.helicsCoreDataLink(cr, "pub1", "fed0/inp1")

    @test_throws h.HELICSErrorInvalidArgument h.helicsCoreMakeConnections(cr, "unknownfile.json")

    # TODO: This test should throw an error
    # @test_throws h.HELICSErrorInvalidArgument h.helicsCoreDataLink(cr, "pub1", "")
    @test_broken false

    cr2 = h.helicsCoreClone(cr)
    @test h.helicsCoreGetAddress(cr2) == h.helicsCoreGetAddress(cr)

    # TODO: this should error as well
    h.helicsFederateEnterExecutingMode(vFed1)
    @test_broken false

    h.helicsFederateFinalize(vFed1)

    destroyFederate(vFed1, fedinfo)
    destroyBroker(broker)
end

@testset "Bad Inputs broker link" begin

    broker = createBroker(1)
    vFed1, fedinfo = createValueFederate(1, "fed0")

    # register the publications

    h.helicsFederateRegisterGlobalTypePublication(vFed1, "pub1", "custom1", "")

    h.helicsFederateRegisterTypeInput(vFed1, "inp1", "custom2", "")

    br = h.helicsBrokerClone(broker)

    h.helicsBrokerDataLink(br, "pub1", "Testfed0/inp1")

    # TODO: This test should throw an error
    # @test_throws h.HELICSErrorInvalidArgument h.helicsBrokerDataLink(br, "pub1", "")
    @test_broken false

    @test_throws h.HELICSErrorInvalidArgument h.helicsBrokerMakeConnections(br, "unknownfile.json")

    @test_throws h.HELICSErrorConnectionFailure h.helicsFederateEnterExecutingMode(vFed1)

    h.helicsFederateFinalize(vFed1)

    destroyFederate(vFed1, fedinfo)

    h.helicsBrokerWaitForDisconnect(broker, 200)

    destroyBroker(broker)

end

@testset "Bad Inputs frees" begin
    broker = createBroker(1)
    vFed1, fedinfo = createValueFederate(1, "fed0")

    fi = h.helicsCreateFederateInfo()
    h.helicsFederateInfoSetBroker(fi, "broker test")
    h.helicsFederateEnterInitializingMode(vFed1)
    h.helicsFederateFinalize(vFed1)

    h.helicsFederateInfoSetBrokerPort(fi, 8929)
    h.helicsFederateInfoSetLocalPort(fi, "8229")

    h.helicsFederateInfoFree(fi)
    h.helicsFederateFree(vFed1)

    @test_throws h.HELICSErrorInvalidObject destroyFederate(vFed1, fedinfo)
    destroyBroker(broker)
end

@testset "Bad Inputs init error 5" begin

    broker = createBroker(1)
    vFed1, fedinfo = createValueFederate(1, "fed0")
    # register the publications

    h.helicsFederateInfoSetSeparator(fedinfo, '-')
    h.helicsFederateSetSeparator(vFed1, '-')

    h.helicsFederateRegisterGlobalTypePublication(vFed1, "pub1", "custom1", "");

    subid = h.helicsFederateRegisterTypeInput(vFed1, "inp1", "custom2", "");

    h.helicsInputAddTarget(subid, "pub1");

    h.helicsFederateSetTimeProperty(vFed1, h.HELICS_PROPERTY_TIME_PERIOD, 1.0);

    @test_throws h.HELICSErrorConnectionFailure resIt = h.helicsFederateEnterExecutingModeIterative(vFed1, h.HELICS_ITERATION_REQUEST_NO_ITERATION);

    @test_throws h.HELICSErrorInvalidFunctionCall h.helicsFederateRequestTimeIterativeAsync(vFed1, 1.0, h.HELICS_ITERATION_REQUEST_NO_ITERATION);

    @test_throws h.HELICSErrorInvalidFunctionCall res = h.helicsFederateRequestTimeIterativeComplete(vFed1)

    h.helicsFederateFinalize(vFed1);

    destroyFederate(vFed1, fedinfo)

    destroyBroker(broker)

end

@testset "Bad Inputs misc tests" begin

    @test h.helicsGetPropertyIndex("") == -1
    @test h.helicsGetPropertyIndex("not_a_property") == -1

    @test h.helicsGetOptionIndex("") == -1
    @test h.helicsGetOptionIndex("not_a_property") == -1

end


text:
# This code is taken from DataFrames.jl/src/other/utils.jl

if VERSION >= v"1.3"
    using Base.Threads: @spawn
else
    # This is the definition of @async in Base
    macro spawn(expr)
        thunk = esc(:(()->($expr)))
        var = esc(Base.sync_varname)
        quote
            local task = Task($thunk)
            if $(Expr(:isdefined, var))
                push!($var, task)
            end
            schedule(task)
        end
    end
end

# Compute chunks of indices, each with at least `basesize` entries
# This method ensures balanced sizes by avoiding a small last chunk
function split_indices(len::Integer, basesize::Integer)
    len′ = Int64(len) # Avoid overflow on 32-bit machines
    @assert len′ > 0
    @assert basesize > 0
    np = Int64(max(1, len ÷ basesize))
    return split_to_chunks(len′, np)
end

function split_to_chunks(len::Integer, np::Integer)
    len′ = Int64(len) # Avoid overflow on 32-bit machines
    np′ = Int64(np)
    @assert len′ > 0
    @assert 0 < np′ <= len′
    return (Int(1 + ((i - 1) * len′) ÷ np):Int((i * len′) ÷ np) for i in 1:np)
end

if VERSION >= v"1.4"
    function _spawn_for_chunks_helper(iter, lbody, basesize)
        lidx = iter.args[1]
        range = iter.args[2]
        quote
            let x = $(esc(range)), basesize = $(esc(basesize))
                @assert firstindex(x) == 1

                nt = Threads.nthreads()
                len = length(x)
                if nt > 1 && len > basesize
                    tasks = [Threads.@spawn begin
                                 for i in p
                                     local $(esc(lidx)) = @inbounds x[i]
                                     $(esc(lbody))
                                 end
                             end
                             for p in split_indices(len, basesize)]
                    foreach(wait, tasks)
                else
                    for i in eachindex(x)
                        local $(esc(lidx)) = @inbounds x[i]
                        $(esc(lbody))
                    end
                end
            end
            nothing
        end
    end
else
    function _spawn_for_chunks_helper(iter, lbody, basesize)
        lidx = iter.args[1]
        range = iter.args[2]
        quote
            let x = $(esc(range))
                for i in eachindex(x)
                    local $(esc(lidx)) = @inbounds x[i]
                    $(esc(lbody))
                end
            end
            nothing
        end
    end
end

"""
    @spawn_for_chunks basesize for i in range ... end
Parallelize a `for` loop by spawning separate tasks
iterating each over a chunk of at least `basesize` elements
in `range`.
A number of task higher than `Threads.nthreads()` may be spawned,
since that can allow for a more efficient load balancing in case
some threads are busy (nested parallelism).
"""
macro spawn_for_chunks(basesize, ex)
    if !(isa(ex, Expr) && ex.head === :for)
        throw(ArgumentError("@spawn_for_chunks requires a `for` loop expression"))
    end
    if !(ex.args[1] isa Expr && ex.args[1].head === :(=))
        throw(ArgumentError("nested outer loops are not currently supported by @spawn_for_chunks"))
    end
    return _spawn_for_chunks_helper(ex.args[1], ex.args[2], basesize)
end

max_stars_repo_name: matthieugomez/GroupedArrays

max_issues_repo_name: matthieugomez/GroupedArrays


text:
"""
ARP fileformat

Defines the geographical location of airports used in the assignment process

Example:
EGSX 3103.267 9.350 EGTT_FIR
EGSY 3203.650 -83.300 EGTT_FIR
EGSZ 3620.000 -102.000 EGPX_FIR
EGTA 3106.883 -56.550 EGTT_FIR
EGTB 3096.700 -48.483 EGTT_FIR

See EUROCONTROL NEST Manual sectino 9.7.7 for Arp fileformat description
"""
module Arp

export read

using CSV

const fileformat = Dict(1=>String, 2=>Float64, 3=>Float64, 4=>String)

const header_format = ["AIRPORT", "LAT_DEG", "LON_DEG", "FIR"]

function read(file)
    df = CSV.read(file, types=fileformat, header=header_format,
    copycols=true)
    reformat!(df)
    return df
end

function reformat!(df)
    df[:,:LAT_DEG] = df[:,:LAT_DEG] / 60.0
    df[:,:LON_DEG] = df[:,:LON_DEG] / 60.0
end

end # module


text:
# Quadrilateral in 3D defined by its 4 vertices.

# NOTE: Quadrilaterals are assumed to be convex and planar (all points in some plane, not
# necessarily xy, yx, etc.)!
# Quadrilaterals must be convex to be valid finite elements. See link below
# https://math.stackexchange.com/questions/2430691/jacobian-determinant-for-bi-linear-quadrilaterals
struct Quadrilateral_3D{T <: AbstractFloat}
    # Counter clockwise order
    points::NTuple{4, Point_3D{T}}
end

# Constructors
# -------------------------------------------------------------------------------------------------
Quadrilateral_3D(p₁::Point_3D{T},
                 p₂::Point_3D{T},
                 p₃::Point_3D{T},
                 p₄::Point_3D{T}) where {T <: AbstractFloat} = Quadrilateral_3D((p₁, p₂, p₃, p₄))

# Methods
# -------------------------------------------------------------------------------------------------
# Interpolation
function (quad::Quadrilateral_3D{T})(r::R, s::S) where {T <: AbstractFloat,
                                                        R <: Real,
                                                        S <: Real}
    # See The Visualization Toolkit: An Object-Oriented Approach to 3D Graphics, 4th Edition
    # Chapter 8, Advanced Data Representation, in the interpolation functions section
    rₜ = T(r)
    sₜ = T(s)
    return (1 - rₜ)*(1 - sₜ)*quad.points[1] +
                 rₜ*(1 - sₜ)*quad.points[2] +
                       rₜ*sₜ*quad.points[3] +
                 (1 - rₜ)*sₜ*quad.points[4]
end

function triangulate(quad::Quadrilateral_3D{T}) where {T <: AbstractFloat}
    # Return the two triangles that partition the domain
    A, B, C, D = quad.points
    return (Triangle_3D(A, B, C), Triangle_3D(C, D, A))
end

function area(quad::Quadrilateral_3D{T}) where {T <: AbstractFloat}
    # Using the convex quadrilateral assumption, just return the sum of the areas of the two
    # triangles that partition the quadrilateral. If the convex assumption ever changes, you
    # need to verify that the triangle pair partitions the quadrilateral. Choosing the wrong
    # pair overestimates the area, so just get the areas of both pairs of valid triangles and use
    # the smaller area.
    return sum(area.(triangulate(quad)))
end

#function intersect(l::LineSegment_3D{T}, quad::Quadrilateral_3D{T}) where {T <: AbstractFloat}
#    # Triangulate the quadrilateral, intersect the triangles.
#    tri = triangulate(quad)
#    intersection1 = l ∩ tri[1]
#    if intersection1[1]
#        return true, intersection1[2]
#    end
#    intersection2 = l ∩ tri[2]
#    if intersection2[1]
#        return true, intersection2[2]
#    end
#    return false, Point_3D(T, 0)
#end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
include("lattices.jl")
include("FFGmain.jl")

using LinearAlgebra
using GraphPlot
using LightGraphs
using Plots
using JLD2, FileIO
using Statistics

# define a model
dims = (12,12) # dimension of the lattice
lattice = Lattices.get_SquareLattice(dims,pbc=true) # make a square lattice
filling = 1 # setting filling ≢ 1 means there are holes.

# gplot(SimpleGraph(adjacency_matrix(lattice.graph)),nodelabel=collect(1:nv(lattice.graph)))

# F = Lattices.get_2DCSL_Wavefunctions(lattice,0.25)
# histogram(F.values,bins=50)
# plot(F.values)

Δ = 0.5
model = Dict(
    "dims" => dims,
    "lattice" => lattice,
    "filling" => filling,
    "hamiltonian" => Lattices.get_2DCSL_Wavefunctions(lattice,Δ),
    "fermi_energy" => 0,
    )

# function observable(chain::FFG.GutzwillerChain)
#     state = FFG.get_State(chain)
#     Sz1 = state.spin_config.sc[1]
#     data = Dict("Sz1" => Sz1)
#     return data
# end

n_runs = 300

policy = FFG.SwapNeighborsPolicy()
mc_spec = Dict("mc_steps" => 3E3*prod(dims),
                "mc_warmup_steps" => 300*prod(dims),
                "sample_interval" => 10*prod(dims))

x = collect(1:2:dims[1])
data = zeros(ComplexF64,n_runs,length(x))
for j = 1:n_runs
    println("On run ",j)
    doublegutz = FFG.DoubleGutzwillerChain()
    FFG.init_Chain!(doublegutz,
            model=model,observable=FFG.MeasureSwap,policy=policy,mc_spec=mc_spec
            )
    FFG.runMC!(doublegutz)
    z = FFG.get_Data(doublegutz)["swap"]
    data[j,:] = z
end

real(1+2im)
print(mean(data,dims=1))

d = -log.(abs.(real.(data)))
y = -log.(abs.(mean(data,dims=1)))'

err = std(d,dims=1)'
plotly()
plot(x,y)

z = 2
# sites = unique(neighborhood(lattice.graph,39,1))
# print(FFG.compute_swapregion(doublegutz,sites),"\n")
# #

#
#
# FFG.runMC!(doublegutz)
# # FFG.aggregate_measurements!(doublegutz)
#
# z = FFG.get_Data(doublegutz)["swap"]
# l = collect(1:2:dims[1])
#
# plot(l,-log.(z))
# # print(-log.(FFG.MeasureSwap(doublegutz)["swap"]))
# print(-log.(z))
# print(FFG.get_Diagnostics(doublegutz))
#
# z = [1im 0.5im; 2 2im]
#
# z = zeros((2,2))
#
# typeof(z)
#
# z[1,1] = 1im
#
# z = 1.04 + 0.23im
#
# t = typeof(z)
#
# z = zeros(ComplexF64,(3,3))
#
# typeof(z) == Complex{Float64}
#
# z = zeros(t,(3,3))
#


text:
using DataStreams, Missings, WeakRefStrings, Dates, Test

import Base: ==
==(a::DataStreams.Data.Schema, b::DataStreams.Data.Schema) = DataStreams.Data.types(a) == DataStreams.Data.types(b) && DataStreams.Data.header(a) == DataStreams.Data.header(b) && isequal(size(a), size(b))

mutable struct Source{T}
    sch::DataStreams.Data.Schema
    nt::T
end

DataStreams.Data.schema(s::Source) = s.sch
DataStreams.Data.isdone(s::Source, row, col, rows, cols) = col > cols || (ismissing(rows) ? row > length(s.nt[col]) : row > rows)
DataStreams.Data.streamfrom(s::Source, ::Type{DataStreams.Data.Column}, ::Type{T}, row, col) where {T} = s.nt[col]
DataStreams.Data.streamfrom(s::Source, ::Type{DataStreams.Data.Field}, ::Type{T}, row, col) where {T} = s.nt[col][row]

mutable struct Sink{T}
    nt::T
end

I = (id = Int64[1, 2, 3, 4, 5],
firstname = (Union{String, Missing})["Benjamin", "Wayne", "Sean", "Charles", missing],
lastname = String["Chavez", "Burke", "Richards", "Long", "Rose"],
salary = (Union{Float64, Missing})[missing, 46134.1, 45046.2, 30555.6, 88894.1],
rate = Float64[39.44, 33.8, 15.64, 17.67, 34.6],
hired = (Union{Date, Missing})[Date("2011-07-07"), Date("2016-02-19"), missing, Date("2002-01-05"), Date("2008-05-15")],
fired = DateTime[DateTime("2016-04-07T14:07:00"), DateTime("2015-03-19T15:01:00"), DateTime("2006-11-18T05:07:00"), DateTime("2002-07-18T06:24:00"), DateTime("2007-09-29T12:09:00")]
)
J = NamedTuple{(:_0, (Symbol("_$i") for i = 1:50)...,)}((["0"], ([i] for i =1:50)...,));
K = NamedTuple{((Symbol("_$i") for i = 1:50)...,)}((([i] for i = 1:50)...,));

nms(::NamedTuple{names}) where {names} = names

I_L = Source(DataStreams.Data.Schema(collect(map(eltype, I)), nms(I), 5), I);
I_M = Source(DataStreams.Data.Schema(collect(map(eltype, I)), nms(I), missing), I);

J_L = Source(DataStreams.Data.Schema(collect(map(eltype, J)), nms(J), 1), J);
J_M = Source(DataStreams.Data.Schema(collect(map(eltype, J)), nms(J), missing), J);

K_L = Source(DataStreams.Data.Schema(collect(map(eltype, K)), nms(K), 1), K);
K_M = Source(DataStreams.Data.Schema(collect(map(eltype, K)), nms(K), missing), K);

Sink(sch::DataStreams.Data.Schema, S, append, args...; reference::Vector{UInt8}=UInt8[]) = Sink(Data.NamedTuple(sch, S, append, args...; reference=reference))
Sink(sink, sch::DataStreams.Data.Schema, S, append; reference::Vector{UInt8}=UInt8[]) = Sink(Data.NamedTuple(sink.nt, sch, S, append; reference=reference))
DataStreams.Data.streamtypes(::Type{Sink}) = [DataStreams.Data.Column, DataStreams.Data.Field]
DataStreams.Data.weakrefstrings(::Type{Sink}) = true
DataStreams.Data.streamto!(sink::Sink, ::Type{DataStreams.Data.Field}, val, row, col) =
    (C = getfield(sink.nt, col); row > length(C) ? push!(C, val) : setindex!(C, val, row); return)
DataStreams.Data.streamto!(sink::Sink, ::Type{DataStreams.Data.Column}, column, col) =
    append!(getfield(sink.nt, col), column)

@testset "DataStreams" begin

@testset "Data.Schema" begin

sch = DataStreams.Data.Schema()
@test size(sch) == (0, 0)
@test DataStreams.Data.header(sch) == String[]
@test DataStreams.Data.types(sch) == ()

sch = DataStreams.Data.Schema((Int,))
@test size(sch) == (0, 1)
@test DataStreams.Data.header(sch) == String["Column1"]
@test DataStreams.Data.types(sch) == (Int,)

sch = DataStreams.Data.Schema([Int])
@test size(sch) == (0, 1)
@test DataStreams.Data.header(sch) == String["Column1"]
@test DataStreams.Data.types(sch) == (Int,)

sch = DataStreams.Data.Schema((Int,), ["col1"])
@test size(sch) == (0, 1)
@test DataStreams.Data.header(sch) == String["col1"]
@test DataStreams.Data.types(sch) == (Int,)

@test_throws ArgumentError DataStreams.Data.Schema((String,), ["col1", "col2"])

getR(::DataStreams.Data.Schema{R}) where {R} = R
getT(::DataStreams.Data.Schema{R, T}) where {R, T} = T

sch = DataStreams.Data.Schema((Int,), ["col1"], 1)
@test size(sch) == (1, 1)
@test DataStreams.Data.header(sch) == String["col1"]
@test DataStreams.Data.types(sch) == (Int,)
@test getR(sch) == true
@test getT(sch) == Tuple{Int}

sch = DataStreams.Data.Schema((String,), ["col1"], missing)
@test size(sch) === (missing, 1)
@test DataStreams.Data.header(sch) == String["col1"]
@test DataStreams.Data.types(sch) == (String,)
@test getR(sch) == false
@test getT(sch) == Tuple{String}

@test_throws ArgumentError DataStreams.Data.Schema((String,), ["col1"], -1)

sch = DataStreams.Data.Schema((Int,), ["col1"], 0, Dict("hey"=>"ho"))
@test size(sch) == (0, 1)
@test DataStreams.Data.header(sch) == String["col1"]
@test DataStreams.Data.types(sch) == (Int,)
@test DataStreams.Data.metadata(sch) == Dict("hey"=>"ho")

sch = DataStreams.Data.Schema((Int,), ["col1"])
@test sch["col1"] == 1

end # @testset "Data.Schema"

@testset "Data.stream!" begin

### ALL PAIRS TESTING

# constructed Sink vs. constructed on-the-fly: A, B
# append vs. not; C, D
# w/ transform functions vs. not: E, F
# Data.Field vs. Data.Column: G, H
# <500, >500, homogenous types code path: I, J, K
# knownrows vs. not: L, M

# incompatible source => constructed Sink
@test_throws ArgumentError DataStreams.Data.stream!(I_L, Sink(I))

# incompatible source => otf Sink
@test_throws ArgumentError DataStreams.Data.stream!(I_L, Sink, I)

## Data.Field
DataStreams.Data.streamtype(::Type{Source}, ::Type{DataStreams.Data.Field}) = true
# A, C, E, G, I, L: append to constructed Sink w/ transforms via Data.Field w/ Source=I_J
source = I_L
sink = Sink(deepcopy(I))
transforms = Dict(1=>x->x+1)
DataStreams.Data.stream!(source, sink; append=true, transforms=transforms)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (10, 7)
@test DataStreams.Data.header(sch) == ["id","firstname","lastname","salary","rate","hired","fired"]
@test DataStreams.Data.types(sch) == (Int64, Union{String, Missing}, String, Union{Float64, Missing}, Float64, Union{Date, Missing}, DateTime)
@test sink.nt.id == [1,2,3,4,5,2,3,4,5,6]

# A, C, F, G, J, L: append to constructed Sink w/o transforms via Data.Field w/ Source=J_L
source = J_L;
sink = Sink(deepcopy(J));
Data.stream!(source, sink; append=true);

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (2, 51)
@test DataStreams.Data.header(sch) == collect("_$i" for i = 0:50)
@test DataStreams.Data.types(sch) == (String, (Int for i = 1:50)...,)
@test sink.nt._0 == ["0", "0"]
@test sink.nt._1 == [1, 1]

# A, D, F, G, J, M: replace constructed Sink w/o transforms via Data.Field w/ Source=J_M
source = J_M;
sink = Sink(deepcopy(J));
Data.stream!(source, sink)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (1, 51)
@test DataStreams.Data.header(sch) == collect("_$i" for i = 0:50)
@test DataStreams.Data.types(sch) == (String, (Int for i = 1:50)...,)
@test sink.nt._0 == ["0"]
@test sink.nt._1 == [1]

# B, D, E, G, J, L: replace otf Sink w/ transforms via Data.Field w/ Source=J_L
source = J_L
sink = Sink(deepcopy(J))
transforms = Dict(2=>x->x+1)
Data.stream!(source, Sink, sink.nt; transforms=transforms)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (1, 51)
@test DataStreams.Data.header(sch) == collect("_$i" for i = 0:50)
@test DataStreams.Data.types(sch) == (String, (Int for i = 1:50)...,)
@test sink.nt._0 == ["0"]
@test sink.nt._1 == [2]

# B, D, E, G, K, L: replace otf Sink w/ transforms via Data.Field w/ Source=K_L
source = K_L
sink = Sink(deepcopy(K))
transforms = Dict(2=>x->x+1)
Data.stream!(source, Sink, sink.nt; transforms=transforms)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (1, 50)
@test DataStreams.Data.header(sch) == collect("_$i" for i = 1:50)
@test DataStreams.Data.types(sch) == ((Int for i = 1:50)...,)
@test sink.nt._1 == [1]
@test sink.nt._2 == [3]

# Test unparameterized NamedTuple sink directly
# B, D, F, G, J, M
source = J_M
sink = Data.stream!(source, Data.NamedTuple)

sch = DataStreams.Data.schema(sink)
@test size(sch) == (1, 51)
@test DataStreams.Data.header(sch) == collect("_$i" for i = 0:50)
@test DataStreams.Data.types(sch) == (String, (Int for i = 1:50)...,)
@test sink._0 == ["0"]
@test sink._1 == [1]

## Data.Column
# A, D, E, H, K, M: replace constructed Sink w/ transforms via Data.Column w/ Source=K_M
DataStreams.Data.streamtype(::Type{<:Source}, ::Type{DataStreams.Data.Column}) = true
source = K_M
sink = Sink(deepcopy(K))
transforms = Dict(2=>x->x.+1)
Data.stream!(source, sink; transforms=transforms)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (1, 50)
@test DataStreams.Data.header(sch) == collect("_$i" for i = 1:50)
@test DataStreams.Data.types(sch) == ((Int for i = 1:50)...,)
@test sink.nt._1 == [1]
@test sink.nt._2 == [3]

# B, C, F, H, K, M: append to otf Sink w/o transforms via Data.Column w/ Source=K_M
source = K_M
sink = Sink(deepcopy(K))
Data.stream!(source, Sink, sink.nt; append=true)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (2, 50)
@test DataStreams.Data.header(sch) == collect("_$i" for i = 1:50)
@test DataStreams.Data.types(sch) == ((Int for i = 1:50)...,)
@test sink.nt._1 == [1, 1]
@test sink.nt._2 == [2, 2]

# B, D, F, H, I, M: replace otf Sink w/o transforms via Data.Column w/ Source=I_M
source = I_M
sink = Sink(deepcopy(I))
transforms = Dict()
Data.stream!(source, Sink, sink.nt)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (5, 7)
@test DataStreams.Data.header(sch) == ["id","firstname","lastname","salary","rate","hired","fired"]
@test DataStreams.Data.types(sch) == (Int64, Union{String, Missing}, String, Union{Float64, Missing}, Float64, Union{Date, Missing}, DateTime)
@test sink.nt.id == [1,2,3,4,5]

# B, D, F, H, I, L: replace otf Sink w/o transforms via Data.Column w/ Source=I_L
source = I_L
sink = Sink(deepcopy(I))
Data.stream!(source, Sink, sink.nt)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (5, 7)
@test DataStreams.Data.header(sch) == ["id","firstname","lastname","salary","rate","hired","fired"]
@test DataStreams.Data.types(sch) == (Int64, Union{String, Missing}, String, Union{Float64, Missing}, Float64, Union{Date, Missing}, DateTime)
@test sink.nt.id == [1,2,3,4,5]

# B, D, E, H, J, L: replace otf Sink w/ transforms via Data.Field w/ Source=J_L
source = J_L
sink = Sink(deepcopy(J))
transforms = Dict(2=>x->x.+1)
Data.stream!(source, Sink, sink.nt; transforms=transforms)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (1, 51)
@test DataStreams.Data.header(sch) == collect("_$i" for i = 0:50)
@test DataStreams.Data.types(sch) == (String, (Int for i = 1:50)...,)
@test sink.nt._0 == ["0"]
@test sink.nt._1 == [2]

# Test unparameterized NamedTuple sink directly
# B, D, F, H, I, M
source = I_M
sink = Data.stream!(source, Data.NamedTuple)

sch = DataStreams.Data.schema(sink)
@test size(sch) == (5, 7)
@test DataStreams.Data.header(sch) == ["id","firstname","lastname","salary","rate","hired","fired"]
@test DataStreams.Data.types(sch) == (Int64, Union{String, Missing}, String, Union{Float64, Missing}, Float64, Union{Date, Missing}, DateTime)
@test sink.id == [1,2,3,4,5]

# Test transforms do not raise a BoundsError
source = J_L
sink = Sink(deepcopy(J))
transforms = Dict(i => x -> x for i in 2:50)
sink = Data.stream!(source, Sink, sink.nt; transforms=transforms)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (1, 51)

# Issue #84
source = (a = String["b"],)
sink = Sink(deepcopy(source))
streamed = Data.stream!(source, sink, transforms=Dict("a" => x->Symbol(x)))
@test streamed.nt == (a = Symbol[:b],)

end # @testset "Data.stream!"

@testset "DataStreams NamedTuple" begin

source = (a=[1,2,3], b=[4.0, 5.0, 6.0], c=["hey", "ho", "neighbor"])

sink = Data.stream!(source, Data.RowTable)
@test Data.header(Data.schema(sink)) == ["a", "b", "c"]
@test size(Data.schema(sink)) == (3, 3)

sink = Data.stream!(source, sink)
@test Data.header(Data.schema(sink)) == ["a", "b", "c"]
@test size(Data.schema(sink)) == (3, 3)

sink = Data.stream!(source, sink, append=true)
@test Data.header(Data.schema(sink)) == ["a", "b", "c"]
@test size(Data.schema(sink)) == (6, 3)

sink2 = Data.stream!(sink, Data.RowTable)
@test Data.header(Data.schema(sink)) == ["a", "b", "c"]
@test size(Data.schema(sink)) == (6, 3)

rows = Data.rows(source)
for (i, row) in enumerate(rows)
    @test row[1] == source[1][i]
end

end

@testset "DataStreams with WeakRefStrings" begin
# Test streaming from a source that has a WeakRefStrings column
# to a sink that converts WeakRefStrings to Strings (like Data.RowTable)
wk = WeakRefStringArray(Any[], Union{Missing, WeakRefString{UInt8}}[])
push!(wk, "hey")
push!(wk, "ho")
push!(wk, "neighbor")

source = (a=[1,2,3], b=[4.0, 5.0, 6.0], c=wk)

sink = Data.stream!(source, Data.RowTable)
@test Data.header(Data.schema(sink)) == ["a", "b", "c"]
@test size(Data.schema(sink)) == (3, 3)

sink = Data.stream!(source, sink)
@test Data.header(Data.schema(sink)) == ["a", "b", "c"]
@test size(Data.schema(sink)) == (3, 3)

sink = Data.stream!(source, sink, append=true)
@test Data.header(Data.schema(sink)) == ["a", "b", "c"]
@test size(Data.schema(sink)) == (6, 3)

sink2 = Data.stream!(sink, Data.RowTable)
@test Data.header(Data.schema(sink)) == ["a", "b", "c"]
@test size(Data.schema(sink)) == (6, 3)

rows = Data.rows(source)
for (i, row) in enumerate(rows)
    @test row[1] == source[1][i]
    @test row[2] == source[2][i]
    @test row[3] == source[3][i]
end

end

include("query.jl")

end # @testset "DataStreams"


text:
module FASTQ

import Automa
import Automa.RegExp: @re_str
import BioCore: BioCore, isfilled
import BioSymbols
import BioSequences
import BufferedStreams
import BufferedStreams: BufferedInputStream

include("quality.jl")
include("record.jl")
include("reader.jl")
include("writer.jl")

end


text:
@testset "I473" begin
    s = """
    r = [3,2,1]
    r[r]
    """
    msgs = lintstr(s)
    @test messageset(msgs) == Set([:I473])
    @test contains(msgs[1].message, "value at position #1 is the referenced r")
    @test contains(msgs[1].message, "OK if it represents permutations")
end

@testset "E434" begin
    msgs = lintstr("""
    r = [3,2,1]
    r[1;r]
    """)
    @test messageset(msgs) == Set([:E434])
    @test contains(msgs[1].message, "value at position #2 is the referenced r. Possible typo?")
end


text:
module _Machin

using Brainstorm._Math._Series._Arctan: arctansum

euler(digits::Integer)           = arctansum(digits, [2,   3], 4*[1,  1])
hermann(digits::Integer)         = arctansum(digits, [2,   7], 4*[2, -1])
hutton(digits::Integer)          = arctansum(digits, [3,   7], 4*[2,  1])
machin(digits::Integer)          = arctansum(digits, [5, 239], 4*[4, -1])

gauss(digits::Integer)           = arctansum(digits, [18,  57,  239], 4*[12,  8, -5])
strassnitzky(digits::Integer)    = arctansum(digits, [ 2,   5,    8], 4*[ 1,  1,  1])
klingenstierna(digits::Integer)  = arctansum(digits, [10, 239,  515], 4*[ 8, -1, -4])
euler3(digits::Integer)          = arctansum(digits, [ 5,  70,   99], 4*[ 4, -1,  1])
loney(digits::Integer)           = arctansum(digits, [ 4,  20, 1985], 4*[ 3,  1,  1])
stoermer(digits::Integer)        = arctansum(digits, [ 8,  57,  239], 4*[ 6,  2,  1])

takano1982(digits::Integer)      = arctansum(digits, Int64[49, 57, 239, 110443], 4*[12, 32,  -5, 12])
stoermer1896(digits::Integer)    = arctansum(digits, Int64[57, 239, 682, 12943], 4*[44,  7, -12, 24])

hwang1997(digits::Integer)       = arctansum(digits, Int64[239, 1023, 5832, 110443, 4841182, 6826318], 4*[183, 32, -68, 12, -12, -100])
hwang2003(digits::Integer)       = arctansum(digits, BigInt[239, 1023, 5832, 113021, 6826318, 33366019650, 43599522992503626068], 4*[183, 32, -68, 12, -100, -12, 12])
wetherfield2004(digits::Integer) = arctansum(digits, BigInt[107, 1710, 103697, 2513489, 18280007883, 7939642926390344818, 3054211727257704725384731479018], 4*[83, 17, -22, -24, -44, 12, 22])

end


text:
"""
Although the BSE exposes a special SOAP-based API for programmatic access to its database, I had too much trouble getting the desired results from it. Therefore the solution in this module is based around regexing of the
human readable HTML output.
"""
module BasisSetExchangeModule
using ..BasisSetModule

export downloadBasisSetExchange, computeBasisSetsAvailable
import ..BasisSetModule.BasisSet
using HTTP
using Printf
using ZipFile

bseZip    = homedir()*"/pqs.zip"
bseFormat = "pqs"

function Base.getindex(dir::ZipFile.Reader,fn::String)
  for f in dir.files
    if f.name == fn
      return f
    end
  end
  @error "$fn not found."
end

function downloadBasisSetExchange(;bseFormat=bseFormat, bseZip=bseZip)
  response = HTTP.request("GET","http://www.basissetexchange.org/download/current/$bseFormat/zip")
  write(bseZip,response.body)
  @assert isfile(bseZip)
  println("Downloaded $bseFormat formatted basis sets to $bseZip.")
end

function computeBasisSetsAvailable(;bseZip=bseZip,bseFormat=bseFormat)
  list = Dict{String,String}()
  zip = ZipFile.Reader(bseZip)
  for f in zip.files
    if occursin(Regex("\\.$bseFormat\$"),f.name)
      shortname = match(Regex(".*/([^\\.]*).*$bseFormat"), f.name).captures[1]
      push!(list, shortname => f.name)
    end
  end
  return list
end

function computeBasisSetExchangeEntry(name::AbstractString; bseZip=bseZip)
  dict = computeBasisSetsAvailable(bseZip=bseZip)
  if haskey(dict,name)
    return dict[name]
  else
    result = String[]
    for (shortname,entry) in dict
      if (occursin(Regex(name,"i"),shortname))
        push!(result,entry)
      end
    end
  end
  if length(result) == 0
    error("No matching basis set entry found (searching for $name in $bseZip)")
  elseif length(result)>1
    display(result)
    error("More than one matching entry found.")
  end
  return result[1]
end

"""
    BasisSet(name)

creates a BasisSet object from the specified file (possibly with a corresponding file extension).
If no such file can be found the argument is interpreted as the name of a basis set and that basis
set obtained from Basis Set Exchange.
"""
function BasisSet(name::AbstractString; bseZip=bseZip)

  if isfile(name)
    return readBasisSetTX93(name)
  end

  for ext in ("tx93", "TX93", "Tx93")
    if isfile("$name.$ext")
      return readBasisSetTX93("$name.$ext")
    end
  end

  @info("BasisSet file $name not found on disk. Obtaining from BasisSetExchange...")
  if !isfile(bseZip)
    downloadBasisSetExchange(bseZip=bseZip)
  end
  location   = computeBasisSetExchangeEntry(name)
  str = ""
  open(bseZip) do bse # trying to get around ZipFile Issue #14
    str      = read(ZipFile.Reader(bse)[location], String)
  end
  return BasisSetModule.BasisSetTX93(str)
end

end # module


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing

    VecBitSites = Sites{Bit, Int, 1}
    for L = 4:10, BC in [Lattices.Periodic, Lattices.Fixed], DType in [Float32, Float64, ComplexF32, ComplexF64]
        precompile(Base.getindex, (SumSingleSiteType{DType, Matrix{DType}, Lattices.Chain{L, BC}}, VecBitSites, VecBitSites))
    end
end

max_stars_repo_name: Roger-luo/QuHamiltonian

max_issues_repo_name: Roger-luo/QuHamiltonian

max_forks_repo_name: Roger-luo/QuHamiltonian


text:
"""
    inverse_axes(a::AbstractMatrix) = inverse_axes(axes(a))
    inverse_axes(a::Tuple{I1,I2}) -> Tuple{I2,I1}

Returns the inverted axes of `a`, corresponding to the `inv` method from the
`LinearAlgebra` package in the standard library.

## Examples
```jldoctest
julia> inverse_axes((Index{:a}(1:4), Index{:b}(1:4)))
(Index{b}(1:4 => Base.OneTo(4)), Index{a}(1:4 => Base.OneTo(4)))
```
"""
inverse_axes(x::AbstractMatrix) = inverse_axes(axes(x))
inverse_axes(x::Tuple{I1,I2}) where {I1,I2} = (last(x), first(x))


text:
using DataFrames, TidyStanza
using Test

include("relocate-test.jl")

df = DataFrame(x = repeat(1:3, inner = 2, outer = 2),
       a = repeat(4:6, inner = 2, outer = 2),
       b = repeat(7:9, inner = 2, outer = 2),
       val1 = ["ce_val1_1","cf_val1_1","ce_val1_2","cf_val1_2","ce_val1_3","cf_val1_3","de_val1_1","df_val1_1","de_val1_2","df_val1_2","de_val1_3","df_val1_3"],
       val2 = ["ce_val2_1","cf_val2_1","ce_val2_2","cf_val2_2","ce_val2_3","cf_val2_3","de_val2_1","df_val2_1","de_val2_2","df_val2_2","de_val2_3","df_val2_3"],
       cname1 = repeat(["c", "d"], inner = 6),
       cname2 = repeat(["e", "f"], 6)
)

wider_df1 = DataFrame(
    x = repeat(1:3, inner=2),
    a =  repeat(4:6, inner=2),
    b =  repeat(7:9, inner=2),
    cname2 = repeat(["e", "f"], 3),
    val1_c = ["ce_val1_1","cf_val1_1", "ce_val1_2", "cf_val1_2", "ce_val1_3", "cf_val1_3"],
    val1_d = ["de_val1_1","df_val1_1", "de_val1_2", "df_val1_2", "de_val1_3", "df_val1_3"],
    val2_c = ["ce_val2_1","cf_val2_1", "ce_val2_2", "cf_val2_2", "ce_val2_3", "cf_val2_3"],
    val2_d = ["de_val2_1","df_val2_1", "de_val2_2", "df_val2_2", "de_val2_3", "df_val2_3"]
)

wider_df2 = DataFrame(
    x = repeat(1:3, inner = 4),
    a = repeat(4:6, inner = 4),
    b = repeat(7:9, inner = 4),
    val2 = ["ce_val2_1","cf_val2_1","de_val2_1","df_val2_1","ce_val2_2","cf_val2_2","de_val2_2","df_val2_2","ce_val2_3","cf_val2_3","de_val2_3","df_val2_3"],
    val1_c_e = vcat("ce_val1_1", repeat([missing],3), "ce_val1_2",repeat([missing],3),  "ce_val1_3", repeat([missing],3) ),
    val1_c_f = vcat(repeat([missing],1), "cf_val1_1", repeat([missing],3), "cf_val1_2",repeat([missing],3),  "cf_val1_3", repeat([missing],2) ),
    val1_d_e = vcat(repeat([missing],2), "de_val1_1", repeat([missing],3), "de_val1_2",repeat([missing],3),  "de_val1_3", repeat([missing],1) ),
    val1_d_f = vcat(repeat([missing],3), "df_val1_1", repeat([missing],3), "df_val1_2",repeat([missing],3),  "df_val1_3")
)

wider_df3 = DataFrame(
    x = 1:3,
    a = 4:6,
    b = 7:9,
    val1_c_e = ["ce_val1_1", "ce_val1_2", "ce_val1_3"],
    val1_c_f = ["cf_val1_1", "cf_val1_2", "cf_val1_3"],
    val1_d_e = ["de_val1_1", "de_val1_2", "de_val1_3"],
    val1_d_f = ["df_val1_1", "df_val1_2", "df_val1_3"],
    val2_c_e = ["ce_val2_1", "ce_val2_2", "ce_val2_3"],
    val2_c_f = ["cf_val2_1", "cf_val2_2", "cf_val2_3"],
    val2_d_e = ["de_val2_1", "de_val2_2", "de_val2_3"],
    val2_d_f = ["df_val2_1", "df_val2_2", "df_val2_3"]
)

@testset "Tidyverse.jl" begin
    # Write your tests here.

    @test pivot_wider(df, names_from = :cname1, values_from = [:val1,:val2]) == wider_df1
    tmp1 = pivot_wider(df, names_from = [:cname1, :cname2], values_from = :val1)
    sort!(tmp1, [:x, :a, :b, :val2])
    @test isequal(tmp1, wider_df2)
    @test pivot_wider(df, names_from = [:cname1,:cname2], values_from = [:val1,:val2]) == wider_df3
    @test pivot_wider(df, names_from = "cname1", values_from = ["val1","val2"]) == wider_df1
    tmp2 = pivot_wider(df, names_from = ["cname1", "cname2"], values_from = "val1")
    sort!(tmp2, [:x, :a, :b, :val2])
    @test isequal(tmp2, wider_df2)
    @test pivot_wider(df, names_from = ["cname1", "cname2"], values_from = ["val1","val2"]) == wider_df3
end


text:
# util/tokenize.jl
# ================
#
# A module defining types and methods for tokenizing strings.
#
# Useful in creating parsers.
#
# This file is a part of BioJulia. License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module Tokenize

using Compat
import Compat.String

export
    Tokenizer, TokenizerResult, tokenize

"""
Tokenizer type that is responsible for splitting a string into token according to a
regex specification. The regex specification is stored in `dict`.

**Fields:**

* `dict`:      A Dictionary of String keys and Regex values.
* `tokenizer`: A regex that is generated from the specification.
"""
type Tokenizer
    dict::Dict{String, Regex}
    tokenizer::Regex
    """
    Inner constructor for Tokenizers. Accepts a vector of String, Regex
    tuples that makes the token specifications.

    **Parameters:**

    * `x`: An array of `Compat.String`, `Regex` tuples.

    **Returns:** Instance of `Tokenizer`
    """
    function Tokenizer(x::Vector{Tuple{String, Regex}})
        dictvalues = [i[2] for i in x]
        stringvalues = Compat.String[i.pattern for i in dictvalues]
        combinedstring = join(stringvalues, "|")
        finalstring = "($combinedstring)"
        return new([i => j for (i, j) in x], Regex(finalstring))
    end
end

immutable TokenizerResult
    tokens::Vector{String}
    failiures::Vector{String}

    function TokenizerResult(tokens = String[], failiures = String[])
        return new(tokens, failiures)
    end
end

function tokenize(s::String, t::Tokenizer)
    lastpos = 0
    out = TokenizerResult()
    for mat in eachmatch(t.tokenizer, s)
        if mat.offset > lastpos + 1
            push!(failiures, str[lastpos + 1 : mat.offset - 1])
        end
        push!(out.tokens, mat.match)
        lastpos = mat.offset + length(mat.match) - 1
    end
    return out
end

end


text:
module mlmlTest

  using FactCheck
  using Distributions
  include(Pkg.dir("MethylUtils", "src","mlml-bs-ox-typed.jl"))
   include(Pkg.dir("MethylUtils", "src","Features.jl"))
  MLML_PATH = Pkg.dir("MethylUtils", "bin","mlml")

  function test_mlml(M_B, u_B, m_O, U_O, p_m,p_h)
      ratio_bs = M_B / (M_B+u_B)
      bs_total = M_B+u_B
      ratio_ox = m_O / (m_O+U_O)
      ox_total = m_O+U_O

      (em_p_m, em_p_h, iter) = EM(M_B, u_B, m_O, U_O, p_m,p_h)
      p_u = 1 - em_p_m - em_p_h
      println("EM: p_m: $em_p_m, p_h: $em_p_h, p_u: $p_u")

      lineA = ["chr1" "2" "3" "CpG:$bs_total" ratio_bs "+"]
      lineB = ["chr1" "2" "3" "CpG:$ox_total" ratio_ox "+"]
      println(typeof(lineA))
      writedlm( "bs.bed", lineA, "\t")
      writedlm( "ox.bed", lineB, "\t")

      cmd = `$MLML_PATH -tolerance 1e-15 -bsseq bs.bed -oxbsseq ox.bed`
      output = readall(cmd)
      (chr,start,stop,mc_txt,hmc,strand) = split(output, "\t")
      (text,mc) = split(mc_txt,":")
      ( mc,hmc) = map( parsefloat, [ mc,hmc])
      (em_p_m,em_p_h,mc,hmc) = map( x -> round(x,4),  [em_p_m,em_p_h,mc,hmc])
      return (em_p_m,em_p_h,mc,hmc)
  end

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

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

  # M_B = Counts of C's in Bis-Seq - including 5mc (m) and 5hmc (h)
  # u_B = counts of T's in Bis-Seq
  # m_O = Counts of C's in Oxidative Bis-Seq - 5mc
  # U_O = Counts of T's in Oxidative Bis-Seq - including 5hmc (h) and umc (u)

  (M_B,u_B,m_O,U_O) = map( int16, [5,5,7,3 ] )

  ##  M ratio  = 0.5 (5/10),  m ratio  = 0.7 (7/10) - overshoot example

  starting_estimate_p_m = M_B / (M_B + u_B)
  starting_estimate_p_h = 1 - starting_estimate_p_m
  betas = Array(Float64,M_B+1)
  betas = expectation(  M_B, u_B, starting_estimate_p_m, starting_estimate_p_h, betas )

  # beta has length M_B+1, as it is calculated over each possible combination of counts
  # for the latent variables (m_B, h_B) - i.e.  ( (0,5), (1,4), (2,3), ..., (5,0) )

  expected_betas = [0.03125,0.15624999999999994,0.3125,0.3125,0.15624999999999994,0.03125]

  # maximisation calculates p_h
  new_p_h = maximisation(M_B, u_B, starting_estimate_p_m, betas)
  new_p_m = update_pm(M_B, u_B, m_O, U_O, betas)
  (expected_new_p_h, expected_new_p_m) = (0.16666666666666666, 0.475)

  ( em_p_m, em_p_h, mlml_p_m, mlml_p_h) = test_mlml(M_B, u_B, m_O, U_O, starting_estimate_p_m,starting_estimate_p_h)

  println( "$em_p_m, $em_p_h, $mlml_p_m, $mlml_p_h")

  facts("expectation") do
   @fact length(betas) => M_B+1
   @fact betas => expected_betas
  end

  facts("maximisation") do
    @fact new_p_h => expected_new_p_h
    @fact new_p_m => expected_new_p_m
  end

  facts("EM1") do
    @fact em_p_m => mlml_p_m
    @fact em_p_h => mlml_p_h
  end

  # 10/15 (0.666)  10/20 (0.5) - No overshoot
  # M_B - m_O = p_h = 0.166
  (M_B,u_B,m_O,U_O) = map( int16, [10,5,10,10 ] )
  starting_estimate_p_m = M_B / (M_B + u_B)
  starting_estimate_p_h = 1 - starting_estimate_p_m

  #
  #
  expected_p_m = ( m_O / ( m_O + U_O ) )
  expected_p_h = starting_estimate_p_m - expected_p_m
  expected_p_h = round(expected_p_h,4)
  expected_p_m = round(expected_p_m,4)

  ( em_p_m, em_p_h, mlml_p_m, mlml_p_h) = test_mlml(M_B, u_B, m_O, U_O, starting_estimate_p_m,starting_estimate_p_h)
   println( "$em_p_m, $em_p_h, $mlml_p_m, $mlml_p_h")

 facts("EM2") do
    @fact em_p_m => expected_p_m
    @fact em_p_h => expected_p_h
  end

  # M_B = 90/190 (0.4736842)
  # m_O = 85/185 (0.4594595)
  # p_h = 0.01422475
  # No overshoot - but only by a small amount
  (M_B,u_B,m_O,U_O) = map( int16, [90,100,85,100 ] )
  starting_estimate_p_m = M_B / (M_B + u_B)
  starting_estimate_p_h = 1 - starting_estimate_p_m

  expected_p_m = ( m_O / ( m_O + U_O ) )
  expected_p_h = starting_estimate_p_m - expected_p_m
  expected_p_h = round(expected_p_h,4)
  expected_p_m = round(expected_p_m,4)

  ( em_p_m, em_p_h, mlml_p_m, mlml_p_h) = test_mlml(M_B, u_B, m_O, U_O, starting_estimate_p_m,starting_estimate_p_h)
   println( "$em_p_m, $em_p_h, $mlml_p_m, $mlml_p_h")

 facts("EM3") do
    @fact em_p_m => expected_p_m
    @fact em_p_h => expected_p_h
  end

  #### OTHER WAY AROUND ####
  # M_B = 85/185
  # m_O = 90/190
  # Overshoot - but only by a small amount
  (M_B,u_B,m_O,U_O) = map( int16, [85,100,90,100 ] )
  starting_estimate_p_m = M_B / (M_B + u_B)
  starting_estimate_p_h = 1 - starting_estimate_p_m

  expected_p_m = ( m_O / ( m_O + U_O ) )
  expected_p_h = starting_estimate_p_m - expected_p_m
  expected_p_h = round(expected_p_h,4)
  expected_p_m = round(expected_p_m,4)

  ( em_p_m, em_p_h, mlml_p_m, mlml_p_h) = test_mlml(M_B, u_B, m_O, U_O, starting_estimate_p_m,starting_estimate_p_h)
   println( "$em_p_m, $em_p_h, $mlml_p_m, $mlml_p_h")

  facts("EM3") do
    @fact em_p_m => mlml_p_m
    @fact em_p_h => mlml_p_h
  end

  # TEST MLML
  bs_counts = rand_CytosineCounts(10,100,5)
  ox_counts = rand_CytosineCounts(10,100,5)
  cgs = run_mlml( bs_counts , ox_counts )
  length(cgs)

  facts("mlml") do
    @fact length(cgs) => 10
  end

  # exitstatus()
end

# put in the other tests
# only calculate em when need to
# parallelise what I can
# put up to travis


text:
import Base: hash,
    string, show, tryparse,
    Float16, Float32, Float64, BigFloat,
    Bool, size, axes

import Base: precision, exponent_mask, significand_mask, eps, exponent, significand,
    zero, one, floatmax, floatmin, nextfloat, prevfloat, frexp, ldexp, modf, decompose,
    ==, <, <=, >=, >, !=, isless, isequal, cmp,
    iszero, isone, isfinite, isinteger, isinf, isnan, issubnormal,
    signbit, sign, flipsign, copysign,
    +, -, *, /, \, abs, inv, sqrt, cbrt,
    mod, mod1, rem, fld, fld1, div, cld, fldmod, fldmod1, divrem,
    ceil, floor, trunc, round,
    muladd, fma, hypot, clamp,
    sum, prod, cumsum

import Base.Math: log, log1p, log10, log2, exp, expm1, exp10, exp2,
    sin, cos, tan, csc, sec, cot, sinpi, cospi,
    asin, acos, atan, acsc, asec, acot,
    sinh, cosh, tanh, csch, sech, coth,
    asinh, acosh, atanh, acsch, asech, acoth

import LinearAlgebra: norm, normalize, dot,
    *, det, tr, inv, lu, qr, factorize,
    svdvals, eigvals, eigvecs

#=
import SpecialFunctions: export
    airyai,
    airyaiprime,
    airybi,
    airybiprime,
    airyaix,
    airyaiprimex,
    airybix,
    airybiprimex,
    besselh,
    besselhx,
    besseli,
    besselix,
    besselj,
    besselj0,
    besselj1,
    besseljx,
    besselk,
    besselkx,
    bessely,
    bessely0,
    bessely1,
    besselyx,
    dawson,
    ellipk,
    ellipe,
    erf,
    erfc,
    erfcinv,
    erfcx,
    erfi,
    erfinv,
    eta,
    digamma,
    invdigamma,
    polygamma,
    trigamma,
    gamma_inc,
    beta_inc,
    beta_inc_inv,
    gamma_inc_inv,
    hankelh1,
    hankelh1x,
    hankelh2,
    hankelh2x,
    zeta,
    sinint,
    cosint,
    lbinomial
=#


text:
using ElectricFields
using Documenter

isdefined(Main, :NOPLOTS) && NOPLOTS || include("plots.jl")

makedocs(;
    modules=[ElectricFields],
    authors="Stefanos Carlström <stefanos.carlstrom@gmail.com> and contributors",
    repo="https://github.com/jagot/ElectricFields.jl/blob/{commit}{path}#L{line}",
    sitename="ElectricFields.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://www.tipota.org/ElectricFields.jl",
        assets=String[],
        mathengine = MathJax(Dict(
            :TeX => Dict(
                :equationNumbers => Dict(:autoNumber => "AMS"),
                :Macros => Dict(
                    :defd => "≝",
                    :abs => ["\\left|#1\\right|", 1],
                    :vec => ["\\mathbf{#1}", 1],
                    :mat => ["\\mathsf{#1}", 1],
                    :conj => ["#1^*", 1],
                    :im => "\\mathrm{i}",
                    :diff => ["\\mathrm{d}#1\\,", 1],
                    :bmat => ["\\begin{bmatrix}#1\\end{bmatrix}", 1]
                )
            )
        ))
    ),
    pages=[
        "Home" => "index.md",
        "Field types" => "field_types.md",
        "Envelopes" => "envelopes.md",
        "Carriers" => "carriers.md",
        "Field properties" => "properties.md",
        "Reference" => "reference.md",
    ],
    doctest=false
)

deploydocs(;
    repo="github.com/jagot/ElectricFields.jl",
)


text:
"""
Basic continuous univariate parameter state type

# Constructors

## BasicContUnvParameterState{N<:Real}(value::N, <optional arguments>)

Construct a basic continuous univariate parameter state with some ``value``.

### Optional arguments:

* ``diagnostickeys::Vector{Symbol}=Symbol[]``: the diagnostic keys of the state.
* ``diagnosticvalues::Vector=Array(Any, length(diagnostickeys))``: the diagnostic values of the state.

## BasicContUnvParameterState{N<:Real}(<optional arguments>)

Construct a basic continuous univariate parameter state with an uninitialized ``value`` (``NaN``).

### Optional arguments:

* ``diagnostickeys::Vector{Symbol}=Symbol[]``: the diagnostic keys of the state.
* ``::Type{N}=Float64``: the element type of the state value.
* ``diagnosticvalues::Vector=Array(Any, length(diagnostickeys))``: the diagnostic values of the state.

# Examples

```julia
julia> state = BasicContUnvParameterState(-1.25, [:accept], [false])
Lora.BasicContUnvParameterState{Float64}(-1.25,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,Bool[false],[:accept])

julia> state.value
-1.25

julia> diagnostics(state)
Dict{Symbol,Bool} with 1 entry:
  :accept => false
```
"""
type BasicContUnvParameterState{N<:Real} <: ParameterState{Continuous, Univariate}
  "Vector value of basic continuous univariate parameter state"
  value::N
  "Value of log-likelihood at the state's value"
  loglikelihood::N
  "Value of log-prior at the state's value"
  logprior::N
  "Value of log-target at the state's value"
  logtarget::N
  "Value of gradient of log-likelihood at the state's value"
  gradloglikelihood::N
  "Value of gradient of log-prior at the state's value"
  gradlogprior::N
  "Value of gradient of log-target at the state's value"
  gradlogtarget::N
  "Value of metric tensor of log-likelihood at the state's value"
  tensorloglikelihood::N
  "Value of metric tensor of log-prior at the state's value"
  tensorlogprior::N
  "Value of metric tensor of log-target at the state's value"
  tensorlogtarget::N
  "Value of derivatives of metric tensor of log-likelihood at the state's value"
  dtensorloglikelihood::N
  "Value of derivatives of metric tensor of log-prior at the state's value"
  dtensorlogprior::N
  "Value of derivatives of metric tensor of log-target at the state's value"
  dtensorlogtarget::N
  "Diagnostic values associated with the sampling of the state"
  diagnosticvalues::Vector
  "Diagnostic keys associated with the sampling of the state"
  diagnostickeys::Vector{Symbol}
end

function BasicContUnvParameterState{N<:Real}(
  value::N,
  diagnostickeys::Vector{Symbol}=Symbol[],
  diagnosticvalues::Vector=Array(Any, length(diagnostickeys))
)
  v = convert(N, NaN)
  BasicContUnvParameterState{N}(value, v, v, v, v, v, v, v, v, v, v, v, v, diagnosticvalues, diagnostickeys)
end

BasicContUnvParameterState{N<:Real}(
  diagnostickeys::Vector{Symbol}=Symbol[],
  ::Type{N}=Float64,
  diagnosticvalues::Vector=Array(Any, length(diagnostickeys))
) =
  BasicContUnvParameterState(convert(N, NaN), diagnostickeys, diagnosticvalues)

value_support{N<:Real}(::Type{BasicContUnvParameterState{N}}) = Continuous
value_support(::BasicContUnvParameterState) = Continuous

variate_form{N<:Real}(::Type{BasicContUnvParameterState{N}}) = Univariate
variate_form(::BasicContUnvParameterState) = Univariate

Base.eltype{N<:Real}(::Type{BasicContUnvParameterState{N}}) = N
Base.eltype{N<:Real}(::BasicContUnvParameterState{N}) = N

generate_empty(state::BasicContUnvParameterState) = BasicContUnvParameterState(state.diagnostickeys, eltype(state))


text:
using Distributed
using Dates
import REPL
using Printf: @sprintf

# parse some command-line arguments
function extract_flag!(args, flag, default=nothing)
    for f in args
        if startswith(f, flag)
            # Check if it's just `--flag` or if it's `--flag=foo`
            if f != flag
                val = split(f, '=')[2]
                if default !== nothing && !(typeof(default) <: AbstractString)
                  val = parse(typeof(default), val)
                end
            else
                val = default
            end

            # Drop this value from our args
            filter!(x -> x != f, args)
            return (true, val)
        end
    end
    return (false, default)
end
do_help, _ = extract_flag!(ARGS, "--help")
if do_help
    println("""
        Usage: runtests.jl [--help] [--list] [--jobs=N] [TESTS...]

               --help             Show this text.
               --list             List all available tests.
               --thorough         Don't allow skipping tests that are not supported.
               --quickfail        Fail the entire run as soon as a single test errored.
               --jobs=N           Launch `N` processes to perform tests (default: Threads.nthreads()).
               --gpus=N           Expose `N` GPUs to test processes (default: 1).
               --sanitize[=tool]  Run the tests under `compute-sanitizer`.
               --snoop=FILE       Snoop on compiled methods and save to `FILE`.

               Remaining arguments filter the tests that will be executed.""")
    exit(0)
end
_, jobs = extract_flag!(ARGS, "--jobs", Threads.nthreads())
_, gpus = extract_flag!(ARGS, "--gpus", 1)
do_sanitize, sanitize_tool = extract_flag!(ARGS, "--sanitize", "memcheck")
do_snoop, snoop_path = extract_flag!(ARGS, "--snoop")
do_thorough, _ = extract_flag!(ARGS, "--thorough")
do_quickfail, _ = extract_flag!(ARGS, "--quickfail")

include("setup.jl")     # make sure everything is precompiled

# choose tests
const tests = ["initialization"]    # needs to run first
const test_runners = Dict()
## files in the test folder
for (rootpath, dirs, files) in walkdir(@__DIR__)
  # find Julia files
  filter!(files) do file
    endswith(file, ".jl") && file !== "setup.jl" && file !== "runtests.jl"
  end
  isempty(files) && continue

  # strip extension
  files = map(files) do file
    file[1:end-3]
  end

  # prepend subdir
  subdir = relpath(rootpath, @__DIR__)
  if subdir != "."
    files = map(files) do file
      joinpath(subdir, file)
    end
  end

  append!(tests, files)
  for file in files
    test_runners[file] = ()->include("$(@__DIR__)/$file.jl")
  end
end
## GPUArrays testsuite
for name in keys(TestSuite.tests)
    push!(tests, "gpuarrays$(Base.Filesystem.path_separator)$name")
    test_runners["gpuarrays$(Base.Filesystem.path_separator)$name"] = ()->TestSuite.tests[name](CuArray)
end
unique!(tests)

# parse some more command-line arguments
## --list to list all available tests
do_list, _ = extract_flag!(ARGS, "--list")
if do_list
    println("Available tests:")
    for test in sort(tests)
        println(" - $test")
    end
    exit(0)
end
## no options should remain
optlike_args = filter(startswith("-"), ARGS)
if !isempty(optlike_args)
    error("Unknown test options `$(join(optlike_args, " "))` (try `--help` for usage instructions)")
end
## the remaining args filter tests
if !isempty(ARGS)
  filter!(tests) do test
    any(arg->startswith(test, arg), ARGS)
  end
end

# check that CI is using the requested toolkit
toolkit_release = CUDA.toolkit_release() # ensure artifacts are downloaded
if parse(Bool, get(ENV, "CI", "false")) && haskey(ENV, "JULIA_CUDA_VERSION")
  @test toolkit_release == VersionNumber(ENV["JULIA_CUDA_VERSION"])
end

# find suitable devices
@info "System information:\n" * sprint(io->CUDA.versioninfo(io))
candidates, driver_version, cuda_driver_version = if has_nvml()
    [(index=i,
      uuid=NVML.uuid(dev),
      name=NVML.name(dev),
      cap=NVML.compute_capability(dev),
      mem=NVML.memory_info(dev).free)
     for (i,dev) in enumerate(NVML.devices())],
    NVML.driver_version(),
    NVML.cuda_driver_version()
else
    # using CUDA to query this information requires initializing a context,
    # which might fail if the device is heavily loaded.
    [(device!(dev);
     (index=i,
      uuid=uuid(dev),
      name=CUDA.name(dev),
      cap=capability(dev),
      mem=CUDA.available_memory()))
     for (i,dev) in enumerate(devices())],
    "(unknown)",
    CUDA.version()
end
## only consider devices that are fully supported by our CUDA toolkit, or tools can fail.
## NOTE: we don't reuse target_support which is also bounded by LLVM support,
#        and is used to pick a codegen target regardless of the actual device.
cuda_support = CUDA.cuda_compat()
filter!(x->x.cap in cuda_support.cap, candidates)
## only consider recent devices if we want testing to be thorough
if do_thorough
    filter!(x->x.cap >= v"7.0", candidates)
end
isempty(candidates) && error("Could not find any suitable device for this configuration")
## order by available memory, but also by capability if testing needs to be thorough
sort!(candidates, by=x->x.mem)
## apply
picks = reverse(candidates[end-gpus+1:end])   # best GPU first
ENV["CUDA_VISIBLE_DEVICES"] = join(map(pick->"GPU-$(pick.uuid)", picks), ",")
@info "Testing using $(length(picks)) device(s): " * join(map(pick->"$(pick.index). $(pick.name) (UUID $(pick.uuid))", picks), ", ")

# determine tests to skip
skip_tests = []
has_cudnn() || push!(skip_tests, "cudnn")
has_nvml() || push!(skip_tests, "nvml")
if !has_cutensor() || CUDA.version() < v"10.1" || first(picks).cap < v"7.0"
    push!(skip_tests, "cutensor")
end
is_debug = ccall(:jl_is_debugbuild, Cint, ()) != 0
if VERSION < v"1.5-" || first(picks).cap < v"7.0"
    push!(skip_tests, "device/wmma")
end
if Sys.ARCH == :aarch64
    # CUFFT segfaults on ARM
    push!(skip_tests, "cufft")
end
for (i, test) in enumerate(skip_tests)
    # we find tests by scanning the file system, so make sure the path separator matches
    skip_tests[i] = replace(test, '/'=>Base.Filesystem.path_separator)
end
# skip_tests is a list of patterns, expand it to actual tests we were going to run
skip_tests = filter(test->any(skip->occursin(skip,test), skip_tests), tests)
if do_thorough
    # we're not allowed to skip tests, so make sure we will mark them as such
    all_tests = copy(tests)
    if !isempty(skip_tests)
        @error "Skipping the following tests: $(join(skip_tests, ", "))"
        filter!(!in(skip_tests), tests)
    end
else
    if !isempty(skip_tests)
        @info "Skipping the following tests: $(join(skip_tests, ", "))"
        filter!(!in(skip_tests), tests)
    end
    all_tests = copy(tests)
end

# add workers
const test_exeflags = Base.julia_cmd()
filter!(test_exeflags.exec) do c
    return !(startswith(c, "--depwarn") || startswith(c, "--check-bounds"))
end
push!(test_exeflags.exec, "--check-bounds=yes")
push!(test_exeflags.exec, "--startup-file=no")
push!(test_exeflags.exec, "--depwarn=yes")
if Base.JLOptions().project != C_NULL
    push!(test_exeflags.exec, "--project=$(unsafe_string(Base.JLOptions().project))")
end
const test_exename = popfirst!(test_exeflags.exec)
function addworker(X; kwargs...)
    exename = if do_sanitize
        sanitizer = CUDA.compute_sanitizer()
        @info "Running under $(readchomp(`$sanitizer --version`))"
        # NVIDIA bug 3263616: compute-sanitizer crashes when generating host backtraces
        `$sanitizer --tool $sanitize_tool --launch-timeout=0 --show-backtrace=no --target-processes=all --report-api-errors=no $test_exename`
    else
        test_exename
    end

    withenv("JULIA_NUM_THREADS" => 1, "OPENBLAS_NUM_THREADS" => 1) do
        procs = addprocs(X; exename=exename, exeflags=test_exeflags, kwargs...)
        @everywhere procs include($(joinpath(@__DIR__, "setup.jl")))
        procs
    end
end
addworker(min(jobs, length(tests)))

# prepare to snoop on the compiler
if do_snoop
    @info "Writing trace of compiled methods to '$snoop_path'"
    snoop_io = open(snoop_path, "w")
end

# pretty print information about gc and mem usage
testgroupheader = "Test"
workerheader = "(Worker)"
name_align        = maximum([textwidth(testgroupheader) + textwidth(" ") +
                             textwidth(workerheader); map(x -> textwidth(x) +
                             3 + ndigits(nworkers()), tests)])
elapsed_align     = textwidth("Time (s)")
gc_align      = textwidth("GC (s)")
percent_align = textwidth("GC %")
alloc_align   = textwidth("Alloc (MB)")
rss_align     = textwidth("RSS (MB)")
printstyled(" "^(name_align + textwidth(testgroupheader) - 3), " | ")
printstyled("         | ---------------- GPU ---------------- | ---------------- CPU ---------------- |\n", color=:white)
printstyled(testgroupheader, color=:white)
printstyled(lpad(workerheader, name_align - textwidth(testgroupheader) + 1), " | ", color=:white)
printstyled("Time (s) | GC (s) | GC % | Alloc (MB) | RSS (MB) | GC (s) | GC % | Alloc (MB) | RSS (MB) |\n", color=:white)
print_lock = stdout isa Base.LibuvStream ? stdout.lock : ReentrantLock()
if stderr isa Base.LibuvStream
    stderr.lock = print_lock
end
function print_testworker_stats(test, wrkr, resp)
    @nospecialize resp
    lock(print_lock)
    try
        printstyled(test, color=:white)
        printstyled(lpad("($wrkr)", name_align - textwidth(test) + 1, " "), " | ", color=:white)
        time_str = @sprintf("%7.2f",resp[2])
        printstyled(lpad(time_str, elapsed_align, " "), " | ", color=:white)

        gpu_gc_str = @sprintf("%5.2f", resp[7])
        printstyled(lpad(gpu_gc_str, gc_align, " "), " | ", color=:white)
        # since there may be quite a few digits in the percentage,
        # the left-padding here is less to make sure everything fits
        gpu_percent_str = @sprintf("%4.1f", 100 * resp[7] / resp[2])
        printstyled(lpad(gpu_percent_str, percent_align, " "), " | ", color=:white)
        gpu_alloc_str = @sprintf("%5.2f", resp[6] / 2^20)
        printstyled(lpad(gpu_alloc_str, alloc_align, " "), " | ", color=:white)

        gpu_rss_str = ismissing(resp[10]) ? "N/A" : @sprintf("%5.2f", resp[10] / 2^20)
        printstyled(lpad(gpu_rss_str, rss_align, " "), " | ", color=:white)

        cpu_gc_str = @sprintf("%5.2f", resp[4])
        printstyled(lpad(cpu_gc_str, gc_align, " "), " | ", color=:white)
        cpu_percent_str = @sprintf("%4.1f", 100 * resp[4] / resp[2])
        printstyled(lpad(cpu_percent_str, percent_align, " "), " | ", color=:white)
        cpu_alloc_str = @sprintf("%5.2f", resp[3] / 2^20)
        printstyled(lpad(cpu_alloc_str, alloc_align, " "), " | ", color=:white)

        cpu_rss_str = @sprintf("%5.2f", resp[9] / 2^20)
        printstyled(lpad(cpu_rss_str, rss_align, " "), " |\n", color=:white)
    finally
        unlock(print_lock)
    end
end
global print_testworker_started = (name, wrkr)->begin
    if do_sanitize
        lock(print_lock)
        try
            printstyled(name, color=:white)
            printstyled(lpad("($wrkr)", name_align - textwidth(name) + 1, " "), " |",
                " "^elapsed_align, "started at $(now())\n", color=:white)
        finally
            unlock(print_lock)
        end
    end
end
function print_testworker_errored(name, wrkr)
    lock(print_lock)
    try
        printstyled(name, color=:red)
        printstyled(lpad("($wrkr)", name_align - textwidth(name) + 1, " "), " |",
            " "^elapsed_align, " failed at $(now())\n", color=:red)
    finally
        unlock(print_lock)
    end
end

# run tasks
results = []
all_tasks = Task[]
try
    # Monitor stdin and kill this task on ^C
    # but don't do this on Windows, because it may deadlock in the kernel
    t = current_task()
    running_tests = Dict{String, DateTime}()
    if !Sys.iswindows() && isa(stdin, Base.TTY)
        stdin_monitor = @async begin
            term = REPL.Terminals.TTYTerminal("xterm", stdin, stdout, stderr)
            try
                REPL.Terminals.raw!(term, true)
                while true
                    c = read(term, Char)
                    if c == '\x3'
                        Base.throwto(t, InterruptException())
                        break
                    elseif c == '?'
                        println("Currently running: ")
                        tests = sort(collect(running_tests), by=x->x[2])
                        foreach(tests) do (test, date)
                            println(test, " (running for ", round(now()-date, Minute), ")")
                        end
                    end
                end
            catch e
                isa(e, InterruptException) || rethrow()
            finally
                REPL.Terminals.raw!(term, false)
            end
        end
    end
    @sync begin
        for p in workers()
            @async begin
                push!(all_tasks, current_task())
                while length(tests) > 0
                    test = popfirst!(tests)
                    local resp
                    wrkr = p
                    snoop = do_snoop ? mktemp() : (nothing, nothing)

                    # tests that muck with the context should not be timed with CUDA events,
                    # since they won't be valid at the end of the test anymore.
                    time_source = in(test, ["initialization", "examples", "exceptions"]) ? :julia : :cuda

                    # run the test
                    running_tests[test] = now()
                    try
                        resp = remotecall_fetch(runtests, wrkr, test_runners[test], test, time_source, snoop[1])
                    catch e
                        isa(e, InterruptException) && return
                        resp = Any[e]
                    end
                    delete!(running_tests, test)
                    push!(results, (test, resp))

                    # act on the results
                    if resp[1] isa Exception
                        print_testworker_errored(test, wrkr)
                        do_quickfail && Base.throwto(t, InterruptException())

                        # the worker encountered some failure, recycle it
                        # so future tests get a fresh environment
                        rmprocs(wrkr, waitfor=30)
                        p = addworker(1)[1]
                    else
                        print_testworker_stats(test, wrkr, resp)
                    end

                    # aggregate the snooped compiler invocations
                    if do_snoop
                        for line in eachline(snoop[2])
                            println(snoop_io, line)
                        end
                        close(snoop[2])
                        rm(snoop[1])
                    end
                end
                if p != 1
                    # Free up memory =)
                    rmprocs(p, waitfor=30)
                end
            end
        end
    end
catch e
    isa(e, InterruptException) || rethrow()
    # If the test suite was merely interrupted, still print the
    # summary, which can be useful to diagnose what's going on
    foreach(task -> begin
            istaskstarted(task) || return
            istaskdone(task) && return
            try
                schedule(task, InterruptException(); error=true)
            catch ex
                @error "InterruptException" exception=ex,catch_backtrace()
            end
        end, all_tasks)
    for t in all_tasks
        # NOTE: we can't just wait, but need to discard the exception,
        #       because the throwto for --quickfail also kills the worker.
        try
            wait(t)
        catch e
            showerror(stderr, e)
        end
    end
finally
    if @isdefined stdin_monitor
        schedule(stdin_monitor, InterruptException(); error=true)
    end
end

# construct a testset to render the test results
o_ts = Test.DefaultTestSet("Overall")
Test.push_testset(o_ts)
completed_tests = Set{String}()
for (testname, (resp,)) in results
    push!(completed_tests, testname)
    if isa(resp, Test.DefaultTestSet)
        Test.push_testset(resp)
        Test.record(o_ts, resp)
        Test.pop_testset()
    elseif isa(resp, Tuple{Int,Int})
        fake = Test.DefaultTestSet(testname)
        for i in 1:resp[1]
            Test.record(fake, Test.Pass(:test, nothing, nothing, nothing))
        end
        for i in 1:resp[2]
            Test.record(fake, Test.Broken(:test, nothing))
        end
        Test.push_testset(fake)
        Test.record(o_ts, fake)
        Test.pop_testset()
    elseif isa(resp, RemoteException) && isa(resp.captured.ex, Test.TestSetException)
        println("Worker $(resp.pid) failed running test $(testname):")
        Base.showerror(stdout, resp.captured)
        println()
        fake = Test.DefaultTestSet(testname)
        for i in 1:resp.captured.ex.pass
            Test.record(fake, Test.Pass(:test, nothing, nothing, nothing))
        end
        for i in 1:resp.captured.ex.broken
            Test.record(fake, Test.Broken(:test, nothing))
        end
        for t in resp.captured.ex.errors_and_fails
            Test.record(fake, t)
        end
        Test.push_testset(fake)
        Test.record(o_ts, fake)
        Test.pop_testset()
    else
        if !isa(resp, Exception)
            resp = ErrorException(string("Unknown result type : ", typeof(resp)))
        end
        # If this test raised an exception that is not a remote testset exception,
        # i.e. not a RemoteException capturing a TestSetException that means
        # the test runner itself had some problem, so we may have hit a segfault,
        # deserialization errors or something similar.  Record this testset as Errored.
        fake = Test.DefaultTestSet(testname)
        Test.record(fake, Test.Error(:test_error, testname, nothing, Any[(resp, [])], LineNumberNode(1)))
        Test.push_testset(fake)
        Test.record(o_ts, fake)
        Test.pop_testset()
    end
end
for test in all_tests
    (test in completed_tests) && continue
    fake = Test.DefaultTestSet(test)
    Test.record(fake, Test.Error(:test_interrupted, test, nothing,
                                    [("skipped", [])], LineNumberNode(1)))
    Test.push_testset(fake)
    Test.record(o_ts, fake)
    Test.pop_testset()
end
println()
Test.print_test_results(o_ts, 1)
if !o_ts.anynonpass
    println("    \033[32;1mSUCCESS\033[0m")
else
    println("    \033[31;1mFAILURE\033[0m\n")
    Test.print_test_errors(o_ts)
    throw(Test.FallbackTestSetException("Test run finished with errors"))
end


text:
export SeisHdr

"""
SeisHdr: header information for seismic events

    S = SeisHdr()

Initialize an empty SeisHdr object. Fields can be initialized at creation with
keywords, e.g., SeisHdr(ot=DateTime("2012-01-03T03:49:45"), int=(0x02, "MMI")).

| Field | Default       | Type               | Meaning                        |
|:----  |:-----         |:-----              |:--------                       |
| id    | ""            | String             | Event ID                       |
| int   | (0x00, "")    | Tuple              | (Intensity, Intensity Scale)   |
| loc   | ()            | EQLoc              | Hypocenter data                |
| mag   | ()            | EQMag              | Magnitude data                 |
| misc  | ()            | Dict{String,Any}() | Non-essential info             |
| ot    | (unix epoch)  | DateTime           | Origin time                    |
| notes | []            | Array{String,1}    | Timestamped notes, logging     |
| src   | ""            | String             | Data source (URL/filename)     |
| typ   | ""            | String             | Event type                     |

See also: EQLoc, EQMag
"""
mutable struct SeisHdr
  id    ::String
  int   ::Tuple{UInt8,String}
  loc   ::EQLoc
  mag   ::EQMag
  misc  ::Dict{String,Any}
  notes ::Array{String,1}
  ot    ::DateTime
  src   ::String
  typ   ::String

  function SeisHdr(
    id    ::String,
    int   ::Tuple{UInt8,String},
    loc   ::EQLoc,
    mag   ::EQMag,
    misc  ::Dict{String,Any},
    notes ::Array{String,1},
    ot    ::DateTime,
    src   ::String,
    typ   ::String
    )

    return new(id, int, loc, mag, misc, notes, ot, src, typ)
  end
end
SeisHdr(;
          id    ::String                      = "",
          int   ::Tuple{UInt8,String}         = (0x00, ""),
          loc   ::EQLoc                       = EQLoc(),
          mag   ::EQMag                       = EQMag(),
          misc  ::Dict{String,Any}            = Dict{String,Any}(),
          notes ::Array{String,1}             = Array{String,1}(undef, 0),
          ot    ::DateTime                    = u2d(0),
          src   ::String                      = "",
          typ   ::String                      = "",
          ) = SeisHdr(id, int, loc, mag, misc, notes, ot, src, typ)

# =============================================================================
# Methods from Base
sizeof(H::SeisHdr) = sum([sizeof(getfield(H,i)) for i in fieldnames(SeisHdr)])

function isempty(H::SeisHdr)
  q = min(getfield(H, :ot)  == u2d(0),
          getfield(H, :int) == (0x00, ""))
  if q == true
    for f in (:id, :loc, :mag, :misc, :notes, :src, :typ)
      q = min(q, isempty(getfield(H, f)))
    end
  end
  return q
end

function isequal(H::SeisHdr, K::SeisHdr)
  q::Bool = true
  for i in fieldnames(SeisHdr)
    if i != :notes
      q = min(q, isequal(getfield(H,i), getfield(K,i)))
    end
  end
  return q
end
==(H1::SeisHdr, H2::SeisHdr) = isequal(H1, H2)

function write(io::IO, H::SeisHdr)
  write(io, Int64(sizeof(H.id)))
  write(io, H.id)
  write(io, H.int[1])
  write(io, Int64(sizeof(H.int[2])))
  write(io, H.int[2])
  write(io, H.loc)
  write(io, H.mag)
  write_misc(io, H.misc)
  write_string_vec(io, H.notes)
  write(io, round(Int64, d2u(getfield(H, :ot))*1.0e6))
  write(io, Int64(sizeof(H.src)))
  write(io, H.src)
  write(io, Int64(sizeof(H.typ)))
  write(io, H.typ)
  return nothing
end

read(io::IO, ::Type{SeisHdr}) = SeisHdr(
  String(read(io, read(io, Int64))),
  (read(io, UInt8), String(read(io, read(io, Int64)))),
  read(io, EQLoc),
  read(io, EQMag),
  read_misc(io, BUF.buf),
  read_string_vec(io, BUF.buf),
  u2d(read(io, Int64)*μs),
  String(read(io, read(io, Int64))),
  String(read(io, read(io, Int64)))
  )

summary(H::SeisHdr) = string(typeof(H), ", ",
  repr("text/plain", H.loc, context=:compact=>true), ", ",
  repr("text/plain", H.mag, context=:compact=>true), ", ",
  H.int[2], " ", H.int[1])

function show(io::IO, H::SeisHdr)
  W = max(80, displaysize(io)[2]-2)-show_os
  println(io, "    ID: ", H.id)
  println(io, "   INT: ", H.int[2], " ", H.int[1])
  println(io, "   LOC: ", repr("text/plain", H.loc, context=:compact=>true))
  println(io, "   MAG: ", repr("text/plain", H.mag, context=:compact=>true))
  println(io, "    OT: ", H.ot)
  println(io, "   SRC: ", str_trunc(H.src, W))
  println(io, "   TYP: ", str_trunc(H.typ, W))
  println(io, "  MISC: ", length(H.misc), " items")
  println(io, " NOTES: ", length(H.notes), " entries")
  return nothing
end
show(S::SeisHdr) = show(stdout, S)


text:
using KitBase, Plots
using KitBase.ProgressMeter: @showprogress
pyplot()

begin
    set = Setup(
        "gas",
        "cylinder",
        "2d2f",
        "kfvs",
        "bgk",
        1, # species
        2, # order of accuracy
        "vanleer", # limiter
        ["maxwell", "extra", "mirrror", "mirror"],
        0.5, # cfl
        10.0, # time
    )
    ps = CSpace2D(1.0, 6.0, 30, 0.0, π, 50, 1, 1)
    vs = VSpace2D(-10.0, 10.0, 48, -10.0, 10.0, 48)
    gas = Gas(5e-2, 3.0, 1.0, 1.0)

    prim0 = [1.0, 0.0, 0.0, 1.0]
    prim1 = [1.0, gas.Ma * sound_speed(1.0, gas.γ), 0.0, 1.0]
    fw = (args...) -> prim_conserve(prim1, gas.γ)
    ff = function (args...)
        prim = conserve_prim(fw(args...), gas.γ)
        h = maxwellian(vs.u, vs.v, prim)
        b = h .* gas.K / 2 / prim[end]
        return h, b
    end
    bc = function (x, y)
        if abs(x^2 + y^2 - 1) < 1e-3
            return prim0
        else
            return prim1
        end
    end
    ib = IB2F(fw, ff, bc)

    ks = SolverSet(set, ps, vs, gas, ib)
end

ctr, a1face, a2face = init_fvm(ks, ks.pSpace)

t = 0.0
dt = timestep(ks, ctr, 0.0)
nt = ks.set.maxTime ÷ dt |> Int
@showprogress for iter = 1:1000#nt
    evolve!(ks, ctr, a1face, a2face, dt; bc = [:maxwell, :extra, :mirror, :mirror])
    update!(
        ks,
        ctr,
        a1face,
        a2face,
        dt,
        zeros(4);
        bc = [:maxwell, :extra, :mirror, :mirror],
    )

    for j = ks.ps.nθ÷2+1:ks.ps.nθ
        ctr[ks.ps.nr+1, j].w .= ks.ib.fw(6, 0)
        ctr[ks.ps.nr+1, j].prim .= conserve_prim(ctr[ks.ps.nr+1, j].w, ks.gas.γ)
        ctr[ks.ps.nr+1, j].sw .= 0.0
        ctr[ks.ps.nr+1, j].h .= maxwellian(ks.vs.u, ks.vs.v, ctr[ks.ps.nr+1, j].prim)
        ctr[ks.ps.nr+1, j].b .=
            ctr[ks.ps.nr+1, j].h .* ks.gas.K ./ 2 ./ ctr[ks.ps.nr+1, j].prim[end]
    end

    global t += dt
end

begin
    sol = zeros(ks.ps.nr, ks.ps.nθ, 4)
    for i in axes(sol, 1), j in axes(sol, 2)
        sol[i, j, :] .= ctr[i, j].prim
        sol[i, j, end] = 1 / sol[i, j, end]
    end
    contourf(
        ps.x[1:ks.ps.nr, 1:ks.ps.nθ],
        ps.y[1:ks.ps.nr, 1:ks.ps.nθ],
        sol[:, :, 1],
        ratio = 1,
    )
end