text:
using Random, Plots, StatsPlots, Distributions, DataFrames, Optim
# Bayesian Inference for Bernoulli parameter (inference success probability)
# Define HPD (highest posterior density) interval
# HPD are difficult to calculate and not be used nowadays.
# Maybe this? https://github.com/tpapp/HighestDensityRegions.jl
"""
HPD interval for beta dist
Optim.jl
http://julianlsolvers.github.io/Optim.jl/v0.9.3/user/tipsandtricks/#dealing-with-constant-parameters

    Inputs
    ------
    ci0:  initial value for calculating HPD interval numerically
    α:    parameter 1 of Beta dist
    β:    parameter 2 of Beta dist
    prob: probabilty of HPD interval

    Outputs
    -------
    HPD interval

"""
function beta_hpdi(ci0, α, β, prob)
    """
    Continuous equation for HPD interval

        Inputs
        ------
        v: HPD interval (2-element Vector)
        a: parameter 1 of Beta dist
        b: parameter 2 of Beta dist
        p: probabilty of HPD interval

        Outputs
        -------
        HPD interval condition equations (2-element Vector)

    """
    function hpdi_conditions(v, a, b, p)
        # HPD interval probability = p
        eq1 = cdf(Beta(a, b), v[2]) - cdf(Beta(a, b), v[1]) - p
        # p(lower|D) = p(upper|D)
        eq2 = pdf(Beta(a, b), v[2]) - pdf(Beta(a, b), v[1])
        return [eq1, eq2]
    end
    return nothing
end

function hpdi_conditions(v, a, b, p)
    # HPD interval probability = p
    eq1 = cdf(Beta(a, b), v[2]) - cdf(Beta(a, b), v[1]) - p
    # p(lower|D) = p(upper|D)
    eq2 = pdf(Beta(a, b), v[2]) - pdf(Beta(a, b), v[1])
    return [eq1, eq2]
end

function hpdi_conditions(v)
    # HPD interval probability = p
    eq1 = cdf(Beta(14, 38), v[2]) - cdf(Beta(14, 38), v[1]) - 0.95
    # p(lower|D) = p(upper|D)
    eq2 = pdf(Beta(14, 38), v[2]) - pdf(Beta(14, 38), v[1])
    return [eq1, eq2]
end

# Optimization
hpdi_conditions([0.17, 0.41])
Optim.optimize(hpdi_conditions, [0.17, 0.41])
hpdi_conditions([0, 0])
Optim.optimize(hpdi_conditions, [0, 0])
NLsolve.nlsolve(hpdi_conditions)

# Define calculation of posterior for Bernoulli parameter
"""
Posterior Distribution (Beta dist) for Bernoulli parameter π

    Inputs
    ------
    data: data (0 or 1)
    a0: parameter 1 of prior (Beta dist)
    b0: parameter 2 of prior (Beta dist)
    prob: credible interval probability (0 < prob < 1)

    Outputs
    -------
    results: posterior summary stats table (DataFrames)
    a: parameter 1 of posterior (Beta dist)
    b: parameter 2 of posterior (Beta dist)

"""
function bernoulli_stats(data, a0, b0, prob)
    # posterior
    n = size(data)[1]
    sum_data = sum(data)
    a = sum_data + a0
    b = n - sum_data + b0
    d_posterior = Beta(a, b)
    # mean
    mean_π = mean(d_posterior)
    # median
    median_π = median(d_posterior)
    # mode
    mode_π = mode(d_posterior)
    # Analytic Mode for Beta Dist
    # (a - 1.0) / (a + b - 2.0)
    # std
    sd_π = std(d_posterior)
    # credible interval
    # ci upper boud
    ci_π_upper = quantile(d_posterior, prob - (1-prob)/2)
    # ci lower boud
    ci_π_lower = quantile(d_posterior, (1-prob)/2)
    ci_π = [ci_π_lower, ci_π_upper]

    posterior_stats = reshape([mean_π, median_π, mode_π, sd_π, ci_π], 1, 5)
    stats_string = ["Mean", "Median", "Mode", "SD", "CI"]
    results = DataFrame(posterior_stats, stats_string)
    return results, a, b
end

# Generate moc data from Bernoulli
p = 0.25
n = 50
Random.seed!(99)
data = rand(Bernoulli(p), n)

# Set prior
a0 = 1.0
b0 = 1.0

# Calculate posterior statistics
prob = 0.95
results, a, b = bernoulli_stats(data, a0, b0, prob)
print(results)

# Visualize Posterior
plot(Beta(a0, b0), xlims=(0, 1), ylims=(0, 7), xlabel="Success Probability q", ylabel="Probability Density", label="prior: Beta(1, 1)", linestyle=:dash)
plot!(Beta(a, b), label="posterior: Beta($a, $b)", linestyle=:solid)

# Plot posterior by different prior
Random.seed!(99)
data = rand(Bernoulli(p), 250) # more data
value_size = [10, 50, 250]
value_a0 = [1.0, 6.0]
value_b0 = [1.0, 4.0]
styles = [:dot, :dashdot, :dash, :solid]

l = @layout [grid(1, 2)]
plots = []
for index = 1:2
    style_index = 1
    a0_i = value_a0[index]
    b0_i = value_b0[index]
    pl = plot(Beta(a0_i, b0_i), linestyle=styles[style_index], label="prior: Beta($a0_i, $b0_i)", color="cornflower blue",
             xlims=(0, 1), ylims=(0, 15.5), xlabel="Success Probability q", ylabel="Probability Density")
    for n_j in value_size
        style_index += 1
        sum_data = sum(data[1:n_j])
        a_j = sum_data + a0_i
        b_j = n_j - sum_data + b0_i
        plot!(Beta(a_j, b_j), linestyle=styles[style_index], label="posterior: Beta($a_j, $b_j)", color="cornflower blue")
    end
    push!(plots, pl)
end
plot(plots..., layout=l)

# Animation
# https://stackoverflow.com/questions/55794068/animating-subplots-using-plots-jl-efficiently
Random.seed!(99)
data = rand(Bernoulli(p), 250) # more data
value_size = Vector(1:1:100)
value_a0 = [1.0, 6.0]
value_b0 = [1.0, 4.0]

gr(fmt = :png)
anim = @animate for n_j in value_size
    sum_data = sum(data[1:n_j])
    a_j_1 = sum_data + value_a0[1]
    b_j_1 = n_j - sum_data + value_b0[1]
    # all_p[1][1][:z] = Beta(a_j, b_j)
    p1 = plot(Beta(a_j_1, b_j_1), color="cornflower blue", label="posterior: Beta($a_j_1, $b_j_1)", xlims=(0, 1), ylims=(0, 15.5), xlabel="Success Probability q", ylabel="Probability Density", title="n=$n_j")
    a0_1 = value_a0[1]
    b0_1 = value_b0[1]
    plot!(Beta(a0_1, b0_1), linestyle=:dot, label="prior: Beta($a0_1, $b0_1)")

    a_j_2 = sum_data + value_a0[2]
    b_j_2 = n_j - sum_data + value_b0[2]
    # all_p[2][1][:z] = Beta(a_j, b_j)
    p2 = plot(Beta(a_j_2, b_j_2), color="cornflower blue", label="posterior: Beta($a_j_2, $b_j_2)", xlims=(0, 1), ylims=(0, 15.5), xlabel="Success Probability q", ylabel="Probability Density", title="n=$n_j")
    a0_2 = value_a0[2]
    b0_2 = value_b0[2]
    plot!(Beta(a0_2, b0_2), linestyle=:dot, label="prior: Beta($a0_2, $b0_2)")

    layout = @layout [grid(1, 2)]
    plot(p1, p2; layout)
end
gif(anim, "posterior.gif", fps=10)


text:
using Test
using CompScienceMeshes

for T in [Float32, Float64]
    mesh = meshrectangle(T(1.0), T(1.0), T(1.0))
    faces = skeleton(mesh, 2)
    local verts = vertices(mesh, faces.faces[1])
    p = simplex(verts)

    @test p.vertices == [
        point(T, 0.0, 0.0, 0.0),
        point(T, 0.0, 1.0, 0.0),
        point(T, 1.0, 0.0, 0.0)]
    @test p.tangents == [
        point(T, -1.0, 0.0, 0.0),
        point(T, -1.0, 1.0, 0.0)]
    @test p.normals == [
        point(T, 0.0, 0.0, -1.0)]
    @test p.volume == T(0.5)
end

ext: jl

max_stars_repo_name: krcools/CompScienceMeshes

max_issues_repo_name: krcools/CompScienceMeshes

max_forks_repo_name: krcools/CompScienceMeshes


text:
using StaticOptim, Optim
using Test

# All of these tests are from OptimTestProblems.jl

using StaticArrays
sx = @SVector ones(2)
sx = 3.2 * sx
rosenbrock(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
res = soptimize(rosenbrock, sx)
@test res.g_converged == true
@test rosenbrock(res.minimizer) == res.minimum
@test show(res) === nothing
res = soptimize(rosenbrock, sx)
@test res.g_converged == true
res = soptimize(rosenbrock, sx, updating = true)
@test res.g_converged == true
res = soptimize(rosenbrock, sx, bto = StaticOptim.Order3())
@test res.g_converged == true
@test rosenbrock(res.minimizer) == res.minimum

res = soptimize(rosenbrock, sx/2)
@test res.g_converged == true
@test rosenbrock(res.minimizer) == res.minimum

# Constrained
@testset "Contrained" begin
function U(h, w)
    h1, h2 = h[1], h[2]
    h1 >= 1 && return -Inf*one(h1)
    h2 >= 1 && return -Inf*one(h1)
    c = 2w*h1 + 1.5w*h2 + 0.1
    c <= 0 && return -Inf*one(h1)
    log(c) + log(1 - h1) + log(1 - h2)
end

l = SVector{2}(0., 0.4)
sx = SVector{2}(0.7, 0.7)
res = constrained_soptimize(x -> -U(x, 1), sx, lower = l)
@test res.g_converged == true
lbs = [SVector{2}(x, y) for x in 0:0.1:0.4, y in 0:0.1:0.4]
ubs = [SVector{2}(x, y) for x in 0.6:.1:0.9, y in 0.6:0.1:0.9]
ws = 0.7:0.1:1.01
sx = SVector{2}(0.5, 0.5)
constrained_soptimize(x -> -U(x, 1.), sx, lower = SVector{2}(0., 0.1))
sols = [(w = w, l = l, res = constrained_soptimize(x -> -U(x, w), sx, lower = l)) for l in lbs, u in ubs, w in ws]
@test all(sol.res.g_converged == true for sol in sols)

## Constrained, 4 people

function U(h, w)
    h1, h2, h3, h4 = h[1], h[2], h[3], h[4]
    h1 >= 1 && return -Inf*one(h1)
    h2 >= 1 && return -Inf*one(h1)
    h3 >= 1 && return -Inf*one(h1)
    h4 >= 1 && return -Inf*one(h1)
    c = 2w*h1 + 1.5w*h2 + 1.2w*h3 + 0.5w*h4 + 0.1
    c <= 0 && return -Inf*one(h1)
    log(c) + log(1 - h1) + log(1 - h2) + log(1 - h3) + log(1 - h4)
end

l = SVector{4}(0., 0., 0., 0.)
sx = SVector{4}(0.7, 0.7, 0.7, 0.7)
res = constrained_soptimize(x -> -U(x, 1), sx, lower = l)
@test res.g_converged == true
lbs = [SVector{4}(x, y, z, t) for x in 0:0.1:0.4, y in 0:0.1:0.4, z in 0:0.1:0.4, t in 0:0.1:0.4]
ubs = [SVector{4}(x, y, z, t) for x in 0.6:.1:0.9, y in 0.6:0.1:0.9, z in 0.6:0.1:0.9, t in 0.6:0.1:0.9]
ws = 0.5:0.5:2.
sx = SVector{4}(0.5, 0.5, 0.5, 0.5)
constrained_soptimize(x -> -U(x, 1), sx, lower = SVector{4}(0.3, 0., 0.2, 0.2))
#sols = [(w = w, l = l, res = constrained_soptimize(x -> -U(x, w), sx, lower = l, upper = u)) for l in lbs, u in ubs, w in ws];
@test all(constrained_soptimize(x -> -U(x, w), sx, lower = l, upper = u).g_converged for l in lbs, u in ubs, w in ws)
#sols = [(w = 1, l = l, res = constrained_soptimize(x -> -U(x, 1), sx, lower = l)) for l in lbs]
end

function fletcher_powell(x::AbstractVector)
    function theta(x::AbstractVector)
        if x[1] > 0
            return atan(x[2] / x[1]) / (2.0 * pi)
        else
            return (pi + atan(x[2] / x[1])) / (2.0 * pi)
        end
    end

    return 100.0 * ((x[3] - 10.0 * theta(x))^2 +
        (sqrt(x[1]^2 + x[2]^2) - 1.0)^2) + x[3]^2
end

sx = @SVector ones(3)
sx = 3.2 * sx
res = soptimize(fletcher_powell, sx)
@test res.g_converged == true
@test fletcher_powell(res.minimizer) == res.minimum
res = soptimize(fletcher_powell, sx/2)
@test res.g_converged == true
@test fletcher_powell(res.minimizer) == res.minimum

function himmelblau(x::AbstractVector)
    return (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
end
sx = @SVector ones(2)
sx = 3.2 * sx
res = soptimize(himmelblau, sx)
@test res.g_converged == true
@test himmelblau(res.minimizer) == res.minimum
res = soptimize(himmelblau, sx/2)
@test res.g_converged == true
@test himmelblau(res.minimizer) == res.minimum

function powell(x::AbstractVector)
    return (x[1] + 10.0 * x[2])^2 + 5.0 * (x[3] - x[4])^2 +
        (x[2] - 2.0 * x[3])^4 + 10.0 * (x[1] - x[4])^4
end
sx = @SVector ones(4)
sx = 3.2 * sx
res = soptimize(powell, sx)
@test res.g_converged == true
@test powell(res.minimizer) == res.minimum
res = soptimize(powell, sx/2)
@test res.g_converged == true
@test powell(res.minimizer) == res.minimum

### Univariate tests

f(x) = x^2 + 2*x
res = soptimize(f, 1.)
@test res.g_converged == true
@test f(res.minimizer) == res.minimum

function U(x)
    x = max(0., x)
    log(x)
end
f(h) = -(U(0.2 + h) + U(1 - h))
res = soptimize(f, 0.9)
@test res.g_converged == true
@test f(res.minimizer) == res.minimum

### Root finding
up(c) = c <= 0 ? Inf*c : 1/c
f(a) = up(2 - a) - .96up(2 + a)
out = sroot(f, 0.5)
@test out.fx < 1e-8
@test f(out.x) == out.fx

out = sroot(f, (-0.5, 0.5))
@test  out.isroot == true
@test f(out.x) == out.fx

const w = 3.
const beta = 0.96
const R = 1.01
const alpha_h = 1.5
uc(c) = 1/c
uh(h) = alpha_h/(1 - h)
function eulerfun(x)
    a, h1, h2 = x[1], x[2], x[3]
    c1 = w*h1 - a
    c2 = w*h2 + R*a
    (h1 >= 1 || h2 >= 1) && return Inf*x
    c1 <= 0 && return Inf*x
    out1 = uc(c1) - beta*R*uc(c2)
    out2 = w*uc(c1) - uh(h1)
    out3 = w*uc(c2) - uh(h2)
    SVector(out1, out2, out3)
end
x = SVector(0., 0.5, 0.5)
res = sroot(eulerfun, x)
@test res.g_converged == true
# Test regular arrays (Nonlinear least squares)
using Random
Random.seed!(1234)
realparam = rand(50)
const data = rand(5000, 50)
data[:, 1] = ones(5000)
const y = data*realparam .+ data*exp.(realparam) .+ randn(5000)
paramg = rand(50)
function obj(param)
yhat = data * param .+ data*exp.(param)
sum((yy - yh)^2 for (yy, yh) in zip(y, yhat))
end
sres = soptimize(obj, paramg);
@test res.g_converged == true
res = optimize(obj, paramg, method = Newton(), autodiff = :forward)
@test all(x -> abs(x) < 1e-6, res.minimizer - sres.minimizer)
l = rand(50)
u = rand(50)
u = [uu < ll ? Inf : uu for (uu, ll) in zip(u, l)]
m = (u + l)/2
cparamg = [isfinite(mm) ? mm : 2ll for (mm, ll) in zip(m, l)]
df = TwiceDifferentiable(obj, cparamg, autodiff = :forward)
dfc = TwiceDifferentiableConstraints(l, u)
resoc = optimize(df, dfc, cparamg, IPNewton())
sresoc = constrained_soptimize(obj, cparamg, lower = l, upper = u)

@test all(x -> abs(x) < 1e-6, resoc.minimizer - sresoc.minimizer)


text:
using Test
using MIPVerify
using MIPVerify: UnrestrictedPerturbationFamily, LInfNormBoundedPerturbationFamily
@isdefined(TestHelpers) || include("../../../TestHelpers.jl")

@timed_testset "mfc+mfc+softmax.jl" begin
    input = gen_array((1, 4, 8, 1), 0, 1)

    l1_height = 16
    l1_width = length(input)
    l1_kernel = gen_array((l1_width, l1_height), -1, 1)
    l1_bias = gen_array((l1_height,), -1, 1)

    m1 = [0, 0, 0, 0, -1, -1, -1, 1, 1, -1, 1, 0, 1, 0, -1, 1]

    l2_height = 8
    l2_width = l1_height
    l2_kernel = gen_array((l2_width, l2_height), -1, 1)
    l2_bias = gen_array((l2_height,), -1, 1)

    m2 = [0, -1, 1, 1, -1, -1, 0, 1]

    l3_height = 4
    l3_width = l2_height
    l3_kernel = gen_array((l3_width, l3_height), -1, 1)
    l3_bias = gen_array((l3_height,), -1, 1)

    nn = Sequential(
        [
            Flatten(4),
            Linear(l1_kernel, l1_bias),
            MaskedReLU(m1),
            Linear(l2_kernel, l2_bias),
            MaskedReLU(m2),
            Linear(l3_kernel, l3_bias),
        ],
        "tests.integration.generated_weights.mfc+mfc+softmax"
    )

    pp_unrestricted = UnrestrictedPerturbationFamily()

    @testset "Basic integration test for MaskedReLU layer." begin
        test_cases = [
            ((1, pp_unrestricted, Inf, 0), 0.08112308),
            ((2, pp_unrestricted, Inf, 0), 0.3622567),
            ((3, pp_unrestricted, Inf, 0), 0),
            ((1, LInfNormBoundedPerturbationFamily(0.08), Inf, 0), NaN),
            ((1, LInfNormBoundedPerturbationFamily(0.085), Inf, 0) => 0.08112308),
        ]

        TestHelpers.batch_test_adversarial_example(nn, input, test_cases)
    end
end


text:
struct ComplexExtensionDerivation{T<:FieldElement, P<:PolyElem{T}} <: Derivation
    domain::AbstractAlgebra.ResField{P}
    D::Derivation
    function ComplexExtensionDerivation(domain::AbstractAlgebra.ResField{P}, D::Derivation) where {T<:FieldElement, P<:PolyElem{T}}
        base_ring(base_ring(base_ring(domain)))==D.domain || error("base ring of domain must be domain of D")
        m = modulus(domain)
        degree(m)==2 && isone(coeff(m, 0)) && iszero(coeff(m, 1)) && isone(coeff(m,2)) ||
            error("domain must be residue field modulo X^2+1.")
        new{T,P}(domain, D)
    end
end

function Base.real(f::AbstractAlgebra.ResFieldElem{P}) where {T<:FieldElement, P<:PolyElem{T}}
    #m = modulus(f)
    #degree(m)==2 && isone(coeff(m, 0)) && iszero(coeff(m, 1)) && isone(coeff(m,2)) ||
    #    error("f must be element of residue field modulo X^2+1.")
    coeff(data(f), 0)
end

function Base.imag(f::AbstractAlgebra.ResFieldElem{P}) where {T<:FieldElement, P<:PolyElem{T}}
    #m = modulus(f)
    #degree(m)==2 && isone(coeff(m, 0)) && iszero(coeff(m, 1)) && isone(coeff(m,2)) ||
    #    error("f must be element of residue field modulo X^2+1.")
    coeff(data(f), 1)
end

import Base: (*)

function (*)(c::F, f::K) where
    {F<:AbstractAlgebra.ResFieldElem, T<:AbstractAlgebra.ResFieldElem, P<:PolyElem{T}, K<:FracElem{P}}
    I = get_I(parent(f))
    (real(c) + imag(c)*I)*f
end

function (*)(c::F, f::P) where
    {F<:AbstractAlgebra.ResFieldElem, T<:AbstractAlgebra.ResFieldElem, P<:PolyElem{T}}
    I = get_I(parent(f))
    (real(c) + imag(c)*I)*f
end

Base.iszero(D::ComplexExtensionDerivation) = iszero(D.D)
isbasic(D::ComplexExtensionDerivation) = isbasic(D.D)
MonomialDerivative(D::ComplexExtensionDerivation) = MonomialDerivative(D.D)
BaseDerivation(D::ComplexExtensionDerivation) = D.D

function constant_field(D::ComplexExtensionDerivation)
    C = constant_field(D.D)
    Cz, I = PolynomialRing(C, :I)
    ResidueField(Cz, I^2+1)
end

isconstant(f::AbstractAlgebra.ResFieldElem{P}, D::ComplexExtensionDerivation) where {T<:FieldElement, P<:PolyElem{T}} =
    isconstant(real(f), D.D) && isconstant(imag(f), D.D)

function constantize(f::AbstractAlgebra.ResFieldElem{P}, D::ComplexExtensionDerivation) where {T<:FieldElement, P<:PolyElem{T}}
    u = constantize(real(f), D.D)
    v = constantize(imag(f), D.D)
    C = parent(u)
    Cz, I = PolynomialRing(C, :I)
    CI = ResidueField(Cz, I^2+1)
    CI(u+v*I)
end

isrational(f::AbstractAlgebra.ResFieldElem{P}) where {T<:FieldElement, P<:PolyElem{T}} =
    isrational(real(f)) && iszero(imag(f))

function rationalize(f::AbstractAlgebra.ResFieldElem{P}) where {T<:FieldElement, P<:PolyElem{T}}
    iszero(imag(f)) || error("not rational")
    rationalize(real(f))
end

contains_I(F::AbstractAlgebra.Ring) = false

contains_I(P::PolyRing{T}) where T<:RingElement = contains_I(base_ring(P))

contains_I(F::FracField{T}) where T<:RingElement = contains_I(base_ring(F))

function contains_I(F::AbstractAlgebra.ResField{P}) where {T<:FieldElement, P<:PolyElem{T}}
    m = modulus(F)
    if degree(m)==2 && isone(coeff(m, 0)) && iszero(coeff(m, 1)) && isone(coeff(m,2))
        return true
    else
        error("contains_I not appplicable for residue field with modulus != X^2+1.")
    end
end

function get_I(F::AbstractAlgebra.Ring)
    error("does not contain I")
end

get_I(P::PolyRing{T}) where T<:RingElement = get_I(base_ring(P)) + zero(P)

get_I(F::FracField{T}) where T<:RingElement = get_I(base_ring(F)) + zero(F)

function get_I(F::AbstractAlgebra.ResField{P}) where {T<:FieldElement, P<:PolyElem{T}}
    m = modulus(F)
    if degree(m)==2 && isone(coeff(m, 0)) && iszero(coeff(m, 1)) && isone(coeff(m,2))
        return gen(base_ring(F)) + zero(F)
    else
        error("get_I not appplicable for residue field with modulus != X^2+1.")
    end
end

function conj(f::AbstractAlgebra.ResFieldElem{P}) where {T<:FieldElement, P<:PolyElem{T}}
    I = get_I(parent(f))
    real(f) - imag(f)*I
end

function (D::ComplexExtensionDerivation)(f::AbstractAlgebra.ResFieldElem{P}) where {T<:FieldElement, P<:PolyElem{T}}
    iscompatible(f, D) || error("f not in domain of D")
    #m = modulus(f)
    #degree(m)==2 && isone(coeff(m, 0)) && iszero(coeff(m, 1)) && isone(coeff(m,2)) ||
    #    error("f must be element of residue field modulo X^2+1.")
    u = coeff(data(f), 0)
    v = coeff(data(f), 1)
    kI = parent(f)
    I = gen(base_ring(kI))
    kI(D.D(u)+I*D.D(v))
end

"""
    Complexify(k, D) -> k1, I, D1

Given a field `k` not containing `√-1` and a derivation `D` on `k`, return the complexified field
`k1=k(√-1)`, the generator `I≈√-1` of this field and the unique extension derivation `D1` on `k1` of `D` that
satisfies `D1(√-1)=0`.
"""
function Complexify(k::AbstractAlgebra.Field, D::Derivation) # where {T <:FieldElement, F<: AbstractAlgebra.Field{T}}
    !contains_I(k) || error("k already contains I=sqrt(-1)")
    kz, I = PolynomialRing(k, :I)
    kI = ResidueField(kz, I^2+1)
    DI = ComplexExtensionDerivation(kI, D)
    kI, kI(I), DI
end

"""
    switch_t_i(K, D) -> (K1, t, I, D1)

Switch generators in `k(t)(√-1)``

Given a complexified field of the form `K=k(t)(√-1)` and a derivation `D` on `K`,
return the field `K1=k(√-1)(t)`, the generators `t` and `I≈√-1` for `K1`, and
the derivation `D1` on `K1` corresponding to `D`
such that the differential fields `(K,D)` and `(K1,D1)` are isomorphic.
"""
function switch_t_i(K::AbstractAlgebra.ResField{P}, D::Derivation) where
    {T<:FieldElement, R<:PolyElem{T}, F<:FracElem{R}, P<:PolyElem{F}}
    domain(D)==K || error("K must be the domain of D")
    k = base_ring(base_ring(base_ring(base_ring(K))))
    D0 = BaseDerivation(BaseDerivation(D))
    kI, I, D0I = Complexify(k, D0)
    v = var(base_ring(base_ring(base_ring(K))))
    kIt, t = PolynomialRing(kI, v)
    K1 = FractionField(kIt)
    if iszero(D)
        D1 = NullDerivation(kIt)
    elseif isbasic(D)
        D1 = BasicDerivation(kIt)
    else
        H = MonomialDerivative(D)(t)
        D1 = ExtensionDerivation(kIt, D0I, H)
    end
    K1, t, I, D1
end

"""
    transform(f, t, I) -> f1

transform elements of `k(t)(√-1)` to elements of `k(√-1)(t)`.

Given an element `f` of a complexified field of the form `K=k(t)(√-1)`
and the generators `t` and `I≈√-1` for the field `K1=k(√-1)(t)`
(as returned by `switch_t_i`), return the corresponding element
`f1` in `K1`
"""
function transform(f::K, t, I) where
    {T<:FieldElement, R<:PolyElem{T}, F<:FracElem{R}, P<:PolyElem{F}, K<:AbstractAlgebra.ResFieldElem{P}}
    a = numerator(real(f))(t)
    b = denominator(real(f))(t)
    c = numerator(imag(f))(t)
    d = denominator(imag(f))(t)
    (a*d+c*b*I)//(b*d)
end

"""
    backtransform(f, t, I) -> f1

transform elements of `k(√-1)(t)` to elements of `k(t)(√-1)`.

Given an element `f` of a field of the form `K=k(√-1)(t)`
and the generators `t` and `I≈√-1` for the isomorphic
field `K1=k(t)(√-1)`, return the corresponding element
`f1` in `K1`
"""
function backtransform(f::K, t, I) where
    {T<:AbstractAlgebra.ResFieldElem, P<:PolyElem{T}, K<:FracElem{P}}
    u = map_coefficients(c->real(c), numerator(f))(t)
    v = map_coefficients(c->imag(c), numerator(f))(t)
    z = map_coefficients(c->real(c), denominator(f))(t)
    w = map_coefficients(c->imag(c), denominator(f))(t)
    d = z^2 + w^2
    (u*z + v*w)//d  + ((v*z - u*w)//d)*I
end

function InFieldDerivative(f::K, D::ComplexExtensionDerivation) where
    {T<:FieldElement, R<:PolyElem{T}, F<:FracElem{R}, P<:PolyElem{F}, K<:AbstractAlgebra.ResFieldElem{P}}
    ktI = parent(f)
    I0 = ktI(gen(base_ring(ktI)))
    t0 = gen(base_ring(base_ring(base_ring(ktI))))
    _, t, I, D =  switch_t_i(ktI, D)
    f = transform(f, t, I)
    u, ρ = InFieldDerivative(f, D)
    backtransform(u, t0, I0), ρ
end

#Note: InFieldLogarithmicDerivative is merely a wrapper for InFieldLogarithmicDerivativeOfRadical

function InFieldLogarithmicDerivativeOfRadical(f::K, D::ComplexExtensionDerivation; expect_one::Bool=false) where
    {T<:FieldElement, R<:PolyElem{T}, F<:FracElem{R}, P<:PolyElem{F}, K<:AbstractAlgebra.ResFieldElem{P}}
    ktI = parent(f)
    I0 = ktI(gen(base_ring(ktI)))
    t0 = gen(base_ring(base_ring(base_ring(ktI))))
    _, t, I, D =  switch_t_i(ktI, D)
    f = transform(f, t, I)
    n, u, ρ = InFieldLogarithmicDerivativeOfRadical(f, D, expect_one=expect_one)
    n, backtransform(u, t0, I0), ρ
end

function RischDE(f::K, g::K, D::ComplexExtensionDerivation) where
    {T<:FieldElement, R<:PolyElem{T}, F<:FracElem{R}, P<:PolyElem{F}, K<:AbstractAlgebra.ResFieldElem{P}}
    ktI = parent(f)
    I0 = ktI(gen(base_ring(ktI)))
    t0 = gen(base_ring(base_ring(base_ring(ktI))))
    _, t, I, D =  switch_t_i(ktI, D)
    f = transform(f, t, I)
    g = transform(g, t, I)
    y, ρ = RischDE(f, g, D)
    backtransform(y, t0, I0), ρ
end

function ParamRischDE(f::K, gs::Vector{K}, D::ComplexExtensionDerivation) where
    {T<:FieldElement, R<:PolyElem{T}, F<:FracElem{R}, P<:PolyElem{F}, K<:AbstractAlgebra.ResFieldElem{P}}
    ktI = parent(f)
    I0 = ktI(gen(base_ring(ktI)))
    t0 = gen(base_ring(base_ring(base_ring(ktI))))
    _, t, I, D =  switch_t_i(ktI, D)
    f = transform(f, t, I)
    gs = [transform(g, t, I) for g in gs]
    hs, A = ParamRischDE(f, gs, D) # could A have complex entries ????
    [backtransform(h, t0, I0) for h in hs], A
end

function LimitedIntegrate(f::K, w::K, D::ComplexExtensionDerivation) where
    {T<:FieldElement, R<:PolyElem{T}, F<:FracElem{R}, P<:PolyElem{F}, K<:AbstractAlgebra.ResFieldElem{P}}
    ktI = parent(f)
    I0 = ktI(gen(base_ring(ktI)))
    t0 = gen(base_ring(base_ring(base_ring(ktI))))
    _, t, I, D =  switch_t_i(ktI, D)
    f = transform(f, t, I)
    w = transform(w, t, I)
    v, cs, ρ = LimitedIntegrate(f, w, D) # could cs have complex elements ???
    backtransform(v, t0, I0), cs, ρ
end

function ParametricLogarithmicDerivative(f::K, w::K, D::ComplexExtensionDerivation) where
    {T<:FieldElement, R<:PolyElem{T}, F<:FracElem{R}, P<:PolyElem{F}, K<:AbstractAlgebra.ResFieldElem{P}}
    ktI = parent(f)
    I0 = ktI(gen(base_ring(ktI)))
    t0 = gen(base_ring(base_ring(base_ring(ktI))))
    _, t, I, D =  switch_t_i(ktI, D)
    f = transform(f, t, I)
    w = transform(w, t, I)
    n, m, v, ρ = ParametricLogarithmicDerivative(f, w, D)
    n, m, backtransform(v, t0, I0), ρ
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct AuthenticationV1Api <: SwaggerApi
    client::Swagger.Client
end

"""

create a TokenReview
Param: body::IoK8sApiAuthenticationV1TokenReview (required)
Param: dryRun::String
Param: fieldManager::String
Param: pretty::String
Return: IoK8sApiAuthenticationV1TokenReview
"""
function createAuthenticationV1TokenReview(_api::AuthenticationV1Api, body; dryRun=nothing, fieldManager=nothing, pretty=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", IoK8sApiAuthenticationV1TokenReview, "/apis/authentication.k8s.io/v1/tokenreviews", ["BearerToken"], body)
    Swagger.set_param(_ctx.query, "dryRun", dryRun)  # type String
    Swagger.set_param(_ctx.query, "fieldManager", fieldManager)  # type String
    Swagger.set_param(_ctx.query, "pretty", pretty)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "application/yaml", "application/vnd.kubernetes.protobuf"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["*/*"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

get available resources
Return: IoK8sApimachineryPkgApisMetaV1APIResourceList
"""
function getAuthenticationV1APIResources(_api::AuthenticationV1Api; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", IoK8sApimachineryPkgApisMetaV1APIResourceList, "/apis/authentication.k8s.io/v1/", ["BearerToken"])
    Swagger.set_header_accept(_ctx, ["application/json", "application/yaml", "application/vnd.kubernetes.protobuf"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "application/yaml", "application/vnd.kubernetes.protobuf"] : [_mediaType])
    Swagger.exec(_ctx)
end

export createAuthenticationV1TokenReview, getAuthenticationV1APIResources


text:
module TestStrayCopses
using Test

for file in sort([file for file in readdir(@__DIR__) if
                                   occursin(r"^test[_0-9]+.*\.jl$", file)])
    m = match(r"test[_0-9]+(.*).jl", file)

    @testset "$(m[1])" begin
        # Here you can optionally exclude some test files
        # VERSION < v"1.1" && file == "test_xxx.jl" && continue

        include(file)
    end
end

end  # module


text:
@testset "heteroscedastic" begin
    N, d = 500, 1
    k = SqExponentialKernel() ∘ ScaleTransform(10.0)
    X, f = generate_f(N, d, k)
    X, g = generate_f(N, d, k; X=X)
    g .-= 3.0
    λ = 2.0
    σ = λ * AGP.logistic.(g)

    y = f + rand.(Normal.(0.0, sqrt.(inv.(σ))))
    floattypes = [Float64]
    tests_likelihood(
        HeteroscedasticLikelihood(λ),
        AGP.GPLikelihoods.HeteroscedasticGaussianLikelihood,
        Dict(
            "VGP" => Dict("AVI" => true, "QVI" => false, "MCVI" => false),
            "SVGP" => Dict("AVI" => true, "QVI" => false, "MCVI" => false),
            "OSVGP" => Dict("AVI" => true, "QVI" => false, "MCVI" => false),
            "MCGP" => Dict("Gibbs" => true, "HMC" => false),
        ),
        floattypes,
        "Regression",
        2,
        X,
        f,
        y,
        k,
    )
end


text:
export rl_unplugged_dm_dataset

using TFRecord

function make_batch_array(type::Type, feature_dims::Int, size::Tuple, batch_size::Int)
    Array{type, feature_dims+1}(undef, size..., batch_size)
end

function dm_buffer_dict(feature_size::Dict{String, Tuple}, batch_size::Int)
    obs_buffer = Dict{Symbol, AbstractArray}()

    buffer_dict = Dict{Symbol, Any}()

    for feature in keys(feature_size)
        feature_dims = length(feature_size[feature])
        if split(feature, "/")[1] == "observation"
            ob_key = Symbol(chop(feature, head=length("observation")+1, tail=0))
            if split(feature, "/")[end] == "egocentric_camera"
                obs_buffer[ob_key] = make_batch_array(UInt8, feature_dims, feature_size[feature], batch_size)
            else
                obs_buffer[ob_key] = make_batch_array(Float32, feature_dims, feature_size[feature], batch_size)
            end
        elseif feature == "action"
            buffer_dict[:action] = make_batch_array(Float32, feature_dims, feature_size[feature], batch_size)
            buffer_dict[:next_action] = make_batch_array(Float32, feature_dims, feature_size[feature], batch_size)
        elseif feature == "step_type"
            buffer_dict[:terminal] = make_batch_array(Bool, feature_dims, feature_size[feature], batch_size)
        else
            ob_key = Symbol(feature)
            buffer_dict[ob_key] = make_batch_array(Float32, feature_dims, feature_size[feature], batch_size)
        end
    end

    buffer_dict[:state] = deepcopy(NamedTuple(obs_buffer))
    buffer_dict[:next_state] = deepcopy(NamedTuple(obs_buffer))

    buffer_dict
end

function size_dict(game::String, type::String)
    if type == "dm_locomotion_humanoid"
        size_dict = DM_LOCOMOTION_HUMANOID_SIZE
    elseif type == "dm_locomotion_rodent"
        size_dict = DM_LOCOMOTION_RODENT_SIZE
    elseif type == "dm_control_suite"
        size_dict = DM_CONTROL_SUITE_SIZE[game]
    else
        error("given game type does not exist")
    end
end

function batch_named_tuple!(dest::NamedTuple, src::NamedTuple, i::Int)
    for fn in fieldnames(typeof(dest))
        xs = getfield(dest, fn)
        x = getfield(src, fn)
        if typeof(xs) <: NamedTuple
            batch_named_tuple!(xs, x, i)
        else
            selectdim(xs, ndims(xs), i) .= x
        end
    end
end

function make_transition(example::TFRecord.Example, feature_size::Dict{String, Tuple})
    f = example.features.feature

    observation_dict = Dict{Symbol, AbstractArray}()
    next_observation_dict = Dict{Symbol, AbstractArray}()
    transition_dict = Dict{Symbol, Any}()

    for feature in keys(feature_size)
        if split(feature, "/")[1] == "observation"
            ob_key = Symbol(chop(feature, head = length("observation")+1, tail=0))
            if split(feature, "/")[end] == "egocentric_camera"
                cam_feature_size = feature_size[feature]
                ob_size = prod(cam_feature_size)
                observation_dict[ob_key] = reshape(f[feature].bytes_list.value[1][1:ob_size], cam_feature_size...)
                next_observation_dict[ob_key] = reshape(f[feature].bytes_list.value[1][ob_size+1:end], cam_feature_size...)
            else
                if feature_size[feature] == ()
                    observation_dict[ob_key] = f[feature].float_list.value
                else
                    ob_size = feature_size[feature][1]
                    observation_dict[ob_key] = f[feature].float_list.value[1:ob_size]
                    next_observation_dict[ob_key] = f[feature].float_list.value[ob_size+1:end]
                end
            end
        elseif feature == "action"
            ob_size = feature_size[feature][1]
            action = f[feature].float_list.value
            transition_dict[:action] = action[1:ob_size]
            transition_dict[:next_action] = action[ob_size+1:end]
        elseif feature == "step_type"
            transition_dict[:terminal] = f[feature].float_list.value[1] == 2
        else
            ob_key = Symbol(feature)
            transition_dict[ob_key] = f[feature].float_list.value[1]
        end
    end
    state_nt = (state = NamedTuple(observation_dict),)
    next_state_nt = (next_state = NamedTuple(next_observation_dict),)
    transition = NamedTuple(transition_dict)

    merge(transition, state_nt, next_state_nt)
end
"""
    rl_unplugged_dm_dataset(game, shards; <keyword arguments>)

Returns a `RingBuffer`(@ref) of `NamedTuple` containing SARTS batches which supports
multi threaded loading. Also contains additional data. The data enclosed within `:state` and
`next_state` is a NamedTuple consisting of all observations that are provided.
Check out keys in `DM_LOCOMOTION_HUMANOID`, `DM_LOCOMOTION_RODENT`, `DM_CONTROL_SUITE_SIZE` for supported
datasets. Also check out `dm_params()` for more info on arguments.

# Arguments

- `game::String`: name of the dataset.
- `shards::Vector{Int}`: the shards that are to be loaded.
- `type::String`: type of the dm_env. can be `dm_control_suite`, `dm_locomotion_humanoid`, `dm_locomotion_rodent`.
- `is_shuffle::Bool`
- `stochasticity::Float32`: represents the stochasticity of the dataset. can be
in the range: `0.0:0.1:0.5`.
- `shuffle_buffer_size::Int=10_000`: size of the shuffle_buffer used in loading AtariRLTransitions.
- `tf_reader_bufsize::Int=10_000`: the size of the buffer `bufsize` that is used internally
in `TFRecord.read`.
- `tf_reader_sz::Int=10_000`: the size of the `Channel`, `channel_size` that is returned by
`TFRecord.read`.
- `batch_size::Int=256`: The number of samples within the batches that are returned by the `Channel`.
- `n_preallocations::Int=nthreads()*12`: the size of the buffer in the `Channel` that is returned.

!!! note

    To enable reading records from multiple files concurrently, remember to set the number of
    threads correctly (See [JULIA_NUM_THREADS](https://docs.julialang.org/en/v1/manual/environment-variables/#JULIA_NUM_THREADS)).
"""
function rl_unplugged_dm_dataset(
    game,
    shards;
    type = "dm_control_suite",
    is_shuffle = true,
    shuffle_buffer_size=10_000,
    tf_reader_bufsize=10_000,
    tf_reader_sz=10_000,
    batch_size=256,
    n_preallocations=nthreads()*12
)
    n = nthreads()

    repo = "rl-unplugged-dm"

    folders= [
        @datadep_str "$repo-$game-$shard"
        for shard in shards
    ]

    ch_files = Channel{String}(length(folders)) do ch
        for folder in cycle(folders)
            file = folder * "/$(readdir(folder)[1])"
            put!(ch, file)
        end
    end

    if is_shuffle
        files = buffered_shuffle(ch_files, length(folders))
    else
        files = ch_files
    end

    feature_size = size_dict(game, type)

    ch_src = Channel{NamedTuple}(n * tf_reader_sz) do ch
        for fs in partition(files, n)
            Threads.foreach(
                TFRecord.read(
                    fs;
                    compression=:gzip,
                    bufsize=tf_reader_bufsize,
                    channel_size=tf_reader_sz,
                );
                schedule=Threads.StaticSchedule()
            ) do x
                put!(ch, make_transition(x, feature_size))
            end
        end
    end

    if is_shuffle
        transitions = buffered_shuffle(
        ch_src,
        shuffle_buffer_size
        )
    else
        transitions = ch_src
    end

    taskref = Ref{Task}()

    buffer_dict = dm_buffer_dict(feature_size, batch_size)

    buffer = NamedTuple(buffer_dict)

    res = RingBuffer(buffer;taskref=taskref, sz=n_preallocations) do buff
        Threads.@threads for i in 1:batch_size
            batch_named_tuple!(buff, take!(transitions), i)
        end
    end

    bind(ch_src, taskref[])
    bind(ch_files, taskref[])
    res
end


text:
function __init__()
    @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
        gpu_or_cpu(x::CuArrays.CuArray) = CuArrays.CuArray
        gpu_or_cpu(x::TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}) = CuArrays.CuArray
        gpu_or_cpu(x::Transpose{<:Any,<:CuArrays.CuArray}) = CuArrays.CuArray
        gpu_or_cpu(x::Adjoint{<:Any,<:CuArrays.CuArray}) = CuArrays.CuArray
        gpu_or_cpu(x::Adjoint{<:Any,TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}}) = CuArrays.CuArray
        gpu_or_cpu(x::Transpose{<:Any,TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}}) = CuArrays.CuArray
    end

    @require NLopt="76087f3c-5699-56af-9a33-bf431cd00edd" begin
        function sciml_train(loss, θ, opt::NLopt.Opt, data = DEFAULT_DATA; maxeval=100)
            local x, cur, state
            cur,state = iterate(data)

            function nlopt_grad!(θ,grad)
              _x,lambda = Flux.Zygote.pullback(θ) do θ
                x = loss(θ)
                first(x)
              end

              if length(grad) > 0
                grad .= first(lambda(1))
                @printf("fval:%2.2e  norm:%2.2e\n", _x[1], norm(θ))
              end

              return _x
            end

            NLopt.min_objective!(opt, nlopt_grad!)
            NLopt.maxeval!(opt, maxeval)

            t0= time()
            (minf,minx,ret) = NLopt.optimize(opt, θ)
            _time = time()

            Optim.MultivariateOptimizationResults(opt,
                                                    θ,# initial_x,
                                                    minx, #pick_best_x(f_incr_pick, state),
                                                    minf, # pick_best_f(f_incr_pick, state, d),
                                                    maxeval, #iteration,
                                                    maxeval >= maxeval, #iteration == options.iterations,
                                                    false, # x_converged,
                                                    0.0,#T(options.x_tol),
                                                    0.0,#T(options.x_tol),
                                                    NaN,# x_abschange(state),
                                                    NaN,# x_abschange(state),
                                                    false,# f_converged,
                                                    0.0,#T(options.f_tol),
                                                    0.0,#T(options.f_tol),
                                                    NaN,#f_abschange(d, state),
                                                    NaN,#f_abschange(d, state),
                                                    false,#g_converged,
                                                    0.0,#T(options.g_tol),
                                                    NaN,#g_residual(d),
                                                    false, #f_increased,
                                                    nothing,
                                                    maxeval,
                                                    maxeval,
                                                    0,
                                                    ret,
                                                    NaN,
                                                    _time-t0,)
        end

        function sciml_train(loss, θ, opt::NLopt.Opt, lower_bounds, upper_bounds, data = DEFAULT_DATA; maxeval=100)
          local x, cur, state
          cur,state = iterate(data)

          function nlopt_grad!(θ,grad)
            _x,lambda = Flux.Zygote.pullback(θ) do θ
              x = loss(θ)
              first(x)
            end

            if length(grad) > 0
              grad .= first(lambda(1))
            end

            return _x
          end

          NLopt.min_objective!(opt, nlopt_grad!)
          NLopt.lower_bounds!(opt, lower_bounds)
          NLopt.upper_bounds!(opt, upper_bounds)
          NLopt.maxeval!(opt, maxeval)

          t0 = time()
          NLopt.optimize(opt, θ)
          _time = time()

          Optim.MultivariateOptimizationResults(opt,
                                                    θ,# initial_x,
                                                    minx, #pick_best_x(f_incr_pick, state),
                                                    minf, # pick_best_f(f_incr_pick, state, d),
                                                    maxeval, #iteration,
                                                    maxeval >= maxeval, #iteration == options.iterations,
                                                    false, # x_converged,
                                                    0.0,#T(options.x_tol),
                                                    0.0,#T(options.x_tol),
                                                    NaN,# x_abschange(state),
                                                    NaN,# x_abschange(state),
                                                    false,# f_converged,
                                                    0.0,#T(options.f_tol),
                                                    0.0,#T(options.f_tol),
                                                    NaN,#f_abschange(d, state),
                                                    NaN,#f_abschange(d, state),
                                                    false,#g_converged,
                                                    0.0,#T(options.g_tol),
                                                    NaN,#g_residual(d),
                                                    false, #f_increased,
                                                    nothing,
                                                    maxeval,
                                                    maxeval,
                                                    0,
                                                    ret,
                                                    NaN,
                                                    _time-t0,)
        end
    end
end


text:
function recompile_everything()
    compile_channel = Channel(100)
    max_title_width = 0

    @schedule begin
        # Find all installed packages (except ourselves), and which packages depend on them
        sorted_pkgs = nothing
        cd(Pkg.dir()) do
            global required, available, installed, fixed, dependencies, conflicts, deps
            put!(compile_channel, ("Reading package lists...", 0,))
            required  = Pkg.Reqs.parse("REQUIRE")
            available = Pkg.Read.available()
            installed = Pkg.Read.installed(available)
            fixed = Pkg.Read.fixed(available, installed)

            put!(compile_channel, ("Computing dependencies...", 0))

            dependencies, conflicts = Pkg.Query.dependencies(available, fixed)
            deps = Dict{AbstractString,Vector{AbstractString}}()
            for (name, (version, fixed)) in installed
                if haskey(dependencies, name)
                    deps[name] = collect(keys(dependencies[name][version].requires))
                else
                    deps[name] = []
                end
            end

            function flatten_deps(pkg_deps)
                if isempty(pkg_deps)
                    return pkg_deps
                end
                flat_deps = AbstractString[]
                for dep in pkg_deps
                    push!(flat_deps, dep)
                    for flat_dep in flatten_deps(deps[dep])
                        push!(flat_deps, flat_dep)
                    end
                end
                return unique(flat_deps)
            end

            # Flatten dependency trees
            for (pkg_name, pkg_deps) in deps
                deps[pkg_name] = flatten_deps(pkg_deps)
            end

            # Now sort packages such that dependencies are always in front of their
            # dependents, so that we are always compiling up the tree, never having
            # to compile a dependent package implicitly

            # We get an initial sort by just sorting by the length of dependencies
            sorted_pkgs = sort(collect(keys(deps)), lt=(a,b) -> length(deps[a]) < length(deps[b]))

            # Remove our own package name because we have obviously already been
            # precompiled since we're running!  Also our name is ridiculously long
            sorted_pkgs = filter(x -> x != "MakePkgUpdatePrecompileInTheBackground", sorted_pkgs)

            # Next, swap around every package until it's in behind all of its dependencies
            for idx in 1:length(sorted_pkgs)
                for dep in deps[sorted_pkgs[idx]]
                    # Is this dependency listed _after_ us?  If so, move it closer
                    dep_idx = find(sorted_pkgs .== dep)[1]
                    if dep_idx > idx
                        splice!(sorted_pkgs, dep_idx)
                        insert!(sorted_pkgs, idx, dep)
                    end
                end
            end
        end

        # Calculate maximum title width just so that the printing algorithm can use it
        max_title_width = maximum([strwidth(x) for x in sorted_pkgs])

        # For each package, check to see if it's been cached and if not, cache it
        for idx in 1:length(sorted_pkgs)
            modname = sorted_pkgs[idx]
            put!(compile_channel, ("$(modname).jl...", idx*1.0/length(sorted_pkgs)))
            cache_paths = Base.find_all_in_cache_path(modname)
            if isempty(cache_paths)
                Base.compilecache(modname, true)
            else
                for cachepath in cache_paths
                    try
                        # Essentiallty recreating recompile_stale() from loading.jl
                        # but without all the info()'s as we want to be silent
                        srcpath = Base.find_in_path(modname, nothing)

                        # Can we just auto-delete orphaned cache files here?
                        if srcpath === nothing
                            rm(cachepath)
                            continue
                        end

                        # If it's stale, recompile!
                        if Base.stale_cachefile(srcpath, cachepath)
                            # Simulate compilation.  lol
                            Base.compilecache(modname, true)
                        end
                    end
                end
            end
        end

        put!(compile_channel, ("Done!", 1))
    end

    last_status = "Pkg:"
    report_progress("Pkg") do width
        if isready(compile_channel)
            title, progress = take!(compile_channel)

            if progress == 1
                return nothing
            end

            titlepad = max(max_title_width + 3 - length(title), 0)
            pbar = progress_bar(Float64(progress), width - length(title) - titlepad - 16)
            last_status = @sprintf("Pkg: %s%s [%s] (%.1f%%)", title, " "^titlepad, pbar, 100.0*progress)
        end
        return last_status
    end
    return nothing
end

function supplant_pkgupdate()
    eval(Base.Pkg, :(using MakePkgUpdatePrecompileInTheBackground))
    eval(Base.Pkg, :(update() = (cd(Entry.update,Dir.getmetabranch()); recompile_everything())))
    return nothing
end


text:
using DynamicHMCModels

Random.seed!(1)

cd(@__DIR__)
include("simulateGaussian.jl")

struct GaussianProb{TY <: AbstractVector}
   "Observations."
   y::TY
end

function (problem::GaussianProb)(θ)
   @unpack y = problem   # extract the data
   @unpack mu, sigma = θ
   loglikelihood(Normal(mu, sigma), y) + logpdf(Normal(0,1), mu) +
   logpdf(Truncated(Cauchy(0,5),0,Inf), sigma)
end

# Define problem with data and inits.

data = simulateGaussian(;Nd=100)
p = GaussianProb(data.y);
p((mu = 0.0, sigma = 1.0))

# Write a function to return properly dimensioned transformation.

problem_transformation(p::GaussianProb) =
  as((mu  = as(Real, -25, 25), sigma = asℝ₊), )

# Use Flux for the gradient.
P = TransformedLogDensity(problem_transformation(p), p)

∇P = ADgradient(:ForwardDiff, P)

#import Zygote
#∇P = ADgradient(:Zygote, P)

# Sample from the posterior.
chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 4000);

# Undo the transformation to obtain the posterior from the chain.
posterior = TransformVariables.transform.(Ref(problem_transformation(p)),
get_position.(chain));

chns = nptochain(posterior, NUTS_tuned)

describe(chns)


text:
module MPFR

export
    BigFloat,
    get_bigfloat_precision,
    set_bigfloat_precision,
    with_bigfloat_precision

import
    Base: (*), +, -, /, <, <=, ==, >, >=, ^, besselj, besselj0, besselj1, bessely,
        bessely0, bessely1, ceil, convert, copysign, deg2rad,
        exp, exp2, exponent, factorial, floor, hypot, isinteger, iround,
        isfinite, isinf, isnan, ldexp, log, log2, log10, max, min, mod, modf,
        nextfloat, prevfloat, promote_rule, rad2deg, rem, round, show,
        showcompact, sum, sqrt, string, print, trunc, precision, exp10, expm1,
        gamma, lgamma, digamma, erf, erfc, zeta, log1p, airyai, iceil, ifloor,
        itrunc, eps, signbit, sin, cos, tan, sec, csc, cot, acos, asin, atan,
        cosh, sinh, tanh, sech, csch, coth, acosh, asinh, atanh, atan2,
        serialize, deserialize, inf, nan, cbrt, typemax, typemin,
        realmin, realmax, get_rounding, set_rounding, maxintfloat, widen

import Base.GMP: ClongMax, CulongMax

import Base.Math.lgamma_r

const ROUNDING_MODE = [0]
const DEFAULT_PRECISION = [256]

# Basic type and initialization definitions

type BigFloat <: FloatingPoint
    prec::Clong
    sign::Cint
    exp::Clong
    d::Ptr{Culong}
    function BigFloat()
        N = get_bigfloat_precision()
        z = new(zero(Clong), zero(Cint), zero(Clong), C_NULL)
        ccall((:mpfr_init2,:libmpfr), Void, (Ptr{BigFloat}, Clong), &z, N)
        finalizer(z, Base.GMP._mpfr_clear_func)
        return z
    end
    # Not recommended for general use
    function BigFloat(prec::Clong, sign::Cint, exp::Clong, d::Ptr{Void})
        new(prec, sign, exp, d)
    end
end

widen(::Type{Float64}) = BigFloat
widen(::Type{BigFloat}) = BigFloat

BigFloat(x::BigFloat) = x

for (fJ, fC) in ((:si,:Clong), (:ui,:Culong), (:d,:Float64))
    @eval begin
        function BigFloat(x::($fC))
            z = BigFloat()
            ccall(($(string(:mpfr_set_,fJ)), :libmpfr), Int32, (Ptr{BigFloat}, ($fC), Int32), &z, x, ROUNDING_MODE[end])
            return z
        end
    end
end

function BigFloat(x::BigInt)
    z = BigFloat()
    ccall((:mpfr_set_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigInt}, Int32), &z, &x, ROUNDING_MODE[end])
    return z
end

function BigFloat(x::String, base::Int)
    z = BigFloat()
    err = ccall((:mpfr_set_str, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{Uint8}, Int32, Int32), &z, x, base, ROUNDING_MODE[end])
    if err != 0; error("incorrectly formatted number"); end
    return z
end
BigFloat(x::String) = BigFloat(x, 10)

BigFloat(x::Integer) = BigFloat(BigInt(x))

BigFloat(x::Union(Bool,Int8,Int16,Int32)) = BigFloat(convert(Clong,x))
BigFloat(x::Union(Uint8,Uint16,Uint32)) = BigFloat(convert(Culong,x))

BigFloat(x::Union(Float16,Float32)) = BigFloat(float64(x))
BigFloat(x::Rational) = BigFloat(num(x)) / BigFloat(den(x))

convert(::Type{Rational}, x::BigFloat) = convert(Rational{BigInt}, x)
convert{S}(::Type{BigFloat}, x::Rational{S}) = BigFloat(x) # to resolve ambiguity
convert(::Type{BigFloat}, x::Real) = BigFloat(x)
convert(::Type{FloatingPoint}, x::BigInt) = BigFloat(x)

for to in (Int8, Int16, Int32, Int64)
    @eval begin
        function convert(::Type{$to}, x::BigFloat)
            (isinteger(x) && (typemin($to) <= x <= typemax($to))) || throw(InexactError())
            convert($to, ccall((:mpfr_get_si,:libmpfr),
                               Clong, (Ptr{BigFloat}, Int32), &x, 0))
        end
    end
end

for to in (Uint8, Uint16, Uint32, Uint64)
    @eval begin
        function convert(::Type{$to}, x::BigFloat)
            (isinteger(x) && (typemin($to) <= x <= typemax($to))) || throw(InexactError())
            convert($to, ccall((:mpfr_get_ui,:libmpfr),
                               Culong, (Ptr{BigFloat}, Int32), &x, 0))
        end
    end
end

function Base.BigInt(x::BigFloat)
    !isinteger(x) && throw(InexactError())
    return itrunc(x)
end
convert(::Type{Float64}, x::BigFloat) =
    ccall((:mpfr_get_d,:libmpfr), Float64, (Ptr{BigFloat},Int32), &x, ROUNDING_MODE[end])
convert(::Type{Float32}, x::BigFloat) =
    ccall((:mpfr_get_flt,:libmpfr), Float32, (Ptr{BigFloat},Int32), &x, ROUNDING_MODE[end])

convert(::Type{Integer}, x::BigFloat) = convert(BigInt, x)

promote_rule{T<:Real}(::Type{BigFloat}, ::Type{T}) = BigFloat
promote_rule{T<:FloatingPoint}(::Type{BigInt},::Type{T}) = BigFloat
promote_rule{T<:FloatingPoint}(::Type{BigFloat},::Type{T}) = BigFloat

rationalize(x::BigFloat; tol::Real=eps(x)) = rationalize(BigInt, x, tol=tol)

# serialization

function serialize(s, n::BigFloat)
    Base.serialize_type(s, BigFloat)
    serialize(s, string(n))
end

deserialize(s, ::Type{BigFloat}) = BigFloat(deserialize(s))

# Basic arithmetic without promotion
# Unsigned addition
function +(x::BigFloat, c::Culong)
    z = BigFloat()
    ccall((:mpfr_add_ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, c, ROUNDING_MODE[end])
    return z
end
+(c::Culong, x::BigFloat) = x + c
+(c::CulongMax, x::BigFloat) = x + convert(Culong, c)
+(x::BigFloat, c::CulongMax) = x + convert(Culong, c)

# Signed addition
function +(x::BigFloat, c::Clong)
    z = BigFloat()
    ccall((:mpfr_add_si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, c, ROUNDING_MODE[end])
    return z
end
+(c::Clong, x::BigFloat) = x + c
+(x::BigFloat, c::ClongMax) = x + convert(Clong, c)
+(c::ClongMax, x::BigFloat) = x + convert(Clong, c)

# Float64 addition
function +(x::BigFloat, c::Float64)
    z = BigFloat()
    ccall((:mpfr_add_d, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Float64, Int32), &z, &x, c, ROUNDING_MODE[end])
    return z
end
+(c::Float64, x::BigFloat) = x + c
+(c::Float32, x::BigFloat) = x + convert(Float64, c)
+(x::BigFloat, c::Float32) = x + convert(Float64, c)

# BigInt addition
function +(x::BigFloat, c::BigInt)
    z = BigFloat()
    ccall((:mpfr_add_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigInt}, Int32), &z, &x, &c, ROUNDING_MODE[end])
    return z
end
+(c::BigInt, x::BigFloat) = x + c

# Unsigned subtraction
function -(x::BigFloat, c::Culong)
    z = BigFloat()
    ccall((:mpfr_sub_ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, c, ROUNDING_MODE[end])
    return z
end
function -(c::Culong, x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_ui_sub, :libmpfr), Int32, (Ptr{BigFloat}, Culong, Ptr{BigFloat}, Int32), &z, c, &x, ROUNDING_MODE[end])
    return z
end
-(x::BigFloat, c::CulongMax) = -(x, convert(Culong, c))
-(c::CulongMax, x::BigFloat) = -(convert(Culong, c), x)

# Signed subtraction
function -(x::BigFloat, c::Clong)
    z = BigFloat()
    ccall((:mpfr_sub_si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, c, ROUNDING_MODE[end])
    return z
end
function -(c::Clong, x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_si_sub, :libmpfr), Int32, (Ptr{BigFloat}, Clong, Ptr{BigFloat}, Int32), &z, c, &x, ROUNDING_MODE[end])
    return z
end
-(x::BigFloat, c::ClongMax) = -(x, convert(Clong, c))
-(c::ClongMax, x::BigFloat) = -(convert(Clong, c), x)

# Float64 subtraction
function -(x::BigFloat, c::Float64)
    z = BigFloat()
    ccall((:mpfr_sub_d, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Float64, Int32), &z, &x, c, ROUNDING_MODE[end])
    return z
end
function -(c::Float64, x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_d_sub, :libmpfr), Int32, (Ptr{BigFloat}, Float64, Ptr{BigFloat}, Int32), &z, c, &x, ROUNDING_MODE[end])
    return z
end
-(x::BigFloat, c::Float32) = -(x, convert(Float64, c))
-(c::Float32, x::BigFloat) = -(convert(Float64, c), x)

# BigInt subtraction
function -(x::BigFloat, c::BigInt)
    z = BigFloat()
    ccall((:mpfr_sub_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigInt}, Int32), &z, &x, &c, ROUNDING_MODE[end])
    return z
end
function -(c::BigInt, x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_z_sub, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigInt}, Ptr{BigFloat}, Int32), &z, &c, &x, ROUNDING_MODE[end])
    return z
end

# Unsigned multiplication
function *(x::BigFloat, c::Culong)
    z = BigFloat()
    ccall((:mpfr_mul_ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, c, ROUNDING_MODE[end])
    return z
end
*(c::Culong, x::BigFloat) = x * c
*(c::CulongMax, x::BigFloat) = x * convert(Culong, c)
*(x::BigFloat, c::CulongMax) = x * convert(Culong, c)

# Signed multiplication
function *(x::BigFloat, c::Clong)
    z = BigFloat()
    ccall((:mpfr_mul_si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, c, ROUNDING_MODE[end])
    return z
end
*(c::Clong, x::BigFloat) = x * c
*(x::BigFloat, c::ClongMax) = x * convert(Clong, c)
*(c::ClongMax, x::BigFloat) = x * convert(Clong, c)

# Float64 multiplication
function *(x::BigFloat, c::Float64)
    z = BigFloat()
    ccall((:mpfr_mul_d, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Float64, Int32), &z, &x, c, ROUNDING_MODE[end])
    return z
end
*(c::Float64, x::BigFloat) = x * c
*(c::Float32, x::BigFloat) = x * convert(Float64, c)
*(x::BigFloat, c::Float32) = x * convert(Float64, c)

# BigInt multiplication
function *(x::BigFloat, c::BigInt)
    z = BigFloat()
    ccall((:mpfr_mul_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigInt}, Int32), &z, &x, &c, ROUNDING_MODE[end])
    return z
end
*(c::BigInt, x::BigFloat) = x * c

# Unsigned division
function /(x::BigFloat, c::Culong)
    z = BigFloat()
    ccall((:mpfr_div_ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, c, ROUNDING_MODE[end])
    return z
end
function /(c::Culong, x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_ui_div, :libmpfr), Int32, (Ptr{BigFloat}, Culong, Ptr{BigFloat}, Int32), &z, c, &x, ROUNDING_MODE[end])
    return z
end
/(x::BigFloat, c::CulongMax) = /(x, convert(Culong, c))
/(c::CulongMax, x::BigFloat) = /(convert(Culong, c), x)

# Signed division
function /(x::BigFloat, c::Clong)
    z = BigFloat()
    ccall((:mpfr_div_si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, c, ROUNDING_MODE[end])
    return z
end
function /(c::Clong, x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_si_div, :libmpfr), Int32, (Ptr{BigFloat}, Clong, Ptr{BigFloat}, Int32), &z, c, &x, ROUNDING_MODE[end])
    return z
end
/(x::BigFloat, c::ClongMax) = /(x, convert(Clong, c))
/(c::ClongMax, x::BigFloat) = /(convert(Clong, c), x)

# Float64 division
function /(x::BigFloat, c::Float64)
    z = BigFloat()
    ccall((:mpfr_div_d, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Float64, Int32), &z, &x, c, ROUNDING_MODE[end])
    return z
end
function /(c::Float64, x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_d_div, :libmpfr), Int32, (Ptr{BigFloat}, Float64, Ptr{BigFloat}, Int32), &z, c, &x, ROUNDING_MODE[end])
    return z
end
/(x::BigFloat, c::Float32) = /(x, convert(Float64, c))
/(c::Float32, x::BigFloat) = /(convert(Float64, c), x)

# BigInt division
function /(x::BigFloat, c::BigInt)
    z = BigFloat()
    ccall((:mpfr_div_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigInt}, Int32), &z, &x, &c, ROUNDING_MODE[end])
    return z
end

# Basic operations
for (fJ, fC) in ((:+,:add), (:-,:sub), (:*,:mul), (:/,:div), (:^, :pow))
    @eval begin
        function ($fJ)(x::BigFloat, y::BigFloat)
            z = BigFloat()
            ccall(($(string(:mpfr_,fC)),:libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[end])
            return z
        end
    end
end

# More efficient commutative operations
for (fJ, fC, fI) in ((:+, :add, 0), (:*, :mul, 1))
    @eval begin
        function ($fJ)(a::BigFloat, b::BigFloat, c::BigFloat)
            z = BigFloat()
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &a, &b, ROUNDING_MODE[end])
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &z, &c, ROUNDING_MODE[end])
            return z
        end
        function ($fJ)(a::BigFloat, b::BigFloat, c::BigFloat, d::BigFloat)
            z = BigFloat()
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &a, &b, ROUNDING_MODE[end])
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &z, &c, ROUNDING_MODE[end])
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &z, &d, ROUNDING_MODE[end])
            return z
        end
        function ($fJ)(a::BigFloat, b::BigFloat, c::BigFloat, d::BigFloat, e::BigFloat)
            z = BigFloat()
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &a, &b, ROUNDING_MODE[end])
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &z, &c, ROUNDING_MODE[end])
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &z, &d, ROUNDING_MODE[end])
            ccall(($(string(:mpfr_,fC)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &z, &e, ROUNDING_MODE[end])
            return z
        end
    end
end

function -(x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_neg, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
    return z
end

function sqrt(x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_sqrt, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
    if isnan(z)
        throw(DomainError())
    end
    return z
end

sqrt(x::BigInt) = sqrt(BigFloat(x))

rad2deg(z::BigFloat) = 180/big(pi)*z
deg2rad(z::BigFloat) = big(pi)/180*z

function ^(x::BigFloat, y::CulongMax)
    z = BigFloat()
    ccall((:mpfr_pow_ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, y, ROUNDING_MODE[end])
    return z
end

function ^(x::BigFloat, y::ClongMax)
    z = BigFloat()
    ccall((:mpfr_pow_si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, y, ROUNDING_MODE[end])
    return z
end

function ^(x::BigFloat, y::BigInt)
    z = BigFloat()
    ccall((:mpfr_pow_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigInt}, Int32), &z, &x, &y, ROUNDING_MODE[end])
    return z
end

for f in (:exp, :exp2, :exp10, :expm1, :digamma, :erf, :erfc, :zeta,
          :cosh,:sinh,:tanh,:sech,:csch,:coth, :cbrt)
    @eval function $f(x::BigFloat)
        z = BigFloat()
        ccall(($(string(:mpfr_,f)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
        return z
    end
end

function airyai(x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_ai, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
    return z
end

function ldexp(x::BigFloat, n::Clong)
    z = BigFloat()
    ccall((:mpfr_mul_2si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, n, ROUNDING_MODE[end])
    return z
end
function ldexp(x::BigFloat, n::Culong)
    z = BigFloat()
    ccall((:mpfr_mul_2ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, n, ROUNDING_MODE[end])
    return z
end
ldexp(x::BigFloat, n::ClongMax) = ldexp(x, convert(Clong, n))
ldexp(x::BigFloat, n::CulongMax) = ldexp(x, convert(Culong, n))
ldexp(x::BigFloat, n::Integer) = x*exp2(BigFloat(n))

function besselj0(x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_j0, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
    return z
end

function besselj1(x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_j1, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
    return z
end

function besselj(n::Integer, x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_jn, :libmpfr), Int32, (Ptr{BigFloat}, Clong, Ptr{BigFloat}, Int32), &z, n, &x, ROUNDING_MODE[end])
    return z
end

function bessely0(x::BigFloat)
    if x < 0
        throw(DomainError())
    end
    z = BigFloat()
    ccall((:mpfr_y0, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
    return z
end

function bessely1(x::BigFloat)
    if x < 0
        throw(DomainError())
    end
    z = BigFloat()
    ccall((:mpfr_y1, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
    return z
end

function bessely(n::Integer, x::BigFloat)
    if x < 0
        throw(DomainError())
    end
    z = BigFloat()
    ccall((:mpfr_yn, :libmpfr), Int32, (Ptr{BigFloat}, Clong, Ptr{BigFloat}, Int32), &z, n, &x, ROUNDING_MODE[end])
    return z
end

function factorial(x::BigFloat)
    if x < 0 || !isinteger(x)
        throw(DomainError())
    end
    ui = convert(Culong, x)
    z = BigFloat()
    ccall((:mpfr_fac_ui, :libmpfr), Int32, (Ptr{BigFloat}, Culong, Int32), &z, ui, ROUNDING_MODE[end])
    return z
end

function hypot(x::BigFloat, y::BigFloat)
    z = BigFloat()
    ccall((:mpfr_hypot, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[end])
    return z
end

for f in (:log, :log2, :log10)
    @eval function $f(x::BigFloat)
        if x < 0
            throw(DomainError())
        end
        z = BigFloat()
        ccall(($(string(:mpfr_,f)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
        return z
    end
end

function log1p(x::BigFloat)
    if x < -1
        throw(DomainError())
    end
    z = BigFloat()
    ccall((:mpfr_log1p, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
    return z
end

function max(x::BigFloat, y::BigFloat)
    z = BigFloat()
    ccall((:mpfr_max, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[end])
    return z
end

function min(x::BigFloat, y::BigFloat)
    z = BigFloat()
    ccall((:mpfr_min, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[end])
    return z
end

function modf(x::BigFloat)
    if isinf(x)
        return (BigFloat(NaN), x)
    end
    zint = BigFloat()
    zfloat = BigFloat()
    ccall((:mpfr_modf, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &zint, &zfloat, &x, ROUNDING_MODE[end])
    return (zfloat, zint)
end

function rem(x::BigFloat, y::BigFloat)
    z = BigFloat()
    ccall((:mpfr_fmod, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[end])
    return z
end

function sum(arr::AbstractArray{BigFloat})
    z = BigFloat(0)
    for i in arr
        ccall((:mpfr_add, :libmpfr), Int32,
            (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Cint),
            &z, &z, &i, 0)
    end
    return z
end

# Functions for which NaN results are converted to DomainError, following Base
for f in (:sin,:cos,:tan,:sec,:csc,
          :acos,:asin,:atan,:acosh,:asinh,:atanh, :gamma)
    @eval begin
        function ($f)(x::BigFloat)
            if isnan(x)
                return x
            end
            z = BigFloat()
            ccall(($(string(:mpfr_,f)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
            if isnan(z)
                throw(DomainError())
            end
            return z
        end
    end
end

# log of absolute value of gamma function
const lgamma_signp = Array(Cint, 1)
function lgamma(x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_lgamma,:libmpfr), Cint, (Ptr{BigFloat}, Ptr{Cint}, Ptr{BigFloat}, Int32), &z, lgamma_signp, &x, ROUNDING_MODE[end])
    return z
end

lgamma_r(x::BigFloat) = (lgamma(x), lgamma_signp[1])

function atan2(y::BigFloat, x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_atan2, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &y, &x, ROUNDING_MODE[end])
    return z
end

# Utility functions
==(x::BigFloat, y::BigFloat) = ccall((:mpfr_equal_p, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &x, &y) != 0
<=(x::BigFloat, y::BigFloat) = ccall((:mpfr_lessequal_p, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &x, &y) != 0
>=(x::BigFloat, y::BigFloat) = ccall((:mpfr_greaterequal_p, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &x, &y) != 0
<(x::BigFloat, y::BigFloat) = ccall((:mpfr_less_p, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &x, &y) != 0
>(x::BigFloat, y::BigFloat) = ccall((:mpfr_greater_p, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &x, &y) != 0

function ==(i::BigInt, f::BigFloat)
    !isnan(f) && ccall((:mpfr_cmp_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigInt}), &f, &i) == 0
end

==(f::BigFloat, i::BigInt) = i == f

function <(i::BigInt, f::BigFloat)
    # note: mpfr_cmp_z returns 0 if isnan(f)
    ccall((:mpfr_cmp_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigInt}), &f, &i) > 0
end

function <(f::BigFloat, i::BigInt)
    # note: mpfr_cmp_z returns 0 if isnan(f)
    ccall((:mpfr_cmp_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigInt}), &f, &i) < 0
end

<=(i::BigInt, f::BigFloat) = !isnan(f) && !(f < i)
<=(f::BigFloat, i::BigInt) = !isnan(f) && !(i < f)

signbit(x::BigFloat) = ccall((:mpfr_signbit, :libmpfr), Int32, (Ptr{BigFloat},), &x) != 0

function precision(x::BigFloat)
    return ccall((:mpfr_get_prec, :libmpfr), Clong, (Ptr{BigFloat},), &x)
end

get_bigfloat_precision() = DEFAULT_PRECISION[end]
function set_bigfloat_precision(x::Int)
    if x < 2
        throw(DomainError())
    end
    DEFAULT_PRECISION[end] = x
end

maxintfloat(x::BigFloat) = BigFloat(2)^precision(x)
maxintfloat(::Type{BigFloat}) = BigFloat(2)^get_bigfloat_precision()

function to_mpfr(r::RoundingMode)
    c = r.code
    if !(0 <= c <= 4)
        error("invalid BigFloat rounding mode")
    end
    c
end

function from_mpfr(c::Integer)
    if !(0 <= c <= 4)
        error("invalid MPFR rounding mode code")
    end
    RoundingMode(c)
end

get_rounding(::Type{BigFloat}) = from_mpfr(ROUNDING_MODE[end])
set_rounding(::Type{BigFloat},r::RoundingMode) = ROUNDING_MODE[end] = to_mpfr(r)

function copysign(x::BigFloat, y::BigFloat)
    z = BigFloat()
    ccall((:mpfr_copysign, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[end])
    return z
end

function exponent(x::BigFloat)
    if x == 0 || !isfinite(x)
        throw(DomainError())
    end
    # The '- 1' is to make it work as Base.exponent
    return ccall((:mpfr_get_exp, :libmpfr), Clong, (Ptr{BigFloat},), &x) - 1
end

function isinteger(x::BigFloat)
    return ccall((:mpfr_integer_p, :libmpfr), Int32, (Ptr{BigFloat},), &x) != 0
end

for f in (:ceil, :floor, :trunc, :round)
    @eval begin
        function ($f)(x::BigFloat)
            z = BigFloat()
            ccall(($(string(:mpfr_,f)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &z, &x)
            return z
        end
    end
end

function itrunc(x::BigFloat)
    z = BigInt()
    ccall((:mpfr_get_z, :libmpfr), Int32, (Ptr{BigInt}, Ptr{BigFloat}, Int32), &z, &x, 0)
    return z
end

iround(x::BigFloat) = itrunc(round(x))

function isinf(x::BigFloat)
    return ccall((:mpfr_inf_p, :libmpfr), Int32, (Ptr{BigFloat},), &x) != 0
end

function isnan(x::BigFloat)
    return ccall((:mpfr_nan_p, :libmpfr), Int32, (Ptr{BigFloat},), &x) != 0
end

isfinite(x::BigFloat) = !isinf(x) && !isnan(x)

@eval inf(::Type{BigFloat}) = $(BigFloat(Inf))
@eval nan(::Type{BigFloat}) = $(BigFloat(NaN))

typemax(::Type{BigFloat}) = inf(BigFloat)
@eval typemin(::Type{BigFloat}) = $(BigFloat(-Inf))

function nextfloat(x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_set, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32),
          &z, &x, ROUNDING_MODE[end])
    ccall((:mpfr_nextabove, :libmpfr), Int32, (Ptr{BigFloat},), &z) != 0
    return z
end

function prevfloat(x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_set, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32),
          &z, &x, ROUNDING_MODE[end])
    ccall((:mpfr_nextbelow, :libmpfr), Int32, (Ptr{BigFloat},), &z) != 0
    return z
end

eps(::Type{BigFloat}) = nextfloat(BigFloat(1)) - BigFloat(1)

realmin(::Type{BigFloat}) = nextfloat(zero(BigFloat))
realmax(::Type{BigFloat}) = prevfloat(inf(BigFloat))

function with_bigfloat_precision(f::Function, precision::Integer)
    old_precision = get_bigfloat_precision()
    set_bigfloat_precision(precision)
    try
        return f()
    finally
        set_bigfloat_precision(old_precision)
    end
end

function string(x::BigFloat)
    lng = 128
    for i = 1:2
        z = Array(Uint8, lng + 1)
        lng = ccall((:mpfr_snprintf,:libmpfr), Int32, (Ptr{Uint8}, Culong, Ptr{Uint8}, Ptr{BigFloat}...), z, lng + 1, "%.Re", &x)
        if lng < 128 || i == 2
            return bytestring(z[1:lng])
        end
    end
end

print(io::IO, b::BigFloat) = print(io, string(b))
show(io::IO, b::BigFloat) = print(io, string(b), " with $(precision(b)) bits of precision")
showcompact(io::IO, b::BigFloat) = print(io, string(b))

end #module

max_stars_repo_name: vchuravy/julia

max_issues_repo_name: vchuravy/julia

max_forks_repo_name: vchuravy/julia


text:
using LightGraphs

function getcolor(color,i)
  if isempty(color)
    return ""
  else
    return "#$(hex(color[i]))"
  end
end

function getlabel(label,i)
  if isempty(label)
    return ""
  else
    return label[i]
  end
end

function getgroup(group,i)
  if isempty(group)
    return ""
  else
    return group[i]
  end
end

function drawGraph(g;label=[],color=[],group=[])
  n = [NodeType(i,"$(getlabel(label,i))","$(getcolor(color,i))","$(getgroup(group,i))") for i in 1:nv(g)]
  e = edges(g)
  VisJSGraph(n,e)
end


text:
include(find_in_path("PatternDispatch.jl"))

module TestDispatch
using PatternDispatch

@show @qpat (x::Int)::Real

@show (@qpat (x::Int)) == (@qpat (x::Int)::Real)
@show (@qpat (x::Int)) == (@qpat (y::Int)::Real)

@show (@qpat x) > (@qpat x::Int)
@show (@qpat x) <= (@qpat x::Int)

@show (@qpat x::Int::String)

end # module


text:
using Test
using CrypticCrosswords
using CrypticCrosswords: definition

@testset "Known clues" begin
    known_clues = [
        ("Ach, Cole wrecked something in the ear", 7, "something in the ear", "cochlea"),
        ("aerial worker anne on the way up", 7, "aerial", "antenna"),
        ("at first congoers like us eschew solving hints", 5, "hints", "clues"),
        ("attractive female engraving", 8, "attractive", "fetching"),
        ("canoe wrecked in large sea", 5, "large sea", "ocean"),
        ("Carryall's gotta be upset", 7, "carryalls", "tote bag"),
        ("couch is unfinished until now", 4, "couch", "sofa"),
        ("cuts up curtains differently for those who use needles", 14, "those who use needles", "acupuncturists"),
        ("Desire bawdy slut", 4, "desire", "lust"),
        ("Dotty, Sue, Pearl, Joy", 8, "joy", "pleasure"),
        ("Endlessly long months and months", 4, "months and months", "year"),
        ("excitedly print Camus document", 10, "document", "manuscript"),
        ("Father returning ring with charm", 6, "charm", "appeal"),
        ("form of licit sea salt", 8, "salt", "silicate"),
        ("improve meal or I eat nuts", 10, "improve", "ameliorate"),
        ("initial meetings disappoint rosemary internally", 6, "initial meetings", "intros"),
        ("Initially congoers like us eschew solving hints", 5, "hints", "clues"),
        ("initially babies are naked", 4, "naked", "bare"),
        ("it's lunacy for dam to back onto ness", 7, "its lunacy", "madness"),
        ("hungary's leader, stuffy and bald", 8, "bald", "hairless"),
        ("male done mixing drink", 8, "drink", "lemonade"),
        ("measuring exotic flowers", 9, "flowers", "geraniums"),
        ("model unusually creepy hat", 9, "model", "archetype"),
        ("mollify with fried sausage", 7, "mollify", "assuage"),
        ("M's Rob Titon pitching slider?", 10, "slider", "trombonist"),
        ("Orchestra: I'm reorganizing conductor", 11, "conductor", "choirmaster"),
        ("Partially misconstrue fulminations; sorry", 6, "sorry", "rueful"),
        ("Propane explodes, theoretically", 7, "theoretically", "on paper"),
        ("Reap pleasure holding fruit", 5, "fruit", "apple"),
        ("Recover via fantastic miracle", 7, "recover", "reclaim"),
        ("returning regal drink", 5, "drink", "lager"),
        ("she literally describes high society", 5, "high society", "elite"),
        ("Significant ataxia overshadows choral piece", 7, "piece", "cantata"), # definition should actually be "choral piece"
        ("signore redefined districts", 7, "districts", "regions"),
        ("Sing gist of laudatory ode loudly", 5, "sing", "yodel"),
        ("singers in special tosca production", 5, "singers", "altos"),
        ("sink graduate with sin", 5, "sink", "basin"),
        ("spin broken shingle", 7, "spin", "english"),
        ("St. Michael transforms metal transformer", 9, "transformer", "alchemist"), # should be "metal transformer"
        ("stirs, spilling soda", 4, "stirs", "ados"),
        ("surprisingly rank height as important", 12, "important", "earthshaking"),
        ("they primarily play Diplomacy", 4, "diplomacy", "tact"),
        ("trimmed complicated test", 7, "test", "midterm"),
    ]

    badly_ranked_clues = [
        ("in glee over unusual color", 10, "color", "olive green"), # TODO: this is solvable, but we get "green olive" equally highly ranked
        ("anagram marvellously conceals structure of language", 7, "language", "grammar"),
        ("clean oneself, but in reverse", 3, "clean oneself", "tub"),
        ("Damaged credential tied together", 10, "tied together", "interlaced"),
        ("during exam I diagrammed viscera", 4, "during", "amid"),
        ("fish or insect for captain", 7, "fish or insect", "skipper"),
        ("figure out price he'd restructured", 8, "figure out", "decipher"),
        ("Inherently helps students over here", 4, "over here", "psst"),
        ("made mistake in deer reduction", 5, "made mistake", "erred"),
        # ("join trio of astronomers in marsh", 6, "join", "fasten"), # TODO: fix these clues. Weirdly low grammar scores.
        # ("sat up, interrupting sibling's balance", 6, "balance", "stasis"),
        ("setting for a cello composition", 6, "setting", "locale"),
        ("small bricks included among durable goods", 4, "small bricks", "lego"),
        ("waste pores vent exhausted resources", 9, "exhausted resources", "overspent"),

    ]

    @time for (clue, length, expected_definition, expected_wordplay) in known_clues
        @show clue
        solutions, state = @time solve(clue, length=length)
        arc = first(solutions)
        @test definition(arc) == expected_definition || endswith(expected_definition, definition(arc))
        @test arc.output == expected_wordplay
        derivations = Iterators.flatten([derive!(state, s) for s in Iterators.take(solutions, 10)])
    end

    @time for (clue, length, expected_definition, expected_wordplay) in badly_ranked_clues
        @show clue
        solutions, state = @time solve(clue, length=length)
        @test any(solutions) do arc
            definition(arc) == expected_definition && arc.output == expected_wordplay
        end
        derivations = Iterators.flatten([derive!(state, s) for s in Iterators.take(solutions, 10)])
    end
end


text:
"""
The {2, 3, 5}-wheel is efficient because it compresses memory
perfectly (1 byte = 30 numbers) and removes 4/15th of all the
multiples already. We don't get that memory efficiency when
extending the wheel to {2, 3, 5, 7}, since we need to store
48 bits per 210 numbers, which is could be done with one 64-bit
integer per 210 numbers, which in fact compresses worse
(64 bits / 210 numbers) than the {2, 3, 5}-wheel
(8 bits / 30 numbers).

What we can do however, is compute the repeating pattern that
the first `n` primes create, and copy that pattern over. That
is, we look at a the numbers modulo p₁ * p₂ * ⋯ * pₙ.

For instance, when presieving all multiples of {2, 3, ..., 19}
we allocate a buffer for the range 1 : 2 * 3 * ... * 19 =
1:9_699_690. In a {2, 3, 5} wheel this means a buffer of
9_699_690 ÷ 30 = 323_323 bytes.
"""
function create_presieve_buffer()
    n_bytes = 7 * 11 * 13 * 17
    xs = fill(0xFF, n_bytes)

    @inbounds for p in (7, 11, 13, 17)
        p²        = p * p
        byte_idx  = p² ÷ 30 + 1
        wheel     = to_idx(p)
        wheel_idx = 8 * (wheel - 1) + wheel
        increment = 0
        @sieve_loop :unroll
    end

    @inbounds xs[1] = 0b11100001 # remove 7, 11, 13 and 17
    return xs
end

"""
When applying the presieve buffer, we have to compute the offset in
"""
function apply_presieve_buffer!(xs::Vector{UInt8}, buffer::Vector{UInt8}, byte_start, byte_stop)

    len = byte_stop - byte_start + 1

    # todo, clean this up a bit.
    from_idx = (byte_start - 1) % length(buffer) + 1
    to = min(len, length(buffer) - from_idx + 1)

    # First copy the remainder of buffer at the front
    copyto!(view(xs, Base.OneTo(to)), view(buffer, from_idx:from_idx + to - 1))
    from = to + 1

    # Then copy buffer multiple times
    while from + length(buffer) - 1 <= len
        copyto!(view(xs, from : from + length(buffer) - 1), buffer)
        from += length(buffer)
    end

    # And finally copy the remainder of buffer again
    copyto!(view(xs, from:len), view(buffer, Base.OneTo(length(from:len))))

    xs
end


text:
import nctools.Ipaper: str_extract, str_extract_all

# match(Regex(), basename(file)).match
"""
    get_model(file; prefix = "day_", postfix = "_hist|_ssp|_piControl")
"""
function get_model(file, prefix = "day_", postfix = "_hist|_ssp|_piControl")
    str_extract(basename(file), "(?<=$prefix).*(?=$postfix)")
end

function get_ensemble(file, pattern::AbstractString = "(?<=_)r\\d.*(?=_\\d{4,8})")
    str_extract(basename(file), pattern)
end

function get_scenario(file, pattern::AbstractString = "[a-z,A-Z,0-9]*(?=_r\\d)")
    str_extract(basename(file), pattern)
end

# get date_begin and date_end from the file name
function get_date(file::AbstractString, pattern::AbstractString = "[0-9]{4,8}")
    str_extract_all(basename(file), pattern)
end

function get_date(files::Vector{<:AbstractString}, pattern::AbstractString = "[0-9]{4,8}")
    dates = map(x -> get_date(x, pattern), files)
    date_begin = map(x -> x[1], dates)
    date_end = map(x -> x[2], dates)
    date_begin, date_end
end

function get_date_nmiss(file)
    dates = nc_date(file)
    dates_nmiss(dates)
end

"""
    CMIPFiles_info(files)

> Note: currently, only works for daily scale

# Return
- `model`:
- `ensemble`:
- `date_begin`, `date_end`:
- `file`:
"""
function CMIPFiles_info(files)
    date_begin, date_end = get_date(files)
    calender = nc_calendar.(files)
    cell_x, cell_y, regular = nc_cellsize(files)

    DataFrame(
        model = get_model.(files),
        ensemble = get_ensemble.(files),
        scenario = get_scenario.(files),
        date_begin = date_begin, date_end = date_end,
        calender = calender,
        nmiss = get_date_nmiss.(files), # v0.1.2
        cell_x = cell_x, cell_y = cell_y, grid_regular = regular, file = files)
end

export get_model, get_ensemble, get_scenario, get_date, get_date_nmiss, CMIPFiles_info


text:
module ImageClipboard

using ImageIO, FileIO, ColorTypes

export clipboard_img

# linux
include("_xclip.jl")
# # mac
# include("_osascript.jl")
# # windows
include("_powershell.jl")

"""
    clipboard_img() -> Matrix{<:Colorant}

Paste an image from clipboard
"""
function clipboard_img()
    if Sys.islinux()
        img = _xclip()
    elseif Sys.iswindows()
        img = _powershell()
    else
        error("Currently, only linux and windows are supported")
    end
    return img
end

"""
    clipboard_img(img::Matrix{<:Colorant})

Copy an image to clipboard
"""
function clipboard_img(img::Matrix{<:Colorant})
    if Sys.islinux()
        _xclip(img)
    elseif Sys.iswindows()
        _powershell(img)
    else
        error("Currently, only linux and windows are supported")
    end
    nothing
end

end # module


text:
"""
**`Helium.getmat`** -*Function*.

    `Helium.getmat(source, DataType, numRows::Int64, numCols::Int64,
    isHeader::Bool, isRowNames::BoolnumRows)` => `matrix, column names, row names`

Read the data matrix and extract column and row names if they exist.

"""
function getmat(he::HeAttributes)

    # Adjust number of columns and rows according to skip option and if there
    # exist column and/or row names.

    #######################
    # Indices adjustement #
    #######################

    if he.hasColNames
        he.numRows -= 1
    end

    if he.hasRowNames
        he.numCols -= 1
    end

    ##################
    # Pre-allocation #
    ##################

    mat = Array{he.matType}(undef, he.numRows, he.numCols)
    if he.hasRowNames
        rowNames = Array{String}(undef, he.numRows, 1)
    else
        rowNames = [""]
    end

    ######################
    # Fill-up the matrix #
    ######################

    global i = 0
    open(he.fileName) do io

        if he.hasColNames
            global colNames = split(readline(io), he.sep)
        else
            global colNames = [""]
        end

        if he.hasRowNames
            while !eof(io)
                tmp = [readline(io)]
                newLine = eval(Meta.parse(string("""rslt=split(replace($tmp[1], r"$(he.strMiss)|NA|MISSING"i => "NaN"), "$(he.sep)")""")))
                global i += 1
                global rowNames[i] = newLine[1]
                global mat[i, :] = parse.(he.matType, newLine[2:end])
            end
            if he.hasColNames
                colNames = colNames[2:end]
            end
        else
            while !eof(io)
                tmp = [readline(io)]
                newLine = eval(Meta.parse(string("""rslt=split(replace($tmp[1], r"$(he.strMiss)|NA|MISSING"i => "NaN"), "$(he.sep)")""")))
                global i += 1
                global mat[i, :] = parse.(he.matType, newLine[1:end])
            end
        end
    end

    return mat, string.(colNames)[:], string.(rowNames)[:], ["" ""]
end

"""
**`Helium.getskipmat`** -*Function*.

    `Helium.getskipmat(source, DataType, numRows::Int64, numCols::Int64,
    isHeader::Bool, isRowNames::BoolnumRows)` => `matrix, column names, row names`

Read the data matrix and extract column and row names if they exist.

"""
function getskipmat(he::HeAttributes)

    # Adjust number of columns and rows according to skip option and if there
    # exist column and/or row names.

    #######################
    # Indices adjustement #
    #######################

    if he.hasColNames
        he.numRows -= 1
    end

    if he.hasRowNames
        he.numCols -= 1 + he.skipCol
        global idxStartMat = 2 + he.skipCol
        global idxEndXtra = idxStartMat - 1
    else
        he.numCols -= he.skipCol
        global idxStartMat = 1 + he.skipCol
        global idxEndXtra = idxStartMat - 1
    end

    ##################
    # Pre-allocation #
    ##################

    mat = Array{he.matType}(undef, he.numRows, he.numCols)

    if he.hasRowNames
        rowNames = Array{String}(undef, he.numRows, 1)
    else
        rowNames = [""]
    end

    global matXtra = Array{String}(undef, he.numRows+he.hasColNames, he.skipCol)

    ######################
    # Fill-up the matrix #
    ######################

    global i = 0
    open(he.fileName) do io

        if he.hasColNames
            global colNames = string.(split(readline(io), he.sep)[:])
            global matXtra[1, :] = colNames[(1+he.hasRowNames):idxEndXtra]
            colNames = colNames[idxStartMat:end]
        else
            global colNames = [""]
        end

        if he.hasRowNames
            while !eof(io)
                global i += 1

                tmp = [readline(io)]
                tmp = split(tmp[1], he.sep)

                global rowNames[i] = tmp[1]

                global matXtra[i+he.hasColNames, :] = tmp[2:idxEndXtra]

                newLine = eval(Meta.parse(string("""replace.($tmp[$idxStartMat:end], r"$(he.strMiss)|NA|MISSING"i => "NaN")""")))
                global mat[i, :] = parse.(he.matType, newLine)
            end
        else
            while !eof(io)
                global i += 1

                tmp = [readline(io)]
                tmp = split(tmp[1], he.sep)

                global matXtra[i+he.hasColNames, :] = tmp[1:idxEndXtra]

                newLine = eval(Meta.parse(string("""replace.($tmp[idxStartMat:end], r"$(he.strMiss)|NA|MISSING"i => "NaN")""")))
                global mat[i, :] = parse.(he.matType, newLine)
            end
        end

    end

    return mat, colNames, string.(rowNames)[:], matXtra
end


text:
using Plots

t = range(0, 4π, length = 100)
r = range(1, 0, length = 100)

x = cos.(t) .* r
y = sin.(t) .* r

@gif for i in eachindex(x)
    scatter((x[i], y[i]), lims = (-1, 1), label = "")
end

@gif for i in eachindex(x)
    plot(x[1:i], y[1:i], lims = (-1, 1), label = "")
    scatter!((x[i], y[i]), color = 1, label = "")
end

@gif for i in eachindex(x)
    plot(x[1:i], y[1:i], alpha = max.((1:i) .+ 10 .- i, 0) / 10, lims = (-1, 1), label = "")
    scatter!((x[i], y[i]), color = 1, label = "")
end


text:
using PySide

w = Widget()
raise(w)


text:
type Tape
  tape::Array{Int, 1}
  pos::Int

  function Tape()
    new(Array(Int, 1), 1)
  end
end

inc(this::Tape) = (this.tape[this.pos] += 1)
dec(this::Tape) = (this.tape[this.pos] -= 1)
get(this::Tape) = this.tape[this.pos]

function advance(this::Tape)
  this.pos += 1
  if this.pos > length(this.tape)
    push!(this.tape, 0)
  end
end

function devance(this::Tape)
  if this.pos > 1
    this.pos -= 1
  end
end

validbfsymbol(x) = in(x, ['>', '<', '+', '-', '.', ',', '[', ']'])

type Program
  code::ASCIIString
  bracket_map::Dict{Int, Int}

  function Program(text)
    code = filter(validbfsymbol, text)
    bracket_map = Dict{Int, Int}()
    stack = Int[]
    pc = 1
    for ch in code
      if ch == '['
        push!(stack, pc)
      elseif ch == ']' && length(stack) > 0
        right = pop!(stack)
        bracket_map[pc] = right
        bracket_map[right] = pc
      end
      pc += 1
    end
    return new(code, bracket_map)
  end
end

function run(this::Program)
  code = this.code
  bracket_map = this.bracket_map
  tape = Tape()
  pc = 1
  while pc <= length(code)
    ch = code[pc]
    if ch == '+'
      inc(tape)
    elseif ch == '-'
      dec(tape)
    elseif ch == '>'
      advance(tape)
    elseif ch == '<'
      devance(tape)
    elseif ch == '['
      if get(tape) == 0
        pc = bracket_map[pc]
      end
    elseif ch == ']'
      if get(tape) != 0
        pc = bracket_map[pc]
      end
    elseif ch == '.'
      print(char(get(tape)))
    end
    pc += 1
  end
end

function main()
  text = open(readall, ARGS[1])
  p = Program(text)
  run(p)
end

main()

max_stars_repo_name: japgolly/kostya-benchmarks

max_issues_repo_name: japgolly/kostya-benchmarks

max_forks_repo_name: japgolly/kostya-benchmarks


text:
"""
    _read_rmf(path_rmf::String)

Reads an RMF calibration file (from HEASARC caldb), loads in energy bands and PI channels
for use when filtering events out of a good energy range

Returns the PI channels, and the min/max good energy ranges
"""
function _read_rmf(path_rmf::String)
    if !isfile(path_rmf)
        throw(SystemError("opening file RMF file `$path_rmf`", 2))
    end

    fits_file = FITSIO.FITS(path_rmf)

    fits_ebounds = fits_file["EBOUNDS"]

    pis = read(fits_ebounds, "CHANNEL")

    e_mins = read(fits_ebounds, "E_MIN")
    e_maxs = read(fits_ebounds, "E_MAX")

    close(fits_file)

    return pis, e_mins, e_maxs
end

function _read_rmf(mission_name::Symbol)
    path_rmf = config(mission_name).path_rmf

    return _read_rmf(path_rmf)
end

"""
    _calibrate_pis(pis::Union{Array,Arrow.Primitive{Int16}}, path_rmf::String)

Loads the RMF calibration data, creates PI channels for energy conversion

Channel bounds are the average of the min and max energy range
"""
function _calibrate_pis(pis::Union{Array,Arrow.Primitive{Int16},Arrow.Primitive{Int64}}, path_rmf::String)
    calp, calEmin, calEmax = _read_rmf(path_rmf)

    return map(x -> (calEmin[x+1] + calEmax[x+1])/2, pis)
end

function _calibrate_pis(pis::Union{Array,Arrow.Primitive{Int16},Arrow.Primitive{Int64}}, mission::Mission)
    path_rmf = mission_paths(mission)[:rmf]

    return _calibrate_pis(pis, path_rmf)
end

function _save_calibrated(mission::Mission, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame}, instrument::Symbol,
        feather_dir::String, calibrated_energy::DataFrames.DataFrame;
        log=true)

    calibrated_file_path = joinpath(feather_dir, "$(instrument)_calib.feather")

    Feather.write(calibrated_file_path, calibrated_energy)

    if log
        _log_add(mission, obs_row,
            Dict("data" =>
                Dict(:feather_cl =>
                    Dict(instrument =>
                        Dict(
                            :path_calib => calibrated_file_path
                        )
                    )
                )
            )
        )
    end

    @info "Saved '$instrument' calib files to $feather_dir"
end

function _read_calibration(path_calib::String)
    calib = Feather.read(path_calib)
end

"""
    calibrate(mission_name::Symbol, obs_row::DataFrames.DataFrame)

Loads in the calibrated event data, as well as the mission calibration RMF file,
then filters the events by the energy ranges/PI channels in the RMF file

Saves the calibrated files as a `calib.feather` if none exists

Loads `calib.feater` file if it does exist

Returns a calibrated (filtered to contain only good energies) `InstrumentData` type
"""
function calibrate(mission::Mission, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame};
        instrument_data::Dict{Symbol,JAXTAM.InstrumentData}=Dict{Symbol,InstrumentData}(), overwrite=false)
    obsid       = obs_row[:obsid]
    instruments = _mission_instruments(mission)
    cl_files    = _log_query(mission, obs_row, "data", :feather_cl)

    if !haskey(instrument_data, instruments)
        instrument_data = read_cl(mission, obs_row)
        # Reload log after read_cl finishes
        cl_files        = _log_query(mission, obs_row, "data", :feather_cl)
    end

    calibration_instrument_data = Dict{Symbol,DataFrames.DataFrame}()
    for instrument in instruments
        if haskey(cl_files[instrument], :path_calib) && !overwrite # Check if log has calibration file path
            @info "Loading feather calib files for $instrument"
            calibration_instrument_data[instrument] = _read_calibration(cl_files[instrument][:path_calib])
        else
            @info "Generating calib files for $instrument"
            calibration_instrument_data[instrument] = DataFrame(E=_calibrate_pis(instrument_data[instrument].events[:PI], mission))

            cl_path = abspath(_obs_path_local(mission, obs_row; kind=:jaxtam), "JAXTAM", _log_query_path(; category=:data, kind=:feather_cl))
            _save_calibrated(mission, obs_row, instrument, cl_path, calibration_instrument_data[instrument])
        end

        instrument_data[instrument].events[:E] = calibration_instrument_data[instrument][:E]
    end

    return instrument_data
end

function calibrate(mission::Mission, append_df::DataFrames.DataFrame, obsid::String; overwrite=false)
    obs_row = master_query(append_df, :obsid, obsid)

    return calibrate(mission, obs_row; overwrite=overwrite)
end

function calibrate(mission::Mission, obsid::String; overwrite=false)
    return calibrate(mission, master(mission), obsid; overwrite=overwrite)
end


text:
type PivotMetadata
  rowNames::Vector{String}
  colNames::Vector{String}
  rowLevels::Vector{Vector{String}}
  colLevels::Vector{Vector{String}}
  statisticLevels::Vector
  statisticInCols::Bool

  function PivotMetadata(rowNames::Vector{String}, colNames::Vector{String},
                         rowLevels, colLevels, statisticNames, statisticInCols)
    nRowLevel, nColLevel = length(rowLevels), length(colLevels)
    nRows,     nCols     = length(rowNames),  length(colNames)

    nRowLevel == nRows ? nothing :
      error("dim mismatch. nrowLevels != nrows, $mRowLevel != $nRows")
    nColLevel == nCols ? nothing :
      error("dim mismatch. ncolLevels != ncols, $nColLevel != $mCols")

    return new(rowNames, colNames, rowLevels, colLevels, statisticNames, statisticInCols)
  end
end
incols(pm::PivotMetadata) = pm.statisticInCols
function rowNames(pm::PivotMetadata)
  if incols(pm)
    return pm.rowNames
  else
    res = copy(pm.rowNames)
    append!(res, ["stats"])
    return res
  end
end
function colNames(pm::PivotMetadata)
  if incols(pm)
    res = copy(pm.colNames)
    append!(res, ["stats"])
    return res
  else
    return pm.colNames
  end
end
function rowLevels(pm::PivotMetadata)
  if incols(pm)
    return pm.rowLevels
  else
    res = copy(pm.rowLevels)
    append!(res, Vector{String}[pm.statisticLevels])
    return res
  end
end
function colLevels(pm::PivotMetadata)
  if incols(pm)
    res = copy(pm.colLevels)
    append!(res, Vector{String}[pm.statisticLevels])
    return res
  else
    return pm.colLevels
  end
end

type Pivot
  data::Vector{Matrix}
  metadata::PivotMetadata

  function Pivot(data, metadata)
#     @assert size(data, 1) == size(metadata.rowLevels, 1)
#     @assert size(data, 2) == size(metadata.colLevels, 2)

    return new(data, metadata)
  end
end
metadata(p::Pivot) = p.metadata

function pivotstyles(p::Pivot)
  empty = [:html => [:class => "empty", :style => "border:0;"]]
  head = [:html => [:class => "head", :style => "text-align:center;"]]
  colhead = [:html => [:class => "subhead", :style => "text-align:center;"]]
  rowhead = [:html => [:class => "subhead", :style => "text-align:right;"]]
  datael = [:html => [:class => "datael", :style => "text-align:right;"]]

  n = length(p.data)
  styleelements = Array(AbstractPresentationStyleElement, n+5)
  styleelements[1:5] = [
    ExpandStyleElement(empty),
    ExpandStyleElement(head),
    UniformStyleElement(head),
    CartesianStyleElement(colhead),
    CartesianStyleElement(rowhead)
    ]
  for i=6:(n+5)
    styleelements[i] = UniformStyleElement(datael)
  end
  return styleelements
end

max_stars_repo_name: cpritcha/Pivots

max_issues_repo_name: cpritcha/Pivots

max_forks_repo_name: cpritcha/Pivots