text:
include("init.jl")

@testset "Filter Type Tests registration" begin

    broker = createBroker(2)
    fFed, fedinfo1 = createMessageFederate(1, "filter")
    mFed, fedinfo2 = createMessageFederate(1, "message")

    h.helicsFederateRegisterGlobalEndpoint(mFed, "port1", "")
    h.helicsFederateRegisterGlobalEndpoint(mFed, "port2", "")

    f1 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_CUSTOM, "filter1")
    h.helicsFilterAddSourceTarget(f1, "port1")

    f2 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_CUSTOM, "filter2")
    h.helicsFilterAddDestinationTarget(f2, "port2")

    @test f1 != f2

    ep1 = h.helicsFederateRegisterEndpoint(fFed, "fout", "")

    f3 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_CUSTOM, "c4")
    h.helicsFilterAddSourceTarget(f3, "Testfilter/fout")

    f1_b = h.helicsFederateGetFilter(fFed, "filter1")
    tmp = h.helicsFilterGetName(f1_b)
    @test tmp == "Testfilter/filter1"

    f1_c = h.helicsFederateGetFilterByIndex(fFed, 2)
    tmp = h.helicsFilterGetName(f1_c)
    @test tmp == "Testfilter/c4"

    @test_throws h.HELICSErrorInvalidArgument f1_n = h.helicsFederateGetFilterByIndex(fFed, -2)

    h.helicsFederateEnterExecutingModeAsync(fFed)
    h.helicsFederateEnterExecutingMode(mFed)
    h.helicsFederateEnterExecutingModeComplete(fFed)
    h.helicsFederateFinalizeAsync(mFed)
    h.helicsFederateFinalize(fFed)
    h.helicsFederateFinalizeComplete(mFed)

    state = h.helicsFederateGetState(fFed)
    @test state == h.HELICS_STATE_FINALIZE

    destroyFederate(fFed, fedinfo1)
    destroyFederate(mFed, fedinfo2)
    destroyBroker(broker)
end

@testset "Filter Type Tests info" begin

    broker = createBroker(2)
    fFed, fedinfo1 = createMessageFederate(1, "filter")
    mFed, fedinfo2 = createMessageFederate(1, "message")

    p1 = h.helicsFederateRegisterGlobalEndpoint(mFed, "port1", "")
    p2 = h.helicsFederateRegisterGlobalEndpoint(mFed, "port2", "")

    h.helicsEndpointSetInfo(p1, "p1_test")
    h.helicsEndpointSetInfo(p2, "p2_test")

    f1 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_CUSTOM, "filter1")
    h.helicsFilterAddSourceTarget(f1, "port1")
    h.helicsFilterSetInfo(f1, "f1_test")

    f2 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_CUSTOM, "filter2")
    h.helicsFilterAddDestinationTarget(f2, "port2")
    h.helicsFilterSetInfo(f2, "f2_test")

    ep1 = h.helicsFederateRegisterEndpoint(fFed, "fout", "")
    h.helicsEndpointSetInfo(ep1, "ep1_test")
    f3 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_CUSTOM, "c4")
    h.helicsFilterAddSourceTarget(f3, "filter0/fout")
    h.helicsFilterSetInfo(f3, "f3_test")

    @test h.helicsEndpointGetInfo(p1) == "p1_test"
    @test h.helicsEndpointGetInfo(p2) == "p2_test"
    @test h.helicsEndpointGetInfo(ep1) == "ep1_test"

    @test h.helicsFilterGetInfo(f1) == "f1_test"
    @test h.helicsFilterGetInfo(f2) == "f2_test"
    @test h.helicsFilterGetInfo(f3) == "f3_test"

    h.helicsFederateEnterExecutingModeAsync(fFed)
    h.helicsFederateEnterExecutingMode(mFed)
    h.helicsFederateEnterExecutingModeComplete(fFed)

    h.helicsFederateFinalizeAsync(mFed)
    h.helicsFederateFinalize(fFed)
    h.helicsFederateFinalizeComplete(mFed)

    destroyFederate(fFed, fedinfo1)
    destroyFederate(mFed, fedinfo2)
    destroyBroker(broker)

end

@testset "Filter Types tests core fitler registration" begin

    core1 = h.helicsCreateCore("inproc", "core1", "--autobroker")

    core2 = h.helicsCoreClone(core1)

    core1IdentifierString = h.helicsCoreGetIdentifier(core1)

    @test core1IdentifierString == "core1"

    sourceFilter1 = h.helicsCoreRegisterFilter(core1, h.HELICS_FILTER_TYPE_DELAY, "core1SourceFilter")

    h.helicsFilterAddSourceTarget(sourceFilter1, "ep1")
    destinationFilter1 = h.helicsCoreRegisterFilter(core1, h.HELICS_FILTER_TYPE_DELAY, "core1DestinationFilter")

    h.helicsFilterAddDestinationTarget(destinationFilter1, "ep2")
    cloningFilter1 = h.helicsCoreRegisterCloningFilter(core1, "ep3")

    h.helicsFilterRemoveDeliveryEndpoint(cloningFilter1, "ep3")
    core1IsConnected = h.helicsCoreIsConnected(core1)
    @test core1IsConnected != 0
    h.helicsCoreSetReadyToInit(core1)
    h.helicsCoreDisconnect(core1)
    h.helicsCoreDisconnect(core2)
    h.helicsCoreFree(core1)
    h.helicsCoreFree(core2)
    h.helicsCloseLibrary()

end

@testset "Filter Type Tests message filter function" begin

    broker = createBroker(2)
    fFed, fedinfo1 = createMessageFederate(1, "filter")
    mFed, fedinfo2 = createMessageFederate(1, "message")

    h.helicsFederateSetFlagOption(mFed, h.HELICS_FLAG_IGNORE_TIME_MISMATCH_WARNINGS, true)
    p1 = h.helicsFederateRegisterGlobalEndpoint(mFed, "port1", "")
    p2 = h.helicsFederateRegisterGlobalEndpoint(mFed, "port2", "")

    f1 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_DELAY, "filter1")
    h.helicsFilterAddSourceTarget(f1, "port1")
    h.helicsFilterSet(f1, "delay", 2.5)

    h.helicsFederateEnterExecutingModeAsync(fFed)
    h.helicsFederateEnterExecutingMode(mFed)
    h.helicsFederateEnterExecutingModeComplete(fFed)

    state = h.helicsFederateGetState(fFed)
    @test state == h.HELICS_STATE_EXECUTION
    data = repeat('a', 500)
    h.helicsEndpointSendMessageRaw(p1, "port2", data)

    h.helicsFederateRequestTimeAsync(mFed, 1.0)
    h.helicsFederateRequestTime(fFed, 1.0)
    h.helicsFederateRequestTimeComplete(mFed)

    @test h.helicsFederateHasMessage(mFed) == false

    h.helicsFederateRequestTimeAsync(mFed, 2.0)
    h.helicsFederateRequestTime(fFed, 2.0)
    h.helicsFederateRequestTimeComplete(mFed)
    @test h.helicsEndpointHasMessage(p2) == false

    h.helicsFederateRequestTimeAsync(fFed, 3.0)
    h.helicsFederateRequestTime(mFed, 3.0)

    @test h.helicsEndpointHasMessage(p2) == true

    m2 = h.helicsEndpointGetMessageObject(p2)
    @test h.helicsMessageGetSource(m2) == "port1"
    @test h.helicsMessageGetOriginalSource(m2) == "port1"
    @test h.helicsMessageGetDestination(m2) == "port2"
    @test h.helicsMessageGetRawDataSize(m2) == length(data)
    @test h.helicsMessageGetTime(m2) == 2.5

    h.helicsFederateRequestTime(mFed, 3.0)
    h.helicsFederateRequestTimeComplete(fFed)
    h.helicsFederateFinalizeAsync(mFed)
    h.helicsFederateFinalize(fFed)
    h.helicsFederateFinalizeComplete(mFed)
    state = h.helicsFederateGetState(fFed)
    @test state == h.HELICS_STATE_FINALIZE

    destroyFederate(fFed, fedinfo1)
    destroyFederate(mFed, fedinfo2)
    destroyBroker(broker)

end

@testset "Filter test types function mObj" begin

    broker = createBroker(2)
    fFed, fedinfo1 = createMessageFederate(1, "filter")
    mFed, fedinfo2 = createMessageFederate(1, "message")

    h.helicsFederateSetFlagOption(mFed, h.HELICS_FLAG_IGNORE_TIME_MISMATCH_WARNINGS, true)
    p1 = h.helicsFederateRegisterGlobalEndpoint(mFed, "port1", "")
    p2 = h.helicsFederateRegisterGlobalEndpoint(mFed, "port2", "")

    f1 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_DELAY, "filter1")
    h.helicsFilterAddSourceTarget(f1, "port1")
    h.helicsFilterSet(f1, "delay", 2.5,)

    h.helicsFederateEnterExecutingModeAsync(fFed)
    h.helicsFederateEnterExecutingMode(mFed)
    h.helicsFederateEnterExecutingModeComplete(fFed)

    state = h.helicsFederateGetState(fFed)
    @test state == h.HELICS_STATE_EXECUTION

    data = repeat('a', 500)
    h.helicsEndpointSendMessageRaw(p1, "port2", data)

    h.helicsFederateRequestTimeAsync(mFed, 1.0)
    h.helicsFederateRequestTime(fFed, 1.0)
    h.helicsFederateRequestTimeComplete(mFed)

    res = h.helicsFederateHasMessage(mFed)
    @test res == false

    h.helicsFederateRequestTimeAsync(mFed, 2.0)
    h.helicsFederateRequestTime(fFed, 2.0)
    h.helicsFederateRequestTimeComplete(mFed)
    @test h.helicsEndpointHasMessage(p2) == false

    h.helicsFederateRequestTimeAsync(fFed, 3.0)
    h.helicsFederateRequestTime(mFed, 3.0)

    @test h.helicsEndpointHasMessage(p2) == true

    m2 = h.helicsEndpointGetMessageObject(p2)
    @test h.helicsMessageGetSource(m2) == "port1"
    @test h.helicsMessageGetOriginalSource(m2) == "port1"
    @test h.helicsMessageGetDestination(m2) == "port2"
    @test h.helicsMessageGetRawDataSize(m2) == length(data)
    @test h.helicsMessageGetTime(m2) == 2.5

    h.helicsFederateRequestTime(mFed, 3.0)
    h.helicsFederateRequestTimeComplete(fFed)
    h.helicsFederateFinalizeAsync(mFed)
    h.helicsFederateFinalize(fFed)
    h.helicsFederateFinalizeComplete(mFed)
    state = h.helicsFederateGetState(fFed)
    @test state == h.HELICS_STATE_FINALIZE

    destroyFederate(fFed, fedinfo1)
    destroyFederate(mFed, fedinfo2)
    destroyBroker(broker)

end

@testset "Filter test types function two stage" begin

    broker = createBroker(3)
    fFed, fedinfo1 = createMessageFederate(1, "filter")
    fFed2, fedinfo2 = createMessageFederate(1, "filter2")
    mFed, fedinfo3 = createMessageFederate(1, "message")

    p1 = h.helicsFederateRegisterGlobalEndpoint(mFed, "port1", "")
    p2 = h.helicsFederateRegisterGlobalEndpoint(mFed, "port2", "")

    f1 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_DELAY, "filter1")
    h.helicsFilterAddSourceTarget(f1, "port1")
    h.helicsFilterSet(f1, "delay", 1.25)

    f2 = h.helicsFederateRegisterFilter(fFed2, h.HELICS_FILTER_TYPE_DELAY, "filter2")
    h.helicsFilterAddSourceTarget(f2, "port1")
    h.helicsFilterSet(f2, "delay", 1.25)

    h.helicsFederateEnterExecutingModeAsync(fFed)
    h.helicsFederateEnterExecutingModeAsync(fFed2)
    h.helicsFederateEnterExecutingMode(mFed)
    h.helicsFederateEnterExecutingModeComplete(fFed)
    h.helicsFederateEnterExecutingModeComplete(fFed2)

    state = h.helicsFederateGetState(fFed)
    @test state == h.HELICS_STATE_EXECUTION
    data = repeat('a', 500)
    h.helicsEndpointSendMessageRaw(p1, "port2", data)

    h.helicsFederateRequestTimeAsync(mFed, .0)
    h.helicsFederateRequestTimeAsync(fFed, 1.0)
    h.helicsFederateRequestTime(fFed2, 1.0)
    h.helicsFederateRequestTimeComplete(mFed)
    h.helicsFederateRequestTimeComplete(fFed)
    @test h.helicsFederateHasMessage(mFed) == false

    h.helicsFederateRequestTimeAsync(mFed, .0)
    h.helicsFederateRequestTimeAsync(fFed2, 2.0)
    h.helicsFederateRequestTime(fFed, 2.0)
    h.helicsFederateRequestTimeComplete(mFed)
    h.helicsFederateRequestTimeComplete(fFed2)
    @test h.helicsEndpointHasMessage(p2) == false

    h.helicsFederateRequestTimeAsync(fFed, 3.0)
    h.helicsFederateRequestTimeAsync(fFed2, 3.0)
    h.helicsFederateRequestTime(mFed, 3.0)
    @test h.helicsEndpointHasMessage(p2) == true

    m2 = h.helicsEndpointGetMessageObject(p2)
    @test h.helicsMessageGetSource(m2) == "port1"
    @test h.helicsMessageGetOriginalSource(m2) == "port1"
    @test h.helicsMessageGetDestination(m2) == "port2"
    @test h.helicsMessageGetRawDataSize(m2) == length(data)
    @test h.helicsMessageGetTime(m2) == 2.5

    h.helicsFederateRequestTimeComplete(fFed)
    h.helicsFederateRequestTimeComplete(fFed2)
    h.helicsFederateFinalizeAsync(mFed)
    h.helicsFederateFinalizeAsync(fFed)
    h.helicsFederateFinalize(fFed2)
    h.helicsFederateFinalizeComplete(mFed)
    h.helicsFederateFinalizeComplete(fFed)
    state = h.helicsFederateGetState(fFed)
    @test state == h.HELICS_STATE_FINALIZE

    destroyFederate(fFed, fedinfo1)
    destroyFederate(fFed2, fedinfo2)
    destroyFederate(mFed, fedinfo3)
    destroyBroker(broker)

end

@testset "Filter test types function2" begin

    broker = createBroker(2)
    fFed, fedinfo1 = createMessageFederate(1, "filter")
    mFed, fedinfo2 = createMessageFederate(1, "message")

    p1 = h.helicsFederateRegisterGlobalEndpoint(mFed, "port1", "")
    p2 = h.helicsFederateRegisterGlobalEndpoint(mFed, "port2", "")

    f1 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_DELAY, "filter1")
    h.helicsFilterAddSourceTarget(f1, "port1")
    h.helicsFilterSet(f1, "delay", 2.5)

    f2 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_DELAY, "filter2")
    h.helicsFilterAddSourceTarget(f2, "port2")
    h.helicsFilterSet(f2, "delay", 2.5)
    # this is expected to fail since a regular filter doesn't have a delivery endpoint
    @test_throws h.HELICSErrorInvalidObject h.helicsFilterAddDeliveryEndpoint(f2, "port1")

    h.helicsFederateEnterExecutingModeAsync(fFed)
    h.helicsFederateEnterExecutingMode(mFed)
    h.helicsFederateEnterExecutingModeComplete(fFed)

    state = h.helicsFederateGetState(fFed)
    @test state == h.HELICS_STATE_EXECUTION

    data = repeat('a', 500)
    h.helicsEndpointSendMessageRaw(p1, "port2", data)

    h.helicsFederateRequestTimeAsync(mFed, 1.0)
    h.helicsFederateRequestTime(fFed, 1.0)
    h.helicsFederateRequestTimeComplete(mFed)

    res = h.helicsFederateHasMessage(mFed)
    @test res == false

    h.helicsEndpointSendMessageRaw(p2, "port1", data)
    h.helicsFederateRequestTimeAsync(mFed, 2.0)
    h.helicsFederateRequestTime(fFed, 2.0)
    h.helicsFederateRequestTimeComplete(mFed)
    @test h.helicsEndpointHasMessage(p2) == false
    # there may be something wrong here yet but this test isn't the one to find it and
    # this may prevent spurious errors for now.
    # std::this_thread::yield()
    h.helicsFederateRequestTime(mFed, 3.0)

    @test h.helicsEndpointHasMessage(p2) == true

    m2 = h.helicsEndpointGetMessageObject(p2)
    @test h.helicsMessageGetSource(m2) == "port1"
    @test h.helicsMessageGetOriginalSource(m2) == "port1"
    @test h.helicsMessageGetDestination(m2) == "port2"
    @test h.helicsMessageGetRawDataSize(m2) == length(data)
    @test h.helicsMessageGetTime(m2) == 2.5

    @test h.helicsEndpointHasMessage(p1) == false

    h.helicsFederateRequestTime(mFed, 4.0)
    @test h.helicsEndpointHasMessage(p1) == true
    h.helicsFederateFinalizeAsync(mFed)
    h.helicsFederateFinalize(fFed)
    h.helicsFederateFinalizeComplete(mFed)
    state = h.helicsFederateGetState(fFed)
    @test state == h.HELICS_STATE_FINALIZE

    destroyFederate(fFed, fedinfo1)
    destroyFederate(mFed, fedinfo2)
    destroyBroker(broker)
end

@testset "Filter test types message filter function3" begin

    broker = createBroker(2)
    fFed, fedinfo1 = createMessageFederate(1, "filter", 1)
    mFed, fedinfo2 = createMessageFederate(1, "message", 1)

    p1 = h.helicsFederateRegisterGlobalEndpoint(mFed, "port1", "")
    p2 = h.helicsFederateRegisterGlobalEndpoint(mFed, "port2", "random")

    f1 = h.helicsFederateRegisterGlobalFilter(fFed, h.HELICS_FILTER_TYPE_CUSTOM, "filter1")
    h.helicsFilterAddSourceTarget(f1, "port1")
    f2 = h.helicsFederateRegisterGlobalFilter(fFed, h.HELICS_FILTER_TYPE_DELAY, "filter2")
    h.helicsFilterAddSourceTarget(f2, "port1")

    h.helicsFederateRegisterEndpoint(fFed, "fout", "")
    f3 = h.helicsFederateRegisterFilter(fFed, h.HELICS_FILTER_TYPE_RANDOM_DELAY, "filter3")
    h.helicsFilterAddSourceTarget(f3, "filter0/fout")

    h.helicsFilterSet(f2, "delay", 2.5)

    h.helicsFederateEnterExecutingModeAsync(fFed)
    h.helicsFederateEnterExecutingMode(mFed)
    h.helicsFederateEnterExecutingModeComplete(fFed)

    state = h.helicsFederateGetState(fFed)
    @test state == h.HELICS_STATE_EXECUTION

    data = "hello world"
    h.helicsEndpointSendMessageRaw(p1, "port2", data)

    h.helicsFederateRequestTimeAsync(mFed, 1.0)
    h.helicsFederateRequestTime(fFed, 1.0)
    h.helicsFederateRequestTimeComplete(mFed)

    @test h.helicsFederateHasMessage(mFed) == false

    h.helicsEndpointSendMessageRaw(p2, "port1", data)
    h.helicsFederateRequestTimeAsync(mFed, 2.0)
    h.helicsFederateRequestTime(fFed, 2.0)
    h.helicsFederateRequestTimeComplete(mFed)
    @test h.helicsEndpointHasMessage(p2) == false
    # there may be something wrong here yet but this test isn't the one to find it and
    # this may prevent spurious errors for now.
    # std::this_thread::yield()
    h.helicsFederateRequestTimeAsync(mFed, 3.0)
    h.helicsFederateRequestTime(fFed, 3.0)
    h.helicsFederateRequestTimeComplete(mFed)

    @test h.helicsEndpointHasMessage(p2)

    m2 = h.helicsEndpointGetMessageObject(p2)
    @test h.helicsMessageGetSource(m2) == "port1"
    @test h.helicsMessageGetOriginalSource(m2) == "port1"
    @test h.helicsMessageGetDestination(m2) == "port2"
    @test h.helicsMessageGetRawDataSize(m2) == length(data)
    @test h.helicsMessageGetTime(m2) == 2.5

    @test h.helicsEndpointHasMessage(p1) == true
    h.helicsFederateFinalize(mFed)
    h.helicsFederateFinalize(fFed)
    state = h.helicsFederateGetState(fFed)
    @test state == h.HELICS_STATE_FINALIZE

    destroyFederate(mFed, fedinfo1)
    destroyFederate(fFed, fedinfo2)
    destroyBroker(broker)

end

@testset "Filter test types clone test" begin

    broker = createBroker(3)
    sFed, fedinfo1 = createMessageFederate(1, "source", 1)
    dFed, fedinfo2 = createMessageFederate(1, "dest", 1)
    dcFed, fedinfo3 = createMessageFederate(1, "dest_clone", 1)

    p1 = h.helicsFederateRegisterGlobalEndpoint(sFed, "src", "")
    p2 = h.helicsFederateRegisterGlobalEndpoint(dFed, "dest", "")
    p3 = h.helicsFederateRegisterGlobalEndpoint(dcFed, "cm", "")

    f1 = h.helicsFederateRegisterCloningFilter(dcFed, "")
    h.helicsFilterAddDeliveryEndpoint(f1, "cm")
    h.helicsFilterAddSourceTarget(f1, "src")

    h.helicsFederateEnterExecutingModeAsync(sFed)
    h.helicsFederateEnterExecutingModeAsync(dcFed)
    h.helicsFederateEnterExecutingMode(dFed)
    h.helicsFederateEnterExecutingModeComplete(sFed)
    h.helicsFederateEnterExecutingModeComplete(dcFed)

    state = h.helicsFederateGetState(sFed)
    @test state == h.HELICS_STATE_EXECUTION
    state = h.helicsFederateGetState(dcFed)
    @test state == h.HELICS_STATE_EXECUTION
    state = h.helicsFederateGetState(dFed)
    @test state == h.HELICS_STATE_EXECUTION

    data = repeat('a', 500)
    h.helicsEndpointSendMessageRaw(p1, "dest", data)

    h.helicsFederateRequestTimeAsync(sFed, 1.0)
    h.helicsFederateRequestTimeAsync(dcFed, 1.0)
    h.helicsFederateRequestTime(dFed, 1.0)
    h.helicsFederateRequestTimeComplete(sFed)
    h.helicsFederateRequestTimeComplete(dcFed)

    @test h.helicsFederateHasMessage(dFed)

    m2 = h.helicsEndpointGetMessageObject(p2)
    @test h.helicsMessageGetSource(m2) == "src"
    @test h.helicsMessageGetOriginalSource(m2) == "src"
    @test h.helicsMessageGetDestination(m2) == "dest"
    @test h.helicsMessageGetRawDataSize(m2) == length(data)

    @test h.helicsFederateHasMessage(dcFed)

    m2 = h.helicsEndpointGetMessageObject(p3)
    @test h.helicsMessageGetSource(m2) == "src"
    @test h.helicsMessageGetOriginalSource(m2) == "src"
    @test h.helicsMessageGetDestination(m2) == "cm"
    @test h.helicsMessageGetOriginalDestination(m2) == "dest"
    @test h.helicsMessageGetRawDataSize(m2) == length(data)

    h.helicsFederateFinalizeAsync(sFed)
    h.helicsFederateFinalizeAsync(dFed)
    h.helicsFederateFinalize(dcFed)
    h.helicsFederateFinalizeComplete(sFed)
    h.helicsFederateFinalizeComplete(dFed)
    state = h.helicsFederateGetState(sFed)
    @test state == h.HELICS_STATE_FINALIZE

    destroyFederate(sFed, fedinfo1)
    destroyFederate(dFed, fedinfo2)
    destroyFederate(dcFed, fedinfo3)
    destroyBroker(broker)

end

@testset "Filter test types clone test connections" begin

    broker = createBroker(3)
    sFed, fedinfo1 = createMessageFederate(1, "source", 1.0, )
    dFed, fedinfo2 = createMessageFederate(1, "dest", 1.0)
    dcFed, fedinfo3 = createMessageFederate(1, "dest_clone", 1.0)

    p1 = h.helicsFederateRegisterGlobalEndpoint(sFed, "src", "")
    p2 = h.helicsFederateRegisterGlobalEndpoint(dFed, "dest", "")
    p3 = h.helicsFederateRegisterGlobalEndpoint(dcFed, "cm", "")

    f1 = h.helicsFederateRegisterGlobalCloningFilter(dcFed, "filt1")
    h.helicsFilterAddDeliveryEndpoint(f1, "cm")
    cr = h.helicsFederateGetCoreObject(sFed)

    h.helicsCoreAddSourceFilterToEndpoint(cr, "filt1", "src")
    h.helicsCoreAddSourceFilterToEndpoint(cr, "", "src")

    h.helicsFederateEnterExecutingModeAsync(sFed)
    h.helicsFederateEnterExecutingModeAsync(dcFed)
    h.helicsFederateEnterExecutingMode(dFed)
    h.helicsFederateEnterExecutingModeComplete(sFed)
    h.helicsFederateEnterExecutingModeComplete(dcFed)

    q = h.helicsCreateQuery("", "filtered_endpoints")
    filteredEndpoints = h.helicsQueryExecute(q, sFed)
    occursin("srcFilters", filteredEndpoints)
    occursin("(cloning)", filteredEndpoints)
    h.helicsQueryFree(q)

    state = h.helicsFederateGetState(sFed)
    @test state == h.HELICS_STATE_EXECUTION
    data = repeat('a', 500)
    h.helicsEndpointSendMessageRaw(p1, "dest", data)

    h.helicsFederateRequestTimeAsync(sFed, 1.0)
    h.helicsFederateRequestTimeAsync(dcFed, 1.0)
    h.helicsFederateRequestTime(dFed, 1.0)
    h.helicsFederateRequestTimeComplete(sFed)
    h.helicsFederateRequestTimeComplete(dcFed)

    @test h.helicsFederateHasMessage(dFed) == true

    m2 = h.helicsEndpointGetMessageObject(p2)
    @test h.helicsMessageGetSource(m2) == "src"
    @test h.helicsMessageGetOriginalSource(m2) == "src"
    @test h.helicsMessageGetDestination(m2) == "dest"
    @test h.helicsMessageGetRawDataSize(m2) == length(data)

    @test h.helicsFederateHasMessage(dcFed) == true

    m2 = h.helicsEndpointGetMessageObject(p3)
    @test h.helicsMessageGetSource(m2) == "src"
    @test h.helicsMessageGetOriginalSource(m2) == "src"
    @test h.helicsMessageGetDestination(m2) == "cm"
    @test h.helicsMessageGetOriginalDestination(m2) == "dest"
    @test h.helicsMessageGetRawDataSize(m2) == length(data)

    h.helicsFederateFinalizeAsync(sFed)
    h.helicsFederateFinalizeAsync(dFed)
    h.helicsFederateFinalize(dcFed)
    h.helicsFederateFinalizeComplete(sFed)
    h.helicsFederateFinalizeComplete(dFed)
    state = h.helicsFederateGetState(sFed)
    @test state == h.HELICS_STATE_FINALIZE

    destroyFederate(sFed, fedinfo1)
    destroyFederate(dFed, fedinfo2)
    destroyFederate(dcFed, fedinfo3)
    destroyBroker(broker)

end

@testset "Filter test types clone test broker connections" begin
    broker = createBroker(3)
    sFed, fedinfo1 = createMessageFederate(1, "source", 1.0)
    dFed, fedinfo2 = createMessageFederate(1, "dest", 1.0)
    dcFed, fedinfo3 = createMessageFederate(1, "dest_clone", 1.0)

    p1 = h.helicsFederateRegisterGlobalEndpoint(sFed, "src", "")
    p2 = h.helicsFederateRegisterGlobalEndpoint(dFed, "dest", "")
    p3 = h.helicsFederateRegisterGlobalEndpoint(dcFed, "cm", "")

    f1 = h.helicsFederateRegisterGlobalCloningFilter(dcFed, "filt1")
    h.helicsFilterAddDeliveryEndpoint(f1, "cm")

    h.helicsBrokerAddSourceFilterToEndpoint(broker, "filt1", "src")
    h.helicsBrokerAddSourceFilterToEndpoint(broker, "", "src")

    h.helicsFederateEnterExecutingModeAsync(sFed)
    h.helicsFederateEnterExecutingModeAsync(dcFed)
    h.helicsFederateEnterExecutingMode(dFed)
    h.helicsFederateEnterExecutingModeComplete(sFed)
    h.helicsFederateEnterExecutingModeComplete(dcFed)

    state = h.helicsFederateGetState(sFed)
    @test state == h.HELICS_STATE_EXECUTION
    data = repeat('a', 500)
    h.helicsEndpointSendMessageRaw(p1, "dest", data)

    h.helicsFederateRequestTimeAsync(sFed, 1.0)
    h.helicsFederateRequestTimeAsync(dcFed, 1.0)
    h.helicsFederateRequestTime(dFed, 1.0)
    h.helicsFederateRequestTimeComplete(sFed)
    h.helicsFederateRequestTimeComplete(dcFed)

    @test h.helicsFederateHasMessage(dFed) == true

    m2 = h.helicsEndpointGetMessageObject(p2)
    @test h.helicsMessageGetSource(m2) == "src"
    @test h.helicsMessageGetOriginalSource(m2) == "src"
    @test h.helicsMessageGetDestination(m2) == "dest"
    @test h.helicsMessageGetRawDataSize(m2) == length(data)

    @test h.helicsFederateHasMessage(dcFed) == true

    m2 = h.helicsEndpointGetMessageObject(p3)
    @test h.helicsMessageGetSource(m2) == "src"
    @test h.helicsMessageGetOriginalSource(m2) == "src"
    @test h.helicsMessageGetDestination(m2) == "cm"
    @test h.helicsMessageGetOriginalDestination(m2) == "dest"
    @test h.helicsMessageGetRawDataSize(m2) == length(data)

    h.helicsFederateFinalizeAsync(sFed)
    h.helicsFederateFinalizeAsync(dFed)
    h.helicsFederateFinalize(dcFed)
    h.helicsFederateFinalizeComplete(sFed)
    h.helicsFederateFinalizeComplete(dFed)
    state = h.helicsFederateGetState(sFed)
    @test state == h.HELICS_STATE_FINALIZE

    destroyFederate(sFed, fedinfo1)
    destroyFederate(dFed, fedinfo2)
    destroyFederate(dcFed, fedinfo3)
    destroyBroker(broker)

end

@testset "Filter test types clone test dest connections" begin
    broker = createBroker(3)
    sFed, fedinfo1 = createMessageFederate(1, "source", 1.0)
    dFed, fedinfo2 = createMessageFederate(1, "dest", 1.0)
    dcFed, fedinfo3 = createMessageFederate(1, "dest_clone", 2.0)

    p1 = h.helicsFederateRegisterGlobalEndpoint(sFed, "src", "")
    p2 = h.helicsFederateRegisterGlobalEndpoint(dFed, "dest", "")
    p3 = h.helicsFederateRegisterGlobalEndpoint(dcFed, "cm", "")

    f1 = h.helicsFederateRegisterGlobalCloningFilter(dcFed, "filt1")
    h.helicsFilterAddDeliveryEndpoint(f1, "cm")

    cr = h.helicsFederateGetCoreObject(sFed)

    h.helicsCoreAddDestinationFilterToEndpoint(cr, "filt1", "dest")

    h.helicsCoreAddDestinationFilterToEndpoint(cr, "", "dest")

    h.helicsCoreFree(cr)

    h.helicsFederateEnterExecutingModeAsync(sFed)
    h.helicsFederateEnterExecutingModeAsync(dcFed)
    h.helicsFederateEnterExecutingMode(dFed)
    h.helicsFederateEnterExecutingModeComplete(sFed)
    h.helicsFederateEnterExecutingModeComplete(dcFed)

    q = h.helicsCreateQuery("", "filtered_endpoints")
    filteredEndpoints = h.helicsQueryExecute(q, dFed)
    @test occursin("cloningdestFilter", filteredEndpoints)
    h.helicsQueryFree(q)

    state = h.helicsFederateGetState(sFed)
    @test state == h.HELICS_STATE_EXECUTION
    data = repeat('a', 500)
    h.helicsEndpointSendMessageRaw(p1, "dest", data)

    h.helicsFederateFinalize(sFed)

    # TODO: implement threading

    # auto dFedExec = [&]() {
    #     h.helicsFederateRequestTime(dFed, 1.0)
    #     m2 = h.helicsEndpointGetMessageObject(p2)
    #     h.helicsFederateFinalize(dFed, "")
    # }

    # h.helics_message m3
    # auto dcFedExec = [&]() {
    #     h.helicsFederateRequestTime(dcFed, 2.0)
    #     auto res = h.helicsFederateHasMessage(dcFed)
    #     if (res == h.helics_false) {
    #         std::this_thread::sleep_for(std::chrono::milliseconds(50))
    #         h.helicsFederateRequestTime(dcFed, 4.0)
    #     }
    #     m3 = h.helicsEndpointGetMessageObject(p3)
    #     h.helicsFederateFinalize(dcFed)
    # }

    # auto threaddFed = std::thread(dFedExec)
    # auto threaddcFed = std::thread(dcFedExec)

    # threaddFed.join()
    # (m2.source, "src")
    # (m2.original_source, "src")
    # (m2.dest, "dest")
    # (m2.length, static_cast<int64_t>(data.size()))

    # threaddcFed.join()

    # (m3.source, "src")
    # (m3.original_source, "src")
    # (m3.dest, "cm")
    # (m3.original_dest, "dest")
    # (m3.length, static_cast<int64_t>(data.size()))

    # (state = h.helicsFederateGetState(sFed))
    # (state == h.helics_state_finalize)

    destroyFederate(sFed, fedinfo1)
    destroyFederate(dFed, fedinfo2)
    destroyFederate(dcFed, fedinfo3)
    destroyBroker(broker)

end

@testset "Filter test types clone test broker dest connections" begin

    broker = createBroker(3)
    sFed, fedinfo1 = createMessageFederate(1, "source", 1.0)
    dFed, fedinfo2 = createMessageFederate(1, "dest", 1.0)
    dcFed, fedinfo3 = createMessageFederate(1, "dest_clone", 1.0)

    p1 = h.helicsFederateRegisterGlobalEndpoint(sFed, "src", "")
    p2 = h.helicsFederateRegisterGlobalEndpoint(dFed, "dest", "")
    p3 = h.helicsFederateRegisterGlobalEndpoint(dcFed, "cm", "")

    f1 = h.helicsFederateRegisterGlobalCloningFilter(dcFed, "filt1")
    h.helicsFilterAddDeliveryEndpoint(f1, "cm")
    h.helicsBrokerAddDestinationFilterToEndpoint(broker, "filt1", "dest")

    h.helicsBrokerAddDestinationFilterToEndpoint(broker, "", "dest")

    h.helicsFederateEnterExecutingModeAsync(sFed)
    h.helicsFederateEnterExecutingModeAsync(dcFed)
    h.helicsFederateEnterExecutingMode(dFed)
    h.helicsFederateEnterExecutingModeComplete(sFed)
    h.helicsFederateEnterExecutingModeComplete(dcFed)

    state = h.helicsFederateGetState(sFed)
    @test state == h.HELICS_STATE_EXECUTION
    data = repeat('a', 500)
    h.helicsEndpointSendMessageRaw(p1, "dest", data)

    h.helicsFederateRequestTimeAsync(sFed, 1.0)
    h.helicsFederateRequestTimeAsync(dcFed, 1.0)
    h.helicsFederateRequestTime(dFed, 1.0)
    h.helicsFederateRequestTimeComplete(sFed)
    h.helicsFederateRequestTimeComplete(dcFed)

    @test h.helicsFederateHasMessage(dFed) == true

    m2 = h.helicsEndpointGetMessageObject(p2)

    @test h.helicsMessageGetSource(m2) == "src"
    @test h.helicsMessageGetOriginalSource(m2) == "src"
    @test h.helicsMessageGetDestination(m2) == "dest"
    @test h.helicsMessageGetRawDataSize(m2) == length(data)

    h.helicsFederateFinalizeAsync(sFed)
    h.helicsFederateFinalizeAsync(dFed)

    # TODO: figure out why this test fails on CI
    @test_broken false
    @show h.helicsFederateHasMessage(dcFed) == false

    h.helicsFederateRequestTime(dcFed, 2.0)

    @test h.helicsFederateHasMessage(dcFed) == true

    m2 = h.helicsEndpointGetMessageObject(p3)
    @test h.helicsMessageGetSource(m2) == "src"
    @test h.helicsMessageGetOriginalSource(m2) == "src"
    @test h.helicsMessageGetDestination(m2) == "cm"
    @test h.helicsMessageGetOriginalDestination(m2) == "dest"
    @test h.helicsMessageGetRawDataSize(m2) == length(data)

    res2 = h.helicsFederateHasMessage(dcFed)

    h.helicsFederateFinalize(dcFed)
    h.helicsFederateFinalizeComplete(sFed)
    h.helicsFederateFinalizeComplete(dFed)
    state = h.helicsFederateGetState(sFed)
    state == h.HELICS_STATE_FINALIZE

    destroyFederate(sFed, fedinfo1)
    destroyFederate(dFed, fedinfo2)
    destroyFederate(dcFed, fedinfo3)
    destroyBroker(broker)

end

@testset "Filter test file load" begin

    filename = joinpath(@__DIR__, "filters.json")
    mFed = h.helicsCreateMessageFederateFromConfig(filename)

    name = h.helicsFederateGetName(mFed)
    @test name == "filterFed"

    @test h.helicsFederateGetEndpointCount(mFed) == 3
    h.helicsFederateFinalize(mFed)
    h.helicsFederateFree(mFed)

end


text:
using Test

import Cn2An:cn2an
import Cn2An:an2cn

@test cn2an("九千八百七十六万五千四百三十二亿九千八百七十六万五千四百三十二") == 9876543298765432
@test an2cn(9876543298765432) == "九千八百七十六万五千四百三十二亿九千八百七十六万五千四百三十二"

ext: jl


text:
using DifferentialEquations
# using CoolDicts

struct NamedViewVector{T} <: AbstractVector{T}
    _namedtuple
    _vector::AbstractVector{T}
end

struct InternalVector{T} <: AbstractVector{T}
    _parent
    _vector::AbstractVector{T}
end
InternalVector(parent, ur::UnitRange) = InternalVector{eltype(ur)}(parent, collect(ur))

function NamedViewVector(nt::NamedTuple)
    _namedtuple = []
    _vector = []
    idx = 1
    for (key, val) in zip(keys(nt), values(nt))
        if val isa Number
            push!(_vector, val)
            push!(_namedtuple, key => idx)
            idx += 1
        elseif val isa AbstractVector
            push!(_vector, val...)
            push!(_namedtuple, key => InternalVector(_vector, idx .+ (0:length(val)-1)))
        # elseif val isa NamedTuple || val isa AbstractDict
        #     val = NamedViewVector(val)
        #     _nt = CoolDict()
        #     for (nvvkey, nvvval) in zip(keys(val), values(val))
        #         len = length(nvvval)
        #         push!(_vector, nvvval...)
        #         _nt[nvvkey] = @view _vector[(end-len+1):end]
        #     end
        #     _namedtuple[key] = NamedTuple(_nt)
        else
            error("$val must be a subtype of of Number or AbstractVector")
        end

    end
    _nt = (; _namedtuple...)
    return NamedViewVector{eltype([_vector...])}(_nt, _vector)
end

ViewableVector = Union{NamedViewVector, InternalVector}

namedtuple(nvv::NamedViewVector) = getfield(nvv, :_namedtuple)

vector(nvv::ViewableVector) = getfield(nvv, :_vector)

Base.size(nvv::ViewableVector) = size(vector(nvv))

Base.length(nvv::ViewableVector) = length(vector(nvv))

Base.iterate(nvv::ViewableVector, args...; kwargs...) =
    iterate(vector(nvv), args...; kwargs...)

Base.getindex(nvv::ViewableVector, idx) = vector(nvv)[idx]
Base.getindex(nvv::NamedViewVector, key::Symbol) = nvv[namedtuple(nvv)[key]]

Base.setindex!(iv::InternalVector, value, idx) = setindex!(iv._parent, value, iv._vector[idx])
Base.setindex!(nvv::NamedViewVector, value, idx) = (vector(nvv)[idx] = value)
Base.setindex!(nvv::NamedViewVector, value, key::Symbol) = (nvv[namedtuple(nvv)[key]] = value)

function Base.getproperty(nvv::NamedViewVector, key::Symbol)
    if key in [:_namedtuple, :_vector]
        error("Hey, the $key field is sorta private. If you want to access it, you can use getfield(myobj, :$key) though. Sorry.")
    else
        return nvv[key]
    end
end
Base.setproperty!(nvv::NamedViewVector, key::Symbol, value) = setindex!(nvv, value, key)

Base.keys(nvv::NamedViewVector) = keys(namedtuple(nvv))

Base.values(nvv::NamedViewVector) = values(namedtuple(nvv))

Base.view(nvv::NamedViewVector, inds...) = view(vector(nvv), namedtuple(nvv)[inds...])

max_stars_repo_name: jonniedie/ComposableDiffEqs

max_issues_repo_name: jonniedie/ComposableDiffEqs

max_forks_repo_name: jonniedie/ComposableDiffEqs


text:
__precompile__()

module EMP

using JuMP
using ReSHOP
using MathOptInterface
MOI = MathOptInterface

using Compat
using Base.Iterators

export @vipair, @NLvipair, EMPmaster, solveEMP, _solveEMP, MathPrgm, EquilibriumProblem, BilevelProblem, addovf!, getsolution, status, termination_status, primal_status, dual_status, get_model_result_num, solve, getobjval

include("JuMPExtension.jl")

" Optimal Value Function (OVF) representation "
mutable struct OVF
    vidx::Int
    args::Vector{Int}
    name::String
    params::Dict{String,Any}
    # No idea how to enforce params to be a Dict{String,Any}
    function OVF(v::JuMP.VariableRef, args::Vector{JuMP.VariableRef}, name::String, params::Dict)
        args_idx = [arg.index.value for arg in args]
        new(v.index.value, args_idx, name, params)
    end
end

" EMP master object "
mutable struct EMPmaster
    backend
    mps::Vector{MathPrgm}
    equils::Vector{Vector{MathPrgm}}
    root::Union{Nothing,MathPrgm,Vector{MathPrgm},Vector{Vector{MathPrgm}}}
    ovfs::Vector{OVF}
    NLoffset::Int
end

"""
    Model([modeling_pkg = "JuMP", solver = "jams"])

Create an EMP master object, using the modeling package given as argument to store the variables and equations.
The solver argument is used to construct the ReSHOPSolver object.
"""
function EMPmaster(; modeling_pkg = "JuMP", solver = "jams")
    if modeling_pkg == "JuMP"
        backend = direct_model(ReSHOP.Optimizer(;solver=solver))
        emp = EMPmaster(backend)
    elseif modeling_pkg == "Convex"
        error("Convex.jl is WIP")
    else
        error("No valid modeling_ds value passed. It should be either ``JuMP'' or ``Convex''")
    end

    return emp
end

"""
    Model(backend)

Create an EMP master object and use the argument as the model data storage object
"""
function EMPmaster(backend)
    if !(backend.moi_backend isa ReSHOP.Optimizer)
        error("The JuMP model must be direct solver, use JuMP.direct_model")
    end
    emp = EMPmaster(backend, Vector{MathPrgm}(), Vector{Vector{MathPrgm}}(), nothing, Vector{OVF}(), -1)

    return emp
end

"""
    MathPrgm(m::EMPmaster)

Create a Mathematical Programm in the EMP master object
"""

" Create an Equilibrium object (not implemented) "
function Equilibrium(emp::EMPmaster, mps::Vector{MathPrgm})
    push!(emp.equils, mps)
end

"""
    EquilibriumProblem(emp::EMPmaster, mps::Vector{MathPrgm})

Define the EMP problem as an equilibrium problem
"""
function EquilibriumProblem(emp::EMPmaster, mps::Vector{MathPrgm})
    emp.root = mps
end

function BilevelProblem(emp::EMPmaster, upper::MathPrgm, lower::MathPrgm)
    emp.root = upper
    push!(upper.mps, lower)
end

include("helpers.jl")

#"""
#    addvar!(mp, var)
#
#Add a variable to a Mathematical Programm
#"""
#function addvar!(mp::MathPrgm, var::JuMP.VariableRef)
#    push!(mp.vars, var.index.value)
#end
#
#function addvar!(mp::MathPrgm, var::Array{JuMP.VariableRef})
#    map(x -> addvar!(mp, x), flatten(var))
#end
#
#function addvar!(mp::MathPrgm, var::JuMP.JuMPArray)
#    map(x-> addvar!(mp, x), values(var))
#end
#
#function addvar!(mp::MathPrgm, var::JuMP.JuMPDict)
#    map(x-> addvar!(mp, x), values(var))
#end

#function addvar!(mp::MathPrgm, var::AbstractArray)
#    map(x -> addvar!(mp, x), flatten(var))
#end

#"""
#    addequ!(mp, eqn)
#
#Add an equation to a Mathematical Programm. The argument eqn is a single or list of JuMP.ConstraintRef, not an expression
#"""
#function addequ!(mp::MathPrgm, eqn::JuMP.ConstraintRef, isNL::Bool)
#    push!(mp.equs, (eqn.index.value, isNL))
#end

#function addequ!(mp::MathPrgm, eqn::Vector{ConstraintRef{M,C}}) where {M <: JuMP.AbstractModel, C <: JuMP.AbstractConstraint}
#    map(x -> addequ!(mp, x), eqn)
#end

#function addequ!(mp::MathPrgm, eqn::Vector{JuMP.ConstraintRef})
#    map(x -> addequ!(mp, x), eqn)
#end

#function addequ!(mp::MathPrgm, eqn::AbstractArray, isNL::Bool)
#    map(x -> addequ!(mp, x, isNL), flatten(eqn))
#
#end
#
#function addequ!(mp::MathPrgm, eqn::Vector, isNL::Bool)
#    map(x-> addequ!(mp, x, isNL), values(eqn))
#end

"""
     addovf!(emp, v, args, name, params)

Tag a variable as a supported OVF

# Arguments
- `emp::EMPmaster`: the EMP object
- `v::JuMP.VariableRef`: the variable
- `args::Vector{JuMP.VariableRef}`: the arguments for the OVF problem
- `name::String`: the name of the OVF
- `params::Dict`: the parameters for defining this OVF
"""
function addovf!(emp::EMPmaster, v::JuMP.VariableRef, args::Vector{JuMP.VariableRef}, name::String, params::Dict)
    push!(emp.ovfs, OVF(v, args, name, params))
end

"""
    solveEMP(emp)

Solve an EMP model
"""
function solveEMP(emp::EMPmaster)
    m = getReSHOPModel(emp)

    if (m.mdl == C_NULL)
        m.mdl = ReSHOP.reshop_alloc(m.ctx)
    end

    ReSHOP.emp_init(m.mdl)

    # This is needed to get the metadata initialized
    if (length(emp.mps) > 0)
        ReSHOP.reshop_set_modeltype(m.ctx, ReSHOP.emp)
    end

    # This is needed to get the NL part loaded, before emptree is called
    if emp.backend.nlp_data !== nothing
        MOI.set(emp.backend, MOI.NLPBlock(), JuMP._create_nlp_block_data(emp.backend))
        empty!(emp.backend.nlp_data.nlconstr_duals)
    end

    emp.NLoffset = m.start_nl_cons

    # Deal with the emptree
    emptree(emp, m.mdl)

    # Deal with the OVF functions
    for ovf in emp.ovfs
        ReSHOP.reshop_ovf(m.mdl, ovf)
    end

    optimize!(emp.backend)
end

#########################################################################################################
# End of public API
#########################################################################################################

include("emptree.jl")

function _solveEMP(emp::EMPmaster)

    # Steps:
    # 1. Define the full models and agents
    # 2. Solve
    # 3. Report values

    reshop_model = getReSHOPModel(emp)

    ReSHOP.make_var_index!(reshop_model)
#    ReSHOP.make_con_index!(reshop_model)

    ########
    #
    # We have to mirror the behavior of ReSHOP, where the nonlinear equation are put before
    # the linear ones.
    #
    # TODO(xhub) fix this mess ...
    #
    #####

    if (length(emp.mps) > 0)
        make_equ_index!(emp.equs, emp.backend.internalModel)
    else
        ReSHOP.make_con_index!(reshop_model)
    end

    ctx = ReSHOP.create_reshop_ctx(reshop_model)
    reshop_model.reshop_ctx = ctx

    reshop_emp = ReSHOP.emp_create(ctx)

    if (length(emp.mps) > 0)
        ReSHOP.reshop_set_modeltype(ctx, ReSHOP.emp)
    else
        ReSHOP.reshop_set_modeltype(reshop_model)
    end

    # Define the EMP tree
    emptree(emp, reshop_emp, ctx)

    for ovf in emp.ovfs
        ReSHOP.reshop_ovf(reshop_emp, ovf)
    end

    setvarnames(ctx, emp.backend)

    gams_ctx, gams_dir = ReSHOP.reshop_setup_gams()
    reshop_model.reshop_ctx_dest = gams_ctx
    reshop_model.gams_dir = gams_dir

    reshop_model.solve_exitcode = ReSHOP.reshop_solve(ctx, gams_ctx, reshop_model.solver_name, reshop_emp)

    ###############################################################################################
    #
    #  Report values and codes
    #
    ###############################################################################################

    if reshop_model.solve_exitcode == 0
        ReSHOP.emp_report_values(reshop_emp, ctx)
        ReSHOP.report_results_common(reshop_model)
        _setvalues(emp.backend, getsolution(emp))
    else
        println("DEBUG: solver failed with status $(reshop_model.solve_exitcode)")
        reshop_model.status = :Error
        reshop_model.solution = fill(NaN,reshop_model.nvar)
        reshop_model.solve_result = "failure"
        reshop_model.solve_result_num = 999
    end

    #    TODO(xhub) move that elsewhere
    # This is left here for now, otherwise, we can't remove the test folder
    # This also plugs a memleak on a context
    ReSHOP.emp_delete(reshop_emp)
end

###################################################################################################
#
#   Accessor and such
#
###################################################################################################

getsolution(emp::EMPmaster) = copy(getReSHOPModel(emp).solution)
getsolvetime(emp::EMPmaster) = getReSHOPModel(emp).solve_time
status(emp::EMPmaster) = status(getReSHOPModel(emp))

function getobjval(mp::MathPrgm)
    reshop_model = getReSHOPModel(mp.emp)
    objvar = ReSHOP.emp_mp_getobjvar(mp.solverobj)
    ReSHOP.ctx_getvarval(reshop_model.reshop_ctx, objvar)
end

# Access to solve results
termination_status(emp::EMPmaster) = JuMP.termination_status(emp.backend)
primal_status(emp::EMPmaster) = JuMP.primal_status(emp.backend)
dual_status(emp::EMPmaster) = JuMP.dual_status(emp.backend)

function Base.show(io::IO, m::EMPmaster)
    println(io, "EMP data structure")
end

end


text:
using Pkg
Pkg.activate("Project.toml")
using Revise # reduce recompilation time
# using Plots
# using Documenter
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using UnPack
using StaticArrays
using DelimitedFiles
using Polyester
using MuladdMacro
using DataFrames
using JLD2
using FileIO

push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D
using Basis2DQuad
using UniformQuadMesh

using SetupDG

# push!(LOAD_PATH, "./examples/EntropyStableEuler.jl/src")
# include("../EntropyStableEuler.jl/src/logmean.jl")
# using EntropyStableEuler
# using EntropyStableEuler.Fluxes2D

@muladd begin

@inline function pfun(rho,rhou,E)
    return (γ-1)*(E-.5*rhou^2/rho)
end

@inline function pfun(rho,rhou,rhov,E)
    return (γ-1)*(E-.5*(rhou^2+rhov^2)/rho)
end

@inline function Efun(rho,u,v,p)
    return p/(γ-1) + .5*rho*(u^2+v^2)
end

@inline function wavespeed_1D(rho,rhou,E)
    p = pfun(rho,rhou,E)
    return abs(rhou/rho) + sqrt(γ*p/rho)
end

@inline function logmean(aL,aR,logL,logR)

    # "from: Entropy stable num. approx. for the isothermal and polytropic Euler"

    da = aR-aL;
    aavg = .5*(aR+aL);
    f = da/aavg;
    v = f^2;
    if abs(f)<1e-4
        # numbers assume the specific value γ = 1.4
        return aavg*(1 + v*(-.2-v*(.0512 - v*0.026038857142857)))
    else
        return -da/(logL-logR)
    end
end

@inline function euler_fluxes_2D(rhoL,uL,vL,betaL,rhologL,betalogL,
                                 rhoR,uR,vR,betaR,rhologR,betalogR)

    rholog  = logmean(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = .5*(rhoL+rhoR)
    uavg   = .5*(uL+uR)
    vavg   = .5*(vL+vR)

    unorm = uL*uR + vL*vR
    pa    = rhoavg/(betaL+betaR)
    f4aux = rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm

    FxS1 = rholog*uavg
    FxS2 = FxS1*uavg + pa
    FxS3 = FxS1*vavg
    FxS4 = f4aux*uavg

    FyS1 = rholog*vavg
    FyS2 = FxS3
    FyS3 = FyS1*vavg + pa
    FyS4 = f4aux*vavg

    return FxS1,FxS2,FxS3,FxS4,FyS1,FyS2,FyS3,FyS4
end

@inline function euler_fluxes_2D_x(rhoL,uL,vL,betaL,rhologL,betalogL,
                                   rhoR,uR,vR,betaR,rhologR,betalogR)

    rholog  = logmean(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = .5*(rhoL+rhoR)
    uavg   = .5*(uL+uR)
    vavg   = .5*(vL+vR)

    unorm = uL*uR + vL*vR
    pa    = rhoavg/(betaL+betaR)
    f4aux = rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm

    FxS1 = rholog*uavg
    FxS2 = FxS1*uavg + pa
    FxS3 = FxS1*vavg
    FxS4 = f4aux*uavg

    return FxS1,FxS2,FxS3,FxS4
end

@inline function euler_fluxes_2D_y(rhoL,uL,vL,betaL,rhologL,betalogL,
                                   rhoR,uR,vR,betaR,rhologR,betalogR)

    rholog  = logmean(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = .5*(rhoL+rhoR)
    uavg   = .5*(uL+uR)
    vavg   = .5*(vL+vR)

    unorm = uL*uR + vL*vR
    pa    = rhoavg/(betaL+betaR)
    f4aux = rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm

    FyS1 = rholog*vavg
    FyS2 = FyS1*uavg
    FyS3 = FyS1*vavg + pa
    FyS4 = f4aux*vavg

    return FyS1,FyS2,FyS3,FyS4
end

@inline function inviscid_flux_prim(rho,u,v,p)
    E = Efun(rho,u,v,p)

    rhou  = rho*u
    rhov  = rho*v
    rhouv = rho*u*v
    Ep    = E+p

    fx1 = rhou
    fx2 = rhou*u+p
    fx3 = rhouv
    fx4 = u*Ep

    fy1 = rhov
    fy2 = rhouv
    fy3 = rhov*v+p
    fy4 = v*Ep

    return fx1,fx2,fx3,fx4,fy1,fy2,fy3,fy4
end

@inline function limiting_param(rhoL,rhouL,rhovL,EL,rhoP,rhouP,rhovP,EP,Lrho,Lrhoe)
    # L - low order, P - P_ij
    l = 1.0
    # Limit density
    if rhoL + rhoP < Lrho
        l = max((Lrho-rhoL)/rhoP, 0.0)
    end

    p = pfun(rhoL+l*rhoP,rhouL+l*rhouP,rhovL+l*rhovP,EL+l*EP)
    if p/(γ-1) > Lrhoe
        return l
    end

    # limiting internal energy (via quadratic function)
    a = rhoP*EP-(rhouP^2+rhovP^2)/2.0
    b = rhoP*EL+rhoL*EP-rhouL*rhouP-rhovL*rhovP-rhoP*Lrhoe
    c = rhoL*EL-(rhouL^2+rhovL^2)/2.0-rhoL*Lrhoe

    d = 1.0/(2.0*a)
    e = b^2-4.0*a*c
    g = -b*d

    l_eps_ij = 1.0
    if e >= 0
        f = sqrt(e)
        h = f*d
        r1 = g+h
        r2 = g-h
        if r1 > TOL && r2 > TOL
            l_eps_ij = min(r1,r2)
        elseif r1 > TOL && r2 < -TOL
            l_eps_ij = r1
        elseif r2 > TOL && r1 < -TOL
            l_eps_ij = r2
        end
    end

    l = min(l,l_eps_ij)
    return l
end

const LIMITOPT = 2 # 1 if elementwise limiting lij, 2 if elementwise limiting li
const POSDETECT = 0 # 1 if turn on detection, 0 otherwise
const LBOUNDTYPE = 0.1 # 0 if use POSTOL as lower bound, if > 0, use LBOUNDTYPE*loworder
const TOL = 5e-16
const POSTOL = 1e-14
const WALLPT = 1.0/6.0
const Nc = 4 # number of components

"Approximation parameters"
const N = parse(Int,ARGS[1])    # N = 2,3,4
const K1D = parse(Int,ARGS[2])  # K = 10,20,40,80,160
const T = 2.0
const dt0 = 1e+1
const XLENGTH = 2.0
const CFL = 0.9
const NUM_THREADS = Threads.nthreads()
const BOTTOMRIGHT = N+1
const TOPRIGHT    = 2*(N+1)
const TOPLEFT     = 3*(N+1)

# Initial condition 2D shocktube
const γ = 1.4

"Mesh related variables"
VX, VY, EToV = uniform_quad_mesh(Int(round(XLENGTH*K1D)),K1D)
@. VX = (VX+1)/2*XLENGTH*5*2
@. VY = (VY+1)/2*5*2

rd = init_reference_quad(N,gauss_lobatto_quad(0,0,N))
"Initialize reference element"
r,_ = gauss_lobatto_quad(0,0,N)   # Reference nodes
VDM = vandermonde_1D(N,r)         # modal to nodal
Dr = grad_vandermonde_1D(N,r)/VDM # nodal differentiation
V1 = vandermonde_1D(1,r)/vandermonde_1D(1,[-1;1]) # nodal linear interpolation

Nq = N
rq,wq = gauss_lobatto_quad(0,0,Nq)
M = diagm(wq)
Mlump = zeros(size(M))
Mlump_inv = zeros(size(M))
for i = 1:Nq+1
    Mlump[i,i] = sum(M[i,:])
    Mlump_inv[i,i] = 1.0/Mlump[i,i]
end

# 1D operators
M1D = Mlump
M1D_inv = Mlump_inv
Q1D = M*Dr
B1D = zeros(N+1,N+1)
B1D[1,1] = -1
B1D[end,end] = 1
L = Array(spdiagm(0=>-2*ones(N+1), 1=>ones(N), -1=>ones(N)))
L[1,1] = -1
L[end,end] = -1
psi = pinv(L)*-1/2*B1D*ones(N+1)
S0 = zeros(N+1,N+1)
for i = 1:N+1
    for j = 1:N+1
        if L[i,j] != 0
            S0[i,j] = psi[j] - psi[i]
        end
    end
end
Q01D = S0+1/2*B1D

# Drop zeros
Q1D = Matrix(droptol!(sparse(Q1D),TOL))
Q01D = Matrix(droptol!(sparse(Q01D),TOL))
B1D = Matrix(droptol!(sparse(B1D),TOL))
S01D = Matrix(droptol!(sparse(S0),TOL))
S1D = Matrix(droptol!(sparse((Q1D-Q1D')/2),TOL))

# Tensor product operators
M = droptol!(sparse(kron(M1D,M1D)),TOL)
Minv = spdiagm(0 => 1 ./ diag(M))
Qr = droptol!(sparse(kron(M1D,Q1D)),TOL)
Qs = droptol!(sparse(kron(Q1D,M1D)),TOL)
Qr0 = droptol!(sparse(kron(M1D,Q01D)),TOL)
Qs0 = droptol!(sparse(kron(Q01D,M1D)),TOL)
Sr = droptol!(sparse((Qr-Qr')/2),TOL)
Ss = droptol!(sparse((Qs-Qs')/2),TOL)
S0r = droptol!(sparse(kron(M1D,S01D)),TOL)
S0s = droptol!(sparse(kron(S01D,M1D)),TOL)
Br = droptol!(sparse(Qr+Qr'),TOL)
Bs = droptol!(sparse(Qs+Qs'),TOL)

@inline function get_infoP(mapP,Fmask,i,k)
    gP = mapP[i,k]           # exterior global face node number
    kP = fld1(gP,Nfp)        # exterior element number
    iP = Fmask[mod1(gP,Nfp)] # exterior node number
    return iP,kP
end

@inline function get_consP(U,f_x,f_y,t,mapP,Fmask,i,iM,xM,yM,uM,vM,k)
    iP,kP = get_infoP(mapP,Fmask,i,k)

    rhoP  = U[1,iP,kP]
    rhouP = U[2,iP,kP]
    rhovP = U[3,iP,kP]
    EP    = U[4,iP,kP]

    return rhoP,rhouP,rhovP,EP
end

@inline function get_fluxP(U,f_x,f_y,t,mapP,Fmask,i,iM,xM,yM,uM,vM,k,rhoP,rhouP,rhovP,EP)
    iP,kP = get_infoP(mapP,Fmask,i,k)

    fx_1_P = f_x[1,iP,kP]
    fx_2_P = f_x[2,iP,kP]
    fx_3_P = f_x[3,iP,kP]
    fx_4_P = f_x[4,iP,kP]
    fy_1_P = f_y[1,iP,kP]
    fy_2_P = f_y[2,iP,kP]
    fy_3_P = f_y[3,iP,kP]
    fy_4_P = f_y[4,iP,kP]

    return fx_1_P,fx_2_P,fx_3_P,fx_4_P,fy_1_P,fy_2_P,fy_3_P,fy_4_P
end

@inline function get_valP(U,f_x,f_y,t,mapP,Fmask,i,iM,xM,yM,uM,vM,k)
    rhoP,rhouP,rhovP,EP = get_consP(U,f_x,f_y,t,mapP,Fmask,i,iM,xM,yM,uM,vM,k)
    fx_1_P,fx_2_P,fx_3_P,fx_4_P,fy_1_P,fy_2_P,fy_3_P,fy_4_P = get_fluxP(U,f_x,f_y,t,mapP,Fmask,i,iM,xM,yM,uM,vM,k,rhoP,rhouP,rhovP,EP)
    return rhoP,rhouP,rhovP,EP,fx_1_P,fx_2_P,fx_3_P,fx_4_P,fy_1_P,fy_2_P,fy_3_P,fy_4_P
end

@inline function is_face_x(i)
    return (((i > BOTTOMRIGHT) & (i <= TOPRIGHT)) | (i > TOPLEFT))
end

@inline function is_face_y(i)
    return (((i > TOPRIGHT) & (i <= TOPLEFT)) | (i <= BOTTOMRIGHT))
end

@inline function update_F_low!(F_low,k,tid,i,j,λ,S0J_ij,U,f)
    # Tensor product elements: f is f_x or f_y
    # S0J_ij - S0xJ_ij or S0yJ_ij
    rho_j  = U[1,j,k]
    rhou_j = U[2,j,k]
    rhov_j = U[3,j,k]
    E_j    = U[4,j,k]
    f_1_j  = f[1,j,k]
    f_2_j  = f[2,j,k]
    f_3_j  = f[3,j,k]
    f_4_j  = f[4,j,k]
    rho_i  = U[1,i,k]
    rhou_i = U[2,i,k]
    rhov_i = U[3,i,k]
    E_i    = U[4,i,k]
    f_1_i  = f[1,i,k]
    f_2_i  = f[2,i,k]
    f_3_i  = f[3,i,k]
    f_4_i  = f[4,i,k]

    FL1 = (S0J_ij*(f_1_i+f_1_j) - λ*(rho_j-rho_i))
    FL2 = (S0J_ij*(f_2_i+f_2_j) - λ*(rhou_j-rhou_i))
    FL3 = (S0J_ij*(f_3_i+f_3_j) - λ*(rhov_j-rhov_i))
    FL4 = (S0J_ij*(f_4_i+f_4_j) - λ*(E_j-E_i))

    F_low[1,i,j,tid] = FL1
    F_low[2,i,j,tid] = FL2
    F_low[3,i,j,tid] = FL3
    F_low[4,i,j,tid] = FL4

    F_low[1,j,i,tid] = -FL1
    F_low[2,j,i,tid] = -FL2
    F_low[3,j,i,tid] = -FL3
    F_low[4,j,i,tid] = -FL4
end

@inline function update_F_high!(F_high,k,tid,i,j,SJ_ij_db,U,rholog,betalog,direction)
    rho_j     = U[1,j,k]
    rhou_j    = U[2,j,k]
    rhov_j    = U[3,j,k]
    E_j       = U[4,j,k]
    p_j       = pfun(rho_j,rhou_j,rhov_j,E_j)
    u_j       = rhou_j/rho_j
    v_j       = rhov_j/rho_j
    beta_j    = rho_j/(2*p_j)
    rholog_j  = rholog[j,k]
    betalog_j = betalog[j,k]
    rho_i     = U[1,i,k]
    rhou_i    = U[2,i,k]
    rhov_i    = U[3,i,k]
    E_i       = U[4,i,k]
    p_i       = pfun(rho_i,rhou_i,rhov_i,E_i)
    u_i       = rhou_i/rho_i
    v_i       = rhov_i/rho_i
    beta_i    = rho_i/(2*p_i)
    rholog_i  = rholog[i,k]
    betalog_i = betalog[i,k]

    if (direction == 0)
        # x direction
        F1,F2,F3,F4 = euler_fluxes_2D_x(rho_i,u_i,v_i,beta_i,rholog_i,betalog_i,
                                        rho_j,u_j,v_j,beta_j,rholog_j,betalog_j)
    else
        # y direction
        F1,F2,F3,F4 = euler_fluxes_2D_y(rho_i,u_i,v_i,beta_i,rholog_i,betalog_i,
                                        rho_j,u_j,v_j,beta_j,rholog_j,betalog_j)
    end
    FH1 = SJ_ij_db*F1
    FH2 = SJ_ij_db*F2
    FH3 = SJ_ij_db*F3
    FH4 = SJ_ij_db*F4

    F_high[1,i,j,tid] = FH1
    F_high[2,i,j,tid] = FH2
    F_high[3,i,j,tid] = FH3
    F_high[4,i,j,tid] = FH4

    F_high[1,j,i,tid] = -FH1
    F_high[2,j,i,tid] = -FH2
    F_high[3,j,i,tid] = -FH3
    F_high[4,j,i,tid] = -FH4
end

function rhs_IDP!(U,rhsU,t,dt,prealloc,ops,geom,in_s1)
    f_x,f_y,rholog,betalog,U_low,U_high,F_low,F_high,F_P,L,wspd_arr,λ_arr,λf_arr,dii_arr,L_plot = prealloc
    S0xJ_vec,S0yJ_vec,S0r_nnzi,S0r_nnzj,S0s_nnzi,S0s_nnzj,SxJ_db_vec,SyJ_db_vec,Sr_nnzi,Sr_nnzj,Ss_nnzi,Ss_nnzj,MJ_inv,BrJ_halved,BsJ_halved,coeff_arr,S_nnzi,S_nnzj = ops
    mapP,Fmask,x,y = geom

    fill!(rhsU,0.0)
    @batch for k = 1:K
        for i = 1:Np
            rho  = U[1,i,k]
            rhou = U[2,i,k]
            rhov = U[3,i,k]
            E    = U[4,i,k]
            p          = pfun(rho,rhou,rhov,E)
            f_x[1,i,k] = rhou
            f_x[2,i,k] = rhou^2/rho+p
            f_x[3,i,k] = rhou*rhov/rho
            f_x[4,i,k] = E*rhou/rho+p*rhou/rho
            f_y[1,i,k] = rhov
            f_y[2,i,k] = rhou*rhov/rho
            f_y[3,i,k] = rhov^2/rho+p
            f_y[4,i,k] = E*rhov/rho+p*rhov/rho
            rholog[i,k]  = log(rho)
            betalog[i,k] = log(rho/(2*p))
        end
    end

    # Precompute wavespeeds
    @batch for k = 1:K
        tid = Threads.threadid()

        # Interior wavespd, leading 2 - x and y directions
        for i = 1:Np
            rho_i  = U[1,i,k]
            rhou_i = U[2,i,k]
            rhov_i = U[3,i,k]
            E_i    = U[4,i,k]
            wspd_arr[i,1,k] = wavespeed_1D(rho_i,rhou_i,E_i)
            wspd_arr[i,2,k] = wavespeed_1D(rho_i,rhov_i,E_i)
        end

        # Interior dissipation coeff
        for c_r = 1:S0r_nnz_hv
            i = S0r_nnzi[c_r]
            j = S0r_nnzj[c_r]
            λ = abs(S0xJ_vec[c_r])*max(wspd_arr[i,1,k],wspd_arr[j,1,k])
            λ_arr[c_r,1,k] = λ
            if in_s1
                dii_arr[i,k] = dii_arr[i,k] + λ
            end
        end

        for c_s = 1:S0s_nnz_hv
            i = S0s_nnzi[c_s]
            j = S0s_nnzj[c_s]
            λ = abs(S0yJ_vec[c_s])*max(wspd_arr[i,2,k],wspd_arr[j,2,k])
            λ_arr[c_s,2,k] = λ
            if in_s1
                dii_arr[i,k] = dii_arr[i,k] + λ
            end
        end
    end

    # Interface dissipation coeff
    @batch for k = 1:K
        for i = 1:Nfp
            iM = Fmask[i]
            BrJ_ii_halved_abs = abs(BrJ_halved[iM])
            BsJ_ii_halved_abs = abs(BsJ_halved[iM])
            xM    = x[iM,k]
            yM    = y[iM,k]

            iP,kP = get_infoP(mapP,Fmask,i,k)

            if is_face_x(i)
                λM = wspd_arr[iM,1,k]
                λP = wspd_arr[iP,1,kP]
                λf = max(λM,λP)*BrJ_ii_halved_abs
                λf_arr[i,k] = λf
                if in_s1
                    dii_arr[iM,k] = dii_arr[iM,k] + λf
                end
            end

            if is_face_y(i)
                λM = wspd_arr[iM,2,k]
                λP = wspd_arr[iP,2,kP]
                λf = max(λM,λP)*BsJ_ii_halved_abs
                λf_arr[i,k] = λf
                if in_s1
                    dii_arr[iM,k] = dii_arr[iM,k] + λf
                end
            end
        end
    end

    # If at the first stage, calculate the time step
    if in_s1
        for k = 1:K
            for i = 1:Np
                dt = min(dt,1.0/MJ_inv[i]/2.0/dii_arr[i,k])
            end
        end
        dt = min(CFL*dt,T-t)
    end

    # =====================
    # Loop through elements
    # =====================
    @batch for k = 1:K
        tid = Threads.threadid()
        for i = 1:Np
            for c = 1:Nc
                U_low[c,i,tid]  = 0.0
                U_high[c,i,tid] = 0.0
            end
        end

        # Calculate low order algebraic flux
        for c_r = 1:S0r_nnz_hv
            i = S0r_nnzi[c_r]
            j = S0r_nnzj[c_r]
            λ = λ_arr[c_r,1,k]
            S0xJ_ij = S0xJ_vec[c_r]
            update_F_low!(F_low,k,tid,i,j,λ,S0xJ_ij,U,f_x)
        end

        for c_s = 1:S0s_nnz_hv
            i = S0s_nnzi[c_s]
            j = S0s_nnzj[c_s]
            λ = λ_arr[c_s,2,k]
            S0yJ_ij = S0yJ_vec[c_s]
            update_F_low!(F_low,k,tid,i,j,λ,S0yJ_ij,U,f_y)
        end

        # Calculate high order algebraic flux
        for c_r = 1:Sr_nnz_hv
            i         = Sr_nnzi[c_r]
            j         = Sr_nnzj[c_r]
            SxJ_ij_db = SxJ_db_vec[c_r]
            update_F_high!(F_high,k,tid,i,j,SxJ_ij_db,U,rholog,betalog,0)
        end

        for c_r = 1:Sr_nnz_hv
            i         = Ss_nnzi[c_r]
            j         = Ss_nnzj[c_r]
            SyJ_ij_db = SyJ_db_vec[c_r]
            update_F_high!(F_high,k,tid,i,j,SyJ_ij_db,U,rholog,betalog,1)
        end

        # Calculate interface fluxes
        for i = 1:Nfp
            iM    = Fmask[i]
            BrJ_ii_halved = BrJ_halved[iM]
            BsJ_ii_halved = BsJ_halved[iM]
            xM    = x[iM,k]
            yM    = y[iM,k]
            rhoM  = U[1,iM,k]
            rhouM = U[2,iM,k]
            rhovM = U[3,iM,k]
            EM    = U[4,iM,k]
            uM    = rhouM/rhoM
            vM    = rhovM/rhoM
            fx_1_M = f_x[1,iM,k]
            fx_2_M = f_x[2,iM,k]
            fx_3_M = f_x[3,iM,k]
            fx_4_M = f_x[4,iM,k]
            fy_1_M = f_y[1,iM,k]
            fy_2_M = f_y[2,iM,k]
            fy_3_M = f_y[3,iM,k]
            fy_4_M = f_y[4,iM,k]

            rhoP,rhouP,rhovP,EP,fx_1_P,fx_2_P,fx_3_P,fx_4_P,fy_1_P,fy_2_P,fy_3_P,fy_4_P = get_valP(U,f_x,f_y,t,mapP,Fmask,i,iM,xM,yM,uM,vM,k)
            λ = λf_arr[i,k]

            # flux in x direction
            if is_face_x(i)
                F_P[1,i,tid] = (BrJ_ii_halved*(fx_1_M+fx_1_P)
                               -λ*(rhoP-rhoM) )
                F_P[2,i,tid] = (BrJ_ii_halved*(fx_2_M+fx_2_P)
                               -λ*(rhouP-rhouM) )
                F_P[3,i,tid] = (BrJ_ii_halved*(fx_3_M+fx_3_P)
                               -λ*(rhovP-rhovM) )
                F_P[4,i,tid] = (BrJ_ii_halved*(fx_4_M+fx_4_P)
                               -λ*(EP-EM) )
            end

            # flux in y direction
            if is_face_y(i)
                F_P[1,i,tid] = (BsJ_ii_halved*(fy_1_M+fy_1_P)
                               -λ*(rhoP-rhoM) )
                F_P[2,i,tid] = (BsJ_ii_halved*(fy_2_M+fy_2_P)
                               -λ*(rhouP-rhouM) )
                F_P[3,i,tid] = (BsJ_ii_halved*(fy_3_M+fy_3_P)
                               -λ*(rhovP-rhovM) )
                F_P[4,i,tid] = (BsJ_ii_halved*(fy_4_M+fy_4_P)
                               -λ*(EP-EM) )
            end
        end

        # Calculate low order solution
        for j = 1:Np
            for i = 1:Np
                for c = 1:Nc
                    U_low[c,i,tid]  = U_low[c,i,tid]  + F_low[c,i,j,tid]
                    U_high[c,i,tid] = U_high[c,i,tid] + F_high[c,i,j,tid]
                end
            end
        end

        for i = 1:Nfp
            iM = Fmask[i]
            for c = 1:Nc
                U_low[c,iM,tid]  = U_low[c,iM,tid]  + F_P[c,i,tid]
                U_high[c,iM,tid] = U_high[c,iM,tid] + F_P[c,i,tid]
            end
        end

        for i = 1:Np
            for c = 1:Nc
                U_low[c,i,tid]  = U[c,i,k] - dt*MJ_inv[i]*U_low[c,i,tid]
                U_high[c,i,tid] = U[c,i,k] - dt*MJ_inv[i]*U_high[c,i,tid]
            end
        end

        is_H_positive = true
        for i = 1:Np
            rhoH_i  = U_high[1,i,tid]
            rhouH_i = U_high[2,i,tid]
            rhovH_i = U_high[3,i,tid]
            EH_i    = U_high[4,i,tid]
            pH_i    = pfun(rhoH_i,rhouH_i,rhovH_i,EH_i)
            if pH_i < POSTOL || rhoH_i < POSTOL
                is_H_positive = false
            end
        end

        if POSDETECT == 0
            is_H_positive = false
        end

        if (LIMITOPT == 1)
            # Calculate limiting parameters
            if !is_H_positive
                for ni = 1:S_nnz_hv
                    i = S_nnzi[ni]
                    j = S_nnzj[ni]
                    coeff_i = dt*coeff_arr[i]
                    coeff_j = dt*coeff_arr[j]
                    rhoi  = U_low[1,i,tid]
                    rhoui = U_low[2,i,tid]
                    rhovi = U_low[3,i,tid]
                    Ei    = U_low[4,i,tid]
                    rhoj  = U_low[1,j,tid]
                    rhouj = U_low[2,j,tid]
                    rhovj = U_low[3,j,tid]
                    Ej    = U_low[4,j,tid]
                    rhoP  = (F_low[1,i,j,tid]-F_high[1,i,j,tid])
                    rhouP = (F_low[2,i,j,tid]-F_high[2,i,j,tid])
                    rhovP = (F_low[3,i,j,tid]-F_high[3,i,j,tid])
                    EP    = (F_low[4,i,j,tid]-F_high[4,i,j,tid])
                    rhoP_i  = coeff_i*rhoP
                    rhouP_i = coeff_i*rhouP
                    rhovP_i = coeff_i*rhovP
                    EP_i    = coeff_i*EP
                    rhoP_j  = -coeff_j*rhoP
                    rhouP_j = -coeff_j*rhouP
                    rhovP_j = -coeff_j*rhovP
                    EP_j    = -coeff_j*EP
                    Lrho      = POSTOL
                    Lrhoe     = POSTOL
                    L[ni,tid] = min(limiting_param(rhoi,rhoui,rhovi,Ei,rhoP_i,rhouP_i,rhovP_i,EP_i,Lrho,Lrhoe),
                                    limiting_param(rhoj,rhouj,rhovj,Ej,rhoP_j,rhouP_j,rhovP_j,EP_j,Lrho,Lrhoe))
                end
            end
        elseif (LIMITOPT == 2)
            # li limiting
            li_min = 1.0 # limiting parameter for element
            for i = 1:Np
                rhoi  = U_low[1,i,tid]
                rhoui = U_low[2,i,tid]
                rhovi = U_low[3,i,tid]
                Ei    = U_low[4,i,tid]
                rhoP_i  = 0.0
                rhouP_i = 0.0
                rhovP_i = 0.0
                EP_i    = 0.0
                coeff   = dt*MJ_inv[i]
                for j = 1:Np
                    rhoP   = (F_low[1,i,j,tid]-F_high[1,i,j,tid])
                    rhouP  = (F_low[2,i,j,tid]-F_high[2,i,j,tid])
                    rhovP  = (F_low[3,i,j,tid]-F_high[3,i,j,tid])
                    EP     = (F_low[4,i,j,tid]-F_high[4,i,j,tid])
                    rhoP_i  = rhoP_i  + coeff*rhoP
                    rhouP_i = rhouP_i + coeff*rhouP
                    rhovP_i = rhovP_i + coeff*rhovP
                    EP_i    = EP_i    + coeff*EP
                end
                if (LBOUNDTYPE == 0)
                    Lrho  = POSTOL
                    Lrhoe = POSTOL
                elseif (LBOUNDTYPE > 0.0)
                    Lrho  = LBOUNDTYPE*rhoi
                    Lrhoe = LBOUNDTYPE*pfun(rhoi,rhoui,rhovi,Ei)/(γ-1)
                end
                l = limiting_param(rhoi,rhoui,rhovi,Ei,rhoP_i,rhouP_i,rhovP_i,EP_i,Lrho,Lrhoe)
                li_min = min(li_min,l)
            end

            for ni = 1:S_nnz_hv
                L[ni,tid] = li_min
            end
        end

        # Elementwise limiting
        l_e = 1.0
        for i = 1:S_nnz_hv
            l_e = min(l_e,L[i,tid])
        end
        l_em1 = l_e-1.0

        if is_H_positive
            l_e = 1.0
            l_em1 = 0.0
        end

        if in_s1
            L_plot[k] = l_e
        end

        for ni = 1:S_nnz_hv
            i     = S_nnzi[ni]
            j     = S_nnzj[ni]
            for c = 1:Nc
                FL_ij = F_low[c,i,j,tid]
                FH_ij = F_high[c,i,j,tid]
                rhsU[c,i,k] = rhsU[c,i,k] + l_em1*FL_ij - l_e*FH_ij
                rhsU[c,j,k] = rhsU[c,j,k] - l_em1*FL_ij + l_e*FH_ij
                # l_e   = L[ni,tid]
                # l_em1 = L[ni,tid]-1.0
                # rhsU[c,i,k] = rhsU[c,i,k] + l_em1*FL_ij - l_e*FH_ij
                # rhsU[c,j,k] = rhsU[c,j,k] - l_em1*FL_ij + l_e*FH_ij
            end
        end

        for i = 1:Nfp
            iM = Fmask[i]
            for c = 1:Nc
                rhsU[c,iM,k] = rhsU[c,iM,k] - F_P[c,i,tid]
            end
        end
    end

    @batch for k = 1:K
        for i = 1:Np
            for c = 1:Nc
                rhsU[c,i,k] = MJ_inv[i]*rhsU[c,i,k]
            end
        end
    end

    @show sum(L_plot)/K

    return dt
end

@unpack Vf,Dr,Ds,LIFT = rd
md = init_mesh((VX,VY),EToV,rd)
@unpack xf,yf,mapM,mapP,mapB,nxJ,nyJ,x,y = md
xb,yb = (x->x[mapB]).((xf,yf))

const K  = size(x,2)
const Nfaces = 4

# Make domain periodic
@unpack Vf = rd
@unpack xf,yf,mapM,mapP,mapB = md
LX,LY = (x->maximum(x)-minimum(x)).((VX,VY)) # find lengths of domain
mapPB = build_periodic_boundary_maps(xf,yf,LX,LY,Nfaces*K,mapM,mapP,mapB)
mapP[mapB] = mapPB
@pack! md = mapP

const Np = (N+1)*(N+1)
const Nfp    = Nfaces*(N+1)
const J   = 25*1.0/K1D/K1D/4*4
const Jf  = 5*1.0/K1D/2*2
const rxJ = 2*K1D*J/5/2
const syJ = 2*K1D*J/5/2
const sJ  = Jf
const rxJ_sq = rxJ^2
const syJ_sq = syJ^2
const rxJ_db = 2*rxJ
const syJ_db = 2*syJ

# Convert diagonal matrices to vectors
Br   = Array(diag(Br))
Bs   = Array(diag(Bs))
M    = Array(diag(M))
Minv = Array(diag(Minv))

MJ_inv    = Minv./J
Br_halved = -sum(S0r,dims=2)
Bs_halved = -sum(S0s,dims=2)
coeff_arr = 2*N*MJ_inv

Fmask  = [1:N+1; (N+1):(N+1):Np; Np:-1:Np-N; Np-N:-(N+1):1]
Fxmask = [(N+2):(2*N+2); (3*N+4):(4*N+4)]
Fymask = [1:(N+1); (2*N+3):(3*N+3)]

S0r_nnz = length(nonzeros(S0r))
S0s_nnz = length(nonzeros(S0s))
Sr_nnz  = length(nonzeros(Sr))
Ss_nnz  = length(nonzeros(Ss))
const S0r_nnz_hv = div(S0r_nnz,2)
const S0s_nnz_hv = div(S0s_nnz,2)
const S0_nnz_hv  = S0r_nnz_hv+S0s_nnz_hv
const Sr_nnz_hv  = div(Sr_nnz,2)
const Ss_nnz_hv  = div(Ss_nnz,2)
const S_nnz_hv   = Sr_nnz_hv+Ss_nnz_hv
S0r_vec  = zeros(S0r_nnz_hv)
S0s_vec  = zeros(S0s_nnz_hv)
S0r_nnzi = zeros(Int32,S0r_nnz_hv)
S0s_nnzi = zeros(Int32,S0s_nnz_hv)
S0r_nnzj = zeros(Int32,S0r_nnz_hv)
S0s_nnzj = zeros(Int32,S0s_nnz_hv)
Sr_vec   = zeros(Sr_nnz_hv)
Ss_vec   = zeros(Ss_nnz_hv)
Sr_nnzi  = zeros(Int32,Sr_nnz_hv)
Ss_nnzi  = zeros(Int32,Ss_nnz_hv)
Sr_nnzj  = zeros(Int32,Sr_nnz_hv)
Ss_nnzj  = zeros(Int32,Ss_nnz_hv)
S_nnzi   = zeros(Int32,S_nnz_hv)
S_nnzj   = zeros(Int32,S_nnz_hv)
global count_r0 = 1
global count_s0 = 1
global count_r  = 1
global count_s  = 1
global count    = 1
for j = 2:Np
    for i = 1:j-1
        S0r_ij = S0r[i,j]
        S0s_ij = S0s[i,j]
        Sr_ij  = Sr[i,j]
        Ss_ij  = Ss[i,j]
        if S0r_ij != 0
            global S0r_vec[count_r0]  = S0r_ij
            global S0r_nnzi[count_r0] = i
            global S0r_nnzj[count_r0] = j
            global count_r0 = count_r0+1
        end
        if S0s_ij != 0
            global S0s_vec[count_s0]  = S0s_ij
            global S0s_nnzi[count_s0] = i
            global S0s_nnzj[count_s0] = j
            global count_s0 = count_s0+1
        end
        if Sr_ij != 0
            global Sr_vec[count_r]  = Sr_ij
            global Sr_nnzi[count_r] = i
            global Sr_nnzj[count_r] = j
            global S_nnzi[count]    = i
            global S_nnzj[count]    = j
            global count_r = count_r+1
            global count   = count+1
        end
        if Ss_ij != 0
            global Ss_vec[count_s]  = Ss_ij
            global Ss_nnzi[count_s] = i
            global Ss_nnzj[count_s] = j
            global S_nnzi[count]    = i
            global S_nnzj[count]    = j
            global count_s = count_s+1
            global count   = count+1
        end
    end
end

# Scale by Jacobian
S0xJ_vec   = rxJ*S0r_vec
S0yJ_vec   = syJ*S0s_vec
SxJ_db_vec = 2*rxJ*Sr_vec
SyJ_db_vec = 2*syJ*Ss_vec
BrJ_halved = Jf*Br_halved
BsJ_halved = Jf*Bs_halved

# 2D shocktube
const TOP_INIT = (1+sqrt(3)/6)/sqrt(3)

# Initial condition 2D shocktube
# at_left(x,y) = y-sqrt(3)*x+sqrt(3)/6 > 0.0
function vortex_sol(x,y,t)

    x0 = 9.0
    y0 = 5.0
    beta = 8.5
    r2 = @. (x-x0-t)^2 + (y-y0)^2

    u = @. 1 - beta*exp(1-r2)*(y-y0)/(2*pi)
    v = @. beta*exp(1-r2)*(x-x0-t)/(2*pi)
    rho = @. 1 - (1/(8*γ*pi^2))*(γ-1)/2*(beta*exp(1-r2))^2
    rho = @. rho^(1/(γ-1))
    p = @. rho^γ

    return (rho, u, v, p)
end

U = zeros(Nc,Np,K)
for k = 1:K
    for i = 1:Np
        # if x[i,k] < 1.0#at_left(x[i,k],y[i,k])
        #     U[1,i,k] = 1.0
        #     U[2,i,k] = 1.0
        #     U[3,i,k] = 0.0
        #     U[4,i,k] = Efun(1,1,0,1)
        # else
        #     U[1,i,k] = 0.5
        #     U[2,i,k] = 1.0
        #     U[3,i,k] = 0.0
        #     U[4,i,k] = Efun(0.5,1,0,1)
        # end
        rho,u,v,p = vortex_sol(x[i,k],y[i,k],0.0)
        U[1,i,k] = rho
        U[2,i,k] = rho*u
        U[3,i,k] = rho*v
        U[4,i,k] = Efun(rho,u,v,p)
    end
end

# Preallocation
rhsU   = zeros(Float64,size(U))
f_x    = zeros(Float64,size(U))
f_y    = zeros(Float64,size(U))
rholog  = zeros(Float64,Np,K)
betalog = zeros(Float64,Np,K)
U_low   = zeros(Float64,Nc,Np,NUM_THREADS)
U_high  = zeros(Float64,Nc,Np,NUM_THREADS)
F_low   = zeros(Float64,Nc,Np,Np,NUM_THREADS)
F_high  = zeros(Float64,Nc,Np,Np,NUM_THREADS)
F_P     = zeros(Float64,Nc,Nfp,NUM_THREADS)
L       =  ones(Float64,S_nnz_hv,NUM_THREADS)
wspd_arr = zeros(Float64,Np,2,K)
λ_arr    = zeros(Float64,S0r_nnz_hv,2,K) # Assume S0r and S0s has same number of nonzero entries
λf_arr   = zeros(Float64,Nfp,K)
dii_arr  = zeros(Float64,Np,K)
L_plot   = zeros(Float64,K)

prealloc = (f_x,f_y,rholog,betalog,U_low,U_high,F_low,F_high,F_P,L,wspd_arr,λ_arr,λf_arr,dii_arr,L_plot)
ops      = (S0xJ_vec,  S0yJ_vec,  S0r_nnzi,S0r_nnzj,S0s_nnzi,S0s_nnzj,
            SxJ_db_vec,SyJ_db_vec,Sr_nnzi, Sr_nnzj, Ss_nnzi, Ss_nnzj,
            MJ_inv,BrJ_halved,BsJ_halved,coeff_arr,S_nnzi,S_nnzj)
geom     = (mapP,Fmask,x,y)

# Time stepping
"Time integration"
t = 0.0
U = collect(U)
resW = zeros(size(U))

# #plotting nodes
# @unpack VDM = rd
# rp,sp = equi_nodes_2D(10)
# Vp = vandermonde_2D(N,rp,sp)/VDM
# gr(aspect_ratio=:equal,legend=false,
#    markerstrokewidth=0,markersize=2)

#=
#Timing
dt = dt0
#rhs_IDP!(U,rhsU,t,dt,prealloc,ops,geom,true);
@btime rhs_IDP!($U,$rhsU,$t,$dt,$prealloc,$ops,$geom,true);
# @profiler rhs_IDP!(U,rhsU,t,dt,prealloc,ops,geom,true);
=#

dt_hist = []
i = 1

# open("/data/yl184/N=$N,K1D=$K1D,t=$t,CFL=$CFL,x,convquad.txt","w") do io
#     writedlm(io,x)
# end
# open("/data/yl184/N=$N,K1D=$K1D,t=$t,CFL=$CFL,y,convquad.txt","w") do io
#     writedlm(io,y)
# end

@time while t < T
#while i < 2
    # SSPRK(3,3)
    fill!(dii_arr,0.0)
    # dt = min(dt0,T-t)
    dt = dt0
    dt = rhs_IDP!(U,rhsU,t,dt,prealloc,ops,geom,true);
    @. resW = U + dt*rhsU
    rhs_IDP!(resW,rhsU,t+dt,dt,prealloc,ops,geom,false);
    @. resW = resW+dt*rhsU
    @. resW = 3/4*U+1/4*resW
    rhs_IDP!(resW,rhsU,t+dt/2,dt,prealloc,ops,geom,false);
    @. resW = resW+dt*rhsU
    @. U = 1/3*U+2/3*resW

    push!(dt_hist,dt)
    global t = t + dt
    println("Current time $t with time step size $dt, and final time $T, at step $i")
    flush(stdout)
    global i = i + 1
    # if ((mod(i,100) == 1) | (i >= 55))
    if (mod(i,100) == 1)
        # open("/data/yl184/N=$N,K1D=$K1D,t=$t,CFL=$CFL,rho,convquad.txt","w") do io
        #     writedlm(io,U[1,:,:])
        # end
        # open("/data/yl184/N=$N,K1D=$K1D,t=$t,CFL=$CFL,rhou,convquad.txt","w") do io
        #     writedlm(io,U[2,:,:])
        # end
        # open("/data/yl184/N=$N,K1D=$K1D,t=$t,CFL=$CFL,rhov,convquad.txt","w") do io
        #     writedlm(io,U[3,:,:])
        # end
        # open("/data/yl184/N=$N,K1D=$K1D,t=$t,CFL=$CFL,E,convquad.txt","w") do io
        #     writedlm(io,U[4,:,:])
        # end
        # open("/data/yl184/N=$N,K1D=$K1D,t=$t,CFL=$CFL,L_plot,convquad.txt","w") do io
        #     writedlm(io,L_plot)
        # end

    end
end

# open("/data/yl184/N=$N,K1D=$K1D,t=$t,CFL=$CFL,rho,convquad.txt","w") do io
#     writedlm(io,U[1,:,:])
# end
# open("/data/yl184/N=$N,K1D=$K1D,t=$t,CFL=$CFL,rhou,convquad.txt","w") do io
#     writedlm(io,U[2,:,:])
# end
# open("/data/yl184/N=$N,K1D=$K1D,t=$t,CFL=$CFL,rhov,convquad.txt","w") do io
#     writedlm(io,U[3,:,:])
# end
# open("/data/yl184/N=$N,K1D=$K1D,t=$t,CFL=$CFL,E,convquad.txt","w") do io
#     writedlm(io,U[4,:,:])
# end
# open("/data/yl184/N=$N,K1D=$K1D,t=$t,CFL=$CFL,L_plot,convquad.txt","w") do io
#     writedlm(io,L_plot)
# end

exact_U = @. vortex_sol.(x,y,T)
exact_rho = [x[1] for x in exact_U]
exact_u = [x[2] for x in exact_U]
exact_v = [x[3] for x in exact_U]
exact_p = [x[4] for x in exact_U]
exact_rhou = exact_rho .* exact_u
exact_rhov = exact_rho .* exact_v
exact_E = Efun.(exact_rho,exact_u,exact_v,exact_p)

rho = U[1,:,:]
rhou = U[2,:,:]
rhov = U[3,:,:]
u = U[2,:,:]./U[1,:,:]
v = U[3,:,:]./U[1,:,:]
E = U[4,:,:]

Linferr = maximum(abs.(exact_rho-rho))/maximum(abs.(exact_rho)) +
          maximum(abs.(exact_rhou-rhou))/maximum(abs.(exact_rhou)) +
          maximum(abs.(exact_rhov-rhov))/maximum(abs.(exact_rhov)) +
          maximum(abs.(exact_E-E))/maximum(abs.(exact_E))

L1err = sum(J*M.*abs.(exact_rho-rho))/sum(J*M.*abs.(exact_rho)) +
        sum(J*M.*abs.(exact_rhou-rhou))/sum(J*M.*abs.(exact_rhou)) +
        sum(J*M.*abs.(exact_rhov-rhov))/sum(J*M.*abs.(exact_rhov)) +
        sum(J*M.*abs.(exact_E-E))/sum(J*M.*abs.(exact_E))

L2err = sqrt(sum(J*M.*abs.(exact_rho-rho).^2))/sqrt(sum(J*M.*abs.(exact_rho).^2)) +
        sqrt(sum(J*M.*abs.(exact_rhou-rhou).^2))/sqrt(sum(J*M.*abs.(exact_rhou).^2)) +
        sqrt(sum(J*M.*abs.(exact_rhov-rhov).^2))/sqrt(sum(J*M.*abs.(exact_rhov).^2)) +
        sqrt(sum(J*M.*abs.(exact_E-E).^2))/sqrt(sum(J*M.*abs.(exact_E).^2))

println("N = $N, K = $K")
println("L1 error is $L1err")
println("L2 error is $L2err")
println("Linf error is $Linferr")

# df = DataFrame(N = Int64[], K = Int64[], T = Float64[], CFL = Float64[], LIMITOPT = Int64[], POSDETECT = Int64[], LBOUNDTYPE = Float64[], L1err = Float64[], L2err = Float64[], Linferr = Float64[])
df = load("dg2D_euler_quad_convergence.jld2","convergence_data")
push!(df,(N,K,T,CFL,LIMITOPT,POSDETECT,LBOUNDTYPE,L1err,L2err,Linferr))
save("dg2D_euler_quad_convergence.jld2","convergence_data",df)

end #muladd

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
using Pkg
# Pkg.develop(path=joinpath(@__DIR__, "../../DojoEnvironments"))
# Pkg.develop(path=joinpath(@__DIR__, "../.."))
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

# ## visualizer
vis = Visualizer()
open(vis)

# ## setup
using Dojo
using IterativeLQR
using LinearAlgebra
using FiniteDiff
using DojoEnvironments

################################################################################
# Continuation
################################################################################
function reset!(env::Environment; rtol=1e-4, btol=1e-3, undercut=5.0)
    env.opts_step.rtol = rtol
    env.opts_step.btol = btol
    env.opts_step.undercut = undercut
    env.opts_grad.rtol = rtol
    env.opts_grad.btol = btol
    env.opts_grad.undercut = undercut
    return nothing
end

function continuation_callback!(solver::Solver, env::Environment; ρ=1.5)
    # contact smoothness continuation
    env.opts_step.rtol = max(1e-6, env.opts_step.rtol/ρ)
    env.opts_step.btol = max(1e-4, env.opts_step.btol/ρ)
    env.opts_grad.rtol = max(1e-6, env.opts_grad.rtol/ρ)
    env.opts_grad.btol = max(1e-4, env.opts_grad.btol/ρ)

    # visualize current policy
    ū = solver.problem.actions
    x̄ = IterativeLQR.rollout(model, x1, ū)
    DojoEnvironments.visualize(env, x̄)

    println("r_tol $(scn(env.opts_grad.rtol))  " *
        "κ_tol $(scn(env.opts_grad.btol))")
    return nothing
end

################################################################################
# ## system
################################################################################
gravity = -9.81
timestep = 0.02
friction_coefficient = 0.8
damper = 0.5
spring = 5.0
env = get_environment(:quadruped,
    representation=:minimal,
    timestep=timestep,
    contact_body=false,
    gravity=gravity,
    friction_coefficient=friction_coefficient,
    damper=damper,
    spring=spring,
    infeasible_control=true,
    vis=vis)

# ## dimensions
n = env.num_states
m = env.num_inputs
nu_infeasible = 6

# ## template
include(joinpath(@__DIR__, "../../DojoEnvironments/src",
    "quadruped/methods/template.jl"))

################################################################################
# ## simulation test
################################################################################
mech = get_mechanism(:quadruped,
    contact_body=false,
    timestep=timestep,
    gravity=gravity,
    friction_coefficient=friction_coefficient,
    damper=damper,
    spring=spring)

initialize!(mech, :quadruped, body_position=[0,0,0.0])
# u_hover = [0.05;0;1.13; 0;-0.01;0; zeros(12)]
u_hover = [0.02;0;0.6; 0;0;0; zeros(12)]
function ctrl!(m, k; u=u_hover)
    nu = input_dimension(m)
    set_input!(m, SVector{nu}(u))
end

storage = simulate!(mech, 1.0, ctrl!,
    record=true,
    verbose=true,
    opts=SolverOptions(rtol=1e-5, btol=1e-4, undercut=5.0, verbose=false),
    )
Dojo.visualize(mech, storage, vis=env.vis)

################################################################################
# ## reference trajectory
################################################################################
N = 1
initialize!(env.mechanism, :quadruped)
xref = quadruped_trajectory(env.mechanism,
    r=0.08,
    z=0.29;
    Δx=-0.04,
    Δfront=0.10,
    width_scale=0.0,
    height_scale=1.0,
    N=18,
    Ncycles=N)
zref = [minimal_to_maximal(env.mechanism, x) for x in xref]
DojoEnvironments.visualize(env, xref)

# ## horizon
T = length(zref)

################################################################################
# ## ILQR problem
################################################################################
# ## model
dyn = IterativeLQR.Dynamics(
    (y, x, u, w) -> dynamics(y, env, x, u, w),
    (dx, x, u, w) -> dynamics_jacobian_state(dx, env, x, u, w),
    (du, x, u, w) -> dynamics_jacobian_input(du, env, x, u, w),
    n, n, m)

model = [dyn for t = 1:T-1]

# ## rollout
x1 = deepcopy(xref[1])
ū = [u_hover for t = 1:T-1]

x̄ = IterativeLQR.rollout(model, x1, ū)
DojoEnvironments.visualize(env, x̄)

# ## objective
############################################################################
qt = [0.3; 0.05; 0.05;
    5e-2 * ones(3);
    1e-3 * ones(3);
    1e-3 * ones(3);
    fill([2, 1e-3], 12)...]
ots = [(x, u, w) -> transpose(x - xref[t]) * Diagonal(timestep * qt) * (x - xref[t]) +
# transpose(u - u_hover) * Diagonal(timestep * 0.01 * ones(m)) * (u - u_hover) for t = 1:T-1]
    transpose(u) * Diagonal(timestep * 0.01 * ones(m)) * u for t = 1:T-1]
oT = (x, u, w) -> transpose(x - xref[end]) * Diagonal(timestep * qt) * (x - xref[end])

cts = [IterativeLQR.Cost(ot, n, m) for ot in ots]
cT = IterativeLQR.Cost(oT, n, 0)
obj = [cts..., cT]

# ## constraints
############################################################################
ul = -1.0 * 1e-3*ones(nu_infeasible)
uu = +1.0 * 1e-3*ones(nu_infeasible)

function contt(x, u, w)
    [
        1e-1 * (ul - u[1:nu_infeasible]);
        1e-1 * (u[1:nu_infeasible] - uu);
    ]
end

function goal(x, u, w)
    Δ = 1e-2 * (x - xref[end])[[1:6;13:2:36]]
    return Δ
end

con_policyt = IterativeLQR.Constraint(contt, n, m, indices_inequality=collect(1:2nu_infeasible))
con_policyT = IterativeLQR.Constraint(goal, n, 0)

cons = [[con_policyt for t = 1:T-1]..., con_policyT]

# ## solver
options = Options(line_search=:armijo,
        max_iterations=50,
        max_dual_updates=12,
        min_step_size=1e-2,
        objective_tolerance=1e-3,
        lagrangian_gradient_tolerance=1e-3,
        constraint_tolerance=1e-4,
        initial_constraint_penalty=1e-1,
        scaling_penalty=10.0,
        max_penalty=1e4,
        verbose=true)

s = IterativeLQR.Solver(model, obj, cons, options=options)

IterativeLQR.initialize_controls!(s, ū)
IterativeLQR.initialize_states!(s, x̄)

# ## solve
local_callback!(solver::IterativeLQR.Solver) = continuation_callback!(solver, env)
reset!(env)
@time IterativeLQR.constrained_ilqr_solve!(s, augmented_lagrangian_callback! = local_callback!)

# ## solution
x_sol, u_sol = IterativeLQR.get_trajectory(s)

# ## visualize
x_view = [[x_sol[1] for t = 1:15]..., x_sol..., [x_sol[end] for t = 1:15]...]
DojoEnvironments.visualize(env, x_view)

# x1 = deepcopy(xref[1])
# ū = vcat(fill(u_sol, 5)...)
# x̄ = IterativeLQR.rollout(fill(dyn, length(ū)), x1, ū)
# DojoEnvironments.visualize(env, x̄)

mech = get_mechanism(:quadruped,
    timestep=timestep,
    contact_body=false,
    gravity=gravity,
    friction_coefficient=friction_coefficient,
    damper=damper,
    spring=spring)

z_sim = get_maximal_state(storage)
x_sim = [maximal_to_minimal(mech, z) for z in z_sim]
H = length(z_sim)
u_sim = [zeros(m) for i=1:H-1]

for i = 1:H-1
    @show i
    set_minimal_state!(mech, x_sim[i])
    set_input!(mech, u_sim[i])
    z = minimal_to_maximal(mech, x_sim[i])
    u = u_sim[i]
    J = minimal_to_maximal_jacobian(mech, x_sim[i])
    # step!(mech, z, u, opts=SolverOptions(rtol=1e-6, btol=2e-4))
end

full_vector(mech.system)

z = z_sim[1]
u = u_sim[1]
step!(mech, z, u, opts=SolverOptions(rtol=1e-8, btol=1e-6))
z20 = get_maximal_state(mech)
z30 = get_next_state(mech)
x20 = maximal_to_minimal(mech, z20)
x30 = maximal_to_minimal(mech, z30)

function get_initial_configurations(mechanism::Mechanism, z::Vector{T}) where T
    configuration_indices = [1,2,3,4,5,6,13,15,17,19,21,23,25,27,29,31,33,35]
    set_maximal_state!(mechanism, z)

    # current configuration
    x = maximal_to_minimal(mechanism, z)
    q2 = x[configuration_indices]

    # previous configuration
    for body in mechanism.bodies
        x, q = previous_configuration(body.state)
        set_maximal_configurations!(body, x=x, q=q)
    end
    z = get_maximal_state(mechanism)
    x = maximal_to_minimal(mechanism, z)
    q1 = x[configuration_indices]
    return q1, q2
end

configuration_indices = [1,2,3,4,5,6,13,15,17,19,21,23,25,27,29,31,33,35]
q10, q20 = get_initial_configurations(mech, z20)
q30 = x30[configuration_indices]

u0 = zeros(m)
w0 = zeros(0)
h0 = 0.0
μ0 = 0.0

θ0 = [q10, q20, u0, w0, μ0, h0]
norm(full_vector(mech.system))

mech.bodies[1]

z = minimal_to_maximal(mech, x_sol[1])
u = u_sol[1]
step!(mech, z, u, opts=SolverOptions(rtol=1e-6, btol=2e-4))

get_next_state(mech)
function action(mechanism::Mechanism)
    a = 0.0

    return a
end

a = action(mech)

function constraint_jacobian_configuration(relative::Symbol, joint::Joint{T,Nλ,Nb},
        xa::AbstractVector, qa::Quaternion,
        xb::AbstractVector, qb::Quaternion,
        η) where {T,Nλ,Nb}

    ∇comp = szeros(T,Nb,7)
    ∇mincoord = minimal_coordinates_jacobian_configuration(relative, joint, xa, qa, xb, qb, attjac=false)
    ∇unlim = joint_constraint_jacobian_configuration(relative, joint, xa, qa, xb, qb, η)

    return [∇comp; ∇mincoord; -∇mincoord; ∇unlim]
end

relative = :parent
joint = mech.joints[4].rotational
xa = srand(3)
xb = srand(3)
qa = Quaternion(normalize(srand(4))...)
qb = Quaternion(normalize(srand(4))...)
η = 0.0

constraint_jacobian_configuration(relative, joint, xa, qa, xb, qb, η)
@benchmark $constraint_jacobian_configuration($relative, $joint, $xa, $qa, $xb, $qb, $η)

function impulse_map(relative::Symbol, joint::Joint{T,Nλ,0},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion,
    η) where {T,Nλ}
    J = constraint_jacobian_configuration(relative, joint, xa, qa, xb, qb, η)
    G = [[Diagonal(sones(T,3)) szeros(T,3,3)]; [szeros(4,3) LVᵀmat(relative == :parent ? qa : qb)]]
    return Diagonal([sones(3); 0.5 * sones(3)]) * transpose(J * G)
end

relative = :parent
joint = mech.joints[4].translational
xa = srand(3)
xb = srand(3)
qa = Quaternion(normalize(srand(4))...)
qb = Quaternion(normalize(srand(4))...)
η = 0.0

impulse_map(relative, joint, xa, qa, xb, qb, η)
@benchmark $impulse_map($relative, $joint, $xa, $qa, $xb, $qb, $η)

function joint_constraint_jacobian_configuration(relative::Symbol, joint::Joint{T},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion, η) where {T}
    X, Q = displacement_jacobian_configuration(relative, joint, xa, qa, xb, qb, attjac=false)
    return constraint_mask(joint) * [X Q]
end

relative = :parent
joint = mech.joints[4].translational
xa = srand(3)
xb = srand(3)
qa = Quaternion(normalize(srand(4))...)
qb = Quaternion(normalize(srand(4))...)
η = 0.0
joint_constraint_jacobian_configuration(relative, joint, xa, qa, xb, qb, η)
@benchmark $joint_constraint_jacobian_configuration($relative, $joint, $xa, $qa, $xb, $qb, $η)

function displacement_jacobian_configuration(relative::Symbol, joint::Rotational{T},
        xa::AbstractVector, qa::Quaternion,
        xb::AbstractVector, qb::Quaternion;
        attjac::Bool=true, vmat=true) where T
    X = szeros(T, 3, 3)
    if relative == :parent
		Q = Lᵀmat(joint.orientation_offset) * Rmat(qb) * Tmat()
		attjac && (Q *= LVᵀmat(qa))
    elseif relative == :child
		Q = Lᵀmat(joint.orientation_offset) * Lᵀmat(qa)
		attjac && (Q *= LVᵀmat(qb))
	end
	# vmat && (Q = Vmat() * Q)
	return X, Q
end

relative = :parent
joint = mech.joints[4].rotational
xa = srand(3)
xb = srand(3)
qa = Quaternion(normalize(srand(4))...)
qb = Quaternion(normalize(srand(4))...)
η = 0.0
displacement_jacobian_configuration(relative, joint, xa, qa, xb, qb)
Main.@code_warntype displacement_jacobian_configuration(relative, joint, xa, qa, xb, qb)
@benchmark $displacement_jacobian_configuration($relative, $joint, $xa, $qa, $xb, $qb)

function displacement_jacobian_configuration(relative::Symbol, joint::Translational{T},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion;
    attjac=true) where T

    vertices = joint.vertices

    if relative == :parent
        d = xb + vector_rotate(vertices[2], qb) - (xa + vector_rotate(vertices[1], qa)) # in the world frame
        X = -rotation_matrix(inv(qa))
        Q = -rotation_matrix(inv(qa)) * ∂rotation_matrix∂q(qa, vertices[1])
        Q += ∂rotation_matrix_inv∂q(qa, d)
        attjac && (Q *= LVᵀmat(qa))
    elseif relative == :child
        X = rotation_matrix(inv(qa))
        Q = rotation_matrix(inv(qa)) * ∂rotation_matrix∂q(qb, vertices[2])
        attjac && (Q *= LVᵀmat(qb))
    end

    return X, Q
end

displacement_jacobian_configuration

joint_constraint_jacobian_configuration
minimal_coordinates_jacobian_configuration
minimal_coordinates_jacobian_configuration
minimal_velocities_jacobian_configuration
minimal_velocities_jacobian_velocity


text:
# ---
# title: 164. Maximum Gap
# id: problem164
# author: AquaIndigo
# date: 2020-11-26
# difficulty: Hard
# categories: Sort
# link: <https://leetcode.com/problems/maximum-gap/description/>
# hidden: true
# ---
#
# Given an unsorted array, find the maximum difference between the successive
# elements in its sorted form.
#
# Return 0 if the array contains less than 2 elements.
#
# **Example 1:**
#
#
#
#     Input: [3,6,9,1]
#     Output: 3
#     Explanation: The sorted form of the array is [1,3,6,9], either
#                  (3,6) or (6,9) has the maximum difference 3.
#
# **Example 2:**
#
#
#
#     Input: [10]
#     Output: 0
#     Explanation: The array contains less than 2 elements, therefore return 0.
#
# **Note:**
#
#   * You may assume all elements in the array are non-negative integers and fit in the 32-bit signed integer range.
#   * Try to solve it in linear time/space.
#
#
## @lc code=start
using LeetCode

function maximum_gap(nums::Vector{Int})::Int
    sort!(nums)
    return length(nums) < 2 ? 0 : maximum(@view(nums[2:end]) - @view(nums[1:(end - 1)]))
end
## @lc code=end


text:
using DataFrames
using Distributions
using Random

function simple_run(individuals_df::DataFrame;
        tracking_prob::Float64=0.0,
        constant_kernel_param::Float64=1.0,
        tracking_delay::Float64=1.0,
        mild_detection_prob=0.0,
        history::Union{Nothing, Vector{Simulation.Event}}=nothing,
        execution_history::Union{Nothing, BitVector}=nothing,
        state_history::Union{Nothing, Vector{Simulation.IndividualState}}=nothing,
        debug_return=false,
        seed=123,
        initial_infections::Integer=100
    )
    rng = MersenneTwister(seed)
    params = Simulation.load_params(
        rng,
        population=individuals_df,

        constant_kernel_param=constant_kernel_param,
        household_kernel_param=1.0,

        backward_tracking_prob=tracking_prob,
        backward_detection_delay=tracking_delay/2,

        forward_tracking_prob=tracking_prob,
        forward_detection_delay=tracking_delay/2,

        testing_time=tracking_delay/2,

        mild_detection_prob = mild_detection_prob
    );
    state = Simulation.SimState(params.progressions |> length, seed=seed)

    sample(rng, 1:length(params.progressions), initial_infections) .|> person_id ->
        push!(state.queue, Simulation.Event(Val(Simulation.OutsideInfectionEvent), 0.0, person_id))

    @time Simulation.simulate!(
        state,
        params,
        history=history,
        execution_history=execution_history,
        state_history=state_history
    )

    return state
end

function simple_run!(state::Simulation.SimState, individuals_df::DataFrame;
        tracking_prob::Float64=0.0,
        constant_kernel_param::Float64=1.0,
        household_kernel_param::Real=0.3,

        tracking_delay::Float64=2.5,
        forward_detection_delay::Union{Nothing,Float64}=nothing,
        backward_detection_delay::Union{Nothing, Float64}=nothing,
        testing_time::Union{Nothing, Float64}=nothing,

        mild_detection_prob=0.0,
        phone_tracking_usage=0.0,
        phone_detection_delay=0.25,
        hospital_kernel_param::Float64=0.0,
        history::Union{Nothing, Vector{Simulation.Event}}=nothing,
        execution_history::Union{Nothing, BitVector}=nothing,
        state_history::Union{Nothing, Vector{Simulation.IndividualState}}=nothing,
        debug_return=false,
        seed=123,
        initial_infections::Integer=100
    )
    Simulation.reset!(state)

    state.rng = MersenneTwister(seed)

    params = Simulation.load_params(
        state.rng,
        population=individuals_df,

        mild_detection_prob = mild_detection_prob,

        constant_kernel_param=constant_kernel_param,
        household_kernel_param=household_kernel_param,
        hospital_kernel_param=hospital_kernel_param,

        backward_tracking_prob=tracking_prob,
        backward_detection_delay= nothing==backward_detection_delay ? tracking_delay/2 : backward_detection_delay,

        forward_tracking_prob=tracking_prob,
        forward_detection_delay= nothing==forward_detection_delay ? tracking_delay/2 : forward_detection_delay,

        testing_time= nothing==testing_time ? tracking_delay/2 : testing_time,
        phone_tracking_usage=phone_tracking_usage,
        phone_detection_delay=phone_detection_delay
    );

    Simulation.initialfeed!(state, initial_infections)

    Simulation.simulate!(
        state,
        params,
        history=history,
        execution_history=execution_history,
        state_history=state_history
    )
    return state
end

function merge_history(history::Vector{Simulation.Event}, execution_history::BitVector, state_history::Vector{Simulation.IndividualState}; pop_last_event=false)
    history = copy(history)
    if pop_last_event
        last_event = pop!(history)
    end
    zip(history, state_history, execution_history)
end

function state2plot(state::Simulation.SimState)
    transmissions = vcat(state.forest.infections...)
    sort!(transmissions, lt=(x,y)->x.time<y.time)
    times = getproperty.(transmissions, :time)
    points = 1:length(transmissions)
    return times, points
end


text:
using FactoredValueMCTS
using Test

using POMDPs
using MultiAgentSysAdmin
using MultiUAVDelivery

@testset "FactoredValueMCTS.jl" begin

    @testset "varel" begin
        @testset "sysadmin" begin
            @testset "local" begin
                mdp = BiSysAdmin{false}()
                solver = FVMCTSSolver()
                planner = solve(solver, mdp)
                s = rand(initialstate(mdp))
                a = action(planner, s)
                @test a isa actiontype(mdp)
            end

            @testset "global" begin
                mdp = BiSysAdmin{true}()
                solver = FVMCTSSolver()
                planner = solve(solver, mdp)
                s = rand(initialstate(mdp))
                a = action(planner, s)
                @test a isa actiontype(mdp)
            end
        end

    end

    @testset "maxplus" begin
         @testset "sysadmin" begin
            @testset "local" begin
                mdp = BiSysAdmin{false}()
                solver = FVMCTSSolver(;coordination_strategy=MaxPlus())
                planner = solve(solver, mdp)
                s = rand(initialstate(mdp))
                a = action(planner, s)
                @test a isa actiontype(mdp)
            end

            @testset "global" begin
                mdp = BiSysAdmin{true}()
                solver = FVMCTSSolver(;coordination_strategy=MaxPlus())
                planner = solve(solver, mdp)
                s = rand(initialstate(mdp))
                a = action(planner, s)
                @test a isa actiontype(mdp)
            end
        end

        @testset "uav" begin
            mdp = FirstOrderMultiUAVDelivery()
            solver = FVMCTSSolver(;coordination_strategy=MaxPlus())
            planner = solve(solver, mdp)
            s = rand(initialstate(mdp))
            a = action(planner, s)
            @test a isa actiontype(mdp)
        end

    end

end


text:
function luhn(numb)
    i = 1
    total = 0
    numba = replace(numb, " " => "")
    numb = replace(numba, !isdigit => "")

    if length(numba) == length(numb) && length(numba)!=1
        for c in reverse(numb)
            tmp = Int(c) - Int('0')
            if i % 2 == 0
                db = tmp * 2
                if db > 9
                    db = db - 9
                end
                tmp = db
            end
            i += 1
            total = total + tmp

        end
    else
        #we force return false as there were nodigits
        total = 1
    end

    ret = false
    if total % 10 == 0
        ret = true
    end
    return ret
end

max_stars_repo_name: stepinski/julia-init

max_issues_repo_name: stepinski/julia-init

max_forks_repo_name: stepinski/julia-init


text:
"""
Financial time series type that includes the instrument (stock, currency pair) metadata.
"""
typealias FinancialTimeSeries{T<:Float64,N,M<:Union(AbstractInstrument, AbstractCurrency)} TimeArray{T,N,M}

function show(io::IO, ft::FinancialTimeSeries)
  # variables
  nrow          = size(ft.values, 1)
  ncol          = size(ft.values, 2)
  intcatcher    = falses(ncol)
  for c in 1:ncol
      rowcheck =  trunc(ft.values[:,c]) - ft.values[:,c] .== 0
      if sum(rowcheck) == length(rowcheck)
          intcatcher[c] = true
      end
  end
  spacetime     = strwidth(string(ft.timestamp[1])) + 3
  firstcolwidth = strwidth(ft.colnames[1])
  colwidth      = Int[]
      for m in 1:ncol
          push!(colwidth, max(strwidth(ft.colnames[m]), strwidth(@sprintf("%.2f", maximum(ft.values[:,m])))))
      end

  # summary line
  #println(io,@sprintf("%dx%d FinancialTimeSeries: %s %s to %s", nrow, ncol, ft.meta.ticker, string(ft.timestamp[1]), string(ft.timestamp[nrow])))
  println(io,@sprintf("%dx%d FinancialTimeSeries for %s, %s to %s", nrow, ncol, string(typeof(ft.meta)), string(ft.timestamp[1]), string(ft.timestamp[nrow])))
  println(io,"")

  # instrument
  Base.show(io, ft.meta)
  println(io, "")

  # row label line
   print(io, ^(" ", spacetime), ft.colnames[1], ^(" ", colwidth[1] + 2 -firstcolwidth))

   for p in 2:length(colwidth)
     print(io, ft.colnames[p], ^(" ", colwidth[p] - strwidth(ft.colnames[p]) + 2))
   end
   println(io,"")

  # timesftmp and values line
    if nrow > 7
        for i in 1:4
            print(io, ft.timestamp[i], " | ")
        for j in 1:ncol
            intcatcher[j] ?
            print(io, rpad(iround(ft.values[i,j]), colwidth[j] + 2, " ")) :
            print(io, rpad(round(ft.values[i,j], 2), colwidth[j] + 2, " "))
        end
        println(io,"")
        end
        println(io,'\u22EE')
        for i in nrow-3:nrow
            print(io, ft.timestamp[i], " | ")
        for j in 1:ncol
            intcatcher[j] ?
            print(io, rpad(iround(ft.values[i,j]), colwidth[j] + 2, " ")) :
            print(io, rpad(round(ft.values[i,j], 2), colwidth[j] + 2, " "))
        end
        println(io,"")
        end
    else
        for i in 1:nrow
            print(io, ft.timestamp[i], " | ")
        for j in 1:ncol
            intcatcher[j] ?
            print(io, rpad(iround(ft.values[i,j]), colwidth[j] + 2, " ")) :
            print(io, rpad(round(ft.values[i,j], 2), colwidth[j] + 2, " "))
        end
        println(io,"")
        end
    end
end


text:
@testset "Composition" begin
  # basic tests
  c = Composition(1,1,1)
  @test parts(c) == (:part1,:part2,:part3)
  @test components(c) == [1,1,1]
  c = Composition(a=1, b=missing)
  @test parts(c) == (:a,:b)
  @test isequal(components(c), [1,missing])
  for c in [Composition((part1=1,part2=2))
            Composition((:part1,:part2), (1,2))
            Composition((1,2))
            Composition(1,2)]
    @test parts(c) == (:part1,:part2)
    @test components(c) == [1,2]
  end

  # equality
  c₀ = Composition(1,1,1)
  c₁ = Composition(1,2,3)
  c₂ = Composition(10,20,30)
  @test c₁ == c₂
  @test c₁ - c₂ == c₀
  @test norm(c₀) ≈ 0.
  @test isapprox(distance(c₁, c₁), 0., atol=1e-6)

  # make sure names are preserved
  c = Composition(a=1.0, b=2.0)
  @test parts(c + c) == parts(c)
  @test parts(-c) == parts(c)
  @test parts(c - c) == parts(c)
  @test parts(2c) == parts(c)

  # unicode names work fine
  c = Composition(CO₂=1.0, CH₄=0.1, N₂O=0.1)
  @test c.CO₂ == 1.0
  @test c.CH₄ == 0.1
  @test c.N₂O == 0.1

  # get part by name
  c = Composition(a=3,b=2,c=1)
  @test c.a == 3
  @test c.b == 2
  @test c.c == 1

  # missing parts
  c = Composition(a=1.0, b=2.0, c=missing)
  @test c.a == 1
  @test c.b == 2
  @test ismissing(c.c)
end


text:
using Documenter
using LuxCat

makedocs(;
    modules = [LuxCat],
    authors = "Jacob Zelko <jacobszelko@gmail.com> and contributors",
    repo = "https://github.com/TheCedarPrince/LuxCat.jl/blob/{commit}{path}#L{line}",
    sitename = "LuxCat.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        assets = String[],
    ),
    pages = [
        "Home" => "index.md",
    ],
)

deploydocs(; devbranch = "main", repo = "github.com/TheCedarPrince/LuxCat.jl.git", push_preview = false, target = "build")