text:
# Not included or used because there is no speedup whatsoever

function cuda_kernel_shared!(previous_ptr::Ptr{T}, current_ptr::Ptr{T}, q_ptr::Ptr, h::Int, w::Int, λ::T, γ::T) where T
    y = (blockIdx().y-1) * blockDim().y + threadIdx().y
    x = (blockIdx().x-1) * blockDim().x + threadIdx().x
    Ψₜ₋₁ = CuDeviceArray((h,w), previous_ptr)
    Ψₜ   = CuDeviceArray((h,w), current_ptr)
    Ψₜ₊₁ = CuDeviceArray((h,w), previous_ptr)
    q    = CuDeviceArray((h,w), q_ptr)
    λsq = λ*λ
    λhalf = λ/2
    eins = one(T)
    zwei = eins + eins

    if 1 < y < h && 1 < x < w
        S = @cuDynamicSharedMem(T, (2, blockDim().y+2, blockDim().x+2))
        if threadIdx().y == 1 || y == 2
            S[1, threadIdx().y,   threadIdx().x+1] = T(q[y-1,x])
            S[2, threadIdx().y,   threadIdx().x+1] = Ψₜ[y-1,x]
        elseif threadIdx().y == blockDim().y || y == h-1
            S[1, threadIdx().y+2, threadIdx().x+1] = T(q[y+1,x])
            S[2, threadIdx().y+2, threadIdx().x+1] = Ψₜ[y+1,x]
        end
        if threadIdx().x == 1 || x == 2
            S[1, threadIdx().y+1, threadIdx().x] = T(q[y,x-1])
            S[2, threadIdx().y+1, threadIdx().x] = Ψₜ[y,x-1]
        elseif threadIdx().x == blockDim().x || x == w-1
            S[1, threadIdx().y+1, threadIdx().x+2] = T(q[y,x+1])
            S[2, threadIdx().y+1, threadIdx().x+2] = Ψₜ[y,x+1]
        end
        S[1, threadIdx().y+1, threadIdx().x+1] = T(q[y,x])
        S[2, threadIdx().y+1, threadIdx().x+1] = Ψₜ[y,x]

        sync_threads()

        q₂₋  = S[1, threadIdx().y+1, threadIdx().x  ]
        Ψₜ₂₋ = S[2, threadIdx().y+1, threadIdx().x  ]
        q₁₋  = S[1, threadIdx().y,   threadIdx().x+1]
        Ψₜ₁₋ = S[2, threadIdx().y,   threadIdx().x+1]
        Ψₜₛ  = S[2, threadIdx().y+1, threadIdx().x+1]
        q₁₊  = S[1, threadIdx().y+2, threadIdx().x+1]
        Ψₜ₁₊ = S[2, threadIdx().y+2, threadIdx().x+1]
        q₂₊  = S[1, threadIdx().y+1, threadIdx().x+2]
        Ψₜ₂₊ = S[2, threadIdx().y+1, threadIdx().x+2]

        K = q₁₋ + q₁₊ + q₂₋ + q₂₊

        Q = q₁₋ * Ψₜ₁₋ + q₁₊ * Ψₜ₁₊ +
            q₂₋ * Ψₜ₂₋ + q₂₊ * Ψₜ₂₊

        γ̄ = γ * ((eins-q₁₋) + (eins-q₁₊) + (eins-q₂₋) + (eins-q₂₊))

        Ψₜ₊₁[y,x] = eins/(eins+λhalf*γ̄) * ((zwei-λsq*K) * Ψₜₛ + (λhalf*γ̄-eins) * Ψₜ₋₁[y,x] + λsq*Q)
    end

    return nothing
end

function cuda_kernel_shared!(previous_ptr::Ptr{T}, current_ptr::Ptr{T}, q_ptr::Ptr, h::Int, w::Int, v::Int, λ::T, γ::T) where T
    y = (blockIdx().y-1) * blockDim().y + threadIdx().y
    x = (blockIdx().x-1) * blockDim().x + threadIdx().x
    z = (blockIdx().z-1) * blockDim().z + threadIdx().z
    Ψₜ₋₁ = CuDeviceArray((h,w,v), previous_ptr)
    Ψₜ   = CuDeviceArray((h,w,v), current_ptr)
    Ψₜ₊₁ = CuDeviceArray((h,w,v), previous_ptr)
    q    = CuDeviceArray((h,w,v), q_ptr)
    λsq = λ*λ
    λhalf = λ/2
    eins = one(T)
    zwei = eins + eins

    if 1 < y < h && 1 < x < w && 1 < z < v
        S = @cuDynamicSharedMem(T, (2, blockDim().y+2, blockDim().x+2, blockDim().z+2))
        sy = threadIdx().y+1
        sx = threadIdx().x+1
        sz = threadIdx().z+1
        if threadIdx().y == 1 || y == 2
            S[1, sy-1, sx, sz] = T(q[y-1, x, z])
            S[2, sy-1, sx, sz] =  Ψₜ[y-1, x, z]
        elseif threadIdx().y == blockDim().y || y == h-1
            S[1, sy+1, sx, sz] = T(q[y+1, x, z])
            S[2, sy+1, sx, sz] =  Ψₜ[y+1, x, z]
        end
        if threadIdx().x == 1 || x == 2
            S[1, sy, sx-1, sz] = T(q[y, x-1, z])
            S[2, sy, sx-1, sz] =  Ψₜ[y, x-1, z]
        elseif threadIdx().x == blockDim().x || x == w-1
            S[1, sy, sx+1, sz] = T(q[y, x+1, z])
            S[2, sy, sx+1, sz] =  Ψₜ[y, x+1, z]
        end
        if threadIdx().z == 1 || z == 2
            S[1, sy, sx, sz-1] = T(q[y, x, z-1])
            S[2, sy, sx, sz-1] =  Ψₜ[y, x, z-1]
        elseif threadIdx().z == blockDim().z || z == v-1
            S[1, sy, sx, sz+1] = T(q[y, x, z+1])
            S[2, sy, sx, sz+1] =  Ψₜ[y, x, z+1]
        end
        S[1, sy, sx, sz] = T(q[y, x, z])
        S[2, sy, sx, sz] =  Ψₜ[y, x, z]

        sync_threads()

        q₃₋  = S[1, sy,   sx,   sz-1]
        Ψₜ₃₋ = S[2, sy,   sx,   sz-1]
        q₂₋  = S[1, sy,   sx-1, sz]
        Ψₜ₂₋ = S[2, sy,   sx-1, sz]
        q₁₋  = S[1, sy-1, sx,   sz]
        Ψₜ₁₋ = S[2, sy-1, sx,   sz]
        Ψₜₛ  = S[2, sy,   sx,   sz]
        q₁₊  = S[1, sy+1, sx,   sz]
        Ψₜ₁₊ = S[2, sy+1, sx,   sz]
        q₂₊  = S[1, sy,   sx+1, sz]
        Ψₜ₂₊ = S[2, sy,   sx+1, sz]
        q₃₊  = S[1, sy,   sx,   sz+1]
        Ψₜ₃₊ = S[2, sy,   sx,   sz+1]

        K = q₁₋ + q₁₊ + q₂₋ + q₂₊ + q₃₋ + q₃₊

        Q = q₁₋ * Ψₜ₁₋ + q₁₊ * Ψₜ₁₊ +
            q₂₋ * Ψₜ₂₋ + q₂₊ * Ψₜ₂₊ +
            q₃₋ * Ψₜ₃₋ + q₃₊ * Ψₜ₃₊

        γ̄ = γ * ((eins-q₁₋) + (eins-q₁₊) + (eins-q₂₋) + (eins-q₂₊) + (eins-q₃₋) + (eins-q₃₊))

        Ψₜ₊₁[y,x,z] = eins/(eins+λhalf*γ̄) * ((zwei-λsq*K) * Ψₜₛ + (λhalf*γ̄-eins) * Ψₜ₋₁[y,x,z] + λsq*Q)
    end

    return nothing
end


text:
"""
    ∂regularize!(∂, rbm; l2_fields = 0, l1_weights = 0, l2_weights = 0, l2l1_weights = 0)

Updates RBM gradients `∂`, with the regularization gradient.
"""
function ∂regularize!(
    ∂::NamedTuple, # unregularized gradient
    rbm::RBM;
    l2_fields::Real = 0, # L2 regularization of visible unit fields
    l1_weights::Real = 0, # L1 regularization of weights
    l2_weights::Real = 0, # L2 regularization of weights
    l2l1_weights::Real = 0, # L2/L1 regularziation of weights (10.7554/eLife.39397, Eq. 8)
    zerosum::Bool = false # whether to zerosum gradients
)
    if !iszero(l2_fields)
        ∂regularize_fields!(∂.visible, rbm.visible; l2_fields)
    end
    if !iszero(l1_weights)
        ∂.w .+= l1_weights * sign.(rbm.w)
    end
    if !iszero(l2_weights)
        ∂.w .+= l2_weights * rbm.w
    end
    if !iszero(l2l1_weights)
        dims = ntuple(identity, ndims(rbm.visible))
        ∂.w .+= l2l1_weights * sign.(rbm.w) .* mean(abs, rbm.w; dims)
    end
    zerosum && zerosum!(∂, rbm)
    return ∂
end

function ∂regularize_fields!(
    ∂::NamedTuple, layer::Union{Binary,Spin,Potts,Gaussian,ReLU,xReLU,pReLU}; l2_fields::Real = 0
)
    if !iszero(l2_fields)
        ∂.θ .+= l2_fields * layer.θ
    end
    return ∂
end

function ∂regularize_fields!(∂::NamedTuple, layer::dReLU; l2_fields::Real = 0)
    if !iszero(l2_fields)
        ∂.θp .+= l2_fields * layer.θp
        ∂.θn .+= l2_fields * layer.θn
    end
    return ∂
end

function ∂regularize(
    rbm::RBM;
    l2_fields::Real = 0, # L2 regularization of visible unit fields
    kw... # weights penalties
)
    visible = ∂regularize_fields(rbm.visible; l2_fields)
    w = ∂regularize_weights(rbm; kw...)
    return (; visible, hidden = nothing, w)
end

∂regularize_fields(layer::Union{Binary,Spin,Potts}; l2_fields::Real = 0) = (
    θ = l2_fields * layer.θ,
)
∂regularize_fields(layer::Union{Gaussian,ReLU}; l2_fields::Real = 0) = (
    θ = l2_fields * layer.θ, γ = zero(layer.γ)
)
∂regularize_fields(layer::dReLU; l2_fields::Real = 0) = (
    θp = l2_fields * layer.θp, γp = zero(layer.γp),
    θn = l2_fields * layer.θn, γn = zero(layer.γn)
)
∂regularize_fields(layer::pReLU; l2_fields::Real = 0) = (
    θ = l2_fields * layer.θ, γ = zero(layer.γ), Δ = zero(layer.Δ), η = zero(layer.η)
)
∂regularize_fields(layer::xReLU; l2_fields::Real = 0) = (
    θ = l2_fields * layer.θ, γ = zero(layer.γ), Δ = zero(layer.Δ), ξ = zero(layer.ξ)
)

function ∂regularize_weights(
    rbm::RBM;
    l1_weights::Real = 0, # L1 regularization of weights
    l2_weights::Real = 0, # L2 regularization of weights
    l2l1_weights::Real = 0 # L2/L1 regularziation of weights (10.7554/eLife.39397, Eq. 8)
)
    dims = ntuple(identity, ndims(rbm.visible))
    ∂l2l1 = l2l1_weights * sign.(rbm.w) .* mean(abs, rbm.w; dims)
    ∂l1 = l1_weights * sign.(rbm.w)
    ∂l2 = l2_weights * rbm.w
    return ∂l2l1 + ∂l1 + ∂l2
end

ext: jl


text:
module WJFKNN
export WjfKNN, knnSearch
using WJFParallelTask

struct CoverSet
    mapping::Dict{UInt32,UInt32}
    subCenters::Vector{UInt32}
    subCoverSetInfos::Vector{Tuple{UInt32,UInt32,Float64,Float64}}
    subCoverSets::Vector{CoverSet}
end

function findCenterSequentially(
    data::Vector,
    items::Vector{UInt32},
    maxNumSamples::UInt32,
    distFun::Function,
)::UInt32
    numItems = length(items)
    numSamples = min(UInt32(ceil(sqrt(numItems))), maxNumSamples)
    samples::Set{UInt32} = Set{UInt32}()
    for idxSample = 1:numSamples
        sample = (rand(UInt32) % numItems) + 1
        while sample in samples
            sample = (rand(UInt32) % numItems) + 1
        end
        push!(samples, sample)
    end
    sampleLists::Vector{UInt32} = [sample for sample in samples]

    mySmallestSumD = Inf64
    myBestSample = UInt32(0)
    for u::UInt32 = 1:UInt32(numSamples)
        sumD::Float64 = 0.0
        for v = 1:numSamples
            sumD += distFun(data[sampleLists[u]], data[sampleLists[v]])
        end
        if sumD < mySmallestSumD
            mySmallestSumD = sumD
            myBestSample = sampleLists[u]
        end
    end
    return items[myBestSample]
end

function findCenter(
    data::Vector,
    items::Vector{UInt32},
    maxNumSamples::UInt32,
    distFun::Function,
)::UInt32
    numItems = length(items)
    numSamples = min(UInt32(ceil(sqrt(numItems))), maxNumSamples)
    samples::Set{UInt32} = Set{UInt32}()
    for idxSample = 1:numSamples
        sample = (rand(UInt32) % numItems) + 1
        while sample in samples
            sample = (rand(UInt32) % numItems) + 1
        end
        push!(samples, sample)
    end
    sampleLists::Vector{UInt32} = [sample for sample in samples]

    function mapFun(u0::UInt32, u1::UInt32)::Tuple{Float64, UInt32}
        mySmallestSumD = Inf64
        myBestSample = UInt32(0)
        for u = u0:u1
            sumD::Float64 = 0.0
            for v = 1:numSamples
                sumD += distFun(data[sampleLists[u]], data[sampleLists[v]])
            end
            if sumD < mySmallestSumD
                mySmallestSumD = sumD
                myBestSample = sampleLists[u]
            end
        end # u
        return (mySmallestSumD, myBestSample)
    end

    function reduceFun(values::Vector{Tuple{Float64, UInt32}}
        )::Tuple{Float64, UInt32}
        mySmallestSumD = Inf64
        myBestSample = UInt32(0)
        for (sumD, sample) in values
            if sumD < mySmallestSumD
                mySmallestSumD = sumD
                myBestSample = sample
            end
        end
        return (mySmallestSumD, myBestSample)
    end

    bestSample = 0
    while bestSample == 0
        (smallestSumD, bestSample) = mapReduce(
            UInt32(1), numSamples, UInt32(100),
            mapFun, reduceFun, (Inf64, UInt32(0)))
    end

    return items[bestSample]
end

function findFarthestItemSequentially(
    data::Vector,
    items::Vector{UInt32},
    center,
    subCenters::Vector{UInt32},
    distFun::Function,
)::UInt32
    numItems = length(items)
    myFarthestdist = 0.0
    myFarthestItem = 0
    for u::UInt32 = 1:UInt32(numItems)
        item = items[u]
        mydist = distFun(data[item], center)
        for j in subCenters
            mydist = min(mydist, distFun(data[item], data[j]))
        end
        if mydist > myFarthestdist
            myFarthestItem = items[u]
            myFarthestdist = mydist
        end
    end # u
    return myFarthestItem
end

function findFarthestItem(
    data::Vector,
    items::Vector{UInt32},
    center,
    subCenters::Vector{UInt32},
    distFun::Function,
)::UInt32
    function mapFun(u0::UInt32, u1::UInt32)::Tuple{Float64,UInt32}
        myFarthestdist = 0.0
        myFarthestItem = 0
        for u = u0:u1
            item = items[u]
            mydist = distFun(data[item], center)
            for j in subCenters
                mydist = min(mydist, distFun(data[item], data[j]))
            end
            if mydist > myFarthestdist
                myFarthestItem = items[u]
                myFarthestdist = mydist
            end
        end # u
        return (myFarthestdist, myFarthestItem)
    end

    function reduceFun(values::Vector{Tuple{Float64,UInt32}}
        )::Tuple{Float64,UInt32}
        myFarthestdist = 0.0
        myFarthestItem = 0
        for (dist, item) in values
            if dist > myFarthestdist
                myFarthestdist = dist
                myFarthestItem = item
            end
        end
        return (myFarthestdist, myFarthestItem)
    end

    numItems = length(items)
    farthestItem = 0
    while farthestItem == 0
        (farthestdist, farthestItem) = mapReduce(
            UInt32(1), UInt32(numItems), UInt32(100),
            mapFun, reduceFun, (0.0, UInt32(0))
        )
    end
    return farthestItem
end

function fillSubCenters!(
    subCenters::Vector{UInt32},
    data::Vector,
    items::Vector{UInt32},
    maxNumSamples::UInt32,
    numPartitionsPerNode::UInt32,
    distFun::Function
)
    center = data[findCenterSequentially(data, items, maxNumSamples, distFun)]
    while length(subCenters) < numPartitionsPerNode
        push!(subCenters, findFarthestItemSequentially(
            data, items, center, subCenters, distFun))
    end
end

function findSubCenters(
    data::Vector,
    items::Vector{UInt32},
    maxNumSamples::UInt32,
    numPartitionsPerNode::UInt32,
    distFun::Function
)::Vector{UInt32}
    center = data[findCenter(data, items, maxNumSamples, distFun)]
    subCenters = Vector{UInt32}()
    while length(subCenters) < numPartitionsPerNode
        push!(subCenters, findFarthestItem(data, items, center, subCenters, distFun))
    end
    return subCenters
end

function fillSubCoverSetInfos!(
    subCoverSetInfos::Vector{Tuple{UInt32,UInt32,Float64,Float64}},
    data::Vector,
    items::Vector{UInt32},
    subCenters::Vector{UInt32},
    distFun::Function,
)
    numItems = length(items)
    mySubCenterDists = fill(0.0, (length(subCenters),))
    mySubCenterRanks = Vector{UInt32}(undef, length(subCenters))
    for u::UInt32 = 1:UInt32(numItems)
        item = items[u]
        for i::UInt32 = 1:length(subCenters)
            mySubCenterRanks[i] = i
            mySubCenterDists[i] =
                distFun(data[item], data[subCenters[i]])
        end
        sort!(mySubCenterRanks, by = x -> mySubCenterDists[x])
        subCoverSetInfos[u] = (
            mySubCenterRanks[1],
            mySubCenterRanks[2],
            mySubCenterDists[mySubCenterRanks[1]],
            mySubCenterDists[mySubCenterRanks[2]],
        )
    end
end

function computeSubCoverSetInfos(
    data::Vector,
    items::Vector{UInt32},
    subCenters::Vector{UInt32},
    distFun::Function,
)::Vector{Tuple{UInt32,UInt32,Float64,Float64}}
    numItems = length(items)
    subCoverSetInfos::Vector{Tuple{Int32,Int32,Float64,Float64}} =
        Vector{Tuple{Int32,Int32,Float64,Float64}}(undef, numItems)
    function mapFun(u0::UInt32, u1::UInt32)
        mySubCenterDists = fill(0.0, (length(subCenters),))
        mySubCenterRanks = Vector{UInt32}(undef, length(subCenters))
        for u = u0:u1
            item = items[u]
            for i::UInt32 = 1:length(subCenters)
                mySubCenterRanks[i] = i
                mySubCenterDists[i] =
                    distFun(data[item], data[subCenters[i]])
            end
            sort!(mySubCenterRanks, by = x -> mySubCenterDists[x])
            subCoverSetInfos[u] = (
                mySubCenterRanks[1],
                mySubCenterRanks[2],
                mySubCenterDists[mySubCenterRanks[1]],
                mySubCenterDists[mySubCenterRanks[2]],
            )
        end # u
    end
    mapOnly(UInt32(1), UInt32(numItems), UInt32(100), mapFun)
    return subCoverSetInfos
end

function getSubCoverItems(
    items::Vector{UInt32},
    subCoverSetInfos::Vector{Tuple{UInt32,UInt32,Float64,Float64}},
    numPartitionsPerNode::UInt32,
)::Vector{Vector{UInt32}}
    result::Vector{Vector{UInt32}} = Vector{Vector{UInt32}}()
    for i = 1:numPartitionsPerNode
        push!(result, Vector{UInt32}())
    end
    for i = 1:length(items)
        (j1, j2, d1, d2) = subCoverSetInfos[i]
        push!(result[j1], items[i])
        push!(result[j2], items[i])
    end
    return result
end

function fillCoverSet!(
    coverSet::CoverSet,
    data::Vector,
    items::Vector{UInt32},
    maxNumSamples::UInt32,
    numPartitionsPerNode::UInt32,
    distFun::Function,
    smallSubCoverSets::Vector{Tuple{CoverSet, Vector{UInt32}}},
)
    fillSubCenters!(
        coverSet.subCenters,
        data,
        items,
        maxNumSamples,
        numPartitionsPerNode,
        distFun
    )
    fillSubCoverSetInfos!(
        coverSet.subCoverSetInfos,
        data,
        items,
        coverSet.subCenters,
        distFun,
    )
    subCoverItems::Vector{Vector{UInt32}} =
        getSubCoverItems(items, coverSet.subCoverSetInfos, numPartitionsPerNode)
    for i = 1:numPartitionsPerNode
        coverSet.subCoverSets[i] = CoverSet(
            data,
            subCoverItems[i],
            maxNumSamples,
            numPartitionsPerNode,
            distFun,
            smallSubCoverSets
        )
    end
end

function CoverSet(
    data::Vector,
    items::Vector{UInt32},
    maxNumSamples::UInt32,
    numPartitionsPerNode::UInt32,
    distFun::Function,
    smallSubCoverSets::Vector{Tuple{CoverSet, Vector{UInt32}}},
)
    numItems = length(items)
    numData = length(data)
    mapping::Dict{UInt32,UInt32} = Dict{UInt32,UInt32}()
    for i = 1:numItems
        mapping[items[i]] = i
    end
    smallSubCoverSetSize = maxNumSamples * UInt32(sqrt(numPartitionsPerNode))
    if length(items) <= maxNumSamples
        return CoverSet(
        mapping,
        Vector{UInt32}(),
        Vector{Tuple{UInt32,UInt32,Float64,Float64}}(),
        Vector{CoverSet}())
    elseif numData > numItems && numItems <= smallSubCoverSetSize
        coverSet = CoverSet(
        mapping,
        Vector{UInt32}(),
        Vector{Tuple{UInt32,UInt32,Float64,Float64}}(undef, numItems),
        Vector{CoverSet}(undef, numPartitionsPerNode))
        push!(smallSubCoverSets, (coverSet,items))
        return coverSet
    else
        subCenters = findSubCenters(
            data,
            items,
            maxNumSamples,
            numPartitionsPerNode,
            distFun
        )
        subCoverSetInfos =
            computeSubCoverSetInfos(data, items, subCenters, distFun)
        subCoverItems::Vector{Vector{UInt32}} =
            getSubCoverItems(items, subCoverSetInfos, numPartitionsPerNode)
        subCoverSets = [
            CoverSet(
                data,
                subCoverItems[i],
                maxNumSamples,
                numPartitionsPerNode,
                distFun,
                smallSubCoverSets
            ) for i = 1:numPartitionsPerNode
        ]
        return CoverSet(mapping, subCenters, subCoverSetInfos, subCoverSets)
    end
end

struct WjfKNN
    data::Vector
    distFun::Function
    maxNumSamples::UInt32
    numPartitionsPerNode::UInt32
    entry::CoverSet
end

function finishSmallSubCoverSets(
    smallSubCoverSets::Vector{Tuple{CoverSet, Vector{UInt32}}},
    data::Vector,
    maxNumSamples::UInt32,
    numPartitionsPerNode::UInt32,
    distFun::Function,
    )::Vector{Tuple{CoverSet, Vector{UInt32}}}
    function mapFun(u0::UInt32,u1::UInt32)
        myResult = Vector{Tuple{CoverSet, Vector{UInt32}}}()
        for u in u0:u1
            subCoverSet = smallSubCoverSets[u][1]
            items = smallSubCoverSets[u][2]
            fillCoverSet!(
                subCoverSet,
                data,
                items,
                maxNumSamples,
                numPartitionsPerNode,
                distFun,
                myResult,
                )
        end
        return myResult
    end
    function blockPrefixFun(x0,xs)
        return xs
    end
    result = mapPrefix(
        UInt32(1),
        UInt32(length(smallSubCoverSets)),
        UInt32(1),
        mapFun,
        blockPrefixFun,
        Vector{Tuple{CoverSet, Vector{UInt32}}}(),
    )
    return result
end

function WjfKNN(
    data::Vector,
    distFun::Function,
    maxNumSamples::UInt32 = UInt32(256),
    numPartitionsPerNode::UInt32 = UInt32(16),
)
    smallSubCoverSets = Vector{Tuple{CoverSet, Vector{UInt32}}}()
    entry = CoverSet(
        data,
        [i for i::UInt32 = 1:length(data)],
        maxNumSamples,
        numPartitionsPerNode,
        distFun,
        smallSubCoverSets,
    )
    while length(smallSubCoverSets) > 0
        smallSubCoverSets = finishSmallSubCoverSets(
            smallSubCoverSets,
            data,
            maxNumSamples,
            numPartitionsPerNode,
            distFun,
        )
    end
    return WjfKNN(
        data,
        distFun,
        maxNumSamples,
        numPartitionsPerNode,
        entry,
    )
end

function getNeighbors(coverSet::CoverSet, i::UInt32,
    bisearchThres::Float64)::Set{UInt32}
    if length(coverSet.subCoverSets) != 0
        (j1, j2, d1, d2) = coverSet.subCoverSetInfos[coverSet.mapping[i]]
        res::Set{UInt32} = getNeighbors(coverSet.subCoverSets[j1],
            i, bisearchThres)
        if d2 * bisearchThres <= d1
            res = union(res,
                getNeighbors(coverSet.subCoverSets[j2], i, bisearchThres))
        end
        return res
    else
        res = Set{UInt32}(keys(coverSet.mapping))
        return res
    end
end

function knnSearch(wjfknn::WjfKNN, i::UInt32, k::UInt32,
    bisearchThres::Float64)::Vector{Tuple{Float64, UInt32}}
    neighbors::Set{UInt32} = getNeighbors(wjfknn.entry, i, bisearchThres)
    results::Vector{Tuple{Float64, UInt32}} = Vector{Tuple{Float64, UInt32}}()
    for neighbor in neighbors
        if neighbor == i
            continue
        end
        dist = wjfknn.distFun(wjfknn.data[i], wjfknn.data[neighbor])
        result = (dist, neighbor)
        insert!(results, searchsortedfirst(results, result), result)
        if length(results) > k
            pop!(results)
        end
    end
    return results
end

function knnSearch(
    wjfknn::WjfKNN, k::UInt32, bisearchThres::Float64
    )::Vector{Vector{Tuple{Float64, UInt32}}}
    numData = length(wjfknn.data)
    result::Vector{Vector{Tuple{Float64, UInt32}}} =
        Vector{Vector{Tuple{Float64, UInt32}}}(undef, numData)
    function mapFun(u1::UInt32, u2::UInt32)
        for u::UInt32 = u1:u2
            result[u] = knnSearch(wjfknn, u, k, bisearchThres)
        end # u
    end
    mapOnly(UInt32(1), UInt32(numData), UInt32(100), mapFun)
    return result
end

end # module


text:
module distorted_block_infsup_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using FinEtools.MeshExportModule
using FinEtoolsDeforLinear.FEMMDeforLinearBaseModule: infsup_gh, infsup_sh
using Test
import LinearAlgebra: norm, cholesky, I, eigen
using UnicodePlots

function distorted_block_infsup_T10()
	lambdatol = sqrt(1e8*eps(1.0));
	E=1000.0;
	nu=0.24;
	parshiftmult= 0.002;
	A = [1.44 -0.741 -0.53; -0.626 1.589 -0.913; -0.55 0.43 1.756] + 1.0I;

	lambdamin = Float64[]
	h = Float64[]
	for ne = [2, 4, 6]
		Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt, orientation::Symbol = (6.0, 6.0, 6.0, ne, ne, ne, :a)

		fens, fes = T4block(Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt, orientation::Symbol)
		fens, fes = T4toT10(fens, fes)
		# @show connasarray(fes)

		for i = 1:count(fens)
			fens.xyz[i,:] = fens.xyz[i,:] + vec(reshape(fens.xyz[i,:], 1, 3)*A);
		end
		# @show fens.xyz

		# File =  "minfsuptest1.vtk"
		# vtkexportmesh(File, fens, fes)
		# try rm(File); catch end

		MR  =  DeforModelRed3D

		material = MatDeforElastIso(MR, E, nu)

		geom = NodalField(fens.xyz)
		u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field
		bfes = meshboundary(fes)
		l1 = connectednodes(bfes)
		setebc!(u, l1, true, 1, 0.0)
		setebc!(u, l1, true, 2, 0.0)
		setebc!(u, l1, true, 3, 0.0)
		numberdofs!(u)

		femm  =  FEMMDeforLinear(MR, IntegDomain(fes, TetRule(1)), material)
		Gh = infsup_gh(femm, geom, u);
		femm  =  FEMMDeforLinear(MR, IntegDomain(fes, TetRule(1)), material)
		Sh = infsup_sh(femm, geom, u);

		lambda, modes = eigen(Matrix(Gh), Matrix(Sh));

		# @show lambda
		abslambda = real.(filter(y -> !isnan(y), lambda));
		ix = findall(y  -> y < 0.0, abslambda);
		if !isempty(ix)
			abslambda[ix] .= 0;
		end

		abslambda = sqrt.(sort(abslambda));
		ix = findall(y  -> y > 0.0, abslambda);
		a = lineplot(1:length(abslambda[ix]), log.(abslambda[ix]), name = "infsup", xlabel = "eigenvalue", ylabel = "log(eigenvalue)", canvas = DotCanvas)
		display(a)

		ix = findall(y  -> y >= lambdatol, abslambda);
		if isempty(ix)
			@error "Bad guess of the number of eigenvalues"
		end
		push!(lambdamin, abslambda[ix[1]])
		push!(h, 1.0/(count(fens))^(1/3))
	end

	@show lambdamin
	a = lineplot(log.(h), log.(lambdamin), name = "infsup", xlabel = "log(Element Size)", ylabel = "log(minimum eigenvalue)", canvas = DotCanvas)
	display(a)

	# @test norm(lambdamin - [0.262065, 0.1709, 0.126159, 0.100228, 0.0828139]) / norm(lambdamin) <= 1.0e-4

end

function distorted_block_infsup_T4()
	lambdatol = sqrt(1e8*eps(1.0));
	E=1000.0;
	nu=0.24;
	parshiftmult= 0.002;
	A = [1.44 -0.741 -0.53; -0.626 1.589 -0.913; -0.55 0.43 1.756] + 1.0I;

	lambdamin = Float64[]
	h = Float64[]
	for ne = [2, 3, 4]
		Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt, orientation::Symbol = (6.0, 6.0, 6.0, ne, ne, ne, :a)

		fens, fes = T4block(Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt, orientation::Symbol)
		# fens, fes = T4toT10(fens, fes)
		# @show connasarray(fes)

		for i = 1:count(fens)
			fens.xyz[i,:] = fens.xyz[i,:] + vec(reshape(fens.xyz[i,:], 1, 3)*A);
		end
		# @show fens.xyz

		# File =  "minfsuptest1.vtk"
		# vtkexportmesh(File, fens, fes)
		# try rm(File); catch end

		MR  =  DeforModelRed3D

		material = MatDeforElastIso(MR, E, nu)

		geom = NodalField(fens.xyz)
		u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field
		bfes = meshboundary(fes)
		l1 = connectednodes(bfes)
		setebc!(u, l1, true, 1, 0.0)
		setebc!(u, l1, true, 2, 0.0)
		setebc!(u, l1, true, 3, 0.0)
		numberdofs!(u)

		femm  =  FEMMDeforLinear(MR, IntegDomain(fes, TetRule(1)), material)
		Gh = infsup_gh(femm, geom, u);
		femm  =  FEMMDeforLinear(MR, IntegDomain(fes, TetRule(1)), material)
		Sh = infsup_sh(femm, geom, u);

		lambda, modes = eigen(Matrix(Gh), Matrix(Sh));

		# @show lambda
		abslambda = real.(filter(y -> !isnan(y), lambda));
		ix = findall(y  -> y < 0.0, abslambda);
		if !isempty(ix)
			abslambda[ix] .= 0;
		end

		abslambda = sqrt.(sort(abslambda));
		ix = findall(y  -> y > 0.0, abslambda);
		a = lineplot(1:length(abslambda[ix]), log.(abslambda[ix]), name = "infsup", xlabel = "eigenvalue", ylabel = "log(eigenvalue)", canvas = DotCanvas)
		display(a)

		ix = findall(y  -> y >= lambdatol, abslambda);
		if isempty(ix)
			@error "Bad guess of the number of eigenvalues"
		end
		push!(lambdamin, abslambda[ix[1]])
		push!(h, 1.0/(count(fens))^(1/3))
	end

	@show lambdamin
	a = lineplot(log.(h), log.(lambdamin), name = "infsup", xlabel = "log(Element Size)", ylabel = "log(minimum eigenvalue)", canvas = DotCanvas)
	display(a)

	# @test norm(lambdamin - [0.262065, 0.1709, 0.126159, 0.100228, 0.0828139]) / norm(lambdamin) <= 1.0e-4

end

function distorted_block_infsup_H8()
	lambdatol = sqrt(1e8*eps(1.0));
	E=1000.0;
	nu=0.24;
	parshiftmult= 0.002;
	A = [1.44 -0.741 -0.53; -0.626 1.589 -0.913; -0.55 0.43 1.756] + 1.0I;

	lambdamin = Float64[]
	h = Float64[]
	for ne = [2, 3, 4]
		Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt, orientation::Symbol = (6.0, 6.0, 6.0, ne, ne, ne, :a)

		fens, fes = H8block(Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt)
		# fens, fes = T4toT10(fens, fes)
		# @show connasarray(fes)

		for i = 1:count(fens)
			fens.xyz[i,:] = fens.xyz[i,:] + vec(reshape(fens.xyz[i,:], 1, 3)*A);
		end
		# @show fens.xyz

		# File =  "minfsuptest1.vtk"
		# vtkexportmesh(File, fens, fes)
		# try rm(File); catch end

		MR  =  DeforModelRed3D

		material = MatDeforElastIso(MR, E, nu)

		geom = NodalField(fens.xyz)
		u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field
		bfes = meshboundary(fes)
		l1 = connectednodes(bfes)
		setebc!(u, l1, true, 1, 0.0)
		setebc!(u, l1, true, 2, 0.0)
		setebc!(u, l1, true, 3, 0.0)
		numberdofs!(u)

		femm  =  FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 1)), material)
		Gh = infsup_gh(femm, geom, u);
		femm  =  FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 1)), material)
		Sh = infsup_sh(femm, geom, u);

		lambda, modes = eigen(Matrix(Gh), Matrix(Sh));

		# @show lambda
		abslambda = real.(filter(y -> !isnan(y), lambda));
		ix = findall(y  -> y < 0.0, abslambda);
		if !isempty(ix)
			abslambda[ix] .= 0;
		end

		abslambda = sqrt.(sort(abslambda));
		ix = findall(y  -> y > 0.0, abslambda);
		a = lineplot(1:length(abslambda[ix]), log.(abslambda[ix]), name = "infsup", xlabel = "eigenvalue", ylabel = "log(eigenvalue)", canvas = DotCanvas)
		display(a)

		ix = findall(y  -> y >= lambdatol, abslambda);
		if isempty(ix)
			@error "Bad guess of the number of eigenvalues"
		end
		push!(lambdamin, abslambda[ix[1]])
		push!(h, 1.0/(count(fens))^(1/3))
	end

	@show lambdamin
	a = lineplot(log.(h), log.(lambdamin), name = "infsup", xlabel = "log(Element Size)", ylabel = "log(minimum eigenvalue)", canvas = DotCanvas)
	display(a)

	# @test norm(lambdamin - [0.262065, 0.1709, 0.126159, 0.100228, 0.0828139]) / norm(lambdamin) <= 1.0e-4

end

function distorted_block_infsup_H20()
	lambdatol = sqrt(1e8*eps(1.0));
	E=1000.0;
	nu=0.24;
	parshiftmult= 0.002;
	A = [1.44 -0.741 -0.53; -0.626 1.589 -0.913; -0.55 0.43 1.756] + 1.0I;

	lambdamin = Float64[]
	h = Float64[]
	for ne = [2, 3, 4]
		Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt, orientation::Symbol = (6.0, 6.0, 6.0, ne, ne, ne, :a)

		fens, fes = H20block(Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt)
		# fens, fes = T4toT10(fens, fes)
		# @show connasarray(fes)

		for i = 1:count(fens)
			fens.xyz[i,:] = fens.xyz[i,:] + vec(reshape(fens.xyz[i,:], 1, 3)*A);
		end
		# @show fens.xyz

		# File =  "minfsuptest1.vtk"
		# vtkexportmesh(File, fens, fes)
		# try rm(File); catch end

		MR  =  DeforModelRed3D

		material = MatDeforElastIso(MR, E, nu)

		geom = NodalField(fens.xyz)
		u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field
		bfes = meshboundary(fes)
		l1 = connectednodes(bfes)
		setebc!(u, l1, true, 1, 0.0)
		setebc!(u, l1, true, 2, 0.0)
		setebc!(u, l1, true, 3, 0.0)
		numberdofs!(u)

		femm  =  FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 2)), material)
		Gh = infsup_gh(femm, geom, u);
		femm  =  FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 2)), material)
		Sh = infsup_sh(femm, geom, u);

		lambda, modes = eigen(Matrix(Gh), Matrix(Sh));

		# @show lambda
		abslambda = real.(filter(y -> !isnan(y), lambda));
		ix = findall(y  -> y < 0.0, abslambda);
		if !isempty(ix)
			abslambda[ix] .= 0;
		end

		abslambda = sqrt.(sort(abslambda));
		ix = findall(y  -> y > 0.0, abslambda);
		a = lineplot(1:length(abslambda[ix]), log.(abslambda[ix]), name = "infsup", xlabel = "eigenvalue", ylabel = "log(eigenvalue)", canvas = DotCanvas)
		display(a)

		ix = findall(y  -> y >= lambdatol, abslambda);
		if isempty(ix)
			@error "Bad guess of the number of eigenvalues"
		end
		push!(lambdamin, abslambda[ix[1]])
		push!(h, 1.0/(count(fens))^(1/3))
	end

	@show lambdamin
	a = lineplot(log.(h), log.(lambdamin), name = "infsup", xlabel = "log(Element Size)", ylabel = "log(minimum eigenvalue)", canvas = DotCanvas)
	display(a)

	# @test norm(lambdamin - [0.262065, 0.1709, 0.126159, 0.100228, 0.0828139]) / norm(lambdamin) <= 1.0e-4

end

function allrun()
    println("#####################################################")
    println("# distorted_block_infsup_T4 ")
    distorted_block_infsup_T4()
    println("#####################################################")
    println("# distorted_block_infsup_H8 ")
    distorted_block_infsup_H8()
    println("#####################################################")
    println("# distorted_block_infsup_T10 ")
    distorted_block_infsup_T10()
    println("#####################################################")
    println("# distorted_block_infsup_H20 ")
    distorted_block_infsup_H20()
    return true
end # function allrun

end # module distorted_block_infsup_examples


text:
#=
Author: Stephan Scholz
Year: 2021

This file contains a Simulation of the Heat equation with a proportional controller
Boundary conditions:
- left side: heat induction (proportional control)
- right side: heat transfer and heat radiation (nonlinear)
Dimension: 1D
Spatial Approximation: Finite Differences
Temporal Approximation: Forward-Euler method
=#

using LinearAlgebra, Plots, SparseArrays, DifferentialEquations

# Define the constants for the PDE
const λ = 45.0          # Conductivity
const ρ = 7800.0;       # Density
const cap = 480.0;      # Specific heat capacitivity
const α = λ/(cap * ρ)   # Diffusivity
const L = 0.1 # 10*10^(-3)    # Length of rod

#Boundary conditions
const h = 10        # Heat transfer coefficient
const ϵ = 0.6;      # Emissivity
const sb = 5.67*10^(-8) # Stefan-Boltzmann constant
const k = ϵ * sb # Radiation coefficient

# Discretization

const N = 101;      # Number of elements
const Δx = L/(N-1)  # Spatial sampling
const Tf = 3000.0 # 15.0;
const dt = 10^(-2) # 1.0# 10^(-4); # Time step width
# o--o--o--o--o--o

#Laplace operator
const M = spdiagm(-1 => ones(N-1), 0 => -2*ones(N), 1 => ones(N-1))
M[1,2] = 2
M[end,end-1] = 2

const Mx = zeros(N)

# x'(t) = A x(t) + B u(t) + E w(t)
const A = α/(Δx^2) * M
const b = 1
const B = spzeros(N); # Input vector
const B[1] = 2 * (α/(λ*Δx)) * b
const E = spzeros(N); # Disturbance vector
const E[end] = 2 * (α/(λ*Δx)) * b

const q = α*dt/(Δx^2);
if q > 0.5
  error("Numerical stability is NOT guaranteed! STOP program!")
end

# Temperatures
const θ₀   = 300.0  # Initial temperature
const θamb = 298.0  # Ambient temperature
const θref = 400.0  # Reference temperature

# Flux towards environment
ϕout(θ) = -h * (θ - θamb)  - k*(θ^4 - θamb^4)

### Controller ###
const Kp = 10^3 # Proportional gain
const yref = 400.0 # Reference temperature

u_in(err) = Kp * err # input signals

global u_hist
u_hist = zeros(1,0) # History of input signals

function heat_eq!(dx,x,p,t)
    global u_hist

    err = yref - x[end]
    u = u_in(err)

    u_hist = hcat(u_hist, u)

    dx .= A * x + B * u + E * ϕout( x[end] )    # Integration of temperatures

  end

x₀ = θ₀ * ones(N)
tspan = (0.0, Tf)

prob = ODEProblem( heat_eq!, x₀, tspan ) # ODE Problem
sol = solve(prob,Euler(),dt=dt ,progress=true, saveat=1.0, save_start=true);

# 1-dimensional grid
xspan = 0 : Δx : L

using Plots
# 2D Heatmap plot: xaxis = time ; yaxis = position
heatmap(sol.t, xspan, sol[1:end,1:end], xaxis="Time [s]", yaxis="Position [m]", title="Evolution of temperature")

# Temperature at left and right side
# plot(sol.t,[sol[1,:], sol[end-1,:]], label=["Left" "Right"], title="Temperature at the left/right end", legend=:bottomright)


text:
module MixingRules

export bruggemanspheres, looyengacylinders, looyengaspheres, lorentzlorenz, maxwellgarnettspheres, monecke, gedf, gem

"""
Calculates the complex refractive index of a porous material using Bruggeman model, for nonmagnetic and isotropic materials. Source: PHYSICAL REVIEW B VOLUME 61, NUMBER 15 15 APRIL 2000-I. It considers spheres (n1) hosted in another material (n2).
Usage: neff = bruggemanspheres(n1,n2,p)
Input:
       n1:    refractive index of material 1
       n2:    refractive index of material 2
       p:     porosity parameter (real number, 0<p<1), proportion of n1 over total
Output:
       neff:  effective refractive index
"""
function bruggemanspheres(n1::AbstractArray{T,N}, n2::AbstractArray{T,N}, p::S) where {T<:Number, N, S<:Float64}
    # dielectric function of each media
    df1 = n1.^2
    df2 = n2.^2
    # effective dielectric function: PHYSICAL REVIEW B VOLUME 61, NUMBER 15 15 APRIL 2000-I. solved with mathematica for dfeff.
    dfeff = 0.25 .* (-df1 + 2 .* df2 + 3 .* df1.*p - 3 .* df2 .* p + sqrt.( 8 .* df1 .* df2 + (-df1 + 2 .* df2 + 3 .* df1.*p - 3 .* df2 .* p).^2 ) )
    # compute effective refractive index
    neff = dfeff.^0.5
    return neff
end # EOF bruggemanspheres(...)

"""
Calculates the complex refractive index of the porous material using Looyenga-Landau-Lifshitz model, for nonmagnetic and isotropic materials. Source: Langmuir 2013, 29, 2784−2789. It considers cylinders (n1) hosted in another material (n2).
Usage:
       neff = looyengacylinders(n1,n2,p)
Input:
       n1:    refractive index of material 1
       n2:    refractive index of material 2
       p:     porosity parameter (real number, 0<p<1), proportion of n1 over total
Output:
       neff:  effective refractive index
"""
function looyengacylinders(n1::AbstractArray{T,N}, n2::AbstractArray{T,N}, p::S) where {T<:Number, N, S<:Float64}
    # compute effective reractive index
    neff = (1-p) .* n2.^2  + n1.^2 .* p
end # function looyengacylinders(...)

"""
Calculates the complex refractive index of the porous material using Looyenga-Landau-Lifshitz model, for nonmagnetic and isotropic materials. Source: IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 38, NO. 3, MAY 2000. It considers spheres (n1) hosted in another material (n2).
Usage:
       neff = looyengaspheres(n1,n2,p)
Input:
       n1:    refractive index of material 1
       n2:    refractive index of material 2
       p:     porosity parameter (real number, 0<p<1), proportion of n1 over total
Output:
       neff:  effective refractive index
"""
function looyengaspheres(n1::AbstractArray{T,N}, n2::AbstractArray{T,N}, p::S) where {T<:Number, N, S<:Float64}
    # dielectric function of each media
    df1 = n1.^2
    df2 = n2.^2
    # compute effective reractive index
    neff = ( ( ( (1-p) .* (df2.^(1/3)) ) + ( (df1.^(1/3)) .* p) ).^3 ).^0.5
    return neff
end # function looyengaspheres(...)

"""
Returns the effective dielectric function of a binary liquid mixture. p is the proportion of component n1.
Usage:
       neff = lorentzlorenz(n1,n2,p)
Input:
       n1:    refractive index of material 1
       n2:    refractive index of material 2
       p:     proportion of n1 over total, parameter (real number, 0<p<1)
Output:
       neff:  effective refractive index
"""
function lorentzlorenz(n1::AbstractArray{T,N}, n2::AbstractArray{T,N}, p::S) where {T<:Number, N, S<:Float64}
    neff = sqrt.( ( n2.^2 .* ( -n1.^2 + 2 * p^2 - 2 - 2 .* p .* n1^2) ) ./ ( n1.^2 .* (p-1) - 2 - p .* n2.^2 )  )
end

"""
Calculates the complex refractive index of the porous material using Maxwell-Garnett model, for nonmagnetic and isotropic materials. Source: IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 38, NO. 3, MAY 2000. It considers spheres (n1) hosted in another material (n2).
Usage:
       neff = maxwellgarnettspheres(n1,n2,p)
Input:
       n1:    refractive index of material 1
       n2:    refractive index of material 2
       p:     porosity parameter (real number, 0<p<1), proportion of n1 over total
Output:
       neff:  effective refractive index
"""
function maxwellgarnettspheres(n1::AbstractArray{T,N}, n2::AbstractArray{T,N}, p::S) where {T<:Number, N, S<:Float64}
    # dielectric function of each media
    df1 = n1.^2
    df2 = n2.^2
    # compute effective reractive index
    nfeff = ( df1 + ( (2 .* p .* df1) .* (df2-df1) ./ ( df2 + df1 - p .* (df2-df1) ) ) ).^0.5
    return neff
end # EOF maxwellgarnettspheres(n1,n2,p)

"""
Calculates the complex refractive index of the porous material using Monecke model, for nonmagnetic and isotropic materials. Source: Phys. Rev. B, Vol. 61, Num. 15, 15 April, 2000-I.
Usage:
        neff = monecke(n1, n2, p)
Input:
        n1:    refractive index of material 1
        n2:    refractive index of material 2
        p:     porosity parameter (real number, 0<p<1), proportion of n1 over total
Output:
        neff:  effective refractive index
"""
function monecke(n1::AbstractArray{T,N}, n2::AbstractArray{T,N}, p::S) where {T<:Number, N, S<:Float64}
    # dielectric function of each media
    df1 = n1.^2
    df2 = n2.^2
    # compute effective dielectric function
    dfeff = ( 2 * ( p * df1 + (1-p) .* df2 ).^2 ) + ( df1 .* df2) ./ ( ( 1 + p) .* df1 ) + ( (2 - p) .* df2 )
    # compute effective reractive index
    neff = dfeff.^0.5
    return neff
end # EOF monecke(...)

"""
Calculates the complex refractive index of a porous material using a general dielectric model, for nonmagnetic and isotropic materials. Source: Appl. Phys. Lett. 108, 102902 (2016); doi: 10.1063/1.4943639. It considers n1) hosted in another material (n2).
Usage: neff = gedf(n1, n2, p, β)
Input:
       n1:    refractive index of filling material
       n2:    refractive index of host material
       p:     porosity parameter (real number, 0<p<1), proportion of n1 over total
       β:     either 1 or 3
Output:
       neff:  effective refractive index
"""
function gedf(n1::AbstractArray{T,N}, n2::AbstractArray{T,N}, p::S, β::S) where {T<:Number, N, S<:Number}
    # dielectric function of each media
    df1 = n1.^2
    df2 = n2.^2
    deff = df2 .* ( df1 + (β * df2) + (β * p * df1) - (β * p * df2) ) ./ ( df1 + (β * df2) - (p * df1) + (p * df2) )
    # compute effective refractive index
    neff = dfeff.^0.5
    return neff
end # EOF gedf(...)

"""
Source: Carbon 40 (2002) 2801–2815. Here we used equation (3) and the parameters ϕ, ϕc and tp as fitted. For a more detailed analysis and relation to the depolarization factors see the reference. It considers (n1) hosted in another material (n2).
Usage: neff = gedf(n1, n2, ϕ, ϕc, tp)
Input:
    n1:    refractive index of filling material (lower index)
    n2:    refractive index of host material (higher index)
    ϕ:     porosity parameter (real number, 0<ϕ<1), proportion of n1 over total
    ϕc:    critical volume fraction at which the high-index phase first percolates
    tp:    exponent related both to ϕc and to the shape of the grains
Output:
    neff:  effective refractive index
"""
function gem(n1::AbstractArray{T,N}, n2::AbstractArray{T,N}, ϕ::S, ϕc::S, tp::S) where {T<:Number, N, S<:Number}
    # dielectric function of each media
    α = (1-ϕc)/ϕc
    β = 1/tp
    df1 = (n1.^2).^β
    df2 = (n2.^2).^β
    deff = ((α*df1 - ϕ*df1 - α*ϕ*df1 - df2 + ϕ*df2 + α*ϕ*df2 + sqrt.(4*α*df1*df2 + ((α*(ϕ-1) + ϕ)*df1 + df2 - (1 + α)*ϕ*df2).^2))/(2*α)).^tp
    # compute effective refractive index
    neff = dfeff.^0.5
    return neff
end # EOF gem(...)

end # module MixingRules


text:
"""
帯グラフは，グループごとの相対度数分布を積み上げ棒グラフにしたものなので，以下のような関数を作っておけばよい。
"""

function obi(table; title="", ylabel="", barwidth=0.9, fontsize=8, width=600, height=300)
    g, n = size(table)
    table2 = table ./ sum(table, dims=2)
    row_names, col_names = names(table)
    p2 = groupedbar(table2[:, end:-1:1],
        bar_width=barwidth,
        direction=:horizontal,
        bar_position = :stack,
        yticks = (1:length(row_names), row_names),
        xticks = (0:0.1:1, string.(0:10:100) .* "%"),
        showaxis=:x,
        title = title,
        ylabel = ylabel,
        size=(width, height)
        )
    for row in 1:g
        cum = cumsum(table2[row, :])
        cum0 = vcat(0, cum[1:end-1]...)
        for col in n:-1:1
            annotate!((cum0[col] + cum[col])/2, row,
                text(col_names[col] * "\n" * string(table[row, col]), fontsize))
        end
    end
    p2
end

"""
使用法は，freqtable() の第 1 引数にグループ分けの変数，第 2 引数に分布を表示する変数として得た二重クロス表を barplot2() に与える。

julia> using RDatasets

julia> haireyecolor = RDatasets.dataset("datasets", "HairEyeColor");

julia> haireyecolor[!, 1:3] = string.(haireyecolor[!, 1:3]); # String7 が使えないパッケージに対する過渡的な対処法

julia> HEC = DataFrame();

julia> for row in eachrow(haireyecolor), j in 1:row.Freq
           push!(HEC, row)
       end

julia> select!(HEC, Not(:Freq)); # Freq 列を除く

julia> using FreqTables, DataFrames, Plots, StatsPlots

julia> freq2 = freqtable(HEC.Hair, HEC.Sex)
4×2 Named Matrix{Int64}
Dim1 ╲ Dim2 │ Female    Male
────────────┼───────────────
Black       │     52      56
Blond       │     81      46
Brown       │    143     143
Red         │     37      34

julia> obi(freq2, ylabel="Hair", title="Distribution of hair color in students")

"""

max_stars_repo_name: r-de-r/stats

max_issues_repo_name: r-de-r/stats

max_forks_repo_name: r-de-r/stats


text:
function GAIL_D_loss(gan_loss)
    (D, 𝒫, 𝒟_ex, 𝒟_π; info = Dict()) ->begin
        Lᴰ(gan_loss, D, 𝒟_ex[:a], 𝒟_π[:a], yD = (𝒟_ex[:s],), yG = (𝒟_π[:s],))
    end
end

function OnPolicyGAIL(;π,
                       S,
                       γ,
                       λ_gae::Float32 = 0.95f0,
                       𝒟_demo,
                       αr::Float32 = 0.5f0,
                       normalize_demo::Bool=true,
                       D::ContinuousNetwork,
                       solver=PPO,
                       gan_loss::GANLoss=GAN_BCELoss(),
                       d_opt::NamedTuple=(;),
                       log::NamedTuple=(;),
                       kwargs...)

    d_opt = TrainingParams(;loss = GAIL_D_loss(gan_loss), name="discriminator_", d_opt...)
    normalize_demo && (𝒟_demo = normalize!(deepcopy(𝒟_demo), S, action_space(π)))
    𝒟_demo = 𝒟_demo |> device(π)

    function GAIL_callback(𝒟; info=Dict(), 𝒮)
        batch_train!(D, d_opt, (;), 𝒟_demo, deepcopy(𝒟), info=info)

        D_out = value(D, 𝒟[:a], 𝒟[:s]) # This is swapped because a->x and s->y and the convention for GANs is D(x,y)
        r = αr * logσ.(D_out) .- (1f0 - αr) * logcompσ.(D_out)
        ignore() do
            info["disc_reward"] = mean(r)
        end

        𝒟[:r] .= r

        eps = episodes(𝒟)
        for ep in eps
            eprange = ep[1]:ep[2]
            fill_gae!(𝒟, eprange, 𝒮.agent.π, λ_gae, γ)
            fill_returns!(𝒟, eprange, γ)
        end
        𝒟[:advantage] .= whiten(𝒟[:advantage])

    end
    solver(;π=π, S=S, post_batch_callback=GAIL_callback, log=(dir="log/onpolicygail", period=500, log...), λ_gae=λ_gae, kwargs...)
end


text:
using Test

using Seis
using Beamforming, Beamforming.Synth

@testset "Vespagrams" begin
    let slow = 10rand(), baz = 360rand(), maxima = rand(1:5), Δ = 100,
            wavefront = :plane, envelope = false
        sx, sy = slow.*(sind(baz), cosd(baz))
        s = Synth.synthetic_arrival(sx, sy, distance=Δ, noise=0.001)
        nsta = length(s)
        for (method, n) in zip((:linear, :nthroot, :phaseweight), (nothing, 3, 3))
            vespa = vespagram(s, 40, 60, 0, 10, 0.1, maxima=maxima, method=method, n=n)
            @test vespa isa Beamforming.VespaGrid
            @test length(vespa.x) == length(vespa.y) == nsta
            @test vespa.slow ≈ collect(0:0.1:10)
            @test vespa.time ≈ collect(40:s[1].delta:60)
            @test size(vespa.power) == (length(vespa.time), length(vespa.slow))
            @test vespa.nmax == maxima
            @test length(vespa.t_max) == length(vespa.s_max) == vespa.nmax
            @test vespa.method == method
            @test vespa.n == n
            @test vespa.envelope == envelope
            @test vespa.wavefront == wavefront
            @test vespa.s_max[1] ≈ slow atol=0.2
            @test vespa.t_max[1] ≈ 50 atol=15s[1].delta
        end
    end
end

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming


text:
struct VarBCRecord
    pdate::String
    class::String
    # key::String
    label::String
    ndata::Int64
    npred::Int64
    predcs::Vector{Int64}
    param0::Vector{Float64}
    params::Vector{Float64}
    hstgrm::Vector{Int64}
    predxcnt::Vector{Int64}
    predmean::Vector{Float64}
    predxcov::Vector{Float64}
 end

 show(io::IO,a::VarBCRecord) = print(io,a.label)

 function show(io::IO,::MIME"text/plain", a::VarBCRecord)
    println(io, "VarBC Record")
    println(io, "label  = $(a.label)")
    println(io, "predcs = $(a.predcs)")
    println(io, "ndata  = $(a.ndata)")
    println(io, "param0 = $(a.param0)")
    println(io, "params = $(a.params)")
 end

 function ==(a::VarBCRecord,b::VarBCRecord)
     a.pdate == b.pdate &&
     a.class == b.class &&
     a.label == b.label &&
     a.ndata == b.ndata &&
     a.npred == b.npred &&
     b.predcs == b.predcs &&
     a.param0 == b.param0 &&
     a.params == b.params &&
     a.hstgrm == b.hstgrm &&
     a.predxcnt == b.predxcnt &&
     a.predmean == b.predmean &&
     a.predxcov == a.predxcov
 end

 # increment(a::VarBCRecord) = a.params - a.param0
 # label(a::VarBCRecord) = a.label

 function read(io::IO,::Type{VarBCRecord})
       rmdi2nan(x) = x == -2.147e+09 ? NaN : x
       parsev(::Type{T},x) where T = split(x) |>  x -> parse.(T,x)  |> x -> rmdi2nan.(x)

       ix       = readline(io)  # |> x -> replace(x, r"^ix="       => "")  |> x -> parse(Int,x)   # ix not used
       pdate    = readline(io) |> x -> replace(x, r"^pdate="     => "")  # |> x -> parse(String,x)
       class    = readline(io) |> x -> replace(x, r"^class="    => "")  #|> x -> parse(String,x)
       key      = readline(io) |> x -> replace(x, r"^key="      => "")  #|> x -> parse(String,x)
       label    = readline(io) |> x -> replace(x, r"^label="    => "")  #|> x -> parse(String,x)
       ndata    = readline(io) |> x -> replace(x, r"^ndata="    => "")  |> x -> parse(Int64,x)
       npred    = readline(io) |> x -> replace(x, r"^npred="    => "")  |> x -> parse(Int64,x)
       predcs   = readline(io) |> x -> replace(x, r"^predcs="   => "")  |> x -> parsev(Int64,x)
       param0   = readline(io) |> x -> replace(x, r"^param0="   => "")  |> x -> parsev(Float64,x)
       params   = readline(io) |> x -> replace(x, r"^params="   => "")  |> x -> parsev(Float64,x)
       hstgrm   = readline(io) |> x -> replace(x, r"^hstgrm="   => "")  |> x -> parsev(Int64,x)
       predxcnt = readline(io) |> x -> replace(x, r"^predxcnt=" => "")  |> x -> parsev(Int64,x)
       predmean = readline(io) |> x -> replace(x, r"^predmean=" => "")  |> x -> parsev(Float64,x)
       predxcov = readline(io) |> x -> replace(x, r"^predxcov=" => "")  |> x -> parsev(Float64,x)
       key => VarBCRecord(pdate,class,label,ndata,npred,predcs,param0,params,hstgrm,predxcnt,predmean,predxcov)
   end

function write(io::IO,a::VarBCRecord, key::String, ix::Int)

   # Print Float  Int  Array{Float,1} etc.
   myprint(f::Float64) = isnan(f) ? "-2.147E+09" : @sprintf("% -.3E",f)
   myprint(i::Int64) = string(i)
   myprint(s::String) = s
   myprint(a::Array{T,1}) where T = join(myprint.(a)," ")

   println(io,"ix=$ix")
   println(io,"pdate=$(myprint(a.pdate))")
   println(io,"class=$(myprint(a.class))")
   println(io,"key=$key")
   println(io,"label=$(myprint(a.label))")
   println(io,"ndata=$(myprint(a.ndata))")
   println(io,"npred=$(myprint(a.npred))")
   println(io,"predcs=$(myprint(a.predcs))")
   println(io,"param0=$(myprint(a.param0))")
   println(io,"params=$(myprint(a.params))")
   println(io,"hstgrm=$(myprint(a.hstgrm))")
   println(io,"predxcnt=$(myprint(a.predxcnt))")
   println(io,"predmean=$(myprint(a.predmean))")
   println(io,"predxcov=$(myprint(a.predxcov))")
end


text:
module TestWDICountries

using Test
using WorldBankData

country_data = Any[ Dict("total"=>4,"per_page"=>"25000","pages"=>1,"page"=>1),
                 Any[
                    Dict("latitude"=>"42.5075","id"=>"AND","iso2Code"=>"AD","incomeLevel"=>Dict("id"=>"NOC","value"=>"High income: nonOECD"),"adminregion"=>Dict("id"=>"","value"=>""),"lendingType"=>Dict("id"=>"LNX","value"=>"Not classified"),"region"=>Dict("id"=>"ECS","value"=>"Europe & Central Asia (all income levels)"),"capitalCity"=>"Andorra la Vella","name"=>"Andorra","longitude"=>"1.5218"),
                    Dict("latitude"=>"","id"=>"ARB","iso2Code"=>"1A","incomeLevel"=>Dict("id"=>"NA","value"=>"Aggregates"),"adminregion"=>Dict("id"=>"","value"=>""),"lendingType"=>Dict("id"=>"","value"=>"Aggregates"),"region"=>Dict("id"=>"NA","value"=>"Aggregates"),"capitalCity"=>"","name"=>"Arab World","longitude"=>""),
                    Dict("latitude"=>"24.4764","id"=>"ARE","iso2Code"=>"AE","incomeLevel"=>Dict("id"=>"NOC","value"=>"High income: nonOECD"),"adminregion"=>Dict("id"=>"","value"=>""),"lendingType"=>Dict("id"=>"LNX","value"=>"Not classified"),"region"=>Dict("id"=>"MEA","value"=>"Middle East & North Africa (all income levels)"),"capitalCity"=>"Abu Dhabi","name"=>"United Arab Emirates","longitude"=>"54.3705"),
                    Dict("latitude"=>"-34.6118","id"=>"ARG","iso2Code"=>"AR","incomeLevel"=>Dict("id"=>"UMC","value"=>"Upper middle income"),"adminregion"=>Dict("id"=>"LAC","value"=>"Latin America & Caribbean (developing only)"),"lendingType"=>Dict("id"=>"IBD","value"=>"IBRD"),"region"=>Dict("id"=>"LCN","value"=>"Latin America & Caribbean (all income levels)"),"capitalCity"=>"Buenos Aires","name"=>"Argentina","longitude"=>"-58.4173"),
                 ]
               ]

df_country = WorldBankData.parse_country(country_data)

@test df_country[:name] == String["Andorra", "Arab World", "United Arab Emirates", "Argentina"]
@test df_country[:iso2c] == String["AD", "1A", "AE", "AR"]

end


text:
using DemoCards: parse_markdown

@testset "regexes" begin
    root = joinpath("assets", "regexes")

    @testset "title" begin
        ref_config_1 = Dict("title" => "This is a title",
                        "id" => "This-is-a-title-1")
        ref_config_2 = Dict("title" => "This is a title",
                        "id" => "custom_id")

        config = parse_markdown(joinpath(root, "title_1.md"))
        @test config == ref_config_1

        config = parse_markdown(joinpath(root, "title_2.md"))
        @test config == ref_config_2

        config = parse_markdown(joinpath(root, "title_3.md"))
        @test config == ref_config_2

        config = parse_markdown(joinpath(root, "title_4.md"))
        @test config == ref_config_1
    end
end


text:
workspace()

using SymPy

include("js_code.jl")

# Define free variable symbols
@syms ln2, l0, l1, apexDistance, ascensionFactor, descensionFactor, panAngle, windowSize, desiredPixelFlow, tanHalfFov, averageLookAtRadius

# Zoom-pan ratio used during ascension/descension
lFp(l0, l1, sl, k) = ((l1 / l0) - 1) / (-(1 / (sl * windowSize) * ln2 * panAngle) / (l0 * tanHalfFov / (windowSize * averageLookAtRadius)) * k)

# eye/lookAt distance while zoom/panning
fL(l0, l1, sl, k) = lFp(l0, l1, sl, k) / (1 + lFp(l0, l1, sl, k)) * desiredPixelFlow

# Ascension/descension time
tADpl(l0, l1, sl, k) = (sl * (log(l0 / l1) / ln2) * windowSize) / fL(l0, l1, sl, k)

# Ascension/descension time in simplified case without panning
tADl(l0, l1, sl, k) = (sl * (log(l0 / l1) / ln2) * windowSize) / desiredPixelFlow

# Panion time
tPp(l, k) = (panAngle * k * averageLookAtRadius * windowSize / (tanHalfFov * l)) / desiredPixelFlow

# Specific time functions for all cases
tApl(h) = simplify(tADpl(l0, h, -1, ascensionFactor))
tAl(h) = simplify(tADl(l0, h, -1, ascensionFactor))
tP(h) = simplify(tPp(h, 1 - ascensionFactor - descensionFactor))
tDpl(h) = simplify(tADpl(h, l1, 1, descensionFactor))
tDl(h) = simplify(tADl(h, l1, 1, descensionFactor))

tBaseLine(l0, l1) = simplify(tADpl(l0, l1, Piecewise((-1, LessThan(l0, l1)), (1, true)), 1))

timefuncs = ((
  ("AscensionZoomPan", tApl),
  ("AscensionZoomOnly", tAl),
  ("Panion", tP),
  ("DescensionZoomPan", tDpl),
  ("DescensionZoomOnly", tDl))

#
# # Derivatives
# dtApl(h) = simplify(diff(tApl(h), h))
# dtAl(h) = simplify(diff(tAl(h), h))
# dtP(h) = simplify(diff(tP(h), h))
# dtDpl(h) = simplify(diff(tDpl(h), h))
# dtDl(h) = simplify(diff(tDl(h), h))
#
# # Second Derivatives
# ddtApl(h) = simplify(diff(tApl(h), h, h))
# ddtAl(h) = simplify(diff(tAl(h), h, h))
# ddtP(h) = simplify(diff(tP(h), h, h))
# ddtDpl(h) = simplify(diff(tDpl(h), h, h))
# ddtDl(h) = simplify(diff(tDl(h), h, h))

function outputFunctions()
  funcbodies = []
  paramvars = []

  for f in timefuncs
    push!(funcbodies, (string("t", f[1]), js_code(cse(f[2](apexDistance)), paramvars)))
    push!(funcbodies, (string("dt", f[1]), js_code(cse(diff(f[2](apexDistance), apexDistance)), paramvars)))
    push!(funcbodies, (string("ddt", f[1]), js_code(cse(diff(f[2](apexDistance), apexDistance, apexDistance)), paramvars)))
  end

  blinecode = js_code(cse(tBaseLine(l0, l1)), paramvars)

  ret = []

  # Output params interface
  sort!(paramvars)

  push!(ret, "import { PathDefinition } from \"./interfaces\";\n")

  push!(ret, "export interface Params {")
  push!(ret, "  ascensionFactor: number;")
  push!(ret, "  descensionFactor: number;")
  push!(ret, "}\n")

  for f in funcbodies
    push!(ret, string("export function ", f[1], "(d: PathDefinition, p: Params, apexDistance: number): number {"))
    push!(ret, f[2])
    push!(ret, "}\n")
  end

  push!(ret, string("export function tBaseLine(d: PathDefinition): number {"))
  push!(ret, blinecode)
  push!(ret, "}\n")

  display(join(ret, "\n"))
end

outputFunctions()

max_stars_repo_name: rihorn/ky-libaries

max_issues_repo_name: rihorn/ky-libaries

max_forks_repo_name: rihorn/ky-libaries


text:
# global cache
# ============

# XXX `JET_REPORT_CACHE` isn't synced with `JET_CODE_CACHE`, and this may lead to a problem ?

"""
    JET_REPORT_CACHE::$(typeof(JET_REPORT_CACHE))

Keeps JET report cache for a `MethodInstance`.
Reports are cached when `JETInterpreter` exits from `_typeinf`.
"""
const JET_REPORT_CACHE = IdDict{UInt, IdDict{MethodInstance,Vector{InferenceErrorReportCache}}}()

"""
    JET_CODE_CACHE::$(typeof(JET_CODE_CACHE))

Keeps `CodeInstance` cache associated with `mi::MethodInstace` that represent the result of
  an inference on `mi` performed by `JETInterpreter`.
This cache is completely separated from the `NativeInterpreter`'s global cache, so that
  JET analysis never interacts with actual code execution.
"""
const JET_CODE_CACHE = IdDict{UInt, IdDict{MethodInstance,CodeInstance}}()

function CC.code_cache(interp::JETInterpreter)
    cache  = JETGlobalCache(interp)
    worlds = WorldRange(get_world_counter(interp))
    return WorldView(cache, worlds)
end

struct JETGlobalCache
    interp::JETInterpreter
end

# cache existence for this `interp` is ensured on its construction
jet_report_cache(interp::JETInterpreter)         = JET_REPORT_CACHE[cache_key(interp)]
jet_report_cache(wvc::WorldView{JETGlobalCache}) = jet_report_cache(wvc.cache.interp)
jet_code_cache(interp::JETInterpreter)           = JET_CODE_CACHE[cache_key(interp)]
jet_code_cache(wvc::WorldView{JETGlobalCache})   = jet_code_cache(wvc.cache.interp)

CC.haskey(wvc::WorldView{JETGlobalCache}, mi::MethodInstance) = haskey(jet_code_cache(wvc), mi)

function CC.get(wvc::WorldView{JETGlobalCache}, mi::MethodInstance, default)
    # ignore code cache for a `MethodInstance` that is not yet analyzed by JET;
    ret = get(jet_code_cache(wvc), mi, default)
    if isa(ret, CodeInstance)
        # cache hit, now we need to append cached reports associated with this `MethodInstance`
        global_cache = get(jet_report_cache(wvc), mi, nothing)
        if isa(global_cache, Vector{InferenceErrorReportCache})
            interp = wvc.cache.interp
            for cached in global_cache
                restored = restore_cached_report!(cached, interp)
                push!(interp.to_be_updated, restored) # should be updated in `abstract_call` (after exiting `typeinf_edge`)
                # # TODO make this hold
                # @assert first(cached.st).linfo === mi "invalid global restoring"
            end
        end
    end
    return ret
end

function CC.getindex(wvc::WorldView{JETGlobalCache}, mi::MethodInstance)
    r = CC.get(wvc, mi, nothing)
    r === nothing && throw(KeyError(mi))
    return r::CodeInstance
end

function CC.setindex!(wvc::WorldView{JETGlobalCache}, ci::CodeInstance, mi::MethodInstance)
    setindex!(jet_code_cache(wvc), ci, mi)
    add_jet_callback!(mi) # register the callback on invalidation
    return nothing
end

function add_jet_callback!(linfo)
    @static BACKEDGE_CALLBACK_ENABLED || return nothing

    if !isdefined(linfo, :callbacks)
        linfo.callbacks = Any[invalidate_jet_cache!]
    else
        if !any(@nospecialize(cb)->cb===invalidate_jet_cache!, linfo.callbacks)
            push!(linfo.callbacks, invalidate_jet_cache!)
        end
    end
    return nothing
end

function invalidate_jet_cache!(replaced, max_world, depth = 0)
    for cache in values(JET_REPORT_CACHE); delete!(cache, replaced); end
    for cache in values(JET_CODE_CACHE); delete!(cache, replaced); end

    if isdefined(replaced, :backedges)
        for mi in replaced.backedges
            mi = mi::MethodInstance
            if !any(cache->haskey(cache, mi), values(JET_REPORT_CACHE)) && !any(cache->haskey(cache, mi), values(JET_CODE_CACHE))
                continue # otherwise infinite loop)
            end
            invalidate_jet_cache!(mi, max_world, depth+1)
        end
    end
    return nothing
end

# local cache
# ===========

struct JETLocalCache
    interp::JETInterpreter
    cache::Vector{InferenceResult}
end

CC.get_inference_cache(interp::JETInterpreter) = JETLocalCache(interp, get_inference_cache(interp.native))

function CC.cache_lookup(linfo::MethodInstance, given_argtypes::Vector{Any}, cache::JETLocalCache)
    inf_result = cache_lookup(linfo, given_argtypes, cache.cache)

    isa(inf_result, InferenceResult) || return nothing

    # constant prop' hits a cycle (recur into same non-constant analysis), we just bail out
    isa(inf_result.result, InferenceState) && return inf_result

    # cache hit, try to restore local report caches
    interp = cache.interp
    sv = interp.current_frame::InferenceState

    analysis_result = jet_cache_lookup(linfo, given_argtypes, interp.cache)

    # corresponds to report throw away logic in `_typeinf(interp::JETInterpreter, frame::InferenceState)`
    filter!(!is_from_same_frame(sv.linfo, linfo), interp.reports)

    if isa(analysis_result, AnalysisResult)
        for cached in analysis_result.cache
            restored = restore_cached_report!(cached, interp)
            push!(interp.to_be_updated, restored) # should be updated in `abstract_call_method_with_const_args`
            # # TODO make this hold
            # @assert first(cached.st).linfo === linfo "invalid local restoring"
        end
    end

    return inf_result
end

# should sync with the implementation of `cache_lookup(linfo::MethodInstance, given_argtypes::Vector{Any}, cache::Vector{InferenceResult})`
function jet_cache_lookup(linfo::MethodInstance, given_argtypes::Vector{Any}, cache::Vector{AnalysisResult})
    method = linfo.def::Method
    nargs::Int = method.nargs
    method.isva && (nargs -= 1)
    length(given_argtypes) >= nargs || return nothing
    for cached_result in cache
        cached_result.linfo === linfo || continue
        cache_match = true
        cache_argtypes = cached_result.argtypes
        cache_overridden_by_const = cached_result.overridden_by_const
        for i in 1:nargs
            if !is_argtype_match(given_argtypes[i],
                                 cache_argtypes[i],
                                 CC.getindex(cache_overridden_by_const, i))
                cache_match = false
                break
            end
        end
        if method.isva && cache_match
            cache_match = is_argtype_match(tuple_tfunc(given_argtypes[(nargs + 1):end]),
                                           cache_argtypes[end],
                                           CC.getindex(cache_overridden_by_const, CC.length(cache_overridden_by_const)))
        end
        cache_match || continue
        return cached_result
    end
    return nothing
end

CC.push!(cache::JETLocalCache, inf_result::InferenceResult) = CC.push!(cache.cache, inf_result)


text:
"""
    function SPR(original_root<:AbstractNode)

Performs SPR on tree. Takes a copy of root of the tree;
Returns a copy of root of altered tree. Throws error if tree is improperly formatted.
"""
function SPR(original_root::T)::T where T <: AbstractNode
    root = deepcopy(original_root)
    SPR!(root)
    return root
end #func

"""
    SPR!(root::T)::T where T <:AbstractNode

Performs SPR on tree in place. Takes reference to root of tree;
Returns reference to root of altered tree. Throws error if tree is improperly formatted.
"""
function SPR!(root::T)::T where T <:AbstractNode
    if length(post_order(root)) <= 2
        error("The tree is too small for SPR")
    end #if
    check_binary(root) || throw(ArgumentError("Not yet implemented for not binary trees"))
    perform_spr(root)
    return root
end #function

"""
    function risky_SPR(original_root::T)::T where T<:AbstractNode

Performs SPR on tree in place. Takes reference to root of tree
Returns copy of root of altered tree. Does not check for correct formatting of tree.
"""
function risky_SPR(original_root::T)::T where T<:AbstractNode
    root = deepcopy(original_root)
    return risky_SPR!(root)
end #function

"""
        risky_SPR!(root<:AbstractNode)

Performs SPR on tree in place.
Returns reference to root of altered tree. Does not check for correct formatting of tree.

* `root` : root node of tree.
"""
function risky_SPR!(root::T)::T where T<:AbstractNode
    return perform_spr(root)
end #func

"""
    perform_spr(root::T) where T <: GeneralNode

performs SPR on binary tree.
Returns root of tree post-SPR.

* `root` : Node of tree on which to perform SPR.
"""
function perform_spr(root::T) where T <: GeneralNode
    # find node to move
    available = [n.num for n in post_order(root)]
    n = rand(available)
    tn::T = find_num(root, n) #this is the root of the subtree which will be moved
    while tn.root || get_mother(tn).root
        n = rand(available)
        tn = find_num(root, n) #this is the root of the subtree which will be moved
    end # while
    tn_mother = get_mother(tn)
    tn_sister = get_sister(tn)
    tn_gm = get_mother(tn_mother)
    remove_child!(tn_gm, tn_mother)
    remove_child!(tn_mother, tn_sister)
    add_child!(tn_gm, tn_sister)
    # find target
    available = [n.num for n in post_order(root)]
    n = rand(available)
    target::T = find_num(root, n) #this is the target of the movement
    while target.root
        n = rand(available)
        target = find_num(root, n)
    end # while
    target_mother = get_mother(target)
    remove_child!(target_mother, target)
    add_child!(target_mother, tn_mother)
    add_child!(tn_mother, target)
    set_binary!(root)
    return root
end #func

"""
    perform_spr(root::T,subtree::T) where T <: GeneralNode

performs SPR on binary tree.
Returns root of tree post-SPR.

* `root` : Node of tree on which to perform SPR.
* `subtree`   : Root of subtree to be moved.
"""
function perform_spr(root::T,subtree::T) where T <: GeneralNode
    tn_mother = get_mother(subtree)
    tn_sister = get_sister(subtree)
    tn_gm = get_mother(tn_mother)
    remove_child!(tn_gm, tn_mother)
    remove_child!(tn_mother, tn_sister)
    add_child!(tn_gm, tn_sister)
    # find target
    available = [n.num for n in post_order(root)]
    n = rand(available)
    target::T = find_num(root, n) #this is the target of the movement
    while target.root
        n = rand(available)
        target = find_num(root, n)
    end # while
    target_mother = get_mother(target)
    remove_child!(target_mother, target)
    add_child!(target_mother, tn_mother)
    add_child!(tn_mother, target)
    set_binary!(root)
    return root
end #func

"""
    perform_spr(root::T,subtree::T,target::T) where T <: GeneralNode

performs SPR on binary tree.
Returns root of tree post-SPR.

* `root` : Node of tree on which to perform SPR.
* `subtree`   : Root of subtree to be moved.
* `target`    : Target of SPR.
"""
function perform_spr(root::T,subtree::T,target::T) where T <: GeneralNode
    @assert !subtree.root ["subtree cannot be the root node!"]
    @assert !get_mother(subtree).root ["subtree cannot be a child of the root node!"]
    tn_mother = get_mother(subtree)
    tn_sister = get_sister(subtree)
    tn_gm = get_mother(tn_mother)
    remove_child!(tn_gm, tn_mother)
    remove_child!(tn_mother, tn_sister)
    add_child!(tn_gm, tn_sister)
    @assert !target.root ["target cannot be the root!"]
    target_mother = get_mother(target)
    remove_child!(target_mother, target)
    add_child!(target_mother, tn_mother)
    add_child!(tn_mother, target)
    set_binary!(root)
    return root
end #func


text:
let
    # Should this be Array{Tainted,1} ?
    @test testval((Y,X) -> typeof(X), [1], [2]) == Array{Int, 1}
    @test testval((Y,X) -> typeof(Y), [1], [2]) == Array{Int, 1}

    @test testval((Y,X) -> eltype(X), [1], [2]) == Int
    @test testval((Y,X) -> typeof(X[1]), [1], [2]) == Int
    @test testval((Y,X) -> typeof(X[1]/1), [1], [2]) == Float64

    f(y,x) = x .+ 1

    @test testval((Y,X) -> broadcasted(+, X, 1)[1], [1], [2]) == 3
    @test testval(f, [1], [2]) == [3]
    @test sparse(testmeta(f, [1], [2])[1]) == sparse([], [], true, 1, 1)

    g(y,x) = y[:] .= x .+ 1
    #g(y,x) = y .= x .+ 1 -- memove

    @test sparse(testmeta(g, [1], [2])[1]) == sparse([1], [1], true)
end


text:
  cd(@__DIR__)
  using Pkg
  pkg"activate ."

  function main()
    include(joinpath("src", "MyFirstGenieApp.jl"))
  end; main()


text:
module day20

include("inputs.jl")

struct Tile1
    id::Int
    tile::Array{Int,2}
    edges::Dict{String,Int}
end

ex = inputs.ex_d20();
input = inputs.input_d20();

tiles_raw = split(ex, "\n\n")
tile_raw = tiles_raw[1]

function parse_tile(tile_raw)
    lines = split(tile_raw, "\n")
    id = parse(Int, lines[1][6:9])

    tile_s = replace.(lines[2:end], "."=>"0") .|>x->replace(x, "#"=>"1")
    tile = permutedims(vcat(hcat([[parse(Int, c) for c in line] for line in tile_s]...)))

    edges = Dict{}()
    edges["t"] = parse(Int, tile_s[1,:][1], base=2)
    edges["b"] = parse(Int, tile_s[end,:][1], base=2)
    edges["l"] = parse(Int, String([c[1] for c in tile_s]), base=2)
    edges["r"] = parse(Int, String([c[end] for c in tile_s]), base=2)

    return Tile1(id, tile, edges)
end

function rotate(tile)
    edges = Dict{}()
    edges["t"] = tile.edges["l"]
    edges["b"] = tile.edges["r"]
    edges["l"] = tile.edges["b"]
    edges["r"] = tile.edges["t"]
    return Tile1(tile.id, rotr90(tile.tile), edges)
end

function flip(tile)
    edges = Dict{}()
    edges["t"] = parse(Int, reverse(bitstring(tile.edges["t"])[end-9:end]), base=2)
    edges["b"] = parse(Int, reverse(bitstring(tile.edges["b"])[end-9:end]), base=2)
    edges["l"] = tile.edges["r"]
    edges["r"] = tile.edges["l"]
    return Tile1(tile.id, reverse(tile.tile, dims=2), edges)
end

apply(f, tile) = f(tile)

# parse(Int, reverse(bitstring(6)[end-11:end]), base=2)
end


text:
verbose = true
 using Dates
 if verbose println(now()) end
 if verbose && !isdefined(Main,:RvSpectML)  println("# Loading RvSpecML")    end
 using RvSpectMLBase
 using EchelleInstruments#, EchelleInstruments.NEID
 using RvSpectML
 if verbose println("# Loading NeidSolarScripts")    end
 using SunAsAStar
 using NeidSolarScripts
 println("# Loading other packages 1/2")
 using ArgParse
 using Markdown

println("# Parsing arguments...")
   #lsf_width = 3.0e3
   function parse_commandline()
     s = ArgParseSettings( description = "Calculate RVs from JLD2 file with NEID's daily order CCFs.")
     #import_settings!(s, s_files_only, args_only=false)
     add_arg_group!(s, "Files", :argg_files)
     @add_arg_table! s begin
         "input"
             help = "Filename for input with daily order CCFs (jld2)"
             arg_type = String
             default = "daily_ccfs.jld2"
         "output"
             help = "Filename for output dailiy RVs (csv)"
             arg_type = String
             default = "daily_rvs.csv"
         "--overwrite"
            help = "Specify it's ok to overwrite the output file."
            #default = true
            action = :store_true
      end
      #=
      add_arg_group!(s, "CCF parameters", :argg_ccf_param)
      @add_arg_table! s begin
        "--orders_to_use"
           help = "First and last order _index_ to compute CCFs for"
           nargs = 2
           arg_type = Int64
           #default = [ min_order(NEID2D()), max_order(NEID2D()) ]
           #default = [ first(orders_to_use_default(NEID2D())), last(orders_to_use_default(NEID2D())) ]
        "--mask_scale_factor"
            help = "Specify CCF mask width scale as multiple of NEID default v width " * string(default_ccf_mask_v_width(NEID2D()))
            arg_type = Float64
            #default = round(lsf_width/default_ccf_mask_v_width(NEID2D()), sigdigits=3)
        "--line_width_50_default"
            help = "Specify default line full width half maximum"
            arg_type = Float64
            #default = 7.9e3 # m/2
      end
      =#
      add_arg_group!(s, "Filter manifest for ", :argg_filter_param)
      @add_arg_table! s begin
         "--target"
            help = "Target field"
            arg_type = String
            default = "Sun"
         "--datestr"
            help = "Filenames that contain date string."
            arg_type = String
        "--max_solar_hour_angle"
           help = "Maximum absolute value of solar hour angle."
           arg_type = Float64
           default = 3.12
         "--max_airmass"
            help = "Maximum airmass."
            arg_type = Float64
            default = 2.0
         "--min_expmeter"
            help = "Minimum exposure meter flux relative to model flux."
            arg_type = Float64
            default = 6e4
#=
         "--min_expmeter_factor"
            help = "Minimum exposure meter flux relative to model flux."
            arg_type = Float64
            default = 0.9
         "--max_expmeter_factor"
            help = "Maximum exposure meter flux relative to model flux."
            arg_type = Float64
            default = 1.1
=#
         "--max_expmeter_frac_rms"
            help = "Maximum fractional RMS of exposure meter flux."
            arg_type = Float64
            default = 0.01
#=
         "--min_pyrhelio_factor"
            help = "Minimum pyrheliometer flux relative to model flux."
            arg_type = Float64
            default = 0.9
         "--max_pyrhelio_factor"
            help = "Maximum pyrheliometer flux relative to model flux."
            arg_type = Float64
            default = 1.1
=#
         "--max_pyrhelio_frac_rms"
            help = "Maximum fractional RMS of pyrheliometer meter flux."
            arg_type = Float64
            default = 0.0035
         #=
         "--min_drp_snr"
            help = "Minimum extracted SNR reported by NEID DRP."
            arg_type = Float64
            #default = 0.0
         =#
            #=
         "--start_time"
            help = "Specify daily start time for CCFs to be processed (HH MM)"
            nargs = 2
            arg_type = Int64
            default = [0, 0]
            #default = [18, 30]
        "--stop_time"
           help = "Specify daily stop time for CCFs to be processed (HH MM)"
           nargs = 2
           arg_type = Int64
           #default = [ min_order(NEID2D()), max_order(NEID2D()) ]
           #default = [22, 30]
           default = [23, 59]
           =#
         "--max_num_spectra"
            help = "Maximum number of spectra to process."
            arg_type = Int
            default = 300   # TODO: Increase  after finish testing
      end

     return parse_args(s)
 end
 args = parse_commandline()

if verbose   println("# Loading other packages 2/2")    end
 using CSV, DataFrames, Query, Dates
 using JLD2, FileIO, MD5 #, SHA
 using StatsBase, Statistics, NaNMath

 # Filename arguments
 @assert isfile(args["input"]) || islink(args["input"])
 @assert match(r"\.jld2$",args["input"]) != nothing
 daily_ccf_filename = args["input"]

 args["overwrite"] = true   # TODO: Comment after done testing
 if isfile(args["output"]) && !args["overwrite"]
    error("# Output file " * args["output"] * " already exists (size = " * string(filesize(args["output"])) * " ).")
 end
 @assert !isfile(args["output"]) || args["overwrite"]
 @assert match(r"\.csv$",args["output"]) != nothing
 daily_rvs_filename = args["output"]
 touch(daily_rvs_filename)

#=
 if isfile(args["report"]) && !args["overwrite"]
    error("# Report file " * args["report"] * " already exists (size = " * string(filesize(args["report"])) * " ).")
 end
 @assert !isfile(args["report"]) || args["overwrite"] == true
 @assert match(r"\.md$",args["report"]) != nothing
 daily_report_filename = args["report"]
 touch(daily_report_filename)
=#

 file_hashes = Dict{String,String}()
 start_processing_time = now()

if filesize(daily_ccf_filename) > 0
   println("# Reading daily CCFs from ", daily_ccf_filename)
   input_data = load(daily_ccf_filename)
else
   println("# Empy daily CCF file.  Creating empty ", daily_rvs_filename)
   touch(daily_rvs_filename)
   exit(0)
end

println("# Filtering for usable observations...")
manifest = input_data["manifest"]

manifest_use = manifest |>
 @filter( args["target"] == nothing || _.target == args["target"] ) |>
 @filter( args["datestr"] == nothing || occursin(args["datestr"],_.target) ) |>
 @filter( _.driftfun == "dailymodel0" ) |>
 @filter( args["max_airmass"] == nothing || _.airmass <= args["max_airmass"] ) |>
 @filter( args["max_solar_hour_angle"] == nothing || abs(_.hour_angle) <= args["max_solar_hour_angle"] ) |>
 #@filter( args["start_time"] == nothing || Time(julian2datetime(_.bjd)) >= start_time ) |>
 #@filter( args["stop_time"] == nothing || Time(julian2datetime(_.bjd)) <= stop_time ) |> # TODO for other instruments may need to deal wtih cross end of 24 UTC
 @filter( args["min_expmeter"] == nothing || _.expmeter_mean >= args["min_expmeter"] ) |>
 DataFrame

if hasproperty(manifest_use,:mean_pyroflux) && hasproperty(manifest_use,:rms_pyroflux)
   manifest_use = manifest_use |>
      @filter( args["max_pyrhelio_frac_rms"] == nothing || _.rms_pyroflux <= args["max_pyrhelio_frac_rms"] * _.mean_pyroflux ) |>
      DataFrame
else
   println("# Warning:  mean_pyroflux and rms_pyroflux not found.  Reverting to expmeter_mean and expmeter_rms.")
   manifest_use = manifest_use |>
      @filter( args["max_expmeter_frac_rms"] == nothing || _.expmeter_rms <= args["max_expmeter_frac_rms"] * _.expmeter_mean ) |>
      DataFrame
end

manifest_use = manifest_use |>
   @take(args["max_num_spectra"] ) |> @orderby(_.bjd) |>
   DataFrame

#@assert size(manifest_use,1) >= 1
if !(size(manifest_use,1) >= 1)
   println("# No usable observations.  Creating empty ", daily_rvs_filename)
   touch(daily_rvs_filename)
   exit(0)
end

println("# Found ", size(manifest_use,1), " files of ",  size(manifest,1), " to use for RVs.")
#   df_out = select(manifest_use,[:drp_ccfjdmod=>:jd_drp,:drp_ccfrvmod => :rv_drp,:drp_dvrmsmod => :σrv_drp ], :Δv_diff_ext, :Δfwhm², :solar_hour_angle, :airmass, :stel,:staz,:alt_sun,:az_sun, :sol_dist, :expmeter_mean, :expmeter_rms, [:mean_pyroflux=>:pyrflux_mean, :rms_pyroflux=>:pyrflux_rms], :exptime, :mean_Δt, :Filename)
   df_out = manifest_use |> @rename(:drp_ccfjdmod=>:jd_drp,:drp_ccfrvmod => :rv_drp,:drp_dvrmsmod => :σrv_drp,
                                    :mean_pyroflux=>:pyrflux_mean, :rms_pyroflux=>:pyrflux_rms ) |> DataFrame

println("# Writing daily RVs.")
 CSV.write(daily_rvs_filename,df_out)


text:
"""
    triangle_area(v1, v2, v3)

Area of a triangle defined by the vertices `v1`, `v2` and `v3`.
"""
function triangle_area(v1, v2, v3)
    u1 = v2 - v1
    u2 = v3 - v1

    return abs(u1[1]*u2[2] - u1[2]*u2[1])/2

end

"""
    isinsideconvex(p, vertices; tol=sqrt(eps()))

Determines if a point `p` is inside a convex region defined by the vertices in the array `vertices`.
"""
function isinsideconvex(p, vertices; tol=sqrt(eps()))
    n = length(vertices)

    center = [0.0, 0.0]
    for i=1:n
        center += vertices[i]
    end
    center = center/n

    convex_area = 0.0
    for i=1:n
        convex_area += triangle_area(center, vertices[i], vertices[mod1(i+1,n)])
    end

    test_area = 0.0
    for i=1:n
        test_area += triangle_area(p, vertices[i], vertices[mod1(i+1,n)])
    end

    return isapprox(test_area, convex_area, rtol = tol)
end


text:
### A Pluto.jl notebook ###
# v0.12.4

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 2360bd70-1a15-11eb-3b28-75fdf85c0706
using NeXLDatabase

# ╔═╡ 290a7a40-1a15-11eb-338e-1d2c5e513ab3
using Gadfly, DataFrames, PlutoUI

# ╔═╡ 30233060-1a15-11eb-2d1f-8b7537f66d45
db=openNeXLDatabase("C:\\Users\\nritchie\\Documents\\DrWatson\\K-Ratio Project\\data\\exp_pro\\kratio.db")

# ╔═╡ 7cfccdb0-1a15-11eb-37b7-4f1f1aff8f2f
begin
	labs = findall(db, DBLaboratory)
	md"""
	Select a lab: $(@bind lab Select([ repr(i) => repr(lab) for (i, lab) in enumerate(labs) ]))
	"""
end

# ╔═╡ ee29aa10-1a17-11eb-3fac-5166acfccd86
md"""
#### Lab Members
"""

# ╔═╡ a391aa40-1a15-11eb-3bd4-45eeff75d6e3
begin
	thelab = labs[parse(Int,lab)]
	members = findall(db, DBPerson, thelab)
	function NeXLUncertainties.asa(::Type{DataFrame}, people::AbstractArray{DBPerson})
		return DataFrame(
			PKey = [ person.pkey for person in people ],
			Name = [ person.name for person in people ],
			EMail = [ person.email for person in people ]
		)
	end
	asa(DataFrame, members)
end

# ╔═╡ fc72e40e-1a17-11eb-2f3c-038a1f5b98d7
md"""
#### Lab Instruments
"""

# ╔═╡ c5eb0000-1a15-11eb-0f46-99d1ca2092f6
begin
	instruments = findall(db, DBInstrument, thelab)
	md"""
	Select a instrument: $(@bind instidx Select([ repr(i) => repr(inst) for (i, inst) in enumerate(instruments) ]))
	"""
end

# ╔═╡ 063a1862-1a18-11eb-0424-ebdcc6ccc526
md"""
#### Detectors
"""

# ╔═╡ eb7faa9e-1a15-11eb-23af-0ba535efb5f6
begin
	theinstrument = instruments[parse(Int,instidx)]
	detectors = findall(db, DBDetector, theinstrument)
	md"""
	Select a detector: $(@bind detidx Select([ repr(i) => repr(det) for (i, det) in enumerate(detectors) ]))
	"""
end

# ╔═╡ 150e5450-1a18-11eb-1f6d-8fe22b72fa94
begin
	thedetector = detectors[parse(Int,detidx)]
	asa(DataFrame,thedetector)
end

# ╔═╡ Cell order:
# ╠═2360bd70-1a15-11eb-3b28-75fdf85c0706
# ╠═290a7a40-1a15-11eb-338e-1d2c5e513ab3
# ╟─30233060-1a15-11eb-2d1f-8b7537f66d45
# ╟─7cfccdb0-1a15-11eb-37b7-4f1f1aff8f2f
# ╟─ee29aa10-1a17-11eb-3fac-5166acfccd86
# ╟─a391aa40-1a15-11eb-3bd4-45eeff75d6e3
# ╟─fc72e40e-1a17-11eb-2f3c-038a1f5b98d7
# ╟─c5eb0000-1a15-11eb-0f46-99d1ca2092f6
# ╟─063a1862-1a18-11eb-0424-ebdcc6ccc526
# ╟─eb7faa9e-1a15-11eb-23af-0ba535efb5f6
# ╟─150e5450-1a18-11eb-1f6d-8fe22b72fa94


text:
__precompile__()

module Bifurcations

    using Luxor, Colors, FixedPointNumbers

    include("utils.jl")
    include("populations.jl")
    include("gradients.jl")

    function plot(pop::Population,
            colmap::Gradient;
            iters::Int64=600,
            steps::Int64=1000,
            xmin::Float64=0.0, xmax::Float64=1.0,
            ymin::Float64=0.0, ymax::Float64=4.0)

        plot(pop, colmap, colmap;
            iters=iters,
            steps=steps,
            xmin=xmin, xmax=xmax,
            ymin=ymin, ymax=ymax)
    end

    function plot(pop::Population,
            colmap_r::Gradient,
            colmap_i::Gradient;
            iters::Int64=600,
            steps::Int64=1000,
            xmin::Float64=0.0, xmax::Float64=1.0,
            ymin::Float64=0.0, ymax::Float64=4.0)

        # We precompute colors for each iteration for performance
        colors_r::Array = [grad(colmap_r, n, 1, iters) for n = 1:iters]
        colors_i::Array = [grad(colmap_i, n, 1, iters) for n = 1:iters]

        coords_r = [Array{Point, 1}() for _ = 1:iters]
        coords_i = [Array{Point, 1}() for _ = 1:iters]

        for x = 1:steps
            p::ComplexF64 = pop.initial
            x::Float64 = interp(x, 1, steps, ymin, ymax)

            for n = 1:iters
                if abs2(p) > maximum(abs.((xmin, xmax, ymin, ymax)))^2
                    break
                end

                # n > card(f) to avoid initial horizontal lines on the plot
                if n > length(pop.fertilities)
                    y_r = 1 - real(p)
                    y_i = 1 - imag(p)
                    if ymin < y_r < ymax
                        push!(coords_r[n], Point(960x, 960y_r))
                    end
                    if ymin < y_i < ymax
                        push!(coords_i[n], Point(960x, 960y_i))
                    end
                end
                f = pop.fertilities[n % end + 1](x)
                p = pop.law(f, p)
            end
        end

        @svg begin
            setline(1)
            setlinecap("round")
            setlinejoin("round")
            for (coords, colors) = zip((coords_r, coords_i), (colors_r, colors_i))
                for (line_coords, color) = Iterators.reverse(zip(coords, colors))
                    setcolor(color)
                    try
                        # poly(line_coords, :stroke)
                        drawbezierpath(makebezierpath(simplify(line_coords, .5)), :stroke)
                        # drawbezierpath(makebezierpath(line_coords), :stroke)
                    catch e
                        if !isa(e, BoundsError)
                            throw(e)
                        end
                    end
                end
            end
        end
    end

    function __init__()
        pop = Population(logistic, [x -> x-.5im*x], .5)
        pop = Population(mandelbrot, [x -> x*im - 0.2], .5)
        pop = Population(logistic, [x->.5x+im], .5) # pour L.I

        black   = RGBA(0.0, 0.0, 0.0, 1.0)
        black_a = RGBA(1.0, 1.0, 1.0, 0.0)
        red     = RGBA(1.0, 0.0, 0.0, 1.0)
        red_a   = RGBA(1.0, 0.0, 0.0, 0.0)
        blue_a  = RGBA(0.0, 0.0, 1.0, 0.0)
        white   = RGBA(1.0, 1.0, 1.0, 1.0)

        black_to_red       = Gradient(black, red)
        black_to_red_fade  = Gradient(black, red_a)
        black_to_blue_fade = Gradient(black, blue_a)
        red_to_black       = Gradient(red,   black)
        red_to_white       = Gradient(red,   white)
        black              = Gradient(black, black)

        plot(pop, black_to_red_fade, black_to_blue_fade;
            # name = string(rand(0000:9999)),
            iters = 100,
            steps = 8000,
            xmin = -4.0, xmax = 4.0,
            ymin = -4.0, ymax = 4.0)
    end
end


text:
@testset "Parameter" begin
    π = Parameter(3.14, :π)
    x = DataFrame(; π=[3.14, 3.14, 3.14, 3.14, 3.14])

    @test π isa Parameter
    @test π.value == 3.14
    @test π.name == :π

    @test sample(π, 5) == x

    @test mean(π) == DataFrame(; π=3.14)

    to_standard_normal_space!(π, x)
    @test x == DataFrame(; π=[3.14, 3.14, 3.14, 3.14, 3.14])

    to_physical_space!(π, x)
    @test x == DataFrame(; π=[3.14, 3.14, 3.14, 3.14, 3.14])
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
#jl #! format: off
# # Hydropower Simulations with [PowerSimulations.jl](https://github.com/NREL-SIIP/PowerSimulations.jl)

# **Originally Contributed by**: Clayton Barrows and Sourabh Dalvi

# ## Introduction

# PowerSimulations.jl supports simulations that consist of sequential optimization problems
# where results from previous problems inform subsequent problems in a variety of ways.
# This example demonstrates a few of the options for modeling hydropower generation.

# ## Dependencies
using SIIPExamples

# ### Modeling Packages
using PowerSystems
using PowerSimulations
const PSI = PowerSimulations
using PowerSystemCaseBuilder

# ### Data management packages
using Dates
using DataFrames

# ### Optimization packages
using HiGHS # solver
solver = optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.05)
odir = mktempdir(".", cleanup = true) #tmpdir for build steps

# ### Data
# PowerSystemCaseBuilder links to some meaningless test data that is suitable for this example.
# We can load three systems of different resolution for the examples here:
c_sys5_hy_wk = build_system(SIIPExampleSystems, "5_bus_hydro_wk_sys")
c_sys5_hy_uc = build_system(SIIPExampleSystems, "5_bus_hydro_uc_sys")
c_sys5_hy_ed = build_system(SIIPExampleSystems, "5_bus_hydro_ed_sys")

c_sys5_hy_wk_targets = build_system(SIIPExampleSystems, "5_bus_hydro_wk_sys_with_targets")
c_sys5_hy_uc_targets = build_system(SIIPExampleSystems, "5_bus_hydro_uc_sys_with_targets")
c_sys5_hy_ed_targets = build_system(SIIPExampleSystems, "5_bus_hydro_ed_sys_with_targets")

# This line just overloads JuMP printing to remove double underscores added by PowerSimulations.jl
PSI.JuMP._wrap_in_math_mode(str) = "\$\$ $(replace(str, "__"=>"")) \$\$"

# ## Two PowerSimulations features determine hydropower representation.
# There are two principal ways that we can customize hydropower representation in
# PowerSimulations. First, we can play with the formulation applied to hydropower generators
# using the `DeviceModel`. We can also adjust how simulations are configured to represent
# different decision making processes and the information flow between those processes.

# ### Hydropower `DeviceModel`s

# First, the assignment of device formulations to particular device types gives us control
# over the representation of devices. This is accomplished by defining `DeviceModel`
# instances. For hydro power representations, we have two available generator types in
# PowerSystems:

print_tree(HydroGen)

# And in PowerSimulations, we have several available formulations that can be applied to
# the hydropower generation devices:

print_tree(PSI.AbstractHydroFormulation)

# Let's see what some of the different combinations create. First, let's apply the
# `HydroDispatchRunOfRiver` formulation to the `HydroEnergyReservoir` generators, and the
# `FixedOutput` formulation to `HydroDispatch` generators.
#  - The `FixedOutput` formulation just acts
# like a load subtractor, forcing the system to accept it's generation.
#  - The `HydroDispatchRunOfRiver` formulation represents the the energy flowing out of
# a reservoir. The model can choose to produce power with that energy or just let it spill by.
template = ProblemTemplate()
set_device_model!(template, HydroEnergyReservoir, HydroDispatchRunOfRiver)
set_device_model!(template, HydroDispatch, FixedOutput)
set_device_model!(template, PowerLoad, StaticPowerLoad)

prob = DecisionModel(template, c_sys5_hy_uc, horizon = 2)
build!(prob, output_dir = odir)

# Now we can see the resulting JuMP model:
PSI.get_jump_model(prob)

# The first two constraints are the power balance constraints that require the generation
# from the controllable `HydroEnergyReservoir` generators to be equal to the load (flat 10.0 for all time periods)
# minus the generation from the `HydroDispatch` generators [1.97, 1.983, ...]. The 3rd and 4th
# constraints limit the output of the `HydroEnergyReservoir` generator to the limit defined by the
# `max_activepwoer` time series. And the last 4 constraints are the lower and upper bounds of
# the `HydroEnergyReservoir` operating range.

#-

# Next, let's apply the `HydroDispatchReservoirBudget` formulation to the `HydroEnergyReservoir` generators.
template = ProblemTemplate()
set_device_model!(template, HydroEnergyReservoir, HydroDispatchReservoirBudget)
set_device_model!(template, PowerLoad, StaticPowerLoad)

prob = DecisionModel(template, c_sys5_hy_uc, horizon = 2)
build!(prob, output_dir = odir)

# And, the resulting JuMP model:
PSI.get_jump_model(prob)

# Finally, let's apply the `HydroDispatchReservoirStorage` formulation to the `HydroEnergyReservoir` generators.
template = ProblemTemplate()
set_device_model!(template, HydroEnergyReservoir, HydroDispatchReservoirStorage)
set_device_model!(template, PowerLoad, StaticPowerLoad)

prob = DecisionModel(template, c_sys5_hy_uc_targets, horizon = 24)
build!(prob, output_dir = odir)

#-
PSI.get_jump_model(prob)

# ### Multi-Stage `SimulationSequence`
# The purpose of a multi-stage simulation is to represent scheduling decisions consistently
# with the time scales that govern different elements of power systems.

# #### Multi-Day to Daily Simulation:
# In the multi-day model, we'll use a really simple representation of all system devices
# so that we can maintain computational tractability while getting an estimate of system
# requirements/capabilities.
template_md = ProblemTemplate()
set_device_model!(template_md, ThermalStandard, ThermalDispatchNoMin)
set_device_model!(template_md, PowerLoad, StaticPowerLoad)
set_device_model!(template_md, HydroEnergyReservoir, HydroDispatchReservoirStorage)

# For the daily model, we can increase the modeling detail since we'll be solving shorter
# problems.
template_da = ProblemTemplate()
set_device_model!(template_da, ThermalStandard, ThermalBasicUnitCommitment)
set_device_model!(template_da, PowerLoad, StaticPowerLoad)
set_device_model!(template_da, HydroEnergyReservoir, HydroDispatchReservoirStorage)
#-

problems = SimulationModels(
    decision_models = [
        DecisionModel(
            template_md,
            c_sys5_hy_wk_targets,
            name = "MD",
            optimizer = solver,
            system_to_file = false,
            initialize_model = false,
        ),
        DecisionModel(
            template_da,
            c_sys5_hy_uc_targets,
            name = "DA",
            optimizer = solver,
            system_to_file = false,
            initialize_model = false,
        ),
    ],
)

# This builds the sequence and passes the the energy dispatch schedule for the `HydroEnergyReservoir`
# generator from the "MD" problem to the "DA" problem in the form of an energy limit over the
# synchronized periods.

sequence = SimulationSequence(
    models = problems,
    feedforwards = Dict(
        "DA" => [
            EnergyLimitFeedforward(
                component_type = HydroEnergyReservoir,
                source = ActivePowerVariable,
                affected_values = [ActivePowerVariable],
                number_of_periods = get_forecast_horizon(c_sys5_hy_uc_targets),
            ),
        ],
    ),
    ini_cond_chronology = InterProblemChronology(),
);

#-

sim = Simulation(
    name = "hydro",
    steps = 1,
    models = problems,
    sequence = sequence,
    simulation_folder = odir,
)

build!(sim)

# We can look at the "MD" Model
PSI.get_jump_model(sim.models.decision_models[1])

# And we can look at the "DA" model
PSI.get_jump_model(sim.models.decision_models[2])

# And we can execute the simulation by running the following command
# ```julia
# execute!(sim)
# ```
#-

# #### 3-Stage Simulation:
transform_single_time_series!(c_sys5_hy_wk, 2, Hour(24)) # TODO fix PSI to enable longer intervals of stage 1

# For the real time model, we can increase the modeling detail since we'll be solving shorter
# problems.
template_ed = ProblemTemplate()
set_device_model!(template_ed, ThermalStandard, ThermalDispatchNoMin)
set_device_model!(template_ed, PowerLoad, StaticPowerLoad)
set_device_model!(template_ed, HydroEnergyReservoir, HydroDispatchReservoirStorage)

problems = SimulationModels(
    decision_models = [
        DecisionModel(
            template_md,
            c_sys5_hy_wk_targets,
            name = "MD",
            optimizer = solver,
            system_to_file = false,
            initialize_model = false,
        ),
        DecisionModel(
            template_da,
            c_sys5_hy_uc_targets,
            name = "DA",
            optimizer = solver,
            system_to_file = false,
            initialize_model = false,
        ),
        DecisionModel(
            template_ed,
            c_sys5_hy_ed_targets,
            name = "ED",
            optimizer = solver,
            system_to_file = false,
            initialize_model = false,
        ),
    ],
)

sequence = SimulationSequence(
    models = problems,
    feedforwards = Dict(
        "DA" => [
            EnergyLimitFeedforward(
                component_type = HydroEnergyReservoir,
                source = ActivePowerVariable,
                affected_values = [ActivePowerVariable],
                number_of_periods = get_forecast_horizon(c_sys5_hy_uc_targets),
            ),
        ],
        "ED" => [
            EnergyLimitFeedforward(
                component_type = HydroEnergyReservoir,
                source = ActivePowerVariable,
                affected_values = [ActivePowerVariable],
                number_of_periods = get_forecast_horizon(c_sys5_hy_ed_targets),
            ),
        ],
    ),
    ini_cond_chronology = InterProblemChronology(),
);

#-

sim = Simulation(
    name = "hydro",
    steps = 1,
    models = problems,
    sequence = sequence,
    simulation_folder = odir,
)

#-

build!(sim)

# We can look at the "MD" Model
PSI.get_jump_model(sim.models.decision_models[1])

# And we can look at the "DA" model
PSI.get_jump_model(sim.models.decision_models[2])

# And we can look at the "ED" model
PSI.get_jump_model(sim.models.decision_models[3])

# And we can execute the simulation by running the following command
# ```julia
# execute!(sim)
# ```
#-

max_stars_repo_name: NREL-SIIP/Examples

max_issues_repo_name: NREL-SIIP/Examples

max_forks_repo_name: NREL-SIIP/Examples