text:
function _notnan!(res, col)
    @inbounds for (i, el) in enumerate(col)
        if typeof(el) <:Real
            res[i] &= !isnan(el)
        end
    end
    return nothing
end

function complete_vals(df::AbstractDataFrame, col::Colon=:)
    if ncol(df) == 0
        throw(ArgumentError("Unable to compute complete vals of a data frame with no columns"))
    end
    res = trues(size(df, 1))
    for i in 1:size(df, 2)
        _notnan!(res, df[!, i])
    end
    res
end

function complete_vals(df::AbstractDataFrame, col::DataFrames.ColumnIndex)
    res = trues(size(df, 1))
    _notnan!(res, df[!, col])
    res
end

complete_vals(df::AbstractDataFrame, cols::Union{AbstractVector, Regex, Not, Between, All}) =
    complete_vals(df[!, cols])

function dropnan(df::AbstractDataFrame, cols =:)
    newdf = df[complete_vals(df, cols), :]
    # disallowmissing && disallowmissing!(newdf, cols)
    newdf
end

function dropnan!(df::AbstractDataFrame,
                      cols=:)
    delete!(df, (!).(complete_vals(df, cols)))
    df
end

function Protocol_names!(df)
    df[!,:env] = map(x -> begin
                        if x == "0.5"
                            return "Low"
                        elseif x == "0.75"
                            return "Med."
                        elseif x == "1.0"
                            return "High"
                        else
                            return "weird"
                        end
                    end,
                    df.Protocol)
    return df
end

function Protocol_colors!(df)
    Protocol_names!(df)
    df[!,:color] = [get(protocol_colors,x,RGB(0,0,0)) for x in df.env]
    return df
end

function Drug_colors!(df)
    if :Treatment in propertynames(df)
        c = :Treatment
    elseif :Phase in propertynames(df)
        c = :Phase
    else
        error("didn't find any color to map")
    end
    df[!,:color] = [get(drug_colors,x,RGB(0,0,0)) for x in df[:,c]]
    return df
end

## calculate the extrema to filter the data out of the 95% of the distribution
conf95(v) = quantile(v,[0.25,0.975])

function conf_ints(vec; mode = :extrema, percent = 95)
    v = jump_missing(vec)
    p = (100 - percent) / 100
    if mode == :extrema
        p_half = p/2
        lower, upper = round.(quantile(v,[p_half,1-p_half]), digits = 3)
    elseif mode == :left
        lower = round.(quantile(v,p), digits = 3)
        upper = round(maximum(v), digits = 3) +1
    elseif mode == :right
        lower = round(minimum(v), digits = 3) -1
        upper = round.(quantile(v,1-p), digits = 3)
    else
        error("unknown confidence intervals mode")
    end
    return [lower,upper]
end

function trim_conf_ints(v; mode = :extrema, percent = 95)
    lower, upper = conf_ints(v; mode = mode, percent = percent)
    [ismissing(x) ? true : lower <= x <= upper for x in v]
end

function trim_conf_ints(df::AbstractDataFrame,x::Symbol; mode = :extrema, percent = 95)
    v = df[:,x]
    idxs = trim_conf_ints(v; mode = mode, percent = percent)
    df[idxs,:]
end

function trim_conf_ints!(df::AbstractDataFrame,x::Symbol; mode = :extrema, percent = 95)
    v = df[:,x]
    idxs = trim_conf_ints(v; mode = mode, percent = percent)
    delete!(df,Not(idxs))
end

function jump_missing(v::AbstractVector)
    res = v[map(!,ismissing.(v))]
    disallowmissing(res)
end

function jump_NaN(v::AbstractVector)
    res = v[map(!,isnan.(v))]
end

skipmean(v) = mean(skipmissing(v))

function binquantile(v, bins = 15)
    qs = nquantile(v,bins)
    nq = [findfirst(q -> o < q, qs) for o in v]
    [isnothing(q) ? 15 : q - 1 for q in nq]
end

"""
    calculate the elapsed time from the last poke of the previous trial assuming it receives a pokes dataframe
"""
function TimeFromLeaving(dfpokes)
    start = dfpokes[1,:PokeIn]
    E = Float64[]
    for (i,o,l) in zip(dfpokes.PokeIn, dfpokes.PokeOut, dfpokes.LastPoke)
        ela = i - start
        push!(E,ela)
        if l
            start = o
        end
    end
    return E
end


text:
# Autogenerated wrapper script for Minuit2_jll for x86_64-apple-darwin
export libMinuit2, libMinuit2Math

JLLWrappers.@generate_wrapper_header("Minuit2")
JLLWrappers.@declare_library_product(libMinuit2, "@rpath/libMinuit2.dylib")
JLLWrappers.@declare_library_product(libMinuit2Math, "@rpath/libMinuit2Math.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libMinuit2,
        "lib/libMinuit2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libMinuit2Math,
        "lib/libMinuit2Math.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
using LinearAlgebra
using SparseArrays
using Plots
using ForwardDiff
using FastGaussQuadrature # more accurate quadrature

"This routine solves Poisson's equation -u''(x) = f(x) in weak form using piecewise
linear finite element methods. It's used to illustrate convergence in different norms."

m = 200 # number of elements

# Manufactured solution
uexact(x) = log(2+sin(4*pi*x))
dudx_exact(x) = ForwardDiff.derivative(uexact,x)
f(x) = -ForwardDiff.derivative(dudx_exact,x) # -d/dx (κ(x)*du/dx) = f
α = uexact(-1)
β = uexact(1)

# define spatial grid
x = LinRange(-1,1,m+2) # x_0, x_1, ..., x_m, x_{m+1} = x_0
# x = @. x + randn() / (2*m) # randomly perturb nodes
# x[1] = -1; x[end] = 1 # reset endpoints

# define local FEM basis
rq,wq = gausslegendre(100) # define overkill-accurate Gauss quadrature
λ(r) = [(1 .-r)./2 (1 .+r)./2]
#dλ(r) = [-.5*ones(length(r)) .5*ones(length(r))]
dλ(r) = [-.5 .5]
map_point(x,a,b) = a + (b-a) * (1+x)/2 # maps x ∈ [-1,1] to interval [a,b]

"∫u' * ϕ_i' = ∫f(x)*ϕ_i(x)"
A = spzeros(m+2,m+2)
A_local = [1 -1;-1 1]
b = zeros(m+2)
for e = 1:m+1 # loop over the intervals
    ids = e:e+1
    h_e = x[e+1]-x[e]

    x_midpt = .5*(x[e+1]+x[e])
    A[ids,ids] .+= A_local * 1/h_e

    # # one point midpoint quadrature
    # b[ids] .+= h_e / 2.0 * f(.5*(x[e]+x[e+1])) # f(.5*(x_i+x_i+1))

    # higher accuracy Gauss quadrature
    xq = map_point.(rq,x[e],x[e+1])
    b[ids] .+= h_e / 2.0 * λ(rq)'*(wq.*f.(xq)) # ∑_k f(x_k) * w_k * ϕ_i(x_k)
end

# impose Dirichlet BCs on both sides
b̃ = b - vec(A[:,1]*α)
b̃[1] = α
A[1,:] .= 0
A[:,1] .= 0
A[1,1] = 1

b̃ = b̃ - vec(A[:,m+2]*β)
b̃[m+2] = β
A[m+2,:] .= 0
A[:,m+2] .= 0
A[m+2,m+2] = 1

# convert b̃ from a SparseVector to a Vector
# (Julia doesn't currently allow "\" to be applied to SparseVectors)
u = A\Vector(b̃)

nodal_err = abs.(u - uexact.(x))
plot(x,u,mark=:dot,ms=3,legend=false,title="Max error at nodes = $(maximum(nodal_err))")
xfine = LinRange(-1,1,1000)
plot!(xfine,uexact.(xfine),legend=false)

# compute H1 and L2 errors
function compute_error(u,m,x,rq,wq)
    H1_err2 = 0.0
    L2_err2 = 0.0
    for e = 1:m+1
        ids = e:e+1
        h_e = x[e+1]-x[e]
        rx = 2 / h_e
        xq = map_point.(rq,x[e],x[e+1])
        dudx_err = rx * dλ(rq)*u[ids] .- dudx_exact.(xq)
        H1_err2 += h_e / 2 * dot(wq, dudx_err.^2)

        u_err = λ(rq)*u[ids] - uexact.(xq) # ∫(u-uexact)^2
        L2_err2 += h_e / 2 * dot(wq, u_err.^2)
    end
    return sqrt(H1_err2), sqrt(L2_err2)
end
H1_err,L2_err = compute_error(u,m,x,rq,wq)
@show H1_err # will be proportional to O(h)
@show L2_err


text:
######################################################################
# Example of how to create a window in borderless fullscreen mode.
# -----
# Licensed under MIT License
using Pkg
Pkg.activate(".")

using KirUtil
using GLFWAbstraction

let vidmode = videomode(monitor(1))
    # Create window in borderless fullscreen mode
    @windowhint RedBits     vidmode.bits.red
    @windowhint GreenBits   vidmode.bits.green
    @windowhint BlueBits    vidmode.bits.blue
    @windowhint RefreshRate vidmode.refresh_rate

    let wnd = window(:mywnd, "Hello, world!", vidmode.width, vidmode.height, monitor(1))
        @assert !wnd.shouldclose
        use(wnd)

        while !wnd.shouldclose
            sleep(0.1)
            swapbuffers(wnd)
            pollevents()
        end
    end
end


text:
using DecFP, Test, Printf, Random, Base.MathConstants, SpecialFunctions

@test DecFP.flags[Threads.threadid()] == 0

import DecFP.isnanstr
@test isnanstr("nan") && isnanstr("  +NAN") && isnanstr("-NaN") && !isnanstr("nano")

function hist(X, n)
    v = zeros(Int, n)
    for x in X
        v[floor(Int, x*n) + 1] += 1
    end
    v
end

function testthreads(T, i, mode)
    @test @sprintf("%.0f", T(i)) == string(i)
    @test rounding(T) == RoundNearest
    setrounding(T, mode) do
        @test rounding(T) == mode
        if mode in (RoundNearest, RoundToZero, RoundDown)
            @test T(1) + eps(T(1)) / 2 == T(1)
        else
            @test T(1) + eps(T(1)) / 2 == T(1) + eps(T(1))
        end
    end
    @test rounding(T) == RoundNearest
    @test_throws InexactError convert(Int64, T("1.5"))
    @test_throws DomainError sqrt(T(-2))
end

for T in (Dec32, Dec64, Dec128)
    @info "TESTING $T    nthreads = $(Threads.nthreads()) ..."

    if T == Dec32
        @test d32"3.2" * d32"4.5" == d32"14.4"
        @test T(1, Int32(maxintfloat(T) - 1), 90) == floatmax(T)
        @test_throws InexactError T(1, Int32(maxintfloat(T)), 90)
        @test_throws InexactError T(1, Int32(maxintfloat(T)) + 1, 0)
        @test T(1, 1, 96) == T("1e96")
        @test T(1, 1, -101) == nextfloat(T(0))
        @test T(1, 10, -102) == nextfloat(T(0))
        @test_throws InexactError T(1, 1, -102)
        @test_throws InexactError T(1, 11, -102)
        @test eps(Dec32) == parse(Dec32, "1e-6")
        @test floatmax(Dec32) == parse(Dec32, "+9999999E+90")
        @test bswap(Dec32(1.5)) == reinterpret(Dec32, bswap(reinterpret(UInt32, Dec32(1.5))))
    elseif T == Dec64
        @test d"3.2" * d"4.5" == d"14.4"
        @test d64"3.2" * d64"4.5" == d64"14.4"
        @test T(1, Int64(maxintfloat(T) - 1), 369) == floatmax(T)
        @test_throws InexactError T(1, Int64(maxintfloat(T)), 369)
        @test_throws InexactError T(1, Int64(maxintfloat(T)) + 1, 0)
        @test T(1, 1, 384) == T("1e384")
        @test T(1, 1, -398) == nextfloat(T(0))
        @test T(1, 10, -399) == nextfloat(T(0))
        @test_throws InexactError T(1, 1, -399)
        @test_throws InexactError T(1, 11, -399)
        @test eps(Dec64) == parse(Dec64, "1e-15")
        @test floatmax(Dec64) == parse(Dec64, "+9999999999999999E+369")
        @test bswap(Dec64(1.5)) == reinterpret(Dec64, bswap(reinterpret(UInt64, Dec64(1.5))))
    else
        @test d128"3.2" * d128"4.5" == d128"14.4"
        @test T(1, Int128(maxintfloat(T)) - 1, 6111) == floatmax(T)
        @test_throws InexactError T(1, Int128(maxintfloat(T)), 6111)
        @test_throws InexactError T(1, Int128(maxintfloat(T)) + 1, 0)
        @test T(1, 1, 6144) == T("1e6144")
        @test T(1, 1, -6176) == nextfloat(T(0))
        @test T(1, 10, -6177) == nextfloat(T(0))
        @test_throws InexactError T(1, 1, -6177)
        @test_throws InexactError T(1, 11, -6177)
        @test eps(Dec128) == parse(Dec128, "1e-33")
        @test floatmax(Dec128) == parse(Dec128, "+9999999999999999999999999999999999E+6111")
        @test bswap(Dec128(1.5)) == reinterpret(Dec128, bswap(reinterpret(UInt128, Dec128(1.5))))
    end

    @test parse(T, "1.5")::T == 1.5
    @test isnan(parse(T, "NaN")::T) && isnan(T(NaN))
    @test isinf(parse(T, "inf")::T) && !isfinite(parse(T, "inf")::T)
    @test parse(T,"inf")::T == typemax(T)::T
    @test parse(T,"-inf")::T == typemin(T)::T
    @test_throws ArgumentError parse(T, "1.0.1")
    @test_throws ArgumentError parse(T, "NaNny")
    @test tryparse(T, "1.5")::T == 1.5
    @test isnan(tryparse(T, "NaN"))
    @test tryparse(T, "NaNny") === nothing
    @test tryparse(T, "1.0.1") === nothing

    @test parse(T, "0.1")::T == T(1//10)
    @test T("0.1")::T == T(1//10)

    if T != Dec128
        @test T(0.1) == T("0.1")
        @test T(0.1, RoundNearest) == T("0.1")
        @test T(0.1, RoundDown) == T("0.1")
        @test T(0.1, RoundUp) == nextfloat(T("0.1"))
    end
    @test T("1", RoundDown) == T(1)
    @test T("1", RoundUp) == T(1)
    @test T("1.0000000000000000000000000000000000000001") == T(1)
    @test T("1.0000000000000000000000000000000000000001", RoundNearest) == T(1)
    @test T("1.0000000000000000000000000000000000000001", RoundDown) == T(1)
    @test T("1.0000000000000000000000000000000000000001", RoundUp) == nextfloat(T(1))
    @test T("0.9999999999999999999999999999999999999999") == T(1)
    @test T("0.9999999999999999999999999999999999999999", RoundNearest) == T(1)
    @test T("0.9999999999999999999999999999999999999999", RoundUp) == T(1)
    @test T("0.9999999999999999999999999999999999999999", RoundDown) == prevfloat(T(1))

    @test T(0, 0, 0) == T(0)
    @test T(1, 0, 10) == T(0)
    @test isequal(T(-1, 0, 0), T(-0.0))
    @test T(1, 125, -2) == T(1.25)
    @test T(-1, 500, -2) == T(-5)
    @test T(1, 5, 2) == T(500)
    @test_throws DomainError T(0, 1, 0)
    @test_throws DomainError T(2, 1, 0)
    @test T(125, -2) == T(1.25)
    @test T(-125, -2) == T(-1.25)
    @test T(1, Int128(maxintfloat(T)), 0) == maxintfloat(T)
    @test T(1, Int128(maxintfloat(T)) * 100, -2) == maxintfloat(T)

    io = IOBuffer()
    show(io, T("NaN")); @test String(take!(io)) == "NaN"
    show(io, T("Inf")); @test String(take!(io)) == "Inf"
    show(io, T("-Inf")); @test String(take!(io)) == "-Inf"
    show(io, T("0")); @test String(take!(io)) == "0.0"
    show(io, T("-0")); @test String(take!(io)) == "-0.0"
    show(io, T("1")); @test String(take!(io)) == "1.0"
    show(io, T("-1")); @test String(take!(io)) == "-1.0"
    show(io, T("1.000")); @test String(take!(io)) == "1.0"
    show(io, T("1e5")); @test String(take!(io)) == "100000.0"
    show(io, T("1e6")); @test String(take!(io)) == "1.0e6"
    show(io, T("1.23456e6")); @test String(take!(io)) == "1.23456e6"
    show(io, T("1e-1")); @test String(take!(io)) == "0.1"
    show(io, T("1e-4")); @test String(take!(io)) == "0.0001"
    show(io, T("1e-5")); @test String(take!(io)) == "1.0e-5"
    show(io, T("1.20e3")); @test String(take!(io)) == "1200.0"
    show(io, T("123.456")); @test String(take!(io)) == "123.456"
    show(io, T("0.00123456")); @test String(take!(io)) == "0.00123456"

    @test @sprintf("%7.2f", T("1.2345")) == "   1.23"
    @test @sprintf("%-7.2f", T("1.2345")) == "1.23   "
    @test @sprintf("%07.2f", T("1.2345")) == "0001.23"
    @test @sprintf("%.0f", T("1.2345")) == "1"
    @test @sprintf("%#.0f", T("1.2345")) == "1."
    @test @sprintf("%.4e", T("1.2345")) == "1.2345e+00"
    @test @sprintf("%.4E", T("1.2345")) == "1.2345E+00"

    @test @sprintf("%6.2f", T("9.999")) == " 10.00"
    @test @sprintf("%9.2e", T("9.999e5")) == " 1.00e+06"

    @test @sprintf("%f", T("Inf")) == "Inf"
    @test @sprintf("%f", T("NaN")) == "NaN"

    @test @sprintf("%.0e", T("3e42")) == "3e+42"
    @test @sprintf("%#.0e", T("3e42")) == "3.e+42"

    @test @sprintf("%e", T("3e42")) == "3.000000e+42"
    @test @sprintf("%E", T("3e42")) == "3.000000E+42"
    @test @sprintf("%e", T("3e-42")) == "3.000000e-42"
    @test @sprintf("%E", T("3e-42")) == "3.000000E-42"

    @sprintf("%.6g", T("12345670.")) == "1.23457e+07"
    @sprintf("%.6g", T("1234567.")) == "1.23457e+06"
    @sprintf("%.6g", T("123456.7")) == "123457"
    @sprintf("%.6g", T("12345.67")) == "12345.7"
    @sprintf("%.6g", T("12340000.0")) == "1.234e+07"

    @test @sprintf("%10.5g", T("123.4")) == "     123.4"
    @test @sprintf("%+10.5g", T("123.4")) == "    +123.4"
    @test @sprintf("% 10.5g", T("123.4")) == "     123.4"
    @test @sprintf("%#10.5g", T("123.4")) == "    123.40"
    @test @sprintf("%-10.5g", T("123.4")) == "123.4     "
    @test @sprintf("%-+10.5g", T("123.4")) == "+123.4    "
    @test @sprintf("%010.5g", T("123.4")) == "00000123.4"
    @test @sprintf("%10.5g", T("-123.4")) == "    -123.4"
    @test @sprintf("%010.5g", T("-123.4")) == "-0000123.4"
    @test @sprintf("%.6g", T("12340000.0")) == "1.234e+07"
    @test @sprintf("%#.6g", T("12340000.0")) == "1.23400e+07"

    @test @sprintf("%.2f %.4f", T("12.34567"), 9.87) == "12.35 9.8700"

    x,y,z = 1.5, -3.25, 0.0625 # exactly represented in binary
    xd = T(x); yd = T(y); zd = T(z)

    @test_throws DomainError sigexp(T(NaN))
    @test_throws DomainError sigexp(T(Inf))
    @test_throws DomainError sigexp(T(-Inf))
    @test sigexp(T(0)) == (1, 0, 0)
    @test sigexp(T("-0")) == (-1, 0, 0)
    @test sigexp(T(1)) == (1, 1, 0)
    @test sigexp(T(-1)) == (-1, 1, 0)
    @test sigexp(T(-1.25)) == (-1, 125, -2)
    @test sigexp(maxintfloat(T) - 1) == (1, Int128(maxintfloat(T) - 1), 0)

    @test fma(xd,yd,zd)::T == muladd(xd,yd,zd)::T == xd*yd+zd == fma(x,y,z)

    @test one(T)::T == 1
    @test zero(T)::T == 0
    @test signbit(xd) == signbit(x) && signbit(yd) == signbit(y) && signbit(parse(T,"-Inf"))
    @test isnan(sign(T(NaN)))
    @test isequal(sign(T(0.0)), T(0.0))
    @test isequal(sign(T(-0.0)), T(-0.0))
    @test sign(xd)::T == sign(x) && sign(yd)::T == sign(y)

    @test nextfloat(xd)::T == xd + eps(T)
    @test prevfloat(xd)::T == xd - eps(T)
    @test nextfloat(T(1.5e10)) == 1.5e10 + eps(T(1.5e10))
    @test prevfloat(T(1.5e10)) == 1.5e10 - eps(T(1.5e10))

    @test floatmin(T) > 0
    @test !issubnormal(floatmin(T))
    @test issubnormal(prevfloat(floatmin(T)))
    @test floatmin(T) == nextfloat(zero(T)) / eps(T)

    @test eps(maxintfloat(T) - 1) == 1
    @test eps(maxintfloat(T)) == 10
    @test maxintfloat(T) == maxintfloat(T)+1 > maxintfloat(T)-1 > maxintfloat(T)-2

    for f in (isnan,isinf,isfinite,issubnormal,abs)
        @test f(xd) == f(x)
    end

    for f in (+,-,*,copysign,==,>,>=,<,<=)
        @test f(xd,yd) == f(x,y)
    end

    for f in (/,hypot,atan,^)
        @test f(xd,yd) ≈ f(x,y)
        if f != ^
            @test f(yd,xd) ≈ f(y,x)
            @test f(yd,yd) ≈ f(y,y)
        end
    end

    for f in (exp,log,sin,cos,tan,asin,acos,atan,sinh,cosh,tanh,asinh,acosh,atanh,log1p,expm1,log10,log2,exp2,exp10,sqrt,cbrt)
        v = f == acosh ? x : z
        @test f(T(v)) ≈ f(v)
    end

    # issue #47
    @test exp10(T(0)) == T(1)

    for x = -3.0:0.25:3.0
        @test all(logabsgamma(T(x)) .≈ logabsgamma(x))
    end
    for x in (NaN, -Inf, -0.0, Inf)
        @test isequal(logabsgamma(T(x)), logabsgamma(x))
    end

    @test_throws DomainError gamma(T(-1))
    @test isnan(gamma(T(NaN)))
    for x in (-1.5, -0.5, -0.0, 0.0, 0.5, 1.0, 1.5, 2.0, Inf)
        @test gamma(T(x)) ≈ gamma(x)
    end

    Random.seed!(1234)
    @test rand(T) != rand(T)
    @test 0 <= rand(T) < 1
    r = rand(convert(T, 97):convert(T, 122))
    @test typeof(r) == T
    @test 97 <= r <= 122
    r = rand(convert(T, 97):convert(T, 2):convert(T, 122), 2)[1]
    @test typeof(r) == T
    @test 97 <= r <= 122
    @test mod(r, 2)==1
    A = Vector{T}(undef, 16)
    rand!(A)
    @test all(0 .<= A .< 1)
    # test uniform distribution
    # array version
    counts = hist(rand(T, 2000), 4)
    @test minimum(counts) > 300 # should fail with proba < 1e-26
    # scalar version
    counts = hist([rand(T) for i in 1:2000], 4)
    @test minimum(counts) > 300

    for c in (π, e, γ, catalan, φ)
        @test T(c) ≈ Float64(c)
    end

    for Tf in (Float64, Float32, Float16)
        @test xd == Tf(x) == T(Tf(x)) == Tf(xd)
    end

    # issue #92
    p = -2
    @test T(2)^-2 == parse(T, "0.25") == T(2)^p

    # exercise literal_pow
    @test T(2)^0 === T(1)
    @test T(2)^1 === T(2)
    @test T(2)^2 === T(4)
    @test T(2)^3 === T(8)

    @test T("0.1") == 1//10
    @test 1//10 == T("0.1")
    @test T("-0.1") == -1//10
    @test T("0.2") == 1//5
    @test T("0.5") == 1//2
    @test T(1) == 1//1
    @test T(0) == 0//1
    @test T(Inf) == 1//0
    @test T(-Inf) == -1//0
    @test T(1) / T(3) != 1//3
    @test T(1) / T(3) < 1//3
    @test T(1) / T(3) <= 1//3
    @test T(2) / T(3) > 2//3
    @test T(2) / T(3) >= 2//3
    @test T(1) / T(6) != 1//6
    @test T(7) / T(2) == 7//2
    @test T(7) / T(5) == 7//5
    @test T(7) / T(100) == 7//100
    @test T(7) / T(300) != 7//300

    for Tf in (Float64, Float32, Float16)
        @test parse(T, "0.1") != parse(Tf, "0.1")
        @test parse(Tf, "0.1") != parse(T, "0.1")
        @test parse(T, "0.7") != parse(Tf, "0.7")
        @test parse(Tf, "0.7") != parse(T, "0.7")
        if Tf == Float16
            @test parse(T, "0.1") > parse(Tf, "0.1")
            @test parse(T, "0.1") >= parse(Tf, "0.1")
            @test parse(Tf, "0.1") < parse(T, "0.1")
            @test parse(Tf, "0.1") <= parse(T, "0.1")
            @test parse(T, "0.7") < parse(Tf, "0.7")
            @test parse(T, "0.7") <= parse(Tf, "0.7")
            @test parse(Tf, "0.7") > parse(T, "0.7")
            @test parse(Tf, "0.7") >= parse(T, "0.7")
        else
            @test parse(T, "0.1") < parse(Tf, "0.1")
            @test parse(T, "0.1") <= parse(Tf, "0.1")
            @test parse(Tf, "0.1") > parse(T, "0.1")
            @test parse(Tf, "0.1") >= parse(T, "0.1")
            @test parse(T, "0.7") > parse(Tf, "0.7")
            @test parse(T, "0.7") >= parse(Tf, "0.7")
            @test parse(Tf, "0.7") < parse(T, "0.7")
            @test parse(Tf, "0.7") <= parse(T, "0.7")
        end
        @test parse(T, "0.5") == parse(Tf, "0.5")
        @test parse(T, "0.5") <= parse(Tf, "0.5")
        @test parse(T, "0.5") >= parse(Tf, "0.5")
        @test parse(Tf, "0.5") == parse(T, "0.5")
        @test parse(Tf, "0.5") <= parse(T, "0.5")
        @test parse(Tf, "0.5") >= parse(T, "0.5")
    end

    @test hash(T(NaN)) == hash(NaN)
    @test hash(T(Inf)) == hash(Inf)
    @test hash(T(-Inf)) == hash(-Inf)
    @test hash(T(0)) == hash(0.0) == hash(0)
    @test hash(T(-0.0)) == hash(-0.0)
    @test hash(T(1)) == hash(1.0)
    @test hash(T(-1)) == hash(-1.0)
    @test hash(T(25)) == hash(25.0)
    @test hash(T(-25)) == hash(-25.0)
    @test hash(T(77)) == hash(77.0)
    @test hash(T(-77)) == hash(-77.0)
    @test hash(T(3500)) == hash(3500.0)
    @test hash(T(-3500)) == hash(-3500.0)
    @test hash(T(1234.5)) == hash(1234.5)
    @test hash(T(-1234.5)) == hash(-1234.5)
    @test hash(T(123.0625)) == hash(123.0625)
    @test hash(T(-123.0625)) == hash(-123.0625)
    @test hash(T(0.0009765625)) == hash(0.0009765625)
    @test hash(T(-0.0009765625)) == hash(-0.0009765625)
    @test hash(T(10)^70) == hash(BigInt(10)^70)

    for x = -5.0:0.25:5.0, y = -5.0:0.25:5.0
        @test isequal(rem(T(x), T(y)), rem(x, y))
        @test isequal(rem(T(x), T(y), RoundNearest), rem(x, y, RoundNearest))
        @test isequal(rem(T(x), T(y), RoundToZero), rem(x, y, RoundToZero))
        @test isequal(rem(T(x), T(y), RoundDown), rem(x, y, RoundDown))
        @test isequal(rem(T(x), T(y), RoundUp), rem(x, y, RoundUp))
        @test isequal(mod(T(x), T(y)), mod(x, y))
        @test isequal(T(x) % T(y), x % y)
    end
    for x in (NaN, -Inf, -0.0, 0.0, 1.0, Inf), y in (NaN, -Inf, -0.0, 0.0, 1.0, Inf)
        @test isequal(rem(T(x), T(y)), rem(x, y))
        @test isequal(rem(T(x), T(y), RoundNearest), rem(x, y, RoundNearest))
        @test isequal(rem(T(x), T(y), RoundToZero), rem(x, y, RoundToZero))
        @test isequal(rem(T(x), T(y), RoundDown), rem(x, y, RoundDown))
        @test isequal(rem(T(x), T(y), RoundUp), rem(x, y, RoundUp))
        @test isequal(mod(T(x), T(y)), mod(x, y))
        @test isequal(T(x) % T(y), x % y)
    end
    for x = -5.0:0.25:5.0
        @test isequal(modf(T(x)), modf(x))
    end
    for x in (NaN, -Inf, -0.0, Inf)
        @test isequal(modf(T(x)), modf(x))
    end

    @test trunc(T(2.7)) === floor(T(2.7)) === round(T(2.7), RoundDown) === round(T(2.7), RoundToZero) === T(2)
    @test ceil(T(2.3)) === round(T(2.3), RoundUp) === round(T(2.3), RoundFromZero) === T(3)
    @test round(T(1.5)) === round(T(2.5)) === round(T(1.5), RoundNearest) === round(T(2.5), RoundNearest) === T(2)
    @test round(T(2.5), RoundNearestTiesAway) === round(T(3.3), RoundNearestTiesAway) === T(3)
    @test round(T(2.5), RoundNearestTiesUp) === round(T(3.3), RoundNearestTiesUp) === T(3)

    for Ti in (Integer,Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128)
        if Ti != Integer
            @test parse(T, "17") == T(Ti(17)) == Ti(17) == Ti(T(17))
        end
        @test trunc(Ti, T(2.7)) === floor(Ti, T(2.7)) === round(Ti, T(2.7), RoundDown) === round(Ti, T(2.7), RoundToZero) === Ti(2)
        @test ceil(Ti, T(2.3)) === round(Ti, T(2.3), RoundUp) === round(Ti, T(2.3), RoundFromZero) === Ti(3)
        @test round(Ti, T(1.5)) === round(Ti, T(2.5)) === round(Ti, T(1.5), RoundNearest) === round(Ti, T(2.5), RoundNearest) === Ti(2)
        @test round(Ti, T(2.5), RoundNearestTiesAway) === round(Ti, T(3.3), RoundNearestTiesAway) === Ti(3)
        @test round(Ti, T(2.5), RoundNearestTiesUp) === round(Ti, T(3.3), RoundNearestTiesUp) === Ti(3)
        @test_throws InexactError convert(Ti, xd)
        @test_throws InexactError trunc(Ti, floatmax(T))
        @test_throws InexactError floor(Ti, floatmax(T))
        @test_throws InexactError ceil(Ti, floatmax(T))
        if Ti <: Signed
            @test parse(T, "-17") == T(Ti(-17)) == Ti(-17) == Ti(T(-17))
        end
        if Ti <: Signed || Ti === Integer
            @test trunc(Ti, T(-2.7)) === round(Ti, T(-2.7), RoundToZero) === Ti(-2)
            @test floor(Ti, T(-2.3)) === round(Ti, T(-2.3), RoundDown) === round(Ti, T(-2.3), RoundFromZero) === Ti(-3)
            @test ceil(Ti, T(-2.7)) === round(Ti, T(-2.7), RoundUp) === Ti(-2)
            @test round(Ti, T(-1.5)) === round(Ti, T(-2.5)) === round(Ti, T(-1.5), RoundNearest) === round(Ti, T(-2.5), RoundNearest) === Ti(-2)
            @test round(Ti, T(-2.5), RoundNearestTiesAway) === round(Ti, T(-3.3), RoundNearestTiesAway) === Ti(-3)
            @test round(Ti, T(-1.5), RoundNearestTiesUp) === round(Ti, T(-0.7), RoundNearestTiesUp) === Ti(-1)
            @test_throws InexactError convert(Ti, yd)
        end
    end

    @test rounding(T) == RoundNearest
    for mode in (RoundNearest, RoundToZero, RoundFromZero, RoundUp, RoundDown)
        setrounding(T, mode) do
            @test rounding(T) == mode
        end
    end
    @test rounding(T) == RoundNearest

    TI = eval(Symbol(string("UInt", sizeof(T)*8)))
    @test bswap(xd) == reinterpret(T, bswap(reinterpret(TI, xd)))

    @test parse(T, "1e10000") == T(Inf)
    @test_throws DomainError asin(xd)
    @test_throws DomainError sqrt(yd)
    @test_throws DomainError acosh(zd)

    @test ldexp10(parse(T, "1"), 3) == 1000
    @test exponent10(parse(T, "1000")) == 3
    @test sqrt(complex(yd)) ≈ sqrt(complex(y))

    @test typeof(xd * pi) == T
    @test typeof((xd+yd*im)*pi) == Complex{T}

    Threads.@threads for (i, mode) in collect(enumerate((RoundNearest, RoundToZero, RoundFromZero, RoundUp, RoundDown)))
        testthreads(T, i, mode)
    end

    # issue #85
    @test T(1.5) == T(T(1.5))
end

@test DecFP.flags[Threads.threadid()] == 0

@test widen(Dec32) == Dec64
@test widen(Dec64) == Dec128
@test widen(Dec32(1)) === Dec64(1)
@test widen(Dec64(1)) === Dec128(1)

@test hash(Dec32(9999999)) == hash(9999999)
@test hash(Dec32(-9999999)) == hash(-9999999)
@test hash(Dec64(9999999999999999)) == hash(9999999999999999)
@test hash(Dec64(-9999999999999999)) == hash(-9999999999999999)
@test hash(Dec128("9999999999999999999999999999999999")) == hash(9999999999999999999999999999999999)
@test hash(Dec128("-9999999999999999999999999999999999")) == hash(-9999999999999999999999999999999999)
@test hash(Dec128(2)^-30) == hash(BigFloat(2)^-30)
@test hash(-Dec128(2)^-30) == hash(-BigFloat(2)^-30)

# issue #37
@test reinterpret(UInt128, Dec128(1.5)) == 0x303e000000000000000000000000000f
# issue #38
@test collect(v for i in 1:1, v in zeros(Dec128, 1)) ==  zeros(Dec128, 1, 1)

@test setrounding(Dec64, RoundDown) do
    Float64(d64"1e100") < 1e100
end

@test Float64(d64"1e100") == 1e100

# issue #79
@test d64"1e100" != 1e100

# issue #93
@test parse(Dec64, "3935767060.093896713") == d64"3.935767060093897e9" ==
      Dec64(d128"3935767060.093896713")

# issue 124
@test parse(Complex{Dec64}, "1.0+2.0im") == Complex(d"1.0", d"2.0")
@test parse(Dec64, SubString("1.3x", 1,3)) == d"1.3"

# issue #139
@test convert(Dec64, big"2.5") == d"2.5"

@testset "precision" begin
    @test precision(Dec32, base=10) == 7
    @test precision(Dec64, base=10) == 16
    @test precision(Dec128, base=10) == 34
    @test precision(Dec64, base=100) == 8
    @test precision(d32"1.2345", base=10) == 7
end


text:
using RobotOS
@rosimport std_msgs.msg: Bool
rostypegen(@__MODULE__)
import .std_msgs.msg: BoolMsg

using ROSTransforms.tf
using Test

init_node("broadcaster", anonymous=true) # delete
tb = TransformBroadcaster()
tf_ref = Transform([1, 2, 3.], [0, 0, 0, 1.])

# Build the success callback function
success = false
function cb(msg::BoolMsg)
    global success = msg.data
end
sub1 = Subscriber{BoolMsg}("/success", cb, queue_size=10)

r = Rate(20.0)
for i ∈ 1:1000
    sendTransform(tb, tf_ref.trans, tf_ref.rot, RobotOS.now(), "child_link",
                  "parent_link")

    success && break  # Break if the listener heard the broadcaster
    rossleep(r)
end

@testset "TransformBroadcaster" begin
    @test success  # Success if listener heard
end


text:
ENV["MPLBACKEND"]="agg"

using Test, ExtendableGrids, GridVisualize
import PyPlot

@testset "Geomspace" begin
    function test_geomspace()
        X1=geomspace(2.0,3.0,0.2,0.2)
        X2=collect(2:0.2:3)
        @assert X1 ≈  X2
        X1=geomspace(2.0,3.0,0.1,0.2)
        X2=geomspace(2.0,3.0,0.2,0.1)

        X1=geomspace(2.0,3.0,1.0e-5,0.2)
        X2=geomspace(2.0,3.0,0.2,1.0e-5)
        (X1[2]-X1[1])  ≈ (X2[end]-X2[end-1]) &&
            (X2[2]-X2[1])  ≈ (X1[end]-X1[end-1])
    end

    function test_geomspace1(n,h0,h1)
        for i=1:n
            X=geomspace(0,1,h0*rand(),h1*rand())
        end
        true
    end

    @test test_geomspace()
    @test test_geomspace1(100,0.01,0.01)
    @test test_geomspace1(100,0.01,0.1)
    @test test_geomspace1(100,0.1,0.01)
    @test test_geomspace1(100,0.001,0.1)

    @test length(geomspace(0,1,0.01,0.1))==27
    @test length(geomspace(0,1,0.001,0.1))==47
    @test length(geomspace(0,1,0.0001,0.1))==68
    @test length(geomspace(0,1,0.00001,0.1))==90

end

@testset "Basic" begin

    function test_prepare_edges()
        ## Compared with pdelib; have more of these
        nodes=[
            0.0 1;
            1 0;
            0 1;
            1 1]'

        cells=[
            1 2 3;
            2 3 4
        ]'

        cellmat=[1,1]
        bfaces=[ 1 2;
                 1 3;
                 2 4]'

        bfacemat=[1,1]

        grid=simplexgrid(nodes,cells,cellmat,bfaces,bfacemat)

        grid[CellEdges]==Int32[1 2; 2 4; 3 5] &&
        grid[FaceNodes]==Int32[1 2 3 3 4; 2 3 1 4 2] &&
            grid[EdgeCells]==Int32[1 1 1 2 2; 0 2 0 0 0]

    end
    @test test_prepare_edges()

    @test let
        g=simplexgrid(0:10)
        bfc=g[BFaceCells]
        bfc[1,1]==1 && bfc[1,2]==10
    end

    @test let
        g=simplexgrid(0:10)
        bfn=g[BFaceNormals]
        bfn==[-1.0 1.0]
    end

    @test let
        g=simplexgrid(0:2, 0:2)
        bfn=g[BFaceNormals]
        bfn==[0.0  -1.0;
              -1.0   0.0;
              0.0  -1.0;
              1.0   0.0;
              0.0   1.0;
              -1.0   0.0;
              1.0   0.0;
              0.0   1.0]'
    end

    @test let
        g=simplexgrid(0:2, 0:2)
        bfc=g[BFaceCells]
        tryfix(bfc)==[1  2  3  3  6  6  7  8;]
    end

    @test let
        X=[0 0; 2 0 ; 1 2;  1 0.5]'
        C=[ 1 2 4 ; 2 3 4 ; 3 1 4]'
        CR=ones(Int,3)
        F=[ 1 2 ; 2 3 ; 3 1]'
        FR=[1,2,3]
        g=simplexgrid(X,C,CR,F,FR)
        bfc=g[BFaceCells]
        tryfix(bfc)==[1  2  3;]
    end

end

@testset "Read/Write" begin
    function testrw(grid)
        ftmp=tempname()
        write(ftmp,grid,format="sg")
        grid1=simplexgrid(ftmp,format="sg")
        seemingly_equal(grid1,grid)
    end
    X=collect(0:0.05:1)
    @test testrw(simplexgrid(X))
    @test testrw(simplexgrid(X,X))
    @test testrw(simplexgrid(X,X,X))
end

function testgrid(grid,testdata)
    (num_nodes(grid),num_cells(grid), num_bfaces(grid))==testdata
end

examples1d=joinpath(@__DIR__,"..","examples","examples1d.jl")
include(examples1d)
examples2d=joinpath(@__DIR__,"..","examples","examples2d.jl")
include(examples2d)
examples3d=joinpath(@__DIR__,"..","examples","examples3d.jl")
include(examples3d)

@testset "1D" begin
    @test testgrid(interval_from_vector(),(21,20,2))
    @test testgrid(interval_localref(),   (27,26,2))
    @test testgrid(interval_multiregion(),(21,20,3))
    @test testgrid(interval_subgrid(),(51,50,2))
end

@testset "2D" begin
    @test testgrid(rectangle(),(441,800,80))
    @test testgrid(rectangle_localref(),(729, 1352, 104))
    @test testgrid(rectangle_multiregion(),(441,800,100))
    @test testgrid(rectangle_subgrid(),(360, 600, 120))
end

@testset "3D" begin
    @test testgrid(quadrilateral(),(330,1200,440))
    @test mask_bedges()

    X=collect(0:0.25:1)
    gxy=simplexgrid(X,X)
    gxyz=simplexgrid(gxy,X)
    g=simplexgrid(X,X,X)
    @test testgrid(gxyz,(125,384,192))
    @test g[Coordinates]≈gxyz[Coordinates]
end

if !Sys.isapple()
    @testset "plotting examples" begin
        include("../docs/makeplots.jl")
        picdir=mktempdir()

        @test makeplot("interval_from_vector",picdir)
        @test makeplot("interval_localref",picdir)
        @test makeplot("interval_multiregion",picdir)
        @test makeplot("interval_subgrid",picdir)
        @test makeplot("rectangle",picdir)
        @test makeplot("rectangle_localref",picdir)
        @test makeplot("rectangle_multiregion",picdir)
        @test makeplot("rectangle_subgrid",picdir)
        @test makeplot("quadrilateral",picdir)
    end
end

function tglue(;dim=2,breg=0)
    X1=linspace(0,1,5)
    Y1=linspace(0,1,5)
    Z1=linspace(0,1,5)

    X2=linspace(1,2,5)
    Y2=linspace(0,0.5,3)
    Z2=linspace(0,0.5,3)

    if dim==1
        g1=simplexgrid(X1)
        g2=simplexgrid(X2)
    end

    if dim==2
        g1=simplexgrid(X1,Y1)
        g2=simplexgrid(X2,Y2)
    end

    if dim==3
        g1=simplexgrid(X1,Y1,Z1)
        g2=simplexgrid(X2,Y2,Z2)
    end

    glue(g1,g2,breg=breg)
end

@testset "Glue" begin
    @test testgrid(tglue(;dim=1,breg=0),(9,8,2))
    @test testgrid(tglue(;dim=1,breg=1),(9,8,3))
    @test testgrid(tglue(;dim=2,breg=0),(37, 48, 24))
    @test testgrid(tglue(;dim=2,breg=1),(37, 48, 26))
    @test testgrid(tglue(;dim=3,breg=0),(161, 480, 256))
    @test testgrid(tglue(;dim=3,breg=1),(161, 480, 264))
end

function voronoitest()
    g=simplexgrid(0:0.1:1)
    @test g[VoronoiFaceCenters]≈[0.05  0.15  0.25  0.35  0.45  0.55  0.65  0.75  0.85  0.95]
    g=simplexgrid(0:0.5:1,0:0.5:1)
    @test g[VoronoiFaceCenters]≈[0.25  0.5   0.25  0.25  0.0   0.75  1.0   0.75  0.75  0.5   0.25  0.25  0.0   1.0   0.75  0.75;
                                  0.0   0.25  0.25  0.5   0.25  0.0   0.25  0.25  0.5   0.75  0.75  1.0   0.75  0.75  0.75  1.0]
end

@testset "Voronoi" begin
    voronoitest()
end

@testset "Grid Stuff" begin
    include("test_gridstuff.jl")
    run_grid_tests()
end


text:
function solve(case, H)
    c = count(enumerate(H)) do (i, h)
        i != firstindex(H) && i != lastindex(H) && H[i - 1] < H[i] > H[i + 1]
    end

    println("Case #$(case): $c")
end

function main()
    T = parse(Int, readline(stdin))
    for case in 1:T
        _ = parse(Int, readline(stdin))
        H = parse.(Int, split(readline(stdin)))
        solve(case, H)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
using Plots
using Random
using LaTeXStrings
using JLD
include("./cmd.jl")
include("./simmd.jl")
include("./extramd.jl")
include("./pextramd.jl")
include("./AutoDiffUtils.jl")
include("./Potentials.jl")

@show α = 10
@show β = 1000

#reset past font scaling
Plots.scalefontsizes()
#scale fonts
Plots.scalefontsizes(1.75)

Random.seed!(123)

lightblue = colorant"rgb(63%,74%,78%)"
orange = colorant"rgb(85%,55%,13%)"
silver = colorant"rgb(69%,67%,66%)"
rust = colorant"rgb(72%,26%,6%)"
brown = colorant"rgb(57%,34%,16%)"

max_iter = 3000

m = 5000
n_terms = 50
n = 1

A = randn(n_terms, m)
b = vec(A[:, 1] + A[:, 2]) / 2 + 0.1 * randn(n_terms)

f(x, y) = sum((A * x - b).^2) + y[1] * (sum(x) - 1)
f_eval(x, y) = sum(abs.(A * x / sum(x) - b))
g(x, y) = - y[1] * (sum(x) - 1)

x0 = ones(m) / m
y0 = zeros(1)

let
    global scatter_cmd = plot()
    global xHist_cmd = zeros(m, max_iter + 1)
    global yHist_cmd = zeros(n, max_iter + 1)
    global iters_cmd = zeros(Int, max_iter + 1)
    global m
    global n
    global A
    global f
    global g
    global x0
    global y0

    ψ_x = shannon_entropy(α)
    ψ_y = squared_distance(β)

    x = ψ_x[2](x0)
    y = ψ_y[2](y0)

    opt = CMD(x, ψ_x, f, y, ψ_y, g, cg!)

    for k = 1 : max_iter
      xHist_cmd[:, k], yHist_cmd[:, k] = iterate_primal(opt)
      iters_cmd[k+1] = iters_cmd[k] + step!(opt)

      opt.x̄ .= ψ_x[2](max.(iterate_primal(opt)[1], 1e-30))
      opt.x̄ .= ψ_x[2](min.(iterate_primal(opt)[1], 1e30))
    end
    xHist_cmd[:, end], yHist_cmd[:, end]= iterate_primal(opt)
end

let
    global xHist_extragd = zeros(m, max_iter + 1)
    global yHist_extragd = zeros(n, max_iter + 1)
    global iters_extragd = zeros(Int, max_iter + 1)
    global m
    global n
    global A
    global f
    global g
    global x0
    global y0

    ψ_x = positive_projection(α)
    ψ_y = squared_distance(β)

    x = ψ_x[2](x0)
    y = ψ_y[2](y0)

    opt = ExtraMD(x, ψ_x, f, y, ψ_y, g)

    for k = 1 : max_iter
      xHist_extragd[:, k], yHist_extragd[:, k] = iterate_primal(opt)
      iters_extragd[k+1] = iters_extragd[k] + step!(opt)

      opt.x̄ .= ψ_x[2](max.(iterate_primal(opt)[1], 1e-30))
      opt.x̄ .= ψ_x[2](min.(iterate_primal(opt)[1], 1e30))
    end
    xHist_extragd[:, end], yHist_extragd[:, end]= iterate_primal(opt)
end

let
    global xHist_extramd = zeros(m, max_iter + 1)
    global yHist_extramd = zeros(n, max_iter + 1)
    global iters_extramd = zeros(Int, max_iter + 1)
    global m
    global n
    global A
    global f
    global g
    global x0
    global y0

    ψ_x = shannon_entropy(α)
    ψ_y = squared_distance(β)

    x = ψ_x[2](x0)
    y = ψ_y[2](y0)

    opt = ExtraMD(x, ψ_x, f, y, ψ_y, g)

    for k = 1 : max_iter
      xHist_extramd[:, k], yHist_extramd[:, k] = iterate_primal(opt)
      iters_extramd[k+1] = iters_extramd[k] + step!(opt)

      opt.x̄ .= ψ_x[2](max.(iterate_primal(opt)[1], 1e-30))
      opt.x̄ .= ψ_x[2](min.(iterate_primal(opt)[1], 1e30))
    end
    xHist_extramd[:, end], yHist_extramd[:, end]= iterate_primal(opt)
end

evals_cmd = [f_eval(xHist_cmd[:, k], yHist_cmd[:, k]) for k in 1 : size(xHist_cmd, 2)]
evals_extramd = [f_eval(xHist_extramd[:, k], yHist_extramd[:, k]) for k in 1 : size(xHist_extramd, 2)]
evals_extragd = [f_eval(xHist_extragd[:, k], yHist_extragd[:, k]) for k in 1 : size(xHist_extragd, 2)]

plot(evals_cmd, label="cmd", color=lightblue)
plot!(evals_extramd, label="extramd", color=orange)
plot!(evals_extragd, label="extragd", color=rust)

save("./out/jld/regression_alpha_$(α)_beta_$(β).jld",
#      "f", f,
#      "f_eval", f_eval,
#      "g", g,
     "A", A,
     "b", b,
     "xHist_cmd", xHist_cmd,
     "yHist_cmd", yHist_cmd,
     "iters_cmd", iters_cmd,
     "evals_cmd", evals_cmd,
     "xHist_extramd", xHist_extramd,
     "yHist_extramd", yHist_extramd,
     "iters_extramd", iters_extramd,
     "evals_extramd", evals_extramd,
     "xHist_extragd", xHist_extragd,
     "yHist_extragd", yHist_extragd,
     "iters_extragd", iters_extragd,
     "evals_extragd", evals_extragd )

max_stars_repo_name: f-t-s/CMD

max_issues_repo_name: f-t-s/CMD

max_forks_repo_name: f-t-s/CMD


text:
export Camera, get_primary_rays

# ------ #
# Camera #
# ------ #

"""
    FixedCameraParams

The Parameters of the [`Camera`](@ref) which should not be updated in inverse rendering
are wrapped into this object.

### Fields:

* `vup`    - Stores the World UP Vector (In most cases this is the Y-Axis Vec3(0.0, 1.0, 0.0))
* `width`  - Width of the image to be rendered
* `height` - Height of the image to be rendered
"""
struct FixedCameraParams{T} <: FixedParams
    vup::Vec3{T}
    width::Int
    height::Int
end

Base.show(io::IO, fcp::FixedCameraParams) =
    print(io, "    Fixed Parameters:\n        World UP - ", fcp.vup,
          "\n        Screen Dimensions - ", fcp.height, " × ", fcp.width)

Base.zero(fcp::FixedCameraParams) = FixedCameraParams(zero(fcp.vup), zero(fcp.width),
                                                      zero(fcp.height))

"""
    Camera

The Perspective View Camera Model

### Fields:

* `lookfrom`    - Position of the Camera in 3D World Space
* `lookat`      - Point in the 3D World Space where the Camera is Pointing
* `vfov`        - Field of View of the Camera
* `focus`       - The focal length of the Camera
* `fixedparams` - An instance of [`FixedCameraParams`](@ref)

### Available Constructors:

* `Camera(;lookfrom = Vec3(0.0f0), lookat = Vec3(0.0f0),
           vup = Vec3(0.0f0, 1.0f0, 0.0f0), vfov = 90.0f0,
           focus = 1.0f0, width = 128, height = 128)`
* `Camera(lookfrom::Vec3{T}, lookat::Vec3{T}, vup::Vec3{T}, vfov::R,
          focus::R, width::Int, height::Int) where {T<:AbstractArray, R<:Real}`
"""
struct Camera{T}
    lookfrom::Vec3{T}
    lookat::Vec3{T}
    vfov::T
    focus::T
    fixedparams::FixedCameraParams{T}
end

Base.show(io::IO, cam::Camera) =
    print(io, "CAMERA Configuration:\n    Lookfrom - ", cam.lookfrom,
          "\n    Lookat - ", cam.lookat, "\n    Field of View - ", cam.vfov[],
          "\n    Focus - ", cam.focus[], "\n", cam.fixedparams)

@diffops Camera

Camera(;lookfrom = Vec3(0.0f0), lookat = Vec3(0.0f0),
        vup = Vec3(0.0f0, 1.0f0, 0.0f0), vfov = 90.0f0,
        focus = 1.0f0, width = 128, height = 128) =
    Camera(lookfrom, lookat, vup, vfov, focus, width, height)

function Camera(lookfrom::Vec3{T}, lookat::Vec3{T}, vup::Vec3{T}, vfov::R,
                focus::R, width::Int, height::Int) where {T<:AbstractArray, R<:Real}
    fixedparams = FixedCameraParams(vup, width, height)
    return Camera(lookfrom, lookat, [vfov], [focus], fixedparams)
end

"""
    get_primary_rays(c::Camera)

Takes the configuration of the camera and returns the origin and the direction
of the primary rays.
"""
function get_primary_rays(c::Camera)
    width = c.fixedparams.width
    height = c.fixedparams.height
    vup = c.fixedparams.vup
    vfov = c.vfov[]
    focus = c.focus[]

    aspect_ratio = width / height
    half_height = tan(deg2rad(vfov / 2))
    half_width = typeof(half_height)(aspect_ratio * half_height)

    origin = c.lookfrom
    w = normalize(c.lookfrom - c.lookat)
    u = normalize(cross(vup, w))
    v = normalize(cross(w, u))

    # Lower Left Corner
    llc = origin - half_width * focus * u - half_height * focus * v - w
    hori = 2 * half_width * focus * u
    vert = 2 * half_height * focus * v

    s = repeat((collect(0:(width - 1)) .+ 0.5f0) ./ width, outer = height)
    t = repeat((collect((height - 1):-1:0) .+ 0.5f0) ./ height, inner = width)

    direction = normalize(llc + s * hori + t * vert - origin)

    return (origin, direction)
end

"""
    get_transformation_matrix(c::Camera)

Returns the `camera_to_world` transformation matrix. This is used for transforming
the coordinates of a Point in 3D Camera Space to 3D World Space using the
[`camera2world`](@ref) function.
"""
function get_transformation_matrix(c::Camera{T}) where {T}
    forward = normalize(c.lookfrom - c.lookat)
    right = normalize(cross(c.fixedparams.vup, forward))
    up = normalize(cross(forward, right))

    return [     right.x[]      right.y[]      right.z[] zero(eltype(T));
                    up.x[]         up.y[]         up.z[] zero(eltype(T));
               forward.x[]    forward.y[]    forward.z[] zero(eltype(T));
            c.lookfrom.x[] c.lookfrom.y[] c.lookfrom.z[]  one(eltype(T))]
end

"""
    compute_screen_coordinates(c::Camera, film_aperture::Tuple,
                               inch_to_mm::Real = 25.4)

Computes the coordinates of the 4 corners of the screen and returns
`top`, `right`, `bottom`, and `left`.
"""
function compute_screen_coordinates(c::Camera, film_aperture::Tuple,
                                    inch_to_mm::Real = 25.4)
    width = c.fixedparams.width
    height = c.fixedparams.height
    vfov = c.vfov[]
    focus = c.focus[]

    film_aspect_ratio = film_aperture[1] / film_aperture[2]
    device_aspect_ratio = width / height

    top = ((film_aperture[2] * inch_to_mm / 2) / focus)
    right = ((film_aperture[1] * inch_to_mm / 2) / focus)

    xscale = 1
    yscale = 1

    if film_aspect_ratio > device_aspect_ratio
        xscale = device_aspect_ratio / film_aspect_ratio
    else
        yscale = film_aspect_ratio / device_aspect_ratio
    end

    right *= xscale
    top *= yscale

    bottom = -top
    left = -right

    return top, right, bottom, left
end


text:
#     using MAT

    struct NCSSeg
        qwTimeStamp::UInt64
        dwChannelNumber::UInt32
        dwSampleFreq::UInt32
        dwNumValidSamples::UInt32
        snSample::Array{Int16,1}
    end

    struct NCSData
        filename::String
        header::String
        seg_num::UInt32
        data::Array{NCSSeg, 1}
        sample::Array{Int16, 1}
    end

    # read neuralynx file
    function readNCSFile(filename)

        f = open(filename)
        header = read(f, UInt8, 16*1024) .|> Char |> String
        bdata = reinterpret(Int16, read(f))
        close(f)

        records = length(bdata) / 522 |> Int
        data = Array{NCSSeg, 1}(records)

        for i = 1:records
            segdata = bdata[1+(i-1)*522:522*i]

            timestamp = reinterpret(UInt64, segdata[1:4])
            chnum = reinterpret(UInt32, segdata[5:6])
            samplefreq = reinterpret(UInt32, segdata[7:8])
            samplenum = reinterpret(UInt32, segdata[9:10])

            data[i] = NCSSeg(timestamp[1], chnum[1], samplefreq[1], samplenum[1], segdata[11:522])
        end

        sample = zeros(Int16, records*512)
        for idx = 1:records
            sample[1+512*(idx-1):512*idx] = data[idx].snSample
        end

        return NCSData(filename, header, records, data, sample)
    end

    # # save into mat-v7.3 file
    # function save(filename::String, ncs::NCSData; mode=:sampleonly)
    #
    #     if mode == :sampleonly
    #         matwrite(filename, Dict("header"=>ncs.header, "sample"=>ncs.sample))
    #     elseif mode == :all
    #         matwrite(filename,
    #             Dict("header"=>ncs.header,
    #                 "seg_num"=>ncs.seg_num,
    #                 "data"=>ncs.data,
    #                 "sample"=>ncs.sample))
    #     else
    #         matwrite(filename, Dict("header"=>ncs.header, "sample"=>ncs.sample))
    #     end
    # end


text:
using
  Documenter,
  Literate,
  Plots,  # so that Literate.jl does not capture precompilation output
  Glob,
  FourierFlows

# Gotta set this environment variable when using the GR run-time on CI machines.
# This happens as examples will use Plots.jl to make plots and movies.
# See: https://github.com/jheinen/GR.jl/issues/278
ENV["GKSwstype"] = "100"

#####
##### Generate examples
#####

const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
const OUTPUT_DIR   = joinpath(@__DIR__, "src/literated")

examples = [
    "OneDShallowWaterGeostrophicAdjustment.jl",
]

for example in examples
  example_filepath = joinpath(EXAMPLES_DIR, example)
  withenv("GITHUB_REPOSITORY" => "FourierFlows/FourierFlowsDocumentation") do
    Literate.markdown(example_filepath, OUTPUT_DIR; flavor = Literate.DocumenterFlavor())
    Literate.notebook(example_filepath, OUTPUT_DIR)
    Literate.script(example_filepath, OUTPUT_DIR)
  end
end

#####
##### Build and deploy docs
#####

# Set up a timer to print a space ' ' every 240 seconds. This is to avoid CI
# timing out when building demanding Literate.jl examples.
Timer(t -> println(" "), 0, interval=240)

format = Documenter.HTML(
    collapselevel = 2,
       prettyurls = get(ENV, "CI", nothing) == "true",
        canonical = "https://fourierflows.github.io/FourierFlowsDocumentation/dev/",
)

pages = [
    "Home" => "index.md",
    "Installation Instructions" => "installation_instructions.md",
    "Code Basics" => "basics.md",
    "Grids" => "grids.md",
    "Aliasing" => "aliasing.md",
    "Problem" => "problem.md",
    "Diagnostics" => "diagnostics.md",
    "Output" => "output.md",
    "GPU" => "gpu.md",
    "Examples" => [
        "literated/OneDShallowWaterGeostrophicAdjustment.md",
        ],
    "Contributor's guide" => "contributing.md",
    "Library" => [
        "Contents" => "library/outline.md",
        "Public" => "library/public.md",
        "Private" => "library/internals.md",
        "Function index" => "library/function_index.md",
        ],
]

makedocs(
   sitename = "FourierFlows.jl",
    authors = "Gregory L. Wagner and Navid C. Constantinou",
    modules = [FourierFlows],
     format = format,
      pages = pages,
    doctest = true,
     strict = :doctest,
      clean = true,
  checkdocs = :exports
)

withenv("GITHUB_REPOSITORY" => "FourierFlows/FourierFlowsDocumentation") do
  deploydocs(        repo = "github.com/FourierFlows/FourierFlowsDocumentation.git",
                 versions = ["stable" => "v^", "v#.#.#", "dev" => "dev"],
             push_preview = false,
                devbranch = "main"
  )
end

@info "Cleaning up temporary .jld2 and .nc files created by doctests..."
for file in vcat(glob("docs/*.jld2"))
    rm(file)
end


text:
"""
Stacked layer is a abstract type and defines common behavior for containers
layers that enables to access them as arrays.
"""
abstract type StackedLayer <: Layer end

add(model::StackedLayer, layers...) = push!(model.layers, layers...)

"""
Sequential layer allows to chain together a number of other layers.

# Usage

```julia
model = Sequential(Conv2D(100),MaxPool(),Dense(10))
```

"""
struct Sequential <:StackedLayer
    layers::Vector

	Sequential(blocks...) = new([blocks...])
end

function call(model::Sequential, X::Tensor)
	for layer in model.layers
		X = layer(X)
	end
	return X
end

"""
Concurrrent layer allows for stacking a number of other layers in parallel and
combining their results before returning it.

This layer will stack on the second last dimension.
So with 2D and 3D convolution this will be the channel layer (WxHxCxN). As a result
other dimensions have to the same.
"""
struct Concurrent <:StackedLayer
    layers::Vector

	Concurrent(blocks...) = new([blocks...])
end

function call(model::Concurrent, X::Tensor)
	out = []
	for layer in model.layers
		push!(out, layer(X))
	end
	cat(out..., dims=ndims(out[1])-1)
end

"""
Residual Layer works like a Sequential layer, however before returning the result
it will be combined with the orginal input (residual). This is a popular techique
in modern neural networds since it allows for better backpropagation.

This will stack on the second last dimension.
So with 2D and 3D convolution this will be the channel layer (WxHxCxN)
"""
struct Residual <:StackedLayer
    layers::Vector

	Residual(blocks...) = new([blocks...])
end

function call(model::Residual, X::Tensor)
	res = X
	for layer in model.layers
		X = layer(X)
	end
	cat(res, X, dims=ndims(res[1])-1)
end

@debug "Loaded Container modules"


text:
using BlackBoxOptim, Optim

cost_function = build_loss_objective(prob1,Tsit5(),L2Loss(t,data),
                                     maxiters=10000,verbose=false)
bound1 = Tuple{Float64, Float64}[(1, 2)]
result = bboptimize(cost_function;SearchRange = bound1, MaxSteps = 11e3)
@test result.archive_output.best_candidate[1] ≈ 1.5 atol=3e-1

cost_function = build_loss_objective(prob2,Tsit5(),L2Loss(t,data,differ_weight=nothing,data_weight=1.0),
                                     maxiters=10000,verbose=false)
bound2 = Tuple{Float64, Float64}[(1, 2),(1, 4)]
result = bboptimize(cost_function;SearchRange = bound2, MaxSteps = 11e3)
@test result.archive_output.best_candidate ≈ [1.5;3.0] atol=3e-1

cost_function = build_loss_objective(prob3,Tsit5(),L2Loss(t,data,differ_weight=10),
                                     maxiters=10000,verbose=false)
bound3 = Tuple{Float64, Float64}[(1, 2),(0, 2), (2, 4), (0, 2)]
result = bboptimize(cost_function;SearchRange = bound3, MaxSteps = 11e3)
@test result.archive_output.best_candidate ≈ [1.5;1.0;3.0;1.0] atol=5e-1

cost_function = build_loss_objective(prob3,Tsit5(),L2Loss(t,data,differ_weight=0.3,data_weight=0.7),
                                     maxiters=10000,verbose=false)
bound3 = Tuple{Float64, Float64}[(1, 2),(0, 2), (1, 4), (0, 2)]
result = bboptimize(cost_function;SearchRange = bound3, MaxSteps = 11e3)
@test result.archive_output.best_candidate ≈ [1.5;1.0;3.0;1.0] atol=5e-1

println("LSA for gradient")
cost_function = build_loss_objective(prob1,Tsit5(),L2Loss(t,data),
                                     maxiters=10000,verbose=false,flsa_gradient=true)
bound1 = Tuple{Float64, Float64}[(1, 2)]
function g_!(stor,x)
    cost_function.cost_function2(x,stor)
end
result = Optim.optimize(cost_function.cost_function,g_!,[1.0], Optim.ConjugateGradient())
@test result.minimizer[1] ≈ 1.5 atol=3e-1

println("Adjoint SA for gradient")
cost_function = build_loss_objective(prob2,Tsit5(),L2Loss(t,data),
                                     maxiters=10000,verbose=false,adjsa_gradient=true)
function g_!(stor,x)
    cost_function.cost_function2(x,stor)
end
result = Optim.optimize(cost_function.cost_function,g_!,[1.3,2.0], Optim.ConjugateGradient())
@test result.minimizer ≈ [1.5;3.0] atol=3e-1


text:
module DumbCompleter

using JSON
using Pkg
using Serialization

export ioserver

# A completion.
struct Leaf
    name::Symbol
    type::Type
    mod::Module
    Leaf(n::Symbol, ::T, m::Module) where T = new(n, T, m)
end

# A completion tree.
struct Tree
    lf::Union{Leaf, Nothing}
    tr::Dict{Char, Tree}
    Tree(lf::Union{Leaf, Nothing}=nothing) = new(lf, Dict{Char, Tree}())
end

Base.show(io::IO, tr::Tree) =
    print(io, "Completion tree with ", length(leaves(tr)), " entries")

# Stored completions.
struct Deps
    exports::Tree
    modules::Dict{Symbol, Tree}
end

const EXPORTS = Ref(Tree())
const MODULES = Ref(Dict{Symbol, Tree}())
const DEPS = joinpath(dirname(@__DIR__), "deps", "completions-$VERSION")
const Command = Dict{Symbol, Union{String, Nothing}}

__init__() = try loaddeps!() catch; loadstoredeps!() end

# Load the completions file.
function loaddeps!()
    deps = open(deserialize, DEPS)
    EXPORTS[] = deps.exports
    MODULES[] = deps.modules
end

# Load default completions and store the completions file.
function loadstoredeps!()
    @info "Loading + storing default completions"
    foreach(loadmodule!, (Core, Base))
    mkpath(dirname(DEPS))
    deps = Deps(EXPORTS[], MODULES[])
    open(io -> serialize(io, deps), DEPS, "w")
end

# Store a single completion.
function putleaf!(tr::Tree, lf::Leaf)
    name = string(lf.name)
    ind = 1
    while ind <= length(name)
        tr = get!(tr.tr, name[ind], Tree(ind == length(name) ? lf : nothing))
        ind = nextind(name, ind)
    end
end

# Compute the normalized module name.
modkey(m::Module) = Symbol(replace(string(m), "Main." => ""))

# Determines whether a symbol is worth using as a completion.
cancomplete(m::Module) = s -> cancomplete(s, m)
cancomplete(s::Symbol, m::Module) = isdefined(m, s) && !startswith(string(s), "#")

# Load a module's names and store them as completions.
function loadmodule!(m::Module; exports::Bool=true)
    m === Main && return
    k = modkey(m)
    haskey(MODULES[], k) && return

    # Populate a tree for this module.
    MODULES[][k] = tr = Tree()
    ns = filter(cancomplete(m), names(m; all=true, imported=true))
    foreach(ns) do n
        lf = Leaf(n, getfield(m, n), m)
        putleaf!(tr, lf)
        exports && Base.isexported(m, n) && putleaf!(EXPORTS[], lf)
    end

    # Do the same for submodules, but don't fill in exports.
    children = map(n -> getfield(m, n), filter(n -> getfield(m, n) isa Module, ns))
    foreach(m -> loadmodule!(m; exports=false), children)
end

function loadmodule!(m::Symbol)
    @eval import $m
    return loadmodule!(eval(m))
end

loadmodule!(m::AbstractString) = loadmodule!(Symbol(m))

# Get all of a tree's completions.
function leaves(tr::Tree)
    lvs = map(tr -> tr.lf, Iterators.filter(tr -> tr.lf !== nothing, values(tr.tr)))
    return vcat(lvs, mapreduce(leaves, vcat, values(tr.tr); init=Leaf[])...)
end

"""
    completions(::Tree, s::AbstractString) -> Vector{Leaf}
    completions(s::AbstractString)
    completions(s::AbstractString, m::Union{Module, Symbol, AbstractString})

Get completions that begin with `s`.
"""
function completions(tr::Tree, s::AbstractString)::Vector{Leaf}
    for c in s
        haskey(tr.tr, c) || return Leaf[]
        tr = tr.tr[c]
    end
    lvs = sort!(leaves(tr); by=lf -> lf.name)
    return tr.lf === nothing ? lvs : [tr.lf; lvs]
end

completions(s::AbstractString, ::Nothing=nothing) = completions(EXPORTS[], s)
completions(s::AbstractString, m::AbstractString) = completions(s, Symbol(m))
completions(s::AbstractString, m::Symbol) =
    completions(s, try eval(m) catch; return Leaf[] end)
completions(s::AbstractString, m::Module) =
    completions(get(MODULES[], modkey(m), Tree()), s)

"""
    activate!(path::AbstractString)

Activate a project and load all of its modules' completions.
"""
function activate!(path::AbstractString)
    toml = joinpath(path, "Project.toml")
    project = Pkg.Types.read_project(toml)
    current = Base.current_project()
    Pkg.activate(path)
    foreach(loadmodule!, keys(projectdeps(project)))
    loadmodule!(projectname(project))
    current === nothing ? Pkg.activate() : Pkg.activate(current)
end

# Pkg's APIs differ between versions.
if VERSION >= v"1.1"
    projectdeps(project::Pkg.Types.Project) = project.deps
    projectname(project::Pkg.Types.Project) = project.name
else
    projectdeps(project::Dict) = project["deps"]
    projectname(project::Dict) = project["name"]
end

JSON.lower(lf::Leaf) = Dict(
    :name => string(lf.name),
    :type => lf.type <: Function ? "Function" : string(lf.type),
    :module => string(modkey(lf.mod)),
)

# Print out some JSON with a newline.
jsonprintln(io::IO, x) = println(io, sprint(JSON.print, x))

"""
    ioserver(input::IO=stdin, output::IO=stdout)

Run a server that listens to some input and prints to some output.
This function does not return.
"""
function ioserver(input::IO=stdin, output::IO=stdout)
    while isopen(input)
        try
            c = JSON.parse(readline(input); dicttype=Command)
            jsonprintln(output, docmd(Val(Symbol(c[:type])), c))
        catch e
            isopen(input) &&
                jsonprintln(output, (; error=sprint(showerror, e), completions=[]))
        end
    end
end

# Do a client command.
docmd(::Val{T}, ::Command) where T = (; error="unknown command type $T", completions=[])
docmd(::Val{nothing}, ::Command) = (; error="type cannot be null", completions=[])
docmd(::Val{:activate}, c::Command) = (@async activate!(c[:path]); (; error=nothing))
docmd(::Val{:completions}, c::Command) = (;
    error=nothing,
    completions=completions(c[:text], c[:module]),
    exports=c[:module] === nothing,
)

end


text:
"""
    theme_aps(width=3.375, colors=COLORS, markers=MARDERS, cycle=CYCLE)

Generate Makie theme for producing figures for APS (American Physical Society). Save the figure using `savefig("path_to_figure.pdf", fig)` or `save("path_to_figure.pdf", fig, pt_per_unit=1.0)`.
"""
function theme_aps(;
                   width=3.25,
                   colors=COLORS,
                   markers=MARKERS,
                   cycle=CYCLE,
                   )
    axis_theme = (
        xlabelsize=10,
        ylabelsize=10,
        spinewidth=1.1,
        xticklabelsize=8,
        yticklabelsize=8,
        # xgridstyle=:dash, ygridstyle=:dash,
        xgridvisible=true,
        ygridvisible=true,
        xtickalign=1,
        ytickalign=1,
        xticksize=5,
        yticksize=5,
        xtickwidth=0.8,
        ytickwidth=0.8,
        xminorticksvisible=true,
        yminorticksvisible=true,
        xminortickalign=1,
        yminortickalign=1,
        xminorticks=IntervalsBetween(5),
        yminorticks=IntervalsBetween(5),
        xminorticksize=3,
        yminorticksize=3,
        xminortickwidth=0.75,
        yminortickwidth=0.75,
        xlabelpadding=-2,
        ylabelpadding=2,
    )

    scatter_theme = (
        cycle=cycle,
        markersize=7,
    )

    legend_theme = (
        nbanks=1,
        framecolor=(:grey, 0.5),
        framevisible=false,
        labelsize=7.5,
        padding=(2, 2, 2, 2),
        margin=(0, 0, 0, 0),
        # position=:rt, # l=left, r=right, c=center; b=bottom, t=top, c=center
        rowgap=-10,
        colgap=4,
    )

    pal = isnothing(markers) ? (color=colors,) : (color=colors, marker=markers)

    return Theme(figure_padding=0,
                 resolution=figsize(width),
                 # font="Helvetica",
                 palette=pal,
                 Axis=axis_theme,
                 Scatter=scatter_theme,
                 Legend=legend_theme,)
end

theme_aps_1col(;
               colors=COLORS,
               markers=MARKERS,
               cycle=CYCLE,
    ) = theme_aps(width=3.375, colors=colors, markers=markers, cycle=cycle)

theme_aps_2col(;
               colors=COLORS,
               markers=MARKERS,
               cycle=CYCLE,
    ) = theme_aps(width=6.75, colors=colors, markers=markers, cycle=cycle)


text:
using Parameters, StaticArrays, BenchmarkTools
import TrajectoryOptimization: dynamics, discrete_dynamics, RK3, AbstractModel, KnotPoint,
    jacobian, discrete_jacobian

@with_kw mutable struct BicycleCar{T} <: AbstractModel
    k::T = 0.1
    a::T = 1.  # dist to front axle
    b::T = 1.  # dist to rear axle
    h::T = 1.
    mass::T = 1000.  # mass
    Iz::T = 10.
    Cf::T = 1.
    Cr::T = 1.
    μ::T = 0.1
    Cd0::T = 0.1
    Cd1_Npmps::T = 0.1
    g::T = 9.81
end

Base.size(::BicycleCar) = 8,2

function dynamics(car::BicycleCar, x, u)
    δ_dot  = u[1]
    fx_dot = u[2]
    δ  = x[1]
    fx = x[2]
    r  = x[3]
    Uy = x[4]
    Ux = x[5]
    Δψ = x[6]
    e  = x[7]
    t  = x[8]

    # Drag Force
    Fx_drag = 0.0

    # Spatial derivative
    s_dot = (Ux * cos(Δψ) - Uy * sin(Δψ)) / (1 - car.k*e)

    # Get tire forces
    Fxf, Fxr, Fyf, Fyr = logit_lateral_force_model(car, x)

    # State Derivatives
    r_dot  = (car.a * (Fxf * cos(δ) + Fxf*sin(δ)) - car.b * Fyr) / car.Iz
    Ux_dot =  r * Uy + (Fxf * cos(δ) - Fyf * sin(δ) + Fxr - Fx_drag) / car.mass
    Uy_dot = -r * Ux + (Fyf * cos(δ) + Fxf * sin(δ)) / car.mass
    Δψ_dot = r - car.k * s_dot
    e_dot = Ux * sin(Δψ) + Uy * cos(Δψ)
    t_dot = 1 / s_dot

    return @SVector [δ_dot, fx_dot, r_dot, Uy_dot, Ux_dot, Δψ_dot, e_dot, t_dot]
end

function FWD_force_model(car::BicycleCar, fx)
    l = car.a + car.b
    Fxf = max(fx * car.b / l, fx)
    Fxr = min(fx * car.a / l, 0.0)
    Fzf = car.mass * car.g * car.b / l
    Fzr = car.mass * car.g * car.a / l
    return Fxf, Fxr, Fzf, Fzr
end

function logit_lateral_force_model(car::BicycleCar, x)
    fx = x[2]
    r  = x[3]
    Uy = x[4]
    Ux = x[5]
    Cα = 1.0

    α_f = atan(Uy + car.a * r, Ux)
    α_r = atan(Uy - car.b * r, Ux)

    Fxf, Fxr, Fzf, Fzr = FWD_force_model(car, fx)
    Fyf = lateral_tire_model(car.μ, Cα, α_f, Fxf, Fzf)
    Fyr = lateral_tire_model(car.μ, Cα, α_r, Fxr, Fzr)
    return Fxf, Fxr, Fyf, Fyr
end

function lateral_tire_model(μ, Cα, α, Fx, Fz)
    Fy_max = sqrt((μ*Fz)^2 - Fx^2)
    slope = 2*Cα / Fy_max
    Fy = Fy_max * (1-2*exp(slope*α)) / (1 + exp(slope*α))
    return Fy
end

struct ContingencyCar{T} <: AbstractModel
    cars::Vector{BicycleCar{T}}
    ix::Vector{SVector{8,Int}}
    iu::Vector{SVector{2,Int}}
end

function ContingencyCar(num_cars=2)
    cars = [BicycleCar() for i = 1:num_cars]
    n,m = 8,2
    ix = [SVector{n}((1:n) .+ (i-1)*n) for i = 1:num_cars]
    iu = [SVector{m}((1:m) .+ (i-1)*m) for i = 1:num_cars]
    ContingencyCar(cars, ix, iu)
end

Base.size(::ContingencyCar) = 16,4

function dynamics(car::ContingencyCar, x, u)
    ix, iu = car.ix, car.iu
    cars = car.cars
    x1,u1 = x[ix[1]], u[iu[1]]
    ẋ1 = dynamics(cars[1], x1, u1)
    x2,u2 = x[ix[2]], u[iu[2]]
    ẋ2 = dynamics(cars[2], x1, u1)
    return [ẋ1; ẋ2]
end

car = ContingencyCar()
n,m = size(car)
x = @SVector rand(n)
u = @SVector rand(m)
z = KnotPoint(x,u,0.1)

@btime dynamics($car, $x, $u)

@btime discrete_dynamics(RK3, $car, $x, $u, 0.1)

@btime discrete_jacobian(RK3, $car, $z)


text: include(joinpath(@__DIR__(), "..", "..", "graphidx", "julia", "src", "GraphIdx.jl"))

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
@testset "Overloads" begin
    @testset "Parent navigation" begin
        test_ex(Parent(HandleDefinition(handle_by_name(:VkDevice))), :(
            parent(device::Device) = device.physical_device
        ))

        test_ex(Parent(HandleDefinition(handle_by_name(:VkSurfaceKHR))), :(
            parent(surface::SurfaceKHR) = surface.instance
        ))
    end
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
"""
    foldxt(step, xf, reducible; [init, simd, basesize, stoppable, nestlevel]) :: T

e**X**tended **t**hreaded fold (reduce).  This is a multi-threaded
`reduce` based on extended fold protocol defined in Transducers.jl.

The "bottom" reduction function `step(::T, ::T) :: T` must be
associative and `init` must be its identity element.

Transducers composing `xf` must be stateless (e.g., [`Map`](@ref),
[`Filter`](@ref), [`Cat`](@ref), etc.) except for [`ScanEmit`](@ref).
Note that [`Scan`](@ref) is not supported (although possible in
theory).  Early termination requires Julia ≥ 1.3.

Use [`tcollect`](@ref) or [`tcopy`](@ref) to collect results into a
container.

See also: [Parallel processing tutorial](@ref tutorial-parallel),
[`foldxl`](@ref), [`foldxd`](@ref).

# Keyword Arguments
- `basesize::Integer = amount(reducible) ÷ nthreads()`: A size of
  chunk in `reducible` that is processed by each worker.  A smaller
  size may be required when:
    * computation time for processing each item fluctuates a lot
    * computation can be terminated by [`reduced`](@ref) or
      transducers using it, such as [`ReduceIf`](@ref)
- `stoppable::Bool`: [This option usually does not have to be set
  manually.]  Transducers.jl's `foldxt` executed in the "stoppable"
  mode used for optimizing reduction with [`reduced`](@ref) has a
  slight overhead if `reduced` is not used.  This mode can be disabled
  by passing `stoppable = false`.  It is usually automatically
  detected and set appropriately.  Note that this option is purely for
  optimization and does not affect the result value.
- `nestlevel::Union{Integer,Val}`: Specify how many inner `Cat`
  (flatten) transducers to be multi-threaded (using [`TCat`](@ref)).
  It must be a positive integer, `Val` of positive integer, or
  `Val(:inf)`.  `Val(:inf)` means to use multi-threading for all `Cat`
  transducers.  Note that `Cat` transducer should be statically known.
  That is to say, `foldxt` sees two `Cat`s in `... |> Map(f) |> Cat()
  |> Cat()` but only one `Cat` in `... |> Map(x -> f(x) |> Cat()) |>
  Cat()` even though they are semantically identical.
- For other keyword arguments, see [`foldl`](@ref).

!!! compat "Transducers.jl 0.4.23"

    Keyword option `stoppable` requires at least Transducers.jl 0.4.23.

# Examples
```jldoctest
julia> using Transducers

julia> foldxt(+, 1:3 |> Map(exp) |> Map(log))
6.0

julia> using BangBang: append!!

julia> foldxt(append!!, Map(x -> 1:x), 1:2; basesize=1, init=Union{}[])
3-element Vector{Int64}:
 1
 1
 2

julia> 1:5 |> Filter(isodd) |> foldxt(+)
9

julia> foldxt(TeeRF(min, max), [5, 2, 6, 8, 3])
(2, 8)
```
"""
foldxt

const _MAPREDUCE_DEPWARN = (
    "`mapreduce(::Transducer, rf, itr)` is deprecated. " *
    " Use `foldxt(rf, ::Transducer, itr)` if you do not need to call single-argument" *
    " `rf` on `complete`." *
    " Use `foldxt(whencomplete(rf, rf), ::Transducer, itr)` to call the" *
    " single-argument method of `rf` on complete."
)

"""
    mapreduce(xf, step, reducible; init, simd)

!!! warning

    $_MAPREDUCE_DEPWARN

Like [`foldxt`](@ref) but `step` is _not_ automatically wrapped by
[`Completing`](@ref).
"""
Base.mapreduce

"""
    Transducers.issmall(reducible, basesize) :: Bool

Check if `reducible` collection is considered small compared to
`basesize` (an integer).  Fold functions such as [`foldxt`](@ref)
switches to sequential `__foldl__` when `issmall` returns `true`.

Default implementation is `amount(reducible) <= basesize`.
"""
issmall

function transduce_assoc(
    xform::Transducer,
    step::F,
    init,
    coll0;
    simd::SIMDFlag = Val(false),
    basesize::Union{Integer,Nothing} = nothing,
    stoppable::Union{Bool,Nothing} = nothing,
    nestlevel::Union{Val,Integer,Nothing} = nothing,
) where {F}
    rf0 = _reducingfunction(xform, step; init = init)
    rf, coll = retransform(rf0, coll0)
    if nestlevel !== nothing
        if basesize === nothing
            throw(ArgumentError("`nestlevel` requires `basesize`"))
        end
        if has(rf, Union{Cat,TCat})
            rf = use_threads_for_inner_cats(rf, basesize, nestlevel)
            basesize = 1
        end
    end
    if stoppable === nothing
        stoppable = _might_return_reduced(rf, init, coll)
    end
    acc = @return_if_reduced _transduce_assoc_nocomplete(
        maybe_usesimd(rf, simd),
        init,
        coll,
        basesize === nothing ? amount(coll) ÷ Threads.nthreads() : basesize,
        stoppable ? CancellableDACContext() : NoopDACContext(),
    )
    result = complete(rf, acc)
    if unreduced(result) isa DefaultInitOf
        throw(EmptyResultError(rf))
        # See how `transduce(rf, init, coll)` is implemented in ./processes.jl
    end
    return result
end

if VERSION >= v"1.3-alpha"
    maybe_collect(coll) = coll
else
    maybe_collect(coll::AbstractArray) = coll
    maybe_collect(coll) = collect(coll)
end

function _transduce_assoc_nocomplete(
    rf::F,
    init,
    coll,
    basesize,
    ctx::DACContext = NoopDACContext(),
) where {F}
    reducible = SizedReducible(maybe_collect(coll), basesize)
    @static if VERSION >= v"1.3-alpha"
        return _reduce(ctx, rf, init, reducible)
    else
        return _reduce_threads_for(rf, init, reducible)
    end
end

function _reduce(ctx, rf::R, init::I, reducible::Reducible) where {R,I}
    if should_abort(ctx)
        return init
    end
    if issmall(reducible)
        acc = _reduce_basecase(rf, init, reducible)
        if acc isa Reduced
            cancel!(ctx)
        end
        return acc
    else
        left, right = _halve(reducible)
        fg, bg = splitcontext(ctx)
        task = @spawn _reduce(bg, rf, init, right)
        a0 = _reduce(fg, rf, init, left)
        b0 = fetch(task)
        a = @return_if_reduced a0
        should_abort(ctx) && return a  # slight optimization
        b0 isa Reduced && return combine_right_reduced(rf, a, b0)
        return combine(rf, a, b0)
    end
end

combine_right_reduced(rf, a, b0::Reduced) =
    reduced(combine(_realbottomrf(rf), a, unreduced(b0)))

# TODO: make this the default behavior of `combine`
function combine_reduced(rf, a, b)
    a isa Reduced && return a
    b isa Reduced && return combine_right_reduced(rf, a, b)
    return combine(rf, a, b)
end

function _reduce_threads_for(rf, init, reducible::SizedReducible{<:AbstractArray})
    arr = reducible.reducible
    basesize = reducible.basesize
    nthreads = max(
        1,
        basesize <= 1 ? length(arr) : length(arr) ÷ basesize
    )
    if nthreads == 1
        return foldl_nocomplete(rf, start(rf, init), arr)
    else
        w = length(arr) ÷ nthreads
        results = Vector{Any}(undef, nthreads)
        Threads.@threads for i in 1:nthreads
            if i == nthreads
                chunk = @view arr[(i - 1) * w + 1:end]
            else
                chunk = @view arr[(i - 1) * w + 1:i * w]
            end
            results[i] = foldl_nocomplete(rf, start(rf, init), chunk)
        end
        # It can be done in `log2(n)` for loops but it's not clear if
        # `combine` is compute-intensive enough so that launching
        # threads is worth enough.  Let's merge the `results`
        # sequentially for now.
        return combine_all(rf, results)
    end
end

function combine_all(rf, results)
    step = combine_step(rf)
    return transduce(ensurerf(Completing(step)), Init(step), results)
end

combine_step(rf) =
    asmonoid() do a0, b0
        a = @return_if_reduced a0
        b0 isa Reduced && return combine_right_reduced(rf, a, b0)
        return combine(rf, a, b0)
    end

# The output of `foldxt` is correct regardless of the value of
# `stoppable`.  Thus, we can use `return_type` here purely for
# optimization.
_might_return_reduced(rf, init, coll) =
    Base.typeintersect(
        Core.Compiler.return_type(
            _reduce_dummy,  # simulate the output type of `_reduce`
            typeof((rf, init, coll)),
        ),
        Reduced,
    ) !== Union{}

_reduce_dummy(rf, init, coll) =
    __reduce_dummy(rf, init, SizedReducible(maybe_collect(coll), 1))

function __reduce_dummy(rf, init, reducible)
    if issmall(reducible)
        return _reduce_basecase(rf, init, reducible)
    else
        left, right = halve(reducible)
        a = _reduce_dummy(rf, init, left)
        b = _reduce_dummy(rf, init, right)
        a isa Reduced && return a
        b isa Reduced && return combine_right_reduced(rf, a, b)
        return combine(rf, a, b)
    end
end

foldxt(step::F, xform::Transducer, itr; init = DefaultInit, kwargs...) where {F} =
    unreduced(transduce_assoc(xform, Completing(step), init, itr; kwargs...))

foldxt(step::F, foldable; kwargs...) where {F} =
    foldxt(step, extract_transducer(foldable)...; kwargs...)

foldxt(rf; kw...) = itr -> foldxt(rf, itr; kw...)

"""
    tcopy(xf::Transducer, T, reducible; basesize) :: Union{T, Empty{T}}
    tcopy(xf::Transducer, reducible::T; basesize) :: Union{T, Empty{T}}
    tcopy([T,] itr; basesize) :: Union{T, Empty{T}}

Thread-based parallel version of [`copy`](@ref).
Keyword arguments are passed to [`foldxt`](@ref).

See also: [Parallel processing tutorial](@ref tutorial-parallel)
(especially [Example: parallel `collect`](@ref tutorial-parallel-collect)).

!!! compat "Transducers.jl 0.4.5"

    New in version 0.4.5.

!!! compat "Transducers.jl 0.4.8"

    `tcopy` now accepts iterator comprehensions and eductions.

# Examples
```jldoctest
julia> using Transducers

julia> tcopy(Map(x -> x => x^2), Dict, 2:2)
Dict{Int64, Int64} with 1 entry:
  2 => 4

julia> using TypedTables

julia> @assert tcopy(Map(x -> (a=x,)), Table, 1:1) == Table(a=[1])

julia> using StructArrays

julia> @assert tcopy(Map(x -> (a=x,)), StructVector, 1:1) == StructVector(a=[1])
```

`tcopy` works with iterator comprehensions and eductions (unlike
[`copy`](@ref), there is no need for manual conversion with
[`eduction`](@ref)):

```jldoctest; setup = :(using Transducers, StructArrays, DataFrames)
julia> table = StructVector(a = [1, 2, 3], b = [5, 6, 7]);

julia> @assert tcopy(
           (A = row.a + 1, B = row.b - 1) for row in table if isodd(row.a)
       ) == StructVector(A = [2, 4], B = [4, 6])

julia> @assert tcopy(
           DataFrame,
           (A = row.a + 1, B = row.b - 1) for row in table if isodd(row.a)
       ) == DataFrame(A = [2, 4], B = [4, 6])

julia> @assert table |>
           Filter(row -> isodd(row.a)) |> Map(row -> (A = row.a + 1, B = row.b - 1)) |>
           tcopy == StructVector(A = [2, 4], B = [4, 6])
```

If you have [`Cat`](@ref) or [`MapCat`](@ref) at the end of the
transducer, consider using [`foldxt`](@ref) directly:

```jldoctest
julia> using Transducers
       using DataFrames

julia> @assert tcopy(
           DataFrame,
           1:2 |> Map(x -> DataFrame(a = [x])) |> MapCat(eachrow);
           basesize = 1,
       ) == DataFrame(a = [1, 2])

julia> using BangBang: Empty, append!!

julia> @assert foldxt(
           append!!,
           Map(x -> DataFrame(a = [x])),
           1:2;
           basesize = 1,
           # init = Empty(DataFrame),
       ) == DataFrame(a = [1, 2])
```

Note that above snippet assumes that it is OK to mutate the dataframe
returned by the transducer.  Use `init = Empty(DataFrame)` if this is
not the case.

This approach of using `foldxt` works with other containers; e.g.,
with `TypedTables.Table`:

```jldoctest; setup = :(using Transducers)
julia> using TypedTables

julia> @assert foldxt(
           append!!,
           Map(x -> Table(a = [x])),
           1:2;
           basesize = 1,
           # init = Empty(Table),
       ) == Table(a = [1, 2])
```
"""
tcopy(xf, T, reducible; kwargs...) =
    foldxt(append!!, Map(SingletonVector) ∘ xf, reducible; init = Empty(T), kwargs...)
tcopy(xf, reducible; kwargs...) = tcopy(xf, _materializer(reducible), reducible; kwargs...)

function tcopy(::Type{T}, itr; kwargs...) where {T}
    xf, foldable = extract_transducer(itr)
    return tcopy(xf, T, foldable; kwargs...)
end

function tcopy(itr; kwargs...)
    xf, foldable = extract_transducer(itr)
    return tcopy(xf, foldable; kwargs...)
end

tcopy(xf, T::Type{<:AbstractSet}, reducible; kwargs...) =
    foldxt(union!!, Map(SingletonVector) ∘ xf, reducible; init = Empty(T), kwargs...)

function tcopy(
    ::typeof(Map(identity)),
    T::Type{<:AbstractSet},
    array::PartitionableArray;
    basesize::Integer = max(1, length(array) ÷ Threads.nthreads()),
    kwargs...,
)
    @argcheck basesize >= 1
    return foldxt(
        union!!,
        Map(identity),
        Iterators.partition(array, basesize);
        init = Empty(T),
        basesize = 1,
        kwargs...,
    )
end

"""
    tcollect(xf::Transducer, reducible; basesize) :: Union{Vector, Empty{Vector}}
    tcollect(itr; basesize) :: Union{Vector, Empty{Vector}}

Thread-based parallel version of [`collect`](@ref).
This is just a short-hand notation of `tcopy(xf, Vector, reducible)`.
Use [`tcopy`](@ref) to get a container other than a `Vector`.

See also: [Parallel processing tutorial](@ref tutorial-parallel)
(especially [Example: parallel `collect`](@ref tutorial-parallel-collect)).

!!! compat "Transducers.jl 0.4.5"

    New in version 0.4.5.

!!! compat "Transducers.jl 0.4.8"

    `tcollect` now accepts iterator comprehensions and eductions.

# Examples
```jldoctest
julia> using Transducers

julia> tcollect(Map(x -> x^2), 1:2)
2-element Vector{Int64}:
 1
 4

julia> tcollect(x^2 for x in 1:2)
2-element Vector{Int64}:
 1
 4
```
"""
tcollect(xf, reducible; kwargs...) = tcopy(xf, Vector, reducible; kwargs...)
tcollect(itr; kwargs...) = tcollect(extract_transducer(itr)...; kwargs...)

verify_nestlevel(lvl::Val{:inf}) = lvl
verify_nestlevel(lvl::Integer) = verify_nestlevel(Val(Int(lvl)))
function verify_nestlevel(::Val{n}) where {n}
    n isa Integer ||
        throw(ArgumentError("`nestlevel` must be an integer, `Val` of `Int`, or `Val(:inf)`"))
    lvl = Int(n)
    lvl > 0 || throw(ArgumentError("`nestlevel` must be positive"))
    return Val(lvl)
end

_dec_lvl(lvl::Val{:inf}) = lvl
_dec_lvl(::Val{n}) where {n} = Val(n - 1)

use_threads_for_inner_cats(rf, basesize, nestlevel) =
    cats_to_tcats(rf, TCat(basesize), verify_nestlevel(nestlevel))

# TODO: handle `TeeRF` etc?
cats_to_tcats(rf::R_, innermost_tcat, lvl::Val) =
    Reduction(xform(rf), cats_to_tcats(inner(rf), innermost_tcat, lvl))
cats_to_tcats(rf::R_{Union{Cat,TCat}}, innermost_tcat, lvl::Val) =
    if has(inner(rf), Union{Cat,TCat})
        if lvl isa Val{1}
            setxform(rf, innermost_tcat)
        else
            Reduction(TCat(1), cats_to_tcats(inner(rf), innermost_tcat, _dec_lvl(lvl)))
        end
    else
        setxform(rf, innermost_tcat)
    end


text:
"""
Fetch the minimum value for this band.

For file formats that don't know this intrinsically, the minimum supported
value for the data type will generally be returned.

### Parameters
* `success`   pointer to a boolean to use to indicate if the returned value is
a tight minimum or not. May be `NULL` (default).

### Returns
the minimum raster value (excluding no data pixels)
"""
_getrasterminimum(band::GDALRasterBandH, success::Ptr{Cint}) =
    GDALGetRasterMinimum(band, success)::Cdouble

function getrasterminimum(band::GDALRasterBandH)
    success = Ref{Cint}()
    raster_min = _getrasterminimum(band, success)
    raster_min, Bool(success[])
end

"""
Fetch the maximum value for this band.

For file formats that don't know this intrinsically, the maximum supported
value for the data type will generally be returned.

### Parameters
* `success`   pointer to a boolean to use to indicate if the returned value is
a tight maximum or not. May be `NULL` (default).

### Returns
the maximum raster value (excluding no data pixels)
"""
_getrastermaximum(band::GDALRasterBandH, success::Ptr{Cint}) =
    GDALGetRasterMaximum(band, success)::Cdouble

function getrastermaximum(band::GDALRasterBandH)
    success = Ref{Cint}()
    raster_max = _getrastermaximum(band, success)
    raster_max, Bool(success[])
end

"""
Fetch image statistics.

Returns the minimum, maximum, mean and standard deviation of all pixel values
in this band. If approximate statistics are sufficient, the `bApproxOK` flag
can be set to `true` in which case overviews, or a subset of image tiles may be
used in computing the statistics.

If `bForce` is `FALSE` results will only be returned if it can be done quickly
(ie. without scanning the data). If bForce is FALSE and results cannot be
returned efficiently, the method will return `CE_Warning` but no warning will
have been issued. This is a non-standard use of the CE_Warning return value
to indicate "nothing done".

Note that file formats using PAM (Persistent Auxiliary Metadata) services will
generally cache statistics in the .pam file allowing fast fetch after the first
request.

### Parameters
* `bApproxOK`   If `TRUE` statistics may be computed based on overviews or a
subset of all tiles.
* `bForce`      If `FALSE` statistics will only be returned if it can be done
without rescanning the image.
* `pdfMin`      Location to load image minimum (may be `NULL`).
* `pdfMax`      Location to load image maximum (may be `NULL`).
* `pdfMean`     Location to load image mean (may be `NULL`).
* `pdfStdDev`   Location to load image standard deviation (may be `NULL`).

### Returns
`CE_None` on success, `CE_Warning` if no values returned, `CE_Failure` if error
"""
_getrasterstatistics(band::GDALRasterBandH,
                     bApproxOK::Integer,
                     bForce::Integer,
                     pdfMin::Ptr{Cdouble},
                     pdfMax::Ptr{Cdouble},
                     pdfMean::Ptr{Cdouble},
                     pdfStdDev::Ptr{Cdouble}) =
    GDALGetRasterStatistics(band, bApproxOK, bForce, pdfMin, pdfMax,
                            pdfMean, pdfStdDev)::CPLErr

"""
Compute image statistics.

Returns the minimum, maximum, mean and standard deviation of all pixel values
in this band. If approximate statistics are sufficient, the `bApproxOK` flag
can be set to `true` in which case overviews, or a subset of image tiles may be
used in computing the statistics.

Once computed, the statistics will generally be "set" back on the raster band
using `SetStatistics()`.

### Parameters
* `bApproxOK`     If `TRUE` statistics may be computed based on overviews or a
subset of all tiles.
* `pdfMin`        Location to load image minimum (may be `NULL`).
* `pdfMax`        Location to load image maximum (may be `NULL`).
* `pdfMean`       Location to load image mean (may be `NULL`).
* `pdfStdDev`     Location to load image standard deviation (may be `NULL`).
* `pfnProgress`     a function to call to report progress, or `NULL`.
* `pProgressData`   application data to pass to the progress function.

### Returns
`CE_None` on success, or `CE_Failure` if an error occurs or processing is
terminated by the user.
"""
_computerasterstatistics(band::GDALRasterBandH,
                         bApproxOK::Integer,
                         pdfMin::Ptr{Cdouble},
                         pdfMax::Ptr{Cdouble},
                         pdfMean::Ptr{Cdouble},
                         pdfStdDev::Ptr{Cdouble},
                         pfnProgress::GDALProgressFunc,
                         pProgressData::Ptr{Void}) =
    GDALComputeRasterStatistics(band, bApproxOK, pdfMin, pdfMax, pdfMean,
                                pdfStdDev, pfnProgress, pProgressData)::CPLErr

"""
Set statistics on band.

This method can be used to store min/max/mean/standard deviation statistics on
a raster band.

The default implementation stores them as metadata, and will only work on
formats that can save arbitrary metadata. This method cannot detect whether
metadata will be properly saved and so may return CE_None even if the
statistics will never be saved.

### Parameters
* `dfMin`       minimum pixel value.
* `dfMax`       maximum pixel value.
* `dfMean`      mean (average) of all pixel values.
* `dfStdDev`    Standard deviation of all pixel values.

### Returns
`CE_None` on success or `CE_Failure` on failure.
"""
_setrasterstatistics(band::GDALRasterBandH,
                     dfMin::Cdouble,
                     dfMax::Cdouble,
                     dfMean::Cdouble,
                     dfStdDev::Cdouble) =
    GDALSetRasterStatistics(band, dfMin, dfMax, dfMean, dfStdDev)::CPLErr

"""
Compute the min/max values for a band.

Parameters
* `bApproxOK`   If `TRUE` statistics may be computed based on overviews or a
subset of all tiles.
* `adfMinMax`   pointer to an Array (of size 2) for storing the min/max
"""
_computerasterminmax(band::GDALRasterBandH,
                     bApproxOK::Integer,
                     adfMinMax::Ptr{Cdouble}) =
    GDALComputeRasterMinMax(band, bApproxOK, adfMinMax)

"""
Compute raster histogram.

Note that the bucket size is (dfMax-dfMin) / nBuckets.

For example to compute a simple 256 entry histogram of eight bit data, the
following would be suitable. The unusual bounds are to ensure that bucket
boundaries don't fall right on integer values causing possible errors due to
rounding after scaling.

```C

    GUIntBig anHistogram[256];
    poBand->GetHistogram( -0.5, 255.5, 256, anHistogram, FALSE, FALSE,
                          GDALDummyProgress, NULL );
```

Note that setting bApproxOK will generally result in a subsampling of the file,
and will utilize overviews if available. It should generally produce a
representative histogram for the data that is suitable for use in generating
histogram based luts for instance. Generally bApproxOK is much faster than an
exactly computed histogram.

### Parameters
* `dfMin`               the lower bound of the histogram.
* `dfMax`               the upper bound of the histogram.
* `nBuckets`            the number of buckets in panHistogram.
* `panHistogram`        array into which the histogram totals are placed.
* `bIncludeOutOfRange`  if `TRUE` values below the histogram range will mapped
into `panHistogram[0]`, and values above will be mapped into
`panHistogram[nBuckets-1]`. Otherwise out of range values are discarded.
* `bApproxOK`           `TRUE` if an approximate, or incomplete histogram OK.
* `pfnProgress`         function to report progress to completion.
* `pProgressData`       application data to pass to pfnProgress.

### Returns
`CE_None` on success, or `CE_Failure` if something goes wrong.
"""
_getrasterhistogram(hBand::GDALRasterBandH,
                    dfMin::Cdouble,
                    dfMax::Cdouble,
                    nBuckets::Integer,
                    panHistogram::Ptr{Cint},
                    bIncludeOutOfRange::Integer,
                    bApproxOK::Integer,
                    pfnProgress::GDALProgressFunc,
                    pProgressData::Ptr{Void}) =
    GDALGetRasterHistogram(hBand, dfMin, dfMax, nBuckets, panHistogram,
                           bIncludeOutOfRange, bApproxOK, pfnProgress,
                           pProgressData)::CPLErr

"""
Fetch default raster histogram.

The default method in GDALRasterBand will compute a default histogram. This
method is overriden by derived classes (such as GDALPamRasterBand, VRTDataset,
HFADataset...) that may be able to fetch efficiently an existing histogram.

### Parameters
* `pdfMin`          pointer to contain the lower bound of the histogram.
* `pdfMax`          pointer to contain the upper bound of the histogram.
* `pnBuckets`       pointer to contain the number of buckets in *ppanHistogram.
* `ppanHistogram`   array into which the histogram totals are placed. To be
freed with VSIFree
* `bForce`          `TRUE` to force the computation. If `FALSE` and no default
histogram is available, the method will return `CE_Warning`
* `pfnProgress`     function to report progress to completion.
* `pProgressData`   application data to pass to pfnProgress.

### Returns
`CE_None` on success, `CE_Failure` if something goes wrong, or `CE_Warning` if
no default histogram is available.
"""
_getdefaulthistogram(hBand::GDALRasterBandH,
                     pdfMin::Ptr{Cdouble},
                     pdfMax::Ptr{Cdouble},
                     pnBuckets::Ptr{Cint},
                     ppanHistogram::Ptr{Ptr{Cint}},
                     bForce::Integer,
                     pfnProgress::GDALProgressFunc,
                     pProgressData::Ptr{Void}) =
    GDALGetDefaultHistogram(hBand, pdfMin, pdfMax, pnBuckets, ppanHistogram,
                            bForce, pfnProgress, pProgressData)::CPLErr

"""
Set default histogram.

Use `GDALSetRasterHistogramEx()` instead for counts exceeding 2 billion.
"""
_setdefaulthistogram(hBand::GDALRasterBandH,
                     dfMin::Cdouble,
                     dfMax::Cdouble,
                     nBuckets::Integer,
                     panHistogram::Ptr{Cint}) =
    GDALSetDefaultHistogram(hBand, dfMin, dfMax, nBuckets,panHistogram)::CPLErr


text:
using ADatasets
using Test
using ADatasets: sampleindices, subsampletrainvalidation

@testset "test functions for construction of datasets" begin
	@test length(sampleindices([1,2,3,4,5], 0.2)) == 1
	@test length(sampleindices([1,2,3,4,5], 1.2)) == 5
	@test length(sampleindices([1,2,3,4,5], 10)) == 5
	@test length(sampleindices([1,2,3,4,5], 0)) == 0
	@test length(sampleindices([1,2,3,4,5], 1)) == 1
end

@testset "test functions for construction of datasets" begin
	x = (randn(2,10), fill(0,10))
	xtrn, xval = subsampletrainvalidation(x, 0.1,0.8)
	labels = xtrn[2]
	@test length(xtrn[2]) == 8
	@test size(xtrn[1],2) == 8
	@test length(xval[2]) == 2
	@test size(xval[1],2) == 2

	x = (randn(2,20), vcat(fill(0,10),fill(1,10)))
	xtrn, xval = subsampletrainvalidation(x, 0.1,0.8)
	@test length(xtrn[2]) == 9
	@test size(xtrn[1],2) == 9
	@test sum(xtrn[2] .== 1) == 1
	@test length(xval[2]) == 11
	@test size(xval[1],2) == 11
	@test sum(xval[2] .== 1) == 9

	@test length(sampleindices([1,2,3,4,5], 1.2)) == 5
	@test length(sampleindices([1,2,3,4,5], 10)) == 5
	@test length(sampleindices([1,2,3,4,5], 0)) == 0
	@test length(sampleindices([1,2,3,4,5], 1)) == 1
end

@testset "degenerative case" begin
	x = ones(3)
	for y in [[1, 0, 1], [0, 1, 1], [1, 1, 0]]
		fpr, tpr = ADatasets.roccurve(x, y)
		@test all(fpr .== [0.0, 1.0]) && all(tpr .== [0.0, 1.0])
		@test ADatasets.auc(fpr,tpr) == 0.5
	end
end

# Compare the these functions to SciKit Learn functions
using PyCall
@pyimport sklearn.metrics as sm

function compareskauc(labels, ascores)
	pyfpr, pytpr, _ = sm.roc_curve(labels, ascores, drop_intermediate = false)
	pyauc = sm.auc(pyfpr, pytpr)

	fpr, tpr = ADatasets.roccurve(ascores, labels)
	auc = ADatasets.auc(fpr, tpr)

	@test pyauc ≈ auc
end

@testset "Scikit learn auc comparison" begin
	# rand dataset
	for i in 1:10
		counts = rand(1:1000, 2)
		labels = vcat(zeros(counts[1]), ones(counts[2]))
		ascores = rand(size(labels))
		compareskauc(labels, ascores)
	end
end


text:
@testset "minimal.atan2_test" begin
    @test atan(emptyinterval(), emptyinterval()) ⊇ emptyinterval()
    @test atan(emptyinterval(), entireinterval()) ⊇ emptyinterval()
    @test atan(entireinterval(), emptyinterval()) ⊇ emptyinterval()
    @test atan(interval(0.0, 0.0), interval(0.0, 0.0)) ⊇ emptyinterval()
    @test atan(entireinterval(), entireinterval()) ⊇ interval(-0x1.921FB54442D19p1, +0x1.921FB54442D19p1)
    @test atan(interval(0.0, 0.0), interval(-Inf, 0.0)) ⊇ interval(0x1.921FB54442D18p1, 0x1.921FB54442D19p1)
    @test atan(interval(0.0, 0.0), interval(0.0, Inf)) ⊇ interval(0.0, 0.0)
    @test atan(interval(0.0, Inf), interval(0.0, 0.0)) ⊇ interval(0x1.921FB54442D18p0, 0x1.921FB54442D19p0)
    @test atan(interval(-Inf, 0.0), interval(0.0, 0.0)) ⊇ interval(-0x1.921FB54442D19p0, -0x1.921FB54442D18p0)
    @test atan(interval(-0x1p-1022, 0.0), interval(-0x1p-1022, -0x1p-1022)) ⊇ interval(-0x1.921FB54442D19p1, +0x1.921FB54442D19p1)
    @test atan(interval(1.0, 1.0), interval(-1.0, -1.0)) ⊇ interval(0x1.2D97C7F3321D2p1, 0x1.2D97C7F3321D3p1)
    @test atan(interval(1.0, 1.0), interval(1.0, 1.0)) ⊇ interval(0x1.921FB54442D18p-1, 0x1.921FB54442D19p-1)
    @test atan(interval(-1.0, -1.0), interval(1.0, 1.0)) ⊇ interval(-0x1.921FB54442D19p-1, -0x1.921FB54442D18p-1)
    @test atan(interval(-1.0, -1.0), interval(-1.0, -1.0)) ⊇ interval(-0x1.2D97C7F3321D3p1, -0x1.2D97C7F3321D2p1)
    @test atan(interval(-0x1p-1022, 0x1p-1022), interval(-0x1p-1022, -0x1p-1022)) ⊇ interval(-0x1.921FB54442D19p1, +0x1.921FB54442D19p1)
    @test atan(interval(-0x1p-1022, 0x1p-1022), interval(0x1p-1022, 0x1p-1022)) ⊇ interval(-0x1.921FB54442D19p-1, +0x1.921FB54442D19p-1)
    @test atan(interval(-0x1p-1022, -0x1p-1022), interval(-0x1p-1022, 0x1p-1022)) ⊇ interval(-0x1.2D97C7F3321D3p1, -0x1.921FB54442D18p-1)
    @test atan(interval(0x1p-1022, 0x1p-1022), interval(-0x1p-1022, 0x1p-1022)) ⊇ interval(0x1.921FB54442D18p-1, 0x1.2D97C7F3321D3p1)
    @test atan(interval(-2.0, 2.0), interval(-3.0, -1.0)) ⊇ interval(-0x1.921FB54442D19p1, +0x1.921FB54442D19p1)
    @test atan(interval(0.0, 2.0), interval(-3.0, -1.0)) ⊇ interval(0x1.0468A8ACE4DF6p1, 0x1.921FB54442D19p1)
    @test atan(interval(1.0, 3.0), interval(-3.0, -1.0)) ⊇ interval(0x1.E47DF3D0DD4Dp0, 0x1.68F095FDF593Dp1)
    @test atan(interval(1.0, 3.0), interval(-2.0, 0.0)) ⊇ interval(0x1.921FB54442D18p0, 0x1.56C6E7397F5AFp1)
    @test atan(interval(1.0, 3.0), interval(-2.0, 2.0)) ⊇ interval(0x1.DAC670561BB4Fp-2, 0x1.56C6E7397F5AFp1)
    @test atan(interval(1.0, 3.0), interval(0.0, 2.0)) ⊇ interval(0x1.DAC670561BB4Fp-2, 0x1.921FB54442D19p0)
    @test atan(interval(1.0, 3.0), interval(1.0, 3.0)) ⊇ interval(0x1.4978FA3269EE1p-2, 0x1.3FC176B7A856p0)
    @test atan(interval(0.0, 2.0), interval(1.0, 3.0)) ⊇ interval(0x0p0, 0x1.1B6E192EBBE45p0)
    @test atan(interval(-2.0, 2.0), interval(1.0, 3.0)) ⊇ interval(-0x1.1B6E192EBBE45p0, +0x1.1B6E192EBBE45p0)
    @test atan(interval(-2.0, 0.0), interval(1.0, 3.0)) ⊇ interval(-0x1.1B6E192EBBE45p0, 0x0p0)
    @test atan(interval(-3.0, -1.0), interval(1.0, 3.0)) ⊇ interval(-0x1.3FC176B7A856p0, -0x1.4978FA3269EE1p-2)
    @test atan(interval(-3.0, -1.0), interval(0.0, 2.0)) ⊇ interval(-0x1.921FB54442D19p0, -0x1.DAC670561BB4Fp-2)
    @test atan(interval(-3.0, -1.0), interval(-2.0, 2.0)) ⊇ interval(-0x1.56C6E7397F5AFp1, -0x1.DAC670561BB4Fp-2)
    @test atan(interval(-3.0, -1.0), interval(-2.0, 0.0)) ⊇ interval(-0x1.56C6E7397F5AFp1, -0x1.921FB54442D18p0)
    @test atan(interval(-3.0, -1.0), interval(-3.0, -1.0)) ⊇ interval(-0x1.68F095FDF593Dp1, -0x1.E47DF3D0DD4Dp0)
    @test atan(interval(-2.0, 0.0), interval(-3.0, -1.0)) ⊇ interval(-0x1.921FB54442D19p1, +0x1.921FB54442D19p1)
    @test atan(interval(-5.0, 0.0), interval(-5.0, 0.0)) ⊇ interval(-0x1.921FB54442D19p1, +0x1.921FB54442D19p1)
    @test atan(interval(0.0, 5.0), interval(-5.0, 0.0)) ⊇ interval(0x1.921FB54442D18p0, 0x1.921FB54442D19p1)
    @test atan(interval(0.0, 5.0), interval(0.0, 5.0)) ⊇ interval(0x0p0, 0x1.921FB54442D19p0)
    @test atan(interval(-5.0, 0.0), interval(0.0, 5.0)) ⊇ interval(-0x1.921FB54442D19p0, 0x0p0)
end


text:
@with_kw struct rtgM1{A,B,C} <: RTGmodel
    model::A = rtgM1_model()
    protein_lookup = get_protein_lookup(model)
    p::C = init_p(model; idx_hill_coefs=[])
    u::B = init_u(model, protein_lookup)
end

"""
Rtg Model: version 1
"""
function rtgM1_model()
    rtgM1 = @reaction_network  begin
        #=
        Input
        =#
        (1), s→s
        #=
        Modulation Layer
        =#
        (k1md*s, kn1md), Rtg2_ina_c ↔ Rtg2_act_c
        (k2md, kn2md), Rtg2_act_c + Mks ↔ Rtg2Mks_c
        (k3md, kn3md), Bmh + Mks ↔ BmhMks
        #=
        Activation Model:
        =#
        # NLSact: Nuclear Localization Signal (NLS) activation
        (k2n + k1*(BmhMks)/(k1m+BmhMks)), Rtg13_a_c → Rtg13_i_c
        (k2, k2n), Rtg3_i_c ↔ Rtg3_a_c
        # Rtg1Binding: Rtg13 formation
        (k3, kn3), Rtg1_c + Rtg3_a_c ↔ Rtg13_a_c
        (k3, kn3), Rtg1_c + Rtg3_i_c ↔ Rtg13_i_c
        # Translocation
        (k4, kn4), Rtg13_a_c ↔ Rtg13_a_n
        (k4, kn4), Rtg3_a_c ↔ Rtg3_a_n

    end k1md kn1md k2md kn2md k3md kn3md k1 k1m k2 k2n k3 kn3 k4 kn4
    return rtgM1
end


text:
using MambaModels

# Data

globe_toss = Dict{Symbol, Any}(
  :w => [6, 7, 5, 6, 6],
  :n => [9, 9, 9, 9, 9]
)
globe_toss[:N] = length(globe_toss[:w]);

# Model Specification

model = Model(
  w = Stochastic(1,
    (n, p, N) ->
      UnivariateDistribution[Binomial(n[i], p) for i in 1:N],
    false
  ),
  p = Stochastic(() -> Beta(1, 1))
);

# Initial Values

inits = [
  Dict(:w => globe_toss[:w], :n => globe_toss[:n], :p => 0.5),
  Dict(:w => globe_toss[:w], :n => globe_toss[:n], :p => rand(Beta(1, 1)))
];

# Sampling Scheme

scheme = [NUTS(:p)]
setsamplers!(model, scheme);

# MCMC Simulations

chn = mcmc(model, globe_toss, inits, 10000, burnin=2500, thin=1, chains=2);

# Describe draws

describe(chn)

# Convert to MCMCChains.Chains object

chn2 = MCMCChains.Chains(chn.value, String.(chn.names))

# Describe the MCMCChains

MCMCChains.describe(chn2)

# End of `02/m2.1m.jl`