text:
"""
    struct Quat{T} <: Rotation{3,T}
    Quat(w, x, y, z)

The `Quat` type is a 3×3 matrix representation of a normalized quaternion.
They allow you to transparently use (fast) quaternion algebra to store, compose
and invert 3D rotations, while at the same time letting you apply rotations
through matrix-vector multiplication.

Note: by default, the constructor will renormalize the input so that the quaternion
has length 1 (w² + x² + y² + z² = 1), and the rotation matrix is orthogonal.

Renormalization can be skipped by passing `false` as an additional constructor
argument, in which case the user provides the guarantee that the input arguments
represent a unit quaternion. Operations on an unnormalized `Quat`, created by
skipping renormalization in this fashion, are not guaranteed to do anything sensible.
"""
struct Quat{T} <: Rotation{3,T}
    w::T
    x::T
    y::T
    z::T

    @inline function Quat{T}(w, x, y, z, normalize::Bool = true) where {T}
        if normalize
            norm = sqrt(w*w + x*x + y*y + z*z)
            new(w/norm, x/norm, y/norm, z/norm)
        else
            new(w, x, y, z)
        end
    end

    Quat{T}(q::Quat) where {T} = new{T}(q.w, q.x, q.y, q.z)
end

@inline function Quat(w::W, x::X, y::Y, z::Z, normalize::Bool = true) where {W, X, Y, Z}
    Quat{promote_type(promote_type(promote_type(W, X), Y), Z)}(w, x, y, z, normalize)
end

# These 3 functions are enough to satisfy the entire StaticArrays interface:
function (::Type{Q})(t::NTuple{9}) where Q<:Quat
    #=
    This function solves the system of equations in Section 3.1
    of https://arxiv.org/pdf/math/0701759.pdf. This cheap method
    only works for matrices that are already orthonormal (orthogonal
    and unit length columns). The nearest orthonormal matrix can
    be found by solving Wahba's problem:
    https://en.wikipedia.org/wiki/Wahba%27s_problem as shown below.

    not_orthogonal = randn(3,3)
    u,s,v = svd(not_orthogonal)
    is_orthogoral = u * diagm([1, 1, sign(det(u * transpose(v)))]) * transpose(v)
    =#

    a = 1 + t[1] + t[5] + t[9]
    b = 1 + t[1] - t[5] - t[9]
    c = 1 - t[1] + t[5] - t[9]
    d = 1 - t[1] - t[5] + t[9]
    max_abcd = max(a, b, c, d)
    if a == max_abcd
        b = t[6] - t[8]
        c = t[7] - t[3]
        d = t[2] - t[4]
    elseif b == max_abcd
        a = t[6] - t[8]
        c = t[2] + t[4]
        d = t[7] + t[3]
    elseif c == max_abcd
        a = t[7] - t[3]
        b = t[2] + t[4]
        d = t[6] + t[8]
    else
        a = t[2] - t[4]
        b = t[7] + t[3]
        c = t[6] + t[8]
    end
    return Q(a, b, c, d)
end

function Base.getindex(q::Quat, i::Int)
    if i == 1
        ww = (q.w * q.w)
        xx = (q.x * q.x)
        yy = (q.y * q.y)
        zz = (q.z * q.z)

        ww + xx - yy - zz
    elseif i == 2
        xy = (q.x * q.y)
        zw = (q.w * q.z)

        2 * (xy + zw)
    elseif i == 3
        xz = (q.x * q.z)
        yw = (q.y * q.w)

        2 * (xz - yw)
    elseif i == 4
        xy = (q.x * q.y)
        zw = (q.w * q.z)

        2 * (xy - zw)
    elseif i == 5
        ww = (q.w * q.w)
        xx = (q.x * q.x)
        yy = (q.y * q.y)
        zz = (q.z * q.z)

        ww - xx + yy - zz
    elseif i == 6
        yz = (q.y * q.z)
        xw = (q.w * q.x)

        2 * (yz + xw)
    elseif i == 7
        xz = (q.x * q.z)
        yw = (q.y * q.w)

        2 * (xz + yw)
    elseif i == 8
        yz = (q.y * q.z)
        xw = (q.w * q.x)

        2 * (yz - xw)
    elseif i == 9
        ww = (q.w * q.w)
        xx = (q.x * q.x)
        yy = (q.y * q.y)
        zz = (q.z * q.z)

        ww - xx - yy + zz
    else
        throw(BoundsError(r,i))
    end
end

function Base.Tuple(q::Quat)
    ww = (q.w * q.w)
    xx = (q.x * q.x)
    yy = (q.y * q.y)
    zz = (q.z * q.z)
    xy = (q.x * q.y)
    zw = (q.w * q.z)
    xz = (q.x * q.z)
    yw = (q.y * q.w)
    yz = (q.y * q.z)
    xw = (q.w * q.x)

    # initialize rotation part
    return (ww + xx - yy - zz,
            2 * (xy + zw),
            2 * (xz - yw),
            2 * (xy - zw),
            ww - xx + yy - zz,
            2 * (yz + xw),
            2 * (xz + yw),
            2 * (yz - xw),
            ww - xx - yy + zz)
end

function Base.:*(q::Quat, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    qo = (-q.x * v[1] - q.y * v[2] - q.z * v[3],
           q.w * v[1] + q.y * v[3] - q.z * v[2],
           q.w * v[2] - q.x * v[3] + q.z * v[1],
           q.w * v[3] + q.x * v[2] - q.y * v[1])

    T = promote_type(eltype(q), eltype(v))

    return similar_type(v, T)(-qo[1] * q.x + qo[2] * q.w - qo[3] * q.z + qo[4] * q.y,
                              -qo[1] * q.y + qo[2] * q.z + qo[3] * q.w - qo[4] * q.x,
                              -qo[1] * q.z - qo[2] * q.y + qo[3] * q.x + qo[4] * q.w)
end

# TODO consider Ac_mul_B, etc...
function Base.:*(q1::Quat, q2::Quat)
    Quat(q1.w*q2.w - q1.x*q2.x - q1.y*q2.y - q1.z*q2.z,
            q1.w*q2.x + q1.x*q2.w + q1.y*q2.z - q1.z*q2.y,
            q1.w*q2.y - q1.x*q2.z + q1.y*q2.w + q1.z*q2.x,
            q1.w*q2.z + q1.x*q2.y - q1.y*q2.x + q1.z*q2.w)
end

function Base.inv(q::Quat)
    Quat(q.w, -q.x, -q.y, -q.z)
end

@inline Base.one(::Type{Quat}) = Quat(1.0, 0.0, 0.0, 0.0, false)
@inline Base.one(::Type{Quat{T}}) where {T} = Quat{T}(one(T), zero(T), zero(T), zero(T), false)

"""
    rotation_between(from, to)

Compute the quaternion that rotates vector `from` so that it aligns with vector
`to`, along the geodesic (shortest path).
"""
rotation_between(from::AbstractVector, to::AbstractVector) = rotation_between(SVector{3}(from), SVector{3}(to))
function rotation_between(from::SVector{3}, to::SVector{3})
    # Robustified version of implementation from https://www.gamedev.net/topic/429507-finding-the-quaternion-betwee-two-vectors/#entry3856228
    normprod = sqrt(dot(from, from) * dot(to, to))
    T = typeof(normprod)
    normprod < eps(T) && throw(ArgumentError("Input vectors must be nonzero."))
    w = normprod + dot(from, to)
    v = abs(w) < 100 * eps(T) ? perpendicular_vector(from) : cross(from, to)
    @inbounds return Quat(w, v[1], v[2], v[3]) # relies on normalization in constructor
end

################################################################################
################################################################################
"""
    struct SPQuat{T} <: Rotation{3,T}
    SPQuat(x, y, z)

An `SPQuat` is a 3D rotation matrix represented by the "stereographic projection" of a normalized quaternion (shortened to "SPQuat"), which is
a 3-element parametrization of a unit quaternion Q formed by the intersection of a line from [-1,0,0,0] to Q, with a plane containing the origin and with normal direction [1,0,0,0]. This
is a compact representation of rotations where the derivatives of the rotation matrix's elements w.r.t. the SPQuat parameters are rational functions (making them useful for optimization).

See:

    [1] "A Recipe on the Parameterization of Rotation Matrices for Non-Linear Optimization using Quaternions",
    Terzakis, Culverhouse, Bugmann, Sharma, Sutton,
    MIDAS technical report, 2011
    http://www.tech.plymouth.ac.uk/sme/springerusv/2011/publications_files/Terzakis%20et%20al%202012,%20A%20Recipe%20on%20the%20Parameterization%20of%20Rotation%20Matrices...MIDAS.SME.2012.TR.004.pdf

    Note 1: the singularity (origin) has been moved from [0,0,0,-1] in Ref[1] to [-1,0,0,0], so the 0 rotation quaternion [1,0,0,0] maps to [0,0,0] as opposed of to [1,0,0].
    Note 2: the quaternion rotation ambiguity q = -q means there will be a rotation with ||SPQuat|| <= 1 and an equivilent rotation with ||SPQuat|| > 1.  This package
            uses the solution with ||SPQuat|| <= 1
    Note 3: it is safe to assume that the corresponding matrix is orthogonal/unitary for any input x, y, z.

"""
struct SPQuat{T} <: Rotation{3,T}
    x::T
    y::T
    z::T

    # TODO should we enforce norm <= 1?
    SPQuat{T}(x, y, z) where {T} = new{T}(x, y, z)
    SPQuat{T}(spq::SPQuat) where {T} = new{T}(spq.x, spq.y, spq.z)
end

@inline SPQuat(x::X, y::Y, z::Z) where {X,Y,Z} = SPQuat{promote_type(promote_type(X, Y), Z)}(x, y, z)

# These functions are enough to satisfy the entire StaticArrays interface:
@inline (::Type{SPQ})(t::NTuple{9}) where {SPQ <: SPQuat} = convert(SPQ, Quat(t))
@inline Base.getindex(spq::SPQuat, i::Int) = convert(Quat, spq)[i]
@inline Base.Tuple(spq::SPQuat) = Tuple(convert(Quat, spq))

@inline function (::Type{Q})(spq::SPQuat) where Q <: Quat
    # Equation (45) in
    # Terzakis et al., "A Recipe on the Parameterization of Rotation Matrices
    # for Non-Linear Optimization using Quaternions":
    alpha2 = spq.x * spq.x + spq.y * spq.y + spq.z * spq.z
    scale = 2 / (alpha2 + 1)
    Q((1 - alpha2) / (alpha2 + 1), scale * spq.x, scale * spq.y, scale * spq.z, false)
end

@inline function (::Type{SPQ})(q::Quat) where SPQ <: SPQuat
    # Simplification of (46) and (47) in
    # Terzakis et al., "A Recipe on the Parameterization of Rotation Matrices
    # for Non-Linear Optimization using Quaternions":
    # α² = (1 - q.w) / (1 + q.w)
    # scale = (α² + 1) / 2
    # scale = ((1 - q.w) / (1 + q.w) + 1) / 2
    # scale = ((1 - q.w + (1 + q.w)) / (1 + q.w) / 2
    # scale = (2 / (1 + q.w)) / 2
    # scale = 1 / (1 + q.w)
    scale = 1 / (1 + q.w)
    SPQ(q.x * scale,  q.y * scale, q.z * scale)
end

@inline Base.:*(spq::SPQuat, x::StaticVector) = Quat(spq) * x

@inline Base.:*(spq::SPQuat, r::Rotation) = Quat(spq) * r
@inline Base.:*(spq::SPQuat, r::RotMatrix) = Quat(spq) * r
@inline Base.:*(r::Rotation, spq::SPQuat) = r * Quat(spq)
@inline Base.:*(r::RotMatrix, spq::SPQuat) = r * Quat(spq)
@inline Base.:*(spq1::SPQuat, spq2::SPQuat) = Quat(spq1) * Quat(spq2)

@inline Base.inv(spq::SPQuat) = SPQuat(-spq.x, -spq.y, -spq.z)

@inline Base.one(::Type{SPQuat}) = SPQuat(0.0, 0.0, 0.0)
@inline Base.one(::Type{SPQuat{T}}) where {T} = SPQuat{T}(zero(T), zero(T), zero(T))


text:
include("comlineoption.jl")
include("iteration.jl")
include("saveresult.jl")
using .Comlineoption
using .Iteration
using .Saveresult

function main(args)
    opt = Comlineoption.construct(args)
    data, solve_tf_param, solve_tf_val, yarray = Iteration.construct(opt)
    xarray = Iteration.iteration!(data, solve_tf_param, solve_tf_val, yarray)
    Saveresult.saveresult(data, xarray, yarray)
end

@time main(ARGS)

ext: jl


text:
function parse_ctfile(ctfilename::String)
    ctstring = read(ctfilename, String)
    parse_ctstring(ctstring)
end

function parse_ctstring(ctstring::String)
    pg_rgx = r"PG:(\w*)"
    ops_rgx = r"Ops:\s*\[([\w,]*)\]"
    irreps_rgx = r"Irreps:\s*\[([\w,]*)\]"
    chars_rgx = r"#Ordered Characters\s*([-\d\s,]*)\s*"
    pg = match(pg_rgx, ctstring).captures[1]
    classes = split(match(ops_rgx, ctstring).captures[1], ",")
    irreps = split(match(irreps_rgx, ctstring).captures[1], ",")
    chars = split(match(chars_rgx, ctstring).captures[1], r"[,\s]+")
    len = size(ops,1)
    char_arr = zeros(len, len)
    for i = 1:len
        for j = 1:len
            char_arr[i,j] = eval(Meta.parse(chars[len*(j-1)+i]))
        end
    end
    symels = ops_to_symels(ops)
    ct = Chartable(pg, symels, irreps, char_arr)
    return ct
end

function pg_to_symels(PG)
    pg = parse_pg_str(PG)
    symels = [Symel("E", [1 0 0; 0 1 0; 0 0 1])]
    σh = [1 0 0; 0 1 0; 0 0 -1]
    if pg.family == "C"
        if pg.subfamily == "h"
            push!(symels, Symel("sigmah", σh)) # sigmah
            if pg.n % 2 == 0
                push!(symels, Symel("i", i()))
            end
            cns = generate_Cn(pg.n)
            sns = generate_Sn(pg.n)
            symels = vcat(symels, cns, sns)
        elseif pg.subfamily == "v"
            cns = generate_Cn(pg.n)
            if pg.n % 2 == 0
                n = div(pg.n, 2)
                σds = generate_σd(n)
            else
                n = pg.n
                σds = []
            end
            σvs = generate_σv(n)
            symels = vcat(symels, cns, σvs, σds)
        elseif pg.subfamily == "s"
            push!(symels, Symel("sigmah", σh))
        elseif pg.subfamily == "i"
            push!(symels, Symel("i", i()))
        elseif isnothing(pg.subfamily)
            cns = generate_Cn(pg.n)
            symels = vcat(symels, cns)
        else
            throw(ArgumentError("Unidentified Point Group!"))
        end
    elseif pg.family == "D"
        if pg.subfamily == "h"
            push!(symels, Symel("sigmah", σh))
            if pg.n % 2 == 0
                push!(symels, Symel("i", i()))
                n = div(pg.n, 2)
                σds = generate_σd(n)
                c2s = generate_C2(pg.n)
            else
                n = pg.n
                σds = []
                c2s = generate_C2(pg.n)
            end
            cns = generate_Cn(pg.n)
            sns = generate_Sn(pg.n)
            σvs = generate_σv(n)
            #c2s = generate_C2(pg.n)
            symels = vcat(symels, cns, sns, σvs, σds, c2s)
        elseif pg.subfamily == "d"
            if pg.n % 2 != 0
                push!(symels, Symel("i", i()))
            end
            cns = generate_Cn(pg.n)
            sns = generate_Sn(pg.n * 2, true)
            σds = generate_σd(pg.n)
            c2s = generate_C2(pg.n)
            symels = vcat(symels, cns, sns, σds, c2s)
        elseif isnothing(pg.subfamily)
            cns = generate_Cn(pg.n)
            c2s = generate_C2(pg.n)
            symels = vcat(symels, cns, c2s)
        else
            throw(ArgumentError("Oh shit, the trout population"))
        end
    elseif pg.family == "S"
        if isnothing(pg.subfamily) & (pg.n % 2 == 0)
            n = div(pg.n, 2)
            if n % 2 != 0
                push!(symels, Symel("i", i()))
            end
            cns = generate_Cn(n)
            sns = generate_Sn(pg.n, true)
            symels = vcat(symels, cns, sns)
        else
            throw(ArgumentError("Oh shit, the trout population"))
        end
    else
        if pg.family == "T"
            if pg.subfamily == "h"
                Ths = generate_Th()
                symels = vcat(symels, Ths)
            elseif pg.subfamily == "d"
                Tds = generate_Td()
                symels = vcat(symels,Tds)
            else
                Ts = generate_T()
                symels = vcat(symels,Ts)
            end
        elseif pg.family == "O"
            if pg.subfamily == "h"
                Ohs = generate_Oh()
                symels = vcat(symels, Ohs)
            else
                Os = generate_O()
                symels = vcat(symels, Os)
            end
        elseif pg.family == "I"
            if pg.subfamily == "h"
                Ihs = generate_Ih()
                symels = vcat(symels, Ihs)
            else
                Is = generate_I()
                symels = vcat(symels, Is)
            end
        else
            throw(ArgumentError("Unidentified Point Group!"))
        end
    end
    return symels
end

function parse_pg_str(s)
    re = r"([A-Z]+)(\d+)?([a-z]+)?"
    m = match(re, s)
    family, n, subfamily = m.captures
    if !isnothing(n)
        n = parse(Int, n)
    end
    if !isnothing(subfamily)
        subfamily = string(subfamily)
    end
    family = string(family)
    return PG(s, family, n, subfamily)
end

function pg_to_chartab(PG)
    pg = parse_pg_str(PG)
    irreps = []
    if pg.family == "C"
        if pg.subfamily == "v"
            Cnv_irr(pg.n)
        elseif pg.subfamily == "h"
            Cnh_irr(pg.n)
        else
            Cn_irrmat(pg.n)
        end
    elseif pg.family == "D"
        if pg.subfamily == "d"
            Dnd_irr(pg.n)
        elseif pg.subfamily == "h"
            Dnh_irr(pg.n)
        else
            Dn_irr(pg.n)
        end
    elseif pg.family == "S"
        Sn_irr(pg.n)
    else
        if pg.family == "T"
            if pg.subfamily == "h"
                println("Th")
            elseif pg.subfamily == "d"
                println("Td")
            else
                println("T")
            end
        elseif pg.family == "O"
            if pg.subfamily == "h"
                println("Oh")
            else
                println("O")
            end
        elseif pg.family == "I"
            if pg.subfamily == "h"
                println("Ih")
            else
                println("I")
            end
        else
            throw(ArgumentError("Unrecognized Point Group"))
        end
    end
end

function char_table_from_str(PG)
    symels = pg_to_symels(PG)
    mtab = build_mult_table(symels)
    println(mtab)
    classes = find_classes(mtab)
    println(classes)
    M = build_Mk(mtab, classes)
    #F = ajd(M; algorithm=:JADE)
    #println(F.ev)
    #G = zeros(length(classes), length(classes))
    #for i = 1:length(classes)
    #    G[:,i] += normalize(F.F[:,i])
    #end
    #println(G)
    ##dixons_alg(classes, M, length(symels))
end


text:
## Old pressure solver

function volume(model::EoSModel, p, T,  z=SA[1.]; phase = "unknown")
    N = length(p)

    ub = [Inf]
    lb = lb_volume(model,z; phase = phase)

    x0 = x0_volume(model,p,T,z; phase = phase)
    f = v -> eos(model, exp10(v[1]), T,z) + exp10(v[1])*p
    #looking at best phase using tunneling
    if phase == "unknown"
        (f_best,v_best) = Solvers.tunneling(f,lb,ub,x0)
        #@show eos.(Ref(model),exp10.(v_best),T,Ref(z)) .- p.*exp10.(v_best)

        return exp10(v_best[1])
    else
        _ub = 100*T/p
        _lb = exp10(only(lb))
        f0(vx) = pressure(model,vx,T,z,phase=phase) - p
        _v0 = exp10(only(x0))
        #try direct newton solving
        vobj = Solvers.ad_newton(f0,_v0)
        if (_lb <= vobj <= _ub)
            return vobj
        else
            vobj = Roots.find_zero(f0,(_lb,_ub),FalsePosition())
        end
        #=
        opt_min = NLopt.Opt(:LD_MMA, length(ub))
        opt_min.lower_bounds = lb
        opt_min.upper_bounds = ub
        opt_min.xtol_rel     = 1e-8
        obj_f0 = x -> f(x)
        obj_f  = (x,g) -> Solvers.NLopt_obj(obj_f0,x,g)
        opt_min.min_objective =  obj_f
        (f_min,v_min) = NLopt.optimize(opt_min, x0)
        #@show eos.(Ref(model),exp10.(v_min),T,Ref(z)) .- p.*exp10.(v_min)
        return exp10(v_min[1])
        =#
    end
end

function Jac_Sat(model::EoSModel, J, T, v_l, v_v)
    #components = model.components
    fun(x) = eos(model, x[2], T,SA[x[1]])
    d2f(x) = ForwardDiff.hessian(fun,x)
    d2f_l = d2f(SA[one(v_l),v_l])
    d2f_v = d2f(SA[one(v_v),v_v])
    (p_scale,μ_scale) = scale_sat_pure(model)
    J[1,1] =  v_l[1]*d2f_l[2,2]*log(10)*p_scale
    J[1,2] = -v_v[1]*d2f_v[2,2]*log(10)*p_scale
    J[2,1] =  v_l[1]*d2f_l[1,2]*log(10)*μ_scale
    J[2,2] = -v_v[1]*d2f_v[1,2]*log(10)*μ_scale
end

function Obj_Jac_sat(model::EoSModel, F, J, T, v_l, v_v)
    Obj_Sat(model, F, T, v_l, v_v)
    Jac_Sat(model, J, T, v_l, v_v)
    F, J
end

function Jvop_sat(x,model::EoSModel,T)
    function Jac_satV(Fv, v)
        #components = model.components
        fun(x) = eos(model, x[2], T,[x[1]])
        d2f(z) = ForwardDiff.hessian(fun,z)
        d2f_l = d2f([1,x[1]])
        d2f_v = d2f([1,x[2]])
        (p_scale,μ_scale) = scale_sat_pure(model)
        Fv[1,] =  x[1]*d2f_l[2,2]*log(10)*p_scale*v[1]-x[2]*d2f_v[2,2]*log(10)*p_scale*v[2]
        Fv[2,] =  x[1]*d2f_l[1,2]*log(10)*μ_scale*v[1]-x[2]*d2f_v[1,2]*log(10)*μ_scale*v[2]
    end
    LinearMap(Jac_satV, length(x))
end

## Mixture critical point solver
# function crit_mix(model::SAFT,x_c)
#     components = model.components
#     z  = create_z(model,x_c)
#     f! = (F,x) -> Obj_Crit_mix(model, F, exp10(x[2]), x[1]*prod(model.params.epsilon[i]^z[i] for i in components), x_c)
#     x0 = [1.5,log10(π/6*N_A*sum(z[i]*model.params.segment[i]*model.params.sigma[i]^3 for i in components)/0.15)]
#     r  = Solvers.nlsolve(f!,x0)
#     T_c = r.zero[1]*prod(model.params.epsilon[i]^z[i] for i in components)
#     v_c = exp10(r.zero[2])
#     p_c = pressure(model, v_c, T_c, x_c)
#     return (T_c, p_c, v_c)
# end
# #
# function Obj_Crit_mix(model::SAFT, F, v_c,T_c,x_c)
#     fun(x)  = eos(model, create_z(model, [x[1],1-x[1]]), v_c, T_c)
#     df(x)   = ForwardDiff.derivative(fun,x)
#     d2f(x)  = ForwardDiff.derivative(df,x)
#     d3f(x)  = ForwardDiff.derivative(d2f,x)
#     F[1] = d2f(x_c[1])
#     F[2] = d3f(x_c[1])
#     println(F)
#     println(v_c)
#     println(T_c)
# end

function Jac_bubble_pressure(model, J, T, v_l, v_v, x, y)
    components = model.components
    append!(y,1-sum(y[i] for i in 1:(length(components)-1)))

    fun(z) = eos(model, z[1:end-1], z[end], T)
    df(z)  = ForwardDiff.gradient(fun,z)
    d2f(z) = ForwardDiff.hessian(fun,z)
    X = deepcopy(x)
    Y = deepcopy(y)
    d2f_l = d2f(append!(X,v_l))
    d2f_v = d2f(append!(Y,v_v))
    for i in 1:(length(components))
        J[i,1] =  log(10)*v_l*d2f_l[i,end]/R̄/model.params.epsilon[components[i]]
        J[i,2] = -log(10)*v_v*d2f_v[i,end]/R̄/model.params.epsilon[components[i]]
    end
    J[end,1] = log(10)*v_l*d2f_l[end,end]*model.params.sigma[components[1]]^3*N_A/R̄/model.params.epsilon[components[1]]
    J[end,2] =-log(10)*v_v*d2f_v[end,end]*model.params.sigma[components[1]]^3*N_A/R̄/model.params.epsilon[components[1]]

    for j in 1:(length(components)-1)
        J[end,j+2] = (d2f_v[end,end-1]-d2f_v[end,j])*model.params.sigma[components[1]]^3*N_A/R̄/model.params.epsilon[components[1]]
    end

    for i in 1:(length(components))
        for j in 1:(length(components)-1)
                J[i,j+2]= -(d2f_v[i,j]-d2f_v[i,end-1])/R̄/model.params.epsilon[components[i]]
        end
    end
end


text:
module Misc
using Gumbo, Cascadia

greet() = print("Hello World!")

include("grm.jl")
include("xys.jl")
include("ping-fan-de-shi-jie.jl")
include("muslin-de-zangli.jl")

end # module


text:
using CoupledConditionalSMC
using RNGPool
using Test

function testapproxequal(a::Vector{Float64}, b::Vector{Float64}, tol::Float64,
  verbose::Bool)
  v::Float64 = maximum(abs.(a-b))
  verbose && (println(a) ; println(b))
  verbose && println("$v < $tol ?")
  @test v < tol
end

function testapproxequal(a::Float64, b::Float64, tol::Float64, verbose::Bool)
  v::Float64 = abs(a-b)
  verbose && println("$v < $tol ?")
  @test v < tol
end

include("ffkJLS_test.jl")
include("ffk_test.jl")
include("estimate_test.jl")
include("couplingTime_test.jl")


text:
# This file was generated, do not modify it. # hide
pdf.(yhat, "virginica")

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
"""
Tools for the analysis of noise data, including computation of the noise
autocorrelation and power-spectral density.

Power spectra require a window function. We suggest `NoiseAnalysis.bartlett` as
a default. Package DSP.jl has many others in module `DSP.Windows`.
"""
module NoiseAnalysis

export compute_autocorr, compute_psd, NoiseResult

using ARMA
using FFTW
using HDF5
using Printf
using Statistics

struct NoiseResult
    autocorr::Vector{Float64}
    psd::Vector{Float64}

    samplesused::Int
    freqstep::Float64
    datasource::String
    model::ARMAModel
end

"""    hdf5save(h5file::HDF5.HDF5File, channum::Integer, nr::NoiseResult)"""
function hdf5save(h5file::HDF5.HDF5File, channum::Integer, nr::NoiseResult)
    chanstring = string(channum)
    g1 = g_create(h5file, chanstring)
    g = g_create(g1, "noise")
    hdf5save(g, nr)
end

"""    hdf5save(g::HDF5.HDF5Group, nr::NoiseResult)"""
function hdf5save(g::HDF5.HDF5Group, nr::NoiseResult)
    g["samplesused"] = nr.samplesused
    g["freqstep"] = nr.freqstep
    g["autocorr"] = nr.autocorr
    g["powerspectrum"] = nr.psd
    g["source"] = nr.datasource
    ARMA.hdf5save(g, nr.model)
end

"""
    model = hdf5load(input, [channumber])

Load and return a `NoiseResult` object from `input`. The argument `input` can be
either the "noise" HDF5 group, or the name of an HDF5 file along
(group "5/noise" would need to be at the root of the file).
"""
function hdf5load(hdf5filename::AbstractString, channum::Integer)
    chanstring = string(channum)
    h5open(hdf5filename, "r") do h5file
        if !(chanstring in names(h5file))
            message = @sprintf("cannot find a NoiseResult for chan %d in HDF5 file '%s'",
                channum, hdf5filename)
            error(message)
        end
        g = h5file[chanstring]["noise"]
        return hdf5load(g)
    end
end

function hdf5load(g::HDF5.DataFile)
    acorr = g["autocorr"][:]
    psd = g["powerspectrum"][:]
    sampused = read(g["samplesused"])
    freqstep = read(g["freqstep"])
    source = read(g["source"])
    model = ARMA.hdf5load(g["ARMAModel"])
    NoiseResult(acorr, psd, sampused, freqstep, source, model)
end

"""
    round_up_dft_length(n)

A length at least as large as `n` convenient for a DFT: (1, 3, or 5) x any power
of 2. This ensures that we never do DFTs on an especially bad vector size.
"""
function round_up_dft_length(n::Integer)
    pow2 = 2^ceil(log2(n))
    if n > 0.75*pow2
        return round(Int, pow2)
    elseif n > 0.625*pow2
        return round(Int, 0.75*pow2)
    else
        return round(Int, 0.625*pow2)
    end
end

"""
    compute_autocorr(data, nlags; chunk_multiple, max_exc)

Compute the autocorrelation function of the data, for lags `0:nlags-1`.

It is assumed that the length of the data vector is much greater than the
number of lags desired and that the data values are continuous (no gaps).

Group the data into segments of lengths `nlags*chunk_multiple` before computing
autocorrelation, so as to have significant number of samples at the maximum lag.
Typically any `chunk_multiple` of at least 3 should be fine (default = 7). In
any chunk, if the extreme value deviates by more than `max_exc` from the mean
value in either direction, the chunk will not be used. By default, `max_exc` is
large (1e99), but it is smart to set it to a more reasonable value for your
application.
"""
compute_autocorr(data::AbstractMatrix, nlags::Integer; kwargs...) =
    compute_autocorr(vcat(data...), nlags; kwargs...)

function compute_autocorr(data::AbstractVector, nlags::Integer;
        chunk_multiple=7, max_exc=1e99)
    m = nlags*chunk_multiple
    nseg = div(length(data), m)
    # Pad the raw data with at least nlags zeros, but also additional ones to ensure
    # that the DFT isn't taken on an unfortunately long and inefficient size.
    m_padded = round_up_dft_length(m+nlags)
    padded_data = zeros(Float64, m_padded)
    ac = zeros(Float64, nlags)

    # This loop uses all data except the last length(data) % m values.
    if nseg<1
        error("nseg=$nseg, should be >=1")
    end
    seg_skipped = 0
    for i=1:nseg
        seg = float(data[(i-1)*m+1:i*m])
        padded_data[1:m] = seg .- mean(seg)
        if maximum(abs.(padded_data[1:m])) > max_exc
            seg_skipped += 1
            continue
        end
        r = rfft(padded_data)
        ac += irfft(abs2.(r), m_padded)[1:nlags]
    end
    seg_used = nseg-seg_skipped
    if seg_used == 0
        error("all segments excluded by max excursion")
    end
    ac /= seg_used
    ac ./ (m:-1:m-nlags+1)
end

# Power spectra and window functions:
"""
    bartlett(n)

Bartlett window of length `n`.
"""
function bartlett(n::Integer)
    [(1 - abs(k*2/(n-1) - 1)) for k=0:(n-1)]
end

"""
    welch(n)

Welch window of length `n`.
"""
function welch(n::Integer)
    [(1 - ((k*2/(n-1)) - 1)^2) for k=0:(n-1)]
end

"""
    hann(n)

Hann window of length `n`.
"""
function hann(n::Integer)
    [1-cos(2π*k/(n-1)) for k=0:(n-1)]
end

# See package DSP.jl with module DSP.Windows for more window functions.

"""
    compute_psd(data, nfreq, dt)

Compute the Power-Spectral Density of `data`, which was sampled at
equal time steps of size `dt`, at `nfreq` frequencies equal to
`fsamp = linspace(0, stop=0.5/dt, length=nfreq)`. Generally, you want `nfreq` to be
1 more than a power of 2, and you need `nfreq ≤ div(1+length(data), 2)`
(ideally, much less than).

The `data` are assumed to be a continuous sequence of noise samples.
Use overlapping segments of the exact needed length (`2(nfreq-1)`), offset by
approximately half their length.
"""
compute_psd(data::AbstractArray, nfreq::Integer, dt::Real; kwargs...) =
    compute_psd(vcat(data...), nfreq, dt; kwargs...)

function compute_psd(data::AbstractVector, nfreq::Integer, dt::Real; max_exc=1e99)

    nsamp = 2(nfreq-1)
    if length(data) < nsamp
        error("data must be at least 2(nfreq-1) in length.")
    end

    # If you have any appreciable power at DC and wish to understand the low-freq
    # power spectrum, you want the hann window or no window, because only these
    # have zero leakage for bins 3:end (zero leakage into bin 2, also, if using
    # no window at all).
    window = hann(nsamp)
    window = window / sqrt(sum(window.^2)) # proper normalization

    # Break data up into half-overlapping segments of the right separation
    nseg = ceil(Int, length(data)/nsamp)
    seg_step = 1
    if nseg>1
        seg_step = Int(floor((length(data)+1-nsamp)/(nseg-1)))
    end

    r = zeros(Float64, nfreq)
    for i=1:nseg
        idx0 = (i-1)*seg_step+1
        dseg = data[idx0:idx0+nsamp-1]
        # Subtract this segment's mean (DC value) to eliminate DC leakage into 3rd bin.
        seg = window .* (dseg .- mean(dseg))
        r += abs2.(rfft(seg))
    end
    r * 2dt / nseg
end

psd_freq(nfreq::Integer, dt::Real) = range(0, stop=0.5/dt, length=nfreq)

end # module


text:
include(joinpath("src", "module.jl"))

config_list = [
    Dict(
        0 => [(1, 1), (1, 2)],
        1 => [(10, 10), (10, 9)]
    ),
    Dict(
        0 => [(5, 5), (5, 6)],
        1 => [(6, 5), (6, 6)]
    )
]

AxelFraud.run_batch(
    config_list=config_list,
    batch_name="competitive",
    grid_dims=(10, 10),
    levels=2,
    steps=3000,
    replicates=300,
    when=10,
    write=true
)

AxelFraud.run_random(
    batch_name="test",
    grid_dims=(10, 10),
    stubborn_count=10,
    steps=100,
    replicates=2,
    when=10,
    rnd_seed=0,
    write=true
)


text:
# used in beamsearch
import Base.(<)
(<)(x::Var, y::Var) = (<)(x.data[1], y.data[1])

type StructuredFeedForward
    feedforward::FeedForward
    logits::Var
    beamid::Int
end

function initmodel!(parser::DepParser, model::Type{StructuredFeedForward})
    parser.model = StructuredFeedForward(parser.model, Var([0f0]), 1)
end

@compat function (m::StructuredFeedForward){T}(ss::AbstractVector{State{T}})
    m.logits = m.feedforward(ss)
    m.beamid = 1
end

@compat function (m::StructuredFeedForward){T}(s::State{T}, act::Int)
    s.score + m.logits[[act], [m.beamid]]
end

typealias Chart Vector{Vector{State{Labeled}}}

function earlystep(golds::Chart, preds::Chart, beamsize::Int; alive=length(golds))
    zerov = Var([0f0])
    step = 0
    for i = 1:min(length(golds), alive)
        step = i
        goldscore = golds[i][1].score
        predscore = beamsize > length(preds[i]) ? zerov : preds[i][beamsize].score
        goldscore < predscore && break
    end
    step
end

function earlyupdate!(opt::SGD, golds::Chart, preds::Chart, beamsize::Int)
    k = earlystep(golds, preds, beamsize)
    k == 1 && return 0.0
    predscores = [s.score for s in preds[k][1:min(length(preds[k]), beamsize)]]
    mat = concat(1, golds[k][1].score, predscores...)
    mat = reshape(mat, length(mat), 1)
    update!(opt, [1], mat)
end

function earlyupdate!(opt::SGD, golds::Vector{Chart}, preds::Vector{Chart}, alive::Vector{Int}, beamsize::Int, batchsize::Int)
    logits = Array{Var}(0)
    for (i, (gold, pred)) in enumerate(zip(golds, preds))
        k = earlystep(gold, pred, beamsize, alive=alive[i])
        k == 1 && continue
        predscores = [s.score for s in pred[k][1:min(length(pred[k]), beamsize)]]
        while length(predscores) < beamsize
            push!(predscores, Param([0f0]))
        end
        mat = concat(1, gold[k][1].score, predscores...)
        mat = reshape(mat, length(mat), 1)
        push!(logits, mat)
    end
    update!(opt, ones(Int, length(logits)), concat(2, logits))
end

import TransitionParser.beamsearch
function beamsearch{T}(parser::DepParser, s::T, beamsize::Int, expand::Function)
    chart = Vector{T}[]
    push!(chart, [s])
    lessthan{T}(x::T, y::T) = x.score > y.score

    k = 1
    while k <= length(chart)
        states = chart[k]
        length(states) > beamsize && sort!(states, lt=lessthan)
        beam = 1:min(beamsize, length(states))
        parser.model(states[beam])
        for i in beam
            parser.model.beamid = i
            for s in expand(states[i])
                while s.step > length(chart)
                    push!(chart, T[])
                end
                push!(chart[s.step], s)
            end
        end
        k += 1
    end
    sort!(chart[end], lt=lessthan)
    chart
end

function beamsearch{T}(parser::DepParser, ss::AbstractVector{T}, beamsize::Int, expand::Function)
    lessthan{T}(x::T, y::T) = x.score > y.score
    charts = Vector{Vector{T}}[]
    for (i, s) in enumerate(ss)
        push!(charts, [])
        push!(charts[i], [s])
    end

    alive = zeros(Int, length(charts))
    k = 1
    while any(v -> v == k-1, alive)
        offset = 0
        parser.model(vcat([chart[k] for chart in charts]...))
        for cid = 1:length(charts)
            cid > 1 && ( offset += length(charts[cid-1][k]) )
            chart = charts[cid]
            states = chart[k]
            if isfinal(states[1])
                push!(chart, states)
            else
                alive[cid] = k
                nexts = map(1:length(states)) do sid
                         parser.model.beamid = offset + sid
                         expand(states[sid])
                end
                nexts = vcat(nexts...)
                sort!(nexts, lt=lessthan)
                push!(chart, nexts[1:min(length(nexts), beamsize)])
            end
        end
        k += 1
    end
    charts, alive
end

function train!{T}(::Type{StructuredFeedForward}, parser::DepParser{T}, trainsents::Doc,
    testsents::Doc=Vector{Token}[]; embed="", beamsize=31, batchsize=32, iter=20,
    progbar=true, opt=SGD(0.0001, 0.9), evaliter=100, outfile="parser.dat")
    info("WILL RUN $iter ITERATIONS")

    saver = ModelSaver(outfile)
    initmodel!(parser, StructuredFeedForward)
    labelsize = targetsize(parser.model.feedforward)

    info("OPTIMIZER: ", typeof(opt))
    info("LABEL SIZE: $labelsize")
    info("BEAM SIZE: $beamsize")
    info("#TRAINSENTS: $(length(trainsents))")

    for i = 1:iter
        info("ITER $i TRAINING")
        batch = sub(trainsents, rand(1:length(trainsents), batchsize))
        loss = 0.0
        # progbar && ( p = Progress(batchsize, 1, "", 50) )
        # res = map(batch) do s
        #     progbar && next!(p)
        #     s = State(s, parser)
        #     golds = beamsearch(parser, s, 1, expandgold)
        #     preds = beamsearch(parser, s, beamsize, expandpred)
        #     loss += earlyupdate!(opt, golds, preds, beamsize)
        #     preds[end][1]
        # end
        ss = map(s -> State(s, parser), batch)
        golds, _ = beamsearch(parser, ss, 1 ,expandgold)
        preds, alive = beamsearch(parser, ss, beamsize, expandpred)
        loss = earlyupdate!(opt, golds, preds, alive, beamsize, batchsize)
        res = map(chart -> chart[end][1], preds)

        info(now())
        info("LOSS: ", loss)
        evaluate(parser, res)

        if i % evaliter == 0
            info("SAVING PARSER TO $outfile")
            open(io -> serialize(io, parser), outfile * "." * string(i), "w")
            # println()
            # info("**ITER $i TESTING**")
            # res = decode(StructuredFeedForward, parser, testsents,
            #                  beamsize=beamsize, batchsize=batchsize)
            # uas, las = evaluate(parser, res)
            # saver(parser, uas)
        end
        println()
    end
end

function decode{T}(::Type{StructuredFeedForward}, parser::DepParser{T}, sents::Doc;
    beamsize=16, batchsize=32, progbar=true)
    res = []
    batches = 1:batchsize:length(sents)
    progbar && ( p = Progress(length(batches), 1, "", 50) )
    for i in batches
        progbar && next!(p)
        batch = sub(sents, i:min(i+batchsize-1, length(sents)))
        ss = map(s -> State(s, parser), batch)
        preds, _ = beamsearch(parser, ss, beamsize, expandpred)
        append!(res, map(chart -> chart[end][1], preds))
        preds = 0 # gc
    end
    res
    # map(sents) do s
    #     progbar && next!(p)
    #     s = State(s, parser)
    #     beamsearch(parser, s, beamsize, expandpred)[end][1]
    # end
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
using OrdinaryDiffEq, Test, Statistics
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear

@testset "Scalar Discrete Problem" begin
  prob = DiscreteProblem(0.5,(0.0,1.0))
  sol =solve(prob,FunctionMap())
  @test sol[1] == sol[end]
  @test sol(0.5:0.1:0.7) == [sol[1],sol[1],sol[1]]

  sol =solve(prob_ode_linear,FunctionMap())
  @test sol[end] ≈ .505
  sol =solve(prob_ode_linear,FunctionMap(scale_by_time=true),dt=1/4)
  sol2 =solve(prob_ode_linear,Euler(),dt=1/4)
  @test sol[end] == sol2[end]
  @test sol(0.53) != sol2(0.53)
end

@testset "Array Discrete Problem" begin
  prob2 = DiscreteProblem(rand(4,2),(0.0,1.0))
  sol =solve(prob2,FunctionMap())
  @test sol[1] == sol[end]
  @test sol(0.5) == sol[1]

  @test_nowarn sol =solve(prob_ode_2Dlinear,FunctionMap())
  @test_nowarn sol2 =solve(prob_ode_2Dlinear,Euler(),dt=1)
  sol =solve(prob_ode_2Dlinear,FunctionMap(scale_by_time=true),dt=1/4)
  sol2 =solve(prob_ode_2Dlinear,Euler(),dt=1/4)
  @test sol[end] == sol2[end]
  @test sol(0.35) != sol2(0.53)
end

@testset "Discrete Problem (Henon Map)" begin
  function henon_map!(u_next, u, _p, t)
    u_next[1] = 1 + u[2] - 1.4 * u[1]^2
    u_next[2] = 0.3 * u[1]
  end

  prob = DiscreteProblem(henon_map!, [0.5, 0.5], (0, 100)) # Integer time
  sol = solve(prob, FunctionMap())
  @test eltype(sol.t) <: Integer
end

@testset "DiscreteProblem in-place time definition" begin
  function f(un,p,t_np1)
    un + t_np1
  end

  prob2 = DiscreteProblem(f,0,(0,5))
  sol2 = solve(prob2,FunctionMap())

  @test sol2.u == [0,1,3,6,10,15]

  function f!(u_np1,un,p,t_np1)
    u_np1[1] = un[1] + 1
    u_np1[2] = un[2] + t_np1
  end

  prob3 = DiscreteProblem(f!,[0,0],(0,5))
  sol3 = solve(prob3,FunctionMap())

  @test sol3.u == [[0,0],[1,1],[2,3],[3,6],[4,10],[5,15]]
end


text:
@enum TeamPreviewArea begin
    POKEMON_ICON
    ITEM_ICON
    GENDER_ICON
    LEVEL_TEXT
    HP_TEXT
    PLAYER_NAME
end

const PokemonIconSheetRect = Rect((13, 1), (100, 132))
const TeamPreviewPlayerBoxRects = [
    # Team A
    Rect((125, 301), (784, 422), resolution = DesignResolution),
    # Team B
    Rect((125, 1201), (784, 422), resolution = DesignResolution),
    # Team A when selecting
    Rect((45, 656), (848, 650), resolution = DesignResolution),
    # Team B when selecting
    Rect((191, 103), (784, 422), resolution = DesignResolution),
]
const TeamPreviewBoxHeight = 120
const TeamPreviewSelectingBoxHeight = 144

# (top, left), (height, width), (top_selecting, left_selecting), (height_selecting, width_selecting)
_get_team_preview_rect_starting(area_type, is_small_box) = if area_type ≡ POKEMON_ICON
    if is_small_box
        (64, 8), (100, 132)
    else
        (14, 2), (100, 132)
    end
elseif area_type ≡ ITEM_ICON
    if is_small_box
        (117, 355), (48, 48)
    else
        (11, 596), (48, 48)
    end
elseif area_type ≡ GENDER_ICON
    if is_small_box
        (116, 268), (50, 50)
    else
        (10, 419), (50, 50)
    end
elseif area_type ≡ LEVEL_TEXT
    if is_small_box
        (119, 143), (46, 120)
    else
        (73, 396), (56, 168)
    end
elseif area_type ≡ HP_TEXT  # Only available for large player box
    (78, 134), (40, 200)
else  # NAME
    (1, 7), (56, 408)
end

function get_area(img, i, player, area_type::TeamPreviewArea)
    box = img[TeamPreviewPlayerBoxRects[player]]
    is_small_box = player ≠ 3
    (top, left), size = _get_team_preview_rect_starting(area_type, is_small_box)
    box_height = is_small_box ? TeamPreviewBoxHeight : TeamPreviewSelectingBoxHeight
    rect = Rect((top + (i - 1) * box_height, left), size, resolution = DesignResolution)
    box[rect]
end

struct TeamPreview <: AbstractPokemonScene
    frame_time::Float64
    team_a::Vector{String}
    uncertain::Int
    team_b::Vector{String}
    genders_a::Vector{Gender}
    genders_b::Vector{Gender}
    levels_a::Vector{Int64}
    levels_b::Vector{Int64}
    items_a::Vector{String}
    hps_a::Vector{Int}
    name_a::String
    name_b::String
end

struct TeamPreviewSelecting <: AbstractPokemonScene
    # Parse this will return a `TeamPreview` struct.
end

feature_image_name(::Type{TeamPreview}, ::PokemonBattle) = "team-preview"
feature_image_name(::Type{TeamPreviewSelecting}, ::PokemonBattle) = "team-preview-selecting"

function parse_pokemon_icon(img, i, player; preprocess=(player == 3))
    icons = pokemon_icons()
    icon = get_area(img, i, player, POKEMON_ICON)
    uncertain = if preprocess
        h, w = size(icon)
        rect = icon.rect
        rh, rw = rect.height, rect.width
        icon = floodfill(icon, (1, 1), 1, 0.1)
        uncertain = Gray(icon[1, w]) < 0.02
        floodfill!(icon, (100 * h ÷ rh, 10 * w ÷ rw), 1, 0.1)
        floodfill!(icon, (1, w), 1, 0.1)
        uncertain
    else
        false
    end
    table_search(
        icon, icons,
        rect = PokemonIconSheetRect
    ), uncertain
end

function parse_item(img, i, player; preprocess=true)
    icons = item_icons()
    icon = get_area(img, i, player, ITEM_ICON)
    if preprocess
        icon = floodfill(icon, (1, 1), 1, 0.1)
    end
    table_search(
        icon, icons
    )
end

function parse_gender(img, i, player, selected=false)
    icons = gender_icons()
    icon = get_area(img, i, player, GENDER_ICON)
    if selected
        icon = floodfill(icon, (1, 1), 1, 0.1)
    end
    r = table_search(
        icon, icons
    )
    r == "null" ? GENDER_NULL : r == "male" ? GENDER_MALE : GENDER_FEMALE
end

function _parse_number(img, i, player, ctx, selected, area_type)
    text_area = get_area(img, i, player, area_type)
    if selected
        text_area = complement.(text_area)
    end
    ocr(Int, text_area, ctx; language = "eng")
end

parse_level(img, i, player, ctx, selected=false) =
    _parse_number(img, i, player, ctx, selected, LEVEL_TEXT)

parse_hp(img, i, ctx, selected=false) =
    _parse_number(img, i, 3, ctx, selected, HP_TEXT)

function parse_name(img, player, ctx; threshold = 0.35f0)
    player == 3 && return ""
    text_area = get_area(img, 1, player, PLAYER_NAME)
    text_area = prepare_text_for_ocr(text_area; threshold = threshold)
    ocr_multiple_lang(text_area, ctx)[1]
end

# player == 1,2,3,4 for A,B,A-Selecting,B-Selecting
function parse_player(img, player, ctx)
    config = ctx.config

    parse_type = config.source.parse_type
    is_player_a = player == 1 || player == 3
    if is_player_a && parse_type ≢ PARSE_BOTH_PLAYERS
        return String[], 0, Gender[], Int64[], String[], Int64[], ""
    end

    team = Vector{String}(undef, 6)
    uncertain = Threads.Atomic{Int64}(0)
    genders = parse_type > PARSE_MINIMAL ? Vector{Gender}(undef, 6) : Gender[]
    items = is_player_a ? Vector{String}(undef, 6) : String[]
    name = ""

    Threads.@threads for i in 1:6
        poke, u = parse_pokemon_icon(img, i, player)
        if u
            uncertain[] = i
        end
        team[i] = poke
        if parse_type > PARSE_MINIMAL
            genders[i] = parse_gender(img, i, player, u)
        end
        if is_player_a
            items[i] = parse_item(img, i, player)
        end
    end
    levels = if parse_type > PARSE_MINIMAL
        Int64[
            parse_level(img, i, player, ctx, uncertain[] ≡ i)
            for i in 1:6
        ]
    else
        Int64[]
    end
    hps = if player == 3
        Int64[
            parse_hp(img, i, ctx, uncertain[] ≡ i)
            for i in 1:6
        ]
    else
        Int64[]
    end
    name = parse_name(img, player, ctx)
    return team, uncertain[], genders, levels, items, hps, name
end

function _parse_scene(::Type{TeamPreview}, frame::VsFrame, ctx::PokemonContext, a = 1, b = 2)
    data = get_current_context!(ctx).data
    if !ctx.data.force && !isnothing(data.team_preview)
        parse_type = ctx.config.source.parse_type
        if parse_type ≢ PARSE_BOTH_PLAYERS || a == 1 || data.uncertain_poke == 0
            # Already parsed everything we need
            return nothing
        end
    end
    img = image(frame)
    team_a, uncertain, genders_a, levels_a, items_a, hps_a, name_a = parse_player(img, a, ctx)
    team_b, _, genders_b, levels_b, _, _, name_b = parse_player(img, b, ctx)

    data.team_preview = true
    TeamPreview(
        frame.time,
        team_a, uncertain, team_b,
        genders_a, genders_b,
        levels_a, levels_b,
        items_a, hps_a,
        name_a, name_b
    )
end

_parse_scene(::Type{TeamPreviewSelecting}, frame::VsFrame, ctx::PokemonContext) =
    _parse_scene(TeamPreview, frame, ctx, 3, 4)

function _vs_update!(ctx::PokemonContext, scene::TeamPreview)
    parsing = get_current_context!(ctx)
    battle = parsing.battle
    parsed_battle = parsing.parsed_battle
    parse_type = ctx.config.source.parse_type

    data = parsing.data
    if parse_type ≡ PARSE_BOTH_PLAYERS
        # TODO: set uncertain
        if isnothing(data.uncertain_poke) && scene.uncertain > 0
            data.uncertain_poke = scene.uncertain
        elseif data.uncertain_poke != scene.uncertain
            data.uncertain_poke = 0
        end
        if scene.name_a != ""
            battle.player.name = scene.name_a
        end
        for (i, value) in enumerate(scene.team_a)
            i == scene.uncertain && continue
            update_team!(parsed_battle, true, i, :id, value)
        end
        for (i, value) in enumerate(scene.genders_a)
            update_team!(parsed_battle, true, i, :gender, value)
        end
        for (i, value) in enumerate(scene.levels_a)
            update_team!(parsed_battle, true, i, :level, value)
        end
        for (i, value) in enumerate(scene.items_a)
            update_team!(parsed_battle, true, i, :item, value)
        end
        for (i, value) in enumerate(scene.hps_a)
            update_team!(parsed_battle, true, i, :hp, value)
        end
    end
    if scene.name_b != ""
        battle.opponent.name = scene.name_b
    end
    battle.opponent_team = [PokemonID(x) for x in scene.team_b]
    for (i, value) in enumerate(scene.team_b)
        update_team!(parsed_battle, false, i, :id, value)
    end
    for (i, value) in enumerate(scene.genders_b)
        update_team!(parsed_battle, false, i, :gender, value)
    end
    for (i, value) in enumerate(scene.levels_b)
        update_team!(parsed_battle, false, i, :level, value)
    end
    ctx
end


text:
@cache mutable struct ImplicitEMCache{uType,rateType,N,noiseRateType,dWType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  gtmp::noiseRateType
  gtmp2::rateType
  dW_cache::dWType
  nlsolver::N
end

function alg_cache(alg::ImplicitEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  gtmp = zero(noise_rate_prototype)
  if is_diagonal_noise(prob)
    gtmp2 = gtmp
    dW_cache = nothing
  else
    gtmp2 = zero(rate_prototype)
    dW_cache = zero(ΔW)
  end

  γ, c = alg.theta,zero(t)
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  fsalfirst = zero(rate_prototype)
  ImplicitEMCache(u,uprev,fsalfirst,gtmp,gtmp2,dW_cache,nlsolver)
end

mutable struct ImplicitEMConstantCache{N} <: StochasticDiffEqConstantCache
  nlsolver::N
end

function alg_cache(alg::ImplicitEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  γ, c = alg.theta,zero(t)
  J, W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  ImplicitEMConstantCache(nlsolver)
end

@cache mutable struct ImplicitEulerHeunCache{uType,rateType,N,noiseRateType,dWType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  gtmp::noiseRateType
  gtmp2::rateType
  gtmp3::noiseRateType
  nlsolver::N
  dW_cache::dWType
end

u_cache(c::ImplicitEulerHeunCache)    = (c.uprev2,c.nlsolver.z,c.nlsolver.dz)
du_cache(c::ImplicitEulerHeunCache)   = (c.nlsolver.k,c.fsalfirst)

function alg_cache(alg::ImplicitEulerHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  gtmp = zero(noise_rate_prototype)
  gtmp2 = zero(rate_prototype)

  if is_diagonal_noise(prob)
      gtmp3 = gtmp2
      dW_cache = nothing
  else
      gtmp3 = zero(noise_rate_prototype)
      dW_cache = zero(ΔW)
  end

  γ, c = alg.theta,zero(t)
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  fsalfirst = zero(rate_prototype)
  ImplicitEulerHeunCache(u,uprev,fsalfirst,gtmp,gtmp2,gtmp3,nlsolver,dW_cache)
end

mutable struct ImplicitEulerHeunConstantCache{N} <: StochasticDiffEqConstantCache
  nlsolver::N
end

function alg_cache(alg::ImplicitEulerHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  γ, c = alg.theta,zero(t)
  J, W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  ImplicitEulerHeunConstantCache(nlsolver)
end

@cache mutable struct ImplicitRKMilCache{uType,rateType,N,noiseRateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  gtmp::noiseRateType
  gtmp2::noiseRateType
  gtmp3::noiseRateType
  nlsolver::N
end

function alg_cache(alg::ImplicitRKMil,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  gtmp = zero(noise_rate_prototype)
  gtmp2 = zero(rate_prototype)
  gtmp3 = zero(rate_prototype)

  γ, c = alg.theta,zero(t)
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  fsalfirst = zero(rate_prototype)
  ImplicitRKMilCache(u,uprev,fsalfirst,gtmp,gtmp2,gtmp3,nlsolver)
end

mutable struct ImplicitRKMilConstantCache{N} <: StochasticDiffEqConstantCache
  nlsolver::N
end

function alg_cache(alg::ImplicitRKMil,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  γ, c = alg.theta,zero(t)
  J, W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  ImplicitRKMilConstantCache(nlsolver)
end


text:
cufftGetVersion() = ccall((:cufftGetVersion,libcufft), Cint, ())

function cufftGetProperty(property::CUDAapi.libraryPropertyType)
  value_ref = Ref{Cint}()
  @check ccall((:cufftGetProperty, libcufft),
               cufftStatus_t,
               (Cint, Ptr{Cint}),
               property, value_ref)
  value_ref[]
end


text:
module InstaRound

# Import Base round
import Base: round

include("custom.jl")
include("utils.jl")

####### Globals Exports #######
export IGRound

end


text:
# Cyber Test Problem
# POMDP, Robust POMDP, and Robust ρPOMDP formulations

struct CyberTestIPOMDP <: IPOMDP{Vector{Int}, Vector{Int}, Vector{Int}}
    discount::Float64
    inforeward::Vector{Vector{Float64}}
    b0::Vector{Float64}
end

# Discount
const discount_cyber_test = 0.75

# State space
const nmoe_test = 3 # number of measures of effectiveness (MOE)
const nlevels_test = 3 # number of levels at which the MOE is measured
const lmin_test = 1 # min level
const lmax_test = 3 # max level
const nS_test = nlevels_test ^ nmoe_test # number of states
S = Vector{Vector{Int}}(nS_test) # states, as [moe_1_level, ...]
ind = 1
for i = 1:nlevels_test, j = 1:nlevels_test, k = 1:nlevels_test
  S[ind] = [i,j,k]
  ind += 1
end
const states_cyber_test = copy(S)

# Action space
const na_test = 2 # number of assessment assets
const nA_test = 6 # total number of actions
A = Vector{Vector{Int}}() # actions, as [asset_1_assignment, ...]
for i = 1:nmoe_test, j = 1:nmoe_test
  (i != j) && (push!(A, [i,j]))
end
const actions_cyber_test = copy(A)

# Observation space
const nZ_test = nlevels_test ^ na_test
Z = Vector{Vector{Int}}(nZ_test)
ind = 1
for i = 1:nlevels_test, j = 1:nlevels_test
  Z[ind] = [i,j]
  ind += 1
end
const observations_cyber_test = copy(Z)

# Belief reward (PWLC vectors: simple)
br = Vector{Vector{Float64}}(nS_test)
br[1] = vcat(1.0, fill(-1/(nS_test-1), nS_test - 1))
br[nS_test] = vcat(fill(-1/(nS_test-1), nS_test - 1), 1.0)
for i = 2:(nS_test-1)
  br[i] = vcat(fill(-1/(nS_test-1), i - 1), 1.0, fill(-1/(nS_test-1), nS_test - i))
end
const inforeward_cyber_test = copy(br)

# Initial state and belief
e11 = zeros(nS_test); e11[1] = 1.0;
const e1_test = copy(e11)        # intial state
const b0_cyber_test = copy(e11)  # initial belief

CyberTestIPOMDP(alphas::Vector{Vector{Float64}}) = CyberTestIPOMDP(discount_cyber_test, alphas, b0_cyber_test)
CyberTestIPOMDP(disc::Float64) = CyberTestIPOMDP(disc, inforeward_cyber_test, b0_cyber_test)
CyberTestIPOMDP() = CyberTestIPOMDP(discount_cyber_test, inforeward_cyber_test, b0_cyber_test)

states(::CyberTestIPOMDP) = states_cyber_test
actions(::CyberTestIPOMDP) = actions_cyber_test
observations(::CyberTestIPOMDP) = observations_cyber_test

n_states(::CyberTestIPOMDP) = nS_test
n_actions(::CyberTestIPOMDP) = nA_test
n_observations(::CyberTestIPOMDP) = nZ_test

state_index(::CyberTestIPOMDP, s::Vector{Int}) = find(x -> x == s, states_cyber_test)[1]
action_index(::CyberTestIPOMDP, a::Vector{Int}) = find(x -> x == a, actions_cyber_test)[1]
observation_index(::CyberTestIPOMDP, z::Vector{Int}) = find(x -> x == z, observations_cyber_test)[1]
obs_index(prob::CyberTestIPOMDP, z::Vector{Int}) = observation_index(prob, z)

initial_state_distribution(::CyberTestIPOMDP) = SparseCat(states_cyber_test, e1_test)
initial_belief(prob::CyberTestIPOMDP) = prob.b0
initial_belief_distribution(prob::CyberTestIPOMDP) = SparseCat(states_cyber_test, initial_belief(prob))

# Transitions
const pϵ_cyber_test = 1e-6 # near-saturation bound
const delt_test = 0.01 # imprecision in transition distribution
const pd_test = 0.15 # prob of decrease
const ps_test = 0.55 # prob of stay
const pi_test = 0.3 # prob of increase
const psd_test = ps_test + pd_test # prob of decline or stay at lower border
const psi_test = ps_test + pi_test # prob of stay or improve at upper border

# calculate trasnsition matrix
function calcTArray_test(S::Vector{Vector{Int}}, A::Vector{Vector{Int}}, par, delt)
  ns = length(S)
  na = length(A)
  T = zeros(ns, na, ns)
  Tl = zeros(ns, na, ns)
  Tu = zeros(ns, na, ns)
  for (si,s) in enumerate(S), (ji,j) in enumerate(S)
      T[si,:,ji] = psj(s,j)
      tl = psj(s,j,par,-delt)
      tu = psj(s,j,par,delt)
      # tl = max(tl, 0.0 + pϵ_cyber)
      # tu = min(tu, 1.0 - pϵ_cyber)
      (tu < tl) && (tu = tl + pϵ_cyber_test / 2)
      Tl[si,:,ji] = tl
      Tu[si,:,ji] = tu
  end
  for (si,s) in enumerate(S), (ai,a) in enumerate(A)
    T[si,ai,:] = T[si,ai,:] ./ sum(T[si,ai,:])
  end
  T, Tl, Tu
end

# Nominal transition function array tarray_cyber[s,a,s'] = Pr(s' | s ,a)
const tarray_cyber_test, tlarray_cyber_test, tuarray_cyber_test = calcTArray_test(states_cyber_test, actions_cyber_test, (pd_test,ps_test,pi_test), delt_test)

# Nominal transition function distributions
const tdist_cyber_test = calcTDist(nS_test, nA_test, tarray_cyber_test)

# Nominal transitions
transition(prob::CyberTestIPOMDP, s::Vector{Int}, a::Vector{Int}) = tdist_cyber_test[state_index(prob, s), action_index(prob, a)]
transition(prob::CyberTestIPOMDP, s::Vector{Int}, a::Vector{Int}, sp::Vector{Int}) = tarray_cyber_test[state_index(prob, s), action_index(prob, a), state_index(prob, sp)]

# Observation function
const delo1_test = pϵ_cyber_test # sensor 1 immprecision
const delo2_test = 0.4 # sensor 2 imprecision
const sensor_imprecision_test = [delo1_test, delo2_test]
const acc1_test = 0.8 - delo1_test # likelihood of correct observation error for sensor 1 (off-nominal)
const acc2_test = 0.9 - delo2_test # likelihood of correct observation error for sensor 2 (off-nominal)
const sensor_accuracy_test = [acc1_test, acc2_test]

function calcOArray_test(S::Vector{Vector{Int}}, A::Vector{Vector{Int}}, Z::Vector{Vector{Int}})
  O = zeros(nA,nS,nZ)
  Ol = zeros(nA,nS,nZ)
  Ou = zeros(nA,nS,nZ)
  for (spi,sp) in enumerate(S), (zi,z) in enumerate(Z), (ai,a) in enumerate(A)
    O[ai, spi, zi] = o(a, sp, z, sensor_accuracy_test)
    neg = o(a, sp, z, sensor_accuracy_test - sensor_imprecision_test)
    pos = o(a, sp, z, sensor_accuracy_test + sensor_imprecision_test)
    ol = min(neg, pos)
    ou = max(neg, pos)
    (ol < 0.0) && (ol = max(ol, 0.0 + pϵ_cyber_test))
    (ou > 1.0) && (ou = min(ou, 1.0 + pϵ_cyber_test))
    (ou < ol) && (ou = ol + pϵ_cyber_test / 2)
    Ol[ai, spi, zi] = ol
    Ou[ai, spi, zi] = ou
  end
  for (spi,sp) in enumerate(S), (ai,a) in enumerate(A)
    O[ai, spi,:] = O[ai,spi,:] ./ sum(O[ai,spi,:])
  end
  O, Ol, Ou
end

const oarray_cyber_test, olarray_cyber_test, ouarray_cyber_test = calcOArray_test(states_cyber_test, actions_cyber_test, observations_cyber_test)

# Nominal transition function distributions
const odist_cyber_test = calcODist(nS_test, nA_test, oarray_cyber_test)

# Nominal transitions
observation(prob::CyberTestIPOMDP, a::Vector{Int}, sp::Vector{Int}) = odist_cyber_test[action_index(prob, a), state_index(prob, sp)]
observation(prob::CyberTestIPOMDP, a::Vector{Int}, sp::Vector{Int}, z::Vector{Int}) = oarray_cyber_test[action_index(prob, a), state_index(prob, sp), observation_index(prob, z)]

function dynamics(prob::CyberTestIPOMDP)
    ns = n_states(prob)
    nz = n_observations(prob)
    na_test = n_actions(prob)
    tarr = states(prob)
    sarr = states(prob)
    aarr = actions(prob)
    zarr = observations(prob)
    p = Array{Float64}(ns,nz,ns,na_test)
    plower = Array{Float64}(ns,nz,ns,na_test)
    pupper = Array{Float64}(ns,nz,ns,na_test)
    for tind in 1:ns, zind in 1:nz, sind in 1:ns, aind in 1:na_test
        pt = transition(prob, sarr[sind], aarr[aind], tarr[tind])
        pz = observation(prob, aarr[aind], tarr[tind], zarr[zind])
        p[tind, zind, sind, aind] = pz * pt
    end
    p
end

function reward(prob::CyberTestIPOMDP, b::Vector{Float64}, a::Vector{Int})
    rmax = -Inf
    for α in prob.inforeward
        rmax = max(rmax, dot(α, b))
    end
    rmax
end
reward(prob::CyberTestIPOMDP, b::Vector{Float64}, s::Vector{Int}, a::Vector{Int}, sp::Vector{Int}) = reward(prob,b,a)

function rewardalpha(prob::CyberTestIPOMDP, b::Vector{Float64}, a::Vector{Int})
    rmax = -Inf
    ralpha = nothing
    for α in prob.inforeward
        rnext = dot(α, b)
        if rnext > rmax
            rmax = rnext
            ralpha = copy(α)
        end
    end
    ralpha
end

discount(p::CyberTestIPOMDP) = p.discount

# Simulation Functions
function initial_state(prob::CyberTestIPOMDP, rng::AbstractRNG)
    d = initial_state_distribution(prob)
    return rand(rng, d)
end

function generate_sor(prob::CyberTestIPOMDP, b, s, a, rng::AbstractRNG)
    sp = rand(rng, transition(prob, s, a))
    o = rand(rng, observation(prob, a, sp))
    r = reward(prob, b, s, a, sp)
    sp, o, r
end


text:
module Bases

import Iterators: groupby, imap
import ..Utils: Interval

export Basis, Basis1D, domain, deriv, order,
       BSplineBasis,
       NURBSBasis

macro assert_ex(condition, error)
    :($condition ? nothing : throw($error))
end

# Abstract Basis types
# ========================================================================

abstract Basis
abstract Basis1D <: Basis

type BasisFunction1D{B<:Basis1D}
    basis::B
    index::Int
    deriv::Int

    function BasisFunction1D(basis, index, deriv)
        @assert_ex(1 <= index <= length(basis), BoundsError())
        @assert_ex(0 <= deriv <= nderivs(basis),
                   ArgumentError("Differentiation order not supported"))
        new(basis, index, deriv)
    end

    BasisFunction1D(basis, index) = BasisFunction1D(basis, index, 0)
end

deriv(b::BasisFunction1D) = typeof(b)(b.basis, b.index, b.deriv+1)
deriv(b::BasisFunction1D, order) = typeof(b)(b.basis, b.index, b.deriv+order)

function Base.call{B<:Basis1D, T<:Real}(b::B, pt::T)
    rng = supported(b, pt)
    (squeeze(evaluate_raw(b, [pt], b.deriv, rng), 2), rng)
end

function Base.call{B<:BasisFunction1D, T<:Real}(b::B, pt::T)
    rng = supported(b.basis, pt)
    if b.index ∉ rng return 0.0 end
    evaluate_raw(b.basis, [pt], b.deriv, rng)[1 + b.index - rng.start, 1]
end

function Base.call{B<:Basis1D, T<:Real}(b::B, pts::Vector{T})
    tp = (Vector{Float64}, UnitRange{Int})
    res = Array(tp, length(pts))

    j = 0
    for (subpts, rng) in supported(b, pts)
        out = evaluate_raw(b, subpts, b.deriv, rng)
        for i in 1:length(subpts)
            res[j+=1] = (out[:,i], rng)
        end
    end

    res
end

function Base.call{B<:BasisFunction1D, T<:Real}(b::B, pts::Vector{T})
    res = zeros(Float64, length(pts))

    i = 1
    for (subpts, rng) in supported(b.basis, pts)
        i += length(subpts)
        if b.index ∉ rng continue end

        out = evaluate_raw(b.basis, subpts, b.deriv, rng)
        res[i-length(subpts):i-1] = out[findin(rng, b.index), :]
    end

    res
end

function Base.call{B<:Basis1D, S<:Real, T<:Real}(b::B, pt::S, coeffs::Vector{T})
    (vals, idxs) = b(pt)
    dot(vals, coeffs[idxs])
end

function Base.call{B<:Basis1D, S<:Real, T<:Real}(b::B, pt::S, coeffs::Matrix{T})
    (vals, idxs) = b(pt)
    vals' * coeffs[idxs,:]
end

Base.call{B<:Basis1D, S<:Real, T<:Real}(b::B, pts::Vector{S}, coeffs::Vector{T}) =
    Float64[dot(vals, coeffs[idxs]) for (vals, idxs) in b(pts)]

function Base.call{B<:Basis1D, S<:Real, T<:Real}(b::B, pts::Vector{S}, coeffs::Matrix{T})
    res = zeros(Float64, length(pts), size(coeffs, 2))

    for (i, (vals, idxs)) in enumerate(b(pts))
        res[i,:] = vals' * coeffs[idxs,:]
    end

    res
end

# Include specific bases
# ========================================================================

include("BSplineBasis.jl")
include("NURBSBasis.jl")

end  # module Bases


text:
# The next line is a deliberate error
sqrt(-3)


text:
"""
Encapsulates all data needed to transform a vector of values into the interval [0.0, 1.0]
by a linear and monotonous transformation.

See `intensities_encode`, `intensities_decode` for the usage.
"""
struct IntensityTransformation
   bottomval::Float64
   topval::Float64

   function IntensityTransformation(v::AbstractVector{Float64},
         q1::Float64 = 0.0, q2::Float64 = 1.0 - q1)

      bottomval = quantile(v, q1)
      topval = quantile(v, q2)
      new(bottomval, topval)
   end
end

"""
    intensities(x)
    intensities(x, q1)
    intensities(x, q1, q2)
Performs a linear and monotonous transformation on the data set `x`
to fit it the values into the interval [0.0, 1.0].
For more information see `intensities_encode`, `intensities_decode`.
"""
function intensities(x::Matrix{Float64}, q1::Float64 = 0.0, q2::Float64 = 1.0 - q1)
   intensities_encode(x, q1, q2)[1]
end

"""
    intensities_decode(x, its)
Backtransforms the intensity values in the data set `x`
(values in the interval [0.0, 1.0])` to the range of the original values
and returns the new data set or vector.
The `its` argument contains the information about the transformation, as it is returned by
`intensities_encode`.

Note that the range is truncated if the original transformation used other quantiles
than 0.0 or 1.0 (minimum and maximum).

# Example:
    x = randn(5, 4)
    xint, its = intensities_encode(x, 0.05)
    dbm = fitdbm(xint)
    xgen = samples(dbm, 5)
    intensities_decode(xgen, its)
"""
function intensities_decode(v::AbstractVector{Float64}, it::IntensityTransformation)
   v .*= it.topval - it.bottomval
   v .+= it.bottomval
end

function intensities_decode(x::AbstractMatrix{Float64},
      its::AbstractVector{IntensityTransformation})

   if length(its) != size(x, 2)
      error("Number of transformations is not equal to the number of columns " *
            "that are to be transformed.")
   end
   mapreduce(i -> intensities_decode(x[:, i], its[i]), hcat, 1:length(its))
end

"""
    intensities_encode(x)
    intensities_encode(x, q1)
    intensities_encode(x, q1, q2)
Performs a linear and monotonous transformation on the data set `x`
to fit it into the interval [0.0, 1.0].
It returns the transformed data set as a first result and
the information to reverse the tranformation as a second result.
If you are only interested in the transformed values, you can use the function `intensities`.

If `q1` is specified, all values below or equal to the quantile specified  by `q1`
are mapped to 0.0.
All values above or equal to the quantile specified by `q2` are mapped to 1.0.
`q2` defaults to `1 - q1`.

The quantiles are calculated per column/variable.

See also `intensities_decode` for the reverse transformation.
"""
function intensities_encode(x::AbstractMatrix{Float64},
      q1::Float64 = 0.0, q2::Float64 = 1.0 - q1)

   intensitytransformations = mapreduce(
         i -> IntensityTransformation(x[:, i], q1, q2),
         vcat, 1:size(x, 2))

   function transform(it::IntensityTransformation, v::AbstractVector{Float64})
      v[v .< it.bottomval] .= it.bottomval
      v[v .> it.topval] .= it.topval
      v .-= it.bottomval
      v ./= it.topval - it.bottomval
   end

   transformed = mapreduce(i -> transform(intensitytransformations[i], x[:, i]),
         hcat, 1:length(intensitytransformations))

   transformed, intensitytransformations
end

"""
    oneornone_decode(x, categories)
Returns a dataset such that
`x .== oneornone_decode(oneornone_encode(x, categories), categories)`.

For more, see `oneornone_encode`.
"""
function oneornone_decode(x::Matrix{Float64}, ncategories::Int)
   ncategoricalvariables, r = divrem(size(x, 2), ncategories - 1)
   if r != 0
      error("Wrong number of categories specified.")
   end
   oneornone_decode(x, fill(ncategories, ncategoricalvariables))
end

function oneornone_decode(x::Matrix{Float64}, nscategories::Vector{Int})
   varranges = ranges(nscategories .- 1)
   nsamples = size(x, 1)
   ncategoricalvariables = length(nscategories)
   xcategorical = zeros(nsamples, ncategoricalvariables)
   for k in 1:ncategoricalvariables
      for i in 1:nsamples
         category = findfirst(x[i, varranges[k]] .== 1.0)
         if category == nothing
            xcategorical[i, k] = 0.0
         else
            xcategorical[i, k] = category
         end
      end
   end
   xcategorical
end

"""
    oneornone_encode(x, categories)
Expects a data set `x` containing values 0.0, 1.0, 2.0 ... encoding the categories.
Returns a data set that encodes the variables/columns in `x` in multiple columns
with only values 0.0 and 1.0, similiar to the one-hot encoding with the deviation that
a zero is encoded as all-zeros.

The `categories` can be specified as
 * integer number if all variables have the same number of categories or as
 * integer vector, containing for each variable the number of categories encoded.

See also `oneornone_decode` for the reverse transformation.
"""
function oneornone_encode(x::M, nscategories::Vector{Int}
      ) where {N <: Number, M <:AbstractArray{N,2}}

   nsamples, norigvariables = size(x)
   n10variables = sum(nscategories .- 1)
   xout = zeros(nsamples, n10variables)
   varoffset = 0
   for k in 1:norigvariables
      for i in 1:nsamples
         if x[i, k] > 0.0
            xout[i, varoffset + Int(x[i,k])] = 1.0
         end
      end
      varoffset += nscategories[k] - 1
   end
   xout
end

function oneornone_encode(x, ncategories::Int
      ) where {N <: Number, M <:AbstractArray{N,2}}

   oneornone_encode(x, fill(ncategories, size(x, 2)))
end

"""
    splitdata(x, ratio)
Splits the data set `x` randomly in two data sets `x1` and `x2`, such that
the fraction of samples in `x2` is equal to (or as close as possible to) the
given `ratio`.

# Example:
    trainingdata, testdata = splitdata(data, 0.1) # Use 10 % as test data
"""
function splitdata(x::Matrix{Float64}, ratio::Float64)
   nsamples = size(x,1)
   testidxs = randperm(nsamples)[1:(round(Int, nsamples*ratio))]
   xtest = x[testidxs,:]
   xtraining = x[setdiff(1:nsamples, testidxs),:]
   xtraining, xtest
end


text:
@testset "multiple alleles" begin
global tree, df, d, net, currT

@testset "test: map alleles to species" begin
    tree = readTopology("(6,(5,(7,(3,4))));");
    PhyloNetworks.expandLeaves!(["7"],tree)
    @test writeTopology(tree) == "(6,(5,((7:0.0,7__2:0.0):1.0,(3,4))));"
    PhyloNetworks.mergeLeaves!(tree)
    @test writeTopology(tree) == "(6,(5,(7:1.0,(3,4))));"
    alleleDF=DataFrame(allele=["1","2"], species=["7","7"])
    # df = CSV.read(joinpath(@__DIR__, "..", "examples", "tableCFCI.csv"))
    # PhyloNetworks.mapAllelesCFtable!(df,alleleDF,[1,2,3,4],true,"CFmapped.csv")
    CSV.write("tmp.csv", alleleDF);
    df = (@test_logs (:warn, r"^not all alleles were mapped") mapAllelesCFtable("tmp.csv",
      joinpath(@__DIR__, "..", "examples", "tableCFCI.csv"),
      # joinpath(dirname(pathof(PhyloNetworks)), "..", "examples", "tableCFCI.csv"),
      filename="CFmapped.csv"))
    rm("CFmapped.csv")
    rm("tmp.csv")
    @test df[!,:t4] == ["4","7","3","7","3","3","7","3","3","3","7","3","3","3","3"]
end

#----------------------------------------------------------#
#   testing sorting of taxa and CFs                        #
#----------------------------------------------------------#
@testset "sorttaxa!" begin

letters = ["a","b","c","d"]; cfvalues = [0.6, 0.39, 0.01] # for ab_cd, ac_bd, ad_bc
d = DataFrame(t1=Array{String}(undef,24),t2=Array{String}(undef,24),t3=Array{String}(undef,24),t4=Array{String}(undef,24),
              CF12_34=Array{Float64}(undef,24), CF13_24=Array{Float64}(undef,24), CF14_23=Array{Float64}(undef,24));
irow=1        # d will contain 6!=24 rows: for all permutations on 4 letters
for i1 in 1:4
  ind234 = deleteat!(collect(1:4),i1)
  for i2 in ind234
    ind34 = deepcopy(ind234)
    deleteat!(ind34, findfirst(isequal(i2), ind34))
    for j in 1:2
      i3=ind34[j]; i4=ind34[3-j]
      d[irow,:t1]=letters[i1]; d[irow,:t2]=letters[i2]; d[irow,:t3]=letters[i3]; d[irow,:t4]=letters[i4]
      # CF12_34 corresponds to CFi1i2_i3i4
      if     (i1,i2)∈[(1,2),(2,1),(3,4),(4,3)] d[irow,:CF12_34] = cfvalues[1]
      elseif (i1,i2)∈[(1,3),(3,1),(2,4),(4,2)] d[irow,:CF12_34] = cfvalues[2]
      elseif (i1,i2)∈[(1,4),(4,1),(2,3),(3,2)] d[irow,:CF12_34] = cfvalues[3]
      end # next: set CF13_24
      if     (i1,i3)∈[(1,2),(2,1),(3,4),(4,3)] d[irow,:CF13_24] = cfvalues[1]
      elseif (i1,i3)∈[(1,3),(3,1),(2,4),(4,2)] d[irow,:CF13_24] = cfvalues[2]
      elseif (i1,i3)∈[(1,4),(4,1),(2,3),(3,2)] d[irow,:CF13_24] = cfvalues[3]
      end # nest: set CF14_23
      if     (i1,i4)∈[(1,2),(2,1),(3,4),(4,3)] d[irow,:CF14_23] = cfvalues[1]
      elseif (i1,i4)∈[(1,3),(3,1),(2,4),(4,2)] d[irow,:CF14_23] = cfvalues[2]
      elseif (i1,i4)∈[(1,4),(4,1),(2,3),(3,2)] d[irow,:CF14_23] = cfvalues[3]
      end
      irow += 1
    end
  end
end
# d
d2 = deepcopy(d);
sorttaxa!(d2);
d3 = DataFrame(t1=repeat([letters[1]],outer=[24]),t2=repeat([letters[2]],outer=[24]),
               t3=repeat([letters[3]],outer=[24]),t4=repeat([letters[4]],outer=[24]),
               CF12_34=repeat([cfvalues[1]],outer=[24]),CF13_24=repeat([cfvalues[2]],outer=[24]),CF14_23=repeat([cfvalues[3]],outer=[24]));
@test d2==d3

dat = readTableCF(d);
net = (@test_logs readTopologyLevel1("(a,((b)#H1,((#H1,c),d)));"));
# earlier warning: "net does not have identifiable branch lengths"
@test_logs topologyQPseudolik!(net, dat);
sorttaxa!(dat)

@test [q.obsCF for q in dat.quartet] == [[0.6,0.39,0.01] for i in 1:24]
@test [q.qnet.expCF for q in dat.quartet] == [[0.6915349833361827,0.12262648039048075,0.1858385362733365] for i in 1:24]
@test [q.taxon for q in dat.quartet] == [letters for i in 1:24]
@test [q.qnet.quartetTaxon for q in dat.quartet] == [letters for i in 1:24]

end # of testset: sorttaxa!

@testset "snaq on multiple alleles" begin
df=DataFrame(t1=["6","6","10","6","6","7","7","7","7","7","7"],
             t2=["7","7","7","10","7","7","7","7","7","7","7"],
             t3=["4","10","4","4","4","8","8","8","10","10","6"],
             t4=["8","8","8","8","10","10","4","6","4","6","4"],
             CF1234=[0.2729102510259939, 0.3967750546426937, 0.30161247267865315, 0.24693940689390592, 0.2729102510259939, 0.155181,  0.792153,  0.486702,  0.962734,  0.202531,  0.486886],
             CF1324=[0.45417949794801216, 0.30161247267865315, 0.30161247267865315, 0.5061211862121882, 0.45417949794801216, 0.673426 ,0.145408,  0.391103,  0.023078,  0.714826,  0.419015],
             CF1423=[0.2729102510259939, 0.30161247267865315, 0.3967750546426937, 0.24693940689390592, 0.2729102510259939, 0.171393,  0.062439,  0.122195,  0.014188,  0.082643,  0.094099])
d = readTableCF(df)
@test !isempty(d.repSpecies)
@test d.repSpecies == ["7"]

tree = "((6,4),(7,8),10);"
currT = readTopology(tree);

originalstdout = stdout
redirect_stdout(open("/dev/null", "w")) # not portable to Windows
estNet = snaq!(currT,d,hmax=1,seed=1010, runs=1, filename="", Nfail=10)
redirect_stdout(originalstdout)
@test 185.27 < estNet.loglik < 185.29 # or: wrong loglik
@test estNet.hybrid[1].k == 4 # or: wrong k
@test estNet.numTaxa == 5 # or: wrong number of taxa

redirect_stdout(open("/dev/null", "w")) # not portable to Windows
estNet = snaq!(currT,d,hmax=1,seed=8378, runs=1, filename="", Nfail=10,
               ftolAbs=1e-6,ftolRel=1e-5,xtolAbs=1e-4,xtolRel=1e-3)
redirect_stdout(originalstdout)
@test 174.58 < estNet.loglik < 174.59 # or: loglik wrong
@test estNet.hybrid[1].k == 5 # or: wrong k in hybrid
@test estNet.numTaxa == 5 # or: wrong # taxa

# net = snaq!(currT,d,hmax=1,seed=8378,filename="")
net = readTopology("(((4,#H7:::0.47411636966376686):0.6360197250223204,10):0.09464128563363322,(7:0.0,(6)#H7:::0.5258836303362331):0.36355727108454877,8);")
@test topologyQPseudolik!(net, d) ≈ 174.58674796123705
@test net.loglik ≈ 174.58674796123705
net = readTopology("(((4,#H1),10),(7,(6)#H1),8);")
net = topologyMaxQPseudolik!(net,d,  # loose tolerance for faster test
        ftolRel=1e-2,ftolAbs=1e-2,xtolAbs=1e-2,xtolRel=1e-2)
@test net.loglik > 174.5
end # test of snaq on multiple alleles

end # overall multiple allele sets of testests


text:
@testset "discriminator_loss" begin
    d_real = ones(Float32, 32, 32, 1, 1)
    d_g_fake_adv = zeros(Float32, 32, 32, 1, 1)
    @test SinGAN.discriminator_loss(d_real, d_g_fake_adv) ≈ 0f0
    @inferred Float32 SinGAN.discriminator_loss(d_real, d_g_fake_adv)
    @code_warntype SinGAN.discriminator_loss(d_real, d_g_fake_adv)

    d_real = zeros(Float32, 32, 32, 1, 1)
    d_g_fake_adv = ones(Float32, 32, 32, 1, 1)
    @test SinGAN.discriminator_loss(d_real, d_g_fake_adv) ≈ 2f0

    d_real = fill(0.5f0, 32, 32, 1, 1)
    d_g_fake_adv = fill(0.5f0, 32, 32, 1, 1)
    @test SinGAN.discriminator_loss(d_real, d_g_fake_adv) ≈ 0.5f0
end

@testset "generator_adv_loss" begin
    d_g_fake_adv = zeros(Float32, 32, 32, 1, 1)
    @test SinGAN.generator_adv_loss(d_g_fake_adv) ≈ 1f0
    @inferred Float32 SinGAN.generator_adv_loss(d_g_fake_adv)

    d_g_fake_adv = ones(Float32, 32, 32, 1, 1)
    @test SinGAN.generator_adv_loss(d_g_fake_adv) ≈ 0f0
end

@testset "update_discriminator!" begin
    opt = ADAM(0.0005f0)
    dscr = SinGAN.build_single_discriminator(3, 32)
    real_img = ones(Float32, 2, 2, 3, 1)
    g_fake_adv = randn(Float32, 2, 2, 3, 1)
    @info SinGAN.update_discriminator!(opt, dscr, real_img, g_fake_adv)
    # @code_warntype(SinGAN.update_discriminator!(opt, dscr, real_img, g_fake_adv))
end

# @testset "update_generator_adv!" begin
#     opt = ADAM(0.0005f0)
#     dscr = SinGAN.build_single_discriminator(3, 4)
#     gen = SinGAN.build_single_generator(3, 4, 1)
#     prev_adv = SinGAN.resize_and_padding(randn(Float32, 2, 2, 3, 1), (2, 2), (4, 4))
#     noise_adv = randn(Float32, 4, 4, 3, 1)
#     @info SinGAN.update_generator_adv!(opt, dscr, gen, prev_adv, noise_adv)
#     @code_warntype SinGAN.update_generator_adv!(opt, dscr, dscr, prev_adv, noise_adv)
# end

# @testset "update_generator_rec!" begin
#     opt = ADAM(0.0005f0)
#     gen = SinGAN.build_single_generator(3, 4, 1)
#     real_img = ones(Float32, 2, 2, 3, 1)
#     prev_rec = randn(Float32, 4, 4, 3, 1) # padded
#     noise_rec = zero(prev_rec)
#     alpha = 10f0
#     @info SinGAN.update_generator_rec!(opt, gen, real_img, prev_rec, noise_rec, alpha)
#     @code_warntype SinGAN.update_generator_rec!(opt, gen, real_img, prev_rec, noise_rec, alpha)
# end

max_stars_repo_name: matsueushi/SinGAN

max_issues_repo_name: matsueushi/SinGAN

max_forks_repo_name: matsueushi/SinGAN


text:
using MicrochipRTDM
using FTD2XX
#=
38400
57600
115200
128000
224000
256000
460800
921600
=#
const projectdirectory = "c:/Users/Chris/MPLABXProjects/lbcmcp/PWM_RTDM.x"
const mapdict = microchip_parsemap(projectdirectory,"production")
const uartconfig1 = UARTConfiguration(921600,8,1,"n",readtimeout = 1000,writetimeout = 1000)
const r = 1
const address = mapdict[:offsetP3]
const buffer = Array(UInt16,1)
const buffer2 = Array(UInt8,128)
const buffer3 = Array(UInt8,9)

function f()
  io = open(FT_DeviceIndex(0),uartconfig1)
  setlatencytimer(io,2); #Should speed things up a little
  rtdminterface = RTDMInterface(mapdict,io)
  for i::UInt16 in 1:1000
    isrtdmok(rtdminterface, retry = r)
    rtdm_read!(rtdminterface, buffer, address, retry = r)
    x = rtdm_read(rtdminterface, UInt16, address, retry = r)
    rtdm_read!(rtdminterface, buffer2, :SPI1CH1_BUFFER, retry = r)
    rtdm_write(rtdminterface, i, :SPI1CH1_BUFFER, retry = r)
    rtdm_write(rtdminterface, i, :PWM_PERIOD, retry = r)
    rtdm_write(rtdminterface, i, :PHASE_UI, retry = r)
    rtdm_write(rtdminterface, i, :DTRA, retry = r)
    rtdm_write(rtdminterface, i, :DTRB, retry = r)
  end
  close(rtdminterface)
  return nothing
end

f()

Profile.clear_malloc_data()
f()
quit()


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =================================================================================================
# Test writing a BMP image.
#
# A 32bit PNM is written as a PAM file so we can't compare those.
@testset "Sanity Check" begin
    local pix = pix_with()
    local formats = [
        pix_write_bmp(pix),
        pix_write_gif(pix),
        pix_write_jpeg(pix),
        pix_write_jp2k(pix),
        pix_write_pdf(pix),
        pix_write_png(pix),
        pix_write_pnm(pix),
        pix_write_ps(pix),
        pix_write_spix(pix)
    ]

    for i = 1:length(formats)
        for j = 1:length(formats)
            if i != j
                @test formats[i] != formats[j]
            end
        end
    end
end


text:
include("../src/RandomField.jl")
using .RandomField

# choose parameters:
## remember, we need that 4|n_i for i=1,...,d
N = [400,400]		# number pixels per dimension
L = [20.0,80.0]		# real space box is [-l,l]^d
f(V::Vector{T} where T<:Real) = gaus_rf(V)
					# spectral density

# derived parameters
d = length(N) 		# dimension
H = (pi/2) .* N ./L
var_y(V::Vector{T} where T<:Real) = 2*(2pi)^(d/2)*prod(L)/pi^d*f(V)

# allocate memory
V = Cloud(d,N)
X = Cloud(d,N)
Y = Cscalar(V)
Z = Cscalar(X)

# initialize grids
init!(X,L)
init!(V,H)

draw_realization!(Y,var_y)

# Y ↦ Z
ift!(Z,Y)

# print:
tovtk(Y, "../data/v.y.vtk")
tocsv(slice(Y,1), "../data/v.y1.csv")
tovtk(Z, "../data/x.z.vtk")
tocsv(slice(Z,1), "../data/x.z1.csv")


text:
"""
    VecColumnTable <: AbstractColumns

A `Tables.jl`-compatible column table that stores data as `Vector{AbstractVector}`
and column names as `Vector{Symbol}`.
Retrieving columns by column names is achieved with a `Dict{Symbol,Int}`
that maps names to indices.

This table type is designed for retrieving and iterating
dynamically generated columns for which
specialization on names and order of columns are not desired.
It is not intended to be directly constructed for interactive usage.
"""
struct VecColumnTable <: AbstractColumns
    columns::Vector{AbstractVector}
    names::Vector{Symbol}
    lookup::Dict{Symbol,Int}
    function VecColumnTable(columns::Vector{AbstractVector}, names::Vector{Symbol},
            lookup::Dict{Symbol,Int})
        # Assume Symbols in names and lookup match
        length(columns) == length(names) == length(lookup) ||
            throw(ArgumentError("arguments do not share the same length"))
        return new(columns, names, lookup)
    end
end

function VecColumnTable(columns::Vector{AbstractVector}, names::Vector{Symbol})
    lookup = Dict{Symbol,Int}(names.=>keys(names))
    return VecColumnTable(columns, names, lookup)
end

function VecColumnTable(data, subset=nothing)
    Tables.istable(data) || throw(ArgumentError("input data is not `Tables.jl`-compatible"))
    names = collect(Tables.columnnames(data))
    ncol = length(names)
    columns = Vector{AbstractVector}(undef, ncol)
    cols = Tables.columns(data)
    if subset === nothing
        @inbounds for i in keys(names)
            columns[i] = Tables.getcolumn(cols, i)
        end
    else
        @inbounds for i in keys(names)
            columns[i] = view(Tables.getcolumn(cols, i), subset)
        end
    end
    return VecColumnTable(columns, names)
end

function VecColumnTable(cols::VecColumnTable, subset=nothing)
    subset === nothing && return cols
    columns = similar(_columns(cols))
    @inbounds for i in keys(columns)
        columns[i] = view(cols[i], subset)
    end
    return VecColumnTable(columns, _names(cols), _lookup(cols))
end

_columns(cols::VecColumnTable) = getfield(cols, :columns)
_names(cols::VecColumnTable) = getfield(cols, :names)
_lookup(cols::VecColumnTable) = getfield(cols, :lookup)

ncol(cols::VecColumnTable) = length(_names(cols))
nrow(cols::VecColumnTable) = ncol(cols) > 0 ? length(_columns(cols)[1])::Int : 0

Base.size(cols::VecColumnTable) = (nrow(cols), ncol(cols))
function Base.size(cols::VecColumnTable, i::Integer)
    if i == 1
        nrow(cols)
    elseif i == 2
        ncol(cols)
    else
        throw(ArgumentError("VecColumnTable only have two dimensions"))
    end
end

Base.length(cols::VecColumnTable) = ncol(cols)
Base.isempty(cols::VecColumnTable) = size(cols, 1) == 0 || size(cols, 2) == 0

@inline function Base.getindex(cols::VecColumnTable, i::Int)
    @boundscheck if !checkindex(Bool, axes(_columns(cols), 1), i)
        throw(BoundsError(cols, i))
    end
    @inbounds return _columns(cols)[i]
end

@inline Base.getindex(cols::VecColumnTable, n::Symbol) = cols[_lookup(cols)[n]]
@inline Base.getindex(cols::VecColumnTable, ns::AbstractArray{Symbol}) = map(n->cols[n], ns)
@inline Base.getindex(cols::VecColumnTable, ::Colon) = _columns(cols)[:]
@inline function Base.getindex(cols::VecColumnTable, I)
    @boundscheck if !checkindex(Bool, axes(_columns(cols), 1), I)
        throw(BoundsError(cols, I))
    end
    @inbounds return _columns(cols)[I]
end

Base.values(cols::VecColumnTable) = _columns(cols)
Base.haskey(cols::VecColumnTable, key::Symbol) = haskey(_lookup(cols), key)
Base.haskey(cols::VecColumnTable, i::Int) = 0 < i <= length(_names(cols))

function Base.:(==)(x::VecColumnTable, y::VecColumnTable)
    size(x) == size(y) || return false
    _names(x) == _names(y) || return false
    eq = true
    for i in 1:size(x, 2)
        # missing could arise
        coleq = x[i] == y[i]
        isequal(coleq, false) && return false
        eq &= coleq
    end
    return eq
end

Base.view(cols::VecColumnTable, I) = VecColumnTable(cols, I)

Base.show(io::IO, cols::VecColumnTable) = print(io, nrow(cols), '×', ncol(cols), " VecColumnTable")

function Base.show(io::IO, ::MIME"text/plain", cols::VecColumnTable)
    show(io, cols)
    if ncol(cols) > 0 && nrow(cols) > 0
        println(io, ':')
        nr, nc = displaysize(io)
        pretty_table(IOContext(io, :limit=>true, :displaysize=>(nr-3, nc)),
            cols, vlines=1:1, hlines=1:1,
            show_row_number=true, show_omitted_cell_summary=false,
            newline_at_end=false, vcrop_mode=:middle)
    end
end

Base.summary(cols::VecColumnTable) =
    string(nrow(cols), '×', ncol(cols), ' ', nameof(typeof(cols)))
Base.summary(io::IO, cols::VecColumnTable) = print(io, summary(cols))

Tables.getcolumn(cols::VecColumnTable, i::Int) = cols[i]
Tables.getcolumn(cols::VecColumnTable, n::Symbol) = cols[n]
Tables.columnnames(cols::VecColumnTable) = _names(cols)

Tables.schema(cols::VecColumnTable) =
    Tables.Schema{(_names(cols)...,), Tuple{(eltype(col) for col in _columns(cols))...}}()
Tables.materializer(::VecColumnTable) = VecColumnTable

Tables.columnindex(cols::VecColumnTable, n::Symbol) = _lookup(cols)[n]
Tables.columntype(cols::VecColumnTable, n::Symbol) = eltype(cols[n])
Tables.rowcount(cols::VecColumnTable) = nrow(cols)

"""
    subcolumns(data, names, rows=Colon(); nomissing=true)

Construct a [`VecColumnTable`](@ref) from `data`
using columns specified with `names` over selected `rows`.

By default, columns are converted to drop support for missing values.
When possible, resulting columns share memory with original columns.
"""
function subcolumns(data, names, rows=Colon(); nomissing=true)
    Tables.istable(data) || throw(ArgumentError("input data is not `Tables.jl`-compatible"))
    names = names isa Vector{Symbol} ? names : Symbol[names...]
    ncol = length(names)
    columns = Vector{AbstractVector}(undef, ncol)
    lookup = Dict{Symbol,Int}()
    @inbounds for i in keys(names)
        col = view(Tables.getcolumn(data, names[i]), rows)
        nomissing && (col = disallowmissing(col))
        columns[i] = col
        lookup[names[i]] = i
    end
    return VecColumnTable(columns, names, lookup)
end

subcolumns(data, names::Symbol, rows=Colon(); nomissing=true) =
    subcolumns(data, [names], rows, nomissing=nomissing)

const VecColsRow = Tables.ColumnsRow{VecColumnTable}

ncol(r::VecColsRow) = length(r)

_rowhash(cols::Tuple{AbstractVector}, r::Int, h::UInt=zero(UInt))::UInt =
    hash(cols[1][r], h)

function _rowhash(cols::Tuple{Vararg{AbstractVector}}, r::Int, h::UInt=zero(UInt))::UInt
    h = hash(cols[1][r], h)
    _rowhash(Base.tail(cols), r, h)
end

# hash is implemented following DataFrames.DataFrameRow for getting unique row values
# Column names are not taken into account
Base.hash(r::VecColsRow, h::UInt=zero(UInt)) =
    _rowhash(ntuple(c->Tables.getcolumns(r)[c], ncol(r)), Tables.getrow(r), h)

# Column names are not taken into account
function Base.isequal(r1::VecColsRow, r2::VecColsRow)
    length(r1) == length(r2) || return false
    return all(((a, b),) -> isequal(a, b), zip(r1, r2))
end

# Column names are not taken into account
function Base.isless(r1::VecColsRow, r2::VecColsRow)
    length(r1) == length(r2) ||
        throw(ArgumentError("compared VecColsRow do not have the same length"))
    for (a, b) in zip(r1, r2)
        isequal(a, b) || return isless(a, b)
    end
    return false
end

Base.sortperm(cols::VecColumnTable; @nospecialize(kwargs...)) =
    sortperm(collect(Tables.rows(cols)); kwargs...)

# names and lookup are not copied
function Base.sort(cols::VecColumnTable; @nospecialize(kwargs...))
    p = sortperm(cols; kwargs...)
    columns = similar(_columns(cols))
    i = 0
    @inbounds for col in cols
        i += 1
        columns[i] = col[p]
    end
    return VecColumnTable(columns, _names(cols), _lookup(cols))
end

function Base.sort!(cols::VecColumnTable; @nospecialize(kwargs...))
    p = sortperm(cols; kwargs...)
    @inbounds for col in cols
        col .= col[p]
    end
end

"""
    apply(d, by::Pair)

Apply a function elementwise to the specified column(s)
in a `Tables.jl`-compatible table `d` and return the result.

Depending on the argument(s) accepted by a function `f`,
it is specified with argument `by` as either `column_index => f` or `column_indices => f`
where `column_index` is either a `Symbol` or `Int` for a column in `d`
and `column_indices` is an iterable collection of such indices for multiple columns.
`f` is applied elementwise to each specified column
to obtain an array of returned values.
"""
apply(d, by::Pair{Symbol,<:Function}) =
    map(by[2], Tables.getcolumn(d, by[1]))

apply(d, @nospecialize(by::Pair)) =
    map(by[2], (Tables.getcolumn(d, c) for c in by[1])...)

"""
    apply_and!(inds::BitVector, d, by::Pair)
    apply_and!(inds::BitVector, d, bys::Pair...)

Apply a function that returns `true` or `false` elementwise
to the specified column(s) in a `Tables.jl`-compatible table `d`
and then update the elements in `inds` through bitwise `and` with the returned array.
If an array instead of a function is provided,
elementwise equality (`==`) comparison is applied between the column and the array.
See also [`apply_and`](@ref).

The way a function is specified is the same as how it is done with [`apply`](@ref).
If multiple `Pair`s are provided,
`inds` are updated for each returned array through bitwise `and`.
"""
function apply_and!(inds::BitVector, d, by::Pair{Symbol,<:Function})
    aux = map(by[2], Tables.getcolumn(d, by[1]))
    inds .&= aux
    return inds
end

function apply_and!(inds::BitVector, d, @nospecialize(by::Pair))
    if by[2] isa Function
        aux = map(by[2], (Tables.getcolumn(d, c) for c in by[1])...)
    else
        aux = Tables.getcolumn(d, by[1]).==by[2]
    end
    inds .&= aux
    return inds
end

function apply_and!(inds::BitVector, d, @nospecialize(bys::Pair...))
    for by in bys
        apply_and!(inds, d, by)
    end
    return inds
end

"""
    apply_and(d, by::Pair)
    apply_and(d, bys::Pair...)

Apply a function that returns `true` or `false` elementwise
to the specified column(s) in a `Tables.jl`-compatible table `d` and return the result.
If an array instead of a function is provided,
elementwise equality (`==`) comparison is applied between the column and the array.
See also [`apply_and!`](@ref).

The way a function is specified is the same as how it is done with [`apply`](@ref).
If multiple `Pair`s are provided,
the returned array is obtained by combining
arrays returned by each function through bitwise `and`.
"""
apply_and(d, by::Pair{Symbol,<:Any}) = Tables.getcolumn(d, by[1]).==by[2]

function apply_and(d, by::Pair{Symbol,<:Function})
    src = Tables.getcolumn(d, by[1])
    out = similar(BitArray, size(src))
    map!(by[2], out, src)
    return out
end

function apply_and(d, @nospecialize(by::Pair))
    if by[2] isa Function
        src = ((Tables.getcolumn(d, c) for c in by[1])...,)
        out = similar(BitArray, size(src[1]))
        map!(by[2], out, src...)
        return out
    else
        return Tables.getcolumn(d, by[1]).==by[2]
    end
end

function apply_and(d, @nospecialize(bys::Pair...))
    inds = apply_and(d, bys[1])
    apply_and!(inds, d, bys[2:end]...)
    return inds
end

_parse_subset(cols::VecColumnTable, by::Pair) = (inds = apply_and(cols, by); return inds)

function _parse_subset(cols::VecColumnTable, inds)
    eltype(inds) <: Pair || return inds
    inds = apply_and(cols, inds...)
    return inds
end

_parse_subset(::VecColumnTable, ::Colon) = Colon()

"""
    TableIndexedMatrix{T,M,R,C} <: AbstractMatrix{T}

Matrix with row and column indices that can be selected
based on row values in a `Tables.jl`-compatible table respectively.
This is useful when how elements are stored into the matrix
are determined by the rows of the tables.

# Parameters
- `T`: element type of the matrix.
- `M`: type of the matrix.
- `R`: type of the table paired with the row indices.
- `C`: type of the table paired with the column indices.

# Fields
- `m::M`: the matrix that stores the elements.
- `r::R`: a table with the same number of rows with `m`.
- `c::C`: a table of which the number of rows is equal to the number of columns of `m`.
"""
struct TableIndexedMatrix{T,M,R,C} <: AbstractMatrix{T}
    m::M
    r::R
    c::C
    function TableIndexedMatrix(m::AbstractMatrix, r, c)
        Tables.istable(r) ||
            throw(ArgumentError("r is not `Tables.jl`-compatible"))
        Tables.istable(c) ||
            throw(ArgumentError("c is not `Tables.jl`-compatible"))
        size(m, 1) == Tables.rowcount(r) ||
            throw(DimensionMismatch("m and r do not have the same number of rows"))
        size(m, 2) == Tables.rowcount(c) || throw(DimensionMismatch(
            "the number of columns of m does not match the number of rows of c"))
        return new{eltype(m), typeof(m), typeof(r), typeof(c)}(m, r, c)
    end
end

Base.size(tm::TableIndexedMatrix) = size(tm.m)
Base.getindex(tm::TableIndexedMatrix, i) = tm.m[i]
Base.getindex(tm::TableIndexedMatrix, i, j) = tm.m[i,j]
Base.IndexStyle(::Type{<:TableIndexedMatrix{T,M}}) where {T,M} = IndexStyle(M)


text:
#
# Auxiliary Functions
#
# Not exactly related, but not exactly not

# make sure all dims are valid (in the Factor)
@inline function _check_dims_valid(dims::NodeNames, ϕ::Factor)
    isempty(dims) && return

    dim = first(dims)
    (dim in ϕ) || not_in_factor_error(dim)

    return _check_dims_valid(dims[2:end], ϕ)
end

# dims are unique
_ckeck_dims_unique(dims::NodeNames) = allunique(dims) || non_unique_dims_error()

"""
    duplicate(A, dims)

Repeates an array only through higer dimensions `dims`.

Custom version of repeate, but only outer repetition, and only duplicates
the array for the number of times specified in `dims` for dimensions greater
than `ndims(A)`. If `dims` is empty, returns a copy of `A`.

```jldoctest
julia> duplicate(collect(1:3), 2)
3×2 Array{Int64,2}:
 1  1
 2  2
 3  3

julia> duplicate([1 3; 2 4], 3)
2×2×3 Array{Int64,3}:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 1  3
 2  4

[:, :, 3] =
 1  3
 2  4
```
"""
function duplicate(A::Array, dims::Dims)
    size_in = size(A)

    length_in = length(A)
    size_out = (size_in..., dims...)::Dims

    B = similar(A, size_out)

    # zero in dims means no looping
    @simd for index in 1:prod(dims)
        unsafe_copy!(B, (index - 1) * length_in + 1, A, 1, length_in)
    end

    return B
end


text:
export EmbedThumbnail,
    EmbedVideo,
    EmbedImage,
    EmbedProvider,
    EmbedAuthor,
    EmbedFooter,
    EmbedField,
    Embed

"""
An [`Embed`](@ref)'s thumbnail image information.
More details [here](https://discordapp.com/developers/docs/resources/channel#embed-object-embed-thumbnail-structure).
"""
struct EmbedThumbnail
    url::Optional{String}
    proxy_url::Optional{String}
    height::Optional{Int}
    width::Optional{Int}
end
@boilerplate EmbedThumbnail :constructors :docs :lower :merge :mock

"""
An [`Embed`](@ref)'s video information.
More details [here](https://discordapp.com/developers/docs/resources/channel#embed-object-embed-video-structure).
"""
struct EmbedVideo
    url::Optional{String}
    height::Optional{Int}
    width::Optional{Int}
end
@boilerplate EmbedVideo :constructors :docs :lower :merge :mock

"""
An [`Embed`](@ref)'s image information.
More details [here](https://discordapp.com/developers/docs/resources/channel#embed-object-embed-image-structure).
"""
struct EmbedImage
    url::Optional{String}
    proxy_url::Optional{String}
    height::Optional{Int}
    width::Optional{Int}
end
@boilerplate EmbedImage :constructors :docs :lower :merge :mock

"""
An [`Embed`](@ref)'s provider information.
More details [here](https://discordapp.com/developers/docs/resources/channel#embed-object-embed-provider-structure).
"""
struct EmbedProvider
    name::Optional{String}
    url::OptionalNullable{String}  # Not supposed to be nullable.
end
@boilerplate EmbedProvider :constructors :docs :lower :merge :mock

"""
An [`Embed`](@ref)'s author information.
More details [here](https://discordapp.com/developers/docs/resources/channel#embed-object-embed-author-structure).
"""
struct EmbedAuthor
    name::Optional{String}
    url::Optional{String}
    icon_url::Optional{String}
    proxy_icon_url::Optional{String}
end
@boilerplate EmbedAuthor :constructors :docs :lower :merge :mock

"""
An [`Embed`](@ref)'s footer information.
More details [here](https://discordapp.com/developers/docs/resources/channel#embed-object-embed-footer-structure).
"""
struct EmbedFooter
    text::String
    icon_url::Optional{String}
    proxy_icon_url::Optional{String}
end
@boilerplate EmbedFooter :constructors :docs :lower :merge :mock

"""
An [`Embed`](@ref) field.
More details [here](https://discordapp.com/developers/docs/resources/channel#embed-object-embed-field-structure).
"""
struct EmbedField
    name::String
    value::String
    inline::Optional{Bool}
end
@boilerplate EmbedField :constructors :docs :lower :merge :mock

"""
A [`Message`](@ref) embed.
More details [here](https://discordapp.com/developers/docs/resources/channel#embed-object).
"""
struct Embed
    title::Optional{String}
    type::Optional{String}
    description::Optional{String}
    url::Optional{String}
    timestamp::Optional{DateTime}
    color::Optional{Int}
    footer::Optional{EmbedFooter}
    image::Optional{EmbedImage}
    thumbnail::Optional{EmbedThumbnail}
    video::Optional{EmbedVideo}
    provider::Optional{EmbedProvider}
    author::Optional{EmbedAuthor}
    fields::Optional{Vector{EmbedField}}
end
@boilerplate Embed :constructors :docs :lower :merge :mock