"""
Clifford Operator specified by the mapping of the basis generators.

```jldoctest
julia> tCNOT
X₁ ⟼ + XX
X₂ ⟼ + _X
Z₁ ⟼ + Z_
Z₂ ⟼ + ZZ

julia> phase_gate = C"Y
                      Z"
X₁ ⟼ + Y
Z₁ ⟼ + Z

julia> stab = S"XI
                IZ";


julia> entangled = tCNOT*stab
+ XX
+ ZZ

julia> CliffordOperator(T"YY")
ERROR: DimensionMismatch: Input tableau should be of size 2n×n (top half is the X mappings and the bottom half are the Z mappings).
[...]
```

[`Destabilizer`](@ref) can also be converted.
```jldoctest
julia> d = Destabilizer(S"Y")
𝒟ℯ𝓈𝓉𝒶𝒷
+ X
𝒮𝓉𝒶𝒷
+ Y

julia> CliffordOperator(d)
X₁ ⟼ + X
Z₁ ⟼ + Y
```
"""
struct CliffordOperator{T<:Tableau,P<:PauliOperator} <: AbstractCliffordOperator
    tab::T
    buffer::P
    function CliffordOperator(tab::Tableau)
        if size(tab,1)==2*size(tab,2)
            p = zero!(tab[1])
            new{typeof(tab),typeof(p)}(tab,p)
        #elseif size(stab,1)==size(stab,2) # TODO be able to work with squara tableaux (by reversing all row operations)
        #    destab = tab(Destabilizer(stab))
        #    new{typeof(destab.phases),typeof(destab.xzs)}(destab) # TODO be smarter about type signatures here... there should be a better way
        else
            throw(DimensionMismatch("Input tableau should be of size 2n×n (top half is the X mappings and the bottom half are the Z mappings)."))
        end
    end
end

macro C_str(a)
    tab = _T_str(a)
    quote CliffordOperator($tab) end
end

CliffordOperator(op::CliffordOperator) = op
CliffordOperator(paulis::AbstractVector{<:PauliOperator}) = CliffordOperator(Tableau(paulis))
CliffordOperator(destab::Destabilizer) = CliffordOperator(tab(destab))
function CliffordOperator(pauli::PauliOperator)
    res = one(CliffordOperator, nqubits(pauli))
    comm!(phases(res), pauli, tab(res))
    tab(res).phases .*= 0x02
    return res
end

Base.:(==)(l::CliffordOperator, r::CliffordOperator) = tab(l) == tab(r)
Base.hash(c::T, h::UInt) where {T<:CliffordOperator} = hash(T, hash(tab(c), h))

Base.getindex(c::CliffordOperator, args...) = getindex(tab(c), args...)
Base.setindex!(c::CliffordOperator, args...) = setindex!(tab(c), args...)

tab(c::CliffordOperator) = c.tab

Base.size(c::CliffordOperator,args...) = size(tab(c),args...)

function row_limit(str, limit=50)
    n = length(str)
    if (n <= limit || limit==-1)
        return str
    end
    padding = Int64(floor(limit/2))
    return SubString(str, 1, padding) * " ... " * SubString(str, n - padding, n - 1)
end

digits_subchars = collect("₀₁₂₃₄₅₆₇₈₉")
digits_substr(n::Int,nwidth::Int) = join(([digits_subchars[d+1] for d in reverse(digits(n, pad=nwidth))]))

function Base.copy(c::CliffordOperator)
    CliffordOperator(copy(tab(c)))
end

@inline nqubits(c::CliffordOperator) = nqubits(tab(c))
@inline phases(c::CliffordOperator) = phases(tab(c))

Base.zero(c::CliffordOperator) = CliffordOperator(zero(tab(c)))
Base.zero(::Type{<:CliffordOperator}, n) = CliffordOperator(zero(Tableau, 2n, n))

function Base.:(*)(l::AbstractCliffordOperator, r::CliffordOperator)
    tab = copy(QuantumClifford.tab(r))
    apply!(Stabilizer(tab),l) # TODO maybe not the most elegant way to perform apply!(::Tableau, gate)
    CliffordOperator(tab)
end

function apply!(r::CliffordOperator, l::AbstractCliffordOperator; phases=false)
    @valbooldispatch _apply!(Stabilizer(tab(r)),l;phases=Val(phases)) phases # TODO maybe not the most elegant way to perform apply!(::Tableau, gate)
    r
end

function apply_inv!(r::CliffordOperator, l::AbstractCliffordOperator; phases=false)
    @valbooldispatch _apply_inv!(Stabilizer(tab(r)),l,phases=Val(phases)) phases
    r
end

"""Nonvectorized version of `apply!` used for unit tests."""
function _apply_nonthread!(stab::AbstractStabilizer, c::CliffordOperator; phases::Bool=true)
    nqubits(stab)==nqubits(c) || throw(DimensionMismatch("The tableau and the Clifford operator need to act on the same number of qubits. Consider calling `apply!(state, indices, operation)` to apply the operation to selected subsystems."))
    s_tab = tab(stab)
    c_tab = tab(c)
    new_stabrow = c.buffer
    for row_stab in eachindex(s_tab)
        zero!(new_stabrow)
        apply_row_kernel!(new_stabrow, row_stab, s_tab, c_tab, phases=Val(phases))
    end
    stab
end

# TODO no need to track phases outside of stabview
function _apply!(stab::AbstractStabilizer, c::CliffordOperator; phases::Val{B}=Val(true)) where B
    nqubits(stab)==nqubits(c) || throw(DimensionMismatch("The tableau and the Clifford operator need to act on the same number of qubits. Consider calling `apply!(state, indices, operation)` to apply the operation to selected subsystems."))
    s_tab = tab(stab)
    c_tab = tab(c)
    threadlocal = c.buffer
    @inbounds for row_stab in eachindex(s_tab)
        zero!(threadlocal) # a new stabrow for temporary storage
        apply_row_kernel!(threadlocal, row_stab, s_tab, c_tab, phases=phases)
    end
    stab
end

function _apply_inv!(stab::AbstractStabilizer, c::CliffordOperator; phases::Val{B}=Val(true)) where B
    _apply!(stab, inv(c); phases=phases)
end

# TODO Added a lot of type assertions to help Julia infer types, but they are much too strict for cases where bitpacking varies (check tests)
@inline function apply_row_kernel!(new_stabrow, row, s_tab, c_tab; phases::Val{B}=Val(true)) where B
    B && (new_stabrow.phase[] = s_tab.phases[row])
    n = nqubits(c_tab)
    for qubit in 1:n
        x,z = s_tab[row,qubit]
        if B&&x&&z
            new_stabrow.phase[] -= 0x1
        end
        if x
            mul_left!(new_stabrow, c_tab, qubit, phases=phases)
        end
        if z
            mul_left!(new_stabrow, c_tab, qubit+n, phases=phases)
        end
    end
    s_tab[row] = new_stabrow
    new_stabrow
end

"""Nonvectorized version of `apply!` used for unit tests."""
function _apply_nonthread!(stab::AbstractStabilizer, c::CliffordOperator, indices_of_application::Base.AbstractVecOrTuple{Int}; phases::Bool=true)
    s_tab = tab(stab)
    c_tab = tab(c)
    new_stabrow = c.buffer
    for row in eachindex(s_tab)
        zero!(new_stabrow)
        apply_row_kernel!(new_stabrow, row, s_tab, c_tab, indices_of_application; phases=Val(phases))
    end
    stab
end

#TODO a lot of code repetition with apply!(stab::AbstractStabilizer, c::CliffordOperator; phases::Bool=true) and apply_row_kernel!
function _apply!(stab::AbstractStabilizer, c::CliffordOperator, indices_of_application::Base.AbstractVecOrTuple{Int}; phases::Val{B}=Val(true)) where B
    #max(indices_of_application)<=nqubits(s) || throw(DimensionMismatch("")) # Too expensive to check every time
    s_tab = tab(stab)
    c_tab = tab(c)
    threadlocal= c.buffer
    @inbounds for row_stab in eachindex(s_tab)
        zero!(threadlocal) # a new stabrow for temporary storage
        apply_row_kernel!(threadlocal, row_stab, s_tab, c_tab, indices_of_application, phases=phases)
    end
    stab
end

function _apply_inv!(stab::AbstractStabilizer, c::CliffordOperator, indices_of_application::Base.AbstractVecOrTuple{Int}; phases::Val{B}=Val(true)) where B
    _apply!(stab, inv(c), indices_of_application; phases=phases)
end

@inline function apply_row_kernel!(new_stabrow, row, s_tab, c_tab, indices_of_application; phases::Val{B}=Val(true)) where B
    B && (new_stabrow.phase[] = s_tab.phases[row])
    n = nqubits(c_tab)
    for (qubit_i, qubit) in enumerate(indices_of_application)
        x,z = s_tab[row,qubit]
        if B&&x&&z
            new_stabrow.phase[] -= 0x1
        end
        if x
            mul_left!(new_stabrow, c_tab, qubit_i, phases=phases)
        end
        if z
            mul_left!(new_stabrow, c_tab, qubit_i+n, phases=phases)
        end
    end
    for (qubit_i, qubit) in enumerate(indices_of_application)
        s_tab[row,qubit] = new_stabrow[qubit_i]
    end
    B && (s_tab.phases[row] = new_stabrow.phase[])
    new_stabrow
end

const tCNOT = C"XX
                IX
                ZI
                ZZ"

const tCPHASE = C"XZ
                  ZX
                  ZI
                  IZ"

const tSWAP = C"IX
                XI
                IZ
                ZI"

const tHadamard = C"Z
                    X"

const tPhase = C"Y
                 Z"

const tId1 = C"X
               Z"
