text:
import Base: summary, show

abstract type AbstractNsolt{T,D} <: PolyphaseFB{T,D} end
abstract type Rnsolt{T,D} <: AbstractNsolt{T,D} end

Rnsolt(df::Integer, ppo::Integer, nchs; kwargs...) = Rnsolt(Float64, df, ppo, nchs; kwargs...)
Rnsolt(::Type{T}, df::Integer, ppo::Integer, nchs::Integer; kwargs...) where {T} = Rnsolt(T, df, ppo, (cld(nchs,2),fld(nchs,2)); kwargs...)
function Rnsolt(::Type{T}, df::Integer, ppo::Integer, nchs::Tuple{Int,Int}; dims::Integer=1, kwargs...) where {T}
    Rnsolt(T, (fill(df,dims)...,), (fill(ppo,dims)...,), nchs; kwargs...)
end

Rnsolt(df::NTuple{D,Int}, ppo::NTuple{D,Int}, nchs::Union{Tuple{Int,Int}, Integer}; kwargs...) where {D} = Rnsolt(Float64, df, ppo, nchs; kwargs...)

Rnsolt(::Type{T}, df::NTuple{D,Int}, ppo::NTuple{D,Int}, nchs::Integer; kwargs...) where {D,T} = Rnsolt(T, df, ppo, (cld(nchs,2), fld(nchs,2)); kwargs...)

function Rnsolt(::Type{T}, df::NTuple{D,Int}, ppo::NTuple{D,Int}, nchs::Tuple{Int,Int}; kwargs...) where {T,D}
    return if nchs[1] == nchs[2]
        RnsoltTypeI(T, df, ppo, nchs; kwargs...)
    else
        RnsoltTypeII(T, df, ppo, nchs; kwargs...)
    end
end

struct RnsoltTypeI{T,D} <: Rnsolt{T,D}
    decimationFactor::NTuple{D, Int}
    nStages::NTuple{D, Int}
    nChannels::Tuple{Int, Int}

    CJ::AbstractMatrix # reverse of permutated DCT

    W0::AbstractMatrix{T} # parameter matrix of initial matrix
    U0::AbstractMatrix{T}

    Udks::NTuple{D, Vector{AbstractMatrix{T}}} # parameter matrices of propagation matrices

    perm::NTuple{D, Int} # permutation of dimensions

    function RnsoltTypeI(::Type{T}, df::NTuple{D, Int}, ppo::NTuple{D, Int}, nchs::Tuple{Int, Int}; perm = (collect(1:D)...,), kwargs...) where {D,T}
        @assert isperm(perm) "invalid permutations"
        @assert (nchs[1] == nchs[2]) "channel size mismatch!"
        @assert (cld(prod(df),2) <= nchs[1] <= sum(nchs) - fld(prod(df),2)) && (fld(prod(df),2) <= nchs[2] <= sum(nchs) - cld(prod(df),2)) "invalid number of channels"
        # CJ = reverse(permdctmtx(T, df...), dims=2) |> Matrix
        CJ = permdctmtx(T, df...)

        W0 = Matrix{T}(I, nchs[1], nchs[1])
        U0 = Matrix{T}(I, nchs[2], nchs[2])
        Udks = ([
            [ (iseven(k) ? 1 : -1) * Matrix{T}(I, nchs[2], nchs[2]) for k in 1:ppo[d] ]
        for d in 1:D ]...,)

        new{T,D}(df, ppo, nchs, CJ, W0, U0, Udks, perm)
    end
end

RnsoltTypeI(df::NTuple{D,Int}, ppo::NTuple{D,Int}, nchs; kwargs...) where {D} = RnsoltTypeI(Float64, df, ppo, nchs; kwargs...)
RnsoltTypeI(::Type{T}, df::NTuple{D,Int}, ppo::NTuple{D,Int}, nchs::Integer; kwargs...) where {T,D} = RnsoltTypeI(T, df, ppo, (cld(nchs,2), fld(nchs,2)); kwargs...)

struct RnsoltTypeII{T,D} <: Rnsolt{T,D}
    decimationFactor::NTuple{D, Int}
    nStages::NTuple{D, Int}
    nChannels::Tuple{Int, Int}

    CJ::AbstractMatrix # reverse of permutated DCT

    # matricesW0U0Odd::AbstractArray{}

    W0::AbstractMatrix{T} # parameter matrix of initial matrix
    U0::AbstractMatrix{T} # parameter matrix of initial matrix

    Wdks::NTuple{D, Vector{AbstractMatrix{T}}} # parameter matrices of propagation matrices
    Udks::NTuple{D, Vector{AbstractMatrix{T}}} # parameter matrices of propagation matrices

    perm::NTuple{D, Int}

    function RnsoltTypeII(::Type{T}, df::NTuple{D, Int}, ppo::NTuple{D, Int}, nchs::Tuple{Int, Int}; perm = (collect(1:D)...,), kwargs...) where {D,T}
        @assert isperm(perm) "invalid permutations"
        @assert (cld(prod(df),2) <= nchs[1] <= sum(nchs) - fld(prod(df),2)) && (fld(prod(df),2) <= nchs[2] <= sum(nchs) - cld(prod(df),2)) "invalid number of channels"
        # CJ = reverse(permdctmtx(T, df...), dims=2) |> Matrix
        CJ = permdctmtx(T, df...)

        W0 = Matrix{T}(I, nchs[1], nchs[1])
        U0 = Matrix{T}(I, nchs[2], nchs[2])

        @assert all(iseven.(ppo)) "polyphase order of each dimension must be odd"
        nStages = fld.(ppo, 2)

        signW, signU = if nchs[1] < nchs[2]; (-1, 1) else (1, -1) end
        Wdks = ([
            [ signW * Matrix{T}(I, nchs[1], nchs[1]) for k in 1:nStages[d] ]
        for d in 1:D ]...,)

        Udks = ([
            [ signU * Matrix{T}(I, nchs[2], nchs[2]) for k in 1:nStages[d] ]
        for d in 1:D ]...,)

        new{T,D}(df, nStages, nchs, CJ, W0, U0, Wdks, Udks, perm)
    end
end

RnsoltTypeII(df::NTuple{D,Int}, ppo::NTuple{D,Int}, nchs; kwargs...) where {D} = RnsoltTypeII(Float64, df, ppo, nchs; kwargs...)
RnsoltTypeII(::Type{T}, df::NTuple{D,Int}, ppo::NTuple{D,Int}, nchs::Integer; kwargs...) where {T,D} = RnsoltTypeII(T, df, ppo, (cld(nchs,2), fld(nchs,2)); kwargs...)

abstract type Cnsolt{T,D} <: AbstractNsolt{T,D} end

Cnsolt(df::NTuple{D,Int}, ppo::NTuple{D,Int}, nchs::Integer; kwargs...) where {D} = Cnsolt(Float64, df, ppo, nchs; kwargs...)
Cnsolt(df::Integer, ppo::Integer, nchs::Integer; kwargs...) = Cnsolt(Float64, df, ppo, nchs; kwargs...)

function Cnsolt(::Type{T}, df::Integer, ppo::Integer, nchs::Integer; dims=1, kwargs...) where {T}
    Cnsolt(T, (fill(df, dims)...,), (fill(ppo, dims)...,), nchs; kwargs...)
end

function Cnsolt(::Type{T}, df::NTuple{D,Int}, ppo::NTuple{D,Int}, nchs::Integer; kwargs...) where {T,D}
    return if iseven(nchs)
        CnsoltTypeI(T, df, ppo, nchs; kwargs...)
    else
        CnsoltTypeII(T, df, ppo, nchs; kwargs...)
    end
end

struct CnsoltTypeI{T,D} <: Cnsolt{T,D}
    decimationFactor::NTuple{D,Int}
    nStages::NTuple{D,Int}
    nChannels::Integer

    FJ::AbstractMatrix

    V0::AbstractMatrix{T}

    Wdks::NTuple{D, Vector{AbstractMatrix{T}}}
    Udks::NTuple{D, Vector{AbstractMatrix{T}}}
    θdks::NTuple{D, Vector{AbstractVector{T}}}

    Φ::Diagonal

    perm::NTuple{D, Int}

    function CnsoltTypeI(::Type{T}, df::NTuple{D,Int}, ppo::NTuple{D,Int}, nchs::Integer; perm = (collect(1:D)...,), kwargs...) where {T,D}
        @assert isperm(perm) "invalid permutations"
        @assert (prod(df) <= sum(nchs)) "number of channels must be greater or equal to decimation factor"
        FJ = cdftmtx(T, df...)

        V0 = Matrix{T}(I, nchs, nchs)

        fP = fld(nchs, 2)
        Wdks = ([
            [ Matrix{T}(I, fP, fP) for k in 1:ppo[d] ]
        for d in 1:D ]...,)

        Udks = ([
            [ -Matrix{T}(I, fP, fP) for k in 1:ppo[d] ]
        for d in 1:D ]...,)

        θdks = ([
            [ zeros(T, fld(nchs, 4)) for k in 1:ppo[d] ]
        for d in 1:D ]...,)

        Φ = Diagonal(cis.(zeros(nchs))) |> complex

        new{T,D}(df, ppo, nchs, FJ, V0, Wdks, Udks, θdks, Φ, perm)
    end
end

CnsoltTypeI(df::NTuple{D,Int}, ppo::NTuple{D,Int}, nchs::Integer; kwargs...) where {D} = CnsoltTypeI(Float64, df, ppo, nchs; kwargs...)

struct CnsoltTypeII{T,D} <: Cnsolt{T,D}
    decimationFactor::NTuple{D,Int}
    nStages::NTuple{D,Int}
    nChannels::Integer

    FJ::AbstractMatrix

    V0::AbstractMatrix{T}

    Wdks::NTuple{D, Vector{AbstractMatrix{T}}}
    Udks::NTuple{D, Vector{AbstractMatrix{T}}}
    θ1dks::NTuple{D, Vector{AbstractVector{T}}}
    Ŵdks::NTuple{D, Vector{AbstractMatrix{T}}}
    Ûdks::NTuple{D, Vector{AbstractMatrix{T}}}
    θ2dks::NTuple{D, Vector{AbstractVector{T}}}

    Φ::Diagonal

    perm::NTuple{D, Int}

    function CnsoltTypeII(::Type{T}, df::NTuple{D,Int}, ppo::NTuple{D,Int}, nchs::Integer; perm = (collect(1:D)...,), kwargs...) where {T,D}
        @assert isperm(perm) "invalid permutations"
        @assert (prod(df) <= sum(nchs)) "number of channels must be greater or equal to decimation factor"
        @assert all(iseven.(ppo)) "polyphase order of each dimension must be odd"
        nStages = fld.(ppo, 2)
        FJ = cdftmtx(T, df...)

        V0 = Matrix{T}(I, nchs, nchs)

        fP, cP = fld(nchs,2), cld(nchs,2)
        Wdks = ([
            [ Matrix{T}(I, fP, fP) for k in 1:nStages[d] ]
        for d in 1:D ]...,)

        Udks = ([
            [ -Matrix{T}(I, fP, fP) for k in 1:nStages[d] ]
        for d in 1:D ]...,)

        θ1dks = ([
            [ zeros(T, fld(nchs, 4)) for k in 1:nStages[d] ]
        for d in 1:D ]...,)

        Ŵdks = ([
            [ Matrix{T}(I, cP, cP) for k in 1:nStages[d] ]
        for d in 1:D ]...,)

        signÛ = Diagonal([fill(-1,cP-1)..., 1])
        Ûdks = ([
            [ signÛ * Matrix{T}(I, cP, cP) for k in 1:nStages[d] ]
        for d in 1:D ]...,)

        θ2dks = ([
            [ zeros(T, fld(nchs, 4)) for k in 1:nStages[d] ]
        for d in 1:D ]...,)

        Φ = Diagonal(cis.(zeros(nchs))) |> complex

        new{T,D}(df, nStages, nchs, FJ, V0, Wdks, Udks, θ1dks, Ŵdks, Ûdks, θ2dks, Φ, perm)
    end
end

CnsoltTypeII(df::NTuple{D,Int}, ppo::NTuple{D,Int}, nchs::Integer; kwargs...) where {D} = CnsoltTypeII(Float64, df, ppo, nchs; kwargs...)

promote_rule(::Type{RnsoltTypeI{TA,D}}, ::Type{RnsoltTypeI{TB,D}}) where {D,TA,TB} = RnsoltTypeI{promote_rule(TA,TB),D}

promote_rule(::Type{RnsoltTypeII{TA,D}}, ::Type{RnsoltTypeII{TB,D}}) where {D,TA,TB} = RnsoltTypeII{promote_rule(TA,TB),D}

promote_rule(::Type{CnsoltTypeI{TA,D}}, ::Type{CnsoltTypeI{TB,D}}) where {D,TA,TB} = CnsoltTypeI{promote_rule(TA,TB),D}

promote_rule(::Type{CnsoltTypeII{TA,D}}, ::Type{CnsoltTypeII{TB,D}}) where {D,TA,TB} = CnsoltTypeII{promote_type(TA,TB),D}

similar(nsolt::Rnsolt{T,DS}, element_type::Type=T, df::NTuple{DD}=decimations(nsolt), ord::NTuple{DD}=orders(nsolt), nch=nchannels_sa(nsolt)) where {T,DS,DD} = Rnsolt(element_type, df, ord, nch)

similar(nsolt::Cnsolt{T,DS}, element_type::Type=T, df::NTuple{DD}=decimations(nsolt), ord::NTuple{DD}=orders(nsolt), nch=nchannels(nsolt)) where {T,DS,DD} = Cnsolt(element_type, df, ord, nch)

nchannels_sa(nsolt::AbstractNsolt) = nchannels(nsolt)
nchannels_sa(nsolt::Rnsolt) = nsolt.nChannels

orders(fb::RnsoltTypeI) = fb.nStages
orders(fb::CnsoltTypeI) = fb.nStages
orders(fb::RnsoltTypeII) = 2 .* fb.nStages
orders(fb::CnsoltTypeII) = 2 .* fb.nStages

nstages(fb::AbstractNsolt) = fb.nStages

istype1(::Type{NS}) where {NS<:RnsoltTypeI} = true
istype1(::Type{NS}) where {NS<:CnsoltTypeI} = true
istype1(::Type{NS}) where {NS<:RnsoltTypeII} = false
istype1(::Type{NS}) where {NS<:CnsoltTypeII} = false

istype1(::T) where {T<:AbstractNsolt} = istype1(T)

istype2(::Type{T}) where {T<:AbstractNsolt} = !istype1(T)
istype2(::T) where {T<:AbstractNsolt} = !istype1(T)

function supertype_nsolt(::Type{NS}) where {NS<:AbstractNsolt}
    RNS = typeintersect(promote_type(NS, Rnsolt), AbstractNsolt)
    CNS = typeintersect(promote_type(NS, Cnsolt), AbstractNsolt)
    typeintersect(RNS, CNS)
end

function copy_params!(dst::NS, src::NS) where {NS<:RnsoltTypeI}
    mycp(a,b) = foreach(a, b) do ad, bd
        foreach((adk, bdk) -> adk .= bdk, ad, bd)
    end

    dst.CJ .= src.CJ

    dst.W0 .= src.W0
    dst.U0 .= src.U0

    mycp(dst.Udks, src.Udks)

    return dst
end

function copy_params!(dst::NS, src::NS) where {NS<:RnsoltTypeII}
    mycp(a,b) = foreach(a, b) do ad, bd
        foreach((adk, bdk) -> adk .= bdk, ad, bd)
    end

    dst.CJ .= src.CJ

    dst.W0 .= src.W0
    dst.U0 .= src.U0

    mycp(dst.Wdks, src.Wdks)
    mycp(dst.Udks, src.Udks)

    return dst
end

function copy_params!(dst::NS, src::NS) where {NS<:CnsoltTypeI}
    mycp(a,b) = foreach(a, b) do ad, bd
        foreach((adk, bdk) -> adk .= bdk, ad, bd)
    end

    dst.FJ .= src.FJ

    dst.V0 .= src.V0

    mycp(dst.Wdks, src.Wdks)
    mycp(dst.Udks, src.Udks)
    mycp(dst.θdks, src.θdks)

    dst.Φ .= src.Φ

    return dst
end

function copy_params!(dst::NS, src::NS) where {NS<:CnsoltTypeII}
    mycp(a,b) = foreach(a, b) do ad, bd
        foreach((adk, bdk) -> adk .= bdk, ad, bd)
    end

    dst.FJ .= src.FJ

    dst.W0 .= src.W0
    dst.U0 .= src.U0

    mycp(dst.Wdks, src.Wdks)
    mycp(dst.Udks, src.Udks)
    mycp(dst.θ1dks, src.θ1dks)
    mycp(dst.Ŵdks, src.Ŵdks)
    mycp(dst.Ûdks, src.Ûdks)
    mycp(dst.θ2dks, src.θ2dks)

    dst.Φ .= src.Φ

    return dst
end

Cnsolt(rn::RnsoltTypeI) = CnsoltTypeI(rn)

function CnsoltTypeI(rn::RnsoltTypeI{T}) where {T}
    Pw = rn.nChannels[1]
    Pu = rn.nChannels[2]
    M = prod(rn.decimationFactor)
    cM, fM = cld(M,2), fld(M,2)

    cn = CnsoltTypeI(T, decimations(rn), orders(rn), nchannels(rn))

    # cn.FJ .= diagm( 0 => [ ones(cM); 1im * ones(fM) ] ) * rn.CJ

    cn.V0 .= begin
        pms = [ collect(1:cM)...,
                  collect((1:Pw-cM) .+ M)...,
                  collect((1:fM) .+ cM)...,
                  collect((1:Pu-fM) .+ (Pw+fM))... ]

        pmtx = foldl(enumerate(pms), init=zero(cn.V0)) do mtx, (idx, pmi)
            setindex!(mtx, 1, idx, pmi)
        end

        # cat(rn.W0, rn.U0, dims=[1,2]) * pmtx
        sgnmtx = diagm( 0 => [ ones(cM); 1im * ones(fM) ] )
        C0 = cat(real( sgnmtx * rn.CJ * cn.FJ'), Matrix(I, (size(cn.V0) .- M)...), dims=[1,2])
        cat(rn.W0, rn.U0, dims=[1,2]) * pmtx * C0
    end

    foreach(cn.Wdks) do cWs
        foreach(cWs) do cW
            cW .= Matrix{T}(I, size(cW)...)
        end
    end

    foreach(cn.Udks, rn.Udks) do cUs, rUs
        foreach(cUs, rUs) do cU, rU
            cU .= rU
        end
    end

    foreach(cn.θdks) do cθs
        foreach(cθs) do cθ
            cθ .= zeros(T, size(cθ)...)
        end
    end

    cn.Φ .= Diagonal([ ones(Pw); -1im * ones(Pu) ])

    return cn
end

function show(io::IO, ::MIME"text/plain", nsolt::Rnsolt)
    print(io, "$(nsolt.nChannels)-channels $(typeof(nsolt)) with Decimation factor=$(nsolt.decimationFactor), Polyphase order=$(orders(nsolt))")
end

function show(io::IO, ::MIME"text/plain", nsolt::Cnsolt)
    print(io, "$(nsolt.nChannels)-channels $(typeof(nsolt)) with Decimation factor=$(nsolt.decimationFactor), Polyphase order=$(orders(nsolt))")
end

# permutedims_prop(n::AbstractNsolt, perm::AbstractVector) = permutedims_prop(n, (perm...,))
#
# function permutedims_prop(nsolt::NS, perm::NTuple{D,N}) where {T,D,NS<:AbstractNsolt{T,D},N<:Integer}
#     @assert isperm(perm) "invalid permutation"
#
#     Nsolt = supertype_nsolt(NS)
#     fp(v) = ([ v[perm[idx]] for idx = 1:D ]...,)
#     out = Nsolt(T, fp(decimations(nsolt)), fp(orders(nsolt)), nchannels(nsolt))
#
#     return setPermutatedParams!(out, nsolt, perm)
# end
#
# function setPermutatedParams!(dst::NS, src::NS, perm) where {NS<:RnsoltTypeI}
#     dst.CJ .= copy(src.CJ)
#     dst.U0 .= copy(src.U0)
#     dst.W0 .= copy(src.W0)
#     for (idxdst, idxsrc) in enumerate(perm)
#         dst.Udks[idxdst] .= copy(src.Udks[idxsrc])
#     end
#     dst
# end
#
# function setPermutatedParams!(dst::NS, src::NS, perm) where {NS<:RnsoltTypeII}
#     dst.CJ .= copy(src.CJ)
#     dst.U0 .= copy(src.U0)
#     dst.W0 .= copy(src.W0)
#     for (idxdst, idxsrc) in enumerate(perm)
#         dst.Udks[idxdst] .= copy(src.Udks[idxsrc])
#         dst.Wdks[idxdst] .= copy(src.Wdks[idxsrc])
#     end
#     dst
# end
#
# function setPermutatedParams!(dst::NS, src::NS, perm) where {NS<:CnsoltTypeI}
#     dst.FJ .= copy(src.FJ)
#     dst.V0 .= copy(src.V0)
#     for (idxdst, idxsrc) in enumerate(perm)
#         dst.Udks[idxdst] .= copy(src.Udks[idxsrc])
#         dst.Wdks[idxdst] .= copy(src.Wdks[idxsrc])
#         dst.θdks[idxdst] .= copy(src.θdks[idxsrc])
#     end
#     dst.Φ .= copy(src.Φ)
#     dst
# end
#
# function setPermutatedParams!(dst::NS, src::NS, perm) where {NS<:CnsoltTypeII}
#     dst.FJ .= copy(src.FJ)
#     dst.V0 .= copy(src.V0)
#     for (idxdst, idxsrc) in enumerate(perm)
#         dst.Udks[idxdst] .= copy(src.Udks[idxsrc])
#         dst.Wdks[idxdst] .= copy(src.Wdks[idxsrc])
#         dst.θ1dks[idxdst] .= copy(src.θ1dks[idxsrc])
#         dst.Ûdks[idxdst] .= copy(src.Ûdks[idxsrc])
#         dst.Ŵdks[idxdst] .= copy(src.Ŵdks[idxsrc])
#         dst.θ2dks[idxdst] .= copy(src.θ2dks[idxsrc])
#     end
#     dst.Φ .= copy(src.Φ)
#     dst
# end

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
# Testing functions
function testAllCases()
    @testset "All tests      " begin
        smallTest()
        #mediumTest()
        #bigTest()
    end
end
function smallTest()
    blockLength = 10
    Nblocks = 5
    NbootstrapReplicates = 2
    Ndata = 500
    statisticHandle = mean
    @testset "Small test     " begin
        testAllData(Ndata, statisticHandle, blockLength, Nblocks, NbootstrapReplicates)
    end
end
function mediumTest()
    blockLength = 40
    Nblocks = 10
    NbootstrapReplicates = 400
    Ndata = 5000
    statisticHandle = mean
    @testset "Medium test    " begin
        testAllData(Ndata, statisticHandle, blockLength, Nblocks, NbootstrapReplicates)
    end
end
function bigTest()
    blockLength = 500
    Nblocks = 100
    NbootstrapReplicates = 10000
    Ndata = 50000
    statisticHandle = mean
    @testset "Big test       " begin
        testAllData(Ndata, statisticHandle, blockLength, Nblocks, NbootstrapReplicates)
    end
end
function testAllData(Ndata, statisticHandle, blockLength, Nblocks, NbootstrapReplicates)
    # Testing different data sources
    @testset "Linear integers" begin
        testDataLinearInts = getTestDataLinearInts(Ndata)
        allTestTypes(testDataLinearInts, statisticHandle, blockLength, Nblocks, NbootstrapReplicates)
    end
    @testset "Random integers" begin
        testDataRandomInts = getTestDataRandomInts(Ndata)
        allTestTypes(testDataRandomInts, statisticHandle, blockLength, Nblocks, NbootstrapReplicates)
    end
    @testset "Linear floats  " begin
        testDataLinearFloats = getTestDataLinearFloats(Ndata)
        allTestTypes(testDataLinearFloats, statisticHandle, blockLength, Nblocks, NbootstrapReplicates)
    end
end
function allTestTypes(testData, statisticHandle, blockLength, Nblocks, NbootstrapReplicates)
    @testset "Sampling tests" begin
        testAllSamplingMethods(testData, blockLength, Nblocks, NbootstrapReplicates)
    end
    @testset "Statistic tests" begin
        testAllStatisticMethods(testData, statisticHandle, blockLength, NbootstrapReplicates)
    end
end
function testAllSamplingMethods(testData, blockLength, Nblocks, NbootstrapReplicates)
    @testset "Sampling method type" begin
        samplingTest(MBBsampler, testData, blockLength, Nblocks, NbootstrapReplicates)
        samplingTest(NBBsampler, testData, blockLength, Nblocks, NbootstrapReplicates)
        samplingTest(CBBsampler, testData, blockLength, Nblocks, NbootstrapReplicates)
    end
end
function testAllStatisticMethods(testData, statisticHandle, blockLength, NbootstrapReplicates)
    @testset "Statistic method type" begin
        statisticTest(MBBsampler, testData, statisticHandle, blockLength, NbootstrapReplicates)
        statisticTest(NBBsampler, testData, statisticHandle, blockLength, NbootstrapReplicates)
        statisticTest(CBBsampler, testData, statisticHandle, blockLength, NbootstrapReplicates)
    end
end
function statisticTest(sampler, testData, statisticHandle, blockLength, NbootstrapReplicates)
    # Testing statistics functions
    dataLength = length(testData)

    # Define bootstrap tests
    samplingFunctionHandle = sampler.handle
    handleStringIn = string(typeof(samplingFunctionHandle).name.mt.name)
    stringToPrint = "Sampling function " * handleStringIn

    @testset "$stringToPrint" begin
        # Set the seed
        Random.seed!(1234)

        # Bootstrap the data
        bootstrapResult = bootstrapStatistic(testData, statisticHandle, sampler, blockLength, NbootstrapReplicates)
        fullDataEstimate = bootstrapResult.fullDataEstimate
        replicateEstimate = bootstrapResult.replicateEstimate
        resampleIndexBB = bootstrapResult.resampleIndexBB
        resampleDataBB = bootstrapResult.resampleDataBB

        # Tests
        @testset "Output size" begin
            @test size(replicateEstimate) == (NbootstrapReplicates,)
            @test size(resampleIndexBB) == (NbootstrapReplicates,)
            @test all(x->x==(dataLength,), size.(resampleIndexBB))
            @test size(resampleDataBB) == (NbootstrapReplicates,)
            @test all(x->x==(dataLength,), size.(resampleDataBB))
        end
        @testset "Output type" begin
            @test typeof(fullDataEstimate) <: Number
            @test typeof(replicateEstimate) <: Array
            @test all(x->x==typeof(fullDataEstimate), typeof.(replicateEstimate))
        end
    end

end
function samplingTest(sampler, testData, blockLength, Nblocks, NbootstrapReplicates)
    # Testing sampling functions
    testDataType = typeof(testData[1])
    bootDataLength = blockLength * Nblocks

    # Define bootstrap tests
    samplingFunctionHandle = sampler.handle
    handleStringIn = string(typeof(samplingFunctionHandle).name.mt.name)
    stringToPrint = "Sampling function " * handleStringIn

    @testset "$stringToPrint" begin
        # Set the seed
        Random.seed!(1234)

        # Resample
        resampleIndexBB, resampleDataBB = samplingFunctionHandle(testData, blockLength, Nblocks, NbootstrapReplicates)

        # Get previous value
        savedIndex = savedTestOutputs(NbootstrapReplicates, handleStringIn)

        # Tests
        @testset "Output size" begin
            @test size(resampleIndexBB) == (NbootstrapReplicates,)
            @test size(resampleDataBB) == (NbootstrapReplicates,)
            @test all(x->x==(bootDataLength,), size.(resampleIndexBB))
            @test all(x->x==(bootDataLength,), size.(resampleDataBB))
        end
        @testset "Output type" begin
            @test all(x->x<:Array, typeof.(resampleIndexBB))
            @test all(x->x<:Array, typeof.(resampleDataBB))
            @test typeof(resampleIndexBB[1][1]) <: Integer
            @test typeof(resampleDataBB[1][1]) <: Real
            @test resampleDataBB[1][1] isa testDataType
        end
        @testset "Output value" begin
            @test resampleIndexBB[1][1] == savedIndex
        end
    end
end
function savedTestOutputs(testSize,methodType)
    # Hardcoded test results for validation
    if methodType == "MBBsample" || methodType == "CBBsample"
        out = 237;
    else
        if testSize == 2
            out = 41;
        elseif testSize == 400
            out = 281;
        end
    end
end

ext: jl


text:
using Base: @kwdef
import Base: isequal, hash

const stdcomment = "// Autogenerated by Bitstream DSL compiler"

struct Variable
    name::Symbol
    type::Symbol
end

getname(x::Variable) = string(x.name)
gettype(x::Variable) = x.type

"""
    Operation

An operation is represented by a list of input types,
a list of output types, and an operator.

# Examples:
```julia
# x represents the addition operation on a pair of
#  stochastic bits (which returns a stochastic bit)
x = Operation([:SBit, :SBit], [:SBit], :+)
```
"""
struct Operation
    inputs::Vector{Symbol}
    outputs::Vector{Symbol}
    operator::Symbol
end
isequal(x::Operation, y::Operation) =
    all(x.inputs .== y.inputs) && all(x.outputs .== y.outputs) && (x.operator == y.operator)
hash(x::Operation) = hash(hash(hash(x.inputs), hash(x.outputs)), hash(x.operator))

opalias(x::Operation) = Operation(x.inputs, x.outputs, opalias(x.operator))
opalias(x) = (x == :÷) ? :div : x

"""
    AbstractHandler

Inherit from this type to create a Verilog generation handler.

This object is responsible for managing all instantiations
of its associated Verilog module.
It should also update the netlist as it adds nets.
"""
abstract type AbstractHandler end

_handlermap = Dict{Operation, Type}()

"""
    register(op::Operation, handler::AbstractHandler)

Register a mapping from an operation to a handler.

# Examples:
```julia
# register a stochastic add to an add handler
register(Operation([:SBit, :SBit], [:SBit], :+), SAddHandler)
```
"""
register(op::Operation, handler::Type{T}) where T <: AbstractHandler =
    _handlermap[opalias(op)] = handler

gethandler(op::Operation) = haskey(_handlermap, op) ? _handlermap[op] :
    error("""
        Cannot find handler for $(op.operator) ($(op.inputs)) -> $(op.outputs).
        Perhaps you forgot to register a handler for your custom operation?
        """)

"""
    HW.allowconstreplace(::Type{T}) where T <: AbstractHandler

Override this function to block constant replacement for the inputs to your handler.

For example, `x + 2` is a node with inputs: `[:x, :2]`.
The constant input could replaced with a Verilog string `Symbol("2'b10")`.
Overriding the return value of this function to `false` will prevent this replacement from happening.
"""
allowconstreplace(::Type{T}) where T <: AbstractHandler = true

"""
    HW.extractrtinfo!(handler::AbstractHandler, innames, outname, opname, inputs, output, op)

Override this function to extract runtime information.

# Fields:
- `handler::AbstractHandler`: an instance of your handler
- `innames::Vector{Symbol}`: the names of the arguments into the node
- `outname::Symbol`: the name of the return out of the node
- `opname::Symbol`: the name of the operator stored in the DFG at this node
- `inputs`: a vector of the actual input argument objects into the node
- `output`: the actual return object from the node
- `op`: the actual function called in Julia (for callable structs this is the struct instance itself)
"""
extractrtinfo!(handler::AbstractHandler, innames, outname, opname, inputs, output, op) = nothing

maptype(x) = @capture(x, Matrix{T_}) ? [Symbol("Array{$T,2}")] :
             @capture(x, Vector{T_}) ? [Symbol("Array{$T,1}")] :
             @capture(x, Number) ? vcat(maptype(:Integer), maptype(:Real)) :
             @capture(x, Real) ? [:Float16, :Float32, :Float64] :
             @capture(x, Integer) ? [:UInt8, :Int8, :UInt16, :Int16, :UInt32, :Int32, :UInt64, :Int64, :UInt128, :Int128] :
             [Symbol(x)]

function register(intypes, outtypes, op, handler::Symbol)
    nintypes = length(intypes)
    typematrix = maptype.(vcat(intypes, outtypes))

    for types in Base.product(typematrix...)
        register(Operation([types[1:nintypes]...], [types[(nintypes + 1):end]...], op), @eval($handler))
    end
end

"""
    HW.@register

A convenience macro for registering a specific handler to a certain operation type.

# Examples

The following code will map any call to `+` with the specified argument/return types
to the `SAddHandler` object.
```julia
@register(SAddHandler, +, begin
    [SBit, SBit] => [SBit]
    [SBit, Vector{SBit}] => [Vector{SBit}]
    [Vector{SBit}, SBit] => [Vector{SBit}]
    [Vector{SBit}, Vector{SBit}] => [Vector{SBit}]
    [SBit, Matrix{SBit}] => [Matrix{SBit}]
    [Matrix{SBit}, SBit] => [Matrix{SBit}]
    [Matrix{SBit}, Matrix{SBit}] => [Matrix{SBit}]
end)
```
"""
macro register(handler, op, ex)
    @capture(unblock(rmlines(ex)), rules__) || error("Could not parse registration.")
    for rule in rules
        @capture(rule, ins_ => outs_) ||
            error("Could not parse registration $rule (use format '[:SBit, :SBit] => [:SBit]')")
        (@capture(ins, [intypes__]) && @capture(outs, [outtypes__])) ||
            error("Specify registration as array (use format '[:SBit, :SBit] => [:SBit]'")
        register(intypes, outtypes, op, handler)
    end

    return :(nothing)
end

include("handlers/utils.jl")
include("handlers/saddhandler.jl")
include("handlers/ssubhandler.jl")
include("handlers/smulthandler.jl")
include("handlers/sdivhandler.jl")
include("handlers/sfdivhandler.jl")
include("handlers/ssqrthandler.jl")
include("handlers/sl2normhandler.jl")
include("handlers/transposehandler.jl")
include("handlers/daddhandler.jl")
include("handlers/dsubhandler.jl")
include("handlers/dmulthandler.jl")
include("handlers/fxpaddhandler.jl")
include("handlers/fxpsubhandler.jl")
include("handlers/fxpmulthandler.jl")
include("handlers/sdmhandler.jl")
include("handlers/delaybufferhandler.jl")


text:
using LibYAML
using Test

@test LibYAML.get_version_string() === "0.2.1"
@test LibYAML.get_version() === (0, 2, 1)


text:
struct CyclicIndex{T <: Integer}
    stop::T
end

CyclicIndex(idx::Base.OneTo) = CyclicIndex(idx.stop)

Base.first(::CyclicIndex) = (1,2)
Base.last(c::CyclicIndex) = (c.stop, 1)
Base.eltype(::Type{CyclicIndex{T}}) where T = Tuple{T,T}

Base.length(c::CyclicIndex) = c.stop
Base.size(c::CyclicIndex) = (length(c), )

function Base.getindex(c::CyclicIndex, i)
    0 < i < c.stop && (i, i+1)
    return (i, 1)
end

Base.iterate(c::CyclicIndex) = c.stop <= 0 ? nothing : first(c), 1

function Base.iterate(c::CyclicIndex, state)
    v = state + 1
    v == c.stop + 1 && return nothing
    v == c.stop && return (v, 1), v
    return (v, v + 1), v
end

struct EdgeIter{T,N, V <: AbstractArray{T,N}}
    v::V
end

Base.length(e::EdgeIter) = length(e.v)
Base.eltype(::Type{EdgeIter{T, N, V}}) where {N,T,V} = Tuple{T,T}

function Base.first(e::EdgeIter)
    if isempty(e.v)
        return nothing
    elseif length(e.v) == 1
        return (e.v[1], e.v[1])
    else
        return (e.v[1],e.v[2])
    end
end

Base.last(e::EdgeIter) = (e.v[end], e.v[1])

Base.iterate(e::EdgeIter) = first(e), 1

function Base.iterate(e::EdgeIter, state)
    sn = state + 1
    sn == length(e.v) + 1 && return nothing
    sn == length(e.v) && return last(e), sn
    return (e.v[sn], e.v[sn + 1]), sn
end

max_stars_repo_name: Michiel-VL/PolyaViz

max_issues_repo_name: Michiel-VL/PolyaViz

max_forks_repo_name: Michiel-VL/PolyaViz


text:
#!/usr/bin/env julia

#
# Hello World server in Julia
# Binds REP socket to tcp://*:5555
# Expects "Hello" from client, replies "World"
#

using ZMQ

context = Context()
socket = Socket(context, REP)
ZMQ.bind(socket, "tcp://*:5555")

while true
    # Wait for next request from client
    message = bytestring(ZMQ.recv(socket))
    println("Received request: $message")

    # Do some 'work'
    sleep(1)

    # Send reply back to client
    ZMQ.send(socket, "World")
end

# classy hit men always clean up when finish the job.
ZMQ.close(socket)
ZMQ.close(context)

max_stars_repo_name: JailbreakFox/LightWeightRepository

max_issues_repo_name: JailbreakFox/LightWeightRepository

max_forks_repo_name: JailbreakFox/LightWeightRepository


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# for converting OpenStreetMap files with road networks to format usable for JEMSS

using JEMSS
import OpenStreetMapX # v0.1.1
const OSM = OpenStreetMapX

include("graph_tools.jl")

# Read OSM file with road network, return nodes and arcs
# kwargs:
# - levels: highway types that will be kept, see ROAD_CLASSES in OSM for more details.
# - boundsLLA: crop map within given bounds.
# Output may need further processing, such as only keeping the largest strongly connected component, and removing duplicate arcs.
# Note that distance calculations done here are different (more accurate) than those in JEMSS.
function readOsmNetworkFile(osmFilename::String;
	levels::Union{Set{Int},Nothing} = nothing, boundsLLA::Union{OSM.Bounds{OSM.LLA},Nothing} = nothing)

	if levels == nothing levels = OSM.ROAD_CLASSES |> values |> Set end

	# read osm file
	@assert(isfile(osmFilename))
	# get osm data; this fails if the osm file does not contain bounds, e.g. <bounds minlat="1.0" minlon="2.0" maxlat="3.0" maxlon="4.0"/>
	osmData = OSM.get_map_data(osmFilename; road_levels = levels, use_cache = false) # 'use_cache = true' ignores other kwargs (e.g. road_levels)
	(nodesENU, roadways) = (osmData.nodes, osmData.roadways)
	@assert(length(nodesENU) > 0)
	@assert(length(roadways) > 0)
	@assert(typeof(first(nodesENU)[2]) == OSM.ENU)

	# crop map within bounds
	if boundsLLA == nothing
		boundsLLA = osmData.bounds
		@assert(boundsLLA != OSM.Bounds{OSM.LLA}(0.0, 0.0, 0.0, 0.0), "missing bounds")
	end
	nodesLLA = OSM.LLA(nodesENU, boundsLLA)
	OSM.crop!(nodesLLA, boundsLLA, roadways)

	# get network data from nodes and ways
	edges, class = OSM.get_edges(nodesENU, roadways) # edges = vector of tuple of vertex keys; class[i] is class (in levels) of edges[i]
	@assert(issubset(Set(class), levels))
	vertices = OSM.get_vertices(edges) # Dict{Int,Int}; vertices[k] gives index (from 1:n) of vertex key k
	weights = OSM.distance(nodesENU, edges) # weights[i] is weight of edges[i]; appears to be equivalent to Geodesy.distance, based on output

	# create nodes from vertices
	numNodes = length(vertices)
	nodes = Vector{Node}(undef, numNodes)
	for (key, i) in vertices
		nodes[i] = Node()
		nodes[i].index = i
		nodes[i].location.x = nodesLLA[key].lon
		nodes[i].location.y = nodesLLA[key].lat
		nodes[i].fields["osm_key"] = key
	end

	# create arcs from edges
	numArcs = length(edges)
	arcs = Vector{Arc}(undef, numArcs)
	for i = 1:numArcs
		v1, v2 = edges[i] # vertex keys (to, from) for edge
		arcs[i] = Arc()
		arcs[i].index = i
		arcs[i].fromNodeIndex = vertices[v1]
		arcs[i].toNodeIndex = vertices[v2]
		arcs[i].fields["osm_class"] = class[i]
		arcs[i].fields["osm_weight"] = weights[i] # edge weight = length
		arcs[i].distance = weights[i] / 1000 # convert metres to km
	end

	return nodes, arcs
end

# merge arc2 into arc1, keeping minimum class (should have higher speed) in arc1
# arc2 will need to be removed in separate step
function mergeDuplicateOsmArcs!(arc1::Arc, arc2::Arc; classFieldName::String = "osm_class")
	@assert(arc1.fromNodeIndex == arc2.fromNodeIndex)
	@assert(arc1.toNodeIndex == arc2.toNodeIndex)
	arc1.fields[classFieldName] = min(arc1.fields[classFieldName], arc2.fields[classFieldName])
	# arcs should already have same weight (length)
end
mergeDuplicateOsmArcs!(arcs::Vector{Arc}, i::Int, j::Int) = mergeDuplicateOsmArcs!(arcs[i], arcs[j])

# Process nodes and arcs from OSM, processing performed:
# - only keep specified highway types (levels)
# - keep only the largest strongly connected component
# - merge duplicate arcs (same start and end node) according to kwarg mergeArcsFunction
# - set offRoadAccess field of nodes, according to kwarg classOffRoadAccess
# - calculate travel time for each arc
# kwargs:
# - levels: highway types that will be kept, see ROAD_CLASSES in OSM for more details.
# - classSpeeds: classSpeeds[i] has dict of the speed (km/hr) of each class of arc, for travel mode i
# - classOffRoadAccess: dict indicating which nodes of different arc classes can be used to get on and off road; node will not have off-road access if all connecting arcs have class returning false in dict
# - mergeArcsFunction: function to merge any duplicate arcs found, given the arcs and two arc indices
function convertOsmNetwork!(nodes::Vector{Node}, arcs::Vector{Arc};
	levels::Union{Set{Int},Nothing} = nothing, classSpeeds::Vector{Dict{Int,Float}} = Dict{Int,Float}[], classOffRoadAccess::Dict{Int,Bool} = Dict{Int,Bool}(), mergeArcsFunction::Function = mergeDuplicateOsmArcs!)

	numTravelModes = length(classSpeeds)

	if levels == nothing levels = OSM.ROAD_CLASSES |> values |> Set end

	# if classOffRoadAccess not given, assume it is true for all levels
	if classOffRoadAccess == Dict()
		classOffRoadAccess = Dict([level => true for level in levels])
	end

	# check that input dicts contain all levels
	for level in levels
		@assert(all(c -> haskey(c, level), classSpeeds))
		@assert(haskey(classOffRoadAccess, level))
	end

	# convert classSpeeds (km/hr) to classInvSpeeds (days/metre)
	# this is needed as the arc lengths (arc field "osm_weight") are in metres
	classInvSpeeds = deepcopy(classSpeeds) # will be days / metre
	for i = 1:numTravelModes
		for (class, speed) in classSpeeds[i]
			classInvSpeeds[i][class] = 1 / (speed * 1000 * 24) # days / metre
		end
	end

	# edit nodes and arcs as needed
	graphRemoveElts!(nodes, arcs, arcFilter = arc -> in(arc.fields["osm_class"], levels)) # keep only specified levels
	graphRemoveDisconnectedArcs!(nodes, arcs) # probably not needed, but use just in case
	graphKeepLargestComponent!(nodes, arcs) # keep only the largest strongly connected component
	graphMergeDuplicateArcs!(nodes, arcs; mergeArcsFunction = mergeArcsFunction)

	# set node offRoadAccess values
	for node in nodes
		node.offRoadAccess = false
	end
	for arc in arcs
		if classOffRoadAccess[arc.fields["osm_class"]]
			nodes[arc.fromNodeIndex].offRoadAccess = true
			nodes[arc.toNodeIndex].offRoadAccess = true
		end
	end

	# set arc travel times
	for i = 1:numTravelModes
		modeFieldName = travelModeString(i)
		for arc in arcs
			class = arc.fields["osm_class"]
			arc.fields[modeFieldName] = arc.fields["osm_weight"] * classInvSpeeds[i][class]
		end
	end
end

function convertOsmNetworkExample!(nodes::Vector{Node}, arcs::Vector{Arc})
	# the 'levels' of highways that will be kept, see ROAD_CLASSES in OSM for more details
	# e.g., levels = Set(1:5) are motorway to tertiary, including links
	# levels = Set(1:6) # roads 1:6 are motorway to residential and unclassified, including links
	levels = Set(1:6)

	# arc speed (km/hr) for each travel mode and each arc class
	v = [
		[110, 110, 60, 60, 60, 60, 60, 60], # travel mode 1, classes 1:8
		[ 90, 90, 45, 45, 45, 45, 45, 45] # travel mode 2, classes 1:8
	]
	classSpeeds = [Dict([j => Float(v[i][j]) for j = 1:length(v[i])]) for i = 1:length(v)]

	# whether nodes along road can be used to get on-road and off-road
	# where set to true, will set both ends of arc to be accessible
	classOffRoadAccess = Dict(
		1 => false, # motorway + link
		2 => false, # trunk + link
		3 => true, # primary + link
		4 => true, # secondary + link
		5 => true, # tertiary + link
		6 => true, # residential + link
	)

	convertOsmNetwork!(nodes, arcs; levels=levels, classSpeeds=classSpeeds, classOffRoadAccess=classOffRoadAccess)
end

function convertOsmNetworkFile(osmFilename::String;
	levels::Union{Set{Int},Nothing} = nothing, boundsLLA::Union{OSM.Bounds{OSM.LLA},Nothing} = nothing, classSpeeds::Vector{Dict{Int,Float}} = Dict{Int,Float}[], classOffRoadAccess::Dict{Int,Bool} = Dict{Int,Bool}(), mergeArcsFunction::Function = mergeDuplicateOsmArcs!)

	(nodes, arcs) = readOsmNetworkFile(osmFilename; levels = levels, boundsLLA = boundsLLA)
	convertOsmNetwork!(nodes, arcs; levels = levels, classSpeeds = classSpeeds, classOffRoadAccess = classOffRoadAccess, mergeArcsFunction = mergeArcsFunction)
	return nodes, arcs
end

# example use of convertOsmNetwork function
function convertOsmNetworkFileExample(osmFilename::String)
	levels = Set(1:6) # roads 1:6 are motorway to residential and unclassified, including links
	(nodes, arcs) = readOsmNetworkFile(osmFilename; levels = levels)
	convertOsmNetworkExample!(nodes, arcs)
	return nodes, arcs
end


text:
"""
    check_filters(node; scale = nothing, symbol = nothing, link = nothing)

Check if the filters are consistant with the mtg onto which they are applied

# Examples

```julia
check_filters(mtg, scale = 1)
check_filters(mtg, scale = (1,2))
check_filters(mtg, scale = (1,2), symbol = "Leaf", link = "<")
```
"""
function check_filters(node; scale = nothing, symbol = nothing, link = nothing) where {T}

    root_node = get_root(node)

    nodeMTG_type = typeof(node.MTG)
    if root_node[:scales] !== nothing
        check_filter(nodeMTG_type, :scale, scale, unique(root_node.attributes[:scales]))
    end

    if root_node[:symbols] !== nothing
        check_filter(nodeMTG_type, :symbol, symbol, unique(root_node.attributes[:symbols]))
    end

    if root_node[:link] !== nothing
        check_filter(nodeMTG_type, :link, link, ("/", "<", "+"))
    end

    return nothing
end

function check_filter(nodetype, type::Symbol, filter, filters)
    if !isnothing(filter)
        filter_type = fieldtype(nodetype, type)
        !(typeof(filter) <: filter_type) &&
            @warn "The $type argument should be of type $filter_type"
        if !(filter in filters)
            @warn "The $type argument should be one of: $filters, and you provided $filter."
        end
    end
end

function check_filter(nodetype, type::Symbol, filter::T, filters) where {T<:Union{Tuple,AbstractArray}}
    for i in filter
        check_filter(nodetype, type, i, filters)
    end
end

"""
    is_filtered(node, scale, symbol, link, filter_fun)

Is a node filtered in ? Returns `true` if the node is kept, `false` if it is filtered-out.
"""
function is_filtered(node, scale, symbol, link, filter_fun)

    link_keep = is_filtered(link, node.MTG.link)
    symbol_keep = is_filtered(symbol, node.MTG.symbol)
    scale_keep = is_filtered(scale, node.MTG.scale)
    filter_fun_keep = isnothing(filter_fun) || filter_fun(node)

    scale_keep && symbol_keep && link_keep && filter_fun_keep
end

function is_filtered(filter, value)
    isnothing(filter) || value in filter
end

function is_filtered(filter::String, value)
    isnothing(filter) || value in (filter,)
end

function is_filtered(filter, value::T) where {T<:Union{Tuple,Array}}
    all(map(x -> is_filtered(filter, x), value))
end

"""
    parse_macro_args(args)

Parse filters and arguments given as a collection of expressions. This function is used to
get the filters as keyword arguments in macros.

# Examples

```julia
args = (:(x = length(node.name)), :(y = node.x + 2), :(scale = 2))
MultiScaleTreeGraph.parse_macro_args(args)
```
"""
function parse_macro_args(args)
    filters = Dict{Symbol,Any}(
        :scale => nothing,
        :symbol => nothing,
        :link => nothing,
        :filter_fun => nothing
    )

    kwargs = Dict{Symbol,Any}(
        :all => true,
        :traversal => AbstractTrees.PreOrderDFS,
    )

    args_array = []

    for i in args
        if i.head == :(=) && i.args[1] in keys(filters)
            filters[i.args[1]] = i.args[2]
        elseif i.head == :(=) && i.args[1] in keys(kwargs)
            kwargs[i.args[1]] = i.args[2]
        else
            push!(args_array, i)
        end
    end
    return (; filters...), (; kwargs...), (args_array...,)
end


text:
# OPM - Load / Create a dataset from flat text files or r datasource

# Set this to true to force a refresh from source datafiles
refresh = false
# Set this to true to install packages from the internet
install = false
# Hardcode a location to the sctfile here.
sctfile = raw"C:\Users\ME\Desktop\data\1973-09-to-2014-06\SCTFILE.TXT"
# This is the location of the root of the files
fileroot = raw"C:\Users\ME\Desktop\data\1973-09-to-2014-06"

# Location of files of interest within our root
subdir = "non-dod/status"
path = joinpath(fileroot, subdir)
# Filename to write our results to disk or read from disk.
result_filename = joinpath(path, "../non-dod-summary.rds")
filtered_filename = joinpath(path, "../non-dod-filtered.rds")
model_filename = joinpath(path, "../non-dod-model.rds")

# List of fields to pull from the files.
fields_of_interest = vcat(:PSEUDO_ID, :FILE_DATE, :AGENCY, :OCCUPATION, :LOS_LEVEL, :PAY_PLAN, :GRADE, :ADJUSTED_BASIC_PAY, :WORK_SCHEDULE)

### Executed code from here on ###
# Get our filelist and parse sctdata
files = [files for (root, dirs, files) in walkdir(path)]
files = [joinpath(path, file) for file in files[1]]
# This loads all files and includes depencenies
include(raw"C:\Users\ME\OneDrive\Documents\julia\FWF\src\FWF.jl")
include("OPMInstallDeps.jl")
sctdata = opm_parse_sctfile_fwf(sctfile)

# For now, load ad-hoc here - julia
files = files[endswith.(files, "_12.txt")]
data = opm_load_filelist(files, sctdata, columns = fields_of_interest, silent=false)
# If refresh...delete our files.
#if (refresh) {
#    file.remove(c(result_filename, filtered_filename, model_filename))
#}

# Load and prepare data...
# otherwise load from source files and sleep for a few hours.
##if (!file.exists(result_filename)) {
#    print(paste(Sys.time(), "No Basic Result File, Reparsing sources."))
#    result <- opm_generate_dataset_from_fwfs(files, sctdata = sctdata, columns = fields_of_interest)
#    write_rds(result, result_filename)
#    # If we wrote this file remove the filtered file, just to be sure
#    file.remove(filtered_filename)
#    print(paste(Sys.time(), "Result file generated"))
#}

#if (!file.exists(filtered_filename)) {
#    print(paste(Sys.time(), "No Filtered File, Performing filtering and grouping"))
#    rm(result)
#    gc()
    # This just always reads from the rds to hopefully allow the above to cleanup memory.
#    result <- read_rds(result_filename) %>% ungroup() %>% opm_filter_unusable_rows() %>% opm_add_years() %>% group_by(`PSEUDO-ID`)
    # Sneak in a gc to eek out a bit more memory.
#    gc()
#    result <- result %>% opm_filter_unusable_ids() %>% opm_number_year_runs() %>% group_by(`PSEUDO-ID`, `group-id`)
#    write_rds(result, filtered_filename)
#    print(paste(Sys.time(), "Filtered file generated"))
#} else if (!file.exists(model_filename)) {
#    print(paste(Sys.time(), "Filtered File found.  Loading for model building."))
#    result <- read_rds(filtered_filename)
#}

#if (!file.exists(model_filename)) {
#    print(paste(Sys.time(), "Model file not found, building model"))
#    model <- opm_generate_model_tibble(result)
#    write_rds(model, model_filename)
#    rm(result)
#    gc()
#   print(paste(Sys.time(), "Model file generated"))
#} else {
#    print(paste(Sys.time(), "Loading model file."))
#    model <- read_rds(model_filename)
#}
# Interesting cases to look at...
# 8,289, 21, 14, 10357722, 166
# 21 is kelli
# 8 departs after 3 years of no increases
# 166, 289, 14 fragmented employment
#
#72738 - termination + long service


text:
using Soss, TupleVectors, MeasureTheory
import SampleChainsDynamicHMC: DynamicHMCChain

m = @model begin
    μ ~ Normal(0,1)
    σ ~ HalfCauchy()
    x ~ Normal(μ,σ)
end;

post = sample(DynamicHMCChain, m() | (;x =3.0)) |> TupleVector

@with post (t_stat=μ/σ,)

import Sobol
ω = SobolHypercube(1)
@with (;ω) post (; x=rand(ω,Normal(μ,σ)))


text:
# # A composite model example with column splits and merges

# The following composite model description comes from a
# discussion at [MLJ issue #166](https://github.com/alan-turing-institute/MLJ.jl/issues/166#issuecomment-533934909):

# > Regress y from x, and classify c from a and b. Then classify w
# > from y and c.

# Below we show how to use MLJ to define a new supervised model type
# `MyComposite` with input (a, b, x) to learn a target (c, y, w)
# according to this prescription. The fields (hyperparmeters) of the
# new composite model will be the two classifiers and regresssor.

# The new model type is obtained by "protyping" the composite model
# using a learning network, and then exporting the network as a
# stand-alone model type.

# Select the relevant MLJ version in the [manual
# entry](https://alan-turing-institute.github.io/MLJ.jl/stable/composing_models/)
# for more on this general procedure.

# To run without issues, this notebook/script should lie in a copy of
# [this
# directory](https://github.com/alan-turing-institute/MLJ.jl/tree/master/examples/composite_example),
# in some tagged release of the [MLJ
# package](https://github.com/alan-turing-institute/MLJ.jl).

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

#-

using MLJ
using Random
Random.seed!(12);

# ### ASSUMPTIONS

# We will make the following assumptions regarding the scientific types
# of the data:

# - `a`, `b`, `x` have scitype  `AbstractMatrix{Continuous}`

# - `c` and `w` have scitype `AbstractVector{<:Finite}`

# - `y` has scitype `AbstractVector{Continuous}`

# All data share the same number of rows (corresponding to
# observations).

# For example,

N = 2
a = fill(1.0, (3N, 3)) + rand(3N, 3)
b = fill(2.0, (3N, 2)) + rand(3N, 2);
x = fill(3.0, (3N, 2)) + rand(3N, 2);
c = categorical(rand("pqr", 3N)); levels!(c, ['p', 'q', 'r'])
w = categorical(rand("PQ", 3N)); levels!(w, ['P', 'Q'])
y = fill(4.0, 3N) + rand(3N);

# I'll suppose the input to our supervised composite model is to be
# presented as a matrix of the form `X = hcat(a, b, x)` where `a`,
# `b`, `x` are of the form above.  For example,

X = hcat(a, b, x);

scitype(X)

# Since the three target variables `c`, `y`, `z` for the composite
# have different types, I'll suppose that these are presented as the
# three columns of a table, with names :c, :y, and :z. For example,

Y = (c=c, y=y, w=w);

scitype(Y)

# We are assuming the learners are:

# - A probabilisitic classifier for learning c from a, b
# - A deterministic regressor for learning y from x
# - A deterministic classifier for learning w from c and y

# Here "classifier" means `AbstractVector{<:Finite}` target scitype,
# and "regressor" means `AbstractVector{<:Continuous}` target scitype.

# We restrict to component models that have `Table(Continous)` input
# scitype and so will need to one-hot encode (c, y), before learning
# w.

# ### PROTYPING THE COMPOSITE MODEL

# Now we define a learning network using component models that will
# become default values for the fields (hyperparameters) of our final
# composite model type.

# We first define ordinary functions to do splitting and merging. The
# functions return a table or vector depending on what the component
# models will be requiring (in this case, tables for inputs, vectors
# for targets):

# Splits:
get_ab(X) = MLJ.table(X[:,1:5], names=[:a1, :a2, :a3, :b1, :b2])
get_x(X)  = MLJ.table(X[:,6:7], names=[:x1, :x2])
get_c(Y)  = Y.c
get_y(Y)  = Y.y
get_w(Y)  = Y.w;

# Merges:
put_cy(c, y) = (c=c, y=y)
put_cyw(c, y, w) = (c=c, y=y, w=w);

#-

get_ab(X) |> pretty
get_x(X) |> pretty
put_cy(c, y) |> pretty
put_cyw(c, y, w) |> pretty

# We now define source nodes. These nodes could simply wrap `nothing`
# instead of concrete data, and the network could still be exported.
# However, to enable testing of the learning network as we build it,
# we will wrap the data defined above. (The author discovered several
# errors in earlier attempts this way.)

X_ = source(X)
Y_ = source(Y, kind=:target)

# Now for the rest of the network.

# Initial splits:
ab_ = node(get_ab, X_)
x_ = node(get_x, X_)
c_ = node(get_c, Y_)
y_ = node(get_y, Y_)
w_ = node(get_w, Y_)

# Node to predict c:
clf1 = @load DecisionTreeClassifier # a model instance
m1 = machine(clf1, ab_, c_)
ĉ_ = predict_mode(m1, ab_)

# Node to predict y:
rgs = @load RidgeRegressor pkg=MultivariateStats
rgs.lambda = 0.1
m = machine(rgs, x_, y_)
ŷ_ = predict(m, x_)

# Merge c and y:
cy_ = node(put_cy, ĉ_, ŷ_)

# Node to do the one-hot-encoding:
hot = OneHotEncoder(drop_last=true)
cy__ = transform(machine(hot, cy_), cy_)

# Node to predict w:
clf2 = @load SVC
m2 = machine(clf2, cy__, w_)
ŵ_ = predict(m2, cy__)

# Final merge:
Ŷ_ = node(put_cyw, ĉ_, ŷ_, ŵ_)

# As a test of functionality, we can fit the final node, which trains
# the whole network...

fit!(Ŷ_, rows=1:2N)

# ... and make a prediction:

Ŷ_(rows=(2N-1):3N)

# ### EXPORT THE LEARNING NETWORK AS STAND-ALONE MODEL

# The next code simultaneously creates a new model type `MyComposite`
# and defines `comp` as an instance, using deep copies of the
# specified learning network component models as default field values:

comp = @from_network MyComposite(classifier1=clf1,
                                 classifier2=clf2,
                                 regressor=rgs) <= Ŷ_

# As a model, this object has no data attached to it. We fit it to
# data, as we do any other model:

X = rand(100, 7);
Y = (c=categorical(rand("abc", 100)),
     y=rand(100),
     w=categorical(rand("AB", 100)));

m = machine(comp, X, Y)
fit!(m, rows=1:80)
Ŷ = predict(m, rows=81:100)
error = sum(Ŷ.w .!= Y.w[81:100])/20

# We can select new component models, for example ...

comp.classifier1 = @load KNNClassifier

# ... and retrain:

fit!(m, rows=1:80)
Ŷ = predict(m, rows=81:100)
error = sum(Ŷ.w .!= Y.w[81:100])/20

using Literate #src
Literate.notebook(@__FILE__, @__DIR__) #src


text:
module QuantumControl

using QuantumPropagators
export propagate, propstep!, propstep

using QuantumControlBase
export ControlProblem, Objective, WeightedObjective, liouvillian
export discretize, discretize_on_midpoints, getcontrols, get_control_parameters
export get_tlist_midpoints
export propagate_objective
export optimize
export chain_infohooks

module Shapes
    using QuantumControlBase.Shapes
    export flattop, box, blackman
end

module Functionals
    using QuantumControlBase.Functionals
    export F_ss, J_T_ss, chi_ss!, F_sm, J_T_sm, chi_sm!, F_re, J_T_re, chi_re!
    export grad_J_T_sm!
end

using Krotov
using GRAPE

using Pkg
using UUIDs

"""Print the versions of the packages constituting QuantumControl."""
function print_versions()
    project_toml = Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))
    version = project_toml["version"]
    direct_deps = project_toml["deps"]
    deps = Pkg.dependencies()
    pkg_names = [name for name in keys(direct_deps) if name ∉ ["Pkg", "UUIDs"]]
    col_width = maximum([length(name) for name in pkg_names])
    for name in reverse(pkg_names)
        pkginfo = deps[UUIDs.UUID(direct_deps[name])]
        if pkginfo.is_tracking_path
            println("$(rpad(name, col_width)): $(pkginfo.version) ($(pkginfo.source))")
        else
            println("$(rpad(name, col_width)): $(pkginfo.version)")
        end
    end
    println("$(rpad("QuantumControl", col_width)): $version")
end

end


text:
"""
    ozone_data

Returns a vector of length 336 = 14*24 with the hourly ozone levels
at Azusa, California, during the first 14 days of July 2014
(discussed on page 319).
"""
function ozone_data()
# Hourly Ozone measurements in Azusa
# Values in parts per million;
# Hourly 12am - 11pm; some measurements missing

k=31;
V = zeros(24, k);
# Tuesday, July 1, 2014
V[:,1] = [ 0.011;  0.011;  0.012;  NaN;    NaN;    0.011;  0.013;  0.014;
           0.022;  0.029;  0.037;  0.050;  0.067;  0.076;  0.086;  0.076;
           NaN;    0.044;  0.035;  0.027;  0.021;  0.018;  0.016;  0.017 ];
# Wednesday, July 2, 2014
V[:,2] = [ 0.018;  0.016;  0.015;  0.016;  NaN;    0.015;  0.013;  0.014;
           NaN;    0.029;  0.042;  0.056;  0.070;  0.083;  0.080;  0.086;
           0.074;  0.067;  0.045;  0.034;  0.024;  0.018;  0.014;  0.012 ];
# Thursday, July 3, 2014
V[:,3] = [ 0.009;  0.009;  0.007;  0.006;  NaN;    0.007;  0.007;  0.014;
           0.018;  0.028;  0.046;  0.069;  0.087;  0.096;  0.105;  0.099;
           0.077;  0.060;  0.046;  0.033;  0.029;  0.018;  0.011;  0.012 ];
# Friday, July 4, 2014
V[:,4] = [ 0.009;  0.008;  0.009;  0.006;  NaN;    0.009;  0.015;  0.020;
           0.028;  0.040;  0.057;  0.073;  0.079;  0.081;  0.070;  0.066;
           0.056;  0.047;  0.040;  0.035;  0.032;  0.026;  0.019;  0.011 ];
# Saturday, July 5, 2014
V[:,5] = [ 0.006;  0.003;  0.004;  0.004;  NaN;    0.004;  0.006;  0.015;
           0.030;  0.047;  0.056;  0.057;  NaN;    NaN;    NaN;    NaN;
           NaN;    NaN;    NaN;    NaN;    0.014;  0.010;  0.008;  0.009 ];
# Sunday, July 6, 2014
V[:,6] = [ 0.008;  0.008;  0.005;  0.005;  NaN;    0.008;  0.013;  0.021;
           0.030;  0.038;  0.050;  NaN;    NaN;    NaN;    NaN;    NaN;
           NaN;    NaN;    NaN;    NaN;    NaN;    NaN;    0.011;  0.009 ];
# Monday, July 7, 2014
V[:,7] = [ 0.007;  0.008;  0.006;  0.003;  NaN;    0.003;  0.006;  0.011;
           0.018;  0.028;  NaN;    NaN;    NaN;    NaN;    NaN;    NaN;
           NaN;    NaN;    NaN;    NaN;    NaN;    NaN;    0.009;  0.014 ];
# Tuesday, July 8, 2014
V[:,8] = [ 0.014;  0.012;  0.011;  NaN;    NaN;    0.007;  0.008;  0.015;
           NaN;    NaN;    NaN;    0.053;  0.070;  NaN;    NaN;    NaN;
           NaN;    NaN;    NaN;    0.031;  0.026;  0.021;  0.019;  0.017 ];
# Wednesday, July 9, 2014
V[:,9] = [ 0.023;  0.020;  0.019;  0.014;  NaN;    0.013;  0.019;  0.016;
           0.023;  0.032;  0.045;  0.057;  0.070;  0.078;  0.082;  0.072;
           0.059;  0.049;  0.039;  0.027;  0.022;  0.019;  0.017;  0.013 ];
# Thursday, July 10, 2014
V[:,10] = [ 0.012;  0.011;  0.011;  NaN;    NaN;    0.013;  0.015;  0.014;
            0.024;  0.033;  0.048;  0.055;  0.069;  0.074;  0.069;  0.062;
            0.052;  0.046;  0.037;  0.028;  0.023;  0.016;  0.011;  0.011 ];
# Friday, July 11, 2014
V[:,11] = [ 0.011;  0.010;  0.011;  0.010;  NaN;    0.004;  0.010;  0.012;
            0.017;  0.022;  0.031;  0.046;  0.051;  0.057;  0.056;  0.041;
            0.034;  0.028;  0.024;  0.020;  0.018;  0.017;  0.017;  0.014 ];
# Saturday, July 12, 2014
V[:,12] = [ 0.011;  0.008;  0.011;  0.009;  NaN;    0.003;  0.007;  0.016;
            0.027;  0.033;  0.039;  0.054;  0.062;  0.070;  0.065;  0.064;
            0.058;  0.046;  0.038;  0.024;  0.023;  0.022;  0.018;  0.014 ];
# Sunday, July 13, 2014
V[:,13] = [ 0.013;  0.009;  0.009;  NaN;    NaN;    0.007;  0.012;  0.021;
            0.033;  0.049;  0.065;  0.079;  0.084;  0.078;  0.077;  0.068;
            0.059;  0.050;  0.038;  0.030;  0.030;  0.032;  0.024;  0.021 ];
# Monday, July 14, 2014
V[:,14] = [ 0.020;  0.017;  0.016;  NaN;    NaN;    0.008;  0.010;  0.013;
            0.013;  0.014;  0.014;  0.019;  0.026;  0.032;  0.038;  0.038;
            0.040;  0.031;  0.030;  0.026;  0.026;  0.023;  0.016;  0.014 ];
# Tuesday, July 15, 2014
V[:,15] = [ 0.013;  0.013;  0.010;  NaN;    NaN;    0.006;  0.010;  0.014;
            0.018;  0.024;  0.030;  0.038;  0.048;  0.052;  0.048;  0.043;
            0.036;  0.026;  0.019;  0.016;  0.015;  0.014;  0.015;  NaN ];
# Wednesday, July 16, 2014
V[:,16] = [ 0.017;  0.015;  0.019;  NaN;    NaN;    0.012;  0.012;  0.012;
            0.014;  0.018;  0.028;  0.036;  0.044;  0.052;  0.058;  NaN;
            NaN;    0.038;  0.028;  0.022;  0.019;  0.016;  0.018;  0.017 ];
# Thursday, July 17, 2014
V[:,17] = [ 0.016;  0.018;  0.017;  0.014;  NaN;    0.017;  0.019;  NaN;
            NaN;    0.015;  0.020;  0.027;  0.037;  0.047;  0.053;  0.038;
            0.035;  0.031;  0.025;  0.020;  0.016;  0.014;  0.013;  0.015 ];
# Friday, July 18, 2014
V[:,18] = [ 0.015;  0.013;  0.012;  NaN;    NaN;    0.008;  0.009;  NaN;
            NaN;    NaN;    NaN;    NaN;    0.044;  NaN;    NaN;    NaN;
            0.037;  0.030;  0.020;  0.013;  0.012;  0.010;  0.008;  0.009 ];
# Saturday, July 19, 2014
V[:,19] = [ 0.010;  0.010;  0.008;  NaN;    NaN;    0.006;  0.007;  0.010;
            0.013;  0.014;  0.016;  0.023;  0.026;  0.031;  0.040;  0.037;
            0.027;  0.024;  0.021;  0.018;  0.015;  0.011;  0.008;  0.007 ];
# Sunday, July 20, 2014
V[:,20] = [ 0.007;  0.005;  0.007;  0.007;  NaN;    0.010;  0.012;  0.016;
            0.019;  0.024;  0.030;  0.040;  0.053;  0.061;  0.068;  0.063;
            0.060;  0.048;  0.028;  0.023;  0.021;  0.018;  0.015;  0.012 ];
# Monday, July 21, 2014
V[:,21] = [ 0.013;  0.011;  0.011;  0.008;  NaN;    0.007;  0.011;  0.012;
            0.016;  0.024;  0.031;  0.038;  0.047;  0.053;  0.048;  0.043;
            0.036;  0.033;  0.027;  0.020;  0.016;  0.012;  0.008;  0.007 ];
# Tuesday, July 22, 2014
V[:,22] = [ 0.008;  0.007;  0.005;  NaN;    NaN;    0.004;  0.006;  NaN;
            NaN;    NaN;    NaN;    0.047;  0.056;  0.058;  0.050;  0.042;
            0.036;  0.031;  0.026;  0.015;  0.010;  0.008;  0.006;  0.007 ];
# Wednesday, July 23, 2014
V[:,23] = [ 0.010;  0.003;  0.003;  NaN;    NaN;    0.004;  0.005;  NaN;
            NaN;    NaN;    NaN;    NaN;    NaN;    NaN;    NaN;    NaN;
            0.042;  0.038;  0.026;  0.019;  0.014;  0.009;  0.005;  0.006 ];
# Thursday, July 24, 2014
V[:,24] = [ 0.004;  0.006;  0.006;  NaN;    NaN;    0.003;  0.005;  0.010;
            0.020;  0.026;  0.045;  0.070;  0.078;  0.094;  0.081;  0.082;
            0.064;  0.056;  0.044;  0.035;  0.024;  0.031;  0.016;  0.006 ];
# Friday, July 25, 2014
V[:,25] = [ 0.009;  0.009;  0.014;  0.018;  NaN;    0.007;  0.008;  0.014;
            0.020;  0.032;  0.041;  0.061;  0.075;  0.081;  0.080;  0.065;
            0.053;  0.046;  0.037;  0.028;  0.027;  0.028;  0.023;  0.020 ];
# Saturday, July 26, 2014
V[:,26] = [ 0.023;  0.021;  0.013;  0.011;  NaN;    0.004;  0.008;  0.014;
            0.023;  0.035;  0.048;  0.062;  0.073;  0.079;  0.074;  0.070;
            0.063;  0.049;  0.041;  0.026;  0.029;  0.025;  0.022;  0.020 ];
# Sunday, July 27, 2014
V[:,27] = [ 0.015;  0.010;  0.013;  NaN;    NaN;    0.008;  0.011;  0.016;
            0.025;  0.033;  0.038;  0.042;  0.058;  0.077;  0.090;  0.086;
            0.064;  0.062;  0.044;  0.039;  0.033;  0.028;  0.030;  0.021 ];
# Monday, July 28, 2014
V[:,28] = [ 0.028;  0.024;  0.017;  NaN;    NaN;    0.012;  0.014;  0.010;
            0.013;  0.024;  0.032;  0.043;  0.051;  0.063;  0.070;  0.064;
            0.055;  0.043;  0.037;  0.033;  0.024;  0.018;  0.013;  0.012 ];
# Tuesday, July 29, 2014
V[:,29] = [ 0.011;  0.013;  0.018;  NaN;    NaN;    0.004;  0.005;  0.009;
            0.017;  0.031;  0.037;  0.044;  0.052;  0.060;  0.057;  0.058;
            0.055;  0.045;  0.037;  0.029;  0.020;  0.013;  0.009;  NaN ];
# Wednesday, July 30, 2014
V[:,30] = [ NaN;    0.006;  0.004;  0.003;  NaN;    0.004;  0.005;  0.008;
            0.019;  0.034;  0.045;  0.064;  0.082;  0.093;  0.081;  0.079;
            0.071;  0.056;  0.035;  0.029;  0.026;  0.014;  0.006;  0.005 ];
# Thursday, July 31, 2014
V[:,31] = [ 0.007;  0.004;  0.004;  0.003;  NaN;    0.003;  0.005;  0.011;
            0.019;  0.032;  0.045;  0.056;  0.069;  0.079;  0.076;  0.075;
            0.063;  0.055;  0.041;  0.026;  0.021;  0.012;  0.008;  0.011 ];
k =14;
return reshape(V[:,1:k], k*24, 1)
end

max_issues_repo_name: lynnlangit/julia-linear-algebra


text:
module CheckConstraintsModule

import ..CoreModule: Node, Options
import ..EquationUtilsModule: compute_complexity

# Check if any binary operator are overly complex
function flag_bin_operator_complexity(
    tree::Node, ::Val{op}, options::Options
)::Bool where {op}
    if tree.degree == 0
        return false
    elseif tree.degree == 1
        return flag_bin_operator_complexity(tree.l, Val(op), options)
    else
        if tree.op == op
            overly_complex::Bool = (
                (
                    (options.bin_constraints[op][1]::Int > -1) && (
                        compute_complexity(tree.l, options) >
                        options.bin_constraints[op][1]::Int
                    )
                ) || (
                    (options.bin_constraints[op][2]::Int > -1) && (
                        compute_complexity(tree.r, options) >
                        options.bin_constraints[op][2]::Int
                    )
                )
            )
            if overly_complex
                return true
            end
        end
        return (
            flag_bin_operator_complexity(tree.l, Val(op), options) ||
            flag_bin_operator_complexity(tree.r, Val(op), options)
        )
    end
end

# Check if any unary operators are overly complex
function flag_una_operator_complexity(
    tree::Node, ::Val{op}, options::Options
)::Bool where {op}
    if tree.degree == 0
        return false
    elseif tree.degree == 1
        if tree.op == op
            overly_complex::Bool = (
                (options.una_constraints[op]::Int > -1) &&
                (compute_complexity(tree.l, options) > options.una_constraints[op]::Int)
            )
            if overly_complex
                return true
            end
        end
        return flag_una_operator_complexity(tree.l, Val(op), options)
    else
        return (
            flag_una_operator_complexity(tree.l, Val(op), options) ||
            flag_una_operator_complexity(tree.r, Val(op), options)
        )
    end
end

"""Count the max number of times an operator of a given degree is nested"""
function count_max_nestedness(tree::Node, degree::Int, op::Int, options::Options)::Int
    if tree.degree == 0
        return 0
    elseif tree.degree == 1
        count = (degree == 1 && tree.op == op) ? 1 : 0
        return count + count_max_nestedness(tree.l, degree, op, options)
    else  # tree.degree == 2
        count = (degree == 2 && tree.op == op) ? 1 : 0
        return count + max(
            count_max_nestedness(tree.l, degree, op, options),
            count_max_nestedness(tree.r, degree, op, options),
        )
    end
end

# function fast_max_nestedness(tree::Node, degree::Int, op_idx::Int, nested_degree::Int, nested_op_idx::Int, options::Options)::Int
function fast_max_nestedness(
    tree::Node,
    degree::Int,
    op_idx::Int,
    nested_degree::Int,
    nested_op_idx::Int,
    options::Options,
)::Int
    # Don't need to branch - once you find operator, run
    # count_max once, then return. Don't need to go deeper!
    if tree.degree == 0
        return 0
    elseif tree.degree == 1
        if degree != tree.degree || tree.op != op_idx
            return fast_max_nestedness(
                tree.l, degree, op_idx, nested_degree, nested_op_idx, options
            )
        end
        return count_max_nestedness(tree.l, nested_degree, nested_op_idx, options)
    else
        if degree != tree.degree || tree.op != op_idx
            return max(
                fast_max_nestedness(
                    tree.l, degree, op_idx, nested_degree, nested_op_idx, options
                ),
                fast_max_nestedness(
                    tree.r, degree, op_idx, nested_degree, nested_op_idx, options
                ),
            )
        end
        return max(
            count_max_nestedness(tree.l, nested_degree, nested_op_idx, options),
            count_max_nestedness(tree.r, nested_degree, nested_op_idx, options),
        )
    end
end

"""Check if there are any illegal combinations of operators"""
function flag_illegal_nests(tree::Node, options::Options)::Bool
    # We search from the top first, then from child nodes at end.
    nested_constraints = options.nested_constraints
    if nested_constraints === nothing
        return false
    end
    for (degree, op_idx, op_constraint) in nested_constraints
        for (nested_degree, nested_op_idx, max_nestedness) in op_constraint
            nestedness = fast_max_nestedness(
                tree, degree, op_idx, nested_degree, nested_op_idx, options
            )
            if nestedness > max_nestedness
                return true
            end
        end
    end
    return false
end

"""Check if user-passed constraints are violated or not"""
function check_constraints(tree::Node, options::Options, maxsize::Int)::Bool
    if compute_complexity(tree, options) > maxsize
        return false
    end
    for i in 1:(options.nbin)
        if options.bin_constraints[i] == (-1, -1)
            continue
        elseif flag_bin_operator_complexity(tree, Val(i), options)
            return false
        end
    end
    for i in 1:(options.nuna)
        if options.una_constraints[i] == -1
            continue
        elseif flag_una_operator_complexity(tree, Val(i), options)
            return false
        end
    end
    if flag_illegal_nests(tree, options)
        return false
    end

    return true
end

function check_constraints(tree::Node, options::Options)::Bool
    return check_constraints(tree, options, options.maxsize)
end

end


text:
module TestProgressRecord

using ProgressLogging: Progress, ProgressString
using Random: shuffle!
using Test
using UUIDs: uuid4

@testset "fraction :: Int" begin
    @test Progress(id = uuid4(), fraction = 0).fraction === 0.0
    @test Progress(id = uuid4(), fraction = 1).fraction === 1.0
    @test Progress(uuid4(), 0).fraction === 0.0
    @test Progress(uuid4(), 1).fraction === 1.0
end

@testset "cmp: ProgressString" begin
    strings = Any[
        let s = string(i)
            rand(Bool) ? s : ProgressString(Progress(uuid4(), 0.0; name = s))
        end
        for i in 10:99
    ]
    shuffle!(strings)
    @test string.(sort!(strings)) == string.(10:99)
end

end  # module


text:
num_vars_to_next_mutation_point(probMutation) = ceil(Int, (-log(rand())) / probMutation)

"""
Provides `apply()` operator that mutates one specified dimension of a parameter
vector.
"""
abstract type GibbsMutationOperator <: MutationOperator end

# apply operator to each dimension
function apply!(m::GibbsMutationOperator, params::AbstractVector{<:Real}, target_index::Int)
    @inbounds @simd for i in eachindex(params)
        params[i] = apply(m, params[i], i, target_index)
    end
    return params
end

"""
Uniform mutation of a parameter vector.
"""
struct UniformMutation{SS<:SearchSpace} <: GibbsMutationOperator
    search_space::SS

    UniformMutation(ss::SS) where {SS<:SearchSpace} = new{SS}(ss)
end

search_space(m::UniformMutation) = m.search_space

apply(m::UniformMutation, v::Number, dim::Int, target_index::Int) =
    dimmin(search_space(m), dim) + rand() * dimdelta(search_space(m), dim)

"""
Mutation clock operator is a more efficient way to mutate vectors than to generate
a random value per variable in the vectors. It instead generates the number of variables
to skip until the next mutation. Then it uses a sub-mutation operator to do the actual
mutation. This is based on the paper:
    Deb and Deb (2012), "Analyzing Mutation Schemes for Real-Parameter Genetic Algorithms"
but we use a Poisson distribution.
"""
mutable struct MutationClock{S<:GibbsMutationOperator} <: MutationOperator
    inner::S
    rate::Float64           # Probability to mutate a dimension
    nextVarToMutate::Int    # dimension index - 1

    MutationClock(inner::S, rate::Float64 = 0.05) where {S<:GibbsMutationOperator} =
        new{S}(inner, rate, 1 + num_vars_to_next_mutation_point(rate))
end

function apply!(mc::MutationClock, v::AbstractVector{<:Real}, target_index::Int)
    n = length(v)
    i = mc.nextVarToMutate
    @inbounds while i <= n
        v[i] = apply(mc.inner, v[i], i, target_index)
        i += num_vars_to_next_mutation_point(mc.rate)
    end
    @assert i > n
    mc.nextVarToMutate = i - n
    return v
end


text:
using RomeoApp
@time msg = unwrapping_main(ARGS)
println(msg)


text:
function calculate_transcriptional_control_variable(W::Array{Float64,1}, binding_rate::Float64, order_param::Float64, frac_bound::Float64)::Float64

    # initialize
    v_variable = 0.0
    v_num = 0.0
    Z = 0.0

    # define states involved in activity
    all_states = [0,1,2]
    activity_states = [1,2]
    f_dependent_states = [2]

    # rename parameters
    ni = order_param
    Ki = binding_rate
    x = frac_bound

    # v-variable function
    for (index, value) ∈ enumerate(all_states)
        Wf = W[index]

        if value ∈ f_dependent_states
            Wf *= (x/Ki)^ni / (1 + (x/Ki)^ni)
        end

        if value ∈ activity_states
            v_num += Wf
        end

        Z += Wf
    end

    v_variable = v_num / Z

    # return
    return v_variable
end


text:
@inline initial_condition_generator(rng,n,model) = (φ=[4], X=[2.0], Y=[0.0])

const tau = first_exit_y(0.0,model.dq.model.b+1.0)

const x = range(0,1.6,length=10_001)


text:
using Plots

cList = [complex(-0.4,-0.6), complex(0,-1),
            complex(-0.12,-0.75), complex(-0.6,0),
                complex(-0.8,0.16), complex(-0.4,0.6)]

function JuliaSet(c)
    abst = zeros((601,601))
    k = (1 + sqrt(1 + (4 * abs(c)))) / 2
    inum = 0
    for i in -1.5:0.005:1.5
        inum += 1
        jnum = 0
        for j in -1.5:0.005:1.5
            jnum += 1
            z = complex(i,j)
            for n in 1:70
                z = (z^2) + c
                if abs(z) > k
                    abst[inum,jnum] = abs(z)
                    break
                end
            end
        end
    end
    return abst
end

theme(:dark)
heatmap(abst, c = cgrad(:thermal, rev = true), legend = false, border=:none)

max_stars_repo_name: shahmari/A-few-fractals-in-Julia

max_issues_repo_name: shahmari/A-few-fractals-in-Julia

max_forks_repo_name: shahmari/A-few-fractals-in-Julia


text:
"""

```
cutoffdistances!(cutoff::Float64,
                 x_solute::AbstractVector{T},
                 x_solvent::Vector{T},
                 lc_solvent::LinkedCells,
                 box::Box,
                 dc::CutoffDistances) where T
```

This routine that returns a list of the distances
between atoms that are smaller than a specified cutoff,
for a given set of coordinates.

Returns `nd`, the number of distances smaller than the cutoff, and modifies `dc`.

"""
function cutoffdistances!(
    cutoff::Float64,
    x_solute::AbstractVector{T},
    x_solvent::Vector{T},
    lc_solvent::LinkedCells,
    box::Box,
    dc::CutoffDistances,
) where {T}

    # Reset the dc structure
    reset!(dc)

    # Loop over solute atoms
    for iat in eachindex(x_solute)
        xat = x_solute[iat]
        # Check the cell of this atom
        i, j, k = icell3D(xat, box)
        # Loop over vicinal cells to compute distances to solvent atoms, and
        # add data to dc structure (includes current cell)
        for ic = i-box.lcell:i+box.lcell
            for jc = j-box.lcell:j+box.lcell
                for kc = k-box.lcell:k+box.lcell
                    cutoffdcell!(
                        cutoff,
                        iat,
                        xat,
                        x_solvent,
                        lc_solvent,
                        box,
                        ic,
                        jc,
                        kc,
                        dc,
                    )
                end
            end
        end
    end

    return nothing
end


text:
module DescriptorZoo

using ASE
using JuLIP

using Reexport

@reexport using ACSF

@reexport using AMP

@reexport using CHSF

@reexport using SHIPs

#@reexport using QUIP

@reexport using SOAPlite

@reexport using DScribe

@reexport using MBTR

end # module


text:
export RemoteEstObj,
    LocalEstObj,
    EstimationWrapper,
    update_reo!,
    reset_reo!,
    theta0, theta1, hess, invhess, grad,
    getworkers,
    parallel_simloglik!,
    update!,
    start_up_workers

"workers() but excluding master"
getworkers() = filter(i -> i != 1, workers())

function start_up_workers(ENV::Base.EnvDict; nprocs = Sys.CPU_THREADS)
    # send primitives to workers
    oldworkers = getworkers()
    println_time_flush("removing workers $oldworkers")
    rmprocs(oldworkers)
    flush(stdout)
    if "SLURM_JOBID" in keys(ENV)
        num_cpus_to_request = parse(Int, ENV["SLURM_TASKS_PER_NODE"])
        println("requesting $(num_cpus_to_request) cpus from slurm.")
        flush(stdout)
        pids = addprocs_slurm(num_cpus_to_request)
    else
        cputhrds = Sys.CPU_THREADS
        cputhrds < nprocs && @warn "using nprocs = $cputhrds < $nprocs specified"
        pids = addprocs(min(nprocs, cputhrds))
    end
    println_time_flush("Workers added: $pids")
    return pids
end

# ----------------------------------

@GenGlobal g_RemoteEstObj

# ----------------------------------

abstract type AbstractEstimationObjects end

struct RemoteEstObj{D<:DataSetofSets} <: AbstractEstimationObjects
    data::D
    llvec::SharedVector{Float64}
    scoremat::SharedMatrix{Float64}
    sim::SimulationDrawsMatrix{Float64,Matrix{Float64}}

    function RemoteEstObj(data::D,llvec,scoremat,sim) where {D}
        k, n = size(scoremat)
        n == num_i(data) == length(llvec) == num_i(sim) || throw(DimensionMismatch("num_i"))
        k == _nparm(data) || throw(DimensionMismatch("k"))
        _nparm(drill(data)) == _nparm(sim) || throw(DimensionMismatch("sim k"))
        return new{D}(data, llvec, scoremat, sim)
    end
end

struct LocalEstObj{D<:DataSetofSets} <: AbstractEstimationObjects
    data::D
    theta0::Vector{Float64}
    theta1::Vector{Float64}
    grad::Vector{Float64}
    hess::Matrix{Float64}
    invhess::Matrix{Float64}

    function LocalEstObj(data::D,theta0,theta1,grad,hess,invhess) where {D}
        k = length(theta0)
        k == length(theta1) == length(grad) || throw(DimensionMismatch())
        k == checksquare(hess) == checksquare(invhess) || throw(DimensionMismatch())
        k == _nparm(data) || throw(DimensionMismatch())

        return new{D}(data,theta0,theta1,grad,hess,invhess)
    end
end

data(x::AbstractEstimationObjects) = x.data

llvec(x::RemoteEstObj) = x.llvec
scoremat(x::RemoteEstObj) = x.scoremat
sim(x::RemoteEstObj) = x.sim
size(x::RemoteEstObj) = size(scoremat(x))
num_i(x::RemoteEstObj) = length(llvec(eo))
_nparm(x::RemoteEstObj) = size(scoremat(x),1)
LL(x::RemoteEstObj) = sum(llvec(x))
negLL(x::RemoteEstObj) = -LL(x)
loglikelihood(x::RemoteEstObj) = LL(x)

theta0(x::LocalEstObj) = x.theta0
theta1(x::LocalEstObj) = x.theta1
grad(  x::LocalEstObj) = x.grad
hess(  x::LocalEstObj) = x.hess
invhess(x::LocalEstObj) = x.invhess
num_i( x::LocalEstObj) = num_i(data(x))
_nparm(x::LocalEstObj) = length(theta0(x))

function reset!(x::LocalEstObj)
    fill!(grad(x), 0)
    fill!(hess(x), 0)
    fill!(invhess(x),0)
    fill!(theta1(x),0)
end

function update!(data::DataSetofSets, theta::Vector)
    _nparm(data) == length(theta) || throw(DimensionMismatch())
    thetas = split_thetas(data, theta)
    for (d, θ) in zip(data, thetas)
        update!(d, θ)
    end
end

@noinline function update_reo!(reo::RemoteEstObj, theta::Vector)
    update!(data(reo), theta)
    return nothing
end

update_reo!(theta::Vector) = update_reo!(get_g_RemoteEstObj(), theta)

@noinline function reset_reo!(reo::RemoteEstObj)
    data_drill = drill(data(reo))
    ddm = _model(data_drill)
    vf = value_function(ddm)
    fill!(vf, 0)
end

reset_reo!() = reset_reo!(get_g_RemoteEstObj())

function invhess!(x::LocalEstObj)
    invhess(x) .= inv(hess(x))
    return invhess(x)
end

function RemoteEstObj(leo::LocalEstObj, M, pids = workers())
    n = num_i(leo)
    k = _nparm(leo)
    llvec    = SharedVector{Float64}(n;    pids=pids)
    scoremat = SharedMatrix{Float64}(k, n; pids=pids)
    sim      = SimulationDraws(M, data(leo))
    return RemoteEstObj(data(leo), llvec, scoremat, sim)
end

function LocalEstObj(data, theta)
    k = length(theta)
    k == _nparm(data) || throw(DimensionMismatch("length(theta) = $(length(theta)) but _nparm(data) = $(_nparm(data)) θ = $theta"))

    theta0 = Vector{Float64}(undef, k)
    theta1 = similar(theta0)
    grad   = similar(theta0)
    hess   = Matrix{Float64}(undef, k, k)
    invhess = similar(hess)
    return LocalEstObj(data, theta0, theta1, grad, hess, invhess)
end

struct EstimationWrapper{D} <: AbstractEstimationObjects
    leo::LocalEstObj{D}
    reo::RemoteEstObj{D}
    function EstimationWrapper(leo::LocalEstObj{D},reo::RemoteEstObj{D}) where {D}
        @assert data(leo) === data(reo)
        return new{D}(leo,reo)
    end
end

LocalEstObj(ew::EstimationWrapper) = ew.leo
RemoteEstObj(ew::EstimationWrapper) = ew.reo
data(ew::EstimationWrapper) = data(RemoteEstObj(ew))
OneTo(ew::EstimationWrapper) = OneTo(num_i(LocalEstObj(ew)))
_nparm(ew::EstimationWrapper) = _nparm(LocalEstObj(ew))

check_theta(ew::EstimationWrapper, theta) = check_theta(data(ew),theta)

function check_theta(data::DataSetofSets, theta)
    thetasvw = split_thetas(data, theta)
    for (dat,θ) in zip(data,thetasvw)
        update!(dat, θ)
    end
end

function update!(ew::EstimationWrapper, theta, dograd)
    l = LocalEstObj(ew)
    r = RemoteEstObj(ew)
    theta1(l) .= theta
    if dograd
        score = scoremat(r)
        h = hess(l)
        g = grad(l)
        mul!(h, score, score')
        sum!(reshape(g, :, 1), score)
        # println("gradient = $g")
        g .*= -1
    end
    nll = negLL(r)
    if dograd
        countplus!(-nll,theta,g)
    else
        countplus!(-nll,theta)
    end
    flush(stdout)
    return nll
end

@noinline function simloglik!(i, theta, dograd, reo::RemoteEstObj; kwargs...)

    check_finite(theta)

    gradi = view(scoremat(reo), :, i)
    grptup = getindex.(data(reo), i)
    simi = view(sim(reo), i)
    θρ = theta_ρ(data(reo), theta)

    fill!(_qm(sim(reo)), 0)
    fill!(gradi, 0)

    # do updates
    @assert isfinite(θρ)
    update!(simi, θρ)

    thetasvw = split_thetas(data(reo), theta)
    idxs = theta_indexes(data(reo))

    llv = llvec(reo)
    llv[i] = simloglik!(gradi, grptup, thetasvw, idxs, simi, dograd; kwargs...)

end

simloglik!(i, theta, dograd; kwargs...) =
    simloglik!(i, theta, dograd, get_g_RemoteEstObj(); kwargs...)

function serial_simloglik!(ew, theta, dograd; n::Integer=num_i(LocalEstObj(ew)), kwargs...)
    check_theta(ew,theta)
    reo = RemoteEstObj(ew)
    update_reo!(reo, theta)
    map(i -> simloglik!(i, theta, dograd, reo; kwargs...), OneTo(n))
    return update!(ew, theta, dograd)
end

function parallel_simloglik!(ew, theta, dograd; kwargs...)
    check_finite(theta)
    check_theta(ew,theta)
    wp = CachingPool(workers())
    @eval @everywhere update_reo!($theta)
    let x=theta, dograd=dograd, kwargs=kwargs
        pmap(i -> simloglik!(i, x, dograd; kwargs...), wp, OneTo(ew))
    end
    return update!(ew, theta, dograd)
end

function test_parallel_simloglik!(ew, theta, dograd; kwargs...)
    @everywhere reset_reo!()
    parallel_simloglik!(ew, theta, dograd; kwargs...)
end

function test_serial_simloglik!(ew, theta, dograd; kwargs...)
    reset_reo!(RemoteEstObj(ew))
    serial_simloglik!(ew, theta, dograd; kwargs...)
end


text:
module ProximalAlgorithms

using ProximalOperators

const RealOrComplex{R} = Union{R,Complex{R}}
const Maybe{T} = Union{T,Nothing}

# utilities

include("utilities/ad.jl")
include("utilities/conjugate.jl")
include("utilities/fb_tools.jl")
include("utilities/iteration_tools.jl")

# acceleration operators

include("accel/lbfgs.jl")
include("accel/anderson.jl")
include("accel/nesterov.jl")
include("accel/broyden.jl")
include("accel/noaccel.jl")

# algorithms

include("algorithms/forward_backward.jl")
include("algorithms/fast_forward_backward.jl")
include("algorithms/zerofpr.jl")
include("algorithms/panoc.jl")
include("algorithms/douglas_rachford.jl")
include("algorithms/drls.jl")
include("algorithms/primal_dual.jl")
include("algorithms/davis_yin.jl")
include("algorithms/li_lin.jl")
include("algorithms/fista.jl")

end # module


text:
## Copyright (c) 2015, Júlio Hoffimann Mendes <juliohm@stanford.edu>
##
## Permission to use, copy, modify, and/or distribute this software for any
## purpose with or without fee is hereby granted, provided that the above
## copyright notice and this permission notice appear in all copies.
##
## THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
## WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
## MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
## ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
## WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
## ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
## OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

function imfilter_gpu{T<:Real,K<:Real,N}(img::AbstractArray{T,N}, kern::AbstractArray{K,N}, border::AbstractString)
  # GPU metadata
  ctx = GPU.ctx; queue = GPU.queue
  mult_kernel = GPU.mult_kernel

  # operations with complex type
  img = map(Complex64, img)
  kern = map(Complex64, kern)

  # Int prefix is a workaround for julia #15276
  prepad  = Int[div(size(kern,i)-1, 2) for i = 1:N]
  postpad = Int[div(size(kern,i),   2) for i = 1:N]
  fullpad = Int[nextprod([2,3], size(img,i) + prepad[i] + postpad[i]) - size(img, i) - prepad[i] for i = 1:N]

  A = border == "inner" ? img : padarray(img, prepad, fullpad, border, Complex64(0))

  # OpenCL FFT expects powers of 2, 3, 5, 7, 11 or 13
  paddings = clfftpad(A)
  A = padarray(A, zeros(Int, ndims(A)), paddings, "symmetric")

  krn = zeros(Complex64, size(A))
  indexesK = ntuple(d->[size(krn,d)-prepad[d]+1:size(krn,d);1:size(kern,d)-prepad[d]], N)
  krn[indexesK...] = reflect(kern)

  # plan FFT
  p = clfft.Plan(Complex64, ctx, size(A))
  clfft.set_layout!(p, :interleaved, :interleaved)
  clfft.set_result!(p, :inplace)
  clfft.bake!(p, queue)

  # populate GPU memory
  bufA = cl.Buffer(Complex64, ctx, :copy, hostbuf=A)
  bufkrn = cl.Buffer(Complex64, ctx, :copy, hostbuf=krn)
  bufRES = cl.Buffer(Complex64, ctx, length(A))

  # compute ifft(fft(A).*fft(kern))
  clfft.enqueue_transform(p, :forward, [queue], bufA, nothing)
  clfft.enqueue_transform(p, :forward, [queue], bufkrn, nothing)
  cl.call(queue, mult_kernel, length(A), nothing, bufA, bufkrn, bufRES)
  clfft.enqueue_transform(p, :backward, [queue], bufRES, nothing)

  # get result back
  AF = reshape(cl.read(queue, bufRES), size(A))

  # undo OpenCL FFT paddings
  AF = padarray(AF, zeros(Int, ndims(AF)), -paddings, "symmetric")

  out = []
  if border == "inner"
    out = Array(realtype(eltype(AF)), ([size(img)...] - prepad - postpad)...)
    indexesA = ntuple(d->postpad[d]+1:size(img,d)-prepad[d], N)
    copyreal!(out, AF, indexesA)
  else
    out = Array(realtype(eltype(AF)), size(img))
    indexesA = ntuple(d->postpad[d]+1:size(img,d)+postpad[d], N)
    copyreal!(out, AF, indexesA)
  end

  out
end

@generated function reflect{T,N}(A::AbstractArray{T,N})
    quote
        B = Array(T, size(A))
        @nexprs $N d->(n_d = size(A, d)+1)
        @nloops $N i A d->(j_d = n_d - i_d) begin
            @nref($N, B, j) = @nref($N, A, i)
        end
        B
    end
end

for N = 1:5
    @eval begin
        function copyreal!{T<:Real}(dst::Array{T,$N}, src, I::Tuple{Vararg{UnitRange{Int}}})
            @nexprs $N d->(I_d = I[d])
            @nloops $N i dst d->(j_d = first(I_d)+i_d-1) begin
                (@nref $N dst i) = real(@nref $N src j)
            end
            dst
        end
        function copyreal!{T<:Complex}(dst::Array{T,$N}, src, I::Tuple{Vararg{UnitRange{Int}}})
            @nexprs $N d->I_d = I[d]
            @nloops $N i dst d->(j_d = first(I_d)+i_d-1) begin
                (@nref $N dst i) = @nref $N src j
            end
            dst
        end
    end
end

realtype{R<:Real}(::Type{R}) = R
realtype{R<:Real}(::Type{Complex{R}}) = R


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

function calc_error_norms(func, u, t, analyzer,
                          mesh::AbstractMeshData{NDIMS}, equations, initial_condition,
                          dg::DGMulti{NDIMS}, cache, cache_analysis) where {NDIMS}
  rd = dg.basis
  md = mesh.md
  @unpack u_values = cache

  # interpolate u to quadrature points
  apply_to_each_field(mul_by!(rd.Vq), u_values, u)

  component_l2_errors = zero(eltype(u_values))
  component_linf_errors = zero(eltype(u_values))
  for i in each_quad_node_global(mesh, dg, cache)
    u_exact = initial_condition(getindex.(md.xyzq, i), t, equations)
    error_at_node = func(u_values[i], equations) - func(u_exact, equations)
    component_l2_errors += md.wJq[i] * error_at_node.^2
    component_linf_errors = max.(component_linf_errors, abs.(error_at_node))
  end
  return sqrt.(component_l2_errors), component_linf_errors
end

function integrate(func::Func, u,
                   mesh::AbstractMeshData,
                   equations, dg::DGMulti, cache; normalize=true) where {Func}
  rd = dg.basis
  md = mesh.md
  @unpack u_values = cache

  # interpolate u to quadrature points
  apply_to_each_field(mul_by!(rd.Vq), u_values, u)

  integral = sum(md.wJq .* func.(u_values, equations))
  if normalize == true
    integral /= sum(md.wJq)
  end
  return integral
end

function analyze(::typeof(entropy_timederivative), du, u, t,
                 mesh::AbstractMeshData, equations, dg::DGMulti, cache)

  rd = dg.basis
  md = mesh.md
  @unpack u_values = cache

  # interpolate u, du to quadrature points
  du_values = similar(u_values) # Todo: simplices. Can we move this to the analysis cache somehow?
  apply_to_each_field(mul_by!(rd.Vq), du_values, du)
  apply_to_each_field(mul_by!(rd.Vq), u_values, u)

  # compute ∫v(u) * du/dt = ∫dS/dt. We can directly compute v(u) instead of computing the entropy
  # projection here, since the RHS will be projected to polynomials of degree N and testing with
  # the L2 projection of v(u) would be equivalent to testing with v(u) due to the moment-preserving
  # property of the L2 projection.
  dS_dt = zero(eltype(first(du)))
  for i in Base.OneTo(length(md.wJq))
    dS_dt += dot(cons2entropy(u_values[i],equations), du_values[i]) * md.wJq[i]
  end
  return dS_dt
end

function create_cache_analysis(analyzer, mesh::AbstractMeshData,
                               equations, dg::DGMulti, cache,
                               RealT, uEltype)
  md = mesh.md

  return (; )
end

SolutionAnalyzer(rd::RefElemData) = rd

nelements(mesh::AbstractMeshData, solver::DGMulti, cache) = mesh.md.num_elements

end # @muladd


text:
using Test
using LatticeTools

@testset "iscompatible" begin
    # hypercube and operation
    @testset "hypercube-operation" begin
        #hypercube = HypercubicLattice([2 0; 0 2])
        hypercube = Hypercube([2 0; 0 2])
        @test_throws DimensionMismatch iscompatible(hypercube, TranslationOperation([1, 0, 0]))
        @test iscompatible(hypercube, TranslationOperation([1, 0]))

        @test_throws DimensionMismatch iscompatible(hypercube, PointOperation([1 0 0; 0 1 0; 0 0 0]))
        @test iscompatible(hypercube, PointOperation([0 1; 1 0]))
        @test !iscompatible(hypercube, PointOperation([0 2; 1 0]))
        @test !iscompatible(Hypercube([3 0; 0 2]), PointOperation([0 1; 1 0]))
    end

    # hypercube and symmetry
    @testset "hypercube-symmetry" begin
        hypercube = Hypercube([2 0; 0 2])
        @test  iscompatible(hypercube, FiniteTranslationSymmetry([2 0; 0 2]))
        @test  iscompatible(hypercube, FiniteTranslationSymmetry([2 2; 0 2]))
        @test !iscompatible(hypercube, FiniteTranslationSymmetry([1 0; 0 1]))
        @test !iscompatible(hypercube, FiniteTranslationSymmetry([2 0; 0 1]))

        psym = project(PointSymmetryDatabase.find("4mm"), [1 0 0; 0 1 0])
        @test  iscompatible(hypercube, psym)
        @test  iscompatible(Hypercube([2 2; 0 2]), psym)
        @test !iscompatible(Hypercube([3 0; 0 2]), psym)
    end

    # translation symmetry and point operation
    @testset "translation-point" begin
        @test  iscompatible(FiniteTranslationSymmetry([2 0; 0 2]), PointOperation([0 1; 1 0]))
        @test  iscompatible(FiniteTranslationSymmetry([2 2; 0 2]), PointOperation([0 1; 1 0]))
        @test !iscompatible(FiniteTranslationSymmetry([3 0; 0 2]), PointOperation([0 1; 1 0]))

        psym = project(PointSymmetryDatabase.find("4mm"), [1 0 0; 0 1 0])
        @test  iscompatible(FiniteTranslationSymmetry([2 0; 0 2]), psym)
        @test  iscompatible(FiniteTranslationSymmetry([2 2; 0 2]), psym)
        @test !iscompatible(FiniteTranslationSymmetry([3 0; 0 2]), psym)

        @test  iscompatible(psym, FiniteTranslationSymmetry([2 0; 0 2]))
        @test  iscompatible(psym, FiniteTranslationSymmetry([2 2; 0 2]))
        @test !iscompatible(psym, FiniteTranslationSymmetry([3 0; 0 2]))
    end

    @testset "translation-irrep-point-operation" begin
        # shape matching, momentum check
        tsym = FiniteTranslationSymmetry([3 0; 0 3])
        for irrep_index in 1:num_irreps(tsym)
            k = tsym.fractional_momenta[irrep_index]
            @test  iscompatible(tsym, irrep_index, PointOperation([0 1; 1  0])) == (k[1] == k[2])
            @test  iscompatible(tsym, irrep_index, PointOperation([1 0; 0 -1])) == (k[2] == 0 || k[2] == 1//2)
        end
        tsym = FiniteTranslationSymmetry([3 3; 0 3])
        for irrep_index in 1:num_irreps(tsym)
            k = tsym.fractional_momenta[irrep_index]
            @test  iscompatible(tsym, irrep_index, PointOperation([0 1; 1  0])) == (k[1] == k[2])
            @test  iscompatible(tsym, irrep_index, PointOperation([1 0; 0 -1])) == (k[2] == 0 || k[2] == 1//2)
        end

        # shape mismatch
        tsym = FiniteTranslationSymmetry([3 0; 0 2])
        for irrep_index in 1:num_irreps(tsym)
            k = tsym.fractional_momenta[irrep_index]
            @test !iscompatible(tsym, irrep_index, PointOperation([0 1; 1 0])) # not compatible with tsym
            @test  iscompatible(tsym, irrep_index, PointOperation([1 0; 0 -1])) == (k[2] == 0 || k[2] == 1//2)
        end

        psym = project(PointSymmetryDatabase.find("4mm"), [1 0 0; 0 1 0])
        @test  iscompatible(FiniteTranslationSymmetry([2 0; 0 2]), psym)
        @test  iscompatible(FiniteTranslationSymmetry([2 2; 0 2]), psym)
        @test !iscompatible(FiniteTranslationSymmetry([3 0; 0 2]), psym)
    end

    @testset "translation-irrep-point-symmetry" begin
        tsym = FiniteTranslationSymmetry([3 0; 0 3])
        psym = project(PointSymmetryDatabase.find("4mm"), [1 0 0; 0 1 0])
        for irrep_index in 1:num_irreps(tsym)
            k = tsym.fractional_momenta[irrep_index]
            @test iscompatible(tsym, irrep_index, psym) == (k in [[0//1, 0//1], [1//2, 1//2]])
            # 4mm only at Gamma and M points
        end

        tsym = FiniteTranslationSymmetry([4 0; 0 4])
        psym = project(PointSymmetryDatabase.find("4mm"), [1 0 0; 0 1 0])
        for irrep_index in 1:num_irreps(tsym)
            k = tsym.fractional_momenta[irrep_index]
            @test iscompatible(tsym, irrep_index, psym) == (k in [[0//1, 0//1], [1//2, 1//2]])
            # 4mm only at Gamma and M points
        end

        tsym = FiniteTranslationSymmetry([4 4; 0 4])
        psym = project(PointSymmetryDatabase.find("4mm"), [1 0 0; 0 1 0])
        for irrep_index in 1:num_irreps(tsym)
            k = tsym.fractional_momenta[irrep_index]
            @test iscompatible(tsym, irrep_index, psym) == (k in [[0//1, 0//1], [1//2, 1//2]])
            # 4mm only at Gamma and M points
        end

        # shape mismatch
        tsym = FiniteTranslationSymmetry([4 2; 0 4])
        psym = project(PointSymmetryDatabase.find("4mm"), [1 0 0; 0 1 0])
        for irrep_index in 1:num_irreps(tsym)
            k = tsym.fractional_momenta[irrep_index]
            @test !iscompatible(tsym, irrep_index, psym)
            # 4mm already not compatible with lattice shape
        end
    end

    @testset "lattice-symmetry" begin
        unitcell = makeunitcell([1.0 0.0; 0.0 1.0]; SiteType=String)
        addsite!(unitcell, "X", FractCoord([0, 0], [0.5, 0.0]))

        lattice = make_lattice(unitcell, [2 0; 0 2])
        @test_throws DimensionMismatch iscompatible(lattice, TranslationOperation([1,2,3]))
        for i in -5:5, j in -5:5
            @test iscompatible(lattice, TranslationOperation([i, j]))
        end
        @test iscompatible(lattice, PointOperation([1 0; 0 -1]))
        @test !iscompatible(lattice, PointOperation([0 1; 1 0]))

        @test iscompatible(lattice, FiniteTranslationSymmetry([2 0; 0 2]))
        @test iscompatible(lattice, FiniteTranslationSymmetry([2 2; 0 2]))
        # ^- This is a choice. Whether to require equality, or just compatibility

        psym = project(PointSymmetryDatabase.get(7), [1 0 0; 0 1 0])
        @test iscompatible(lattice, psym)

        psym = project(PointSymmetryDatabase.get(13), [1 0 0; 0 1 0])
        @test !iscompatible(lattice, psym)
    end
end


text:
# module Bukdu.Octo

export singularize, pluralize

function singularize(s::String)::String
    word = lowercase(s)
    m = match(r"(\w*)s", word)
    if isa(m, RegexMatch)
        return first(m.captures)
    else
        return word
    end
end

function pluralize(s::String)::String
    word = lowercase(s)
    string(word, "s")
end

# irregular


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "IVerilog"
version = v"12.0.0"

# Collection of sources required to complete build
sources = [
   GitSource("https://github.com/steveicarus/iverilog.git", "2693dd32b075243cca20400cf3a808cef119477e")
]

dependencies = [
    HostBuildDependency("Bison_jll"),
    HostBuildDependency("gperf_jll"),
    Dependency("Readline_jll"; compat="8.1.1"),
    Dependency("Zlib_jll"; compat="1.2.12"),
]

# Bash recipe for building across all platforms
script = raw"""
cd iverilog
export CFLAGS="-O2"
export CXXFLAGS="-O2"
sh ./autoconf.sh
# The generated configure script is replacing a variable incorrectly
sed -i 's/EXEEXT=$ac_cv_build_exeext/BUILD_EXEEXT=$ac_cv_build_exeext/' configure
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("iverilog", :iverilog),
    ExecutableProduct("iverilog-vpi", :iverilog_vpi),
    ExecutableProduct("vvp", :vvp)
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: staticfloat/BinaryBuilderRepos

max_issues_repo_name: staticfloat/BinaryBuilderRepos

max_forks_repo_name: staticfloat/BinaryBuilderRepos