text:
include("LabeledPolynomial.jl")

TOTAL_DEGREES = []

TOPREDUCTIONSF5 = 0
NORMALFORMS = 0

RR_F5_DETECTED_USELESS_NF = 0
RR_F5_USELESS_NF   = 0
RR_F5_NF        = 0
RR_F5_BASIS_SIZE      = 0

function rr_zerof5()
    global RR_F5_DETECTED_USELESS_NF
    global RR_F5_USELESS_NF
    global RR_F5_NF
    global RR_F5_BASIS_SIZE
    RR_F5_DETECTED_USELESS_NF = 0
    RR_F5_USELESS_NF        = 0
    RR_F5_NF      = 0
    RR_F5_BASIS_SIZE   = 0
end

function katsuran(n; ground=QQ)
    R, x = PolynomialRing(ground, ["x$i" for i in 0:n])

    return [
        (sum(x[abs(l)+1]*x[abs(m-l)+1] for l=-n:n if abs(m-l)<=n) -
        x[m+1] for m=0:n-1)...,
        x[1] + 2sum(x[i+1] for i=1:n) - 1
    ]
end

using OffsetArrays

debug() = false

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

without_k(arr, k) = arr[1:length(arr) .!= k]

function interreduce_basis(F::Vector{AbstractAlgebra.Generic.MPoly{T}}) where {T}

    F = Vector{AbstractAlgebra.Generic.MPoly{T}}(F)

    @label start

    newF = empty(F)

    for (i, f) in enumerate(F)
        reducers = union(newF, F[i+1:end])
        _, nf = divrem(f, reducers)
        if !iszero(nf)
            push!(newF, nf)
        end
    end

    if F != newF
        F = newF
        @goto start
    end

    newF
end

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

function homogenize(F)
    base = parent(F[1])
    s = copy(symbols(base))
    push!(s, :h)
    R, xsu = PolynomialRing(base_ring(base), s, ordering=ordering(base))
    xs = xsu[1:end-1]
    u = xsu[end]

    homF = Vector{eltype(F)}(undef, 0)
    for f in F
        tdeg = total_degree(leading_monomial(f))
        f = evaluate(f, xs)
        homf = zero(R)
        for t in terms(f)
            homf += t * u^(tdeg - total_degree(t))
        end
        push!(homF, homf)
    end

    homF
end

function dehomogenize(homF)
    R = parent(homF[1])
    xs = collect(gens(R))
    xs[end] = one(R)

    F = Vector{eltype(homF)}(undef, 0)
    for homf in homF
        f = evaluate(homf, xs)
        push!(F, f)
    end

    F
end

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

function reduceby(f, g)
    newf = f
    leadg = leading_monomial(g)
    leadf = leading_monomial(f)
    success, u = divides(leadf, leadg)
    if success
        u = divexact(leadf, leadg)
        newf = f - u*g
    end
    return success, newf
end

function mynormalform(f, G, r)

    global NORMALFORMS
    NORMALFORMS += 1

    a, b = divrem(f, [r[g].ev for g in G])
    return b

    newf = f

    @label start
    updated = false
    for g in G
        success, newf = reduceby(newf, r[g].ev)
        updated = updated || success
    end

    if updated
        @goto start
    end

    newf
end

#-----------------------------------------------------------------------------
# Groebner basis things

function f5_total_degree_lead_cmp(x, y)
    t1, t2 = total_degree(x), total_degree(y)
    if t1 < t2
        return true
    elseif t1 == t2
        return x < y
    else
        return false
    end
end

# S-polynomial

# returns u, v,  such that
# u*a - v*b is spoly(a, b)
function mults(a, b)
    u = lcm(leading_monomial(a), leading_monomial(b))
    u = divexact(u, leading_term(a))

    v = lcm(leading_monomial(a), leading_monomial(b))
    v = divexact(v, leading_term(b))

    u, v
end

# S-polynomial of a and b
function spoly(a, b)
    u, v = mults(a, b)
    u*a - v*b
end

struct RewriteRule{Poly}
    index::Int
    monom::Poly
end

function addrule!(Rule::Vector{Vector{RewriteRule{Poly}}}, sgn, k) where {Poly}
    push!(Rule[sgn.index], RewriteRule{Poly}(k, sgn.monom))
end

function findrewriting(u, k, r, Rule)
    sgn = r[k].sgn

    rules = Rule[sgn.index]
    ctr = length(rules)

    debug() && @info "search for rewriter $rules" u k sgn*u

    while ctr >= 1
        rule = rules[ctr]
        flag, t = divides(u * sgn.monom, rule.monom)
        if flag
            debug() && @info "found rewriter of $(u * sgn) by $rule"
            return rule.index
        end
        ctr -= 1
    end

    return k
end

# changed
function findrewriting2(u, k, r, Rule)
    sgn = r[k].sgn

    rules = Rule[sgn.index]
    ctr = length(rules)

    @info "search for rewriter $rules" u k sgn*u

    println(rules)
    println("EEEE")

    while ctr >= 1
        rule = rules[ctr]
        flag, t = divides(u * sgn.monom, rule.monom)
        if flag
            debug() && @info "found rewriter of $(u * sgn) by $rule"
            if rule.index != k
                return rule.index
            end
        end
        ctr -= 1
    end

    return k
end

function isrewritable(u, k, r, Rule)
    j = findrewriting(u, k, r, Rule)
    j != k
end

function computespolys(P, r::Vector{Poly}, Rule) where {Poly}
    S = Int[]

    sort!(P, by= x -> x.t)

    debug() && @info "computespolys" P r Rule

    for p in P
        debug() && @info "spairs $p"

        # ADDED
        sgnk = r[p.k].sgn
        sgnl = r[p.l].sgn
        debug() && @info "is rewritable?? $(p.u) or $(p.v)"
        if isrewritable(p.u, p.k, r, Rule) || isrewritable(p.v, p.l, r, Rule)
            global RR_F5_DETECTED_USELESS_NF
            RR_F5_DETECTED_USELESS_NF += 1

            debug() && @info "rewritten $p"
            continue
        end

        s = spoly(r[p.k].ev, r[p.l].ev)
        sgns = r[p.k].sgn * p.u
        newlp = LabeledPolynomial(s, sgns)

        push!(r, newlp)

        addrule!(Rule, sgns, length(r))

        # @warn "Added" newlp
        # println(r)

        if !iszero(s)
            push!(S, length(r))
        end
    end

    sort!(S, by= x -> r[x].sgn)

    S
end

struct CriticalPair{Poly}
    t::Poly
    k::Int
    u::Poly
    l::Int
    v::Poly
end

# some magic to prin2t ModuleElement nicely
Base.show(io::IO, ::MIME"text/plain", m::CriticalPair) = println(io, "$((m.t, m.k, m.u, m.l, m.v))")
Base.show(io::IO, m::CriticalPair) = print(io, "$((m.t, m.k, m.u, m.l, m.v))")

function istopreducible(f, Gprev, r)
    t = leading_monomial(f)
    debug() && @info "istopreducible $f by " Gprev r[Gprev]
    for gi in Gprev
        g = r[gi].ev
        flag, _ = divides(t, leading_monomial(g))
        if flag
            debug() && @info "found top reducer" f gi r[gi]
            return true
        end
    end
    return false
end

function criticalpair(i, k, l, Gprev, r, Rule)
    tk = leading_monomial(r[k].ev)
    tl = leading_monomial(r[l].ev)

    t = lcm(tk, tl)

    u1 = divexact(t, tk)
    u2 = divexact(t, tl)

    sgn1 = r[k].sgn
    sgn2 = r[l].sgn

    # debug() && @info "for lps" i k l r[k] r[l]
    # debug() && @info "criticalpair reducers" sgn1 sgn2 Gprev

    debug() && @info "for lps" i k l r[k] r[l]
    debug() && @info "criticalpair reducers" sgn1 sgn2 Gprev

    usgn1 = u1*sgn1
    usgn2 = u2*sgn2

    if usgn1 == usgn2
        global RR_F5_DETECTED_USELESS_NF
        RR_F5_DETECTED_USELESS_NF += 1

        return CriticalPair(t, 0, u1, 0, u2)
    end

    @assert usgn2 != usgn1

    debug() && @info "multiplied sgns" usgn1 usgn2

    if sgn1.index == i && istopreducible(usgn1.monom, Gprev, r)
        global RR_F5_DETECTED_USELESS_NF
        RR_F5_DETECTED_USELESS_NF += 1

        debug() && @info "discarded $((k, l))" u1*sgn1
        return CriticalPair(t, 0, u1, 0, u2)
    end
    if sgn2.index == i && istopreducible(usgn2.monom, Gprev, r)
        global RR_F5_DETECTED_USELESS_NF
        RR_F5_DETECTED_USELESS_NF += 1

        debug() && @info "discarded $((k, l))" u2*sgn2
        return CriticalPair(t, 0, u1, 0, u2)
    end

    # @error "is??" usgn1.monom, l, r[l]

    # Added
    # if istopreducible(usgn1.monom, [l], r)
    #    @error "NEW discarded $((k, l))" u2*sgn2
    #    return CriticalPair(t, 0, u1, 0, u2)
    # end

    if isrewritable(u1, k, r, Rule) || isrewritable(u2, l, r, Rule)
        global RR_F5_DETECTED_USELESS_NF
        RR_F5_DETECTED_USELESS_NF += 1

        debug() && @info "rewritten $u1 or $u2"
        return CriticalPair(t, 0, u1, 0, u2)
    end

    if u1*sgn1 < u2*sgn2
        u1, u2 = u2, u1
        k, l = l, k
    end

    # debug() && @info "CriticalPair" t k u1 l u2
    debug() && @info "CriticalPair" t k u1 l u2

    CriticalPair(t, k, u1, l, u2)
end

function findreductor(k, Gprev, newGcurr, r, Rule)
    t = leading_monomial(r[k].ev)
    sgnk = r[k].sgn
    for j in newGcurr
        tj = leading_monomial(r[j].ev)
        flag, u = divides(t, tj)
        if flag
            sgnj = r[j].sgn
            sgnju = sgnj*u
            if sgnju.index != sgnk.index || sgnju.monom != sgnk.monom
                if !isrewritable(u, j, r, Rule)
                    if !istopreducible(sgnju.monom, Gprev, r)
                        return j
                    end
                end
            end
        end
    end
    return 0
end

function topreduction(k, Gprev, newGcurr, r, Rule)
    if iszero(r[k].ev)
        @error "Reduction to zero!"

        global RR_F5_USELESS_NF
        RR_F5_USELESS_NF += 1

        return 0, ()
    end

    j = findreductor(k, Gprev, newGcurr, r, Rule)

    (debug() || false) && @info "for k=$k reductor j=$j"

    if iszero(j)
        normalize!(r[k])
        return k, ()
    end

    p = r[k]
    q = r[j]

    (debug() || false) && println("reductor poly $q")

    u = divexact(leading_monomial(p.ev), leading_monomial(q.ev))
    c = leading_coefficient(p.ev) // leading_coefficient(q.ev)

    global TOPREDUCTIONSF5
    TOPREDUCTIONSF5 += 1

    newpev = p.ev - c*u*q.ev

    (debug() || false) && println("reductor res $p")

    if !iszero(newpev)
        newpev = normalize(newpev)
    end

    if u*r[j].sgn < r[k].sgn
        r[k].ev = newpev
        return 0, (k,)
    else
        newr = LabeledPolynomial(newpev, u*r[j].sgn)
        push!(r, newr)
        addrule!(Rule, newr.sgn, length(r))
        return 0, (k, length(r))
    end
end

function reduction(S, Gprev, Gcurr, r, Rule)
    todo = S
    completed = Int[]

    newGcurr = copy(Gcurr)

    debug() && @info "reducing $S with"
    if debug()
        println(Gprev, " // ", Gcurr)
    end

    if debug() || false
        println("r = $r")
    end

    while !isempty(todo)
        sort!(todo, by= x -> r[x].sgn)

        k = todo[1]
        todo = todo[2:end]

        (debug() || false) && @info "reduction, handle $k ($(r[k])) by $(r[Gprev])"

        global RR_F5_NF
        RR_F5_NF += 1
        r[k].ev = mynormalform(r[k].ev, Gprev, r)

        (debug() || false) && @info "from normal form" r[k]

        if debug() || false
            println("GOING TO TOPREDUCTION with Gprev = $Gprev and Gcurr=$newGcurr")
        end

        newly_computed, redo = topreduction(k, Gprev, newGcurr, r, Rule)

        (debug() || false) && @info "After top" r[k]

        (debug() || false) && @info "from topreduction" newly_computed redo
        (debug() || false) && @info "todo = $todo"

        if !iszero(newly_computed)
            push!(completed, newly_computed)
            push!(newGcurr, newly_computed)
        end

        for j in redo
            push!(todo, j)
        end

    end

    completed
end

function incremental_basis(i, Gprev, r::Vector{LabeledPolynomial{Poly}}, Rule) where {Poly}
    curr_idx = length(r)

    Gcurr = copy(Gprev)
    push!(Gcurr, curr_idx)

    push!(Rule, Vector{RewriteRule{Poly}}(undef, 0) )

    P = Vector{CriticalPair{Poly}}(undef, 0)
    for j in Gprev
        p = criticalpair(i, curr_idx, j, Gprev, r, Rule)
        if !iszero(p.k)
            push!(P, p)
        end
    end

    # global TOTAL_DEGREES

    debug() && @info "env" Gcurr
    @warn "$i iter, generated pairs" length(P) # P

    if debug() || false
        println("Pairs are:")
        for p in P
            println("($(p.k), $(p.l))")
        end
    end

    while !isempty(P)
        sort!(P, by= p -> total_degree(p.t))

        d = total_degree(P[1].t)
        Pd = Vector{CriticalPair{Poly}}(undef, 0)

        j  = 1
        while j <= length(P) && total_degree(P[j].t) == d
            push!(Pd, P[j])
            j += 1
        end
        # j = max(1, div(j, 5))
        # Pd = Pd[1:j]
        P = P[j:end]

        @warn "SELECTED $(length(Pd)) pairs of degree $d out of $(length(Pd) + length(P))"

        # println("Spolys: $Pd")
        # println("$(Pd[1].k)-th poly $(r[Pd[1].k])")
        # println("$(Pd[1].l)-th poly $(r[Pd[1].l])")

        if length(Pd) > 1
            # println("$(Pd[2].k)-th poly $(r[Pd[2].k])")
            # println("$(Pd[2].l)-th poly $(r[Pd[2].l])")
        end

        debug() && @info "selected pairs" d j Pd

        S = computespolys(Pd, r, Rule)

        if (debug() || false) && (length(S) > 0)
            println("# SPOLYS $(length(S))")
            println(S)
            println("First spoly: $(r[S[1]])")
        end

        R = reduction(S, Gprev, Gcurr, r, Rule)

        if (debug() || false) && (length(R) > 0)
            println("AFTER reduction $(length(R))")
            println(R)
            println("First reduced: $(r[R[1]])")
        end

        for k in R
            for j in Gcurr
                p = criticalpair(i, j, k, Gprev, r, Rule)
                if p.l != 0
                    push!(P, p)
                end
            end
            (debug() || true) && println("New # pairs for poly $k:", length(P))
            push!(Gcurr, k)
            (debug() || true) &&  println("New # Gcurr:", length(Gcurr))
        end

        if debug()
            println("###############")
            println("###############")
            println("In degree $d:")
            println("P = ", P)
            println("Gcurr = ", Gcurr)
            println("r = ", r)
            println("Rule =", Rule)
            println("###############")
            println("###############")
        end

        # push!(TOTAL_DEGREES[i], d)
    end

    if debug() || false
        println("ENDED with $Gcurr polynomials")
    end

    Gcurr
end

function getreducers(i, Gprev, r)
    ans = Gprev[1:(i-1)]
    ans = filter!(x -> !iszero(r[x].ev), append!(ans, Gprev[(i+1):end]))
    ans
end

function interreduction(Gprev, r)

    @label start
    updated = false

    for (i, idx) in enumerate(Gprev)
        reducers = getreducers(i, Gprev, r)
        newf = mynormalform(r[idx].ev, reducers, r)
        # @warn "" i idx newf r[idx].ev
        updated = updated || (newf != r[idx].ev)
        r[idx].ev = newf
    end

    if updated
        @goto start
    end

    res = []
    for i in Gprev
        if !iszero(r[i].ev)
            push!(res, r[i].ev)
        end
    end

    return res
end

function setup_reduced_basis(Gprev, r, Rule)
    B = interreduction(Gprev, r)

    Gcurr = collect(1:length(B))

    resize!(r, length(B))
    r[1:length(B)] = [
        LabeledPolynomial(B[i], Signature2(i, one(B[i])))
        for i in 1:length(B)
    ]

    resize!(Rule, length(B)+1)
    for i in 0:length(Rule)-1
        Rule[i] = []
    end

    for j in 1:length(B)-1
        t = leading_monomial(B[j])
        for k in j + 1:length(B)
            u = divexact(lcm(t, leading_monomial(B[k])), leading_monomial(B[k]))
            addrule!(Rule, Signature2(k, u), 0)
        end
    end

    return Gcurr
end

function signature_groebner_basis2(F::Vector{Poly}) where {Poly}
    sort!(F, by=leading_monomial, lt=f5_total_degree_lead_cmp)

    debug() && @warn "after initial sort" F

    # F = interreduce_basis(F)

    debug() && @warn "after interreduction" F

    m = length(F)

    # normalize f1
    lp1 = LabeledPolynomial(F[1], Signature2(1, one(F[1])))
    normalize!(lp1)

    r = [lp1]

    Gprev = [1]

    i = 2
    Rule = Vector{Vector{RewriteRule{Poly}}}(undef, 0)
    push!(Rule, Vector{RewriteRule{Poly}}(undef, 0))
    push!(Rule, Vector{RewriteRule{Poly}}(undef, 0))
    # Rule = OffsetArray(Rule, 0:1)

    # global TOTAL_DEGREES
    # TOTAL_DEGREES = []
    # push!(TOTAL_DEGREES, [])

    @debug "at the start" Gprev Rule

    while i <= m
        # push!(TOTAL_DEGREES, [])

        lpi = LabeledPolynomial(F[i], Signature2(i, one(F[i])))
        normalize!(lpi)
        push!(r, lpi)

        @debug "before incremental" r
        Gprev = incremental_basis(i, Gprev, r, Rule)

        # Gprev = setup_reduced_basis(Gprev, r, Rule)

        #println(Gprev)
        #println(r)
        #println(Rule)

        if debug()
            println("Gprev =", Gprev, r, Rule)
        end

        i += 1
    end

    if debug()
        println("##############")
        println("##############")
    end

    global RR_F5_BASIS_SIZE
    RR_F5_BASIS_SIZE = length(Gprev)

    # println(Gprev)
    # println(r[Gprev])
    B = interreduction(Gprev, r)

    sort(B, by=leading_monomial, rev=true)
end

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

using Logging
Logging.global_logger(ConsoleLogger(Logging.Info))

#=
R, (x,y,z) = PolynomialRing(QQ, ["x","y", "z"], ordering=:degrevlex)

F = [x*z + 1, y*z + 1]

# evaluation
G = signature_groebner_basis(F)
Gev = [g.ev for g in G]
println("############")
println(G)
println(Gev)
=#


text:
mutable struct SPPointMassConstraintOut <: SumProductRule{PointMassConstraint} end
outboundType(::Type{SPPointMassConstraintOut}) = Message{PointMass}
isApplicable(::Type{SPPointMassConstraintOut}, input_types::Vector{<:Type}) = (length(input_types) == 1) && (input_types[1] == Nothing)

ext: jl


text:
struct Config
    n_nodes::Int64
    n0::Int64
    k::Int64
    agentcount::Int64
    initial_trash::Int64
    initial_litterers::Int64
    garbage_collection::Int64
    n_iter::Int64
end

Config(;
    n_nodes::Int64, n0::Int64, k::Int64,
    agentcount::Int64, initial_trash::Int64,
    initial_litterers::Int64, garbage_collection::Int64,
    n_iter::Int64
) = Config(
    n_nodes, n0, k,
    agentcount, initial_trash,
    initial_litterers, garbage_collection,
    n_iter
)

Config() = Config(
    n_nodes=1_000, n0=3, k=3,
    agentcount=100, initial_trash=500,
    initial_litterers=20, garbage_collection=0,
    n_iter=1_000
)


text:
module CompositeDataFramesTimings

using DataArrays, DataFrames
using DataFramesMeta

srand(1)
const n = 5_000_000
a = rand(n)
b = rand(n)
cdf = CompositeDataFrame(a = a, b = b)
df = DataFrame(cdf)

function dot1(df::AbstractDataFrame)
    x = 0.0
    for i in 1:size(df, 1)
        x += df[:a][i] * df[:a][i]
    end
    return x
end

function dot2(df::AbstractDataFrame)
    x = 0.0
    for i in 1:size(df, 1)
        x += df[i,:a] * df[i,:a]
    end
    return x
end

function dot3(df::AbstractDataFrame)
    @with df begin
        x = 0.0
        for i in 1:length(:a)
            x += (:a)[i] * (:b)[i]
        end
        x
    end
end

@show t1 = @elapsed dot1(df)
@show t2 = @elapsed dot2(df)
@show t3 = @elapsed dot3(df)
@show t1c = @elapsed dot1(cdf)
@show t2c = @elapsed dot2(cdf)
@show t3c = @elapsed dot3(cdf)

end # module


text:
using StatsBase: mad

###############################################################################
# Location Estimators

@testset "Trivial $E"  for E in [MMMBackground, SourceExtractorBackground, BiweightLocationBackground]
    estimator = E()

    @test estimator(ones(1)) == 1

    data = ones(10, 10)

    @test estimator(data) ≈ 1.0
    @test estimator(data, dims = 1) ≈ ones(1, 10)
    @test estimator(data, dims = 2) ≈ ones(10)

    data = zeros(10, 10)

    @test estimator(data) ≈ 0.0
    @test estimator(data, dims = 1) ≈ zeros(1, 10)
    @test estimator(data, dims = 2) ≈ zeros(10)

    data = randn(100, 100)
end

@testset "SourceExtractorBackground" begin
    # test skewed distribution
    data = float(collect(1:100))
    data[71:end] .= 1e7

    @test SourceExtractorBackground()(data) ≈ median(data)
end

@testset "BiweightLocationBackground" begin
    b = BiweightLocationBackground()
    @test b([1, 3, 5, 500, 2]) ≈ 2.745 atol = 1e-3
end

###############################################################################
# RMS Estimators

@testset "Trivial $E"  for E in [StdRMS, MADStdRMS, BiweightScaleRMS]
    estimator = E()

    @test estimator(ones(1)) == 0

    data = ones(10, 10)
    @test estimator(data) ≈ 0.0
    @test estimator(data, dims = 1) ≈ zeros(1, 10)
    @test estimator(data, dims = 2) ≈ zeros(10)

    data = randn(10000, 10000)
    @test estimator(data) ≈ 1 atol = 3e-2
end

@testset "StdRMS" begin
    s = StdRMS()
    data = randn(100)
    @test s(data) == std(data, corrected = false)
end

@testset "MADStdRMS" begin
    s = MADStdRMS()
    data = randn(100)
    @test s(data) == mad(data, normalize = true)
end

@testset "BiweightScaleRMS" begin
    s = BiweightScaleRMS()
    @test s([1, 3, 5, 500, 2]) ≈ 1.70992562072
end


text:
#************************************************************************
# Global definitions
#************************************************************************

# Key and button actions
@enum Action::Cint begin
	RELEASE                = 0
	PRESS                  = 1
	REPEAT                 = 2
end

@enum Key::Cint begin
	# Unknown key
	KEY_UNKNOWN            = -1

	# Printable keys
	KEY_SPACE              = 32
	KEY_APOSTROPHE         = 39  # '
	KEY_COMMA              = 44  # ,
	KEY_MINUS              = 45  # -
	KEY_PERIOD             = 46  # .
	KEY_SLASH              = 47  # /
	KEY_0                  = 48
	KEY_1                  = 49
	KEY_2                  = 50
	KEY_3                  = 51
	KEY_4                  = 52
	KEY_5                  = 53
	KEY_6                  = 54
	KEY_7                  = 55
	KEY_8                  = 56
	KEY_9                  = 57
	KEY_SEMICOLON          = 59  # ;
	KEY_EQUAL              = 61  # =
	KEY_A                  = 65
	KEY_B                  = 66
	KEY_C                  = 67
	KEY_D                  = 68
	KEY_E                  = 69
	KEY_F                  = 70
	KEY_G                  = 71
	KEY_H                  = 72
	KEY_I                  = 73
	KEY_J                  = 74
	KEY_K                  = 75
	KEY_L                  = 76
	KEY_M                  = 77
	KEY_N                  = 78
	KEY_O                  = 79
	KEY_P                  = 80
	KEY_Q                  = 81
	KEY_R                  = 82
	KEY_S                  = 83
	KEY_T                  = 84
	KEY_U                  = 85
	KEY_V                  = 86
	KEY_W                  = 87
	KEY_X                  = 88
	KEY_Y                  = 89
	KEY_Z                  = 90
	KEY_LEFT_BRACKET       = 91  # [
	KEY_BACKSLASH          = 92  # \
	KEY_RIGHT_BRACKET      = 93  # ]
	KEY_GRAVE_ACCENT       = 96  # `
	KEY_WORLD_1            = 161 # non-US #1
	KEY_WORLD_2            = 162 # non-US #2

	# Function keys
	KEY_ESCAPE             = 256
	KEY_ENTER              = 257
	KEY_TAB                = 258
	KEY_BACKSPACE          = 259
	KEY_INSERT             = 260
	KEY_DELETE             = 261
	KEY_RIGHT              = 262
	KEY_LEFT               = 263
	KEY_DOWN               = 264
	KEY_UP                 = 265
	KEY_PAGE_UP            = 266
	KEY_PAGE_DOWN          = 267
	KEY_HOME               = 268
	KEY_END                = 269
	KEY_CAPS_LOCK          = 280
	KEY_SCROLL_LOCK        = 281
	KEY_NUM_LOCK           = 282
	KEY_PRINT_SCREEN       = 283
	KEY_PAUSE              = 284
	KEY_F1                 = 290
	KEY_F2                 = 291
	KEY_F3                 = 292
	KEY_F4                 = 293
	KEY_F5                 = 294
	KEY_F6                 = 295
	KEY_F7                 = 296
	KEY_F8                 = 297
	KEY_F9                 = 298
	KEY_F10                = 299
	KEY_F11                = 300
	KEY_F12                = 301
	KEY_F13                = 302
	KEY_F14                = 303
	KEY_F15                = 304
	KEY_F16                = 305
	KEY_F17                = 306
	KEY_F18                = 307
	KEY_F19                = 308
	KEY_F20                = 309
	KEY_F21                = 310
	KEY_F22                = 311
	KEY_F23                = 312
	KEY_F24                = 313
	KEY_F25                = 314
	KEY_KP_0               = 320
	KEY_KP_1               = 321
	KEY_KP_2               = 322
	KEY_KP_3               = 323
	KEY_KP_4               = 324
	KEY_KP_5               = 325
	KEY_KP_6               = 326
	KEY_KP_7               = 327
	KEY_KP_8               = 328
	KEY_KP_9               = 329
	KEY_KP_DECIMAL         = 330
	KEY_KP_DIVIDE          = 331
	KEY_KP_MULTIPLY        = 332
	KEY_KP_SUBTRACT        = 333
	KEY_KP_ADD             = 334
	KEY_KP_ENTER           = 335
	KEY_KP_EQUAL           = 336
	KEY_LEFT_SHIFT         = 340
	KEY_LEFT_CONTROL       = 341
	KEY_LEFT_ALT           = 342
	KEY_LEFT_SUPER         = 343
	KEY_RIGHT_SHIFT        = 344
	KEY_RIGHT_CONTROL      = 345
	KEY_RIGHT_ALT          = 346
	KEY_RIGHT_SUPER        = 347
	KEY_MENU               = 348
end

# Modifier key bit flags
const MOD_SHIFT              = 0x0001
const MOD_CONTROL            = 0x0002
const MOD_ALT                = 0x0004
const MOD_SUPER              = 0x0008

# Mouse buttons
@enum MouseButton::Cint begin
	MOUSE_BUTTON_1         = 0
	MOUSE_BUTTON_2         = 1
	MOUSE_BUTTON_3         = 2
	MOUSE_BUTTON_4         = 3
	MOUSE_BUTTON_5         = 4
	MOUSE_BUTTON_6         = 5
	MOUSE_BUTTON_7         = 6
	MOUSE_BUTTON_8         = 7
end

const MOUSE_BUTTON_LEFT      = MOUSE_BUTTON_1
const MOUSE_BUTTON_RIGHT     = MOUSE_BUTTON_2
const MOUSE_BUTTON_MIDDLE    = MOUSE_BUTTON_3

# Joysticks
@enum Joystick::Cint begin
	JOYSTICK_1             = 0
	JOYSTICK_2             = 1
	JOYSTICK_3             = 2
	JOYSTICK_4             = 3
	JOYSTICK_5             = 4
	JOYSTICK_6             = 5
	JOYSTICK_7             = 6
	JOYSTICK_8             = 7
	JOYSTICK_9             = 8
	JOYSTICK_10            = 9
	JOYSTICK_11            = 10
	JOYSTICK_12            = 11
	JOYSTICK_13            = 12
	JOYSTICK_14            = 13
	JOYSTICK_15            = 14
	JOYSTICK_16            = 15
end

# Error codes
@enum ErrorCode::Cint begin
	NOT_INITIALIZED        = 0x00010001  # GLFW has not been initialized.
	NO_CURRENT_CONTEXT     = 0x00010002  # No context is current for this thread.
	INVALID_ENUM           = 0x00010003  # One of the arguments to the function was an invalid enum value.
	INVALID_VALUE          = 0x00010004  # One of the arguments to the function was an invalid value.
	OUT_OF_MEMORY          = 0x00010005  # A memory allocation failed.
	API_UNAVAILABLE        = 0x00010006  # GLFW could not find support for the requested API on the system.
	VERSION_UNAVAILABLE    = 0x00010007  # The requested OpenGL or OpenGL ES version is not available.
	PLATFORM_ERROR         = 0x00010008  # A platform-specific error occurred that does not match any of the more specific categories.
	FORMAT_UNAVAILABLE     = 0x00010009  # The requested format is not supported or available.
	NO_WINDOW_CONTEXT      = 0x0001000A  # The specified window does not have an OpenGL or OpenGL ES context.
end

const FOCUSED                = 0x00020001
const ICONIFIED              = 0x00020002
const RESIZABLE              = 0x00020003
const VISIBLE                = 0x00020004
const DECORATED              = 0x00020005
const AUTO_ICONIFY           = 0x00020006
const FLOATING               = 0x00020007
const MAXIMIZED              = 0x00020008

const RED_BITS               = 0x00021001
const GREEN_BITS             = 0x00021002
const BLUE_BITS              = 0x00021003
const ALPHA_BITS             = 0x00021004
const DEPTH_BITS             = 0x00021005
const STENCIL_BITS           = 0x00021006
const ACCUM_RED_BITS         = 0x00021007
const ACCUM_GREEN_BITS       = 0x00021008
const ACCUM_BLUE_BITS        = 0x00021009
const ACCUM_ALPHA_BITS       = 0x0002100A
const AUX_BUFFERS            = 0x0002100B
const STEREO                 = 0x0002100C
const SAMPLES                = 0x0002100D
const SRGB_CAPABLE           = 0x0002100E
const REFRESH_RATE           = 0x0002100F
const DOUBLEBUFFER           = 0x00021010

const CLIENT_API             = 0x00022001
const CONTEXT_VERSION_MAJOR  = 0x00022002
const CONTEXT_VERSION_MINOR  = 0x00022003
const CONTEXT_REVISION       = 0x00022004
const CONTEXT_ROBUSTNESS     = 0x00022005
const OPENGL_FORWARD_COMPAT  = 0x00022006
const OPENGL_DEBUG_CONTEXT   = 0x00022007
const OPENGL_PROFILE         = 0x00022008
const CONTEXT_RELEASE_BEHAVIOR = 0x00022009
const CONTEXT_NO_ERROR       = 0x0002200A
const CONTEXT_CREATION_API   = 0x0002200B

const OPENGL_API             = 0x00030001
const OPENGL_ES_API          = 0x00030002
const NO_API                 =          0

const NO_ROBUSTNESS          =          0
const NO_RESET_NOTIFICATION  = 0x00031001
const LOSE_CONTEXT_ON_RESET  = 0x00031002

const OPENGL_ANY_PROFILE     =          0
const OPENGL_CORE_PROFILE    = 0x00032001
const OPENGL_COMPAT_PROFILE  = 0x00032002

const CURSOR                 = 0x00033001
const STICKY_KEYS            = 0x00033002
const STICKY_MOUSE_BUTTONS   = 0x00033003

const CURSOR_NORMAL          = 0x00034001
const CURSOR_HIDDEN          = 0x00034002
const CURSOR_DISABLED        = 0x00034003

const ANY_RELEASE_BEHAVIOR   =          0
const RELEASE_BEHAVIOR_FLUSH = 0x00035001
const RELEASE_BEHAVIOR_NONE  = 0x00035002

const NATIVE_CONTEXT_API     = 0x00036001
const EGL_CONTEXT_API        = 0x00036002

# Standard cursor shapes
@enum StandardCursorShape::Cint begin
	ARROW_CURSOR           = 0x00036001
	IBEAM_CURSOR           = 0x00036002
	CROSSHAIR_CURSOR       = 0x00036003
	HAND_CURSOR            = 0x00036004
	HRESIZE_CURSOR         = 0x00036005
	VRESIZE_CURSOR         = 0x00036006
end

@enum DeviceConfigEvent::Cint begin
	CONNECTED              = 0x00040001
	DISCONNECTED           = 0x00040002
end

const DONT_CARE              = -1

#************************************************************************
# GLFW API types
#************************************************************************

struct Monitor
	handle::Ptr{Cvoid}
end
Base.show(io::IO, m::Monitor) = write(io, m.handle != C_NULL ? "$(m.name) Monitor ($(m.videomode))" : "No Monitor")

mutable struct Window
	handle::Ptr{Cvoid}
end

#Came from GLWindow.jl/screen.jl
"""
Function to create a pure GLFW OpenGL window
"""
function Window(
		;name = "GLWindow",
		resolution = standard_screen_resolution(),
		debugging = false,
		major = 3,
		minor = 3,# this is what GLVisualize needs to offer all features
		windowhints = standard_window_hints(),
		contexthints = standard_context_hints(major, minor),
		visible = true,
		focus = false,
		fullscreen = false,
		monitor = nothing,
		share = Window(C_NULL)
	)
	WindowHint(VISIBLE, visible)
	WindowHint(FOCUSED, focus)
	for ch in contexthints
		WindowHint(ch[1], ch[2])
	end
	for wh in windowhints
		WindowHint(wh[1], wh[2])
	end

	@static if Sys.isapple()
		if debugging
			@warn "OpenGL debug message callback is not available in macOS"
			debugging = false
		end
	end

	WindowHint(OPENGL_DEBUG_CONTEXT, Cint(debugging))

	monitor = if monitor == nothing
		GetPrimaryMonitor()
	elseif isa(monitor, Integer)
		GetMonitors()[monitor]
	elseif isa(monitor, Monitor)
		monitor
	else
		error("Monitor needs to be nothing, int, or GLFW.Monitor. Found: $monitor")
	end

	window = CreateWindow(resolution..., String(name), Monitor(C_NULL), share)

	if fullscreen
		SetKeyCallback(window, (_1, button, _2, _3, _4) -> begin
			button == KEY_ESCAPE && make_windowed!(window)
		end)
		make_fullscreen!(window, monitor)
	end
	MakeContextCurrent(window)
	window
end

function make_windowed!(window::Window)
	width, height = standard_screen_resolution()
	SetWindowMonitor(window, Monitor(C_NULL), 0, 0, width, height, DONT_CARE)
	return
end

function make_fullscreen!(window::Window, monitor::Monitor = GetPrimaryMonitor())
	vidmodes = GetVideoModes(monitor)[end]
	SetWindowMonitor(window, monitor, 0, 0, vidmodes.width, vidmodes.height, GLFW.DONT_CARE)
	return
end

"""
Sets visibility of OpenGL window. Will still render if not visible.
Only applies to the root screen holding the opengl context.
"""
function set_visibility!(screen::Window, visible::Bool)
	if visible
		ShowWindow(screen)
	else !visible
		HideWindow(screen)
	end
	return
end

struct Cursor
	handle::Ptr{Cvoid}
end

struct VidMode
	width::Cint         # The width, in screen coordinates, of the video mode.
	height::Cint        # The height, in screen coordinates, of the video mode.
	redbits::Cint       # The bit depth of the red channel of the video mode.
	greenbits::Cint     # The bit depth of the green channel of the video mode.
	bluebits::Cint      # The bit depth of the blue channel of the video mode.
	refreshrate::Cint   # The refresh rate, in Hz, of the video mode.
end
Base.show(io::IO, vm::VidMode) = write(io, "$(vm.width)x$(vm.height)@$(vm.refreshrate)Hz")

struct GLFWImage
	width::Cint
	height::Cint
	pixels::Ptr{UInt8}
end

function Base.cconvert(::Type{Ref{GLFWImage}}, image::AbstractMatrix{NTuple{4,UInt8}})
	imaget = permutedims(image) # c libglfw expects row major matrix
	return Base.RefValue{GLFWImage}(GLFWImage(size(imaget)..., pointer(imaget)))
end

function Base.cconvert(::Type{Ref{GLFWImage}}, images::Vector{<:AbstractMatrix{NTuple{4,UInt8}}})
	imagest = permutedims.(images)
	out = Vector{GLFWImage}(undef, length(imagest))
	@inbounds for i in 1:length(imagest)
		out[i] = GLFWImage(size(imagest[i])..., pointer(imagest[i]))
	end
	return (out, imagest)
end

function Base.unsafe_convert(::Type{Ref{GLFWImage}}, data::Tuple{Vector{GLFWImage},Vector{<:AbstractMatrix{NTuple{4,UInt8}}}})
	Base.unsafe_convert(Ref{GLFWImage}, data[1])
end

struct GLFWError <: Exception
	code::Union{ErrorCode, Integer}
	description::String

	GLFWError(code::ErrorCode, desc) = new(code, desc)
	GLFWError(code::Integer, desc) = new(try ErrorCode(code) catch; code end, desc)
end
Base.showerror(io::IO, e::GLFWError) = print(io, "GLFWError ($(e.code)): ", e.description)

#************************************************************************
# GLFW API functions
#************************************************************************

const INITIALIZED = Ref(false)

is_initialized() = INITIALIZED[]

# Initialization and version information
function Init()
	INITIALIZED[] = Bool(ccall((:glfwInit, libglfw), Cint, ())) || error("glfwInit failed")
end

function Terminate()
	INITIALIZED[] = false
	ccall((:glfwTerminate, libglfw), Cvoid, ())
end

GetVersionString() = unsafe_string(ccall((:glfwGetVersionString, libglfw), Cstring, ()))

# Error handling
@callback Error(code::Cint, description::Cstring) -> (GLFWError(code, unsafe_string(description)),)

# Monitor handling
function GetMonitors()
	count = Ref{Cint}()
	ptr = ccall((:glfwGetMonitors, libglfw), Ptr{Monitor}, (Ref{Cint},), count)
	unsafe_wrap(Array, ptr, count[])
end

GetPrimaryMonitor() = ccall((:glfwGetPrimaryMonitor, libglfw), Monitor, ())

function GetMonitorPos(monitor::Monitor)
	x, y = Ref{Cint}(), Ref{Cint}()
	ccall((:glfwGetMonitorPos, libglfw), Cvoid, (Monitor, Ref{Cint}, Ref{Cint}), monitor, x, y)
	(x=x[], y=y[])
end

function GetMonitorPhysicalSize(monitor::Monitor)
	width, height = Ref{Cint}(), Ref{Cint}()
	ccall((:glfwGetMonitorPhysicalSize, libglfw), Cvoid, (Monitor, Ref{Cint}, Ref{Cint}), monitor, width, height)
	(width=width[], height=height[])
end

GetMonitorName(monitor::Monitor) = unsafe_string(ccall((:glfwGetMonitorName, libglfw), Cstring, (Monitor,), monitor))
@callback Monitor(monitor::Monitor, event::DeviceConfigEvent)

function GetVideoModes(monitor::Monitor)
	count = Ref{Cint}()
	ptr = ccall((:glfwGetVideoModes, libglfw), Ptr{VidMode}, (Monitor, Ref{Cint}), monitor, count)
	unsafe_wrap(Array, ptr, count[])
end

GetVideoMode(monitor::Monitor) = unsafe_load(ccall((:glfwGetVideoMode, libglfw), Ptr{VidMode}, (Monitor,), monitor))
SetGamma(monitor::Monitor, gamma::Real) = ccall((:glfwSetGamma, libglfw), Cvoid, (Monitor, Cfloat), monitor, gamma)

# Window handling
DefaultWindowHints() = ccall((:glfwDefaultWindowHints, libglfw), Cvoid, ())
WindowHint(target::Integer, hint::Integer) = ccall((:glfwWindowHint, libglfw), Cvoid, (Cint, Cint), target, hint)

# Pair window handles with a callback function list
# to prevent them from being garbage-collected
const _window_callbacks = Dict{Window, Ref{Vector{Callback}}}()

function CreateWindow(width::Integer, height::Integer, title::AbstractString, monitor::Monitor=Monitor(C_NULL), share::Window=Window(C_NULL))
	window = ccall((:glfwCreateWindow, libglfw), Window, (Cint, Cint, Cstring, Monitor, Window), width, height, title, monitor, share)
	callbacks = Ref{Vector{Callback}}(fill(nothing, _window_callback_num[]))
	_window_callbacks[window] = callbacks
	ccall((:glfwSetWindowUserPointer, libglfw), Cvoid, (Window, Ref{Vector{Callback}}), window, callbacks)
	window
end

callbacks(window) = ccall((:glfwGetWindowUserPointer, libglfw), Ref{Vector{Callback}}, (Window,), window)

function DestroyWindow(window::Window)
	ccall((:glfwDestroyWindow, libglfw), Cvoid, (Window,), window)
	delete!(_window_callbacks, window)
end

WindowShouldClose(window::Window) = ccall((:glfwWindowShouldClose, libglfw), Cint, (Window,), window) != 0
SetWindowShouldClose(window::Window, value::Bool) = ccall((:glfwSetWindowShouldClose, libglfw), Cvoid, (Window, Cint), window, value)
SetWindowTitle(window::Window, title::AbstractString) = ccall((:glfwSetWindowTitle, libglfw), Cvoid, (Window, Cstring), window, title)

"""
    GLFW.SetWindowIcon(window::Window, image::Matrix{NTuple{4, UInt8}})

	GLFW.SetWindowIcon(window::Window, images::Vector{<:AbstractMatrix{NTuple{4,UInt8}}})

Set the window icon, where a single image may be passed or a vector of images with different icon sizes.
The images must be of RGBA format. Before calling this function it might be necessary to reinterpret the image
as a matrix of element type NTuple{4, UInt8}  if the icons are loaded with type RGBA{N0f8}

# Examples
```julia-repl
using FileIO
icons = load.(["icon-16.png", "icon-32.png", "icon-128.png"])
buffs = reinterpret.(NTuple{4, UInt8}, icons)
GLFW.SetWindowIcon(win, buffs)
GLFW.PollEvents() # needs a poll events to become active
```
"""
SetWindowIcon

function SetWindowIcon(window::Window, images::Vector{<:AbstractMatrix{NTuple{4,UInt8}}})
	ccall((:glfwSetWindowIcon, libglfw), Cvoid, (Window, Cint, Ref{GLFWImage}), window, length(images), images)
end

function SetWindowIcon(window::Window, image::AbstractMatrix{NTuple{4,UInt8}})
	ccall((:glfwSetWindowIcon, libglfw), Cvoid, (Window, Cint, Ref{GLFWImage}), window, 1, image)
end

function GetWindowPos(window::Window)
	x, y = Ref{Cint}(), Ref{Cint}()
	ccall((:glfwGetWindowPos, libglfw), Cvoid, (Window, Ref{Cint}, Ref{Cint}), window, x, y)
	(x=x[], y=y[])
end

SetWindowPos(window::Window, x::Integer, y::Integer) = ccall((:glfwSetWindowPos, libglfw), Cvoid, (Window, Cint, Cint), window, x, y)

function GetWindowSize(window::Window)
	width, height = Ref{Cint}(), Ref{Cint}()
	ccall((:glfwGetWindowSize, libglfw), Cvoid, (Window, Ref{Cint}, Ref{Cint}), window, width, height)
	(width=width[], height=height[])
end

SetWindowSizeLimits(window, minwidth, minheight, maxwidth, maxheight) = ccall((:glfwSetWindowSizeLimits, libglfw),
	Cvoid, (Window, Cint, Cint, Cint, Cint), window, minwidth, minheight, maxwidth, maxheight)
SetWindowAspectRatio(window, numer, denom) = ccall((:glfwSetWindowAspectRatio, libglfw),
	Cvoid, (Window, Cint, Cint), window, numer, denom)
SetWindowAspectRatio(window, r::Rational) = SetWindowAspectRatio(window, r.num, r.den)
SetWindowSize(window::Window, width::Integer, height::Integer) = ccall((:glfwSetWindowSize, libglfw), Cvoid, (Window, Cint, Cint), window, width, height)

function GetFramebufferSize(window::Window)
	width, height = Ref{Cint}(), Ref{Cint}()
	ccall((:glfwGetFramebufferSize, libglfw), Cvoid, (Window, Ref{Cint}, Ref{Cint}), window, width, height)
	(width=width[], height=height[])
end

function GetWindowFrameSize(window::Window)
	left, top, right, bottom = Ref{Cint}(), Ref{Cint}(), Ref{Cint}(), Ref{Cint}()
	ccall((:glfwGetWindowFrameSize, libglfw), Cvoid, (Window, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}), window, left, top, right, bottom)
	(left=left[], top=top[], right=right[], bottom=bottom[])
end

IconifyWindow(window::Window) = ccall((:glfwIconifyWindow, libglfw), Cvoid, (Window,), window)
RestoreWindow(window::Window) = ccall((:glfwRestoreWindow, libglfw), Cvoid, (Window,), window)
MaximizeWindow(window) = ccall((:glfwMaximizeWindow, libglfw), Cvoid, (Window,), window)
ShowWindow(window::Window) = ccall((:glfwShowWindow, libglfw), Cvoid, (Window,), window)
HideWindow(window::Window) = ccall((:glfwHideWindow, libglfw), Cvoid, (Window,), window)
GetWindowMonitor(window::Window) = ccall((:glfwGetWindowMonitor, libglfw), Monitor, (Window,), window)
# TODO: Add SetWindowMonitor variants:
# - Monitor with video mode
# - Nothing with size and position
SetWindowMonitor(window, monitor, xpos, ypos, width, height, refreshRate) = ccall((:glfwSetWindowMonitor, libglfw),
	Cvoid, (Window, Monitor, Cint, Cint, Cint, Cint, Cint), window, monitor, xpos, ypos, width, height, refreshRate)
GetWindowAttrib(window::Window, attrib::Integer) = ccall((:glfwGetWindowAttrib, libglfw), Cint, (Window, Cint), window, attrib)
@windowcallback WindowPos(window::Window, x::Cint, y::Cint)
@windowcallback WindowSize(window::Window, width::Cint, height::Cint)
@windowcallback WindowClose(window::Window)
@windowcallback WindowRefresh(window::Window)
@windowcallback WindowFocus(window::Window, focused::Cint) -> (window, Bool(focused))
@windowcallback WindowIconify(window::Window, iconified::Cint) -> (window, Bool(iconified))
@windowcallback FramebufferSize(window::Window, width::Cint, height::Cint)
PollEvents() = ccall((:glfwPollEvents, libglfw), Cvoid, ())
WaitEvents() = ccall((:glfwWaitEvents, libglfw), Cvoid, ())
WaitEvents(timeout) = ccall((:glfwWaitEventsTimeout, libglfw), Cvoid, (Cdouble,), timeout)
PostEmptyEvent() = ccall((:glfwPostEmptyEvent, libglfw), Cvoid, ())

# Input handling
function GetInputMode(window::Window, mode::Integer)
	value = ccall((:glfwGetInputMode, libglfw), Cint, (Window, Cint), window, mode)
	if mode in (STICKY_KEYS, STICKY_MOUSE_BUTTONS)
		value = Bool(value)
	end
	value
end

SetInputMode(window::Window, mode::Integer, value::Integer) = ccall((:glfwSetInputMode, libglfw), Cvoid, (Window, Cint, Cint), window, mode, value)

function GetKeyName(key, scancode=0)
	ptr = ccall((:glfwGetKeyName, libglfw), Cstring, (Cint, Cint), key, scancode)
	if ptr != C_NULL
		unsafe_string(ptr)
	end
end

GetKey(window::Window, key) = Bool(ccall((:glfwGetKey, libglfw), Cint, (Window, Cint), window, key))
GetMouseButton(window::Window, button::MouseButton) = Bool(ccall((:glfwGetMouseButton, libglfw), Cint, (Window, Cint), window, button))

function GetCursorPos(window::Window)
	x, y = Ref{Cdouble}(), Ref{Cdouble}()
	ccall((:glfwGetCursorPos, libglfw), Cvoid, (Window, Ref{Cdouble}, Ref{Cdouble}), window, x, y)
	(x=x[], y=y[])
end

SetCursorPos(window::Window, x::Real, y::Real) = ccall((:glfwSetCursorPos, libglfw), Cvoid, (Window, Cdouble, Cdouble), window, x, y)
# TODO: Add glfwCreateCursor
CreateStandardCursor(shape::StandardCursorShape) = ccall((:glfwCreateStandardCursor, libglfw), Cursor, (Cint,), shape)
DestroyCursor(cursor::Cursor) = ccall((:glfwDestroyCursor, libglfw), Cvoid, (Cursor,), cursor)
SetCursor(window::Window, cursor::Cursor) = ccall((:glfwSetCursor, libglfw), Cvoid, (Window, Cursor), window, cursor)
SetCursor(window::Window, ::Nothing) = SetCursor(window, Cursor(C_NULL))
@windowcallback Key(window::Window, key::Key, scancode::Cint, action::Action, mods::Cint)
@windowcallback Char(window::Window, codepoint::Cuint) -> (window, convert(Char, codepoint))
@windowcallback CharMods(window::Window, codepoint::Cuint, mods::Cint) -> (window, convert(Char, codepoint), mods)
@windowcallback MouseButton(window::Window, button::MouseButton, action::Action, mods::Cint)
@windowcallback CursorPos(window::Window, x::Cdouble, y::Cdouble)
@windowcallback CursorEnter(window::Window, entered::Cint) -> (window, Bool(entered))
@windowcallback Scroll(window::Window, xoffset::Cdouble, yoffset::Cdouble)
@windowcallback Drop(window::Window, count::Cint, paths::Ptr{Cstring}) -> (window, map(unsafe_string, unsafe_wrap(Array, paths, count)))
JoystickPresent(joy::Joystick) = Bool(ccall((:glfwJoystickPresent, libglfw), Cint, (Cint,), joy))

function GetJoystickAxes(joy::Joystick)
	count = Ref{Cint}()
	ptr = ccall((:glfwGetJoystickAxes, libglfw), Ptr{Cfloat}, (Cint, Ref{Cint}), joy, count)
	if ptr != C_NULL
		unsafe_wrap(Array, ptr, count[])
	end
end

function GetJoystickButtons(joy::Joystick)
	count = Ref{Cint}()
	ptr = ccall((:glfwGetJoystickButtons, libglfw), Ptr{Cuchar}, (Cint, Ref{Cint}), joy, count)
	if ptr != C_NULL
		unsafe_wrap(Array, ptr, count[])
	end
end

function GetJoystickName(joy::Joystick)
	ptr = ccall((:glfwGetJoystickName, libglfw), Cstring, (Cint,), joy)
	if ptr != C_NULL
		unsafe_string(ptr)
	end
end

@callback Joystick(joy::Joystick, event::DeviceConfigEvent)

# Context handling
MakeContextCurrent(window::Window) = ccall((:glfwMakeContextCurrent, libglfw), Cvoid, (Window,), window)
GetCurrentContext() = Base.cconvert(Window, ccall((:glfwGetCurrentContext, libglfw), Window, ()))
SwapBuffers(window::Window) = ccall((:glfwSwapBuffers, libglfw), Cvoid, (Window,), window)
SwapInterval(interval::Integer) = ccall((:glfwSwapInterval, libglfw), Cvoid, (Cint,), interval)
ExtensionSupported(extension::AbstractString) = Bool(ccall((:glfwExtensionSupported, libglfw), Cint, (Cstring,), extension))
GetProcAddress(procname::AbstractString) = ccall((:glfwGetProcAddress, libglfw), Ptr{Cvoid}, (Cstring,), procname)

#came from GLWindow/core.jl
"""
Returns the monitor resolution of the primary monitor.
"""
function primarymonitorresolution()
	props = MonitorProperties(GetPrimaryMonitor())
	w,h = props.videomode.width, props.videomode.height
	# Vec(Int(w),Int(h))
	(Int(w),Int(h))
end

#Came from GLWindow.jl/screen.jl

"""
Takes half the resolution of the primary monitor.
This should make for sensible defaults!
"""
function standard_screen_resolution()
	w, h = primarymonitorresolution()
	(div(w,2), div(h,2)) # half of total resolution seems like a good fit!
end

"""
Tries to create sensible context hints!
Taken from lessons learned at:
[GLFW](http://www.glfw.org/docs/latest/window.html)
"""
function standard_context_hints(major, minor)
	# this is spaar...Modern OpenGL !!!!
	major < 3 && error("OpenGL major needs to be at least 3.0. Given: $major")
	# core profile is only supported for OpenGL 3.2+ (and a must for OSX, so
	# for the sake of homogenity, we try to default to it for everyone!)
	if (major > 3 || (major == 3 && minor >= 2 ))
		profile = OPENGL_CORE_PROFILE
	else
		profile = OPENGL_ANY_PROFILE
	end
	[
		(CONTEXT_VERSION_MAJOR, major),
		(CONTEXT_VERSION_MINOR, minor),
		(OPENGL_FORWARD_COMPAT, Cint(1)),
		(OPENGL_PROFILE, profile)
	]
end

"""
Standard window hints for creating a plain context without any multisampling
or extra buffers beside the color buffer
"""
function standard_window_hints()
	[
		(SAMPLES,      0),
		(DEPTH_BITS,   0),

		(ALPHA_BITS,   8),
		(RED_BITS,     8),
		(GREEN_BITS,   8),
		(BLUE_BITS,    8),

		(STENCIL_BITS, 0),
		(AUX_BUFFERS,  0)
	]
end


text:
@doc raw"""
Nonlocal term coming from norm-conserving pseudopotentials in Kleinmann-Bylander form.
``\text{Energy} = \sum_a \sum_{ij} \sum_{n} f_n <ψ_n|p_{ai}> D_{ij} <p_{aj}|ψ_n>.``
"""
struct AtomicNonlocal end
function (::AtomicNonlocal)(basis::PlaneWaveBasis{T}) where {T}
    model = basis.model

    # keep only pseudopotential atoms and positions
    psp_groups = [group for group in model.atom_groups
                  if model.atoms[first(group)] isa ElementPsp]
    psps          = [model.atoms[first(group)].psp      for group in psp_groups]
    psp_positions = [model.positions[group] for group in psp_groups]

    isempty(psp_groups) && return TermNoop()
    ops = map(basis.kpoints) do kpt
        P = build_projection_vectors_(basis, kpt, psps, psp_positions)
        D = build_projection_coefficients_(T, psps, psp_positions)
        NonlocalOperator(basis, kpt, P, D)
    end
    TermAtomicNonlocal(ops)
end

struct TermAtomicNonlocal <: Term
    ops::Vector{NonlocalOperator}
end

@timing "ene_ops: nonlocal" function ene_ops(term::TermAtomicNonlocal,
                                             basis::PlaneWaveBasis{T},
                                             ψ, occ; kwargs...) where {T}
    isnothing(ψ) && return (E=T(Inf), ops=term.ops)

    E = zero(T)
    for (ik, kpt) in enumerate(basis.kpoints)
        Pψ = term.ops[ik].P' * ψ[ik]  # nproj x nband
        band_enes = dropdims(sum(real.(conj.(Pψ) .* (term.ops[ik].D * Pψ)), dims=1), dims=1)
        E += basis.kweights[ik] * sum(band_enes .* occ[ik])
    end
    E = mpi_sum(E, basis.comm_kpts)

    (E=E, ops=term.ops)
end

@timing "forces: nonlocal" function compute_forces(::TermAtomicNonlocal,
                                                   basis::PlaneWaveBasis{TT},
                                                   ψ, occ; kwargs...) where TT
    T = promote_type(TT, real(eltype(ψ[1])))
    model = basis.model
    unit_cell_volume = model.unit_cell_volume
    psp_groups = [group for group in model.atom_groups
                  if model.atoms[first(group)] isa ElementPsp]

    # early return if no pseudopotential atoms
    isempty(psp_groups) && return nothing

    # energy terms are of the form <psi, P C P' psi>, where P(G) = form_factor(G) * structure_factor(G)
    forces = [zero(Vec3{T}) for _ in 1:length(model.positions)]
    for group in psp_groups
        element = model.atoms[first(group)]

        C = build_projection_coefficients_(T, element.psp)
        for (ik, kpt) in enumerate(basis.kpoints)
            # we compute the forces from the irreductible BZ; they are symmetrized later
            qs_cart = Gplusk_vectors_cart(basis, kpt)
            qs = Gplusk_vectors(basis, kpt)
            form_factors = build_form_factors(element.psp, qs_cart)
            for idx in group
                r = model.positions[idx]
                structure_factors = [cis2pi(-dot(q, r)) for q in qs]
                P = structure_factors .* form_factors ./ sqrt(unit_cell_volume)

                forces[idx] += map(1:3) do α
                    dPdR = [-2T(π)*im*q[α] for q in qs] .* P
                    ψk = ψ[ik]
                    dHψk = P * (C * (dPdR' * ψk))
                    -sum(occ[ik][iband] * basis.kweights[ik] *
                         2real(dot(ψk[:, iband], dHψk[:, iband]))
                         for iband=1:size(ψk, 2))
                end  # α
            end  # r
        end  # kpt
    end  # group

    forces = mpi_sum!(forces, basis.comm_kpts)
    symmetrize_forces(basis, forces)
end

# TODO possibly move over to pseudo/NormConservingPsp.jl ?
# Build projection coefficients for a atoms array generated by term_nonlocal
# The ordering of the projector indices is (A,l,m,i), where A is running over all
# atoms, l, m are AM quantum numbers and i is running over all projectors for a
# given l. The matrix is block-diagonal with non-zeros only if A, l and m agree.
function build_projection_coefficients_(T, psps, psp_positions)
    # TODO In the current version the proj_coeffs still has a lot of zeros.
    #      One could improve this by storing the blocks as a list or in a
    #      BlockDiagonal data structure
    n_proj = count_n_proj(psps, psp_positions)
    proj_coeffs = zeros(T, n_proj, n_proj)

    count = 0
    for (psp, positions) in zip(psps, psp_positions), _ in positions
        n_proj_psp = count_n_proj(psp)
        block = count+1:count+n_proj_psp
        proj_coeffs[block, block] = build_projection_coefficients_(T, psp)
        count += n_proj_psp
    end # psp, r
    @assert count == n_proj

    proj_coeffs
end

# Builds the projection coefficient matrix for a single atom
# The ordering of the projector indices is (l,m,i), where l, m are the
# AM quantum numbers and i is running over all projectors for a given l.
# The matrix is block-diagonal with non-zeros only if l and m agree.
function build_projection_coefficients_(T, psp::NormConservingPsp)
    n_proj = count_n_proj(psp)
    proj_coeffs = zeros(T, n_proj, n_proj)
    count = 0
    for l in 0:psp.lmax, m in -l:l
        n_proj_l = size(psp.h[l + 1], 1)  # Number of i's
        range = count .+ (1:n_proj_l)
        proj_coeffs[range, range] = psp.h[l + 1]
        count += n_proj_l
    end # l, m
    proj_coeffs
end

"""
Build projection vectors for a atoms array generated by term_nonlocal

H_at  = sum_ij Cij |pi> <pj|
H_per = sum_R sum_ij Cij |pi(x-R)> <pj(x-R)|
      = sum_R sum_ij Cij |pi(x-R)> <pj(x-R)|

<e_kG'|H_per|e_kG> = ...
                   = 1/Ω sum_ij Cij pihat(k+G') pjhat(k+G)^*

where pihat(q) = ∫_R^3 pi(r) e^{-iqr} dr

We store 1/√Ω pihat(k+G) in proj_vectors.
"""
function build_projection_vectors_(basis::PlaneWaveBasis{T}, kpt::Kpoint,
                                   psps, psp_positions) where {T}
    unit_cell_volume = basis.model.unit_cell_volume
    n_proj = count_n_proj(psps, psp_positions)
    n_G    = length(G_vectors(basis, kpt))
    proj_vectors = zeros(Complex{T}, n_G, n_proj)

    # Compute the columns of proj_vectors = 1/√Ω pihat(k+G)
    # Since the pi are translates of each others, pihat(k+G) decouples as
    # pihat(q) = ∫ p(r-R) e^{-iqr} dr = e^{-iqR} phat(q).
    # The first term is the structure factor, the second the form factor.
    offset = 0  # offset into proj_vectors
    for (psp, positions) in zip(psps, psp_positions)
        # Compute position-independent form factors
        form_factors = build_form_factors(psp, Gplusk_vectors_cart(basis, kpt))

        # Combine with structure factors
        for r in positions
            # k+G in this formula can also be G, this only changes an unimportant phase factor
            structure_factors = map(q -> cis2pi(-dot(q, r)), Gplusk_vectors(basis, kpt))
            @views for iproj = 1:count_n_proj(psp)
                proj_vectors[:, offset+iproj] .= (
                    structure_factors .* form_factors[:, iproj] ./ sqrt(unit_cell_volume)
                )
            end
            offset += count_n_proj(psp)
        end
    end
    @assert offset == n_proj
    proj_vectors
end

"""
Build form factors (Fourier transforms of projectors) for an atom centered at 0.
"""
function build_form_factors(psp, qs)
    qnorms = norm.(qs)
    T = real(eltype(qnorms))
    # Compute position-independent form factors
    form_factors = zeros(Complex{T}, length(qs), count_n_proj(psp))
    count = 1
    for l in 0:psp.lmax, m in -l:l
        prefac_lm = im^l .* ylm_real.(l, m, qs)
        n_proj_l = size(psp.h[l + 1], 1)

        for iproj in 1:n_proj_l
            radial_il = eval_psp_projector_fourier.(psp, iproj, l, qnorms)
            form_factors[:, count] = prefac_lm .* radial_il
            count += 1
        end
    end
    @assert count == count_n_proj(psp) + 1
    form_factors
end


text:
using Taro, DataFrames

file = joinpath(Pkg.dir("GlassCatalog"),"assets","schott","schott_optical_glass_catalogue_excel_may_2016.xls");

firstrow = 4
lastrow = 20
firstcol = "A"
lastcol = "FJ"
range = "$firstcol$firstrow:$lastcol$lastrow";

table = Taro.readxl(file,"Datatable",range)

typeof(table)

#describe(table)


text:
using Test
using AcceleratedDubins
using Dubins

for run in 1:1
    tol = 1e-6    # for isapprox comparison
    params = [5, 15, 2.6, 4]

    @testset "original vs. accelerated" begin
        # compare lengths of the original path (pkg Dubins) and
        # accelerated with minimum turning radius (pkg AcceleratedDubins)
        start = [0, 0, 0.5]
        stop = [3, 0, pi]
        r_min = 1.

        # same positions, but shifted and rotated (check is rotation works correctly)
        start_r = [5, 5, 0.5+pi/2]
        stop_r = [5, 8, pi+pi/2]

        _, path_d = dubins_shortest_path(start, stop, r_min)
        path_a, _ = retrieve_path(start, stop, [r_min, r_min], params) # shortest path

        _, path_d_r = dubins_shortest_path(start_r, stop_r, r_min)
        path_a_r, _ = retrieve_path(start_r, stop_r, [r_min, r_min], params) # shortest path

        @test isapprox(dubins_path_length(path_d), path_len(path_a), atol=tol)
        @test isapprox(dubins_path_length(path_d_r), path_len(path_a), atol=tol)

    end

    @testset "straight path" begin
        start = [0, 0, 0.]
        stop = [6, 0, 0.]

        _, path_d = dubins_shortest_path(start, stop, 1.)
        len_d = dubins_path_length(path_d)

        path_s, _ = retrieve_path(start, stop, [1., 1.], params) # shortest path
        path_f, _ = fastest_path(start, stop, [1., 2., 3.], params)
        path_o, _ = optimized_path(start, stop, [1., 3.], params)

        len_s = path_len(path_s)
        len_f = path_len(path_f)
        len_o = path_len(path_o)

        @test isapprox(len_d, path_len(path_s), atol=tol)
        @test isapprox(len_s, len_f, atol=tol)
        @test isapprox(len_s, len_o, atol=tol)
    end

    @testset "proposed paths better or equal" begin
        for testrun in 1:20
            start = [10*rand(), 10*rand(), 2*pi*rand()]
            stop = [10*rand()+10, 10*rand()+5, 2*pi*rand()]

            _, path_d = dubins_shortest_path(start, stop, 1.)
            len_d = dubins_path_length(path_d)

            path_s, _ = retrieve_path(start, stop, [1., 1.], params) # shortest path
            path_f, _ = fastest_path(start, stop, [1., 2., 3.], params)

            time_s = path_time(path_s, params)
            time_f = path_time(path_f, params)

            @test time_s >= time_f
        end
    end
end


text:
using Plots
using Distributions
pyplot()

function gen_rand_uniform_mult(min1, max1, min2, max2, n)
    rand(Uniform(min1, max1), n) .* rand(Uniform(min2, max2), n)
end

histogram(gen_rand_uniform_mult(0, 10, 35, 40, 10000))

max_stars_repo_name: JakobWyatt/code-playground

max_issues_repo_name: JakobWyatt/code-playground

max_forks_repo_name: JakobWyatt/code-playground


text:
using Distributed
# add processes to workspace
addprocs(2)

@everywhere using ModelingToolkit, OrdinaryDiffEq

# create the Lorenz system
@everywhere @parameters t σ ρ β
@everywhere @variables x(t) y(t) z(t)
@everywhere D = Differential(t)

@everywhere eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

@everywhere @named de = ODESystem(eqs)
@everywhere ode_func = ODEFunction(de, [x,y,z], [σ, ρ, β])

@everywhere u0 = [19.,20.,50.]
@everywhere params = [16.,45.92,4]

@everywhere ode_prob = ODEProblem(ode_func, u0, (0., 10.),params)

@everywhere begin

    using OrdinaryDiffEq
    using ModelingToolkit

    function solve_lorenz(ode_problem)
        print(solve(ode_problem,Tsit5()))
    end
end

solve_lorenz(ode_prob)

future = @spawn solve_lorenz(ode_prob)
@test_broken fetch(future)

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
import Base.show
import Base64.stringmime

# Methods for mglGraph
function show(io::IO, m::MIME"text/html", gr::mglGraph)
    htm = """\n"
    show(io, m, HTML(htm))
end

function show(io::IO, m::MIME"text/latex", gr::mglGraph)
    name = "MathGL-"*randstring()*".png"
    WriteFrame(gr, name)
    write(io, "\\includegraphics[width=\\textwidth]{$name}")
end

function show(io::IO, m::MIME"application/x-latex", gr::mglGraph)
    name = "MathGL-"*randstring()*".png"
    WriteFrame(gr, name)
    show(io, m, "\\includegraphics[width=\\textwidth]{$name}")
end

function show(io::IO, m::MIME"image/png", gr::mglGraph)
    name = tempname()*".png"
    WriteFrame(gr, name)
    file = open(name, "r")
    write(io, read(file))
    close(file)
    rm(name)
end

# Methods for a plotOpStack
function show(io::IO, m::MIME"text/html", ops::MathGL.plotOpStack)
    htm = """\n"
    show(io, m, HTML(htm))
end

function show(io::IO, m::MIME"text/latex", ops::MathGL.plotOpStack)
    name = "MathGL-"*randstring()*".png"
    WriteFrame(draw(ops), name)
    write(io, "\\includegraphics[width=\\textwidth]{$name}")
end

function show(io::IO, m::MIME"application/x-latex", ops::MathGL.plotOpStack)
    name = "MathGL-"*randstring()*".png"
    WriteFrame(draw(ops), name)
    show(io, m, "\\includegraphics[width=\\textwidth]{$name}")
end

function show(io::IO, m::MIME"image/png", ops::MathGL.plotOpStack)
    name = tempname()*".png"
    WriteFrame(draw(ops), name)
    file = open(name, "r")
    write(io, read(file))
    close(file)
    rm(name)
end


text:
module TestMLBaseWrapper

include(joinpath("..", "fixture_learners.jl"))
using .FixtureLearners
fcp = FeatureClassification()

using FactCheck

importall Orchestra.Transformers.MLBaseWrapper

facts("MLBase transformers") do
  context("StandardScaler transforms features") do
    instances = [
      5 10;
      -5 0;
      0 5;
    ]
    labels = [
      "x";
      "y";
      "z";
    ]
    expected_transformed = [
      1.0 1.0;
      -1.0 -1.0;
      0.0 0.0;
    ]
    standard_scaler = StandardScaler()
    fit!(standard_scaler, instances, labels)
    transformed = transform!(standard_scaler, instances)

    @fact transformed => expected_transformed
  end
end

end # module


text:
using DiscreteFiltering
using LinearAlgebra
using Plots
using SparseArrays

## Domain
a = 0.0
b = 2π
# domain = PeriodicIntervalDomain(a, b)
domain = ClosedIntervalDomain(a, b)

## Discretization
N = 200
M = N
x = discretize(domain, M)
ξ = discretize(domain, N)
Δx = (b - a) / N

## Filter
# h₀ = Δx
# h(x) = h₀ # * (1 - 1 / 2 * cos(x))
# dh(x) = 0.0 # h₀ / 2 * sin(x)
# filter = TopHatFilter(h)

h₀ = 5.1Δx
h(x) = h₀ # * (1 - 1 / 2 * cos(x))
dh(x) = 0.0 # h₀ / 2 * sin(x)
filter = GaussianFilter(h, Δx / 2)
# filter = ConvolutionalFilter(h, x -> (-Δx / 2 ≤ x ≤ Δx / 2) / Δx)

# Filter matrix
W = filter_matrix(filter, domain, M, N)

## Equations
f = (ξ, t) -> 0.0
g_a = t -> 0.0
g_b = t -> 0.0
equation = DiffusionEquation(domain, IdentityFilter(), f, g_a, g_b)
equation_filtered = DiffusionEquation(domain, filter, f, g_a, g_b)

## Time
T = 0.04
t = T

# ODE solver tolerances
tols = (; abstol = 1e-6, reltol = 1e-4)
degmax = 50
λ = 1e-4

## Plot filter
plot(ξ, h)
ylims!((0, ylims()[2]))

## Exact solutions (force zero BC)
_u₀(x) = sin(x) + 0.6cos(5x) + 0.04sin(20(x - 1))
u₀(x) = _u₀(x) - (_u₀(a) + (_u₀(b) - _u₀(a)) * (x - a) / (b - a))
u₀_int(x) =
    -cos(x) + 0.6 / 5 * sin(5x) - 0.04 / 20 * cos(20(x - 1)) -
    (_u₀(a)x + (_u₀(b) - _u₀(a)) * (x - a)^2 / 2(b - a))
if domain isa PeriodicIntervalDomain
    ū₀(x) = 1 / 2h(x) * (u₀_int(x + h(x)) - u₀_int(x - h(x)))
else
    ū₀(x) = begin
        α = max(a, x - h(x))
        β = min(b, x + h(x))
        1 / (β - α) * (u₀_int(β) - u₀_int(α))
    end
end

## Discrete initial conditions
uₕ = u₀.(ξ)
ūₕ = ū₀.(x)

plot(ξ, uₕ, label = "Unfiltered")
plot!(x, ūₕ, label = "Filtered")
title!("Initial conditions")

## Solve discretized problem
sol = solve(
    equation,
    u₀,
    (0.0, T),
    M,
    N;
    method = "discretizefirst",
    boundary_conditions = "derivative",
    tols...,
    degmax,
    λ,
)
plot(x, uₕ, label = "Initial conditions")
plot!(x, sol(t), label = "Discretized")
title!("Solution")

## Solve filtered-and-then-discretized problem with ADBC
ū_adbc = solve_adbc(equation_filtered, u₀, (0.0, T), M, T / 100_000)
plot(x, uₕ, label = "Initial")
# plot!(x, ūₕ, label = "Initial filtered")
plot!(x, ū_adbc, label = "Filtered-then-discretized (ADBC)")
title!("Solution")

## Solve discretized-and-then-filtered problem
sol_allbar = solve(
    equation_filtered,
    u₀,
    (0.0, T),
    M,
    N;
    method = "discretizefirst",
    boundary_conditions = "derivative",
    tols...,
    degmax,
    λ,
)

plot(ξ, uₕ, label = "Initial")
# plot!(x, ūₕ, label = "Initial discretized-then-filtered")
plot!(x, sol_allbar(t), label = "Discretized-then-filtered")
title!("Solution")

## Comparison
plot(ξ, uₕ, label = "Initial")
# plot!(x, ūₕ, label = "Initial filtered")
plot!(ξ, sol(t), label = "Discretized")
# plot!(x, sol_bar(t), label = "Filtered-then-discretized")
plot!(x, ū_adbc, label = "Filtered-then-discretized")
plot!(x, sol_allbar(t), label = "Discretized-then-filtered")
# plot!(ξ, [u.(ξ, t), ū.(ξ, t)], label = "Exact")
ylims!(minimum(uₕ), maximum(uₕ))
title!("Solution")

## Relative error
u_exact = sol(t)
ū_exact = W * sol(t)
# err = abs.(sol(t) - u_exact) ./ maximum(abs.(u_exact))
err_bar = abs.(ū_adbc - ū_exact) ./ maximum(abs.(ū_exact))
err_allbar = abs.(sol_allbar(t) - ū_exact) ./ maximum(abs.(ū_exact))

##
plot(yaxis = :log)
# plot!(ξ, err, label = "Unfiltered discretized")
plot!(x, err_bar, label = "Filtered-then-discretized (ADBC)")
plot!(x, err_allbar, label = "Discretized-then-filtered")
title!("Relative error")


text:
"""
A Julia package to iteratively solve linear problems
with custom matrix-vector operations.

This package contains some codes from IterativeSolvers.jl: Copyright (c) 2013--2016 The Julia Language.
For the copyright notice, please check the LICENSE file.
"""
module LinearSolvers

export bicgstab!, qmrcgstab!

using LinearAlgebra
using Random
using Printf

import Base: iterate

abstract type VectorLike end
abstract type MatrixLike end
abstract type ScalarLike end

mul!(x::VectorLike, A::MatrixLike, y::VectorLike) = begin x .= A * y end
mul!(x::VectorLike, y::ScalarLike, z::VectorLike) = begin x .= y .* z end
sub!(x::VectorLike, y::VectorLike, z::VectorLike) = begin x .= y .- z end
add!(x::VectorLike, y::VectorLike, z::VectorLike) = begin x .= y .+ z end
dot(x::Nothing, y::Nothing) = nothing # LinearAlgebra.dot(x,y)
norm(x::Nothing) = nothing # LinearAlgebra.norm(x)
add_scmul!(x::VectorLike, scalar::ScalarLike, y::VectorLike) = begin x .+= scalar .* y end
sub_scmul!(x::VectorLike, scalar::ScalarLike, y::VectorLike) = begin x .-= scalar .* y end
eltype(x::Nothing) = nothing # Base.eltype(x)
size(x::Nothing) = nothing # Base.size(x)
fill!(x::Nothing, scalar::ScalarLike) = nothing # Base.fill!(x, scalar)
rand!(x::Nothing) = nothing # Random.rand!(x)
copyto!(x::Nothing, y::Nothing) = nothing # Base.copyto!(x, y)

include("bicgstab.jl")
include("qmrcgstab.jl")

end


text:
"""
  isInside(polygon, n, p)
Returns true if the point p lies inside the polygon with n vertices
"""
function is_inside(polygon::Vector{Point2D{T}}, n::Int, p::Point2D{T}; extreme = Point2D(100000.0, p[2])) where T<:Real
  # There must be at least 3 vertices in polygon
  if (n < 3) return false end

  # Count intersections of the above line with sides of polygon
  count = 0
  for i=1:n
    next = mod(i,n)+1
    # Check if the line segment from 'p' to 'extreme' intersects
    # with the line segment from 'polygon[i]' to 'polygon[next]'
    if (do_intersect(polygon[i], polygon[next], p, extreme))
      # If the point 'p' is colinear with line segment 'i-next',
      # then check if it lies on segment. If it lies, return true,
      # otherwise false
      if (get_orientation(polygon[i], p, polygon[next]) == 0)
         return on_segment(polygon[i], p, polygon[next])
      end
      count += 1
    end
  end

  # Return true if count is odd, false otherwise
  if isodd(count)
    return true
  else
    return false
  end
end

"""
do_intersect(p1,q1,p2,q2) -> bool

Checks whether two lines, defined by (p1,q1) and (p2,q2) where pi,qi are 2-arrays.

Algorithm explained in https://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/.

Note we have NOT implemented the special case for overlapping line segments, as that will never be the case in finite element applications
"""
function do_intersect(p1::Point2D{T},q1::Point2D{T},p2::Point2D{T},q2::Point2D{T}) where T<:Real
  # Find the four orientations needed for general case
  o1 = get_orientation(p1, q1, p2)
  o2 = get_orientation(p1, q1, q2)
  o3 = get_orientation(p2, q2, p1)
  o4 = get_orientation(p2, q2, q1)

  # General case
  if (o1 != o2 && o3 != o4)
    return true
  else
    return false
  end
end

"""
  orientation(p,q,r) -> Int

Computes orientation of 3 (ordered) points p,q,r. where each is a 2-array
2 -> counterclockwise (CCW)
1 -> clockwise        (CW)
0 -> collinear

Algorithm is explained in https://www.geeksforgeeks.org/orientation-3-ordered-points/
"""
function get_orientation(p::Point2D{T},q::Point2D{T},r::Point2D{T}) where T<:Real
  val = (q.y-p.y)*(r.x-q.x) - (q.x-p.x)*(r.y-q.y);

  if (val == 0) return 0 end  # colinear

  return (val > 0) ? 1 : 2 # clock- or counterclock-wise
end

"""
  onSegment(p,q,r)

Checks whether the point `q` is on the line segment from `p` to `r`.
Returns a bool.
"""
function on_segment(p::Point2D{T},q::Point2D{T},r::Point2D{T}) where T<:Real
  if -eps() < (distance(p, q) + distance(q, r) - distance(p, r)) < eps()
    return true
  else
    return false
  end
end

function distance(a::Point2D{T},b::Point2D{T}) where T<:Real
  return sqrt((a.x - b.x)^2 + (a.y - b.y)^2)
end

#    NearestPoint!(point,node,wall)
#
#Compute nearest point to particle on LineWall, CircleWall, and ArcWall
function get_nearest_point!(point::Point2D{T},node::Particle2D{T},wall::LineWall) where T<:Real
  px=node.pos.x; py=node.pos.y

  Ax=wall.nodes[1].x; Ay=wall.nodes[1].y
  Bx=wall.nodes[2].x; By=wall.nodes[2].y

  bx=px-Ax; by=py-Ay
  ax=Bx-Ax; ay=By-Ay

  ℓ2 = ax^2+ay^2

  dotprod = ax*bx + ay*by

  point.x = dotprod*ax/ℓ2 + Ax
  point.y = dotprod*ay/ℓ2 + Ay

  nothing
end
function get_nearest_point!(point::Point2D{T},node::Particle2D{T},wall::CircleWall) where T<:Real
  px=node.pos.x; py=node.pos.y
  cx=wall.center.x; cy=wall.center.y
  r = wall.radius
  θ = atan(py-cy,px-cx)

  point.x = cx + r*cos(θ)
  point.y = cy + r*sin(θ)

  nothing
end
function get_nearest_point!(point::Point2D{T},node::Particle2D{T},wall::ArcWall) where T<:Real
  px= node.pos.x; py=node.pos.y
  cx= wall.nodes[2].x; cy=wall.nodes[2].y
  r = sqrt((cx-wall.nodes[1].x)^2 + (cy-wall.nodes[1].y)^2)
  θ = atan(py-cy,px-cx)

  point.x = cx + r*cos(θ)
  point.y = cy + r*sin(θ)

  nothing
end
function get_nearest_point(node::Particle2D{T},wall::W) where {T<:Real,W<:AbstractWall}
  point = Point2D(0.0,0.0)
  get_nearest_point!(point,node,wall)
  return point
end

# function to determine whether quadrature node (sx,sy) is within line
function is_in_line(wall::LineWall{T},point::Point2D{T}) where T<:Real
  return on_segment(wall.nodes[1],point,wall.nodes[2])
end
function is_in_line(wall::CircleWall{T},point::Point2D{T}) where T<:Real
  cx = wall.center.x; cy=wall.center.y
  val = abs(wall.radius - sqrt((cx-point.x)^2 + (cy-point.y)^2))
  TOL = 1e-12
  return (val < TOL ? true : false)
end
function is_in_line(wall::ArcWall{T},point::Point2D{T}) where T<:Real
  cx = wall.nodes[2].x; cy=wall.nodes[2].y
  radius = sqrt((cx-wall.nodes[1].x)^2 + (cy-wall.nodes[1].y)^2)
  val = abs(radius - sqrt((cx-point.x)^2 + (cy-point.y)^2))

  θ1 = atan(wall.nodes[1].y-cy,wall.nodes[1].x-cx) - pi/24
  θ2 = atan(wall.nodes[3].y-cy,wall.nodes[3].x-cx) + pi/24
  θs = atan(point.y-cy,point.x-cx)

  TOL = 1e-12
  isOn = false
  if val < TOL
      if θ1 > θ2
          if θs > θ1 || θs < θ2
              isOn = true
          end
      else
          if θs > θ1 && θs < θ2
              return true
          end
      end
  end

  return isOn
end

# determines whether or not point is inside a rectangle
# also includes being on bottom or left boundary
function is_inside_rect(rect::Vector{T},point::P)  where {T<:Real,P<:AbstractPoint}
  x0=rect[1]; x1=rect[2]; y0=rect[3]; y1=rect[4]

  if point.x<=x1 && point.x>=x0 && point.y<=y1 && point.y>=y0
    return true
  else
    return false
  end
end


text:
ENV["JULIA_ASSETREGISTRY_BASEURL"] = "/mybasedir"

using AssetRegistry
using Distributed
ps = addprocs(1)

using Test

using JSON

@testset "register" begin
    key = AssetRegistry.register(pwd())
    @test startswith(key, ENV["JULIA_ASSETREGISTRY_BASEURL"])
    @test key == AssetRegistry.getkey(pwd())
    @test AssetRegistry.isregistered(pwd())
    @test JSON.parse(String(read(joinpath(homedir(), ".jlassetregistry.json")))) == Dict(AssetRegistry.filekey(key) => [pwd(), 1])

    # test that a new proc doesn't leave any residue
    run(```
        $(Base.julia_cmd())
        -e 'using AssetRegistry; AssetRegistry.register(pwd());'
        ```)

    @test JSON.parse(String(read(joinpath(homedir(), ".jlassetregistry.json")))) == Dict(AssetRegistry.filekey(key) => [pwd(), 1])

    AssetRegistry.deregister(pwd())
    @test !AssetRegistry.isregistered(pwd())
    @test JSON.parse(String(read(joinpath(homedir(), ".jlassetregistry.json")))) == Dict()
end

@testset "register for not normalized paths" begin
    path = joinpath(pwd(), "..", "src")
    key = AssetRegistry.register(path)
    @test startswith(key, ENV["JULIA_ASSETREGISTRY_BASEURL"])
    @test key == AssetRegistry.getkey(path)
    @test AssetRegistry.isregistered(path)
    @test JSON.parse(String(read(joinpath(homedir(), ".jlassetregistry.json")))) == Dict(AssetRegistry.filekey(key) => [normpath(path), 1])
    AssetRegistry.deregister(path)
    @test !AssetRegistry.isregistered(path)
    @test JSON.parse(String(read(joinpath(homedir(), ".jlassetregistry.json")))) == Dict()
end


text:
"Abstract Type for canopy scattering"
abstract type AbstractCanopyScatteringType{FT<:AbstractFloat} end

"Model for bi-lambertian canopy leaf scattering"
Base.@kwdef struct BiLambertianCanopyScattering{FT<:AbstractFloat} <: AbstractCanopyScatteringType{FT}
    "Lambertian Reflectance"
    R::FT = FT(0.3)
    "Lambertian Transmission"
    T::FT  = FT(0.1)
    "Number of quadrature points in inclination angle"
    nQuad::Int = 10
end

"Model for specular canopy leaf scattering"
Base.@kwdef struct SpecularCanopyScattering{FT<:AbstractFloat} <: AbstractCanopyScatteringType{FT}
    "Refractive index"
    nᵣ::FT = FT(1.5)
    "Roughness parameter"
    κ::FT  = FT(0.5)
    "Number of quadrature points in azimuth"
    nQuad::Int = 20
end

"Abstract Type for leaf distributions"
abstract type AbstractLeafDistribution end

"""
    struct LeafDistribution{FT<:AbstractFloat}
A struct that defines the leaf angular distribution in radians (from 0->π/2; here scaled to [0,1])
# Fields
$(DocStringExtensions.FIELDS)
"""
struct LeafDistribution{FT<:AbstractFloat} <: AbstractLeafDistribution
    "Julia Univariate Distribution from Distributions.js"
    LD::UnivariateDistribution
    "Scaling factor to normalize distribution (here mostly 2/π as Beta distribution is from [0,1])"
    scaling::FT
end

####################################
"Abstract Type for leaf composition"
abstract type AbstractLeafProperties end

"""
    struct LeafProperties{FT}
A struct which stores important variables of leaf chemistry and structure
# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct LeafProspectProProperties{FT} <: AbstractLeafProperties
    ### Prospect related parameters
    "Leaf structure parameter [0-3]"
    N::FT       = FT(1.4  )
    "Chlorophyll a+b content `[µg cm⁻²]`"
    Ccab::FT    = FT(40.0)#u"µg/cm^2"
    "Carotenoid content `[µg cm⁻²]`"
    Ccar::FT    = FT(10.0)#u"µg/cm^2"
    "Anthocynanin content `[nmol cm⁻²]`"
    Canth::FT   = FT(0.5)#u"nmol/cm^2"
    "Brown pigments content in arbitrary units"
    Cbrown::FT = FT(0.0)
    "Equivalent water thickness `[cm]`, typical about 0.002-0.015"
    Cw::FT    = FT(0.012)#u"cm"
    "Dry matter content (dry leaf mass per unit area) `[g cm⁻²]`, typical about 0.003-0.016"
    Cm::FT   = FT(0.0)#u"g/cm^2"
    "protein content `[g/cm]`"
    Cprot::FT  = FT(0.001)#u"g/cm^2"
    "Carbone-based constituents content in `[g/cm⁻²]` (cellulose, lignin, sugars...)"
    Ccbc::FT   =  FT(0.009)#u"g/cm^2"
end

####################################
"Abstract type for Prospect model versions"
abstract type AbstractProspectProperties end

"""
    struct PigmentOpticalProperties{FT}
A struct which stores important absorption cross sections of pigments, water, etc
# Fields
$(DocStringExtensions.FIELDS)
"""
struct LeafOpticalProperties{FT,FT2} <: AbstractProspectProperties
    "Wavelength `[length]`"
    λ::FT2 #typeof(([1.0])u"nm")
    ### Prospect-PRO related parameters
    "Refractive index of leaf material "
    nᵣ::Array{FT, 1}
    "specific absorption coefficient of chlorophyll (a+b) `[cm² μg⁻¹]`"
    Kcab::Array{FT, 1}
    "specific absorption coefficient of carotenoids `[cm² μg⁻¹]`"
    Kcar::Array{FT, 1}
    "specific absorption coefficient of Anthocyanins `[cm² nmol⁻¹]`"
    Kant::Array{FT, 1}
    "specific absorption coefficient of brown pigments (arbitrary units)"
    Kb::Array{FT, 1}
    "specific absorption coefficient of water `[cm⁻¹]`"
    Kw::Array{FT, 1}
    "specific absorption coefficient of dry matter `[cm² g⁻¹]`"
    Km::Array{FT, 1}
    "specific absorption coefficient of proteins `[cm² g⁻¹]`"
    Kp::Array{FT, 1}
    "specific absorption coefficient of carbon based constituents `[cm² g⁻¹]`"
    Kcbc::Array{FT, 1}
end


text:
function __init__smsspam()

    DEPNAME = "SMSSpamCollection"
    LINK = "https://raw.githubusercontent.com/mohitgupta-omg/Kaggle-SMS-Spam-Collection-Dataset-/master/"
    LINK = "https://archive.ics.uci.edu/ml/machine-learning-databases/00228/"
    DOCS = "https://archive.ics.uci.edu/ml/datasets/SMS+Spam+Collection#"
    DATA = "smsspamcollection.zip"

    register(DataDep(
        DEPNAME,
        """
        Dataset: The SMS Spam Collection v.1
        Website: $DOCS
        """,
        LINK .* [DATA],
        "1587ea43e58e82b14ff1f5425c88e17f8496bfcdb67a583dbff9eefaf9963ce3",
        post_fetch_method = unpack
    ))
end

"""
    SMSSpamCollection(; dir=nothing)

The SMS Spam Collection v.1 (hereafter the corpus) is a set of SMS tagged messages
that have been collected for SMS Spam research. It contains one set of SMS messages
in English of 5,574 messages, tagged acording being ham (legitimate) or spam.
The corpus has a total of 4,827 SMS legitimate messages (86.6%) and a total of 747 (13.4%) spam messages.

The corpus has been collected by Tiago Agostinho de Almeida (http://www.dt.fee.unicamp.br/~tiago)
and José María Gómez Hidalgo (http://www.esp.uem.es/jmgomez).

```julia-repl
julia> using MLDatasets: SMSSpamCollection

julia> targets = SMSSpamCollection.targets();

julia> summary(targets)
"5574-element Vector{Any}"

julia> targets[1]
"ham"

julia> summary(features)
"5574-element Vector{Any}"
"""
struct SMSSpamCollection <: SupervisedDataset
    metadata::Dict{String,Any}
    features::Vector{String}
    targets::Vector{String}
end

function SMSSpamCollection(; dir = nothing)
    DEPNAME = "SMSSpamCollection"
    path = datafile(DEPNAME, "SMSSpamCollection", dir)
    spam_data = open(readlines, path)
    spam_data = [split(str, "\t") for str in spam_data]
    @assert all(x -> length(x)==2, spam_data)
    targets = [s[1] for s in spam_data]
    features = [s[2] for s in spam_data]

    metadata = Dict{String,Any}()
    metadata["n_observations"] = length(features)
    SMSSpamCollection(metadata, features, targets)
end

# DEPRECATED in V0.6
function Base.getproperty(::Type{SMSSpamCollection}, s::Symbol)
    if s == :features
        @warn "SMSSpamCollection.features() is deprecated, use `SMSSpamCollection().features` instead."
        return () -> SMSSpamCollection().features
    elseif s == :targets
        @warn "SMSSpamCollection.targets() is deprecated, use `SMSSpamCollection().targets` instead."
        return () -> SMSSpamCollection().targets
    else
        return getfield(Titanic, s)
    end
end


text:
"""
## module Potentials

### Summary

This module implements some basic interatomic potentials in pure Julia, as well
as provides building blocks and prototypes for further implementations
The implementation is done in such a way that they can be used either in "raw"
form or withinabstractframeworks.

### Types

### `evaluate`, `evaluate_d`, `evaluate_dd`, `grad`

### The `@D`, `@DD`, `@GRAD` macros

TODO: write documentation

"""
module Potentials

using JuLIP: Atoms, AbstractAtoms, AbstractCalculator,
      JVec, mat, vec, JMat, SVec, vecs, SMat,
      positions, set_positions!, fltype_intersect

using JuLIP.Chemistry: atomic_number

using StaticArrays: @SMatrix, SVector

using NeighbourLists

using LinearAlgebra: norm

using SparseArrays: sparse

import JuLIP: energy, forces, cutoff, virial, hessian_pos, hessian,
              site_energies, r_sum,
              site_energy, site_energy_d,
              energy!, forces!, virial!,
              alloc_temp, alloc_temp_d, alloc_temp_dd,
              read_dict, write_dict, fltype, rfltype

export PairPotential, SitePotential, ZeroSitePotential

# the following are prototypes for internal functions around which IPs are
# defined

function evaluate end
function evaluate_d end
function evaluate_ed end
function evaluate_dd end
function evaluate! end
function evaluate_d! end
function evaluate_dd! end
function precon! end

include("potentials_base.jl")
# * @pot, @D, @DD
# *

"""
`SitePotential`:abstractsupertype for generic site potentials
"""
abstract type SitePotential <: AbstractCalculator end

"""
`SimpleSitePotential`:abstractsupertype for generic site potentials,
ignoring species
"""
abstract type SimpleSitePotential <: SitePotential end

"""
`PairPotential`:abstractsupertype for pair potentials
"""
abstract type PairPotential <: SitePotential end

"""
`SimplePairPotential`:abstractsupertype for pair potentials,
ignoring species
"""
abstract type SimplePairPotential <: PairPotential end

abstract type ExplicitPairPotential <: SimplePairPotential end

# ---- redirect with some convenience functions ----

# redirect allocating -> non-allocating calls

evaluate(V::SitePotential, R, args...) =
      evaluate!(alloc_temp(V, length(R)), V, R, args...)

evaluate_d(V::SitePotential, R::AbstractVector{JVec{T}}, args...) where {T} =
      evaluate_d!(zeros(JVec{fltype_intersect(V, T)}, length(R)),
                  alloc_temp_d(V, length(R)),
                  V, R, args...)

evaluate_dd(V::SitePotential, R::AbstractVector{JVec{T}}, args...) where {T} =
      evaluate_dd!(zeros(JMat{fltype_intersect(V, T)}, length(R), length(R)),
                   alloc_temp_dd(V, length(R)),
                   V, R, args...)

# ----- interface for SimpleSitePotential

evaluate!(tmp, V::SimpleSitePotential, R, Z, z0) =
      evaluate!(tmp, V, R)
evaluate_d!(dEs, tmp, V::SimpleSitePotential, R, Z, z0) =
      evaluate_d!(dEs, tmp, V, R)
evaluate_dd!(hEs, tmp, V::SimpleSitePotential, R, Z, z0) =
      evaluate_dd!(hEs, tmp, V, R)
precon!(hEs, tmp, V::SimpleSitePotential, R, Z, z0, innerstab) =
      precon!(hEs, tmp, V, R, innerstab)

# ------- Neighbourlist related business -------------------

NeighbourLists.sites(at::AbstractAtoms, rcut::AbstractFloat) =
      sites(neighbourlist(at, rcut))

NeighbourLists.pairs(at::AbstractAtoms, rcut::AbstractFloat) =
      pairs(neighbourlist(at, rcut))

"""
`neigsz!(tmp, nlist::PairList, at::Atoms, i::Integer) -> j, R Z`

requires a temporary storage array `tmp` with fields
`tmp.R, tmp.Z`.
"""
function neigsz!(tmp, nlist::PairList, at::Atoms, i::Integer)
   j, R = neigs!(tmp.R, nlist, i)
   Z = tmp.Z
   for n = 1:length(j)
      Z[n] = at.Z[j[n]]
   end
   return j, R, (@view Z[1:length(j)])
end

function neigsz(nlist::PairList, at::Atoms, i::Integer)
   j, R = NeighbourLists.neigs(nlist, i)
   return j, R, at.Z[j]
end

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

"a site potential that just returns zero"
mutable struct ZeroSitePotential <: SitePotential
end

@pot ZeroSitePotential

cutoff(::ZeroSitePotential) = Bool(0)
energy(V::ZeroSitePotential, at::AbstractAtoms{T}; kwargs...) where T = zero(T)
forces(V::ZeroSitePotential, at::AbstractAtoms{T}; kwargs...) where T = zeros(JVec{T}, length(at))
evaluate!(tmp, p::ZeroSitePotential, args...) = Bool(0)
evaluate_d!(dEs, tmp, V::ZeroSitePotential, args...) = fill!(dEs, zero(eltype(dEs)))
evaluate_dd!(hEs, tmp, V::ZeroSitePotential, args...) = fill!(hEs, zero(eltype(hEs)))

# Implementation of a generic site potential
# ================================================

alloc_temp_site(N::Integer, T=Float64) =
      ( R = zeros(JVec{T}, N),
        Z = zeros(AtomicNumber, N), )

alloc_temp(V::SitePotential, at::AbstractAtoms) =
      alloc_temp(V, maxneigs(neighbourlist(at, cutoff(V))))

alloc_temp(V::SitePotential, N::Integer) =
      ( R = zeros(JVecF, N),
        Z = zeros(AtomicNumber, N), )

alloc_temp_d(V::SitePotential, at::AbstractAtoms) =
      alloc_temp_d(V, maxneigs(neighbourlist(at, cutoff(V))))

alloc_temp_d(V::SitePotential, N::Integer) =
      (dV = zeros(JVec{fltype(V)}, N),
        R = zeros(JVecF, N),
        Z = zeros(AtomicNumber, N), )

alloc_temp_dd(V::SitePotential, args...) = nothing

# -------------- Implementations of energy, forces, virials
#                for a generic site potential

energy(V::SitePotential, at::AbstractAtoms; kwargs...) =
      energy!(alloc_temp(V, at), V, at; kwargs...)

virial(V::SitePotential, at::AbstractAtoms; kwargs...) =
      virial!(alloc_temp_d(V, at), V, at; kwargs...)

forces(V::SitePotential, at::AbstractAtoms; kwargs...) =
      forces!(zeros(JVec{fltype_intersect(V, at)}, length(at)),
              alloc_temp_d(V, at), V, at; kwargs...)

function energy!(tmp, calc::SitePotential, at::Atoms;
                 domain=1:length(at))
   TFL = fltype_intersect(calc, at)
   E = zero(TFL)
   nlist = neighbourlist(at, cutoff(calc))
   for i in domain
      j, R, Z = neigsz!(tmp, nlist, at, i)
      E += evaluate!(tmp, calc, R, Z, at.Z[i])
   end
   return E
end

function forces!(frc, tmp, calc::SitePotential, at::Atoms;
                 domain=1:length(at), reset=true)
   TFL = fltype_intersect(calc, at)
   if reset; fill!(frc, zero(eltype(frc))); end
   nlist = neighbourlist(at, cutoff(calc))
   for i in domain
      j, R, Z = neigsz!(tmp, nlist, at, i)
      if length(j) > 0
         evaluate_d!(tmp.dV, tmp, calc, R, Z, at.Z[i])
         for a = 1:length(j)
            frc[j[a]] -= tmp.dV[a]
            frc[i]    += tmp.dV[a]
         end
      end
   end
   return frc
end

site_virial(dV::AbstractVector{JVec{T1}}, R::AbstractVector{JVec{T2}}
            ) where {T1, T2} =  (
      length(R) > 0 ? (- sum( dVi * Ri' for (dVi, Ri) in zip(dV, R) ))
                    : zero(JMat{fltype_intersect(T1, T2)})
      )

function virial!(tmp, calc::SitePotential, at::Atoms; domain=1:length(at))
   TFL = fltype_intersect(calc, at)
   nlist = neighbourlist(at, cutoff(calc))
   vir = zero(JMat{TFL})
   for i in domain
      j, R, Z = neigsz!(tmp, nlist, at, i)
      if length(j) > 0
         evaluate_d!(tmp.dV, tmp, calc, R, Z, at.Z[i])
         vir += site_virial(tmp.dV, R)
      end
   end
   return vir
end

function site_energies(V::SitePotential, at::AbstractAtoms; kwargs...)
   TFL = fltype_intersect(V, at)
   return site_energies!(zeros(TFL, length(at)),
                         alloc_temp(V, at), V, at; kwargs...)
end

function site_energies!(Es, tmp, V::SitePotential, at::AbstractAtoms;
                        domain = 1:length(at))
   nlist = neighbourlist(at, cutoff(V))
   for i in domain
      _j, R, Z = neigsz!(tmp, nlist, at, i)
      Es[i] = evaluate!(tmp, V, R, Z, at.Z[i])
   end
   return Es
end

site_energy(V::SitePotential, at::AbstractAtoms, i0::Integer) =
      energy(V, at; domain = (i0,))

site_energy_d(V::SitePotential, at::AbstractAtoms, i0::Integer) =
      rmul!(forces(V, at; domain = (i0,)), -one(fltype(at)))

# ------------------------------------------------
#  specialisation for Pair potentials

include("analyticpotential.jl")

include("cutoffs.jl")

include("pairpotentials.jl")

include("adsite.jl")

include("stillingerweber.jl")

include("splines.jl")
include("eam.jl")

include("onebody.jl")

include("hessians.jl")

include("emt.jl")

end


text:
@testset verbose=true "Math" begin
    components = [:x, :y, :z]

    basic_functions = [conj, abs2, abs, inv, log, exp, sqrt]

    ### Test equivalence with complex math
    ℍ_to_ℂ(q::Real, component) = q + 0im
    function ℍ_to_ℂ(q::AbstractQuaternion, component)
        c = q.w + im * getproperty(q, component)
    end
    function ℂ_to_ℍ(c::Complex, component)
        Quaternion(c.re, [comp_i == component ? c.im : zero(real(c)) for comp_i in components]...)
    end
    @testset "Complex equivalence $T" for T in FloatTypes
        ϵ = 10eps(T)
        scalars = [zero(T), one(T), -one(T)]
        for c in [a+b*im for a in scalars for b in scalars]
            # The following constructs a quaternion from the complex number `c` by equating `im`
            # with one of the three quaternionic imaginaries, looping over each option.  We then
            # apply the quaternionic math function to that quaternion, and map back to ℂ; the result
            # should be identical to just applying the complex math function directly.  The only
            # exceptions are for `log` and `sqrt` when the complex number is a negative real number,
            # for which we have to arbitrarily choose the "vector" part of the quaternion, which is
            # chosen to be `imz` in the code.
            for component in components
                q = ℂ_to_ℍ(c, component)
                for unary_function in basic_functions
                    if unary_function ∈ [log, sqrt] && component!=:z && c.re < zero(T) && c.im == zero(T)
                        continue  # We arbitrarily chose the z component for these return values
                    end
                    cval = unary_function(c)
                    qval = ℍ_to_ℂ(unary_function(q), component)
                    @test cval ≈ qval rtol=ϵ nans=true
                    # Repeat the test on QuatVec for `exp` for pure-imaginary input
                    if unary_function ∈ [exp] && iszero(c.re)
                        qval = ℍ_to_ℂ(unary_function(QuatVec(q)), component)
                        @test cval ≈ qval rtol=ϵ nans=true
                    end
                    # Repeat the test on Rotor for `log` and `sqrt` for nonzero input
                    if unary_function ∈ [log, sqrt] && !iszero(c)
                        cval = unary_function(c/abs(c))
                        qval = ℍ_to_ℂ(unary_function(Rotor(q)), component)
                        @test cval ≈ qval rtol=ϵ nans=true
                    end
                end
                for unary_function in [angle]  # Not identical, but related
                    cval = unary_function(c)
                    qval = ℍ_to_ℂ(unary_function(q), component)
                    @test 2*abs(cval) ≈ qval rtol=ϵ nans=true
                    if !iszero(c)
                        cval = unary_function(c/abs(c))
                        qval = ℍ_to_ℂ(unary_function(Rotor(q)), component)
                        @test 2*abs(cval) ≈ qval rtol=ϵ nans=true
                    end
                end
                for binary_function in [(^)]
                    if binary_function ∈ [(^)] && ((component!=:z && c.re < zero(T) && c.im == zero(T)) || iszero(c))
                        continue  # We arbitrarily chose the z component for these return values
                    end
                    if T !== BigFloat
                        S = 10randn(T, 100)
                        for s ∈ S
                            cval = binary_function(c, float(s))
                            qval = ℍ_to_ℂ(binary_function(q, s), component)
                            @test cval ≈ qval rtol=6ϵ nans=true
                            if !iszero(c)
                                cval = binary_function(c/abs(c), s)
                                qval = ℍ_to_ℂ(binary_function(Rotor(q), s), component)
                                @test cval ≈ qval rtol=6ϵ nans=true
                            end
                        end
                    end
                    S = [-3, -2, -1, 0, 1, 2, 3]
                    for s ∈ S
                        cval = binary_function(c, s)
                        qval = ℍ_to_ℂ(binary_function(q, s), component)
                        @test cval ≈ qval rtol=ϵ nans=true
                        if iszero(c.re)
                            qval = ℍ_to_ℂ(binary_function(QuatVec(q), s), component)
                            @test cval ≈ qval rtol=ϵ nans=true
                        end
                        if !iszero(c)
                            cval = binary_function(c/abs(c), s)
                            qval = ℍ_to_ℂ(binary_function(Rotor(q), s), component)
                            @test cval ≈ qval rtol=ϵ nans=true
                        end
                    end
                end

                @test exp(log(q)) ≈ q rtol=ϵ nans=true
                @test exp(zero(q)) ≈ one(q) rtol=ϵ nans=true
                @test exp(log(Rotor(q))) ≈ Rotor(q) rtol=ϵ nans=true
                @test exp(zero(QuatVec(q))) ≈ one(q) rtol=ϵ nans=true

            end
        end
    end

end


text:
#
# Factors Broadcasting
#

# THIS IS VERY POORLY WRITTEN

"""
    broadcast(f, ϕ::Factor, dims, values)

Broadcast a vector (or array of vectors) across the dimension(s) `dims`

See `Base.broadcast` for more info.
"""
Base.broadcast(f, ϕ::Factor, dims, values) = broadcast!(f, deepcopy(ϕ), dims, values)

"""
    broadcast!(f, ϕ::Factor, dims, values)

Broadcast a vector (or array of vectors) across the dimension(s) `dims`

See `Base.broadcast!` for more info.
"""
Base.broadcast!(f, ϕ::Factor, dim::Symbol, values) = Base.broadcast!(f, ϕ, [dim], [values])

function Base.broadcast!(f, ϕ::Factor, dims::Vector{Symbol}, values)
    _check_dims_valid(dims, ϕ)

    (length(dims) != length(values)) &&
        throw(DimensionMismatch("Number of dimensions does not match number of values to broadcast"))

    # broadcast will check if the dimensions of each value are valid
    # check all the dimensions match or its a singleton value
    #any( (length(ϕ, dims) .!= map(length, values)) & map(length, values) .!= 1) ) &&
    #    throw(ArgumentError("Length of dimensions don't match lengths of broadcast array"))

    inds = indexin(dims, ϕ)

    # each vector gets reshaped to be 1x1... except for the index of its dimension
    reshape_dims = ones(Int, ndims(ϕ))
    new_values = Vector{Array}(length(values))

    for (i, val) in enumerate(values)
        if isa(val, Vector)
            # reshape to the proper dimension
            dim_loc = inds[i]
            reshape_dims[dim_loc] = length(val)
            new_values[i] = reshape(val, reshape_dims...)
            reshape_dims[dim_loc] = 1
        else # its a number (?)
            new_values[i] = [val]
        end
    end

    broadcast!(f, ϕ.potential, ϕ.potential, new_values...)

    return ϕ
end


text:
cd(@__DIR__)
using Pkg
pkg"activate ."

using Revise

const ROOT_PATH = pwd()

"""
    bootstrap_genie() :: Nothing

Bootstraps the Genie framework setting up paths and workers. Invoked automatically.
"""
function bootstrap() :: Nothing
  printstyled("""
   _____         _
  |   __|___ ___|_|___
  |  |  | -_|   | | -_|
  |_____|___|_|_|_|___|

  """, color = :magenta)

  DEFAULT_NWORKERS_REPL = 1
  DEFAULT_NWORKERS_SERVER = 1

  isfile("env.jl") && include("env.jl")
  haskey(ENV, "GENIE_ENV") || (ENV["GENIE_ENV"] = "dev")
  @info "\nStarting Genie in >> $(ENV["GENIE_ENV"] |> uppercase) << mode \n\n"

  push!(LOAD_PATH, pwd(), "src")

  nothing
end

bootstrap()

using Genie, Genie.App, Genie.Toolbox

Genie.load(context = @__MODULE__)

Genie.run()


text:
using FractionalDiffEq
h=0.005; tf=(0, 30)
alpha = [0.99, 0.99, 0.99]
x0 = [1, 0, 1]

function f(du, u)
    a, b, c = 10, 28, 8/3
    du[1] = a*(u[2]-u[1])
    du[2] = u[1]*(b-u[3])-u[2]
    du[3] = u[1]*u[2]-c*u[3]
end
prob = FODESystem(f, alpha, x0, tf)
result = solve(prob, h, PIEX())

using Plots
plot(result[:, 1], result[:, 2])