text:
using FunctionalTables, Test
using FunctionalTables:
    # utilities
    cancontain, narrow, append1, split_namedtuple, is_prefix, merge_default_types,
    # column ordering building blocks
    ColumnOrdering, mask_ordering, table_ordering, cmp_ordering, retained_ordering,
    ordering_repr, split_compatible_ordering,
    # column collection building blocks
    SINKCONFIG, collect_column, collect_columns, RLEVector, TrustLength
import Tables
using Statistics: mean

include("utilities.jl")         # utilities for tests

####
#### Utilities and low-level building blocks
####

@testset "narrow" begin
    @test narrow(true) ≡ true
    @test narrow(false) ≡ false
    @test narrow(1) ≡ Int8(1)
    @test narrow(0) ≡ Int8(0)
    @test narrow(-9) ≡ Int8(-9)
    @test narrow(128) ≡ Int16(128)
    @test narrow(2^17) ≡ Int32(2^17)
    @test narrow(2^32) ≡ Int(2^32)
end

@testset "cancontain" begin
    @test cancontain(Int8, -7)
    @test !cancontain(Int8, typemax(Int8) + 1)
    @test !cancontain(String, -7)
    @test cancontain(Float64, Int8(9))
    @test cancontain(Float64, maxintfloat(Float64))
    @test !cancontain(Float64, Int(maxintfloat(Float64)) + 1)
    for T in (Float32, Float64)
        # Test around `maxintfloat`, as it will fail if the number can be contained by
        # accident.
        for elt in vcat(Int(maxintfloat(T)) .+ (-1:1), -Int(maxintfloat(T)) .+ (-1:1))
            v = T[]
            push!(v, elt)
            (oftype(elt, v[1]) ≡ elt) ≡ cancontain(T, elt)
        end
    end
end

@testset "append1" begin
    a1 = append1(ones(Int8, 2), missing)
    @test eltype(a1) ≡ Union{Missing, Int8}
    @test [1, 1, missing] ≅ a1
end

@testset "splitting named tuples" begin
    s = NamedTuple{(:a, :c)}
    @test @inferred(split_namedtuple(s, (a = 1, b = 2, c = 3, d = 4))) ≡
        ((a = 1, c = 3), (b = 2, d = 4))
    @test @inferred(split_namedtuple(s ,(c = 1, b = 2, a = 3, d = 4))) ≡
        ((a = 3, c = 1), (b = 2, d = 4))
    @test_throws ErrorException split_namedtuple(s, (a = 1, b = 2))
end

@testset "merge_default_types" begin
    T = NamedTuple{(:a, :b), Tuple{Union{}, Union{}}}
    @test merge_default_types(T, typeof(NamedTuple())) ≡ T
    @test merge_default_types(T, NamedTuple{(:c, :d), Tuple{Missing, String}}) ≡ T
    @test merge_default_types(T, NamedTuple{(:a, ), Tuple{Missing}}) ≡
        NamedTuple{(:a, :b), Tuple{Missing, Union{}}}
end

@testset "picking and wrapping" begin
    x = 1:5
    ax = wrapping(:a).(x)
    @test ax == collect((a = i,) for i in x)
    @test picking(:a).(ax) == x
end

####
#### columns
####

@testset "collect by names" begin
    itr = [(a = i, b = Float64(i), c = 'a' + i - 1) for i in 1:10]
    ordering = table_ordering((:a, :b, :c))
    len, result, rule = collect_columns(SinkConfig(;useRLE = false), itr, TrustOrdering(ordering))
    @test len ≡ TrustLength(10)
    @test result isa NamedTuple{(:a, :b, :c), Tuple{Vector{Int8}, Vector{Float64}, Vector{Char}}}
    @test result.a ≅ 1:10
    @test result.b ≅ Float64.(1:10)
    @test result.c ≅ (0:9) .+ 'a'
    @test rule ≡ TrustOrdering(ordering)
end

@testset "RLE creation" begin
    r = RLEVector{Missing}(Int8, Missing)
    @test eltype(typeof(r)) ≡ Missing
    @test length(r) == 0

    r = RLEVector{Missing}(Int8, Float64)
    @test eltype(typeof(r)) ≡ Float64
    @test length(r) == 0

    r = RLEVector{Missing}(Int8, Union{})
    @test eltype(typeof(r)) ≡ Union{}
    @test length(r) == 0

    r = RLEVector{Missing}(Int8, Union{Missing, Float64})
    @test eltype(typeof(r)) ≡ Union{Missing, Float64}
    @test length(r) == 0
end

@testset "simple RLE" begin
    v = vcat(fill(1, 10), fill(missing, 5), fill(2, 20))
    s = collect_column(SINKCONFIG, v)
    @test length(s) == 35
    @test eltype(typeof(s)) ≡ Union{Int8, Missing}
    @test s.data == [1, 2]
    @test s.counts == [10, -5, 20]
    @test collect(s) ≅ v
end

@testset "RLE corner cases" begin
    s = collect_column(SINKCONFIG, [missing])
    @test length(s) == 1
    @test collect(s)[1] ≡ missing
    @test eltype(typeof(s)) ≡ Missing
    s = collect_column(SINKCONFIG, [130])
    @test length(s) == 1
    @test eltype(typeof(s)) ≡ Int16
    @test collect(s) == [130]
end

@testset "overrun RLE" begin
    v = vcat(fill(1, 300), fill(missing, 5), fill(2, 20), fill(missing, 200))
    s = collect_column(SINKCONFIG, v)
    @test length(s) == length(v)
    @test eltype(typeof(s)) ≡ Union{Int8, Missing}
    @test s.data == [1, 1, 1, 2]
    @test s.counts == [127, 127, 300-(2*127), -5, 20, -128, -200+128]
    @test collect(s) ≅ v
end

@testset "large collection" begin
    v = randvector(1000)
    len, columns, ordering = collect_columns(SINKCONFIG, [(a = a, ) for a in v], TrustOrdering())
    @test len ≡ TrustLength(length(v))
    @test collect(columns.a) ≅ v
    @test ordering ≡ TrustOrdering()
end

@testset "is_prefix" begin
    @test is_prefix((:a, ), (:a, :b, :c))
    @test is_prefix((:a, :b, :c), (:a, :b, :c))
    @test is_prefix((), ())
    @test !is_prefix((:b, :c), (:a, :b, :c))
    @test !is_prefix((:b, ), ())
    @test !is_prefix((:b, :c, :a), (:a, :b, :c))
end

####
#### Ordering building blocks
####

@testset "column ordering specifications" begin
    @test table_ordering((:a, :b => reverse, ColumnOrdering(:c, false))) ==
        (ColumnOrdering(:a, false), ColumnOrdering(:b, true), ColumnOrdering(:c, false))
    @test_throws ArgumentError table_ordering(("foobar", "baz")) # invalid
    @test_throws ArgumentError TrustOrdering(:a, :a)             # duplicate

    @test ordering_repr(()) == "no ordering"
    @test ordering_repr(table_ordering((:a, :b => reverse))) == "ordering ↑a ↓b"
end

@testset "mask_ordering" begin
    t = table_ordering((:a, :b => reverse, :c))
    @test mask_ordering(t, (:a, )) ≡ t[[1]]
    @test mask_ordering(t, (:b, )) ≡ ()
    @test mask_ordering(t, (:c, )) ≡ ()
    @test mask_ordering(t, (:b, :a, )) ≡ t[[1,2]]
    @test mask_ordering(t, (:c, :a, )) ≡ t[[1]]

    o = table_ordering((:a, :b, :c))
    @test mask_ordering(o, (:d, :e), true) ≡ o
    @test mask_ordering(o, (:c, :b), true) ≡ table_ordering((:a, ))
    @test mask_ordering(o, (:a, :b, :c), true) ≡ table_ordering(())
end

@testset "retained ordering" begin
    o = table_ordering((:a, :b => reverse))
    row = (a = 1, b = 2, c = 3)
    @test retained_ordering(o, row, row) ≡ o
    @test retained_ordering(o, row, (a = 2, b = 1, c = -1)) ≡ o
    @test retained_ordering(o, row, (a = 2, b = 3, c = -1)) ≡ table_ordering((:a, ))
    @test @inferred(retained_ordering(table_ordering(()), row, row)) ≡ table_ordering(())
end

@testset "split compatible ordering" begin
    ordering = table_ordering((:a, :b => reverse, :c))
    @test split_compatible_ordering(ordering, (:a, )) ≡ ordering[1:1]
    @test split_compatible_ordering(ordering, (:a, :b)) ≡ ordering[1:2]
    @test split_compatible_ordering(ordering, (:b, :a)) ≡ ordering[[2, 1]]
    @test split_compatible_ordering(ordering, (:b, :a, :d)) ≡
        table_ordering((ordering[[2, 1]]..., :d))
    @test split_compatible_ordering((), (:b, :a, :d)) ≡ table_ordering((:b, :a, :d))
    @test split_compatible_ordering((), ()) ≡ ()
end

####
#### FunctionalTable API
####

@testset "FunctionalTable basics and column operations" begin
    A = 1:10
    B = 'a':('a'+9)
    C = Float64.(21:30)
    ft = FunctionalTable((a = A, b = B, c = C))
    @test Base.IteratorEltype(ft) ≡ Base.HasEltype()
    @test eltype(typeof(ft)) ≡ typeof((a = first(A), b = first(B), c = first(C)))
    @test Base.IteratorSize(ft) ≡ Base.HasLength()
    @test length(ft) ≡ length(A)
    @test @inferred(keys(ft)) ≡ (:a, :b, :c)
    @test @inferred(propertynames(ft)) ≡ (:a, :b, :c)
    @test @inferred(columns(ft)) ≡ getfield(ft, :columns)
    column_a(ft) = ft.a         # wrap for inference to kick in
    @test @inferred(column_a(ft)) ≡ A
    @test_throws ErrorException ft.nonexistent
    @test @inferred(values(ft)) ≡ values(columns(ft))
    @test ft[(:a, :b)] ≅ FunctionalTable((a = A, b = B))
    @test ft[drop = (:a, :b)] ≅ FunctionalTable((c = C,))
    @test ft[:a] == A
    @test FunctionalTable(ft) ≅ ft
    cols = @inferred map(collect, columns(ft))
    @test all(isa.(values(cols), AbstractVector))
    @test cols.a == A && cols.a ≢ A
    @test cols.b == B && cols.b ≢ B
    @test cols.c == C && cols.c ≢ C
    @test @inferred(FunctionalTable(ft)) ≡ ft # same object
end

@testset "merging" begin
    A = 1:10
    B = 'a':('a'+9)
    C = Float64.(21:30)
    A2 = .-A
    ft = FunctionalTable((a = A, b = B), VerifyOrdering(:a, :b))
    @test @inferred(merge(ft, FunctionalTable((c = C, )))) ≅
        FunctionalTable((a = A, b = B, c = C), VerifyOrdering(:a, :b))
    @test_throws ArgumentError merge(ft, FunctionalTable((c = C, a = A2)))
    @test @inferred(merge(ft, FunctionalTable((c = C, a = A2)); replace = true)) ≅
        FunctionalTable((a = A2, b = B, c = C))
    @test @inferred(merge(ft, FunctionalTable((c = C, a = A2)); replace = true)) ≅
        merge(ft, (c = C, a = A2); replace = true)
end

@testset "map (direct)" begin
    A = 1:10
    B = 'a':('a'+9)
    ft = @inferred FunctionalTable((a = A, b = B), VerifyOrdering(:a, :b))
    f(row) = (b = row.a + 1, c = row.b + 2)
    B2 = A .+ 1
    C = collect(B .+ 2)
    ft2 = map(f, ft)
    # NOTE map removes ordering
    @test ft2 ≅ FunctionalTable((b = B2, c = C))
    ft3 = merge(f, ft; replace = true)
    # NOTE as :b is replaced, its ordering is removed
    @test ft3 ≅ FunctionalTable((a = A, b = B2, c = C), VerifyOrdering((:a, )))
    # overlap, without replacement
    @test_throws ArgumentError merge(f, ft)
end

@testset "filter" begin
    A = 1:5
    B = 'a':'e'
    o = VerifyOrdering(:a, :b)
    ft = FunctionalTable((a = A, b = B), o)
    @test filter(row -> isodd(row.a), ft) ≅
        FunctionalTable((a = [1, 3, 5], b = ['a', 'c', 'e']), o)
end

@testset "RepeatRow" begin
    rr = RepeatRow((a = 1, b = 2))
    ft = FunctionalTable((a = fill(1, 3), b = fill(2, 3)))
    @test @inferred(FunctionalTable(3, rr)) ≅ ft
    ft2 = FunctionalTable((c = 4:6, d = 7:9))
    @test @inferred(merge(rr, ft2)) ≅ @inferred(merge(ft, ft2))
    @test @inferred(merge(ft2, rr)) ≅ @inferred(merge(ft2, ft))
end

@testset "split by 1" begin
    keycounts = [:a => 10, :b => 17, :c => 19]
    ft = FunctionalTable(mapreduce(((k, c), ) -> [(sym = k, val = i)
                                                  for i in 1:c], vcat, keycounts),
                         TrustOrdering(:sym, :val))
    g = by(ft, (:sym, ))        # FIXME eventually this should be @inferred
    @test eltype(typeof(g)) ≡ Tuple{NamedTuple{(:sym, ), Tuple{Symbol}},
                                    FunctionalTable{NamedTuple{(:val, ), Tuple{Vector{Int8}}},
                                                    Tuple{}}} # ordering, NOTE will change to :val
    cg = @inferred collect(g)
    for (i, (s, c)) in enumerate(keycounts)
        @test cg[i] ≅ ((sym = s, ), FunctionalTable((val = 1:c, )))
    end

    # empty split keys
    @test collect(by(ft, ())) ≅ [(NamedTuple(), FunctionalTable(ft, VerifyOrdering()))]
end

@testset "split by 2" begin
    A = [1, 1, 1, 2, 2]
    B = 'a':'e'
    ft = FunctionalTable((a = A, b = B), VerifyOrdering(:a))
    g = by(ft, (:a, ))
    @test Base.IteratorSize(g) ≡ Base.SizeUnknown()
    result = collect(g)
    @test result ≅ [((a = 1, ), FunctionalTable((b = ['a', 'b', 'c'],))),
                    ((a = 2, ), FunctionalTable((b = ['d', 'e'],)))]
end

@testset "tables interface" begin
    cols = (a = 1:10, b = 'a' .+ (1:10))
    ft = FunctionalTable(cols)
    # general sanity checks
    @test Tables.columntable(ft) == cols
    @test Tables.rowtable(ft) == Tables.rowtable(cols)
    # detailed API testing
    @test Tables.istable(ft)
    @test Tables.rowaccess(ft)
    @test Tables.schema(Tables.rows(ft)) == Tables.schema(Tables.rows(cols))
    @test Tables.columnaccess(ft)
    @test Tables.schema(Tables.columns(ft)) == Tables.schema(Tables.columns(cols))
end

@testset "ordering comparisons" begin
    ordering = table_ordering((:a, :b => reverse))
    @test @inferred cmp_ordering(ordering, (a = 1, b = 2), (a = 1, b = 2)) == 0
    @test @inferred cmp_ordering(ordering, (a = 1, b = 3), (a = 1, b = 2)) == -1
    @test @inferred cmp_ordering(ordering, (a = 0, b = 3), (a = 1, b = 2)) == -1
    @test @inferred cmp_ordering(ordering, (a = 1, b = 2), (b = 2, a = 1)) == 0 # order irrelevant
    @test_throws ErrorException cmp_ordering(ordering, (c = 1, ), (c = 1, ))    # no such field
end

@testset "sort" begin
    ft = FunctionalTable((a = [1, -1, 3, 1, 2],
                          b = [2, 2, 1, 2, 2],
                          c = 1:5))
    sft = sort(ft, (:b, :a => reverse))
    @test sft ≅ FunctionalTable((a = [3, 2, 1, 1, -1],
                                 b = [1, 2, 2, 2, 2],
                                 c = [3, 5, 1, 4, 2]),
                                VerifyOrdering(:b, :a => reverse))
    @test columns(sort(sft, (:b, ))) ≡ columns(sft) # prefix ordering

    let ordering = (:c => reverse, ) # test convenience syntax
        @test sort(ft, ordering) ≅ sort(ft, ordering...)
    end
end

@testset "map by" begin
    a = [1, 1, 1, 2, 2]
    b = 1:5
    ft = FunctionalTable((a = a, b = b), VerifyOrdering(:a, :b))
    f(_, ft) = map(sum, columns(ft))
    @test map(Ref ∘ f, by(ft, (:a, ))) ≅
        FunctionalTable((a = [1, 2], b = [6, 9]), TrustOrdering(:a, ))
end

@testset "map by with filtering" begin
    ft = FunctionalTable((a = 1:5, b = 1:5), VerifyOrdering(:a))
    @test map((_, ft) -> isodd(first(ft.b)) ? ft : (), by(ft, :a)) ≅
        FunctionalTable((a = [1, 3, 5], b = [1, 3, 5]), TrustOrdering(:a))
end

@testset "map by empty subtables" begin
    ft = FunctionalTable((a = [1, 1, 1, 2, 2], ), VerifyOrdering(:a))
    @test map((_, ft) -> Ref((len = length(ft), )), by(ft, (:a, ))) ≅
        FunctionalTable((a = [1, 2], len = [3, 2]), TrustOrdering(:a))
end

@testset "corner cases for collecting and ordering" begin
    A = (a = 1, )
    AA = [A, A]

    # different keys by row
    @test_throws ArgumentError FunctionalTable([A, (a = 1, b = 2)])

    # field specified by ordering is missing
    @test_throws ErrorException FunctionalTable(AA, VerifyOrdering(:b, ))
    @test_throws ArgumentError FunctionalTable(AA, TrustOrdering(:b, ))

    # prefix narrows ordering silently
    @test FunctionalTable(AA, TryOrdering(:b)) ≅ FunctionalTable(AA)

    # FIXME not implemented yet
    @test_skip FunctionalTable((a = [2, 1], ), (:a, ), TryOrdering()) ≅
        FunctionalTable((a = [1, 2], ), (), TrustOrdering())

    # nothing in the rows, also special constructor
    emptyrows = (NamedTuple() for _ in 1:1000)
    @test FunctionalTable(emptyrows) ≅ FunctionalTable(1000)

    # no rows
    @test FunctionalTable(()) ≅ FunctionalTable(0, NamedTuple(), TrustOrdering())
end

@testset "printing" begin
    ft = FunctionalTable((a = [1, 2], b = [3, 4]), TrustOrdering(:a, :b => reverse))
    reprft = """
    FunctionalTable of 2 rows, ordering ↑a ↓b
        a = Int64[1, 2]
        b = Int64[3, 4]"""
    @test repr(ft) == reprft
end

@testset "first" begin
    a = 1:1000
    b = fill(9, length(a))
    rt = Tables.rowtable((a = a, b = b))
    s = VerifyOrdering(:a, :b)
    ft = FunctionalTable(rt, s)
    @test first(ft, 100) ≅ FunctionalTable(Iterators.take(rt, 100), s)
end

@testset "rename" begin
    ft = FunctionalTable((a = 1:3, b = 4:6), VerifyOrdering(:b, :a => reverse))
    @test rename(ft, (a = :aa, b = :bb)) ≅
        FunctionalTable((aa = 1:3, bb = 4:6), VerifyOrdering(:bb, :aa => reverse))
    @test rename(ft, (a = :aa, )) ≅
        FunctionalTable((aa = 1:3, b = 4:6), VerifyOrdering(:b, :aa => reverse))

    # extras - error
    @test_throws ArgumentError rename(ft, (a = :aa, b = :bb, c = :cc))

    # duplicate keys
    @test_throws ArgumentError rename(ft, (a = :aa, b = :aa))

    # exchange names
    @test rename(ft, (a = :b, b = :a)) ≅
        FunctionalTable((b = 1:3, a = 4:6), VerifyOrdering(:a, :b => reverse))
end

@testset "eltype split tables" begin
    ft = FunctionalTable((a = [1, 1, 1, 2, 2, 2],
                          b = [1, 1, 1, 4, 4, 4]))
    ft2 = map(by(ft, :a)) do _, ft
        ft2 = filter(r -> r.b ≤ 1, ft) # here, ft2 is empty, but inference b has Int64 eltype
        Ref((bsum = sum(ft2.b), ))
    end
    @test ft2 ≅ FunctionalTable((a = Int8[1, 2], bsum = Int8[3, 0]), TrustOrdering(:a))
end

@testset "aggregator and ignoreindex" begin
    ft = FunctionalTable((a = [1, 1, 1, 2, 2, 2],
                          b = [1, 1, 1, 4, 4, 4]))
    ft_mean = FunctionalTable((a = Int8[1, 2], b = Int8[1, 4]), TrustOrdering(:a))
    @test map(aggregator(mean), by(ft, :a)) ≅ ft_mean # same as below
    @test map(ignoreindex(t -> Ref(map(mean, columns(t)))), by(ft, :a)) ≅ ft_mean
end


text:
module TCX
using EzXML, Dates, DataFrames, Geodesy
import Base.show

export parse_tcx_dir, parse_tcx_file, getActivityType, getDataFrame, getDistance, getDistance2, getDuration, getAverageSpeed, getAveragePace

struct TrackPoint
    Time::DateTime
    Latitude::Float64
    Longtitude::Float64
    HeartRateBpm::Int32
    AltitueMeter::Float64
    DistanceMeter::Float64
end

struct TCXRecord
    Id::DateTime
    Name::String
    ActivityType::String
    DistanceStatic::Float64
    DurationStatic::Float64
    HeartRate::Int32
    TrackPoints::Array{TrackPoint}
end

DictDateFormats = Dict(
    :24 => "yyyy-mm-ddTHH:MM:SS.sssZ",
    :20 => "yyyy-mm-ddTHH:MM:SSZ",
    )

function parse_tcx_file(file::String)
    file_path = abspath(file)
    if isfile(file_path) == false
        return 404, nothing
    end

    xmldoc = try readxml(file_path)
    catch e
       if isa(e, EzXML.XMLError)
           # Not a valid XML document
           @warn "Invalid XML document: $file_path"
           return 400, nothing
       else
           return 500, nothing
       end
    end

    root_element = root(xmldoc)
    # Check if TCX
    if nodename(root_element) != "TrainingCenterDatabase"
        @warn "Invalid TCX document: $file_path"
        return 400, nothing
    end

    # Type - "/*/*[1]/*[1]/@Sport"
    aType = nodecontent(findfirst("/*/*[1]/*[1]/@Sport", xmldoc))
    # Id - "/*/*[1]/*/*[1]"
    xid =nodecontent(findfirst("/*/*[1]/*/*[1]", xmldoc))
    aId = DateTime(xid,DictDateFormats[length(xid)])
    # Name = "/*/*[1]/*[1]/*[3]"
    aName = nodecontent(findfirst("/*/*[1]/*/*[2]", xmldoc))
    # Lap - "/*/*[1]/*/*[2]"
    # TotalSeconds - "/*/*[1]/*/*[2]/*[1]"
    aTime = parse(Float64, nodecontent(findfirst("/*/*[1]/*/*[2]/*[1]", xmldoc)))
    aDistance = parse(Float64, nodecontent(findfirst("/*/*[1]/*/*[2]/*[2]", xmldoc)))
    # DistanceMeters - "/*/*[1]/*/*[2]/*[2]"
    # AverageHeartRateBpm - "/*/*[1]/*/*[2]/*[5]/*[1]"
    xbpm = findfirst("/*/*[1]/*/*[2]/*[5]/*[1]", xmldoc)
    if xbpm === nothing
        aHeartRateBpm = 0
    else
        aHeartRateBpm = parse(Int32, nodecontent(xbpm))
    end
    # TrackPoints - "/*/*[1]/*/*[2]/*[9]/*"
    tp_Points = findall("/*/*[1]/*/*[2]/*[9]/*", xmldoc)
    aTrackPoints = Array{TrackPoint, size(tp_Points, 1)}[]
    for tp in tp_Points
        xtime = nodecontent(findfirst("./*[1]", tp))
        tp_time = DateTime(xtime, DictDateFormats[length(xtime)])
        xlat = findfirst("./*[2]/*[1]", tp)
        if xlat !== nothing
            tp_lat = parse(Float64, nodecontent(xlat))
        else
            continue
        end
        tp_lont = parse(Float64, nodecontent(findfirst("./*[2]/*[2]", tp)))
        xbpm = findfirst("./*[5]/*[1]", tp)
        if xbpm !== nothing
            tp_bpm = parse(Int32, nodecontent(findfirst("./*[5]/*[1]", tp)))
        else
            tp_bm = 0
        end
        tp_dist = parse(Float64, nodecontent(findfirst("./*[3]", tp)))
        tp_alt = parse(Float64, nodecontent(findfirst("./*[4]", tp)))

        aTrackPoints = vcat(aTrackPoints, TrackPoint(tp_time, tp_lat, tp_lont, tp_bpm, tp_dist, tp_alt))
    end

    return 200, TCXRecord(aId, aName, aType, aDistance, aTime, aHeartRateBpm, aTrackPoints)
end

function parse_tcx_dir(path::String)
    if ispath(path) == false
        @warn "Invalid path: $path"
        return 500, nothing
    end

    tcxArray = Array{TCXRecord}[]
    searchdir(path, key) = filter(x->occursin(key, x), readdir(path))

    for f in searchdir(path, ".tcx")
        err, tcx = parse_tcx_file(joinpath(path, f))
        if err == 200
            tcxArray = vcat(tcxArray, tcx)
        end
    end

    if length(tcxArray) > 0
        return 200, tcxArray
    else
        return 404, nothing
    end
end

function getActivityType(record::TCXRecord)
    return record.ActivityType
end

function getDataFrame(record::TCXRecord)
    return DataFrame(record.TrackPoints)
end

function getDataFrame(tcxArray::Array{Any, 1})
    aTP = Array{TrackPoint}[]
    for t in tcxArray
        aTP = vcat(aTP, t.TrackPoints)
    end
    return DataFrame(aTP)
end

function getDistance(record::TCXRecord)
    return record.DistanceStatic
end

function getDistance2(record::TCXRecord)
    total_distance = 0
    df = getDataFrame(record)
    num_of_rows = size(df, 1)
    for i in 1:num_of_rows
        if i < num_of_rows
            total_distance += distance(
                                       LLA(df[i, :Latitude], df[i, :Longtitude], df[i, :AltitueMeter]),
                                       LLA(df[i+1, :Latitude], df[i+1, :Longtitude], df[i+1, :AltitueMeter])
               )
        end
    end
    return total_distance
end

function getAverageSpeed(record::TCXRecord)
    return (record.DistanceStatic /1000) / (record.DurationStatic / 3600)  # km/h
end

function getAveragePace(record::TCXRecord)
    return (record.DurationStatic / 60) / (record.DistanceStatic / 1000) # min/km
end

function getDuration(record::TCXRecord)
    return record.DurationStatic
end

Base.show(io::IO, tcx::TCXRecord) = print(io, "$(tcx.ActivityType) $(tcx.DistanceStatic/1000) km at $(tcx.Id) for $(tcx.DurationStatic) seconds.")
end #module_end

ext: jl


text:
function pretrain(epochs, L, ps, data, opt; print_freq = 1000, lr_init = 1e-2,
	cyclic_lr = false)
	training_loss = 0.0
	new_lr = 0.0
	for ep in 1:epochs
		if cyclic_lr
			new_lr = cyclic_LR(ep, epochs, lr_init=lr_init, lr_ratio=0.05)
			opt.eta = new_lr
		end
		for d in data
			gs = gradient(ps) do
				training_loss = L(d...)
				return training_loss
		    end
		    Flux.update!(opt, ps, gs)
		end
		if (mod(ep,print_freq) == 0 )|| (ep == epochs - 1)
			println("epochs: ",ep," LR: ",new_lr, " Loss: ",training_loss)
	    end

	end
	return ps
end
function weight_uncertainty(model, cost, data, opt; callback =()->(return 0),
	prior_dist = Normal(0.0,10.0), σ_z = 10.0,	σ_m = 10.0, σ_p = 10.0,
	itr = 100, T=10, c=1, M=3, print_freq=1)

	W_swa, P, chn = subspace_inference(model, cost, data, opt,
	callback =()->(return 0),
	prior_dist = prior_dist, σ_z = σ_z,	σ_m = σ_m, σ_p = σ_p,
	itr = itr, T=T, c=c, M=M, print_freq=print_freq)
	n_samples = length(chn["z[1]"]);
	z_samples = Array{Float64}(undef,M,n_samples)
	for j in 1:M
		z_samples[j,:] = chn["z[$j]"].data
	end
	chn_weights =  P*z_samples .+ W_swa

	return chn_weights
end


text:
using SoapySDR, SoapyRTLSDR_jll

# Here we want to test the behavior in a loop, to ensure
# that we can block on buffer overflow conditions, and
# handle partial reads, measure latnecy, etc

function rapid_read()
    dev = Devices()[1]
    rx_chan = dev.rx
    rx_stream = SoapySDR.Stream(rx_chan)
    @show SoapySDR.mtu(rx_stream)
    SoapySDR.activate!(rx_stream)
    bufs = [SoapySDR.SampleBuffer(rx_stream, 10^6) for i = 1:2]
    @show bufs[1].packet_count
    @show bufs[2].packet_count
    flip = true
    while true
        # double buffer
        flip = !flip
        current_buff = bufs[Int(flip)+1]
        prev_buff = bufs[Int(!flip)+1]
        @assert length(current_buff.bufs[1])%rx_stream.mtu == 0

        read!(rx_stream, current_buff)

        # sanity checks?
        #nequal = 0
        #for i in eachindex(current_buff.bufs)
        #    nequal += Int(current_buff.bufs[1][i] == prev_buff.bufs[1][i])
        #end
        #@show current_buff.timens
        #@show nequal, current_buff.timens, delta_t
    end

end


text:
# Autogenerated wrapper script for tree_sitter_highlight_jll for x86_64-apple-darwin
export libtree_sitter_highlight

JLLWrappers.@generate_wrapper_header("tree_sitter_highlight")
JLLWrappers.@declare_library_product(libtree_sitter_highlight, "@rpath/libtree_sitter_highlight-5f224cd64cd26dec.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libtree_sitter_highlight,
        "lib/libtree_sitter_highlight.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
## Simulation Input: Initial Condition / Design parameters
using Plots
α       = 1
dim     = 2
N       = 3
s_Bias  = ComponentArray(δ = 0.01, n = 1, r_ref = 10E3, k = 9, m = 10)

V_M_0   = 300
γ_M_0 	= deg2rad(45)
χ_M_0_deg = 30
χ_f_d_deg = 180
σ_M_lim_deg = 60
χ_M_0 	= deg2rad(χ_M_0_deg)
χ_f_d 	= deg2rad(χ_f_d_deg)
σ_M_lim = deg2rad(σ_M_lim_deg)

function sim_plot(N, dim, α, χ_M_0, χ_f_d, σ_M_lim, s_Bias)
    if dim == 2
        (p_M_0, v_M_0)  = (zeros(dim), V_M_0*[sin(χ_M_0); cos(χ_M_0)])
        (p_T_0, v_T_0)  = ([5E3; 1E3], zeros(dim))
        γ_f_d = 0
    elseif dim == 3
        (p_M_0, v_M_0)  = (zeros(dim), V_M_0*[cos(γ_M_0)*sin(χ_M_0); cos(γ_M_0)*cos(χ_M_0); sin(γ_M_0)])
        (p_T_0, v_T_0)  = ([10E3; 5E3; 5E3], 100*[-1; 0; 0])
        γ_f_d = -0.01
    end
    v̂_f_d   = [cos(γ_f_d)*sin(χ_f_d); cos(γ_f_d)*cos(χ_f_d); sin(γ_f_d)]

    s_BPNG = BPNG(N, dim, α, σ_M_lim, v̂_f_d, Bias_zero, s_Bias)
    # Bias = Bias_zero, Bias_IACG_StationaryTarget, Bias_IACG_StationaryTarget_2D

    # Execute simulation
    df = main(p_M_0, v_M_0, p_T_0, v_T_0, s_BPNG; A_M_max = 100)

    ts = df.time
    p_Ms = df.sol |> Map(datum -> datum.pursuer.p) |> collect
    p_Ts = df.sol |> Map(datum -> datum.evador.p)  |> collect

    p_Ms = hcat(p_Ms...)'
    p_Ts = hcat(p_Ts...)';

    ## Plotting
    legend_string = ["Missile" "Target"]

    if dim == 2
        f = fig_print([p_Ms[:,1] p_Ts[:,1]], [p_Ms[:,2] p_Ts[:,2]], "", legend_string, "x [m]", "y [m]"; ar_val = :equal, save_file = 0, N_markers = 10)
        plot!(f, xlims=(0, 6E3), ylims=(-3E3, 3E3))

    elseif dim == 3
        f = plot([p_Ms[:,1] p_Ts[:,1]], [p_Ms[:,2] p_Ts[:,2]], [p_Ms[:,3] p_Ts[:,3]], label = legend_string)
    end

    f = fig_print([], [], "filename"; fig_handle = f, ar_val = :equal) 	# To just save the result
end

# Parallel execution
# Threads.nthreads()

# n   = 100 # number of scenarios
# Ns  = 1:n |> Map(i -> rand(1)[1]) |> collect
# Ns= 1:0.1:3
# sim_results = Ns |> Map(N -> main(N)) |> tcollect


text:
# ---
# title: 890. Find and Replace Pattern
# id: problem890
# author: Indigo
# date: 2021-02-01
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/find-and-replace-pattern/description/>
# hidden: true
# ---
#
# You have a list of `words` and a `pattern`, and you want to know which words
# in `words` matches the pattern.
#
# A word matches the pattern if there exists a permutation of letters `p` so
# that after replacing every letter `x` in the pattern with `p(x)`, we get the
# desired word.
#
# ( _Recall that a permutation of letters is a bijection from letters to
# letters: every letter maps to another letter, and no two letters map to the
# same letter._ )
#
# Return a list of the words in `words` that match the given pattern.
#
# You may return the answer in any order.
#
#
#
# **Example 1:**
#
#
#
#     Input: words = ["abc","deq","mee","aqq","dkd","ccc"], pattern = "abb"
#     Output: ["mee","aqq"]
#     Explanation: "mee" matches the pattern because there is a permutation {a -> m, b -> e, ...}.
#     "ccc" does not match the pattern because {a -> c, b -> c, ...} is not a permutation,
#     since a and b map to the same letter.
#
#
#
# **Note:**
#
#   * `1 <= words.length <= 50`
#   * `1 <= pattern.length = words[i].length <= 20`
#
#
## @lc code=start
using LeetCode

function find_and_replace_pattern(words::Vector{String}, pattern::String)
    function arr_pattern(s::String)
        arr_pt = Char[]
        dct = Dict{Char, Char}()
        cur_ch = 'a'
        for ch in s
            if !haskey(dct, ch)
                dct[ch] = cur_ch
                cur_ch += 1
            end
            push!(arr_pt, dct[ch])
        end
        arr_pt
    end
    res = String[]
    a_pt = arr_pattern(pattern)
    for w in words
        if a_pt == arr_pattern(w)
            push!(res, w)
        end
    end
    res
end
## @lc code=end


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

new_sis3316_data_table() = Table(
    daqevtno = Vector{Int32}(),
    bufferno = Vector{Int32}(),
    channel = Vector{Int32}(),
    timestamp = Vector{Int64}(),
    energy = Vector{Int32}(),
    trigmax = Vector{Int32}(),
    samples = VectorOfVectors{UInt16}(),
    mawvalues = VectorOfVectors{Int32}(),
)

const SIS3316Data = typeof(new_sis3316_data_table())
export SIS3316Data

SIS3316Data() = new_sis3316_data_table()

Base.@kwdef struct EvtFormat
    save_mawvals::Bool = false              # Save (Energy or Trigger) MAW Test Data
    save_energy::Bool = false               # Save energy values (at trigger and maximum during trigger window)
    save_ft_maw::Bool = false               # Save fast trigger MAW values (max value, value before Trigger, value with Trigger)
    save_acc_78::Bool = false               # Save accumulator values for gates 7 and 8
    save_ph_acc16::Bool = false             # Save peak height and accumulator values for gates 1,2, ..,6
    n_samples::Int = 0                      # Number of raw samples to save
    n_mawvals::Int = 0                      # Number of MAW test data words to save

    # require n_samples % 2 == 0
end
export EvtFormat

function raw_event_data_size(fmt::EvtFormat)
    @argcheck fmt.n_samples % 2 == 0
    nEvtWords =
        2 +
        (fmt.save_ph_acc16  ? 7 : 0) +
        (fmt.save_acc_78    ? 2 : 0) +
        (fmt.save_ft_maw    ? 3 : 0) +
        (fmt.save_energy    ? 2 : 0) +
        1 +
        div(fmt.n_samples, 2) +
        (fmt.save_mawvals ? fmt.n_mawvals : 0)

    nEvtWords * sizeof(Int32)
end
export raw_event_data_size

Base.@kwdef struct ReadoutConfig
    channels::Vector{Int} = Vector{Int}(all_channels)
    buf_fill_threshold::Float64 = 0.5
    readout_interval::Float64 = 2.0
end
export ReadoutConfig

"""
    SIS3316Digitizer

Struck SIS 3316 VME Digitizer.

Constructors:

```julia
SIS3316Digitizer(mem::SIS3316Memory)
SIS3316Digitizer(hostname::AbstractString, port::Integer = 0xE000)
```
"""
mutable struct SIS3316Digitizer{Mem<:SIS3316Memory}
    mem::Mem
    fifo_mem::NTuple{4,Mem}
    fifo_mem_lock::NTuple{4,ReentrantLock}
    events::Observable{SIS3316Data}
    evtrate::Observable{Vector{Float64}}
    readout_config::ReadoutConfig
    readout_reqs::Channel{Bool}
end

export SIS3316Digitizer

function SIS3316Digitizer(hostname::AbstractString, port::Integer = 0xE000)
    gw = SIS3316Gateway(hostname, port)
    mem = SIS3316Memory(gw)
    fifo_mem = (mem, mem, mem, mem)
    fifo_mem_lock = (ReentrantLock(), ReentrantLock(), ReentrantLock(), ReentrantLock())
    events = Observable(SIS3316Data())
    evtrate = Observable(fill(0.0, length(all_channels)))
    readout_config = ReadoutConfig()
    readout_reqs = Channel{Bool}(); close(readout_reqs)

    dev = SIS3316Digitizer(mem, fifo_mem, fifo_mem_lock, events, evtrate, readout_config, readout_reqs)
    dev.readout_reqs = readout_task(dev)

    finalizer(x -> close(x), dev)

    dev
end

Base.show(io::IO, dev::SIS3316Digitizer) = print(io, "SIS3316Digitizer($(repr(dev.mem)))")

@inline Base.Broadcast.broadcastable(dev::SIS3316Digitizer) = Ref(dev)

Base.isopen(dev::SIS3316Digitizer) = isopen(dev.mem)

function Base.close(dev::SIS3316Digitizer)
    close(dev.mem)
    for m in dev.fifo_mem
        close(m)
    end
end

function readout_task(dev::SIS3316Digitizer; sleep_interval_ns::Integer = 50 * 10^6)
    Channel{Bool}(0; taskref = nothing, spawn = true) do readout_reqs
        try
            @debug "SIS3316Digitizer readout task started"
            while isopen(dev)
                if take!(readout_reqs) == true
                    @info "Sampling starting"
                    sampling_loop(dev, readout_reqs; sleep_interval_ns = sleep_interval_ns)
                end
                @info "Sampling stopped"
            end
            @debug "SIS3316Digitizer readout task terminating"
        catch err
            @error err
            close(dev)
        end
    end
end

function sampling_loop(dev::SIS3316Digitizer, readout_reqs::Channel{Bool}; sleep_interval_ns::Integer = 50 * 10^6)
    nbuf::Int = 1
    last_check::UInt64 = time_ns()
    last_readout::UInt64 = last_check
    nbytes_filled_last = fill(0.0, length(all_channels))
    start_capture(dev)

    active::Bool = true
    while active && isopen(dev)
        if isready(readout_reqs) && take!(readout_reqs) == false
            stop_capture(dev)
            active = false
            #!!! TODO: Read out last buffer if necessary
        else
            bnkinfo = get_bank_info(dev, all_channels)
            curr_time = time_ns()

            delta_t_check = (curr_time - last_check) * 1e-9
            delta_t_readout = (curr_time - last_readout) * 1e-9
            nbytes_filled = map(x -> Int(x.sampling_nbytes), bnkinfo)
            evt_nbytes = map(x -> raw_event_data_size(x.event_format), bnkinfo)

            dev.evtrate[] = (nbytes_filled .- nbytes_filled_last) ./ evt_nbytes ./ delta_t_check

            readout_config = dev.readout_config
            buf_over_thresh = any(nbytes_filled ./ ch_buffer_size .>= readout_config.buf_fill_threshold)

            if buf_over_thresh || delta_t_readout > readout_config.readout_interval
                @debug "Swapping buffers and reading out"
                swap_banks(dev)
                events = readout_parsed(dev, nbuf, readout_config.channels)
                dev.events[] = events
                last_readout = curr_time
                nbuf += 1
            end
            last_check = curr_time
            nbytes_filled_last .= nbytes_filled
        end
    end
    nothing
end

function set_capture!(dev::SIS3316Digitizer, state::Bool)
    push!(dev.readout_reqs, state)
end
export set_capture!

function arm_bank(dev::SIS3316Digitizer, bank::Integer)
    @argcheck bank in all_banks
    # sync #!!!
    dev.mem[key_disarm_and_arm_bank_addr[bank]] = 1
end
export arm_bank

function swap_banks(dev::SIS3316Digitizer)
    # sync #!!!
    sampling_bank = dev.mem[actual_sample_address_regs[1].bank]
    arm_bank(dev, sampling_bank == 1 ? 2 : 1)
end
export swap_banks

function clear_and_arm(dev::SIS3316Digitizer)
    arm_bank(dev, 2)
    arm_bank(dev, 1)
end
export clear_and_arm

function disarm(dev)
    # sync #!!!
    dev.mem[key_disarm_addr] = 1
    nothing
end
export disarm

function start_capture(dev)
    @info "Starting capture"
    disarm(dev)
    clear_and_arm(dev)
end
export start_capture

function stop_capture(dev)
    @info "Stopping capture"
    disarm(dev)
end
export stop_capture

function event_format(dev::SIS3316Digitizer, ch::Integer)
    @argcheck ch in all_channels

    # TODO: Add support for extended sample length
    # TODO: Add support for averaging value data format

    mem = dev.mem
    group, grp_ch = fpga_num(ch), fpga_ch(ch)

    @mt_out_of_order begin
        cfg = mem[dataformat_config_regs[group].ch[grp_ch]]
        n_samples = mem[raw_data_buffer_config_regs[group].sample_length]
        n_mawvals = mem[maw_test_buffer_config_regs[group].buffer_len]
    end

    EvtFormat(
        save_mawvals = cfg.save_mawvals,
        save_energy = cfg.save_energy,
        save_ft_maw = cfg.save_ft_maw,
        save_acc_78 = cfg.save_acc_78,
        save_ph_acc16 = cfg.save_ph_acc16,
        n_samples = n_samples,
        n_mawvals = n_mawvals
    )
end
export event_format

Base.@kwdef struct BankInfo
    sampling_bank::Int = 0
    sampling_nbytes::Int = 0
    readout_bank::Int = 0
    readout_nbytes::Int = 0
    event_format::EvtFormat = EvtFormat()
end

export BankInfo

function BankInfo(dev::SIS3316Digitizer, ch::Integer)
    @mt_out_of_order begin
        sampling_bank::Int = dev.mem[actual_sample_address_regs[ch].bank]
        sampling_nbytes::Int = dev.mem[actual_sample_address_regs[ch].sample_addr]
        readout_bank::Int = dev.mem[previous_bank_sample_address_regs[ch].bank]
        readout_nbytes::Int = dev.mem[previous_bank_sample_address_regs[ch].sample_addr]
        evtfmt::EvtFormat = event_format(dev, ch)
    end

    BankInfo(sampling_bank, sampling_nbytes, readout_bank, readout_nbytes, evtfmt)
end

function get_bank_info(dev::SIS3316Digitizer, channels::AbstractVector{<:Integer} = all_channels)
    result = fetch.(map(ch -> @mt_async(BankInfo(dev, ch)), all_channels))
    convert(Vector{BankInfo}, result)::Vector{BankInfo}
end
function get_bank_info end

function print_bank_info(dev::SIS3316Digitizer)
    bankinfo = get_bank_info(dev, all_channels)
    for ch in all_channels
        bi = bankinfo[ch]
        raw_event_size = raw_event_data_size(bi.event_format)

        sampling_n_events = bi.sampling_nbytes / raw_event_size
        readout_n_events = bi.readout_nbytes / raw_event_size

        @info "Channel $ch: Bank $(bi.sampling_bank) armed: bank $(bi.sampling_nbytes) bytes (sampling_n_events), bank $(bi.readout_bank) readout: $(bi.readout_nbytes) bytes (bi.readout_n_events events)"
    end
end
export print_bank_info

function read_ch_bank_data(dev::SIS3316Digitizer, ch::Integer, bank::Integer, from::Integer = 0, n_bytes::Integer = ch_buffer_size)
    group = fpga_num(ch)
    lock(dev.fifo_mem_lock[group]) do
        read_ch_bank_data(dev.fifo_mem[group], ch::Integer, bank::Integer, from, n_bytes)
    end
end
export read_ch_bank_data

function readout_raw(dev::SIS3316Digitizer, channels::AbstractVector{<:Integer})
    binfo = get_bank_info(dev, channels)

    futures = map(channels) do ch
        @mt_async begin
            bi = binfo[ch]
            read_ch_bank_data(dev, ch, bi.readout_bank, 0, bi.readout_nbytes)
        end
    end

    fetch.(futures)
end
export readout_raw

function readout_parsed(dev::SIS3316Digitizer, bufno::Integer, channels::AbstractVector{<:Integer})
    binfo = get_bank_info(dev, channels)

    futures = map(channels) do ch
        @mt_async begin
            bi = binfo[ch]
            data = read_ch_bank_data(dev, ch, bi.readout_bank, 0, bi.readout_nbytes)
            parse_bank_data(data, ch, bufno, bi.event_format)
        end
    end

    vcat(fetch.(futures)...)::SIS3316Data
end
export readout_parsed

function parse_bank_data!(result::Table, data::AbstractVector{UInt8}, channel::Integer, bufno::Integer, fmt::EvtFormat)
    evt_nbytes = raw_event_data_size(fmt)
    total_nbytes = length(eachindex(data))
    mod(total_nbytes, evt_nbytes) == 0 || throw(ArgumentError("Length of data buffer not a multiple of event size"))
    nevents = div(total_nbytes, evt_nbytes)

    nmawvalues = fmt.save_mawvals ? fmt.n_mawvals : 0
    firmware_type = FW250 # ToDo: auto-detect

    tmpbuffer = Vector{Int32}()
    evtdatabuf = IOBuffer(data)
    evtno = 0

    for i in 1:nevents
        evtno += 1
        evt = read(evtdatabuf, RawChEvent, nmawvalues, firmware_type, tmpbuffer)

        push!(result.daqevtno, evtno)
        push!(result.channel, channel)
        push!(result.bufferno, bufno)
        push!(result.timestamp, evt.timestamp)
        push!(result.energy, isnothing(evt.energy) ? Int32(0) : evt.energy.maximum)
        push!(result.trigmax, isnothing(evt.trig_maw) ? Int32(0) : evt.trig_maw.maximum)
        push!(result.samples, evt.samples)
        push!(result.mawvalues, evt.mawvalues)
    end

    result
end
parse_bank_data!

function parse_bank_data(data::AbstractVector{UInt8}, channel::Integer, bufno::Integer, fmt::EvtFormat)
    parse_bank_data!(SIS3316Data(), data, channel, bufno, fmt)
end
export parse_bank_data


text:
using Transformers
using Transformers.Basic
using Test

import Flux

const tests = [
    "transformer",
    "nntopo",
    "embed",
    "basic",
    "gpt",
    "bert",
    "pretrain",
]

if v"1.0.0" <= VERSION  <= v"1.2.0"
    @test isempty(detect_ambiguities(Transformers, Transformers.GenerativePreTrain, Transformers.Basic, Transformers.Datasets))
end

@testset "Transformers" begin
    for t in tests
        fp = joinpath(dirname(@__FILE__), "test_$t.jl")
        @info "Test $(uppercase(t))"
        include(fp)
    end

    if Base.find_package("CuArrays") != nothing
        include("test_cuda.jl")
    end
end


text:
module SynapseTools

using Pkg

if !haskey(Pkg.installed(),"DISSEQT")
    warn("Module DISSEQT not installed. Please refer to SynapseTools installation instructions at https://github.com/rasmushenningsson/SynapseTools.jl")
elseif !haskey(Pkg.installed(),"SynapseClient")
    warn("Module SynapseClient not installed. Please refer to SynapseTools installation instructions at https://github.com/rasmushenningsson/SynapseTools.jl")
else

    using SynapseClient
    using DISSEQT
    using DataFrames
    using Missings

    import SynapseClient: AbstractEntity, Entity, File, Folder, Project, Activity

    export
        retrystore,
        getchildbyname,
        childpath,
        createchildfolder,
        uploadiflocal,
        listfiles,
        localpath,
        UploadList,
        markforupload!,
        uploadfiles,
        listaligned,
        referencefromlog,
        wrongreferencefromlog,
        getreferencegenomes,
        metadatatemplate,
        filtermetadata,
        getsamplemetadata,
        appendsynapseids!,
        appendswarmids!,
        appendalignedids!,
        downloadbymeta!,
        downloadswarms!,
        downloadconsensuses!,
        downloadbams!,
        downloadbais!,
        downloadaligned!,
        getsamplefitness,
        appendfitness!

    include("pattern.jl")
    include("paths.jl")
    include("samples.jl")
    include("referencegenomes.jl")
    include("metadata.jl")
    include("fitness.jl")
end

end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using Combinatorics
import JuLIP: evaluate

struct PureBasis{TB}
   basis::TB
   Ibas::Vector{Int}
end

Base.length(b::PureBasis) = length(b.Ibas)

(b::PureBasis)(args...) = evaluate(b, args...)

PureBasis(ship::SHIPBasis) = PureBasis(ship, collect(1:length(ship)))

function PureBasis(ship::SHIPBasis, N::Integer)
   @assert length(ship.zlist) == 1
   Ibas = ACE.Utils.findall_basis_N(ship, N)
   return PureBasis(ship, Ibas)
end

function evaluate(b::PureBasis, Rs::AbstractVector)
   N = length(Rs)
   ship = b.basis
   @assert N <= 5
   @assert length(ship.zlist) == 1
   Zs = zeros(Int16, N)
   tmp = alloc_temp(ship, Rs)

   # one-particle basis
   phi = precompute_phi(tmp, Rs, Zs, ship)

   # permutation symmetrised tensor products
   Phi = precompute_prod_phi(tmp, phi, ship)

   # get B via the same rotation-symmetrisation as in the density trick case
   _my_mul!(tmp.Bc[1], ship.A2B[1], Phi)
   return real.(tmp.Bc[1])[b.Ibas]
end

function precompute_phi(tmp, Rs, Zs, ship)
   alist = ship.alists[1]
   phi = zeros(Complex{Float64}, length(alist), length(Rs))
   for (iR, (R, Z)) in enumerate(zip(Rs, Zs))
      # evaluate the r-basis and the R̂-basis for the current neighbour at R
      evaluate!(tmp.J, tmp.tmpJ, ship.J, norm(R))
      evaluate!(tmp.Y, tmp.tmpY, ship.SH, R)
      # write the products into phi
      iz = z2i(ship, Z)
      for i = alist.firstz[iz]:(alist.firstz[iz+1]-1)
         zklm = alist[i]
         phi[i, iR] = tmp.J[zklm.k+1] * tmp.Y[index_y(zklm.l, zklm.m)]
      end
   end
   return phi
end

function precompute_prod_phi(tmp, phi, ship)
   aalist = ship.aalists[1]
   Phi = fill!(tmp.AA[1], 0)
   for i = 1:length(aalist)
      Phi[i] = _VN(phi, (@view aalist.i2Aidx[i, 1:aalist.len[i]]))
   end
   return Phi
end

function _VN(phi, phiidx)
   nR = size(phi, 2)     # number of neighbours
   N  = length(phiidx)   # order of the basis function
   Phi = zero(ComplexF64)
   for A in combinations(1:nR, N), σA in permutations(A)
      Phi += prod( phi[phiidx[α], σα]  for (α, σα) in enumerate(σA) )
   end
   return Phi / factorial(N)
end


text:
struct CacheKey{C <: AbstractCache, D <: PSY.Device}
    cache_type::Type{C}
    device_type::Type{D}
end

function CacheKey(cache::C) where {C <: AbstractCache}
    return CacheKey(C, cache.device_type)
end
# The current implementation will require all custom caches to have the same data template
"""
Tracks the last time status of a device changed in a simulation
"""
mutable struct TimeStatusChange <: AbstractCache
    device_type::Type{<:PSY.Device}
    value::JuMP.Containers.DenseAxisArray{Dict{Symbol, Float64}}
    units::Dates.TimePeriod
    ref::UpdateRef

    function TimeStatusChange(
        device_type::Type{<:PSY.Device},
        value::JuMP.Containers.DenseAxisArray{Dict{Symbol, Float64}},
        ref::UpdateRef,
        units::Dates.TimePeriod = Dates.Hour(1),
    )
        units = IS.time_period_conversion(units)
        new(device_type, value, units, ref)
    end
end

function TimeStatusChange(::Type{T}, name::AbstractString) where {T <: PSY.Device}
    value_array = JuMP.Containers.DenseAxisArray{Dict{Symbol, Float64}}(undef, 1)
    return TimeStatusChange(T, value_array, UpdateRef{JuMP.VariableRef}(T, name))
end

mutable struct StoredEnergy <: AbstractCache
    device_type::Type{<:PSY.Device}
    value::JuMP.Containers.DenseAxisArray{Float64}
    ref::UpdateRef
end

function StoredEnergy(::Type{T}, name::AbstractString) where {T <: PSY.Device}
    value_array = JuMP.Containers.DenseAxisArray{Float64}(undef, 1)
    return StoredEnergy(T, value_array, UpdateRef{JuMP.VariableRef}(T, name))
end

cache_value(cache::AbstractCache, key) = cache.value[key]


text:
# # Polynomial Optimization

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Polynomial Optimization/polynomial_optimization.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Polynomial Optimization/polynomial_optimization.ipynb)
# **Contributed by**: Benoît Legat

# ## Introduction

# Consider the polynomial optimization problem [L09, Example 2.2] of
# minimizing the polynomial $x^3 - x^2 + 2xy - y^2 + y^3$
# over the polyhedron defined by the inequalities $x \ge 0, y \ge 0$ and $x + y \geq 1$.

# [L09] Lasserre, J. B.
# *Moments, positive polynomials and their applications*.
# World Scientific, **2009**.

using Test #src
using DynamicPolynomials
@polyvar x y
p = x^3 - x^2 + 2x*y -y^2 + y^3
using SumOfSquares
S = @set x >= 0 && y >= 0 && x + y >= 1
p(x=>1, y=>0), p(x=>1//2, y=>1//2), p(x=>0, y=>1)

# The optimal solutions are $(x, y) = (1, 0)$ and $(x, y) = (0, 1)$ with objective value $0$ but [Ipopt](https://github.com/jump-dev/Ipopt.jl/) only finds the local minimum $(1/2, 1/2)$ with objective value $1/4$.

import Ipopt
model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
@variable(model, a >= 0)
@variable(model, b >= 0)
@constraint(model, a + b >= 1)
@NLobjective(model, Min, a^3 - a^2 + 2a*b - b^2 + b^3)
optimize!(model)
@test termination_status(model) == MOI.LOCALLY_SOLVED #src
@show termination_status(model)
@test value(a) ≈ 0.5 rtol=1e-5 #src
@show value(a)
@test value(b) ≈ 0.5 rtol=1e-5 #src
@show value(b)
@test objective_value(model) ≈ 0.25 rtol=1e-5 #src
@show objective_value(model)

# Note that the problem can be written equivalently as follows using [registered functions](https://jump.dev/JuMP.jl/v0.21/nlp/#User-defined-Functions-1).

using Ipopt
model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
@variable(model, a >= 0)
@variable(model, b >= 0)
@constraint(model, a + b >= 1)
peval(a, b) = p(x=>a, y=>b)
register(model, :peval, 2, peval, autodiff=true)
@NLobjective(model, Min, peval(a, b))
optimize!(model)
@test termination_status(model) == MOI.LOCALLY_SOLVED #src
@show termination_status(model)
@test value(a) ≈ 0.5 rtol=1e-5 #src
@show value(a)
@test value(b) ≈ 0.5 rtol=1e-5 #src
@show value(b)
@test objective_value(model) ≈ 0.25 rtol=1e-5 #src
@show objective_value(model)

## Sum-of-Squares approach

# We will now see how to find the optimal solution using Sum of Squares Programming.
# We first need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v0.21.6/installation/#Supported-solvers) for a list of the available choices.

import CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)

# A Sum-of-Squares certificate that $p \ge \alpha$ over the domain `S`, ensures that $\alpha$ is a lower bound to the polynomial optimization problem.
# The following program searches for the largest upper bound and finds zero.

model = SOSModel(solver)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c3, p >= α, domain = S)
optimize!(model)
@test termination_status(model) == MOI.OPTIMAL #src
@show termination_status(model)
@test objective_value(model) ≈ 0.0 atol=1e-5 #src
@show objective_value(model)

# Using the solution $(1/2, 1/2)$ found by Ipopt of objective value $1/4$
# and this certificate of lower bound $0$ we know that the optimal objective value is in the interval $[0, 1/4]$
# but we still do not know what it is (if we consider that we did not try the solutions $(1, 0)$ and $(0, 1)$ as done in the introduction).
# If the dual of the constraint `c3` was atomic, its atoms would have given optimal solutions of objective value $0$ but that is not the case.

ν3 = moment_matrix(c3)
@test extractatoms(ν3, 1e-3) === nothing #src
extractatoms(ν3, 1e-3) # Returns nothing as the dual is not atomic

# Fortunately, there is a hierarchy of programs with increasingly better bounds that can be solved until we get one with atom dual variables.
# This comes from the way the Sum-of-Squares constraint with domain `S` is formulated.
# The polynomial $p - \alpha$ is guaranteed to be nonnegative over the domain `S` if there exists Sum-of-Squares polynomials $s_0$, $s_1$, $s_2$, $s_3$ such that
# ```math
# p - \alpha = s_0 + s_1 x + s_2 y + s_3 (x + y - 1).
# ```
# Indeed, in the domain `S`, $x$, $y$ and $x + y - 1$ are nonnegative so the right-hand side is a sum of squares hence is nonnegative.
# Once the degrees of $s_1$, $s_2$ and $s_3$ have been decided, the degree needed for $s_0$ will be determined but we have a freedom in choosing the degrees of $s_1$, $s_2$ and $s_3$.
# By default, they are chosen so that the degrees of $s_1 x$, $s_2 y$ and $s_3 (x + y - 1)$ match those of $p - \alpha$ but this can be overwritten using the `maxdegree` keyword argument.

# ### The maxdegree keyword argument

# The maximum total degree (i.e. maximum sum of the exponents of $x$ and $y$) of the monomials of $p$ is 3 so the constraint in the program above is equivalent to `@constraint(model, p >= α, domain = S, maxdegree = 3)`.
# That is, since $x$, $y$ and $x + y - 1$ have total degree 1, the sum of squares polynomials $s_1$, $s_2$ and $s_3$ have been chosen with maximum total degree $2$.
# Since these polynomials are sums of squares, their degree must be even so the next maximum total degree to try is 4.
# For this reason, the keywords `maxdegree = 4` and `maxdegree = 5` have the same effect in this example.
# In general, if the polynomials in the domain are not all odd or all even, each value of `maxdegree` has a different effect in the choice of the maximum total degree of some $s_i$.

model = SOSModel(solver)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c5, p >= α, domain = S, maxdegree = 5)
optimize!(model)
@test termination_status(model) == MOI.OPTIMAL #src
@show termination_status(model)
@test objective_value(model) ≈ 0.0 atol=1e-5 #src
@show objective_value(model)

# As shown below, for `maxdegree = 5`, the dual variable is atomic as it is the moments of the measure
# $$0.5 \delta(x-1, y) + 0.5 \delta(x, y-1)$$
# where $\delta(x, y)$ is the dirac measure centered at $(0, 0)$.
# Therefore the program provides both a certificate that $0$ is a lower bound and a certificate that it is also an upper bound since it is attained at the global minimizers $(1, 0)$ and $(0, 1)$.

ν5 = moment_matrix(c5)
atoms5 = extractatoms(ν5, 1e-3) #src
@test atoms5.atoms[1].weight ≈ 0.5 rtol=1e-2 #src
@test atoms5.atoms[2].weight ≈ 0.5 rtol=1e-2 #src
@test atoms5.atoms[1].center[2:-1:1] ≈ atoms5.atoms[2].center[1:2] rtol=1e-2 #src
extractatoms(ν5, 1e-3)

# ## A deeper look into atom extraction

# The `extractatoms` function requires a `ranktol` argument that we have set to `1e-3` in the preceding section.
# This argument is used to transform the dual variable into a system of polynomials equations whose solutions give the atoms.
# This transformation uses the SVD decomposition of the moment matrix and discards the equations corresponding to a singular value lower than `ranktol`.
# When this system of equation has an infinite number of solutions, `extractatoms` concludes that the measure is not atomic.
# For instance, with `maxdegree = 3`, we obtain the system
# $$x + y = 1$$
# which contains a whole line of solution.
# This explains `extractatoms` returned `nothing`.

ν3 = moment_matrix(c3)
SumOfSquares.MultivariateMoments.computesupport!(ν3, 1e-3)
@test length(ν3.support.I.p) == 1 #src

# With `maxdegree = 5`, we obtain the system
# ```math
# \begin{aligned}
#   x + y & = 1\\
#   y^2 & = y\\
#   xy & = 0\\
#   x^2 + y & = 1
# \end{aligned}
# ```

ν5 = moment_matrix(c5)
SumOfSquares.MultivariateMoments.computesupport!(ν5, 1e-3)

# This system can be reduced to the equivalent system
# ```math
# \begin{aligned}
#   x + y & = 1\\
#   y^2 & = y
# \end{aligned}
# ```
# which has the solutions $(0, 1)$ and $(1, 0)$.

SemialgebraicSets.computegröbnerbasis!(ideal(ν5.support))
ν5.support
@test length(ν5.support.I.p) == 2 #src

# The function `extractatoms` then reuses the matrix of moments to find the weights $1/2$, $1/2$ corresponding to the diracs centered respectively at $(0, 1)$ and $(1, 0)$.
# This details how the function obtained the result
# $$0.5 \delta(x-1, y) + 0.5 \delta(x, y-1)$$
# given in the previous section.

# ## HomotopyContinuation

# As discussed in the previous section, the atom extraction relies on the solution
# of a system of algebraic equations. The `extractatoms` function takes an optional
# `solver` argument that is used to solve this system of equation.
# If no solver is provided, the default solver of SemialgebraicSets.jl is used which
# currently computes the Gröbner basis, then the multiplication matrices and
# then the Schur decomposition of a random combination of these matrices.
# As the system of equations is obtained from a numerical solution and is represented
# using floating point coefficients, homotopy continuation is recommended as it is
# more numerically robust than Gröbner basis computation.
# The following uses homotopy continuation to solve the system of equations.

using HomotopyContinuation
solver = SemialgebraicSetsHCSolver(; excess_residual_tol = 2e-2, real_tol = 2e-2, compile = false)
atoms5 = extractatoms(ν5, 1e-3, solver) #src
@test length(atoms5.atoms) == 2 #src
@test atoms5.atoms[1].weight + atoms5.atoms[2].weight ≈ 1.0 rtol=1e-2 #src
@test atoms5.atoms[1].center[2:-1:1] ≈ atoms5.atoms[2].center[1:2] rtol=1e-2 #src
extractatoms(ν5, 1e-3, solver)

# As the system has 3 equations for 2 variables and the coefficients of the equations
# are to be treated with tolerance since they originate from the solution of an SDP,
# we need to set `excess_residual_tol` and `real_tol` to a high tolerance otherwise,
# HomotopyContinuation would consider that there is no solution.
# Indeed, as the system is overdetermined (it has more equations than variables)
# HomotopyContinuation expects to have excess solution hence it filters out
# excess solution among the solution found. It determines which solution are in excess
# by comparing the infinity norm of the residuals of the equations at the solution with `excess_residual_tol`.
# It also filters out solution for which the absolute value of the imaginary part of one of the entry
# is larger than `real_tol` and strips out the imaginary part.
# The raw solutions obtained by HomotopyContinuation can be obtained as follows:

F = HomotopyContinuation.System(ν5.support)
res = HomotopyContinuation.solve(F, solver.options...)
r = path_results(res) #src
@test length(r) == 4 #src
@test all(HomotopyContinuation.is_excess_solution, r) #src
path_results(res)

# The printed `residual` above shows why `2e-2` allows to filter how the 2 actual
# solutions from the 2 excess solutions.

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
#=
Create a tiny image as test data
=#

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

function createTestMaskedImage()
    imagePath = "/home/bootch/git/resynthesizerJulia/test/data/in/tinylighthouse.png"
    image = load(imagePath)

    # The mask is more or less in the middle, with a frisket
    mask = falses(size(image))
    mask[3:4, 2:3] .= true

    # @printf("Image length %d\n", length(image))
    testMaskedImage = MaskedImage(image, mask)
    return testMaskedImage
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
#### Precipitation component in atmosphere model
abstract type PrecipitationModel end

export NoPrecipitation, RainModel, RainSnowModel

eq_tends(pv::PV, m::PrecipitationModel, ::Flux{O}) where {PV, O} = ()

using ..Microphysics

vars_state(::PrecipitationModel, ::AbstractStateType, FT) = @vars()

function atmos_nodal_update_auxiliary_state!(
    ::PrecipitationModel,
    m::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function flux_first_order!(
    ::PrecipitationModel,
    atmos::AtmosModel,
    flux::Grad,
    args,
) end
function compute_gradient_flux!(
    ::PrecipitationModel,
    diffusive,
    ∇transform,
    state,
    aux,
    t,
) end
function flux_second_order!(
    ::PrecipitationModel,
    flux::Grad,
    ::AtmosModel,
    args,
) end
function flux_first_order!(::PrecipitationModel, _...) end
function compute_gradient_argument!(
    ::PrecipitationModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end

source!(::PrecipitationModel, args...) = nothing

"""
    NoPrecipitation <: PrecipitationModel

No precipitation.
"""
struct NoPrecipitation <: PrecipitationModel end

"""
    RainModel <: PrecipitationModel

Precipitation model with rain.
"""
struct RainModel <: PrecipitationModel end

vars_state(::RainModel, ::Prognostic, FT) = @vars(ρq_rai::FT)
vars_state(::RainModel, ::Gradient, FT) = @vars(q_rai::FT)
vars_state(::RainModel, ::GradientFlux, FT) = @vars(∇q_rai::SVector{3, FT})

function atmos_nodal_update_auxiliary_state!(
    precip::RainModel,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end

function compute_gradient_argument!(
    precip::RainModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    transform.precipitation.q_rai = state.precipitation.ρq_rai * ρinv
end

function compute_gradient_flux!(
    precip::RainModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.precipitation.∇q_rai = ∇transform.precipitation.q_rai
end

function flux_first_order!(
    precip::RainModel,
    atmos::AtmosModel,
    flux::Grad,
    args,
)
    tend = Flux{FirstOrder}()
    flux.precipitation.ρq_rai =
        Σfluxes(eq_tends(Rain(), atmos, tend), atmos, args)
end

function flux_second_order!(
    precip::RainModel,
    flux::Grad,
    atmos::AtmosModel,
    args,
)
    tend = Flux{SecondOrder}()
    flux.precipitation.ρq_rai =
        Σfluxes(eq_tends(Rain(), atmos, tend), atmos, args)
end

function source!(m::RainModel, source::Vars, atmos::AtmosModel, args)
    tend = Source()
    source.precipitation.ρq_rai =
        Σsources(eq_tends(Rain(), atmos, tend), atmos, args)
end

"""
    RainSnowModel <: PrecipitationModel

Precipitation model with rain and snow.
"""
struct RainSnowModel <: PrecipitationModel end

vars_state(::RainSnowModel, ::Prognostic, FT) = @vars(ρq_rai::FT, ρq_sno::FT)
vars_state(::RainSnowModel, ::Gradient, FT) = @vars(q_rai::FT, q_sno::FT)
vars_state(::RainSnowModel, ::GradientFlux, FT) =
    @vars(∇q_rai::SVector{3, FT}, ∇q_sno::SVector{3, FT})

function atmos_nodal_update_auxiliary_state!(
    precip::RainSnowModel,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end

function compute_gradient_argument!(
    precip::RainSnowModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    transform.precipitation.q_rai = state.precipitation.ρq_rai * ρinv
    transform.precipitation.q_sno = state.precipitation.ρq_sno * ρinv
end

function compute_gradient_flux!(
    precip::RainSnowModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.precipitation.∇q_rai = ∇transform.precipitation.q_rai
    diffusive.precipitation.∇q_sno = ∇transform.precipitation.q_sno
end

function flux_first_order!(
    precip::RainSnowModel,
    atmos::AtmosModel,
    flux::Grad,
    args,
)
    tend = Flux{FirstOrder}()
    flux.precipitation.ρq_rai =
        Σfluxes(eq_tends(Rain(), atmos, tend), atmos, args)
    flux.precipitation.ρq_sno =
        Σfluxes(eq_tends(Snow(), atmos, tend), atmos, args)
end

function flux_second_order!(
    precip::RainSnowModel,
    flux::Grad,
    atmos::AtmosModel,
    args,
)
    tend = Flux{SecondOrder}()
    flux.precipitation.ρq_rai =
        Σfluxes(eq_tends(Rain(), atmos, tend), atmos, args)
    flux.precipitation.ρq_sno =
        Σfluxes(eq_tends(Snow(), atmos, tend), atmos, args)
end

function source!(m::RainSnowModel, source::Vars, atmos::AtmosModel, args)
    tend = Source()
    source.precipitation.ρq_rai =
        Σsources(eq_tends(Rain(), atmos, tend), atmos, args)
    source.precipitation.ρq_sno =
        Σsources(eq_tends(Snow(), atmos, tend), atmos, args)
end

#####
##### Tendency specifications
#####

eq_tends(pv::PV, ::RainModel, ::Flux{FirstOrder}) where {PV <: Rain} =
    (PrecipitationFlux{PV}(),)
eq_tends(pv::PV, ::RainModel, ::Flux{SecondOrder}) where {PV <: Rain} =
    (Diffusion{PV}(),)

eq_tends(
    pv::PV,
    ::RainSnowModel,
    ::Flux{FirstOrder},
) where {PV <: Precipitation} = (PrecipitationFlux{PV}(),)
eq_tends(
    pv::PV,
    ::RainSnowModel,
    ::Flux{SecondOrder},
) where {PV <: Precipitation} = (Diffusion{PV}(),)

max_stars_repo_name: jm-c/CLIMA

max_issues_repo_name: jm-c/CLIMA

max_forks_repo_name: jm-c/CLIMA


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

+(a::IntervalBox, b::IntervalBox) = IntervalBox( a.v .+ b.v )
+(a::IntervalBox, b::Real) = IntervalBox( a.v .+ b )
+(a::Real, b::IntervalBox) = IntervalBox( a .+ b.v )

-(a::IntervalBox, b::IntervalBox) = IntervalBox( a.v .- b.v )
-(a::IntervalBox, b::Real) = IntervalBox( a.v .- b )
-(a::Real, b::IntervalBox) = IntervalBox( a .- b.v )
-(a::IntervalBox) = IntervalBox( .- a.v )

*(a::IntervalBox, b::Real) = IntervalBox( a.v .* b )
*(a::Real, b::IntervalBox) = IntervalBox( a .* b.v )

/(a::IntervalBox, b::Real) = IntervalBox( a.v ./ b )

# broadcasting:

# wrap decides whether to wrap the result in an IntervalBox or not, based on the return type
wrap(v::SVector{N,T} where {N,T<:Interval}) = IntervalBox(v)
wrap(v) = v

@inline Base.broadcast(f, X::IntervalBox) = wrap(f.(X.v))
@inline Base.broadcast(f, X::IntervalBox, Y::IntervalBox) = wrap(f.(X.v, Y.v))
@inline Base.broadcast(f, X::IntervalBox, y) = wrap(f.(X.v, y))

for op in (:+, :-, :∩, :∪, :⊆, :isinterior, :dot, :setdiff, :×)
    @eval $(op)(a::SVector, b::IntervalBox) = $(op)(IntervalBox(a), b)
end

for op in (:+, :-, :∩, :∪, :⊆, :isinterior, :dot, :setdiff, :×)
    @eval $(op)(a::IntervalBox, b::SVector) = $(op)(a, IntervalBox(b))
end


text:
s1 = "I am a string."
s2 = """I am also a "string". """

# Single quotes are for characters
typeof('a')

name = "Jane"
num_fingers = 10
num_toes = 10

# String interpolation
println("Hello, my name is $name.")
println("I have $num_fingers fingers and $num_toes toes, which is $(num_fingers + num_toes) total digits.")

# String concatenation
string(s1, " ", s2)
string(s1, " ", s2, 210)

s3 = s1*s2
s4 = "$s1$s2"

max_stars_repo_name: AustinTSchaffer/DailyProgrammer

max_issues_repo_name: AustinTSchaffer/DailyProgrammer

max_forks_repo_name: AustinTSchaffer/DailyProgrammer


text:
using AdaptiveSampling
using Test

@testset "AdaptiveSampling.jl" begin
    fs = 5e6
    dt = 1/fs
    t = collect(0:dt:0.00004)
    f = 1e5
    y = sin.(2*pi*f*t)
    y[10:12]   .= 5;
    detail_th  = 0.1
    @time tc,yc, idx2save = AdaptiveSampling.anchor_point(t,y,detail_th)
    @test length(tc) < length(t) # compressed signal should be smaller than original
    @info "compression_ratio" compression_ratio = length(t)/length(tc)

    @time tc,yc, idx2save = AdaptiveSampling.anchor_point(t,y,detail_th,true)
    @test length(tc) < length(t) # compressed signal should be smaller than original
    @info "compression_ratio" compression_ratio = length(t)/length(tc)

    @time tc,yc, idx2save = AdaptiveSampling.anchor_point(t,y,detail_th,false)
    @test length(tc) < length(t) # compressed signal should be smaller than original
    @info "compression_ratio" compression_ratio = length(t)/length(tc)

    # @time tc,yc, idx2save, anim = AdaptiveSampling.anchor_point_gif(t,y,thr)
    # @test length(tc) < length(t) # compressed signal should be smaller than original
    # @info "compression_ratio" compression_ratio = length(t)/length(tc)
end


text:
function add_re!(sys)
    re = RenewableDispatch(
        "WindBusA",
        true,
        get_component(Bus, sys, "bus5"),
        0.0,
        0.0,
        1.200,
        PrimeMovers.WT,
        (min = 0.0, max = 0.0),
        1.0,
        TwoPartCost(0.220, 0.0),
        100.0,
    )
    add_component!(sys, re)
    copy_time_series!(re, get_component(PowerLoad, sys, "bus2"))

    for g in get_components(HydroEnergyReservoir, sys)
        tpc = get_operation_cost(g)
        smc = StorageManagementCost(
            variable = get_variable(tpc),
            fixed = get_fixed(tpc),
            start_up = 0.0,
            shut_down = 0.0,
            energy_shortage_cost = 10.0,
            energy_surplus_cost = 10.0,
        )
        set_operation_cost!(g, smc)
    end
end

function run_test_sim(result_dir::String)
    sim_name = "results_sim"
    sim_path = joinpath(result_dir, sim_name)

    if ispath(sim_path)
        @info "Reading UC system from" sim_path
        c_sys5_hy_uc = System(
            joinpath(sim_path, "..", "c_sys5_hy_uc.json"),
            time_series_read_only = true,
        )
        @info "Reading ED system from" sim_path
        c_sys5_hy_ed = System(
            joinpath(sim_path, "..", "c_sys5_hy_ed.json"),
            time_series_read_only = true,
        )
        executions = tryparse.(Int, readdir(sim_path))
        sim = joinpath(sim_path, string(maximum(executions)))
        @info "Reading results from last execution" sim
    else
        @info "Building UC system from"
        c_sys5_hy_uc = PSB.build_system(PSB.SIIPExampleSystems, "5_bus_hydro_uc_sys")
        @info "Building ED system from"
        c_sys5_hy_ed = PSB.build_system(PSB.SIIPExampleSystems, "5_bus_hydro_ed_sys")

        @info "Adding extra RE"
        add_re!(c_sys5_hy_uc)
        add_re!(c_sys5_hy_ed)
        to_json(c_sys5_hy_uc, joinpath(sim_path, "..", "c_sys5_hy_uc.json"))
        to_json(c_sys5_hy_ed, joinpath(sim_path, "..", "c_sys5_hy_ed.json"))

        mkpath(result_dir)
        GLPK_optimizer =
            optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => GLPK.GLP_MSG_OFF)

        template_hydro_st_uc = template_unit_commitment()
        set_device_model!(template_hydro_st_uc, HydroDispatch, FixedOutput)
        set_device_model!(
            template_hydro_st_uc,
            HydroEnergyReservoir,
            HydroDispatchReservoirStorage,
        )

        template_hydro_st_ed = template_economic_dispatch()
        set_device_model!(template_hydro_st_ed, HydroDispatch, FixedOutput)
        set_device_model!(
            template_hydro_st_ed,
            HydroEnergyReservoir,
            HydroDispatchReservoirStorage,
        )
        template_hydro_st_ed.services = Dict() #remove ed services

        problems = SimulationProblems(
            UC = OperationsProblem(
                template_hydro_st_uc,
                c_sys5_hy_uc,
                optimizer = GLPK_optimizer,
                system_to_file = false,
            ),
            ED = OperationsProblem(
                template_hydro_st_ed,
                c_sys5_hy_ed,
                optimizer = GLPK_optimizer,
                constraint_duals = [:CopperPlateBalance],
                system_to_file = false,
                balance_slack_variables = true,
            ),
        )

        sequence = SimulationSequence(
            problems = problems,
            feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
            intervals = Dict(
                "UC" => (Hour(24), Consecutive()),
                "ED" => (Hour(1), Consecutive()),
            ),
            feedforward = Dict(
                ("ED", :devices, :ThermalStandard) => SemiContinuousFF(
                    binary_source_problem = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
                ("ED", :devices, :HydroEnergyReservoir) => IntegralLimitFF(
                    variable_source_problem = PSI.ACTIVE_POWER,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            cache = Dict(
                ("UC",) => TimeStatusChange(PSY.ThermalStandard, PSI.ON),
                ("UC", "ED") => StoredEnergy(PSY.HydroEnergyReservoir, PSI.ENERGY),
            ),
            ini_cond_chronology = InterProblemChronology(),
        )
        sim = Simulation(
            name = "results_sim",
            steps = 2,
            problems = problems,
            sequence = sequence,
            simulation_folder = result_dir,
        )
        build!(sim)
        execute!(sim)
    end

    results = SimulationResults(sim)
    results_uc = get_problem_results(results, "UC")
    set_system!(results_uc, c_sys5_hy_uc)
    results_ed = get_problem_results(results, "ED")
    set_system!(results_ed, c_sys5_hy_ed)

    return results_uc, results_ed
end

function run_test_prob()
    c_sys5_hy_uc = PSB.build_system(PSB.SIIPExampleSystems, "5_bus_hydro_uc_sys")
    add_re!(c_sys5_hy_uc)
    GLPK_optimizer =
        optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => GLPK.GLP_MSG_OFF)

    template_hydro_st_uc = template_unit_commitment()
    set_device_model!(template_hydro_st_uc, HydroDispatch, FixedOutput)
    set_device_model!(
        template_hydro_st_uc,
        HydroEnergyReservoir,
        HydroDispatchReservoirStorage,
    )

    prob = OperationsProblem(
        template_hydro_st_uc,
        c_sys5_hy_uc,
        optimizer = GLPK_optimizer,
        horizon = 12,
        use_parameters = true,
    )
    build!(prob, output_dir = mktempdir())
    solve!(prob)
    res = ProblemResults(prob)
    return res
end


text:
import PWS: DirectMCEstimate, chemotaxis_system, ConditionalEnsemble, MarginalEnsemble, marginal_configuration, generate_configuration, collect_samples, mutual_information
using Test
using Statistics

system = chemotaxis_system()

algorithm = DirectMCEstimate(1000)
dtimes = collect(range(0.0, 2.0, length=51)[2:end])
num_responses = 10

result = mutual_information(system, algorithm, num_responses=num_responses, dtimes=dtimes)

mi = hcat(result.MutualInformation...)

@test all(isfinite.(mi))

mean(mi, dims=2)


text:
using FromFile
using Test

@from "test/file.jl" using A

path = abspath("test/file.jl")
m = getfield(From.__toplevel__, Symbol(path))


text:
import Base: -, +, *, filter, real, complex
import LinearAlgebra: norm, promote_leaf_eltypes

@inline +(φ1::T, φ2::T) where {T <: AbstractProperty} = T( φ1.val + φ2.val )
@inline -(φ1::T, φ2::T) where {T <: AbstractProperty} = T( φ1.val - φ2.val )
@inline -(φ::T) where {T <: AbstractProperty} = T( -φ.val)
@inline *(a::Union{Number, AbstractMatrix}, φ::T) where {T <: AbstractProperty} =
      T(a * φ.val)
@inline *(φ::T, a::Union{Number, AbstractMatrix}) where {T <: AbstractProperty} =
      T(φ.val * a)
@inline norm(φ::AbstractProperty) = norm(φ.val)
@inline Base.length(φ::AbstractProperty) = length(φ.val)
@inline Base.size(φ::AbstractProperty) = size(φ.val)
@inline Base.zero(φ::T) where {T <: AbstractProperty} = T(zero(φ.val))
@inline Base.zero(::Type{T}) where {T <: AbstractProperty} = zero(T())

promote_leaf_eltypes(φ::T) where {T <: AbstractProperty} =
      promote_leaf_eltypes(φ.val)

Base.convert(T::Type{TP}, φ::TP) where {TP <: AbstractProperty} = φ
Base.convert(T::Type, φ::AbstractProperty) = convert(T, φ.val)
Base.convert(T::Type{Any}, φ::AbstractProperty) = φ

# some type piracy ...
# TODO: hack like this make #27 important!!!
# *(a::SArray{Tuple{L1,L2,L3}}, b::SVector{L3}) where {L1, L2, L3} =
#       reshape( reshape(a, L1*L2, L3) * b, L1, L2)

coco_o_daa(φ::AbstractProperty, b::SVector) = coco_o_daa(φ.val, b)
coco_o_daa(cc::Number, b::SVector) = cc * b
coco_o_daa(cc::SVector, b::SVector) = cc * transpose(b)
coco_o_daa(cc::SMatrix{N1,N2}, b::SVector{N3}) where {N1,N2,N3} =
		reshape(cc[:] * transpose(b), Size(N1, N2, N3))
coco_o_daa(cc::SArray{Tuple{N1,N2,N3}}, b::SVector{N4}) where {N1,N2,N3,N4} =
		reshape(cc[:] * transpose(b), Size(N1, N2, N3, N4))

Base.isapprox(φ1::T, φ2::T) where {T <: AbstractProperty} =
      isapprox(φ1.val, φ2.val)

"""
`struct Invariant{D}` : specifies that the output of an ACE is
an invariant scalar.
"""
struct Invariant{T} <: AbstractProperty
   val::T
end

Invariant{T}() where {T <: Number} = Invariant{T}(zero(T))

Invariant(T::DataType = Float64) = Invariant{T}()

real(φ::Invariant) = Invariant(real(φ.val))
complex(φ::Invariant) = Invariant(complex(φ.val))
complex(::Type{Invariant{T}}) where {T} = Invariant{complex(T)}

filter(φ::Invariant, b::Array) = ( length(b) <= 1 ? true :
     iseven(sum(bi.l for bi in b)) && iszero(sum(bi.m for bi in b))  )

rot3Dcoeffs(::Invariant, T=Float64) = Rot3DCoeffs(T)

# TODO: this is a naive implementation of differentiation.
#       cf https://github.com/ACEsuit/ACE.jl/issues/27
#       for further discussion
*(φ::Invariant, dAA::SVector) = φ.val * dAA

coco_init(::Invariant, l, m, μ, T, A) = (
      l == m == μ == 0 ? Invariant(T(1)) : Invariant(T(0))  )

coco_zeros(::Invariant, ll, mm, kk, T, A) = Invariant(zero(T))

coco_filter(::Invariant, ll, mm) =
            iseven(sum(ll)) && (sum(mm) == 0)

coco_filter(::Invariant, ll, mm, kk) =
            iseven(sum(ll)) && (sum(mm) == sum(kk) == 0)

coco_dot(u1::Invariant, u2::Invariant) = u1.val * u2.val

# --------------------- EuclideanVector

@doc raw"""
`struct EuclideanVector{D, T}` : specifies that the output $\varphi$ of an
ACE is an equivariant $\varphi \in \mathbb{R}^{3}$, i.e., it transforms under
$O(3)$ as
```math
      \varphi \circ Q = Q \cdot \varphi,
```
where $\cdot$ denotes the standard matrix-vector product.
"""
struct EuclideanVector{T} <: AbstractProperty
   val::SVector{3, T}
end

EuclideanVector{T}() where {T <: Number} = EuclideanVector{T}(zero(SVector{3, T}))

EuclideanVector(T::DataType=Float64) = EuclideanVector{T}()

filter(φ::EuclideanVector, b::Array) = ( length(b) <= 1 ? true :
             isodd( sum(bi.l for bi in b)) &&
            (abs(sum(bi.m for bi in b)) <= 1) )

rot3Dcoeffs(::EuclideanVector,T=Float64) = Rot3DCoeffsEquiv{T,1}(Dict[], ClebschGordan(T))

# differentiation - cf #27
# *(φ::EuclideanVector, dAA::SVector) = φ.val * dAA'

coco_init(phi::EuclideanVector{CT}, l, m, μ, T, A) where {CT} = (
      (l == 1 && abs(m) <= 1 && abs(μ) <= 1)
         ? [EuclideanVector{CT}(rmatrices[m,μ][:,k]) for k=1:3]
         : coco_zeros(phi, l, m, μ, T, A)  )

coco_zeros(φ::EuclideanVector, ll, mm, kk, T, A) = zeros(typeof(φ), 3)

coco_filter(::EuclideanVector, ll, mm) =
            isodd(sum(ll)) && (abs(sum(mm)) <= 1)

coco_filter(::EuclideanVector, ll, mm, kk) =
      abs(sum(mm)) <= 1 &&
      abs(sum(kk)) <= 1 &&
      isodd(sum(ll))

coco_dot(u1::EuclideanVector, u2::EuclideanVector) = dot(u1.val, u2.val)

rmatrices = Dict(
  (-1,-1) => SMatrix{3, 3, ComplexF64, 9}(1/6, 1im/6, 0, -1im/6, 1/6, 0, 0, 0, 0),
  (-1,0) => SMatrix{3, 3, ComplexF64, 9}(0, 0, 0, 0, 0, 0, 1/(3*sqrt(2)), 1im/(3*sqrt(2)), 0),
  (-1,1) => SMatrix{3, 3, ComplexF64, 9}(-1/6, -1im/6, 0, -1im/6, 1/6, 0, 0, 0, 0),
  (0,-1) => SMatrix{3, 3, ComplexF64, 9}(0, 0, 1/(3*sqrt(2)), 0, 0, -1im/(3*sqrt(2)), 0, 0, 0),
  (0,0) => SMatrix{3, 3, ComplexF64, 9}(0, 0, 0, 0, 0, 0, 0, 0, 1/3),
  (0,1) => SMatrix{3, 3, ComplexF64, 9}(0, 0, -1/(3*sqrt(2)), 0, 0, -1im/(3*sqrt(2)), 0, 0, 0),
  (1,-1) => SMatrix{3, 3, ComplexF64, 9}(-1/6, 1im/6, 0, 1im/6, 1/6, 0, 0, 0, 0),
  (1,0) => SMatrix{3, 3, ComplexF64, 9}(0, 0, 0, 0, 0, 0, -1/(3*sqrt(2)), 1im/(3*sqrt(2)), 0),
  (1,1) => SMatrix{3, 3, ComplexF64, 9}(1/6, -1im/6, 0, 1im/6, 1/6, 0, 0, 0, 0)
  )

# --------------------- SphericalVector

struct SphericalVector{L, LEN, T} <: AbstractProperty
   val::SVector{LEN, T}
   _valL::Val{L}
end

# # differentiation - cf #27
# *(φ::SphericalVector, dAA::SVector) = φ.val * dAA'

getL(φ::SphericalVector{L}) where {L} = L

# L = 0 -> (0,0)
# L = 1 -> (0,0), (1,-1), (1,0), (1,1)  -> 4
# L = 3 ->  ... + 5 -> 9
# 1 + 3 + 5 + ... + 2*L+1
# = L + 2 * (1 + ... + L) = L+1 + 2 * L * (L+1) / 2 = (L+1)^2
function SphericalVector(L::Integer; T = Float64)
   LEN = 2L+1   # length of SH basis up to L
   return SphericalVector( zero(SVector{LEN, T}), Val{L}() )
end

function SphericalVector{L, LEN, T}(x::AbstractArray) where {L, LEN, T}
   @assert length(x) == LEN
   SphericalVector{L, LEN, T}( SVector{LEN, T}(x), Val(L) )
end

SphericalVector{L, LEN, T}()  where {L, LEN, T} =
      SphericalVector( zero(SVector{LEN, T}), Val{L}() )

filter(φ::SphericalVector, b::Array) = ( length(b) <= 1 ? true :
        ( ( iseven(sum(bi.l for bi in b)) == iseven(getL(φ)) ) &&
         ( abs(sum(bi.m for bi in b)) <= getL(φ) )  ) )

rot3Dcoeffs(::SphericalVector, T::DataType=Float64) = Rot3DCoeffs(T)

const __rotcoeff_inv = Rotations3D.Rot3DCoeffs(Invariant())

using ACE.Wigner: rotation_D_matrix_ast, rotation_D_matrix

# Equation (1.2) - vector value coupling coefficients
# ∫_{SO3} D^{ll}_{μμmm} D^*(Q) e^t dQ -> 2L+1 column vector
function vec_cou_coe(rotc::Rot3DCoeffs{T},
					      l::Integer, m::Integer, μ::Integer,
					      L::Integer, t::Integer) where {T,N}
	@assert 0 < t <= 2L+1
	D = rotation_D_matrix_ast(L)   # Dt = D[:,t]  -->  # D^* ⋅ e^t
	LL = SA[l, L]
	Z = ntuple(i -> begin
			cc = (rotc(LL, SA[μ, D[i, t].m], SA[m, D[i, t].μ]).val)::T
			D[i, t].sign * cc
		end, 2*L+1)
	return SphericalVector{L, 2L+1, Complex{T}}(SVector(Z))
end

function _select_t(φ::SphericalVector{L}, l, M, K) where {L}
	D = rotation_D_matrix_ast(L)
	tret = -1; numt = 0
	for t = 1:2L+1
		prodμt = prod( (D[i, t].μ + M) for i in 1:2L+1)  # avoid more allocations
		prodmt = prod( (D[i, t].m + K) for i in 1:2L+1)
		if prodμt == prodmt == 0
			tret = t; numt += 1
		end
	end
	# We assumed that there is only one coefficient; this will warn us if it fails
	@assert numt == 1
	return tret
end

coco_zeros(φ::TP, ll, mm, kk, T, A)  where {TP <: SphericalVector} = zero(TP)

coco_dot(u1::SphericalVector, u2::SphericalVector) =
		dot(u1.val, u2.val)

coco_filter(φ::SphericalVector{L}, ll, mm) where {L} =
		iseven(sum(ll) + L) && (abs(sum(mm)) <=  L)

coco_filter(φ::SphericalVector{L}, ll, mm, kk) where {L} =
      iseven(sum(ll) + L) && (abs(sum(mm)) <=  L) && (abs(sum(kk)) <= L)

coco_init(φ::SphericalVector{L}, l, m, μ, T, A) where {L} =
			vec_cou_coe(__rotcoeff_inv, l, m, μ, L, _select_t(φ, l, m, μ))

# --------------- SphericalMatrix

struct SphericalMatrix{L1, L2, LEN1, LEN2, T, LL} <: AbstractProperty
   val::SMatrix{LEN1, LEN2, T, LL}
   _valL1::Val{L1}
   _valL2::Val{L2}
end

# differentiation - cf #27
# actually this here appears to be the generic form how to do the
# differentiation for arbtirary order tensors.
# *(φ::SphericalMatrix{L1, L2, LEN1, LEN2}, dAA::SVector{N}
#       ) where {L1, L2, LEN1, LEN2, N} =
#       reshape(φ.val[:] * dAA', Size(LEN1, LEN2, N))

getL(φ::SphericalMatrix{L1,L2}) where {L1,L2} = L1, L2

# L = 0 -> (0,0)
# L = 1 -> (0,0), (1,-1), (1,0), (1,1)  -> 4
# L = 3 ->  ... + 5 -> 9
# 1 + 3 + 5 + ... + 2*L+1
# = L + 2 * (1 + ... + L) = L+1 + 2 * L * (L+1) / 2 = (L+1)^2
function SphericalMatrix(L1::Integer, L2::Integer; T = Float64)
   LEN1 = 2L1+1   # length of SH basis up to L
   LEN2 = 2L2+1
   return SphericalMatrix( zero(SMatrix{LEN1, LEN2, T}), Val{L1}(), Val{L2}() )
end

# this is a mess - need to fix it, we should need just one constructor???

function SphericalMatrix{L1, L2, LEN1, LEN2, T, LL}(x::AbstractMatrix) where {L1, L2, LEN1, LEN2, T, LL}
   @assert size(x) == (LEN1, LEN2)
   SphericalMatrix{L1, L2, LEN1, LEN2, T, LL}( SMatrix{LEN1, LEN2, T}(x), Val(L1), Val(L2) )
end

SphericalMatrix{L1, L2, LEN1, LEN2, T, LL}()  where {L1, L2, LEN1, LEN2, T, LL} =
      SphericalMatrix( zero(SMatrix{LEN1, LEN2, T}), Val{L1}(), Val{L2}() )

SphericalMatrix{L1, L2, LEN1, LEN2, T}()  where {L1, L2, LEN1, LEN2, T, LL} =
		SphericalMatrix( zero(SMatrix{LEN1, LEN2, T}), Val{L1}(), Val{L2}() )

filter(φ::SphericalMatrix, b::Array) = ( length(b) < 1 ? true :
        ( ( iseven(sum(bi.l for bi in b)) == iseven(sum(getL(φ))) ) &&
         ( abs(sum(bi.m for bi in b)) <= sum(getL(φ)) )  ) )

rot3Dcoeffs(::SphericalMatrix, T::DataType=Float64) = Rot3DCoeffs(T)

function mat_cou_coe(rotc::Rot3DCoeffs{T},
				   		l::Integer, m::Integer, μ::Integer,
					   	a::Integer, b::Integer,
							::Val{L1}, ::Val{L2}) where {T, L1, L2}
	@assert (0 < a <= 2L1 + 1) && (0 < b <= 2L2 + 1)
	Z = zero(MMatrix{2L1+1, 2L2+1, Complex{T}})  # zeros(2 * L1 + 1, 2 * L2 + 1)
	Dp = rotation_D_matrix_ast(L1)
	Dq = rotation_D_matrix(L2)
	LL = SA[l, L1, L2]
	for i = 1:(2 * L1 + 1)
		for j = 1:(2 * L2 + 1)
			MM = SA[μ, Dp[i,a].m, Dq[b,j].m]
			KK = SA[m, Dp[i,a].μ, Dq[b,j].μ]
			cc = (rotc(LL, MM, KK).val)::T
			Z[i,j] = Dp[i,a].sign * Dq[b,j].sign * cc
		end
	end
	return SphericalMatrix(SMatrix(Z), Val{L1}(), Val{L2}())
	# return SphericalMatrix(SMatrix{2L1+1,2L2+1,Complex{T}}(Z), Val{L1}(), Val{L2}())
end

function _select_ab(φ::SphericalMatrix{L1,L2}, M, K) where {L1,L2}
   Dp = rotation_D_matrix_ast(L1)
   Dq = rotation_D_matrix(L2)
   list_ab = Tuple{Int, Int}[]
   for a = 1:2L1+1
      for b = 1:2L2+1
			# pm = prod( ma[i] + mb[j] + K for i = 1:2L1+1, j = 1:2L2+1)
			pm = prod( Dp[i,a].m + Dq[b,j].m + K for i = 1:2L1+1, j = 1:2L2+1)
			# pμ = prod(μa[i] + μb[j] + M for i = 1:2L1+1, j = 1:2L2+1)
			pμ = prod( Dp[i,a].μ + Dq[b,j].μ + M for i = 1:2L1+1, j = 1:2L2+1)
         if pμ == pm ==0
				push!(list_ab, (a,b))
         end
      end
   end
	return list_ab
end

function coco_init(φ::SphericalMatrix{L1,L2}, l, m, μ, T, A) where{L1,L2}
   list = _select_ab(φ, m, μ)
	@assert length(list) > 0
	# MAIN CASE
   if iseven(l + L1 + L2) && abs(m) <= L1+L2 && abs(μ) <= L1+L2
		return [ mat_cou_coe(__rotcoeff_inv, l, m, μ, a, b, Val(L1), Val(L2))
				   for (a,b) in list ]
	end

	# @warn("SHOULDN'T BE HERE!!")
   return fill( zero(typeof(φ)), length(list) )
end

coco_zeros(φ::TP, ll, mm, kk, T, A) where{TP <: SphericalMatrix} =
            zeros(TP, length(_select_ab(φ, sum(mm), sum(kk))))

coco_filter(φ::SphericalMatrix{L1,L2}, ll, mm) where {L1,L2} =
            iseven(sum(ll)) == iseven(L1+L2) && (abs(sum(mm)) <=  L1+L2)

coco_filter(φ::SphericalMatrix{L1,L2}, ll, mm, kk) where {L1,L2} =
            iseven(sum(ll) + L1 + L2) &&
            (abs(sum(mm)) <=  L1+L2) &&
            (abs(sum(kk)) <=  L1+L2)

coco_dot(u1::SphericalMatrix, u2::SphericalMatrix) =
		dot(u1.val, u2.val)


text:
function tryretry(f, retries)
    for i in 1:retries - 1
        try
           return f()
        catch
            @info "Left $(retries - i) trials"
        end
        sleep(rand(1:10))
    end

    f()
end

struct SOError{T} <: Exception
    msg::T
end

struct SOClient
    ep::String
end

function SOClient(; ep = "https://api.stackexchange.com", api_version = "2.2", use_globally = true)
    ep = ep * "/" * api_version * "/"

    client = SOClient(ep)
    if use_globally
        SOGlobal[] = client
    end

    return client
end

const SOGlobal = Ref(SOClient(; use_globally = false))

########################################
# JSON structures
########################################
mutable struct TagQuestionId
    tags::Vector{String}
    question_id::Int
end
TagQuestionId() = TagQuestionId(String[], -1)
StructTypes.StructType(::Type{TagQuestionId}) = StructTypes.Mutable()

mutable struct SOResponse{T}
    items::Vector{T}
    has_more::Bool
    quota_max::Int
    quota_remaining::Int
end
SOResponse{T}() where T = SOResponse(T[], false, 0, 0)
StructTypes.StructType(::Type{<:SOResponse}) = StructTypes.Mutable()

mutable struct Owner
    display_name::String
    link::String
end
Owner() = Owner("Unknown", "")
StructTypes.StructType(::Type{Owner}) = StructTypes.Mutable()

mutable struct Question
    tags::Vector{String}
    question_id::Int
    body::String
    owner::Owner
    title::String
    link::String
    creation_date::Int
    last_activity_date::Int
    score::Int
    answer_count::Int
    is_answered::Bool
end
Question() = Question([], 0, "", Owner(), "", "", 0, 0, 0, 0, false)
StructTypes.StructType(::Type{Question}) = StructTypes.Mutable()

mutable struct Answer
    is_accepted::Bool
    body::String
    owner::Owner
    question_id::Int
    answer_id::Int
    creation_date::Int
    last_activity_date::Int
    score::Int
end
Answer() = Answer(false, "", Owner(), 0, 0, 0, 0, 0)
StructTypes.StructType(::Type{Answer}) = StructTypes.Mutable()

########################################
# SO processing utils
########################################
function process_response(response, T = Nothing)
    headers = Dict(response.headers)
    if get(headers, "content-encoding", "") == "gzip"
        stream = response.body |> IOBuffer |> GzipDecompressorStream
        if T === Nothing
            return JSON3.read(stream)
        else
            return JSON3.read(stream, T)
        end
    else
        # Do not know, what to do, let's throw an error
        @error "In stackoverflow response unknown content-encoding " * get(headers, "content-encoding", "\"empty field\"")
        throw(SOError("Not gzip encoded header: " * JSON3.write(headers)))
    end
end

function searchtag(client::SOClient = SOGlobal[]; params...)
    params = Dict(params) |> HTTP.URIs.escapeuri
    url = client.ep * "search?"
    response = HTTP.get(url * params)

    return process_response(response, SOResponse{TagQuestionId})
end

getquestions(qids, retries = 5; params...) = getquestions(SOGlobal[], qids, retries; params...)
function getquestions(client::SOClient, qids, retries; params...)
    url = client.ep * "questions/"
    qids = join(qids, ";") |> HTTP.URIs.escapeuri
    params = Dict(params) |> HTTP.URIs.escapeuri

    uri = url * qids * "?" * params
    println("============================================")
    println(uri)
    println("============================================")
    tryretry(retries) do
        response = HTTP.get(uri)
        process_response(response, SOResponse{Question})
    end
end

function getqanswers(client::SOClient, qids, retries; params...)
    url = client.ep * "questions/"
    qids = join(qids, ";") |> HTTP.URIs.escapeuri
    params = Dict(params) |> HTTP.URIs.escapeuri

    uri = url * qids * "/answers?" * params
    println("============================================")
    println(uri)
    println("============================================")

    tryretry(retries) do
        response = HTTP.get(uri)
        process_response(response, SOResponse{Answer})
    end
end

getallqanswers(qids, retries = 5; params...) = getallqanswers(SOGlobal[], qids, retries; params...)
function getallqanswers(client::SOClient, qids, retries; params...)
    answers = []
    hasmore = true
    page = 1
    while hasmore
        res = getqanswers(client, qids, retries; page = page, params...)
        push!(answers, res)
        page += 1
        hasmore = res.has_more
    end

    return answers
end


text:
"""
    Node

Abstract type for nodes in the AST tree.
"""
abstract type Node end

max_stars_repo_name: Sohva/oolong

max_issues_repo_name: Sohva/oolong

max_forks_repo_name: oolongtensor/oolong


text:
# function to aggregate the ressource flows of the simulation into
# average weekly data over all locations
function ressource_status(ressource_flow::Array{Int64,3},
                          incidents::DataFrame,
                          simulation_capacity::Array{Int64,2})
    # Dataframe to be filled
    capacity_status = DataFrame(sum(ressource_flow, dims = 2)[:,1,:],
                                        [:at_location,
                                        :to_incident,
                                        :at_incident,
                                        :to_location,
                                        :at_backlog,
                                        :at_patrol,
                                        :backlog_minutes])
    # append the weekhour to the ressource flow DataFrame
    current_weekhour = incidents[1,:weekhour]
    current_minute   = minute(unix2datetime(incidents[1,:epoch]))
    capacity_status[:,:weekhour] .= 0
    for i = 1:size(capacity_status,1)
        if current_minute == 60
            current_minute = 0
            if current_weekhour >= size(simulation_capacity,1)
                current_weekhour = 0
            end
            current_weekhour += 1
        end
        capacity_status[i,:weekhour]  = current_weekhour
        current_minute +=1
    end
    # group the DataFrame by weekhour
    capacity_status = groupby(capacity_status, :weekhour)
    capacity_status = combine(capacity_status,
                        :at_location => mean => :at_location,
                        :to_incident => mean => :to_incident,
                        :at_incident => mean => :at_incident,
                        :to_location => mean => :to_location,
                        :at_backlog  => mean => :at_backlog,
                        :at_patrol   => mean => :at_patrol,
                        :backlog_minutes => mean => :backlog_minutes)
    capacity_status = sort!(capacity_status, :weekhour)

    # create an additional DataFrame to hold the average backlog for
    # each weekhour over all locations
    capacity_status = select!(capacity_status, Not([:backlog_minutes]))
    main_capacity_status = combine(capacity_status, [n => mean => n for n in names(capacity_status)])
    capacity_status = sort!(stack(capacity_status, 2:7), [:weekhour])

    return capacity_status, main_capacity_status
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta
import ._ProtoBuf_Top_.opentelemetry

mutable struct ExportTraceServiceRequest <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ExportTraceServiceRequest(; kwargs...)
        obj = new(meta(ExportTraceServiceRequest), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct ExportTraceServiceRequest
const __meta_ExportTraceServiceRequest = Ref{ProtoMeta}()
function meta(::Type{ExportTraceServiceRequest})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ExportTraceServiceRequest)
            __meta_ExportTraceServiceRequest[] =
                target = ProtoMeta(ExportTraceServiceRequest)
            allflds = Pair{Symbol,Union{Type,String}}[:resource_spans=>Base.Vector{
                opentelemetry.proto.trace.v1.ResourceSpans,
            }]
            meta(
                target,
                ExportTraceServiceRequest,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_ExportTraceServiceRequest[]
    end
end
function Base.getproperty(obj::ExportTraceServiceRequest, name::Symbol)
    if name === :resource_spans
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{
            opentelemetry.proto.trace.v1.ResourceSpans,
        }
    else
        getfield(obj, name)
    end
end

mutable struct ExportTraceServiceResponse <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ExportTraceServiceResponse(; kwargs...)
        obj = new(meta(ExportTraceServiceResponse), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct ExportTraceServiceResponse
const __meta_ExportTraceServiceResponse = Ref{ProtoMeta}()
function meta(::Type{ExportTraceServiceResponse})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ExportTraceServiceResponse)
            __meta_ExportTraceServiceResponse[] =
                target = ProtoMeta(ExportTraceServiceResponse)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(
                target,
                ExportTraceServiceResponse,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_ExportTraceServiceResponse[]
    end
end

# service methods for TraceService
const _TraceService_methods = MethodDescriptor[MethodDescriptor(
    "Export",
    1,
    ExportTraceServiceRequest,
    ExportTraceServiceResponse,
)] # const _TraceService_methods
const _TraceService_desc = ServiceDescriptor(
    "opentelemetry.proto.collector.trace.v1.TraceService",
    1,
    _TraceService_methods,
)

TraceService(impl::Module) = ProtoService(_TraceService_desc, impl)

mutable struct TraceServiceStub <: AbstractProtoServiceStub{false}
    impl::ProtoServiceStub
    TraceServiceStub(channel::ProtoRpcChannel) =
        new(ProtoServiceStub(_TraceService_desc, channel))
end # mutable struct TraceServiceStub

mutable struct TraceServiceBlockingStub <: AbstractProtoServiceStub{true}
    impl::ProtoServiceBlockingStub
    TraceServiceBlockingStub(channel::ProtoRpcChannel) =
        new(ProtoServiceBlockingStub(_TraceService_desc, channel))
end # mutable struct TraceServiceBlockingStub

Export(
    stub::TraceServiceStub,
    controller::ProtoRpcController,
    inp::ExportTraceServiceRequest,
    done::Function,
) = call_method(stub.impl, _TraceService_methods[1], controller, inp, done)
Export(
    stub::TraceServiceBlockingStub,
    controller::ProtoRpcController,
    inp::ExportTraceServiceRequest,
) = call_method(stub.impl, _TraceService_methods[1], controller, inp)

export ExportTraceServiceRequest,
    ExportTraceServiceResponse,
    TraceService,
    TraceServiceStub,
    TraceServiceBlockingStub,
    Export


text:
const Assignment = Dict{NodeName, Any}

"""
    nodenames(a::Assignment)
Return a vector of NodeNames (aka symbols) for the assignment
"""
nodenames(a::Assignment) = collect(keys(a))

"""
    consistent(a::Assignment, b::Assignment)
True if all shared NodeNames have the same value
"""
function consistent(a::Assignment, b::Assignment)

    for key in keys(a)
        if haskey(b,key) && b[key] != a[key]
            return false
        end
    end

    return true
end