text:
# kernels for reduction

### facilities for debugging

function safe_sum(x)
    r = zero(eltype(x))
    for i in 1 : length(x)
        r = r + x[i]
    end
    r
end

function safe_max(x)
    r = typemin(eltype(x))
    for i in 1 : length(x)
        r = max(r, x[i])
    end
    r
end

function safe_min(x)
    r = typemax(eltype(x))
    for i in 1 : length(x)
        r = min(r, x[i])
    end
    r
end

# rename @inbounds to @inbounds (when want to disable @inbounds effect)
macro inboundss(x)
    esc(x)
end

#################################################
#
#   Core skeleton
#
#################################################

function saccum{Op<:Functor{2}}(op::Op, n::Int, a::Array, ia::Int)
    ia_end = ia + n
    @inbounds if n <= 4
        v = getvalue(a, ia)
        ia += 1
        while ia < ia_end
            v = evaluate(op, v, getvalue(a, ia))
            ia += 1
        end
    else
        v1 = getvalue(a, ia)
        v2 = getvalue(a, ia+1)
        v3 = getvalue(a, ia+2)
        v4 = getvalue(a, ia+3)
        ia += 4
        ia_ = ia_end - 3
        while ia < ia_
            v1 = evaluate(op, v1, getvalue(a, ia))
            v2 = evaluate(op, v2, getvalue(a, ia+1))
            v3 = evaluate(op, v3, getvalue(a, ia+2))
            v4 = evaluate(op, v4, getvalue(a, ia+3))
            ia += 4
        end
        v = evaluate(op, evaluate(op, v1, v2), evaluate(op, v3, v4))
        while ia < ia_end
            v = evaluate(op, v, getvalue(a, ia))
            ia += 1
        end
    end
    return v
end

function saccum{Op<:Functor{2}}(op::Op, n::Int, fun::Functor{1}, a::Array, ia::Int)
    ia_end = ia + n
    @inbounds if n <= 4
        v = evaluate(fun, getvalue(a, ia))
        ia += 1
        while ia < ia_end
            v = evaluate(op, v, evaluate(fun, getvalue(a, ia)))
            ia += 1
        end
    else
        v1 = evaluate(fun, getvalue(a, ia))
        v2 = evaluate(fun, getvalue(a, ia+1))
        v3 = evaluate(fun, getvalue(a, ia+2))
        v4 = evaluate(fun, getvalue(a, ia+3))
        ia += 4
        ia_ = ia_end - 3
        while ia < ia_
            v1 = evaluate(op, v1, evaluate(fun, getvalue(a, ia)))
            v2 = evaluate(op, v2, evaluate(fun, getvalue(a, ia+1)))
            v3 = evaluate(op, v3, evaluate(fun, getvalue(a, ia+2)))
            v4 = evaluate(op, v4, evaluate(fun, getvalue(a, ia+3)))
            ia += 4
        end
        v = evaluate(op, evaluate(op, v1, v2), evaluate(op, v3, v4))
        while ia < ia_end
            v = evaluate(op, v, evaluate(fun, getvalue(a, ia)))
            ia += 1
        end
    end
    return v
end

function saccum{Op<:Functor{2}}(op::Op, n::Int, fun::Functor{2}, a::ArrayOrNum, ia::Int, b::ArrayOrNum, ib::Int)
    ia_end = ia + n
    @inbounds if n <= 4
        v = evaluate(fun, getvalue(a, ia), getvalue(b, ib))
        ia += 1
        ib += 1
        while ia < ia_end
            v = evaluate(op, v, evaluate(fun, getvalue(a, ia), getvalue(b, ib)))
            ia += 1
            ib += 1
        end
    else
        v1 = evaluate(fun, getvalue(a, ia), getvalue(b, ib))
        v2 = evaluate(fun, getvalue(a, ia+1), getvalue(b, ib+1))
        ia += 2
        ib += 2
        ia_ = ia_end - 1
        while ia < ia_
            v1 = evaluate(op, v1, evaluate(fun, getvalue(a, ia), getvalue(b, ib)))
            v2 = evaluate(op, v2, evaluate(fun, getvalue(a, ia+1), getvalue(b, ib+1)))
            ia += 2
            ib += 2
        end
        v = evaluate(op, v1, v2)
        if ia < ia_end
            v = evaluate(op, v, evaluate(fun, getvalue(a, ia), getvalue(b, ib)))
        end
    end
    return v
end

function saccum_fdiff{Op<:Functor{2}}(op::Op, n::Int, fun::Functor{1}, a::ArrayOrNum, ia::Int, b::ArrayOrNum, ib::Int)
    ia_end = ia + n
    @inbounds if n <= 4
        v = evaluate(fun, getvalue(a, ia) - getvalue(b, ib))
        ia += 1
        ib += 1
        while ia < ia_end
            v = evaluate(op, v, evaluate(fun, getvalue(a, ia) - getvalue(b, ib)))
            ia += 1
            ib += 1
        end
    else
        v1 = evaluate(fun, getvalue(a, ia) - getvalue(b, ib))
        v2 = evaluate(fun, getvalue(a, ia+1) - getvalue(b, ib+1))
        ia += 2
        ib += 2
        ia_ = ia_end - 1
        while ia < ia_
            v1 = evaluate(op, v1, evaluate(fun, getvalue(a, ia) - getvalue(b, ib)))
            v2 = evaluate(op, v2, evaluate(fun, getvalue(a, ia+1) - getvalue(b, ib+1)))
            ia += 2
            ib += 2
        end
        v = evaluate(op, v1, v2)
        if ia < ia_end
            v = evaluate(op, v, evaluate(fun, getvalue(a, ia) - getvalue(b, ib)))
        end
    end
    return v
end

function saccum{Op<:Functor{2}}(op::Op, n::Int, fun::Functor{3}, a::ArrayOrNum, ia::Int, b::ArrayOrNum, ib::Int, c::ArrayOrNum, ic::Int)
    ia_end = ia + n
    @inbounds if n <= 4
        v = evaluate(fun, getvalue(a, ia), getvalue(b, ib), getvalue(c, ic))
        ia += 1
        ib += 1
        ic += 1
        while ia < ia_end
            v = evaluate(op, v, evaluate(fun, getvalue(a, ia), getvalue(b, ib), getvalue(c, ic)))
            ia += 1
            ib += 1
            ic += 1
        end
    else
        v1 = evaluate(fun, getvalue(a, ia), getvalue(b, ib), getvalue(c, ic))
        v2 = evaluate(fun, getvalue(a, ia+1), getvalue(b, ib+1), getvalue(c, ic+1))
        ia += 2
        ib += 2
        ic += 2
        ia_ = ia_end - 1
        while ia < ia_
            v1 = evaluate(op, v1, evaluate(fun, getvalue(a, ia), getvalue(b, ib), getvalue(c, ic)))
            v2 = evaluate(op, v2, evaluate(fun, getvalue(a, ia+1), getvalue(b, ib+1), getvalue(c, ic+1)))
            ia += 2
            ib += 2
            ic += 2
        end
        v = evaluate(op, v1, v2)
        if ia < ia_end
            v = evaluate(op, v, evaluate(fun, getvalue(a, ia), getvalue(b, ib), getvalue(c, ic)))
        end
    end
    return v
end

function saccum{Op<:Functor{2}}(op::Op, n::Int, a::Array, ia::Int, sa::Int)
    ia_end = ia + sa * n
    @inbounds if n <= 4
        v = getvalue(a, ia)
        ia += sa
        while ia < ia_end
            v = evaluate(op, v, getvalue(a, ia))
            ia += sa
        end
    else
        sa2 = 2sa
        v1 = getvalue(a, ia)
        v2 = getvalue(a, ia+sa)
        ia += sa2
        ia_ = ia_end - sa
        while ia < ia_
            v1 = evaluate(op, v1, getvalue(a, ia))
            v2 = evaluate(op, v2, getvalue(a, ia+sa))
            ia += sa2
        end
        v = evaluate(op, v1, v2)
        if ia < ia_end
            v = evaluate(op, v, getvalue(a, ia))
        end
    end
    return v
end

function saccum{Op<:Functor{2}}(op::Op, n::Int, fun::Functor{1}, a::Array, ia::Int, sa::Int)
    ia_end = ia + sa * n
    @inbounds if n <= 4
        v = evaluate(fun, getvalue(a, ia))
        ia += sa
        while ia < ia_end
            v = evaluate(op, v, evaluate(fun, getvalue(a, ia)))
            ia += sa
        end
    else
        sa2 = 2sa
        v1 = evaluate(fun, getvalue(a, ia))
        v2 = evaluate(fun, getvalue(a, ia+sa))
        ia += sa2
        ia_ = ia_end - sa
        while ia < ia_
            v1 = evaluate(op, v1, evaluate(fun, getvalue(a, ia)))
            v2 = evaluate(op, v2, evaluate(fun, getvalue(a, ia+sa)))
            ia += sa2
        end
        v = evaluate(op, v1, v2)
        if ia < ia_end
            v = evaluate(op, v, evaluate(fun, getvalue(a, ia)))
        end
    end
    return v
end

function saccum{Op<:Functor{2}}(op::Op, n::Int, fun::Functor{2}, a::ArrayOrNum, ia::Int, sa::Int, b::ArrayOrNum, ib::Int, sb::Int)
    ia_end = ia + sa * n
    @inbounds if n <= 4
        v = evaluate(fun, getvalue(a, ia), getvalue(b, ib))
        ia += sa
        ib += sb
        while ia < ia_end
            v = evaluate(op, v, evaluate(fun, getvalue(a, ia), getvalue(b, ib)))
            ia += sa
            ib += sb
        end
    else
        sa2 = 2sa
        sb2 = 2sb
        v1 = evaluate(fun, getvalue(a, ia), getvalue(b, ib))
        v2 = evaluate(fun, getvalue(a, ia+sa), getvalue(b, ib+sb))
        ia += sa2
        ib += sb2
        ia_ = ia_end - sa
        while ia < ia_
            v1 = evaluate(op, v1, evaluate(fun, getvalue(a, ia), getvalue(b, ib)))
            v2 = evaluate(op, v2, evaluate(fun, getvalue(a, ia+sa), getvalue(b, ib+sb)))
            ia += sa2
            ib += sb2
        end
        v = evaluate(op, v1, v2)
        if ia < ia_end
            v = evaluate(op, v, evaluate(fun, getvalue(a, ia), getvalue(b, ib)))
        end
    end
    return v
end

function saccum_fdiff{Op<:Functor{2}}(op::Op, n::Int, fun::Functor{1}, a::ArrayOrNum, ia::Int, sa::Int, b::ArrayOrNum, ib::Int, sb::Int)
    ia_end = ia + sa * n
    @inbounds if n <= 4
        v = evaluate(fun, getvalue(a, ia) - getvalue(b, ib))
        ia += sa
        ib += sb
        while ia < ia_end
            v = evaluate(op, v, evaluate(fun, getvalue(a, ia) - getvalue(b, ib)))
            ia += sa
            ib += sb
        end
    else
        sa2 = 2sa
        sb2 = 2sb
        v1 = evaluate(fun, getvalue(a, ia) - getvalue(b, ib))
        v2 = evaluate(fun, getvalue(a, ia+sa) - getvalue(b, ib+sb))
        ia += sa2
        ib += sb2
        ia_ = ia_end - sa
        while ia < ia_
            v1 = evaluate(op, v1, evaluate(fun, getvalue(a, ia) - getvalue(b, ib)))
            v2 = evaluate(op, v2, evaluate(fun, getvalue(a, ia+sa) - getvalue(b, ib+sb)))
            ia += sa2
            ib += sb2
        end
        v = evaluate(op, v1, v2)
        if ia < ia_end
            v = evaluate(op, v, evaluate(fun, getvalue(a, ia) - getvalue(b, ib)))
        end
    end
    return v
end

function saccum{Op<:Functor{2}}(op::Op, n::Int, fun::Functor{3}, a::ArrayOrNum, ia::Int, sa::Int,
                                                                b::ArrayOrNum, ib::Int, sb::Int,
                                                                c::ArrayOrNum, ic::Int, sc::Int)
    ia_end = ia + sa * n
    @inbounds if n <= 4
        v = evaluate(fun, getvalue(a, ia), getvalue(b, ib), getvalue(c, ic))
        ia += sa
        ib += sb
        ic += sc
        while ia < ia_end
            v = evaluate(op, v, evaluate(fun, getvalue(a, ia), getvalue(b, ib), getvalue(c, ic)))
            ia += sa
            ib += sb
            ic += sc
        end
    else
        sa2 = 2sa
        sb2 = 2sb
        sc2 = 2sc
        v1 = evaluate(fun, getvalue(a, ia), getvalue(b, ib), getvalue(c, ic))
        v2 = evaluate(fun, getvalue(a, ia+sa), getvalue(b, ib+sb), getvalue(c, ic+sc))
        ia += sa2
        ib += sb2
        ic += sc2
        ia_ = ia_end - sa
        while ia < ia_
            v1 = evaluate(op, v1, evaluate(fun, getvalue(a, ia), getvalue(b, ib), getvalue(c, ic)))
            v2 = evaluate(op, v2, evaluate(fun, getvalue(a, ia+sa), getvalue(b, ib+sb), getvalue(c, ic+sc)))
            ia += sa2
            ib += sb2
            ic += sc2
        end
        v = evaluate(op, v1, v2)
        if ia < ia_end
            v = evaluate(op, v, evaluate(fun, getvalue(a, ia), getvalue(b, ib), getvalue(c, ic)))
        end
    end
    return v
end

function paccum!{Op<:Functor{2}}(op::Op, n::Int, d::Array, id::Int, a::Array, ia::Int)
    ia_end = ia + n
    ia_ = ia_end - 3
    @inbounds while ia < ia_
        d[id] = evaluate(op, d[id], getvalue(a, ia))
        d[id+1] = evaluate(op, d[id+1], getvalue(a, ia+1))
        d[id+2] = evaluate(op, d[id+2], getvalue(a, ia+2))
        d[id+3] = evaluate(op, d[id+3], getvalue(a, ia+3))
        ia += 4
        id += 4
    end

    @inbounds while ia < ia_end
        d[id] = evaluate(op, d[id], getvalue(a, ia))
        id += 1
        ia += 1
    end
end

function paccum!{Op<:Functor{2}}(op::Op, n::Int, d::Array, id::Int, fun::Functor{1}, a::Array, ia::Int)
    ia_end = ia + n
    ia_ = ia_end - 3
    @inbounds while ia < ia_
        d[id] = evaluate(op, d[id], evaluate(fun, getvalue(a, ia)))
        d[id+1] = evaluate(op, d[id+1], evaluate(fun, getvalue(a, ia+1)))
        d[id+2] = evaluate(op, d[id+2], evaluate(fun, getvalue(a, ia+2)))
        d[id+3] = evaluate(op, d[id+3], evaluate(fun, getvalue(a, ia+3)))
        ia += 4
        id += 4
    end

    @inbounds while ia < ia_end
        d[id] = evaluate(op, d[id], evaluate(fun, getvalue(a, ia)))
        id += 1
        ia += 1
    end
end

function paccum!{Op<:Functor{2}}(op::Op, n::Int, d::Array, id::Int, fun::Functor{2}, a::ArrayOrNum, ia::Int, b::ArrayOrNum, ib::Int)
    ia_end = ia + n
    ia_ = ia_end - 1
    @inbounds while ia < ia_
        d[id] = evaluate(op, d[id], evaluate(fun, getvalue(a, ia), getvalue(b, ib)))
        d[id+1] = evaluate(op, d[id+1], evaluate(fun, getvalue(a, ia+1), getvalue(b, ib+1)))
        ia += 2
        ib += 2
        id += 2
    end

    @inbounds if ia < ia_end
        d[id] = evaluate(op, d[id], evaluate(fun, getvalue(a, ia), getvalue(b, ib)))
    end
end

function paccum_fdiff!{Op<:Functor{2}}(op::Op, n::Int, d::Array, id::Int, fun::Functor{1}, a::ArrayOrNum, ia::Int, b::ArrayOrNum, ib::Int)
    ia_end = ia + n
    ia_ = ia_end - 1
    @inbounds while ia < ia_
        d[id] = evaluate(op, d[id], evaluate(fun, getvalue(a, ia) - getvalue(b, ib)))
        d[id+1] = evaluate(op, d[id+1], evaluate(fun, getvalue(a, ia+1) - getvalue(b, ib+1)))
        ia += 2
        ib += 2
        id += 2
    end

    @inbounds if ia < ia_end
        d[id] = evaluate(op, d[id], evaluate(fun, getvalue(a, ia) - getvalue(b, ib)))
    end
end

function paccum!{Op<:Functor{2}}(op::Op, n::Int, d::Array, id::Int, fun::Functor{3}, a::ArrayOrNum, ia::Int,
                                                                                b::ArrayOrNum, ib::Int,
                                                                                c::ArrayOrNum, ic::Int)
    ia_end = ia + n
    ia_ = ia_end - 1
    @inbounds while ia < ia_
        d[id] = evaluate(op, d[id], evaluate(fun, getvalue(a, ia), getvalue(b, ib), getvalue(c, ic)))
        d[id+1] = evaluate(op, d[id+1], evaluate(fun, getvalue(a, ia+1), getvalue(b, ib+1), getvalue(c, ic+1)))
        ia += 2
        ib += 2
        ic += 2
        id += 2
    end

    @inbounds if ia < ia_end
        d[id] = evaluate(op, d[id], evaluate(fun, getvalue(a, ia), getvalue(b, ib), getvalue(c, ic)))
    end
end

function paccum!{Op<:Functor{2}}(op::Op, n::Int, d::Array, id::Int, sd::Int, a::Array, ia::Int, sa::Int)
    ia_end = ia + sa * n
    sa2 = 2sa
    sd2 = 2sd
    ia_ = ia_end - sa
    @inbounds while ia < ia_
        d[id] = evaluate(op, d[id], getvalue(a, ia))
        d[id + sd] = evaluate(op, d[id+sd], getvalue(a, ia+sa))
        ia += sa2
        id += sd2
    end

    @inbounds if ia < ia_end
        d[id] = evaluate(op, d[id], getvalue(a, ia))
    end
end

function paccum!{Op<:Functor{2}}(op::Op, n::Int, d::Array, id::Int, sd::Int, fun::Functor{1}, a::Array, ia::Int, sa::Int)
    ia_end = ia + sa * n
    sa2 = 2sa
    sd2 = 2sd
    ia_ = ia_end - sa
    @inbounds while ia < ia_
        d[id] = evaluate(op, d[id], evaluate(fun, getvalue(a, ia)))
        d[id + sd] = evaluate(op, d[id+sd], evaluate(fun, getvalue(a, ia+sa)))
        ia += sa2
        id += sd2
    end

    @inbounds if ia < ia_end
        d[id] = evaluate(op, d[id], evaluate(fun, getvalue(a, ia)))
    end
end

function paccum!{Op<:Functor{2}}(op::Op, n::Int, d::Array, id::Int, sd::Int, fun::Functor{2}, a::ArrayOrNum, ia::Int, sa::Int,
                                                                                b::ArrayOrNum, ib::Int, sb::Int)
    ia_end = ia + sa * n
    sa2 = 2sa
    sb2 = 2sb
    sd2 = 2sd
    ia_ = ia_end - sa
    @inbounds while ia < ia_
        d[id] = evaluate(op, d[id], evaluate(fun, getvalue(a, ia), getvalue(b, ib)))
        d[id + sd] = evaluate(op, d[id+sd], evaluate(fun, getvalue(a, ia+sa), getvalue(b, ib+sb)))
        ia += sa2
        ib += sb2
        id += sd2
    end

    @inbounds if ia < ia_end
        d[id] = evaluate(op, d[id], evaluate(fun, getvalue(a, ia), getvalue(b, ib)))
    end
end

function paccum_fdiff!{Op<:Functor{2}}(op::Op, n::Int, d::Array, id::Int, sd::Int, fun::Functor{1}, a::ArrayOrNum, ia::Int, sa::Int,
                                                                                b::ArrayOrNum, ib::Int, sb::Int)
    ia_end = ia + sa * n
    sa2 = 2sa
    sb2 = 2sb
    sd2 = 2sd
    ia_ = ia_end - sa
    @inbounds while ia < ia_
        d[id] = evaluate(op, d[id], evaluate(fun, getvalue(a, ia) - getvalue(b, ib)))
        d[id + sd] = evaluate(op, d[id+sd], evaluate(fun, getvalue(a, ia+sa) - getvalue(b, ib+sb)))
        ia += sa2
        ib += sb2
        id += sd2
    end

    @inbounds if ia < ia_end
        d[id] = evaluate(op, d[id], evaluate(fun, getvalue(a, ia) - getvalue(b, ib)))
    end
end

function paccum!{Op<:Functor{2}}(op::Op, n::Int, d::Array, id::Int, sd::Int, fun::Functor{3}, a::ArrayOrNum, ia::Int, sa::Int,
                                                                                b::ArrayOrNum, ib::Int, sb::Int,
                                                                                c::ArrayOrNum, ic::Int, sc::Int)
    ia_end = ia + sa * n
    sa2 = 2sa
    sb2 = 2sb
    sc2 = 2sc
    sd2 = 2sd
    ia_ = ia_end - sa
    @inbounds while ia < ia_
        d[id] = evaluate(op, d[id], evaluate(fun, getvalue(a, ia), getvalue(b, ib), getvalue(c, ic)))
        d[id + sd] = evaluate(op, d[id+sd], evaluate(fun, getvalue(a, ia+sa), getvalue(b, ib+sb), getvalue(c, ic+sc)))
        ia += sa2
        ib += sb2
        ic += sc2
        id += sd2
    end

    @inbounds if ia < ia_end
        d[id] = evaluate(op, d[id], evaluate(fun, getvalue(a, ia), getvalue(b, ib), getvalue(c, ic)))
    end
end

#################################################
#
#   BLAS acceleration
#
#################################################

saccum{T<:BlasFP}(::Add, n::Int, f::AbsFun, a::Array{T}, ia::Int) =
    BLAS.asum(n, pointer(a, ia), 1)

saccum{T<:BlasFP}(::Add, n::Int, f::Abs2Fun, a::Array{T}, ia::Int) =
    (p = pointer(a, ia); BLAS.dot(n, p, 1, p, 1))

saccum{T<:BlasFP}(::Add, n::Int, f::Multiply, a::Array{T}, ia::Int, b::Array{T}, ib::Int) =
    (pa = pointer(a, ia); pb = pointer(b, ib); BLAS.dot(n, pa, 1, pb, 1))


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, StatsBase, Distributions, Plots, LaTeXStrings
	Random.seed!(0)
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing5.16"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	mu0, mu1a, mu1b, mu1c, sd = 15, 16, 18, 20, 2
	tauGrid = 5:0.1:25

	dist0 = Normal(mu0,sd)
	dist1a, dist1b, dist1c  = Normal(mu1a,sd), Normal(mu1b,sd), Normal(mu1c,sd)

	falsePositive = ccdf.(dist0,tauGrid)
	truePositiveA, truePositiveB, truePositiveC =
	    ccdf.(dist1a,tauGrid), ccdf.(dist1b,tauGrid), ccdf.(dist1c,tauGrid)
end

# ╔═╡ a80f7e58-3af9-11eb-04b3-57fcb9ee857c
begin
	plot(falsePositive, [truePositiveA truePositiveB truePositiveC],
		c=[:blue :red :green],
		label=[L"H1a: \mu_1 = 16" L"H1b: \mu_1 = 18" L"H1c: \mu_1 = 20"])
	plot!([0,1], [0,1], c=:black, ls=:dash, label="H0 = H1 = 15",
		xlims=(0,1), ylims=(0,1), xlabel=L"\alpha", ylabel="Power",
		ratio=:equal, legend=:bottomright)
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing5.16"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═a80f7e58-3af9-11eb-04b3-57fcb9ee857c
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12

ext: jl


text:
function ψp0_BSIM6(Vg,Vfb,Na,Cox,T,phif)
    G = γC(Na,Cox)/sqrt(kb*T)
    vg = Vg/kb/T
    vfb = Vfb/kb/T
    temp = (vg-vfb)*0.5-3*(1+G/sqrt(2))
    return temp + sqrt( temp^2 + 6*(vg-vfb) )
end

function ψp_BSIM6(Vg,Vfb,Na,Cox,T,phif)
    G = γC(Na,Cox)/sqrt(kb*T)
    vg = Vg/kb/T
    vfb = Vfb/kb/T
    δpD = 0.0
    ψ0 = ψp0_BSIM6(Vg,Vfb,Na,Cox,T,phif)
    if vg - vfb > 0
        return 1 -exp(-ψ0) + (sqrt(vg-vfb-1 + exp(-ψ0)+(G/2)^2)-G/2)^2
    else
        return -log(1-ψ0+ ( (vg-vfb-ψ0)/G )^2 )
    end
end

function nq_BSIM6(Vg,Vfb,Na,Cox,T,phif)
    G = γC(Na,Cox)/sqrt(kb*T)
    return 1 + G/2/sqrt(ψp_BSIM6(Vg,Vfb,Na,Cox,T,phif))
end

function q_BSIM6(Vg,Vch,Vfb,Na,Cox,T,phib)
    gam = γC(Na,Cox)/sqrt(kb*T)
    vch = Vch/kb/T
    phib = phib/kb/T/nq_BSIM6(Vg,Vfb,Na,Cox,T,phib)
    psip = ψp_BSIM6(Vg,Vfb,Na,Cox,T,phib)
    T8 = 0.5 * (psip + 1.0 + sqrt((psip - 1.0) * (psip - 1.0) + 0.25 * 2.0 * 2.0))
    sqrtpsip = sqrt(T8)
    T9 = 1.0 + gam / (2.0 * sqrtpsip)
    T0 = (1.0 + (gam / (2.0 * sqrtpsip))) / gam
    T1 = psip - 2.0 * phib - vch
    T2 = T1 - log(4.0 * T0 * sqrtpsip);
    T8 = 0.5 * (T2 - 0.201491 - sqrt(T2 * (T2 + 0.402982) + 2.446562));
    sqrtpsisa = sqrtpsip
    if (T8 <= -68.0)
        T4 = -100.0
        T5 = 20.0
        if (T8 < T4 - 0.5 * T5)
            T3 = exp(T4)
        else
            if (T8 > T4 + 0.5 * T5)
                T3 = exp(T8)
            else
                T2 = (T8 - T4) / T5
                T6 = T2 * T2
                T3 = exp(T4 + T5 * ((5.0 / 64.0) + 0.5 * T2 + T6 * ((15.0 / 16.0) - T6 * (1.25 - T6))))
            end
        end
        q = T3 * (1.0 + T1 - T8 - log(2.0 * T0 * (T3 * 2.0 * T0 + 2.0 * sqrtpsisa)))
    else
        T3 = exp(T8)
        sqrtpsisainv = 1.0 / sqrtpsisa
        T4 = 2.0 * T3 + log(T3 * 2.0 * T0 * (T3 *  2.0 * T0 + 2.0 * sqrtpsisa)) - T1
        T5 = 2.0 + (1.0 / T3) + (T0 + sqrtpsisainv) / (T0 * T3 + sqrtpsisa)
        T3 = T3 - T4 / T5
        T4 = 2.0 * T3 + log(T3 * 2.0 * T0 * (T3 * 2.0 * T0 + 2.0 * sqrtpsisa)) - T1
        T5 = 2.0 + (1.0 / T3) + (T0 + sqrtpsisainv) / (T0 * T3 + sqrtpsisa)
        T6 = ((T0 + sqrtpsisainv) / (T0 * T3 + sqrtpsisa)) * ((T0 + sqrtpsisainv) / (T0 * T3 + sqrtpsisa))
        T7 = -((1.0 / T3) * (1.0 / T3)) - (1.0 / (sqrtpsisa * sqrtpsisa * sqrtpsisa * (T0 * T3 + sqrtpsisa))) - T6
        q  = T3 - (T4 / T5) * (1.0 + T4 * T7 / (2.0 * T5 * T5))
    end
    return q
end

function id_BSIM6(Vg,Vd,Vs,Vfb,Na,Cox,T,phif)
    qs = q_BSIM6(Vg,Vs,Vfb,Na,Cox,T,phif)
    qd = q_BSIM6(Vg,Vd,Vfb,Na,Cox,T,phif)
    nq = nq_BSIM6(Vg,Vfb,Na,Cox,T,phif)
    return 2*nq*Cox*(kb*T)^2*(qs-qd)*(qs+qd+1)
end
function Id_BSIM6(Vg,Vd,Vs,Vfb,Na,Cox,T,phif)
    qs = q_BSIM6(Vg,Vs,Vfb,Na,Cox,T,phif)
    qd = q_BSIM6(Vg,Vd,Vfb,Na,Cox,T,phif)
    nq = nq_BSIM6(Vg,Vfb,Na,Cox,T,phif)
    return μ_BSIM6(Vg,Vs,Vfb,Na,Cox,T,phif)*2*nq*Cox*(kb*T)^2*(qs-qd)*(qs+qd+1)
end

function μ_BSIM6(Vg,Vs,Vfb,Na,Cox,T,phif)
    U0_i = 67e-3
    UTE_i = -1.5
    UA_i = 0.001
    UC_i = 0
    UD_i = 0.001
    UCS_i = 2.0
    UA1_i = 1e-3
    UC1_i = 0.056e-9
    TRatio = T/300.0
    delTemp = T-300.0
    UD1_i = 0
    UCSTE_i = -4.775e-3
    EU_i = 1.5
    U0_t   = U0_i * pow(TRatio, UTE_i);
    UA_t   = UA_i * hypersmooth(1.0 + UA1_i * delTemp - 1.0e-6, 1.0e-3);
    UC_t   = UC_i * hypersmooth(1.0 + UC1_i * delTemp - 1.0e-6, 1.0e-3);
    UD_t   = UD_i * pow(TRatio, UD1_i);
    UCS_t  = UCS_i * pow(TRatio, UCSTE_i);
    nq = nq_BSIM6(Vg,Vfb,Na,Cox,T,phif)
    qs = q_BSIM6(Vg,Vs,Vfb,Na,Cox,T,phif)
    phit = kb*T
    psip = ψp_BSIM6(Vg,Vfb,Na,Cox,T,phif)
    n = n_BSIM6(γC(Na,Cox),phif,Vs,phit)
    qis = qis_BSIM6(nq,n,phit,qs)
    qbs = qbs_BSIM6(nq,phit,Vg,Vfb,qs,n,psip)
    Eeffs = Eeff_BSIM6(qbs,qis,1/Cox/3.9/ϵ₀)
    T2 = pow(0.5 * (1.0 + (qis / qbs)), UCS_t);
    T3 = (UA_t + UC_t * 0.0) * pow(Eeffs, EU_i) + UD_t / T2;
    T4 = 1.0 + T3;
    return U0_t/Smooth(T4,1.0,0.0015)
end

mutable struct BSIM6Model <: TransistorModel
    mos::Union{MOSStructure,Nothing}
    W::Number
    L::Number
    μ_0
    VFB
    Nb
    Cox
    T0
    ϕb
end

function BSIM6Model(mos::MOSStructure,W,L)
    BSIM6Model(mos,W,L,67e-3,VFB(mos),abs(C(mos)),Cox(mos),Temperature(mos),ϕb(mos))
end

pow(a,b) = a^b
Eeff_BSIM6(qba,qia,tox,ϵratio=11.7/3.9,η=0.5) = 1e-8*(qba+η*qia)/(ϵratio*tox)
T0_BSIM6(ϕb,Vs,phit) = hypersmooth(2*ϕb+Vs/phit,1e-3)
n_BSIM6(γ,ϕb,Vs,phit) = 1+ γ/(2*sqrt(T0_BSIM6(ϕb,Vs,phit)))
qbs_BSIM6(nq,phit,Vg,Vfb,qs,n,psip) = Smooth(nq*phit * (Vg/phit-Vfb/phit - psip - 2.0 * qs * (n - 1.0)),0,0.1);
qis_BSIM6(nq,n,phit,qs) = 2*n*nq*phit*qs
hypersmooth(x,c)= 0.5*(x+sqrt(x*x+4.0*c*c))
Smooth(x,x0,deltax) =  0.5 * (x + x0 + sqrt((x - x0) * (x - x0) + 0.25 * deltax * deltax))
id(vg,vd,vs,m::BSIM6Model) = id_BSIM6(vg,vd,vs,m.VFB,m.Nb,m.Cox,m.T0,m.ϕb)
ψp(vg,m::BSIM6Model) = ψp_BSIM6(vg,m.VFB,m.Nb,m.Cox,m.T0,m.ϕb)
qi(Vg,V,m::BSIM6Model) = q_BSIM6(Vg,V,m.VFB,m.Nb,m.Cox,m.T0,m.ϕb)
#id(vg,vd,vs,m::BSIM6Model) = id_BSIM6(vg,vd,vs,m.VFB,m.Nb,m.Cox,m.T0,m.ϕb)
Id(vg,vd,vs,m::BSIM6Model) = m.W/m.L*Id_BSIM6(vg,vd,vs,m.VFB,m.Nb,m.Cox,m.T0,m.ϕb)
gm(vg,vd,vs,m::BSIM6Model) = Zygote.ForwardDiff.derivative(x-> Id(x, vd, vs,m),vg)
gds(vg,vd,vs,m::BSIM6Model) = Zygote.ForwardDiff.derivative(x-> Id(vg, x, vs,m),vd)

## Temperature Model


text:
export draw_curve

"""
`draw_curve(plist)` draws a closed curve through points in the
`plist` which is a 1-dimensional list of complex numbers.
Use `draw_curve(plist,false)` for an open curve.
"""
function draw_curve(pts::Array{Complex{T},1}, closed::Bool = true; opts...) where {T}
    x = real.(pts)
    y = imag.(pts)
    np = length(pts)

    kind = closed ? :closed : :open
    maxt = closed ? np + 1 : np
    Sx = Spline(x, kind)
    Sy = Spline(y, kind)

    fx = funk(Sx)
    fy = funk(Sy)

    plot!(fx, fy, 1, maxt; opts...)
end


text:
token_endpoint = "https://launch.smarthealthit.org/v/r4/auth/token"

client_id = "eyJ0eXAiOiJKV1QiLCJhbGciOiJIUzI1NiJ9.eyJwdWJfa2V5IjoiTFMwdExTMUNSVWRKVGlCUVZVSk1TVU1nUzBWWkxTMHRMUzBLVFVsSlFrbHFRVTVDWjJ0eGFHdHBSemwzTUVKQlVVVkdRVUZQUTBGUk9FRk5TVWxDUTJkTFEwRlJSVUYzWVdvMVoza3hkRXRvVGtOWVdYTjNOV0YzVkFwd1p5OVRja2xuYm1sU1YybElVVmwyZDFseWFrSk5WSHBYUkhkcGQyRTNXbkZLTDNSalRFTk5lR1Y1T0dRMlRHdDRWbkpoYldOb1lqWkdSMnhaZERaUkNtRnZNbkpRWlRSNGJVZ3hkak4zZW1kbVZqaEljbTFUTTI5R2NqbDRjRFJPTm5rNGNtdFdkekZ2Vmtoc2RqZHpVRTV0VlRkell5OHhhU3RJY1ZOUlRFb0thM3BWY1dOQ2FubzFVME14YkhwMlpYaG5jVzkxWjNKNGRUVk5abWwwTmtGd1pHRjFSVGc0U3k5dVNGVk9TM1l2T1ROWmFqTkNaM3BNSzBGV1UwUkpRUW92Ynpsc2VFVlplVmxHV1RBek5HaFJSVmhwVFVFME4yY3ZVRk5ZU20xU2NHWkRXV2hhVUc4MFNtTkdjRXBoU0V4amVGbGhiRmxVZUdSdVZDODVlREJuQ21sQlJETnJjMlZaY20wNFprd3JjRU5EY1V4bFdHbEZXVm94Y0d0R1pqWjFjMkZ5WVZScVMyeGlaSGxMWjJadEwyNWtWemR5V2xkemJVSkZSVVVyUWtVS01GRkpSRUZSUVVJS0xTMHRMUzFGVGtRZ1VGVkNURWxESUV0RldTMHRMUzB0Q2c9PSIsImlzcyI6Imh0dHBzOi8vd2hhdGV2ZXIuc21hcnQvb3VyLXNhbXBsZS1iYWNrZW5kLXNlcnZpY2UiLCJhY2Nlc3NUb2tlbnNFeHBpcmVJbiI6MSwiaWF0IjoxNTEwNzY2MTQzfQ.7YooXIb64Y3_j38n-Gqwa1PqXc-hz-4xJAJF5oqxJVo"

smart_config = BackendServicesConfig(;
    iss = "https://whatever.smart/our-sample-backend-service",
    sub = client_id,
    private_key = JSONWebTokens.RS384(test_private_key),
    scope = "system/*.*",
    token_endpoint = token_endpoint,
)

smart_result = backend_services(smart_config)

@test smart_result isa SMARTBackendServices.BackendServicesResult

access_token = get_fhir_access_token(smart_result)

@test access_token isa AbstractString

@test length(access_token) > 1


text:
using Test
import ClimateMachine
using ClimateMachine.VariableTemplates: @vars, Vars
using ClimateMachine.Mesh.Grids:
    EveryDirection, HorizontalDirection, VerticalDirection
using ClimateMachine.MPIStateArrays: weightedsum
using ClimateMachine.BalanceLaws

import GaussQuadrature
using MPI
using LinearAlgebra

ClimateMachine.init()

@testset "Exponential and Cutoff filter matrix" begin
    let
        # Values computed with:
        #   https://github.com/tcew/nodal-dg/blob/master/Codes1.1/Codes1D/Filter1D.m
        #! format: off
        W = [0x3fe98f3cd0d725e8  0x3fddfd863c6c9a44  0xbfe111110d0fd334  0x3fddbe357bce0b5c  0xbfc970267f929618
             0x3fb608a150f6f927  0x3fe99528b1a1cd8d  0x3fcd41d41f8bae45  0xbfc987d5fabab8d5  0x3fb5da1cd858af87
             0xbfb333332eb1cd92  0x3fc666666826f178  0x3fe999999798faaa  0x3fc666666826f176  0xbfb333332eb1cd94
             0x3fb5da1cd858af84  0xbfc987d5fabab8d4  0x3fcd41d41f8bae46  0x3fe99528b1a1cd8e  0x3fb608a150f6f924
             0xbfc970267f929618  0x3fddbe357bce0b5c  0xbfe111110d0fd333  0x3fddfd863c6c9a44  0x3fe98f3cd0d725e8]
        #! format: on
        W = reinterpret.(Float64, W)

        N = size(W, 1) - 1

        topology = ClimateMachine.Mesh.Topologies.BrickTopology(
            MPI.COMM_SELF,
            -1.0:2.0:1.0,
        )

        grid = ClimateMachine.Mesh.Grids.DiscontinuousSpectralElementGrid(
            topology;
            polynomialorder = N,
            FloatType = Float64,
            DeviceArray = Array,
        )

        filter = ClimateMachine.Mesh.Filters.ExponentialFilter(grid, 0, 32)
        nf = length(filter.filter_matrices)
        @test all(ntuple(i -> filter.filter_matrices[i] ≈ W, nf))
    end

    let
        # Values computed with:
        #   https://github.com/tcew/nodal-dg/blob/master/Codes1.1/Codes1D/Filter1D.m
        #! format: off
        W = [0x3fd822e5f54ecb62   0x3fedd204a0f08ef8   0xbfc7d3aa58fd6968   0xbfbf74682ac4d276
             0x3fc7db36e726d8c1   0x3fe59d16feee478b   0x3fc6745bfbb91e20   0xbfa30fbb7a645448
             0xbfa30fbb7a645455   0x3fc6745bfbb91e26   0x3fe59d16feee478a   0x3fc7db36e726d8c4
             0xbfbf74682ac4d280   0xbfc7d3aa58fd6962   0x3fedd204a0f08ef7   0x3fd822e5f54ecb62]
        #! format: on
        W = reinterpret.(Float64, W)

        N = size(W, 1) - 1

        topology = ClimateMachine.Mesh.Topologies.BrickTopology(
            MPI.COMM_SELF,
            -1.0:2.0:1.0,
        )
        grid = ClimateMachine.Mesh.Grids.DiscontinuousSpectralElementGrid(
            topology;
            polynomialorder = N,
            FloatType = Float64,
            DeviceArray = Array,
        )

        filter = ClimateMachine.Mesh.Filters.ExponentialFilter(grid, 1, 4)
        nf = length(filter.filter_matrices)
        @test all(ntuple(i -> filter.filter_matrices[i] ≈ W, nf))
    end

    let
        T = Float64
        N = (5, 3)
        Nc = (4, 2)

        topology = ClimateMachine.Mesh.Topologies.BrickTopology(
            MPI.COMM_SELF,
            -1.0:2.0:1.0,
        )
        grid = ClimateMachine.Mesh.Grids.DiscontinuousSpectralElementGrid(
            topology;
            polynomialorder = N,
            FloatType = T,
            DeviceArray = Array,
        )

        ξ = ClimateMachine.Mesh.Grids.referencepoints(grid)
        ξ1 = ξ[1]
        ξ2 = ξ[2]
        a1, b1 = GaussQuadrature.legendre_coefs(T, N[1])
        a2, b2 = GaussQuadrature.legendre_coefs(T, N[2])
        V1 = GaussQuadrature.orthonormal_poly(ξ1, a1, b1)
        V2 = GaussQuadrature.orthonormal_poly(ξ2, a2, b2)

        Σ1 = ones(T, N[1] + 1)
        Σ2 = ones(T, N[2] + 1)
        Σ1[(Nc[1]:N[1]) .+ 1] .= 0
        Σ2[(Nc[2]:N[2]) .+ 1] .= 0

        W1 = V1 * Diagonal(Σ1) / V1
        W2 = V2 * Diagonal(Σ2) / V2

        filter = ClimateMachine.Mesh.Filters.CutoffFilter(grid, Nc)
        @test filter.filter_matrices[1] ≈ W1
        @test filter.filter_matrices[2] ≈ W2
    end
end

struct FilterTestModel{N} <: ClimateMachine.BalanceLaws.BalanceLaw end
ClimateMachine.BalanceLaws.vars_state(::FilterTestModel, ::Auxiliary, FT) =
    @vars()
ClimateMachine.BalanceLaws.init_state_auxiliary!(::FilterTestModel, _...) =
    nothing

# Legendre Polynomials
l0(r) = 1
l1(r) = r
l2(r) = (3 * r^2 - 1) / 2
l3(r) = (5 * r^3 - 3r) / 2

low(x, y, z) = l0(x) * l0(y) + 4 * l1(x) * l1(y) + 5 * l1(z) + 6 * l1(z) * l1(x)

high(x, y, z) = l2(x) * l3(y) + l3(x) + l2(y) + l3(z) * l1(y)

filtered(::EveryDirection, dim, x, y, z) = high(x, y, z)
filtered(::VerticalDirection, dim, x, y, z) =
    (dim == 2) ? l2(x) * l3(y) + l2(y) : l3(z) * l1(y)
filtered(::HorizontalDirection, dim, x, y, z) =
    (dim == 2) ? l2(x) * l3(y) + l3(x) : l2(x) * l3(y) + l3(x) + l2(y)

ClimateMachine.BalanceLaws.vars_state(::FilterTestModel{4}, ::Prognostic, FT) =
    @vars(q1::FT, q2::FT, q3::FT, q4::FT)
function ClimateMachine.BalanceLaws.init_state_prognostic!(
    ::FilterTestModel{4},
    state::Vars,
    aux::Vars,
    localgeo,
    filter_direction,
    dim,
)
    (x, y, z) = localgeo.coord
    state.q1 = low(x, y, z) + high(x, y, z)
    state.q2 = low(x, y, z) + high(x, y, z)
    state.q3 = low(x, y, z) + high(x, y, z)
    state.q4 = low(x, y, z) + high(x, y, z)

    if !isnothing(filter_direction)
        state.q1 -= filtered(filter_direction, dim, x, y, z)
        state.q3 -= filtered(filter_direction, dim, x, y, z)
    end
end

@testset "Exponential and Cutoff filter application" begin
    N = 3
    Ne = (1, 1, 1)

    @testset for FT in (Float64, Float32)
        @testset for dim in 2:3
            @testset for direction in (
                EveryDirection,
                HorizontalDirection,
                VerticalDirection,
            )
                brickrange = ntuple(
                    j -> range(FT(-1); length = Ne[j] + 1, stop = 1),
                    dim,
                )
                topl = ClimateMachine.Mesh.Topologies.BrickTopology(
                    MPI.COMM_WORLD,
                    brickrange,
                    periodicity = ntuple(j -> true, dim),
                )

                grid =
                    ClimateMachine.Mesh.Grids.DiscontinuousSpectralElementGrid(
                        topl,
                        FloatType = FT,
                        DeviceArray = ClimateMachine.array_type(),
                        polynomialorder = N,
                    )

                filter = ClimateMachine.Mesh.Filters.CutoffFilter(grid, 2)

                model = FilterTestModel{4}()
                dg = ClimateMachine.DGMethods.DGModel(
                    model,
                    grid,
                    nothing,
                    nothing,
                    nothing;
                    state_gradient_flux = nothing,
                )

                @testset for target in ((1, 3), (:q1, :q3))
                    Q = ClimateMachine.DGMethods.init_ode_state(
                        dg,
                        nothing,
                        dim,
                    )
                    ClimateMachine.Mesh.Filters.apply!(
                        Q,
                        target,
                        grid,
                        filter,
                        direction = direction(),
                    )
                    P = ClimateMachine.DGMethods.init_ode_state(
                        dg,
                        direction(),
                        dim,
                    )
                    @test Array(Q.data) ≈ Array(P.data)
                end
            end
        end
    end
end

ClimateMachine.BalanceLaws.vars_state(
    ::FilterTestModel{1},
    ::Prognostic,
    FT,
) where {N} = @vars(q::FT)
function ClimateMachine.BalanceLaws.init_state_prognostic!(
    ::FilterTestModel{1},
    state::Vars,
    aux::Vars,
    localgeo,
)
    (x, y, z) = localgeo.coord
    state.q = abs(x) - 0.1
end

@testset "TMAR filter application" begin

    N = 4
    Ne = (2, 2, 2)

    @testset for FT in (Float64, Float32)
        @testset for dim in 2:3
            brickrange =
                ntuple(j -> range(FT(-1); length = Ne[j] + 1, stop = 1), dim)
            topl = ClimateMachine.Mesh.Topologies.BrickTopology(
                MPI.COMM_WORLD,
                brickrange,
                periodicity = ntuple(j -> true, dim),
            )

            grid = ClimateMachine.Mesh.Grids.DiscontinuousSpectralElementGrid(
                topl,
                FloatType = FT,
                DeviceArray = ClimateMachine.array_type(),
                polynomialorder = N,
            )

            model = FilterTestModel{1}()
            dg = ClimateMachine.DGMethods.DGModel(
                model,
                grid,
                nothing,
                nothing,
                nothing;
                state_gradient_flux = nothing,
            )

            @testset for target in ((1,), (:q,), :)
                Q = ClimateMachine.DGMethods.init_ode_state(dg)

                initialsumQ = weightedsum(Q)
                @test minimum(Q.realdata) < 0

                ClimateMachine.Mesh.Filters.apply!(
                    Q,
                    target,
                    grid,
                    ClimateMachine.Mesh.Filters.TMARFilter(),
                )

                sumQ = weightedsum(Q)

                @test minimum(Q.realdata) >= 0
                @test isapprox(initialsumQ, sumQ; rtol = 10 * eps(FT))
            end
        end
    end
end

max_forks_repo_name: jm-c/CLIMA


text:
# This file was generated, do not modify it. # hide
println(rpad(" Feature", 11), "| ", "Coefficient")
println("-"^24)
for (i, name) in enumerate(names(X))
    println(rpad("$name", 11), "| ", round(fp.coefs[i], sigdigits=3))
end
println(rpad("Intercept", 11), "| ", round(fp.intercept, sigdigits=3))


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module Modia3D.Shapes

Visual elements used for animation. The visual elements are passed to an
external renderer. Currently, the (free) community edition and the
(commercial) professional editions of the DLR Visualization library are supported
(not all visualization elements of this library are yet interfaced to Modia3D).
Modia3D is designed so that other renderers can be supported as well.

Visualization elements with shape and visualization material:

![VisuElements](../../resources/images/VisuElementsWith.jpg)

Visualization elements that do not have a visualization material:

![VisuElementsWithout](../../resources/images/VisuElementsWithout.jpg)

# Main developers
Andrea Neumayr and Martin Otter, [DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
```
"""
module Shapes

using  StaticArrays
import Modia3D
import Modia3D.Basics
import Modia3D.Frames

import JSON

# Colors and material
export rgb, defaultColor, VisualMaterial, colorPalette, visualMaterialPalette

# Visual elements
export Sphere, Ellipsoid, Box, Cylinder, Capsule, Beam, Cone
export Spring, GearWheel, CoordinateSystem, Grid, FileMesh
export changeDimensionsBox
export Font, TextShape, Screen, XY_Plane, XZ_Plane, YZ_Plane, Left, Right, Center

export Visual

# Concave properties
export ConcaveProperties

# Solid geometries
export Solid
export volume, centroid, bottomArea, topArea, longestEdge, lengthGeo, inertiaMatrix

export supportPoint_Sphere, supportPoint_Ellipsoid, supportPoint_Box, supportPoint_Cylinder, supportPoint_Cone, supportPoint_Capsule, supportPoint_Beam, supportPoint_FileMesh

export supportPoint_i_Box, supportPoint_i_Cylinder, supportPoint_i_Cone, supportPoint_i_Capsule,supportPoint_i_Beam, supportPoint_i_Ellipsoid, supportPoint_i_FileMesh

# Solid materials
export InternalMassProperties, MassProperties, SolidMaterial, ContactMaterialElastic, dummyMassProperties

export MassPropertiesFromShape, MassPropertiesFromShapeAndMass

export solidMaterial, solidMaterialPalette, defaultContactMaterial
export regularize, resultantCoefficientOfRestitution, resultantDampingCoefficient
export solidMaterialPairsPalette, CommonCollisionProperties, getCommonCollisionProperties

# Observer materials

export contactPairMaterialPalette, ElasticContactPairMaterial, TwoNamesKey

export NoContactPairMaterial, ObserverContactPairMaterial, ImpulseContactPairMaterial, WheelRailContactPairMaterial

export getShapeKind

# for identifying shapes at compile time
export ShapeKind, UndefinedShapeKind, SphereKind, EllipsoidKind, BoxKind, CylinderKind, ConeKind, CapsuleKind, BeamKind, FileMeshKind, CoordinateSystemKind, GridKind, SpringKind, GearWheelKind, TextKind

@enum ShapeKind UndefinedShapeKind SphereKind EllipsoidKind BoxKind CylinderKind ConeKind CapsuleKind BeamKind FileMeshKind CoordinateSystemKind GridKind SpringKind GearWheelKind TextKind

include("color.jl")
include("visualMaterial.jl")
include("text.jl")
include("geometry.jl")
include("shape.jl")
include("visual.jl")
include("concaveProperties.jl")

include("inertiaTensorAndVolume.jl")
include("computePropertiesFileMes.jl")
include("boundingBoxes.jl")
include("utilities.jl")
include("solidMaterial.jl")
include("massProperties.jl")
include("contactPairMaterials.jl")
include("setCollisionSmoothingRadius.jl")
include("solid.jl")

end


text:
# The code is written for Julia v0.6 and later.
# https://github.com/JuliaLang/julia

struct GoSlice
    arr::Ptr{Void}
    len::Int64
    cap::Int64
end
GoSlice(a::Vector, cap=10) = GoSlice(pointer(a), length(a), cap)

struct GoStr
    p::Ptr{Cchar}
    len::Int64
end
GoStr(s::String) = GoStr(pointer(s), length(s))

const libawesome = "awesome.so"

Add(x,y) = ccall((:Add, libawesome), Int,(Int,Int), x,y)
Cosine(x) = ccall((:Cosine, libawesome), Float64, (Float64,), x)
function Sort(vals)
    ccall((:Sort, libawesome), Void, (GoSlice,), GoSlice(vals))
    return vals # for convenience
end
Log(msg) = ccall((:Log, libawesome), Int, (GoStr,), GoStr(msg))

for ex in [:(Add(12, 9)),:(Cosine(1)), :(Sort([77,12,5,99,28,23]))]
    println("awesome.$ex = $(eval(ex))")
end
Log("Hello from Julia!")

#=
> julia client.jl
awesome.Add(12, 9) = 21
awesome.Cosine(1) = 0.5403023058681398
awesome.Sort([77, 12, 5, 99, 28, 23]) = [5, 12, 23, 28, 77, 99]
Hello from Julia!
=#

max_stars_repo_name: dpurfield/go-cshared-examples

max_issues_repo_name: dpurfield/go-cshared-examples

max_forks_repo_name: dpurfield/go-cshared-examples


text:
using Turing, Random
using Turing: ParticleContainer, weights, resample!,
    effectiveSampleSize, Trace, current_trace, VarName,
    Sampler, consume, produce, copy, fork
using Turing.Core: logZ
using Test

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

@testset "container.jl" begin
    @turing_testset "copy particle container" begin
        pc = ParticleContainer{Trace}(x -> x * x)
        newpc = copy(pc)

        @test newpc.logE        == pc.logE
        @test newpc.logWs       == pc.logWs
        @test newpc.conditional == pc.conditional
        @test newpc.n_consume   == pc.n_consume
    end
    @turing_testset "particle container" begin
        n = Ref(0)

        alg = PG(5)
        spl = Turing.Sampler(alg, empty_model())
        dist = Normal(0, 1)

        function fpc()
            t = TArray(Float64, 1);
            t[1] = 0;
            while true
                ct = current_trace()
                vn = VarName(gensym(), :x, "[$n]", 1)
                Turing.assume(spl, dist, vn, ct.vi); n[] += 1;
                produce(0)
                vn = VarName(gensym(), :x, "[$n]", 1)
                Turing.assume(spl, dist, vn, ct.vi); n[] += 1;
                t[1] = 1 + t[1]
            end
        end

        pc = ParticleContainer{Trace}(fpc)
        model = Turing.Model{(:x,),()}(fpc, NamedTuple(), NamedTuple())
        tr = Trace(pc.model, model, spl, Turing.VarInfo())
        push!(pc, tr)
        tr = Trace(pc.model, model, spl, Turing.VarInfo())
        push!(pc, tr)
        tr = Trace(pc.model, model, spl, Turing.VarInfo())
        push!(pc, tr)

        @test weights(pc) == [1/3, 1/3, 1/3]
        @test logZ(pc) ≈ log(3)
        @test pc.logE ≈ log(1)

        @test consume(pc) == log(1)

        resample!(pc)
        @test pc.num_particles == length(pc)
        @test weights(pc) == [1/3, 1/3, 1/3]
        @test logZ(pc) ≈ log(3)
        @test pc.logE ≈ log(1)
        @test effectiveSampleSize(pc) == 3

        @test consume(pc) ≈ log(1)
        resample!(pc)
        @test consume(pc) ≈ log(1)
    end
    @turing_testset "trace" begin
        n = Ref(0)

        alg = PG(5)
        spl = Turing.Sampler(alg, empty_model())
        dist = Normal(0, 1)
        function f2()
            t = TArray(Int, 1);
            t[1] = 0;
            while true
                ct = current_trace()
                vn = VarName(gensym(), :x, "[$n]", 1)
                Turing.assume(spl, dist, vn, ct.vi); n[] += 1;
                produce(t[1]);
                vn = VarName(gensym(), :x, "[$n]", 1)
                Turing.assume(spl, dist, vn, ct.vi); n[] += 1;
                t[1] = 1 + t[1]
            end
        end

        # Test task copy version of trace
        model = Turing.Model{(:x,),()}(f2, NamedTuple(), NamedTuple())
        tr = Trace(f2, model, spl, Turing.VarInfo())

        consume(tr); consume(tr)
        a = fork(tr);
        consume(a); consume(a)

        @test consume(tr) == 2
        @test consume(a) == 4
    end
end


text:
import PDESolver: evalJacobian, evalJacobianStrong

"""
  Euler implementation of `evalJacobian`.
"""
function evalJacobian(mesh::AbstractMesh, sbp::AbstractOperator, eqn::EulerData,
                      opts::Dict, assembler::AssembleElementData, t=0.0;
                      start_comm=false)

  #println("\nEntered evalJacobian")
  time = eqn.params.time
  eqn.params.t = t  # record t to params
  myrank = mesh.myrank

  @assert opts["parallel_type"] == 2
#  println("entered evalResidual")
  time.t_send += @elapsed if start_comm
    setParallelData(eqn.shared_data, opts["parallel_data"])
    startSolutionExchange(mesh, sbp, eqn, opts)
  end

  time.t_dataprep_diff += @elapsed dataPrep_diff(mesh, sbp, eqn, opts)
  #println("dataPrep @time printed above")

  time.t_volume_diff += @elapsed if opts["addVolumeIntegrals"]
    evalVolumeIntegrals_diff(mesh, sbp, eqn, opts, assembler)
    #println("volume integral diff @time printed above")
  end

  if opts["use_GLS"]
    error("use_GLS not supported by evalJacobian")
  end

  time.t_bndry_diff += @elapsed if opts["addBoundaryIntegrals"]
    evalBoundaryIntegrals_diff(mesh, sbp, eqn, opts, assembler)
    #println("boundary integral diff @time printed above")
  end

  time.t_stab_diff += @elapsed if opts["addStabilization"]
    addStabilization_diff(mesh, sbp, eqn, opts, assembler)
#    println("stabilizing @time printed above")
  end

  time.t_face_diff += @elapsed if mesh.isDG && opts["addFaceIntegrals"]
    evalFaceIntegrals_diff(mesh, sbp, eqn, opts, assembler)
    #println("face integral diff @time printed above")
  end

  time.t_shock_diff += @elapsed if opts["addShockCapturing"]
    evalShockCapturing_diff(mesh, sbp, eqn, opts, assembler)
  end

  time.t_sharedface_diff += @elapsed if mesh.commsize > 1
    evalSharedFaceIntegrals_diff(mesh, sbp, eqn, opts, assembler)
#    println("evalSharedFaceIntegrals @time printed above")
  end

  time.t_source_diff += @elapsed evalSourceTerm_diff(mesh, sbp, eqn, opts, assembler)
  #println("source integral diff @time printed above")

  # apply inverse mass matrix to eqn.res, necessary for CN
  #TODO: this will have to be done at the element level
  if opts["use_Minv"]
#    error("use_Minv must be applied at the element level")
#    applyMassMatrixInverse3D(mesh, sbp, eqn, opts, eqn.res)
  end

  return nothing
end

function evalJacobianStrong(mesh::AbstractMesh, sbp::AbstractOperator,
                      eqn::EulerData,
                      opts::Dict, assembler::AssembleElementData, t=0.0;
                      start_comm=false)
# currently this function neglects the SAT terms (including boundary conditions)

  time = eqn.params.time
  eqn.params.t = t  # record t to params
  myrank = mesh.myrank

  time.t_volume_diff += @elapsed if opts["addVolumeIntegrals"]
    calcVolumeIntegralsStrong_nopre_diff(mesh, sbp, eqn, opts, assembler)
#    println("volume integral @time printed above")
  end

  return nothing
end

"""
  Differentiated version of [`dataPrep`](@ref).  The Jacobian computation
  does use any of the arrays of precomputed values (options keys
  `precompute_*`).  Those keys can be specified by the user and will affect
  [`evalResidual`](@ref), but not [`evalJacobian`](@ref).

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
"""
function dataPrep_diff(mesh::AbstractMesh{Tmsh},
                      sbp::AbstractOperator,
                      eqn::AbstractEulerData{Tsol, Tres}, opts) where {Tmsh, Tsol, Tres}
# gather up all the data needed to do vectorized operatinos on the mesh
# calculates all mesh wide quantities in eqn
# this is almost the exact list of everything we *shouldn't* be storing, but
# rather recalculating on the fly

#println("Entered dataPrep()")

#  println("typeof(eqn) = ", typeof(eqn))
#  println("typeof(eqn.params) = ", typeof(eqn.params))

  # apply filtering to input
  if eqn.params.use_filter
    applyFilter(mesh, sbp, eqn, eqn.q, opts)
  end

  # zero out res
#  fill!(eqn.res, 0.0)
#  fill!(eqn.res_edge, 0.0)

  getAuxVars(mesh, eqn)
#  println("  getAuxVars @time printed above")

  if opts["check_density"]
    checkDensity(eqn, mesh)
#    println("  checkDensity @time printed above")
  end

  if opts["check_pressure"]
#    throw(ErrorException("I'm done"))
    checkPressure(eqn, mesh)
#    println("  checkPressure @time printed above")
  end

  # calculate fluxes

  if opts["use_staggered_grid"]
    error("staggered grid not supported by calcJacobian()")
  end

  # the jacobian computation doesn't precompute anything

  if eqn.params.use_edgestab
    error("Use_edgestab not supported by calcJacobian()")
    stabscale(mesh, sbp, eqn)
  end
#  println("  stabscale @time printed above")

  return nothing
end # end function dataPrep

"""
  Differentiated version of [`evalVolumeIntegrals`](@ref).  Throws an error
  for unsupported schemes.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * assembler
"""
function evalVolumeIntegrals_diff(mesh::AbstractMesh{Tmsh},
                       sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                       opts, assembler::AssembleElementData
                       ) where {Tmsh,  Tsol, Tres, Tdim}

  integral_type = opts["volume_integral_type"]

  if opts["Q_transpose"] == false
    error("Q_transpose == false not supported by evalJacobian")
  end  # end if Q_transpose

  if integral_type == 1  # regular volume integrals
    calcVolumeIntegrals_nopre_diff(mesh, sbp, eqn, opts, assembler)
  elseif integral_type == 2  # entropy stable formulation
    calcVolumeIntegralsSplitForm_diff(mesh, sbp, eqn, opts,
                                      eqn.volume_flux_func_diff, assembler)
  else
    throw(ErrorException("Unsupported volume integral type = $integral_type"))
  end

end  # end evalVolumeIntegrals

"""
  Differentiated version of [`evalBoundaryIntegrals`](@ref).

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * assembler.
"""
function evalBoundaryIntegrals_diff(mesh::AbstractMesh{Tmsh},
                               sbp::AbstractOperator,
                               eqn::EulerData{Tsol, Tres, Tdim}, opts,
                               assembler::AssembleElementData) where {Tmsh, Tsol, Tres, Tdim}

  if mesh.isDG
    # when precompute_boundary_flux == false, this fuunction does the
    # integration too, updating res
    getBCFluxes_diff(mesh, sbp, eqn, opts, assembler)
  else
    error("CG meshes not supported by evalJacobian()")
  end

  return nothing

end  # end evalBoundaryIntegrals

"""
  Differentiated version of [`addStabilization`](@ref).

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * assembler
"""
function addStabilization_diff(mesh::AbstractMesh{Tmsh},
                          sbp::AbstractOperator, eqn::EulerData{Tsol}, opts,
                          assembler::AssembleElementData) where {Tmsh,  Tsol}

#  println("==== start of addStabilization ====")

  # ----- Edge Stabilization -----#
  if eqn.params.use_edgestab  #TODO: use options dictionary instead
#    println("applying edge stabilization")
    error("edge stabilization not supported by evalJacobian()")
  end

  # ----- Filtering -----
  if eqn.params.use_res_filter
    error("use_res_filter not supported by evalJacobian()")
  end

  # ----- Artificial Dissipation -----
  if eqn.params.use_dissipation
    error("use_dissipation not supported by evalJacobian()")
  end

  if opts["use_GLS2"]
    error("use_GLS2 not supported by evalJacobian()")
     applyGLS3(mesh, sbp, eqn, opts)
#    test_GLS(mesh, sbp, eqn, opts)
  end

  if opts["use_lps"]
    applyLPStab_diff(mesh, sbp, eqn, opts, assembler)
  end
#  println("==== end of addStabilization ====")

  return nothing
end

"""
  Differentiated version of [`evalFaceIntegrals_diff`](@ref).  Throws an
  error for unsupported schemes.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * assembler
"""
function evalFaceIntegrals_diff(mesh::AbstractDGMesh{Tmsh},
                    sbp::AbstractOperator,
                    eqn::EulerData{Tsol}, opts,
                    assembler::AssembleElementData) where {Tmsh, Tsol}

  face_integral_type = opts["face_integral_type"]
  if face_integral_type == 1
      calcFaceIntegral_nopre_diff(mesh, sbp, eqn, opts, eqn.flux_func_diff, mesh.interfaces, assembler)

  elseif face_integral_type == 2

    if opts["use_staggered_grid"]
      error("explicit jacobian calculation does not support staggered grids")
    else
      getFaceElementIntegral_diff(mesh, sbp, eqn,
                  eqn.face_element_integral_func,
                  eqn.flux_func_diff, mesh.sbpface, mesh.interfaces, assembler)
    end

  else
    throw(ErrorException("Unsupported face integral type = $face_integral_type"))
  end

  return nothing
end

"""
  Differentiated version of [`evalSharedFaceIntegrals`](@ref).
  Currently only supports `opts[parallel_data] == element` because it is
  used for Jacobian calculation.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * assembler
"""
function evalSharedFaceIntegrals_diff(mesh::AbstractDGMesh, sbp, eqn, opts,
                                      assembler::AssembleElementData)

  face_integral_type = opts["face_integral_type"]
  eqn.assembler = assembler  # stash this here to get retrieved later
  if face_integral_type == 1

    finishExchangeData(mesh, sbp, eqn, opts, eqn.shared_data, calcSharedFaceIntegrals_element_diff)

  elseif face_integral_type == 2
    finishExchangeData(mesh, sbp, eqn, opts, eqn.shared_data, calcSharedFaceElementIntegrals_element_diff)
  else
    throw(ErrorException("unsupported face integral type = $face_integral_type"))
  end

  # we don't want to keep a reference to the matrix inside the assembler
  # that would prevent it from getting freed
  eqn.assembler = NullAssembleElementData
  return nothing
end

"""
  Differentiated version of [`evalSourceTerm`](@ref).  Currently, source
  terms are not supported.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * assembler
"""
function evalSourceTerm_diff(mesh::AbstractMesh{Tmsh},
                     sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                     opts, assembler::AssembleElementData) where {Tmsh, Tsol, Tres, Tdim}

  # source terms are not a function of q, so they don't affect the Jacobian

  #=
  # placeholder for multiple source term functionality (similar to how
  # boundary conditions are done)
  if opts["use_src_term"]
    error("source terms not supported by evalJacobian()")
  end
  =#
  return nothing
end  # end function

"""
  Differentiated version of [`evalShockCapturing`](@ref).

 **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * assembler
"""
function evalShockCapturing_diff(mesh::AbstractMesh{Tmsh},
                       sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                       opts, assembler::AssembleElementData
                       ) where {Tmsh,  Tsol, Tres, Tdim}

  capture = eqn.shock_capturing
  applyShockCapturing_diff(mesh, sbp, eqn, opts, capture, assembler)

  return nothing
end


text:
include("Utils.jl")
include("Testing.jl")
using Flux
using FFTW
using LinearAlgebra
using BSON

####################
# TRAINING DATA
####################
function hannWindow(window)
    N = length(window)
    return [window[n] * 0.5 * (1 - cos(2π * n / N)) for n = 1:N]
end

function spectrum(window)
    S = fft(window)
    S = S[1:Int(length(S) / 2)]
    scale = 1 / length(window)
    global S = [sqrt(real(n)^2 + imag(n)^2) * scale for n in S]
end

function test_case(f, Fs, W)
    K = Int(floor(Fs / 2 / f)) - 1
    frequencies = [f * i for i = 1:K]
    test = Testing.generate_test(f, W, Fs)

    S = spectrum(hannWindow(test.signal))
    bin_size = Fs / 2 / length(S)

    correct_peaks = frequencies ./ bin_size .+ 1
    assumed_peaks = Int.(round.(correct_peaks))

    inputs = []
    answers = []
    for i = 1:length(assumed_peaks)
        idx = assumed_peaks[i]
        if S[idx-1] < S[idx] && S[idx] > S[idx+1]
            push!(inputs, S[[(idx - 1), (idx + 1)]] ./ S[idx])
            push!(answers, correct_peaks[i] - assumed_peaks[i])
        end
    end
    return inputs, answers
end

function get_training_data(N)
    W_list = [1024, 2048, 4096]
    Fs = 48000
    inputs = [[[] for n = 1:N] for W in W_list]
    answers = [[[] for n = 1:N] for W in W_list]
    frequencies = [[[] for n = 1:N] for W in W_list]

    for w = 1:length(W_list)
        W = W_list[w]
        bin_size = Fs / W
        low_note = bin_size * 4
        high_note = 4800

        note_range = Utils.getRange(low_note, high_note)
        Threads.@threads for n = 1:N
            for f in [note_range[(n-1)%length(note_range)+1], rand(low_note:0.01:high_note)]
                x, y = test_case(f, Fs, W)

                push!(inputs[w][n], x)
                push!(answers[w][n], y)
                push!(frequencies[w][n], f)
            end
        end
    end
    X = reduce(hcat, reduce(vcat, reduce(vcat, reduce(vcat, inputs))))
    Y = float.(reduce(vcat, reduce(vcat, reduce(vcat, answers))))
    frequencies = float.(reduce(vcat, reduce(vcat, reduce(vcat, answers))))
    return X, Y, frequencies
end

###############################
# MODEL
###############################
function get_params(param_sizes)
    W = Matrix{Float64}[]
    b = Vector{Float64}[]
    for layer in param_sizes
        fan_in = layer[1]
        fan_out = layer[2]
        limit = sqrt(6 / (fan_in + fan_out))
        push!(W, rand(-limit:0.001:limit, fan_in, fan_out))
        push!(b, zeros(layer[2]))
    end
    return W, b
end

function params_to_vector(W, b)
    reduce(vcat, [[W[i]..., b[i]...] for i = 1:length(W)])
end

function vector_to_params(β, param_sizes)
    W = Matrix{Float64}[]
    b = Vector{Float64}[]
    idx = 1
    for layer in param_sizes
        len = layer[1] * layer[2]
        push!(W, reshape(β[idx:idx+len-1], layer))
        idx += len
        push!(b, β[idx:idx+layer[2]-1])
        idx += layer[2]
    end
    return W, b
end

function training_model(x, β)
    W, b = vector_to_params(β, param_sizes)
    vec(model(x, W, b))
end

function model(x, W, b)
    W[3]' * tanh.(W[2]' * tanh.(W[1]' * x .+ b[1]) .+ b[2]) .+ b[3]
end

###############################
# Training
###############################
param_sizes = [(2, 3), (3, 4), (4, 1)]
X, Y, frequencies = get_training_data(2500)
W, b = get_params(param_sizes)
β = params_to_vector(W, b)

using LsqFit
fit = curve_fit(training_model, X, Y, β)
β = coef(fit)

W, b = vector_to_params(β, param_sizes)
BSON.@save "Models/Params.bson" W b

###############################
# EVALUATION
###############################
X, Y, frequencies = get_training_data(2500)
sort_order = sortperm(frequencies)
f_sort = frequencies[sort_order]
X_sort = X[:, sort_order]
Y_sort = Y[sort_order]
guesses = vec(model(X_sort, W, b))

err = Y_sort .- guesses
avg_err = sum(err) / length(err)
maximum(err)

using GLMakie
lines(err)


text:
module UpperEnvelope

export PiecewiseLinear
export compute_envelope

include("piecewise_linear.jl")

end # module


text:
push!(Base.LOAD_PATH, joinpath(dirname(@__FILE__),"..",".."))

using JuMP
using GLPK
using ParameterJuMP
using Test
using Profile
using ProfileView

function bench1_p(N::Int, M::Int)
    m_slave = ModelWithParams(GLPK.Optimizer)

    x = add_parameters(m_slave, 4.0*ones(M))
    @variable(m_slave, y[1:N])
    for i in 1:100
        @constraint(m_slave, sum(3.0*y[i] for i in 1:N) >= 2.0 - sum(7.0*x[i] for i in 1:M))
    end
    for i in 1:300
        @constraint(m_slave, sum(3.0*y[i] for i in 1:N) >= 2.0)
    end
end
function bench1_v(N::Int, M::Int)
    m_slave = Model(GLPK.Optimizer)

    @variable(m_slave, x[1:M] == 4.0)
    @variable(m_slave, y[1:N])
    for i in 1:100
        @constraint(m_slave, sum(3.0*y[i] for i in 1:N) >= 2.0 - sum(7.0*x[i] for i in 1:M))
    end
    for i in 1:300
        @constraint(m_slave, sum(3.0*y[i] for i in 1:N) >= 2.0)
    end
end
function bench2_p(N::Int)
    m_slave = ModelWithParams(GLPK.Optimizer)

    x = add_parameters(m_slave, 4.0*ones(N))
    @variable(m_slave, y[1:N])
    for i in 1:100
        @constraint(m_slave, sum(3.0*y[i] + 3.0 + 7.0*x[i] for i in 1:N) >= 2.0)
    end
    for i in 1:300
        @constraint(m_slave, sum(3.0*y[i] for i in 1:N) >= 2.0)
    end
end
function bench2_v(N::Int)
    m_slave = Model(GLPK.Optimizer)

    @variable(m_slave, x[1:N] == 4.0)
    @variable(m_slave, y[1:N])
    for i in 1:100
        @constraint(m_slave, sum(3.0*y[i] + 3.0 + 7.0*x[i] for i in 1:N) >= 2.0)
    end
    for i in 1:300
        @constraint(m_slave, sum(3.0*y[i] for i in 1:N) >= 2.0)
    end
end

bench1_p(10, 3)
bench1_v(10, 3)
bench2_p(10)
bench2_v(10)

N = 1000
M = 1000
GC.gc()
@time bench1_p(N, M)
GC.gc()
@time bench1_v(N, M)
GC.gc()
@time bench2_p(N)
GC.gc()
@time bench2_v(N)
;

function benchmark()
    # Any setup code goes here.

    # Run once, to force compilation.
    println("======================= First run:")
    # srand(666)
    @time bench2_p(10)

    # Run a second time, with profiling.
    println("\n\n======================= Second run:")
    # srand(666)
    Profile.init(delay=0.0001)
    Profile.clear()
    # clear_malloc_data()
    @profile @time bench2_p(1000)

    # Write profile results to profile.bin.
    # r = Profile.retrieve()
    # f = open("profile.bin", "w")
    # serialize(f, r)
    # close(f)

    ProfileView.view()
end
benchmark()


text:
module Coincheck

using JSON
using Nettle
using WebSocketClient
import WebSocketClient: on_text
import Requests: URI

include("common.jl")
include("HttpUtil.jl")

include("http.jl")
include("websocket.jl")

end


text:
include("preamble.jl")

T = fill(0.0,1,1)
c = [1.0]
model = BoundedFluidQueue(T,c,1.0)
hx = 1.0
nodes = 0.0:hx:model.b

x_lims = (-0.1,1.1)
n_evals = 10_001
## First, lets examin closing vectors for QBDRAP
orders = 1:2:21
CDFs = [
    (z,i)->Float64(z>0.5);
    (z,i)->(z<0.5)*2*z+(z>=0.5);
    (z,i)->(z>0.5)*(z^2/2-z/2+1/8)*8;
    (z,i)->z;
    (z,i)->z^2;
    (z,i)->-2*z^3+3*z^2;
    (z,i)->(1-exp(-3*z))/(1-exp(-3));
    (z,i)->2*(z/2 - sin(2*pi*z)/(4*pi));
    (z,i)->z + sin(4*pi*z)/(4*pi);
    (z,i)->cdf(ConcentratedMatrixExponential(5,mean=0.5))(z)
    ]
PDFs = [
    (z,i)->NaN;
    (z,i)->2*Float64(z<0.5);
    (z,i)->(z>0.5)*8*(z-0.5);
    (z,i)->1.0;
    (z,i)->2.0*z;
    (z,i)->-6*z^2+6*z;
    (z,i)->(3*exp(-3*z))/(1-exp(-3));
    (z,i)->cos(2*π*(z+0.5))+1;
    (z,i)->cos(4*π*(z+0.5))+1;
    (z,i)->pdf(ConcentratedMatrixExponential(5,mean=0.5))(z)
    ]
distns = [x->interior_point_mass(0.5,1,x);
    [ x->SFMDistribution(PDFs[i],x,TrapezoidRule;fun_evals=n_evals) for i in 2:length(PDFs)]]
    # x->SFMDistribution(PDFs[2],x,TrapezoidRule;fun_evals=2001);
    # x->SFMDistribution(PDFs[3],x,TrapezoidRule;fun_evals=2001);
    # [ x->SFMDistribution(PDFs[i],x) for i in 4:length(PDFs)]]
# x->SFMDistribution((z,i)->1.0,x),
# x->SFMDistribution((z,i)->2.0*z,x),
# x->SFMDistribution((z,i)->-6*z^2+6*z,x),
# x->SFMDistribution((z,i)->3*exp(-3*z)/(1-exp(-3)),x),
# x->SFMDistribution((z,i)->cos(2*π*(z+0.5))+1,x),
# x->SFMDistribution((z,i)->cos(4*π*(z+0.5))+1,x),
# x->SFMDistribution((z,i)->pdf(ConcentratedMatrixExponential(5,mean=0.5))(z),x),
# ]
for (func_count,init_dist_fun) in enumerate(distns)
    L1_cdf_errors = DataFrame(unnormalised = Float64[], naive = Float64[], normalised = Float64[])
    ks_errors = DataFrame(unnormalised = Float64[], naive = Float64[], normalised = Float64[])
    L1_pdf_errors = DataFrame(unnormalised = Float64[], naive = Float64[], normalised = Float64[])
    L2_pdf_errors = DataFrame(unnormalised = Float64[], naive = Float64[], normalised = Float64[])

    pth = mkpath((@__FILE__)[1:end-3]*"/func_count_"*string(func_count)*"/QDBRAP")
    mkpath(pth*"/data")
    mkpath(pth*"/figs")

    for plt_type in (cdf,pdf)
        p = plot(layout=(3,length(orders)))
        linetypes = [:solid,:solid,:solid]
        c_order = 0
        for order in orders
            c_order += 1
            c_mesh = 0
            rec_types = (eval(Symbol(:unnormalised_closing_operator_,plt_type)),
                eval(Symbol(:naive_normalised_closing_operator_,plt_type)),
                eval(Symbol(:normalised_closing_operator_,plt_type)))
            L1_cdf_error_row = [ 0.0 0.0 0.0]
            ks_error_row = [ 0.0 0.0 0.0]
            L1_pdf_error_row = [ 0.0 0.0 0.0]
            L2_pdf_error_row = [ 0.0 0.0 0.0]
            for reconstruction in rec_types
                c_mesh += 1
                mesh = FRAPMesh(nodes,order)

                dq = DiscretisedFluidQueue(model,mesh)

                d0 = init_dist_fun(dq)
                V = DiscretisedFluidQueues.vandermonde(order)

                rec = plt_type(d0, reconstruction)

                ## errors
                (plt_type==cdf) && (L1_cdf_error_row[c_mesh] = log10(DiscretisedFluidQueues.Lp(x->rec(x,1),x->CDFs[func_count](x,1),range(0,1,length=n_evals))))
                (plt_type==cdf) && (ks_error_row[c_mesh] = log10(DiscretisedFluidQueues.kolmogorov_smirnov(x->rec(x,1),x->CDFs[func_count](x,1),range(0,1,length=n_evals))))
                (plt_type==pdf) && (L1_pdf_error_row[c_mesh] = log10(DiscretisedFluidQueues.Lp(x->rec(x,1),x->PDFs[func_count](x,1),range(0,1,length=n_evals))))
                (plt_type==pdf) && (L2_pdf_error_row[c_mesh] = log10(DiscretisedFluidQueues.Lp(x->rec(x,1),x->PDFs[func_count](x,1),range(0,1,length=n_evals),2)))

                ## plotting...
                y_lim_vals = (plt_type==pdf) ? (-3.0,9.0) : (-0.2,1.4)
                y_tick_vals = ((plt_type==pdf) ? (y_lim_vals[1]:3:y_lim_vals[2]-1) : (0:0.5:1))
                std_plot(args...; kwargs...) = plot!(
                    p.layout.grid[c_mesh,c_order],
                    x->rec(x,1),
                    nodes[1]+eps(), nodes[end]-eps();
                    ylim=y_lim_vals,
                    linestyle=linetypes[c_mesh],
                    xticks=false,
                    yticks=false,
                    grid=false,
                    kwargs...,
                )
                if c_order == 1
                    std_plot(;label=false,ylabel=string(reconstruction)[1:findfirst("_",string(reconstruction))[1]-1],legend=false)
                    yticks!(p.layout.grid[c_mesh,c_order], y_tick_vals)
                elseif c_order == length(orders)
                    std_plot(;label=string(reconstruction)[1:findfirst("_",string(reconstruction))[1]-1],legend=false)#:outertopright)
                else
                    std_plot(;label=false,legend=false)
                end
                if c_mesh == 1
                    plot!(p.layout.grid[c_mesh,c_order]; title=string(order))
                elseif c_mesh == 3
                    xticks!(p.layout.grid[c_mesh,c_order], nodes)
                    plot!(p.layout.grid[c_mesh,c_order]; xlabel="x")
                end
                std_plot2(args... ; kwargs...) = if plt_type==cdf
                    plot!(
                        p.layout.grid[c_mesh,c_order],
                        x->CDFs[func_count](x,1),
                        nodes[1], nodes[end],
                        args...;
                        linestyle=:dot,
                        kwargs...,
                    )
                elseif plt_type==pdf
                    plot!(
                        p.layout.grid[c_mesh,c_order],
                        x->PDFs[func_count](x,1),
                        nodes[1], nodes[end],
                        args...;
                        subplot=c_order,
                        linestyle=:dot,
                        kwargs...
                    )
                end
                if c_order==length(orders)
                    std_plot2(;label="truth",legend=false)#:outertopright)
                else
                    std_plot2(;label=false)
                end
            end
            (plt_type==cdf) && push!(L1_cdf_errors,L1_cdf_error_row)
            (plt_type==cdf) && push!(ks_errors,ks_error_row)
            (plt_type==pdf) && push!(L1_pdf_errors,L1_pdf_error_row)
            (plt_type==pdf) && push!(L2_pdf_errors,L2_pdf_error_row)
        end
        #display(p)
        savefig(p,pth*"/figs/"*string(plt_type)*"_func_count_"*string(func_count)*".svg")
    end
    # #display(L1_cdf_errors)
    file = pth*"/data/l1_cdf_"*"func_count_"*string(func_count)
    CSV.write(file*".csv",L1_cdf_errors)
    q = plot()
    for reconstruction in names(L1_cdf_errors)
        plot!(orders,L1_cdf_errors[:,reconstruction],label=reconstruction)
    end
    plot!(xlabel="Order", ylabel="log₁₀ Error",
        title="L¹ error between the true CDF and approximations",
        legend=:right)
    # #display(q)
    file = pth*"/figs/l1_cdf_"*"func_count_"*string(func_count)
    savefig(q,file*".svg")

    # #display(ks_errors)
    file = pth*"/data/ks_cdf_"*"func_count_"*string(func_count)
    CSV.write(file*".csv",ks_errors)
    q = plot()
    for reconstruction in names(ks_errors)
        plot!(orders,ks_errors[:,reconstruction],label=reconstruction)
    end
    plot!(xlabel="Order", ylabel="log₁₀ Error",
        title="KS error between the true CDF and approximations",
        legend=:right)
    # #display(q)
    file = pth*"/figs/ks_cdf_"*"func_count_"*string(func_count)
    savefig(q,file*".svg")

    # #display(L1_pdf_errors)
    file = pth*"/data/l1_pdf_"*"func_count_"*string(func_count)
    CSV.write(file*".csv",L1_pdf_errors)
    q = plot()
    for reconstruction in names(L1_pdf_errors)
        plot!(orders,L1_pdf_errors[:,reconstruction],label=reconstruction)
    end
    plot!(xlabel="Order", ylabel="log₁₀ Error",
        title="L¹ error between the true PDF and approximations",
        legend=:right)
    # #display(q)
    file = pth*"/figs/l1_pdf_"*"func_count_"*string(func_count)
    savefig(q,file*".svg")

    # #display(L2_pdf_errors)
    file = pth*"/data/l2_pdf_"*"func_count_"*string(func_count)
    CSV.write(file*".csv",L2_pdf_errors)
    q = plot()
    for reconstruction in names(L2_pdf_errors)
        plot!(orders,L2_pdf_errors[:,reconstruction],label=reconstruction)
    end
    plot!(xlabel="Order", ylabel="log₁₀ Error",
        title="L² error between the true PDF and approximations",
        legend=:right)
    # #display(q)
    file = pth*"/figs/l2_pdf_"*"func_count_"*string(func_count)
    savefig(q,file*".svg")

end

for (func_count,init_dist_fun) in enumerate(distns)
    pth = mkpath((@__FILE__)[1:end-3]*"/func_count_"*string(func_count)*"/mesh_comp")
    mkpath(pth*"/data")
    mkpath(pth*"/figs")

    L1_cdf_errors = DataFrame(DG = Float64[], Order_1 = Float64[], QBDRAP = Float64[])
    ks_errors = DataFrame(DG = Float64[], Order_1 = Float64[], QBDRAP = Float64[])
    L1_pdf_errors = DataFrame(DG = Float64[], Order_1 = Float64[], QBDRAP = Float64[])
    L2_pdf_errors = DataFrame(DG = Float64[], Order_1 = Float64[], QBDRAP = Float64[])
    for plt_type in (cdf,pdf)
        p = plot(layout=(3,length(orders)),legend=false)
        linetypes = [:solid,:solid,:solid]
        c_order = 0
        for order in orders
            c_order += 1
            c_mesh = 0

            L1_cdf_error_row = [ 0.0 0.0 0.0]
            ks_error_row = [ 0.0 0.0 0.0]
            L1_pdf_error_row = [ 0.0 0.0 0.0]
            L2_pdf_error_row = [ 0.0 0.0 0.0]
            for mtype in (DGMesh,DGMesh,FRAPMesh)
                c_mesh += 1
                tmp_nodes = c_mesh==2 ? (nodes[1]:(hx/order):nodes[end]) : nodes
                tmp_order = c_mesh==2 ? 1 : order
                mesh = mtype(tmp_nodes,tmp_order)

                dq = DiscretisedFluidQueue(model,mesh)

                d0 = init_dist_fun(dq)
                # V = DiscretisedFluidQueues.vandermonde(order)

                rec = (mtype==FRAPMesh) ? plt_type(d0, eval(Symbol(:normalised_closing_operator_,plt_type))) : plt_type(d0)

                ## errors
                (plt_type==cdf) && (L1_cdf_error_row[c_mesh] = log10(DiscretisedFluidQueues.Lp(x->rec(x,1),x->CDFs[func_count](x,1),range(eps(),1-eps(),length=n_evals))))
                (plt_type==cdf) && (ks_error_row[c_mesh] = log10(DiscretisedFluidQueues.kolmogorov_smirnov(x->rec(x,1),x->CDFs[func_count](x,1),range(eps(),1-eps(),length=n_evals))))
                (plt_type==pdf) && (L1_pdf_error_row[c_mesh] = log10(DiscretisedFluidQueues.Lp(x->rec(x,1),x->PDFs[func_count](x,1),range(eps(),1-eps(),length=n_evals))))
                (plt_type==pdf) && (L2_pdf_error_row[c_mesh] = log10(DiscretisedFluidQueues.Lp(x->rec(x,1),x->PDFs[func_count](x,1),range(eps(),1-eps(),length=n_evals),2)))

                ## plotting...
                y_lim_vals = (plt_type==pdf) ? (-3.0,9.0) : (-0.2,1.4)
                y_tick_vals = ((plt_type==pdf) ? (y_lim_vals[1]:3:y_lim_vals[2]-1) : (0:0.5:1))
                std_plot(args...; kwargs...) = plot!(
                    p.layout.grid[c_mesh,c_order],
                    x->rec(x,1),
                    nodes[1]+eps(), nodes[end]-eps();
                    ylim=y_lim_vals,
                    linestyle=linetypes[c_mesh],
                    xticks=false,
                    yticks=false,
                    grid=false,
                    legend=false,
                    kwargs...,
                )
                if c_order == 1
                    std_plot(;label=false,ylabel=string(mtype)[1:findfirst("M",string(mtype))[1]-1],legend=false)
                    yticks!(p.layout.grid[c_mesh,c_order], y_tick_vals)
                elseif c_order == length(orders)
                    std_plot(;legend=false)#:outertopright)
                else
                    std_plot(;label=false,legend=false)
                end
                if c_mesh == 1
                    plot!(p.layout.grid[c_mesh,c_order]; title=string(order),legend=false)
                elseif c_mesh == 3
                    xticks!(p.layout.grid[c_mesh,c_order], nodes)
                    plot!(p.layout.grid[c_mesh,c_order]; xlabel="x",legend=false)
                end
                std_plot2(args... ; kwargs...) = if plt_type==cdf
                    plot!(
                        p.layout.grid[c_mesh,c_order],
                        x->CDFs[func_count](x,1),
                        nodes[1], nodes[end],
                        args...;
                        linestyle=:dot,
                        legend=false,
                        kwargs...,
                    )
                elseif plt_type==pdf
                    plot!(
                        p.layout.grid[c_mesh,c_order],
                        x->PDFs[func_count](x,1),
                        nodes[1], nodes[end],
                        args...;
                        subplot=c_order,
                        linestyle=:dot,
                        legend=false,
                        kwargs...
                    )
                end
                if c_order==length(orders)
                    std_plot2(;legend=false)#:outertopright)
                else
                    std_plot2(;label=false)
                end
            end
            (plt_type==cdf) && push!(L1_cdf_errors,L1_cdf_error_row)
            (plt_type==cdf) && push!(ks_errors,ks_error_row)
            (plt_type==pdf) && push!(L1_pdf_errors,L1_pdf_error_row)
            (plt_type==pdf) && push!(L2_pdf_errors,L2_pdf_error_row)
        end
        #display(p)
        # display(pth)
        savefig(p,pth*"/figs/mesh_"*string(plt_type)*"_func_count_"*string(func_count)*".svg")
    end
    # #display(L1_cdf_errors)
    file = pth*"/data/meshs_l1_cdf_"*"func_count_"*string(func_count)
    CSV.write(file*".csv",L1_cdf_errors)
    q = plot()
    for reconstruction in names(L1_cdf_errors)
        plot!(orders,L1_cdf_errors[:,reconstruction],label=reconstruction)
    end
    plot!(xlabel="Order", ylabel="log₁₀ Error",
        title="L¹ error between the true CDF and approximations",
        legend=:right)
    #display(q)
    # #display(pth)
    file = pth*"/figs/meshs_l1_cdf_"*"func_count_"*string(func_count)
    savefig(q,file*".svg")

    # #display(ks_errors)
    file = pth*"/data/meshs_ks_"*"func_count_"*string(func_count)
    CSV.write(file*".csv",ks_errors)
    q = plot()
    for reconstruction in names(ks_errors)
        plot!(orders,ks_errors[:,reconstruction],label=reconstruction)
    end
    plot!(xlabel="Order", ylabel="log₁₀ Error",
        title="KS error between the true CDF and approximations",
        legend=:right)
    #display(q)
    file = pth*"/figs/meshs_ks_"*"func_count_"*string(func_count)
    savefig(q,file*".svg")

    # #display(L1_pdf_errors)
    file = pth*"/data/meshs_l1_pdf_"*"func_count_"*string(func_count)
    CSV.write(file*".csv",L1_pdf_errors)
    q = plot()
    for reconstruction in names(L1_pdf_errors)
        plot!(orders,L1_pdf_errors[:,reconstruction],label=reconstruction)
    end
    plot!(xlabel="Order", ylabel="log₁₀ Error",
        title="L¹ error between the true PDF and approximations",
        legend=:right)
    #display(q)
    file = pth*"/figs/meshs_l1_pdf_"*"func_count_"*string(func_count)
    savefig(q,file*".svg")

    # #display(L2_pdf_errors)
    file = pth*"/data/meshs_l2_pdf_"*"func_count_"*string(func_count)
    CSV.write(file*".csv",L2_pdf_errors)
    q = plot()
    for reconstruction in names(L2_pdf_errors)
        plot!(orders,L2_pdf_errors[:,reconstruction],label=reconstruction)
    end
    plot!(xlabel="Order", ylabel="log₁₀ Error",
        title="L² error between the true PDF and approximations",
        legend=:right)
    #display(q)
    file = pth*"/figs/meshs_l2_pdf_"*"func_count_"*string(func_count)
    savefig(q,file*".svg")

end


text:
function score(x, y)
    distance_to_center = √(x^2 + y^2)
    if 0 <= distance_to_center <= 1 return 10 end
    if 1 < distance_to_center <= 5 return 5 end
    if 5 < distance_to_center <= 10 return 1 end
    return 0
end

max_stars_repo_name: vladflore/exercism

max_issues_repo_name: vladflore/exercism

max_forks_repo_name: vladflore/exercism


text:
using Revise
using AdFem
using PyPlot
using LinearAlgebra
using Statistics
using MAT
SMALL_SIZE = 20
MEDIUM_SIZE = 20
BIGGER_SIZE = 20

matplotlib.rc("font", size=SMALL_SIZE)          # controls default text sizes
matplotlib.rc("axes", titlesize=SMALL_SIZE)     # fontsize of the axes title
matplotlib.rc("axes", labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
matplotlib.rc("xtick", labelsize=SMALL_SIZE)    # fontsize of the tick labels
matplotlib.rc("ytick", labelsize=SMALL_SIZE)    # fontsize of the tick labels
matplotlib.rc("legend", fontsize=SMALL_SIZE)    # legend fontsize
matplotlib.rc("figure", titlesize=BIGGER_SIZE)  # fontsize of the figure title

function f(x, y)
    0.0001*(10/(1+x^2) + x * y + 10*y^2)
end

Δt = 0.01
mmesh = Mesh(30, 30, 1/30)

xy = gauss_nodes(mmesh)
x, y = xy[:,1], xy[:, 2]
init_nu = 0.001*(1. .+ 7.  .- 4. *x + 11. *y)
E = eval_f_on_gauss_pts(f, mmesh)

E2 = matread("fenics/bwd2-61.mat")["nu"]
E3 = matread("pixel/bwd3-16.mat")["nu"]

close("all")
figure(figsize=(32, 5))
subplot(141)
title("Exact")
visualize_scalar_on_gauss_points(E, mmesh)
gca().set_rasterized(true)
subplot(142)
title("Difference (DNN)")
visualize_scalar_on_gauss_points(abs.(E-E2), mmesh, vmin=0.0, vmax = 2e-4)
xlabel("x")
subplot(143)
title("Difference (Discretization on Gauss Points)")
visualize_scalar_on_gauss_points(abs.(E-E3), mmesh, vmin=0.0, vmax = 2e-4)
xlabel("x")

subplot(144)
title("Initial Guess for the Discretization Method")
visualize_scalar_on_gauss_points(init_nu, mmesh)
xlabel("x")

tight_layout()
savefig("Burgers.png")
savefig("Burgers.pdf")


text:
maybe_counting_sort!(::Any, ::Any) = false

function maybe_counting_sort!(
    v::AbstractVector{<:Integer},
    order::Union{typeof(Base.Forward),typeof(Base.Reverse)},
)
    (rangelen, minval), ok = preprocess_counting_sort(v)
    (ok && rangelen < div(length(v), 2)) || return false

    # For `Vector{Int64}`, using threads starts to make sense when the
    # length is around 2^17.

    if length(v) >= 2^17
        parallel_sort_int_range!(v, rangelen, minval, order === Base.Reverse)
    else
        sort_int_range!(v, rangelen, minval, order === Base.Reverse ? reverse : identity)
    end

    return true
end

function preprocess_counting_sort(v::AbstractVector{<:Integer})
    minval, maxval = ThreadsX.extrema(v)
    (diff, o1) = Base.sub_with_overflow(maxval, minval)
    (rangelen, o2) = Base.add_with_overflow(diff, oneunit(diff))
    return (rangelen, minval, maxval), (!o1 && !o2)
end

function sort_int_range!(
    x::AbstractVector{<:Integer},
    rangelen,
    minval,
    maybereverse = identity,
)
    offs = 1 - minval

    where = fill(0, rangelen)
    @inbounds for i in eachindex(x)
        where[x[i]+offs] += 1
    end

    idx = firstindex(x)
    @inbounds for i in maybereverse(1:rangelen)
        lastidx = idx + where[i] - 1
        val = i - offs
        for j in idx:lastidx
            x[j] = val
        end
        idx = lastidx + 1
    end

    return x
end

function parallel_sort_int_range!(
    x::AbstractVector{<:Integer},
    rangelen,
    minval,
    rev = false;
    nthreads = Threads.nthreads(),
)

    tasks = Task[]
    chunks = _partition(x, cld(length(x), nthreads))
    @sync for xchunk in chunks
        push!(tasks, @spawn count_ints(xchunk, rangelen, minval))
    end
    where = fetch(tasks[1])
    for t in @view tasks[2:end]
        where .+= fetch(t)
    end
    if rev
        reverse!(where)
    end
    pushfirst!(where, 0)
    acc = 0
    for i in 2:length(where)
        acc = @inbounds where[i] += acc
    end

    indices = _partition(1:rangelen, cld(rangelen, nthreads))
    if rev
        indices = reverse(indices)
    end
    @sync for ichunk in indices
        @spawn scatter_ints!(x, minval, where, ichunk)
    end

    return x
end

function count_ints(xs, rangelen, minval)
    offs = 1 - minval
    where = fill(0, rangelen)
    @inbounds for x in xs
        where[x+offs] += 1
    end
    where
end

function scatter_ints!(x, minval, offsets, indices)
    offs = 1 - minval
    @inbounds for i in indices
        val = i - offs
        for j in offsets[i]:offsets[i+1]-1
            x[firstindex(x)+j] = val
        end
    end
end

function parallel_counting_sort!(
    v::AbstractVector{<:Integer},
    order::Union{typeof(Base.Forward),typeof(Base.Reverse)} = Base.Forward,
)
    (rangelen, minval, maxval), ok = preprocess_counting_sort(v)
    ok || error("Value range too large: min = $(repr(minval)), max = $(repr(maxval))")
    parallel_sort_int_range!(v, rangelen, minval, order === Base.Reverse)
    return v
end


text:
# Set the measure to be α*exp(-α^2*x^2)/sqrt(π), for an α > 0.

function getlogγGaussianmercer(β::T, n::Int)::T where T <: Real
    return 0.5*(log(β) -(n-1)*log(2) -lgamma(n))
end

function setupmercer(θ::GaussianKernel1DType, α::T) where T <: Real
    @assert α > zero(T)
    @assert θ.ϵ_sq > zero(T)

    ϵ = sqrt(θ.ϵ_sq)

    β = (1+(2*ϵ/α)^2)^(0.25)
    δ_sq = α^2/2*(β^2-1)

    return ϵ, β, δ_sq
end

# function gethermitepolynomials(x::T, M::Int)::Vector{T} where T <: Real
#     if M == 0
#         return [one(T)]
#     elseif M == 1
#         return [one(T); 2*x]
#     end
#
#     H = Vector{T}(undef, M+1)
#     H[1] = one(T)
#     H[2] = 2*x
#     for i = 2:M
#         m = i-1
#         H[i+1] = 2*x*H[i] -2*m*H[i-1]
#     end
#
#     return H
# end
# # # test code.
# # x = randn()
# # oracle = [1 2*x 4x^2-2 8*x^3-12*x 16*x^4-48*x^2+12]

function evaleigenfunction( θ::GaussianKernel1DType,
                            α::T,
                            β::T,
                            δ_sq::T,
                            n::Int,
                            x::T) where T <: Real

    return exp( getlogγGaussianmercer(β, n) -δ_sq*x^2 ) *GSL.sf_hermite_phys(n-1, α*β*x)
end

# the measure is μ(x) = α/sqrt(π)*exp(-α^2*x^2).
function evaleigenfunctions(θ::GaussianKernel1DType,
                            N::Int,
                            α::T,
                            x::T) where T <: Real

    ϵ, β, δ_sq = setupmercer(θ, α)
    #H = gethermitepolynomials(α*β*x, N-1)

    # eigenfunction evaluations at x.
    ϕ_evals = Vector{T}(undef,N)
    fill!(ϕ_evals,Inf)
    for n = 1:N
        i = n+1

        #ϕ_evals[n] = exp( getlogγGaussianmercer(β, n) -δ_sq*x^2 ) *GSL.sf_hermite_phys(n-1,α*β*x)
        ϕ_evals[n] = evaleigenfunction(θ, α, β, δ_sq, n, x)
    end

    return ϕ_evals
end

function geteigenvalues(    θ::GaussianKernel1DType,
                            N::Int,
                            α::T) where T <: Real

    ϵ, β, δ_sq = setupmercer(θ, α)
    α_sq = α^2
    ϵ_sq = ϵ^2

    # eigenfunction evaluations at x.
    λn = Vector{T}(undef,N)
    fill!(λn,Inf)
    for n = 1:N
        λn[n] = sqrt(α_sq/(α_sq+δ_sq+ϵ_sq))*(ϵ_sq/(α_sq+δ_sq+ϵ_sq))^(n-1)
    end

    return λn
end

# explore SOS
# M is the truncation limit.
function getQmat(c::Vector{T}, ϕ::Function, X, M::Int) where T <: Real
    Q = Matrix{T}(undef,M,M)

    for i = 1:M
        for j = 1:M
            Q[i,j] = sum( c[l]*ϕ(i,X[l])*ϕ(j,X[l]) for l = 1:length(c) )
        end
    end

    return Q
end

# get the Φ matrix.
function getΦmat(ϕ::Function, X::Vector{Vector{T}}, M::Int) where T <: Real
    N = length(X)

    Φ = Matrix{T}(undef, M, N)
    for m = 1:M
        for i = 1:N
            Φ[m,i] = ϕ(m,X[i])
        end
    end

    return Φ
end

# get the Ψ matrix. ψ_b = ϕ_i * ϕ_j, b is col-major ordering.
function getΨmat(ϕ::Function, X::Vector{Vector{T}}, M::Int) where T <: Real
    N = length(X)
    𝑖 = CartesianIndices((M,M))

    Ψ = Matrix{T}(undef, M*M, N)
    fill!(Ψ, Inf)
    for m = 1:M*M
        for l = 1:N

            Ψ[m,l] = ϕ(𝑖[m][1],X[l])*ϕ(𝑖[m][2],X[l])
        end
    end

    return Ψ
end

# alternative to getΨmat. Does the same thing. Less computationally efficient.
function getΨmat2(ϕ::Function, X::Vector{Vector{T}}, M::Int) where T <: Real
    N = length(X)
    Φ = getΦmat(ϕ, X, M)
    B = kron(Φ,Φ)

    Ψ = Matrix{T}(undef, M*M, N)
    for l = 1:N
        Ψ[:,l] = B[:,(l-1)*N+l]
    end

    return Ψ
end

# alternative to getΨmat. Does the same thing. Less computationally efficient.
function getQfromα(ϕ::Function, X::Vector{Vector{T}}, M::Int, α::Vector{T}) where T <: Real

    Q = zeros(T, M, M)
    for i = 1:M
        for j = 1:M
            for l = 1:N
                Q[i,j] += α[l]*ϕ(i,X[l])*ϕ(j,X[l])
            end
        end
    end

    return Q
end

function getQvec(Q::Matrix{T}, M::Int) where T <: Real
    N = length(X)
    𝑖 = CartesianIndices((M,M))

    Q_vec = Vector{T}(undef, M*M)
    for m = 1:M*M
        Q_vec[m] = Q[𝑖[m]]
    end

    return Q_vec
end

function generatecandQSOS(α::Vector{T}, ϕ::Function, X::Vector{Vector{T}}) where T <: Real
    N = length(X)

    Q = Matrix{T}(undef, M, M)
    for i = 1:M
        for j = 1:M
            Q[i,j] = sum( α[l]*ϕ(i,X[l])*ϕ(j,X[l]) for l = 1:N )
        end
    end

    return Q
end

function evalChebyshevpolynomial(n::Int, x::T)::T where T <: Real
    return cos(n*acos(x))
end

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
using BenchmarkTools
using GeoPhyInv
using Test

# # Intro

#md # ```@docs
#md # SeisForwExpt
#md # Base.getindex(::GeoPhyInv.PFdtd, ::Symbol, ::Int)
#md # ```

# # Methods
#md # ```@docs
#md # update!(::GeoPhyInv.PFdtd)
#md # update!(::GeoPhyInv.PFdtd, ::Medium)
#md # update!(::GeoPhyInv.PFdtd, ::SrcWav, ::Any)
#md # ```
#
#
#
#
#
#
#


text:
# licor_split.jl
#
#   Split Li-7200 TXT files into GHG archives
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 0.7
# 20.05.2014
# Last Edit: 12.04.2019

"# licor_split(source::String,destination::String;verbose::Bool=false)

Split an Li-7200 text file into small files in the standard GHG format. If possible, files will start on the usual four hour marks (0000, 0400, 0800, etc.).

`licor_split(source,destination)` Split Li7200 files in the source directors\n
* **source**::String = Directory of Li-7200 text files with names similar to 2016-03-17T104500.txt
* **destination**::String = Destination directory\n\n

---\n

#### Keywords:\n
* verbose::Bool = Display information as the function runs, TRUE is default
\n\n"
function licor_split(Dr::String,Dest::String;verbose::Bool=true)
    ###################
    ##  Definitions  ##
    ###################
    # Dr = Source Directory
    # Dest = Destination Directory

    ####################
    ##  Parse Inputs  ##
    ####################
    # Check Directory Validity
    if ~isdir(Dr) | ~isdir(Dest)
        error("Source or destination directory does not exist")
    end

    (files,folders) = dirlist(Dr,regex=r"\d{4}-\d{2}-\d{2}T\d{6}\.txt$",recur=1) # List TXT files

    #####################
    ##  Process Files  ##
    #####################
    if isempty(files)
	    println("No Metdata files to split")
    else
    	if verbose
	    	println("=== Splitting Metdata Files ===")
	    end
    end
    time = DateTime[]
    starttime = DateTime[]
    endtime = DateTime[]
    fname = String[]
    header = Array{String}(undef,0) # Initialized header array
    sentinel = 20; # Sentinel value for header while loop
    for i=1:1:length(files)
    	if verbose
	    	println("   " * files[i])
	    end
	    s = stat(files[i]) # Get file info

	    fid = open(files[i],"r")

	    # Load Header
	    j = 1;
	    l = "Bierchen" # Needs any string at least 5 characters long
	    while l[1:5] != "DATAH" && j < sentinel
	    	l = readline(fid,keep=true)
	    	header = [header;l]
	    	j += 1
	    end
	    data_pos = position(fid) # Position of where the data starts

	    #Get Start Time
	    l = readline(fid,keep=true)
		ft = findall(x -> x == '\t',l)
	    starttime = DateTime(l[ft[6]+1:ft[8]-1],"yyyy-mm-dd\tHH:MM:SS:sss")

	    # Get End Time
	    endline = true
	    lastpos = 4 # Characters from the end of the file where the last line begins
	    while endline == true
		    seek(fid,s.size-lastpos)
		    l = readline(fid,keep=true)
		    if eof(fid) == false
			    endline = false
		    else
			    lastpos += 2
		    end
	    end
	    l = readline(fid,keep=true)
		ft = findall(x -> x == '\t',l)
	    endtime = DateTime(l[ft[6]+1:ft[8]-1],"yyyy-mm-dd\tHH:MM:SS:sss")

	    # Split the File
	    seek(fid,data_pos) # Move the cursor to the start of the data
	    fid2 = open(files[i]);close(fid2)
	    next_start = []
	    first_file = true
	    while eof(fid) == false
		    if isopen(fid2) == false
			    # No output stream available
			    # Open new output file from previous line

			    if first_file == true
				    # There is not previously loaded line, load one now
				    l = readline(fid,keep=true)
				    first_file = false
			    end

			    # Calculate last value of file
				ft = findall(x -> x == '\t',l)
			    temp_start = DateTime(l[ft[6]+1:ft[8]-1],"yyyy-mm-dd\tHH:MM:SS:sss")
			    HH2 = Dates.Hour[]
			    DD2 = Dates.Day(temp_start)
			    if Dates.Hour(0) <= Dates.Hour(temp_start) < Dates.Hour(4)
				    HH2 = Dates.Hour(4)
			    elseif Dates.Hour(4) <= Dates.Hour(temp_start) < Dates.Hour(8)
				    HH2 = Dates.Hour(8)
			    elseif Dates.Hour(8) <= Dates.Hour(temp_start) < Dates.Hour(12)
				    HH2 = Dates.Hour(12)
			    elseif Dates.Hour(12) <= Dates.Hour(temp_start) < Dates.Hour(16)
				    HH2 = Dates.Hour(16)
			    elseif Dates.Hour(16) <= Dates.Hour(temp_start) < Dates.Hour(20)
				    HH2 = Dates.Hour(20)
			    elseif Dates.Hour(20) <= Dates.Hour(temp_start) < Dates.Hour(24)
				    HH2 = Dates.Hour(0)
				    DD2 = DD2 + Dates.Day(1)
			    end
			    next_start = DateTime(Dates.Year(temp_start),Dates.Month(temp_start),DD2,HH2)

			    # Generate File Name
			    for j=1:1:length(header)
			    	if occursin(r"^Instrument:",header[j])
			    		fname = joinpath(Dest,Dates.format(temp_start,"yyyy-mm-ddTHHMMSS") * "_" * header[j][13:end-1] * ".data")
			    	end
				end
				if verbose
			    	println("\t\t- " * fname)
			    end
			    fid2 = open(fname,"w+")

			    # Find and replace the Timestamp header line
			    for j=1:1:length(header);
			    	if occursin(r"^Timestamp",header[j])
			    		header[j] = "Timestamp:  " * Dates.format(temp_start,"yyyy-mm-dd HH:MM:SS:sss") * "\n"
			    	end
			    end

			    # Iteratively write the header
			    for j=1:1:length(header)
			    	write(fid2,header[j])
			    end

			    write(fid2,l) # Write the line used to create the file name
		    else
			    # An output stream is available
			    l = readline(fid,keep=true) # Load in another line

			    # Write or Close
				ft = findall(x -> x == '\t',l)
			    temp_end = DateTime(l[ft[6]+1:ft[8]-1],"yyyy-mm-dd\tHH:MM:SS:sss")
			    if temp_end >= next_start
				    # The current line is newer than the start of the next file, close the current file
				    close(fid2)

				    # Zipped File Name (minus extension)
				    fzip = splitext(fname)[1]

				    # Zip File
				    if verbose
				    	println("\t\t\tCompressing...")
				    end
				    fzip = fzip * ".ghg"
					ziptextfiles(fzip,fname) # Zip file
				    rm(fname)
			    else
				    # Still within range, write the current line
				    write(fid2,l)
			    end
		    end
	    end

	    close(fid2)
	    close(fid)

		# Zipped File Name (minus extension)
		fzip = splitext(fname)[1]

		# Zip the final file
		if verbose
			println("\t\t\tCompressing...")
		end
		fzip = fzip * ".ghg"
	    fzip2 = fzip * ".ghg"
		ziptextfiles(fzip,fname) # Zip file
	    rm(fname)
    end

    if verbose
    	println("Complete")
    end
end

max_stars_repo_name: JeremyRueffer/ClimateDataIO

max_issues_repo_name: JeremyRueffer/ClimateDataIO