implement sm⊗sm, tHadamard*sm, pcT⊗P"X" , pcT*tId1 , P"-Y"*sm and copy(sm) from #260

[Fe-r-oz]

This PR aims to address some TODOs to become familiar with @Krastanov's very cool work on non-clifford functionality with the goal to resolve #309 by completing the non-clifford functionalities.

[codecov]

Codecov Report

Attention: Patch coverage is 81.25000% with 12 lines in your changes missing coverage. Please review.

Project coverage is 82.89%. Comparing base (59e399d) to head (b51df40).

Files	Patch %	Lines
src/nonclifford.jl	76.59%	11 Missing ⚠️
src/QuantumClifford.jl	90.90%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #316      +/-   ##
==========================================
- Coverage   82.97%   82.89%   -0.08%     
==========================================
  Files          61       61              
  Lines        4023     4069      +46     
==========================================
+ Hits         3338     3373      +35     
- Misses        685      696      +11     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[Fe-r-oz]

I hope the TODOs are satisfactory. After the TODOs, I will again mark the PR as draft.

Attaching the benchmark results for TODO1

julia> using BenchmarkTools

julia> n = 100000; a = rand(Int32, n);b = rand(Int32, n); c = rand(Int32, n);

julia> function _allthreesumtozero(a, b, c)
         @inbounds @simd for i in 1:length(a)
           iseven(a[i] + b[i] + c[i]) || return false
         end
         true
       end
_allthreesumtozero (generic function with 1 method)

julia> function _allthreesumtozero_optimized(a, b, c)
           @inbounds @simd for i in 1:length(a)
               ((a[i] + b[i] + c[i]) & 1) == 0 || return false
           end
           true
       end
_allthreesumtozero_optimized (generic function with 1 method)

julia> println(@benchmark _allthreesumtozero($a, $b, $c))
Trial(10.601 ns)

julia> println(@benchmark _allthreesumtozero_optimized($a, $b, $c))
Trial(3.505 ns)

[Fe-r-oz]

The failure of the downgrade and build test is due to stuff not related to this PR, It has been pointed out in the issue 299 and 245. As you said back then (i245) that it's a statistical test which fails sometime in a thousand runs but recently, it's been failing quite often.

 Test Summary:                                  | Pass  Fail  Total  Time
noisycircuits                                  |   90     1     91  7.2s
  Noisy Gates                                  |    2            2  1.5s
  Monte Carlo Purification examples            |    9            9  0.5s
  Perturbative expansion Purification examples |   11     1     12  4.3s
    Comparison to MC                           |    8     1      9  2.1s
    Symbolic                                   |    3            3  2.2s
  Measurements                                 |   53           53  0.5s
  Classical Bits                               |   15           15  0.2s

[Fe-r-oz]

What's changed in apply! ? In summary, using filter instead of deleting keys while iterating along with parameterization of prunethreshold. goal: To make it a cleaner and safer code.

I sought some advice on your question/concern: "Is this actually necessary? Is it actually dangerous to modify the dictionary keys while iterating over it?" After discussions, it turned out that this is dangerous and filter! is safe approach according to folks."

Valentin provided good reasons (forwarded them) of why it's dangerous and suggested that filter! is safe. 'why not use filter!(.., result_dict), which handles the edge case of rehashing correctly for you' and 'deleting like you're doing is definitely dangerous'. Hence, So, I have used the filter! and also tried to parameterize the prune threshold. The default threshold is 1e-14 and user can enter their desired threshold. If they don't enter anything, then the default prune threshold is used.

function apply!(state::GeneralizedStabilizer, gate::AbstractPauliChannel; prune_threshold::Union{Nothing, Float64}=nothing)
    if prune_threshold === nothing
        prune_threshold = 1e-14  # Default value
    end
    dict = state.destabweights
    stab = state.stab
    dtype = _dictvaltype(dict)
    tzero = zero(dtype)
    tone = one(dtype)
    newdict = typeof(dict)(tzero) 
    for ((dᵢ,dⱼ), χ) in dict # the state
        for ((Pₗ,Pᵣ), w) in zip(gate.paulis, gate.weights) # the channel
            phaseₗ, dₗ, dₗˢᵗᵃᵇ = rowdecompose(Pₗ, stab)
            phaseᵣ, dᵣ, dᵣˢᵗᵃᵇ = rowdecompose(Pᵣ, stab)
            c = (dot(dₗˢᵗᵃᵇ, dᵢ) + dot(dᵣˢᵗᵃᵇ, dⱼ)) * 2
            dᵢ′ = dₗ .⊻ dᵢ
            dⱼ′ = dᵣ .⊻ dⱼ
            χ′ = χ * w * (-tone)^c * (im)^(-phaseₗ + phaseᵣ + 4)
            if abs(χ′) >= prune_threshold
                if haskey(newdict, (dᵢ′, dⱼ′))
                    newdict[(dᵢ′, dⱼ′)] += χ′
                else
                    newdict[(dᵢ′, dⱼ′)] = χ′
                end
            end
    end
    filter!(x -> abs(x[2]) >= prune_threshold, newdict)
    state.destabweights = newdict
    state
end

I have made refined further by using in-built get! instead of the if else condition.

Julia documentation for get!: https://docs.julialang.org/en/v1/base/collections/#Base.get!

Using newdict[(dᵢ′, dⱼ′)] = get!(newdict, (dᵢ′, dⱼ′), 0)+χ′ is a more idiomatic and concise compared to the if haskey(newdict, (dᵢ′, dⱼ′)) else . It combines the existence check and value update into a single line, improving readability and using built-in functions for more cleaner code.

function apply!(state::GeneralizedStabilizer, gate::AbstractPauliChannel; prune_threshold::Union{Nothing, Float64}=nothing)
    if prune_threshold === nothing
        prune_threshold = 1e-14  # Default value
    end
    dict = state.destabweights
    stab = state.stab
    dtype = _dictvaltype(dict)
    tzero = zero(dtype)
    tone = one(dtype)
    newdict = typeof(dict)(tzero)
    for ((dᵢ,dⱼ), χ) in dict # the state
        for ((Pₗ,Pᵣ), w) in zip(gate.paulis, gate.weights) # the channel
            phaseₗ, dₗ, dₗˢᵗᵃᵇ = rowdecompose(Pₗ,stab)
            phaseᵣ, dᵣ, dᵣˢᵗᵃᵇ = rowdecompose(Pᵣ,stab)
            c = (dot(dₗˢᵗᵃᵇ,dᵢ) + dot(dᵣˢᵗᵃᵇ,dⱼ))*2
            dᵢ′ = dₗ .⊻ dᵢ
            dⱼ′ = dᵣ .⊻ dⱼ
            χ′ = χ * w * (-tone)^c * (im)^(-phaseₗ+phaseᵣ+4)
           if abs(χ′) >= prune_threshold
               newdict[(dᵢ′,dⱼ′)] = get!(newdict,(dᵢ′,dⱼ′),0)+χ′
           end
    end
    filter!(x -> abs(x[2]) >= prune_threshold, newdict)
    state.destabweights = newdict
    state
end

I hope this is satisfactory.

[Fe-r-oz]

Note: Please note that the latest commit is 67165b9 but it appear before the resolve-merge-conflict(8ace9ee) not the bb0baf0 which appears after the resolve-merge-conflict.

[Fe-r-oz]

In order to become more familiar with non-clifford requirements, I went through some TODOs to understand your cool work so that I can ask better questions. The first goal was first to wrap my head around your implementation by doing TODOs. I hope that the latest TODOs improvements are satisfactory.

Besides, I have also read Ted's paper as well (till the 3 algorithms section) so that I can follow the concepts used here. Before answering questions, it is necessary for me to describe all the context to ask reasonable, efficient questions.

Please skip to questions directly if you like. Please take your time answering them at your convenience.

Given that:

1) Full stabilizer `T = (S, D)` where `S` is full Stabilizer and `D` is full Destabilizer.
2) `T = (S, D)` and Pauli operator `P ∈ Gₙ`, `P = ads` in O(n²) where `α ∈ {±1, ±i}`, `d ∈ D`, and `s ∈ S`.
3) `B(T ) = B(S, D)` stabilizers a unique pure stabilizer state, `|ψₛ⟩`, or `ρₛ = |ψₛ⟩⟨ψₛ| `
4) `dₖρₛdₖ = (½)ⁿ( I - sₖ) ∏ (I + sᵢ) `where `∏` is from `ᵢ ≠ ₖ`
5) `B(T) = B(S, D) = {d|ψₛ⟩  d ∈ D}` ; stabilizer basis corresponding to T.
6) Generalized Stabilizer for τ(arbitrary state) is a two-tuple` (χ, B(S, D))`
where `χ` is a density matrix and ` B(S, D) = B(T )` is the basis in which `χ` is expressed.
7) arbitrary state: `τ = ∑∑ χᵢⱼ dᵢ |ψₛ⟩⟨ψₛ| dⱼ = ∑∑ χᵢⱼ dᵢ ρₛ dⱼ`
8 Inverse sparsity: Λ(τ) = min B(S,D) Λ(χ) = min S Λ(χ) where Λ(χ) is # of non zero elements in density matrix `χ`
9) Bounds: rank(τ) ≤ Λ(τ) ≤ 4ⁿ
10) Pauli channel: `E(τ ) = ∑ϕᵢⱼ Bᵢ ρ Bⱼ†`
where `ϕᵢ` is complex number and `Bᵢ `is Pauli operators `Bᵢ ∈ Gₙ.`
11) Λ(E) ≡ Λ(ϕ) is Pauli Channel complexity.
12) Anything a S representation can do is also possible with a (χ, B(S, D)).
13) This follows from 10) that χ is all zeros except for a single 1 along on the diagonal with appropriate stabilizer basis B(S, D)

Algorithms:

1) Clifford Gates: Suppose Clifford circuit C.
CτC† = ∑ χᵢⱼ CdᵢC† CρₛC† CdⱼC† = ∑ χᵢⱼ d'ᵢ ρₛ' dⱼ'
------------------------------------------------------------
input: τ and Clifford gate, G
output: S', D'
Procedure:
Conjugate each generator of  S: S' = ⟨Gs₁G†, Gs₂G†, ... , GsₙG†⟩
Conjugate each generator of D: D' = ⟨Gd₁G†, Gd₂G†, ... , GdₙG†⟩
------------------------------------------------------------

2) Pauli measurements
----------------------------------------------------------------------------------------------------------------
input: τ = (χ, B(S, D)) and α, b, c from rowdecompose where M = a*d_b*s_c , M ∈ Gₙ
output:  τ' = (χ, B(S', D'))
procedure:
Initialize sparse matrix χ'
if b = 0 then
    for all χᵢⱼ ≠ 0 do
        if α(−1)^(i·c) = 1 and α(−1)^(j·c) = 1 then
            χᵢⱼ' ← χᵢⱼ'  + χᵢⱼ
else
k ← vector of zeros with a single 1 in the location of b’s first 1
    for all χᵢⱼ ≠ 0 do
    x ← i
    y ← j
    q ← 1
    if i · k = 1 then
        q ← α(−1)^(i·c) q
        x ← i + b
    if j · k = 1 then  
        q ← α∗(−1)^(j·c) q
        y ← j + b
     χ_xy' ← χ_xy'  + ½qχᵢⱼ
Return χ', S', D'
Note: (S', D') can be found from (S, D) by traditional stabilizer update via Algorithm 1
----------------------------------------------------------------------------------------------------------------

3) Pauli Channels:
--------------------------------------------------------------------------------
input: τ = (χ, B(S, D)), {αₕ, bₕ, cₕ}ₕ for each Bₕ
output:  χ' s.t. τ = (χ', B(S, D))
procedure:
Initialize sparse matrix χ
    for all ϕₖₗ ≠ 0 do
        for all χᵢⱼ ≠ 0 do
            x ← i + bₖ
            y ← j + bₗ
            χ_xy' ← χ_xy'  + [αₕαₗ*(−1)^(i·cₖ+j·cₗ)]ϕₖₗχᵢⱼ
Return χ'
--------------------------------------------------------------------------------

Going on a tangent:

I was wondering why the destabweights are updated following `apply!(GeneralizedStabilizer(S"-X"), pcT)`.
The answer seems to be the following (taken from the paper):

1) An algorithm that uses just a density matrix representation could also use sparsity to improve its runtime for simulation of such non-stabilizer states through stabilizer circuits.
2) However, simulating g Clifford gates will generally result in an exponential (i.e. 2 ^O(g)) decrease in the sparsity of the density matrix representation, such that, when it comes time to measure the state at the conclusion of the simulation, the calculation will be intractable.
3) Our method avoids this by updating the stabilizer basis instead of the density matrix directly.
4) The price we pay is having to interpret the concluding measurements in terms of our modified basis, but this is polynomial in n by the decomposition lemma.

I see the rowdecompose utility rowdecompose does the following: Given a full tableau T = (S, D) and P ∈ G_n, we can find α ∈ {±1, ±i}, d ∈ D, and s ∈ S such that P = αds in O(n^2) time. Say we want to measure the (hermitian) Pauli operator M ∈ G_n. First, we decompose M in terms of the stabilizer and destabilizer, M = αdb where α = ±1, ±i.

I am also going through project_reset_trace.jl to understand the project! was done for Stabilizer, Destabilizer, MixedStabilizer, MixedDestabilizer objects.

Questions:

In my rough estimate, the algorithm 3 in Ted's paper seems similar for apply!(s::genstab, g::PauliChannel) function and for project! is algorithm 2 (updates the generalized stabilizer τ = (χ, B(S, D)) by a measurement of M = αdb).

1. Why is MixedDestabilizer used instead of Destabilizer as it relaxes the constraint that the 'full stabilizer S stabilizes unique pure state ρₛ = |ψₛ⟩⟨ψₛ| ? The user might, unknowingly, enter a mixed state which the MixedDestabilizer permits. This is also related to the #TODO where you raise this exact question. If we are going with MixedDestabilizer, then we should document this stuff and warn the user to not enter mixed states. Filler details: Stabilizer basis corresponding to the full tableau (T) is given as B(T) = B(S, D) = {d|ψₛ⟩ d ∈ D}. Generalized stabilizer uses B(S, D) to represent arbitrary state τ

2. How much of algorithm 2 is expect(p::PauliOperator, s::GeneralizedStabilizer) handling or is it something preliminary before starting the algorithm 2 for project!(sm::GeneralizedStabilizer, p::PauliOperator)? I see that algorithm 2 starts by considering the two cases for b where b is retrieved from rowdecompose such that M = a*d_b*s_c. The two cases are either (b = (0, 0, ... , 0) or b is nonzero.

In algorithm 2, in the else case (second nested else), there is this α∗ which corresponds to conjugate of 'phase', but phase from rowdecompose returns 0 to 3 depending on 1, i, -1,-i and this continues. E.g. phase =5, isiHow to properly take conjugate of phase?

I think we can use like this: im^(phase) , and then conjugate will be conj(im^(phase)), Correct, Right?

if phase = 0 , then im^(0) = 1 
if phase = 1 , then im^(0) = i
if phase = 2, then im^(0) = -1
if phase = 3 , then im^(0) = -i  

3. What the role of _proj₋ and _proj₊ for project and how it fits in the context of algorithm 2. I might be wrong because I briefly looked at project_reset_trace.jl, but I didn't see this clear differentiation of helper functions there?

4. I think that expect to be a portion of algorithm 2 and is implementing equations (14), (15), (16) from Ted's paper? Does the last line return (-1)^(phase÷2) * e takes care of the fact that you consider both ¼(I + M)τ(I + M) and ¼(I - M)τ(I - M) where M is a Pauli operator? But I think currently, expect is not handling (18-26) equations for the case when b is not 0 where b is retrieved from rowdecompose such that M = a*d_b*s_c. So, for algorithm 2, we are currently dealing with the if case only? P.S: (Q2 and Q4 are connected, I separated them so it's easier to read.)

5. There are some test failures of expect upon which you commented that "I have messed up a convention somewhere and we are getting expectation values that are occasionally imaginary." Also, there is a #TODO optimize, what optimization is needed? How can this be fixed?

6. How is _allthreesumtozero presented in Ted's paper (what does it relate to there), I am not sure where this routine was used in his algorithms. The only thing that I found closely similar about this utility was from algorithm 2 preliminary that: In the case where (b = (0, 0, ... , 0) , either M or −M is a stabilizer group element. In what follows, all vector additions and dot products are conducted modulo two.

7. A very trivial question, the .stab in genstab is that tau τ? A bit more context about this two-tuple representation ((dᵢ,, dⱼ), χ) given by .destabweights . In many places, I have observed that destabweights is used to check that whether a specific pair of (dᵢ, dⱼ) along with their associated weightχexists.

8. I find dᵢ: If D has r generators, we will often treat the index of dᵢ as an r-bit binary vector such that the factorization of dᵢ into the generators dᵢ of D is dᵢ =∏ (x=1 to r) (dₓ)^(iₓ). Is dⱼ follows a similar description? (please skip this)

9. In the algorithm 3, which I am currently assuming is apply, dᵢ corresponds to bₖ and dⱼ corresponds to bₗ after rowdecomponse. Also, instead of cₖ and cᵢ, it was combined as c = (dot(dₗˢᵗᵃᵇ,dᵢ) + dot(dᵣˢᵗᵃᵇ,dⱼ))*2, May I know the reason why? Also, I don't understand this -phaseₖ+phaseₗ+4 part of apply? I see that in algorithm 3, the phases are represented as αₕαₗ*?

I hope that I raised reasonable questions that you would be expecting from someone who is working on non-clifford. Maybe some of these questions are bad. I think that maybe addressing all the TODOs would be helpful to detect why is expect is giving warnings, and then proceed with completing functionality. I am looking for your feedback.

Kindly please take your time answering them. In the meantime, I will read the Ted's paper again, and will correspond reading with the implementation.

Thank you!

[Krastanov]

Is this for review or still WIP?

[Fe-r-oz]

Technically, it's WIP. I left it as open because I am looking for your comments on the TODOs. Thanks for your help.

In the meantime, I will turn it into draft.

[Krastanov]

Sorry, I will not be able to provide particularly detailed answers to this set of questions -- to some extent the point of the bounty is for the contributor to work through the paper and code on their own. Happy to attempt to answer the occasional brief question though.

[Fe-r-oz]

No problem. I totally understand. The set of questions helped me keep track of points, and eventually reading the paper frequently helped me refer to them :) It's like keeping a scratch space, although I should note them down locally on text document. I did realize that I needed to answer brief questions specific to one particular thing.... not essay questions!

[Krastanov]

indeed, the occasional brief question is much more probable to get a quick answer (and to not get lost)

[Fe-r-oz]

Kindly please have a look at implementation of Algorithm 2 and conj_destabs!. I hope you find it not soo bad. I tried to follow the paper as closely as possible. I really appreciate your guidance so that we can make it better.

I'm unsure if this is feasible, so I'd value your input before committing.

Algorithm 2 (page 8 of Ted's paper)

It covers both two cases for b where b is retrieved from rowdecompose such that P = a*d_b*s_c. The two cases are either (b = (0, 0, ... , 0) or b is nonzero.

I understand that atm, it's only updating χ to χ' not S to S' and D to D'.

I am unsure whether algorithm 2 is equivalent to your implementation of expect as it seems you considered a different approach for expect that differs from steps outlined of algorithm 2 (pg. 8).

Applied the same test of expect on algorithm 2, namely expect(Operator(p),ρ) == expect(p, gs), it seems worse.

function expect(p::PauliOperator, s::GeneralizedStabilizer)
    dict = s.destabweights
    dtype = _dictvaltype(dict)
    tzero = zero(dtype)
    tone = one(dtype)
    newdict = typeof(dict)(tzero)
    phase, b, c = rowdecompose(p, s.stab)
    if all(x -> x == 0, b)
        for ((dᵢ, dⱼ), χ) in dict
            if ((-tone)^(dot(dᵢ, c)))*im^phase == 1 && ((-tone)^(dot(dⱼ, c)))*im^phase == 1
                newdict[(dᵢ, dⱼ)] = get(newdict, (dᵢ, dⱼ), 0) + χ
            end
        end
        s.destabweights = newdict
        return s
    else
        k = _create_k(b)
        for ((dᵢ, dⱼ), χ) in dict
            q = 1
            dᵢ′ = copy(dᵢ)
            dⱼ′ = copy(dⱼ)
            if dot(dᵢ, k) == 1
                q *= ((-tone)^(dot(dᵢ, c)))*im^phase
                dᵢ′ = b .⊻ dᵢ
            end
            if dot(dⱼ, k) == 1
                q *= (-tone)^(dot(dⱼ, c))*conj(im^phase)
                dⱼ′ = b .⊻ dⱼ
            end
            newdict[(dᵢ′, dⱼ′)] = get(newdict, (dᵢ′,dⱼ′), 0) + χ * q / 2
        end
         s.destabweights = newdict
        return s
    end
end

function _create_k(b)
    k = zeros(Int, length(b))
    idx = findfirst(x -> x == 1, b)
    if idx !== nothing
        k[idx] = true
    end
    return k
end

Some runs:

julia> s = S"X"; p = P"Y";

julia> gs = GeneralizedStabilizer(s)
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
+ X
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + _ | + _

julia> apply!(gs, pcT)
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
+ X
with ϕᵢⱼ | Pᵢ | Pⱼ:
 0.0+0.353553im | + _ | + Z
 0.0-0.353553im | + Z | + _
 0.853553+0.0im | + _ | + _
 0.146447+0.0im | + Z | + Z

julia> expect(p, gs)
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
+ X
with ϕᵢⱼ | Pᵢ | Pⱼ:
 0.853553+0.0im | + _ | + _

Conjugate Destabs (conj_destabs!)

function conj_destabs!(state::GeneralizedStabilizer, gate::AbstractCliffordOperator)
    dict = state.destabweights
    dtype = _dictvaltype(dict)
    tzero = zero(dtype)
    newdict = typeof(dict)(tzero)
    gate_fields = fieldnames(typeof(gate))
    g_field = nothing
    for field in gate_fields
        if getfield(gate, field) isa Vector{UInt64}
            g_field = field
            break
        end
    end
    if g_field === nothing
        error("No Vector{UInt64} field found in the gate object")
    end
    g_part = getfield(gate, g_field)
    C = Vector{Bool}(g_part)
    for ((dᵢ, dⱼ), χ) in dict
        # Compute dᵢ' and dⱼ'
        dᵢ′ = C .* dᵢ .* C
        dⱼ′ = C .* dⱼ .* C
        newdict[(dᵢ′, dⱼ′)] = χ
    end
    state.destabweights = newdict
    state
end

Some Runs

julia> gs = GeneralizedStabilizer(S"-Y")
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
- Y
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + _ | + _

julia> apply!(gs, pcT)
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
- Y
with ϕᵢⱼ | Pᵢ | Pⱼ:
 0.0+0.353553im | + _ | + Z
 0.0-0.353553im | + Z | + _
 0.853553+0.0im | + _ | + _
 0.146447+0.0im | + Z | + Z

julia> QuantumClifford.conj_destabs!(gs, P"-X")
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
- Y
with ϕᵢⱼ | Pᵢ | Pⱼ:
 0.853553+0.0im | + _ | + _
 0.0+0.353553im | + _ | + Z
 0.0-0.353553im | + Z | + _
 0.146447+0.0im | + Z | + Z

julia> QuantumClifford.conj_destabs!(gs, P"X")
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
- Y
with ϕᵢⱼ | Pᵢ | Pⱼ:
 0.853553+0.0im | + _ | + _
 0.0+0.353553im | + _ | + Z
 0.0-0.353553im | + Z | + _
 0.146447+0.0im | + Z | + Z

julia> gs = GeneralizedStabilizer(S"Z")
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ X
𝒮𝓉𝒶𝒷
+ Z
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + _ | + _

julia> apply!(gs, pcT)
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ X
𝒮𝓉𝒶𝒷
+ Z
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + _ | + _

julia> QuantumClifford.conj_destabs!(gs, P"Y")
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ X
𝒮𝓉𝒶𝒷
+ Z
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + _ | + _

[Fe-r-oz]

Apologies for ALL the inconvenience. This PR is rewired to complete the following tasks from #260 , continuing the work from #259. These tasks are as follows:

• sm⊗sm
• tHadamard*sm
• tId1*pcT
• pcT⊗P"X"
• copy(sm)

Hence, dividing the nonclifford functionality into smaller PRs. This helps me not to get lost by taking definite smaller steps. Hopefully, this approach is satisfactory.

I will focus on the doctests, and tests for these, and will request review in a day or two.

Given that:

julia> sm = GeneralizedStabilizer(S"-X")
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
- X
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + _ | + _

Outputs are attached as follows:

sm⊗sm

julia> sm ⊗ sm
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z_
+ _Z
𝒮𝓉𝒶𝒷
- X_
- _X
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + __ | + __

tHadamard*sm

julia> tHadamard*sm
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ X
𝒮𝓉𝒶𝒷
- Z
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + _ | + _

pcT⊗P"X"

julia> pcT⊗P"X"
A unitary Pauli channel P = ∑ ϕᵢ Pᵢ with the following branches:
with ϕᵢ | Pᵢ
 0.853553+0.353553im | + _X
 0.146447-0.353553im | + ZX

copy(sm)

julia> copy(sm)
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
- X
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + _ | + _

tId1*pcT

julia> tId1 * pcT
A unitary Pauli channel P = ∑ ϕᵢ Pᵢ with the following branches:
with ϕᵢ | Pᵢ
 0.853553+0.353553im | X₁ ⟼ + X
Z₁ ⟼ + Z
 0.146447-0.353553im | X₁ ⟼ - X
Z₁ ⟼ + Z

julia> tHadamard * pcT
A unitary Pauli channel P = ∑ ϕᵢ Pᵢ with the following branches:
with ϕᵢ | Pᵢ
0.853553+0.353553im | X₁ ⟼ + Z
Z₁ ⟼ + X
0.146447-0.353553im | X₁ ⟼ + Z
Z₁ ⟼ - X

[Fe-r-oz]

The non-cliff interfaces mentioned in this comment: #316 (comment) are doctested along basic tests.

I have squashed all the commets related to these functionalities in this one commit for review: a6c7b14

[Krastanov]

@Fe-r-oz , you mentioned that the commits are squashed to make it more straightforward to review, but that does not seem to be the case. If you look at the list of commits, there is a lot of them, not organized in any particular way: https://github.com/QuantumSavory/QuantumClifford.jl/pull/316/commits

Could you make a new pull request from this branch towards the nonclif branch (i.e. towards this PR #259) so that it is easy to look only at the changes you have made (and please make sure the commits are organized).

[Fe-r-oz]

Apologies for this frustrating experience.

I only squashed the last 5 commits not all.

I was unable to squash all of them because of some commits Could not apply c64492b... Concatenated quantum code (#289) when try to squash them all.

[Fe-r-oz]

Closing this PR of it's making the review quite unfeasible.