text:
"""
    Gate(x, y, logic)

Constructs a logic gate at position (x, y) and returns type Block. Acceptable
inputs for logic are (:AND, :NAND, :OR, :NOR, :XOR, :XNOR and :NOT)
"""
function Gate(
        x::Real, y::Real, logic::Symbol, inputs=[:in1, :in2],
        output=:out, scale_x=1.0, scale_y=0.8
    )
    if      logic == :AND   return and_gate(x, y, inputs, output, scale_x, scale_y)
    elseif  logic == :NAND  return nand_gate(x, y, inputs, output, scale_x, scale_y)
    elseif  logic == :OR    return or_gate(x, y, inputs, output, scale_x, scale_y)
    elseif  logic == :NOR   return nor_gate(x, y, inputs, output, scale_x, scale_y)
    elseif  logic == :XOR   return xor_gate(x, y, inputs, output, scale_x, scale_y)
    elseif  logic == :XNOR  return xnor_gate(x, y, inputs, output, scale_x, scale_y)
    elseif  logic == :NOT   return not_gate(x, y, [:in], output, scale_x, scale_y)
    else
        println("Failed to recognize Logic")
        return nothing
    end
end

# Basic shape of AND
function and_base()
    compose(
        context(),
        line([(0.001, 1), (-1, 1), (-1, -1), (0.001, -1)]),
        curve((0, 1), (4/3, 1), (4/3, -1), (0, -1))
    )
end

# Basic shape of OR
function or_base()
    compose(
        context(),
        curve((-1, -1), (-0.5, -0.5), (-0.5, 0.5), (-1, 1)),
        curve((-1, 1), (0, 1), (0.5, 1), (1, 0)),
        curve((-1, -1), (0, -1), (0.5, -1), (1, 0))
    )
end

# Basic shape of XOR
function xor_base()
    compose(
        context(),
        curve((-1, -1), (-0.5, -0.5), (-0.5, 0.5), (-1, 1)),
        curve((-1.2, -1), (-0.7, -0.5), (-0.7, 0.5), (-1.2, 1)),
        curve((-1, 1), (0, 1), (0.5, 1), (1, 0)),
        curve((-1, -1), (0, -1), (0.5, -1), (1, 0))
    )
end

function and_gate(x, y, inputs, output=:out, sx=1.0, sy=1.0)
    p = (x, y)
    in_pos = partition(
        p + 0.5 * (-1sx, -1sy),
        p + 0.5 * (-1sx, 1sy),
        length(inputs)
    )
    connections = Dict{Symbol, XYTuple}(
        map((label, pos) -> Pair(label, pos), inputs, in_pos)...,
        :out => p + 0.5 * (1sx, 0sy)
    )
    composition = compose(
        context(
            x-0.55sx, y-0.55sy, 1.1sx, 1.1sy,
            units = UnitBox(-1.1, -1.1, 2.2, 2.2)
        ),
        and_base()
    )
    return Block(composition, connections)
end

function nand_gate(x, y, inputs, output=:out, sx=1.0, sy=1.0)
    p = (x, y)
    in_pos = partition(
        p + 0.5 * (-1sx, -1sy),
        p + 0.5 * (-1sx, 1sy),
        length(inputs)
    )
    connections = Dict{Symbol, XYTuple}(
        map((label, pos) -> Pair(label, pos), inputs, in_pos)...,
        :out => p + 0.5 * (1.4sx, 0.0sy)
    )
    composition = compose(
        context(
            x-0.55sx, y-0.55sy, 1.3sx, 1.1sy,
            units = UnitBox(-1.1, -1.1, 2.6, 2.2)
        ),
        and_base(),
        (context(), circle(1.2sx, 0, 0.2sx), fill(nothing))
    )
    return Block(composition, connections)
end

function or_gate(x, y, inputs, output=:out, sx=1.0, sy=1.0)
    @assert (length(inputs) == 2) "Only two inputs are implemented for OR."
    p = (x, y)
    connections = Dict{Symbol, XYTuple}(
        inputs[1] => p + 0.5 * (-0.7025sx, -0.5sy),
        inputs[2] => p + 0.5 * (-0.7025sx, 0.5sy),
        :out      => p + 0.5 * (1sx, 0sy)
    )
    composition = compose(
        context(
            x-0.55sx, y-0.55sy, 1.1sx, 1.1sy,
            units = UnitBox(-1.1, -1.1, 2.2, 2.2)
        ),
        or_base()
    )
    return Block(composition, connections)
end

function nor_gate(x, y, inputs, output=:out, sx=1.0, sy=1.0)
    @assert (length(inputs) == 2) "Only two inputs are implemented for NOR."
    p = (x, y)
    connections = Dict{Symbol, XYTuple}(
        inputs[1] => p + 0.5 * (-0.7025sx, -0.5sy),
        inputs[2] => p + 0.5 * (-0.7025sx, 0.5sy),
        :out      => p + 0.5 * (1.4sx, 0.0sy)
    )
    composition = compose(
        context(
            x-0.55sx, y-0.55sy, 1.3sx, 1.1sy,
            units = UnitBox(-1.1, -1.1, 2.6, 2.2)
        ),
        or_base(),
        (context(), circle(1.2sx, 0, 0.2sx), fill(nothing))
    )
    return Block(composition, connections)
end

function xor_gate(x, y, inputs, output=:out, sx=1.0, sy=1.0)
    @assert (length(inputs) == 2) "Only two inputs are implemented for OR."
    p = (x, y)
    connections = Dict{Symbol, XYTuple}(
        inputs[1] => p + 0.5 * (-0.7025sx, -0.5sy),
        inputs[2] => p + 0.5 * (-0.7025sx, 0.5sy),
        :out      => p + 0.5 * (1sx, 0sy)
    )
    composition = compose(
        context(
            x-0.6sx, y-0.55sy, 1.15sx, 1.1sy,
            units = UnitBox(-1.2, -1.1, 2.3, 2.2)
        ),
        xor_base()
    )
    return Block(composition, connections)
end

function xnor_gate(x, y, inputs, output=:out, sx=1.0, sy=1.0)
    @assert (length(inputs) == 2) "Only two inputs are implemented for NOR."
    p = (x, y)
    connections = Dict{Symbol, XYTuple}(
        inputs[1] => p + 0.5 * (-0.7025sx, -0.5sy),
        inputs[2] => p + 0.5 * (-0.7025sx, 0.5sy),
        :out      => p + 0.5 * (1.4sx, 0.0sy)
    )
    composition = compose(
        context(
            x-0.6sx, y-0.55sy, 1.35sx, 1.1sy,
            units = UnitBox(-1.2, -1.1, 2.7, 2.2)
        ),
        xor_base(),
        (context(), circle(1.2sx, 0, 0.2sx), fill(nothing))
    )
    return Block(composition, connections)
end

function not_gate(x, y, inputs, output=:out, sx=1.0, sy=1.0)
    p = (x, y)
    in_pos = partition(
        p + 0.5 * (-1.0sx, -1.0sy),
        p + 0.5 * (-1.0sx, 1.0sy),
        length(inputs)
    )
    connections = Dict{Symbol, XYTuple}(
        map((label, pos) -> Pair(label, pos), inputs, in_pos)...,
        :out => p + 0.5 * (1.0sx, 0.0sy)
    )
    composition = compose(
        context(
            x-0.55sx, y-0.55sy, 1.1sx, 1.1sy,
            units = UnitBox(-1.1, -1.1, 2.2, 2.2)
        ),
        (
            context(),
            polygon([(-1, -1), (-1, 1), (0.68, 0.)]),
            circle(0.85, 0, 0.15),
            fill(nothing)
        )
    )
    return Block(composition, connections)
end

# function bezier(anchor1, ctrl1, ctrl2, anchor2)
#     f(t) = (1 - t)^3 * anchor1 +
#         3 * t * (1 - t)^2 * (anchor1 + ctrl1) +
#         3 * t^2 * (1 - t) * (anchor2 - ctrl2) +
#         t^3 * anchor2
# end


text:
using Dates

@testset "Date.TimeTypes" begin

res = Parsers.xparse(Date, "")
x, code = res.val, res.code
@test code == INVALID | EOF
res = Parsers.xparse(Date, "2018-01-01")
x, code = res.val, res.code
@test x === Date(2018, 1, 01)
@test code == OK | EOF
res = Parsers.xparse(DateTime, "2018-01-01")
x, code = res.val, res.code
@test x === DateTime(2018, 1, 01)
@test code == OK | EOF
res = Parsers.xparse(Time, "01:02:03")
x, code = res.val, res.code
@test x === Time(1, 2, 3)
@test code == OK | EOF

res = Parsers.xparse(Date, "")
x, code = res.val, res.code
@test code == INVALID | EOF
res = Parsers.xparse(Date, "\"\"")
x, code = res.val, res.code
@test code == QUOTED | INVALID | EOF
res = Parsers.xparse(Date, "2018-01-01")
x, code = res.val, res.code
@test x === Date(2018, 1, 01)
@test code == OK | EOF
res = Parsers.xparse(Date, "\"2018-01-01\"")
x, code = res.val, res.code
@test x === Date(2018, 1, 01)
@test code == QUOTED | OK | EOF
res = Parsers.xparse(DateTime, "\"2018-01-01")
x, code = res.val, res.code
@test code == OK | QUOTED | EOF | INVALID_QUOTED_FIELD
res = Parsers.xparse(Date, "\"abcd\"")
x, code = res.val, res.code
@test code == QUOTED | INVALID | EOF

res = Parsers.xparse(Date, "NA", sentinel=["NA"])
x, code = res.val, res.code
@test code === SENTINEL | EOF
res = Parsers.xparse(Date, "\\N", sentinel=["\\N"])
x, code = res.val, res.code
@test code === SENTINEL | EOF
res = Parsers.xparse(Date, "NA2", sentinel=["NA"])
x, code = res.val, res.code
@test code === SENTINEL | INVALID_DELIMITER | EOF
res = Parsers.xparse(Date, "-", sentinel=["-"])
x, code = res.val, res.code
@test code === SENTINEL | EOF
res = Parsers.xparse(Date, "£", sentinel=["£"])
x, code = res.val, res.code
@test code === SENTINEL | EOF
res = Parsers.xparse(Date, "null")
x, code = res.val, res.code
@test code === INVALID_DELIMITER | EOF
res = Parsers.xparse(Date, ",")
x, code = res.val, res.code
@test code === INVALID | DELIMITED
res = Parsers.xparse(Date, "1,")
x, code = res.val, res.code
@test x === Date(1)
@test code === OK | DELIMITED

@test Parsers.parse(DateTime, "1996/Feb/15", Parsers.Options(dateformat="yy/uuu/dd")) === DateTime(1996, 2, 15)
@test Parsers.parse(DateTime, "1996, Jan, 15", Parsers.Options(dateformat="yyyy, uuu, dd")) === DateTime(1996, 1, 15)

@test_throws Parsers.Error Parsers.parse(Date, "2020-05-32")
@test_throws Parsers.Error Parsers.parse(DateTime, "2020-05-32")
@test_throws Parsers.Error Parsers.parse(Time, "25:00:00")
@test_throws Parsers.Error Parsers.parse(DateTime, "2020-05-05T00:00:60")

end

@testset "Date.TimeTypes IO" begin

    res = Parsers.xparse(Date, IOBuffer(""))
    x, code = res.val, res.code

    @test code == INVALID | EOF
    res = Parsers.xparse(Date, IOBuffer("2018-01-01"))
    x, code = res.val, res.code
    @test x === Date(2018, 1, 01)
    @test code == OK | EOF
    res = Parsers.xparse(DateTime, IOBuffer("2018-01-01"))
    x, code = res.val, res.code
    @test x === DateTime(2018, 1, 01)
    @test code == OK | EOF
    res = Parsers.xparse(Time, IOBuffer("01:02:03"))
    x, code = res.val, res.code
    @test x === Time(1, 2, 3)
    @test code == OK | EOF

    res = Parsers.xparse(Date, IOBuffer("\"\""))
    x, code = res.val, res.code

    @test code == QUOTED | INVALID | EOF
    res = Parsers.xparse(Date, IOBuffer("\"2018-01-01\""))
    x, code = res.val, res.code
    @test x === Date(2018, 1, 01)
    @test code == QUOTED | OK | EOF
    res = Parsers.xparse(DateTime, IOBuffer("\"2018-01-01"))
    x, code = res.val, res.code
    @test code == OK | QUOTED | EOF | INVALID_QUOTED_FIELD
    res = Parsers.xparse(Date, IOBuffer("\"abcd\""))
    x, code = res.val, res.code

    @test code == QUOTED | INVALID | EOF

    res = Parsers.xparse(Date, IOBuffer("NA"), sentinel=["NA"])
    x, code = res.val, res.code

    @test code === SENTINEL | EOF
    res = Parsers.xparse(Date, IOBuffer("\\N"), sentinel=["\\N"])
    x, code = res.val, res.code

    @test code === SENTINEL | EOF
    res = Parsers.xparse(Date, IOBuffer("NA2"), sentinel=["NA"])
    x, code = res.val, res.code

    @test code === SENTINEL | INVALID_DELIMITER | EOF
    res = Parsers.xparse(Date, IOBuffer("-"), sentinel=["-"])
    x, code = res.val, res.code

    @test code === SENTINEL | EOF
    res = Parsers.xparse(Date, IOBuffer("£"), sentinel=["£"])
    x, code = res.val, res.code

    @test code === SENTINEL | EOF
    res = Parsers.xparse(Date, IOBuffer("null"))
    x, code = res.val, res.code

    @test code === INVALID_DELIMITER | EOF
    res = Parsers.xparse(Date, IOBuffer(","))
    x, code = res.val, res.code

    @test code === INVALID | DELIMITED
    res = Parsers.xparse(Date, IOBuffer("1,"))
    x, code = res.val, res.code
    @test x === Date(1)
    @test code === OK | DELIMITED

    @test Parsers.parse(DateTime, IOBuffer("1996/Feb/15"), Parsers.Options(dateformat="yy/uuu/dd")) === DateTime(1996, 2, 15)
    @test Parsers.parse(DateTime, IOBuffer("1996, Jan, 15"), Parsers.Options(dateformat="yyyy, uuu, dd")) === DateTime(1996, 1, 15)

    @test_throws Parsers.Error Parsers.parse(Date, IOBuffer("2020-05-32"))
    @test_throws Parsers.Error Parsers.parse(DateTime, IOBuffer("2020-05-32"))
    @test_throws Parsers.Error Parsers.parse(Time, IOBuffer("25:00:00"))
    @test_throws Parsers.Error Parsers.parse(DateTime, IOBuffer("2020-05-05T00:00:60"))

end

@test_throws Parsers.Error Parsers.parse(DateTime, "a02/15/1996 25:00", Parsers.Options(dateformat=Parsers.Format("mm/dd/yyyy HH:MM")))
@test_throws Parsers.Error Parsers.parse(DateTime, "02/15/1996 25:00", Parsers.Options(dateformat=Parsers.Format("mm/dd/yyyy HH:MM")))
@test_throws Parsers.Error Parsers.parse(DateTime, "1996-Jan-15", Parsers.Options(dateformat="yy-mm-dd"))
@test_throws Parsers.Error Parsers.parse(DateTime, "96/2/15", Parsers.Options(dateformat="yy/uuu/dd"))
@test_throws Parsers.Error Parsers.parse(DateTime, "2017-Mar-17 00:00:00.1234", Parsers.Options(dateformat="y-u-d H:M:S.s"))
@test_throws Parsers.Error Parsers.parse(DateTime, "96/2/15", Parsers.Options(dateformat="yy𒀱uuu/dd"))

@test_throws Parsers.Error Parsers.parse(Time, "24:00")  # invalid hours
@test_throws Parsers.Error Parsers.parse(Time, "00:60")  # invalid minutes
@test_throws Parsers.Error Parsers.parse(Time, "00:00:60")  # invalid seconds
@test_throws Parsers.Error Parsers.parse(Time, "20:03:20", Parsers.Options(dateformat="HH:MM"))  # too much precision
@test_throws Parsers.Error Parsers.parse(Time, "10:33:51", Parsers.Options(dateformat="yyyy-mm-dd HH:MM:SS"))  # Time can't hold year/month/day
@test Parsers.parse(Time, "2021-06-26 10:33:51", Parsers.Options(dateformat="yyyy-mm-dd HH:MM:SS")) == Time(10, 33, 51)

Dates.LOCALES["french"] = Dates.DateLocale(
    ["janvier", "février", "mars", "avril", "mai", "juin",
        "juillet", "août", "septembre", "octobre", "novembre", "décembre"],
    ["janv", "févr", "mars", "avril", "mai", "juin",
        "juil", "août", "sept", "oct", "nov", "déc"],
    ["lundi", "mardi", "mercredi", "jeudi", "vendredi", "samedi", "dimanche"],
    [""],
)

Dates.LOCALES["gobblygook"] = Dates.DateLocale(
    ["𒀱"],
    ["𒀱"],
    [""],
    [""],
)

testcases = [
    ("yy-mm-dd", "96-01-15", Dates.DateTime(96, 1, 15)),
    ("yy-mm-dd", "96-1-15", Dates.DateTime(96, 1, 15)),
    ("yy-mm-dd", "96-1-1", Dates.DateTime(96, 1, 1)),
    ("yy-mm-dd", "1996-1-15", Dates.DateTime(1996, 1, 15)),
    ("yy/uuu/dd", "96/Feb/15", Dates.DateTime(96, 2, 15)),
    ("yy/uuu/dd", "1996/Feb/15", Dates.DateTime(1996, 2, 15)),
    ("yy/uuu/dd", "96/Feb/1", Dates.DateTime(96, 2, 1)),

    ("yyyy, uuu, dd", "1996, Jan, 15", Dates.DateTime(1996, 1, 15)),
    ("yyyy.U.dd", "1996.February.15", Dates.DateTime(1996, 2, 15)),
    ("yyyymmdd", "19960315", Dates.DateTime(1996, 3, 15)),
    ("yyyy-mm-dd HH:MM:SS", "1996-12-15 10:00:00", Dates.DateTime(1996, 12, 15, 10)),
    ("ymd", "999", Dates.DateTime(9, 9, 9)),
    ("/yyyy/m/d", "/1996/5/15", Dates.DateTime(1996, 5, 15)),
    ("yyyy年mm月dd日", "2009年12月01日", Dates.DateTime(2009, 12, 1)),
    ("yyyy𒀱mm𒀱dd", "2021𒀱6𒀱28", Dates.Date(2021, 6, 28)),

    (Parsers.Format("dd uuuuu YYYY", "french"), "28 mai 2014", Dates.DateTime(2014, 5, 28)),
    (Parsers.Format("dd uuuuu yyyy", "french"), "28 févr 2014", Dates.DateTime(2014, 2, 28)),
    (Parsers.Format("dd uuuuu yyyy", "french"), "28 août 2014", Dates.DateTime(2014, 8, 28)),
    (Parsers.Format("dd u yyyy", "french"), "28 avril 2014", Dates.DateTime(2014, 4, 28)),
    (Parsers.Format("dduuuuyyyy", "french"), "28mai2014", Dates.DateTime(2014, 5, 28)),
    (Parsers.Format("dduuuuyyyy", "french"), "28août2014", Dates.DateTime(2014, 8, 28)),
    (Parsers.Format("dd uuuuu YYYY", "gobblygook"), "28 𒀱 2014", Dates.DateTime(2014, 1, 28)),

    ("[HH:MM:SS.sss]", "[14:51:00.118]", Dates.DateTime(1, 1, 1, 14, 51, 0, 118)),
    ("HH:MM:SS.sss", "14:51:00.118", Dates.DateTime(1, 1, 1, 14, 51, 0, 118)),
    ("[HH:MM:SS.sss?", "[14:51:00.118?", Dates.DateTime(1, 1, 1, 14, 51, 0, 118)),
    ("?HH:MM:SS.sss?", "?14:51:00.118?", Dates.DateTime(1, 1, 1, 14, 51, 0, 118)),
    ("xHH:MM:SS.sss]", "x14:51:00.118]", Dates.DateTime(1, 1, 1, 14, 51, 0, 118)),
    ("HH:MM:SS.sss]", "14:51:00.118]", Dates.DateTime(1, 1, 1, 14, 51, 0, 118)),

    (Dates.RFC1123Format, "Sat, 23 Aug 2014 17:22:15", Dates.DateTime(2014, 8, 23, 17, 22, 15)),
    ("E, dd u yyyy HH:MM:SS", "Saturday, 23 Aug 2014 17:22:15", Dates.DateTime(2014, 8, 23, 17, 22, 15)),
    # milliseconds
    ("y-u-d H:M:S.s", "2017-Mar-17 00:00:00.0000", Dates.DateTime(2017, 3, 17, 0, 0, 0, 0)),
    ("y-u-d H:M:S.s", "2017-Mar-17 00:00:00.1", Dates.DateTime(2017, 3, 17, 0, 0, 0, 100)),
    ("y-u-d H:M:S.s", "2017-Mar-17 00:00:00.12", Dates.DateTime(2017, 3, 17, 0, 0, 0, 120)),
    ("y-u-d H:M:S.s", "2017-Mar-17 00:00:00.123", Dates.DateTime(2017, 3, 17, 0, 0, 0, 123)),
    ("y-u-d H:M:S.s", "2017-Mar-17 00:00:00.1230", Dates.DateTime(2017, 3, 17, 0, 0, 0, 123)),
]

for useio in (true, false)
    for case in testcases
        fmt, str, dt = case
        @test Parsers.parse(typeof(dt), useio ? IOBuffer(str) : str, Parsers.Options(dateformat=fmt)) == dt
    end
end

@static if VERSION >= v"1.3-DEV"
@testset "AM/PM" begin
    for (t12,t24) in (("12:00am","00:00"), ("12:07am","00:07"), ("01:24AM","01:24"),
                    ("12:00pm","12:00"), ("12:15pm","12:15"), ("11:59PM","23:59"))
        d = DateTime("2018-01-01T$t24:00")
        t = Time("$t24:00")
        for HH in ("HH","II")
            @test Parsers.parse(DateTime, "2018-01-01 $t12", Parsers.Options(dateformat="yyyy-mm-dd $HH:MMp")) == d
            @test Parsers.parse(Time, "$t12", Parsers.Options(dateformat="$HH:MMp")) == t
        end
    end
    for bad in ("00:24am", "00:24pm", "13:24pm", "2pm", "12:24p.m.", "12:24 pm", "12:24pµ")
        @test_throws Parsers.Error Parsers.parse(Time, bad, Parsers.Options(dateformat="II:MMp"))
    end
    # if am/pm is missing, defaults to 24-hour clock
    @test Parsers.parse(Time, "13:24", Parsers.Options(dateformat="II:MMp")) == Parsers.parse(Time, "13:24", Parsers.Options(dateformat="HH:MM"))
end
end

ext: jl


text:
@testset "check_calibration_arguments(::AbstractVector{<: Real}, ::AbstractVector{<: Bool}))" begin
    @test isnothing(
        Escape.check_calibration_arguments([0.1, 0.2, 0.3], [true, false, true])
    )
    @test_throws ErrorException Escape.check_calibration_arguments(
        [0.1, 0.2], [true, false, true]
    )
    @test_throws ErrorException Escape.check_calibration_arguments(
        [0.3, 1, 1.1], [true, false, true]
    )
end

@testset "binned_intervals(::AbstractVector{<: Real}, ::AbstractVector{<: Bool})" begin
    theta = collect(0.3:0.1:1)
    y = [0, 0, 0, 1, 0, 1, 1, 1]

    expected_1 = mean(theta[1:4])
    expected_2 = mean(theta[5:8])
    observed_1 = mean(y[1:4])
    observed_2 = mean(y[5:8])

    interval_1 = quantile(rand(Beta(2, 4), 50000), (0.025, 0.975))
    interval_2 = quantile(rand(Beta(4, 2), 50000), (0.025, 0.975))

    @test Escape.binned_intervals(theta, Bool.(y), bins = 2) isa DataFrames.DataFrame
    df = Escape.binned_intervals(theta, Bool.(y), bins = 2)

    theta = range(0, 1, length = 100000)
    y = map(x -> rand(Bernoulli(x)), theta)
    df = Escape.binned_intervals(theta, y, bins = 10, lower = 0.0001, upper = 0.9999)

    for row in eachrow(df)
        @test row[:lower] <= row[:observed] <= row[:upper]
    end
end

@testset "@recipe f(::HLAModelResult)" begin
    result = @suppress Escape.run(
        Escape.HLAModel{1}(), Escape.HLADataset("Test").data[1],
        iter = 10, chains = 1, warmup = 10
    )

    @test @suppress Escape.calibration_plot(result, samples = 10) isa Plots.Plot
end


text:
export CircuitTable

"""
    CircuitTable{I<:Integer}

Holding all circuits describing a mixed cell cone ``C_M`` of mixed cell ``M``.
"""
struct CircuitTable{I<:Integer}
    # circuits is an (2n + 1) × m Matrix
    # The first 2n rows are the entries corresponding to the columns of the mixed cell
    # the last row corresponds to the index of the choosen extra column.
    circuits::Matrix{I}
    # linear indices of the mixed cell which corresponds to this circuit table
    indices::Vector{I} # length 2n
    # The entries of the columns *not* in the mixed cell, i.e.
    # for which we have an inequality.
    entries::BitVector
end

function CircuitTable(A::CayleyConfiguration{I}, cellindices::CellIndices) where {I<:Integer}
    # construct a circuit table freshly
    indices = linearindices(A, cellindices)

    n2, m = size(A)

    circuits = fill(zero(I), n2+1, m) # (2n + 1) × m

    # construct 2n × 2n submatrix indexed by the cell columns
    D = A[:, indices]

    # For each column γ **not** in the cell we need to compute a circuit c
    # which represents a facet of the MixedCell cone. This is a generator
    # of the nullspace of [D A[:,γ]]. We compute this by
    # [I D^-1 A[:,γ]] = 0 and scale the result such that the first entry is det(D)
    LU = lufact(D)
    D_inv = inv(LU)
    γ = round(I, det(LU))

    k = 1
    entries = trues(m)
    for i=1:m
        if k ≤ n2 && i == indices[k]
            k += 1
            entries[i] = false
            continue
        end
        c = D_inv * A[:,i]
        scale!(c, abs(γ))
        for j=1:n2
            circuits[j, i] = round(I, c[j])
        end
        circuits[n2+1, i] = γ
    end

    CircuitTable(circuits, indices, entries)
end

Base.size(T::CircuitTable) = size(T.circuits)
Base.size(T::CircuitTable, I) = size(T.circuits, I)

"""
    dot(circuittable, j, ω)

Compute the dot product of the `j`-th circuit with the vector `ω`.
"""
function Base.LinAlg.dot(table::CircuitTable{I}, j, ω::AbstractVector{I}) where {I}
    N = size(table, 1)
    sum = zero(I)
    for k = 1:N-1
        @inbounds sum = table.circuits[k,j] * ω[table.indices[k]]
    end
    sum + table.circuits[N, j] * ω[j]
end


text:
function SimplexSolver()
    return SimplexSolver(())
end

function SimplexSolver(arg0::jdouble)
    return SimplexSolver((jdouble,), arg0)
end

function SimplexSolver(arg0::jdouble, arg1::jint)
    return SimplexSolver((jdouble, jint), arg0, arg1)
end

function SimplexSolver(arg0::jdouble, arg1::jint, arg2::jdouble)
    return SimplexSolver((jdouble, jint, jdouble), arg0, arg1, arg2)
end

function do_optimize(obj::SimplexSolver)
    return jcall(obj, "doOptimize", PointValuePair, ())
end

function optimize(obj::SimplexSolver, arg0::Vector{OptimizationData})
    return jcall(obj, "optimize", PointValuePair, (Vector{OptimizationData},), arg0)
end


text:
module FundStabABM

runmodel() = include("./src/runmodel.jl")

end


text:
using Documenter, TruncatedDistributions

makedocs(sitename="TruncatedDistributions.jl")

deploydocs(repo = "github.com/yoninazarathy/TruncatedDistributions.jl.git")


text:
#
# this example is an adoption of the chain example
# in Programming Elixir ch. 15
# see https://media.pragprog.com/titles/elixir16/code/spawn/chain.exs
#
using .Threads

function counter(next::Channel{Int})
    return ch -> put!(next, take!(ch)+1)
end

function create_processes(start::Channel{Int}, n::Int; spawn=false)
    next = start
    for i in 1:n
        next = Channel{Int}(counter(next), spawn=spawn)
    end

    put!(next, 0)

    take!(start)
end

function startonthread(id::Int, start::Channel{Int}, n::Int; spawn=false)
    t = Task(nothing)
    @threads for i in 1:nthreads()
        if i == id
            t = @async create_processes(start, n, spawn=spawn)
        end
    end
    fetch(t)
end

const start = Channel{Int}()


text:
using ACE, Test, Printf, LinearAlgebra, StaticArrays, BenchmarkTools

##
@testset "ACE.jl" begin
    # --------------------------------------------
    #   basic polynomial basis building blocks
    @testset "Ylm" begin include("polynomials/test_ylm.jl") end
    @testset "TestWigner" begin include("testing/test_wigner.jl") end
    @testset "Transforms" begin include("transforms/test_transforms.jl") end
    @testset "Lambdas" begin include("transforms/test_lambdas.jl") end
    @testset "OrthogonalPolynomials" begin include("polynomials/test_orthpolys.jl") end

    # --------------------------------------------
    #  states, ..
    @testset "States" begin include("test_states.jl") end

    # --------------------------------------------
    # core permutation-invariant functionality
    @testset "1-Particle Basis"  begin include("test_1pbasis.jl") end
    @testset "Categorical1pBasis" begin include("test_discrete.jl") end
    @testset "PIBasis" begin include("test_pibasis.jl") end

    # --------------------------------------------
    #   O(3) equi-variance
    @testset "Clebsch-Gordan" begin include("test_cg.jl") end
    @testset "SymmetricBasis" begin include("test_symmbasis.jl") end
    @testset "EuclideanVector" begin include("test_euclvec.jl") end
    @testset "EuclideanMatrix" begin include("test_EuclideanMatrix.jl") end
    @testset "Multiple SH Bases" begin include("test_multish.jl") end

    # Model tests
    @testset "LinearACEModel"  begin include("test_linearmodel.jl") end
    @testset "MultipleProperties"  begin include("test_multiprop.jl") end
    @testset "AD-LinearACEModel"  begin include("test_admodel.jl") end

    # Experimental material
    @testset "Sparsification" begin include("test_sparsify.jl") end
    # @testset "Multipliers" begin include("test_multiplier.jl") end
end

    # -----------------------------------------
    #    old tests to be re-introduced - maybe
    # include("test_real.jl")
    # include("test_orth.jl")
    # include("bonds/test_cylindrical.jl")
    # include("bonds/test_fourier.jl")


text:
#///////////////////////////////////////
#// File Name: RATiLQR.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2020/10/28
#// Description: Julia package for RAT iLQR and baseline methods
#///////////////////////////////////////

module RATiLQR

using Distributions
using Distributed
import Future.randjump
using ForwardDiff
using LinearAlgebra
using Printf
using Random
using Statistics

# Optimal Control Problems
export
    OptimalControlProblem,
    FiniteHorizonRiskSensitiveOptimalControlProblem,
    FiniteHorizonGenerativeOptimalControlProblem
include("optimal_control_problems.jl")

# iterative Linear-Exponential-Quadratic-Gaussian
export
    simulate_dynamics,
    integrate_cost,
    ILEQGSolver,
    initialize!,
    ApproximationResult,
    approximate_model,
    DynamicProgrammingResult,
    solve_approximate_dp!,
    solve_approximate_dp,
    increase_μ_and_Δ!,
    decrease_μ_and_Δ!,
    line_search!,
    step!,
    solve!
include("ileqg.jl")

# Cross Entropy Bilevel Optimization (i.e. RATiLQR)
export
    CrossEntropyBilevelOptimizationSolver,
    compute_value_worker,
    compute_cost,
    compute_cost_serial,
    get_positive_samples,
    step!,
    solve!
include("cross_entropy_bilevel_optimization.jl")

# Nelder-Mead Simplex Bilevel Optimization
export
    NelderMeadBilevelOptimizationSolver,
    compute_cost_worker,
    initialize!,
    step!,
    solve!
include("nelder_mead_bilevel_optimization.jl")

# PETS
export
    CrossEntropyDirectOptimizationSolver,
    initialize!,
    compute_cost_worker,
    compute_cost,
    compute_cost_serial,
    get_elite_samples,
    compute_new_distribution,
    step!,
    solve!
include("pets.jl")

end # module


text:
using Random

import MinimalRLCore

module MountainCarConst
const vel_limit = (-0.07, 0.07)
const pos_limit = (-1.2, 0.5)
const pos_initial_range = (-0.2, 0.2)
const vel_initial_range = (-0.04,0.04)

# const pos_initial_range = (-1.2, 0.5)
# const vel_initial_range = (-0.07, 0.07)

const Reverse=1
const Neutral=2
const Accelerate=3
end

"""
MountainCar(pos=0.0, vel=0.0, normalized::Bool=false)
MountainCar(rng::AbstractRNG, normalized::Bool=false)

    The mountain car environment.

"""
mutable struct MountainCar <: MinimalRLCore.AbstractEnvironment
    pos::Float64
    vel::Float64
    normalized::Bool
    function MountainCar(pos=0.0, vel=0.0, normalized::Bool=false)
        mcc = MountainCarConst
        @boundscheck (pos >= mcc.pos_limit[1] && pos <= mcc.pos_limit[2])
        @boundscheck (vel >= mcc.vel_limit[1] && vel <= mcc.vel_limit[2])
        new(pos, vel, normalized)
    end
end

MinimalRLCore.get_actions(env::MountainCar) = [MountainCarConst.Reverse, MountainCarConst.Neutral, MountainCarConst.Accelerate]
valid_action(env::MountainCar, action) = action in [MountainCarConst.Reverse, MountainCarConst.Neutral, MountainCarConst.Accelerate]

function MinimalRLCore.reset!(env::MountainCar, rng::AbstractRNG=Random.GLOBAL_RNG)
    env.pos = (rand(rng)*(MountainCarConst.pos_initial_range[2]
                          - MountainCarConst.pos_initial_range[1])
               + MountainCarConst.pos_initial_range[1])
    # env.vel = 0.0
    env.vel = rand(rng)*(MountainCarConst.vel_initial_range[2] - MountainCarConst.vel_initial_range[1]) + MountainCarConst.vel_initial_range[1]
end

function MinimalRLCore.reset!(env::MountainCar,
                              start_state::T) where {T<:AbstractArray}
    if env.normalized == false
        env.pos = start_state[1]
        env.vel = start_state[2]
    else
        pos_limit = MountainCarConst.pos_limit
        vel_limit = MountainCarConst.vel_limit
        env.pos = (start_state[1]*(pos_limit[2] - pos_limit[1])) + pos_limit[1]
        env.vel = (start_state[2]*(vel_limit[2] - vel_limit[1])) + vel_limit[1]
    end
end

function MinimalRLCore.environment_step!(env::MountainCar, action, rng::AbstractRNG=Random.GLOBAL_RNG)

    @boundscheck valid_action(env, action)
    env.vel =
        clamp(env.vel + (action - 2)*0.001 - 0.0025*cos(3*env.pos),
              MountainCarConst.vel_limit...)
    env.pos = clamp(env.pos + env.vel,
                    MountainCarConst.pos_limit...)
end

function MinimalRLCore.get_reward(env::MountainCar) # -> determines if the agent_state is terminal
    if env.pos >= MountainCarConst.pos_limit[2]
        return -0.01f0
    end
    return -0.01f0
end

function MinimalRLCore.is_terminal(env::MountainCar) # -> determines if the agent_state is terminal
    return env.pos >= MountainCarConst.pos_limit[2]
end

function MinimalRLCore.get_state(env::MountainCar)
    if env.normalized
        return get_normalized_state(env)
    else
        return Float32[env.pos, env.vel]
    end
end

function get_normalized_state(env::MountainCar,state)
    pos_limit = MountainCarConst.pos_limit
    vel_limit = MountainCarConst.vel_limit
    return Float32[(state[1] - pos_limit[1])/(pos_limit[2] - pos_limit[1]),
                   (state[2] - vel_limit[1])/(vel_limit[2] - vel_limit[1])]
end

function get_normalized_state(env::MountainCar)
    pos_limit = MountainCarConst.pos_limit
    vel_limit = MountainCarConst.vel_limit
    return Float32[(env.pos - pos_limit[1])/(pos_limit[2] - pos_limit[1]),
                   (env.vel - vel_limit[1])/(vel_limit[2] - vel_limit[1])]
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using BinaryProvider # requires BinaryProvider 0.5.0 or later
using Libdl

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
@info "####Prefix: $prefix"
mkpath(prefix.path)
products = [
    LibraryProduct(prefix, ["libz"], :libz),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.4"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc) => ("$bin_prefix/Zlib.v1.2.11.aarch64-linux-gnu.tar.gz", "3431c1a7d937cbad379ae5fe78c11dd0d425a6bce02b12c9d7c89e0a3173af97"),
    Linux(:aarch64, libc=:musl) => ("$bin_prefix/Zlib.v1.2.11.aarch64-linux-musl.tar.gz", "7f89328d100f16ab75d37759cdb4df98d8f5b192146eacbe23fa859c6d7bca09"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf) => ("$bin_prefix/Zlib.v1.2.11.arm-linux-gnueabihf.tar.gz", "04abd9100d2d24a56e50d536679c2df317d6c344b3ae06feb77681533377652f"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf) => ("$bin_prefix/Zlib.v1.2.11.arm-linux-musleabihf.tar.gz", "fd4fd26462039de7135b5ad09f0c86ec441afd9bb7feb2689f705bc9b087931f"),
    Linux(:i686, libc=:glibc) => ("$bin_prefix/Zlib.v1.2.11.i686-linux-gnu.tar.gz", "41495344e20619158a05561db8a7ccb9c57192dc5038744419bb9bd58eba7b7a"),
    Linux(:i686, libc=:musl) => ("$bin_prefix/Zlib.v1.2.11.i686-linux-musl.tar.gz", "c5d101f21c440156fdc8616bbe9f13f437dfed634cb17dbc65ed83c6b180533c"),
    Windows(:i686) => ("$bin_prefix/Zlib.v1.2.11.i686-w64-mingw32.tar.gz", "84d4d5037459ec706029c78a94c98b07760943662a4adbe21cd3d353ac3a0484"),
    Linux(:powerpc64le, libc=:glibc) => ("$bin_prefix/Zlib.v1.2.11.powerpc64le-linux-gnu.tar.gz", "c485627191c75411eaa66088427813a182c7ca0549f0dc21971c9452c45d8167"),
    MacOS(:x86_64) => ("$bin_prefix/Zlib.v1.2.11.x86_64-apple-darwin14.tar.gz", "a28fb652f94f1b1548197a3b3c9a71ec56fa50bb7d39454c5e12320335d30934"),
    Linux(:x86_64, libc=:glibc) => ("$bin_prefix/Zlib.v1.2.11.x86_64-linux-gnu.tar.gz", "1c3b1c8520713f98d3f605ee1ca5e2e3656d92ddb6441abeeeff0ae12a11a620"),
    Linux(:x86_64, libc=:musl) => ("$bin_prefix/Zlib.v1.2.11.x86_64-linux-musl.tar.gz", "17da1e44f7a815f5dc6eb66e3c7df8c80846c67d5212148a33fb68ea677e38dc"),
    FreeBSD(:x86_64) => ("$bin_prefix/Zlib.v1.2.11.x86_64-unknown-freebsd11.1.tar.gz", "87522b29c8dfb7681209fea541dcc09bd5863f55df9cb1de1645ec8484aad7b8"),
    Windows(:x86_64) => ("$bin_prefix/Zlib.v1.2.11.x86_64-w64-mingw32.tar.gz", "13934d974c5b1fd99897897b9af4ef7cce1025a1cdf1a57f14a9dd8e0258508a"),
)

# A simple source build fallback for platforms not supported by BinaryBuilder
# Assumes that tar, GNU make, and a C compiler are available
function sourcebuild()
    @info "Trying to install zlib from the source code."
    srcdir = joinpath(@__DIR__, "src")
    libdir = joinpath(@__DIR__, "lib")
    z = "zlib-1.2.11"
    for d = [srcdir, libdir]
        isdir(d) && rm(d, force=true, recursive=true)
        mkpath(d)
    end
    download("https://zlib.net/$(z).tar.gz", joinpath(srcdir, "$(z).tar.gz"))
    cd(srcdir) do
        run(`tar xzf $(z).tar.gz`)
    end
    cd(joinpath(srcdir, z)) do
        run(`./configure --prefix=.`)
        make = Sys.isbsd() && !Sys.isapple() ? `gmake` : `make`
        run(`$make -j$(Sys.CPU_THREADS)`)
    end
    libz = nothing
    for f in readdir(joinpath(srcdir, z))
        if f == "libz." * Libdl.dlext
            libz = joinpath(srcdir, z, f)
            break
        end
    end
    libz === nothing && error("zlib was unable to build properly")
    open(joinpath(@__DIR__, "deps.jl"), "w") do io
        println(io, """
            function check_deps()
                ptr = Libdl.dlopen_e("$libz")
                loaded = ptr != C_NULL
                Libdl.dlclose(ptr)
                if !loaded
                    error("Unable to load zlib from $libz. Please rerun " *
                          "`Pkg.build(\\"CodecZlib\\")` and restart Julia.")
                end
            end
            const libz = "$libz"
            """)
    end
end

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    @warn "ZlibBuilder provides no prebuilt binary for your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\")."
    sourcebuild()
else
    # If we have a download, and we are unsatisfied (or the version we're
    # trying to install is not itself installed) then load it up!
    if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
        # Download and install binaries
        install(dl_info...; prefix=prefix, force=true, verbose=verbose)
    end

    # Write out a deps.jl file that will contain mappings for our products
    write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)
end


text:
using DistanceTransforms
using CUDA
using Test
using TestSetExtensions


text:
import Random, Base.size

BASE_WALLS = [0 0 0 0 0 1 0 0 0 0 0;
              0 0 0 0 0 1 0 0 0 0 0;
              0 0 0 0 0 0 0 0 0 0 0;
              0 0 0 0 0 1 0 0 0 0 0;
              0 0 0 0 0 1 0 0 0 0 0;
              1 0 1 1 1 1 0 0 0 0 0;
              0 0 0 0 0 1 1 1 0 1 1;
              0 0 0 0 0 1 0 0 0 0 0;
              0 0 0 0 0 1 0 0 0 0 0;
              0 0 0 0 0 0 0 0 0 0 0;
              0 0 0 0 0 1 0 0 0 0 0;]

UP = 1
RIGHT = 2
DOWN = 3
LEFT = 4

ROOM_TOP_LEFT = 1
ROOM_TOP_RIGHT = 2
ROOM_BOTTOM_LEFT = 3
ROOM_BOTTOM_RIGHT = 4

mutable struct environment
    state
    walls::Array{Bool, 2}
    environment() = new([0,0], convert(Array{Bool}, BASE_WALLS))
    environment(walls::Array{Int64, 2}) = new([0,0], convert(Array{Bool}, walls))
    function environment(size, wall_list)
        walls = zeros(Bool, size[1], size[2])
        for wall in wall_list
            walls[wall[1], wall[2]] = true
        end
        new((0,0), walls)
    end
end

is_terminal(env::environment) = false
reward(env::environment) = 0
state(env::environment) = env.state
is_wall(env::environment, state) = env.walls[state[1], state[2]] == 1
random_state(env::environment, rng) = [rand(rng, 1:size(env.walls)[1]), rand(rng, 1:size(env.walls)[2])]
size(env::environment) = size(env.walls)

function which_room(state)
    room = -1
    if state[1] < 6
        # LEFT
        if state[2] < 6
            # TOP
            room = ROOM_TOP_LEFT
        else
            # Bottom
            room = ROOM_BOTTOM_LEFT
        end
    else
        # RIGHT
        if state[2] < 7
            # TOP
            room = ROOM_TOP_RIGHT
        else
            # Bottom
            room = ROOM_BOTTOM_RIGHT
        end
    end
    return room
end

function start!(env::environment; rng=Random.GLOBAL_RNG)
    # state = random_state(rng)
    # state = rand(rng, [1:size(env.walls)[1], 1:size(env.walls)[2]])
    state = random_state(env, rng)
    while !env.walls[state[1], state[2]]
        state = random_state(env, rng)
    end
    env.state = state
    return state
end

function step!(env::environment, action)
    (env.state, rew, terminal) = step(env, action)
    return env.state, rew, terminal
end

function step(env::environment, action)
    # print("hello world")
    next_state = copy(env.state)

    if action == UP
        next_state[1] -= 1
    elseif action == DOWN
        next_state[1] += 1
    elseif action == RIGHT
        next_state[2] += 1
    elseif action == LEFT
        next_state[2] -= 1
    end

    next_state[1] = clamp(next_state[1], 1, size(env.walls)[1])
    next_state[2] = clamp(next_state[2], 1, size(env.walls)[2])
    if is_wall(env, next_state)
        next_state = env.state
    end

    return next_state, 0, false
end

function step(env::environment, state, action)
    # print("hello world")
    next_state = copy(state)

    if action == UP
        next_state[1] -= 1
    elseif action == DOWN
        next_state[1] += 1
    elseif action == RIGHT
        next_state[2] += 1
    elseif action == LEFT
        next_state[2] -= 1
    end

    next_state[1] = clamp(next_state[1], 1, size(env.walls)[1])
    next_state[2] = clamp(next_state[2], 1, size(env.walls)[2])
    if is_wall(env, next_state)
        next_state = state
    end

    return next_state, 0, false
end


text:
# ArgMacros
# macros.jl

"Ensure flags is not empty and that no flag overlaps with those used by help"
function _validateflags(local_name, flags)
    @assert !isempty(flags) "Argument $local_name must have at least one flag."
    @assert !("-h" in flags) && !("--help" in flags) "Flags -h and --help are reserved"
end

"""
    @inlinearguments begin ... end

Denote and setup a block with other macros from `ArgMacros`

# Example
```julia
function julia_main()
    @inlinearguments begin
        ...
        @argumentrequired Int foo "-f" "--foo"
        @argumentdefault Int 5 bar "-b" "--bar"
        ...
    end
    ...
end
```
"""
macro inlinearguments(block::Expr)
    _validateorder(block) # Validate ordering of macros

    return quote
        splitargs = _split_arguments($(esc(:ARGS))) # Consumed by the macros in the block
        allowextraargs = false

        _help_check(splitargs, $(Expr(:quote, block))) # End and print help if needed

        $block

        if !allowextraargs && !isempty(splitargs)
            _quit_try_help("Too many or unrecognized arguments provided")
        end
    end
end

"""
    @beginarguments begin ... end

This macro is deprecated beginning in ArgMacros v0.2.0
Please use @inlinearguments, which has the same interface
Or consider @structarguments, @tuplearguments, and @dictarguments depending on your use case
"""
macro beginarguments(block::Expr)
    @warn """
        Macro `@beginarguments` is deprecated in ArgMacros v0.2.0
        Please use @inlinearguments, which has the same interface
        Or consider @structarguments, @tuplearguments, and @dictarguments depending on your use case
    """
    esc(:(@inlinearguments $block))
end

"""
    @structarguments mutable typename begin ... end

Denote and setup a block with other macros from `ArgMacros`
Defines an optionally mutable struct type based on the arguments and a zero-argument constructor
which will generate an instance of the struct based on the parsed arguments.

# Example
```julia
@structarguments false Args begin
    ...
    @argumentrequired Int foo "-f" "--foo"
    @argumentdefault Int 5 bar "-b" "--bar"
    ...
end

function julia_main()
    args = Args()
    ...
end
```
"""
macro structarguments(mutable::Bool, name::Symbol, block::Expr)
    Expr(:block,
        Expr(:struct, mutable, name, Expr(:block, _getargumentpairs(block)...)),
        esc(Expr(:function,
            :($name()::$name),
            Expr(:block,
                :(@inlinearguments $block),
                Expr(:call,
                    name,
                    (:($(pair.args[1])) for pair in _getargumentpairs(block))...
                )
            )
        ))
    )
end

"""
    @tuplearguments begin ... end

Denote and setup a block with other macros from `ArgMacros`
Return a NamedTuple with the arguments instead of dumping them in the enclosing namespace

# Example
```julia
function julia_main()
    args = @tuplearguments begin
        ...
        @argumentrequired Int foo "-f" "--foo"
        @argumentdefault Int 5 bar "-b" "--bar"
        ...
    end
    ...
end
```
"""
macro tuplearguments(block::Expr)
    Expr(:let, Expr(:block), esc(Expr(:block,
        :(@inlinearguments $block),
        Expr(:tuple,
            (:($(pair.args[1]) = $(pair.args[1])) for pair in _getargumentpairs(block))...
        )
    )))
end

"""
    @dictarguments begin ... end

Denote and setup a block with other macros from `ArgMacros`
Return a Dict with the arguments instead of dumping them in the enclosing namespace

# Example
```julia
function julia_main()
    args = @dictarguments begin
        ...
        @argumentrequired Int foo "-f" "--foo"
        @argumentdefault Int 5 bar "-b" "--bar"
        ...
    end
    ...
end
```
"""
macro dictarguments(block::Expr)
    Expr(:let, Expr(:block), esc(Expr(:block,
        :(@inlinearguments $block),
        Expr(:call,
            :(Dict{Symbol, Any}),
            (:($(Meta.quot(pair.args[1])) => $(pair.args[1])) for pair in _getargumentpairs(block))...
        )
    )))
end

#=
Remaining macros require two escapes to get to the level of the calling scope
Because they will be nested in the @beginarguments macro
=#

"""
    @argumentrequired type local_name flags::String...

Get a required argument specified by the given flags and store
in the variable `local_name` with the specified type.

Print a message and exit the program if the value is not found or
cannot be converted to the specified type.

Must be used in `@beginarguments begin ... end` block

# Example
```julia
@beginarguments begin
    ...
    @argumentrequired String output_file "-o" "--output"
    ...
end
```
"""
macro argumentrequired(type::Symbol, local_name::Symbol, flags::String...)
    _validateflags(local_name, flags)
    return quote
        #  _converttype! is completely type safe
        local $(esc(esc(local_name)))::$(esc(esc(type))) = _converttype!($type, _pop_argval!(splitargs, [$flags...]), $(flags[end]))
    end
end

"""
    @argumentdefault type default_value local_name flags::String...

Attempt to get an argument specified by the given flags and store
in the variable `local_name` with the specified type.
Store the default value instead if the flags cannot be found.
Default value automatically converted to specified type.

Must be used in `@beginarguments begin ... end` block

# Example
```julia
@beginarguments begin
    ...
    @argumentdefault String "output.txt" output_file "-o" "--output"
    ...
end
```
"""
macro argumentdefault(type::Symbol, default_value, local_name::Symbol, flags::String...)
    _validateflags(local_name, flags)
    return quote
        potential_val::Union{String, Nothing} = _pop_argval!(splitargs, [$flags...])
        # Convert either potential or the default value, allows default to be specified with wrong literal type
        local $(esc(esc(local_name)))::$(esc(esc(type))) = _converttype!($type, something(potential_val, $default_value), $(flags[end]))
    end
end

"""
    @argumentoptional type local_name flags::String...

Attempt to get an argument specified by the given flags and store
in the variable `local_name` with the type `Union{type, Nothing}`.
Store `nothing` if the flags cannot be found.

Must be used in `@beginarguments begin ... end` block

# Example
```julia
@beginarguments begin
    ...
    @argumentoptional String output_file "-o" "--output"
    ...
end
```
"""
macro argumentoptional(type::Symbol, local_name::Symbol, flags::String...)
    _validateflags(local_name, flags)
    return quote
        potential_val::Union{String, Nothing} = _pop_argval!(splitargs, [$flags...])
        # Return nothing directly without calling _converttype! if arg not found
        local $(esc(esc(local_name)))::$(esc(esc(Union))){$(esc(esc(type))), $(esc(esc(Nothing)))} =
            isnothing(potential_val) ? nothing : _converttype!($type, potential_val, $(flags[end]))
    end
end

"""
    @argumentflag local_name flags::String...

Store `true` in the variable `local_name` with type `Bool` if one or more of the flags is found.
Otherwise, store `false`.

Must be used in `@beginarguments begin ... end` block

# Example
```julia
@beginarguments begin
    ...
    @argumentflag verbose "-v" "--verbose"
    ...
end
```
"""
macro argumentflag(local_name::Symbol, flags::String...)
    _validateflags(local_name, flags)
    return quote
        local $(esc(esc(local_name)))::$(esc(esc(Bool))) = _pop_flag!(splitargs, [$flags...])
    end
end

"""
    @argumentcount local_name flag::String

Store the number of occurrences of `flag` in `local_name` with type `Int`.

Must be used in `@beginarguments begin ... end` block

# Example
```julia
@beginarguments begin
    ...
    @argumentcount verbose "-v"
    ...
end
```
"""
macro argumentcount(local_name::Symbol, flag::String)
    return quote
        local $(esc(esc(local_name)))::$(esc(esc(Int))) = _pop_count!(splitargs, $flag)
    end
end

"""
    @positionalrequired type local_name [help_name::String]

Attempt to get a positional argument and store in variable `local_name` with the specified type.

Positional arguments are read in order after all flag/option arguments have been read.
`help_name` used instead of `local_name` in messages to user if specified.

Print a message and exit the program if a value is not found or cannot be converted to the specified type.

Must be used in `@beginarguments begin ... end` block

# Example
```julia
@beginarguments begin
    ...
    @positionalrequired String output_file "output"
    ...
end
```
"""
macro positionalrequired(type::Symbol, local_name::Symbol, help_name::Union{String, Nothing}=nothing)
    help_name_str::String = something(help_name, String(local_name))
    return quote
        local $(esc(esc(local_name)))::$(esc(esc(type))) = _converttype!(
            $type,
            !isempty(splitargs) ? popfirst!(splitargs) : nothing,
            $help_name_str
        )
    end
end

"""
    @positionaldefault type default_value local_name [help_name::String]

Attempt to get a positional argument and store in variable local_name with the specified type.
Store the default value instead if an argument cannot be found.

Positional arguments are read in order after all flag/option arguments have been read.
`help_name` used instead of `local_name` in messages to user if specified.

Default value automatically converted to specified type.

Must be used in `@beginarguments begin ... end` block

# Example
```julia
@beginarguments begin
    ...
    @positionaldefault String output_file "output"
    ...
end
```
"""
macro positionaldefault(type::Symbol, default_value, local_name::Symbol, help_name::Union{String, Nothing}=nothing)
    help_name_str::String = something(help_name, String(local_name))
    return quote
        local $(esc(esc(local_name)))::$(esc(esc(type))) =_converttype!(
            $type,
            !isempty(splitargs) ? popfirst!(splitargs) : $default_value,
            $help_name_str
        )
    end
end

"""
    @positionaloptional type local_name [help_name]

Attempt to get a positional argument and store in variable `local_name` with the type
`Union{type, Nothing}`. Store `nothing` if an argument is not found.

Positional arguments are read in order after all flag/option arguments have been read.
`help_name` used instead of `local_name` in messages to user if specified.

Must be used in `@beginarguments begin ... end` block

# Example
```julia
@beginarguments begin
    ...
    @positionaloptional String output_file "output"
    ...
end
```
"""
macro positionaloptional(type::Symbol, local_name::Symbol, help_name::Union{String, Nothing}=nothing)
    help_name_str::String = something(help_name, String(local_name))
    return quote
        local $(esc(esc(local_name)))::$(esc(esc(Union))){$(esc(esc(type))), $(esc(esc(Nothing)))} =
            !isempty(splitargs) ? _converttype!($type, popfirst!(splitargs), $help_name_str) : nothing
    end
end

"""
    @positionalleftover type local_name [help_name]

Get any leftover positional arguments after all other arguments have been parsed, and store in
variable `local_name` with the type `Vector{type}`. This vector will be empty if there are no
leftover arguments.

This macro must be used after all other arguments have been declared, as positional arguments are
read in order after all flag/option arguments have been read.
Using this macro means all input will be captured by this argument if not by another, so
`@allowextraarguments` is not necesesary when this macro is used.
`help_name` used instead of `local_name` in messages to user if specified.

Must be used in `@beginarguments begin ... end` block

# Example
```julia
@beginarguments begin
    ...
    @positionalleftover String file_names "files"
    ...
end
```
"""
macro positionalleftover(type::Symbol, local_name::Symbol, help_name::Union{String, Nothing}=nothing)
    help_name_str::String = something(help_name, String(local_name))
    return quote
        local $(esc(esc(local_name)))::$(esc(esc(Vector))){$(esc(esc(type)))} =
            _converttype!($type, splitargs, $help_name_str)
        empty!(splitargs)
    end
end

"""
    @argtest argname func [desc]

Apply `func` to the value stored in `argname`, printing an error message (optionally
specified by `desc`) and the program  if `func` returns `false`.
Test skipped if `argname` has value `nothing` (only possible for optional arguments).
This macro must be used AFTER declaring the arugment with another macro.

Must be used in `@beginarguments begin ... end` block

# Example
```julia
@beginarguments begin
    ...
    @positionalrequired String input_file
    @argtest input_file isfile "Couldn't find the input file."
    ...
end
```
"""
macro argtest(argname::Symbol, func::Union{Symbol, Expr}, desc::Union{String, Nothing}=nothing)
    errstr::String = something(desc, "Tests for argument $argname failed.")
    return quote
        if !(isnothing($(esc(esc(argname)))) || $(esc(esc(func)))($(esc(esc(argname)))))
            _quit_try_help($errstr)
        end
    end
end

"""
    @allowextraarguments

Disables the default behavior of printing a message and exiting
the program when not all values in `ARGS` could be assigned to
specified arguments.

This will not capture any arguments, use `@positionalleftover`
if you want to capture extra arguments that could not otherwise
be assigned.

Must be used in `@beginarguments begin ... end` block

# Example
```julia
@beginarguments begin
    ...
    @allowextraarguments
end
```
"""
macro allowextraarguments()
    return quote
        allowextraargs = true
    end
end


text:
# =========================================================================== #
# Test
# =========================================================================== #

function test_norm()

    a = [10.0 11; 5 2]
    println(a)

    println( linalg.norm_L1(a) )
    println( linalg.norm_Linf(a) )

    println( linalg.norm(a, :1) )
    println( linalg.norm(a, :inf) )

end

# =========================================================================== #
# Benchmark
# =========================================================================== #

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
"""$(SIGNATURES) Element type of a `rv::RandVar`

The element type of a RandVar is the type of value you get when you sample from it.
from it.  Since `RandVar`s are normal julia functions, this may not be a single concrete type.

Warning.  `elemtype` relies on Julia's type inference.
We use it mostly for printing and debugging.
Relying on it as part of an algorithm is inadvisable.
"""
function elemtype(rv::RandVar, ΩT::Type{OT} = defΩ(rv)) where OT
  Base.promote_op(apl, typeof(rv), OT)
end
@spec rand(rv) isa _res


text:
struct TapeFiles{T1 <: Nullable{IOStream}, T2 <: Nullable{IOStream}}
    ts::IOStream
    rs::IOStream
    vs::T1
    ps::T2
end

function TapeFiles(ts_file, rs_file, vs_file = nothing, ps_file = nothing)
    pos = position(rs_file)
    @assert isnothing(vs_file) || pos ≡ position(vs_file)
    @assert isnothing(ps_file) || pos ≡ position(ps_file)
    T1 = isnothing(vs_file) ? Nothing : IOStream
    T2 = isnothing(ps_file) ? Nothing : IOStream
    TapeFiles{T1, T2}(ts_file, rs_file, vs_file, ps_file)
end

function TapeFiles(
    ts_filename::AbstractString,
    rs_filename::AbstractString,
    vs_filename::Nullable{AbstractString} = nothing,
    ps_filename::Nullable{AbstractString} = nothing;
    read = true,
)
    files = [
        if isnothing(filename) || read && !isfile(filename)
            nothing
        else
            open(filename, read ? "r" : "a+")
        end for filename in (ts_filename, rs_filename, vs_filename, ps_filename)
    ]
    TapeFiles(files...)
end

TapeFiles(datadir::AbstractString; filename_prefix = "tape", has_vs = true, has_ps = false, read = true) =
    TapeFiles(
        joinpath(datadir, "$filename_prefix-ts.dat"),
        joinpath(datadir, "$filename_prefix-rs.dat"),
        has_vs ? joinpath(datadir, "$filename_prefix-vs.dat") : nothing,
        has_ps ? joinpath(datadir, "$filename_prefix-ps.dat") : nothing;
        read = read,
    )

has_vs(tf::TapeFiles) = !isnothing(tf.vs)
has_ps(tf::TapeFiles) = !isnothing(tf.ps)

function read_tape(tf::TapeFiles, ::Type{T}, N::Integer) where {T <: MosiVector}
    times = read_vector_all(tf.ts, Float64)
    rs = read_vectors_all(tf.rs, T, N)
    vs = has_vs(tf) ? read_vectors_all(tf.vs, T, N) : Vector{T}[]
    periods = has_ps(tf) ? read_vector3s_all(tf.ps, T, N) : Vector{T}[]
    SimulationTape(times, rs, vs, periods)
end

function update_tape(tf::TapeFiles, s::SimulationState)
    write(tf.ts, time(s))
    write(tf.rs, positions(s))
    if has_vs(tf)
        write(tf.vs, velocities(s))
    end
    if has_ps(tf)
        write(tf.ps, periods(s))
    end
end
all_files(tf::TapeFiles) =
    let a = has_vs(tf), b = has_ps(tf)
        a ? b ? (tf.ts, tf.rs, tf.vs, tf.ps) : (tf.ts, tf.rs, tf.vs) :
        b ? (tf.ts, tf.rs, tf.ps) : (tf.ts, tf.rs)
    end
Base.close(tf::TapeFiles) = close.(all_files(tf))
Base.flush(tf::TapeFiles) = flush.(all_files(tf))
Base.seekend(tf::TapeFiles) = seekend.(all_files(tf))
Base.seekstart(tf::TapeFiles) = seekstart.(all_files(tf))

struct MultiFileMemoryMapTape{TM <: MosiModel, T1 <: Nullable{IOStream}, T2 <: Nullable{IOStream}} <:
       SimulationTape{TM}
    tape_files::TapeFiles{T1, T2}
    len::Int
    model::TM
end

function MultiFileMemoryMapTape(
    datadir::AbstractString,
    model;
    filename_prefix = "tape",
    has_vs = true,
    has_ps = false,
)
    tape_files = TapeFiles(datadir; filename_prefix = filename_prefix, has_vs = has_vs, has_ps = has_ps)
    len = filesize(joinpath(datadir, "$filename_prefix-ts.dat")) ÷ sizeof(Float64)
    MultiFileMemoryMapTape(tape_files, len, model)
end

function fetch_data!(file::IOStream, ::Type{T}, len::Integer, i::Integer) where {T}
    record_size = sizeof(T) * len
    offset = record_size * (i - 1)
    seek(file, offset)
    vec = Vector{T}(undef, len)
    read!(file, vec)
    vec
end

has_ps(tape::MultiFileMemoryMapTape) = has_ps(tape.tape_files)
Base.length(tape::MultiFileMemoryMapTape) = tape.len
natoms(tape::MultiFileMemoryMapTape) = natoms(tape.model)
vectype(tape::MultiFileMemoryMapTape) = vectype(tape.model)
Base.time(tape::MultiFileMemoryMapTape, i) = fetch_data!(tape.tape_files.ts, Float64, 1, i)[1]
times(tape::MultiFileMemoryMapTape) = read_vector_all(tape.tape_files.ts, Float64)
positions(tape::MultiFileMemoryMapTape, i) =
    fetch_data!(tape.tape_files.rs, vectype(tape), natoms(tape), i)
positions(tape::MultiFileMemoryMapTape) =
    read_vectors_all(tape.tape_files.rs, vectype(tape), natoms(tape))
velocities(tape::MultiFileMemoryMapTape, i) =
    fetch_data!(tape.tape_files.vs, vectype(tape), natoms(tape), i)
velocities(tape::MultiFileMemoryMapTape) =
    if has_vs(tape)
        read_vectors_all(tape.tape_files.vs, vectype(tape), natoms(tape))
    else
        Vector{vectype(tape)}[]
    end
periods(tape::MultiFileMemoryMapTape, i) =
    if has_ps(tape)
        fetch_data!(tape.tape_files.ps, vectype(tape), natoms(tape), i)
    else
        vectype(tape)[]
    end
periods(tape::MultiFileMemoryMapTape) =
    if has_ps(tape)
        read_vectors_all(tape.tape_files.ps, vectype(tape), natoms(tape))
    else
        Vector{vectype(tape)}[]
    end
configuration(tape::MultiFileMemoryMapTape, i) =
    ConfigurationSystem(positions(tape, i), periods(tape, i), pbc_box(tape.model))

function get_configuration_func(
    tape::MultiFileMemoryMapTape;
    atom_range = 1:natoms(tape),
    buffer_size = 256,
    return_copy = false,
)
    model = tape.model
    T = vectype(model)
    hasps = has_ps(tape)
    box = pbc_box(model)
    N = natoms(tape)
    a, b = first(atom_range), last(atom_range)
    @assert a > 0 && b ≤ N
    s1 = sizeof(T)
    rsize = s1 * N
    a, b = (a - 1) * s1, b * s1
    buffer_len = buffer_size * rsize
    buffer_rs = zeros(UInt8, buffer_len)
    buffer_ps = hasps ? zeros(UInt8, buffer_len) : nothing
    range = 1:buffer_size
    f_rs = tape.tape_files.rs
    f_ps = tape.tape_files.ps
    seekstart(f_rs)
    readbytes!(f_rs, buffer_rs; all = false)
    if hasps
        seekstart(f_ps)
        readbytes!(f_ps, buffer_ps; all = false)
    end
    (i::Int) -> begin
        if i ∉ range
            t = (i - 1) ÷ buffer_size
            range = t*buffer_size+1:(t+1)*buffer_size
            seek(f_rs, rsize * t * buffer_size)
            readbytes!(f_rs, buffer_rs; all = false)
            if hasps
                seek(f_ps, rsize * t * buffer_size)
                readbytes!(f_ps, buffer_ps; all = false)
            end
        end
        j = (i - 1) % buffer_size
        i_start = j * rsize
        rs = reinterpret(T, view(buffer_rs, i_start+1+a:i_start+b))
        if return_copy
            rs = copy(rs)
        end
        if hasps
            ps = reinterpret(T, view(buffer_ps, i_start+1+a:i_start+b))
            if return_copy
                ps = copy(ps)
            end
            ConfigurationSystem((rs, ps), box = box, update_periods = false)
        else
            ConfigurationSystem(rs, box = box, update_periods = false)
        end
    end
end


text:
module FlowNodeLayersInputLayerTest

using Test
using ReactiveMP
using LinearAlgebra

using ReactiveMP: getdim

@testset "Input Layer" begin

    @testset "Constructor" begin

        # check for placeholder creation
        layer = InputLayer(2)
        @test typeof(layer)   == InputLayer
        @test typeof(layer)   <: ReactiveMP.AbstractLayerPlaceholder

        # check error handling
        @test_throws AssertionError InputLayer(1)
        @test_throws AssertionError InputLayer(-1)

    end

    @testset "Get" begin

        # check for layer
        for k = 2:10
            layer = InputLayer(k)
            @test layer.dim     == getdim(layer)
            @test getdim(layer) == k
        end

    end

    @testset "Integration in flow model" begin

        model = FlowModel((InputLayer(5), AdditiveCouplingLayer(PlanarFlow()), AdditiveCouplingLayer(PlanarFlow(); permute=false)))
        compiled_model = compile(model)
        @test typeof(first(getlayers(compiled_model))) <: AdditiveCouplingLayer
        for layer in getlayers(compiled_model)
            @test getdim(layer) == 5
        end

    end

end

end


text:
#
# utils.jl -
#
# General purpose methods.
#
#-------------------------------------------------------------------------------
#
# This file is part of LazyAlgebra (https://github.com/emmt/LazyAlgebra.jl)
# released under the MIT "Expat" license.
#
# Copyright (c) 2017-2020 Éric Thiébaut.
#

@noinline bad_argument(args...) = bad_argument(string(args...))
bad_argument(mesg::String) = throw(ArgumentError(mesg))

@noinline bad_size(args...) = bad_size(string(args...))
bad_size(mesg::String) = throw(DimensionMismatch(mesg))

arguments_have_incompatible_axes() =
    bad_size("arguments have incompatible dimensions/indices")

operands_have_incompatible_axes() =
    bad_size("operands have incompatible dimensions/indices")

"""
    message([io=stdout,] header, args...; color=:blue)

prints a message on `io` with `header` text in bold followed by a space,
`args...` and a newline.  Keyword `color` can be used to specify the text color
of the message.

"""
message(header::String, args...; kwds...) =
    message(stdout, header, args...; kwds...)

@noinline function message(io::IO, header::String, args...;
                           color::Symbol=:blue)
    printstyled(io, header; color=color, bold=true)
    printstyled(io, " ", args...; color=color, bold=false)
    println(io)
end

"""
    warn([io=stdout,] args...)

prints a warning message in yellow on `io` with `"Warning: "` in bold followed
by `args...` and a newline.

"""
warn(args...) = warn(stderr, args...)
warn(io::IO, args...) = message(io, "Warning:", args...; color=:yellow)

#inform(args...) = inform(stderr, args...)
#inform(io::IO, args...) = message(io, "Info:", args...; color=:blue)

"""
    promote_multiplier(λ, T)

yields multiplier `λ` converted to a suitable floating-point type for
multiplying values or expressions of type `T`.  This method is *type stable*.
The result has the same floating-point precision as `T` and is a real if `λ` is
real or a complex if `λ` is complex.

Multiple arguments can be specified after the multiplier `λ`:

    promote_multiplier(λ, args...)

to have `T` the promoted type of all types in `args...` or all element types of
arrays in `args...`.

See methods [`LazyAlgebra.multiplier_type`](@ref) and
[`LazyAlgebra.multiplier_floatingpoint_type`](@ref).

""" promote_multiplier

# Note taht the only direct sub-types of `Number` are abstract types `Real` and
# `Complex`.  Also see discussion here
# (https://github.com/emmt/LinearInterpolators.jl/issues/7) for details about
# the following implementation.

@inline function promote_multiplier(λ::Real, args...)
    T = multiplier_floatingpoint_type(args...)
    return convert(T, λ)::T
end

@inline function promote_multiplier(λ::Complex{<:Real}, args...)
    T = multiplier_floatingpoint_type(args...)
    return convert(Complex{T}, λ)::Complex{T}
end

"""
    multiplier_floatingpoint_type(args...) -> T::AbstractFloat

yields the multiplier floating-point type for the arguments `args...` of the
multiplier.  Each argument may be anything acceptable for
[`LazyAlgebra.multiplier_type`](@ref).  The result is guaranteed to be a
concrete floating-point type.

See methods [`LazyAlgebra.promote_multiplier`](@ref) and
[`LazyAlgebra.multiplier_type`](@ref).

""" multiplier_floatingpoint_type

multiplier_floatingpoint_type(::Tuple{}) =
    throw(ArgumentError("at least one other argument must be specified"))

@inline function multiplier_floatingpoint_type(args...)
    T = promote_type(map(multiplier_type, args)...)
    (T <: Number && isconcretetype(T)) || error(
        "resulting multiplier type ", T, " is not a concrete real type")
    return float(real(T))
end

"""
    multiplier_type(x) -> T::Number

yields the *element* type to be imposed to multipliers of `x`.  The result must
be a concrete number type.  Argument `x` may be an array, a number, or a data
type.  Other packages are however encouraged to specialize this method for
their needs.

See methods [`LazyAlgebra.promote_multiplier`](@ref) and
[`LazyAlgebra.multiplier_floatingpoint_type`](@ref).

"""
multiplier_type(::Type{T}) where  {T<:Number} = T
multiplier_type(::AbstractArray{T}) where {T<:Number} = T
multiplier_type(::T) where  {T<:Number} = T

"""
    to_tuple(arg)

converts `arg` into an `N`-tuple where `N` is the number of elements of `arg`.
This is equivalent to `Tuple(arg)` or `(arg...,)` for a vector but it is much
faster for small vectors.

""" to_tuple

to_tuple(x::Tuple) = x

# The cutoff at n = 10 below reflects what is used by `ntuple`.  This value is
# somewhat arbitrary, on the machines where I tested the code, the explicit
# unrolled expression for n = 10 is still about 44 times faster than `(x...,)`.
# Calling `ntuple` for n ≤ 10 is about twice slower; for n > 10, `ntuple` is
# slower than `(x...,)`.
function to_tuple(x::AbstractVector)
    n = length(x)
    @inbounds begin
        n == 0 ? () :
        n > 10 || firstindex(x) != 1 ? (x...,) :
        n == 1 ? (x[1],) :
        n == 2 ? (x[1], x[2]) :
        n == 3 ? (x[1], x[2], x[3]) :
        n == 4 ? (x[1], x[2], x[3], x[4]) :
        n == 5 ? (x[1], x[2], x[3], x[4], x[5]) :
        n == 6 ? (x[1], x[2], x[3], x[4], x[5], x[6]) :
        n == 7 ? (x[1], x[2], x[3], x[4], x[5], x[6], x[7]) :
        n == 8 ? (x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8]) :
        n == 9 ? (x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9]) :
        (x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10])
    end
end

"""
    @certify expr [mesg]

asserts that expression `expr` is true; otherwise, throws an `AssertionError`
exception with message `mesg`.  If unspecified, `mesg` is `expr` converted into
a string.  Compared to `@assert`, the assertion made by `@certify` may never be
disabled whatever the optimization level.

"""
macro certify(expr)
    _certify(expr, string(expr))
end
macro certify(expr, mesg::Union{Expr,Symbol})
    _certify(expr, :(string($(esc(mesg)))))
end
macro certify(expr, mesg::AbstractString)
    _certify(expr, mesg)
end
macro certify(expr, mesg)
    _certify(expr, string(mesg))
end
_certify(expr, mesg) = :($(esc(expr)) ? nothing : throw(AssertionError($mesg)))


text:
# Functions:
function fft_len(x::AbstractVector, len::Int)
    if len > length(x)
        return fft([x; zeros(len - length(x))])
    else
        return fft(x[1:len])
    end
end

function conv_same(x::AbstractVector, y::AbstractVector)
    # RoundUp is used to handle scenarios where length(y) is odd
    s = div(length(y), 2, RoundUp)
    e = length(x) + s - 1
    return conv(x, y)[s:e]
end

####################################################################

# Julia code for Example 10.5
using Plots

x = zeros(1000, 20)
t = range(-2, 2, length=1000)

for col in eachcol(x)
    col .= rand() * cos.(2π * t)
end

plot(t, x, linewidth=2, color=:gray80, legend=false)
plot!(t, 0.5 * cos.(2π * t), linewidth=4, color=:darkred)

####################################################################

# Julia code for Example 10.6
using Plots

x = zeros(1000, 20)
t = range(-2, 2, length=1000)

for col in eachcol(x)
    col .= cos.(2π*t .+ 2π*rand())
end

plot(t, x, linewidth=2, color=:gray80, legend=false)
hline!([0], linewidth=4, color=:darkred)

####################################################################

# Julia code for Example 10.7
using Plots

x = zeros(21, 20)
t = 0:20

for col in eachcol(x)
    col .= rand() .^ t
end

sticks(t, x, linewidth=2, shape=:circle, color=:gray80, legend=false)
sticks!(t, 1 ./ (t .+ 1), linewidth=2, shape=:circle, markersize=6, color=:darkred)

####################################################################

# Julia code for Example 10.11: Plot the time function
using Plots

x = zeros(1000, 20)
t = range(-2, 2, length=1000)

for col in eachcol(x)
    col .= rand() * cos.(2π * t)
end

plot(t, x, linewidth=2, color=:gray80, legend=false)
plot!(t, 0.5 * cos.(2π * t), linewidth=4, color=:darkred)
scatter!(zeros(20), x[501, :], markersize=4, color=:orange)
scatter!(0.5 * ones(20), x[626, :], markersize=4, color=:blue)

# Julia code for Example 10.11: Plot the autocorrelation function
t = range(-1, 1, length=1000)
R = 1/3 * cos.(2π * t) * cos.(2π * t')

heatmap(t, t, R, yticks= -1:0.25:1, xticks= -1:0.25:1, legend=false, yflip=true)

####################################################################

# Julia code for Example 10.12: Plot the time function
using ToeplitzMatrices
using Plots

x = zeros(1000, 20)
t = range(-2, 2, length=1000)

for col in eachcol(x)
    col .= cos.(2π*t .+ 2π*rand())
end

plot(t, x, linewidth=2, color=:gray80, legend=false)
plot!(t, 0 * cos.(2π * t), linewidth=4, color=:darkred)
scatter!(zeros(20), x[501, :], markersize=4, color=:orange)
scatter!(0.5 * ones(20), x[626, :], markersize=4, color=:blue)

# Julia code for Example 10.12: Plot the autocorrelation function
t = range(-1, 1, length=1000)
R = SymmetricToeplitz(0.5 * cos.(2π * t))

heatmap(t, t, R, yticks= -1:0.25:1, xticks= -1:0.25:1, legend=false, yflip=true)

####################################################################

# Julia code to demonstrate autocorrelation
using Plots
using DSP

# Figure 1
X₁ = randn(1000)
X₂ = randn(1000)

plot(X₁, linewidth=2, linestyle=:dot, color=:blue, legend=false)
plot!(X₂, linewidth=2, linestyle=:dot, color=:black)

# Figure 2
N  = 1000 # number of sample paths
T  = 1000 # number of time stamps
X  = randn(N, T)
xc = zeros(N, 2T - 1)

for i in 1:N
    xc[i, :] .= xcorr(X[i, :], X[i, :]) / T
end

plot(xc[1, :], linewidth=2, linestyle=:dot, color=:blue, label="correlation of sample 1")
plot!(xc[2, :], linewidth=2, linestyle=:dot, color=:black, label="correlation of sample 2")
####################################################################

# Julia code for Example 10.15
using Plots, LaTeXStrings
using Images
using DSP

# Figure 1
t = -10:0.001:10
L = length(t)
X = randn(L)
h = 10 * max.(1 .- abs.(t), 0) / 1000
Y = imfilter(X, h, "circular")

plot(t, X, linewidth=0.5, color=:gray60, label=L"X(t)")
hline!([0], linewidth=4, color=:gray20, label=L"\mu_X(t)")
plot!(t, Y, linewidth=2, color=:red, label=L"Y(t)")
hline!([0], linewidth=4, linestyle=:dot, color=:yellow, label=L"\mu_Y(t)")
xlims!((-10, 10))

# Figure 2
h₂ = conv(h, h)
Rₓ = zeros(40001)
Rₓ[20001] = 0.2

plot(-20:0.001:20, Rₓ, linewidth=2, color=:gray60, label=L"R_X(t)")
plot!(-20:0.001:20, h₂, linewidth=2, color=:red, label=L"R_Y(t)")
xlims!((-2, 2))
ylims!((-0.05, 0.2))

####################################################################

# Julia code to solve the Yule Walker Equation
using Plots, LaTeXStrings
using ToeplitzMatrices
using DelimitedFiles
using DSP

y = vec(readdlm("ch10_LPC_data.txt"))
K = 10
N = length(y)
y_corr = xcorr(y, y)
R = SymmetricToeplitz(y_corr[N:N+K-1])
lhs = y_corr[N+1:N+K]
h = R \ lhs

# Figure 1
plot(y, linewidth=4, color=:blue, label=L"Y[n]")

# Figure 2
plot(y_corr, linewidth=4, color=:black, label=L"R_Y[k]")

####################################################################

# Julia code to predict the samples
using ToeplitzMatrices
using DelimitedFiles
using Plots
using DSP

y = vec(readdlm("ch10_LPC_data_02.txt"))
K = 10
N = length(y)

y_corr = xcorr(y, y)
R = SymmetricToeplitz(y_corr[N:N+K-1])
lhs = y_corr[N+1:N+K]
h = R \ lhs

z = y[311:320]
ŷ = zeros(340)
ŷ[1:320] .= y

for t in 1:20
    pred = z' * h
    z = [z[2:10]; pred]
    ŷ[320+t] = pred
end

plot(ŷ, linewidth=3, color=:red, label="Prediction", legend=:bottomleft)
plot!(y, linewidth=4, linestyle=:dash, color=:gray60, label="Input")

####################################################################

# Julia code for Wiener filtering
using DSP, FFTW
using Plots

w = 0.05 * randn(320)
x = y + w

Ry = xcorr(y, y)
Rw = xcorr(w, w)
Sy = fft(Ry)
Sw = fft(Rw)
H = Sy ./ (Sy + Sw)

Ŷ = H .* fft_len(x, 639)
ŷ = real.(ifft(Ŷ))

# Figure 1
plot(x, linewidth=4, color=:gray, label="Noise Input X[n]", legend=:bottomleft)
plot!(ŷ[1:320], linewidth=2, color=:red, label="Wiener Filtered Yhat[n]")
plot!(y, linewidth=2, linestyle=:dot, color=:black, label="Ground Truth Y[n]")

# Figure 2
plot(Rw, linewidth=4, color=:blue, label="h[n]")
####################################################################

# Julia code to solve the Wiener deconvolution problem
using DelimitedFiles
using DSP, FFTW
using Plots

y = vec(readdlm("ch10_wiener_deblur_data.txt"))
g = ones(32) / 32
w = 0.02 * randn(320)
x = conv_same(y, g) + w

Ry = xcorr(y, y)
Rw = xcorr(w, w)
Sy = fft(Ry)
Sw = fft(Rw)
G  = fft_len(g, 639)

H = @. (conj(G) * Sy) / (abs(G)^2 * Sy + Sw)
Ŷ = H .* fft_len(x, 639)
ŷ = real.(ifft(Ŷ))

plot(x, linewidth=4, color=:gray, label="Noise Input X[n]")
plot!(16:320+15, ŷ[1:320], linewidth=2, color=:red, label="Wiener Filtered Yhat[n]")
plot!(y, linewidth=2, linestyle=:dot, color=:black, label="Ground Truth Y[n]")


text:
# This file was generated, do not modify it. # hide
X_iris, y_iris = @load_iris;
schema(X_iris)

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
# +begin_src jupyter-julia
using LightGraphs
using GraphPlot
using GraphRecipes
using Plots
using Printf
using LaTeXStrings
println("Loading complete")
# +end_src

# +begin_src jupyter-julia
push!(LOAD_PATH, "../dev")
using Dice
# +end_src

# +begin_src jupyter-julia
G8 = SimpleGraph(8)
add_edge!(G8, 1, 2)
add_edge!(G8, 2, 3)
add_edge!(G8, 3, 4)
add_edge!(G8, 4, 5)
add_edge!(G8, 5, 6)
add_edge!(G8, 6, 7)
add_edge!(G8, 7, 8)
add_edge!(G8, 1, 8)
add_edge!(G8, 1, 5)
add_edge!(G8, 2, 6)
add_edge!(G8, 3, 7)
add_edge!(G8, 4, 8)

p = graphplot(G8, nodesize=0.05, nodeshape=:circle, curves=false, nodecolor=:black)
# +end_src

# +begin_src jupyter-julia
# The parameters below may differ from those used for the paper,
# therefore, the exact profile may differ
G = G8
mtds = sine
scale = 0.15 / nv(G)
tmax = 350

Ksstart = 0
Ksmax = 7
Ksstep = 0.1
Ksran = collect(Ksstart:Ksstep:Ksmax)
MM = length(Ksran)

model = Model(G, mtds, scale, Ksstart)

Nrep = 10000
probs = zeros(MM)
for j in 1:MM
    model.Ks = Ksran[j]
    count_best = 0

    for i in 1:Nrep
        V3 = get_initial(nv(G), (-0.15, 0.15));
        V3t = propagate(model, tmax, V3);
        state = roundup(V3t[:, end])
        curcut = cut(G, extract_configuration(state, 0))
        if curcut == 10
            count_best += 1
        end
    end
    println("For Ks = $(model.Ks):$Ksmax - success $count_best out of $Nrep")
    probs[j] = count_best / Nrep
end
plot(Ksran, probs)
# +end_src

max_stars_repo_name: merement/Dice

max_issues_repo_name: merement/Dice

max_forks_repo_name: merement/Dice


text:
module AstroAngles

export deg2dms,
       deg2hms,
       deg2ha,
       ha2dms,
       ha2hms,
       ha2deg,
       ha2rad,
       rad2dms,
       rad2hms,
       rad2ha,
       hms2deg,
       hms2ha,
       hms2rad,
       dms2deg,
       dms2rad,
       dms2ha,
       parse_dms,
       parse_hms,
       @dms_str,
       @hms_str,
       format_angle

include("conversions.jl")
include("parsing.jl")
include("printing.jl")

end


text:
abstract Separator

doc"""
ConstraintSeparator is a separator that represents a constraint defined directly
using `@constraint`.
"""
# CHANGE TO IMMUTABLE!
type ConstraintSeparator <: Separator
    variables::Vector{Symbol}
    separator::Function
    contractor::Contractor
    expression::Expr
end

doc"""CombinationSeparator is a separator that is a combination (union, intersection,
or complement) of other separators.
"""
type CombinationSeparator <: Separator
    variables::Vector{Symbol}
    separator::Function
    expression::Expr
end

doc"""Create a ConstraintSeparator from a given constraint expression."""
function ConstraintSeparator(ex::Expr)
    expr, constraint = parse_comparison(ex)

    if constraint == entireinterval()
        throw(ArgumentError("Must give explicit constraint"))
    end

    if isa(expr, Symbol)
        expr = :(1 * $expr)  # convert symbol into expression
    end

    C = Contractor(expr)
    variables = C.variables[2:end]

    a = constraint.lo
    b = constraint.hi

    local outer

    f = X -> begin

        inner = C(a..b, X...)  # closure over the function C

        if a == -∞
            outer = C(b..∞, X...)

        elseif b == ∞
            outer = C(-∞..a, X...)

        else

            outer1 = C(-∞..a, X...)
            outer2 = C(b..∞, X...)

            outer = [ hull(x1, x2) for (x1,x2) in zip(outer1, outer2) ]
        end

        return (inner, (outer...))

    end

    expression = :($expr ∈ $constraint)

    return ConstraintSeparator(variables, f, C, expression)

end

macro separator(ex::Expr)  # alternative name for constraint -- remove?
    ex = Meta.quot(ex)
    :(ConstraintSeparator($ex))
end

doc"""Create a separator from a given constraint expression, written as
standard Julia code.

e.g. `C = @constraint x^2 + y^2 <= 1`

The variables (`x` and `y`, in this case) are automatically inferred.
External constants can be used as e.g. `$a`:

```
a = 3
C = @constraint x^2 + y^2 <= $a
```
"""
macro constraint(ex::Expr)
    ex = Meta.quot(ex)
    :(ConstraintSeparator($ex))
end

function show(io::IO, S::Separator)
    println(io, "Separator:")
    print(io, "- variables: ")
    print(io, join(map(string, S.variables), ", "))
    println(io)
    print(io, "- expression: ")
    println(io, S.expression)
end

show_code(S::ConstraintSeparator) = show_code(S.contractor)

@compat (S::ConstraintSeparator)(X) = S.separator(X)
@compat (S::CombinationSeparator)(X) = S.separator(X)

# show_code(S::Separator) = show_code(S.contractor)

doc"Unify the variables of two separators"
function unify_variables(vars1, vars2)

    variables = unique(sort(vcat(vars1, vars2)))
    numvars = length(variables)  # total number of variables

    indices1 = indexin(vars1, variables)
    indices2 = indexin(vars2, variables)

    where1 = zeros(Int, numvars)  # inverse of indices1
    where2 = zeros(Int, numvars)

    for (i, which) in enumerate(indices1)
        where1[which] = i
    end

    for (i, which) in enumerate(indices2)
        where2[which] = i
    end

    return variables, indices1, indices2, where1, where2
end

# TODO: when S1 and S2 have different variables -- amalgamate!

doc"""
    ∩(S1::Separator, S2::Separator)

Separator for the intersection of two sets given by the separators `S1` and `S2`.
Takes an iterator of intervals (`IntervalBox`, tuple, array, etc.), of length
equal to the total number of variables in `S1` and `S2`;
it returns inner and outer tuples of the same length
"""
function ∩(S1::Separator, S2::Separator)

    variables, indices1, indices2, where1, where2 = unify_variables(S1.variables, S2.variables)
    numvars = length(variables)

    f = X -> begin

        inner1, outer1 = S1([X[i] for i in indices1])
        inner2, outer2 = S2([X[i] for i in indices2])

        if any(isempty, inner1)
            inner1 = emptyinterval(IntervalBox(X))
        else
            inner1 = [i ∈ indices1 ? inner1[where1[i]] : X[i] for i in 1:numvars]
        end

        if any(isempty, outer1)
            outer1 = emptyinterval(IntervalBox(X))
        else
            outer1 = [i ∈ indices1 ? outer1[where1[i]] : X[i] for i in 1:numvars]
        end

        if any(isempty, inner2)
            inner2 = emptyinterval(IntervalBox(X))
        else
            inner2 = [i ∈ indices2 ? inner2[where2[i]] : X[i] for i in 1:numvars]
        end

        if any(isempty, outer2)
            outer2 = emptyinterval(IntervalBox(X))
        else
            outer2 = [i ∈ indices2 ? outer2[where2[i]] : X[i] for i in 1:numvars]
        end

        # Treat as if had X[i] in the other directions, except if empty

        inner = tuple( [x ∩ y for (x,y) in zip(inner1, inner2) ]... )
        outer = tuple( [x ∪ y for (x,y) in zip(outer1, outer2) ]... )

        return (inner, outer)

    end

    expression = :($(S1.expression) ∩ $(S2.expression))

    return CombinationSeparator(variables, f, expression)

end

function ∪(S1::Separator, S2::Separator)

    variables, indices1, indices2, where1, where2 = unify_variables(S1.variables, S2.variables)
    numvars = length(variables)

    f = X -> begin
        inner1, outer1 = S1(tuple([X[i] for i in indices1]...))
        inner2, outer2 = S2(tuple([X[i] for i in indices2]...))

        if any(isempty, inner1)
            inner1 = emptyinterval(IntervalBox(X))
        else
            inner1 = [i ∈ indices1 ? inner1[where1[i]] : X[i] for i in 1:numvars]
        end

        if any(isempty, outer1)
            outer1 = emptyinterval(IntervalBox(X))
        else
            outer1 = [i ∈ indices1 ? outer1[where1[i]] : X[i] for i in 1:numvars]
        end

        if any(isempty, inner2)
            inner2 = emptyinterval(IntervalBox(X))
        else
            inner2 = [i ∈ indices2 ? inner2[where2[i]] : X[i] for i in 1:numvars]
        end

        if any(isempty, outer2)
            outer2 = emptyinterval(IntervalBox(X))
        else
            outer2 = [i ∈ indices2 ? outer2[where2[i]] : X[i] for i in 1:numvars]
        end

        inner = tuple( [x ∪ y for (x,y) in zip(inner1, inner2) ]... )
        outer = tuple( [x ∩ y for (x,y) in zip(outer1, outer2) ]... )

        return (inner, outer)
    end

    expression = :($(S1.expression) ∪ $(S2.expression))

    return CombinationSeparator(variables, f, expression)

end

function !(S::Separator)
    f = X -> begin
        inner, outer = S(X)
        return (outer, inner)
    end

    expression = :(!($(S.expression)))

    return CombinationSeparator(S.variables, f, expression)
end


text:
using RCall, JLD

R"
library(gridExtra)
library(cowplot)
library(tidyverse)
library(latex2exp)
library(viridis)
"

DF = load("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_6_2.jld")
EV_NGG__approx = DF["NGG_approx"]

DF_NGGM_approx = DataFrame(Exp = first.(EV_NGG__approx),Var =last.(EV_NGG__approx), beta = first.(grid_vec), sigma = last.(grid_vec))
DF_NGGM_approx.std = sqrt.(DF_NGGM_approx.Var)
DF_NGGM_approx.beta_log= log.(DF_NGGM_approx.beta)

R"
df_r= $DF_NGGM_approx
df_r$beta_scaled= range01(df_r$beta_log)
df_r$color = rgb(df_r$beta_scaled, df_r$sigma,0.5,1)
df_r$color_sigma = rgb(0, df_r$sigma,0.5,1)
df_r$color_beta = rgb(df_r$beta_scaled, 0,0.5,1)
"
R"
p <- df_r %>% ggplot(aes(x=Exp, y = std, group=sigma))  + geom_line(alpha = 0.8,linetype = 'longdash') + geom_point(aes(colour=color), size=4)
p_nggm_approx <- p  + geom_line(data =df_r, aes(x=Exp, y = std, group=beta), alpha = 0.9,linetype = 'dotted') +
         xlim(0, 100)+ ylim(0,15)+labs(y='Std', x='Expectation')+scale_colour_identity()+theme_classic()+ theme(plot.title = element_text(hjust = 0.5,size = 15), axis.text = element_text(size=15),legend.position='none')
"

R"p_nggm_approx
pdf('/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/figures/figure_6_2.pdf', width= 4, height = 4)
plot(p_nggm_approx)
dev.off()"

max_stars_repo_name: konkam/approximatingBNPpriors

max_issues_repo_name: konkam/approximatingBNPpriors

max_forks_repo_name: konkam/approximatingBNPpriors


text:
export ConcatCube, concatenateCubes
export mergeAxes
import ..ESDLTools.getiperm
import ..Cubes: _read, gethandle, caxes, iscompressed, cubechunks, chunkoffset

mutable struct PermCube{T,N,C} <: AbstractCubeData{T,N}
  parent::C
  perm::NTuple{N,Int}
end
Base.size(x::PermCube)=ntuple(i->size(x.parent,x.perm[i]),ndims(x))
Base.size(x::PermCube,i)=size(x.parent,x.perm[i])
caxes(v::PermCube)=caxes(v.parent)[collect(v.perm)]
getCubeDes(v::PermCube)=getCubeDes(v.parent)
permtuple(t,perm)=ntuple(i->t[perm[i]],length(t))
iscompressed(c::PermCube)=iscompressed(c.parent)
cubechunks(c::PermCube)=cubechunks(c.parent)[collect(c.perm)]
chunkoffset(c::PermCube)=chunkoffset(c.parent)[collect(c.perm)]
function subsetcube(p::PermCube; kwargs...)
  subscube = subsetcube(p.parent;kwargs...)
  lostaxes = setdiff(caxes(p.parent), caxes(subscube))
  ilose    = map(ax->findAxis(ax,p.parent), lostaxes)
  losemap  = cumsum([!in(i,ilose) for i=1:ndims(p.parent)])
  losemap[findall([in(i,ilose) for i=1:ndims(p.parent)])] .= 0
  newperm = []
  foreach(i->losemap[i]==0 ? nothing : push!(newperm,losemap[i]),p.perm)
  isempty(newperm) ? subscube : permutedims(subscube,(newperm...,))
end
function _read(x::PermCube{T,N},thedata::AbstractArray,r::CartesianIndices{N}) where {T,N}
  perm=x.perm
  iperm=getiperm(perm)
  r2=CartesianIndices((permtuple(r.indices,iperm)))
  sr = size(r2)
  aout = zeros(T,sr...)
  _read(x.parent,aout,r2)
  permutedims!(thedata,aout,perm)

end
Base.permutedims(x::AbstractCubeData{T,N},perm) where {T,N}=PermCube{T,N,typeof(x)}(x,perm)
function gethandle(c::PermCube,block_size)
  handle = gethandle(c.parent)
  PermutedDimsArray(handle,c.perm)
end

import Base.Iterators.product

mutable struct TransformedCube{T,N,F} <: AbstractCubeData{T,N}
  parents
  op::F
  cubeaxes::Vector{CubeAxis}
  properties::Dict{String,Any}
end

function Base.map(op, incubes::AbstractCubeData...; T::Type=eltype(incubes[1]))
  axlist=copy(caxes(incubes[1]))
  N=ndims(incubes[1])
  for i=2:length(incubes)
    all(caxes(incubes[i]).==axlist) || error("All cubes must have the same axes, cube number $i does not match")
    ndims(incubes[i])==N || error("All cubes must have the same dimension")
  end
  props=merge(cubeproperties.(incubes)...)
  TransformedCube{T,N,typeof(op)}(incubes,op,axlist,props)
end

Base.size(x::TransformedCube)=size(x.parents[1])
Base.size(x::TransformedCube{T,N},i) where {T,N}=size(x.parents[1],i)
caxes(v::TransformedCube)=v.cubeaxes
getCubeDes(v::TransformedCube)="Transformed cube $(getCubeDes(v.parents[1]))"
iscompressed(v::TransformedCube)=any(iscompressed,v.parents)
cubechunks(v::TransformedCube)=cubechunks(v.parents[1])
chunkoffset(v::TransformedCube)=chunkoffset(v.parents[1])
function subsetcube(v::TransformedCube{T}; kwargs...) where T
  newcubes = map(i->subsetcube(i; kwargs...), v.parents)
  map(v.op, newcubes..., T=T)
end

using Base.Cartesian
function _read(x::TransformedCube{T,N},thedata::AbstractArray,r::CartesianIndices{N}) where {T,N}
  ainter=map(x.parents) do c
    aouti=Array{eltype(c)}(undef,size(thedata))
    _read(c,aouti,r)
    aouti
  end
  map!(x.op,thedata,ainter...)
  return thedata
end

ops2 = [:+, :-,:/, :*, :max, :min]
for op in ops2
  eval(:(Base.$(op)(x::AbstractCubeData, y::AbstractCubeData)=map($op, x,y)))
  eval(:(Base.$(op)(x::AbstractCubeData, y::Number)          =map(i->$(op)(i,y),x)))
  eval(:(Base.$(op)(x::Number, y::AbstractCubeData)          =map(i->$(op)(x,i),y)))
end

ops1 = [:sin, :cos, :log, :log10, :exp, :abs]
for op in ops1
  eval(:(Base.$(op)(x::AbstractCubeData)=map($op, x)))
end

"""
    ConcatCube

Concatenate a list of cubes of the same type and axes to a composite cube.
"""
mutable struct ConcatCube{T,N} <: AbstractCubeData{T,N}
  cubelist::Vector
  cataxis::CubeAxis
  cubeaxes::Vector{CubeAxis}
  properties::Dict{String,Any}
end

"""
    function concatenateCubes(cubelist, cataxis::CategoricalAxis)

Concatenates a vector of datacubes that have identical axes to a new single cube along the new
axis `cataxis`
"""
function concatenateCubes(cl,cataxis::CubeAxis)
  length(cataxis.values)==length(cl) || error("cataxis must have same length as cube list")
  axlist=copy(caxes(cl[1]))
  T=eltype(cl[1])
  N=ndims(cl[1])
  for i=2:length(cl)
    all(caxes(cl[i]).==axlist) || error("All cubes must have the same axes, cube number $i does not match")
    eltype(cl[i])==T || error("All cubes must have the same element type, cube number $i does not match")
    ndims(cl[i])==N || error("All cubes must have the same dimension")
  end
  props=mapreduce(cubeproperties,merge,cl,init=Dict())
  ConcatCube{T,N+1}(collect(cl),cataxis,axlist,props)
end
function concatenateCubes(;kwargs...)
  cubenames = String[]
  for (n,c) in kwargs
    push!(cubenames,string(n))
  end
  cubes = map(i->i[2],collect(kwargs))
  findAxis("Variable",cubes[1]) === nothing || error("Input cubes must not contain a variable kwarg concatenation")
  concatenateCubes(cubes, CategoricalAxis("Variable",cubenames))
end
Base.size(x::ConcatCube)=(size(x.cubelist[1])...,length(x.cataxis))
Base.size(x::ConcatCube{T,N},i) where {T,N}=i==N ? length(x.cataxis) : size(x.cubelist[1],i)
caxes(v::ConcatCube)=[v.cubeaxes;v.cataxis]
iscompressed(x::ConcatCube)=any(iscompressed,x.cubelist)
cubechunks(x::ConcatCube)=(cubechunks(x.cubelist[1])...,1)
chunkoffset(x::ConcatCube)=(chunkoffset(x.cubelist[1])...,0)
getCubeDes(v::ConcatCube)="Collection of $(getCubeDes(v.cubelist[1]))"
using Base.Cartesian
function _read(x::ConcatCube{T,N},thedata::AbstractArray,r::CartesianIndices{N}) where {T,N}
  rnew = CartesianIndices(r.indices[1:end-1])
  for (j,i)=enumerate(r.indices[end])
    a=selectdim(thedata,N,j)
    _read(x.cubelist[i],a,rnew)
  end
  return thedata
end
function _write(x::ConcatCube{T,N},thedata::AbstractArray,r::CartesianIndices{N}) where {T,N}
  rnew = CartesianIndices(r.indices[1:end-1])
  for (j,i)=enumerate(r.indices[end])
    a=selectdim(thedata,N,j)
    _write(x.cubelist[i],a,rnew)
  end
  return nothing
end

using RecursiveArrayTools
function gethandle(c::ConcatCube,block_size)
  VectorOfArray(map(gethandle,c.cubelist))
end

struct SplitDimsCube{T,N,C} <: AbstractCubeData{T,N}
    parent::C
    newaxes
    newchunks::Tuple{Int,Int}
    isplit::Int
end

export splitdim
"""
Splits an axis into two, reshaping the data cube into a higher-order cube.
"""
function splitdim(c, dimtosplit, newdims)
    isplit = findAxis(dimtosplit,c)
    isplit ===nothing && error("Could not find axis in cube")
    axesold = caxes(c)
    length(axesold[isplit]) == prod(length.(newdims)) || error("Size of new dimensions don't match the old one")
    all(i->isa(i,CubeAxis), newdims) || error("Newdims must be a list of cube axes")
    splitchunk = cubechunks(c)[isplit]
    newchunks = if splitchunk==length(newdims[1])
        (length(newdims[1]),1)
    elseif splitchunk>length(newdims[1])
        iszero(rem(splitchunk,length(newdims[1]))) || error("You can not do this split since chunks woudl overlap")
        (length(newdims[1]),splitchunk ÷ length(newdims[1]))
    else
        iszero(rem(length(newdims[1]),splitchunk)) || error("You can not do this split since chunks woudl overlap")
        (length(newdims[1]) ÷ splitchunk, 1)
    end
    SplitDimsCube{eltype(c), ndims(c)+1, typeof(c)}(c,(newdims...,),newchunks,isplit)
end

function Base.size(x::SplitDimsCube)
    sp = size(x.parent)
    (sp[1:isplit-1]...,length.(newaxes)...,sp[isplit+1:end]...)
end
Base.size(x::SplitDimsCube{T,N},i) where {T,N} = if i<x.isplit
    size(x.parent,i)
elseif i==x.isplit
    length(x.newaxes[1])
elseif i==x.isplit+1
    length(x.newaxes[2])
else
    size(x.parent,i+1)
end
function caxes(v::SplitDimsCube)
    axold = caxes(v.parent)
    [axold[1:v.isplit-1]...;v.newaxes...;axold[v.isplit+1:end]...]
end
getCubeDes(v::SplitDimsCube)=getCubeDes(v.parent)
iscompressed(v::SplitDimsCube)=iscompressed(v.parent)
function cubechunks(v::SplitDimsCube)
    cc = cubechunks(v.parent)
    (cc[1:v.isplit-1]...,v.newchunks...,cc[v.isplit+1:end]...)
end
function chunkoffset(v::SplitDimsCube)
    co = chunkoffset(v.parent)
    #This can not be determined, so we just assume zero here
    (co[1:v.isplit-1]...,0,0,co[v.isplit+1:end]...)
end
function subsetcube(v::SplitDimsCube; kwargs...)
    newax = collect(v.newaxes)
    foreach(kwargs) do (k,v)
        findAxis(k,newax) !== nothing && error("Subsetting split axes is not yet supported")
    end
    newcube = subsetcube(v.parent; kwargs...)
    splitdim(newcube,caxes(newcube)[isplit],v.newaxes)
end
using Base.Cartesian
function _read(x::SplitDimsCube{T,N},thedata::AbstractArray,r::CartesianIndices{N}) where {T,N}
  r1 = r.indices[x.isplit]
  r2 = r.indices[x.isplit+1]
  if length(r2)==1 && length(r1)==length(x.newaxes[1])
    error("This read is not yet implemented")
  end
  istart = (first(r2)-1)*size(x,x.isplit) + first(r1)
  iend   = (last(r2)-1) *size(x,x.isplit) + last(r1)
  indsold = r.indices
  inew = CartesianIndices((r.indices[1:x.isplit-1]...,istart:iend,r.indices[x.isplit+2:end]...))
  a2 = reshape(thedata,size(inew))
  _read(x.parent,a2,inew)
  return thedata
end


text:
todisplay = VERSION <= VersionNumber("1.2") ? true : false

using CAGD
using SparseArrays
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

if todisplay
	using ViewerGL
	GL = ViewerGL
	include("../examples/3d/views.jl")
end

V = [
    0.0 +1  0  0  0
    0.0  0  0 +1 -1
    0.0  0 +1  0  0
];

W = [
     0.0  0.0  0.0
    +1.0  0.0  0.0
     0.0  0.0 +1.0
     0.0 +1.0  0.0
     0.0 -1.0  0.0
];

EV = SparseArrays.sparse(Int8.([
    -1 +1  0  0  0
    -1  0 +1  0  0
    -1  0  0 +1  0
    -1  0  0  0 +1
     0 -1 +1  0  0
     0 -1  0 +1  0
     0 -1  0  0 +1
     0  0 -1 +1  0
     0  0 -1  0 +1
]));

FE = SparseArrays.sparse(Int8.([
    +1 -1  0  0 +1   0  0  0  0
    +1  0 -1  0  0  +1  0  0  0
    +1  0  0 -1  0   0 +1  0  0
     0 +1 -1  0  0   0  0 +1  0
     0 +1  0 -1  0   0  0  0 +1
     0  0  0  0 +1  -1  0 +1  0
     0  0  0  0 +1   0 -1  0 +1
]));

CF = SparseArrays.sparse(Int8.([
    -1 +1  0 -1  0 +1  0
    +1  0 -1  0 +1  0 -1
     0 -1 +1 +1 -1 -1 +1
]));

model = CAGD.Model(V)
CAGD.addModelCells!(model, 1, EV)
CAGD.addModelCells!(model, 2, FE)
CAGD.addModelCells!(model, 3, CF)

atol = 1e-6;
if todisplay  displayModel(model)  end

split_model = CAGD.facesplitting(model, atol = atol)

if todisplay  displayModel(split_model)  end

congr_model = CAGD.mergeModelVertices(split_model, err=atol, signed_merge=true)
#split_model = CAGD.facesplitting(congr_model, atol = atol)
#congr_model = CAGD.mergeModelVertices(split_model, err=atol, signed_merge=true)

if todisplay  displayModel(congr_model)  end

gift_model = deepcopy(congr_model)
FCcycles, bicon_comps = CAGD.tgw(congr_model, 3, atol = atol)
CAGD.addModelCells!(gift_model, 3, convert(Lar.ChainOp, FC'))

if todisplay  viewExplode(gift_model, complete = false)  end

# Since first face is inner, then retrieving the external is needed

externals = [bc[CAGD.getExternalBoundary(congr_model, 3, FCcycles, bc)] for bc in bicon_comps]
internals = [setdiff(bicon_comps[i], externals[i]) for i = 1 : length(externals)]
assembled_model = deepcopy(congr_model)

FC = FCcycles[:, internals...]
JK = ones(Int8, length(externals))
FC = [FCcycles * SparseArrays.sparse(externals, JK, JK, size(FCcycles, 2), 1) FC]
CAGD.addModelCells!(assembled_model, 3, convert(Lar.ChainOp, FC'))

if todisplay  viewExplode(assembled_model)  end


text:
module NormalMeanVarianceNodeTest

using Test
using ReactiveMP
using Random

import ReactiveMP: make_node

@testset "NormalMeanVarianceNode" begin

    @testset "Creation" begin
        node = make_node(NormalMeanVariance)

        @test functionalform(node)          === NormalMeanVariance
        @test sdtype(node)                  === Stochastic()
        @test name.(interfaces(node))       === (:out, :μ, :v)
        @test factorisation(node)           === ((1, 2, 3), )
        @test localmarginalnames(node)      === (:out_μ_v, )
        @test metadata(node)                === nothing

        node = make_node(NormalMeanVariance, meta = 1)

        @test functionalform(node)          === NormalMeanVariance
        @test sdtype(node)                  === Stochastic()
        @test name.(interfaces(node))       === (:out, :μ, :v)
        @test factorisation(node)           === ((1, 2, 3), )
        @test localmarginalnames(node)      === (:out_μ_v, )
        @test metadata(node)                === 1

        node = make_node(NormalMeanVariance, factorisation = ((1, ), (2, 3)))

        @test functionalform(node)          === NormalMeanVariance
        @test sdtype(node)                  === Stochastic()
        @test name.(interfaces(node))       === (:out, :μ, :v)
        @test factorisation(node)           === ((1, ), (2, 3, ))
        @test localmarginalnames(node)      === (:out, :μ_v, )
        @test metadata(node)                === nothing
    end

    @testset "AverageEnergy" begin

        begin
            q_out = PointMass(0.956629)
            q_μ   = NormalMeanVariance(0.255332, 0.762870)
            q_τ   = GammaShapeRate(0.93037, 0.79312)

            for N in (NormalMeanPrecision, NormalMeanVariance, NormalWeightedMeanPrecision), G in (GammaShapeRate, GammaShapeScale)
                marginals = (Marginal(q_out, false, false), Marginal(convert(N, q_μ), false, false), Marginal(convert(G, q_τ), false, false))
                @test score(AverageEnergy(), NormalMeanPrecision, Val{(:out, :μ, :τ)}, marginals, nothing) ≈ 1.8879401707928936
            end
        end

        begin
            q_out = NormalMeanVariance(0.148725, 0.483501)
            q_μ   = NormalMeanVariance(0.992776, 0.545851)
            q_v   = GammaShapeRate(0.309396, 0.343814)

            for N in (NormalMeanPrecision, NormalMeanVariance, NormalWeightedMeanPrecision), G in (GammaShapeRate, GammaShapeScale)
                marginals = (Marginal(q_out, false, false), Marginal(convert(N, q_μ), false, false), Marginal(convert(G, q_v), false, false))
                @test score(AverageEnergy(), NormalMeanPrecision, Val{(:out, :μ, :τ)}, marginals, nothing) ≈ 2.86416186172411
            end
        end

        begin
            q_out_μ = MvNormalMeanCovariance([0.2818402997601115, 0.0847764277628964], [0.7059042678955475 0.3595204552322394; 0.3595204552322394 0.22068491824258746])
            q_v     = GammaShapeRate(0.49074414, 0.4071772)

            for N in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision), G in (GammaShapeRate, GammaShapeScale)
                marginals = (Marginal(convert(N, q_out_μ), false, false), Marginal(convert(G, q_v), false, false))
                @test score(AverageEnergy(), NormalMeanPrecision, Val{(:out_μ, :τ)}, marginals, nothing) ≈ 1.6231194045861121
            end
        end

        begin
            q_out_μ = MvNormalMeanCovariance([0.8378350736808462, 0.41494396892699026], [1.0074451742986652 0.4369298270351709; 0.4369298270351709 0.19572138039218784])
            q_v     = GammaShapeRate(0.435485, 0.5269575)

            for N in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision), G in (GammaShapeRate, GammaShapeScale)
                marginals = (Marginal(convert(N, q_out_μ), false, false), Marginal(convert(G, q_v), false, false))
                @test score(AverageEnergy(), NormalMeanPrecision, Val{(:out_μ, :τ)}, marginals, nothing) ≈ 1.969539193740776
            end
        end

    end

end
end