text:
# load all cases (necessary)
## the .csv file has to contain 7 columns representing the incidents
### column 1: caseid -> unique id for each incident (Int)
### column 2: priority -> priority of the incident starting with 1 as
###           as the highest priority (Int)
### column 3: cars -> number of cars requested (Int)
### column 4: length -> length of the incident in minutes (Int)
### column 5: backlog -> length of the paperwork/ backlog after the
###           incident in minutes (Int)
### column 6: location -> unique id of the corresponding basic area
###           where the incident took place (Int)
### column 7: epoch -> epoch time of the beginning of the incident
###           in seconds (Int)
    incidents = CSV.read("data/$problem/incidents_$problem.csv", DataFrame)

# load the aerial distances between the basic areas (necessary)
## 2-dimensional array containing the airial distance between all
## combinations of basic areas (scale irrelevant, Float64)
## example: row 3 column 5 corresponds to the airial distance
##          between basic area 3 and basic area 5
    airdist = readdlm("data/$problem/airdistances_$problem.csv", Float64)

# load the driving times between the basic areas(necessary)
## 2-dimensional array containing the driving time between all
## combinations of basic areas in minutes (Float64)
## example: row 3 column 5 corresponds to the driving time
##          between basic area 3 and basic area 5 in minutes
    drivingtime = readdlm("data/$problem/drivingtimes_$problem.csv", Float64)

# load the shift pattern for each weekhour (necessary)
## array with 168 entrys, each entry corresponds to one weekhour
## each shift should be assigned with a number starting from 1
    shifts = readdlm("data/shifts_8hours.csv", Int64)

# load the real capacity plan of the emergency service
## array with 168 rows, where each entry corresponds to one weekhour
## the number of rows depends on the number of districts (number_districts)
    if real_capacity == true
        simulation_capacity = readdlm("data/$problem/capacity_$problem.csv", Int64)
        if size(simulation_capacity,2) != number_districts
            error("Number of districts does not match with the capacity plan!")
        end
        total_capacity = sum(simulation_capacity)/168
    end

# load the adjacency matrix between the basic areas(if available)
## 2-dimensional array containing the adjacency between all
## combinations of basic areas (Int64). Important, if rivers or
## other barriers are contained within the area
## example: row 3 column 5 corresponds to the adjacency
##          between basic area 3 and basic area 5
## if no data is available, a matrix is created based on the
## airial distance between the hexagonal BAs
    if isfile("data/$problem/adjacency_$problem.csv")
        adjacent = readdlm("data/$problem/adjacency_$problem.csv", Bool)
    else
        adjacent = adjacency_matrix(airdist::Array{Float64,2})
    end

# load the array of potential BAs for a district center (if available)
# if no array is available, all BAs are potential locations (Int)
## binary array with a 1 if a BA represents a potential center, else 0
## the first entry in the array corresponds to the hexagon with the id 1
    if isfile("data/$problem/potential_locations_$problem.csv")
        potential_locations = readdlm("data/$problem/potential_locations_$problem.csv", Bool)
    else
        potential_locations = Vector{Bool}(undef,size(airdist,1)) .= 1
    end

# load the array of current departments (if available)
# if no array is available, no BA is assumed as current location (Int)
## binary array with a 1 if a BA represents a current location, else 0
## the first entry in the array corresponds to the hexagon with the id 1
    if isfile("data/$problem/current_departments_$problem.csv")
        current_locations = Vector{Bool}(undef,size(airdist,1)) .= 0
        cl = readdlm("data/$problem/current_departments_$problem.csv", Int64)[:,1]
        for i = 1:length(cl)
            current_locations[cl[i]] = true
        end
    else
        current_locations = Vector{Bool}(undef,size(airdist,1)) .= 0
    end

# load the traffic flow parameter for each weekhour (if available)
# if no array is available, no markup for the traffic is expected
## the .csv file has to contain 2-columns with 168 entrys (Float64)
### column 1: traffic flow parameter q
### column 2: standard deviation of q
### the first entry corresponds to the data for the first weekhour
### the first weekhour starts each monday at 00:00
    if isfile("data/$problem/traffic_flow_$problem.csv")
        traffic = readdlm("data/$problem/traffic_flow_$problem.csv", Float64)
    else
        traffic = Array{Float64,2}(undef,168,2)
        traffic[:,1] .= 1
        traffic[:,2] .= 0
    end

# load the shape file of the problem instance (if available)
## the .shp file can contain multiple columns, 2 columns are important
## for our model:
### column :id hexagon id
### column :geometry shape file of the hexagon
    if isfile("grids/$problem/grid_$problem.shp")
        hexshape = Shapefile.Table("grids/$problem/grid_$problem.shp")
        hexshape = DataFrame(hexshape)
        hexshape = sort!(hexshape, :id)
    else
        hexshape = nothing
    end

# Create a theme for the Plots
    theme(:vibrant,
    titlefontsize = 8,
    legendfontsize = 3,
    legendtitlefontsize = 6,
    size = (800,400), show = true, dpi = 300,
    legend = :outerright)


text:
module HERE
include("XMLiterate.jl")

include("TrafficAPI.jl")
export flow

include("WeatherAPI.jl")
export weather

end

ext: jl


text:
using Test
using MDCDL
using FFTW
using Random

@testset "CNSOLT" begin
    include("testsetGenerator.jl")

    ccsd1D = cnsoltValidConfigSet1D()
    ccsd2D = cnsoltValidConfigSet2D()
    ccsd3D = cnsoltValidConfigSet3D()

    ccsd = [ ccsd1D, ccsd2D, ccsd3D ]

    Random.seed!(392875929)

    @testset "Constructor" begin

        maxDims = 3
        defaultType = Float64

        # cfgset = [ (crdf.I, nch, crord.I .- 1) for crdf in CartesianIndices(tuple(fill(4,d)...)), nch in 2:20, crord in CartesianIndices(tuple(fill(6,d)...)) ]

        for d in 1:maxDims
            allcfgset = [ (crdf.I, crord.I .- 1, nch) for crdf in CartesianIndices(tuple(fill(4,d)...)), crord in CartesianIndices(tuple(fill(6+1,d)...)), nch in 2:20 ]
            cfgset = randsubseq(vec(allcfgset), 100 / length(allcfgset))

            for (df, ord, nch) in cfgset
                if prod(df) > sum(nch)
                    @test_throws AssertionError Cnsolt(df, ord, nch)
                    continue
                end

                if isodd(sum(nch)) && any(isodd.(ord))
                    @test_throws AssertionError Cnsolt(df, ord, nch)
                    continue
                end

                nsolt = Cnsolt(df, ord, nch)

                if iseven(sum(nch))
                    @test nsolt isa CnsoltTypeI
                    @test istype1(nsolt) == true && istype2(nsolt) == false
                else
                    @test nsolt isa CnsoltTypeII
                    @test istype1(nsolt) == false && istype2(nsolt) == true
                end
            end
        end
    end

    @testset "DefaultValues" begin
        @testset "maximally-sampled" begin
            for d in 1:3
                for n = 1:20
                    df = (rand(1:4, d)...,)
                    ord = (fill(0, d)...,)
                    nch = prod(df)

                    x = rand(df...)
                    A = MDCDL.cdftmtx(df...)
                    expctd = A*vec(x)

                    nsolt = Cnsolt(df, ord, nch)
                    pvx = reshape(x, prod(df), :)
                    y = analyze(nsolt, pvx, (fill(1,d)...,))
                    actual = reshape(y, prod(df))

                    @test actual ≈ expctd
                end
            end
        end

        @testset "over-sampled" begin
            for d in 1:3
                for n = 1:20
                    df = (rand(1:4, d)...,)
                    ord = (fill(0, d)...,)
                    nch = prod(df) + rand(0:6)

                    x = rand(df...)
                    A = MDCDL.cdftmtx(df...)
                    expctd = A*vec(x)

                    nsolt = Cnsolt(df, ord, nch)
                    pvx = reshape(x, prod(df), :)
                    y = analyze(nsolt, pvx, (fill(1,d)...,))
                    actual = reshape(y[1:prod(df),:], prod(df))

                    @test actual ≈ expctd
                end
            end
        end
    end

    @testset "FilterSymmetry" begin
        for d in 1:length(ccsd), (df, ord, nch) in ccsd[d]
            nsolt = Cnsolt(df, ord, nch)
            rand!(nsolt)

            afb = analysisbank(nsolt)
            hsafb = nsolt.Φ' * afb

            @test hsafb ≈ conj(reverse(hsafb; dims=2))
        end
    end

    @testset "AnalysisSynthesis" begin
        # output mode options for analyzer
        oms = [ Shapes.Separated, Shapes.Combined, Shapes.Vec ]
        @testset "NormalOrder" begin
            for d in 1:length(ccsd), (df, ord, nch) in ccsd[d]
                szx = df .* (ord .+ 1)
                nsolt = Cnsolt(df, ord, nch)
                rand!(nsolt)

                x = rand(Complex{Float64}, szx...)

                foreach(oms) do om
                    nsop = createTransform(nsolt, om(size(x)))

                    y = analyze(nsop, x)
                    rx = synthesize(nsop, y)

                    @test size(x) == size(rx)
                    @test rx ≈ x
                end
            end
        end
        @testset "DimensionPermutation" begin
            for d in 1:length(ccsd), (df, ord, nch) in ccsd[d]
                szx = df .* (ord .+ 1)
                nsolt = Cnsolt(df, ord, nch, perm=(randperm(d)...,))
                rand!(nsolt)

                x = rand(Complex{Float64}, szx...)

                foreach(oms) do om
                    nsop = createTransform(nsolt, om(size(x)))

                    y = analyze(nsop, x)
                    rx = synthesize(nsop, y)

                    @test size(x) == size(rx)
                    @test rx ≈ x
                end
            end
        end
    end

    @testset "AnalyzerKernel" begin
        for d in 1:length(ccsd), (df, ord, nch) in ccsd[d]
            szx = df .* (ord .+ 1)
            x = rand(Complex{Float64}, szx)

            nsolt = Cnsolt(df, ord, nch)
            rand!(nsolt)

            ana = createTransform(nsolt, Shapes.Separated())
            ya = analyze(ana, x)

            afs = analysiskernels(nsolt)
            myfilter(A, h) = begin
                ha = zero(A)
                ha[UnitRange.(1, size(h))...] = h
                ifft(fft(A).*fft(ha))
            end
            offset = df .- 1

            sys = map(afs) do af
                fx = myfilter(x, af)
                dwfx = downsample(fx, df, offset)
                circshift(dwfx, (-1 .* fld.(ord,2)))
            end
            @test all(ya .≈ sys)
        end
    end

    @testset "SynthesizerKernel" begin
        for d in 1:length(ccsd), (df, ord, nch) in ccsd[d]

            nsolt = Cnsolt(df, ord, nch)
            rand!(nsolt)

            y = [ rand(Complex{Float64},((ord.+1) .* df)...) for p in 1:sum(nch) ]

            syn = createTransform(nsolt, Shapes.Separated())
            x = synthesize(syn, y)

            sfs = synthesiskernels(nsolt)
            myfilter(A, h) = begin
                ha = zero(A)
                ha[UnitRange.(1, size(h))...] = h
                ifft(fft(A).*fft(ha))
            end
            offset = df .- 1

            subrxs = sum(map((yp,fp)->myfilter( MDCDL.upsample(yp, df), fp ), y, sfs))
            subrx = circshift(subrxs, -1 .* df .* cld.(ord,2))
            @test size(x) == size(subrx)
            @test subrx ≈ x
        end
    end

    @testset "Factorization" begin
        for d in 1:length(ccsd), (df, ord, nch) in ccsd[d]
            src = Cnsolt(df, ord, nch)
            dst = Cnsolt(df, ord, nch)
            rand!(src)

            (angs, mus) = getrotations(src)
            setrotations!(dst, angs, mus)

            @test src.V0 ≈ dst.V0
            foreach(src.Wdks, dst.Wdks) do srcWs, dstWs
                @test all(srcWs .≈ dstWs)
            end
            foreach(src.Udks, dst.Udks) do srcUs, dstUs
                @test all(srcUs .≈ dstUs)
            end

            if istype1(src)
                foreach(src.θdks, dst.θdks) do srcθs, dstθs
                    @test all(srcθs .≈ dstθs)
                end
            else
                foreach(src.θ1dks, dst.θ1dks) do srcθ1s, dstθ1s
                    @test all(srcθ1s .≈ dstθ1s)
                end
                foreach(src.Ŵdks, dst.Ŵdks) do srcŴs, dstŴs
                    @test all(srcŴs .≈ dstŴs)
                end
                foreach(src.Ûdks, dst.Ûdks) do srcÛs, dstÛs
                    @test all(srcÛs .≈ dstÛs)
                end
                foreach(src.θ2dks, dst.θ2dks) do srcθ2s, dstθ2s
                    @test all(srcθ2s .≈ dstθ2s)
                end
            end
        end
    end
end

nothing

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
using Pcap
using Base.Test

# write your own tests here
#@test 1 == 1

#println(pcap_lookupdev())
#pcap_findalldevs()

function test_eth_hdr(ethhdr::EthHdr)
    @test "74:de:2b:08:78:09" == ethhdr.src_mac
    @test "00:24:fe:b1:8f:dc" == ethhdr.dest_mac
    @test 2048 == ethhdr.ptype
end # function test_ip_hdr

function test_ip_hdr(iphdr::IpHdr)
    @test 4 == iphdr.version
    @test 20 == iphdr.length
    @test 0 == iphdr.services
    @test 63 == iphdr.totlen
    @test 20831 == iphdr.id
    @test 64    == iphdr.ttl
    @test 17    == iphdr.protocol
    @test "192.168.0.51" == iphdr.src_ip
    @test "192.168.0.1"  == iphdr.dest_ip
    @test true           == iphdr.checksum
end # function test_ip_hdr

function test_udp_hdr(udphdr::UdpHdr)
    @test 34904  == udphdr.src_port
    @test 53     == udphdr.dest_port
    @test 43     == udphdr.length
    @test 0xa24a == udphdr.checksum
end # function test_udp_hdr

# test little endian captures and header dissectors
cap     = PcapOffline("data/dns-query-response.pcap")
rec     = pcap_get_record(cap)
layers  = decode_pkt(rec.payload)

test_eth_hdr(layers.datalink)
test_ip_hdr(layers.network)
test_udp_hdr(layers.protocol)

# test big endian capture
cap = PcapOffline("data/big-endian-icmp-echo-request.pcap")
rec = pcap_get_record(cap)
@test rec.ts_sec == 0x569558b0
@test rec.ts_usec == 0x000d4d10
@test rec.incl_len == 0x00000062
@test rec.orig_len == 0x00000062


text:
"""
    Feather.write(file::String, tbl; description="", metadata="")
    Feather.write(io::IO, tbl; description="", metadata="")

Write any Tables.jl implementation as a feather-formatted file. Optionally, a `description`
and `metadata` can be provided as Strings.
"""
function write end

write(filename::AbstractString; kwargs...) = x->write(filename, x; kwargs...)
function write(io::IO, tbl; description::String="", metadata::String="")
    writepadded(io, FEATHER_MAGIC_BYTES)
    metacols = Metadata.Column[]
    columns = Tables.columns(tbl)
    names = propertynames(columns)
    len = 0
    for nm in names
        column = getarrow(getproperty(columns, nm))
        len = length(column)
        vals = writecontents(Metadata.PrimitiveArray, io, column)
        push!(metacols, Metadata.Column(String(nm), vals,
                                        getmetadata(io, eltype(column), column), ""))
    end
    ctable = Metadata.CTable(description, len, metacols, FEATHER_VERSION, metadata)
    metalen = writemetadata(io, ctable)
    Base.write(io, metalen)
    Base.write(io, FEATHER_MAGIC_BYTES)
    io
end
function write(filename::AbstractString, tbl; description::String="", metadata::String="")
    open(filename, "w+") do io
        write(io, tbl, description=description, metadata=metadata)
    end
    filename
end

getarrow(col::AbstractVector{T}) where {T} = arrowformat(col)
function getarrow(col::AbstractVector{Union{T, Missing}}) where {T}
    hasmissing = any(ismissing, col)
    return arrowformat(hasmissing ? col : convert(AbstractVector{T}, col))
end
function getarrow(col::AbstractVector{Missing})
    throw(ArgumentError("Feather format does not support writing `AbstractVector{Missing}`. "*
                        "Consider converting column to `AbstractVector{Union{T,Missing}}` where "*
                        "`T` is a supported type."))
end

function Metadata.PrimitiveArray(A::ArrowVector{J}, off::Integer, nbytes::Integer) where J
    Metadata.PrimitiveArray(feathertype(J), Metadata.PLAIN, off, length(A), nullcount(A), nbytes)
end
function Metadata.PrimitiveArray(A::DictEncoding, off::Integer, nbytes::Integer)
    Metadata.PrimitiveArray(feathertype(eltype(references(A))), Metadata.PLAIN, off, length(A),
                            nullcount(A), nbytes)
end

writecontents(io::IO, A::Primitive) = writepadded(io, A)
writecontents(io::IO, A::NullablePrimitive) = writepadded(io, A, bitmask, values)
writecontents(io::IO, A::List) = writepadded(io, A, offsets, values)
writecontents(io::IO, A::NullableList) = writepadded(io, A, bitmask, offsets, values)
writecontents(io::IO, A::BitPrimitive) = writepadded(io, A, values)
writecontents(io::IO, A::NullableBitPrimitive) = writepadded(io, A, bitmask, values)
writecontents(io::IO, A::DictEncoding) = writecontents(io, references(A))
function writecontents(::Type{Metadata.PrimitiveArray}, io::IO, A::ArrowVector)
    a = position(io)
    writecontents(io, A)
    b = position(io)
    Metadata.PrimitiveArray(A, a, b-a)
end

function writemetadata(io::IO, ctable::Metadata.CTable)
    meta = FlatBuffers.build!(ctable)
    rng = (meta.head+1):length(meta.bytes)
    writepadded(io, view(meta.bytes, rng))
    Int32(length(rng))
end


text:
module OBC

using Distributions
using Iterators
using SAMC

export
    plot_all,
    plot_trace,
    logsum,
    get_bbox,
    gen_grid,
    gelman_rubin,
    get_bbox,
    gen_grid,
    gen_unit_grid

type BinaryClassifier{T<:Sampler,S<:MCMC}
    cls1 :: T
    cls2 :: T
    mcmc1 :: S
    mcmc2 :: S
end

include("cls_utils.jl")
include("utils.jl")

end # module


text:
SUITE["dot"] = BenchmarkGroup()
SUITE["dcontract"] = BenchmarkGroup()
SUITE["otimes"] = BenchmarkGroup()
SUITE["other"] = BenchmarkGroup()
SUITE["promotion"] = BenchmarkGroup()
SUITE["constructors"] = BenchmarkGroup()
SUITE["basic-operations"] = BenchmarkGroup()

for dim in (ALL_DIMENSIONS ? (1,2,3) : (3,))
    for T in (Float64, Float32, dT)
        v1 = tensor_dict[(dim, 1, T)]
        V2 = tensor_dict[(dim, 2, T)]
        V4 = tensor_dict[(dim, 4, T)]
        V2sym = symtensor_dict[(dim, 2, T)]
        V4sym = symtensor_dict[(dim, 4, T)]

        # dot
        SUITE["dot"]["Vec{$dim, $T} ⋅ Vec{$dim, $T}"]                               = @benchmarkable dot($v1, $v1)
        SUITE["dot"]["Tensor{2, $dim, $T} ⋅ Vec{$dim, $T}"]                         = @benchmarkable dot($V2, $v1)
        SUITE["dot"]["Vec{$dim, $T} ⋅ Tensor{2, $dim, $T}"]                         = @benchmarkable dot($v1, $V2)
        SUITE["dot"]["SymmetricTensor{2, $dim, $T} ⋅ Vec{$dim, $T}"]                = @benchmarkable dot($V2sym, $v1)
        SUITE["dot"]["Vec{$dim, $T} ⋅ SymmetricTensor{2, $dim, $T}"]                = @benchmarkable dot($v1, $V2sym)
        SUITE["dot"]["Tensor{2, $dim, $T} ⋅ Tensor{2, $dim, $T}"]                   = @benchmarkable dot($V2, $V2)
        SUITE["dot"]["SymmetricTensor{2, $dim, $T} ⋅ SymmetricTensor{2, $dim, $T}"] = @benchmarkable dot($V2sym, $V2sym)
        if MIXED_SYM_NONSYM
            SUITE["dot"]["SymmetricTensor{2, $dim, $T} ⋅ Tensor{2, $dim, $T}"] = @benchmarkable dot($V2sym, $V2)
            SUITE["dot"]["Tensor{2, $dim, $T} ⋅ SymmetricTensor{2, $dim, $T}"] = @benchmarkable dot($V2, $V2sym)
        end

        # dcontract
        SUITE["dcontract"]["Tensor{2, $dim, $T} ⊡ Tensor{2, $dim, $T}"]                   = @benchmarkable dcontract($V2, $V2)
        SUITE["dcontract"]["SymmetricTensor{2, $dim, $T} ⊡ SymmetricTensor{2, $dim, $T}"] = @benchmarkable dcontract($V2sym, $V2sym)
        if MIXED_SYM_NONSYM
            SUITE["dcontract"]["Tensor{2, $dim, $T} ⊡ SymmetricTensor{2, $dim, $T}"] = @benchmarkable dcontract($V2, $V2sym)
            SUITE["dcontract"]["SymmetricTensor{2, $dim, $T} ⊡ Tensor{2, $dim, $T}"] = @benchmarkable dcontract($V2sym, $V2)
        end
        SUITE["dcontract"]["Tensor{4, $dim, $T} ⊡ Tensor{2, $dim, $T}"] = @benchmarkable dcontract($V4, $V2)
        SUITE["dcontract"]["Tensor{2, $dim, $T} ⊡ Tensor{4, $dim, $T}"] = @benchmarkable dcontract($V2, $V4)
        SUITE["dcontract"]["Tensor{4, $dim, $T} ⊡ Tensor{4, $dim, $T}"] = @benchmarkable dcontract($V4, $V4)
        if MIXED_SYM_NONSYM
            SUITE["dcontract"]["SymmetricTensor{4, $dim, $T} ⊡ Tensor{2, $dim, $T}"] = @benchmarkable dcontract($V4sym, $V2)
            SUITE["dcontract"]["Tensor{2, $dim, $T} ⊡ SymmetricTensor{4, $dim, $T}"] = @benchmarkable dcontract($V2, $V4sym)
            SUITE["dcontract"]["Tensor{4, $dim, $T} ⊡ SymmetricTensor{2, $dim, $T}"] = @benchmarkable dcontract($V4, $V2sym)
            SUITE["dcontract"]["SymmetricTensor{2, $dim, $T} ⊡ Tensor{4, $dim, $T}"] = @benchmarkable dcontract($V2sym, $V4)
            SUITE["dcontract"]["SymmetricTensor{4, $dim, $T} ⊡ Tensor{4, $dim, $T}"] = @benchmarkable dcontract($V4sym, $V4)
            SUITE["dcontract"]["Tensor{4, $dim, $T} ⊡ SymmetricTensor{4, $dim, $T}"] = @benchmarkable dcontract($V4, $V4sym)
        end
        SUITE["dcontract"]["Tensor{4, $dim, $T} ⊡ Tensor{4, $dim, $T}"]                   = @benchmarkable dcontract($V4, $V4)
        SUITE["dcontract"]["SymmetricTensor{4, $dim, $T} ⊡ SymmetricTensor{4, $dim, $T}"] = @benchmarkable dcontract($V4sym, $V4sym)

        # otimes
        SUITE["otimes"]["Vec{$dim, $T} ⊗ Vec{$dim, $T}"]                               = @benchmarkable otimes($v1, $v1)
        SUITE["otimes"]["Tensor{2, $dim, $T} ⊗ Tensor{2, $dim, $T}"]                   = @benchmarkable otimes($V2, $V2)
        SUITE["otimes"]["SymmetricTensor{2, $dim, $T} ⊗ SymmetricTensor{2, $dim, $T}"] = @benchmarkable otimes($V2sym, $V2sym)
        if MIXED_SYM_NONSYM
            SUITE["otimes"]["Tensor{2, $dim, $T} ⊗ SymmetricTensor{2, $dim, $T}"] = @benchmarkable otimes($V2, $V2sym)
            SUITE["otimes"]["SymmetricTensor{2, $dim, $T} ⊗ Tensor{2, $dim, $T}"] = @benchmarkable otimes($V2sym, $V2)
        end

        # other
        for (i, V2t) in enumerate((V2, V2sym))
            TensorType = i == 2 ? "SymmetricTensor" : "Tensor"
            for f in (norm, tr, vol, det, inv, transpose, symmetric, skew, eig, mean, dev)
                (i == 1 || typeof(V2t) <: Tensor || T == dT) && f == eig && continue
                SUITE["other"]["$f($TensorType{2, $dim, $T})"] = @benchmarkable $f($V2t)
            end
        end

        for (i, V4t) in enumerate((V4, V4sym))
            TensorType = i == 2 ? "SymmetricTensor" : "Tensor"
            for f in (norm, symmetric)
                SUITE["other"]["$f($TensorType{4, $dim, $T})"] = @benchmarkable $f($V4t)
            end
        end

        if T in (Float32, Float64)
            for f in (:+, :-)
                SUITE["basic-operations"]["Vec{$dim, $T} $f Vec{$dim, $T}"]                               = @benchmarkable $f($v1, $v1)
                SUITE["basic-operations"]["Tensor{2, $dim, $T} $f Tensor{2, $dim, $T}"]                   = @benchmarkable $f($V2, $V2)
                SUITE["basic-operations"]["SymmetricTensor{2, $dim, $T} $f SymmetricTensor{2, $dim, $T}"] = @benchmarkable $f($V2sym, $V2sym)
                SUITE["basic-operations"]["Tensor{4, $dim, $T} $f Tensor{4, $dim, $T}"]                   = @benchmarkable $f($V4, $V4)
                SUITE["basic-operations"]["SymmetricTensor{4, $dim, $T} $f SymmetricTensor{4, $dim, $T}"] = @benchmarkable $f($V4sym, $V4sym)
                if MIXED_SYM_NONSYM
                    SUITE["basic-operations"]["Tensor{2, $dim, $T} $f SymmetricTensor{2, $dim, $T}"]      = @benchmarkable $f($V2, $V2sym)
                    SUITE["basic-operations"]["Tensor{4, $dim, $T} $f SymmetricTensor{4, $dim, $T}"]      = @benchmarkable $f($V4, $V4sym)
                end
            end
            for f in (:*, :/)
                n = rand(T)
                SUITE["basic-operations"]["Vec{$dim, $T} $f $T"]                = @benchmarkable $f($v1, $n)
                SUITE["basic-operations"]["Tensor{2, $dim, $T} $f $T"]          = @benchmarkable $f($V2, $n)
                SUITE["basic-operations"]["SymmetricTensor{2, $dim, $T} $f $T"] = @benchmarkable $f($V2sym, $n)
                SUITE["basic-operations"]["Tensor{4, $dim, $T} $f $T"]          = @benchmarkable $f($V4, $n)
                SUITE["basic-operations"]["SymmetricTensor{4, $dim, $T} $f $T"] = @benchmarkable $f($V4sym, $n)
                if f == :*
                    SUITE["basic-operations"]["$T $f Vec{$dim, $T}"]                = @benchmarkable $f($n, $v1)
                    SUITE["basic-operations"]["$T $f Tensor{2, $dim, $T}"]          = @benchmarkable $f($n, $V2)
                    SUITE["basic-operations"]["$T $f SymmetricTensor{2, $dim, $T}"] = @benchmarkable $f($n, $V2sym)
                    SUITE["basic-operations"]["$T $f Tensor{4, $dim, $T}"]          = @benchmarkable $f($n, $V4)
                    SUITE["basic-operations"]["$T $f SymmetricTensor{4, $dim, $T}"] = @benchmarkable $f($n, $V4sym)
                end
            end
        end
    end
end

for dim in (ALL_DIMENSIONS ? (1,2,3) : (3,))
    a32 = tensor_dict[(dim, 1, Float32)]
    a64 = tensor_dict[(dim, 1, Float64)]
    A32 = tensor_dict[(dim, 2, Float32)]
    A64 = tensor_dict[(dim, 2, Float64)]
    AA32 = tensor_dict[(dim, 4, Float32)]
    AA64 = tensor_dict[(dim, 4, Float64)]
    A32s = symtensor_dict[(dim, 2, Float32)]
    A64s = symtensor_dict[(dim, 2, Float64)]
    AA32s = symtensor_dict[(dim, 4, Float32)]
    AA64s = symtensor_dict[(dim, 4, Float64)]

    # promotion between tensortypes
    SUITE["promotion"]["SymmetricTensor{2, $dim} -> Tensor{2, $dim}"] = @benchmarkable promote($A64, $A64s)
    SUITE["promotion"]["SymmetricTensor{4, $dim} -> Tensor{4, $dim}"] = @benchmarkable promote($AA64, $AA64s)

    # element type promotion
    SUITE["promotion"]["Vec{dim, Float32} -> Vec{$dim, Float64}"]                                = @benchmarkable promote($a64, $a32)
    SUITE["promotion"]["Tensor{2, $dim, Float32} -> Tensor{2, $dim, Float64}"]                   = @benchmarkable promote($A64, $A32)
    SUITE["promotion"]["SymmetricTensor{2, $dim, Float32} -> SymmetricTensor{2, $dim, Float64}"] = @benchmarkable promote($A64s, $A32s)
    SUITE["promotion"]["Tensor{4, $dim, Float32} -> Tensor{4, $dim, Float64}"]                   = @benchmarkable promote($A64, $A32)
    SUITE["promotion"]["SymmetricTensor{4, $dim, Float32} -> SymmetricTensor{4, $dim, Float64}"] = @benchmarkable promote($A64s, $A32s)

    # test just some operations with mixed tensortype and eltype and hope it catches things like
    # https://github.com/KristofferC/Tensors.jl/pull/5#issuecomment-282518974
    if MIXED_ELTYPES
        n = 5
        SUITE["promotion"]["Tensor{2, $dim, Float64} * $Int"]                     = @benchmarkable *($A64, $n)
        SUITE["promotion"]["Tensor{2, $dim, Float64} / $Int"]                     = @benchmarkable /($A64, $n)
        SUITE["promotion"]["Tensor{2, $dim, Float32} + Tensor{2, $dim, Float64}"] = @benchmarkable +($A32, $A64)
        SUITE["promotion"]["Tensor{2, $dim, Float32} - Tensor{2, $dim, Float64}"] = @benchmarkable -($A32, $A64)
        SUITE["promotion"]["Tensor{2, $dim, Float32} ⋅ Tensor{2, $dim, Float64}"] = @benchmarkable dot($A32, $A64)
        SUITE["promotion"]["Tensor{2, $dim, Float32} ⊡ Tensor{2, $dim, Float64}"] = @benchmarkable dcontract($A32, $A64)
        SUITE["promotion"]["Tensor{2, $dim, Float32} ⊗ Tensor{2, $dim, Float64}"] = @benchmarkable otimes($A32, $A64)
    end
end

# constructors (only testing for dim = 3)
# could be done cleaner, but https://github.com/JuliaCI/BenchmarkTools.jl/issues/50
for f in (:zero, :one, :ones, :rand)
    SUITE["constructors"]["$f(Tensor{2, 3, Float32})"]          = @benchmarkable $(f)(Tensor{2, 3, Float32})
    SUITE["constructors"]["$f(Tensor{4, 3, Float32})"]          = @benchmarkable $(f)(Tensor{4, 3, Float32})
    SUITE["constructors"]["$f(SymmetricTensor{2, 3, Float32})"] = @benchmarkable $(f)(SymmetricTensor{2, 3, Float32})
    SUITE["constructors"]["$f(SymmetricTensor{4, 3, Float32})"] = @benchmarkable $(f)(SymmetricTensor{4, 3, Float32})
    SUITE["constructors"]["$f(Tensor{2, 3, Float64})"]          = @benchmarkable $(f)(Tensor{2, 3, Float64})
    SUITE["constructors"]["$f(Tensor{4, 3, Float64})"]          = @benchmarkable $(f)(Tensor{4, 3, Float64})
    SUITE["constructors"]["$f(SymmetricTensor{2, 3, Float64})"] = @benchmarkable $(f)(SymmetricTensor{2, 3, Float64})
    SUITE["constructors"]["$f(SymmetricTensor{4, 3, Float64})"] = @benchmarkable $(f)(SymmetricTensor{4, 3, Float64})
    if f != :one
        SUITE["constructors"]["$f(Vec{3, Float32})"] = @benchmarkable $(f)(Vec{3, Float32})
        SUITE["constructors"]["$f(Vec{3, Float64})"] = @benchmarkable $(f)(Vec{3, Float64})
    end
end

# create from a Julia array
for order in (1, 2, 4)
    dim = 3
    A = rand(Tensors.n_components(Tensor{order, dim}))
    SUITE["constructors"]["Tensor{$order, $dim}(A::Array)"] = @benchmarkable Tensor{$order, $dim}($A)
    if order != 1
        As = rand(Tensors.n_components(SymmetricTensor{order, dim}))
        SUITE["constructors"]["SymmetricTensor{$order, $dim}(A::Array)"] = @benchmarkable SymmetricTensor{$order, $dim}($As)
    end
end

begin
    dim = 3
    f1 = (i) -> float(i)
    f2 = (i, j) -> float(i)
    f4 = (i, j, k, l) -> float(i)
    SUITE["constructors"]["Vec{3}(f::Function)"] =       @benchmarkable Vec{3}($f1)
    SUITE["constructors"]["Tensor{2, 3}(f::Function)"] = @benchmarkable Tensor{2, 3}($f2)
    SUITE["constructors"]["Tensor{4, 3}(f::Function)"] = @benchmarkable Tensor{4, 3}($f4)
    SUITE["constructors"]["SymmetricTensor{2, 3}(f::Function)"] = @benchmarkable SymmetricTensor{2, 3}($f2)
    SUITE["constructors"]["SymmetricTensor{4, 3}(f::Function)"] = @benchmarkable SymmetricTensor{4, 3}($f4)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module RouteFilterClient

using Random
using Dates
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_RouteFilterRulesApi.jl")
include("api_RouteFiltersApi.jl")

# export models
export convert, Access
export convert, CircuitConnectionStatus
export convert, CloudError
export convert, CloudErrorBody
export convert, ExpressRouteCircuitConnection
export convert, ExpressRouteCircuitConnectionPropertiesFormat
export convert, ExpressRouteCircuitPeering
export convert, ExpressRouteCircuitPeeringConfig
export convert, ExpressRouteCircuitPeeringPropertiesFormat
export convert, ExpressRouteCircuitStats
export convert, ExpressRouteConnectionId
export convert, ExpressRoutePeeringState
export convert, ExpressRoutePeeringType
export convert, Ipv6CircuitConnectionConfig
export convert, Ipv6ExpressRouteCircuitPeeringConfig
export convert, PatchRouteFilter
export convert, PatchRouteFilterRule
export convert, PeerExpressRouteCircuitConnection
export convert, PeerExpressRouteCircuitConnectionPropertiesFormat
export convert, ProvisioningState
export convert, Resource
export convert, RouteFilter
export convert, RouteFilterListResult
export convert, RouteFilterPropertiesFormat
export convert, RouteFilterRule
export convert, RouteFilterRuleListResult
export convert, RouteFilterRulePropertiesFormat
export convert, SubResource
export convert, TagsObject

# export operations
export convert, RouteFilterRulesApi, RouteFiltersApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
#=
Support same input format as asap2 sorting:
    - 2 records per core
    - Sort high to low: pass low keys.
    - Record is 5 words key, 45 word payload
    - 0x2A followed by 0 is the flush signal, which arrives between records.
    -- Treat any non-0 as flush.
    -- Post-flush value is the number of records to pass; 0 for the first core,
       but will generally be >0 within the chain.
    -- The final core will supress the flush and count signals, and just transmits
       output records.
    - 0 is the non-flush signal, indicates a new record to process.
AG layout:
    - AG0 points to low record key's low word, 5 words, stride backwards.
    - AG1 points to low record, 50 words, stride forward.
    - AG2 points to high record, 50 words, stride forward.
Storage:
    - 0: 1 if this is the first core in the chain, read single input.
    - 1: 1 if this is the last core in the chain, send to single output, dont pass count.
Set signed vs. unsigned key to configurable.
=#

function initialize_snakesort(core)
    # Initialize address generators
    AsapSim.set_stride!(core.address_generators[1], 0xFF)
    AsapSim.set_stride!(core.address_generators[2], 1)
    AsapSim.set_stride!(core.address_generators[3], 1)

    # Initialize CXMASK0
    AsapSim.set!(core.condexec[1], 1 << 2)
end

# Helper function for testing. Make a 50 word record with the given key.
function makerecord(key; zerofirst = false)
    if zerofirst
        return vcat(Int16.([0, 0,0,0,0,key]), rand(Int16, 45))
    else
        return vcat(Int16.([0,0,0,0,key]), rand(Int16, 45))
    end
end

#Function(Snakesort)
@asap4asm function snakesort(firstcore, lastcore, lastcore_flush)
    new_key_4 = Loc(AsapSim.DMEM, 1)
    new_key_3 = Loc(AsapSim.DMEM, 2)
    new_key_2 = Loc(AsapSim.DMEM, 3)
    new_key_1 = Loc(AsapSim.DMEM, 4)
    new_key_0 = Loc(AsapSim.DMEM, 5)
    num_records_to_pass = Loc(AsapSim.DMEM, 8)

    # Record 0 stored between 20 and 69
    # Key is 5 word range starting from low word of key and going backward.
    # So 24 down to 20.
    ag_config_full_record_0 = ((69 << 8) + 20)
    ag_config_key_0         = ((20 << 8 + 24))
    local_key_0_4 = Loc(AsapSim.DMEM, 20)
    local_key_0_3 = Loc(AsapSim.DMEM, 21)
    local_key_0_2 = Loc(AsapSim.DMEM, 22)
    local_key_0_1 = Loc(AsapSim.DMEM, 23)
    local_key_0_0 = Loc(AsapSim.DMEM, 24)

    # Record 1 stored between 70 and 119
    ag_config_full_record_1 = ((119 << 8) + 70)
    ag_config_key_1         = ((70 << 8) + 74)
    local_key_1_4 = Loc(AsapSim.DMEM, 70)
    local_key_1_3 = Loc(AsapSim.DMEM, 71)
    local_key_1_2 = Loc(AsapSim.DMEM, 72)
    local_key_1_1 = Loc(AsapSim.DMEM, 73)
    local_key_1_0 = Loc(AsapSim.DMEM, 74)

    low_key_ptr_pi  = :ag0pi
    low_key_ptr     = :ag0
    low_key_ptr_cfg = Loc(AsapSim.AG_START, 0)

    low_record_ptr_pi  = :ag1pi
    low_record_ptr     = :ag1
    low_record_ptr_cfg = Loc(AsapSim.AG_START, 1)

    high_record_ptr_pi  = :ag2pi
    high_record_ptr     = :ag
    high_record_ptr_cfg = Loc(AsapSim.AG_START, 2)

    # Begin main loop.

    # Start a new record batch, after reset.
    @label start
        # Peek at ibuf, check if the end of block code (any non-0) is found,
        #  and go to flush code if so.
        MOVE(null, ibuf0)
        BR(flush_0_local_record_with_count, Z, neg)
        # Otherwise not flushing

        # Reuse the storage to low function.  Cycle low pointer between
        #  memory blocks.
        MOVI(low_record_ptr_cfg, ag_config_full_record_0)
        # Set high to the same location; this is used if a flush signal is
        #  seen before storing another record.
        MOVI(high_record_ptr_cfg, ag_config_full_record_0, nop3) #TODO: nop2 this
        # Load an initial record.
        BRL(store_new_record_at_low, j)

        # Read ibuf, check if the end of block code (any non-0) is found,
        #  and go to flush code if so.
        MOVE(null, ibuf0)
        BR(flush_1_local_record_with_count, Z, neg)

        # Otherwise not flushing, so load a second record.
        # Set low pointer to next block.
        MOVI(low_record_ptr_cfg, ag_config_full_record_1, nop3)
        BRL(store_new_record_at_low, j)

        # Find which record is lower to set proper pointers.
        # This function jumps to get_new_key when done.
        BRL(compare_local_records)

    @label get_new_key
        # Peek at ibuf, check if the end of block code is found.
        # 0 means not end of block; non-zero is end of block.
        # Pass the control word to the next core either way.
        # Last core will not write out if not passing flush signals.
        if(lastcore == 1 && lastcore_flush == 0)
            MOVE(null, ibuf0)
        else
            MOVE(output[0], ibuf0)
        end
        BR(flush_2_local_record, Z, neg)

        # Get the next key.
        MOVE(new_key_4, ibuf0)
        MOVE(new_key_3, ibuf0)
        MOVE(new_key_2, ibuf0)
        MOVE(new_key_1, ibuf0)
        MOVE(new_key_0, ibuf0)

        # Compare to the lowest local key.
        # TODO: fast compare of high words and jump if new key is less, maybe.
        SUBU(null, low_key_ptr_pi, bypass[1])
        SUBCU(null, low_key_ptr_pi, bypass[3])
        SUBCU(null, low_key_ptr_pi, new_key_2)
        SUBCU(null, low_key_ptr_pi, new_key_3)
        #if(Unsigned0_Signed1){
        SUBC(null, low_key_ptr_pi, new_key_4, cxs0) #Set if carry if true.
        #}else{
        #SUBCU(null, low_key_ptr_pi, new_key_4, cxs0) #Set if carry if true.
        #}
        # If carry, new key is higher.

        # If new key lower or equal, pass it with record to output.
        # This function jumps to get_new_key when done.
        BRL(pass_input_to_output, C, neg)

        # If new key higher, send the local low record to output.
        # Store the new key and record at the previous low location.
        # Key.
        MOVE(output[0], low_record_ptr)
        MOVE(low_record_ptr_pi, new_key_4)
        MOVE(output[0], low_record_ptr)
        MOVE(low_record_ptr_pi, new_key_3)
        MOVE(output[0], low_record_ptr)
        MOVE(low_record_ptr_pi, new_key_2)
        MOVE(output[0], low_record_ptr)
        MOVE(low_record_ptr_pi, new_key_1)
        MOVE(output[0], low_record_ptr)
        MOVE(low_record_ptr_pi, new_key_0)

        # Rest of upper half.
        # 2x10=20
        RPT(10)
            MOVE(output[0], low_record_ptr)
            MOVE(low_record_ptr_pi, ibuf0)
            MOVE(output[0], low_record_ptr)
            MOVE(low_record_ptr_pi, ibuf0)
        END_RPT()

        # Lower half.
        # 2x12=24, add one extra at end for 25.
        RPT(12)
            if lastcore
                MOVE(output[0], low_record_ptr)
            else
                MOVE(output[1], low_record_ptr)
            end

            if firstcore
                MOVE(low_record_ptr_pi, ibuf0)
            else
                MOVE(low_record_ptr_pi, ibuf1)
            end

            if lastcore
                MOVE(output[0], low_record_ptr)
            else
                MOVE(output[1], low_record_ptr)
            end

            if firstcore
                MOVE(low_record_ptr_pi, ibuf0)
            else
                MOVE(low_record_ptr_pi, ibuf1)
            end
        END_RPT()
        # Copy of above.
        if lastcore
            MOVE(output[0], low_record_ptr)
        else
            MOVE(output[1], low_record_ptr)
        end

        if firstcore
            MOVE(low_record_ptr_pi, ibuf0)
        else
            MOVE(low_record_ptr_pi, ibuf1)
        end

    # Compare the two stored records to find which is now lower.
    # This block jumps to get_new_key when done.
    # Shared with startup code.
    @label compare_local_records
        # Check the keys against each other.
        SUBU(null, local_key_0_0, local_key_1_0)
        SUBCU(null, local_key_0_1, local_key_1_1)
        SUBCU(null, local_key_0_2, local_key_1_2)
        SUBCU(null, local_key_0_3, local_key_1_3)
        #if(Unsigned0_Signed1){
        SUBC(null, local_key_0_4, local_key_1_4, cxs0) #Set if carry is true.
        #}else{
        #SUBCU(null, local_key_0_4, local_key_1_4, cxs0) #Set if carry is true.
        #}

        # If carry set, key 0 is lower.
        MOVI(low_record_ptr_cfg, ag_config_full_record_0, cxt0)
        MOVI(low_key_ptr_cfg, ag_config_key_0, cxt0)
        MOVI(high_record_ptr_cfg, ag_config_full_record_1, cxt0)

        # Otherwise key 1 is lower.
        MOVI(low_record_ptr_cfg, ag_config_full_record_1, cxf0)
        MOVI(low_key_ptr_cfg, ag_config_key_1, cxf0)
        MOVI(high_record_ptr_cfg, ag_config_full_record_0, cxf0)

        # Ready to get the next record.
        BRL(get_new_key)

    @label pass_input_to_output
        MOVE(output[0], new_key_4)
        MOVE(output[0], new_key_3)
        MOVE(output[0], new_key_2)
        MOVE(output[0], new_key_1)
        MOVE(output[0], new_key_0)
        # 5x4=20
        RPT(5)
            MOVE(output[0], ibuf0)
            MOVE(output[0], ibuf0)
            MOVE(output[0], ibuf0)
            MOVE(output[0], ibuf0)
        END_RPT()
        # 5x5=25
        RPT(5)
            if firstcore && !lastcore
                MOVE(output[1], ibuf0)
                MOVE(output[1], ibuf0)
                MOVE(output[1], ibuf0)
                MOVE(output[1], ibuf0)
                MOVE(output[1], ibuf0)
            elseif (firstcore && lastcore)
                MOVE(output[0], ibuf0)
                MOVE(output[0], ibuf0)
                MOVE(output[0], ibuf0)
                MOVE(output[0], ibuf0)
                MOVE(output[0], ibuf0)
            elseif (!firstcore && lastcore)
                MOVE(output[0], ibuf1)
                MOVE(output[0], ibuf1)
                MOVE(output[0], ibuf1)
                MOVE(output[0], ibuf1)
                MOVE(output[0], ibuf1)
            else
                MOVE(output[1], ibuf1)
                MOVE(output[1], ibuf1)
                MOVE(output[1], ibuf1)
                MOVE(output[1], ibuf1)
                MOVE(output[1], ibuf1)
            end
        END_RPT()
        BRL(get_new_key)

    @label store_new_record_at_low
        # MOVE(low_record_ptr_pi, ibuf0)
        # MOVE(low_record_ptr_pi, ibuf0)
        # MOVE(low_record_ptr_pi, ibuf0)
        # MOVE(low_record_ptr_pi, ibuf0)
        # MOVE(low_record_ptr_pi, ibuf0)
        RPT(25, nop3)
            MOVE(low_record_ptr_pi, ibuf0)
        END_RPT()
        RPT(25, nop3)
            if firstcore
                MOVE(low_record_ptr_pi, ibuf0)
            else
                MOVE(low_record_ptr_pi, ibuf1)
            end
        END_RPT()
        BRL(back)

    @label flush_0_local_record_with_count
        # Get count and add 0 for records being passed.
        # Skip the output for hte last core when not passing flush signals.
        if lastcore == 0 || lastcore_flush == 1
            MOVE(num_records_to_pass, ibuf0_next)
            ADD(output[0], ibuf0, 0)
        else
            MOVE(num_records_to_pass, ibuf0)
        end
        BRL(flush_0_local_record, nop2)

    @label flush_1_local_record_with_count
        # Get count and add 1 for records being passed.
        if lastcore == 0 || lastcore_flush == 1
            MOVE(num_records_to_pass, ibuf0_next)
            ADD(output[0], ibuf0, 1)
        else
            MOVE(num_records_to_pass, ibuf0)
        end
        BRL(flush_1_local_record)

    @label flush_2_local_record
        # Get count and add 2 for records being passed.
        if lastcore == 0 || lastcore_flush == 1
            MOVE(num_records_to_pass, ibuf0_next)
            ADD(output[0], ibuf0, 2)
        else
            MOVE(num_records_to_pass, ibuf0)
        end

        # Send low record.
        # Not expanded due to imem limit.
        RPT(25, nop3)
            MOVE(output[0], low_record_ptr_pi)
        END_RPT()
        RPT(25, nop3)
            if lastcore
                MOVE(output[0], low_record_ptr_pi)
            else
                MOVE(output[1], low_record_ptr_pi)
            end
        END_RPT()

    @label flush_1_local_record
        # Send high record.
        # 5x5=25
        RPT(5)
            MOVE(output[0], high_record_ptr_pi)
            MOVE(output[0], high_record_ptr_pi)
            MOVE(output[0], high_record_ptr_pi)
            MOVE(output[0], high_record_ptr_pi)
            MOVE(output[0], high_record_ptr_pi)
        END_RPT()
        RPT(25, nop3)
            if lastcore
                MOVE(output[0], high_record_ptr_pi)
            else
                MOVE(output[1], high_record_ptr_pi)
            end
            END_RPT()

    @label flush_0_local_record
        # Pass records until count of how many to pass goes negative.
        # Maximum of ~32k records supported.
        # Start with a 0 check.
        MOVE(null, num_records_to_pass)
        BR(start, Z)

    @label pass_record
    # Loop unroll the repeats to remove the need for nops.
    # 5x5=25
    RPT(5)
        MOVE(output[0], ibuf0)
        MOVE(output[0], ibuf0)
        MOVE(output[0], ibuf0)
        MOVE(output[0], ibuf0)
        MOVE(output[0], ibuf0)
    END_RPT()
    # 5x5=25
    RPT(5)
        if firstcore && !lastcore
            MOVE(output[1], ibuf0)
            MOVE(output[1], ibuf0)
            MOVE(output[1], ibuf0)
            MOVE(output[1], ibuf0)
            MOVE(output[1], ibuf0)
        elseif firstcore && lastcore
            MOVE(output[0], ibuf0)
            MOVE(output[0], ibuf0)
            MOVE(output[0], ibuf0)
            MOVE(output[0], ibuf0)
            MOVE(output[0], ibuf0)
        elseif !firstcore && lastcore
            MOVE(output[0], ibuf1)
            MOVE(output[0], ibuf1)
            MOVE(output[0], ibuf1)
            MOVE(output[0], ibuf1)
            MOVE(output[0], ibuf1)
        else
            MOVE(output[1], ibuf1)
            MOVE(output[1], ibuf1)
            MOVE(output[1], ibuf1)
            MOVE(output[1], ibuf1)
            MOVE(output[1], ibuf1)
        end
    END_RPT()

    # Check count and loop.
    SUBU(num_records_to_pass, num_records_to_pass, 1)
    BRL(pass_record, Z, neg)
    BRL(start)

end

# Start_Initialization
#   # Verify storage size.
#   if(m->storage.size() != 3){
#       __debugbreak();
#   }
#
#   # Set up address generator strides.
#   MOVI(AG0STRIDE, 0xFF)
#   MOVI(AG1STRIDE, 1)
#   MOVI(AG2STRIDE, 1)
#
#   # Set condex 0 to look at carry flag.
#   MOVI(CXMASK0, 1<<COND_EXEC_MASK_CARRY_BIT)
#
# End_Initialization


text:
include("flows/normalising_flow.jl")
include("flows/inverse_normalising_flow.jl")


text:
using TableReader

if length(ARGS) != 1
    error("a function name is expected")
end

func = ARGS[1]
if func === "readcsv"
    @time readcsv("test/test.csv")
elseif func === "readtsv"
    @time readtsv("test/test.tsv")
else
    error("invalid function name: '$(func)'")
end


text:
import StatsBase: entropy
entropy(probs::ProbabilityWeights) = -sum(x->(x * log(x)), probs[probs .> 0])

function entropy(counts::FrequencyWeights; method=:Naive)
    method == :Naive && return entropy_naive(counts)
    method == :Shrink && return entropy_shrinkage(counts)
    method == :CS && return entropy_cs(counts)
    throw(ArgumentError("Unknown method $method"))
end

entropy_naive(counts::FrequencyWeights) =
    entropy(ProbabilityWeights(counts / sum(counts)))

"Chao-Shen (2003) entropy estimator."
function entropy_cs(counts::FrequencyWeights)
  n = sum(counts)
  θ_ML = counts / n

  f1 = sum(counts .== 1)
  f1 = (f1 == n) ? n - 1 : f1 # avoid C=0

  # Estimate coverage
  C = (1 - f1 / n)
  p_a = C * θ_ML
  l_a = (1 .- (1 .- p_a) .^ n)

  return - sum(p_a .* log.(p_a) ./ l_a)
end

#=
From http://www.jmlr.org/papers/volume10/hausser09a/hausser09a.pdf
=#
"Shrinkage entropy estimator."
function entropy_shrinkage(counts::FrequencyWeights)
    n = sum(counts)
    θ_ML = counts / n

    # Uniform distribution
    t_k = 1 / length(θ_ML)

    den = (n - 1) * sum((θ_ML .- t_k) .^2)
    if den < 1e-10
      return entropy(ProbabilityWeights(θ_ML))
    else
      # Regularization parameter
      λ = (1 - sum(θ_ML .^ 2)) / den
      return entropy(ProbabilityWeights(λ * t_k .+ (1 - λ) * θ_ML))
    end
end

"""
Estimate the entropy of an array using a naive (frequencies-based),
Chao-Shen, or shrinkage estimator. Chao-Shen and shrinkage estimators reduce
the bias for small samples and a large number of classes.
"""
function estimate_entropy(data::AbstractVector; method::Symbol=:Naive)
  count_values = values(countmap(data))
  freqs = FrequencyWeights(collect(count_values))
  return entropy(freqs; method=method)
end

function estimate_joint_entropy(x::AbstractVector, y::AbstractVector; method::Symbol=:Naive)
  @assert length(x) == length(y) "Vectors must be the same length"

  count = countmap([hash(yᵢ, hash(xᵢ)) for (xᵢ, yᵢ) ∈ zip(x, y)])
  freqs = FrequencyWeights(collect(values(count)))
  entropy(freqs; method=method)
end


text:
using Unitful
using DataFrames, DataFramesMeta
using CSV
# using Plots, StatsPlots

module OrganismUnits
    using Unitful
    # Reference unit for single entity
    @refunit entity "entity" Entity Unitful.𝐍 true

    # Units for specific organism types.  All are equal to one entity, but
    # the different types enforce apples/orages correctness.
    organisms = ["Acorn", "Bacteria", "Beaver", "Beetle", "Beetle", "Bird",
        "Caddisfly", "Caribou", "Caterpillar", "Cell", "Cockle", "Copepod", "Cladoceran", "Crab",
        "Ctenophore", "Elk", "Fish", "Fly", "Gazelle", "Insect", "Locust",
        "Molecule", "Mouse", "Mussel", "Nut", "Patch", "Plant", "Seal", "Seastar", "School",
        "Seed", "Shark", "Squirrel", "Superpatch", "Tree", "Tuna", "Turd", "Whale",
        "Wolf", "Worm"]
    for org in organisms
        abbr = lowercase(org)
        symb = Symbol(abbr)
        name = Symbol(org)
        @eval @unit $symb $abbr $name 1entity true
    end

    # a couple of (approximate) time units that aren't defined in base Unitful
    @unit yr "yr" Year 365.0u"d" false
    @unit mo "mo" Month 30.0u"d" false

    function __init__()
        Unitful.register(OrganismUnits)
    end
end
# using Main.OrganismUnits
string(unit(3.2 * u"elk/yr"))

# helper function to parse units straight from a string
function parseunits(u::AbstractString)
    @eval @u_str($u)
end
parseunits(u::Missing) = missing

raw = CSV.read(joinpath(@__DIR__, "data/PatchinessData_QC.csv"), DataFrame,
    missingstring="NA")
meta_vars = [:timestamp, :contributor, :consumer_resource_pair,
    :consumer_resource_abbrev,:genus, :species, :ecosystem,
    :interaction_system, :consumer_type]
consolidated = raw[!, meta_vars]

measurements = [:speed, :turning_interval, :generation_time, :consumption_rate,
    :mortality_rate, :consumer_body_size, :consumer_body_mass,
    :patch_duration, :patch_length_scale, :resource_body_size, :resource_body_mass,
    :resource_density, :consumer_density, :tsearch, :tconsumption]

for meas in measurements
    meas_units = Symbol(meas, "_units")
    consolidated[!, meas] = raw[!, meas] .* parseunits.(raw[!, meas_units])
end
consolidated[!, :patch_duration] = uconvert.(u"s", consolidated[!, :patch_duration])

consolidated = @transform(consolidated,
    :patch_length_scale = uconvert.(u"m", :patch_length_scale),
    :patch_duration = uconvert.(u"s", :patch_duration),
    :speed = uconvert.(u"m/s", :speed),
    :turning_interval = uconvert.(u"s", :turning_interval),
    :tsearch_diff = uconvert.(u"s", :patch_length_scale.^2 ./ (:speed.^2 .* :turning_interval)),
    :tsearch_dir = uconvert.(u"s", coalesce.(:tsearch, :patch_length_scale ./ :speed)),
    :treprod = uconvert.(u"s", :generation_time),
    :tconsumption = uconvert.(u"s",
        coalesce.(:tconsumption, :resource_density ./ (:consumer_density .* :consumption_rate))),
    :consumer_body_mass = uconvert.(u"g", :consumer_body_mass),
    :consumer_body_size = uconvert.(u"m", :consumer_body_size),
    :resource_body_mass = uconvert.(u"g", :resource_body_mass),
    :resource_body_size = uconvert.(u"m", :resource_body_size))

consolidated = @transform(consolidated,
    :Fr_diff = :patch_duration ./ :tsearch_diff,
    :Fr_dir = :patch_duration ./ :tsearch_dir,
    :Str = :patch_duration ./ :treprod,
    :Le = :patch_duration ./ :tconsumption)

output = consolidated[!, meta_vars]
for meas in [measurements; [:tsearch_diff, :tsearch_dir, :treprod]]
    meas_units = Symbol(meas, "_units")
    output[!, meas] = ustrip.(consolidated[!, meas])
    output[!, meas_units] = unit.(consolidated[!, meas])
end
output = [output consolidated[!, [:Fr_diff, :Fr_dir, :Str, :Le]]]
CSV.write(joinpath(@__DIR__, "data/PatchinessData_processed.csv"), output,
    missingstring="NA")

max_stars_repo_name: allicramer/EcologicalPatchiness

max_issues_repo_name: allicramer/EcologicalPatchiness

max_forks_repo_name: allicramer/EcologicalPatchiness


text:
mutable struct point_cloud_t <: LCMType
  utime::Int64       # unix time associated with this cloud
  id::Int64       # root pose id (sequential, 0-based)

  n::Int32           # number of points in the cloud
  points::Array{Float64,2} # [n][3] points (relative to the root pose)
  colors::Array{UInt8,2}     # [n][3] point colors (rgb)
end

@lcmtypesetup(point_cloud_t,
  points => (n, 3),
  colors => (n, 3)
)


text:
depths = parse.(Int, readlines("input.txt"))

#! Part 1

function countLarger(vec)
    isLarger = a -> a > 0 ? 1 : 0;
    sum([isLarger(vec[j] - vec[i]) for i in 1:length(vec) for j = i+1 if j <= length(vec)])
end

countLarger(depths)

#! Part 2

# A (1) -> 1, 2, 3
# B (2) -> 2, 3, 4
# C (3) -> 3, 4, 5
# ⋮

giveIndices = window -> [window + 1, window + 2, window + 3]

countLarger([sum(depths[giveIndices(i)]) for i in 1:length(depths)-2 if max(giveIndices(i)...) <= length(depths)])

max_stars_repo_path: Day 1/main.jl

max_issues_repo_path: Day 1/main.jl

max_forks_repo_path: Day 1/main.jl


text:
ENV["PYTHON"] = "F:/Anaconda/envs/ML/python.exe"
ENV["JUPYTER"] = "F:/Anaconda/envs/ML/Scripts/jupyter.exe"
ENV["JULIA_PARDISO"] = "F:/Code/PARDISO"
ENV["CMDSTAN"] = "C:/Users/12475/Documents/.cmdstanr/cmdstan-2.29.0"
try
    using Revise
catch e
    @warn "Error initializing Revise" exception = (e, catch_backtrace())
end
atreplinit() do repl
    try
        @eval using OhMyREPL
    catch e
        @warn "error while importing OhMyREPL" e
    end
end

atreplinit() do repl
    try
        @eval begin
            using Logging: global_logger
            using TerminalLoggers: TerminalLogger
            global_logger(TerminalLogger())
        end
    catch
    end
end


text:
"""
    This file specifically is meant for the reading of HDF4 files using Julia's
    Conda wrapper, because Julia does not have inbuilt support for HDF4 files.
    It was found necessary to create this file for accessing older datasets when
    HDF5 had not yet been released.

    The following datasets are in HDF4:
        - Tropical Rainfall Measuring Mission (TRMM)
        - MODerate resolution Imaging Spectroradiometer (MODIS)
"""

function h4read(file::AbstractString,variable::AbstractString)
    fhdf = pyhdf.SD(file,pyhdf.SDC.READ); dataob = fhdf.select(variable);
end


text:
using SpatialGraphs
using Base.Test

@test Point{Float64} <: Geometry{Float64}

@test Point(1, 2) + Point(2, 1) == Point(3, 3)
@test Point(1.0, 2.0) + Point(2, 1) == Point(3.0, 3.0)
@test Point(1, 2) * 2 == Point(2, 4)
@test Point(1, 2) * 2.0 == Point(2.0, 4.0)
@test 2 * Point(1, 2) == Point(2.0, 4.0)

SimpleSpatialVertex(1, Point(1.2, 3.5))
SimpleSpatialVertex(1, 1.2, 3.5)


text:
include("../PlotUtils/PlotUtils.jl")
include("../PlotUtils/imports.jl")

# using Distributed
# rmprocs(filter(w -> w > 1, workers()))
# addprocs(2)

if length(ARGS) == 0
  results_dir = "/scratchdata/alui2/cytof/results/repfam/test-sim-6-7-3"
else
  results_dir = ARGS[1]  # path to results directory
end

"""
Returns a list of monitors, where each element is the result from
one temperature.
"""
function getsamples(path_to_output)
  output = BSON.load(path_to_output)
  samples = output[:samples]
  num_tempers = length(output[:tempers])
  num_monitors = length(samples)

  return [let
     monitor = [let
                  num_samples = length(output[:samples][m])
                  [output[:samples][m][i][t] for i in 1:num_samples]
                end for m in 1:num_monitors]
   end for t in 1:num_tempers]
end

# @everywhere include("../PlotUtils/PlotUtils.jl")
# @everywhere include("../PlotUtils/imports.jl")

# @everywhere function makeplots(path_to_output)
function makeplots(path_to_output)
  output_dir = splitdir(path_to_output)[1]
  path_to_simdat = joinpath(output_dir, "simdat.bson")
  simdat = BSON.load(path_to_simdat)[:simdat]

  println(path_to_output)
  output = BSON.load(path_to_output)

  # Extract samples
  extract(s) = map(o -> o[s], output[:samples][1])
  Zs = extract(:theta__Z)
  Ws = extract(:theta__W)
  lams = extract(:theta__lam)
  alphas = extract(:theta__alpha)
  vs = extract(:theta__v)
  deltas = extract(:theta__delta)
  sig2s = extract(:theta__sig2)
  etas = extract(:theta__eta)

  # Number of samples
  I = length(sig2s[1])

  # Create a directory for images / txt if needed.
  dir_to_output, _ = splitdir(path_to_output)
  imgdir = joinpath(dir_to_output, "img")
  mkpath("$(imgdir)/txt/")

  # Print R
  computeR(w) = vec(sum(w .> 0, dims=2))
  open("$(imgdir)/txt/Rcounts.txt", "w") do io
    I, K = size(Ws[1])
    Rs = map(computeR, Ws)
    Rcounts = zeros(Int, I, K)
    for R in Rs
      for i in 1:I
        Rcounts[i, R[i]] += 1
      end
    end
    writedlm(io, Rcounts)
  end

  # For parallel tempering
  if :swapcounts in keys(output) && :paircounts in keys(output)
    swapcounts = output[:swapcounts]
    paircounts = output[:paircounts]
    tempers = output[:tempers]
    ntempers = length(tempers)
    plt.imshow(swapcounts ./ (paircounts .+ 1e-6))
    plt.colorbar()
    plt.xticks(0:(ntempers-1), round.(tempers, digits=2), rotation=45)
    plt.yticks(0:(ntempers-1), round.(tempers, digits=2))
    plt.savefig("$(imgdir)/swapprops.pdf", bbox_inches="tight")
    plt.close()
  end

  # Number of burn in iterations
  nburn = output[:nburn]

  # Plot loglikelihood
  llkey = if :loglike in keys(output)
    PlotUtils.plot_loglike(output[:loglike], imgdir, fname="loglike.pdf")
    PlotUtils.plot_loglike(output[:loglike][(nburn + 1):end],
                           imgdir, fname="loglike_postburn.pdf")
  elseif :lls in keys(output)
    lls = output[:lls]
    for t in 1:ntempers
      temper = tempers[t]
      title = "temperature: $(temper)"
      ll_path = joinpath(imgdir, "loglike")
      mkpath(ll_path)
      PlotUtils.plot_loglike(lls[t], ll_path, fname="loglike_$(t).pdf",
                             title=title)
      PlotUtils.plot_loglike(lls[t][(nburn + 1):end],
                             ll_path, fname="loglike_postburn_$(t).pdf",
                             title=title)
    end
  end

  # Plot missing mechanism
  for i in 1:I
    PlotUtils.plot_missmech(output[:c].constants.beta, i, xlim=[-5, 0])
    plt.savefig("$(imgdir)/missmech_$(i).pdf", bbox_inches="tight")
    plt.close()
  end

  # Print missing mechanism
  open(joinpath(imgdir, "txt/beta.txt"), "w") do io
    writedlm(io, output[:c].constants.beta, ',')
  end

  # Plot parameters
  PlotUtils.plot_W(Ws, imgdir=imgdir, W_true=simdat[:W])
  PlotUtils.plot_v(vs, imgdir)
  PlotUtils.plot_alpha(alphas, imgdir)
  PlotUtils.plot_mus(deltas, imgdir)
  PlotUtils.plot_sig2(sig2s, imgdir, sig2_true=simdat[:sig2])

  # TODO: Plot data density
  PlotUtils.plot_dden(ddens=output[:dden],
                      etas=etas, Ws=Ws, Zs=Zs, sig2s=sig2s, deltas=deltas,
                      ygrid=output[:c].constants.y_grid,
                      imgdir=imgdir, simdat=simdat)

  # Plot y / Z
  PlotUtils.make_yz(simdat[:y], Zs, Ws, lams, imgdir, vlim=(-4,4),
                    Z_true=simdat[:Z], w_thresh=0.0)
end

### MAIN ###

# Name of output file
OUTPUT_FILE = "output.bson"

# PATH TO ALL OUTPUT FILES
output_paths = [joinpath(root, OUTPUT_FILE)
                for (root, _, files) in walkdir(results_dir)
                if OUTPUT_FILE in files]

# Make y, Z plots.
status = map(makeplots, output_paths)
println(status)

# Remove extra processors
# rmprocs(filter(w -> w > 1, workers()))

println("DONE!")


text:
using Chron
using Test, Statistics, StatsBase

@testset "Utilities" begin include("testUtilities.jl") end
@testset "Eruption / deposition age distributions" begin include("testDist.jl") end
@testset "Strat only" begin include("testStratOnly.jl") end
@testset "Radiocarbon" begin include("testRadiocarbon.jl") end
@testset "Coupled model" begin include("testCoupled.jl") end


text:
function make_MF_vertex(node::DecisionTree.Node)
    return MFNode(node.featid, node.featval, make_MF_vertex(node.left), make_MF_vertex(node.right))
end
function make_MF_vertex(leaf::DecisionTree.Leaf)
    return MFLeaf()
end

function get_MF_leaf(μ::MFMean, p::MFIteration)
    return get_MF_leaf(μ, p.θ, p.y_lo)
end
function get_MF_leaf(μ::MFMean, θ, y_lo)
    return get_MF_leaf(μ, vcat(θ, y_lo))
end
function get_MF_leaf(μ::MFMean, features)
    return get_MF_leaf(μ.root, features)
end
function get_MF_leaf(node::MFNode, features)
    child = features[node.featid] < node.featval ? node.left : node.right
    return get_MF_leaf(child, features)
end
function get_MF_leaf(leaf::MFLeaf, features)
    return leaf
end

Base.length(node::MFNode) = length(node.left) + length(node.right)
Base.length(::MFLeaf) = 1

leaves(node::MFNode) = Iterators.flatten((leaves(node.left), leaves(node.right)))
leaves(leaf::MFLeaf)::Tuple{MFLeaf} = (leaf, )
leaves(μ::MFMean) = leaves(μ.root)

function (μ::MFMean)(p)
    l = get_MF_leaf(μ, p)
    return exp(l.logμ)
end
(μ::MFMean)(θ, y_lo) = μ(vcat(θ, y_lo))

function DecisionTree.print_tree(io::IO, leaf::MFLeaf, _k::Int64, depth=-1, indent=0)
    print(io, "\\State \\Return \$\\nu_{$(_k)} = $(round(exp(leaf.logμ), digits=3))\$.\n")
    return _k+1
end
function DecisionTree.print_tree(io::IO, node::MFNode, _k::Int64, depth=-1, indent=0)
    if depth==indent
        print(io, "\n")
        return
    end

    θ_flag = (node.featid <= 3)
    _n = θ_flag ? node.featid : (node.featid - 3)

    if θ_flag
        pars = ["k_{1}", "k_{-1}", "k_2"]
        if_str = "\\If{\$ $(pars[node.featid]) \\leq $(round(node.featval, digits=3)) \$\n}"
    else
        if_str = "\\If{\$ y_{$(_n)} \\leq $(round(node.featval, digits=3)) \$}\n"
    end
    print(io, if_str)
    _k = print_tree(io, node.left, _k, depth, indent+1)
    print(io, "\\Else\n")
    _k = print_tree(io, node.right, _k, depth, indent+1)
    print(io, "\\EndIf\n")
    return _k
end
function DecisionTree.print_tree(io::IO, μ::MFMean, args...)
    print(io, "\\begin{algorithmic}\n")
    print(io, "\\Require \$\\theta=(k_1, k_{-1}, k_2)\$; \$y_{\\lo} = (y_1, y_2, \\dots, y_{10})\$.\n")
    print_tree(io, μ.root, 1, args...)
    print(io, "\\end{algorithmic}\n")
    return true
end

max_stars_repo_name: tpprescott/mf-lf

max_issues_repo_name: tpprescott/mf-lf

max_forks_repo_name: tpprescott/mf-lf


text:
function _asKs(; Ks=nothing, k=nothing, ρ=1e3, μ=1e-3, g=9.81)
    if !isnothing(Ks)
        @argcheck isnothing(k) "may only assign one of Ks and k, got both"
        @argcheck Ks>zero(Ks)
        return Ks

    elseif !isnothing(k)
        @argcheck k>zero(k); @argcheck ρ>zero(ρ); @argcheck μ>zero(μ); @argcheck g>zero(g)
        return ρ*g*k/μ

    else
        return 1
    end
end


text:
a = 3

# points to current directory
@__DIR__

# Points to where julia gets started
pwd()

module MyModule
function divide(a, b)
    return a*b
end
end

MyModule.divide(6,3)


text:
import IterTools.filter

function sumby_intrange{T<:Integer, S}(by::Vector{T}, val::Vector{S}, bylow::T, byhigh::T)
  lby = byhigh - bylow + 1
  res = zeros(T, lby)

  for (byi,vali) in zip(by, val)
    @inbounds res[(byi-bylow) + 1] += vali
  end
  IterTools.filter(x -> x[2] != 0, zip(bylow:byhigh,res))
end
@time res = sumby_intrange(id6, v1, 1, 2_500_000); #7 seconds


text:
let
    S = CBSEnv.State
    A = CBSEnv.Action
    P = PathNode{S,A}
    vtx_grid = initialize_dense_vtx_grid(4,4)
    # 1  5   9  13
    # 2  6  10  14
    # 3  7  11  15
    # 4  8  12  16
    G = initialize_grid_graph_from_vtx_grid(vtx_grid)
    starts = [CBSEnv.State(1,0),CBSEnv.State(2,0)]
    goals = [CBSEnv.State(vtx=5),CBSEnv.State(vtx=6)]
    heuristic = PerfectHeuristic(G,map(s->s.vtx,starts),map(s->s.vtx,goals))
    env = CBSEnv.LowLevelEnv(graph=G,heuristic=heuristic)
    mapf = MAPF(env,starts,goals)
    let
        root_node = initialize_root_node(mapf)
        paths = [
            Path([P(S(1,1),A(Edge(1,2),2),S(2,3))]),
            Path([P(S(2,1),A(Edge(2,1),2),S(1,3))]),
        ]
        table = detect_conflicts(paths)
        @test CRCBS.is_action_conflict(get_next_conflict(table))
        set_solution_path!(root_node.solution,paths[1],1)
        set_solution_path!(root_node.solution,paths[2],2)
        detect_conflicts!(root_node)
        @test count_conflicts(root_node) == 1
        @test count_conflicts(root_node,1,2) == count_conflicts(root_node,2,1)
        @test count_conflicts(root_node,[1,2],[1,2]) == 1
        c = get_next_conflict(root_node)
        constraints = generate_constraints_from_conflict(c)
        @test length(constraints) == 2
        for c in constraints
            @test is_action_constraint(c)
        end
        add_constraint!(env,root_node,constraints[1])
    end
    let
        root_node = initialize_root_node(mapf)
        paths = [
            Path([P(S(1,0),A(Edge(1,2),0),S(2,1))]),
            Path([P(S(2,0),A(Edge(2,2),0),S(2,1))]),
        ]
        set_solution_path!(root_node.solution,paths[1],1)
        set_solution_path!(root_node.solution,paths[2],2)
        child1 = initialize_child_search_node(root_node)
        child2 = initialize_child_search_node(root_node)
        detect_conflicts!(child1,[1])
        @test count_conflicts(child1) == 1
        @test count_conflicts(child1,1,2) == 1
        @test count_conflicts(child1,2,1) == 1
        @test count_conflicts(child1,[1,2],[1,2]) == 1

        c = get_next_conflict(child1)
        constraints = generate_constraints_from_conflict(c)
        @test length(constraints) == 2
        for c in constraints
            @test is_state_constraint(c)
        end
        add_constraint!(env,child1,constraints[1])

        @test count_conflicts(child2) == 0
        @test count_conflicts(root_node) == 0

        @test convert_to_vertex_lists(
            get_paths(child1.solution)[1]
            ) == convert_to_vertex_lists(
            get_paths(root_node.solution)[1]
            )

        set_solution_path!(child1.solution,Path([P()]),1)
        @test convert_to_vertex_lists(
            get_paths(child1.solution)[1]
            ) != convert_to_vertex_lists(
            get_paths(root_node.solution)[1]
            )

        detect_conflicts!(child1,[1])
        @test count_conflicts(child1) == 0
    end
end
let
    S = CBSEnv.State
    A = CBSEnv.Action
    P = PathNode{S,A}

    vtx_grid = initialize_dense_vtx_grid(4,4)
    # 1  5   9  13
    # 2  6  10  14
    # 3  7  11  15
    # 4  8  12  16
    G = initialize_grid_graph_from_vtx_grid(vtx_grid)
    env = CBSEnv.LowLevelEnv(graph=G,agent_idx=1)
    isa(state_space_trait(env),DiscreteSpace)
    let
        for v in vertices(G)
            s = S(v,1)
            idx,t = serialize(env,s,s.t)
            sp, _ = deserialize(env,S(),idx,t)
            @test get_vtx(s) == get_vtx(sp)
        end
        for e in edges(G)
            a = A(e,1)
            idx,t = serialize(env,a,1)
            ap, _ = deserialize(env,A(),idx,t)
            @test get_e(a) == get_e(ap)
        end
    end
    let
        for constraints in [
            discrete_constraint_table(env,1),
            # ConstraintTable{P}(agent_id = 1)
            ]
            # add a constraint whose agent id does not match (should throw an error)
            s = S(1,1)
            a = A(Edge(1,2),1)
            sp = S(2,2)
            c = state_constraint(get_agent_id(constraints)+1,P(),2)
            @test_throws AssertionError add_constraint!(env,constraints,c) # wrong agent id
            @test_throws AssertionError CRCBS.has_constraint(env,constraints,c) # wrong agent id
            # add a constraint whose agent id DOES match
            c = state_constraint(get_agent_id(constraints),P(s,a,sp),2)
            add_constraint!(env,constraints,c)
            @test CRCBS.has_constraint(env,constraints,c)
            remove_constraint!(env,constraints,c)
            @test !CRCBS.has_constraint(env,constraints,c)
            add_constraint!(env,constraints,c)
            sc,ac = CRCBS.search_constraints(env,constraints,get_path_node(c))
            @test length(sc) == 1
            @test length(ac) == 0
            # # add a constraint whose agent id DOES match
            c = action_constraint(get_agent_id(constraints),P(s,a,sp),2)
            add_constraint!(env,constraints,c)
            @test CRCBS.has_constraint(env,constraints,c)
            remove_constraint!(env,constraints,c)
            @test !CRCBS.has_constraint(env,constraints,c)
            add_constraint!(env,constraints,c)
            sc,ac = CRCBS.search_constraints(env,constraints,get_path_node(c))
            @test length(sc) == 1
            @test length(ac) == 1
            # # query sorted constraints
            c = sorted_state_constraints(env,constraints)[1]
            @test get_vtx(get_sp(c)) == get_vtx(sp)
            c = sorted_action_constraints(env,constraints)[1]
            @test get_e(get_a(c)) == get_e(a)
            # Test interval constraints
            c = state_constraint(get_agent_id(constraints),P(s,a,sp),3,4)
            add_constraint!(env,constraints,c)
            for t in 3:4
                @test CRCBS.has_constraint(env,constraints,state_constraint(get_agent_id(constraints),P(s,a,sp),t))
            end
            remove_constraint!(env,constraints,c)
            for t in 3:4
                @test !CRCBS.has_constraint(env,constraints,state_constraint(get_agent_id(constraints),P(s,a,sp),t))
            end

        end
    end
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((a+a*sin(c+d*x))^7//2, x) == :(-256//35*a^4*d^-1*(a+a*sin(c+d*x))^-1//2*cos(c+d*x)+-64//35*a^3*d^-1*(a+a*sin(c+d*x))^(1/2)*cos(c+d*x)+-24//35*a^2*d^-1*(a+a*sin(c+d*x))^3//2*cos(c+d*x)+-2//7*a*d^-1*(a+a*sin(c+d*x))^5//2*cos(c+d*x))
@test integrate((a+a*sin(c+d*x))^5//2, x) == :(-64//15*a^3*d^-1*(a+a*sin(c+d*x))^-1//2*cos(c+d*x)+-16//15*a^2*d^-1*(a+a*sin(c+d*x))^(1/2)*cos(c+d*x)+-2//5*a*d^-1*(a+a*sin(c+d*x))^3//2*cos(c+d*x))
@test integrate((a+a*sin(c+d*x))^3//2, x) == :(-8//3*a^2*d^-1*(a+a*sin(c+d*x))^-1//2*cos(c+d*x)+-2//3*a*d^-1*(a+a*sin(c+d*x))^(1/2)*cos(c+d*x))
@test integrate((a+a*sin(c+d*x))^(1/2), x) == :(-2*a*d^-1*(a+a*sin(c+d*x))^-1//2*cos(c+d*x))
@test integrate((a+a*sin(c+d*x))^-1//2, x) == :(-1*2^(1/2)*a^-1//2*d^-1*arctanh((1/2)*2^(1/2)*a^(1/2)*(a+a*sin(c+d*x))^-1//2*cos(c+d*x)))
@test integrate((a+a*sin(c+d*x))^-3//2, x) == :(-1//2*d^-1*(a+a*sin(c+d*x))^-3//2*cos(c+d*x)+-1//4*2^(1/2)*a^-3//2*d^-1*arctanh((1/2)*2^(1/2)*a^(1/2)*(a+a*sin(c+d*x))^-1//2*cos(c+d*x)))
@test integrate((a+a*sin(c+d*x))^-5//2, x) == :(-1//4*d^-1*(a+a*sin(c+d*x))^-5//2*cos(c+d*x)+-3//16*a^-1*d^-1*(a+a*sin(c+d*x))^-3//2*cos(c+d*x)+-3//32*2^(1/2)*a^-5//2*d^-1*arctanh((1/2)*2^(1/2)*a^(1/2)*(a+a*sin(c+d*x))^-1//2*cos(c+d*x)))
@test integrate((5+3*sin(c+d*x))^-1, x) == :(1//4*x+(1/2)*d^-1*arctan((3+sin(c+d*x))^-1*cos(c+d*x)))
@test integrate((5+3*sin(c+d*x))^-2, x) == :(5//64*x+5//32*d^-1*arctan((3+sin(c+d*x))^-1*cos(c+d*x))+3//16*d^-1*(5+3*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((5+3*sin(c+d*x))^-3, x) == :(59//2048*x+59//1024*d^-1*arctan((3+sin(c+d*x))^-1*cos(c+d*x))+3//32*d^-1*(5+3*sin(c+d*x))^-2*cos(c+d*x)+45//512*d^-1*(5+3*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((5+3*sin(c+d*x))^-4, x) == :(385//32768*x+385//16384*d^-1*arctan((3+sin(c+d*x))^-1*cos(c+d*x))+1//16*d^-1*(5+3*sin(c+d*x))^-3*cos(c+d*x)+25//512*d^-1*(5+3*sin(c+d*x))^-2*cos(c+d*x)+311//8192*d^-1*(5+3*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((5+-3*sin(c+d*x))^-1, x) == :(1//4*x+-1//2*d^-1*arctan((3+-1*sin(c+d*x))^-1*cos(c+d*x)))
@test integrate((5+-3*sin(c+d*x))^-2, x) == :(5//64*x+-5//32*d^-1*arctan((3+-1*sin(c+d*x))^-1*cos(c+d*x))+-3//16*d^-1*(5+-3*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((5+-3*sin(c+d*x))^-3, x) == :(59//2048*x+-59//1024*d^-1*arctan((3+-1*sin(c+d*x))^-1*cos(c+d*x))+-45//512*d^-1*(5+-3*sin(c+d*x))^-1*cos(c+d*x)+-3//32*d^-1*(5+-3*sin(c+d*x))^-2*cos(c+d*x))
@test integrate((5+-3*sin(c+d*x))^-4, x) == :(385//32768*x+-385//16384*d^-1*arctan((3+-1*sin(c+d*x))^-1*cos(c+d*x))+-311//8192*d^-1*(5+-3*sin(c+d*x))^-1*cos(c+d*x)+-25//512*d^-1*(5+-3*sin(c+d*x))^-2*cos(c+d*x)+-1//16*d^-1*(5+-3*sin(c+d*x))^-3*cos(c+d*x))
@test integrate((-5+3*sin(c+d*x))^-1, x) == :(-1//4*x+(1/2)*d^-1*arctan((3+-1*sin(c+d*x))^-1*cos(c+d*x)))
@test integrate((-5+3*sin(c+d*x))^-2, x) == :(5//64*x+-5//32*d^-1*arctan((3+-1*sin(c+d*x))^-1*cos(c+d*x))+-3//16*d^-1*(5+-3*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((-5+3*sin(c+d*x))^-3, x) == :(-59//2048*x+59//1024*d^-1*arctan((3+-1*sin(c+d*x))^-1*cos(c+d*x))+3//32*d^-1*(5+-3*sin(c+d*x))^-2*cos(c+d*x)+45//512*d^-1*(5+-3*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((-5+3*sin(c+d*x))^-4, x) == :(385//32768*x+-385//16384*d^-1*arctan((3+-1*sin(c+d*x))^-1*cos(c+d*x))+-311//8192*d^-1*(5+-3*sin(c+d*x))^-1*cos(c+d*x)+-25//512*d^-1*(5+-3*sin(c+d*x))^-2*cos(c+d*x)+-1//16*d^-1*(5+-3*sin(c+d*x))^-3*cos(c+d*x))
@test integrate((-5+-3*sin(c+d*x))^-1, x) == :(-1//4*x+-1//2*d^-1*arctan((3+sin(c+d*x))^-1*cos(c+d*x)))
@test integrate((-5+-3*sin(c+d*x))^-2, x) == :(5//64*x+5//32*d^-1*arctan((3+sin(c+d*x))^-1*cos(c+d*x))+3//16*d^-1*(5+3*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((-5+-3*sin(c+d*x))^-3, x) == :(-59//2048*x+-59//1024*d^-1*arctan((3+sin(c+d*x))^-1*cos(c+d*x))+-45//512*d^-1*(5+3*sin(c+d*x))^-1*cos(c+d*x)+-3//32*d^-1*(5+3*sin(c+d*x))^-2*cos(c+d*x))
@test integrate((-5+-3*sin(c+d*x))^-4, x) == :(385//32768*x+385//16384*d^-1*arctan((3+sin(c+d*x))^-1*cos(c+d*x))+1//16*d^-1*(5+3*sin(c+d*x))^-3*cos(c+d*x)+25//512*d^-1*(5+3*sin(c+d*x))^-2*cos(c+d*x)+311//8192*d^-1*(5+3*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((3+5*sin(c+d*x))^-1, x) == :(-1//4*d^-1*log(3*cos((1/2)*c+(1/2)*d*x)+sin((1/2)*c+(1/2)*d*x))+1//4*d^-1*log(3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x)))
@test integrate((3+5*sin(c+d*x))^-2, x) == :(-3//64*d^-1*log(3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x))+3//64*d^-1*log(3*cos((1/2)*c+(1/2)*d*x)+sin((1/2)*c+(1/2)*d*x))+-5//16*d^-1*(3+5*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((3+5*sin(c+d*x))^-3, x) == :(-43//2048*d^-1*log(3*cos((1/2)*c+(1/2)*d*x)+sin((1/2)*c+(1/2)*d*x))+43//2048*d^-1*log(3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x))+-5//32*d^-1*(3+5*sin(c+d*x))^-2*cos(c+d*x)+45//512*d^-1*(3+5*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((3+5*sin(c+d*x))^-4, x) == :(-279//32768*d^-1*log(3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x))+279//32768*d^-1*log(3*cos((1/2)*c+(1/2)*d*x)+sin((1/2)*c+(1/2)*d*x))+-995//24576*d^-1*(3+5*sin(c+d*x))^-1*cos(c+d*x)+-5//48*d^-1*(3+5*sin(c+d*x))^-3*cos(c+d*x)+25//512*d^-1*(3+5*sin(c+d*x))^-2*cos(c+d*x))
@test integrate((3+-5*sin(c+d*x))^-1, x) == :(-1//4*d^-1*log(-3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x))+1//4*d^-1*log(-1*sin((1/2)*c+(1/2)*d*x)+3*cos((1/2)*c+(1/2)*d*x)))
@test integrate((3+-5*sin(c+d*x))^-2, x) == :(-3//64*d^-1*log(-1*sin((1/2)*c+(1/2)*d*x)+3*cos((1/2)*c+(1/2)*d*x))+3//64*d^-1*log(-3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x))+5//16*d^-1*(3+-5*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((3+-5*sin(c+d*x))^-3, x) == :(-43//2048*d^-1*log(-3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x))+43//2048*d^-1*log(-1*sin((1/2)*c+(1/2)*d*x)+3*cos((1/2)*c+(1/2)*d*x))+-45//512*d^-1*(3+-5*sin(c+d*x))^-1*cos(c+d*x)+5//32*d^-1*(3+-5*sin(c+d*x))^-2*cos(c+d*x))
@test integrate((3+-5*sin(c+d*x))^-4, x) == :(-279//32768*d^-1*log(-1*sin((1/2)*c+(1/2)*d*x)+3*cos((1/2)*c+(1/2)*d*x))+279//32768*d^-1*log(-3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x))+-25//512*d^-1*(3+-5*sin(c+d*x))^-2*cos(c+d*x)+5//48*d^-1*(3+-5*sin(c+d*x))^-3*cos(c+d*x)+995//24576*d^-1*(3+-5*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((-3+5*sin(c+d*x))^-1, x) == :(-1//4*d^-1*log(-1*sin((1/2)*c+(1/2)*d*x)+3*cos((1/2)*c+(1/2)*d*x))+1//4*d^-1*log(-3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x)))
@test integrate((-3+5*sin(c+d*x))^-2, x) == :(-3//64*d^-1*log(-1*sin((1/2)*c+(1/2)*d*x)+3*cos((1/2)*c+(1/2)*d*x))+3//64*d^-1*log(-3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x))+5//16*d^-1*(3+-5*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((-3+5*sin(c+d*x))^-3, x) == :(-43//2048*d^-1*log(-1*sin((1/2)*c+(1/2)*d*x)+3*cos((1/2)*c+(1/2)*d*x))+43//2048*d^-1*log(-3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x))+-5//32*d^-1*(3+-5*sin(c+d*x))^-2*cos(c+d*x)+45//512*d^-1*(3+-5*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((-3+5*sin(c+d*x))^-4, x) == :(-279//32768*d^-1*log(-1*sin((1/2)*c+(1/2)*d*x)+3*cos((1/2)*c+(1/2)*d*x))+279//32768*d^-1*log(-3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x))+-25//512*d^-1*(3+-5*sin(c+d*x))^-2*cos(c+d*x)+5//48*d^-1*(3+-5*sin(c+d*x))^-3*cos(c+d*x)+995//24576*d^-1*(3+-5*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((-3+-5*sin(c+d*x))^-1, x) == :(-1//4*d^-1*log(3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x))+1//4*d^-1*log(3*cos((1/2)*c+(1/2)*d*x)+sin((1/2)*c+(1/2)*d*x)))
@test integrate((-3+-5*sin(c+d*x))^-2, x) == :(-3//64*d^-1*log(3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x))+3//64*d^-1*log(3*cos((1/2)*c+(1/2)*d*x)+sin((1/2)*c+(1/2)*d*x))+-5//16*d^-1*(3+5*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((-3+-5*sin(c+d*x))^-3, x) == :(-43//2048*d^-1*log(3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x))+43//2048*d^-1*log(3*cos((1/2)*c+(1/2)*d*x)+sin((1/2)*c+(1/2)*d*x))+-45//512*d^-1*(3+5*sin(c+d*x))^-1*cos(c+d*x)+5//32*d^-1*(3+5*sin(c+d*x))^-2*cos(c+d*x))
@test integrate((-3+-5*sin(c+d*x))^-4, x) == :(-279//32768*d^-1*log(3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x))+279//32768*d^-1*log(3*cos((1/2)*c+(1/2)*d*x)+sin((1/2)*c+(1/2)*d*x))+-995//24576*d^-1*(3+5*sin(c+d*x))^-1*cos(c+d*x)+-5//48*d^-1*(3+5*sin(c+d*x))^-3*cos(c+d*x)+25//512*d^-1*(3+5*sin(c+d*x))^-2*cos(c+d*x))
@test integrate((a+b*sin(c+d*x))^7//2, x) == :(-2//7*b*d^-1*(a+b*sin(c+d*x))^5//2*cos(c+d*x)+-24//35*a*b*d^-1*(a+b*sin(c+d*x))^3//2*cos(c+d*x)+-2//105*b*d^-1*(a+b*sin(c+d*x))^(1/2)*(25*b^2+71*a^2)*cos(c+d*x)+-1//105*d^-1*((a+b)^-1*(a+b*sin(c+d*x)))^(1/2)*(a+b*sin(c+d*x))^-1//2*(-50*b^4+142*a^4+-92*a^2*b^2)*Elliptic.F((1/2)*c+-1//4*pi+(1/2)*d*x,2*b*(a+b)^-1)+32//105*a*d^-1*((a+b)^-1*(a+b*sin(c+d*x)))^-1//2*(a+b*sin(c+d*x))^(1/2)*(11*a^2+13*b^2)*Elliptic.E((1/2)*c+-1//4*pi+(1/2)*d*x,2*b*(a+b)^-1))
@test integrate((a+b*sin(c+d*x))^5//2, x) == :(-2//5*b*d^-1*(a+b*sin(c+d*x))^3//2*cos(c+d*x)+-16//15*a*b*d^-1*(a+b*sin(c+d*x))^(1/2)*cos(c+d*x)+1//15*d^-1*((a+b)^-1*(a+b*sin(c+d*x)))^-1//2*(a+b*sin(c+d*x))^(1/2)*(18*b^2+46*a^2)*Elliptic.E((1/2)*c+-1//4*pi+(1/2)*d*x,2*b*(a+b)^-1)+-16//15*a*d^-1*((a+b)^-1*(a+b*sin(c+d*x)))^(1/2)*(a+b*sin(c+d*x))^-1//2*(a^2+-1*b^2)*Elliptic.F((1/2)*c+-1//4*pi+(1/2)*d*x,2*b*(a+b)^-1))
@test integrate((a+b*sin(c+d*x))^3//2, x) == :(-2//3*b*d^-1*(a+b*sin(c+d*x))^(1/2)*cos(c+d*x)+-1//3*d^-1*((a+b)^-1*(a+b*sin(c+d*x)))^(1/2)*(a+b*sin(c+d*x))^-1//2*(-2*b^2+2*a^2)*Elliptic.F((1/2)*c+-1//4*pi+(1/2)*d*x,2*b*(a+b)^-1)+8//3*a*d^-1*((a+b)^-1*(a+b*sin(c+d*x)))^-1//2*(a+b*sin(c+d*x))^(1/2)*Elliptic.E((1/2)*c+-1//4*pi+(1/2)*d*x,2*b*(a+b)^-1))
@test integrate((a+b*sin(c+d*x))^(1/2), x) == :(2*d^-1*((a+b)^-1*(a+b*sin(c+d*x)))^-1//2*(a+b*sin(c+d*x))^(1/2)*Elliptic.E((1/2)*c+-1//4*pi+(1/2)*d*x,2*b*(a+b)^-1))
@test integrate((a+b*sin(c+d*x))^-1//2, x) == :(2*d^-1*((a+b)^-1*(a+b*sin(c+d*x)))^(1/2)*(a+b*sin(c+d*x))^-1//2*Elliptic.F((1/2)*c+-1//4*pi+(1/2)*d*x,2*b*(a+b)^-1))
@test integrate((a+b*sin(c+d*x))^-3//2, x) == :(2*b*d^-1*(a+b*sin(c+d*x))^-1//2*(a^2+-1*b^2)^-1*cos(c+d*x)+2*d^-1*((a+b)^-1*(a+b*sin(c+d*x)))^-1//2*(a+b*sin(c+d*x))^(1/2)*(a^2+-1*b^2)^-1*Elliptic.E((1/2)*c+-1//4*pi+(1/2)*d*x,2*b*(a+b)^-1))
@test integrate((a+b*sin(c+d*x))^-5//2, x) == :(-2*d^-1*((a+b)^-1*(a+b*sin(c+d*x)))^(1/2)*(a+b*sin(c+d*x))^-1//2*(-3*b^2+3*a^2)^-1*Elliptic.F((1/2)*c+-1//4*pi+(1/2)*d*x,2*b*(a+b)^-1)+2*b*d^-1*(a+b*sin(c+d*x))^-3//2*(-3*b^2+3*a^2)^-1*cos(c+d*x)+8//3*a*b*d^-1*(a+b*sin(c+d*x))^-1//2*(a^2+-1*b^2)^-2*cos(c+d*x)+8//3*a*d^-1*((a+b)^-1*(a+b*sin(c+d*x)))^-1//2*(a+b*sin(c+d*x))^(1/2)*(a^2+-1*b^2)^-2*Elliptic.E((1/2)*c+-1//4*pi+(1/2)*d*x,2*b*(a+b)^-1))
@test integrate((a+b*sin(c+d*x))^-7//2, x) == :(2*b*d^-1*(a+b*sin(c+d*x))^-5//2*(-5*b^2+5*a^2)^-1*cos(c+d*x)+-16//15*a*d^-1*((a+b)^-1*(a+b*sin(c+d*x)))^(1/2)*(a+b*sin(c+d*x))^-1//2*(a^2+-1*b^2)^-2*Elliptic.F((1/2)*c+-1//4*pi+(1/2)*d*x,2*b*(a+b)^-1)+1//15*d^-1*((a+b)^-1*(a+b*sin(c+d*x)))^-1//2*(a+b*sin(c+d*x))^(1/2)*(a^2+-1*b^2)^-3*(18*b^2+46*a^2)*Elliptic.E((1/2)*c+-1//4*pi+(1/2)*d*x,2*b*(a+b)^-1)+2//15*b*d^-1*(a+b*sin(c+d*x))^-1//2*(a^2+-1*b^2)^-3*(9*b^2+23*a^2)*cos(c+d*x)+16//15*a*b*d^-1*(a+b*sin(c+d*x))^-3//2*(a^2+-1*b^2)^-2*cos(c+d*x))

max_stars_repo_path: test/4 Trig functions/4.1 Sine/4.1.1.1 (a+b sin)^n.jl


text:
module IntervalSpecialFunctions

using IntervalArithmetic
using SpecialFunctions
using IntervalRootFinding

import IntervalArithmetic: @round, big53

import SpecialFunctions:
    erf, erfc, erfinv, erfcinv

export
    # functions of 1 arg exported by SpecialFunctions and ArbNumerics
    airyai, airyaiprime, airybi, airybiprime,
    besselj0, besselj1, bessely0, bessely1,
    erf, erfc, erfinv, erfcinv,
    gamma, lgamma, digamma, trigamma,
    eta, zeta,

    # functions of 1 arg exported by both with distinct names
    cosint, ci, sinint, si, ellipk, elliptic_k, ellipe, elliptic_e,

    # functions of 1 arg exported by ArbNumerics only
    chi, shi, ei, xi,  lambertw,

    # functions of a 2 args exported by ArbNumerics only
    elliptic_pi, elliptic_pi2, # the 2 arg versions
    polylog,
    hypergeometric_0F1, regular_hypergeometric_0F1,

    # functions of a 3 args exported by ArbNumerics only
    elliptic_pi, elliptic_pi2, # the 3 arg versions
    hypergeometric_1F1, regular_hypergeometric_1F1,

    # functions of a 4 args exported by ArbNumerics only
    hypergeometric_2F1, regular_hypergeometric_2F1

# Cannot make erfinv and erfcinv work
include("erf.jl")

include("arb.jl")

end # module


text:
"""
    simulate_forward(eventHistory, maxT, startT, bg, kappa, kern)

Simulate a Hawkes process given we have already observed some events. This can be used to forecast future events.

# Arguments
* `eventHistory` The events that have ocured so far.
* `maxT` Maximum time value to simulate forward to.
* `startT` Starting time from which to simulate forward.
* `bg` Background value
* `kappa` Kappa value
* `kern` Kernel function or distribution

"""
function simulate_forward(eventHistory, maxT::Number, startT::Number, bg::Number, kappa::Number, kern::Function)
    bgF(x) = intensity(collect(x .+ startT), eventHistory, bg, kappa, kern)
    events = simulate(bgF, kappa, kern, maxT)
    events .+ startT
end

function simulate_forward(eventHistory, maxT::Number, startT::Number, bg::Number, kappa::Number, kern::Distribution)
    kernD, kern_cdf = generate_kernel_functions(kern)
    simulate_forward(eventHistory, maxT, startT, bg, kappa, kernD)
end


text:
using BridgeDiffEq, StaticArrays
using Test

α=1
β=1
u0=1/2
f(u,p,t) = α*u
g(u,p,t) = β*u
dt = 1//2^(4)
tspan = (0.0,1.0)

prob = ODEProblem(f,u0,(0.0,1.0))
sol = solve(prob,BridgeR3(),dt=dt)
sol = solve(prob,BridgeBS3(),dt=dt)

prob = SDEProblem(f,g,u0,(0.0,1.0))
@time sol = solve(prob,BridgeEuler(),dt=dt)
sol = solve(prob,BridgeHeun(),dt=dt)
sol = solve(prob,BridgeSRK(),dt=dt)

u0 = @SVector [2.0,3.0]
prob = ODEProblem(f,u0,(0.0,1.0))
sol = solve(prob,BridgeR3(),dt=dt)
sol = solve(prob,BridgeBS3(),dt=dt)

prob = SDEProblem(f,g,u0,(0.0,1.0))
sol = solve(prob,BridgeEuler(),dt=dt)
sol = solve(prob,BridgeHeun(),dt=dt)


text:
module Tmp
using Random, BenchmarkTools, Test, InteractiveUtils

include("src/params.jl")
include("src/functions.jl")
import .Params
import .Func

function flipcoin(randval::Float64)
    if randval<0.5
        return "H"
    else
        return "T"
    end
end
@show code_typed(flipcoin(rand()))

end # module


text:
using UUIDs
import Pkg.TOML
using AutoHashEquals

@auto_hash_equals struct RegistryData
    name::String
    uuid::UUID
    repo::Union{String, Nothing}
    description::Union{String, Nothing}
    packages::Dict{String, Dict{String, Any}}
    extra::Dict{String, Any}
end

function RegistryData(
    name::AbstractString,
    uuid::Union{UUID, AbstractString};
    repo::Union{AbstractString, Nothing}=nothing,
    description::Union{AbstractString, Nothing}=nothing,
    packages::Dict=Dict(),
    extras::Dict=Dict(),
)
    RegistryData(name, UUID(uuid), repo, description, packages, extras)
end

function parse_registry(io::IO)
    data = TOML.parse(io)

    name = pop!(data, "name")
    uuid = pop!(data, "uuid")
    repo = pop!(data, "repo", nothing)
    description = pop!(data, "description", nothing)
    packages = pop!(data, "packages", Dict())

    RegistryData(name, UUID(uuid), repo, description, packages, data)
end

parse_registry(str::AbstractString) = open(parse_registry, str)

function TOML.print(io::IO, reg::RegistryData)
    println(io, "name = ", repr(reg.name))
    println(io, "uuid = ", repr(string(reg.uuid)))

    if reg.repo !== nothing
        println(io, "repo = ", repr(reg.repo))
    end

    if reg.description !== nothing
        # print long-form string if there are multiple lines
        if '\n' in reg.description
            print(io, """\n
                description = \"\"\"
                $(reg.description)\"\"\"
                """
            )
        else
            println(io, "description = ", repr(reg.description))
        end
    end

    for (k, v) in pairs(reg.extra)
        TOML.print(io, Dict(k => v), sorted=true)
    end

    println(io, "\n[packages]")
    for (uuid, data) in sort!(collect(reg.packages), by=first)
        print(io, uuid, " = { ")
        print(io, "name = ", repr(data["name"]))
        print(io, ", path = ", repr(data["path"]))
        println(io, " }")
    end

    nothing
end

max_stars_repo_name: sisl/registratura

max_issues_repo_name: sisl/registratura

max_forks_repo_name: sisl/registratura


text:
"""
    PeriodicTransform(f)

Makes a kernel periodic by mapping a scalar input onto the unit circle. Samples from a GP
with a kernel with this transformation applied will produce samples with frequency `f`.
"""
struct PeriodicTransform{Tf<:AbstractVector{<:Real}} <: Transform
    f::Tf
end

@functor PeriodicTransform

PeriodicTransform(f::Real) = PeriodicTransform([f])

dim(t::PeriodicTransform) = 2

(t::PeriodicTransform)(x::Real) = [sinpi(2 * first(t.f) * x), cospi(2 * first(t.f) * x)]

function _map(t::PeriodicTransform, x::AbstractVector{<:Real})
    return RowVecs(hcat(sinpi.((2 * first(t.f)) .* x), cospi.((2 * first(t.f)) .* x)))
end

function Base.isequal(t1::PeriodicTransform, t2::PeriodicTransform)
    return isequal(first(t1.f), first(t2.f))
end

function Base.show(io::IO, t::PeriodicTransform)
    print(io, "Periodic Transform with frequency $(first(t.f))")
end


text:
import QMTK: SubSites, Bit, Spin, Half, data, sitetype
using Compat.Test

@testset "Construction" begin

    @test SubSites(Bit, 1, 0, 1).data == (1, 0, 1)
    @test SubSites(Bit, (1, 0, 1)).data == (1, 0, 1)
    @test eltype(SubSites(Bit, Float32, 1, 0, 1)) == Float32
    @test eltype(SubSites(Bit, Float32, (1, 0, 1))) == Float32

    @test sitetype(SubSites{Bit, Int, 2}) == Bit

    self = SubSites(Bit, 1, 0, 1)
    @test data(self) === self.data
end

@testset "Conversions" begin

    for TYPE in (UInt8, UInt16, UInt32, UInt64, UInt128,
                 Int8, Int16, Int32, Int64, Int128)

        @test TYPE(SubSites(Bit, 1, 0, 0)) === TYPE(0b001)
        @test TYPE(SubSites(Bit, 1, 1, 0)) === TYPE(0b011)

        @test TYPE(SubSites(Spin, 1, -1, -1)) === TYPE(0b001)
        @test TYPE(SubSites(Spin, 1,  1, -1)) === TYPE(0b011)

        @test TYPE(SubSites(Half, 0.5, -0.5, -0.5)) === TYPE(0b001)
        @test TYPE(SubSites(Half, 0.5,  0.5, -0.5)) === TYPE(0b011)

    end

    @test convert(SubSites{Bit, Int, 3}, SubSites(Spin, 1, -1, -1)) == SubSites(Bit, 1, 0, 0)
    @test convert(SubSites{Spin, Int, 3}, SubSites(Bit, 1, 0, 0)) == SubSites(Spin, 1, -1, -1)
    @test convert(SubSites{Half, Float64, 3}, SubSites(Bit, 1, 0, 1)) == SubSites(Half, 0.5, -0.5, 0.5)
end


text:
# This script actually solves the CoeurdacierReyWinant model with a risk-adjusted linearization
# and times the methods, if desired
using BenchmarkTools, RiskAdjustedLinearizations, Test, JLD2
include("bansal_yaron_2004.jl")
include("het_risk_aversion.jl")

# What to do?
solve_rep_agent = false
solve_het_agent = true
check_mat       = false

# Set up
if solve_rep_agent
    m_rep = BansalYaron2004()
    m = bansal_yaron_2004(m_rep; sparse_arrays = true, sparse_jacobian = [:μ, :ξ])
    z0 = copy(m.z)
    y0 = copy(m.y)
    Ψ0 = copy(m.Ψ)

    # Solve!
    solve!(m, m.z, m.y; algorithm = :deterministic)
    solve!(m; algorithm = :relaxation)

    # Now solve with a higher risk aversion, using previous steady state as a guess
    # Note the parameters used in the equations for the RiskAdjustedLinearization
    # are from the field `m_rep.p`, which is a LabelledArray and thus allows us to
    # mutate the parameters without re-compiling the functions contained in `m`.
    m_rep.p[:γ] = 9.5
    solve!(m; algorithm = :relaxation)
    m_rep.p[:γ] = 9. # restore the original parameter so, for example, the check_mat test passes
end

if solve_het_agent
    m_het = HetRiskAversion()
    m = het_risk_aversion(m_het; sparse_arrays = false) # , sparse_jacobian = [:μ, :ξ])
    z0 = copy(m.z)
    y0 = copy(m.y)
    Ψ0 = copy(m.Ψ)

    # Solve!
    # solve!(m; algorithm = :homotopy, step = .1)

#=    # Now solve with a higher risk aversion, using previous steady state as a guess
    # Note the parameters used in the equations for the RiskAdjustedLinearization
    # are from the field `m_rep.p`, which is a LabelledArray and thus allows us to
    # mutate the parameters without re-compiling the functions contained in `m`.
    m_het.p[:γ] = 9.5
    solve!(m; algorithm = :relaxation)
    m_het.p[:γ] = 9. # restore the original parameter so, for example, the check_mat test passes=#
end

if check_mat
    using MAT
    matout = matread("dynare_ss.mat")
    @unpack J, S = m_rep

    solve!(m, z0, y0; algorithm = :deterministic)
    @test log(matout["Q"]) ≈ m.y[J[:q]] atol=1e-6
    @test log(matout["Omega"]) ≈ m.y[J[:ω]] atol=1e-6
    @test log(matout["V"]) ≈ m.y[J[:v]] atol=1e-6
    @test log(matout["CE"]) ≈ m.y[J[:ce]] atol=1e-6
    @test log(matout["PQ"]) ≈ m.y[J[:pq1]] atol=1e-6
    @test log(matout["DQ"]) ≈ m.y[J[:dq1]] atol=1e-6
    @test log(matout["POmega"]) ≈ m.y[J[:pω1]] atol=1e-6
    @test log(matout["DOmega"]) ≈ m.y[J[:dω0]] atol=1e-6
    @test log(matout["X"]) ≈ m.z[S[:x]] atol=1e-6
    @test log(matout["Y"]) ≈ m.z[S[:yy]] atol=1e-6
    @test matout["SigSq"] ≈ m.z[S[:σ²_y]] atol=1e-6
end


text:
function _record_pointers(o::IO)

    n = o.size ÷ 80                             # number of records in IO (36 records/block)

    r = [(i-1)*80 for i=1:n]                    # start-of-record pointers (80 bytes/record)

    return r

end

function _block_pointers(o::IO)

    n = o.size ÷ 2880                           # number of blocks in IO

    b = [(i-1)*2880 for i=1:n]                  # start-of-block pointers (2880 bytes/block)

    return b

end

function _header_pointers(o::IO)

    b = _block_pointers(o::IO)                  # b: start-of-block pointers

    h::Array{Int,1} = []                        # h: init start-of-header pointers

    for i ∈ Base.eachindex(b)                   # i: start-of-block pointer
        Base.seek(o,b[i])
        key = String(Base.read(o,8))
        key ∈ ["SIMPLE  ", "XTENSION"] ? Base.push!(h,b[i]) : false
    end

    return h                                    # return start-of-header pointers

end

function _hdu_pointers(o::IO)

    h = _header_pointers(o::IO)                 # h: start-of-header pointers

    return h                                    # return start-of-HDU (= start-of-header) pointers

end

function _data_pointers(o::IO)

    b = _block_pointers(o::IO)                  # b: start-of-block pointers

    d::Array{Int,1} = []                        # h: init start-of-data pointers

    for i ∈ eachindex(b)                        # i: start-of-block pointer (36 records/block = 2880 bytes)
        Base.seek(o,b[i])
        [(key = String(Base.read(o,8)); key == "END     " ? Base.push!(d,b[i]+2880) : Base.skip(o,72)) for j=0:35]
    end

    return d                                    # return start-of-data pointers

end


text:
"""
	minimumN(kout, kin, shape_function; tol = 1e-9, N_points = 10_000,
        N_start = 400, N_min = 100, N_max = 1_000) -> N, err

Return the minimum `N` necessary (i.e. `2N` nodes) to achieve error of at most
`tol` in the electric field for a `ShapeParams` inclusion created by
`shape_function(N)` which is filled with material of wavenumber `kin` and
surrounded by free space with wavenumber `k0`.
Error is calculated on `N_points` points on the scattering disk (`s.R`), by
assuming a fictitious line source and comparing its field to that produced
by the resulting potential densities.

Since the error scales with \$N^{-3}\$ for moderate wavelengths and errors, we
estimate `N` using the error of `N_start`, then binary search based on that
guess between `N_min` and `N_max`.
"""
function minimumN(kout, kin, shape_function; tol = 1e-9, N_points = 10_000, N_start = 400, N_min = 100, N_max = 1_000)
    # we know that in general err(N) = a/N^3
    #estimate a = err1*N1^3, then N2 = (a/err2)^(1/3) = N1*(err1/err2)^(1/3)
    #assumes unit line source at (0,0), unit plane wave outside

    s = shape_function(N_max) #just for radius
    err_points = [s.R*f(i*2*pi/N_points) for i=0:(N_points-1), f in (cos,sin)]
    E_ana = (0.25im*besselh(0, 1, kout*s.R))*ones(Complex{Float64}, N_points)
    E_comp = Array{Complex{Float64}}(undef, N_points)

    N = N_start
    err_start = minimumN_helper(N, kout, kin, shape_function, err_points, E_comp, E_ana)
    #if starting point was too good, it can be our new max
    err_start < tol && (N_max = min(N_start,N_max))

    N = round(Int, N_start*(err_start/tol)^(1/3))
    N = max(N, N_min)

    err = minimumN_helper(N, kout, kin, shape_function, err_points, E_comp, E_ana)

    if err > tol
        if N > N_max
            @warn("minimumN: Cannot enter binary search mode.")
        else
            @warn("minimumN: Entering binary search mode in [$N,$N_max].")
            N, err = binarySearch(N_ -> -minimumN_helper(N_, kout, kin,
                shape_function, err_points, E_comp, E_ana), -tol, N, N_max)
        end
        err > tol && @warn("minimumN: Failed to find err(N) < tol.")
    elseif err < tol
		if N == N_min
			return N,err
		end
        @warn("minimumN: Entering binary search in [$N_min,$N].")
        N, err = binarySearch(N_ -> -minimumN_helper(N_, kout, kin,
                    shape_function, err_points, E_comp, E_ana), -tol, N_min, N)
    end
    return N, err
end

function minimumN_helper(N, kout, kin, shape_function, err_points, E_comp, E_ana)
    #helper function for computing error for given N
    s = shape_function(N)
    mu_sigma = solvePotential_forError(kin, kout, s, [0.0 0.0], [1], 0.0)

    E_comp[:] = scatteredfield(mu_sigma, kout, s, err_points)
    this_err = norm(E_ana-E_comp)/norm(E_ana)
end

#TODO: similar minimumP for circles, comparing perhaps to collection of line sources

"""
	minimumP(k0, kin, s::ShapeParams; tol = 1e-9, N_points = 10_000, P_min = 1,
        P_max = 60, dist = 2) -> P, errP

Return the minimum `P` necessary to achieve error of at most `tol` in the
electric field, when compared to that obtained with `2N` discretization, for a
`ShapeParams` inclusion filled with material of wavenumber `kin` and surrounded
by free space with wavenumber `k0`.
Error is calculated on `N_points` points on a disk of radius `dist*s.R`.

Uses binary search between `P_min` and `P_max`.
"""
function minimumP(k0, kin, s::ShapeParams; tol = 1e-9, N_points = 10_000, P_min = 1, P_max = 60, dist = 2.0)
	err_points = [dist*s.R*f(i*2*pi/N_points) for i=0:(N_points-1), f in (cos,sin)]
	@warn("""minimumP: Calculating at $dist*R, thus implicitly assuming all scatterers have same radius and we do not care about closer implications.""")

    #compute direct solution for comparison
    inner = get_potentialPW(k0, kin, s, 0.0)
    E_quadrature = scatteredfield(inner, k0, s, err_points)
    E_multipole = Array{Complex{Float64}}(undef, N_points)

    P, errP = binarySearch(P_ -> -minimumP_helper(k0, kin, s, P_,
		N_points, err_points, E_quadrature, E_multipole), -tol, P_min, P_max)
end

function minimumP_helper(k0, kin, s, P, N_points, err_points, E_quadrature, E_multipole)
    #helper function for computing error for given P
    E_multipole[:] .= 0.0
	sp = ScatteringProblem([s],[1],[0.0 0.0],[0.0])
    beta_p,inn = solve_particle_scattering(k0, kin, P, sp, PlaneWave(0.0),
		verbose = false)
    scattered_field_multipole!(E_multipole, k0, beta_p, P, [0.0 0.0], [1],
		err_points, 1:N_points)
    this_err = norm(E_quadrature - E_multipole)/norm(E_quadrature)
end