text:
using
    Wavelets,
    WaveletsExt,
    LinearAlgebra,
    Random,
    DataFrames,
    Statistics,
    CSV

function addnoise(x::AbstractArray{<:Number,1}, s::Real=0.1)
    ϵ = randn(length(x))
    ϵ = ϵ/norm(ϵ)
    y = x + ϵ * s * norm(x)
    return y
end

# compute metrics
function computemetrics(x::Array{T}, x₀::Array{T}) where T<:Number
    # PSNR
    r_psnr = [psnr(x[:,i], x₀[:,i]) for i in axes(x₀,2)]
    # SSIM
    r_ssim = [ssim(x[:,i], x₀[:,i]) for i in axes(x₀,2)]
    return mean([r_psnr r_ssim], dims = 1)
end

# single comparison
function singlecomparison(x::Vector{T}, wt::DiscreteWavelet,
        noise::AbstractFloat, shifts::Integer, samples::Integer) where T<:Number

    n = size(x,1)
    result = Array{AbstractFloat, 2}(undef, (109,3))
    # generate signals
    X₀ = duplicatesignals(x, samples, shifts)
    X = hcat([addnoise(X₀[:,i], noise) for i in 1:samples]...)

    # NOISY
    result[1,:] = [0 computemetrics(X, X₀)]

    # DWT
    L = maxtransformlevels(n)
    y = dwtall(X, wt, L)
    σ₁ = [noisest(y[:,i], false) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,i], false) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :dwt, wt, L=L, dnt=VisuShrink(n), estnoise=σ₁, bestTH=nothing)
    result[2,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :dwt, wt, L=L, dnt=VisuShrink(n), estnoise=σ₁, bestTH=mean)
    result[3,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed denoiseall(y, :dwt, wt, L=L, dnt=VisuShrink(n), estnoise=σ₁, bestTH=median)
    result[4,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :dwt, wt, L=L, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=nothing)
    result[5,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :dwt, wt, L=L, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=mean)
    result[6,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed denoiseall(y, :dwt, wt, L=L, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=median)
    result[7,:] = [time computemetrics(X̂, X₀)]

    # WPT-L4
    tree = maketree(n, 4, :full)
    y = wptall(X, wt, 4)
    σ₁ = [noisest(y[:,i], false, tree) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,i], false, tree) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=nothing)
    result[8,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=mean)
    result[9,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=median)
    result[10,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=nothing)
    result[11,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=mean)
    result[12,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=median)
    result[13,:] = [time computemetrics(X̂, X₀)]

    # WPT-L8
    tree = maketree(n, 8, :full)
    y = wptall(X, wt, 8)
    σ₁ = [noisest(y[:,i], false, tree) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,i], false, tree) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=nothing)
    result[14,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=mean)
    result[15,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=median)
    result[16,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=nothing)
    result[17,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=mean)
    result[18,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=median)
    result[19,:] = [time computemetrics(X̂, X₀)]

    # WPT-BT
    Xw = wpdall(X, wt)
    tree = bestbasistreeall(Xw, BB())
    y = getbasiscoefall(Xw, tree)
    σ₁ = [noisest(y[:,i], false, tree[:,i]) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,i], false, tree[:,i]) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed hcat([
        denoise(y[:,i], :wpt, wt, tree=tree[:,i], dnt=VisuShrink(n), estnoise=σ₁[i]) for i in 1:samples
    ]...)
    result[20,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed hcat([
        denoise(y[:,i], :wpt, wt, tree=tree[:,i], dnt=VisuShrink(n), estnoise=mean(σ₁)) for i in 1:samples
    ]...)
    result[21,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed hcat([
        denoise(y[:,i], :wpt, wt, tree=tree[:,i], dnt=VisuShrink(n), estnoise=median(σ₁)) for i in 1:samples
    ]...)
    result[22,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed hcat([
        denoise(y[:,i], :wpt, wt, tree=tree[:,i], dnt=RelErrorShrink(), estnoise=σ₂[i]) for i in 1:samples
    ]...)
    result[23,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed hcat([
        denoise(y[:,i], :wpt, wt, tree=tree[:,i], dnt=RelErrorShrink(), estnoise=mean(σ₂)) for i in 1:samples
    ]...)
    result[24,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed hcat([
        denoise(y[:,i], :wpt, wt, tree=tree[:,i], dnt=RelErrorShrink(), estnoise=median(σ₂)) for i in 1:samples
    ]...)
    result[25,:] = [time computemetrics(X̂, X₀)]

    # JBB
    tree = bestbasistree(Xw, JBB())
    y = getbasiscoefall(Xw, tree)
    σ₁ = [noisest(y[:,i], false, tree) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,i], false, tree) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=nothing)
    result[26,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=mean)
    result[27,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=median)
    result[28,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=nothing)
    result[29,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=mean)
    result[30,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=median)
    result[31,:] = [time computemetrics(X̂, X₀)]

    # LSDB
    tree = bestbasistree(Xw, LSDB())
    y = getbasiscoefall(Xw, tree)
    σ₁ = [noisest(y[:,i], false, tree) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,i], false, tree) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=nothing)
    result[32,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=mean)
    result[33,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=median)
    result[34,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=nothing)
    result[35,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=mean)
    result[36,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed denoiseall(y, :wpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=median)
    result[37,:] = [time computemetrics(X̂, X₀)]

    # ACWT
    L = maxtransformlevels(n)
    y = acdwtall(X, wt, L)
    σ₁ = [noisest(y[:,:,i], true) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,:,i], true, nothing) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :acwt, wt, L=L, dnt=VisuShrink(n), estnoise=σ₁, bestTH=nothing)
    result[38,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :acwt, wt, L=L, dnt=VisuShrink(n), estnoise=σ₁, bestTH=mean)
    result[39,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed denoiseall(y, :acwt, wt, L=L, dnt=VisuShrink(n), estnoise=σ₁, bestTH=median)
    result[40,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :acwt, wt, L=L, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=nothing)
    result[41,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :acwt, wt, L=L, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=mean)
    result[42,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed denoiseall(y, :acwt, wt, L=L, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=median)
    result[43,:] = [time computemetrics(X̂, X₀)]

    # ACWPT-L4
    tree = maketree(n, 4, :full)
    y = acwpdall(X, wt, 8)
    σ₁ = [noisest(y[:,:,i], true, tree) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,:,i], true, tree) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=nothing)
    result[44,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=mean)
    result[45,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=median)
    result[46,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=nothing)
    result[47,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=mean)
    result[48,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=median)
    result[49,:] = [time computemetrics(X̂, X₀)]

    # ACWPT-L8
    tree = maketree(n, 8, :full)
    σ₁ = [noisest(y[:,:,i], true, tree) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,:,i], true, tree) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=nothing)
    result[50,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=mean)
    result[51,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=median)
    result[52,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=nothing)
    result[53,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=mean)
    result[54,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=median)
    result[55,:] = [time computemetrics(X̂, X₀)]

    # ACWPT-BT
    tree = bestbasistreeall(y, BB(redundant=true))
    σ₁ = [noisest(y[:,:,i], true, tree[:,i]) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,:,i], true, tree[:,i]) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed hcat([
        denoise(y[:,:,i], :acwpt, wt, tree=tree[:,i], dnt=VisuShrink(n), estnoise=σ₁[i]) for i in 1:samples
    ]...)
    result[56,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed hcat([
        denoise(y[:,:,i], :acwpt, wt, tree=tree[:,i], dnt=VisuShrink(n), estnoise=mean(σ₁)) for i in 1:samples
    ]...)
    result[57,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed hcat([
        denoise(y[:,:,i], :acwpt, wt, tree=tree[:,i], dnt=VisuShrink(n), estnoise=median(σ₁)) for i in 1:samples
    ]...)
    result[58,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed hcat([
        denoise(y[:,:,i], :acwpt, wt, tree=tree[:,i], dnt=RelErrorShrink(), estnoise=σ₂[i]) for i in 1:samples
    ]...)
    result[59,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed hcat([
        denoise(y[:,:,i], :acwpt, wt, tree=tree[:,i], dnt=RelErrorShrink(), estnoise=mean(σ₂)) for i in 1:samples
    ]...)
    result[60,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed hcat([
        denoise(y[:,:,i], :acwpt, wt, tree=tree[:,i], dnt=RelErrorShrink(), estnoise=median(σ₂)) for i in 1:samples
    ]...)
    result[61,:] = [time computemetrics(X̂, X₀)]

    # ACWPT-JBB
    tree = bestbasistree(y, JBB(redundant=true))
    σ₁ = [noisest(y[:,:,i], true, tree) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,:,i], true, tree) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=nothing)
    result[62,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=mean)
    result[63,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=median)
    result[64,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=nothing)
    result[65,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=mean)
    result[66,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=median)
    result[67,:] = [time computemetrics(X̂, X₀)]

    # ACWPT-LSDB
    tree = bestbasistree(y, LSDB(redundant=true))
    σ₁ = [noisest(y[:,:,i], true, tree) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,:,i], true, tree) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=nothing)
    result[68,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=mean)
    result[69,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=median)
    result[70,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=nothing)
    result[71,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=mean)
    result[72,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed denoiseall(y, :acwpt, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=median)
    result[73,:] = [time computemetrics(X̂, X₀)]

    # SDWT
    L = maxtransformlevels(n)
    y = sdwtall(X, wt, L)
    σ₁ = [noisest(y[:,:,i], true) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,:,i], true, nothing) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :sdwt, wt, L=L, dnt=VisuShrink(n), estnoise=σ₁, bestTH=nothing)
    result[74,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :sdwt, wt, L=L, dnt=VisuShrink(n), estnoise=σ₁, bestTH=mean)
    result[75,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed denoiseall(y, :sdwt, wt, L=L, dnt=VisuShrink(n), estnoise=σ₁, bestTH=median)
    result[76,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :sdwt, wt, L=L, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=nothing)
    result[77,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :sdwt, wt, L=L, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=mean)
    result[78,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed denoiseall(y, :sdwt, wt, L=L, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=median)
    result[79,:] = [time computemetrics(X̂, X₀)]

    # SWPT-L4
    tree = maketree(n, 4, :full)
    y = swpdall(X, wt, 8)
    σ₁ = [noisest(y[:,:,i], true, tree) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,:,i], true, tree) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=nothing)
    result[80,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=mean)
    result[81,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=median)
    result[82,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=nothing)
    result[83,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=mean)
    result[84,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=median)
    result[85,:] = [time computemetrics(X̂, X₀)]

    # SWPT-L8
    tree = maketree(n, 8, :full)
    σ₁ = [noisest(y[:,:,i], true, tree) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,:,i], true, tree) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=nothing)
    result[86,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=mean)
    result[87,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=median)
    result[88,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=nothing)
    result[89,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=mean)
    result[90,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=median)
    result[91,:] = [time computemetrics(X̂, X₀)]

    # SWPT-BT
    tree = bestbasistreeall(y, BB(redundant=true))
    σ₁ = [noisest(y[:,:,i], true, tree[:,i]) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,:,i], true, tree[:,i]) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed hcat([
        denoise(y[:,:,i], :swpd, wt, tree=tree[:,i], dnt=VisuShrink(n), estnoise=σ₁[i]) for i in 1:samples
    ]...)
    result[92,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed hcat([
        denoise(y[:,:,i], :swpd, wt, tree=tree[:,i], dnt=VisuShrink(n), estnoise=mean(σ₁)) for i in 1:samples
    ]...)
    result[93,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed hcat([
        denoise(y[:,:,i], :swpd, wt, tree=tree[:,i], dnt=VisuShrink(n), estnoise=median(σ₁)) for i in 1:samples
    ]...)
    result[94,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed hcat([
        denoise(y[:,:,i], :swpd, wt, tree=tree[:,i], dnt=RelErrorShrink(), estnoise=σ₂[i]) for i in 1:samples
    ]...)
    result[95,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed hcat([
        denoise(y[:,:,i], :swpd, wt, tree=tree[:,i], dnt=RelErrorShrink(), estnoise=mean(σ₂)) for i in 1:samples
    ]...)
    result[96,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed hcat([
        denoise(y[:,:,i], :swpd, wt, tree=tree[:,i], dnt=RelErrorShrink(), estnoise=median(σ₂)) for i in 1:samples
    ]...)
    result[97,:] = [time computemetrics(X̂, X₀)]

    # SJBB
    tree = bestbasistree(y, JBB(redundant=true))
    σ₁ = [noisest(y[:,:,i], true, tree) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,:,i], true, tree) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=nothing)
    result[98,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=mean)
    result[99,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=median)
    result[100,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=nothing)
    result[101,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=mean)
    result[102,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=median)
    result[103,:] = [time computemetrics(X̂, X₀)]

    # SLSDB
    tree = bestbasistree(y, LSDB(redundant=true))
    σ₁ = [noisest(y[:,:,i], true, tree) for i in 1:samples]
    σ₂ = [relerrorthreshold(y[:,:,i], true, tree) for i in 1:samples]
    ## visushrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=nothing)
    result[104,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=mean)
    result[105,:] = [time computemetrics(X̂, X₀)]
    ## visushrink & bestTH = median
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=VisuShrink(n), estnoise=σ₁, bestTH=median)
    result[106,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = individual
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=nothing)
    result[107,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = mean
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=mean)
    result[108,:] = [time computemetrics(X̂, X₀)]
    ## relerrorshrink & bestTH = median
    X̂, time = @timed denoiseall(y, :swpd, wt, tree=tree, dnt=RelErrorShrink(), estnoise=σ₂, bestTH=median)
    result[109,:] = [time computemetrics(X̂, X₀)]

    return result
end

# repeated comparisons
function repeatedcomparisons(x::Vector{T}, wt::DiscreteWavelet,
        noise::AbstractFloat, shifts::Integer, samples::Integer,
        repeats::Integer) where T<:Number

    results =
        [singlecomparison(x, wt, noise, shifts, samples) for _ in 1:repeats] |>
        y -> cat(y..., dims = 3) |>
        y -> mean(y, dims = 3) |>
        y -> reshape(y, :, 3) |>
        y -> DataFrame(y, [:Time, :PSNR, :SSIM])

    params = DataFrame(
        transform = ["None";
            repeat(["DWT", "WPT-L4", "WPT-L8", "WPT-BT", "JBB", "LSDB"], inner=6);
            repeat(["ACWT", "ACWPT-L4", "ACWPT-L8", "ACWPT-BT", "ACJBB", "ACLSDB"], inner=6);
            repeat(["SDWT", "SWPT-L4", "SWPT-L8", "SWPT-BT", "SJBB", "SLSDB"], inner=6)],
        method = ["None";
            repeat(["VisuShrink", "RelErrorShrink"], inner=3, outer=18)],
        THType = ["None"; repeat(["Hard"], inner=108)],
        bestTH = ["None";
            repeat(["Individual", "Average", "Median"], outer=36)]
    )
    return [params results]
end

wt = wavelet(WT.db4)
L = 8

# compute and save results
samples = 100
repeats = 50
results = Dict{String, DataFrame}()
results["blocks"] = repeatedcomparisons(generatesignals(:blocks,L), wt, 0.5, 2, samples, repeats)
CSV.write("./results/blocks.csv", results["blocks"])
results["bumps"] = repeatedcomparisons(generatesignals(:bumps,L), wt, 0.5, 2, samples, repeats)
CSV.write("./results/bumps.csv", results["bumps"])
results["heavisine"] = repeatedcomparisons(generatesignals(:heavisine,L), wt, 0.5, 2, samples, repeats)
CSV.write("./results/heavisine.csv", results["heavisine"])
results["doppler"] = repeatedcomparisons(generatesignals(:doppler,L), wt, 0.5, 2, samples, repeats)
CSV.write("./results/doppler.csv", results["doppler"])
results["quadchirp"] = repeatedcomparisons(generatesignals(:quadchirp,L), wt, 0.5, 2, samples, repeats)
CSV.write("./results/quadchirp.csv", results["quadchirp"])
results["mishmash"] = repeatedcomparisons(generatesignals(:mishmash,L), wt, 0.5, 2, samples, repeats)
CSV.write("./results/mishmash.csv", results["mishmash"])

# display results
vscodedisplay(results["blocks"])
vscodedisplay(results["bumps"])
vscodedisplay(results["heavysine"])
vscodedisplay(results["doppler"])
vscodedisplay(results["quadchirp"])
vscodedisplay(results["mishmash"])

max_stars_repo_name: zengfung/WaveletsDenoisingExperiment

max_issues_repo_name: zengfung/WaveletsDenoisingExperiment

max_forks_repo_name: zengfung/WaveletsDenoisingExperiment


text:
mutable struct EmptyRoomDirected{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Wall, Goal}}
    agent_pos::CartesianIndex{2}
    agent_dir::AbstractDirection
    reward::T
    rng::R
    terminal_reward::T
    goal_pos::CartesianIndex{2}
    done::Bool
end

@generate_getters(EmptyRoomDirected)
@generate_setters(EmptyRoomDirected)

mutable struct EmptyRoomUndirected{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Wall, Goal}}
    agent_pos::CartesianIndex{2}
    reward::T
    rng::R
    terminal_reward::T
    goal_pos::CartesianIndex{2}
    done::Bool
end

@generate_getters(EmptyRoomUndirected)
@generate_setters(EmptyRoomUndirected)

#####
# Directed
#####

function EmptyRoomDirected(; T = Float32, height = 8, width = 8, rng = Random.GLOBAL_RNG)
    objects = (AGENT, WALL, GOAL)
    world = GridWorldBase(objects, height, width)

    room = Room(CartesianIndex(1, 1), height, width)
    place_room!(world, room)

    goal_pos = CartesianIndex(height - 1, width - 1)
    world[GOAL, goal_pos] = true

    agent_pos = CartesianIndex(2, 2)
    world[AGENT, agent_pos] = true

    agent_dir = RIGHT
    reward = zero(T)
    terminal_reward = one(T)
    done = false

    env = EmptyRoomDirected(world, agent_pos, agent_dir, reward, rng, terminal_reward, goal_pos, done)

    RLBase.reset!(env)

    return env
end

RLBase.state_space(env::EmptyRoomDirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const EMPTY_ROOM_DIRECTED_LAYERS = SA.SVector(2, 3)
RLBase.state(env::EmptyRoomDirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_agent_dir(env), get_half_size(env), EMPTY_ROOM_DIRECTED_LAYERS)

RLBase.state_space(env::EmptyRoomDirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::EmptyRoomDirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = (copy(get_grid(env)), get_agent_dir(env))

RLBase.action_space(env::EmptyRoomDirected, ::RLBase.DefaultPlayer) = DIRECTED_NAVIGATION_ACTIONS
RLBase.reward(env::EmptyRoomDirected, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::EmptyRoomDirected) = get_done(env)

function RLBase.reset!(env::EmptyRoomDirected{T}) where {T}
    world = get_world(env)
    rng = get_rng(env)

    world[AGENT, get_agent_pos(env)] = false
    world[GOAL, get_goal_pos(env)] = false

    new_goal_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_goal_pos!(env, new_goal_pos)
    world[GOAL, new_goal_pos] = true

    new_agent_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_agent_pos!(env, new_agent_pos)
    world[AGENT, new_agent_pos] = true

    set_agent_dir!(env, rand(rng, DIRECTIONS))

    set_reward!(env, zero(T))
    set_done!(env, false)

    return nothing
end

function (env::EmptyRoomDirected{T})(action::AbstractTurnAction) where {T}
    new_dir = turn(action, get_agent_dir(env))
    set_agent_dir!(env, new_dir)
    world = get_world(env)

    if world[GOAL, get_agent_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

function (env::EmptyRoomDirected{T})(action::AbstractMoveAction) where {T}
    world = get_world(env)

    agent_pos = get_agent_pos(env)
    dest = move(action, get_agent_dir(env), agent_pos)
    if !world[WALL, dest]
        world[AGENT, agent_pos] = false
        set_agent_pos!(env, dest)
        world[AGENT, dest] = true
    end

    if world[GOAL, get_agent_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

#####
# Undirected
#####

function EmptyRoomUndirected(; T = Float32, height = 8, width = 8, rng = Random.GLOBAL_RNG)
    objects = (AGENT, WALL, GOAL)
    world = GridWorldBase(objects, height, width)

    room = Room(CartesianIndex(1, 1), height, width)
    place_room!(world, room)

    goal_pos = CartesianIndex(height - 1, width - 1)
    world[GOAL, goal_pos] = true

    agent_pos = CartesianIndex(2, 2)
    world[AGENT, agent_pos] = true

    reward = zero(T)
    terminal_reward = one(T)
    done = false

    env = EmptyRoomUndirected(world, agent_pos, reward, rng, terminal_reward, goal_pos, done)

    RLBase.reset!(env)

    return env
end

RLBase.state_space(env::EmptyRoomUndirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const EMPTY_ROOM_UNDIRECTED_LAYERS = SA.SVector(2, 3)
RLBase.state(env::EmptyRoomUndirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_half_size(env), EMPTY_ROOM_UNDIRECTED_LAYERS)

RLBase.state_space(env::EmptyRoomUndirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::EmptyRoomUndirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = copy(get_grid(env))

RLBase.action_space(env::EmptyRoomUndirected, player::RLBase.DefaultPlayer) = UNDIRECTED_NAVIGATION_ACTIONS
RLBase.reward(env::EmptyRoomUndirected, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::EmptyRoomUndirected) = get_done(env)

function RLBase.reset!(env::EmptyRoomUndirected{T}) where {T}
    world = get_world(env)
    rng = get_rng(env)

    world[AGENT, get_agent_pos(env)] = false
    world[GOAL, get_goal_pos(env)] = false

    new_goal_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_goal_pos!(env, new_goal_pos)
    world[GOAL, new_goal_pos] = true

    new_agent_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_agent_pos!(env, new_agent_pos)
    world[AGENT, new_agent_pos] = true

    set_reward!(env, zero(T))
    set_done!(env, false)

    return nothing
end

function (env::EmptyRoomUndirected{T})(action::AbstractMoveAction) where {T}
    world = get_world(env)

    agent_pos = get_agent_pos(env)
    dest = move(action, agent_pos)

    if !world[WALL, dest]
        world[AGENT, agent_pos] = false
        set_agent_pos!(env, dest)
        world[AGENT, dest] = true
    end

    if world[GOAL, get_agent_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

ext: jl


text:
if basename(pwd()) == "aoc"
    cd("2018/5")
end

part1(filename::AbstractString) = reaction(readchomp(filename))

function reaction(str::AbstractString)
    len = -1
    while length(str) != len
        len = length(str)
        str = replace(
            str,
            Iterators.flatten(
                (c * C => "", C * c => "")
                for (c, C) in zip('a':'z', 'A':'Z')
            )...
        )
    end
    len
end

@assert reaction("aA") == 0
@assert reaction("abBA") == 0
@assert reaction("abAB") == 4
@assert reaction("aabAAB") == 6
@assert part1("example.txt") == 10
part1("input.txt")

function part2(filename::AbstractString)
    str = readchomp(filename)
    mini = minimum(uppercase(str))
    maxi = maximum(lowercase(str))
    ranges = zip(lowercase(mini):maxi, mini:uppercase(maxi))
    minimum(reaction, replace(str, c => "", C => "") for (c, C) in ranges)
end

@assert part2("example.txt") == 4
part2("input.txt")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
function is_file_path_ok(path_to_file::String)::VLResult

    if (isfile(path_to_file) == false)

        # error message -
        error_message = "Ooops! $(path_to_file) does not exist."
        error_object = ArgumentError(error_message)

        # return -
        return VLResult{ArgumentError}(error_object)
    end

    # default: return nothing -
    return VLResult(nothing)
end

function is_dir_path_ok(path_to_file::String)::VLResult

    # get the dir of the path -
    dir_name = dirname(path_to_file)

    # check, is this a legit dir?
    if (isdir(dir_name) == false)

        # error message -
        error_message = "Ooops! $(dir_name) does not exist."
        error_object = ArgumentError(error_message)

        # return -
        return VLResult{ArgumentError}(error_object)
    end

    # default: return nothing -
    return VLResult(nothing)
end

function check(result::VLResult; logger::Union{Nothing,SimpleLogger} = nothing)::(Union{Nothing,T} where T<:Any)

    # ok, so check, do we have an error object?
    # Yes: log the error if we have a logger, then throw the error.
    # No: return the result.value

    # Error case -
    if (isa(result.value, Exception) == true)

        # get the error object -
        error_object = result.value

        # get the error message as a String -
        error_message = sprint(showerror, error_object)

        # log -
        if (isnothing(logger) == false)
            with_logger(logger) do
                @error(error_message)
            end
        end

        # throw -
        throw(result.value)
    end

    # default -
    return result.value
end

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
export AbstractNDSparse

import Dagger: chunktype, domain, tochunk, distribute,
               chunks, Context, compute, gather, free!

import IndexedTables: eltypes, astuple, colnames, ndsparse, pkeynames, valuenames

boundingrect(x::IndexSpace) = x.boundingrect
interval(x::IndexSpace) = x.interval

"""
    DNDSparse{K,V} <: AbstractNDSparse

A distributed [NDSparse](@ref) datastructure. Can be constructed by:
- [`ndsparse`](@ref) from Julia objects
- [`loadndsparse`](@ref) from data on disk
- [`distribute`](@ref) from an [`NDSparse`](@ref) object
"""
mutable struct DNDSparse{K,V} <: AbstractNDSparse
    domains::Vector{IndexSpace{K}}
    chunks::Vector
    freed::Bool
    function DNDSparse{K, V}(domains, chunks) where {K, V}
        x = new(domains, chunks, false)
        Dagger.refcount_chunks(x.chunks)
        finalizer(x, free!)
        x
    end
end
function free!(x::DNDSparse)
    if !x.freed
        @schedule Dagger.free_chunks(x.chunks)
        x.freed = true
    end
    nothing
end

const DDataset = Union{DNextTable, DNDSparse}

function ndsparse(::Val{:distributed}, ks::Tup,
                  vs::Union{Tup, AbstractArray};
                  agg=nothing,
                  allowoverlap=agg === nothing,
                  closed=true, chunks=nothing, kwargs...)

    if chunks === nothing
        # this means the vectors are distributed.
        # pick the first distributed vector and distribute
        # all others similarly
        i = findfirst(x->isa(x, ArrayOp), ks)
        if i != 0
            darr = ks[i]
        elseif vs isa Tup && any(x->isa(x, ArrayOp, vs))
            darr = vs[findfirst(x->isa(x, ArrayOp), vs)]
        elseif isa(vs, ArrayOp)
            darr = vs
        else
            error("Don't know how to distribute. specify `chunks`")
        end

        chunks = domainchunks(compute(darr))
    end
    kdarrays = map(x->distribute(x, chunks), ks)
    vdarrays = isa(vs, Tup) ? map(x->distribute(x, chunks), vs) : distribute(vs, chunks)

    if isempty(kdarrays)
        error("NDSparse must be constructed with at least one index column")
    end

    nchunks = length(kdarrays[1].chunks)
    inames = isa(ks, NamedTuple) ? fieldnames(ks) : nothing
    ndims = length(ks)
    dnames = isa(vs, NamedTuple) ? fieldnames(vs) : nothing
    iscols = isa(vs, Tup)

    function makechunk(args...)
        k = Columns(args[1:ndims]..., names=inames)
        v = iscols ? Columns(args[ndims+1:end]..., names=dnames) : args[end]
        ndsparse(k,v; agg=agg, kwargs...)
    end

    cs = Array{Any}(nchunks)
    for i = 1:nchunks
        args = Any[map(x->x.chunks[i], kdarrays)...]
        append!(args, isa(vs, Tup) ? [map(x->x.chunks[i], vdarrays)...] :
                [vdarrays.chunks[i]])
        cs[i] = delayed(makechunk)(args...)
    end
    fromchunks(cs, closed=closed,
               merge=(x,y)->merge(x,y, agg=agg),
               allowoverlap=allowoverlap)
end

function ndsparse(x::Dagger.DArray{<:Tup}, y; kwargs...)
    ndsparse(columns(x), y; kwargs...)
end

function ndsparse(x::DNDSparse; kwargs...)
    ndsparse(columns(x, pkeynames(x)), values(x); kwargs...)
end
Base.eltype(dt::DNDSparse{K,V}) where {K,V} = V
Base.keytype(dt::DNDSparse{K,V}) where {K,V} = IndexedTables.astuple(K)
IndexedTables.dimlabels(dt::DNDSparse{K}) where {K} = fieldnames(K)
Base.ndims(dt::DNDSparse{K}) where {K} = nfields(K)
keytype(dt::DNDSparse{K}) where {K} = astuple(K)

# TableLike API
Base.@pure function IndexedTables.colnames(t::DNDSparse{K,V}) where {K,V}
    dnames = V<:Tup ? fieldnames(V) : [1]
    if all(x->isa(x, Integer), dnames)
        dnames = map(x->x+nfields(K), dnames)
    end
    vcat(fieldnames(K), dnames)
end

Base.@pure pkeynames(t::DNDSparse{K,V}) where {K, V} = (fieldnames(K)...)
function IndexedTables.valuenames(t::DNDSparse{K,V}) where {K,V}
    if V <: Tup
        if V<:NamedTuple
            (fieldnames(V)...)
        else
            ((ndims(t) + (1:nfields(V)))...)
        end
    else
        ndims(t) + 1
    end
end

const compute_context = Ref{Union{Void, Context}}(nothing)
get_context() = compute_context[] == nothing ? Context() : compute_context[]

"""
    compute(t::DNDSparse; allowoverlap, closed)

Computes any delayed-evaluations in the `DNDSparse`.
The computed data is left on the worker processes.
Subsequent operations on the results will reuse the chunks.

If `allowoverlap` is false then the computed data is re-sorted if required to have no
chunks with overlapping index ranges if necessary.

If `closed` is true then the computed data is re-sorted if required to have no
chunks with overlapping OR continuous boundaries.

See also [`collect`](@ref).

!!! warning
    `compute(t)` requires at least as much memory as the size of the
    result of the computing `t`. You usually don't need to do this for the whole dataset.
    If the result is expected to be big, try `compute(save(t, "output_dir"))` instead.
    See [`save`](@ref) for more.
"""
compute(t::DNDSparse; kwargs...) = compute(get_context(), t; kwargs...)

function compute(ctx, t::DNDSparse; allowoverlap=true, closed=false)
    if any(Dagger.istask, t.chunks)
        # we need to splat `thunks` so that Dagger knows the inputs
        # are thunks and they need to be staged for scheduling
        vec_thunk = delayed((refs...) -> [refs...]; meta=true)(t.chunks...)
        cs = compute(ctx, vec_thunk) # returns a vector of Chunk objects
        t1 = fromchunks(cs, allowoverlap=allowoverlap, closed=closed)
        compute(t1)
    else
        map(Dagger.unrelease, t.chunks) # don't let this be freed
        foreach(Dagger.persist!, t.chunks)
        t
    end
end

"""
    collect(t::DNDSparse)

Gets distributed data in a DNDSparse `t` and merges it into
[NDSparse](#IndexedTables.NDSparse) object

!!! warning
    `collect(t)` requires at least as much memory as the size of the
    result of the computing `t`. If the result is expected to be big,
    try `compute(save(t, "output_dir"))` instead. See [`save`](@ref) for more.
    This data can be loaded later using [`load`](@ref).
"""
collect(t::DNDSparse) = collect(get_context(), t)
function Base.copy(x::DNDSparse{K,V}) where {K, V}
    fromchunks(NDSparse,
               delayedmap(copy, x.chunks),
               domains=x.domains,
               KV=(K,V))
end

function collect(ctx::Context, dt::DNDSparse{K,V}) where {K,V}
    cs = dt.chunks
    if length(cs) > 0
        collect(ctx, treereduce(delayed(_merge), cs))
    else
        ndsparse(similar(IndexedTables.arrayof(K), 0),
                 similar(IndexedTables.arrayof(V), 0))
    end
end

# Fast-path merge if the data don't overlap
function _merge(f, a::NDSparse, b::NDSparse)
    if isempty(a)
        b
    elseif isempty(b)
        a
    elseif last(a.index) < first(b.index)
        # can vcat
        NDSparse(vcat(a.index, b.index), vcat(a.data, b.data),
              presorted=true, copy=false)
    elseif last(b.index) < first(a.index)
        _merge(b, a)
    else
        f(a, b) # Keep equal index elements
    end
end

_merge(f, x::NDSparse) = x
function _merge(f, x::NDSparse, y::NDSparse, ys::NDSparse...)
    treereduce((a,b)->_merge(f, a, b), [x,y,ys...])
end

_merge(x::NDSparse, y::NDSparse...) = _merge((a,b) -> merge(a, b, agg=nothing), x, y...)

"""
    map(f, t::DNDSparse)

Applies a function `f` on every element in the data of table `t`.
"""
Base.map(f, dt::DNDSparse) = mapchunks(c->map(f, c), dt)

struct EmptySpace{T} end

# Teach dagger how to automatically figure out the
# metadata (in dagger parlance "domain") about an NDSparse chunk.
function Dagger.domain(nd::NDSparse)
    T = eltype(keys(nd))

    if isempty(nd)
        return EmptySpace{T}()
    end

    wrap = T<:NamedTuple ? T : tuple

    interval = Interval(wrap(first(nd.index)...), wrap(last(nd.index)...))
 #  cs = astuple(nd.index.columns)
 #  extr = map(extrema, cs[2:end]) # we use first and last value of first column
 #  boundingrect = Interval(wrap(first(cs[1]), map(first, extr)...),
 #                          wrap(last(cs[1]), map(last, extr)...))
    return IndexSpace(interval, interval, Nullable{Int}(length(nd)))
end

function subindexspace(t::Union{NDSparse, NextTable}, r)
    ks = pkeys(t)
    T = eltype(typeof(ks))
    wrap = T<:NamedTuple ? T : tuple

    if isempty(r)
        return EmptySpace{T}()
    end

    interval = Interval(wrap(ks[first(r)]...), wrap(ks[last(r)]...))
    cs = astuple(columns(ks))
    extr = map(c -> extrema_range(c, r), cs[2:end])

    boundingrect = Interval(wrap(cs[1][first(r)], map(first, extr)...),
                            wrap(cs[1][last(r)], map(last, extr)...))
    return IndexSpace(interval, boundingrect, Nullable{Int}(length(r)))
end

Base.eltype(::IndexSpace{T}) where {T} = T
Base.eltype(::EmptySpace{T}) where {T} = T

Base.isempty(::EmptySpace) = true
Base.isempty(::IndexSpace) = false

nrows(td::IndexSpace) = td.nrows
nrows(td::EmptySpace) = Nullable(0)

Base.ndims(::IndexSpace{T}) where {T}  = nfields(T)
Base.ndims(::EmptySpace{T}) where {T}  = nfields(T)

Base.first(td::IndexSpace) = first(td.interval)
Base.last(td::IndexSpace) = last(td.interval)

mins(td::IndexSpace) = first(td.boundingrect)
maxes(td::IndexSpace) = last(td.boundingrect)

function Base.merge(d1::IndexSpace, d2::IndexSpace, collisions=true)
    n = collisions || isnull(d1.nrows) || isnull(d2.nrows) ?
        Nullable{Int}() :
        Nullable(get(d1.nrows) + get(d2.nrows))

    interval = merge(d1.interval, d2.interval)
    boundingrect = boxmerge(d1.boundingrect, d2.boundingrect)
    IndexSpace(interval, boundingrect, n)
end
Base.merge(d1::IndexSpace, d2::EmptySpace) = d1
Base.merge(d1::EmptySpace, d2::Union{IndexSpace, EmptySpace}) = d2

function Base.intersect(d1::IndexSpace, d2::IndexSpace)
    interval = intersect(d1.interval, d2.interval)
    boundingrect = intersect(d1.boundingrect, d2.boundingrect)
    IndexSpace(interval, boundingrect, Nullable{Int}())
end

function Base.intersect(d1::EmptySpace, d2::Union{IndexSpace,EmptySpace})
    d1
end

# given a chunks index constructed above, give an array of
# index spaces spanned by the chunks in the index
function index_spaces(t::NDSparse)
    intervals = map(x-> Interval(map(first, x), map(last, x)), t.index)
    boundingrects = map(x-> Interval(map(first, x), map(last, x)), t.data.columns.boundingrect)
    map(IndexSpace, intervals, boundingrects, t.data.columns.length)
end

"""
The length of the `DNDSparse` if it can be computed. Will throw an error if not.
You can get the length of such tables after calling [`compute`](@ref) on them.
"""
function Base.length(t::DNDSparse)
    l = trylength(t)
    if isnull(l)
        error("The length of the DNDSparse is not yet known since some of its parts are not yet computed. Call `compute` to compute them, and then call `length` on the result of `compute`.")
    else
        get(l)
    end
end

function _has_overlaps(firsts, lasts, closed)
    if !issorted(firsts)
        p = sortperm(firsts)
        firsts = firsts[p]
        lasts = lasts[p]
    end

    for i = 1:length(firsts)
        s_i = firsts[i]
        j = searchsortedfirst(lasts, s_i)

        # allow repeated indices between chunks
        if j != i && j <= length(lasts) && isless(s_i, lasts[j])
            return true
        elseif closed && j != i && j <= length(lasts) && s_i == lasts[j]
            return true
        end
    end
    return false
end

function has_overlaps(domains; closed=true)
    _has_overlaps(first.(domains), last.(domains), closed)
end

function has_overlaps(domains, dims::AbstractVector)
    sub(x) = x[dims]
    fs = sub.(first.(domains))
    ls = sub.(last.(domains))
    _has_overlaps(fs, ls, true)
end

function with_overlaps(f, t::DDataset, closed=false)
    domains = t.domains
    chunks = t.chunks

    if isempty(domains)
        return t
    end

    if !issorted(domains, by=first)
        perm = sortperm(domains, by = first)
        domains = domains[perm]
        chunks = chunks[perm]
    end

    stack = [domains[1]]
    groups = [[1]]
    for i in 2:length(domains)
        sub = domains[i]
        if hasoverlap(stack[end].interval, sub.interval)
            stack[end] = merge(stack[end], sub)
            push!(groups[end], i)
        else
            push!(stack, sub)
            push!(groups, [i])
        end
    end

    cs = collect([f(chunks[group]) for group in groups])
    fromchunks(cs)
end

function fromchunks(::Type{<:NDSparse}, chunks::AbstractArray;
                    domains::AbstractArray = map(domain, chunks),
                    KV = getkvtypes(chunks),
                    merge=_merge,
                    closed = false,
                    allowoverlap = true)

    nzidxs = find(x->!isempty(x), domains)
    domains = domains[nzidxs]

    dt = DNDSparse{KV...}(domains, chunks[nzidxs])

    if !allowoverlap && has_overlaps(domains, closed=closed)
        return rechunk(dt, closed=closed, merge=merge)
    else
        return dt
    end
end

function cache_thunks(t)
    for c in t.chunks
        if isa(c, Dagger.Thunk)
            Dagger.cache_result!(c)
        end
    end
    t
end

function getkvtypes(::Type{N}) where N<:NDSparse
    eltype(N.parameters[3]), N.parameters[1]
end

_promote_type(T,S) = promote_type(T,S)
_promote_type(T::Type{<:IndexTuple}, S::Type{<:IndexTuple}) = map_params(_promote_type, T, S)

function getkvtypes(xs::AbstractArray)
    kvtypes = getkvtypes.(chunktype.(xs))
    K, V = kvtypes[1]
    for (Tk, Tv) in kvtypes[2:end]
        K = _promote_type(Tk, K)
        V = _promote_type(Tv, V)
    end
    (K, V)
end

### Distribute a NDSparse into a DNDSparse

"""
    distribute(itable::NDSparse, rowgroups::AbstractArray)

Distributes an NDSparse object into a DNDSparse by splitting it up into chunks
of `rowgroups` elements. `rowgroups` is a vector specifying the number of
rows in the chunks.

Returns a `DNDSparse`.
"""
function distribute(nds::NDSparse{V}, rowgroups::AbstractArray;
                     allowoverlap = true, closed = false) where V
    splits = cumsum([0, rowgroups;])

    if splits[end] != length(nds)
        throw(ArgumentError("the row groups don't add up to total number of rows"))
    end

    ranges = map(UnitRange, splits[1:end-1].+1, splits[2:end])

    # this works around locality optimizations in Dagger to make
    # sure that the parts get distributed instead of being left on
    # the master process - which would lead to all operations being serial.
    chunks = map(r->delayed(identity)(subtable(nds, r)), ranges)
    domains = map(r->subindexspace(nds, r), ranges)
    realK = eltypes(typeof(nds.index.columns))
    cache_thunks(fromchunks(chunks, domains=domains, KV = (realK, V),
                            allowoverlap=allowoverlap, closed=closed))
end

distribute(t::DNDSparse, cs) = ndsparse(t, chunks=cs)

"""
    distribute(itable::NDSparse, nchunks::Int=nworkers())

Distributes an NDSparse object into a DNDSparse of `nchunks` chunks
of approximately equal size.

Returns a `DNDSparse`.
"""
function distribute(nds::NDSparse, nchunks::Int=nworkers();
                    allowoverlap = true, closed = false)
    N = length(nds)
    q, r = divrem(N, nchunks)
    nrows = vcat(collect(_repeated(q, nchunks)))
    nrows[end] += r
    distribute(nds, nrows; allowoverlap=allowoverlap, closed=closed)
end

"""
    mapchunks(f, t::DNDSparse; keeplengths=true)

Applies a function to each chunk in `t`. Returns a new DNDSparse.
If `keeplength` is false, this means that the lengths of the output
chunks is unknown before [`compute`](@ref). This function is used
internally by many DNDSparse operations.
"""
function mapchunks(f, t::DNDSparse{K,V}; keeplengths=true) where {K,V}
    chunks = map(delayed(f), t.chunks)
    if keeplengths
        DNDSparse{K, V}(t.domains, chunks)
    else
        DNDSparse{K, V}(map(null_length, t.domains), chunks)
    end
end

function null_length(x::IndexSpace)
    IndexSpace(x.interval, x.boundingrect, Nullable{Int}())
end

function subtable(t::DNDSparse{K,V}, idx) where {K, V}
    if isnull(trylength(t))
        t = compute(t)
    end
    if isempty(idx)
        return DNDSparse{K, V}(similar(t.domains, 0), similar(t.chunks, 0))
    end
    ls = map(x->get(nrows(x)), t.domains)
    cumls = cumsum(ls)
    i = searchsortedlast(cumls, first(idx))
    j = searchsortedfirst(cumls, last(idx))

    # clip first and last chunks
    strt = first(idx) - get(cumls, i-1, 0)
    fin = cumls[j] - last(idx)

    ds = t.domains[i:j]
    cs = convert(Vector{Any}, t.chunks[i:j])
    if i==j
        cs[1] = delayed(x->subtable(x, strt:length(x)-fin))(cs[1])
        i = ds[1]
        ds[1] = IndexSpace(i.interval, i.boundingrect, Nullable(length(idx)))
    else
        cs[1] = delayed(x->subtable(x, strt:length(x)))(cs[1])
        cs[end] = delayed(x->subtable(x, 1:length(x)-fin))(cs[end])
        ds[1] = null_length(ds[1])
        ds[end] = null_length(ds[2])
    end

    DNDSparse{K,V}(ds, cs)
end

import IndexedTables: showtable

function Base.show(io::IO, big::DNDSparse)
    h, w = displaysize(io)
    showrows = h - 5 # This will trigger an ellipsis when there's
                     # more to see than the screen fits
    t = first(Iterators.partition(big, showrows))
    if !(values(t) isa Columns)
        cnames = colnames(keys(t))
        eltypeheader = "$(eltype(t))"
    else
        cnames = colnames(t)
        nf = nfields(eltype(t))
        if eltype(t) <: NamedTuple
            eltypeheader = "$(nf) field named tuples"
        else
            eltypeheader = "$(nf)-tuples"
        end
    end
    len = trylength(big)
    vals = isnull(len) ? "of" : "with $(get(len)) values"
    header = "$(ndims(t))-d Distributed NDSparse $vals ($eltypeheader) in $(length(big.chunks)) chunks:"
    showtable(io, t; header=header, cnames=cnames,
              divider=ndims(t), ellipsis=:end)
end

Base.@deprecate_binding DTable DNDSparse


text:
using LinearAlgebra
"""
mag_tor = magneric_torque(i_m, B)

磁気トルカによるトルクを計算

# Arguments
- 'i_m'：各磁気トルカに流れる電流
- 'B'：磁束密度ベクトル＠SCSF

# Return
- 'mag_tor'：磁気トルカによる原点周りのトルク

"""
function magnetic_torque(i_m, B)
    μ = 5000 #比透磁率
    n = 400 #巻き数
	L = 80  #コア長さ
	D = 10  #コア直径
	m_r = [0., 0., 0.]#[0.01, 0.01, 0.01] #残留磁気モーメント
	m = zeros(4,3)
	T = zeros(4,3)
	mag_tor = zeros(3)

	S = π*(D/1000)^2/4  # 断面積 [m]
	p = L/D
	μeff = 1 / (1/μ + (log(p)-1)/(p^2))  # コア実効比透磁率

	#磁気モーメント計算

    m[1,:] = μeff*n*S*i_m[1]*[1,0,0]
    m[2,:] = μeff*n*S*i_m[2]*[0,1,0]
	m[3,:] = μeff*n*S*i_m[3]*[0,0,1]
	m[4,:] = m_r

	#各軸周りのトルク計算
	for i=1:4
		T[i,:] = cross(m[i,:], B)
	end

	#全トルク計算
	for i=1:3
		for j=1:4
			global mag_tor[i] += T[j, i]
		end
	end

	return mag_tor
end

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
@testset "geometry.cartesian" begin

@testset "point2d" begin
    coords = float.([1 2; 3 4; 5 6; 7 8])
    p1 = Point2D(coords, 3)
    p2 = Point2D(1.2, -4.6)
    @test p1.x == 5
    @test p2.y == -4.6
    @test isapprox(toarray(p1 + p2), [6.2, 1.4])
    @test isapprox(toarray(p1 - p2), [3.8, 10.6])
    @test isapprox(toarray(p2 * 2.0), [2.4, -9.2])
    @test isapprox(cross2d(p1, Point2D(coords, 4)), -2)

    v = Point2D(3.0, 0.0)
    @test isapprox(interiorangle(v, Point2D(5.0, 0.0)), 0)
    @test isapprox(interiorangle(v, Point2D(-20.0, 0.0)), pi)
    @test isapprox(interiorangle(v, Point2D(3.0, -3sqrt(3))), pi / 3)
    @test isapprox(interiorangle(v, Point2D(-3.0, 3sqrt(3))), 2pi / 3)
    @test isnan(interiorangle(v, Point2D(0.0, 0.0)))
end

@testset "segment2d" begin
    coords = float.([1 2; 3 4; 5 6; 7 8])
    seg1 = Segment{Point2D}(coords, 2, 4)
    seg2 = Segment(Point2D(1.2, 4.5), Point2D(3.4, 7.9))
    @test toarray(seg1.u) == [3, 4]
    @test toarray(seg2.v) == [3.4, 7.9]
    @test seg2.u.x == 1.2
    @test seg1.u.y == 4
    @test seg1.v.x == 7
    @test seg2.v.y == 7.9

    @test isapprox(toarray(seg2.v - seg2.u), [2.2, 3.4])
    @test isapprox(toarray(midpoint(seg2)), [2.3, 6.2])
    tl = translate(seg1, -pi / 2, sqrt(2) / 2)
    @test isapprox(toarray(tl), [3.5 3.5; 7.5 7.5])

    seg3 = Segment(Point2D(1, 1), Point2D(5, 5))
    @test isapprox(toarray(trim_u(seg3, 1 / 4)), [2 2; 5 5])
    @test isapprox(toarray(trim_v(seg3, 1 / 4)), [1 1; 4 4])
    @test isapprox(toarray(trim_uv(seg3, 1 / 4)), [1.5 1.5; 4.5 4.5])

    @test isapprox(cross2d(seg1.u, seg1.v), -4)
end

@testset "cartesian2d" begin
    coords = float.([1 2; 3 4; 5 6; 7 8])
    @test sum(x_components(coords)) == 16
    @test sum(y_components(coords)) == 20
    seg = toarray(coords, 3, 4)
    @test seg isa SubArray
end

@testset "cartesian3d" begin
    coords = float.([1 2 3; 4 5 6; 7 8 9; 10 11 12])
    @test sum(x_components(coords)) == 22
    @test sum(y_components(coords)) == 26
    @test sum(z_components(coords)) == 30
end

@testset "isclockwise" begin
    @test isclockwise(float.([0 1; 1 0; 0 -1; -1 0]))
    @test !isclockwise(float.([0 1; -1 0; 0 -1; 1 0]))
    @test isclockwise(float.([1 1; 1 -1; -1 1; -1 -1])) === nothing
    @test isclockwise(float.([-1 1; 1 1; 1 2; 2 2; 2 -2; -2 -2; -2 2; -1 2]))
end

end # geometry.cartesian


text:
@cache struct SymplecticEulerCache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
end

function alg_cache(alg::SymplecticEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  SymplecticEulerCache(u,uprev,similar(u),zero(rate_prototype),zero(rate_prototype))
end

struct SymplecticEulerConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::SymplecticEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = SymplecticEulerConstantCache()

@cache struct VelocityVerletCache{uType,rateType,uEltypeNoUnits} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  half::uEltypeNoUnits
end

struct VelocityVerletConstantCache{uEltypeNoUnits} <: OrdinaryDiffEqConstantCache
    half::uEltypeNoUnits
end

function alg_cache(alg::VelocityVerlet,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = zero(rate_prototype)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  half = uEltypeNoUnits(1//2)
  VelocityVerletCache(u,uprev,k,tmp,fsalfirst,half)
end

alg_cache(alg::VelocityVerlet,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = VelocityVerletConstantCache(uEltypeNoUnits(1//2))

@cache struct Symplectic2Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::VerletLeapfrog,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = VerletLeapfrogConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  Symplectic2Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::VerletLeapfrog,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) =
      VerletLeapfrogConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

function alg_cache(alg::PseudoVerletLeapfrog,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = PseudoVerletLeapfrogConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  Symplectic2Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::PseudoVerletLeapfrog,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) =
PseudoVerletLeapfrogConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

function alg_cache(alg::McAte2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = McAte2ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  Symplectic2Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::McAte2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) =
McAte2ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

@cache struct Symplectic3Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::Ruth3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = Ruth3ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  Symplectic3Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::Ruth3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) =
Ruth3ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

function alg_cache(alg::McAte3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = McAte3ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  Symplectic3Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::McAte3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) =
McAte3ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

@cache struct Symplectic4Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::McAte4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = McAte4ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  Symplectic4Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::McAte4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = McAte4ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

function alg_cache(alg::CandyRoz4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = CandyRoz4ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  Symplectic4Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::CandyRoz4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = McAte4ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

@cache struct Symplectic45Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::CalvoSanz4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = CalvoSanz4ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  Symplectic45Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::CalvoSanz4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = CalvoSanz4ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

function alg_cache(alg::McAte42,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = McAte42ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  Symplectic45Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::McAte42,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = McAte42ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

@cache struct Symplectic5Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::McAte5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = McAte5ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  Symplectic5Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::McAte5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = McAte5ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

@cache struct Symplectic6Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::Yoshida6,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = Yoshida6ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  Symplectic6Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::Yoshida6,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Yoshida6ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

@cache struct Symplectic62Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::KahanLi6,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = KahanLi6ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  Symplectic62Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::KahanLi6,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = KahanLi6ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

@cache struct McAte8Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::McAte8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = McAte8ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  McAte8Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::McAte8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = McAte8ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

@cache struct KahanLi8Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::KahanLi8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = KahanLi8ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  KahanLi8Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::KahanLi8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = KahanLi8ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

@cache struct SofSpa10Cache{uType,rateType,tableauType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  tab::tableauType
end

function alg_cache(alg::SofSpa10,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = SofSpa10ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  SofSpa10Cache(u,uprev,k,tmp,fsalfirst,tab)
end

alg_cache(alg::SofSpa10,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = SofSpa10ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))


text:
using Mux
using Test
using HTTP
import HTTP.ExceptionRequest: StatusError

@testset "Mux" begin

@testset "misc" begin
  function f()
    @app foo = (Mux.defaults)
  end

  @test f() == nothing

  @test Mux.notfound()(Dict())[:status] == 404
end

@testset "basic server" begin
  d1 = Dict("one"=> "1", "two"=> "2")
  d2 = Dict("one"=> "1", "two"=> "")
  @app test = (
    Mux.defaults,
    page("/",respond("
Hello World!
")),
    page("/about", respond("
Boo!
")),
    page("/user/:user", req -> "
Hello, $(req[:params][:user])!
"),
    query(d1, respond("
query1
")),
    query(d2, respond("
query2
")),
    Mux.notfound())
  serve(test)

  @testset "page" begin
    @test String(HTTP.get("http://localhost:8000").body) ==
                "
Hello World!
"
    @test String(HTTP.get("http://localhost:8000/about").body) ==
                "
Boo!
"
    @test String(HTTP.get("http://localhost:8000/user/julia").body) ==
                "
Hello, julia!
"
  end

  @testset "query" begin
    @test String(HTTP.get("http://localhost:8000/dum?one=1&two=2").body) ==
                "
query1
"
    @test_throws StatusError String(HTTP.get("http://localhost:8000/dum?one=1").body)
    @test_throws StatusError String(HTTP.get("http://localhost:8000/dum?one=1&two=2&sarv=boo").body)
    @test_throws StatusError String(HTTP.get("http://localhost:8000/dum?one=1").body)

    @test String(HTTP.get("http://localhost:8000/dum?one=1&two=56").body) ==
                "
query2
"
    @test String(HTTP.get("http://localhost:8000/dum?one=1&two=hfjd").body) ==
                "
query2
"
    @test_throws StatusError String(HTTP.get("http://localhost:8000/dum?one=1").body)
    @test_throws StatusError String(HTTP.get("http://localhost:8000/dum?one=1&two=2&sarv=boo").body)
  end
end

@testset "MIME types" begin
  # Issue #68
  @test Mux.fileheaders("foo.css")["Content-Type"] == "text/css"
  @test Mux.fileheaders("foo.html")["Content-Type"] == "text/html"
  @test Mux.fileheaders("foo.js")["Content-Type"] == "application/javascript"
end

# Check that prod_defaults don't completely break things
# And check prod_defaults error handler
@testset "prod defaults" begin
  throwapp() = (_...) -> error("An error!")

  # Used for wrapping stderrcatch so we can check its output and stop it spewing
  # all over the test results.
  path, mock_stderr = mktemp()

  @app test = (
    Mux.prod_defaults,
    (app, req) -> redirect_stderr(() -> Mux.stderrcatch(app, req), mock_stderr),
    page("/",respond("
Hello World!
")),
    page("/about", respond("
Boo!
")),
    page("/user/:user", req -> "
Hello, $(req[:params][:user])!
"),
    throwapp(),
    Mux.notfound())

  serve(test, 8001)

  @test String(HTTP.get("http://localhost:8001").body) ==
              "
Hello World!
"
  @test String(HTTP.get("http://localhost:8001/about").body) ==
              "
Boo!
"
  @test String(HTTP.get("http://localhost:8001/user/julia").body) ==
              "
Hello, julia!
"
  @test String(HTTP.get("http://localhost:8001/badurl";
                        status_exception=false).body) ==
               "Internal server error"

  # Check our error was logged and close fake stderr
  seekstart(mock_stderr)
  @test occursin("An error!", read(mock_stderr, String))
  close(mock_stderr)
  rm(path)
end

# Test page and route are callable without a string argument
# (previously the first two raised StackOverflowError)
@testset "bare page()" begin
  @test page(identity, identity) isa Function
  @test route(identity, identity) isa Function
  @test page(identity) isa Function
  @test route(identity) isa Function

  # Test you can pass the string last if you really want.
  @test page(identity, "") isa Function
  @test route(identity, "") isa Function
end

@testset "WebSockets" begin
  import Mux.WebSockets

  @app h = (
    Mux.defaults,
    page("/", respond("
Hello World!
")),
    Mux.notfound());

  @app w = (
    Mux.wdefaults,
    route("/ws_io", Mux.echo),
    Mux.wclose,
    Mux.notfound());

  serve(h, w, 2333)

  WebSockets.open("ws://localhost:2333/ws_io") do ws_client
    message = "Hello WebSocket!"
    WebSockets.writeguarded(ws_client, message)
    data, success = WebSockets.readguarded(ws_client)
    @test success
    @test String(data) == message
  end
end

@testset "Secure WebSockets" begin
  import Mux.WebSockets

  @app h = (
    Mux.defaults,
    page("/", respond("
Hello World!
")),
    Mux.notfound());

  @app w = (
    Mux.wdefaults,
    route("/ws_io", Mux.echo),
    Mux.wclose,
    Mux.notfound());

  cert = joinpath(@__DIR__, "test.cert")
  key = joinpath(@__DIR__, "test.key")
  #serve(h, w, "127.0.0.1", 2444; sslconfig=WebSockets.SSLConfig(cert, key))
  serve(h, w, 2444; sslconfig=WebSockets.SSLConfig(cert, key))

  WebSockets.open("wss://localhost:2444/ws_io"; sslconfig=WebSockets.SSLConfig(false)) do ws_client
    message = "Hello WebSocket!"
    WebSockets.writeguarded(ws_client, message)
    data, success = WebSockets.readguarded(ws_client)
    @test success
    @test String(data) == message
  end
end

end


text:
module Helium_HF_FEM_Eigen
    include("gausslegendre.jl")
    include("helium_hf_fem_eigen_module.jl")
    include("helium_vh_fem_module.jl")
    using LinearAlgebra
    using Match
    using MKL
    using .GaussLegendre
    using .Helium_HF_FEM_Eigen_module
    using .Helium_Vh_FEM_module

    const MAXR = 50.0
    const MINR = 0.0
    const NODE_TOTAL = 5000

    function construct()
        hfem_param = Helium_HF_FEM_Eigen_module.Helium_HF_FEM_Eigen_param(NODE_TOTAL - 1, NODE_TOTAL, MAXR, MINR)
        hfem_val = Helium_HF_FEM_Eigen_module.Helium_HF_FEM_Eigen_variables(
            Symmetric(zeros(hfem_param.ELE_TOTAL, hfem_param.ELE_TOTAL)),
            Array{Float64}(undef, hfem_param.ELE_TOTAL),
            Array{Float64, 3}(undef, hfem_param.ELE_TOTAL, 2, 2),
            Array{Float64, 3}(undef, hfem_param.ELE_TOTAL, 2, 2),
            Array{Int64, 2}(undef, hfem_param.ELE_TOTAL, 2),
            Array{Float64, 2}(undef, hfem_param.ELE_TOTAL, 2),
            Array{Float64}(undef, hfem_param.NODE_TOTAL),
            Array{Float64}(undef, hfem_param.NODE_TOTAL),
            Symmetric(zeros(hfem_param.ELE_TOTAL, hfem_param.ELE_TOTAL)),
            Array{Float64}(undef, hfem_param.NODE_TOTAL))

        return hfem_param, hfem_val
    end

    function make_wavefunction(iter, hfem_param, hfem_val, vh_val)
        if iter == 0
            # データの生成
            make_data!(hfem_param, hfem_val)

            return nothing
        end

        # 要素行列の生成
        make_element_matrix!(hfem_param, hfem_val, vh_val)

        # 全体行列を生成
        hg_tmp, ug_tmp = make_global_matrix(hfem_param, hfem_val)

        # 境界条件処理を行う
        boundary_conditions!(hfem_param, hfem_val, hg_tmp, ug_tmp)

        # 一般化固有値問題を解く
        eigenval, phi = eigen!(hfem_val.hg, hfem_val.ug)

        # 基底状態の固有ベクトルを取り出す
        hfem_val.phi = @view(phi[:,1])

        # 固有ベクトルの要素数を増やす
        resize!(hfem_val.phi, NODE_TOTAL)

        # 端点rcでの値を0にする
        hfem_val.phi[NODE_TOTAL] = 0.0

        # 固有ベクトル（波動関数）を規格化
        normalize!(hfem_val)

        return eigenval[1]
    end

    rVh(hfem_param, hfem_val, vh_val, r) = let
        klo = 1
        max = hfem_param.NODE_TOTAL
        khi = max

        # 表の中の正しい位置を二分探索で求める
        @inbounds while khi - klo > 1
            k = (khi + klo) >> 1

            if hfem_val.node_r_glo[k] > r
                khi = k
            else
                klo = k
            end
        end

        # yvec_[i] = f(xvec_[i]), yvec_[i + 1] = f(xvec_[i + 1])の二点を通る直線を代入
        return (vh_val.ug[khi] - vh_val.ug[klo]) / (hfem_val.node_r_glo[khi] - hfem_val.node_r_glo[klo]) * (r - hfem_val.node_r_glo[klo]) + vh_val.ug[klo]
    end

    function boundary_conditions!(hfem_param, hfem_val, hg_tmp, ug_tmp)
        hfem_val.hg = Symmetric(zeros(hfem_param.ELE_TOTAL, hfem_param.ELE_TOTAL))
        hfem_val.ug = Symmetric(zeros(hfem_param.ELE_TOTAL, hfem_param.ELE_TOTAL))

        @inbounds for i = 1:hfem_param.ELE_TOTAL
            for j = i - 1:i + 1
                if j != 0 && j != hfem_param.NODE_TOTAL
                    hfem_val.hg.data[j, i] = hg_tmp.data[j, i]
                    hfem_val.ug.data[j, i] = ug_tmp.data[j, i]
                end
            end
        end
    end

    get_A_matrix_element(e, le, p, q, hfem_param, hfem_val, vh_val) = let
        ed = float(e - 1)
        @match p begin
            1 =>
                @match q begin
                    1 => return 0.5 * le * (ed * ed + ed + 1.0 / 3.0) - le * le * (2.0 * ed / 3.0 + 1.0 / 6.0) +
                                GaussLegendre.gl_integ(r -> r * (hfem_val.node_r_ele[e, 2] - r) / hfem_val.length[e] * rVh(hfem_param, hfem_val, vh_val, r) * (hfem_val.node_r_ele[e, 2] - r) / hfem_val.length[e],
                                                       hfem_val.node_r_ele[e, 1],
                                                       hfem_val.node_r_ele[e, 2],
                                                       vh_val)
                    2 => return -0.5 * le * (ed * ed + ed + 1.0 / 3.0) - le * le * (ed / 3.0 + 1.0 / 6.0) +
                                GaussLegendre.gl_integ(r -> r * (hfem_val.node_r_ele[e, 2] - r) / hfem_val.length[e] * rVh(hfem_param, hfem_val, vh_val, r) * (r - hfem_val.node_r_ele[e, 1]) / hfem_val.length[e],
                                                       hfem_val.node_r_ele[e, 1],
                                                       hfem_val.node_r_ele[e, 2],
                                                       vh_val)
                    _ => return 0.0
                end

            2 =>
                @match q begin
                    2 => return 0.5 * le * (ed * ed + ed + 1.0 / 3.0) - le * le * (2.0 * ed / 3.0 + 0.5) +
                                GaussLegendre.gl_integ(r -> r * (r - hfem_val.node_r_ele[e, 1]) / hfem_val.length[e] * rVh(hfem_param, hfem_val, vh_val, r) * (r - hfem_val.node_r_ele[e, 1]) / hfem_val.length[e],
                                                            hfem_val.node_r_ele[e, 1],
                                                            hfem_val.node_r_ele[e, 2],
                                                            vh_val)
                    _ => return 0.0
                end

            _ => return 0.0
        end
    end

    get_B_matrix_element(e, le, p, q) = let
        ed = float(e - 1)
        @match p begin
            1 =>
                @match q begin
                    1 => return le * le * le * (ed * ed / 3.0 + ed / 6.0 + 1.0 / 30.0)
                    2 => return le * le * le * (ed * ed / 6.0 + ed / 6.0 + 0.05)
                    _ => return 0.0
                end

            2 =>
                @match q begin
                    2 => return le * le * le * (ed * ed / 3.0 + ed / 2.0 + 0.2)
                    _ => return 0.0
                end

            _ => return 0.0
        end
    end

    get_totalenergy(eigenenergy, hfem_param, hfem_val, vh_val) = let
        sum = 0.0
        max = length(hfem_val.phi) - 2

        # Simpsonの公式によって数値積分する
        @inbounds @simd for i = 1:2:max
            f0 = rVh(hfem_param, hfem_val, vh_val, hfem_val.node_r_glo[i]) * hfem_val.phi[i] * hfem_val.phi[i] * hfem_val.node_r_glo[i]
            f1 = rVh(hfem_param, hfem_val, vh_val, hfem_val.node_r_glo[i + 1]) * hfem_val.phi[i + 1] * hfem_val.phi[i + 1] * hfem_val.node_r_glo[i + 1]
            f2 = rVh(hfem_param, hfem_val, vh_val, hfem_val.node_r_glo[i + 2]) * hfem_val.phi[i + 2] * hfem_val.phi[i + 2] * hfem_val.node_r_glo[i + 2]

            sum += (f0 + 4.0 * f1 + f2)
        end

        return 2.0 * eigenenergy - sum * hfem_val.length[1] / 3.0
    end

    function make_data!(hfem_param, hfem_val)
        # Global節点のx座標を定義(R_MIN～R_MAX）
        dr = (hfem_param.R_MAX - hfem_param.R_MIN) / float(hfem_param.ELE_TOTAL)
        @inbounds for i = 0:hfem_param.NODE_TOTAL - 1
            # 計算領域を等分割
            hfem_val.node_r_glo[i + 1] = hfem_param.R_MIN + float(i) * dr
        end

        @inbounds for e = 1:hfem_param.ELE_TOTAL
            hfem_val.node_num_seg[e, 1] = e
            hfem_val.node_num_seg[e, 2] = e + 1
        end

        @inbounds for e = 1:hfem_param.ELE_TOTAL
            for i = 1:2
                hfem_val.node_r_ele[e, i] = hfem_val.node_r_glo[hfem_val.node_num_seg[e, i]]
            end
        end

        # 各線分要素の長さを計算
        @inbounds for e = 1:hfem_param.ELE_TOTAL
            hfem_val.length[e] = abs(hfem_val.node_r_ele[e, 2] - hfem_val.node_r_ele[e, 1])
        end
    end

    function make_element_matrix!(hfem_param, hfem_val, vh_val)
        # 要素行列の各成分を計算
        @inbounds for e = 1:hfem_param.ELE_TOTAL
            le = hfem_val.length[e]
            for j = 1:2
                for i = 1:j
                    hfem_val.mat_A_ele[e, i, j] = get_A_matrix_element(e, le, i, j, hfem_param, hfem_val, vh_val)
                    hfem_val.mat_B_ele[e, i, j] = get_B_matrix_element(e, le, i, j)
                end
            end
        end
    end

    function make_global_matrix(hfem_param, hfem_val)
        hg_tmp = Symmetric(zeros(hfem_param.NODE_TOTAL, hfem_param.NODE_TOTAL))
        ug_tmp = Symmetric(zeros(hfem_param.NODE_TOTAL, hfem_param.NODE_TOTAL))

        @inbounds for e = 1:hfem_param.ELE_TOTAL
            for j = 1:2
                for i = 1:j
                    hg_tmp.data[hfem_val.node_num_seg[e, i], hfem_val.node_num_seg[e, j]] += hfem_val.mat_A_ele[e, i, j]
                    ug_tmp.data[hfem_val.node_num_seg[e, i], hfem_val.node_num_seg[e, j]] += hfem_val.mat_B_ele[e, i, j]
                end
            end
        end

        return hg_tmp, ug_tmp
    end

    function normalize!(hfem_val)
        sum = 0.0
        len = length(hfem_val.phi)
        max = len - 2

        # Simpsonの公式によって数値積分する
        @inbounds @simd for i = 1:2:max
            f0 = hfem_val.phi[i] * hfem_val.phi[i] * hfem_val.node_r_glo[i] * hfem_val.node_r_glo[i]
            f1 = hfem_val.phi[i + 1] * hfem_val.phi[i + 1] * hfem_val.node_r_glo[i + 1] * hfem_val.node_r_glo[i + 1]
            f2 = hfem_val.phi[i + 2] * hfem_val.phi[i + 2] * hfem_val.node_r_glo[i + 2] * hfem_val.node_r_glo[i + 2]

            sum += (f0 + 4.0 * f1 + f2)
        end

        hfem_val.phi = abs.(hfem_val.phi * sqrt(sum * hfem_val.length[1] / 3.0))
    end
end


text:
kernels=filter(x->contains(x,".cu")&~contains(x,"~"),readdir("./"))

println("Compiling Kernels:")
for k in kernels
    println(k)
#    run(`nvcc -arch sm_30 -use_fast_math -ptx $k`) # Jetson TK1
    run(`/usr/local/cuda/bin/nvcc -c -use_fast_math -ptx $k`) # Titan GTX
end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
using Documenter, BinaryBuilder

makedocs(
    modules = [BinaryBuilder],
    format = :html,
    sitename = "BinaryBuilder.jl",
    pages = [
        "Home" => "index.md",
        "Building Packages" => "build_tips.md",
        "FAQ" => "FAQ.md",
        "Internals" => [
            "RootFS" => "rootfs.md",
            "Environment Variables" => "environment_variables.md",
            "Tricksy Gotchas" => "tricksy_gotchas.md",
            "Reference" => "reference.md",
        ],
    ],
)

deploydocs(
    deps = nothing,
    make = nothing,
    target = "build",
    repo = "github.com/JuliaPackaging/BinaryBuilder.jl.git",
    julia = "0.6",
)


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie
CairoMakie.activate!() #HIDE
let
    f(x,y) = (x + 2y^2) * abs(sin(y) + cos(x))
    x = y = 1:0.2:20
    z = [f(x,y) for x in x, y in y]
    fig = Figure(resolution=(900,400), fontsize = 16)
    ax1 = Axis(fig, aspect = 1, xlabel = "x", ylabel = "y", tickalign = 1)
    ax2 = Axis(fig, aspect = 1,  xlabel = "x", tickalign = 1)
    ax3 = Axis(fig, aspect = 1,  xlabel = "x", tickalign = 1)
    p1 = heatmap!(ax1, x, y, z, colormap = :plasma)
    contour!(ax2, x, y, z, color = :black, levels = 100:1:101)
    heatmap!(ax3, x, y, z, colormap = (:plasma, 0.5))
    contour!(ax3, x, y, z, color = :white, levels = 100:1:101)
    cbar = Colorbar(fig, p1, width = 10, ticksize=10, ticklabelsize = 10,
        tickalign = 1, height = Relative(0.7))
    limits!(ax1, 1,20,1,20)
    limits!(ax2, 1,20,1,20)
    limits!(ax3, 1,20,1,20)
    hideydecorations!(ax2, grid = false, ticks = false)
    hideydecorations!(ax3, grid = false, ticks = false)
    fig[1,1] = ax1
    fig[1,2] = ax2
    fig[1,3] = ax3
    fig[1,4] = cbar
    fig
    save(joinpath(@__DIR__, "output", "ContourOverHeatmap.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status("CairoMakie") # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using VSL
import Random: rand, rand!
import Test: @test, @test_throws
import Printf: @printf

brng_types = [
    VSL_BRNG_MCG31, VSL_BRNG_R250, VSL_BRNG_MRG32K3A, VSL_BRNG_MCG59, VSL_BRNG_WH, VSL_BRNG_SOBOL, VSL_BRNG_NIEDERR,
    VSL_BRNG_MT19937, VSL_BRNG_MT2203, VSL_BRNG_SFMT19937, VSL_BRNG_NONDETERM, VSL_BRNG_ARS5, VSL_BRNG_PHILOX4X32X10
]

printresults(values::AbstractArray) = printresults.(values)
printresults(value::Number) = println(bitstring(value))

macro test_rand(expr)
    values = @eval $expr
    printresults(values)
end

test_folder = joinpath(dirname(pathof(VSL)), "..", "test")


text:
using Test
using J4Org

function readcode(filename::AbstractString)
    if VERSION < v"0.7"
        code = readstring(filename)
    else
        code = read(filename,String)
    end
    code
end

# write your own tests here
@testset "J4Org" begin

    # initialize a default boxing module
    #
    initialize_boxing_module()

    # include test files
    #
    include("tokenizer.jl")
    include("extract.jl")
    include("find_tag.jl")
    include("links.jl")
    include("evaluate.jl")
    include("documented_item.jl")
    include("real_code.jl")

    import J4Org:
    tokenized,
    find_tag,
    org_string_comment,
    org_string_code,
    identifier,
    create_documented_item,
    org_string_documented_item

    code_1 = String("
    # This is a test file

    #+MyTag
    #
    # My documentedItem
    #
    # line 2
    #
    function myFunction (;x::Float64=sin(1))
        return 0
    end
        ")

    @testset "code_1" begin
        @test !isempty(code_1)
        t=tokenized(code_1)
        r=find_tag(t,1,"MyTag")
        di=create_documented_item(t,tag_idx(r))
        @test org_string_comment(di,[di],[di],"","BoxingModule") == "#+BEGIN_QUOTE\nMy documentedItem\n\nline 2\n#+END_QUOTE\n"
        @test org_string_code(di) == "#+BEGIN_SRC julia :eval never :exports code\nfunction myFunction (;x::Float64=sin(1))\n#+END_SRC\n"
        @test identifier(di) == "myFunction"
        @test org_string_documented_item(di,[di],[di]) ==  "- @@latex:\\phantomsection@@ *=myFunction=* \n#+BEGIN_SRC julia :eval never :exports code\nfunction myFunction (;x::Float64=sin(1))\n#+END_SRC\n#+BEGIN_QUOTE\nMy documentedItem\n\nline 2\n#+END_QUOTE\n\n"
    end

end;
nothing


text:
using OffsetArrays

Flux = OffsetVector
zero_flux(Nx :: Int) = Flux(zeros(Nx+1), 0 : Nx)
zero_flux(U :: Vector) = Flux(zeros(length(U)+1), 0 : length(U))

struct ScalarQuantity
    data :: Vector
    posFlux :: Flux
    negFlux :: Flux

    ScalarQuantity(U :: Vector) = new(copy(U), zero_flux(U), zero_flux(U))
end

struct RelaxPb
    U :: ScalarQuantity
    V :: ScalarQuantity

    fluxP :: Flux
    fluxV :: Flux
    dxP :: Vector
    dxV :: Vector

    f :: Function
    p :: Function

    function RelaxPb(U :: Vector, V :: Vector,
                     f :: Function, p :: Function)

        Nx = length(U)
        pbU, pbV = ScalarQuantity(U), ScalarQuantity(V)
        fluxV, fluxP = zero_flux(Nx), zero_flux(Nx)
        dxV, dxP = zeros(Nx), zeros(Nx)

        new(pbU,pbV,fluxP,fluxV,dxV,dxP,f,p)
    end
end

function Base.copy(pb :: RelaxPb)
    RelaxPb(pb.U.data,pb.V.data,pb.f,pb.p)
end

max_stars_repo_name: tremelow/StiffKinetic

max_issues_repo_name: tremelow/StiffKinetic

max_forks_repo_name: tremelow/StiffKinetic


text:
# Adapt.adapt(T, x::Array) = T.(x)
# Adapt.adapt(T, x::Real) = T(x)

"""
  to32(m)

  convert all numbers, arrays, and parameters in model to Float32

"""
to32(m) = mapleaves(x -> Adapt.adapt(Float32, x), m)

logit_cross_entropy(logit, y) = -sum(y.*(logit.-logsumexp(logit, 1))) / size(logit,2)

weighted_logit_cross_entropy(logit, y, w::AbstractVector) = weighted_logit_cross_entropy(logit, y, w')
weighted_logit_cross_entropy(logit, y::AbstractVector, w::AbstractMatrix) = weighted_logit_cross_entropy(logit, Flux.onehotbatch(y, 1:size(logit,1)), w)
weighted_logit_cross_entropy(logit, y::AbstractMatrix, w::AbstractMatrix) = -sum(y.*w.*(logit.-logsumexp(logit, 1)))

entropy_dim(x,dim::Int=1) = -mean(sum(x .* log.(x),dim))

using StatsBase

function classweightvector(::Type{S},y::Vector{T}) where {S<:AbstractFloat,T<:Integer}
  w=ones(S,size(y,1))
  classsizes = StatsBase.countmap(y)
  for j in 1:size(y,1)
    w[j] = S(1/size(y,2))/classsizes[y[j]]
  end
  return(w./sum(w))
end

function classweightvector(y::AbstractArray{T},classweights::Vector{S}) where {T<:Integer,S<:AbstractFloat}
  classweights = classweights./sum(classweights)
  w=ones(S,size(y,1))
  classsizes = StatsBase.countmap(y)
  for j in 1:size(y,1)
    w[j] = classweights[y[j]]/classsizes[y[j]]
  end
  if sum(isnan.(w))>0 || sum(isinf.(w))>0
    save("error.jld2","y",y,"classweights",classweights)
    error("nans or infs in classweights")
  end
  w
end

function confusion(y,pred,k)
	c = zeros(Int,k,k)
	foreach(i -> c[i[1],i[2]] += 1, zip(y,pred))
	c
end

max_stars_repo_name: pevnak/FluxExtensions

max_issues_repo_name: pevnak/FluxExtensions

max_forks_repo_name: pevnak/FluxExtensions


text:
function gensys_uncertain_zlb(prob_vec::AbstractVector{S}, Th::AbstractMatrix{S}, Ch::AbstractVector{S},
                              Tzlbs::Vector{Matrix{S}}, Rzlbs::Vector{Matrix{S}}, Czlbs::Vector{Vector{S}},
                              Γ0_til::AbstractMatrix{S}, Γ1_til::AbstractMatrix{S}, Γ2_til::AbstractMatrix{S},
                              C_til::AbstractVector{S}, Ψ_til::AbstractMatrix{S}) where {S <: Real}

    Tbars = Vector{Matrix{S}}(undef, length(Tzlbs) + 1)
    Cbars = Vector{Vector{S}}(undef, length(Tzlbs) + 1)
    Tout  = Vector{Matrix{S}}(undef, length(Tzlbs) + 1)
    Rout  = Vector{Matrix{S}}(undef, length(Tzlbs) + 1)
    Cout  = Vector{Vector{S}}(undef, length(Tzlbs) + 1)

    # Calculate "uncertain" ZLB matrices and back out the transition equation
    for i in 1:length(Tzlbs)
        Tbars[i] = prob_vec[1] * Tzlbs[i] + prob_vec[2] * Th # it is assumed Tzlbs is a vector of the T_{t + 1}^{(zlb)} matrices
        Cbars[i] = prob_vec[1] * Czlbs[i] + prob_vec[2] * Ch

        Tout[i]  = (Γ2_til * Tbars[i] + Γ0_til) \ Γ1_til
        Rout[i]  = (Γ2_til * Tbars[i] + Γ0_til) \ Ψ_til
        Cout[i]  = (Γ2_til * Tbars[i] + Γ0_til) \ (C_til  - Γ2_til * Cbars[i])
    end

    # Add boundary condition
    Tout[end] = Tzlbs[end]
    Rout[end] = Rzlbs[end]
    Cout[end] = Czlbs[end]

    return Tout, Rout, Cout
end


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

export my_function  # remember to export the functions you want to use!

# Let's define a function that will be used in the "tutorial_basic_sintax.txt" file
# NOTE: DO NOT MODIFY THIS FUNCTION: IT?S USED FOR TESTS!
function my_function(x::Float64)
     return x+1.0
end

#    PAY ATTENTION: ALL THE NUMBERS DEFINED IN THE SCENEFILE ARE PARSED *ALWAYS*
#    ARE PARSED *ALWAYS* AS FLOAT64.
#
#    => DO NOT DEFINE FUNCTIONS THAT TAKE AS ARGUMENTS INTEGER VALUES,
#         THEY WILL NEVER BE USED!!!

function func(x::T) where {T<:Real}
     return (x, 0.5)
end

function earth_moon_sun(x, P_rev_moon=200, P_rot_earth=50, d=3.0, v=1.0)
     moon_x = d * cos(2 * π * x/P_rev_moon)
     moon_y =  d * sin(2 * π * x/P_rev_moon)
     moon_z =  v * cos(2 * π * x/P_rev_moon) * 0.0
     moon_rotang = 2 * π * x/P_rev_moon
     earth_rotang = 2 * π * x/P_rot_earth
     return (moon_x, moon_y, moon_z, moon_rotang, earth_rotang)
end


text:
using IterativeSolvers

fnames = (:cg,:minres,:gmres,:idrs,:bicgstabl,:chebyshev,:jacobi,:gauss_seidel,:sor,:ssor,:lsmr,:lsqr)
for fname ∈ fnames
    @eval begin
        function IterativeSolvers.$(fname)(mls::MaxwellLS, args...; kwargs...)
            $(Symbol(fname,"!"))(view(mls.solution.tot.values,:,1), mls.maxwell.A, mls.j, args...; kwargs...)
            return mls.solution
        end
        function IterativeSolvers.$(Symbol(fname,"!"))(mls::MaxwellLS, args...; kwargs...)
            _, hist = $(Symbol(fname,"!"))(view(mls.solution.tot.values,:,1), mls.maxwell.A, reshape(Array(mls.j),length(mls.j),:), args...; kwargs..., log=true)
            for i ∈ eachindex(mls.equivalent_source.field)
                k = mod1(i,length(mls.j)÷3)
                mls.solution.incident[i] = mls.equivalent_source.field[i]*mls.equivalent_source.mask(mls.equivalent_source.field.pos[k])
                mls.solution.scattered[i] = mls.solution.total[i] - mls.solution.incident[i]
            end
            mls.solved[] = true
            mls.converged[] = hist.isconverged
            return nothing
        end
    end
end