"""
A tableau representation of the non-commutative canonical form of a set of Paulis,
which is used in [`commutify`](@ref).

They are organized in the same form as [`MixedDestabilizer`](@ref)
with a stabilizer, destabilizer, logical X, and logical Z components.
"""
mutable struct SubsystemCodeTableau{T <: Tableau} <: AbstractStabilizer
    tab::T
    index::Int
    r::Int
    m::Int
    k::Int
end

function SubsystemCodeTableau(t::Tableau)
    index = 1
    for i in range(1, stop=length(t), step=2)
        if i + 1 > length(t)
            break
        end
        if comm(t[i], t[i+1]) == 0x01
            index = i+2 # index to split t into non-commuting pairs and commuting operators
        end
    end
    s = Stabilizer(t[index:length(t)])
    ind = 1
    if length(s)>nqubits(s)#if stabilizer is overdetermined, Destabilizer constructor throws error
        m = length(s)
        tab = zero(Tableau, length(t)+length(s), nqubits(t))
        for i in s
            tab[ind] = zero(PauliOperator, nqubits(s))
            ind+=1
        end
    else
        d = Destabilizer(s)
        m = length(d)
        tab = zero(Tableau, length(t)+length(d), nqubits(t))
        for p in destabilizerview(d)
            tab[ind] = p
            ind+=1
        end
    end
    for i in range(1, stop=index-1, step=2)
        tab[ind] = t[i]
        ind+=1
    end
    for p in s
        tab[ind] = p
        ind+=1
    end
    for i in range(2, stop=index, step=2)
        tab[ind] = t[i]
        ind+=1
    end
    return SubsystemCodeTableau(tab, index, length(s), m, Int((index-1)/2))
end


Base.copy(t::SubsystemCodeTableau) = SubsystemCodeTableau(copy(t.tab), t.index, t.r, t.m, t.k)

function Base.show(io::IO, t::SubsystemCodeTableau)
    r = t.r+t.m+2*t.k
    q = nqubits(t)
    if get(io, :compact, false) | haskey(io, :typeinfo)
        print(io, "MixedDestablizer $r×$q")
    elseif get(io, :limit, false)
        h,w = displaysize(io)
        println(io, "𝒟ℯ𝓈𝓉𝒶𝒷" * "━"^max(min(w-9,size(t.tab,2)-4),0))
        _show(io, destabilizerview(t).tab, w, h÷4)
        if r != q
            println(io)
            println(io, "𝒳" * "━"^max(min(w-5,size(t.tab,2)),0))
            _show(io, logicalxview(t).tab, w, h÷4)
        end
        println(io)
        println(io, "𝒮𝓉𝒶𝒷" * "━"^max(min(w-7,size(t.tab,2)-2),0))
        _show(io, stabilizerview(t).tab, w, h÷4)
        if r != q
            println(io)
            println(io, "𝒵" * "━"^max(min(w-5,size(t.tab,2)),0))
            _show(io, logicalzview(t).tab, w, h÷4)
        end
    else
        println(io, "𝒟ℯ𝓈𝓉𝒶𝒷" * "━"^max(size(t.tab,2)-4,0))
        _show(io, destabilizerview(t).tab, missing, missing)
        if r != q
            println(io)
            println(io, "𝒳ₗ" * "━"^max(size(t.tab,2),0))
            _show(io, logicalxview(t).tab, missing, missing)
        end
        println(io)
        println(io, "𝒮𝓉𝒶𝒷" * "━"^max(size(t.tab,2)-2,0))
        _show(io, stabilizerview(t).tab, missing, missing)
        if r != q
            println(io)
            println(io, "𝒵ₗ" * "━"^max(size(t.tab,2)))
            _show(io, logicalzview(t).tab, missing, missing)
        end
    end
end

@inline stabilizerview(s::SubsystemCodeTableau) = Stabilizer(@view tab(s)[s.m+s.k+1:s.m+s.k+s.r])
@inline destabilizerview(s::SubsystemCodeTableau) = Stabilizer(@view tab(s)[1:s.r])
@inline logicalxview(s::SubsystemCodeTableau) = Stabilizer(tab(s)[s.m+1:s.m+s.k])
@inline logicalzview(s::SubsystemCodeTableau) = Stabilizer(tab(s)[length(s.tab)-s.k+1:length(s.tab)])

"""
Return the full stabilizer group represented by the input generating set (a [`Stabilizer`](@ref)).

The returned object is exponentially long.

```jldoctest
julia> groupify(S"XZ ZX")
+ __
+ XZ
+ ZX
+ YY
```
"""
function groupify(s::Stabilizer; phases=false)
    # Create a `Tableau` of 2ⁿ n-qubit identity Pauli operators(where n is the size of
    # `Stabilizer` s), then multiply each one by a different subset of the elements in s to
    # create all 2ⁿ unique elements in the group generated by s, then return the `Tableau`.
    if phases == true
        throw(ArgumentError("in groupify phases=true functionality not yet implemented"))
    end
    gen_set = minimal_generating_set(s)
    n = length(gen_set)::Int
    group = zero(Tableau, 2^n, nqubits(gen_set))
    for i in 0:2^n-1
        for (digit_order, j) in enumerate(digits(i, base=2, pad=n))
            if j == 1
                mul_left!(group, i+1, tab(gen_set), digit_order, phases=Val(false))
            end
        end
    end
    return group
end


"""
For a not-necessarily-minimal generating set,
return the minimal generating set.

The input has to have only real phases.

```jldoctest
julia> minimal_generating_set(S"__ XZ ZX YY")
+ XZ
+ ZX
```
"""
function minimal_generating_set(s::Stabilizer)
    # Canonicalize `Stabilizer` s, then return a `Stabilizer` with all non-identity Pauli operators
    # in the result. If s consists of only identity operators, return the negative
    # identity operator if one is contained in s, and the positive identity operator otherwise.
    s, _, r = canonicalize!(copy(s), ranks=true)
    if r == 0
        gs = zero(Stabilizer, 1, nqubits(s))
        if 0x02 in phases(s)
            phases(gs)[1] = 0x2
        end
        return gs
    else
        return s[1:r, :]
    end
end

"""
For a not-necessarily commutative set of Paulis, return a generating set of the form
⟨A₁, A₂, ... Aₖ, Aₖ₊₁, ... Aₘ, B₁, B₂, ... Bₖ⟩ where pairs Aₖ, Bₖ anticommute and all other pairings commute. Based on [RevModPhys.87.307](@cite)

Returns the generating set as a data structure of type [`SubsystemCodeTableau`](@ref). The
`logicalxview` function returns the ⟨A₁, A₂,... Aₖ⟩, and the `logicalzview`
function returns ⟨B₁, B₂, ... Bₖ⟩. `stabilizerview` returns ⟨Aₖ₊₁, ... Aₘ⟩
as a Stabilizer, and `destabilizerview` returns the Destabilizer of that Stabilizer.

Phases are zeroed-out in this canonicalization.

```jldoctest
julia> canonicalize_noncomm(T"XX XZ XY")
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z_
𝒳━━
+ XX
𝒮𝓉𝒶𝒷
+ X_
𝒵━━
+ XZ
```
"""
function canonicalize_noncomm(t::Tableau)
    for i in eachindex(t) phases(t)[i] = 0x00 end
    loc = zero(Tableau, 2*length(t), nqubits(t))
    x_index = 0
    z_index = length(loc)
    for i in eachindex(t)
        for j in eachindex(t)
            if comm(t[i], t[j]) == 0x01
                for k in eachindex(t)
                    if k !=i && k != j
                        if comm(t[k], t[i]) == 0x01
                            mul_left!(t, k, j, phases=Val(false))
                        end
                        if comm(t[k], t[j]) == 0x01
                            mul_left!(t, k, i, phases=Val(false))
                        end
                    end
                end
                if !(t[i] in loc)
                    x_index+=1
                    loc[x_index]= t[i]
                end
                if !(t[j] in loc)
                    loc[z_index]= t[j]
                    z_index-=1
                end
            end
        end
    end
    ind = 2*x_index+1 #format tableau for SubsystemCodeTableau
    k = x_index
    m = length(t)-2*k
    for i in range(k, stop =1, step=-1)
        loc[i+m]= loc[i]
    end
    i = m+1
    j = m+k
    while i < j #ensure anticommutative pairs are lined up correctly
        loc[i], loc[j], = loc[j], loc[i]
        i+=1
        j-=1
    end
    s_index = m+k+1
    for i in eachindex(t)
        if !(t[i] in loc)
            loc[s_index]= t[i]
            s_index+=1
        end
    end
    r = s_index-m-k-1
    if r <= nqubits(t)
        d = destabilizerview(Destabilizer(Stabilizer(loc[m+k+1:s_index-1])))
        for i in eachindex(d)
            loc[i] =d[i]
        end
    else for i in 1:m zero!(loc, i) end end
    return SubsystemCodeTableau(loc, ind, r, m, k)
end

canonicalize_noncomm(ps::Base.AbstractVecOrTuple{PauliOperator}) = canonicalize_noncomm(Tableau(ps))

"""
For a not-necessarily commutative set of Paulis S,
computed S', the [non-commutative canonical form](@ref canonicalize_noncomm) of of S.
For each pair Aₖ, Bₖ of anticommutative Paulis in S', add a qubit to each Pauli in the set:
X to Aₖ, Z to Bₖ, and I to each other operator to produce S'', a fully commutative set. Return
S'' as well as a list of the indices of the added qubits.

The returned object is a Stabilizer that is also useful for the [`matroid_parent`](@ref) function.

```jldoctest
julia> commutify(T"XX XZ XY")[1]
+ XXX
+ X__
+ XZZ

julia> commutify(T"XX XZ XY")[2]
3:3
```
"""
function commutify(t)
    loc = canonicalize_noncomm(t)
    commutative = zero(Stabilizer, 2*loc.k+loc.r, nqubits(loc)+loc.k)
    puttableau!(commutative, logicalxview(loc), 0,0)
    puttableau!(commutative, stabilizerview(loc), loc.k,0)
    puttableau!(commutative, logicalzview(loc), loc.k+loc.r,0)
    x= T"X"
    z=T"Z"
    for i in 0:(loc.k-1)
        puttableau!(commutative, x, i, nqubits(loc)+i)
        puttableau!(commutative, z, i+loc.r+loc.k, nqubits(loc)+i)
    end
    to_delete = nqubits(loc)+1:nqubits(loc)+loc.k
    return commutative, to_delete
end

"""
For a given set S of Paulis that does not necessarily represent a state,
return a set of Paulis S' that represents a state.
S' is a superset of [commutified](@ref commutify) S.
Additionally returns two arrays representing deletions needed to produce S.
Based on [goodenough2024bipartiteentanglementnoisystabilizer](@cite)

By deleting the qubits in the first output array from S', taking the [`normalizer`](@ref) of S', then
deleting the qubits in the second returned array from the [`normalizer`](@ref) of S', S is reproduced.

```jldoctest; filter = r"( [(]empty range[)])?\$"
julia> matroid_parent(T"XX")[1]
+ X_X
+ XX_
+ ZZZ

julia> matroid_parent(T"XX")[2]
3:3

julia> matroid_parent(T"XX")[3]
3:2
```
"""
function matroid_parent(t::Tableau)
    com, d1= commutify(t)
    norm = normalizer(com.tab)
    state, d2 = commutify(norm)
    return state, d2, d1
end

"""
Return the full Pauli group of a given length. Phases are ignored by default,
but can be included by setting `phases=true`.

```jldoctest
julia> pauligroup(1)
+ _
+ X
+ Z
+ Y

julia> pauligroup(1, phases=true)
+ _
+ X
+ Z
+ Y
- _
- X
- Z
- Y
+i_
+iX
+iZ
+iY
-i_
-iX
-iZ
-iY
```
"""
function pauligroup(n::Int; phases=false)
    if phases
        s = zero(Tableau, 4^(n + 1), n)
        paulis = ((false, false), (true, false), (false, true), (true, true))
        for (i, P) in enumerate(Iterators.product(Iterators.repeated(paulis, n)...))
            for (j, p) in enumerate(P)
                s[i, j] = p
            end
        end
        for i in 1:4^n
            s[i+4^n] = -1 * s[i]
        end
        for i in 4^n+1:2*4^n
            s[i+4^n] = -1im * s[i]
        end
        for i in 2*4^n+1:3*4^n
            s[i+4^n] = -1 * s[i]
        end
        return s
    else
        s = zero(Tableau, 4^n, n)
        paulis = ((false, false), (true, false), (false, true), (true, true))
        for (i, P) in enumerate(Iterators.product(Iterators.repeated(paulis, n)...))
            for (j, p) in enumerate(P)
                s[i, j] = p
            end
        end
        return s
    end
end

"""
Return all Pauli operators with the same number of qubits as the given `Tableau` `t`
that commute with all operators in `t`.

```jldoctest
julia> normalizer(T"X")
+ _
+ X
```
"""
function normalizer(t::Tableau; phases=false)
    # For each `PauliOperator` p in the with same number of qubits as the `Stabilizer` s, iterate through s and check each
    # operator's commutivity with p. If they all commute, add p a vector of `PauliOperators`. Return the vector
    # converted to `Tableau`.
    n = nqubits(t)
    ptype =typeof(P"I")
    norm = ptype[]

    p = zero(PauliOperator, n)
    paulis = ((false, false), (true, false), (false, true), (true, true))
    for i in Iterators.product(Iterators.repeated(paulis, n)...)
        zero!(p)
        for (j, k) in enumerate(i)
            p[j] = k
        end
        commutes = true
        for q in t
            if comm(p, q) == 0x01
                commutes = false
            end
        end
        if commutes
            push!(norm, copy(p))
        end
        if phases
            for phase in [-1, 1im, -1im]
                push!(norm, phase *p)
            end
        end
    end
    return Tableau(norm)
end

"""
For a given set of Paulis (in the form of a `Tableau`), return the subset of Paulis that commute with all Paulis in set.

```jldoctest
julia> centralizer(T"XX ZZ _Z")
+ ZZ
```
"""
function centralizer(t::Tableau)
    center = typeof(t[1])[]
    for P in t
        commutes = 0
        for Q in t
            if comm(P, Q) == 0x01
                commutes = 1
                break
            end
        end
        if commutes == 0
            push!(center, P)
        end
    end
    if length(center) == 0
        return Tableau(zeros(Bool, 1,1))
    end
    c = Tableau(center)
    return c
end

"""
Return the subset of Paulis in a Stabilizer that have identity operators on all qubits corresponding to
the given subset, without the entries corresponding to subset. Based on [goodenough2024bipartiteentanglementnoisystabilizer](@cite)

```jldoctest
julia> contractor(S"_X X_", [1])
+ X
```
"""
function contractor(s::Stabilizer, subset)
    result = typeof(s[1])[]
    for p in s
        contractable = true
        for i in subset
            if p[i] != (false, false)
                contractable = false
                break
            end
        end
        if contractable push!(result, p[setdiff(1:length(p), subset)]) end
    end
    if length(result) > 0
        return Tableau(result)
    else
        return Tableau(zeros(Bool, 1,1))
    end
end
