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

using MGVI

using LinearAlgebra
using Random
using SparseArrays
using ValueShapes

Test.@testset "test_fisher_values" begin

    include("information_utils.jl")

    Random.seed!(42)
    epsilon = 5E-2
    num_runs = 100000

    function test_univariate(dist, params)
        model = p -> dist(p...)
        res = MGVI.fisher_information(model(params))
        truth = fisher_information_mc(model, params, num_runs)
        Test.@test norm((Matrix(res) - truth)) / norm(truth) < epsilon
    end

    test_univariate(Normal, [0.1, 0.2])
    test_univariate(Exponential, [0.3])
    test_univariate(Poisson, [5.75])

    function test_mvnormal(dim)
        dim = 1
        cov = I*5 + Symmetric(rand(dim, dim))
        mean = rand(dim)
        params = vcat(mean, cov[:])
        model = p -> MvNormal(p[1:dim], reshape(p[dim+1:end], dim, dim))
        res = MGVI.fisher_information(model(params))
        explicit = explicit_mv_normal_fi(cov)
        truth = fisher_information_mc(model, params, num_runs)
        Test.@test norm((Matrix(res) - truth)) / norm(truth) < epsilon
        Test.@test norm((Matrix(res) - explicit)) / norm(explicit) < epsilon
    end

    test_mvnormal(1)  # test mvnormal 1d
    test_mvnormal(2)  # test mvnormal 2d
    test_mvnormal(3)  # test mvnormal 3d

end

Test.@testset "test_fisher_information_combinations" begin

    epsilon = 1E-5

    # test Product(Univariates)
    μ1, σ1 = 0.1, 0.2
    μ2, σ2 = 0.1, 0.3
    dists = [Normal(μ1, σ1), Normal(μ2, σ2)]
    res = MGVI.fisher_information(Product(dists))
    truth = blockdiag(MGVI.fisher_information.(dists)...)
    Test.@test norm(Matrix(res) - Matrix(truth)) < epsilon

    # test NamedTupleDist
    dists = NamedTupleDist(a=Normal(0.1, 0.2),
                           b=Product([Normal(0.1, 0.2), Normal(0.3, 0.1)]),
                           c=MvNormal([0.2, 0.3], [2. 0.1; 0.1 4.5]))
    res = MGVI.fisher_information(dists)
    truth = blockdiag((MGVI.without_chol ∘ MGVI.fisher_information).(values(dists))...)
    Test.@test norm(Matrix(res) - Matrix(truth)) < epsilon

end


text:
module TransformationElasticitySolver
using CutCellDG
using LinearAlgebra
include("analytical-solver.jl")
PS = PlaneStrainSolver

function bulk_modulus(l, m)
    return l + 2m / 3
end

function onboundary(x, L, W)
    return x[2] ≈ 0.0 || x[1] ≈ L || x[2] ≈ W || x[1] ≈ 0.0
end

function nodal_displacement(
    mesh,
    basis,
    cellquads,
    facequads,
    interfacequads,
    stiffness,
    theta0,
    boundarydisplacement,
    penalty;
    eta = 1,
)

    L, W = CutCellDG.mesh_widths(mesh)
    lambda1, mu1 = CutCellDG.lame_coefficients(stiffness, +1)
    transfstress =
        CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)

    sysmatrix = CutCellDG.SystemMatrix()
    sysrhs = CutCellDG.SystemRHS()

    CutCellDG.assemble_displacement_bilinear_forms!(
        sysmatrix,
        basis,
        cellquads,
        stiffness,
        mesh,
    )
    CutCellDG.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        cellquads,
        mesh,
    )
    CutCellDG.assemble_interelement_condition!(
        sysmatrix,
        basis,
        facequads,
        stiffness,
        mesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_interelement_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        mesh,
    )

    CutCellDG.assemble_coherent_interface_condition!(
        sysmatrix,
        basis,
        interfacequads,
        stiffness,
        mesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_coherent_interface_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        interfacequads,
        mesh,
    )

    CutCellDG.assemble_penalty_displacement_bc!(
        sysmatrix,
        sysrhs,
        boundarydisplacement,
        basis,
        facequads,
        stiffness,
        mesh,
        x -> onboundary(x, L, W),
        penalty,
    )
    CutCellDG.assemble_penalty_displacement_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        mesh,
        x -> onboundary(x, L, W),
    )

    matrix = CutCellDG.sparse_displacement_operator(sysmatrix, mesh)
    rhs = CutCellDG.displacement_rhs_vector(sysrhs, mesh)

    solution = matrix \ rhs

    return solution
end

function construct_merged_mesh_and_quadratures(
    cutmesh,
    levelset,
    numqp;
    tinyratio = 0.3,
)
    cellquads = CutCellDG.CellQuadratures(cutmesh, levelset, numqp)
    interfacequads = CutCellDG.InterfaceQuadratures(cutmesh, levelset, numqp)
    facequads = CutCellDG.FaceQuadratures(cutmesh, levelset, numqp)

    mergedmesh = CutCellDG.MergedMesh!(
        cutmesh,
        cellquads,
        facequads,
        interfacequads,
        tinyratio = tinyratio,
    )

    return mergedmesh, cellquads, facequads, interfacequads
end

function parent_stress_strain(
    nodaldisplacement,
    basis,
    refpoints,
    refcellids,
    mesh,
    stiffness,
)
    parentstrain = CutCellDG.parent_strain(
        nodaldisplacement,
        basis,
        refpoints,
        refcellids,
        mesh,
    )
    parentstress = CutCellDG.parent_stress(parentstrain, stiffness)
    # parentstrainenergy =
    #     V0 * CutCellDG.strain_energy(parentstress, parentstrain)
    #
    # parentradialtraction =
    #     CutCellDG.traction_force_at_points(parentstress, normals)
    # parentsrr = CutCellDG.traction_component(parentradialtraction, normals)
    # parentdilatation = CutCellDG.dilatation(parentstrain)
    #
    # parentcompwork = V0 * (1.0 .+ parentdilatation) .* parentsrr
    #
    # parentpotential = parentstrainenergy - parentcompwork
    return parentstress, parentstrain
end

function product_stress_strain(
    nodaldisplacement,
    basis,
    refpoints,
    refcellids,
    mesh,
    stiffness,
    theta0,
)

    productstrain = CutCellDG.product_elastic_strain(
        nodaldisplacement,
        basis,
        theta0,
        refpoints,
        refcellids,
        mesh,
    )
    productstress = CutCellDG.product_stress(productstrain, stiffness, theta0)

    # productstrainenergy =
    #     V0 * CutCellDG.strain_energy(productstress, productstrain)
    # productradialtraction =
    #     CutCellDG.traction_force_at_points(productstress, normals)
    # productsrr = CutCellDG.traction_component(productradialtraction, normals)
    # productdilatation = CutCellDG.dilatation(productstrain)
    # productcompwork = V0 * (1.0 .+ productdilatation) .* productsrr
    #
    # productpotential = productstrainenergy - productcompwork

    return productstress, productstrain
end

function potential_difference_at_closest_points(
    querypoints,
    nodaldisplacement,
    basis,
    spatialseedpoints,
    seedcellids,
    mesh,
    levelset,
    stiffness,
    theta0,
    V01,
    V02,
    tol,
    boundingradius,
)

    closestpoints, closestcellids = CutCellDG.closest_points_on_zero_levelset(
        querypoints,
        spatialseedpoints,
        seedcellids,
        levelset,
        tol,
        boundingradius,
    )

    normals = CutCellDG.collect_normals_at_spatial_points(
        closestpoints,
        closestcellids,
        levelset,
    )

    parentclosestrefpoints = CutCellDG.map_to_reference_on_merged_mesh(
        closestpoints,
        closestcellids,
        -1,
        mesh,
    )
    productclosestrefpoints = CutCellDG.map_to_reference_on_merged_mesh(
        closestpoints,
        closestcellids,
        +1,
        mesh,
    )

    # parentpotential = parent_potential(
    #     nodaldisplacement,
    #     basis,
    #     parentclosestrefpoints,
    #     closestcellids,
    #     normals,
    #     mesh,
    #     stiffness,
    #     V02,
    # )
    parentstress, parentstrain = parent_stress_strain(
        nodaldisplacement,
        basis,
        parentclosestrefpoints,
        closestcellids,
        mesh,
        stiffness,
    )
    # productpotential = product_potential(
    #     nodaldisplacement,
    #     basis,
    #     productclosestrefpoints,
    #     closestcellids,
    #     normals,
    #     mesh,
    #     stiffness,
    #     theta0,
    #     V01,
    # )
    productstress, productstrain = product_stress_strain(
        nodaldisplacement,
        basis,
        productclosestrefpoints,
        closestcellids,
        mesh,
        stiffness,
        theta0,
    )

    parentstrainenergy =
        V02 * CutCellDG.strain_energy(parentstress, parentstrain)
    parentradialtraction =
        CutCellDG.traction_force_at_points(parentstress, normals)
    parentsrr = CutCellDG.traction_component(parentradialtraction, normals)
    parentdilatation = CutCellDG.dilatation(parentstrain)

    productstrainenergy =
        V01 * CutCellDG.strain_energy(productstress, productstrain)
    productradialtraction =
        CutCellDG.traction_force_at_points(productstress, normals)
    productsrr = CutCellDG.traction_component(productradialtraction, normals)
    productdilatation = CutCellDG.dilatation(productstrain)

    meansrr = 0.5 * (parentsrr + productsrr)
    parentsrr = productsrr = meansrr

    strainenergyjump = productstrainenergy - parentstrainenergy
    compworkjump = (
        V01 * (1.0 .+ productdilatation) .* productsrr -
        V02 * (1.0 .+ parentdilatation) .* parentsrr
    )

    pd = strainenergyjump - compworkjump

    return pd
end

function average(v)
    return sum(v) / length(v)
end

function time_step_size(levelsetspeed, dx, CFL)
    s = maximum(abs.(levelsetspeed))
    return CFL * dx / s
end

function potential_difference_at_nodal_coordinates(
    cutmesh,
    basis,
    levelset,
    stiffness,
    theta0,
    boundarydisplacement,
    numqp,
    penalty,
    spatialseedpoints,
    seedcellids,
    V01,
    V02,
    ΔG0,
    tol,
    boundingradius,
)

    mesh, cellquads, facequads, interfacequads =
        construct_merged_mesh_and_quadratures(cutmesh, levelset, numqp)

    nodaldisplacement = nodal_displacement(
        mesh,
        basis,
        cellquads,
        facequads,
        interfacequads,
        stiffness,
        theta0,
        boundarydisplacement,
        penalty,
    )

    querypoints =
        CutCellDG.nodal_coordinates(CutCellDG.background_mesh(levelset))
    potentialdifference =
        (
            ΔG0 .+
            1e9 * potential_difference_at_closest_points(
                querypoints,
                nodaldisplacement,
                basis,
                spatialseedpoints,
                seedcellids,
                mesh,
                levelset,
                stiffness,
                theta0,
                V01,
                V02,
                tol,
                boundingradius,
            )
        ) / abs(ΔG0)

    return potentialdifference
end

function step_levelset(
    levelset,
    levelsetspeed,
    spatialseedpoints,
    seedcellids,
    paddedmesh,
    tol,
    boundingradius,
    CFL,
)

    paddedlevelset = CutCellDG.BoundaryPaddedLevelSet(
        paddedmesh,
        spatialseedpoints,
        seedcellids,
        levelset,
        tol,
        boundingradius,
    )

    dx = minimum(CutCellDG.grid_size(paddedmesh))
    dt = time_step_size(levelsetspeed, dx, CFL)

    newcoeffs =
        CutCellDG.step_first_order_levelset(paddedlevelset, levelsetspeed, dt)

    return newcoeffs
end

function step_interface(
    dgmesh,
    basis,
    levelset,
    stiffness,
    theta0,
    V01,
    V02,
    ΔG0,
    numqp,
    interfacecenter,
    outerradius;
    tol = 1e4eps(),
    boundingradius = 4.5,
    penaltyfactor = 1e3,
    CFL = 0.5,
    maxiter = 50,
)

    cutmesh = CutCellDG.CutMesh(dgmesh, levelset)
    refseedpoints, refseedcellids =
        CutCellDG.seed_zero_levelset(2, levelset, cutmesh)
    spatialseedpoints =
        CutCellDG.map_to_spatial(refseedpoints, refseedcellids, cutmesh)

    lambda1, mu1 = CutCellDG.lame_coefficients(stiffness, +1)
    lambda2, mu2 = CutCellDG.lame_coefficients(stiffness, -1)

    dx = minimum(CutCellDG.element_size(dgmesh))
    penalty = penaltyfactor / dx * 0.5 * (lambda1 + mu1 + lambda2 + mu2)

    interfaceradius =
        average(mapslices(norm, spatialseedpoints .- interfacecenter, dims = 1))
    analyticalsolution = PS.CylindricalSolver(
        interfaceradius,
        outerradius,
        interfacecenter,
        lambda1,
        mu1,
        lambda2,
        mu2,
        theta0,
    )

    potentialdifference = potential_difference_at_nodal_coordinates(
        cutmesh,
        basis,
        levelset,
        refseedpoints,
        refseedcellids,
        spatialseedpoints,
        stiffness,
        theta0,
        V01,
        V02,
        ΔG0,
        numqp,
        penalty,
        analyticalsolution,
        tol,
        boundingradius,
        maxiter,
    )

    newcoeffs = step_levelset(
        levelset,
        potentialdifference,
        refseedpoints,
        refseedcellids,
        spatialseedpoints,
        tol,
        boundingradius,
        CFL,
    )

    return newcoeffs
end

end

ext: jl


text:
# import methods that clash with Base and imported packages
import Base: flush

# server impleneation intentionally empty as this is just the client part


text:
using Test, InlineStrings, Parsers, Serialization, Random
import Parsers: SENTINEL, OK, EOF, OVERFLOW, QUOTED, DELIMITED, INVALID_QUOTED_FIELD, ESCAPED_STRING, NEWLINE, SUCCESS, peekbyte, incr!, checksentinel, checkdelim, checkcmtemptylines

@testset "InlineString basics" begin

y = "abcdef"
x = InlineString(y)
x, overflow = InlineStrings.addcodeunit(x, UInt8('g'))
@test !overflow
@test x == "abcdefg"
x, overflow = InlineStrings.addcodeunit(x, UInt8('g'))
@test overflow

x = InlineString("abc")
@test x == InlineString7(x) == InlineString15(x) == InlineString31(x) == InlineString63(x)
@test x == InlineString127(x) == InlineString255(x)
y = InlineString7(x)
@test InlineString3(y) == x
@test_throws ArgumentError InlineString3(InlineString("abcd"))
@test_throws ArgumentError InlineString1(InlineString("ab"))
x = InlineString("a")
y = InlineString7(x)
@test x == y
@test InlineString1(y) == x
@test InlineString1(x) == x

@test promote_type(InlineString1, InlineString3) === InlineString3
@test promote_type(InlineString1, InlineString255) === InlineString255
@test promote_type(InlineString31, InlineString127) === InlineString127
@test promote_type(InlineString255, InlineString7) === InlineString255
@test promote_type(InlineString63, InlineString15) === InlineString63

# Ensure we haven't caused ambiguity with Base.
# https://discourse.julialang.org/t/having-trouble-implementating-a-tables-jl-row-table-when-using-badukgoweiqitools-dataframe-tbl-no-longer-works/63622/1
@test promote_type(Union{}, String) == String

# make sure we don't read past an end of a buffer
buf = Vector{UInt8}("hey")
x = InlineString7(buf, 1, 3)
@test x == "hey"
@test typeof(x) == InlineString7

# https://github.com/JuliaData/WeakRefStrings.jl/issues/88
@test InlineString(String1("a")) === String1("a")

# https://github.com/JuliaData/InlineStrings.jl/issues/2
x = InlineString("hey")
@test typeof(string(x)) == String

# https://github.com/JuliaData/InlineStrings.jl/issues/8
# construction from pointer
ptrstr = "hey"
@test_throws ArgumentError String1(Ptr{UInt8}(0))
@test_throws ArgumentError String1(pointer(ptrstr), 2)
@test String1(pointer(ptrstr), 1) === String1("h")
ptrstr1 = UInt8['h', 0x00]
@test String1(pointer(ptrstr1)) === String1("h")
ptrstr2 = UInt8['h', 'e', 0x00]
@test_throws ArgumentError String1(pointer(ptrstr2))
@test String3(pointer(ptrstr1)) === String3("h")
ptrstr3 = UInt8['h', 'e', 'y', '1', 0x00]
@test_throws ArgumentError String3(pointer(ptrstr3))
@test_throws ArgumentError String3(pointer(ptrstr3), 4)
@test String3(pointer(ptrstr3), 3) === String3("hey")

end # @testset

@testset "InlineString operations" begin
    for y in ("", "🍕", "a", "a"^3, "a"^7, "a"^15, "a"^31, "a"^63, "a"^127, "a"^255)
        x = InlineString(y)
        @show typeof(x)
        @test codeunits(x) == codeunits(y)
        @test sizeof(x) == sizeof(y)
        @test ncodeunits(x) == ncodeunits(y)
        @test length(x) == length(y)
        @test codeunit(x) == UInt8
        @test lastindex(x) == lastindex(y)
        @test isempty(x) == isempty(y)
        @test String(x) === y
        @test Symbol(x) == Symbol(y)
        @test Vector{UInt8}(x) == Vector{UInt8}(y)
        @test Array{UInt8}(x) == Array{UInt8}(y)
        @test isascii(x) == isascii(y)
        @test x * x == y * y
        @test x^5 == y^5
        @test string(x) == string(y)
        @test join([x, x]) == join([y, y])
        @test reverse(x) == reverse(y)
        y != "" && @test startswith(x, "a") == startswith(y, "a")
        y != "" && @test endswith(x, "a") == endswith(y, "a")
        y != "" && @test findfirst(==(x[1]), x) === findfirst(==(x[1]), y)
        y != "" && @test findlast(==(x[1]), x) === findlast(==(x[1]), y)
        y != "" && @test findnext(==(x[1]), x, 1) === findnext(==(x[1]), y, 1)
        y != "" && @test findprev(==(x[1]), x, length(x)) === findprev(==(x[1]), y, length(x))
        @test lpad(x, 12) == lpad(y, 12)
        @test rpad(x, 12) == rpad(y, 12)
        y != "" && @test replace(x, x[1] => 'a') == replace(y, x[1] => 'a')
        r1 = match(Regex(x), x)
        r2 = match(Regex(y), y)
        @test r1 === r2 || r1.match == r2.match
        for i = 1:ncodeunits(x)
            @test codeunit(x, i) == codeunit(y, i)
            @test isvalid(x, i) == isvalid(y, i)
            @test thisind(x, i) == thisind(y, i)
            @test nextind(x, i) == nextind(y, i)
            @test prevind(x, i) == prevind(y, i)
            @test iterate(x, i) == iterate(y, i)
        end
        for i = 1:length(x)
            if isvalid(x, i)
                @test x[i] == y[i]
            end
        end
        @test x == x
        @test x == y
        @test y == x
        @test hash(x) == hash(y)
        @test cmp(x, x) == 0
        @test cmp(x, y) == 0
        @test cmp(y, x) == 0
        io = IOBuffer()
        write(io, x)
        seekstart(io)
        @test read(io, typeof(x)) === x
        @test chomp(x) == chomp(y)
        if typeof(x) != String255
            @test chomp(InlineString(x * "\n")) == chomp(y * "\n")
            @test chomp(InlineString(x * "\r\n")) == chomp(y * "\r\n")
        end
        @test chop(x) == chop(y)
        @test chop(x; head=1, tail=0) == chop(y; head=1, tail=0)
        @test chop(x; head=1) == chop(y; head=1)
        @test chop(x; head=sizeof(x) + 1) == chop(y; head=sizeof(y) + 1)
        @test chop(x; head=sizeof(x) + 1, tail=sizeof(x) + 1) == chop(y; head=sizeof(x) + 1, tail=sizeof(y) + 1)
        y != "" && @test first(x) == first(y)
        @test first(x, sizeof(x)) == first(y, sizeof(y))
        y != "" && @test first(x, sizeof(x) - 1) == first(y, sizeof(y) - 1)
        @test first(x, sizeof(x) + 1) == first(y, sizeof(y) + 1)
        y != "" && @test last(x) == last(y)
        @test last(x, sizeof(x)) == last(y, sizeof(y))
        y != "" && @test last(x, sizeof(x) - 1) == last(y, sizeof(y) - 1)
        @test last(x, sizeof(x) + 1) == last(y, sizeof(y) + 1)
    end
end

@testset "InlineString parsing" begin
testcases = [
    ("", InlineString7(""), NamedTuple(), OK | EOF),
    (" ", InlineString7(" "), NamedTuple(), OK | EOF),
    (" \"", InlineString7(), NamedTuple(), OK | QUOTED | EOF | INVALID_QUOTED_FIELD), # invalid quoted
    (" \"\" ", InlineString7(), NamedTuple(), OK | QUOTED | EOF), # quoted
    (" \" ", InlineString7(), NamedTuple(), OK | QUOTED | INVALID_QUOTED_FIELD | EOF), # invalid quoted
    (" \" \" ", InlineString7(" "), NamedTuple(), OK | QUOTED | EOF), # quoted
    ("NA", InlineString7(), (; sentinel=["NA"]), EOF | SENTINEL), # sentinel
    ("\"\"", InlineString7(), NamedTuple(), OK | QUOTED | EOF), # same e & cq
    ("\"\",", InlineString7(), NamedTuple(), OK | QUOTED | EOF | DELIMITED), # same e & cq
    ("\"\"\"\"", InlineString7("\""), NamedTuple(), OK | QUOTED | ESCAPED_STRING | EOF), # same e & cq
    ("\"\\", InlineString7(), (; escapechar=UInt8('\\')), OK | QUOTED | INVALID_QUOTED_FIELD | EOF), # \\ e, invalid quoted
    ("\"\\\"\"", InlineString7("\""), (; escapechar=UInt8('\\')), OK | QUOTED | ESCAPED_STRING | EOF), # \\ e, valid
    ("\"\"", InlineString7(), (; escapechar=UInt8('\\')), OK | QUOTED | EOF), # diff e & cq
    ("\"a", InlineString7(), NamedTuple(), OK | QUOTED | INVALID_QUOTED_FIELD | EOF), # invalid quoted
    ("\"a\"", InlineString7("a"), NamedTuple(), OK | QUOTED | EOF), # quoted
    ("\"a\" ", InlineString7("a"), NamedTuple(), OK | QUOTED | EOF), # quoted
    ("\"a\",", InlineString7("a"), NamedTuple(), OK | QUOTED | EOF | DELIMITED), # quoted
    ("a,", InlineString7("a"), NamedTuple(), OK | EOF | DELIMITED),
    ("a__", InlineString7("a"), (; delim="__"), OK | EOF | DELIMITED),
    ("a,", InlineString7("a"), (; ignorerepeated=true), OK | EOF | DELIMITED),
    ("a__", InlineString7("a"), (; delim="__", ignorerepeated=true), OK | EOF | DELIMITED),
    ("a\n", InlineString7("a"), (; ignorerepeated=true), OK | NEWLINE | EOF),
    ("a\r", InlineString7("a"), (; ignorerepeated=true), OK | NEWLINE | EOF),
    ("a\r\n", InlineString7("a"), (; ignorerepeated=true), OK | NEWLINE | EOF),
    ("a", InlineString7("a"), (; ignorerepeated=true), OK | EOF),
    ("a,,\n", InlineString7("a"), (; ignorerepeated=true), OK | DELIMITED | NEWLINE | EOF),
    ("a\n", InlineString7("a"), (; delim="__", ignorerepeated=true), OK | NEWLINE | EOF),
    ("a\r", InlineString7("a"), (; delim="__", ignorerepeated=true), OK | NEWLINE | EOF),
    ("a\r\n", InlineString7("a"), (; delim="__", ignorerepeated=true), OK | NEWLINE | EOF),
    ("a", InlineString7("a"), (; delim="__", ignorerepeated=true), OK | EOF),
    ("a____\n", InlineString7("a"), (; delim="__", ignorerepeated=true), OK | DELIMITED | NEWLINE | EOF),
    ("a\n", InlineString7("a"), NamedTuple(), OK | NEWLINE | EOF),
    ("a\r", InlineString7("a"), NamedTuple(), OK | NEWLINE | EOF),
    ("a\r\n", InlineString7("a"), NamedTuple(), OK | NEWLINE | EOF),
    ("abcdefg", InlineString7("abcdefg"), (; delim=nothing), OK | EOF),
    ("", InlineString7(), (; sentinel=missing), SENTINEL | EOF),
]

for (i, case) in enumerate(testcases)
    println("testing case = $i")
    buf, check, opts, checkcode = case
    res = Parsers.xparse(InlineString7, buf; opts...)
    @test check === res.val
    @test checkcode == res.code
end

res = Parsers.xparse(InlineString1, "")
@test Parsers.overflow(res.code)
res = Parsers.xparse(InlineString1, "ab")
@test Parsers.overflow(res.code)
res = Parsers.xparse(InlineString1, "b")
@test res.val === InlineString("b")

end # @testset

@testset "InlineString Serialization symmetry" begin
    for str in ("",  "🍕", "a", "a"^3, "a"^7, "a"^15, "a"^31, "a"^63, "a"^127, "a"^255)
        buf = IOBuffer()
        i_str = InlineString(str)
        Serialization.serialize(buf, i_str)
        seekstart(buf)
        i_str_copy = Serialization.deserialize(buf)

        @test typeof(i_str_copy) == typeof(i_str)
        @test i_str_copy == i_str
    end
end # @testset

@testset "alias tests" begin
    @test String1 == InlineString1
    @test String3 == InlineString3
    @test String7 == InlineString7
    @test String15 == InlineString15
    @test String31 == InlineString31
    @test String63 == InlineString63
    @test String127 == InlineString127
    @test String255 == InlineString255
end

@testset "sorting tests" begin
    for nelems in (50, 100, 500, 1000, 5000, 100_000)
        for T in (String1, String3, String7, String15, String31, String63, String127, String255)
            x = [randstring(rand(1:(max(1, sizeof(T) - 1)))) for _ = 1:nelems];
            y = map(T, x);
            @test sort(x) == sort(y)
            @test sort(x; rev=true) == sort(y; rev=true)
        end
    end
    x = [missing, String1("b"), String1("a")]
    sort!(x)
    @test isequal(x, ["a", "b", missing])
end

@testset "inlinestrings" begin

    @test inlinestrings([]) == []

    x = inlinestrings("$i" for i in (1, 10, 100))
    @test eltype(x) === String3
    @test x == ["1", "10", "100"]
    @test eltype(inlinestrings([randstring(i) for i = 1:3])) === String3
    @test eltype(inlinestrings([randstring(i) for i = 1:4])) === String7
    @test eltype(inlinestrings(["a", missing, "abcd"])) === Union{String7, Missing}
    @test eltype(inlinestrings([missing, "abcd"])) === Union{String7, Missing}

    # Base.SizeUnknown() case
    t() = true
    x = inlinestrings("$i" for i in (1, 10, 100) if t())
    @test eltype(x) === String3
    @test x == ["1", "10", "100"]
    @test eltype(inlinestrings(randstring(i) for i = 1:3 if t())) === String3
    @test eltype(inlinestrings(randstring(i) for i = 1:4 if t())) === String7
    @test eltype(inlinestrings(x for x in ["a", missing, "abcd"] if t())) === Union{String7, Missing}
    @test eltype(inlinestrings(x for x in [missing, "abcd"] if t())) === Union{String7, Missing}

    x = [randstring(i) for i = 1:31]
    @test InlineString.(x) == map(InlineString, x) == collect(InlineString, x)
    @test eltype(InlineString.(x)) == eltype(map(InlineString, x)) == eltype(collect(InlineString, x)) == InlineString31

    # promote all the way to String
    x = inlinestrings(randstring(i) for i = 1:256)
    @test eltype(x) === String
    x = inlinestrings(i == 1 ? missing : randstring(i) for i = 1:256 if t())
    @test eltype(x) === Union{Missing, String}
end


text:
using LinearAlgebra
using DoubleFloats
# using SemiLagMPI
using MPI
using SemiLagrangian

function printout(
    advd::AdvectionData{T,N,timeopt},
    str,
) where {T,N,timeopt}
    if timeopt != MPIOpt || advd.adv.mpid.ind == 1
        println(str)
    end
end
printout(str) = println(str)

# function trace_energy(
#     advd::AdvectionData{T,N,timeopt},
#     t,
# ) where {T,N,timeopt}
#     compute_charge!(advd)
#     compute_elfield!(advd)
#     # clockend(cl_obs,6)
#     # clockbegin(cl_obs,7)
#     elenergy = compute_ee(advd)
#     # clockend(cl_obs,7)
#     # clockbegin(cl_obs,8)
#     kinenergy = compute_ke(advd)
#     # clockend(cl_obs,8)
#     energyall = elenergy + kinenergy
#     return energyall
# end

function landau(advd::AdvectionData, nbdt)

    # global cl_obs
    # clockreset(cl_obs)
    maxdiff = 0
    dt = advd.adv.dt_base
    refel = getenergyall(advd)
 #   @show typeof(refel)
    #    printall(cl_obs)
    #    clockreset(cl_obs)
    maxel = refel
    minel = refel
    for i = 1:nbdt
        while advection!(advd)
        end
        el = getenergyall(advd)
 #       @show typeof(el)
        maxel = max(maxel ,el)
        minel = min(minel, el)
        maxdiff = maxel - minel
        # @show i,el,diff, maxdiff
        # flush(stdout)
        # printall(cl_obs)
        # clockreset(cl_obs)
    end
    # printall(cl_obs)
    return maxdiff
end

function landau2(
    t_max::T,
    timeopt,
    dt::T,
    sz,
    interp::AbstractInterpolation,
    type::TimeAlgorithm,
    epsilon::T,
    typeadd,
) where {T}
    nbdt = Int(round(t_max/dt))
    spmin, spmax, nsp = T(0), T(4big(pi)), sz[1]
    vmin, vmax, nv = -T(10), T(10), sz[2]

    mesh_sp = UniformMesh(spmin, spmax, nsp)
    mesh_v = UniformMesh(vmin, vmax, nv)

    tabst = [( [1,2], 2, 1, false) ]

    adv = Advection(
        (mesh_sp,mesh_v,),
        [interp,interp],
        dt,
        tabst,
        tab_coef=nosplit(dt),
        timeopt = timeopt,
        timealg=type,
        ordalg=typeadd)

    fct_sp(x) = epsilon * cos(x / 2) + 1
    fct_v(v) = exp(-v^2 / 2) / sqrt(2T(pi))

    lgn_sp = fct_sp.(mesh_sp.points)
    lgn_v = fct_v.(mesh_v.points)

    data = dotprod((lgn_sp, lgn_v))

    resint = sum(data)*(vmax-vmin)/length(data)
    data /= resint

    pvar = getpoissonvar(adv, type=StdPoisson2d, typeadd=typeadd)

    advd = AdvectionData(adv, data, pvar)

    if type == StdABinit
        # tabst_2 = [( [2,1], 1, 1, true, true),( [1,2], 1, 2, true, false) ]
        # adv_2 = Advection(
        # (mesh_sp,mesh_v,),
        # [interp,interp],
        # dt,
        # tabst_2,
        # tab_coef=hamsplit_3_11(dt),
        # timeopt = timeopt)

        # pvar_2 = getpoissonvar(adv_2)

        # advd_2 = AdvectionData(adv_2, data, pvar_2)

        # init_data = map(x->zeros(T,sz), 1:(typeadd-1))
        # for i=1:(typeadd-1)
        #     while advection!(advd_2)
        #     end
        #     pvar.init_data[i] .= advd_2.data
        # end
        for i = 1:(typeadd-1)
            pvar.init_data[i] .= data
        end
    end

    # advdata = Advection1dData(adv, data, pvar)

    return landau(advd, nbdt)
end

function landau1_1(
    t_max::T,
    timeopt,
    dt::T,
    sz,
    interp::AbstractInterpolation,
    tab_coef,
    epsilon::T
) where {T}
    nbdt = Int(round(t_max/dt))
    spmin, spmax, nsp = T(0), T(4big(pi)), sz[1]
    vmin, vmax, nv = -T(10), T(10), sz[2]

    mesh_sp = UniformMesh(spmin, spmax, nsp)
    mesh_v = UniformMesh(vmin, vmax, nv)

    tabst = [( [2,1], 1, 1, true),( [1,2], 1, 2, true) ]
#    tabst = [( [1,2], 1, 1, true),( [2,1], 1, 2, true) ]

    adv = Advection(
        (mesh_sp,mesh_v,),
        [interp,interp],
        dt,
        tabst,
        tab_coef=tab_coef,
        timeopt = timeopt)

    fct_sp(x) = epsilon * cos(x / 2) + 1
    fct_v(v) = exp(-v^2 / 2) / sqrt(2T(pi))

    lgn_sp = fct_sp.(mesh_sp.points)
    lgn_v = fct_v.(mesh_v.points)

    data = dotprod((lgn_sp, lgn_v))

    resint = sum(data)*(vmax-vmin)/length(data)
    data /= resint

    pvar = getpoissonvar(adv)

    advd = AdvectionData(adv, data, pvar)

    # advdata = Advection1dData(adv, data, pvar)

    return landau(advd, nbdt)
end

function run_mesure(
    t_max::T,
    timeopt,
    sz,
    interp,
    epsilon
) where{T}
# tabsplit = [standardsplit, strangsplit, triplejumpsplit, order6split, hamsplit_3_11]
# tabsplit = [standardsplit, strangsplit, triplejumpsplit, table2split]
tabsplit = [standardsplit, strangsplit, hamsplit_3_11]
# tabtype = [StdPoisson2d, StdAB, StdAB, StdAB, StdAB, StdAB, StdAB, StdAB, StdAB, StdAB, StdAB, StdAB, StdAB, StdAB, StdAB, StdAB, StdAB, StdAB]
# tabtype = [NoTimeAlg, ABTimeAlg,ABTimeAlg,ABTimeAlg,ABTimeAlg,ABTimeAlg,ABTimeAlg,]
# tabtype = [StdRK4,]
tabtype=[]
tabtypeadd = [0, 2, 3, 4, 5, 6]
# tabtxtsplit = ["stdsplit", "strangsplit", "triplejumpsplit", "order6split", "fernandosplit"]
# tabtxt = ["stdsplit", "strangsplit", "fernandosplit", "std2d", "stdAB_2", "stdAB_3", "stdAB_4", "stdAB_5", "stdAB_6",]
# tabtxt = ["stdsplit", "strangsplit", "fernandosplit", "stdRK4",]
tabtxt = ["stdsplit", "strangsplit", "fernandosplit"]
tabnbdt = [10,20,50,100,200,500,1000,2000,5000,10000,20000,50000,100000,200000,500000,1000000]

    res = zeros(Float64, length(tabtxt)+1, length(tabnbdt))

    res[1,:] .= Float64(t_max) ./ tabnbdt

    lenitc = length(tabsplit)
    for inbdt=1:length(tabnbdt), itc=1:length(tabtxt)
        tc = itc <= lenitc ? tabsplit[itc] : missing
        tp = itc > lenitc ? tabtype[itc-lenitc] : missing
        typadd = itc > lenitc ? tabtypeadd[itc-lenitc] : missing
        nbdt = tabnbdt[inbdt]
        dt = t_max/nbdt
        res[itc+1, inbdt] = if itc <= lenitc
            landau1_1(t_max, timeopt, dt, sz, interp, tc(dt), epsilon)
        else
            landau2(t_max, timeopt, dt, sz, interp, tp, epsilon, typadd)
        end
        if timeopt != MPIOpt || MPI.Comm_rank(MPI.COMM_WORLD) == 1
            println("# sz=$sz t_max=$t_max interp=$interp")
            for txt in tabtxt
                print("# $txt\t")
            end
            println("")
            for j=1:size(res,2),i=1:size(res,1)
                        print("$(res[i,j])")
                if i == size(res,1)
                    print("\n")
                else
                    print("\t")
                end
            end

	        println("free memory : $(Sys.free_memory()/2^30)")
            flush(stdout)
        end
    end
end
    # to run with mpi here is the command for n process (n must be a number) from SemiLagrangian.jl :
    # ./mpiexec.loc n notebooks/run_mesures_poisson.jl
    # else run like this :
    # julia --project=.  notebooks/run_mesures_poisson.jl

T=Double64
#run_mesure(T(1), NoTimeOpt, (128,128), Lagrange(11,T), T(0.5))
#run_mesure(T(1), MPIOpt, (128,128), Lagrange(11,T), T(0.5))
run_mesure(T(1), SimpleThreadsOpt, (128,128), Lagrange(11,T), T(0.5))


text:
module ProgressMeter

using Compat

export Progress, next!, cancel

type Progress
    n::Int
    dt::Float64
    counter::Int
    tfirst::Float64
    tlast::Float64
    printed::Bool    # true if we have issued at least one status update
    desc::String     # prefix to the percentage, e.g.  "Computing..."
    barlen::Int      # progress bar size (default is available terminal width)
    color::Symbol    # default to green

    function Progress(n::Integer, dt::Real = 1.0, desc::String = "Progress: ", barlen::Int = 0, color::Symbol = :green)
        this = new(convert(Int, n), convert(Float64, dt), 0)
        this.tfirst = time()
        this.tlast = this.tfirst
        this.printed = false
        this.desc = desc
        this.barlen = barlen
        this.color = color
        this
    end

    function Progress(n::Integer, desc::String = "Progress: ")
        this = new(convert(Int, n), convert(Float64, 0.01), 0)
        this.tfirst = time()
        this.tlast = this.tfirst
        this.printed = false
        this.desc = desc
        #...length of percentage and ETA string with days is 29 characters
        this.barlen = max(0, Base.tty_size()[2] - (length(desc)+29))
        this.color = :green
        this
    end
end

function next!(p::Progress)
    t = time()
    p.counter += 1
    if p.counter >= p.n
        if p.printed
            percentage_complete = 100.0 * p.counter / p.n
            bar = barstring(p.barlen, percentage_complete)
            dur = durationstring(t-p.tfirst)
            msg = @sprintf "%s%3u%%%s Time: %s" p.desc round(Int, percentage_complete) bar dur
            printover(msg, p.color)
            println()
        end
        return
    end

    if t > p.tlast+p.dt
        percentage_complete = 100.0 * p.counter / p.n
        bar = barstring(p.barlen, percentage_complete)
        elapsed_time = t - p.tfirst
        est_total_time = 100 * elapsed_time / percentage_complete
        eta_sec = round(Int, est_total_time - elapsed_time )
        eta = durationstring(eta_sec)
        msg = @sprintf "%s%3u%%%s  ETA: %s" p.desc round(Int, percentage_complete) bar eta
        printover(msg, p.color)
        # Compensate for any overhead of printing. This can be especially important
        # if you're running over a slow network connection.
        p.tlast = t + 2*(time()-t)
        p.printed = true
    end
end

function next!(p::Progress, color::Symbol)
  p.color = color
  next!(p)
end

function cancel(p::Progress, msg::String = "Aborted before all tasks were completed", color = :red)
    if p.printed
        printover(msg, color)
        println()
    end
    return
end

function printover(io::IO, s::String, color::Symbol = :color_normal)
    if isdefined(Main, :IJulia)
        print(io, "\r" * s)
    else
        print(io, "\u1b[1G")   # go to first column
        print_with_color(color, io, s)
        print(io, "\u1b[K")    # clear the rest of the line
    end
end

printover(s::String, color::Symbol = :color_normal) = printover(STDOUT, s, color)

function barstring(barlen, percentage_complete)
    bar = ""
    if barlen>0
        nsolid = round(Int, barlen * percentage_complete / 100)
        nempty = barlen - nsolid
        bar = string("|", repeat("#",nsolid), repeat(" ",nempty), "|")
    end
    bar
end

function durationstring(nsec)
    days = div(nsec, 60*60*24)
    r = nsec - 60*60*24*days
    hours = div(r,60*60)
    r = r - 60*60*hours
    minutes = div(r, 60)
    seconds = r - 60*minutes

    hhmmss = @sprintf "%u:%02u:%02u" hours minutes seconds
    if days>0
       return @sprintf "%u days, %s" days hhmmss
    end
    hhmmss
end

end


text:
using Gadfly

##

run_bench = true
addprocs(1)

## Setup

@everywhere begin

    import Base: minimum, string

    using BlackBoxOptimizationBenchmarking
    const BBOB = BlackBoxOptimizationBenchmarking

    include("/Users/jbieler/.julia/v0.6/CMAES/src/CMAES.jl")
    include("/Users/jbieler/.julia/v0.6/CMAES2/src/CMAES2.jl")
    include(joinpath(Pkg.dir(),"BlackBoxOptimizationBenchmarking/scripts/optimizers_interface.jl"))

    function optimize(::Type{CMAES2.CMAESOpt},f,D,run_length)
        try rm("CMAES.jld") end
        (CMAES2_res(), CMAES2.minimize(f, pinit(D), 5, fill(-100,D), fill(100,D); maxfevals = run_length)...)
    end

    struct CMAES2_res end

    xmin(mfit::Tuple{CMAES2_res,Array{Float64,1},Float64,Int64}) = mfit[3]
    minimizer(mfit::Tuple{CMAES2_res,Array{Float64,1},Float64,Int64}) = mfit[2]

    dimensions = [3]
    Ntrials = 3
    Δf = 1e-6
#    run_lengths = round.(Int,linspace(20,60_000,20))
    run_lengths = round.(Int,linspace(20,6_000,6))
    funcs = 1:6#length(enumerate(BBOBFunction))

end

string(opt::CMAES.CMA{T,K}) where {T,K} = string(typeof(opt).name,K)

optimizers = [
    CMAES.CMA(CMAES.FullCovMatrix,5.0,;dimension=3),
#    CMAES.CMA(CMAES.BD_CovMatrix,5.0,;dimension=3),
#    CMAES2.CMAESOpt,
#    CMAES.CMA(12,10.0),
#    PyCMA(),
]

opt_strings = map(string,optimizers)

## run benchmark

mean_succ, mean_dist, mean_fmin, runtime = BBOB.benchmark(
    optimizers, funcs, run_lengths, Ntrials, dimensions, Δf,
)

## make plots

outdir = joinpath(Pkg.dir(),"CMAES","plots")
mkpath(outdir)
cols = Colors.distinguishable_colors(size(mean_succ,1)+1,colorant"white")[2:end]

#all together

idx = sortperm([mean(mean_succ[i,:,end,:]) for i=1:length(optimizers)],rev=true)

p = plot(
    [layer(
        x=collect(run_lengths),y=mean(mean_succ[i,:,:,:],(1,3)),Geom.line,
        Theme(default_color=cols[i], line_style= mod(i,2)==0 ? :solid : :dash, line_width=2pt)
    ) for i in 1:size(mean_succ,1)]...,
    Coord.cartesian(ymax=1),
    Guide.title("All functions"), Guide.xlabel("Run Length"), Guide.ylabel("Success rate"),
    Guide.manual_color_key("", opt_strings[idx], cols[idx])
)

#draw(PNG(joinpath(outdir,"mean_succ.png"),dpi=150,16cm,10cm),p)

##
# each dimension
for (k,D) in zip(1:length(dimensions), dimensions)

    idx = sortperm([mean(mean_succ[i,:,end,k]) for i=1:length(optimizers)],rev=true)

    p = plot(
        [layer(
            x=collect(run_lengths),y=mean(mean_succ[i,:,:,k],1),Geom.line,
            Theme(default_color=cols[i], line_style= mod(i,2)==0 ? :solid : :dash, line_width=2pt)
        ) for i in 1:size(mean_succ,1)]...,
        Coord.cartesian(ymax=1),
        Guide.title("All functions, D: $D"), Guide.xlabel("Run Length"), Guide.ylabel("Success rate"),
        Guide.manual_color_key("", opt_strings[idx], cols[idx])
    )
#    draw(PNG(joinpath(outdir,"mean_succ_$(D).png"),dpi=150,16cm,10cm),p)

end

## runtime

idx = sortperm(mean(runtime,2)[:],rev=true)

p = plot(
    layer(
        y=opt_strings[idx], x=(mean(runtime,2)/Base.minimum(mean(runtime,2)))[idx],
        color = cols[idx],
        Geom.bar(orientation=:horizontal)
    ),
    Guide.title("All functions"), Guide.ylabel(""),
    Guide.xlabel("Relative Run Time",orientation=:horizontal),
    Scale.x_log10,
#    Scale.y_log10,
)

#draw(PNG(joinpath(outdir,"runtime.png"),dpi=150,16cm,10cm),p)
p

##


text:
import Base: >>>

export hasstate,
       clickable,
       selectable,
       draggable,
       resizable,
       leftbutton,
       rightbutton,
       scrollbutton

abstract Behaviour <: Tile

subscribe(t::Behaviour, s::Input; absorb=true) =
    subscribe(t, t.name, s, absorb=absorb)

(>>>)(t::Behaviour, s::Input) = subscribe(t, s)

immutable WithState{attr} <: Behaviour
    name::Symbol
    tile::Tile
    trigger::String
end

render{attr}(t::WithState{attr}) =
    render(t.tile) << Elem("watch-state",
        attributes=[:name=>t.name, :attr=>attr, :trigger=>t.trigger])

hasstate(tile::Tile; name=:_state, attr="value", trigger="change") =
    WithState{symbol(attr)}(name, tile, trigger)

# Sample a bunch of signals upon changes to another bunch of signals
# Returns a signal of dict of signal values
@api samplesignals => SignalSampler <: Behaviour begin
    arg(signals::AbstractArray)
    arg(triggers::AbstractArray)
    curry(tile::Tile)
    typedkwarg(name::Symbol=:_sampler)
end

samplesignals(tosample::Symbol, triggers::Symbol, x...; name=:_sampler) =
    samplesignals([tosample], [triggers], x...; name=name)
samplesignals(tosample::Symbol, triggers, x...; name=:_sampler) =
    samplesignals([tosample], [triggers], x...; name=name)
samplesignals(tosample, triggers::Symbol, x...; name=:_sampler) =
    samplesignals(tosample, [triggers], x...; name=name)

render(sig::SignalSampler) =
    render(sig.tile) <<
        Elem("signal-sampler",
            name=sig.name,
            signals=sig.signals,
            triggers=sig.triggers)

abstract MouseButton

@terms MouseButton begin
    leftbutton => LeftButton
    rightbutton => RightButton
    scrollbutton => ScrollButton
end

@api clickable => Clickable <: Behaviour begin
    typedarg(buttons::AbstractArray=[leftbutton])
    curry(tile::Tile)
    kwarg(name::Symbol=:_clicks)
end

button_number(::LeftButton) = 1
button_number(::RightButton) = 2
button_number(::ScrollButton) = 3

render(c::Clickable) =
    render(c.tile) << Elem("clickable-behaviour", name=c.name,
                        buttons=string(map(button_number, c.buttons)))

@api selectable => Selectable <: Behaviour begin
    curry(tile::Tile)
    kwarg(name::Symbol=:_clicks)
end

render(s::Selectable) =
    render(s.tile) << Elem("selectable-behaviour", name=s.name)

convert(::Type{MouseButton}, x::Int) =
    try [leftbutton, rightbutton, scrollbutton][x]
    catch error("Invalid mouse button code: $x")
    end

abstract MouseState

@terms MouseState begin
    mousedown => MouseDown
    mouseup => MouseUp
end

@api hoverable => Hoverable <: Behaviour begin
    typedarg(get_coords::Bool=false)
    curry(tile::Tile)
    kwarg(name::Symbol=:_hover)
end

immutable Hover
    state::MouseState
    position::(Float64, Float64)
end

immutable Resizable <: Behaviour
    name::Symbol
    tile::Tile
    container::Tile
end

immutable Size
    proportions::(Float64, Float64)
end

resizable(tile, container; name=:_size) = Resizable(name, tile, container)

immutable Draggable <: Behaviour
    name::Symbol
    tile::Tile
    container::Tile
end

draggable(tile, container; name=:_drag) = Draggable(name, tile, container)

immutable Drag
    fraction::(Float64, Float64)
end

immutable Sortable <: Behaviour
    name::Symbol
    tile::AbstractArray
end

immutable PositionSwap
    swap::(Int, Int)
end

immutable Editable <: Behaviour
    name::Symbol
    tile::Tile
end

## UI-side global channels
import Base: send, recv
export send, recv, wire

immutable ChanSend <: Tile
    chan::Symbol
    watch::Symbol
    tile::Tile
end
send(chan::Symbol, b::Behaviour) =
    ChanSend(chan, b.name, b)

render(chan::ChanSend) =
    render(chan.tile) <<
        Elem("chan-send",
            chan=chan.chan, watch=chan.watch)

immutable ChanRecv <: Tile
    chan::Symbol
    attr::Symbol
    tile::Tile
end
recv(chan::Symbol, t, attr) =
    ChanRecv(chan, attr, t)

render(chan::ChanRecv) =
    render(chan.tile) <<
        Elem("chan-recv",
            chan=chan.chan, attr=chan.attr)

wire(a::Behaviour, b, chan, attribute) =
    send(chan, a), recv(chan, b, attribute)


text:
mutable struct Entropy
    data::Ptr{Cvoid}
    sources::Vector{Any}

    function Entropy()
        ctx = new()
        ctx.data = Libc.malloc(100000)  # Exact byte count is 75088; playing it safe with some buffer
        ctx.sources = Any[]
        ccall((:mbedtls_entropy_init, libmbedcrypto), Cvoid, (Ptr{Cvoid},), ctx.data)
        finalizer(ctx->begin
            ccall((:mbedtls_entropy_free, libmbedcrypto), Cvoid, (Ptr{Cvoid},), ctx.data)
            Libc.free(ctx.data)
        end, ctx)
        ctx
    end
end

function add_source!(ctx::Entropy, f_source::Ptr, f, threshold, strong)
    ret = ccall((:mbedtls_entropy_add_source, libmbedcrypto), Cint,
      (Ptr{Cvoid}, Ptr{Cvoid}, Any, Csize_t, Cint),
      ctx.data, f_source, f, threshold, strong)
    Int(ret)
end

function jl_entropy(f, output, len, olen)
    output_jl = unsafe_wrap(Array, convert(Ptr{UInt8}, output), len)
    sz = f(output_jl)
    unsafe_store!(convert(Ptr{Csize_t}, olen), Csize_t(sz))
    return Cint(0)
end

function add_source!(ctx::Entropy, f::Function, threshold, strong)
    push!(ctx.sources, f)
    add_source!(ctx, c_entropy[], f, threshold, strong ? 1 : 0)
end

const c_entropy = Ref{Ptr{Cvoid}}(C_NULL)
function __entropyinit__()
    c_entropy[] = @cfunction(jl_entropy, Cint, (Any, Ptr{Cvoid}, Csize_t, Ptr{Cvoid}))
end

function gather(ctx::Entropy)
    @err_check ccall((:mbedtls_entropy_gather, libmbedcrypto), Cint,
      (Ptr{Cvoid},), ctx.data)
end

function update_manual(ctx::Entropy, data::Vector{UInt8})
    @err_check ccall((:mbedtls_entropy_update_manual, libmbedcrypto), Cint,
      (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), ctx.data, pointer(data), length(data))
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Graphics (Modia3D/Graphics/_module.jl)
#
# Content:
#   This file contains elements to show text in the visualization window

const fontFamilyPalette = Dict{String,String}([("FreeSans"     , "FreeSans"),
                                               ("FreeSerif"    , "FreeSerif"),
                                               ("Arial"        , "ARIAL"),
                                               ("ArialNarrow"  , "ARIALN"),
                                               ("CourierNew"   , "COUR"),
                                               ("TimesNewRoman", "Times"),
                                               ("Verdana"      , "verdana")
                                               ])

function getFontFileName(fontFamily::String, bold::Bool, italic::Bool)::String
   # Handle special cases first

   if fontFamily == "ArialNarrow" && bold && !italic
      fontFileName = "ARIALNB"
   elseif fontFamily == "Verdana"
      fontFileName = bold ? ( italic ? "verdanaz" : "verdanab") :
                            ( italic ? "verdanai" : "verdana" )
   else
      if !haskey(fontFamilyPalette, fontFamily)
         error("\nError from Modia3D.visualizer: fontFamily = \"", fontFamily, "\" is not supported.\n",
               "Supported font families: FreeSans, FreeSerif, Arial, ArialNarrow, CourierNew, TimesNewRoman, Verdana.")
      end
      name = fontFamilyPalette[fontFamily]
      fontFileName = bold ? ( italic ? name * "BI" : name * "BD" ) :
                            ( italic ? name * "I"  : name )
   end
   fontFileName = fontFileName * ".TTF"
end

"""
    font = Modia3D.Font(;fontFamily="FreeSans", bold=false, italic=false, charSize=0.1,
                         color="LightBlue", transparency=0.0)

Return a Font definition.

# Arguments
- `fontFamily::String`: Font family ("FreeSans", "FreeSerif", "Arial", "ArialNarrow", "CourierNew", "TimesNewRoman", or "Verdana").
- `bold::Bool`: = true, if bold font.
- `italic::Bool`: = true, if italic font.
- `charSize::Number`: Character size in [m].
- `color::Modia3D.RGBColor`: Color; Examples: rgb("Blue"), rgb([0,0,255]), rgb(0,0,255).
- `transparency::Number`: 0.0 (opaque) ... 1.0 (transparent)

# Examples
```julia
import Modia3D
font1 = Modia3D.Font()
font2 = Modia3D.Font(fontFamily="Arial", bold=true, charSize=0.2,
                     color="LightBlue", transparency=0.5)
```
"""
struct Font
   fontFamily::String
   bold::Bool
   italic::Bool
   charSize::Float64
   color::RGBColor
   transparency::Float64

   # Constructed by constructor
   fontFileName::String  # File name of font

   function Font(;fontFamily::String    = "FreeSans",
                  bold::Bool            = false,
                  italic::Bool          = false,
                  charSize::Number      = 0.1,
                  color                 = "LightBlue",
                  transparency::Number  = 0.0)
      @assert(charSize >= 0.0)
      @assert(transparency >= 0.0 && transparency <= 1.0)
      fontFileName = getFontFileName(fontFamily, bold, italic)

      new(fontFamily, bold, italic, charSize, rgb(color), transparency, fontFileName)
   end
end

@enum AxisAlignment::Cint Screen=0 XY_Plane=1 XZ_Plane=2 YZ_Plane=3
@enum Alignment::Cint     Left=1 Right=2 Center=3

"""
    textShape = Modia3D.TextShape(text; font=Modia3D.Font(), offset=[0.0,0.0,0.0],
                          axisAlignment=Modia3D.Screen, alignment=Modia3D.Center)

Return a text shape.

# Arguments
- `text::AbstractString`: String of the text.
- `font::`[`Modia3D.Font`](@ref): Font of the text.
- `offset::AbstractVector`: Offset from origin to text alignment point.
- `axisAlignment::Modia3D.AxisAlignment`: Alignment of Text (parallel to screen or in
   planes of frame: = Modia3D.Screen, Modia3D.XY_Plane, Modia3D.XZ_Plane, Modia3D.YZ_Plane).
- `axisAlignment::Modia3D.Alignment`: Alignment of Text relative to its origin
  (= Modia3D.Left, Modia3D.Right or Modia3D.Center).

# Examples
```julia
import Modia3D
font  = Modia3D.Font(fontFamily="Arial", charSize=0.4, color=Modia3D.rgb("Red"))
text1 = Modia3D.TextShape("This is a box")
text2 = Modia3D.TextShape("This is the xy plane";
                          font=font, axisAlignment=Modia3D.XY_Plane,
                          alignment=Modia3D.Left)
```
"""
mutable struct TextShape <: Modia3D.AbstractVisualElement
   text::String                   # String to be displayed
   font::Font                     # Font
   offset::MVector{3,Float64}     # Offset from origin to text alignment point
   axisAlignment::AxisAlignment   # Alignment of Text (parallel to screen or in planes of frame)
   alignment::Alignment           # Alignment of Text relative to its origin

   function TextShape(text::AbstractString;
                      font::Font=Font(),
                      offset::AbstractVector = Basics.zeroMVector(),
                      axisAlignment::AxisAlignment=Screen,
                      alignment::Alignment=Center)

      new(text,font,offset,axisAlignment,alignment)
   end
end


text:
using Mamba

include("utils.jl")

const openbugstests = [
  "asthma",
  "birats",
  "blocker",
  "bones",
  "dogs",
  "dyes",
  "epil",
  "equiv",
  "eyes",
  "inhalers",
  "jaws",
  "leuk",
  "lsat",
  "magnesium",
  "mice",
  "oxford",
  "pumps",
  "rats",
  "salm",
  "seeds",
  "stacks",
  "surgical"
]

const contributedtests = [
  "gk",
  "line_abc",
  "line_amwg_slice",
  "pollution"
]

println("Running examples:")

for t in openbugstests
  @everywhere srand(123)
  @runtest "../doc/examples/" t
  gelmandiag(sim) |> show
end

for t in contributedtests
  @everywhere srand(123)
  @runtest "../doc/examples/" t
end


text:
export Dirac

struct Dirac{X} <: PrimitiveMeasure
    x::X
end

function Pretty.quoteof(d::Dirac)
    q = Pretty.quoteof(d.x)
    :(Dirac($q))
end

sampletype(μ::Dirac{X}) where {X} = X

function (μ::Dirac{X})(s) where {X}
    μ.x ∈ s && return 1
    return 0
end

density(μ::Dirac, x) = x == μ.x

logdensity(μ::Dirac, x) = (x == μ.x) ? 0.0 : -Inf

Base.rand(::Random.AbstractRNG, T::Type, μ::Dirac) = μ.x

export dirac

dirac(d::AbstractMeasure) = Dirac(rand(d))

function logdensity(μ::Dirac{M}, ν::Dirac{M}, x) where {M}
    logdensity(μ, x) - logdensity(ν, x)
end

testvalue(d::Dirac) = d.x


text:
export BrickletVoltageVoltageCallbackThreshold
struct BrickletVoltageVoltageCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletVoltageAnalogValueCallbackThreshold
struct BrickletVoltageAnalogValueCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletVoltageIdentity
struct BrickletVoltageIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletVoltage
"""
Measures DC voltage between 0V and 50V
"""
mutable struct BrickletVoltage <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletVoltage(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_voltage")
        deviceInternal = package.BrickletVoltage(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_voltage
"""
    $(SIGNATURES)

Returns the voltage of the sensor.

If you want to get the voltage periodically, it is recommended to use the
:cb:`Voltage` callback and set the period with
:func:`Set Voltage Callback Period`.
"""
function get_voltage(device::BrickletVoltage)
    return device.deviceInternal.get_voltage()
end

export get_analog_value
"""
    $(SIGNATURES)

Returns the value as read by a 12-bit analog-to-digital converter.

.. note::
 The value returned by :func:`Get Voltage` is averaged over several samples
 to yield less noise, while :func:`Get Analog Value` gives back raw
 unfiltered analog values. The only reason to use :func:`Get Analog Value` is,
 if you need the full resolution of the analog-to-digital converter.

If you want the analog value periodically, it is recommended to use the
:cb:`Analog Value` callback and set the period with
:func:`Set Analog Value Callback Period`.
"""
function get_analog_value(device::BrickletVoltage)
    return device.deviceInternal.get_analog_value()
end

export set_voltage_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Voltage` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Voltage` callback is only triggered if the voltage has changed since
the last triggering.
"""
function set_voltage_callback_period(device::BrickletVoltage, period)
    device.deviceInternal.set_voltage_callback_period(period)
end

export get_voltage_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Voltage Callback Period`.
"""
function get_voltage_callback_period(device::BrickletVoltage)
    return device.deviceInternal.get_voltage_callback_period()
end

export set_analog_value_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Analog Value` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Analog Value` callback is only triggered if the analog value has
changed since the last triggering.
"""
function set_analog_value_callback_period(device::BrickletVoltage, period)
    device.deviceInternal.set_analog_value_callback_period(period)
end

export get_analog_value_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Analog Value Callback Period`.
"""
function get_analog_value_callback_period(device::BrickletVoltage)
    return device.deviceInternal.get_analog_value_callback_period()
end

export set_voltage_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Voltage Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the voltage is *outside* the min and max values"
 "'i'",    "Callback is triggered when the voltage is *inside* the min and max values"
 "'<'",    "Callback is triggered when the voltage is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the voltage is greater than the min value (max is ignored)"
"""
function set_voltage_callback_threshold(device::BrickletVoltage, option, min, max)
    device.deviceInternal.set_voltage_callback_threshold(option, min, max)
end

export get_voltage_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Voltage Callback Threshold`.
"""
function get_voltage_callback_threshold(device::BrickletVoltage)
    return device.deviceInternal.get_voltage_callback_threshold()
end

export set_analog_value_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Analog Value Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the analog value is *outside* the min and max values"
 "'i'",    "Callback is triggered when the analog value is *inside* the min and max values"
 "'<'",    "Callback is triggered when the analog value is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the analog value is greater than the min value (max is ignored)"
"""
function set_analog_value_callback_threshold(device::BrickletVoltage, option, min, max)
    device.deviceInternal.set_analog_value_callback_threshold(option, min, max)
end

export get_analog_value_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Analog Value Callback Threshold`.
"""
function get_analog_value_callback_threshold(device::BrickletVoltage)
    return device.deviceInternal.get_analog_value_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callbacks

* :cb:`Voltage Reached`,
* :cb:`Analog Value Reached`

are triggered, if the thresholds

* :func:`Set Voltage Callback Threshold`,
* :func:`Set Analog Value Callback Threshold`

keep being reached.
"""
function set_debounce_period(device::BrickletVoltage, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletVoltage)
    return device.deviceInternal.get_debounce_period()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletVoltage)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletVoltage, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
"""
    asymptoticdistribution(x, wi, mu, sigmas);

Simulate the asymptotic distribution of test statistic for `kstest`.
`nrep` is the number of random values to generate.
`debuginfo` is whether to show the debug information.
When number of components of null distribution is greater than 1, the test statistic has no closed form asymptotic distribution. When the null distribution is just normal, the asymptotic distribution is just `Chisq(2)`.

"""
function asymptoticdistribution(x::RealVector{Float64}, wi::Vector{Float64}, mu::Vector{Float64}, sigmas::Vector{Float64}; nrep::Int=10000, debuginfo::Bool=false)

    n = length(x)
    C = length(wi)
    m = MixtureModel(map((u, v) -> Normal(u, v), mu, sigmas), wi)
    llC = zeros(n, C)
    S_π = zeros(n, C-1)
    S_μσ = zeros(n, 2*C)
    S_λ = zeros(n, 2*C)
    ll = logpdf(m, x)
    for i in 1:n, kcom in 1:C
        llC[i, kcom] = logpdf(m.components[kcom], x[i])
    end

    for kcom in 1:(C-1)
        S_π[:, kcom] = exp(llC[:, kcom] .- ll) .- exp(llC[:, C] .- ll) #(llC[:, kcom] .- llC[:, C]) ./ ll
    end
    for i in 1:n
        for kcom in 1:C
            llC[i, kcom] = exp(log(wi[kcom]) + llC[i, kcom] - ll[i])
            S_μσ[i, 2*kcom-1] = H1(x[i], mu[kcom], sigmas[kcom]) * llC[i, kcom]
            S_μσ[i, 2*kcom] = H2(x[i], mu[kcom], sigmas[kcom]) * llC[i, kcom]
            S_λ[i, 2*kcom-1] = H3(x[i], mu[kcom], sigmas[kcom]) * llC[i, kcom]
            S_λ[i, 2*kcom] = H4(x[i], mu[kcom], sigmas[kcom]) * llC[i, kcom]
        end
    end
    S_η = hcat(S_π, S_μσ)
    debuginfo && println(round(llC[1:5,:], 6))
    debuginfo && println(sum(S_η, 1))
    debuginfo && println(sum(S_λ, 1))
    I_η = S_η'*S_η./n
    I_λη = S_λ'*S_η./n
    I_λ = S_λ'*S_λ./n
    I_all = vcat(hcat(I_η, I_λη'), hcat(I_λη, I_λ))
    if 1/cond(I_all) < eps(Float64)
        D, V = eig(I_all)
        debuginfo && println(D)
        tol2 = maximum(abs(D)) * 1e-14
        D[D.<tol2] = tol2
        I_all = V*diagm(D)*V'
    end
    debuginfo && println(round(cor(S_η), 6))
    debuginfo && println(round(cor(S_λ), 6))
    I_λ_η = I_all[(3*C):(5*C-1), (3*C):(5*C-1)] - I_all[(3*C):(5*C-1), 1:(3*C-1)] * inv(I_all[1:(3*C-1), 1:(3*C-1)]) * I_all[1:(3*C-1),(3*C):(5*C-1)]
    debuginfo && println(round(I_λ_η, 6))
    #I_λ_η=(I_λ_η .+ I_λ_η')./2
    D, V = eig(I_λ_η)
    D[D.<0.] = 0.
    debuginfo && println(D)
    I_λ_η2 = V * diagm(sqrt(D)) * V'
    u = randn(nrep, 2*C) * I_λ_η2
    EM = zeros(nrep, C)
    T = zeros(nrep)
    for kcom in 1:C
        EM[:, kcom] = sum(u[:, (2*kcom-1):(2*kcom)] * inv(I_λ_η[(2*kcom-1):(2*kcom), (2*kcom-1):(2*kcom)]) .* u[:, (2*kcom-1):(2*kcom)], 2)
    end
    for i in 1:nrep
        T[i] = maximum(EM[i, :])
    end
    debuginfo && println(EM[1:10,:])
    T
end

"""
    gmmrepeat(x, C)

Repeat the `gmm` for `ntrials` randomly generated starting values and pick the one with largest penalized likelihood.

 - `wi_init`, `mu_lb`, `mu_ub`, `sigmas_lb`, `sigmas_ub`: specify how to generate the random starting values
 - `taufixed`, `whichtosplit`, `tau`: whether keep the ratio between `wi[whichtosplit]` and `wi[whichtosplit]+wi[whichtosplit+1]` fixed at a constant `tau`.
 - `sn` and `an`: the penalty
 - `debuginfo` and `tol`: whether print the debug information and the convergence critera
 - `pl`: wheter the penalty on `sigmas` be included in the log likelihood in the final two EM steps. Note that the starting value with largest penalized log likelihood is picked, but the penalty term should not be included in the likelihood ratio
 - `ptau`: whether to add the penalty on `tau` be included in likelihood. Better to be `true` since the more `tau` values we try the larger the test statistic

"""
function gmmrepeat(x::RealVector, C::Int,
    wi_init::Vector{Float64},
    mu_init::Vector{Float64},
    sigmas_init::Vector{Float64};
    ntrials::Int=25,
    taufixed::Bool=false, whichtosplit::Int=1, tau::Real=0.5,
   sn::Vector{Float64}=std(x).*ones(C), an::Real=1/length(x), debuginfo::Bool=false, tol::Real=.001, pl::Bool=false, ptau::Bool=false)

    n = length(x)
    tau = min(tau, 1-tau)
    mu_lb = minimum(x) .* ones(C)
    mu_ub = maximum(x) .* ones(C)
    if taufixed && C>2
        for kcom in 1:(C-1)
            tmp = (mu_init[kcom+1]*(wi_init[kcom+1])^.25 +mu_init[kcom]*(wi_init[kcom])^.25) / ((wi_init[kcom+1])^.25+(wi_init[kcom])^.25)
            mu_lb[kcom+1] = tmp
            mu_ub[kcom] = tmp
        end
        mu_lb[whichtosplit+1]=mu_lb[whichtosplit]
        mu_ub[whichtosplit]=mu_ub[whichtosplit+1]
    end
    debuginfo && println(mu_lb, mu_ub)
    sigmas_lb = 0.25 .* sigmas_init
    sigmas_ub = 2 .* sigmas_init

    if taufixed
        tmp = wi_init[whichtosplit] + wi_init[whichtosplit+1]
        wi_init[whichtosplit] = tmp*tau
        wi_init[whichtosplit+1] = tmp*(1-tau)
        mu_init[whichtosplit] -= 1e-3
    end

    wi = repmat(wi_init, 1, 4*ntrials)
    mu = repmat(mu_init, 1, 4*ntrials)
    sigmas = repmat(sigmas_init, 1, 4*ntrials)
    ml = -Inf .* ones(4*ntrials)
    for i in 1:4*ntrials
        #set the first initial value as mu_init, sigmas_init
        if i >= 2
            mu[:, i] = rand(C) .* (mu_ub .- mu_lb) .+ mu_lb
            sigmas[:, i] = rand(C) .* (sigmas_ub .- sigmas_lb) .+ sigmas_lb
        end

        wi[:, i], mu[:, i], sigmas[:, i], ml[i] =
             gmm(x, C, wi[:, i], mu[:, i], sigmas[:, i],
             taufixed=taufixed, whichtosplit=whichtosplit, tau=tau,
             mu_lb=mu_lb, mu_ub=mu_ub,
             maxiteration=100, sn=sn, an=an,
             tol=tol, pl=true, ptau=false)
    end

    mlperm = sortperm(ml)
    for j in 1:ntrials
        i = mlperm[4*ntrials+1 - j] # start from largest ml
        wi[:, i], mu[:, i], sigmas[:, i], ml[i] =
            gmm(x, C, wi[:, i], mu[:, i], sigmas[:, i],
            taufixed=taufixed, whichtosplit=whichtosplit, tau=tau,
            mu_lb=mu_lb, mu_ub=mu_ub,
            sn=sn, an=an,
            tol=tol, pl=true, ptau=false)
    end

    mlmax, imax = findmax(ml[mlperm[(3*ntrials+1):4*ntrials]])
    imax = mlperm[3*ntrials+imax]

    re=gmm(x, C, wi[:, imax], mu[:, imax], sigmas[:, imax],
         maxiteration=2, an=an, sn=sn, tol=0., pl=pl, ptau=ptau, whichtosplit=whichtosplit)
    debuginfo && println("Trial:", re)
    return(re)
end

"""
    kstest(x, C0)

Do the EM test under null Hypothesis of `C0` components.
If rejected, then it suggest the true number of components is greater than `C0`.
Optional arguments for `kstest`

 - `vtau`: the finite set of `tau` value, default to try 0.5 only
 - `ntrials`: the number of initial values to try, default to be 25
 - `debuginfo`: whether show the debug information

"""
function kstest(x::RealVector{Float64}, C0::Int; vtau::Vector{Float64}=[.5;],
    ntrials::Int=25, debuginfo::Bool=false, tol::Real=0.001,ptau::Bool=true)

    C1 = C0+1
    n = length(x)

    wi_init, mu_init, sigmas_init, ml_C0 = gmm(x, C0)
    wi_init, mu_init, sigmas_init, ml_C0 = gmmrepeat(x, C0, wi_init, mu_init, sigmas_init, pl=false, ptau=false, an=1/n)
    debuginfo && println(wi_init, mu_init, sigmas_init, ml_C0)
    if C0 > 1
        trand=GaussianMixtureTest.asymptoticdistribution(x, wi_init, mu_init, sigmas_init)
    end

    if debuginfo
        println("ml_C0=", ml_C0)
    end
    minx = minimum(x)
    maxx = maximum(x)
    or = sortperm(mu_init)
    wi0 = wi_init[or]
    mu0 = mu_init[or]
    sigmas0 = sigmas_init[or]
    an = decidepenalty(wi0, mu0, sigmas0, n)
    lr = 0.0
    lrv = zeros(length(vtau), C0)
    for whichtosplit in 1:C0, i in 1:length(vtau)
        ind = [1:whichtosplit, whichtosplit:C0;]

        wi_C1 = wi0[ind]
        wi_C1[whichtosplit] = wi_C1[whichtosplit]*vtau[i]
        wi_C1[whichtosplit+1] = wi_C1[whichtosplit+1]*(1-vtau[i])

        lrv[i, whichtosplit] = gmmrepeat(x, C1,
         wi_C1, mu0[ind], sigmas0[ind],
         ntrials=ntrials,
         taufixed=true, whichtosplit=whichtosplit, tau=vtau[i],
         sn=sigmas0[ind], an=an,
         debuginfo=debuginfo, tol=tol, pl=false, ptau=ptau)[4]
       if debuginfo
           println(whichtosplit, " ", vtau[i], "->",
           lrv[i, whichtosplit])
       end
   end
   lr = maximum(lrv)
   if debuginfo
       println("lr=", lr)
   end
   Tvalue = 2*(lr - ml_C0)
   if C0 == 1
       pvalue = 1 - cdf(Chisq(2), Tvalue)
   else
       pvalue = mean(trand .> Tvalue)
   end
   return(Tvalue, pvalue)
end


text:
## Method 1

function get_uncondensed_state(ind1::T, ind2::T) where T <: Tuple
    i,j = ind1
    k,l = ind2

    δ = 1 * (i == j == k == l) +
    2 * ((i == j == k) & (i != l)) +
    3 * ((i == j == l) & (i != k)) +
    4 * ((i == k == l) & (i != j)) +
    5 * ((j == k == l) & (i != j)) +
    6 * ((i == j) & (k == l) & (i != k)) +
    7 * ((i == j) & (k != l) & (i != k) & (i != l)) +
    8 * ((k == l) & (i != j) & (i != k) & (j != k)) +
    9 * ((i == k) & (j == l) & (i != j) & (k != l)) +
    10 * ((i == k) & (j != l) & (i != j) & (k != l)) +
    11 * ((j == l) & (i != k) & (i != j) & (k != l)) +
    12 * ((i == l) & (j == k) & (i != j) & (k != l)) +
    13 * ((i == l) & (j != k) & (i != j) & (k != l)) +
    14 * ((i != l) & (j == k) & (i != j) & (k != l)) +
    15 * ((i != j) & (i != k) & (i != l) & (j != k) & (j != l) & (k != l))

    return δ
end

s1(p) = [p[1], p[1]^2, p[1]^2, p[1]^3, p[1]^2, p[1]^3, p[1]^2, p[1]^3, p[1]^4]

s2(p) = [0.0, 0.0, prod(p[[1,4]]), 2.0 * prod(p[[1,4]]) * p[1], 0.0, 0.0, 0.0, prod(p[[1,4]]) * p[1], 2.0 * prod(p[[1,4]]) * p[1]^2] #D3/S2

s3(p) = [0.0, 0.0, prod(p[[1,3]]), 2.0 * prod(p[[1,3]]) * p[1], 0.0, 0.0, 0.0, prod(p[[1,3]]) * p[1], 2.0 * prod(p[[1,3]]) * p[1]^2] #D3/S3

s4(p) = [0.0, 0.0, 0.0, 0.0, prod(p[[1,2]]), 2.0 * prod(p[[1,2]]) * p[1], 0.0, prod(p[[1,2]]) *p[1], 2.0 * prod(p[[1,2]]) * p[1]^2] #D5/S4

s5(p) = [0.0, 0.0, 0.0, 0.0, prod(p[[1,2]]), 2.0 * prod(p[[1,2]]) * p[2], 0.0, prod(p[[1,2]]) *p[2], 2.0 * prod(p[[1,2]]) * p[2]^2] #D5/S5

s6(p) = [0.0, prod(p[[1,3]]), 0.0, prod(p[[1,3]]) * p[3], 0.0, prod(p[[1,3]]) * p[1], 0.0, 0.0, prod(p[[1,3]])^2] #D2/S6

s7(p) = [0.0, 0.0, 0.0, 2.0 * prod(p[[1,3,4]]), 0.0, 0.0, 0.0, 0.0, 2.0 * prod(p[[1,3,4]]) * p[1]] #D4/S7

s8(p) = [0.0, 0.0, 0.0, 0.0, 0.0, 2.0 * prod(p[[1,2,3]]), 0.0, 0.0, 2.0 * prod(p[[1,2,3]]) * p[3]] #D6/S8

s9(p) =  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0 * prod(p[[1,2]]), prod(p[[1,2]]) * sum(p[[1,2]]), 4.0 * prod(p[[1,2]])^2] #D7/S9

s10(p) = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, prod(p[[1,2,4]]), 4.0 * prod(p[[1,2,4]]) * p[1]] #D8/S10

s11(p) = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, prod(p[[1,2,3]]), 4.0 * prod(p[[1,2,3]]) * p[2]] #D8/S11

s12(p) = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0 * prod(p[[1,2]]), prod(p[[1,2]]) * sum(p[[1,2]]), 4.0 * prod(p[[1,2]])^2] #D7/S12

s13(p) = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, prod(p[[1,2,3]]), 4.0 * prod(p[[1,2,3]]) * p[1]] #D8/S13

s14(p) = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, prod(p[[1,2,4]]), 4.0 * prod(p[[1,2,4]]) * p[2]] #D8/S14

s15(p) = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0 * prod(p)] #D9/S15

δ = (s1, s2, s3, s4, s5, s6, s7, s8, s9, s10, s11, s12, s13, s14, s15)

function probability_states(ind1::T, ind2::T, alleles::Dict) where T <: Tuple
    i,j = ind1
    k,l = ind2
    p = [alleles[i], alleles[j], alleles[k], alleles[l]]

    return δ[get_uncondensed_state(ind1, ind2)](p)
end

s1(p) = (p[1], p[1]^2, p[1]^2, p[1]^3, p[1]^2, p[1]^3, p[1]^2, p[1]^3, p[1]^4)

s2(p) = (0.0, 0.0, prod(p[[1,4]]), 2.0 * prod(p[[1,4]]) * p[1], 0.0, 0.0, 0.0, prod(p[[1,4]]) * p[1], 2.0 * prod(p[[1,4]]) * p[1]^2) #D3/S2

s3(p) = (0.0, 0.0, prod(p[[1,3]]), 2.0 * prod(p[[1,3]]) * p[1], 0.0, 0.0, 0.0, prod(p[[1,3]]) * p[1], 2.0 * prod(p[[1,3]]) * p[1]^2) #D3/S3

s4(p) = (0.0, 0.0, 0.0, 0.0, prod(p[[1,2]]), 2.0 * prod(p[[1,2]]) * p[1], 0.0, prod(p[[1,2]]) *p[1], 2.0 * prod(p[[1,2]]) * p[1]^2) #D5/S4

s5(p) = (0.0, 0.0, 0.0, 0.0, prod(p[[1,2]]), 2.0 * prod(p[[1,2]]) * p[2], 0.0, prod(p[[1,2]]) *p[2], 2.0 * prod(p[[1,2]]) * p[2]^2) #D5/S5

s6(p) = (0.0, prod(p[[1,3]]), 0.0, prod(p[[1,3]]) * p[3], 0.0, prod(p[[1,3]]) * p[1], 0.0, 0.0, prod(p[[1,3]])^2) #D2/S6

s7(p) = (0.0, 0.0, 0.0, 2.0 * prod(p[[1,3,4]]), 0.0, 0.0, 0.0, 0.0, 2.0 * prod(p[[1,3,4]]) * p[1]) #D4/S7

s8(p) = (0.0, 0.0, 0.0, 0.0, 0.0, 2.0 * prod(p[[1,2,3]]), 0.0, 0.0, 2.0 * prod(p[[1,2,3]]) * p[3]) #D6/S8

s9(p) =  (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0 * prod(p[[1,2]]), prod(p[[1,2]]) * sum(p[[1,2]]), 4.0 * prod(p[[1,2]])^2) #D7/S9

s10(p) = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, prod(p[[1,2,4]]), 4.0 * prod(p[[1,2,4]]) * p[1]) #D8/S10

s11(p) = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, prod(p[[1,2,3]]), 4.0 * prod(p[[1,2,3]]) * p[2]) #D8/S11

s12(p) = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0 * prod(p[[1,2]]), prod(p[[1,2]]) * sum(p[[1,2]]), 4.0 * prod(p[[1,2]])^2) #D7/S12

s13(p) = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, prod(p[[1,2,3]]), 4.0 * prod(p[[1,2,3]]) * p[1]) #D8/S13

s14(p) = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, prod(p[[1,2,4]]), 4.0 * prod(p[[1,2,4]]) * p[2]) #D8/S14

s15(p) = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0 * prod(p)) #D9/S15

δ_t = (s1, s2, s3, s4, s5, s6, s7, s8, s9, s10, s11, s12, s13, s14, s15)

function probability_states_tuple(ind1::T, ind2::T, alleles::Dict) where T <: Tuple
    i,j = ind1
    k,l = ind2
    p = [alleles[i], alleles[j], alleles[k], alleles[l]]

    return δ_t[get_uncondensed_state(ind1, ind2)](p)
end

## Method 2
S1 = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4]
S2 = [1 2; 3 4]
S3 = [1 2; 1 3; 2 3]
S4 = [1 2]
S5 = [1 3; 1 4; 3 4]
S6 = [3 4]
S7 = [1 3; 2 4]
S8 = [1 3]
S9 = [1 3][[1 3][:,1] .!= 1, :]
ibd_states = [S1, S2, S3, S4, S5, S6, S7, S8, S9]

function get_jacquard_state(ind1::T, ind2::T, ibd_states = ibd_states) where T <: Tuple
    the_shared = [intersect(ind1, ind2)...,0]
    the_lone = [setdiff(ind1, ind2)..., setdiff(ind2, ind1)...,0]

    ind1_resort = [ind1[(sum(ind1 .∈ the_shared',dims=2) .>= 1)[:,1]]..., ind1[(sum(ind1 .∈ the_lone',dims=2) .>= 1)[:,1]]...]
    ind2_resort = [ind2[(sum(ind2 .∈ the_shared',dims=2) .>= 1)[:,1]]..., ind2[(sum(ind2 .∈ the_lone',dims=2) .>= 1)[:,1]]...]

    adj_mat = Array{Int8}(undef, 6, 2)
    idx = 0
    for i in 1:3, j in i+1:4
        idx += 1
        adj_mat[idx,1:2] = [i,j] .* (1.0 * ([ind1_resort..., ind2_resort...][i] == [ind1_resort..., ind2_resort...][j]))
    end
    adj_mat = adj_mat[adj_mat[:,1] .> 0,:]

    jaq_state = 0
    for i in 1:9
        jaq_state += i * (adj_mat == ibd_states[i])
    end
    return jaq_state
end

## Single type?
struct JacquardPair
    genotype1::Genotype
    genotype2::Genotype
    state::Int
end

function probability_state(genos::JacquardPair, alleles::Dict)::Vector{Float64}
    state = genos.state
    if state == 1
        p = alleles[genos.genotype1[1]]
        [p, p^2, p^2, p^3, p^2, p^3, p^2, p^3, p^4]
    elseif state == 2
        p = (alleles[genos.genotype1[1]], alleles[genos.genotype2[1]])
        [0.0, prod(p), 0.0, prod(p) * p[2], 0.0, prod(p) * p[1], 0.0, 0.0, prod(p)^2]

    elseif state == 3
        p = (alleles[genos.genotype1[1]], alleles[genos.genotype2[2]])
        [0.0, 0.0, prod(p), 2.0 * prod(p) * p[1], 0.0, 0.0, 0.0, prod(p) * p[1], 2.0 * prod(p) * p[1]^2]

    elseif state == 4
        p = (alleles[genos.genotype1[1]], alleles[genos.genotype2[1]], alleles[genos.genotype2[2]])
        [0.0, 0.0, 0.0, 2.0 * prod(p), 0.0, 0.0, 0.0, 0.0, 2.0 * prod(p) * p[1]]

    elseif state == 5
        p = (alleles[genos.genotype1[1]], alleles[genos.genotype1[2]])
        [0.0, 0.0, 0.0, 0.0, prod(p), 2.0 * prod(p) * p[1], 0.0, prod(p) *p[1], 2.0 * prod(p) * p[1]^2]

    elseif state == 6
        p = (alleles[genos.genotype2[1]], alleles[genos.genotype1[1]], alleles[genos.genotype1[2]])
        [0.0, 0.0, 0.0, 0.0, 0.0, 2.0 * prod(p), 0.0, 0.0, 2.0 * prod(p) * p[1]]

    elseif state == 7
        p = (alleles[genos.genotype1[1]], alleles[genos.genotype1[2]])
        [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0 * prod(p), prod(p) * sum(p), 4.0 * prod(p)^2]

    elseif state == 8
        p = (alleles[genos.genotype1[1]], alleles[genos.genotype1[2]], alleles[genos.genotype2[2]])
        [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, prod(p), 4.0 * prod(p) * p[1]]

    elseif state == 9
        p = (alleles[genos.genotype1[1]], alleles[genos.genotype1[2]], alleles[genos.genotype2[1]], alleles[genos.genotype2[2]])
        [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0 * prod(p)]
    else
        [-9.0, -9.0, -9.0, -9.0, -9.0, -9.0, -9.0, -9.0, -9.0]
    end
end

## multiple types?
struct JacquardZero
    genotype1::Genotype
    genotype2::Genotype
end

struct JacquardOne
    genotype1::Genotype
    genotype2::Genotype
end

struct JacquardTwo
    genotype1::Genotype
    genotype2::Genotype
end

struct JacquardThree
    genotype1::Genotype
    genotype2::Genotype
end

struct JacquardFour
    genotype1::Genotype
    genotype2::Genotype
end

struct JacquardFive
    genotype1::Genotype
    genotype2::Genotype
end

struct JacquardSix
    genotype1::Genotype
    genotype2::Genotype
end

struct JacquardSeven
    genotype1::Genotype
    genotype2::Genotype
end

struct JacquardEight
    genotype1::Genotype
    genotype2::Genotype
end

struct JacquardNine
    genotype1::Genotype
    genotype2::Genotype
end

function probability_state(genos::JacquardZero, alleles::Dict)::Vector{Float64}
    [-9.0, -9.0, -9.0, -9.0, -9.0, -9.0, -9.0, -9.0, -9.0]
end

function probability_state(genos::JacquardOne, alleles::Dict)::Vector{Float64}
    p = alleles[genos.genotype1[1]]
    [p, p^2, p^2, p^3, p^2, p^3, p^2, p^3, p^4]
end

function probability_state(genos::JacquardTwo, alleles::Dict)::Vector{Float64}
    p = (alleles[genos.genotype1[1]], alleles[genos.genotype2[1]])
    [0.0, prod(p), 0.0, prod(p) * p[2], 0.0, prod(p) * p[1], 0.0, 0.0, prod(p) * prod(p)]
end

function probability_state(genos::JacquardThree, alleles::Dict)::Vector{Float64}
    p = (alleles[genos.genotype1[1]], alleles[genos.genotype2[1]])
    [0.0, 0.0, prod(p), 2.0 * prod(p) * p[1], 0.0, 0.0, 0.0, prod(p) * p[1], 2.0 * prod(p) * p[1]^2]
end

function probability_state(genos::JacquardFour, alleles::Dict)::Vector{Float64}
    p = (alleles[genos.genotype1[1]], alleles[genos.genotype2[1]], alleles[genos.genotype2[2]])
    [0.0, 0.0, 0.0, 2.0 * prod(p), 0.0, 0.0, 0.0, 0.0, 2.0 * prod(p) * p[1]]
end

function probability_state(genos::JacquardFive, alleles::Dict)::Vector{Float64}
    p = (alleles[genos.genotype1[2]], alleles[genos.genotype1[1]])
    [0.0, 0.0, 0.0, 0.0, prod(p), 2.0 * prod(p) * p[1], 0.0, prod(p) *p[1], 2.0 * prod(p) * p[1]^2]
end

function probability_state(genos::JacquardSix, alleles::Dict)::Vector{Float64}
    p = (alleles[genos.genotype2[1]], alleles[genos.genotype1[1]], alleles[genos.genotype1[2]])
    [0.0, 0.0, 0.0, 0.0, 0.0, 2.0 * prod(p), 0.0, 0.0, 2.0 * prod(p) * p[1]]
end

function probability_state(genos::JacquardSeven, alleles::Dict)::Vector{Float64}
    p = (alleles[genos.genotype1[1]], alleles[genos.genotype1[2]])
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0 * prod(p), prod(p) * sum(p), 4 * prod(p) * prod(p)]
end

function probability_state(genos::JacquardEight, alleles::Dict)::Vector{Float64}
    p = (alleles[genos.genotype1[1]], alleles[genos.genotype1[2]], alleles[genos.genotype2[2]])
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, prod(p), 4 * prod(p) * p[1]]
end

function probability_state(genos::JacquardNine, alleles::Dict)::Vector{Float64}
    p = (alleles[genos.genotype1[1]], alleles[genos.genotype1[2]], alleles[genos.genotype2[1]], alleles[genos.genotype2[2]])
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4 * prod(p)]
end


text:
module DiscreteControl
using Interpolations
using Optim, ForwardDiff
using Distributions

export DynamicSystem1D, DynamicSystemTrajectory1D, OfflineSystemControl1D,
    MPCSystem1D, OLFCSystem1D
export DynamicSystem
export initializestate!, step!, simulatetrajectories
export solvebellman!, onlinedecision, solveolfc!

include("types.jl")
include("bellman.jl")
include("mpc.jl")
include("olfc.jl")
include("api.jl")
include("utils.jl")

end # module


text:
using LinearAlgebra

abstract type AbstractForecasting end

export StateSpaceModel

"""

Space-State model is defined through the following equations

Let ``{x_0, x_1, \\dots, x_K }`` denote the discrete-time state process,
which is typically unknown to the observer. Let ``{y_1,y_2,... ,y_K}``
represent the observation process, which describes the measurements of the
system that are available to the observer.

These state and observation processes are related following a nonlinear
state-space model: The first equation describes how the system state ``x_k``
evolves through a nonlinear dynamical model ``f_{k−1,k}`` between successive
time steps ``t_{k−1}`` and ``t_k``.
``η = {η1,η2,... ,ηK}`` is the model noise process, which accounts for
the imperfections of the model; it is assumed to be independent and
identically distributed. At each time step, ``η_k`` is assumed
to be Gaussian with zero mean and covariance matrix Q. The second equation
models how an observation ``y_k`` at time ``t_k`` is obtained from the state
``x_k`` through a linear operator ``H``. ``ε = {ε1, ε2, . . . , εk}`` is an
independent and identically distributed process representing the observation
errors; ``ε_k`` is assumed
Gaussian with zero mean and covariance matrix ``R``. Finally, the
processes ``η`` and ``ε`` are assumed to be jointly independent and independent
of the initial (background) state x0, which is assumed to be Gaussian
with mean xb and covariance B.

```math
\\left\\{
\\begin{array}{l}
X_t = m(X_{t-1}) + \\eta_t, \\\\
Y_t = H(X_t) + \\varepsilon_t,
\\end{array}
\\right.
```

- X : hidden variables
- Y : observed variables

- `dt_integration`is the numerical time step used to solve the ODE.
- `dt_states` is the number of `dt_integration` between ``X_{t-1}`` and ``X_t``.
- `dt_obs` is the number of `dt_integration` between ``Y_{t-1}`` and ``Y_t``.

"""
struct StateSpaceModel <: AbstractForecasting

    model::Function
    dt_integration::Float64
    dt_states::Int64
    dt_obs::Int64
    params::Vector{Float64}
    var_obs::Vector{Int64}
    nb_loop_train::Int64
    nb_loop_test::Int64
    sigma2_catalog::Float64
    sigma2_obs::Float64

end


text:
# Demonstrate/test morphological functions

using ImageProjectiveGeometry
using PyPlot
using PerceptualColourMaps

img = imread("lena.tif");
figure(10), imshow(img);
set_cmap(ColorMap("gray"))
title("Input image")

img2 = rand(512,512) .> 0.999;
figure(11), imshow(img2);
title("Dots input image")

imgd = imerode(img, "RECT", [9, 15]);
figure(12), imshow(imgd);
title("Rectangular dilate")
keypause()

imgd = imerode(img, "RECT", 21);
figure(12), imshow(imgd);
title("Rectangular erode")
keypause()

imgd = imdilate(img, "OCT", 9);
figure(12), imshow(imgd);
title("Octagonal dilate")
keypause()

imgd = imerode(img, "OCT", 9);
figure(12), imshow(imgd);
title("Octagonal erode")
keypause()

img2d = imdilate(img2, "RECT", [19, 15]);
figure(12), imshow(img2d);
title("Dots: Rectangular dilate")
keypause()

img2d = imdilate(img2, "OCT", 9);
figure(12), imshow(img2d);
title("Dots: Octagonal dilate")
keypause()

se = circularstruct(7);
img2d = imdilate(img2, se);
figure(12), imshow(img2d);
title("Dots: Circular dilate")


text:
"""
This module provides high-level functionality to generate, manage and refine
meshes based on constrained Delaunay triangulations.
"""

module DelaunayMeshes

import DiffEqPDEBase
import GeometricalPredicates

import Base.convert
import Base.push!
import Base.next
import Base.done
import Base.start

include("QuadEdge.jl")
include("Triangulation.jl")
include("DelaunayMesh.jl")
include("PostProcess.jl")

export getvoronoivertices, setboundingbox
export scalepoints, unscalepoints, addconstraint!
export getdelaunaycoordinates
end # module


text:
function r_init!(t::RelaxInterval, g::DAT, b::IntervalCache{T}) where T<:Real
    z = set(b, 1) ∩ g.sink_bnd
    b[1] = z
    return !isempty(z)
end

function rprop!(t::RelaxInterval, v::Variable, g::DAT, c::IntervalCache{T}, k) where T<:Real
    z = z ∩ set(t, c, k)
    c[k] = z
    return !isempty(z)
end

function rprop!(t::RelaxInterval, v::Subexpression, g::DAT, c::IntervalCache{T}, k) where T<:Real
    store_subexpression!(c, set(t, c, k), first_index(g, k))
    return true
end

# needed for O(n) reverse interval propagation of +
# returns q for x = q + y
function hukuhara_diff(x::Interval{T}, y::Interval{T}) where T<:Real
    isempty(x) && return x
    isempty(y) && return y
    l = sub_round(x.lo, y.lo, RoundDown)
    u = sub_round(x.hi, y.hi, RoundUp)
    Interval{T}(l, u)
end

"""
$(FUNCTIONNAME)

Updates storage tapes with reverse evalution of node representing `n = +(x,y,z...)` which updates x, y, z and so on.
"""
function rprop!(t::RelaxInterval, v::Val{PLUS}, g::DAT, c::IntervalCache{T}, k::Int) where T<:Real
    tsum = sum(j -> set(t, c, j), children(g, k))
    for j in children(g, k)
        tmsum = hukuhara_diff(tsum, set(c, j))
        _, w, _ = IntervalContractors.plus_rev(set(t, c, k), set(t, c, i), tmsum)
        isempty(w) && return false
        c[i] = v
    end
    return true
end

# needed for close to O(n) reverse interval propagation of *
# returns q for x = q*y
function hukuhara_div(x::Interval{T}, y::Interval{T}) where T<:Real
    isempty(x) && return x
    isempty(y) && return y
    if y.lo >= zero(T)
        if x.lo >= zero(T)
            l = div_round(x.lo, y.lo, RoundDown)
            u = div_round(x.hi, y.hi, RoundUp)
        elseif x.hi <= zero(T)
            l = div_round(x.lo, y.hi, RoundDown)
            u = div_round(x.hi, y.lo, RoundUp)
        else
            l = div_round(x.lo, y.hi, RoundDown)
            u = div_round(x.hi, y.hi, RoundUp)
        end
        return Interval(l, u), true
    elseif y.hi <= zero(T)
        if x.lo >= zero(T)
            l = div_round(x.hi, y.lo, RoundDown)
            u = div_round(x.lo, y.hi, RoundUp)
        elseif x.hi <= zero(T)
            l = div_round(x.hi, y.hi, RoundDown)
            u = div_round(x.lo, y.lo, RoundUp)
        else
            l = div_round(x.hi, y.lo, RoundDown)
            u = div_round(x.lo, y.lo, RoundUp)
        end
        return Interval(l, u), true
    else
        if x.lo > zero(T)
            l = div_round(x.hi, y.lo, RoundDown)
            u = div_round(x.hi, y.hi, RoundUp)
            return Interval(l, u), true
        elseif x.hi < zero(T)
            l = div_round(x.lo, y.hi, RoundDown)
            u = div_round(x.lo, y.lo, RoundUp)
            return Interval(l, u), true
        end

    end
    empty(Interval{T}), false
end

"""
$(FUNCTIONNAME)

Updates storage tapes with reverse evalution of node representing `n = *(x,y,z...)` which updates x, y, z and so on.
"""
function rprop!(t::RelaxInterval, v::Val{MULT}, g::DAT, c::IntervalCache{T}, k::Int) where T<:Real
    tmul = sum(j -> set(t, c, j), children(g, k))
    for j in children(g, k)
        tmulm, hdiv = hukuhara_div(tmul, set(t, b, j))
        if !hdiv
            tmulmf = one(Interval{T})
            for i in children(g, k)
                if i != j
                    tmulmf *= set(t, b, i)
                end
            end
            tmulm = tmulmf
        end
        _, w, _ = IntervalContractors.mul_rev(set(t, b, k), set(t, b, c), tmulm)
        isempty(w) && return false
        c[i] = w
    end
    return true
end

for (f, fc, F) in ((-, MINUS, IntervalContractors.minus_rev),
                   (^, POW, IntervalContractors.power_rev),
                   (/, DIV, IntervalContractors.div_rev))
    @eval function rprop!(t::RelaxInterval, v::Val{$fc}, g::DAT, b::IntervalCache{T}, k) where T<:Real
        x = child(g, 1, k)
        y = child(g, 2, k)
        z, u, v = ($F)(set(t, b, k), set(t, b, x), set(t, b, y))
        isempty(u) && return false
        isempty(v) && return false
        b[x] = u
        b[y] = v
        return true
    end
end

rprop!(t::RelaxInterval, v::Val{USER}, g::DAT, b::IntervalCache, k::Int) = true
rprop!(t::RelaxInterval, v::Val{USERN}, g::DAT, b::IntervalCache, k::Int) = true

for ft in UNIVARIATE_ATOM_TYPES
    f = UNIVARIATE_ATOM_DICT[ft]
    (f == :user || f == :+ || f == :-) && continue
    @eval rprop!(t::RelaxInterval, v::Val{$ft}, g::DAT, b::IntervalCache, k::Int) = true
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: MatthewStuber/EAGO


text:
using Documenter, AWSBatch

makedocs(;
    modules=[AWSBatch],
    format=Documenter.HTML(
        prettyurls=get(ENV, "CI", nothing) == "true",
        assets=[
            "assets/invenia.css",
        ],
    ),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/JuliaCloud/AWSBatch.jl/blob/{commit}{path}#L{line}",
    sitename="AWSBatch.jl",
    authors="Invenia Technical Computing",
    strict = true,
    checkdocs = :none,
)

deploydocs(;
    repo="github.com/JuliaCloud/AWSBatch.jl",
    push_preview=true,
    devbranch = "main"
)