text:
"""
    struct Configuration{<:AbstractOrbital}

Represents a configuration -- a set of orbitals and their associated
occupation number.  Furthermore, each orbital can be in one of the
following _states_: `:open`, `:closed` or `:inactive`.

# Constructors

    Configuration(orbitals :: Vector{<:AbstractOrbital},
                  occupancy :: Vector{Int},
                  states :: Vector{Symbol}
                  [; sorted=false])

    Configuration(orbitals :: Vector{Tuple{<:AbstractOrbital, Int, Symbol}}
                  [; sorted=false])

In the first case, the parameters of each orbital have to be passed as
separate vectors, and the orbitals and occupancy have to be of the
same length. The `states` vector can be shorter and then the latter
orbitals that were not explicitly specified by `states` are assumed to
be `:open`.

The second constructor allows you to pass a vector of tuples instead,
where each tuple is a triplet `(orbital :: AbstractOrbital, occupancy
:: Int, state :: Symbol)` corresponding to each orbital.

In all cases, all the orbitals have to be distinct. The orbitals in
the configuration will be sorted (if `sorted`) according to the
ordering defined for the particular [`AbstractOrbital`](@ref).
"""
struct Configuration{O<:AbstractOrbital}
    orbitals::Vector{O}
    occupancy::Vector{Int}
    states::Vector{Symbol}
    sorted::Bool
    function Configuration(
        orbitals::Vector{O},
        occupancy::Vector{Int},
        states::Vector{Symbol}=[:open for o in orbitals];
        sorted=false) where {O<:AbstractOrbital}
        length(orbitals) == length(occupancy) ||
            throw(ArgumentError("Need to specify occupation numbers for all orbitals"))
        length(states) ≤ length(orbitals) ||
            throw(ArgumentError("Cannot provide more states than orbitals"))

        length(orbitals) == length(unique(orbitals)) ||
            throw(ArgumentError("Not all orbitals are unique"))

        if length(states) < length(orbitals)
            append!(states, repeat([:open], length(orbitals) - length(states)))
        end
        sd = setdiff(states, [:open, :closed, :inactive])
        isempty(sd) || throw(ArgumentError("Unknown orbital states $(sd)"))

        for i in eachindex(orbitals)
            occ = occupancy[i]
            occ < 1 && throw(ArgumentError("Invalid occupancy $(occ)"))
            orb = orbitals[i]
            degen_orb = degeneracy(orb)
            if occ > degen_orb
                if O <: Orbital || O <: SpinOrbital
                    throw(ArgumentError("Higher occupancy than possible for $(orb) with degeneracy $(degen_orb)"))
                elseif O <: RelativisticOrbital
                    orb_conj = flip_j(orb)
                    degen_orb_conj = degeneracy(orb_conj)
                    if occ == degen_orb + degen_orb_conj && orb_conj ∉ orbitals
                        occupancy[i] = degen_orb
                        push!(orbitals, orb_conj)
                        push!(occupancy, degen_orb_conj)
                        push!(states, states[i])
                    elseif orb_conj ∈ orbitals
                        throw(ArgumentError("Higher occupancy than possible for $(orb) with degeneracy $(degen_orb)"))
                    else
                        throw(ArgumentError("Can only specify higher occupancy for $(orb) if completely filling the $(orb),$(orb_conj) subshell (for total degeneracy of $(degen_orb+degen_orb_conj))"))
                    end
                end
            end
        end

        for i in eachindex(orbitals)
            states[i] == :closed && occupancy[i] < degeneracy(orbitals[i]) &&
                throw(ArgumentError("Can only close filled orbitals"))
        end

        p = sorted ? sortperm(orbitals) : eachindex(orbitals)

        new{O}(orbitals[p], occupancy[p], states[p], sorted)
    end
end

function Configuration(orbs::Vector{Tuple{O,Int,Symbol}}; sorted=false) where {O<:AbstractOrbital}
    orbitals = Vector{O}()
    occupancy = Vector{Int}()
    states = Vector{Symbol}()
    for (orb,occ,state) in orbs
        push!(orbitals, orb)
        push!(occupancy, occ)
        push!(states, state)
    end
    Configuration(orbitals, occupancy, states, sorted=sorted)
end

Configuration(orbital::O, occupancy::Int, state::Symbol=:open; sorted=false) where {O<:AbstractOrbital} =
    Configuration([orbital], [occupancy], [state], sorted=sorted)

Configuration{O}(;sorted=false) where {O<:AbstractOrbital} =
    Configuration(O[], Int[], Symbol[], sorted=sorted)

Base.copy(cfg::Configuration) =
    Configuration(copy(cfg.orbitals), copy(cfg.occupancy), copy(cfg.states),
                  sorted=cfg.sorted)

const RelativisticConfiguration = Configuration{<:RelativisticOrbital}

"""
    issorted(cfg::Configuration)

Tests if the orbitals of `cfg` is sorted.
"""
Base.issorted(cfg::Configuration) = cfg.sorted || issorted(cfg.orbitals)

"""
    sort(cfg::Configuration)

Returns a sorted copy of `cfg`.
"""
Base.sort(cfg::Configuration) =
    Configuration(cfg.orbitals, cfg.occupancy, cfg.states, sorted=true)

"""
    sorted(cfg::Configuration)

Returns `cfg` if it is already sorted or a sorted copy otherwise.
"""
sorted(cfg::Configuration) =
    issorted(cfg) ? cfg : sort(cfg)

"""
    issimilar(a::Configuration, b::Configuration)

Compares the electronic configurations `a` and `b`, only considering
the constituent orbitals and their occupancy, but disregarding their
ordering and states (`:open`, `:closed`, &c).

# Examples

```jldoctest
julia> a = c"1s 2s"
1s 2s

julia> b = c"2si 1s"
2sⁱ 1s

julia> issimilar(a, b)
true

julia> a==b
false
```

"""
function issimilar(a::Configuration{<:O}, b::Configuration{<:O}) where {O<:AbstractOrbital}
    ca = sorted(a)
    cb = sorted(b)
    ca.orbitals == cb.orbitals && ca.occupancy == cb.occupancy
end

"""
    ==(a::Configuration, b::Configuration)

Tests if configurations `a` and `b` are the same, considering orbital occupancy,
ordering, and states.

# Examples

```jldoctest
julia> c"1s 2s" == c"1s 2s"
true

julia> c"1s 2s" == c"1s 2si"
false

julia> c"1s 2s" == c"2s 1s"
false
```

"""
Base.:(==)(a::Configuration{<:O}, b::Configuration{<:O}) where {O<:AbstractOrbital} =
    a.orbitals == b.orbitals &&
    a.occupancy == b.occupancy &&
    a.states == b.states

Base.hash(c::Configuration, h::UInt) = hash(c.orbitals, hash(c.occupancy, hash(c.states, h)))

"""
    fill(c::Configuration)

Returns a corresponding configuration where the orbitals are completely filled (as
determined by [`degeneracy`](@ref)).

See also: [`fill!`](@ref)
"""
Base.fill(config::Configuration) = fill!(deepcopy(config))

"""
    fill!(c::Configuration)

Sets all the occupancies in configuration `c` to maximum, as determined by the
[`degeneracy`](@ref) function.

See also: [`fill`](@ref)
"""
function Base.fill!(config::Configuration)
    config.occupancy .= (degeneracy(o) for o in config.orbitals)
    return config
end

"""
    close(c::Configuration)

Return a corresponding configuration where where all the orbitals are marked `:closed`.

See also: [`close!`](@ref)
"""
Base.close(config::Configuration) = close!(deepcopy(config))

"""
    close!(c::Configuration)

Marks all the orbitals in configuration `c` as closed.

See also: [`close`](@ref)
"""
function close!(config::Configuration)
    if any(occ != degeneracy(o) for (o, occ, _) in config)
        throw(ArgumentError("Can't close $config due to unfilled orbitals."))
    end
    config.states .= :closed
    return config
end

function write_orbitals(io::IO, config::Configuration)
    ascii = get(io, :ascii, false)
    for (i,(orb,occ,state)) in enumerate(config)
        i > 1 && write(io, " ")
        show(io, orb)
        occ > 1 && write(io, ascii ? "$occ" : to_superscript(occ))
        state == :closed && write(io, ascii ? "c" : "ᶜ")
        state == :inactive && write(io, ascii ? "i" : "ⁱ")
    end
end

function Base.show(io::IO, config::Configuration{O}) where O
    ascii = get(io, :ascii, false)
    nc = length(config)
    if nc == 0
        write(io, ascii ? "empty" : "∅")
        return
    end
    noble_core_name = get_noble_core_name(config)
    core_config = core(config)
    ncc = length(core_config)
    if !isnothing(noble_core_name)
        write(io, "[$(noble_core_name)]")
        write(io, ascii ? "c" : "ᶜ")
        ngc = length(get_noble_gas(O, noble_core_name))
        core_config = core(config)
        if ncc > ngc
            write(io, ' ')
            write_orbitals(io, core_config[ngc+1:end])
        end
        nc > ncc && write(io, ' ')
    elseif ncc > 0
        write_orbitals(io, core_config)
    end
    write_orbitals(io, peel(config))
end

function Base.ascii(config::Configuration)
    io = IOBuffer()
    ctx = IOContext(io, :ascii=>true)
    show(ctx, config)
    String(take!(io))
end

"""
    get_noble_core_name(config::Configuration)

Returns the name of the noble gas with the most electrons whose configuration still forms
the first part of the closed part of `config`, or `nothing` if no such element is found.

```jldoctest
julia> AtomicLevels.get_noble_core_name(c"[He] 2s2")
"He"

julia> AtomicLevels.get_noble_core_name(c"1s2c 2s2c 2p6c 3s2c")
"Ne"

julia> AtomicLevels.get_noble_core_name(c"1s2") === nothing
true
```
"""
function get_noble_core_name(config::Configuration{O}) where O
    nc = length(config)
    nc == 0 && return nothing
    core_config = core(config)
    ncc = length(core_config)
    if length(core_config) > 0
        for gas in Iterators.reverse(noble_gases)
            gas_cfg = get_noble_gas(O, gas)
            ngc = length(gas_cfg)
            if ncc ≥ ngc && issimilar(core_config[1:length(gas_cfg)], gas_cfg)
                return gas
            end
        end
    end
    return nothing
end

function state_sym(state::AbstractString)
    if state == "c" || state == "ᶜ"
        :closed
    elseif state == "i" || state == "ⁱ"
        :inactive
    else
        :open
    end
end

function core_configuration(::Type{O}, element::AbstractString, state::AbstractString, sorted) where {O<:AbstractOrbital}
    element ∉ keys(noble_gases_configurations[O]) && throw(ArgumentError("Unknown noble gas $(element)"))
    state = state_sym(state == "" ? "c" : state) # By default, we'd like cores to be frozen
    core_config = noble_gases_configurations[O][element]
    Configuration(core_config.orbitals, core_config.occupancy,
                  [state for o in core_config.orbitals],
                  sorted=sorted)
end

function parse_orbital_occupation(::Type{O}, orb_str) where {O<:AbstractOrbital}
    m = match(r"^(([0-9]+|.)([a-z]|\[[0-9]+\])[-]{0,1})([0-9]*)([*ci]{0,1})$", orb_str)
    m2 = match(r"^(([0-9]+|.)([a-z]|\[[0-9]+\])[-]{0,1})([¹²³⁴⁵⁶⁷⁸⁹⁰]*)([ᶜⁱ]{0,1})$", orb_str)
    orb,occ,state = if !isnothing(m)
        m[1], m[4], m[5]
    elseif !isnothing(m2)
        m2[1], from_superscript(m2[4]), m2[5]
    else
        throw(ArgumentError("Unknown subshell specification $(orb_str)"))
    end
    parse(O, orb) , (occ == "") ? 1 : parse(Int, occ), state_sym(state)
end

function Base.parse(::Type{Configuration{O}}, conf_str; sorted=false) where {O<:AbstractOrbital}
    isempty(conf_str) && return Configuration{O}(sorted=sorted)
    orbs = split(conf_str, r"[\. ]")
    core_m = match(r"\[([a-zA-Z]+)\]([*ciᶜⁱ]{0,1})", first(orbs))
    if !isnothing(core_m)
        core_config = core_configuration(O, core_m[1], core_m[2], sorted)
        if length(orbs) > 1
            peel_config = Configuration(parse_orbital_occupation.(Ref(O), orbs[2:end]))
            Configuration(vcat(core_config.orbitals, peel_config.orbitals),
                          vcat(core_config.occupancy, peel_config.occupancy),
                          vcat(core_config.states, peel_config.states),
                          sorted=sorted)
        else
            core_config
        end
    else
        Configuration(parse_orbital_occupation.(Ref(O), orbs), sorted=sorted)
    end
end

function parse_conf_str(::Type{O}, conf_str, suffix) where O
    suffix ∈ ["", "s"] || throw(ArgumentError("Unknown configuration suffix $suffix"))
    parse(Configuration{O}, conf_str, sorted=suffix=="s")
end

"""
    @c_str -> Configuration{Orbital}

Construct a [`Configuration`](@ref), representing a non-relativistic
configuration, out of a string. With the added string macro suffix
`s`, the configuration is sorted.

# Examples

```jldoctest
julia> c"1s2 2s"
1s² 2s

julia> c"1s² 2s"
1s² 2s

julia> c"1s2.2s"
1s² 2s

julia> c"[Kr] 4d10 5s2 4f2"
[Kr]ᶜ 4d¹⁰ 5s² 4f²

julia> c"[Kr] 4d10 5s2 4f2"s
[Kr]ᶜ 4d¹⁰ 4f² 5s²
```
"""
macro c_str(conf_str, suffix="")
    parse_conf_str(Orbital, conf_str, suffix)
end

"""
    @rc_str -> Configuration{RelativisticOrbital}

Construct a [`Configuration`](@ref) representing a relativistic
configuration out of a string. With the added string macro suffix `s`,
the configuration is sorted.

# Examples

```jldoctest
julia> rc"[Ne] 3s 3p- 3p"
[Ne]ᶜ 3s 3p- 3p

julia> rc"[Ne] 3s 3p-2 3p4"
[Ne]ᶜ 3s 3p-² 3p⁴

julia> rc"[Ne] 3s 3p-² 3p⁴"
[Ne]ᶜ 3s 3p-² 3p⁴

julia> rc"2p- 1s"s
1s 2p-
```
"""
macro rc_str(conf_str, suffix="")
    parse_conf_str(RelativisticOrbital, conf_str, suffix)
end

"""
    @scs_str -> Vector{<:SpinConfiguration{<:Orbital}}

Generate all possible spin-configurations out of a string. With the
added string macro suffix `s`, the configuration is sorted.

# Examples

```jldoctest
julia> scs"1s2 2p2"
15-element Array{Configuration{SpinOrbital{Orbital{Int64},Tuple{Int64,HalfIntegers.Half{Int64}}}},1}:
 1s₀α 1s₀β 2p₋₁α 2p₋₁β
 1s₀α 1s₀β 2p₋₁α 2p₀α
 1s₀α 1s₀β 2p₋₁α 2p₀β
 1s₀α 1s₀β 2p₋₁α 2p₁α
 1s₀α 1s₀β 2p₋₁α 2p₁β
 1s₀α 1s₀β 2p₋₁β 2p₀α
 1s₀α 1s₀β 2p₋₁β 2p₀β
 1s₀α 1s₀β 2p₋₁β 2p₁α
 1s₀α 1s₀β 2p₋₁β 2p₁β
 1s₀α 1s₀β 2p₀α 2p₀β
 1s₀α 1s₀β 2p₀α 2p₁α
 1s₀α 1s₀β 2p₀α 2p₁β
 1s₀α 1s₀β 2p₀β 2p₁α
 1s₀α 1s₀β 2p₀β 2p₁β
 1s₀α 1s₀β 2p₁α 2p₁β
```
"""
macro scs_str(conf_str, suffix="")
    spin_configurations(parse_conf_str(Orbital, conf_str, suffix))
end

"""
    @rscs_str -> Vector{<:SpinConfiguration{<:RelativisticOrbital}}

Generate all possible relativistic spin-configurations out of a
string. With the added string macro suffix `s`, the configuration is
sorted.

# Examples

```jldoctest
julia> rscs"1s2 2p2"
6-element Array{Configuration{SpinOrbital{RelativisticOrbital{Int64},Tuple{HalfIntegers.Half{Int64}}}},1}:
 1s(-1/2) 1s(1/2) 2p(-3/2) 2p(-1/2)
 1s(-1/2) 1s(1/2) 2p(-3/2) 2p(1/2)
 1s(-1/2) 1s(1/2) 2p(-3/2) 2p(3/2)
 1s(-1/2) 1s(1/2) 2p(-1/2) 2p(1/2)
 1s(-1/2) 1s(1/2) 2p(-1/2) 2p(3/2)
 1s(-1/2) 1s(1/2) 2p(1/2) 2p(3/2)
```
"""
macro rscs_str(conf_str, suffix="")
    spin_configurations(parse_conf_str(RelativisticOrbital, conf_str, suffix))
end

Base.getindex(conf::Configuration{O}, i::Integer) where O =
    (conf.orbitals[i], conf.occupancy[i], conf.states[i])
Base.getindex(conf::Configuration{O}, i::Union{<:UnitRange{<:Integer},<:AbstractVector{<:Integer}}) where O =
    Configuration(Tuple{O,Int,Symbol}[conf[ii] for ii in i], sorted=conf.sorted)

Base.iterate(conf::Configuration{O}, (el, i)=(length(conf)>0 ? conf[1] : nothing,1)) where O =
    i > length(conf) ? nothing : (el, (conf[i==length(conf) ? i : i+1],i+1))

Base.length(conf::Configuration) = length(conf.orbitals)
Base.lastindex(conf::Configuration) = length(conf)
Base.eltype(conf::Configuration{O}) where O = (O,Int,Symbol)

function Base.write(io::IO, conf::Configuration{O}) where O
    write(io, 'c')
    if O <: Orbital
        write(io, 'n')
    elseif O <: RelativisticOrbital
        write(io, 'r')
    elseif O <: SpinOrbital
        write(io, 's')
    end
    write(io, length(conf))
    for (orb,occ,state) in conf
        write(io, orb, occ, first(string(state)))
    end
    write(io, conf.sorted)
end

function Base.read(io::IO, ::Type{Configuration})
    read(io, Char) == 'c' || throw(ArgumentError("Failed to read a Configuration from stream"))
    kind = read(io, Char)
    O = if kind == 'n'
        Orbital
    elseif kind == 'r'
        RelativisticOrbital
    elseif kind == 's'
        SpinOrbital
    else
        throw(ArgumentError("Unknown orbital type $(kind)"))
    end
    n = read(io, Int)
    occupations = Vector{Int}(undef, n)
    states = Vector{Symbol}(undef, n)
    orbitals = map(1:n) do i
        orb = read(io, O)
        occupations[i] = read(io, Int)
        s = read(io, Char)
        states[i] = if s == 'o'
            :open
        elseif s == 'c'
            :closed
        elseif s == 'i'
            :inactive
        else
            throw(ArgumentError("Unknown orbital state $(s)"))
        end
        orb
    end
    sorted = read(io, Bool)
    Configuration(orbitals, occupations, states, sorted=sorted)
end

function Base.isless(a::Configuration{<:O}, b::Configuration{<:O}) where {O<:AbstractOrbital}
    a = sorted(a)
    b = sorted(b)
    l = min(length(a),length(b))
    # If they are equal up to orbital l, designate the shorter config
    # as the smaller one.
    a[1:l] == b[1:l] && return length(a) == l
    norm_occ = (orb,w) -> 2w ≥ degeneracy(orb) ? degeneracy(orb) - w : w
    for ((orba,occa,statea),(orbb,occb,stateb)) in zip(a[1:l],b[1:l])
        if orba < orbb
            return true
        elseif orba == orbb
            # This is slightly arbitrary, but we designate the orbital
            # with occupation closest to a filled shell as the smaller
            # one.
            norm_occ(orba,occa) < norm_occ(orbb,occb) && return true
        else
            return false
        end
    end
    false
end

"""
    num_electrons(c::Configuration) -> Int

Return the number of electrons in the configuration.

```jldoctest
julia> num_electrons(c"1s2")
2

julia> num_electrons(rc"[Kr] 5s2 5p-2 5p2")
42
```
"""
num_electrons(c::Configuration) = sum(c.occupancy)

"""
    num_electrons(c::Configuration, o::AbstractOrbital) -> Int

Returns the number of electrons on orbital `o` in configuration `c`. If `o` is not part of
the configuration, returns `0`.

```jldoctest
julia> num_electrons(c"1s 2s2", o"2s")
2

julia> num_electrons(rc"[Rn] Qf-5 Pf3", ro"Qf-")
5

julia> num_electrons(c"[Ne]", o"3s")
0
```
"""
function num_electrons(c::Configuration, o::AbstractOrbital)
    idx = findfirst(isequal(o), c.orbitals)
    isnothing(idx) && return 0
    c.occupancy[idx]
end

"""
    in(o::AbstractOrbital, c::Configuration) -> Bool

Checks if orbital `o` is part of configuration `c`.

```jldoctest
julia> in(o"2s", c"1s2 2s2")
true

julia> o"2p" ∈ c"1s2 2s2"
false
```
"""
Base.in(orb::O, conf::Configuration{O}) where {O<:AbstractOrbital} =
    orb ∈ conf.orbitals

"""
    filter(f, c::Configuration) -> Configuration

Filter out the orbitals from configuration `c` for which the predicate `f` returns `false`.
The predicate `f` needs to take three arguments: `orbital`, `occupancy` and `state`.

```julia
julia> filter((o,occ,s) -> o.ℓ == 1, c"[Kr]")
2p⁶ᶜ 3p⁶ᶜ 4p⁶ᶜ
```
"""
Base.filter(f::Function, conf::Configuration) =
    conf[filter(j -> f(conf[j]...), eachindex(conf.orbitals))]

"""
    core(::Configuration) -> Configuration

Return the core configuration (i.e. the sub-configuration of all the orbitals that are
marked `:closed`).

```jldoctest
julia> core(c"1s2c 2s2c 2p6c 3s2")
[Ne]ᶜ

julia> core(c"1s2 2s2")
∅

julia> core(c"1s2 2s2c 2p6c")
2s²ᶜ 2p⁶ᶜ
```
"""
core(conf::Configuration) = filter((orb,occ,state) -> state == :closed, conf)

"""
    peel(::Configuration) -> Configuration

Return the non-core part of the configuration (i.e. orbitals not marked `:closed`).

```jldoctest
julia> peel(c"1s2c 2s2c 2p3")
2p³

julia> peel(c"[Ne] 3s 3p3")
3s 3p³
```
"""
peel(conf::Configuration) = filter((orb,occ,state) -> state != :closed, conf)

"""
    inactive(::Configuration) -> Configuration

Return the part of the configuration marked `:inactive`.

```jldoctest
julia> inactive(c"1s2c 2s2i 2p3i 3s2")
2s²ⁱ 2p³ⁱ
```
"""
inactive(conf::Configuration) = filter((orb,occ,state) -> state == :inactive, conf)

"""
    active(::Configuration) -> Configuration

Return the part of the configuration marked `:open`.

```jldoctest
julia> active(c"1s2c 2s2i 2p3i 3s2")
3s²
```
"""
active(conf::Configuration) = filter((orb,occ,state) -> state != :inactive, peel(conf))

"""
    bound(::Configuration) -> Configuration

Return the bound part of the configuration (see also [`isbound`](@ref)).

```jldoctest
julia> bound(c"1s2 2s2 2p4 Ks2 Kp1")
1s² 2s² 2p⁴
```
"""
bound(conf::Configuration) = filter((orb,occ,state) -> isbound(orb), conf)

"""
    continuum(::Configuration) -> Configuration

Return the non-bound (continuum) part of the configuration (see also [`isbound`](@ref)).

```jldoctest
julia> continuum(c"1s2 2s2 2p4 Ks2 Kp1")
Ks² Kp
```
"""
continuum(conf::Configuration) = filter((orb,occ,state) -> !isbound(orb), peel(conf))

Base.empty(conf::Configuration{O}) where {O<:AbstractOrbital} =
    Configuration(empty(conf.orbitals), empty(conf.occupancy), empty(conf.states), sorted=conf.sorted)

"""
    parity(::Configuration) -> Parity

Return the parity of the configuration.

```jldoctest
julia> parity(c"1s 2p")
odd

julia> parity(c"1s 2p2")
even
```

See also: [`Parity`](@ref)
"""
parity(conf::Configuration) = mapreduce(o -> parity(o[1])^o[2], *, conf)

"""
    replace(conf, a => b[; append=false])

Substitute one electron in orbital `a` of `conf` by one electron in
orbital `b`. If `conf` is unsorted the substitution is performed
in-place, unless `append`, in which case the new orbital is appended
instead.

# Examples

```jldoctest
julia> replace(c"1s2 2s", o"1s" => o"2p")
1s 2p 2s

julia> replace(c"1s2 2s", o"1s" => o"2p", append=true)
1s 2s 2p

julia> replace(c"1s2 2s"s, o"1s" => o"2p")
1s 2s 2p
```
"""
function Base.replace(conf::Configuration{O₁}, orbs::Pair{O₂,O₃};
                      append=false) where {O<:AbstractOrbital,O₁<:O,O₂<:O,O₃<:O}
    src,dest = orbs
    orbitals = promote_type(O₁,O₂,O₃)[]
    append!(orbitals, conf.orbitals)
    occupancy = copy(conf.occupancy)
    states = copy(conf.states)

    i = findfirst(isequal(src), orbitals)
    isnothing(i) && throw(ArgumentError("$(src) not present in $(conf)"))

    j = findfirst(isequal(dest), orbitals)
    if isnothing(j)
        j = (append ? length(orbitals) : i) + 1
        insert!(orbitals, j, dest)
        insert!(occupancy, j, 1)
        insert!(states, j, :open)
    else
        occupancy[j] == degeneracy(dest) &&
            throw(ArgumentError("$(dest) already maximally occupied in $(conf)"))
        occupancy[j] += 1
    end

    occupancy[i] -= 1
    if occupancy[i] == 0
        deleteat!(orbitals, i)
        deleteat!(occupancy, i)
        deleteat!(states, i)
    end

    Configuration(orbitals, occupancy, states, sorted=conf.sorted)
end

"""
    -(configuration::Configuration, orbital::AbstractOrbital[, n=1])

Remove `n` electrons in the orbital `orbital` from the configuration
`configuration`. If the orbital had previously been `:closed` or
`:inactive`, it will now be `:open`.
"""
function Base.:(-)(configuration::Configuration{O₁}, orbital::O₂, n::Int=1) where {O<:AbstractOrbital,O₁<:O,O₂<:O}
    orbitals = promote_type(O₁,O₂)[]
    append!(orbitals, configuration.orbitals)
    occupancy = copy(configuration.occupancy)
    states = copy(configuration.states)

    i = findfirst(isequal(orbital), orbitals)
    isnothing(i) && throw(ArgumentError("$(orbital) not present in $(configuration)"))

    occupancy[i] ≥ n ||
        throw(ArgumentError("Trying to remove $(n) electrons from orbital $(orbital) with occupancy $(occupancy[i])"))

    occupancy[i] -= n
    states[i] = :open
    if occupancy[i] == 0
        deleteat!(orbitals, i)
        deleteat!(occupancy, i)
        deleteat!(states, i)
    end

    Configuration(orbitals, occupancy, states, sorted=configuration.sorted)
end

"""
    +(::Configuration, ::Configuration)

Add two configurations together. If both configuration have an orbital, the number of
electrons gets added together, but in this case the status of the orbitals must match.

```jldoctest
julia> c"1s" + c"2s"
1s 2s

julia> c"1s" + c"1s"
1s²
```
"""
function Base.:(+)(a::Configuration{O₁}, b::Configuration{O₂}) where {O<:AbstractOrbital,O₁<:O,O₂<:O}
    orbitals = promote_type(O₁,O₂)[]
    append!(orbitals, a.orbitals)
    occupancy = copy(a.occupancy)
    states = copy(a.states)
    for (orb,occ,state) in b
        i = findfirst(isequal(orb), orbitals)
        if isnothing(i)
            push!(orbitals, orb)
            push!(occupancy, occ)
            push!(states, state)
        else
            occupancy[i] += occ
            states[i] == state || throw(ArgumentError("Incompatible states for $(orb): $(states[i]) and $state"))
        end
    end
    Configuration(orbitals, occupancy, states, sorted=a.sorted && b.sorted)
end

"""
    delete!(c::Configuration, o::AbstractOrbital)

Remove the entire subshell corresponding to orbital `o` from configuration `c`.

```jldoctest
julia> delete!(c"[Ar] 4s2 3d10 4p2", o"4s")
[Ar]ᶜ 3d¹⁰ 4p²
```
"""
function Base.delete!(c::Configuration{O}, o::O) where O <: AbstractOrbital
    idx = findfirst(isequal(o), [co[1] for co in c])
    idx === nothing && return c
    deleteat!(c.orbitals, idx)
    deleteat!(c.occupancy, idx)
    deleteat!(c.states, idx)
    return c
end

"""
    ⊗(::Union{Configuration, Vector{Configuration}}, ::Union{Configuration, Vector{Configuration}})

Given two collections of `Configuration`s, it creates an array of `Configuration`s with all
possible juxtapositions of configurations from each collection.

# Examples

```jldoctest
julia> c"1s" ⊗ [c"2s2", c"2s 2p"]
2-element Array{Configuration{Orbital{Int64}},1}:
 1s 2s²
 1s 2s 2p

julia> [rc"1s", rc"2s"] ⊗ [rc"2p-", rc"2p"]
4-element Array{Configuration{RelativisticOrbital{Int64}},1}:
 1s 2p-
 1s 2p
 2s 2p-
 2s 2p
```
"""
⊗(a::Vector{<:Configuration}, b::Vector{<:Configuration}) =
    [x+y for x in a for y in b]
⊗(a::Union{<:Configuration,Vector{<:Configuration}}, b::Configuration) =
    a ⊗ [b]
⊗(a::Configuration, b::Vector{<:Configuration}) =
    [a] ⊗ b

"""
    rconfigurations_from_orbital(n, ℓ, occupancy)

Generate all `Configuration`s with relativistic orbitals corresponding to the
non-relativistic orbital with `n` and `ℓ` quantum numbers, with given occupancy.

# Examples

```jldoctest
julia> AtomicLevels.rconfigurations_from_orbital(3, 1, 2)
3-element Array{Configuration{#s21} where #s21<:RelativisticOrbital,1}:
 3p-²
 3p- 3p
 3p²
```
"""
function rconfigurations_from_orbital(n::N, ℓ::Int, occupancy::Int) where {N<:MQ}
    n isa Integer && ℓ + 1 > n && throw(ArgumentError("ℓ=$ℓ too high for given n=$n"))
    occupancy > 2*(2ℓ + 1) && throw(ArgumentError("occupancy=$occupancy too high for given ℓ=$ℓ"))

    degeneracy_ℓm, degeneracy_ℓp = 2ℓ, 2ℓ + 2 # degeneracies of nℓ- and nℓ orbitals
    nlow_min = max(occupancy - degeneracy_ℓp, 0)
    nlow_max = min(degeneracy_ℓm, occupancy)
    confs = RelativisticConfiguration[]
    for nlow = nlow_max:-1:nlow_min
        nhigh = occupancy - nlow
        conf = if nlow == 0
            Configuration([RelativisticOrbital(n, ℓ, ℓ + 1//2)], [nhigh])
        elseif nhigh == 0
            Configuration([RelativisticOrbital(n, ℓ, ℓ - 1//2)], [nlow])
        else
            Configuration(
                [RelativisticOrbital(n, ℓ, ℓ - 1//2),
                 RelativisticOrbital(n, ℓ, ℓ + 1//2)],
                [nlow, nhigh]
            )
        end
        push!(confs, conf)
    end
    return confs
end

"""
    rconfigurations_from_orbital(orbital::Orbital, occupancy)

Generate all `Configuration`s with relativistic orbitals corresponding to the
non-relativistic version of the `orbital` with a given occupancy.

# Examples

```jldoctest
julia> AtomicLevels.rconfigurations_from_orbital(o"3p", 2)
3-element Array{Configuration{#s21} where #s21<:RelativisticOrbital,1}:
 3p-²
 3p- 3p
 3p²
```
"""
function rconfigurations_from_orbital(orbital::Orbital, occupation::Integer)
    rconfigurations_from_orbital(orbital.n, orbital.ℓ, occupation)
end

function rconfigurations_from_nrstring(orb_str::AbstractString)
    m = match(r"^([0-9]+|.)([a-z]+)([0-9]+)?$", orb_str)
    isnothing(m) && throw(ArgumentError("Invalid orbital string: $(orb_str)"))
    n = parse_orbital_n(m)
    ℓ = parse_orbital_ℓ(m)
    occupancy = isnothing(m[3]) ? 1 : parse(Int, m[3])
    return rconfigurations_from_orbital(n, ℓ, occupancy)
end

"""
    relconfigurations(c::Configuration{<:Orbital}) -> Vector{<:Configuration{<:RelativisticOrbital}}

Generate all relativistic configurations from the non-relativistic
configuration `c`, by applying [`rconfigurations_from_orbital`](@ref)
to each subshell and combining the results.
"""
function relconfigurations(c::Configuration{<:Orbital})
    rcs = map(c) do (orb,occ,state)
        orcs = rconfigurations_from_orbital(orb, occ)
        if state == :closed
            close!.(orcs)
        end
        orcs
    end
    if length(rcs) > 1
        rcs = map(Iterators.product(rcs...)) do subshells
            reduce(⊗, subshells)
        end
    end
    reduce(vcat, rcs)
end

relconfigurations(cs::Vector{<:Configuration{<:Orbital}}) =
    reduce(vcat, map(relconfigurations, cs))

"""
    @rcs_str -> Vector{Configuration{RelativisticOrbital}}

Construct a `Vector` of all `Configuration`s corresponding to the non-relativistic `nℓ`
orbital with the given occupancy from the input string.

The string is assumed to have the following syntax: `\$(n)\$(ℓ)\$(occupancy)`, where `n`
and `occupancy` are integers, and `ℓ` is in spectroscopic notation.

# Examples

```jldoctest
julia> rcs"3p2"
3-element Array{Configuration{#s21} where #s21<:RelativisticOrbital,1}:
 3p-²
 3p- 3p
 3p²
```
"""
macro rcs_str(s)
    rconfigurations_from_nrstring(s)
end

"""
    spin_configurations(configuration)

Generate all possible configurations of spin-orbitals from `configuration`, i.e. all
permissible values for the quantum numbers `n`, `ℓ`, `mℓ`, `ms` for each electron. Example:

```jldoctest
julia> spin_configurations(c"1s2")
1-element Array{Configuration{SpinOrbital{Orbital{Int64},Tuple{Int64,HalfIntegers.Half{Int64}}}},1}:
 1s₀α 1s₀β

julia> spin_configurations(c"1s2"s)
1-element Array{Configuration{SpinOrbital{Orbital{Int64},Tuple{Int64,HalfIntegers.Half{Int64}}}},1}:
 1s₀α 1s₀β

julia> spin_configurations(c"1s ks")
4-element Array{Configuration{SpinOrbital{#s16,Tuple{Int64,HalfIntegers.Half{Int64}}} where #s16<:Orbital},1}:
 1s₀α ks₀α
 1s₀β ks₀α
 1s₀α ks₀β
 1s₀β ks₀β
```
"""
function spin_configurations(cfg::Configuration{O}) where O
    states = Dict{O,Symbol}()
    orbitals = map(cfg) do (orb,occ,state)
        states[orb] = state
        sorbs = spin_orbitals(orb)
        collect(combinations(sorbs, occ))
    end
    map(Iterators.product(orbitals...)) do choice
        c = vcat(choice...)
        s = [states[orb.orb] for orb in c]
        Configuration(c, ones(Int,length(c)), s, sorted=cfg.sorted)
    end |> vec
end

Base.convert(::Type{Configuration{O}}, c::Configuration) where {O<:AbstractOrbital} =
    Configuration(Vector{O}(c.orbitals), c.occupancy, c.states, sorted=c.sorted)

Base.convert(::Type{Configuration{SO}}, c::Configuration{<:SpinOrbital}) where {SO<:SpinOrbital} =
    Configuration(Vector{SO}(c.orbitals), c.occupancy, c.states, sorted=c.sorted)

Base.promote_type(::Type{Configuration{SO}}, ::Type{Configuration{SO}}) where {SO<:SpinOrbital} =
    Configuration{SO}

Base.promote_type(::Type{CA}, ::Type{CB}) where {
    A<:SpinOrbital,CA<:Configuration{A},
    B<:SpinOrbital,CB<:Configuration{B}
} =
    Configuration{promote_type(A,B)}

Base.promote_type(::Type{Cfg}, ::Type{Union{}}) where {SO<:SpinOrbital,Cfg<:Configuration{SO}} = Cfg

"""
    spin_configurations(configurations)

For each configuration in `configurations`, generate all possible configurations of
spin-orbitals.
"""
function spin_configurations(cs::Vector{<:Configuration})
    scs = map(spin_configurations, cs)
    sort(reduce(vcat, scs))
end

"""
    SpinConfiguration

Specialization of [`Configuration`](@ref) for configurations
consisting of [`SpinOrbital`](@ref)s.
"""
const SpinConfiguration{O<:SpinOrbital} = Configuration{O}

function Base.parse(::Type{SpinConfiguration{SO}}, conf_str; sorted=false) where {O<:AbstractOrbital,SO<:SpinOrbital{O}}
    isempty(conf_str) && return SpinConfiguration{SO}(sorted=sorted)
    orbs = split(conf_str, r"[ ]")
    core_m = match(r"\[([a-zA-Z]+)\]([*ci]{0,1})", first(orbs))
    if !isnothing(core_m)
        core_config = first(spin_configurations(core_configuration(O, core_m[1], core_m[2], sorted)))
        if length(orbs) > 1
            peel_orbitals = parse.(Ref(SO), orbs[2:end])
            orbitals = vcat(core_config.orbitals, peel_orbitals)
            Configuration(orbitals,
                          ones(Int, length(orbitals)),
                          vcat(core_config.states, fill(:open, length(peel_orbitals))),
                          sorted=sorted)
        else
            core_config
        end
    else
        Configuration(parse.(Ref(SO), orbs), ones(Int, length(orbs)), sorted=sorted)
    end
end

"""
    @sc_str -> SpinConfiguration{<:SpinOrbital{<:Orbital}}

A string macro to construct a non-relativistic [`SpinConfiguration`](@ref).

# Examples

```jldoctest
julia> sc"1s₀α 2p₋₁β"
1s₀α 2p₋₁β

julia> sc"ks(0,-1/2) l[4](-3,1/2)"
ks₀β lg₋₃α
```
"""
macro sc_str(conf_str, suffix="")
    parse(SpinConfiguration{SpinOrbital{Orbital}}, conf_str, sorted=suffix=="s")
end

"""
    @rsc_str -> SpinConfiguration{<:SpinOrbital{<:RelativisticOrbital}}

A string macro to construct a relativistic [`SpinConfiguration`](@ref).

# Examples

```jldoctest
julia> rsc"1s(1/2) 2p(-1/2)"
1s(1/2) 2p(-1/2)

julia> rsc"ks(-1/2) l[4]-(-5/2)"
ks(-1/2) lg-(-5/2)
```
"""
macro rsc_str(conf_str, suffix="")
    parse(SpinConfiguration{SpinOrbital{RelativisticOrbital}}, conf_str, sorted=suffix=="s")
end

function Base.show(io::IO, c::SpinConfiguration)
    ascii = get(io, :ascii, false)
    nc = length(c)
    if nc == 0
        write(io, ascii ? "empty" : "∅")
        return
    end

    orbitals = Dict{Orbital, Vector{<:SpinOrbital}}()
    core_orbitals = sort(unique(map(o -> o.orb, core(c).orbitals)))
    core_cfg = Configuration(core_orbitals, ones(Int, length(core_orbitals))) |> fill |> close

    if !isempty(core_cfg)
        show(io, core_cfg)
        write(io, " ")
    end
    # if !c.sorted
    # In the unsorted case, we do not yet contract subshells for
    # printing; to be implemented.
    so = string.(peel(c).orbitals)
    write(io, join(so, " "))
    return
    # end
    # for orb in peel(c).orbitals
    #     orbitals[orb.orb] = push!(get(orbitals, orb.orb, SpinOrbital[]), orb)
    # end
    # map(sort(collect(keys(orbitals)))) do orb
    #     ℓ = orb.ℓ
    #     g = degeneracy(orb)
    #     sub_shell = orbitals[orb]
    #     if length(sub_shell) == g
    #         format("{1:s}{2:s}", orb, to_superscript(g))
    #     else
    #         map(Iterators.product()mℓrange(orb)) do mℓ
    #             mℓshell = findall(o -> o.mℓ == mℓ, sub_shell)
    #             if length(mℓshell) == 2
    #                 format("{1:s}{2:s}{3:s}", orb, to_subscript(mℓ), to_superscript(2))
    #             elseif length(mℓshell) == 1
    #                 string(sub_shell[mℓshell[1]])
    #             else
    #                 ""
    #             end
    #         end |> so -> join(filter(s -> !isempty(s), so), " ")
    #     end
    # end |> so -> write(io, join(so, " "))
end

"""
    substitutions(src::SpinConfiguration, dst::SpinConfiguration)

Find all orbital substitutions going from spin-configuration `src` to
configuration `dst`.
"""
function substitutions(src::Configuration{A}, dst::Configuration{B}) where {A<:SpinOrbital,B<:SpinOrbital}
    src == dst && return []
    # This is only valid for spin-configurations, since occupation
    # numbers are not dealt with.
    num_electrons(src) == num_electrons(dst) ||
        throw(ArgumentError("Configurations not of same amount of electrons"))
    r = Vector{Pair{A,B}}()
    missing = A[]
    same = Int[]
    for (i,(orb,occ,state)) in enumerate(src)
        j = findfirst(isequal(orb), dst.orbitals)
        if j === nothing
            push!(missing, orb)
        else
            push!(same, j)
        end
    end
    new = [j for j ∈ 1:num_electrons(dst)
           if j ∉ same]
    [mo => dst.orbitals[j] for (mo,j) in zip(missing, new)]
end

# We need to declare noble_gases first, with empty entries for Orbital and RelativisticOrbital
# since parse(Configuration, ...) uses it.
const noble_gases_configurations = Dict(
    O => Dict{String,Configuration{<:O}}()
    for O in [Orbital, RelativisticOrbital]
)
const noble_gases_configuration_strings = [
    "He" => "1s2",
    "Ne" => "[He] 2s2 2p6",
    "Ar" => "[Ne] 3s2 3p6",
    "Kr" => "[Ar] 3d10 4s2 4p6",
    "Xe" => "[Kr] 4d10 5s2 5p6",
    "Rn" => "[Xe] 4f14 5d10 6s2 6p6",
]
const noble_gases = [gas for (gas, _) in noble_gases_configuration_strings]
for (gas, configuration) in noble_gases_configuration_strings, O in [Orbital, RelativisticOrbital]
    noble_gases_configurations[O][gas] = parse(Configuration{O}, configuration)
end

# This construct is needed since when showing configurations, they
# will be specialized on the Orbital parameterization, which we cannot
# index noble_gases with.
for O in [Orbital, RelativisticOrbital]
    @eval get_noble_gas(::Type{<:$O}, k) = noble_gases_configurations[$O][k]
end

"""
    nonrelconfiguration(c::Configuration{<:RelativisticOrbital}) -> Configuration{<:Orbital}

Reduces a relativistic configuration down to the corresponding non-relativistic configuration.

```jldoctest
julia> c = rc"1s2 2p-2 2s 2p2 3s2 3p-"s
1s² 2s 2p-² 2p² 3s² 3p-

julia> nonrelconfiguration(c)
1s² 2s 2p⁴ 3s² 3p
```
"""
function nonrelconfiguration(c::Configuration{<:RelativisticOrbital})
    mq = Union{mqtype.(c.orbitals)...}
    nrorbitals, nroccupancies, nrstates = Orbital{<:mq}[], Int[], Symbol[]
    for (orbital, occupancy, state) in c
        nrorbital = nonrelorbital(orbital)
        nridx = findfirst(isequal(nrorbital), nrorbitals)
        if isnothing(nridx)
            push!(nrorbitals, nrorbital)
            push!(nroccupancies, occupancy)
            push!(nrstates, state)
        else
            nrstates[nridx] == state || throw(ArgumentError("Mismatching states for $(nrorbital). Previously found $(nrstates[nridx]), but $(orbital) has $(state)"))
            nroccupancies[nridx] += occupancy
        end
    end
    Configuration(nrorbitals, nroccupancies, nrstates)
end

"""
    multiplicity(::Configuration)

Calculates the number of Slater determinants corresponding to the configuration.
"""
multiplicity(c::Configuration) = prod(binomial.(degeneracy.(c.orbitals), c.occupancy))

export Configuration, @c_str, @rc_str, @sc_str, @rsc_str, @scs_str, @rscs_str, issimilar,
    num_electrons, core, peel, active, inactive, bound, continuum, parity, ⊗, @rcs_str,
    SpinConfiguration, spin_configurations, substitutions, close!,
    nonrelconfiguration, relconfigurations


text:
testOnAllContexts("Comparers") do ctx
    println("\tComparer")
    arr = [1,3,2]
    sorted = sort!(arr, lt = Comparer())
    @test sorted == [1,2,3]
    sorted = sort!(arr, lt = Comparer((x,y) -> >(x, y)))
    @test sorted == [3,2,1]
    sorted = sort!(arr, lt = Comparer(Comparer()))
    @test sorted == [1,2,3]

    println("\tCalculateComparer")
    arr = [1,3,2,1]
    calculations = 0
    factory = x -> (calculations += 1; -x)
    sorted = sort!(arr, lt = CalculateComparer(factory))
    @test sorted == [3,2,1,1]
    sorted = sort!(arr, lt = CalculateComparer(factory, (x,y) -> >(x, y)))
    @test sorted == [1,1,2,3]
    sorted = sort!(arr, lt = CalculateComparer(factory, Comparer()))
    @test sorted == [3,2,1,1]
    @test calculations == 3 * 3

    println("\tFastComparer")
    arr = [1,2,2,1,2,1,1,1,1,2,1]
    comparisons = 0
    compare = (x,y) -> (comparisons += 1; <(x, y))
    sorted = sort!(arr, lt = Comparer(compare))
    @test sorted == [1,1,1,1,1,1,1,2,2,2,2]
    @test comparisons == 28

    fc = FastComparer(compare)

    comparisons = 0
    sorted = sort!(arr, lt = fc)
    @test sorted == [1,1,1,1,1,1,1,2,2,2,2]
    @test comparisons == 6

    comparisons = 0
    sorted = sort!(arr, lt = fc)
    @test sorted == [1,1,1,1,1,1,1,2,2,2,2]
    @test comparisons == 0

    comparisons = 0
    reset!(fc)
    sorted = sort!(arr, lt = fc)
    @test sorted == [1,1,1,1,1,1,1,2,2,2,2]
    @test comparisons == 6

    println("\tArrayComparer")
    ac = ArrayComparer()

    println("\t- throws when attempt to compare non arrays")
    @test_throws TypeError ac(1, 2)

    println("\t- compares arrays of same size")
    a1 = [1, 2, 3]
    a2 = [2, 3, 4]
    @test ac(a1, a2) == true
    @test ac(a2, a1) == false
    @test ac(a1, a1) == false
    @test ac(a2, a2) == false

    println("\t- compares arrays of different sizes")
    a1 = [1, 2, 3, 4]
    a2 = [4, 5]
    @test ac(a1, a2) == true
    @test ac(a2, a1) == false
    @test ac(a1, a1) == false
    @test ac(a2, a2) == false

    println("\t- compares same arrays")
    a1 = [1, 2, 3]
    a2 = [1, 2, 3]
    @test ac(a1, a2) == false
    @test ac(a2, a1) == false
    @test ac(a1, a1) == false
    @test ac(a2, a2) == false
end

ext: jl

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
using Debugger: eval_code

# Simple evaling of function argument
function evalfoo1(x,y)
    x+y
end
frame = JuliaInterpreter.enter_call_expr(:($(evalfoo1)(1,2)))
res = eval_code(frame, "x")

@test res == 1

res = eval_code(frame, "y")
@test res == 2

# Evaling with sparams
function evalsparams(x::T) where T
    x
end
frame = JuliaInterpreter.enter_call_expr(:($(evalsparams)(1)))
res = eval_code(frame, "x")
@test res == 1
eval_code(frame, "x = 3")
res = eval_code(frame, "x")
@test res == 3
res = eval_code(frame, "T")
@test res == Int
eval_code(frame, "T = Float32")
res = eval_code(frame, "T")
@test res == Float32

# Evaling with keywords
evalkw(x; bar=true) = x
frame = @make_frame evalkw(2)
res = eval_code(frame, "x")
@test res == 2

# Evaling with symbols
evalsym() = (x = :foo)
frame = @make_frame evalsym()
res = eval_code(frame, "x")
@test res == :foo

frame = Debugger.@make_frame evalsym()
state = dummy_state(frame)
Debugger.add_watch_entry!(state, "x")
watch_str = sprint(Debugger.show_watch_list, state)
@test occursin("1] x: :foo", watch_str)
Debugger.clear_watch_list!(state)


text:
if !isdefined(Main, :CustomStencil)
    include("../module/CustomStencil.jl")
    using Main.CustomStencil
end
include("misc.jl")
include("cpu_stencils.jl")

## Stencil Definition
coefs = Float64.([0.25; 0.125; 0.06])
dense_stencil = @def_stencil_expression(
            @sum(i,-2,2,
                @sum(j,-2,2,
                    @sum(k,-2,2, c[max(abs(i),abs(j),abs(k))]*D[x+i,y+j,z+k]))))
st_def = CreateStencilDefinition(dense_stencil, coefs)
st_inst = NewStencilInstance(st_def)

## Input Data size Definition
radius = 2
nx = 8
ny = 8
nz = 8
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

t_steps = 4
## GPU Implementation
original_data = CreateData(dx,dy,dz)
data = PadData(radius, original_data)
dout = ApplyStencil(st_inst, data, t_steps)
gpu_out = Array(dout)


text:
import LinearAlgebra.dot
import GeometryBasics.Point
import Random.GLOBAL_RNG

export random_vector_matrix, perlin_noise

"""
    random_vector_matrix([rng,] rows, cols)

Produce a matrix with the given number of rows and columns, in which every
entry is a random unit vector. If provided, the given random number generator
is used to determine the vectors.
"""
function random_vector_matrix(rng, rows, cols)
	map(frac_rotation.(rand(rng, rows, cols))) do rot
		rot * unitvec
	end
end

function random_vector_matrix(rows, cols)
    random_vector_matrix(GLOBAL_RNG, rows, cols)
end

function smoothstep(v1, v2, w)
	w = clamp(w, 0., 1.)
	w = 6. * w^5 - 15. * w^4 + 10. * w^3
	(1. - w) * v1 + w * v2
end

"""
    perlin_noise(control_points, point)

Compute the noise gradient for the given point, given a grid of control points.
"""
function perlin_noise(control_points, point::Point)
    rows, cols = size(control_points)
    cell = convert.(Int, (floor.(point)))
    xi, yi = cell

    offset = point - cell
    xf, yf = offset

    n00 = dot(control_points[xi % rows + 1, yi % cols + 1],
        Point(xf, yf))
    n01 = dot(control_points[xi % rows + 1, (yi + 1) % cols + 1],
        Point(xf, yf - 1))
    n10 = dot(control_points[(xi + 1) % rows + 1, yi % cols + 1],
        Point(xf - 1, yf))
    n11 = dot(control_points[(xi + 1) % rows + 1, (yi + 1) % cols + 1],
        Point(xf - 1, yf - 1))

    smoothstep(
        smoothstep(n00, n01, yf),
        smoothstep(n10, n11, yf),
        xf
    )
end


text:
"""
    arrows(points, directions; kwargs...)
    arrows(x, y, u, v)
    arrows(x::AbstractVector, y::AbstractVector, u::AbstractMatrix, v::AbstractMatrix)
    arrows(x, y, z, u, v, w)

Plots arrows at the specified points with the specified components.
`u` and `v` are interpreted as vector components (`u` being the x
and `v` being the y), and the vectors are plotted with the tails at
`x`, `y`.

If `x, y, u, v` are `<: AbstractVector`, then each 'row' is plotted
as a single vector.

If `u, v` are `<: AbstractMatrix`, then `x` and `y` are interpreted as
specifications for a grid, and `u, v` are plotted as arrows along the
grid.

`arrows` can also work in three dimensions.

## Attributes
$(ATTRIBUTES)
"""
@recipe(Arrows, points, directions) do scene
    theme = Attributes(
        arrowhead = automatic,
        arrowtail = nothing,
        linecolor = :black,
        linewidth = 1,
        arrowsize = 0.3,
        linestyle = nothing,
        scale = Vec3f0(1),
        normalize = false,
        lengthscale = 1.0f0,
        colormap = :viridis
    )
    # connect arrow + linecolor by default
    get!(theme, :arrowcolor, theme[:linecolor])
    theme
end

# For the matlab/matplotlib users
const quiver = arrows
const quiver! = arrows!
export quiver, quiver!

arrow_head(N, marker) = marker
arrow_head(N, marker::Automatic) = N == 2 ? '▲' : Pyramid(Point3f0(0, 0, -0.5), 1f0, 1f0)

scatterfun(N) = N == 2 ? scatter! : meshscatter!

convert_arguments(::Type{<: Arrows}, x, y, u, v) = (Point2f0.(x, y), Vec2f0.(u, v))
function convert_arguments(::Type{<: Arrows}, x::AbstractVector, y::AbstractVector, u::AbstractMatrix, v::AbstractMatrix)
    (vec(Point2f0.(x, y')), vec(Vec2f0.(u, v)))
end
convert_arguments(::Type{<: Arrows}, x, y, z, u, v, w) = (Point3f0.(x, y, z), Vec3f0.(u, v, w))

function plot!(arrowplot::Arrows{<: Tuple{AbstractVector{<: Point{N, T}}, V}}) where {N, T, V}
    @extract arrowplot (points, directions, lengthscale, arrowhead, arrowsize, arrowcolor, colormap)
    headstart = lift(points, directions, lengthscale) do points, directions, s
        map(points, directions) do p1, dir
            dir = arrowplot[:normalize][] ? StaticArrays.normalize(dir) : dir
            Point{N, Float32}(p1) => Point{N, Float32}(p1 .+ (dir .* Float32(s)))
        end
    end
    linesegments!(
        arrowplot, headstart,
        color = arrowplot[:linecolor], linewidth = arrowplot[:linewidth],
        linestyle = arrowplot[:linestyle], colormap = colormap,
    )
    scatterfun(N)(
        arrowplot,
        lift(x-> last.(x), headstart),
        marker = lift(x-> arrow_head(N, x), arrowhead), markersize = arrowsize,
        color = arrowcolor, rotations = directions,  strokewidth = 0.0, colormap = colormap,
    )
end


text:
#Run this in Escher to visualize a WheelGraph

using LightGraphs
using NetworkViz
using ThreeJS
using Colors

main(window) = begin
    num = Signal(10)
    width = Signal(1)
    nodesize = Signal(0.2)
    toggle = Signal(false)

    push!(window.assets,("ThreeJS","threejs"))
    push!(window.assets,"widgets")
        vbox(
        h1("WheelGraph Example"),
        vskip(2em),
        vbox(
            "Number of Nodes",
            slider(10:100) >>> num
        ),
        vskip(2em),
        vbox(
            "2D/3D",
            hskip(2em),
            togglebutton() >>> toggle
        ),
        vskip(2em),

        map(toggle,num) do t,n
            g = WheelGraph(n)
            c = Color[parse(Colorant,"#00004d") for i in 1:nv(g)]
            n = NodeProperty(c,0.2,1)
            e = EdgeProperty("#ff3333",1)
            drawGraph(g,node=n,edge=e,z=t)
        end
        ) |> pad(2em)
end


text:
using Documenter
using VectorAutoregressions
using Literate

EXAMPLE = dirname(dirname(pathof(VectorAutoregressions)))*"/examples/lutkepohl_example.jl"
OUTPUT = joinpath(@__DIR__, "src/generated")
codefence = "```@repl lutkepohl_example" => "```"
Literate.markdown(EXAMPLE, OUTPUT, codefence = codefence, documenter=true)

makedocs(
    modules=[VectorAutoregressions],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Example" => "generated/lutkepohl_example.md",
        "API"  => "api.md"
    ],
    sitename="VectorAutoregressions.jl",
    authors="Stephen Nicar",
)

deploydocs(
    repo="github.com/stephenbnicar/VectorAutoregressions.jl",
)


text:
@testset "SVD" begin
    for n in [1,3,10]

m = n + rand(0:5)
A = rand(m, n)

r = polarfact(A, alg = :svd);

# Test unitary matrix U

U = r.U
H = r.H

@test U'*U ≈ Matrix(I,n,n) atol=1e-7

# Test Hermitian positive semifefinite matrix H

@test issymmetric(H)

for i in eigvals(H)
    @test i >= 0.
end

@test A ≈ U*H

    end
end


text:
function forward(X, Y, CH::ConditionalLayerHINT; logdet = nothing, x_lane = false)
    isnothing(logdet) ? logdet = (CH.logdet && ~CH.is_reversed) : logdet = logdet

    # Y-lane
    ~isnothing(CH.C_Y) ? (Yp = CH.C_Y.forward(Y)) : (Yp = copy(Y))
    logdet ? (Zy, logdet2) = CH.CL_Y.forward(Yp) : Zy = CH.CL_Y.forward(Yp)

    # X-lane: coupling layer
    ~isnothing(CH.C_X) ? (Xp = CH.C_X.forward(X)) : (Xp = copy(X))
    logdet ? (X, logdet1) = CH.CL_X.forward(Xp) : X = CH.CL_X.forward(Xp)

    # X-lane: conditional layer
    logdet ? (Zx, logdet3) = CH.CL_YX.forward(Yp, X)[2:3] : Zx = CH.CL_YX.forward(Yp, X)[2]

    if x_lane == false
        logdet ? (return Zx, Zy, logdet1 + logdet2 + logdet3) : (return Zx, Zy)
    else
        logdet ? (return Zx, Zy, logdet1 + logdet3) : (return Zx, Zy)
    end
end

function inverse(Zx, Zy, CH::ConditionalLayerHINT; logdet = nothing, x_lane = false)
    isnothing(logdet) ? logdet = (CH.logdet && CH.is_reversed) : logdet = logdet

    # Y-lane
    logdet ? (Yp, logdet1) = CH.CL_Y.inverse(Zy; logdet = true) :
    Yp = CH.CL_Y.inverse(Zy; logdet = false)
    ~isnothing(CH.C_Y) ? (Y = CH.C_Y.inverse(Yp)) : (Y = copy(Yp))

    # X-lane: conditional layer
    YZ = tensor_cat(Yp, Zx)
    logdet ? (X, logdet2) = CH.CL_YX.inverse(Yp, Zx; logdet = true)[2:3] :
    X = CH.CL_YX.inverse(Yp, Zx)[2]

    # X-lane: coupling layer
    logdet ? (Xp, logdet3) = CH.CL_X.inverse(X; logdet = true) :
    Xp = CH.CL_X.inverse(X; logdet = false)
    ~isnothing(CH.C_X) ? (X = CH.C_X.inverse(Xp)) : (X = copy(Xp))

    if x_lane == false
        logdet ? (return X, Y, logdet1 + logdet2 + logdet3) : (return X, Y)
    else
        logdet ? (return X, Y, logdet2 + logdet3) : (return X, Y)
    end
end

function inverse(Y, H::CouplingLayerHINT; scale = 1, permute = nothing, logdet = nothing)
    isnothing(logdet) ? logdet = (H.logdet && H.is_reversed) : logdet = logdet
    isnothing(permute) ? permute = H.permute : permute = permute

    # Permutation
    permute == "both" && (Y = H.C.forward(Y))
    Ya, Yb = tensor_split(Y)

    # Check for recursion
    recursive = false
    if typeof(Y) <: AbstractArray{Float32,4} && size(Y, 3) > 4
        recursive = true
    elseif typeof(Y) <: AbstractArray{Float32,5} && size(Y, 4) > 4
        recursive = true
    end

    # Coupling layer
    if recursive
        Xa, logdet1 = inverse(Ya, H; scale = scale + 1, permute = "none")
        if logdet
            Y_temp, logdet2 = H.CL[scale].inverse(Xa, Yb; logdet = true)[[2, 3]]
        else
            Y_temp = H.CL[scale].inverse(Xa, Yb)[2]
            logdet2 = 0.0f0
        end
        Xb, logdet3 = inverse(Y_temp, H; scale = scale + 1, permute = "none")
        logdet_full = logdet1 + logdet2 + logdet3
    else
        Xa = copy(Ya)
        if logdet
            Xb, logdet_full = H.CL[scale].inverse(Ya, Yb; logdet = true)[[2, 3]]
        else
            Xb = H.CL[scale].inverse(Ya, Yb)[2]
            logdet_full = 0.0f0
        end
    end

    # Initial permutation
    permute == "lower" && (Xb = H.C.inverse(Xb))
    X = tensor_cat(Xa, Xb)
    if permute == "full" || permute == "both"
        X = H.C.inverse(X)
    end

    if scale == 1
        logdet ? (return X, logdet_full) : (return X)
    else
        return X, logdet_full
    end
end

function ResidualBlock(
    n_in,
    n_hidden;
    k1 = 3,
    k2 = 3,
    p1 = 1,
    p2 = 1,
    s1 = 1,
    s2 = 1,
    fan = false,
    ndims = 2,
)

    k1 = Tuple(k1 for i = 1:ndims)
    k2 = Tuple(k2 for i = 1:ndims)
    # Initialize weights
    W1 = Parameter(glorot_uniform(k1..., n_in, n_hidden))
    W2 = Parameter(glorot_uniform(k2..., n_hidden, n_hidden))
    W3 = Parameter(1.0f-2 * glorot_uniform(k1..., 2 * n_in, n_hidden))
    b1 = Parameter(zeros(Float32, n_hidden))
    b2 = Parameter(zeros(Float32, n_hidden))

    return ResidualBlock(W1, W2, W3, b1, b2, fan, (s1, s2), (p1, p2))
end

function forward(X, Y, CH::NetworkConditionalHINT; logdet = nothing, x_lane = false)
    isnothing(logdet) ? logdet = (CH.logdet && ~CH.is_reversed) : logdet = logdet

    depth = length(CH.CL)
    logdet_ = 0.0f0
    logdet_x = 0.0f0
    for j = 1:depth
        logdet ? (X_, logdet1) = CH.AN_X[j].forward(X) : X_ = CH.AN_X[j].forward(X)
        logdet ? (Y_, logdet2) = CH.AN_Y[j].forward(Y) : Y_ = CH.AN_Y[j].forward(Y)
        logdet ? (X, Y, logdet3) = CH.CL[j].forward(X_, Y_; x_lane = x_lane) :
        (X, Y) = CH.CL[j].forward(X_, Y_)
        logdet && (logdet_ += (logdet1 + logdet2 + logdet3))
        logdet && (logdet_x += (logdet1 + logdet3))
    end
    if x_lane == false
        logdet ? (return X, Y, logdet_) : (return X, Y)
    else
        logdet ? (return X, Y, logdet_x) : (return X, Y)
    end
end

function inverse(Zx, Zy, CH::NetworkConditionalHINT; logdet = nothing, x_lane = false)
    isnothing(logdet) ? logdet = (CH.logdet && CH.is_reversed) : logdet = logdet

    depth = length(CH.CL)
    logdet_ = 0.0f0
    logdet_x = 0.0f0
    for j = depth:-1:1
        logdet ?
        (Zx_, Zy_, logdet1) = CH.CL[j].inverse(Zx, Zy; logdet = true, x_lane = x_lane) :
        (Zx_, Zy_) = CH.CL[j].inverse(Zx, Zy; logdet = false)
        logdet ? (Zy, logdet2) = CH.AN_Y[j].inverse(Zy_; logdet = true) :
        Zy = CH.AN_Y[j].inverse(Zy_; logdet = false)
        logdet ? (Zx, logdet3) = CH.AN_X[j].inverse(Zx_; logdet = true) :
        Zx = CH.AN_X[j].inverse(Zx_; logdet = false)
        logdet && (logdet_ += (logdet1 + logdet2 + logdet3))
        logdet && (logdet_x += (logdet1 + logdet3))
    end
    if x_lane == false
        logdet ? (return Zx, Zy, logdet_) : (return Zx, Zy)
    else
        logdet ? (return Zx, Zy, logdet_x) : (return Zx, Zy)
    end
end


text:
"""
    Perceptron.jl file

Implement the BetaML.Perceptron module

`?BetaML.Perceptron` for documentation

- [Importable source code (most up-to-date version)](https://github.com/sylvaticus/BetaML.jl/blob/master/src/Perceptron.jl) - [Julia Package](https://github.com/sylvaticus/BetaML.jl)
- [Demonstrative static notebook](https://github.com/sylvaticus/lmlj.jl/blob/master/notebooks/Perceptron.ipynb)
- [Demonstrative live notebook](https://mybinder.org/v2/gh/sylvaticus/BetaML.jl/master?filepath=notebooks%2FPerceptron.ipynb) (temporary personal online computational environment on myBinder) - it can takes minutes to start with!
- Theory based on [MITx 6.86x - Machine Learning with Python: from Linear Models to Deep Learning](https://github.com/sylvaticus/MITx_6.86x) ([Unit 3](https://github.com/sylvaticus/MITx_6.86x/blob/master/Unit%2003%20-%20Neural%20networks/Unit%2003%20-%20Neural%20networks.md))
- New to Julia? [A concise Julia tutorial](https://github.com/sylvaticus/juliatutorial) - [Julia Quick Syntax Reference book](https://julia-book.com)

"""

"""
    Perceptron module

Provide linear and kernel classifiers.

See a [runnable example on myBinder](https://mybinder.org/v2/gh/sylvaticus/BetaML.jl/master?filepath=notebooks%2FPerceptron.ipynb)

- [`perceptron`](@ref): Train data using the classical perceptron
- [`kernelPerceptron`](@ref): Train data using the kernel perceptron
- [`pegasos`](@ref): Train data using the pegasos algorithm
- [`predict`](@ref): Predict data using parameters from one of the above algorithms

All algorithms are multiclass, with `perceptron` and `pegasos` employing a one-vs-all strategy, while `kernelPerceptron` employs a _one-vs-one_ approach, and return a "probability" for each class in term of a dictionary for each record. Use `mode(ŷ)` to return a single class prediction per record.

The binary equivalent algorithms, accepting only `{-1,+1}` labels, are available as `peceptronBinary`, `kernelPerceptronBinary` and `pegasosBinary`. They are slighly faster as they don't need to be wrapped in the multi-class equivalent and return a more informative output.

The multi-class versions are available in the MLJ framework as `PerceptronClassifier`,`KernelPerceptronClassifier` and `PegasosClassifier` respectivly.
"""
module Perceptron

using LinearAlgebra, Random, ProgressMeter, Reexport, CategoricalArrays

using ForceImport
@force using ..Api
@force using ..Utils

export perceptron, perceptronBinary, kernelPerceptron, kernelPerceptronBinary, pegasos, pegasosBinary, predict

# Todo (breaking): change the API so that a struct kernelPerceptronModel or linearModel is reported
# instead of a named tuple
# linearModel <: betaSupervisedModel <: betaModel

"""
    perceptron(x,y;θ,θ₀,T,nMsgs,shuffle,forceOrigin,returnMeanHyperplane)

Train the multiclass classifier "perceptron" algorithm  based on x and y (labels).

The perceptron is a _linear_ classifier. Multiclass is supported using a one-vs-all approach.

# Parameters:
* `x`:           Feature matrix of the training data (n × d)
* `y`:           Associated labels of the training data, can be in any format (string, integers..)
* `θ`:           Initial value of the weights (parameter) [def: `zeros(d)`]
* `θ₀`:          Initial value of the weight (parameter) associated to the constant
                 term [def: `0`]
* `T`:           Maximum number of iterations across the whole set (if the set
                 is not fully classified earlier) [def: 1000]
* `nMsg`:        Maximum number of messages to show if all iterations are done [def: `0`]
* `shuffle`:     Whether to randomly shuffle the data at each iteration [def: `false`]
* `forceOrigin`: Whether to force `θ₀` to remain zero [def: `false`]
* `returnMeanHyperplane`: Whether to return the average hyperplane coefficients instead of the final ones  [def: `false`]
* `rng`:         Random Number Generator (see [`FIXEDSEED`](@ref)) [deafult: `Random.GLOBAL_RNG`]

# Return a named tuple with:
* `θ`:          The weights of the classifier
* `θ₀`:         The weight of the classifier associated to the constant term
* `classes`:    The classes (unique values) of y

# Notes:
* The trained parameters can then be used to make predictions using the function `predict()`.
* This model is available in the MLJ framework as the `PerceptronClassifier`

# Example:
```jldoctest
julia> model = perceptron([1.1 2.1; 5.3 4.2; 1.8 1.7], [-1,1,-1])
julia> ŷ     = predict([2.1 3.1; 7.3 5.2], model.θ, model.θ₀, model.classes)
```
"""
function perceptron(x::AbstractMatrix, y::AbstractVector; θ=zeros(size(x,2)),θ₀=0.0, T=1000, nMsgs=0, shuffle=false, forceOrigin=false, returnMeanHyperplane=false, rng = Random.GLOBAL_RNG)
    yclasses = unique(y)
    nCl      = length(yclasses)
    if nCl == 2
        outθ        = Array{Vector{Float64},1}(undef,1)
        outθ₀       = Array{Float64,1}(undef,1)
    else
        outθ        = Array{Vector{Float64},1}(undef,nCl)
        outθ₀       = Array{Float64,1}(undef,nCl)
    end
    for (i,c) in enumerate(yclasses)
        ybin = ((y .== c) .*2 .-1)  # conversion to -1/+1
        outBinary = perceptronBinary(x, ybin; θ=θ,θ₀=θ₀, T=T, nMsgs=nMsgs, shuffle=shuffle, forceOrigin=forceOrigin, rng=rng)
        if returnMeanHyperplane
            outθ[i]  = outBinary.avgθ
            outθ₀[i] = outBinary.avgθ₀
        else
            outθ[i]  = outBinary.θ
            outθ₀[i] = outBinary.θ₀
        end
        if nCl == 2
            break    # if there are only two classes we do compute only one passage, as A vs B would be the same as B vs A
        end
    end
    return (θ=outθ,θ₀=outθ₀,classes=yclasses)
end

"""
    perceptronBinary(x,y;θ,θ₀,T,nMsgs,shuffle,forceOrigin)

Train the binary classifier "perceptron" algorithm based on x and y (labels)

# Parameters:
* `x`:           Feature matrix of the training data (n × d)
* `y`:           Associated labels of the training data, in the format of ⨦ 1
* `θ`:           Initial value of the weights (parameter) [def: `zeros(d)`]
* `θ₀`:          Initial value of the weight (parameter) associated to the constant
                 term [def: `0`]
* `T`:           Maximum number of iterations across the whole set (if the set
                 is not fully classified earlier) [def: 1000]
* `nMsg`:        Maximum number of messages to show if all iterations are done
* `shuffle`:     Whether to randomly shuffle the data at each iteration [def: `false`]
* `forceOrigin`: Whether to force `θ₀` to remain zero [def: `false`]
* `rng`:         Random Number Generator (see [`FIXEDSEED`](@ref)) [deafult: `Random.GLOBAL_RNG`]

# Return a named tuple with:
* `θ`:          The final weights of the classifier
* `θ₀`:         The final weight of the classifier associated to the constant term
* `avgθ`:       The average weights of the classifier
* `avgθ₀`:      The average weight of the classifier associated to the constant term
* `errors`:     The number of errors in the last iteration
* `besterrors`: The minimum number of errors in classifying the data ever reached
* `iterations`: The actual number of iterations performed
* `separated`:  Weather the data has been successfully separated

# Notes:
* The trained parameters can then be used to make predictions using the function `predict()`.

# Example:
```jldoctest
julia> model = perceptronBinary([1.1 2.1; 5.3 4.2; 1.8 1.7], [-1,1,-1])
```
"""
function perceptronBinary(x, y; θ=zeros(size(x,2)),θ₀=0.0, T=1000, nMsgs=10, shuffle=false, forceOrigin=false, rng = Random.GLOBAL_RNG)
   if nMsgs != 0
       @codeLocation
       println("***\n*** Training perceptron for maximum $T iterations. Random shuffle: $shuffle")
   end
   x = makeMatrix(x)
   (n,d) = size(x)
   bestϵ = Inf
   lastϵ = Inf
   if forceOrigin θ₀ = 0.0; end
   sumθ = θ; sumθ₀ = θ₀
   @showprogress 1 "Training Perceptron..." for t in 1:T
       ϵ = 0
       if shuffle
          # random shuffle x and y
          ridx = Random.shuffle(rng, 1:size(x)[1])
          x = x[ridx, :]
          y = y[ridx]
       end
       for i in 1:n
           if y[i]*(θ' * x[i,:] + θ₀) <= eps()
               θ  = θ + y[i] * x[i,:]
               θ₀ = forceOrigin ? 0.0 : θ₀ + y[i]
               sumθ += θ; sumθ₀ += θ₀
               ϵ += 1
           end
       end
       if (ϵ == 0)
           if nMsgs != 0
               println("*** Avg. error after epoch $t : $(ϵ/size(x)[1]) (all elements of the set has been correctly classified")
           end
           return (θ=θ,θ₀=θ₀,avgθ=sumθ/(n*T),avgθ₀=sumθ₀/(n*T),errors=0,besterrors=0,iterations=t,separated=true)
       elseif ϵ < bestϵ
           bestϵ = ϵ
       end
       lastϵ = ϵ
       if nMsgs != 0 && (t % ceil(T/nMsgs) == 0 || t == 1 || t == T)
         println("Avg. error after iteration $t : $(ϵ/size(x)[1])")
       end
   end
   return  (θ=θ,θ₀=θ₀,avgθ=sumθ/(n*T),avgθ₀=sumθ₀/(n*T),errors=lastϵ,besterrors=bestϵ,iterations=T,separated=false)
end

"""
   kernelPerceptron(x,y;K,T,α,nMsgs,shuffle)

Train a multiclass kernel classifier "perceptron" algorithm based on x and y.

`kernelPerceptron` is a (potentially) non-linear perceptron-style classifier employing user-defined kernel funcions. Multiclass is supported using a one-vs-one approach.

# Parameters:
* `x`:        Feature matrix of the training data (n × d)
* `y`:        Associated labels of the training data, in the format of ⨦ 1
* `K`:        Kernel function to employ. See `?radialKernel` or `?polynomialKernel`for details or check `?BetaML.Utils` to verify if other kernels are defined (you can alsways define your own kernel) [def: [`radialKernel`](@ref)]
* `T`:        Maximum number of iterations (aka "epochs") across the whole set (if the set is not fully classified earlier) [def: 100]
* `α`:        Initial distribution of the errors [def: `zeros(length(y))`]
* `nMsg`:     Maximum number of messages to show if all iterations are done [def: `0`]
* `shuffle`:  Whether to randomly shuffle the data at each iteration [def: `false`]
* `rng`:      Random Number Generator (see [`FIXEDSEED`](@ref)) [deafult: `Random.GLOBAL_RNG`]

# Return a named tuple with:
* `x`: The x data (eventually shuffled if `shuffle=true`)
* `y`: The label
* `α`: The errors associated to each record
* `classes`: The labels classes encountered in the training

# Notes:
* The trained model can then be used to make predictions using the function `predict()`.
* This model is available in the MLJ framework as the `KernelPerceptronClassifier`

# Example:
```jldoctest
julia> model  = kernelPerceptron([1.1 2.1; 5.3 4.2; 1.8 1.7], [-1,1,-1])
julia> ŷtrain = Perceptron.predict(xtrain,model.x,model.y,model.α, model.classes,K=model.K)
julia> ϵtrain = error(ytrain, mode(ŷtrain))
```
"""
function kernelPerceptron(x, y; K=radialKernel, T=100, α=zeros(Int64,length(y)), nMsgs=0, shuffle=false, rng = Random.GLOBAL_RNG)
    x         = makeMatrix(x)
    yclasses  = unique(y)
    nCl       = length(yclasses)
    nModels   = Int((nCl  * (nCl - 1)) / 2)
    (n,d) = size(x)
    outX = Array{typeof(x),1}(undef,nModels)
    outY = Array{Array{Int64,1},1}(undef,nModels)
    outα = Array{Array{Int64,1},1}(undef,nModels)

    modelCounter = 1
    for (i,c) in enumerate(yclasses)
        for (i2,c2) in enumerate(yclasses)
            if i2 <= i continue end
            ids = ( (y .== c) .| (y .== c2) )
            thisx = x[ids,:]
            thisy = y[ids]
            thisα = α[ids]
            ybin = ((thisy .== c) .*2 .-1)  # conversion to +1 (if c) or -1 (if c2)
            outBinary = kernelPerceptronBinary(thisx, ybin; K=K, T=T, α=thisα, nMsgs=nMsgs, shuffle=shuffle, rng = rng)
            outX[modelCounter] = outBinary.x
            outY[modelCounter] = outBinary.y
            outα[modelCounter] = outBinary.α
            modelCounter += 1
        end
    end
    return (x=outX,y=outY,α=outα,classes=yclasses,K=K)
end

"""
   kernelPerceptronBinary(x,y;K,T,α,nMsgs,shuffle)

Train a multiclass kernel classifier "perceptron" algorithm based on x and y

# Parameters:
* `x`:        Feature matrix of the training data (n × d)
* `y`:        Associated labels of the training data, in the format of ⨦ 1
* `K`:        Kernel function to employ. See `?radialKernel` or `?polynomialKernel`for details or check `?BetaML.Utils` to verify if other kernels are defined (you can alsways define your own kernel) [def: [`radialKernel`](@ref)]
* `T`:        Maximum number of iterations across the whole set (if the set is not fully classified earlier) [def: 1000]
* `α`:        Initial distribution of the errors [def: `zeros(length(y))`]
* `nMsg`:     Maximum number of messages to show if all iterations are done
* `shuffle`:  Whether to randomly shuffle the data at each iteration [def: `false`]
* `rng`:      Random Number Generator (see [`FIXEDSEED`](@ref)) [deafult: `Random.GLOBAL_RNG`]

# Return a named tuple with:
* `x`: the x data (eventually shuffled if `shuffle=true`)
* `y`: the label
* `α`: the errors associated to each record
* `errors`: the number of errors in the last iteration
* `besterrors`: the minimum number of errors in classifying the data ever reached
* `iterations`: the actual number of iterations performed
* `separated`: a flag if the data has been successfully separated

# Notes:
* The trained data can then be used to make predictions using the function `predict()`. **If the option `shuffle` has been used, it is important to use there the returned (x,y,α) as these would have been shuffle compared with the original (x,y)**.
* Please see @kernelPerceptron for a multi-class version
# Example:
```jldoctest
julia> model = kernelPerceptronBinary([1.1 2.1; 5.3 4.2; 1.8 1.7], [-1,1,-1])
```
"""
function kernelPerceptronBinary(x, y; K=radialKernel, T=1000, α=zeros(Int64,length(y)), nMsgs=10, shuffle=false, rng = Random.GLOBAL_RNG)
    if nMsgs != 0
        @codeLocation
        println("***\n*** Training kernel perceptron for maximum $T iterations. Random shuffle: $shuffle")
    end
    x = makeMatrix(x)
    (n,d) = size(x)
    bestϵ = Inf
    lastϵ = Inf
    @showprogress 1 "Training Kernel Perceptron..." for t in 1:T
        ϵ = 0
        if shuffle
           # random shuffle x, y and alpha
           ridx = Random.shuffle(rng, 1:size(x)[1])
           x = x[ridx, :]
           y = y[ridx]
           α = α[ridx]
        end
        for i in 1:n
            if y[i]*sum([α[j]*y[j]*K(x[j,:],x[i,:]) for j in 1:n]) <= 0 + eps()
                α[i] += 1
                ϵ += 1
            end
        end
        if (ϵ == 0)
            if nMsgs != 0
                println("*** Avg. error after epoch $t : $(ϵ/size(x)[1]) (all elements of the set has been correctly classified")
            end
            return (x=x,y=y,α=α,errors=0,besterrors=0,iterations=t,separated=true,K=K)
        elseif ϵ < bestϵ
            bestϵ = ϵ
        end
        lastϵ = ϵ
        if nMsgs != 0 && (t % ceil(T/nMsgs) == 0 || t == 1 || t == T)
          println("Avg. error after iteration $t : $(ϵ/size(x)[1])")
        end
    end
    return  (x=x,y=y,α=α,errors=lastϵ,besterrors=bestϵ,iterations=T,separated=false)
end

"""
    pegasos(x,y;θ,θ₀,λ,η,T,nMsgs,shuffle,forceOrigin,returnMeanHyperplane)

Train the multiclass classifier "pegasos" algorithm according to x (features) and y (labels)

Pegasos is a _linear_, gradient-based classifier. Multiclass is supported using a one-vs-all approach.

# Parameters:
* `x`:           Feature matrix of the training data (n × d)
* `y`:           Associated labels of the training data, can be in any format (string, integers..)
* `θ`:           Initial value of the weights (parameter) [def: `zeros(d)`]
* `θ₀`:          Initial value of the weight (parameter) associated to the constant term [def: `0`]
* `λ`:           Multiplicative term of the learning rate
* `η`:           Learning rate [def: (t -> 1/sqrt(t))]
* `T`:           Maximum number of iterations across the whole set (if the set is not fully classified earlier) [def: 1000]
* `nMsg`:        Maximum number of messages to show if all iterations are done
* `shuffle`:     Whether to randomly shuffle the data at each iteration [def: `false`]
* `forceOrigin`: Whehter to force `θ₀` to remain zero [def: `false`]
* `returnMeanHyperplane`: Whether to return the average hyperplane coefficients instead of the average ones  [def: `false`]
* `rng`:         Random Number Generator (see [`FIXEDSEED`](@ref)) [deafult: `Random.GLOBAL_RNG`]

# Return a named tuple with:
* `θ`:          The weights of the classifier
* `θ₀`:         The weight of the classifier associated to the constant term
* `classes`:    The classes (unique values) of y

# Notes:
* The trained parameters can then be used to make predictions using the function `predict()`.
* This model is available in the MLJ framework as the `PegasosClassifier`

# Example:
```jldoctest
julia> model = pegasos([1.1 2.1; 5.3 4.2; 1.8 1.7], [-1,1,-1])
julia> ŷ     = predict([2.1 3.1; 7.3 5.2], model.θ, model.θ₀, model.classes)
```
"""
function pegasos(x, y; θ=zeros(size(x,2)),θ₀=0.0, λ=0.5,η= (t -> 1/sqrt(t)), T=1000, nMsgs=0, shuffle=false, forceOrigin=false,returnMeanHyperplane=false, rng = Random.GLOBAL_RNG)
    yclasses = unique(y)
    nCl      = length(yclasses)
    if nCl == 2
        outθ        = Array{Vector{Float64},1}(undef,1)
        outθ₀       = Array{Float64,1}(undef,1)
    else
        outθ        = Array{Vector{Float64},1}(undef,nCl)
        outθ₀       = Array{Float64,1}(undef,nCl)
    end
    for (i,c) in enumerate(yclasses)
        ybin = ((y .== c) .*2 .-1)  # conversion to -1/+1
        outBinary = pegasosBinary(x, ybin; θ=θ,θ₀=θ₀, λ=λ,η=η, T=T, nMsgs=nMsgs, shuffle=shuffle, forceOrigin=forceOrigin, rng=rng)
        if returnMeanHyperplane
            outθ[i]  = outBinary.avgθ
            outθ₀[i] = outBinary.avgθ₀
        else
            outθ[i]  = outBinary.θ
            outθ₀[i] = outBinary.θ₀
        end
        if nCl == 2
            break    # if there are only two classes we do compute only one passage, as A vs B would be the same as B vs A
        end
    end
    return (θ=outθ,θ₀=outθ₀,classes=yclasses)
end

"""
    pegasosBinary(x,y;θ,θ₀,λ,η,T,nMsgs,shuffle,forceOrigin)

Train the peagasos algorithm based on x and y (labels)

# Parameters:
* `x`:           Feature matrix of the training data (n × d)
* `y`:           Associated labels of the training data, in the format of ⨦ 1
* `θ`:           Initial value of the weights (parameter) [def: `zeros(d)`]
* `θ₀`:          Initial value of the weight (parameter) associated to the constant term [def: `0`]
* `λ`:           Multiplicative term of the learning rate
* `η`:           Learning rate [def: (t -> 1/sqrt(t))]
* `T`:           Maximum number of iterations across the whole set (if the set is not fully classified earlier) [def: 1000]
* `nMsg`:        Maximum number of messages to show if all iterations are done
* `shuffle`:    Whether to randomly shuffle the data at each iteration [def: `false`]
* `forceOrigin`: Whether to force `θ₀` to remain zero [def: `false`]

# Return a named tuple with:
* `θ`:          The final weights of the classifier
* `θ₀`:         The final weight of the classifier associated to the constant term
* `avgθ`:       The average weights of the classifier
* `avgθ₀`:      The average weight of the classifier associated to the constant term
* `errors`:     The number of errors in the last iteration
* `besterrors`: The minimum number of errors in classifying the data ever reached
* `iterations`: The actual number of iterations performed
* `separated`:  Weather the data has been successfully separated

# Notes:
* The trained parameters can then be used to make predictions using the function `predict()`.

# Example:
```jldoctest
julia> pegasos([1.1 2.1; 5.3 4.2; 1.8 1.7], [-1,1,-1])
```
"""
function pegasosBinary(x, y; θ=zeros(size(x,2)),θ₀=0.0, λ=0.5,η= (t -> 1/sqrt(t)), T=1000, nMsgs=10, shuffle=false, forceOrigin=false, rng = Random.GLOBAL_RNG)
  if nMsgs != 0
      @codeLocation
      println("***\n*** Training pegasos for maximum $T iterations. Random shuffle: $shuffle")
  end
  x = makeMatrix(x)
  (n,d) = size(x)
  bestϵ = Inf
  lastϵ = Inf
  if forceOrigin θ₀ = 0.0; end
  sumθ = θ; sumθ₀ = θ₀
  @showprogress 1 "Training Pegasos..." for t in 1:T
      ϵ = 0
      ηₜ = η(t)
      if shuffle
         # random shuffle x and y
         ridx = Base.shuffle(rng, 1:size(x)[1])
         x = x[ridx, :]
         y = y[ridx]
      end
      for i in 1:n
          if y[i]*(θ' * x[i,:] + θ₀) <= eps()
              θ  = (1-ηₜ*λ) * θ + ηₜ * y[i] * x[i,:]
              θ₀ = forceOrigin ? 0.0 : θ₀ + ηₜ * y[i]
              sumθ += θ; sumθ₀ += θ₀
              ϵ += 1
          else
              θ  = (1-ηₜ*λ) * θ
          end
      end
      if (ϵ == 0)
          if nMsgs != 0
              println("*** Avg. error after epoch $t : $(ϵ/size(x)[1]) (all elements of the set has been correctly classified")
          end
          return (θ=θ,θ₀=θ₀,avgθ=sumθ/(n*T),avgθ₀=sumθ₀/(n*T),errors=0,besterrors=0,iterations=t,separated=true)
      elseif ϵ < bestϵ
          bestϵ = ϵ
      end
      lastϵ = ϵ
      if nMsgs != 0 && (t % ceil(T/nMsgs) == 0 || t == 1 || t == T)
        println("Avg. error after iteration $t : $(ϵ/size(x)[1])")
      end
  end
  return  (θ=θ,θ₀=θ₀,avgθ=sumθ/(n*T),avgθ₀=sumθ₀/(n*T),errors=lastϵ,besterrors=bestϵ,iterations=T,separated=false)
end

# ------------------------------------------------------------------------------
# Other functions

"""
  predict(x,θ,θ₀)

Predict a binary label {-1,1} given the feature vector and the linear coefficients

# Parameters:
* `x`:        Feature matrix of the training data (n × d)
* `θ`:        The trained parameters
* `θ₀`:       The trained bias barameter [def: `0`]

# Return :
* `y`: Vector of the predicted labels

# Example:
```julia
julia> predict([1.1 2.1; 5.3 4.2; 1.8 1.7], [3.2,1.2])
```
"""
function predict(x,θ,θ₀=0.0)
    x = makeMatrix(x)
    θ = makeColVector(θ)
    (n,d) = size(x)
    d2 = length(θ)
    if (d2 != d) error("x and θ must have the same dimensions."); end
    y = zeros(Int64,n)
    for i in 1:n
        y[i] = (θ' * x[i,:] + θ₀) > eps() ? 1 : -1  # no need to divide by the norm to get the sign!
    end
    return y
end

"""
  predict(x,θ,θ₀,classes)

Predict a multiclass label given the feature vector, the linear coefficients and the classes vector

# Parameters:
* `x`:       Feature matrix of the training data (n × d)
* `θ`:       Vector of the trained parameters for each one-vs-all model (i.e. `model.θ`)
* `θ₀`:      Vector of the trained bias barameter for each one-vs-all model (i.e. `model.θ₀`)
* `classes`: The overal classes encountered in training (i.e. `model.classes`)

# Return :
* `ŷ`: Vector of dictionaries `label=>probability`

 # Notes:
 * Use `mode(ŷ)` if you want a single predicted label per record

# Example:
```julia
julia> model  = perceptron([1.1 2.1; 5.3 4.2; 1.8 1.7], [-1,1,-1])
julia> ŷtrain = predict([10 10; 2.5 2.5],model.θ,model.θ₀, model.classes)
"""
function predict(x,θ::AbstractVector{T},θ₀::AbstractVector{Float64},classes::Vector{Tcl}) where {T<: AbstractVector{Float64},Tcl}
    (n,d) = size(x)
    nCl   = length(classes)
    y     = Array{Dict{Tcl,Float64},1}(undef,n)
    for i in 1:n
        probRaw = Array{Float64,1}(undef,nCl)
        for (c,cl) in enumerate(classes)
            if nCl == 2 && c ==2
                probRaw[2] = - probRaw[1]
            else
                probRaw[c] = (θ[c]' * x[i,:] + θ₀[c])
            end
        end
        prob = softmax(probRaw)
        y[i] = Dict(zip(classes,prob))
    end
    return y
end

"""
  predict(x,xtrain,ytrain,α;K)

Predict a binary label {-1,1} given the feature vector and the training data together with their errors (as trained by a kernel perceptron algorithm)

# Parameters:
* `x`:      Feature matrix of the training data (n × d)
* `xtrain`: The feature vectors used for the training
* `ytrain`: The labels of the training set
* `α`:      The errors associated to each record
* `K`:      The kernel function used for the training and to be used for the prediction [def: [`radialKernel`](@ref)]

# Return :
* `y`: Vector of the predicted labels

# Example:
```julia
julia> predict([1.1 2.1; 5.3 4.2; 1.8 1.7], [3.2,1.2])
```
"""
function predict(x,xtrain,ytrain,α;K=radialKernel)
    x = makeMatrix(x)
    xtrain = makeMatrix(xtrain)
    (n,d) = size(x)
    (ntrain,d2) = size(xtrain)
    if (d2 != d) error("xtrain and x must have the same dimensions."); end
    if ( length(ytrain) != ntrain || length(α) != ntrain) error("xtrain, ytrain and α must all have the same length."); end
    y = zeros(Int64,n)
    for i in 1:n
        y[i] = sum([ α[j] * ytrain[j] * K(x[i,:],xtrain[j,:]) for j in 1:ntrain]) > eps() ? 1 : -1
    end
    return y
 end

 """
   predict(x,xtrain,ytrain,α,classes;K)

 Predict a multiclass label given the new feature vector and a trained kernel perceptron model.

 # Parameters:
 * `x`:      Feature matrix of the training data (n × d)
 * `xtrain`: A vector of the feature matrix used for training each of the one-vs-one class matches (i.e. `model.x`)
 * `ytrain`: A vector of the label vector used for training each of the one-vs-one class matches (i.e. `model.y`)
 * `α`:      A vector of the errors associated to each record (i.e. `model.α`)
 * `classes`: The overal classes encountered in training (i.e. `model.classes`)
 * `K`:      The kernel function used for the training and to be used for the prediction [def: [`radialKernel`](@ref)]

 # Return :
 * `ŷ`: Vector of dictionaries `label=>probability` (warning: it isn't really a probability, it is just the standardized number of matches "won" by this class compared with the other classes)

  # Notes:
  * Use `mode(ŷ)` if you want a single predicted label per record

 # Example:
 ```julia
 julia> model  = kernelPerceptron([1.1 2.1; 5.3 4.2; 1.8 1.7], [-1,1,-1])
 julia> ŷtrain = Perceptron.predict([10 10; 2.2 2.5],model.x,model.y,model.α, model.classes,K=model.K)
 ```
 """
 function predict(x,xtrain,ytrain,α,classes::AbstractVector{Tcl};K=radialKernel) where {Tcl}
     (n,d)   = size(x)
     nCl     = length(classes)
     y       = Array{Dict{Tcl,Float64},1}(undef,n)
     nModels = Int((nCl  * (nCl - 1)) / 2)
     if !(nModels == length(xtrain) == length(ytrain) == length(α)) error("xtrain, ytrain or α have a length not compatible with the number of classes in this model."); end
     x = makeMatrix(x)
     d2 = size(xtrain[1],2)
     if (d2 != d) error("xtrain and x must have the same dimensions."); end
     for i in 1:n
         #countByClass = zeros(Float64,nCl)
         countByClass = zeros(Int64,nCl)
         mCounter = 1
          for (ic,c) in enumerate(classes)
             for (ic2,c2) in enumerate(classes)
                 if ic2 <= ic
                     continue
                 end
                 nThisModel = size(xtrain[mCounter],1)
                 if ( length(ytrain[mCounter]) != nThisModel || length(α[mCounter]) != nThisModel) error("xtrain, ytrain and α must all have the same length."); end
                 # note that we assign "winning" scores between pair of classes matches only based on who win, not by how much he did
                 # todo check
                 #if sum([ α[mCounter][j] * ((ytrain[mCounter][j] .== c) .*2 .-1) * K(x[i,:],xtrain[mCounter][j,:]) for j in 1:nThisModel]) > eps()
                 score = sum([ α[mCounter][j] * (ytrain[mCounter][j]) * K(x[i,:],xtrain[mCounter][j,:]) for j in 1:nThisModel])
                 #println(score)
                 if score > eps()
                     countByClass[ic] += 1
                     #countByClass[ic] += score
                 else
                     countByClass[ic2] += 1
                     #countByClass[ic2] += -score
                 end
                 mCounter += 1
             end
         end
         #println(countByClass)
         prob = softmax(countByClass)
         y[i] = Dict(zip(classes,prob))
    end
    return y
end

# MLJ interface
include("Perceptron_MLJ.jl")

end


text:
# Type aliasing
const vx_char = UInt8
const vx_uint8 = UInt8
const vx_uint16 = UInt16
const vx_uint32 = UInt32
const vx_uint64 = UInt64
const vx_int8 = Int8
const vx_int16 = Int16
const vx_int32 = Int32
const vx_int64 = Int64
const vx_bitfield = UInt32
const vx_float32 = Cfloat
const vx_float64 = Cdouble
const vx_reference = Ptr{Cvoid}
const vx_enum = Int32
const vx_size = Csize_t
const vx_df_image = UInt32
const vx_map_id = Csize_t
const vx_bool = vx_enum
const vx_status = vx_enum
const vx_action = vx_enum

# Opaque pointer types
const vx_scalar = Ptr{Cvoid}
const vx_image = Ptr{Cvoid}
const vx_kernel = Ptr{Cvoid}
const vx_parameter = Ptr{Cvoid}
const vx_node = Ptr{Cvoid}
const vx_graph = Ptr{Cvoid}
const vx_context = Ptr{Cvoid}
const vx_delay = Ptr{Cvoid}
const vx_lut = Ptr{Cvoid}
const vx_distribution = Ptr{Cvoid}
const vx_matrix = Ptr{Cvoid}
const vx_pyramid = Ptr{Cvoid}
const vx_threshold = Ptr{Cvoid}
const vx_convolution = Ptr{Cvoid}
const vx_remap = Ptr{Cvoid}
const vx_array = Ptr{Cvoid}
const vx_object_array = Ptr{Cvoid}
const vx_tensor = Ptr{Cvoid}
const vx_meta_format = Ptr{Cvoid}
const vx_nodecomplete_f = Ptr{Cvoid}
const vx_publish_kernels_f = Ptr{Cvoid}
const vx_unpublish_kernels_f = Ptr{Cvoid}
const vx_kernel_f = Ptr{Cvoid}
const vx_kernel_initialize_f = Ptr{Cvoid}
const vx_kernel_deinitialize_f = Ptr{Cvoid}
const vx_kernel_validate_f = Ptr{Cvoid}
const vx_kernel_image_valid_rectangle_f = Ptr{Cvoid}
const vx_log_callback_f = Ptr{Cvoid}
const vx_import = Ptr{Cvoid}
const vx_platform = Ptr{Cvoid}
const vx_kernel_input_validate_f = Ptr{Cvoid}
const vx_kernel_output_validate_f = Ptr{Cvoid}
const vx_tiling_kernel_f = Ptr{Cvoid}

# Skipped: Function prototype for VX_CALLBACK

# Structs
struct vx_imagepatch_addressing_t
    dim_x::vx_uint32
    dim_y::vx_uint32
    stride_x::vx_int32
    stride_y::vx_int32
    scale_x::vx_uint32
    scale_y::vx_uint32
    step_x::vx_uint32
    step_y::vx_uint32
    function vx_imagepatch_addressing_t(;dim_x = 0, dim_y = 0, stride_x = 0, stride_y = 0, scale_x = 0, scale_y = 0, step_x = 0, step_y = 0)
        new(dim_x, dim_y, stride_x, stride_y, scale_x, scale_y, step_x, step_y)
    end
end

const IMAGEPATCH_ADDR_INIT = vx_imagepatch_addressing_t()

struct vx_perf_t
    tmp::vx_uint64
    beg::vx_uint64
    _end::vx_uint64
    sum::vx_uint64
    avg::vx_uint64
    min::vx_uint64
    num::vx_uint64
    max::vx_uint64
end

const PERF_INIT = vx_perf_t(0, 0, 0, 0, 0, 0, 0, 0)

struct vx_hough_lines_p_t
    rho::vx_float32
    theta::vx_float32
    threshold::vx_int32
    line_length::vx_int32
    line_gap::vx_int32
    theta_max::vx_float32
    theta_min::vx_float32
end

struct vx_line2d_t
    start_x::vx_float32
    start_y::vx_float32
    end_x::vx_float32
    end_y::vx_float32
end

struct vx_tensor_matrix_multiply_params_t
    transpose_input1::vx_bool
    transpose_input2::vx_bool
    transpose_input3::vx_bool
end

struct vx_kernel_info_t
    enumeration::vx_enum
    name::NTuple{256, vx_char}
end

struct vx_keypoint_t
    x::vx_int32
    y::vx_int32
    strength::vx_float32
    scale::vx_float32
    orientation::vx_float32
    tracking_status::vx_int32
    error::vx_float32
end

struct vx_rectangle_t
    start_x::vx_uint32
    start_y::vx_uint32
    end_x::vx_uint32
    end_y::vx_uint32
end

struct vx_coordinates2d_t
    x::vx_uint32
    y::vx_uint32
end

struct vx_coordinates2df_t
    x::vx_float32
    y::vx_float32
end

struct vx_coordinates3d_t
    x::vx_uint32
    y::vx_uint32
    z::vx_uint32
end

struct vx_pixel_value_t
    RGB::NTuple{3, vx_uint8}
    RGBX::NTuple{4, vx_uint8}
    YUV::NTuple{3, vx_uint8}
    U8::vx_uint8
    U16::vx_uint16
    S16::vx_int16
    U32::vx_uint32
    S32::vx_int32
    reserved::NTuple{16, vx_uint8}
    function vx_pixel_value_t(;RGB=(0,0,0), RGBX=(0,0,0,0), YUV=(0,0,0), U8=0, U16=0, S16=0,
                              U32=0, S32=0, reserved=(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))
        new(RGB, RGBX, YUV, U8, U16, S16, U32, S32, reserved)
    end
end

struct vx_hog_t
    cell_width::vx_int32
    cell_height::vx_int32
    block_width::vx_int32
    block_height::vx_int32
    block_stride::vx_int32
    num_bins::vx_int32
    window_width::vx_int32
    window_height::vx_int32
    window_stride::vx_int32
    threshold::vx_float32
end

struct vx_border_t
    mode::vx_enum
    constant_value::vx_pixel_value_t
end

struct vx_delta_rectangle_t
    delta_start_x::Cint
    delta_start_y::Cint
    delta_end_x::Cint
    delta_end_y::Cint
end

struct vx_nn_convolution_params_t
    padding_x::vx_size
    padding_y::vx_size
    overflow_policy::vx_enum
    rounding_policy::vx_enum
    down_scale_size_rounding::vx_enum
    dilation_x::vx_size
    dilation_y::vx_size
end

struct vx_nn_deconvolution_params_t
    padding_x::vx_size
    padding_y::vx_size
    overflow_policy::vx_enum
    rounding_policy::vx_enum
    a_x::vx_size
    a_y::vx_size
end

struct vx_nn_roi_pool_params_t
    pool_type::vx_enum
end
## TODO: perhaps need to redefine these, in C it looks like array of structs per
#field of the specs for find/replace
struct ANONYMOUS1_find
    kernel_id::vx_enum
    arg_spec::NTuple{32, vx_uint32}
end

struct ANONYMOUS2_replace
    kernel_id::vx_enum
    arg_spec::NTuple{32, vx_uint32}
end

struct AgoNodeMergeRule
    find::ANONYMOUS1_find
    replace::ANONYMOUS2_replace
end

struct vx_tile_block_size_t
    width::vx_int32
    height::vx_int32
end

struct vx_neighborhood_size_t
    left::vx_int32
    right::vx_int32
    top::vx_int32
    bottom::vx_int32
end

struct vx_image_description_t
    width::vx_uint32
    height::vx_uint32
    format::vx_df_image
    planes::vx_uint32
    range::vx_enum
    space::vx_enum
end

struct vx_tile_t
    base::NTuple{4, Ptr{vx_uint8}}
    tile_x::vx_uint32
    tile_y::vx_uint32
    addr::NTuple{4, vx_imagepatch_addressing_t}
    tile_block::vx_tile_block_size_t
    neighborhood::vx_neighborhood_size_t
    image::vx_image_description_t
end

struct AgoImageFormatDescription
    components::vx_size
    planes::vx_size
    pixelSizeInBitsNum::vx_size
    colorSpace::vx_color_space_e
    channelRange::vx_channel_range_e
    pixelSizeInBitsDenom::vx_size
end

struct AgoTargetAffinityInfo
    device_type::vx_uint32
    device_info::vx_uint32
    reserved::NTuple{2, vx_uint32}
end

struct AgoContextTextMacroInfo
    macroName::NTuple{256, vx_char}
    text::Cstring
end

struct AgoGraphImportInfo
    text::Cstring
    num_ref::vx_uint32
    ref::Ptr{vx_reference}
    dumpToConsole::vx_int32
    data_registry_callback_f::Ptr{Cvoid}
    data_registry_callback_obj::Ptr{Cvoid}
end

struct AgoGraphExportInfo
    fileName::NTuple{256, vx_char}
    num_ref::vx_uint32
    ref::Ptr{vx_reference}
    comment::NTuple{64, vx_char}
end