text:
module ScalarBenchmarks

include(joinpath(dirname(@__FILE__), "..", "utils", "RandUtils.jl"))

using .RandUtils
using BenchmarkTools

const SUITE = BenchmarkGroup()

const INTS = (UInt, Int, BigInt)
const FLOATS = (Float32, Float64, BigFloat)
const REALS = (INTS..., FLOATS...)
const COMPS = map(R -> Complex{R}, REALS)
const NUMS = (REALS..., COMPS...)
const BIGNUMS = (BigInt, BigFloat, Complex{BigInt}, Complex{BigFloat})

##############
# predicates #
##############

g = addgroup!(SUITE, "predicate", ["isinteger", "isinf", "isnan", "iseven", "isodd"])

for T in NUMS
    x = one(T)
    tstr = string(T)
    tol = in(T, BIGNUMS) ? 0.40 : 0.25
    g["isequal", tstr]   = @benchmarkable isequal($x, $x) time_tolerance=tol
    g["isinteger", tstr] = @benchmarkable isinteger($x) time_tolerance=tol
    g["isinf", tstr]     = @benchmarkable isinf($x) time_tolerance=tol
    g["isfinite", tstr]  = @benchmarkable isfinite($x) time_tolerance=tol
    g["isnan", tstr]     = @benchmarkable isnan($x) time_tolerance=tol
end

for T in REALS
    x = one(T)
    tstr = string(T)
    tol = in(T, BIGNUMS) ? 0.40 : 0.25
    g["isless", tstr] = @benchmarkable isless($x, $x) time_tolerance=tol
end

for T in INTS
    x = one(T)
    tstr = string(T)
    tol = in(T, BIGNUMS) ? 0.40 : 0.25
    g["iseven", tstr] = @benchmarkable iseven($x) time_tolerance=tol
    g["isodd", tstr]  = @benchmarkable isodd($x) time_tolerance=tol
end

##############
# arithmetic #
##############

arith = addgroup!(SUITE, "arithmetic")
fstmth = addgroup!(SUITE, "fastmath", ["arithmetic"])

for X in NUMS
    x = one(X)
    xstr = string(X)
    isbignum = in(X, BIGNUMS)
    fstmth["add", xstr] = @benchmarkable @fastmath($x + $(copy(x))) time_tolerance=0.40
    fstmth["sub", xstr] = @benchmarkable @fastmath($x - $(copy(x))) time_tolerance=0.40
    fstmth["mul", xstr] = @benchmarkable @fastmath($x * $(copy(x))) time_tolerance=0.40
    fstmth["div", xstr] = @benchmarkable @fastmath($x / $(copy(x))) time_tolerance=0.40
    for Y in NUMS
        y = one(Y)
        ystr = string(Y)
        tol = (X != Y || isbignum) ? 0.50 : 0.25
        arith["add", xstr, ystr] = @benchmarkable +($x, $y) time_tolerance=tol
        arith["sub", xstr, ystr] = @benchmarkable -($x, $y) time_tolerance=tol
        arith["mul", xstr, ystr] = @benchmarkable *($x, $y) time_tolerance=tol
        arith["div", xstr, ystr] = @benchmarkable /($x, $y) time_tolerance=tol
    end
end

for X in (INTS..., Char, Bool)
    x = X(1) # one(X) is not valid for X==Char
    xstr = string(X)
    for Y in (INTS..., Bool)
        tol = (X != Y || X == BigInt || Y == BigInt) ? 0.40 : 0.25
        arith["rem type", xstr, string(Y)] = @benchmarkable %($x, $Y) time_tolerance=tol
    end
end

function perf_13786(k::Float64,τ::Float64,twopiN::Float64)
    W = getW(k,twopiN)
    y = 1.00-k*τ
    sqy = y>0.0 ? sqrt(y) : 0.0
    tof::Float64 = sqy * (W*y + τ)
    return tof
end

function getW(k::Float64,twopiN::Float64)
    local sqrt2l= sqrt(2.0)
    local K2c0  = 0.7542472332656508
    local K2c1 = -0.2
    local K2c2 = 0.08081220356417687
    local K2c3 = -0.031746031746031744
    local K2c4  = 0.012244273267299524
    local K2c5 = -2.0/429.0
    local K2c6 = 8.0 * sqrt(2.0)/6435.0
    local K2c7 = -8.0/12155.0
    local K2c8 = 8.0*sqrt(2.0)/46189.0
    local k00 =0.3535533905932738
    local k02  =3.00*k00*0.25
    local k03  =-2.00/3.0
    local k04  =15.0*sqrt(2.0)/128.0
    local k05  =-2.0/5.0
    local k06  =35.0*sqrt(2.0)/512.0
    local k07  =-8/35.0
    local k08  =315.0*sqrt(2.0)/8192.0
    local twopi = 2.0*π

    W = 0.0

    @fastmath begin
        if k<=-0.02
            t2 = (k*k)
            t3 = acos(t2 - 1.0)
            t4 = 2.0 - t2
            t6 = twopi - t3   + twopiN
            t5 = 1.0/t4
            W  = (t6 * sqrt(t5) - k)*t5
        elseif (k>=0.02) & (k < (sqrt2l-0.02+twopiN))
            t2 = (k*k)
            t3 = acos(t2 - 1.00)
            t4 = 2.00 - t2
            t6 = (t3  + twopiN)
            t5 = 1.00/t4
            W  = (t6 * sqrt(t5) - k)*t5
        elseif k<0.02 # then!series to k ~= 0 improve convergence
            xt = twopiN+π
            t2 = k*k
            t3 = t2*k
            t4 = t2*t2
            t5 = t4*k
            t6 = t3*t3
            t7 = t5*t2
            t8 = t4*t4
            W  = xt*k00 - k + (xt*k02)*t2 + k03*t3 + (xt*k04)*t4 + k05*t5 + (xt*k06)*t6  + k07*t7 + (xt*k08)*t8
        elseif (k> (sqrt2l-0.02) && k< (sqrt2l+0.02))
            xt = k-sqrt2l
            t2 = xt*xt
            t3 = t2*xt
            t4 = t3*xt
            t5 = t4*xt
            t6 = t5*xt
            t7 = t6*xt
            t8 = t7*xt
            W  = K2c0 + K2c1*k + K2c2*t2 + K2c3*t3 + K2c4*t4 + K2c5*t5 + K2c6*t6 + K2c7*t7 + K2c8*t8
        elseif k >= (sqrt2l+0.02)
            t7 = (k+1.0) * (k-1.00)
            t3 = log(t7 + sqrt(t7*t7-1.0))
            t4 = t7-1.0
            t5 = 1.0/t4
            W  = (-t3 * sqrt(t5) + k)*t5
        end
    end

    return W
end

fstmth["13786"] = @benchmarkable perf_13786(0.4, 0.5, 0.0)

#############
# iteration #
#############

function perf_iterate_indexed(n, v)
    s = 0
    for i = 1:n
        for j = 1:1
            @inbounds k = v[j]
            s += k
        end
    end
    s
end

function perf_iterate_in(n, v)
    s = 0
    for i = 1:n
        for k in v
            s += k
        end
    end
    s
end

g = addgroup!(SUITE, "iteration", ["indexed", "in"])

g["indexed"] = @benchmarkable perf_iterate_indexed(10^5, 3) time_tolerance=0.25
g["in"]      = @benchmarkable perf_iterate_in(10^5, 3) time_tolerance=0.25

#######
# exp #
#######

g = addgroup!(SUITE, "floatexp")

g["ldexp", "norm -> norm",       "Float64"] = @benchmarkable ldexp($(1.7), $(10))
g["ldexp", "norm -> norm",       "Float32"] = @benchmarkable ldexp($(1.7f0), $(10))
g["ldexp", "norm -> subnorm",    "Float64"] = @benchmarkable ldexp($(1.7), $(-1028))
g["ldexp", "norm -> subnorm",    "Float32"] = @benchmarkable ldexp($(1.7f0), $(-128))
g["ldexp", "norm -> inf",        "Float64"] = @benchmarkable ldexp($(1.7), $(1028))
g["ldexp", "norm -> inf",        "Float32"] = @benchmarkable ldexp($(1.7f0), $(128))
g["ldexp", "subnorm -> norm",    "Float64"] = @benchmarkable ldexp($(1.7e-310), $(100))
g["ldexp", "subnorm -> norm",    "Float32"] = @benchmarkable ldexp($(1.7f-40), $(100))
g["ldexp", "subnorm -> subnorm", "Float64"] = @benchmarkable ldexp($(1.7e-310), $(-3))
g["ldexp", "subnorm -> subnorm", "Float32"] = @benchmarkable ldexp($(1.7f-40), $(-3))
g["ldexp", "inf -> inf",         "Float64"] = @benchmarkable ldexp($(Inf), $(10))
g["ldexp", "inf -> inf",         "Float32"] = @benchmarkable ldexp($(Inf32), $(10))

g["frexp", "norm",    "Float64"] = @benchmarkable frexp($(1.7))
g["frexp", "norm",    "Float32"] = @benchmarkable frexp($(1.7f0))
g["frexp", "subnorm", "Float64"] = @benchmarkable frexp($(1.7e-310))
g["frexp", "subnorm", "Float32"] = @benchmarkable frexp($(1.7f-40))
g["frexp", "inf",     "Float64"] = @benchmarkable frexp($(Inf))
g["frexp", "inf",     "Float32"] = @benchmarkable frexp($(Inf32))

g["exponent", "norm",    "Float64"] = @benchmarkable exponent($(1.7))
g["exponent", "norm",    "Float32"] = @benchmarkable exponent($(1.7f0))
g["exponent", "subnorm", "Float64"] = @benchmarkable exponent($(1.7e-310))
g["exponent", "subnorm", "Float32"] = @benchmarkable exponent($(1.7f-40))

g["significand", "norm",    "Float64"] = @benchmarkable significand($(1.7))
g["significand", "norm",    "Float32"] = @benchmarkable significand($(1.7f0))
g["significand", "subnorm", "Float64"] = @benchmarkable significand($(1.7e-310))
g["significand", "subnorm", "Float32"] = @benchmarkable significand($(1.7f-40))

g["exp", "normal path, k = 2",              "Float64"] = @benchmarkable exp($(1.5))
g["exp", "fast path, k = 1",                "Float64"] = @benchmarkable exp($(0.5))
g["exp", "no agument reduction, k = 9",     "Float64"] = @benchmarkable exp($(0.1))
g["exp", "small argument path",             "Float64"] = @benchmarkable exp($(2.0^-30))
g["exp", "normal path -> small, k = -1045", "Float64"] = @benchmarkable exp($(-724.0))
g["exp", "overflow",                        "Float64"] = @benchmarkable exp($(900.0))
g["exp", "underflow",                       "Float64"] = @benchmarkable exp($(-900.0))

g["exp", "normal path, k = 2",              "Float32"] = @benchmarkable exp($(1f5))
g["exp", "fast path, k = 1",                "Float32"] = @benchmarkable exp($(0f5))
g["exp", "no agument reduction, k = 9",     "Float32"] = @benchmarkable exp($(0f1))
g["exp", "small argument path",             "Float32"] = @benchmarkable exp($(2f0^-15))
g["exp", "normal path -> small, k = -1045", "Float32"] = @benchmarkable exp($(-724f0))
g["exp", "overflow",                        "Float32"] = @benchmarkable exp($(150f0))
g["exp", "underflow",                       "Float32"] = @benchmarkable exp($(-150f0))

g["exp10", "no agument reduction, k = 1",     "Float64"] = @benchmarkable exp10($(0.25))
g["exp10", "direct approx, k = 0",            "Float64"] = @benchmarkable exp10($(0.01))
g["exp10", "taylor expansion",                "Float32"] = @benchmarkable exp10($(2f0^-35))
g["exp10", "agument reduction, k = 2",        "Float64"] = @benchmarkable exp10($(0.5))
g["exp10", "agument reduction, k = 83",       "Float64"] = @benchmarkable exp10($(25.0))
g["exp10", "normal path -> small, k = -1075", "Float64"] = @benchmarkable exp10($(-323.6))
g["exp10", "overflow",                        "Float64"] = @benchmarkable exp10($(400.0))
g["exp10", "underflow",                       "Float64"] = @benchmarkable exp10($(-400.0))

g["exp10", "no agument reduction, k = 1",    "Float32"] = @benchmarkable exp10($(0.25f0))
g["exp10", "direct approx, k = 0",           "Float32"] = @benchmarkable exp10($(0.01f0))
g["exp10", "taylor expansion",               "Float32"] = @benchmarkable exp10($(2f0^-25))
g["exp10", "agument reduction, k = 2",       "Float32"] = @benchmarkable exp10($(0.5f0))
g["exp10", "agument reduction, k = 83",      "Float32"] = @benchmarkable exp10($(25.0f0))
g["exp10", "normal path -> small, k = -150", "Float32"] = @benchmarkable exp10($(-45.5f0))
g["exp10", "overflow",                       "Float32"] = @benchmarkable exp10($(100f0))
g["exp10", "underflow",                      "Float32"] = @benchmarkable exp10($(-100f0))

for b in values(g)
    b.params.time_tolerance = 0.40
end

##############
# intfuncs   #
##############

g = addgroup!(SUITE, "intfuncs", ["prevpow2", "nextpow2"])

for T in INTS
    x = T[1, 2, 3, 4, 10, 100, 1024, 10000, 2^30, 2^30-1]
    if T == BigInt
        push!(x, big(2)^3000, big(2)^3000-1)
    end
    y = similar(x)
    tol = in(T, BIGNUMS) ? 0.40 : 0.25
    tstr = string(T)
    for funpow2 = (x -> prevpow(2, x), x -> nextpow(2, x)), sgn = (+,)
        g[string(funpow2), tstr, string(sgn)] = @benchmarkable map!($funpow2, $y, $(sgn(x))) time_tolerance=tol
    end
end

##########
# mod2pi #
##########
g = addgroup!(SUITE, "mod2pi")

# -π/4 <= x <= π/4
g["no reduction", "zero", "Float64"] = @benchmarkable mod2pi($(0.0))
g["no reduction", "positive argument", "Float64"] = @benchmarkable mod2pi($(pi/6))
g["no reduction", "negative argument", "Float64"] = @benchmarkable mod2pi($(-pi/6))
# -2π/4 <= x <= 2π/4
g["no reduction", "positive argument", "Float64"] = @benchmarkable mod2pi($(2*pi/4-0.1))
g["no reduction", "negative argument", "Float64"] = @benchmarkable mod2pi($(-2*pi/4+0.1))
g["no reduction", "positive argument", "Float64"] = @benchmarkable mod2pi($(2*pi/4))
g["no reduction", "negative argument", "Float64"] = @benchmarkable mod2pi($(-2*pi/4))
# -3π/4 <= x <= 3π/4
g["no reduction", "positive argument", "Float64"] = @benchmarkable mod2pi($(3*pi/4-0.1))
g["no reduction", "negative argument", "Float64"] = @benchmarkable mod2pi($(-3*pi/4+0.1))
# -4π/4 <= x <= 4π/4
g["no reduction", "positive argument", "Float64"] = @benchmarkable mod2pi($(pi-0.1))
g["no reduction", "negative argument", "Float64"] = @benchmarkable mod2pi($(-pi+0.1))
g["no reduction", "positive argument", "Float64"] = @benchmarkable mod2pi($(Float64(pi)))
g["no reduction", "negative argument", "Float64"] = @benchmarkable mod2pi($(Float64(-pi)))
# -5π/4 <= x <= 5π/4
g["argument reduction (easy) abs(x) < 5π/4", "positive argument", "Float64"] = @benchmarkable mod2pi($(5*pi/4-0.1))
g["argument reduction (easy) abs(x) < 5π/4", "negative argument", "Float64"] = @benchmarkable mod2pi($(-5*pi/4+0.1))
# -6π/4 <= x <= 6π/4
g["argument reduction (easy) abs(x) < 6π/4", "positive argument", "Float64"] = @benchmarkable mod2pi($(6*pi/4-0.1))
g["argument reduction (easy) abs(x) < 6π/4", "negative argument", "Float64"] = @benchmarkable mod2pi($(-6*pi/4+0.1))
g["argument reduction (hard) abs(x) < 6π/4", "positive argument", "Float64"] = @benchmarkable mod2pi($(6*pi/4))
g["argument reduction (hard) abs(x) < 6π/4", "negative argument", "Float64"] = @benchmarkable mod2pi($(-6*pi/4))
# -7π/4 <= x <= 7π/4
g["argument reduction (easy) abs(x) < 7π/4", "positive argument", "Float64"] = @benchmarkable mod2pi($(7*pi/4-0.1))
g["argument reduction (easy) abs(x) < 7π/4", "negative argument", "Float64"] = @benchmarkable mod2pi($(-7*pi/4+0.1))
# -8π/4 <= x <= 8π/4
g["argument reduction (easy) abs(x) < 8π/4", "positive argument", "Float64"] = @benchmarkable mod2pi($(2*pi-0.1))
g["argument reduction (easy) abs(x) < 8π/4", "negative argument", "Float64"] = @benchmarkable mod2pi($(-2*pi+0.1))
g["argument reduction (hard) abs(x) < 8π/4", "positive argument", "Float64"] = @benchmarkable mod2pi($(2*pi))
g["argument reduction (hard) abs(x) < 8π/4", "negative argument", "Float64"] = @benchmarkable mod2pi($(-2*pi))
# -9π/4 <= x <= 9π/4
g["argument reduction (easy) abs(x) < 9π/4", "positive argument", "Float64"] = @benchmarkable mod2pi($(9*pi/4-0.1))
g["argument reduction (easy) abs(x) < 9π/4", "negative argument", "Float64"] = @benchmarkable mod2pi($(-9*pi/4+0.1))
# -2.0^20π/2 <= x <= 2.0^20π/2
g["argument reduction (easy) abs(x) < 2.0^20π/4", "positive argument", "Float64"] = @benchmarkable mod2pi($(2.0^10*pi/4-0.1))
g["argument reduction (easy) abs(x) < 2.0^20π/4", "negative argument", "Float64"] = @benchmarkable mod2pi($(-2.0^10*pi/4+0.1))
# abs(x) >= 2.0^20π/2
# idx < 0
g["argument reduction (easy) abs(x) > 2.0^20*π/2", "positive argument", "Float64"] = @benchmarkable mod2pi($(2.0^30*pi/4-0.1))
g["argument reduction (easy) abs(x) > 2.0^20*π/2", "negative argument", "Float64"] = @benchmarkable mod2pi($(-2.0^30*pi/4+0.1))
# idx > 0
g["argument reduction (easy) abs(x) > 2.0^20*π/2", "positive argument", "Float64"] = @benchmarkable mod2pi($(2.0^80*pi/4-1.2))
g["argument reduction (easy) abs(x) > 2.0^20*π/2", "negative argument", "Float64"] = @benchmarkable mod2pi($(-2.0^80*pi/4+1.2))

# trig benchmarks

# The benchmark groups below benchmark trig functions in Base. Numeri-
# cal evaluation of trig functions consists of two steps: argument re-
# duction followed by evaluation of a polynomial on the reduced argu-
# ment. Below, "no reduction" means that the kernel functions are cal-
# led directly, "argument reduction (easy)" means that we are using
# the two coefficient Cody-Waite method, "argument reduction (hard)"
# means that we are using a more precise but more expensive Cody-Waite
# scheme, and "argument reduction (paynehanek)" means that we are us-
# ing the expensive Payne-Hanek scheme for large values. "(hard)"
# values are either around integer multiples of pi/2 or for the medium
# size arguments 9pi/4 <= |x| <= 2.0^20π/2. (paynehanek) vales are for
# |x| >= 2.0^20π/2. The tags "sin_kernel" and "cos_kernel" refer to
# the actual polynomial being used. "z_kernel" evaluates a polynomial
# that approximates z∈{sin, cos} on the interval of x's such that
# |x| <= pi/4.

#######
# sin #
#######
arg_string(::Type{Float32}) = "Float32"
arg_string(::Type{Float64}) = "Float64"
g = addgroup!(SUITE, "sin")
for T in (Float32, Float64)
    _arg_string = arg_string(T)
    # -π/4 <= x <= π/4
    g["no reduction", "zero", _arg_string] = @benchmarkable sin($(T(0.0)))
    g["no reduction", "positive argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(T(pi)/6))
    g["no reduction", "negative argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(-T(pi)/6))
    # -2π/4 <= x <= 2π/4
    g["argument reduction (easy) abs(x) < 2π/4", "positive argument", _arg_string, "cos_kernel"] = @benchmarkable sin($(2*T(pi)/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 2π/4", "negative argument", _arg_string, "cos_kernel"] = @benchmarkable sin($(-2*T(pi)/4+T(0.1)))
    g["argument reduction (hard) abs(x) < 2π/4", "positive argument", _arg_string, "cos_kernel"] = @benchmarkable sin($(2*T(pi)/4))
    g["argument reduction (hard) abs(x) < 2π/4", "negative argument", _arg_string, "cos_kernel"] = @benchmarkable sin($(-2*T(pi)/4))
    # -3π/4 <= x <= 3π/4
    g["argument reduction (easy) abs(x) < 3π/4", "positive argument", _arg_string, "cos_kernel"] = @benchmarkable sin($(3*T(pi)/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 3π/4", "negative argument", _arg_string, "cos_kernel"] = @benchmarkable sin($(-3*T(pi)/4+T(0.1)))
    # -4π/4 <= x <= 4π/4
    g["argument reduction (easy) abs(x) < 4π/4", "positive argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(T(pi)-T(0.1)))
    g["argument reduction (easy) abs(x) < 4π/4", "negative argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(-T(pi)+T(0.1)))
    g["argument reduction (hard) abs(x) < 4π/4", "positive argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(T(pi)))
    g["argument reduction (hard) abs(x) < 4π/4", "negative argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(-T(pi)))
    # -5π/4 <= x <= 5π/4
    g["argument reduction (easy) abs(x) < 5π/4", "positive argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(5*T(pi)/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 5π/4", "negative argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(-5*T(pi)/4+T(0.1)))
    # -6π/4 <= x <= 6π/4
    g["argument reduction (easy) abs(x) < 6π/4", "positive argument", _arg_string, "cos_kernel"] = @benchmarkable sin($(6*T(pi)/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 6π/4", "negative argument", _arg_string, "cos_kernel"] = @benchmarkable sin($(-6*T(pi)/4+T(0.1)))
    g["argument reduction (hard) abs(x) < 6π/4", "positive argument", _arg_string, "cos_kernel"] = @benchmarkable sin($(6*T(pi)/4))
    g["argument reduction (hard) abs(x) < 6π/4", "negative argument", _arg_string, "cos_kernel"] = @benchmarkable sin($(-6*T(pi)/4))
    # -7π/4 <= x <= 7π/4
    g["argument reduction (easy) abs(x) < 7π/4", "positive argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(7*T(pi)/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 7π/4", "negative argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(-7*T(pi)/4+T(0.1)))
    # -8π/4 <= x <= 8π/4
    g["argument reduction (easy) abs(x) < 8π/4", "positive argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(2*T(pi)-T(0.1)))
    g["argument reduction (easy) abs(x) < 8π/4", "negative argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(-2*T(pi)+T(0.1)))
    g["argument reduction (hard) abs(x) < 8π/4", "positive argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(2*T(pi)))
    g["argument reduction (hard) abs(x) < 8π/4", "negative argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(-2*T(pi)))
    # -9π/4 <= x <= 9π/4
    g["argument reduction (easy) abs(x) < 9π/4", "positive argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(9*T(pi)/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 9π/4", "negative argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(-9*T(pi)/4+T(0.1)))
    # -2.0^20π/2 <= x <= 2.0^20π/2
    g["argument reduction (easy) abs(x) < 2.0^20π/4", "positive argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(T(2.0)^10*T(pi)/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 2.0^20π/4", "negative argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(-T(2.0)^10*T(pi)/4+T(0.1)))
    # abs(x) >= 2.0^20π/2
    # idx < 0
    g["argument reduction (paynehanek) abs(x) > 2.0^20*π/2", "positive argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(T(2.0)^30*T(pi)/4-T(0.1)))
    g["argument reduction (paynehanek) abs(x) > 2.0^20*π/2", "negative argument", _arg_string, "sin_kernel"] = @benchmarkable sin($(-T(2.0)^30*T(pi)/4+T(0.1)))
    # idx > 0
    g["argument reduction (paynehanek) abs(x) > 2.0^20*π/2", "positive argument", _arg_string, "cos_kernel"] = @benchmarkable sin($(T(2.0)^80*T(pi)/4-1.2))
    g["argument reduction (paynehanek) abs(x) > 2.0^20*π/2", "negative argument", _arg_string, "cos_kernel"] = @benchmarkable sin($(-T(2.0)^80*T(pi)/4+1.2))
end

#######
# cos #
#######

g = addgroup!(SUITE, "cos")
for T in (Float32, Float64)
    _arg_string = arg_string(T)
    # -π/4 <= x <= π/4
    g["no reduction", "zero", _arg_string] = @benchmarkable cos($(0.0))
    g["no reduction", "positive argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(pi/6))
    g["no reduction", "negative argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(-pi/6))
    # -2π/4 <= x <= 2π/4
    g["argument reduction (easy) abs(x) < 2π/4", "positive argument", _arg_string, "sin_kernel"] = @benchmarkable cos($(2*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 2π/4", "negative argument", _arg_string, "sin_kernel"] = @benchmarkable cos($(-2*pi/4+T(0.1)))
    g["argument reduction (hard) abs(x) < 2π/4", "positive argument", _arg_string, "sin_kernel"] = @benchmarkable cos($(2*pi/4))
    g["argument reduction (hard) abs(x) < 2π/4", "negative argument", _arg_string, "sin_kernel"] = @benchmarkable cos($(-2*pi/4))
    # -3π/4 <= x <= 3π/4
    g["argument reduction (easy) abs(x) < 3π/4", "positive argument", _arg_string, "sin_kernel"] = @benchmarkable cos($(3*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 3π/4", "negative argument", _arg_string, "sin_kernel"] = @benchmarkable cos($(-3*pi/4+T(0.1)))
    # -4π/4 <= x <= 4π/4
    g["argument reduction (easy) abs(x) < 4π/4", "positive argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(pi-T(0.1)))
    g["argument reduction (easy) abs(x) < 4π/4", "negative argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(-pi+T(0.1)))
    g["argument reduction (hard) abs(x) < 4π/4", "positive argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(T(pi)))
    g["argument reduction (hard) abs(x) < 4π/4", "negative argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(T(-pi)))
    # -5π/4 <= x <= 5π/4
    g["argument reduction (easy) abs(x) < 5π/4", "positive argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(5*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 5π/4", "negative argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(-5*pi/4+T(0.1)))
    # -6π/4 <= x <= 6π/4
    g["argument reduction (easy) abs(x) < 6π/4", "positive argument", _arg_string, "sin_kernel"] = @benchmarkable cos($(6*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 6π/4", "negative argument", _arg_string, "sin_kernel"] = @benchmarkable cos($(-6*pi/4+T(0.1)))
    g["argument reduction (hard) abs(x) < 6π/4", "positive argument", _arg_string, "sin_kernel"] = @benchmarkable cos($(6*pi/4))
    g["argument reduction (hard) abs(x) < 6π/4", "negative argument", _arg_string, "sin_kernel"] = @benchmarkable cos($(-6*pi/4))
    # -7π/4 <= x <= 7π/4
    g["argument reduction (easy) abs(x) < 7π/4", "positive argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(7*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 7π/4", "negative argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(-7*pi/4+T(0.1)))
    # -8π/4 <= x <= 8π/4
    g["argument reduction (easy) abs(x) < 8π/4", "positive argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(2*pi-T(0.1)))
    g["argument reduction (easy) abs(x) < 8π/4", "negative argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(-2*pi+T(0.1)))
    g["argument reduction (hard) abs(x) < 8π/4", "positive argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(2*pi))
    g["argument reduction (hard) abs(x) < 8π/4", "negative argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(-2*pi))
    # -9π/4 <= x <= 9π/4
    g["argument reduction (easy) abs(x) < 9π/4", "positive argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(9*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 9π/4", "negative argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(-9*pi/4+T(0.1)))
    # -2.0^20π/2 <= x <= 2.0^20π/2
    g["argument reduction (easy) abs(x) < 2.0^20π/4", "positive argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(T(2.0)^10*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 2.0^20π/4", "negative argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(-T(2.0)^10*pi/4+T(0.1)))
    # abs(x) >= 2.0^20π/2
    # idx < 0
    g["argument reduction (paynehanek) abs(x) > 2.0^20*π/2", "positive argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(T(2.0)^30*pi/4-T(0.1)))
    g["argument reduction (paynehanek) abs(x) > 2.0^20*π/2", "negative argument", _arg_string, "cos_kernel"] = @benchmarkable cos($(-T(2.0)^30*pi/4+T(0.1)))
    # idx > 0
    g["argument reduction (paynehanek) abs(x) > 2.0^20*π/2", "positive argument", _arg_string, "sin_kernel"] = @benchmarkable cos($(T(2.0)^80*pi/4-1.2))
    g["argument reduction (paynehanek) abs(x) > 2.0^20*π/2", "negative argument", _arg_string, "sin_kernel"] = @benchmarkable cos($(-T(2.0)^80*pi/4+1.2))
end

##########
# sincos #
##########
g = addgroup!(SUITE, "sincos")
for T in (Float32, Float64)
    _arg_string = arg_string(T)
    # -π/4 <= x <= π/4
    g["no reduction", "zero", _arg_string] = @benchmarkable sincos($(0.0))
    g["no reduction", "positive argument", _arg_string] = @benchmarkable sincos($(pi/6))
    g["no reduction", "negative argument", _arg_string] = @benchmarkable sincos($(-pi/6))
    # -2π/4 <= x <= 2π/4
    g["argument reduction (easy) abs(x) < 2π/4", "positive argument", _arg_string] = @benchmarkable sincos($(2*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 2π/4", "negative argument", _arg_string] = @benchmarkable sincos($(-2*pi/4+T(0.1)))
    g["argument reduction (hard) abs(x) < 2π/4", "positive argument", _arg_string] = @benchmarkable sincos($(2*pi/4))
    g["argument reduction (hard) abs(x) < 2π/4", "negative argument", _arg_string] = @benchmarkable sincos($(-2*pi/4))
    # -3π/4 <= x <= 3π/4
    g["argument reduction (easy) abs(x) < 3π/4", "positive argument", _arg_string] = @benchmarkable sincos($(3*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 3π/4", "negative argument", _arg_string] = @benchmarkable sincos($(-3*pi/4+T(0.1)))
    # -4π/4 <= x <= 4π/4
    g["argument reduction (easy) abs(x) < 4π/4", "positive argument", _arg_string] = @benchmarkable sincos($(pi-T(0.1)))
    g["argument reduction (easy) abs(x) < 4π/4", "negative argument", _arg_string] = @benchmarkable sincos($(-pi+T(0.1)))
    g["argument reduction (hard) abs(x) < 4π/4", "positive argument", _arg_string] = @benchmarkable sincos($(T(pi)))
    g["argument reduction (hard) abs(x) < 4π/4", "negative argument", _arg_string] = @benchmarkable sincos($(T(-pi)))
    # -5π/4 <= x <= 5π/4
    g["argument reduction (easy) abs(x) < 5π/4", "positive argument", _arg_string] = @benchmarkable sincos($(5*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 5π/4", "negative argument", _arg_string] = @benchmarkable sincos($(-5*pi/4+T(0.1)))
    # -6π/4 <= x <= 6π/4
    g["argument reduction (easy) abs(x) < 6π/4", "positive argument", _arg_string] = @benchmarkable sincos($(6*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 6π/4", "negative argument", _arg_string] = @benchmarkable sincos($(-6*pi/4+T(0.1)))
    g["argument reduction (hard) abs(x) < 6π/4", "positive argument", _arg_string] = @benchmarkable sincos($(6*pi/4))
    g["argument reduction (hard) abs(x) < 6π/4", "negative argument", _arg_string] = @benchmarkable sincos($(-6*pi/4))
    # -7π/4 <= x <= 7π/4
    g["argument reduction (easy) abs(x) < 7π/4", "positive argument", _arg_string] = @benchmarkable sincos($(7*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 7π/4", "negative argument", _arg_string] = @benchmarkable sincos($(-7*pi/4+T(0.1)))
    # -8π/4 <= x <= 8π/4
    g["argument reduction (easy) abs(x) < 8π/4", "positive argument", _arg_string] = @benchmarkable sincos($(2*pi-T(0.1)))
    g["argument reduction (easy) abs(x) < 8π/4", "negative argument", _arg_string] = @benchmarkable sincos($(-2*pi+T(0.1)))
    g["argument reduction (hard) abs(x) < 8π/4", "positive argument", _arg_string] = @benchmarkable sincos($(2*pi))
    g["argument reduction (hard) abs(x) < 8π/4", "negative argument", _arg_string] = @benchmarkable sincos($(-2*pi))
    # -9π/4 <= x <= 9π/4
    g["argument reduction (easy) abs(x) < 9π/4", "positive argument", _arg_string] = @benchmarkable sincos($(9*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 9π/4", "negative argument", _arg_string] = @benchmarkable sincos($(-9*pi/4+T(0.1)))
    # -2.0^20π/2 <= x <= 2.0^20π/2
    g["argument reduction (easy) abs(x) < 2.0^20π/4", "positive argument", _arg_string] = @benchmarkable sincos($(T(2.0)^10*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 2.0^20π/4", "negative argument", _arg_string] = @benchmarkable sincos($(-T(2.0)^10*pi/4+T(0.1)))
    # abs(x) >= 2.0^20π/2
    # idx < 0
    g["argument reduction (paynehanek) abs(x) > 2.0^20*π/2", "positive argument", _arg_string] = @benchmarkable sincos($(T(2.0)^30*pi/4-T(0.1)))
    g["argument reduction (paynehanek) abs(x) > 2.0^20*π/2", "negative argument", _arg_string] = @benchmarkable sincos($(-T(2.0)^30*pi/4+T(0.1)))
    # idx > 0
    g["argument reduction (paynehanek) abs(x) > 2.0^20*π/2", "positive argument", _arg_string] = @benchmarkable sincos($(T(2.0)^80*pi/4-1.2))
    g["argument reduction (paynehanek) abs(x) > 2.0^20*π/2", "negative argument", _arg_string] = @benchmarkable sincos($(-T(2.0)^80*pi/4+1.2))
end

########
# tan #
########

g = addgroup!(SUITE, "tan")
for T in (Float32, Float64)
    _arg_string = arg_string(T)
    g["zero", _arg_string] = @benchmarkable tan($(zero(T)))
    g["very small", "positive argument", _arg_string] = @benchmarkable tan($(nextfloat(zero(T))))
    g["very small", "negative argument", _arg_string] = @benchmarkable tan($(prevfloat(zero(T))))
    g["small", "positive argument", _arg_string] = @benchmarkable tan($(prevfloat(sqrt(eps(T))/2)))
    g["small", "negative argument", _arg_string] = @benchmarkable tan($(nextfloat(-sqrt(eps(T))/2)))
    g["medium", "positive argument", _arg_string] = @benchmarkable tan($T(0.6743))
    g["medium", "negative argument", _arg_string] = @benchmarkable tan($T(-0.6743))
    g["large", "positive argument", _arg_string] = @benchmarkable tan($T(0.6745))
    g["large", "negative argument", _arg_string] = @benchmarkable tan($T(-0.6745))
    g["large", "positive argument", _arg_string] = @benchmarkable tan($T(2.5))
    g["large", "negative argument", _arg_string] = @benchmarkable tan($T(-2.5))
end

############
# rem_pio2 #
############

g = addgroup!(SUITE, "rem_pio2")
const _rem = try
    hasmethod(Base.Math.ieee754_rem_pio2, Tuple{Float64})
    Base.Math.ieee754_rem_pio2
catch
    Base.Math.rem_pio2_kernel
end

for T in (Float64, )# (Float32, Float64) add Float32 later
    _arg_string = arg_string(T)
    # -2π/4 <= x <= 2π/4
    g["argument reduction (easy) abs(x) < 2π/4", "positive argument", _arg_string] = @benchmarkable _rem($(2*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 2π/4", "negative argument", _arg_string] = @benchmarkable _rem($(-2*pi/4+T(0.1)))
    g["argument reduction (hard) abs(x) < 2π/4", "positive argument", _arg_string] = @benchmarkable _rem($(2*pi/4))
    g["argument reduction (hard) abs(x) < 2π/4", "negative argument", _arg_string] = @benchmarkable _rem($(-2*pi/4))
    # -3π/4 <= x <= 3π/4
    g["argument reduction (easy) abs(x) < 3π/4", "positive argument", _arg_string] = @benchmarkable _rem($(3*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 3π/4", "negative argument", _arg_string] = @benchmarkable _rem($(-3*pi/4+T(0.1)))
    # -4π/4 <= x <= 4π/4
    g["argument reduction (easy) abs(x) < 4π/4", "positive argument", _arg_string] = @benchmarkable _rem($(pi-T(0.1)))
    g["argument reduction (easy) abs(x) < 4π/4", "negative argument", _arg_string] = @benchmarkable _rem($(-pi+T(0.1)))
    g["argument reduction (hard) abs(x) < 4π/4", "positive argument", _arg_string] = @benchmarkable _rem($(Float64(pi)))
    g["argument reduction (hard) abs(x) < 4π/4", "negative argument", _arg_string] = @benchmarkable _rem($(Float64(-pi)))
    # -5π/4 <= x <= 5π/4
    g["argument reduction (easy) abs(x) < 5π/4", "positive argument", _arg_string] = @benchmarkable _rem($(5*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 5π/4", "negative argument", _arg_string] = @benchmarkable _rem($(-5*pi/4+T(0.1)))
    # -6π/4 <= x <= 6π/4
    g["argument reduction (easy) abs(x) < 6π/4", "positive argument", _arg_string] = @benchmarkable _rem($(6*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 6π/4", "negative argument", _arg_string] = @benchmarkable _rem($(-6*pi/4+T(0.1)))
    g["argument reduction (hard) abs(x) < 6π/4", "positive argument", _arg_string] = @benchmarkable _rem($(6*pi/4))
    g["argument reduction (hard) abs(x) < 6π/4", "negative argument", _arg_string] = @benchmarkable _rem($(-6*pi/4))
    # -7π/4 <= x <= 7π/4
    g["argument reduction (easy) abs(x) < 7π/4", "positive argument", _arg_string] = @benchmarkable _rem($(7*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 7π/4", "negative argument", _arg_string] = @benchmarkable _rem($(-7*pi/4+T(0.1)))
    # -8π/4 <= x <= 8π/4
    g["argument reduction (easy) abs(x) < 8π/4", "positive argument", _arg_string] = @benchmarkable _rem($(2*pi-T(0.1)))
    g["argument reduction (easy) abs(x) < 8π/4", "negative argument", _arg_string] = @benchmarkable _rem($(-2*pi+T(0.1)))
    g["argument reduction (hard) abs(x) < 8π/4", "positive argument", _arg_string] = @benchmarkable _rem($(2*pi))
    g["argument reduction (hard) abs(x) < 8π/4", "negative argument", _arg_string] = @benchmarkable _rem($(-2*pi))
    # -9π/4 <= x <= 9π/4
    g["argument reduction (easy) abs(x) < 9π/4", "positive argument", _arg_string] = @benchmarkable _rem($(9*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 9π/4", "negative argument", _arg_string] = @benchmarkable _rem($(-9*pi/4+T(0.1)))
    # -2.0^20π/2 <= x <= 2.0^20π/2
    g["argument reduction (easy) abs(x) < 2.0^20π/4", "positive argument", _arg_string] = @benchmarkable _rem($(T(2.0)^10*pi/4-T(0.1)))
    g["argument reduction (easy) abs(x) < 2.0^20π/4", "negative argument", _arg_string] = @benchmarkable _rem($(-T(2.0)^10*pi/4+T(0.1)))
    # abs(x) >= 2.0^20π/2
    # idx < 0
    g["argument reduction (paynehanek) abs(x) > 2.0^20*π/2", "positive argument", _arg_string] = @benchmarkable _rem($(T(2.0)^30*pi/4-T(0.1)))
    g["argument reduction (paynehanek) abs(x) > 2.0^20*π/2", "negative argument", _arg_string] = @benchmarkable _rem($(-T(2.0)^30*pi/4+T(0.1)))
    # idx > 0
    g["argument reduction (paynehanek) abs(x) > 2.0^20*π/2", "positive argument", _arg_string] = @benchmarkable _rem($(T(2.0)^80*pi/4-1.2))
    g["argument reduction (paynehanek) abs(x) > 2.0^20*π/2", "negative argument", _arg_string] = @benchmarkable _rem($(-T(2.0)^80*pi/4+1.2))
end

########
# asin #
########

g = addgroup!(SUITE, "asin")
for T in (Float32, Float64)
    _arg_string = arg_string(T)
    g["zero", _arg_string] = @benchmarkable asin($(zero(T)))
    g["small", "positive argument", _arg_string] = @benchmarkable asin($(nextfloat(zero(T))))
    g["small", "negative argument", _arg_string] = @benchmarkable asin($(prevfloat(zero(T))))
    g["one", "positive argument", _arg_string] = @benchmarkable asin($(one(T)))
    g["one", "negative argument", _arg_string] = @benchmarkable asin($(-one(T)))
    g["abs(x) < 0.5", "positive argument", _arg_string] = @benchmarkable asin($(T(0.45)))
    g["abs(x) < 0.5", "negative argument", _arg_string] = @benchmarkable asin($(T(-0.45)))
    g["0.5 <= abs(x) < 0.975", "positive argument", _arg_string] = @benchmarkable asin($(T(0.6)))
    g["0.5 <= abs(x) < 0.975", "negative argument", _arg_string] = @benchmarkable asin($T(-0.6))
    if T == Float64
        g["0.975 <= abs(x) < 1.0", "positive argument", _arg_string] = @benchmarkable asin($(0.98))
        g["0.975 <= abs(x) < 1.0", "negative argument", _arg_string] = @benchmarkable asin($(-0.98))
    end
end

########
# acos #
########

g = addgroup!(SUITE, "acos")
for T in (Float32, Float64)
    _arg_string = arg_string(T)
    g["zero", _arg_string] = @benchmarkable acos($(zero(T)))
    g["small", "positive argument", _arg_string] = @benchmarkable acos($(nextfloat(zero(T))))
    g["small", "negative argument", _arg_string] = @benchmarkable acos($(prevfloat(zero(T))))
    g["one", "positive argument", _arg_string] = @benchmarkable acos($(one(T)))
    g["one", "negative argument", _arg_string] = @benchmarkable acos($(-one(T)))
    g["abs(x) < 0.5", "positive argument", _arg_string] = @benchmarkable acos($(T(0.45)))
    g["abs(x) < 0.5", "negative argument", _arg_string] = @benchmarkable acos($(T(-0.45)))
    g["0.5 <= abs(x) < 1", "positive argument", _arg_string] = @benchmarkable acos($(T(0.6)))
    g["0.5 <= abs(x) < 1", "negative argument", _arg_string] = @benchmarkable acos($T(-0.6))
end

########
# atan #
########

# Before calculating atan(x) a range reduction is performed. The various inter-
# vals below are chosen such that all branches of the reduction and evaluation
# phases are reached and benchmarked.

g = addgroup!(SUITE, "atan")
for T in (Float32, Float64)
    _arg_string = arg_string(T)
    g["zero", _arg_string] = @benchmarkable atan($(zero(T)))
    g["very small", "positive argument", _arg_string] = @benchmarkable atan($(nextfloat(zero(T))))
    g["very small", "negative argument", _arg_string] = @benchmarkable atan($(prevfloat(zero(T))))
    g["very large", "positive argument", _arg_string] = @benchmarkable atan($(T(2.0)^67))
    g["very large", "negative argument", _arg_string] = @benchmarkable atan($(-T(2.0)^67))
    g["0 <= abs(x) < 7/16", "positive argument", _arg_string] = @benchmarkable atan($(T(6/16)))
    g["0 <= abs(x) < 7/16", "negative argument", _arg_string] = @benchmarkable atan($(-T(6/16)))
    g["7/16 <= abs(x) < 11/16", "positive argument", _arg_string] = @benchmarkable atan($(T(10/16)))
    g["7/16 <= abs(x) < 11/16", "negative argument", _arg_string] = @benchmarkable atan($(-T(10/16)))
    g["11/16 <= abs(x) < 19/16", "positive argument", _arg_string] = @benchmarkable atan($(T(18/16)))
    g["11/16 <= abs(x) < 19/16", "negative argument", _arg_string] = @benchmarkable atan($(-T(18/16)))
    g["19/16 <= abs(x) < 39/16", "positive argument", _arg_string] = @benchmarkable atan($(T(38/16)))
    g["19/16 <= abs(x) < 39/16", "negative argument", _arg_string] = @benchmarkable atan($(-T(38/16)))
    g["39/16 <= abs(x) < 2^66", "positive argument", _arg_string] = @benchmarkable atan($(T(50/16)))
    g["39/16 <= abs(x) < 2^66", "negative argument", _arg_string] = @benchmarkable atan($(-T(50/16)))
end

#########
# atan2 #
#########

g = addgroup!(SUITE, "atan2")
for T in (Float32, Float64)
    # when referring to y and x we refer to y and x in atan(y, x)
    _arg_string = arg_string(T)
    r = T(randn())
    absr = abs(r)
    # x is one
    g["x one", _arg_string] = @benchmarkable atan($(T(r)), $(one(T)))
    # y zero
    g["y zero", "y positive", "x positive", _arg_string] = @benchmarkable atan($(T(r)), $(one(T)))
    g["y zero", "y positive", "x positive", _arg_string] = @benchmarkable atan($(zero(T)), $absr)
    g["y zero", "y negative", "x positive", _arg_string] = @benchmarkable atan($(-zero(T)), $absr)
    g["y zero", "y positive", "x negative", _arg_string] = @benchmarkable atan($(zero(T)), $(-absr))
    g["y zero", "y negative", "x negative", _arg_string] = @benchmarkable atan($(-zero(T)), $(-absr))
    # x zero and y not zero
    g["x zero", "y positive", _arg_string] = @benchmarkable atan($(one(T)), $(zero(T)))
    g["x zero", "y negative", _arg_string] = @benchmarkable atan($(-one(T)), $(zero(T)))
    # isinf(x) == true && isinf(y) == true
    g["y infinite", "y positive", "x infinite", "x positive", _arg_string] = @benchmarkable atan($(T(Inf)), $(T(Inf)))
    g["y infinite", "y negative", "x infinite", "x positive", _arg_string] = @benchmarkable atan($(-T(Inf)), $(T(Inf)))
    g["y infinite", "y positive", "x infinite", "x negative", _arg_string] = @benchmarkable atan($(T(Inf)), $(-T(Inf)))
    g["y infinite", "y negative", "x infinite", "x negative", _arg_string] = @benchmarkable atan($(-T(Inf)), $(-T(Inf)))
    # isinf(x) == true && isinf(y) == false
    # m in 0 through 3 are different cases explained in the atan code
    g["y finite", "y positive", "x infinite", "x positive", _arg_string] = @benchmarkable atan($(absr), $(T(Inf))) # m == 0
    g["y finite", "y negative", "x infinite", "x positive", _arg_string] = @benchmarkable atan($(-absr), $(T(Inf))) # m == 1
    g["y finite", "y positive", "x infinite", "x negative", _arg_string] = @benchmarkable atan($(absr), $(-T(Inf))) # m == 2
    g["y finite", "y negative", "x infinite", "x negative", _arg_string] = @benchmarkable atan($(-absr), $(-T(Inf))) # m == 3
    # isinf(y) == true && isinf(x) == false
    g["y infinite", "y positive", "x finite", "x positive", _arg_string] = @benchmarkable atan($(T(Inf)), $(absr))
    g["y infinite", "y negative", "x finite", "x positive", _arg_string] = @benchmarkable atan($(-T(Inf)), $(absr))
    g["y infinite", "y positive", "x finite", "x negative", _arg_string] = @benchmarkable atan($(T(Inf)), $(-absr))
    g["y infinite", "y negative", "x finite", "x negative", _arg_string] = @benchmarkable atan($(-T(Inf)), $(-absr))
    # |y/x| above high threshold
    atanpi = T(1.5707963267948966)
    g["abs(y/x) high", "y positive", "x positive", _arg_string] = @benchmarkable atan($(T(2.0^61)), $(T(1.0)))
    g["abs(y/x) high", "y negative", "x positive", _arg_string] = @benchmarkable atan($(-T(2.0^61)), $(T(1.0)))
    g["abs(y/x) high", "y positive", "x negative", _arg_string] = @benchmarkable atan($(T(2.0^61)), $(-T(1.0)))
    g["abs(y/x) high", "y negative", "x negative", _arg_string] = @benchmarkable atan($(-T(2.0^61)), $(-T(1.0)))
    g["abs(y/x) high", "y infinite", "y negative", "x finite", "x negative", _arg_string] = @benchmarkable atan($(-T(Inf)), $(-absr))
    # |y|/x between 0 and low threshold
    g["abs(y/x) small", "y positive", "x positive", _arg_string] = @benchmarkable atan($(T(2.0^-61)), $(T(1.0)))
    g["abs(y/x) small", "y positive", "x negative", _arg_string] = @benchmarkable atan($(T(2.0^-61)), $(-T(1.0)))
    # y/x is "safe" ("arbitrary values", just need to hit the branch)
    _ATAN2_PI_LO(::Type{Float32}) = -8.7422776573f-08
    _ATAN2_PI_LO(::Type{Float64}) = 1.2246467991473531772E-16
    g["abs(y/x) safe (small)", "y positive", "x positive", _arg_string] = @benchmarkable atan($(T(5.0)), $(T(2.5)))
    g["abs(y/x) safe (small)", "y negative", "x positive", _arg_string] = @benchmarkable atan($(-T(5.0)), $(T(2.5)))
    g["abs(y/x) safe (small)", "y positive", "x negative", _arg_string] = @benchmarkable atan($(T(5.0)), $(-T(2.5)))
    g["abs(y/x) safe (small)", "y negative", "x negative", _arg_string] = @benchmarkable atan($(-T(5.0)), $(-T(2.5)))
    g["abs(y/x) safe (large)", "y positive", "x positive", _arg_string] = @benchmarkable atan($(T(1235.2341234)), $(T(2.5)))
    g["abs(y/x) safe (large)", "y negative", "x positive", _arg_string] = @benchmarkable atan($(-T(1235.2341234)), $(T(2.5)))
    g["abs(y/x) safe (large)", "y positive", "x negative", _arg_string] = @benchmarkable atan($(T(1235.2341234)), $(-T(2.5)))
    g["abs(y/x) safe (large)", "y negative", "x negative", _arg_string] = @benchmarkable atan($(-T(1235.2341234)), $(-T(2.5)))
end

########
# sinh #
########

g = addgroup!(SUITE, "sinh")
for T in (Float32, Float64)
    _arg_string = arg_string(T)
    g["zero", _arg_string] = @benchmarkable sinh($(zero(T)))
    g["very small", "positive argument", _arg_string] = @benchmarkable sinh($(nextfloat(zero(T))))
    g["very small", "negative argument", _arg_string] = @benchmarkable sinh($(prevfloat(zero(T))))
    g["very large", "positive argument", _arg_string] = @benchmarkable sinh($(T(1000)))
    g["very large", "negative argument", _arg_string] = @benchmarkable sinh($(-T(1000)))
end
g["0 <= abs(x) < 2f-12", "positive argument", "Float32"] = @benchmarkable sinh($(prevfloat(2f-12)))
g["0 <= abs(x) < 2f-12", "negative argument", "Float32"] = @benchmarkable sinh($(nextfloat(-2f-12)))
g["2f-12 <= abs(x) < 9f0", "positive argument", "Float32"] = @benchmarkable sinh($(5f0))
g["2f-12 <= abs(x) < 9f0", "negative argument", "Float32"] = @benchmarkable sinh($(-5f0))
g["9f0 <= abs(x) < 88.72283f0", "positive argument", "Float32"] = @benchmarkable sinh($(22f0))
g["9f0 <= abs(x) < 88.72283f0", "negative argument", "Float32"] = @benchmarkable sinh($(-22f0))
g["0 <= abs(x) < 2.0^-28", "positive argument", "Float64"] = @benchmarkable sinh($(prevfloat(2.0^-28)))
g["0 <= abs(x) < 2.0^-28", "negative argument", "Float64"] = @benchmarkable sinh($(nextfloat(-2.0^-28)))
g["2.0^-28 <= abs(x) < 22.0", "positive argument", "Float64"] = @benchmarkable sinh($(5.0))
g["2.0^-28 <= abs(x) < 22.0", "negative argument", "Float64"] = @benchmarkable sinh($(-5.0))
g["22.0 <= abs(x) < 709.7822265633563", "positive argument", "Float64"] = @benchmarkable sinh($(30.0))
g["22.0 <= abs(x) < 709.7822265633563", "negative argument", "Float64"] = @benchmarkable sinh($(-30.0))

########
# cosh #
########

g = addgroup!(SUITE, "cosh")
for T in (Float32, Float64)
    _arg_string = arg_string(T)
    g["zero", _arg_string] = @benchmarkable cosh($(zero(T)))
    g["very small", "positive argument", _arg_string] = @benchmarkable cosh($(nextfloat(zero(T))))
    g["very small", "negative argument", _arg_string] = @benchmarkable cosh($(prevfloat(zero(T))))
    g["very large", "positive argument", _arg_string] = @benchmarkable cosh($(T(1000)))
    g["very large", "negative argument", _arg_string] = @benchmarkable cosh($(-T(1000)))
end
g["0 <= abs(x) < 0.00024414062f0", "positive argument", "Float32"] = @benchmarkable cosh($(prevfloat(0.00024414062f0)))
g["0 <= abs(x) < 0.00024414062f0", "negative argument", "Float32"] = @benchmarkable cosh($(nextfloat(-0.00024414062f0)))
g["0.00024414062f0 <= abs(x) < 9f0", "positive argument", "Float32"] = @benchmarkable cosh($(5f0))
g["0.00024414062f0 <= abs(x) < 9f0", "negative argument", "Float32"] = @benchmarkable cosh($(-5f0))
g["9f0 <= abs(x) < 88.72283f0", "positive argument", "Float32"] = @benchmarkable cosh($(22f0))
g["9f0 <= abs(x) < 88.72283f0", "negative argument", "Float32"] = @benchmarkable cosh($(-22f0))
g["0 <= abs(x) < 2.7755602085408512e-17", "positive argument", "Float64"] = @benchmarkable cosh($(prevfloat(2.7755602085408512e-17)))
g["0 <= abs(x) < 2.7755602085408512e-17", "negative argument", "Float64"] = @benchmarkable cosh($(nextfloat(-2.7755602085408512e-17)))
g["2.7755602085408512e-17 <= abs(x) < 22.0", "positive argument", "Float64"] = @benchmarkable cosh($(5.0))
g["2.7755602085408512e-17 <= abs(x) < 22.0", "negative argument", "Float64"] = @benchmarkable cosh($(-5.0))
g["22.0 <= abs(x) < 709.7822265633563", "positive argument", "Float64"] = @benchmarkable cosh($(30.0))
g["22.0 <= abs(x) < 709.7822265633563", "negative argument", "Float64"] = @benchmarkable cosh($(-30.0))

########
# tanh #
########

g = addgroup!(SUITE, "tanh")
for T in (Float32, Float64)
    _arg_string = arg_string(T)
    g["zero", _arg_string] = @benchmarkable tanh($(zero(T)))
    g["very small", "positive argument", _arg_string] = @benchmarkable tanh($(nextfloat(zero(T))))
    g["very small", "negative argument", _arg_string] = @benchmarkable tanh($(prevfloat(zero(T))))
    g["very large", "positive argument", _arg_string] = @benchmarkable tanh($(T(2.0)^67))
    g["very large", "negative argument", _arg_string] = @benchmarkable tanh($(-T(2.0)^67))
end
g["0 <= abs(x) < 2f0^-12", "positive argument", "Float32"] = @benchmarkable tanh($(prevfloat(2f0^-12)))
g["0 <= abs(x) < 2f0^-12", "negative argument", "Float32"] = @benchmarkable tanh($(nextfloat(-2f0^-12)))
g["2f0^-12 <= abs(x) < 1f0", "positive argument", "Float32"] = @benchmarkable tanh($(0.5f0))
g["2f0^-12 <= abs(x) < 1f0", "negative argument", "Float32"] = @benchmarkable tanh($(-0.5f0))
g["1f0 <= abs(x) < 9f0", "positive argument", "Float32"] = @benchmarkable tanh($(8f0))
g["1f0 <= abs(x) < 9f0", "negative argument", "Float32"] = @benchmarkable tanh($(-8f0))
g["0 <= abs(x) < 2.0^-28", "positive argument", "Float64"] = @benchmarkable tanh($(prevfloat(2.0^-28)))
g["0 <= abs(x) < 2.0^-28", "negative argument", "Float64"] = @benchmarkable tanh($(nextfloat(-2.0^-28)))
g["2.0^-28 <= abs(x) < 1.0", "positive argument", "Float64"] = @benchmarkable tanh($(0.5))
g["2.0^-28 <= abs(x) < 1.0", "negative argument", "Float64"] = @benchmarkable tanh($(-0.5))
g["1.0 <= abs(x) < 22.0", "positive argument", "Float64"] = @benchmarkable tanh($(14.0))
g["1.0 <= abs(x) < 22.0", "negative argument", "Float64"] = @benchmarkable tanh($(-14.0))

#########
# asinh #
#########

g = addgroup!(SUITE, "asinh")
for T in (Float32, Float64)
    _arg_string = arg_string(T)
    g["zero", _arg_string] = @benchmarkable asinh($(zero(T)))
    g["very small", "positive argument", _arg_string] = @benchmarkable asinh($(nextfloat(zero(T))))
    g["very small", "negative argument", _arg_string] = @benchmarkable asinh($(prevfloat(zero(T))))
    g["very large", "positive argument", _arg_string] = @benchmarkable asinh($(T(2)^28))
    g["very large", "negative argument", _arg_string] = @benchmarkable asinh($(-T(2)^28))
    g["0 <= abs(x) < 2^-28", "positive argument", _arg_string] = @benchmarkable asinh($(T(2)^-28))
    g["0 <= abs(x) < 2^-28", "negative argument", _arg_string] = @benchmarkable asinh($(-T(2)^-28))
    g["2^-28 <= abs(x) < 2", "positive argument", _arg_string] = @benchmarkable asinh($(T(1.5)))
    g["2^-28 <= abs(x) < 2", "negative argument", _arg_string] = @benchmarkable asinh($(-T(1.5)))
    g["2 <= abs(x) < 2^28", "positive argument", _arg_string] = @benchmarkable asinh($(T(1000)))
    g["2 <= abs(x) < 2^28", "negative argument", _arg_string] = @benchmarkable asinh($(-T(1000)))
end

#########
# acosh #
#########

g = addgroup!(SUITE, "acosh")
for T in (Float32, Float64)
    _arg_string = arg_string(T)
    g["one", _arg_string] = @benchmarkable acosh($(one(T)))
    g["very large", "positive argument", _arg_string] = @benchmarkable acosh($(T(2.0)^28))
    g["1 <= abs(x) < 2", "positive argument", _arg_string] = @benchmarkable acosh($(T(1.5)))
    g["2 <= abs(x) < 2^28", "positive argument", _arg_string] = @benchmarkable acosh($(T(1000)))
end

#########
# atanh #
#########

g = addgroup!(SUITE, "atanh")
for T in (Float32, Float64)
    _arg_string = arg_string(T)
    g["zero", _arg_string] = @benchmarkable atanh($(zero(T)))
    g["one", _arg_string] = @benchmarkable atanh($(one(T)))
    g["very small", "positive argument", _arg_string] = @benchmarkable atanh($(nextfloat(zero(T))))
    g["very small", "negative argument", _arg_string] = @benchmarkable atanh($(prevfloat(zero(T))))
    g["2^-28 <= abs(x) < 0.5", "positive argument", _arg_string] = @benchmarkable atanh($(T(0.25)))
    g["2^-28 <= abs(x) < 0.5", "negative argument", _arg_string] = @benchmarkable atanh($(-T(0.25)))
    g["0.5 <= abs(x) < 1", "positive argument", _arg_string] = @benchmarkable atanh($(T(0.75)))
    g["0.5 <= abs(x) < 1", "negative argument", _arg_string] = @benchmarkable atanh($(-T(0.75)))
end

#########
# cbrt #
#########

g = addgroup!(SUITE, "cbrt")
for T in (Float32, Float64)
    _arg_string = arg_string(T)
    g["zero", _arg_string] = @benchmarkable cbrt($(zero(T)))
    g["one", _arg_string] = @benchmarkable cbrt($(one(T)))
    g["very small", "positive argument", _arg_string] = @benchmarkable cbrt($(nextfloat(zero(T))))
    g["very small", "negative argument", _arg_string] = @benchmarkable cbrt($(prevfloat(zero(T))))
    g["small", "positive argument", _arg_string] = @benchmarkable cbrt($(T(0.25)))
    g["small", "negative argument", _arg_string] = @benchmarkable cbrt($(-T(0.25)))
    g["medium", "positive argument", _arg_string] = @benchmarkable cbrt($(T(95.75)))
    g["medium", "negative argument", _arg_string] = @benchmarkable cbrt($(-T(95.75)))
    g["large", "positive argument", _arg_string] = @benchmarkable cbrt($(T(3.0)^84))
    g["large", "negative argument", _arg_string] = @benchmarkable cbrt($(-T(3.0)^84))
end
#########
# exp2 #
#########

g = addgroup!(SUITE, "exp2")
for T in (Float32, Float64)
    _arg_string = arg_string(T)
    g["zero", _arg_string] = @benchmarkable exp2($(zero(T)))
    g["one", _arg_string] = @benchmarkable exp2($(one(T)))
    g["very small", "positive argument", _arg_string] = @benchmarkable exp2($(nextfloat(zero(T))))
    g["very small", "negative argument", _arg_string] = @benchmarkable exp2($(prevfloat(zero(T))))
    g["small", "positive argument", _arg_string] = @benchmarkable exp2($(T(0.25)))
    g["small", "negative argument", _arg_string] = @benchmarkable exp2($(-T(0.25)))
    g["2pow3", "positive argument", _arg_string] = @benchmarkable exp2($(T(2)^3))
    g["2pow3", "negative argument", _arg_string] = @benchmarkable exp2($(-T(2)^3))
    g["2pow35", "positive argument", _arg_string] = @benchmarkable exp2($(T(2.0)^35))
    g["2pow35", "negative argument", _arg_string] = @benchmarkable exp2($(-T(2.0)^35))
end
g["2pow1023", "positive argument", Float64] = @benchmarkable exp2($((2.0)^1023))
g["2pow1023", "negative argument", Float64] = @benchmarkable exp2($(-(2.0)^1023))
g["2pow127", "positive argument", Float32] = @benchmarkable exp2($((2.0f0)^127))
g["2pow127", "negative argument", Float32] = @benchmarkable exp2($(-(2.0f0)^127))

#########
# expm1 #
#########

g = addgroup!(SUITE, "expm1")
for T in (Float32, Float64)
    _arg_string = arg_string(T)
    g["zero", _arg_string] = @benchmarkable expm1($(zero(T)))
    g["one", _arg_string] = @benchmarkable expm1($(one(T)))
    g["very small", "positive argument", _arg_string] = @benchmarkable expm1($(nextfloat(zero(T))))
    g["very small", "negative argument", _arg_string] = @benchmarkable expm1($(prevfloat(zero(T))))
    g["small", "positive argument", _arg_string] = @benchmarkable expm1($(T(0.25)))
    g["small", "negative argument", _arg_string] = @benchmarkable expm1($(-T(0.25)))
    g["medium", "positive argument", _arg_string] = @benchmarkable expm1($(T(95.75)))
    g["medium", "negative argument", _arg_string] = @benchmarkable expm1($(-T(95.75)))
    g["large", "positive argument", _arg_string] = @benchmarkable expm1($(T(2.0)^20))
    g["large", "negative argument", _arg_string] = @benchmarkable expm1($(-T(2.0)^20))
    g["large", "positive argument", _arg_string] = @benchmarkable expm1($(T(2.0)^40))
    g["large", "negative argument", _arg_string] = @benchmarkable expm1($(-T(2.0)^40))
end

g["huge", "positive argument", "Float64"] = @benchmarkable expm1($(-56.0*0.6931471805599453-10))
g["arg reduction I", "positive argument", "Float64"] = @benchmarkable expm1($(0.5*0.6931471805599453+0.1))
g["arg reduction I", "negative argument", "Float64"] = @benchmarkable expm1($(0.5*0.6931471805599453+0.1))
g["arg reduction II", "positive argument", "Float64"] = @benchmarkable expm1($(1.5*0.6931471805599453-0.1))
g["arg reduction II", "negative argument", "Float64"] = @benchmarkable expm1($(1.5*0.6931471805599453-0.1))
g["huge", "positive argument", "Float3"] = @benchmarkable expm1($(-18.714973f0-10f0))
g["arg reduction I", "positive argument", "Float32"] = @benchmarkable expm1($(0.5f0*0.6931472f0+0.1f0))
g["arg reduction I", "negative argument", "Float32"] = @benchmarkable expm1($(0.5f0*0.6931472f0+0.1f0))
g["arg reduction II", "positive argument", "Float32"] = @benchmarkable expm1($(1.5f0*0.6931472f0-0.1f0))
g["arg reduction II", "negative argument", "Float32"] = @benchmarkable expm1($(1.5f0*0.6931472f0-0.1f0))

end # module


text:
module LabNotebook

include("setup.jl")
include("maintenance.jl")

end # module

ext: jl


text:
# ```@meta
# CurrentModule = ProximalAlgorithms
# ```
# # Getting started
#
# The methods implemented in ProximalAlgorithms are commonly referred to as (you've guessed it) *proximal algorithms*,
# in that they rely on the *proximal operator* (or *mapping*) to deal with non-differentiable terms in the objective.
# Loosely speaking, the algorithms in this package can be used to solve problems of the form
#
# ```math
# \operatorname*{minimize}_x\ \sum_{i=1}^N f_i(x)
# ```
#
# where ``N`` depends on the specific algorithm, together with specific assumptions on the terms ``f_i`` (like smoothness, convexity, strong convexity).
# The problem above is solved by iteratively accessing specific *first order information* on the terms ``f_i``,
# like their gradient ``\nabla f_i`` or their proximal mapping ``\operatorname{prox}_{f_i}``:
# ```math
# \mathrm{prox}_{\gamma f_i}(x) = \arg\min_z \left\{ f_i(z) + \tfrac{1}{2\gamma}\|z-x\|^2 \right\}
# ```
# The literature on proximal operators and algorithms is vast: for an overview, one can refer to [Parikh2014](@cite), [Beck2017](@cite).
#
# To evaluate these first-order primitives, in ProximalAlgorithms:
# * ``\nabla f_i`` falls back to using automatic differentiation (as provided by [Zygote](https://github.com/FluxML/Zygote.jl)).
# * ``\operatorname{prox}_{f_i}`` relies on the intereface of [ProximalOperators](https://github.com/JuliaFirstOrder/ProximalOperators.jl).
# Both of the above can be implemented for custom function types, as [documented here](@ref custom_terms).
#
#md # !!! note
#md #
#md #     Each of the [implemented algorithms](@ref problems_algorithms) assumes a different structure of the objective to be optimized
#md #     (e.g., a specific number ``N`` of terms), with specific assumptions (e.g., smoothness, convexity).
#md #     Furthermore, multiple algorithms can often be applied to the same problem (possibly through appropriate reformulation) and are expected
#md #     to perform differently; sometimes, even the *same algorithm* can be applied in multiple ways to the same problem, by grouping (*splitting*)
#md #     the terms in different ways.
#md #
#md #     Because of these reasons, ProximalAlgorithms **does not offer a modeling language** to automagically minimize any objective:
#md #     rather, the user is expected to formulate their problem by providing the right objective terms to the algorithm of choice.
#md #     Please refer to the [this section of the manual](@ref problems_algorithms) for information on what terms can be provided and under which assumptions.
#
# ## [Interface to algorithms](@id algorithm_interface)
#
# At a high level, using algorithms from ProximalAlgorithms amounts to the following.
# 1. Instantiate the algorithm, with options like the termination tolerance, verbosity level, or other algorithm-specific parameters.
# 2. Call the algorithm on the problem description: this amounts to the initial point, the objective terms, and possibly additional required information (e.g. Lipschitz constants).
#
# See [here](@ref problems_algorithms) for the list of available algorithm constructors, for different types of problems.
# In general however, algorithms are instances of the [`IterativeAlgorithm`](@ref) type.
#
# ## [Example: box constrained quadratic](@id box_qp)
#
# As a simple example, consider the minimization of a 2D quadratic function subject to box constraints,
# which we will solve using the fast proximal gradient method (also known as fast forward-backward splitting):

using LinearAlgebra
using ProximalOperators
using ProximalAlgorithms

quadratic_cost(x) = dot([3.4 1.2; 1.2 4.5] * x, x) / 2 + dot([-2.3, 9.9], x)
box_indicator = ProximalOperators.IndBox(0, 1)

ffb = ProximalAlgorithms.FastForwardBackward(maxit=1000, tol=1e-5, verbose=true)

# Here, we defined the cost function `quadratic_cost`, and the constraint indicator `box_indicator`.
# Then we set up the optimization algorithm of choice, [`FastForwardBackward`](@ref),
# with options for the maximum number of iterations, termination tolerance, verbosity.
# Finally, we run the algorithm by providing an initial point and the objective terms defining the problem:

solution, iterations = ffb(x0=ones(2), f=quadratic_cost, g=box_indicator)

# We can verify the correctness of the solution by checking that the negative gradient is orthogonal to the constraints, pointing outwards:

-ProximalAlgorithms.gradient(quadratic_cost, solution)[1]

# Or by plotting the solution against the cost function and constraint:

using Plots

contour(-1:0.1:2, -1:0.1:2, (x,y) -> quadratic_cost([x, y]), fill=true, framestyle=:none, background=nothing)
plot!(Shape([0, 1, 1, 0], [0, 0, 1, 1]), opacity=.5, label="feasible set")
scatter!([solution[1]], [solution[2]], color=:red, markershape=:star5, label="computed solution")

# ## [Iterator interface](@id iterator_interface)
#
# Under the hood, algorithms are implemented in the form of standard Julia iterators:
# constructing such iterator objects directly, and looping over them, allows for more fine-grained control over the termination condition,
# or what information from the iterations get logged.
#
# Each iterator is constructed with the full problem description (objective terms and, if needed, additional information like Lipschitz constats)
# and algorithm options (usually step sizes, and any other parameter or option of the algorithm),
# and produces the sequence of states of the algorithm, so that one can do (almost) anything with it.
#
#md # !!! note
#md #     Iterators only implement the algorithm iteration logic, and not additional details like stopping criteria.
#md #     As such, iterators usually yield an infinite sequence of states:
#md #     when looping over them, be careful to properly guard the loop with a stopping criterion.
#
#md # !!! warning
#md #
#md #     To save on allocations, most (if not all) algorithms re-use state objects when iterating,
#md #     by updating the state *in place* instead of creating a new one. For this reason:
#md #     - one **should not** mutate the state object in any way, as this may corrupt the algorithm's logic;
#md #     - one **should not** `collect` the sequence of states, since this will result in an array of identical objects.
#md #
#
# Iterator types are named after the algorithm they implement, so the relationship should be obvious:
# * the `ForwardBackward` algorithm uses the `ForwardBackwardIteration` iterator type;
# * the `FastForwardBackward` algorithm uses the `FastForwardBackwardIteration` iterator type;
# * the `DouglasRachford` algorithm uses the `DouglasRachfordIteration` iterator type;
# and so on.
#
# Let's see what this means in terms of the previous example.
#
# ## [Example: box constrained quadratic (cont)](@id box_qp_cont)
#
# Let's solve the problem from the [previous example](@ref box_qp) by directly interacting with the underlying iterator:
# the `FastForwardBackward` algorithm internally uses a [`FastForwardBackwardIteration`](@ref) object.

ffbiter = ProximalAlgorithms.FastForwardBackwardIteration(x0=ones(2), f=quadratic_cost, g=box_indicator)

# We can now perform anything we want throughout the iteration, by just looping over the iterator:
# for example, we can store the sequence of iterates from the algorithm, to later plot them,
# and stop whenever two successive iterates are closer than a given tolerance.

xs = []
for state in ffbiter
    push!(xs, copy(state.x))
    if length(xs) > 1 && norm(xs[end] - xs[end - 1]) / (1 + norm(xs[end])) <= 1e-5
        break
    end
end

contour(-1:0.1:2, -1:0.1:2, (x,y) -> quadratic_cost([x, y]), fill=true, framestyle=:none, background=nothing)
plot!(Shape([0, 1, 1, 0], [0, 0, 1, 1]), opacity=.5, label="feasible set")
plot!([x[1] for x in xs], [x[2] for x in xs], markershape=:circle, label="algorithm trajectory")
scatter!([solution[1]], [solution[2]], color=:red, markershape=:star5, label="computed solution")

#md # !!! note
#md #
#md #     Since each algorithm iterator type has its own logic, it will also have its own dedicated state structure.
#md #     Interacting with the state then requires being familiar with its structure, and with the nature of its attributes.
#md #


text:
function reduce!(
    op::MonoidUnion, w::GBVector, A::GBMatOrTranspose;
    mask = nothing, accum = nothing, desc = nothing
)
    mask, accum = _handlenothings(mask, accum)
    desc === nothing && (desc = DEFAULTDESC)
    A, desc, _ = _handletranspose(A, desc, nothing)
    op = getoperator(op, eltype(w))
    accum = getaccum(accum, eltype(w))
    libgb.GrB_Matrix_reduce_Monoid(w, mask, accum, op, A, desc)
    return w
end

function reduce!(
    op::Function, w::GBVector, A::GBMatOrTranspose;
    mask = nothing, accum = nothing, desc = nothing
)
    #try to find an existing monoid, if not error:
    return reduce!(Monoids.Monoid(op), w, A; mask, accum, desc)
end

function Base.reduce(
    op::MonoidUnion,
    A::GBMatOrTranspose;
    dims = :,
    typeout = nothing,
    init = nothing,
    mask = nothing,
    accum = nothing,
    desc = nothing
)
    mask, accum = _handlenothings(mask, accum)
    desc === nothing && (desc = DEFAULTDESC)
    if typeout === nothing
        typeout = eltype(A)
    end
    if dims == 2
        w = GBVector{typeout}(size(A, 1))
        reduce!(op, w, A; desc, accum, mask)
        return w
    elseif dims == 1
        desc = desc + T0
        w = GBVector{typeout}(size(A, 2))
        reduce!(op, w, A; desc, accum, mask)
        return w
    elseif dims == (1,2) || dims == Colon()
        if init === nothing
            c = Ref{typeout}()
            typec = typeout
        else
            c = Ref(init)
            typec = typeof(init)
        end
        op = getoperator(op, typec)
        A, desc, _ = _handletranspose(A, desc, nothing)
        accum = getaccum(accum, typec)
        libgb.scalarmatreduce[typeout](c, accum, op, A, desc)
        return c[]
    end
end

function Base.reduce(
    op::Function, A::GBMatOrTranspose;
    dims = :, typeout = nothing, init = nothing, mask = nothing, accum = nothing, desc = nothing
)
    #try to find an existing monoid, if not error:
    return reduce(Monoids.Monoid(op), A; mask, accum, desc, dims, typeout, init)
end

function Base.reduce(
    op::MonoidUnion,
    v::GBVector;
    typeout = nothing,
    init = nothing,
    accum = nothing,
    desc = nothing
)
    accum, desc = _handlenothings(accum, desc)
    if typeout === nothing
        typeout = eltype(v)
    end
    if init === nothing
        c = Ref{typeout}()
        typec = typeout
    else
        c = Ref(init)
        typec = typeof(init)
    end
    op = getoperator(op, typec)
    accum = getaccum(accum, typec)
    libgb.scalarvecreduce[typeout](c, accum, op, v, desc)
    return c[]
end

function Base.reduce(
    op::Function, v::GBVector;
    typeout = nothing, init = nothing, accum = nothing, desc = nothing
)
    #try to find an existing monoid, if not error:
    return reduce(Monoids.Monoid(op), v; accum, desc, typeout, init)
end

function Base.reduce(
    ::BinaryUnion,
    ::GBArray;
    dims = 2,
    typeout = nothing,
    init = nothing,
    mask = nothing,
    accum = nothing,
    desc = nothing
)
    throw(ArgumentError("reduce requires a Monoid op."))
end

"""
    reduce(op::Monoid, A::GBMatrix, dims=:; kwargs...)
    reduce(op::Monoid, v::GBVector; kwargs...)

Reduce `A` along dimensions of A with monoid `op`.

# Arguments
- `op::MonoidUnion`: the monoid reducer. This may not be a BinaryOp.
- `A::GBArray`: `GBVector` or optionally transposed `GBMatrix`.
- `dims = :`: Optional dimensions for GBMatrix, may be `1`, `2`, or `:`.

# Keywords
- `typeout`: Optional output type specification. Defaults to `eltype(A)`.
- `init`: Optional initial value.
- `mask::Union{Nothing, GBMatrix} = nothing`: optional mask.
- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation
    where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied.
- `desc = nothing`
"""
reduce

Base.maximum(A::GBArray) = reduce(max, A)
Base.maximum(f::Function, A::GBArray) = reduce(max, map(f, A))

Base.minimum(A::GBArray) = reduce(min, A)
Base.minimum(f::Function, A::GBArray) = reduce(min, map(f, A))


text:
f(x) = x ^4
f2(x) = 4.0 * x ^ 3.0

x0 = 1

for i = 1:40
	@printf("%.4e\n",x0)
	x0 = x0 - f(x0) / f2(x0)
end

max_stars_repo_name: Magikis/Uniwersity

max_issues_repo_name: Magikis/Uniwersity

max_forks_repo_name: Magikis/Uniwersity


text:
module Private

using EEG         # General
using TimeModels  # Kalman
using LsqFit      # Detrend
using Loess       # Detrend
using PyCall      # Cross Spectral Density
using SIUnits, SIUnits.ShortUnits
using ProgressMeter
using DataFrames
using Plots

export
    kalman_filter,
    set_nans,
	model_amplitude
include("statistics/kalman/kalman.jl")

export
    hotelling,
    plot_hotelling
include("statistics/hotelling.jl")

export
    subsample,
    blank
include("reshaping/subsample.jl")

export
    detrend,
    trend
include("reshaping/detrend.jl")

export
    readLFT,
    readSRF,
    readLOC,
    readELP
include("read_write/elp.jl")
include("read_write/srf.jl")
include("read_write/loc.jl")
include("read_write/lft.jl")

export
    plotSRF
include("plotting/srf.jl")

export
    import_leadfield,
    calculate_specific_leadfield,
    calculate_specific_leadfield2
include("types/Leadfield/Leadfield.jl")
include("types/Leadfield/processing.jl")

export
    beamformer_lcmv,
    cross_spectral_density,
    reduce_epochs,
    retain_svd,
    correct_midline
include("source_analysis/lcmv.jl")
include("source_analysis/cpsd.jl")
include("source_analysis/misc.jl")

end # module


text:
using PooledArrays
import PooledArrays.PooledArray
pool = [1:Int64(round(N/K))...]

@time id6_pooled = PooledArray(PooledArrays.RefArray(rand(Int32(1):Int32(round(N/K)), N)), pool)
@time sumby(id6_pooled, v1)


text:
"""
From (InfrastructureSystems.jl)[https://github.com/NREL/InfrastructureSystems.jl/blob/master/test/test_generate_structs.jl].
"""
@testset "Test generated data standardization structs" begin

    descriptor_file = joinpath(SRC_DIR, "descriptors", "standardize_data_structs.json")
    existing_dir = joinpath(SRC_DIR, "parse", "generated")

    output_dir = "tmp-test-generated-structs"

    if isdir(output_dir)
        rm(output_dir; recursive=true)
    end

    mkdir(output_dir)

    IS.generate_structs(descriptor_file, output_dir; print_results=false)

    matched = true
    try
        run(`diff $output_dir $existing_dir`)
    catch(err)
        @error "Generated structs do not match the descriptor file."
        matched = false
    finally
        rm(output_dir; recursive=true)
    end

    @test matched
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
# ----------------------------------------------------------------------
# What are the shared memory limits?
#=
$ ipcs -lm --human

------ Shared Memory Limits --------
max number of segments = 4096
max seg size = 16E
max total shared memory = 16E
min seg size = 1B
=#

# another way
#=
$ sysctl kernel.shmmni kernel.shmall kernel.shmmax
kernel.shmmni = 4096
kernel.shmall = 18446744073692774399
kernel.shmmax = 18446744073692774399
=#

# ----------------------------------------------------------------------
# How to adjust?

#=
$ sudo sysctl -w kernel.shmmax=137438953472

$ sysctl kernel.shmmax
kernel.shmmax = 137438953472
=#

# ----------------------------------------------------------------------
# Current size of /dev/shm

#=
$ df -h
Filesystem      Size  Used Avail Use% Mounted on
devtmpfs         16G     0   16G   0% /dev
tmpfs            16G     0   16G   0% /dev/shm
tmpfs            16G  448K   16G   1% /run
tmpfs            16G     0   16G   0% /sys/fs/cgroup
/dev/nvme0n1p1   50G   30G   21G  59% /
tmpfs           3.1G     0  3.1G   0% /run/user/1000
=#

# ----------------------------------------------------------------------
# Start new julia process and create a SharedArray

#=
julia> using SharedArrays

julia> A = SharedArray{Float64}(10000,10000);

julia> A[:] = rand(10000, 10000);

julia> varinfo(Main, r"A")
name        size summary
–––– ––––––––––– ––––––––––––––––––––––––––––––––––
A    762.940 MiB 10000×10000 SharedArray{Float64,2}
=#

# Check system again
#=
$ df -h | egrep '(Used|shm)'
Filesystem      Size  Used Avail Use% Mounted on
shmfs            16G  763M   16G   5% /dev/shm
=#

# ----------------------------------------------------------------------
# How to add more space
#=
$ sudo mount -t tmpfs shmfs -o size=28g /dev/shm

$ df -h | egrep '(Used|shm)'
Filesystem      Size  Used Avail Use% Mounted on
shmfs            28G     0   28G   0% /dev/shm
=#


text:
"""
DateDataFrame add timestamp to following objects:

df::DataFrame
df[1, :] # DataFrameRow
df[!, :symbol] # Vector{Float64} # pandas.series like
df[1:2, :] # DataFrame
"""
const DateTimeOrVec = Union{AbstractVector{DateTime}, DateTime}

struct DateDataFrame{TT <: DateTimeOrVec, DT <: Union{AbstractDataFrame, DataFrameRow, AbstractVector}}
    timestamp::TT
    df::DT
end

const DateDataFrameVecEnd = DateDataFrame{<:DateTimeOrVec, <:AbstractVector}

function show_timestamp(ts::AbstractVector{DateTime})
    return "$time::$(typeof(ts)) span=> $(ts[1])->$(ts[end])"
end

function show_timestamp(ts::DateTime)
    "date_time=> $ts"
end

function Base.show(io::IO, ddf::DateDataFrame)
    print(io, "DateDataFrame: $(show_timestamp(ddf.timestamp)) extended df=")
    sdf = hcat(DataFrame(Dict(:__date => ddf.timestamp)), ddf.df)
    Base.show(io, sdf)
end

type_translate(timestamp, df) = DateDataFrame(timestamp, df)
type_translate(timestamp, df::Number) = df

function Base.getindex(ddf::DateDataFrame, row_idx, col_idx)
    t_idx = translate_idx(ddf.timestamp, row_idx)
    # @show t_idx
    df = ddf.df[t_idx, col_idx]
    # timstamp = ddf.timestamp[t_idx]
    timestamp = translate_timestamp(ddf.timestamp, t_idx)
    return type_translate(timestamp, df)
end

function Base.getindex(ddf::DateDataFrameVecEnd, row_idx)
    t_idx = translate_idx(ddf.timestamp, row_idx)
    return DateDataFrame(ddf.timestamp[t_idx], ddf.df[t_idx])
end

function Base.getindex(ddf::DateDataFrameVecEnd, row_idx::Union{Int, DateTime})
    t_idx = translate_idx(ddf.timestamp, row_idx)
    return ddf.df[t_idx]
end

function Base.setindex!(ddf::DateDataFrame, value, row_idx, col_idx)
    t_idx = translate_idx(ddf.timestamp, row_idx)
    ddf.df[t_idx, col_idx] = value
end

function Base.setindex!(ddf::DateDataFrame{<:Any, <:AbstractVector}, value, idx)
    t_idx = translate_idx(ddf.timestamp, idx)
    ddf.df[t_idx] = value
end

# General fallback
function translate_idx(ts, row_idx)
    return row_idx
end

function translate_idx(ts::AbstractVector{DateTime}, dt::DateTime)
    idx = searchsorted(ts, dt)
    if length(idx) > 0
        return first(idx)
    else length(idx) == 0
        error("$dt not found in index")
    end
end

function translate_idx(ts::StepRange{DateTime, <:Any}, dt::DateTime)
    d = dt - ts.start
    step = typeof(d)(ts.step)
    # @show ts ts.start d ts.step typeof(d)(ts.step) (d % typeof(d)(ts.step))
    if d % step == zero(typeof(d))
        idx = d ÷ step + 1
        return idx
    end
    error("DateTime $dt is not in index")
end

function translate_idx(ts, dt::AbstractVector{DateTime})
    return translate_idx.([ts], dt)
end

# TODO, add a dispatch for row_idx::StepRange{DateTime, <:Any}

function translate_timestamp(ts::AbstractVector{<:DateTime}, idx)
    return ts[idx]
end

function translate_timestamp(ts::AbstractVector{<:DateTime}, ::typeof(!))
    return ts
end

Base.eachcol(ddf::DateDataFrame) = Base.eachcol(ddf.df)
Base.size(ddf::DateDataFrame, args...) = Base.size(ddf.df, args...)

Base.sum(ddf::DateDataFrameVecEnd) = sum(ddf.df)
Base.first(ddf::DateDataFrame, row_idx=1:1) = ddf[row_idx, :]
Base.first(ddf::DateDataFrameVecEnd, row_idx=1:1) = ddf[row_idx]

Base.cumsum(ddf::DateDataFrameVecEnd) = DateDataFrame(ddf.timestamp, cumsum(ddf.df))

# While it's possible to implement `iterate` interface for DateTimeOrVec, I don't do it to capture other errors easily.

#=
Reference:
https://github.com/JuliaData/DataFrames.jl/blob/main/src/abstractdataframe/abstractdataframe.jl#L400
=#
function Base.:(==)(ddf1::DateDataFrame, ddf2::DateDataFrame)
    return ddf1.timestamp == ddf2.timestamp && ddf1.df == ddf2.df
end

Base.:(-)(dff::DateDataFrameVecEnd) = DateDataFrame(ddf.timestamp, -ddf.df)

function Base.hcat(ddf1::DateDataFrame, ddf2::DateDataFrame)
    @assert ddf1.timestamp == ddf2.timestamp
    return DateDataFrame(ddf1.timestamp, hcat(ddf1.df, ddf2.df))
end

function Base.vcat(ddf1::DateDataFrame, ddf2::DateDataFrame)
    @assert ddf1.timestamp[end] <= ddf2.timestamp[begin]
    return DateDataFrame(vcat(ddf1.timestamp, ddf2.timestamp), vcat(ddf1.df, ddf2.df))
end


text:
import SigDigits.use_sigdigits
import SigDigits.count_sigdigits

@testset "Individual significant digits" begin
    @test count_sigdigits(0.0000007) == 1
    @test count_sigdigits(5.034001) == 7
    @test count_sigdigits(7000) == 1
    @test count_sigdigits(1.000) == 1
    @test count_sigdigits(6.983410) == 6
    @test count_sigdigits(14080.093149) == 11
end

@testset "Simple expression significant digits" begin
    @test use_sigdigits("5.4+3.5") == 8.9
    @test use_sigdigits("5.4-3.5") == 1.9
    @test use_sigdigits("5.4*3.5") == 19
    @test use_sigdigits("5.4/3.5") == 1.5
    @test use_sigdigits("5.402+3.5") == 8.9
    @test use_sigdigits("5.402-3.5") == 1.9
    @test use_sigdigits("5.402*3.5") == 19
    @test use_sigdigits("5.402/3.5") == 1.5
end

@testset "Compound expression significant digits" begin
    @test use_sigdigits("(8.958*1398.002/23.95)/(150.005)") == 3.486
    @test use_sigdigits("(8.958*1398.002/23.95)/(150.005/0.01 - 45.29*-3.59/-0.59)") == 0.04
    @test use_sigdigits("(8.958*1398.002+293.45/23.95)/(150.005/0.01 - 45.29*-3.59/-0.59)") == 0.9
end


text: balancedbrackets(str::AbstractString) = foldl((x, y) -> x < 0 ? -1 : x + y, 0, collect((x == '[') - (x == ']') for x in str)) == 0

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# This is a barebones left-right setup, which can quickly be
# changed to see the phase portrait of any function. However,
# there are no cut tests
using ComplexVisual
@ComplexVisual.import_huge

const CVG = try using ComplexVisualGtk; ComplexVisualGtk; catch nothing; end

trafo = z -> z^2
domain   = CV_Math2DCanvas(-2.0 + 2.0im, 2.0 - 2.0im, 100)
codomain = CV_Math2DCanvas(-4.0 + 4.0im, 4.0 - 4.0im, 50)

lr_start_kwargs = Dict(:z_start => 0.2 +0.8im)
scene = cv_scene_lr_std(trafo, domain, codomain; lr_start_kwargs)

cv_save_image(cv_get_can_layout(scene), "./zsquare_lr.png")

if CVG !== nothing
    handler = CVG.cvg_visualize(scene)
    CVG.cvg_wait_for_destroy(handler.window)

    CVG.cvg_close(handler); CVG.cv_destroy(scene)
end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
"""

The `QuickHeaps` module implements versatile binary heaps and priority queues
for Julia.

Wikipedia page https://en.wikipedia.org/wiki/Binary_heap has very clear
explanations about binary heaps.

This code was much inspired by the DataStructures.jl package (see
https://github.com/JuliaCollections/DataStructures.jl) but has a number of
improvements:

- sorting methods are 2 to 4 times faster (in part because NaN are
  ignored but also because of carefully in-lining critical sections);

- methods not having the `unsafe` word in their name are safer and check
  their arguments for correctness (methods with the `unsafe` word in
  their name are in tended to be as fast as possible).

"""
module QuickHeaps

export
    # Form Base.Order:
    Ordering, ForwardOrdering, ReverseOrdering, Forward, Reverse,

    # From this package:
    FastForward, FastReverse,
    FastMin, FastMax, SafeMin, SafeMax,
    AbstractBinaryHeap, BinaryHeap, FastBinaryHeap,
    AbstractPriorityQueue, PriorityQueue, FastPriorityQueue,
    heapify, heapify!, isheap,

    # From Base (in recent versions of Julia):
    peek,

    # From DataStructures:
    enqueue!, dequeue!, dequeue_pair!

import Base:
    IndexStyle,
    IteratorEltype,
    IteratorSize,
    Pair,
    Tuple,
    copy,
    delete!,
    eltype,
    empty!,
    first,
    get,
    getindex,
    haskey,
    isempty,
    iterate,
    keytype,
    keys,
    length,
    peek,
    pop!,
    push!,
    resize!,
    setindex!,
    show,
    size,
    sizehint!,
    valtype,
    values

# The `peek` method appeared in Julia 1.5.
@static if isdefined(Base, :peek)
    import Base: peek
end

using Base: @propagate_inbounds, OneTo, has_offset_axes, HasEltype, HasLength

import DataStructures:
    # heapify!, heapify, isheap
    enqueue!, dequeue!, dequeue_pair!

using ArrayTools

#------------------------------------------------------------------------------
# In order to perform fast sorting (without taking care of NaN's), we
# extend `Base.Order.lt` method for specialized ordering types.  The
# ordering can be an instance or its type.

using Base.Order: Ordering, ForwardOrdering, ReverseOrdering, Forward, Reverse
import Base.Order: lt

"""
    FastForwardOrdering

is the singleton type for fast *forward* ordering without considering NaN's.

"""
struct FastForwardOrdering <: Ordering end

lt(::FastForwardOrdering, a, b) = a < b

const FastForward = FastForwardOrdering()
const FastReverse = ReverseOrdering(FastForward)

const FastMin = FastForward
const FastMax = FastReverse

const SafeMin = Forward
const SafeMax = Reverse

"""
    default_ordering(T)

yields the default ordering for an ordered data structure of type `T`.  This
method shall be specialized for each ordered data structure.

"""
default_ordering(x::Any) = default_ordering(typeof(x))
default_ordering(::Type) = Forward

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

include("utilities.jl")
include("binaryheaps.jl")
include("priorityqueues.jl")

end # module


text:
module ContinuumWorld

# package code goes here
using Random
using LinearAlgebra
using POMDPs
using StaticArrays
using Parameters
using GridInterpolations
using POMDPModelTools
using POMDPModels
using Plots
plotly()

export
    CWorld,
    CWorldVis,
    CircularRegion,
    Vec2,
    CWorldSolver,

    evaluate,
    action_ind

const Vec2 = SVector{2, Float64}

struct CircularRegion
    center::Vec2
    radius::Float64
end

Base.in(v::Vec2, r::CircularRegion) = LinearAlgebra.norm(v-r.center) <= r.radius

const card_and_stay = [Vec2(1.0, 0.0), Vec2(-1.0, 0.0), Vec2(0.0, 1.0), Vec2(0.0, -1.0), Vec2(0.0, 0.0)]
const cardinal = [Vec2(1.0, 0.0), Vec2(-1.0, 0.0), Vec2(0.0, 1.0), Vec2(0.0, -1.0)]
const default_regions = [CircularRegion(Vec2(3.5, 2.5), 0.5),
                         CircularRegion(Vec2(3.5, 5.5), 0.5),
                         CircularRegion(Vec2(8.5, 2.5), 0.5),
                         CircularRegion(Vec2(7.5, 7.5), 0.5)]
const default_rewards = [-10.0, -5.0, 10.0, 3.0]

@with_kw struct CWorld <: MDP{Vec2, Vec2}
    xlim::Tuple{Float64, Float64}                   = (0.0, 10.0)
    ylim::Tuple{Float64, Float64}                   = (0.0, 10.0)
    reward_regions::Vector{CircularRegion}          = default_regions
    rewards::Vector{Float64}                        = default_rewards
    terminal::Vector{CircularRegion}                = default_regions
    stdev::Float64                                  = 0.5
    actions::Vector{Vec2}                           = cardinal
    discount::Float64                               = 0.95
end

POMDPs.actions(w::CWorld) = w.actions
POMDPs.n_actions(w::CWorld) = length(w.actions)
POMDPs.discount(w::CWorld) = w.discount

function POMDPs.gen(::DDNNode{:sp}, w::CWorld, s::AbstractVector, a::AbstractVector, rng::AbstractRNG)
    return s + a + w.stdev*randn(rng, Vec2)
end

function POMDPs.reward(w::CWorld, s::AbstractVector, a::AbstractVector, sp::AbstractVector) # XXX inefficient
    rew = 0.0
    for (i,r) in enumerate(w.reward_regions)
        if sp in r
            rew += w.rewards[i]
        end
    end
    return rew
end

function POMDPs.isterminal(w::CWorld, s::Vec2) # XXX inefficient
    for r in w.terminal
        if s in r
            return true
        end
    end
    return false
end

function POMDPs.initialstate(w::CWorld, rng::AbstractRNG)
    x = w.xlim[1] + (w.xlim[2] - w.xlim[1]) * rand(rng)
    y = w.ylim[1] + (w.ylim[2] - w.ylim[1]) * rand(rng)
    return Vec2(x,y)
end

include("solver.jl")
include("visualization.jl")

end # module


text:
module GalacticSpeedMapping

using SpeedMapping, GalacticOptim, GalacticOptim.SciMLBase

export SpeedMappingOpt

struct SpeedMappingOpt end

function __map_optimizer_args(prob::OptimizationProblem, opt::SpeedMappingOpt;
    callback=nothing,
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    kwargs...)

    # add optimiser options from kwargs
    mapped_args = (; kwargs...)

    if !(isnothing(maxiters))
        @info "maxiters defines maximum gradient calls for $(opt)"
        mapped_args = (; mapped_args..., maps_limit=maxiters)
    end

    if !(isnothing(maxtime))
        mapped_args = (; mapped_args..., time_limit=maxtime)
    end

    if !isnothing(abstol)
        @warn "common abstol is currently not used by $(opt)"
    end

    if !isnothing(reltol)
        @warn "common reltol is currently not used by $(opt)"
    end

    return mapped_args
end

function SciMLBase.__solve(prob::OptimizationProblem, opt::SpeedMappingOpt;
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    progress=false,
    kwargs...)
    local x

    maxiters = GalacticOptim._check_and_convert_maxiters(maxiters)
    maxtime = GalacticOptim._check_and_convert_maxtime(maxtime)

    f = GalacticOptim.instantiate_function(prob.f, prob.u0, prob.f.adtype, prob.p)

    _loss = function (θ)
        x = f.f(θ, prob.p)
        return first(x)
    end

    if isnothing(f.grad)
        @info "SpeedMapping's ForwardDiff AD backend is used to calculate the gradient information."
    end

    opt_args = __map_optimizer_args(prob, opt, maxiters=maxiters, maxtime=maxtime, abstol=abstol, reltol=reltol; kwargs...)

    t0 = time()
    opt_res = SpeedMapping.speedmapping(prob.u0; f=_loss, (g!)=f.grad, lower=prob.lb, upper=prob.ub, opt_args...)
    t1 = time()
    opt_ret = Symbol(opt_res.converged)

    SciMLBase.build_solution(prob, opt, opt_res.minimizer, _loss(opt_res.minimizer); original=opt_res, retcode=opt_ret)
end

end


text:
using Libdl

savepath = normpath(joinpath(@__DIR__, "..", "src", "lib_lightgbm.$(Libdl.dlext)"))

# These precompiled ones are CPU-only, can incorporate stuff later for user-built GPU capable binaries
GH_URL = "https://github.com/microsoft/LightGBM/releases/download/v3.1.0/lib_lightgbm.$(Libdl.dlext)"

download(GH_URL, savepath)

# The end


text:
using NearestNeighbors
include("../types/types.jl")
include("../setting/constants.jl")
include("./fortran_call.jl")

function sem_mesh_locate_kdtree2!(mesh_data::sem_mesh_data, npoint::Int64, xyz::Array{Float64,2}, idoubling::Vector{Int64}, nnearest::Int64, max_search_dist::Float64, max_misloc::Float64, iproc_old::Int64)
    # * init result
    location_result = Vector{sem_mesh_location}(undef, npoint)
    for i in 1:npoint
        location_result[i] = sem_mesh_location()
    end

    # * GLL colocation points and lagrange interpolation weights
    xigll = zeros(Float64, NGLLX)
    wxgll = zeros(Float64, NGLLX)
    yigll = zeros(Float64, NGLLY)
    wygll = zeros(Float64, NGLLY)
    zigll = zeros(Float64, NGLLZ)
    wzgll = zeros(Float64, NGLLZ)
    hlagx = zeros(Float64, NGLLX)
    hlagy = zeros(Float64, NGLLY)
    hlagz = zeros(Float64, NGLLZ)

    # * mesh dimension variable
    nspec = mesh_data.nspec

    # * coordinates of GLL points
    zwgljd!(xigll, wxgll, NGLLX, GAUSSALPHA, GAUSSBETA)
    zwgljd!(yigll, wygll, NGLLY, GAUSSALPHA, GAUSSBETA)
    zwgljd!(zigll, wzgll, NGLLZ, GAUSSALPHA, GAUSSBETA)

    # * anchor points of mesh elements

    # get index of anchor points in the GLL element
    iax = zeros(Int32, NGNOD)
    iay = zeros(Int32, NGNOD)
    iaz = zeros(Int32, NGNOD)
    anchor_point_index!(iax, iay, iaz)

    # get anchor and center points of each GLL element
    xyz_elem = zeros(Float64, 3, nspec)
    xyz_anchor = zeros(Float64, 3, NGNOD, nspec)

    for ispec = 1:nspec
        for ia = 1:NGNOD
            iglob = mesh_data.ibool[iax[ia], iay[ia], iaz[ia], ispec]
            xyz_anchor[:,ia,ispec] = mesh_data.xyz_glob[:,iglob]
        end
        # the last anchor point is the element center
        xyz_elem[:,ispec] = xyz_anchor[:,NGNOD,ispec]
    end

    # * build kdtree
    kdtree = KDTree(xyz_elem)
    # * locate each point
    # initialize location results
    for i in 1:npoint
        location_result[i].stat = -1
        location_result[i].eid = -1
        location_result[i].misloc = HUGEVAL
        location_result[i].lagrange = zeros(Float64, NGLLX, NGLLY, NGLLZ)
        location_result[i].uvw = zeros(Float64, 3)
    end

    # loop points
    for ipoint = 1:npoint
        xyz1 = xyz[:,ipoint]
       # * get the n nearest elements in the mesh
        idxs, _ = knn(kdtree, xyz1, nnearest)

       # * test each neighbour elements to see if target point is located inside
        for inn = 1:nnearest
            ispec = idxs[inn]
            # skip the element a certain distance away
            dist = sqrt(sum((xyz_elem[:,ispec] .- xyz1).^2))
            # if (ipoint == 492835)
            #     @info "@@@@"
            #     @info dist, max_search_dist, inn
            # end
            if dist > max_search_dist
                continue
            end

            # only use element with the same layer ID (i.e. idoubling)
            if (idoubling[ipoint] != IFLAG_DUMMY && mesh_data.idoubling[ispec] != idoubling[ipoint])
                continue
            end

            # locate point to this element
            uvw1 = similar(xyz1)
            misloc1 = HUGEVAL
            flag_inside = false
            misloc1, flag_inside = xyz2cube_bounded!(xyz_anchor[:,:,ispec], xyz1, uvw1, misloc1, flag_inside)
            # if (ipoint == 492835)
            #     @info "@@@@"
            #     @info flag_inside, misloc1, uvw1, ispec
            # end
            if flag_inside == true
                location_result[ipoint].stat = 1
                location_result[ipoint].eid = ispec
                location_result[ipoint].misloc = misloc1
                location_result[ipoint].uvw = uvw1
                break
            else
                if (misloc1 < max_misloc) && (misloc1 < location_result[ipoint].misloc)
                    location_result[ipoint].stat = 0
                    location_result[ipoint].eid = ispec
                    location_result[ipoint].misloc = misloc1
                    location_result[ipoint].uvw = uvw1
                end
            end
        end

        # * set interpolation weights on GLL points if located
        if location_result[ipoint].stat != -1
            lagrange_poly!(location_result[ipoint].uvw[1], NGLLX, xigll, hlagx)
            lagrange_poly!(location_result[ipoint].uvw[2], NGLLY, yigll, hlagy)
            lagrange_poly!(location_result[ipoint].uvw[3], NGLLZ, zigll, hlagz)

            for igllz = 1:NGLLZ
                for iglly = 1:NGLLY
                    for igllx = 1:NGLLX
                        location_result[ipoint].lagrange[igllx,iglly,igllz] = hlagx[igllx] * hlagy[iglly] * hlagz[igllz]
                    end
                end
            end
        end
        # to debug
        # if (ipoint == 492835) && (location_result[ipoint].stat != -1)
        #     @info "@@@@"
        #     @info iproc_old
        #     @info xyz1
        #     @info idxs
        #     @info hlagx
        #     @info hlagy
        #     @info hlagz
        #     @info location_result[ipoint].uvw
        #     @info "@@@@"
        # end
    end
    return location_result
end

function prepare_kdtree(mesh_data::sem_mesh_data)
    # * GLL colocation points and lagrange interpolation weights
    xigll = zeros(Float64, NGLLX)
    wxgll = zeros(Float64, NGLLX)
    yigll = zeros(Float64, NGLLY)
    wygll = zeros(Float64, NGLLY)
    zigll = zeros(Float64, NGLLZ)
    wzgll = zeros(Float64, NGLLZ)
    hlagx = zeros(Float64, NGLLX)
    hlagy = zeros(Float64, NGLLY)
    hlagz = zeros(Float64, NGLLZ)

    # * mesh dimension variable
    nspec = mesh_data.nspec

    # * coordinates of GLL points
    zwgljd!(xigll, wxgll, NGLLX, GAUSSALPHA, GAUSSBETA)
    zwgljd!(yigll, wygll, NGLLY, GAUSSALPHA, GAUSSBETA)
    zwgljd!(zigll, wzgll, NGLLZ, GAUSSALPHA, GAUSSBETA)

    # * anchor points of mesh elements

    # get index of anchor points in the GLL element
    iax = zeros(Int32, NGNOD)
    iay = zeros(Int32, NGNOD)
    iaz = zeros(Int32, NGNOD)
    anchor_point_index!(iax, iay, iaz)

    # get anchor and center points of each GLL element
    xyz_elem = zeros(Float64, 3, nspec)
    xyz_anchor = zeros(Float64, 3, NGNOD, nspec)

    for ispec = 1:nspec
        for ia = 1:NGNOD
            iglob = mesh_data.ibool[iax[ia], iay[ia], iaz[ia], ispec]
            xyz_anchor[:,ia,ispec] = mesh_data.xyz_glob[:,iglob]
        end
        # the last anchor point is the element center
        xyz_elem[:,ispec] = xyz_anchor[:,NGNOD,ispec]
    end
    # * build kdtree
    kdtree = KDTree(xyz_elem)

    return kdtree, xyz_elem, xyz_anchor, xigll, yigll, zigll, hlagx, hlagy, hlagz
end

function sem_mesh_locate_kdtree_series!(ipoint::Int64, nnearest::Int64, max_search_dist::Float64, max_misloc::Float64, nmodel::Int64, kdtree::KDTree, mesh_data::sem_mesh_data, xyz::Array{Float64,2}, idoubling::Vector{Int64}, xyz_elem::Array{Float64,2}, xyz_anchor::Array{Float64,3}, xigll::Array{Float64,1}, yigll::Array{Float64,1}, zigll::Array{Float64,1}, model_gll_old::Array{Float64,5}, hlagx::Array{Float64,1}, hlagy::Array{Float64,1}, hlagz::Array{Float64,1}, stat_final::Array{Int64,1}, misloc_final::Array{Float64,1}, model_interp::Array{Float64,2})
    location_result = sem_mesh_location()
    # init
    location_result.stat = -1
    location_result.eid = -1
    location_result.misloc = HUGEVAL
    location_result.lagrange = zeros(Float64, NGLLX, NGLLY, NGLLZ)
    location_result.uvw = zeros(Float64, 3)

    # calculate location_result
    # * get the n nearest elements in the mesh
    xyz1 = xyz[:,ipoint]
    idxs, _ = knn(kdtree, xyz1, nnearest)

    # * test each neighbour elements to see if target point is located inside
    for inn = 1:nnearest
        ispec = idxs[inn]
        # skip the element a certain distance away
        dist = sqrt(sum((xyz_elem[:,ispec] .- xyz1).^2))
        if dist > max_search_dist
            continue
        end
        # only use element with the same layer ID (i.e. idoubling)
        if (idoubling[ipoint] != IFLAG_DUMMY && mesh_data.idoubling[ispec] != idoubling[ipoint])
            continue
        end

        # locate point to this element
        uvw1 = similar(xyz1)
        misloc1 = HUGEVAL
        flag_inside = false
        misloc1, flag_inside = xyz2cube_bounded!(xyz_anchor[:,:,ispec], xyz1, uvw1, misloc1, flag_inside)

        if flag_inside == true
            location_result.stat = 1
            location_result.eid = ispec
            location_result.misloc = misloc1
            location_result.uvw = uvw1
            break
        else
            if (misloc1 < max_misloc) && (misloc1 < location_result.misloc)
                location_result.stat = 0
                location_result.eid = ispec
                location_result.misloc = misloc1
                location_result.uvw = uvw1
            end
        end
    end

    # * set interpolation weights on GLL points if located
    if location_result.stat != -1
        lagrange_poly!(location_result.uvw[1], NGLLX, xigll, hlagx)
        lagrange_poly!(location_result.uvw[2], NGLLY, yigll, hlagy)
        lagrange_poly!(location_result.uvw[3], NGLLZ, zigll, hlagz)
        for igllz = 1:NGLLZ
            for iglly = 1:NGLLY
                for igllx = 1:NGLLX
                    location_result.lagrange[igllx,iglly,igllz] = hlagx[igllx] * hlagy[iglly] * hlagz[igllz]
                end
            end
        end
    end

    # * perform the model interpolation
    if (stat_final[ipoint] == 1 && location_result.stat == 1)
        return
    end
    # for point located inside one element in the first time or closer to one element than located before
    if location_result.stat == 1 || (location_result.stat == 0 && location_result.misloc < misloc_final[ipoint])
        for imodel = 1:nmodel
            model_interp[imodel,ipoint] = sum(location_result.lagrange .* model_gll_old[imodel,:,:,:,location_result.eid])
        end
        stat_final[ipoint] = location_result.stat
        misloc_final[ipoint] = location_result.misloc
    end

end

max_stars_repo_name: ziyixi/seisflow

max_issues_repo_name: ziyixi/seisflow

max_forks_repo_name: ziyixi/seisflow


text:
typedef struct WrongSpelled *Dab;

struct SpelledWrong {
  int x;
}

int main() {
  return 0;
}

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
type Exp3 <: BanditAlgorithm
  gamma::Float64
  weights::Vector{Float64}
end

function Exp3(n_arms::Int64, gamma::Float64)
  Exp3(gamma, ones(n_arms))
end

function initialize(algo::Exp3, n_arms::Int64)
  algo.weights = ones(n_arms)
end

function initialize(algo::Exp3, arms::Vector)
  initialize(algo, length(arms))
end

function select_arm(algo::Exp3)
  K = length(algo.weights)
  probs = (1 - algo.gamma) * (algo.weights ./ sum(algo.weights)) + (algo.gamma / K)
  return rand(Categorical(probs))
end

function update(algo::Exp3, chosen_arm::Int64, reward::Real)
  K = length(algo.weights)
  probs = (1 - algo.gamma) * (algo.weights ./ sum(algo.weights)) + (algo.gamma / K)

  x = zeros(K)
  x[chosen_arm] = reward / probs[chosen_arm]

  algo.weights = algo.weights .* exp((algo.gamma / K) .* x)
end

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook