Hi

I have a problem calculating the partial trace of a mixed destabilizer. I tried

julia> state = MixedDestabilizer(random_destabilizer(5))
Rank 5 stabilizer
+ _Z_YX
+ _Z_XZ
- ZZ___
- YYX_Y
+ Z_YZY
═══════
- Z_ZYY
+ YXYYY
- _X_XX
- XXXXX
- __XZY
═══════

Then I try to trace out the first qubit

julia> traceout!(copy(state),[1])
Rank 3 stabilizer
- YYX_Y
+ Z_YZY
- ZZ___
━━━━━━━
+ ___ZY
+ YXXXX
━━━━━━━
- ___XX
- __XZY
- _X_XX
━━━━━━━
- Z_ZYY
+ XXX__

Why do I have a rank 3 stabilizer state rather than rank 2.

I think the stabilizers are

julia> stabilizerview(state)
- Z_ZYY
+ YXYYY
- _X_XX
- XXXXX
- __XZY

Only two of them are nonzero after I trace out the first qubit, right?

Thanks

Hi,

I just start learning the use of this package. My attempt is to use Clifford circuit to simulate a quantum circuit with random unitary evolution, projective measurement and amplitude damping. Is it possible to achieve this or at least some approximation of the amplitude damping channel? By the way, if this is not the correct place to ask this kind of question, I apologize for it. (Is there any forum for beginners of this package?)

Here, I have written a code for the implementation of dephasing. I am not sure if I have to make more modifications to my code in order to have a better speed. I am yet to benchmark it. The definition of dephasing is defined in page 4 of supplemental material of the paper https://arxiv.org/pdf/2202.12905.pdf. Please have a look:

function rowswapping(mixdstat,row1,row2)
    tabl1 = tab(mixdstat)
    if row1 == row2
        return tabl1
    else
        tabl1[row1], tabl1[row2] = tabl1[row2], tabl1[row1]
        numq = nqubits(mixdstat)
        tabl1[row1+numq], tabl1[row2+numq] = tabl1[row2+numq], tabl1[row1+numq]
        return tabl1
    end
end

function dephasingch(sysstate, qloc, syssize)
    orgrank = sysstate.rank
    prostate, indx, results = projectZ!(copy(sysstate), qloc) 
    if indx == 0 || indx > orgrank
        return sysstate
    else
        swappedtabl = rowswapping(sysstate, indx, orgrank)
        updatstate = MixedDestabilizer(swappedtabl, orgrank-1)
        return updatstate
    end
end

There are many ways to synthesize a circuit (with nothing bigger than 2-qubit gates) out of given Clifford Operation. They are also equivalent to a canonicalization procedure.

Track and implement them.

Implement sparse versions for each of the "assembler" gates.

Modify show so that too large tableaux do not pollute the CLI

Some julia projects have taken to the style of having each test file be a self contained module that can be copy pasted directly into a REPL. See for instance https://github.com/Krastanov/QuantumSavory.jl/blob/master/test/test_register_interface.jl which uses @safetestset. This style would make it easier to debug failing tests: currently we need to find what constants are defined and what libraries are imported outside of the test file, before rerunning the test interactively; with this change it would be easy to just copy-paste one single file.

I implemented two ways of calculating bipartite entanglement entropy base on this repo. One is by clipping algorithm. (Ref: Phys. Rev. B 100, 134306 (2019).) The other is by calculating the rank of bipartite adjacency matrix obtained by Graphs.Graph(s::AbstractStabilizer). (Ref: arXiv:quant-ph/0602096.) And I also used Qiskit to cross check my results, using "statevector_simulator" backend.

In my 800 test cases, clipping algorithm agrees with Qiskit simulation. But there are 30 cases that entropy values from graphs differ by 1 or 2. So I think this may be a result of wrong graphs with more or less edges between two subsystems.

To help reproduce the issue:

• Codes to calculate entanglement entropy and run tests: https://github.com/royess/QuantumClifford.jl/tree/dev-clipping
• Codes to generate Qiskit data: https://github.com/royess/test-stab-entanglement/blob/main/get_data.py
• My test data in HDF5: stab_ent.zip

I know the issue involves many implementation details. I am willing to do more checks if needed.

https://github.com/quantumlib/Stim

num_threads causes invalidations

using SnoopCompileCore
invalidations = @snoopr begin
    using QuantumClifford
    #QuantumClifford._precompile()
end
using SnoopCompile
##
length(uinvalidated(invalidations))
##
trees = invalidation_trees(invalidations)
##
ftrees = filtermod(QuantumClifford, trees) # ← a lot of num_threads issues
#tinf = @snoopi_deep QuantumClifford._my_precompile_()
#Core.Compiler.Timings.ROOT() with 54 direct children

I've noticed that, when trying to measure an operator which is in the logical group of a MixedDestabilizer, a BoundsError is thrown. I've tried to run the following minimal example:

using QuantumClifford
using QuantumClifford.Experimental.NoisyCircuits
using Test

state = Register(MixedDestabilizer(S"ZZ"), zeros(Bool, 1))
meas = DenseMeasurement(P"ZI", 1)

@test_throws nothing state, flag = applyop!(state, meas)

The result I see is:

Test Failed at REPL[98]:1
  Expression: (state, flag) = applyop!(state, op)
    Expected: nothing
      Thrown: BoundsError(UInt8[0x00, 0x00], (0,))
ERROR: There was an error during testing

I think that this is happening because this case involves a random measurement result, without the measured operator anticommuting with a stabilizer.

Taking a look around the codebase, I've noticed a few things. First, the BoundsError is coming from either line 436 or 438 of NoisyCircuits.jl, which attempts to set stabilizerview(stab).phases[anticom], where stab and anticom come from stab,anticom,r = project!(stab, op.pauli) on line 433.
In project!(d::MixedDestabilizer,pauli::PauliOperator;keep_result::Bool=true,phases::Bool=true) (beginning line 356, project_trace_reset.jl), there's some logic to determine the index of the stabilizer that pauli anticommutes with. There's a separate variable called anticomlog which gets set to j if pauli anticommutes with logical j (or something like this), but this index never gets returned.
This means that, when we get to the if isnothing(r) that determines whether we're going to put a phase on a row of the tableau, we've correctly determined that the measurement result is random, but the index of the new stabilizer that needs a phase update is set to 0, rather than the index of the new stabilizer which is proportional to pauli.
This is a little tricky to deal with, since it's not clear that the index of the anticommuting logical will be the same as the index of the new stabilizer that needs its phase set.

Weirdly, findfirst seems slower...

function f1(stabilizer, qubit)
    anticommutes = 0
    r = length(stabilizer)
    for i in 1:r # The explicit loop is faster than anticommutes = findfirst(row->comm(pauli,stabilizer,row)!=0x0, 1:r); both do not allocate.
        if isZ(stabilizer,i,qubit)
            anticommutes = i
            break
        end
    end
    return anticommutes
end

function f2(stabilizer, qubit)
    r = length(stabilizer)
    findfirst(i->isZ(stabilizer,i,qubit), 1:r)
end
##
N = 1000
s = ghz(N);
##
@benchmark f1(_s, N) setup=(_s=copy(s))
##
@benchmark f2(_s, N) setup=(_s=copy(s))

expose that functionality

Provide the Clifford operation and the circuit that enable a given logical operation for a given code

• by simple conjugation with the encoding circuit
• by implementing https://github.com/nrenga/symplectic-arxiv18a

Two of the Aqua code-quality tests have to be disabled. Why? Can they be fixed?

https://github.com/Krastanov/QuantumClifford.jl/blob/a8c856b438d62d06252f4b495b361b921bff6781/test/runtests.jl#L44

There are a number of broken tests in the unit tests (marked with @test_broken). Clean them up and either fix the issues or implement throw(...) error messages explaining why there is an issue (switching from @test_broken to @test_throws).

As seen in the benchmark log: https://github.com/Krastanov/QuantumCliffordBenchmarksLog
Fix it.

See https://www.nature.com/articles/nature03350 and https://quantum-journal.org/papers/q-2021-07-06-497/

As I mentioned in #57 , I have created a fork https://github.com/royess/QuantumClifford.jl/tree/dev-clipping to add clipping algorithm and entanglement entropy calculation. (Ref: Phys. Rev. B 100, 134306 (2019).) And I think it will be great if I can merge them into the original repo.

My implementations allow us to do following things:

1. Clipping algorithm. Briefly speaking, the algorithm can turn a stabilizer state into clipped gauge, where at each qubit the number of left and right endpoints of stabilizers stay are 2. And the two endpoints at each qubit is different Pauli operator. (Ref: Phys. Rev. X 7, 031016 (2017).)
2. Bigram. Based on clipping algorithm, we can collect a table of locations of endpoints, called the bigram. Bigram is unique for a stabilizer state. And the length distribution from bigram also gives some physical information in the first reference paper.
3. Bipartite entanglement entropy. Based on bigram, we can easily calculate bipartite entanglement entropy of a contiguous subregion.

The implementations are tested by random stabilizer states. For entanglement entropy, @Krastanov and I have cross checked with a pure-Julia and a Qiskit implementation.

Besides, as dicussed in #57 , @Krastanov and I proposed three ways (including clipping algorithm) to calculate entanglement entropy. And I think we need find a good way to arrange them.

If @Krastanov agrees, we may dicuss in details about a PR. Thanks!

In the function canonicalize_clip!, on line entanglement.jl, if you change

 k2 = findfirst(let j=j, k1=k1; k->
                           (|(tab[k,j]...) & # not identity
                           (tab[k,j]!=tab[k1,j])) end, # not same as k1
                           k1+1:columns)

to

 k2 = findfirst(let j=j, k1=k1; k->
                           (|(tab[k,j]...) & # not identity
                           (tab[k,j]!=tab[k1,j])) end, # not same as k1
                           k1+1:rows)

then the clipped gauge construction works for mixed states as well. Then the associated entanglement_entropy function needs to be modified (line 188):

length(subsystem_range) - count(r->(r[1] in subsystem_range && r[2] in subsystem_range, eachrow(bg))

Hi,

Is there a function which can calculate the entropy like Renyi entropy or von Neuman entropy? Do we have a function to calculate negativity? Suppose I want to divide the system into A and B. Then I want to trace out the part of B and use the partial density operator of A to calculate these quantities, how to do it?

Thanks

Using https://journals.aps.org/pra/pdf/10.1103/PhysRevA.73.022334 (check ref 17-19 as well)

good_bell_state = S"XX
                    ZZ"
canonicalize_rref!(good_bell_state)[1]

function state_tensor_pow(state::Stabilizer{Vector{UInt8}, Matrix{UInt64}},power::Integer)
    result=state
    for i in 1:power-1
        result=result⊗state
    end
    return result
end

function append_permutation_op!(circuit,permutation,qubit)
    if permutation==2
        gate=Hadamard⊗Hadamard
    elseif permutation==3
        gate=(Hadamard*Phase*Hadamard)⊗(Phase*Hadamard*Phase)
    elseif permutation==4
        gate=(Phase*Hadamard)⊗(Hadamard*Phase*Hadamard*Phase)
    elseif permutation==5
        gate=(Phase*Hadamard*Phase*Hadamard)⊗(Hadamard*Phase*Hadamard*Phase*Hadamard*Phase*Hadamard*Phase)
    elseif permutation==6
        gate=(Hadamard*Phase*Hadamard*Phase*Hadamard)⊗(Hadamard*Hadamard*Phase*Hadamard*Phase*Hadamard*Phase*Hadamard*Phase)
    end
    if permutation!=1
        circuit=append!(circuit,[SparseGate(gate,[qubit*2-1,qubit*2])])
    end
    return circuit
end 

function append_permutation_cnot!(circuit,permutation,pair)
    a,b=pair
    circuit=append_permutation_op!(circuit,permutation[1],a)
    circuit=append_permutation_op!(circuit,permutation[2],b)
    circuit=append!(circuit,[SparseGate(CNOT,[a*2-1,b*2-1]),SparseGate(CNOT,[a*2,b*2])])
    return circuit
end

function convert(abstract_circuit,netnoise_value,N)
    #convert abstract circuit representation to complete representation
    circuit=[]
    netnoise = UnbiasedUncorrelatedNoise(netnoise_value/N)
    for i in abstract_circuit
        if i[1]=="M"
            if i[2][1]==1
                measurement1=X
                measurement2=X
            elseif i[1][1]==2
                measurement1=Y
                measurement2=-Y
            else
                measurement1=Z
                measurement2=Z
            end
            for q in i[3]
                circuit=append!(circuit,[BellMeasurement([measurement1,measurement2],[q*2-1,q*2])])
                circuit=append!(circuit,[Reset(good_bell_state,[q*2-1,q*2])])
                circuit=append!(circuit,[NoiseOp(netnoise,[q*2-1,q*2])])
            end
        else
            circuit=append_permutation_cnot!(circuit,i[2],i[3])
        end
    end
    circuit=append!(circuit,[VerifyOp(good_bell_state, [1,2])])
    return circuit
end

make_noisy(g::SparseGate, noise) = NoisyGate(g, UnbiasedUncorrelatedNoise(1//3*noise))
make_noisy(m::BellMeasurement, noise) = NoisyBellMeasurement(m, noise)
make_noisy(other_op, noise) = other_op
make_noisy(circuit::AbstractVector, noise) = [make_noisy(op, noise) for op in circuit];

circuit=convert([("CN", [1, 6], [1, 3]), ("CN", [4, 1], [1, 2]), ("CN", [5, 3], [3, 2]), ("CN", [5, 5], [2, 3]), ("M", [1], [3]), ("CN", [5, 1], [1, 3]), ("CN", [6, 6], [3, 2]), ("M", [2], [3]), ("CN", [1, 3], [1, 2]), ("CN", [1, 1], [2, 3]), ("CN", [4, 6], [1, 2]), ("CN", [4, 4], [2, 3]), ("CN", [5, 4], [3, 2]), ("M", [2], [2, 3])],0.1,3)
N=3
netnoise_value=0.1
localnoise_value=0.01
initial_state=state_tensor_pow(good_bell_state,N)
netnoise = UnbiasedUncorrelatedNoise(netnoise_value/N)
netnoise_opall = NoiseOpAll(netnoise);
c = [netnoise_opall,circuit...]
c = make_noisy(c, localnoise_value)
@time pe_allnoise = petrajectories(initial_state, c)

output:

4.772920 seconds (8.90 M allocations: 517.571 MiB, 3.64% gc time, 43.27% compilation time)
Dict{Symbol, Float64} with 3 entries:
  :undetected_failure => 0.0220938
  :detected_failure   => 0.401384
  :true_success       => 0.15239

problem:
petrajectories takes long time when calling long circuits. This step takes the majority of time of the optimization circuits.

using JET; report_package("QuantumClifford") should be used in CI

When I first looked through the documentation, I found this example on the Manual page:
apply!(s, tCNOT, [4,2])

This became the basis for how I thought to use apply!(). For a while, I could not figure out how to use apply!() to apply a simple X gate to qubit (let's say qubit 2). The documentation says that sX is supported by apply!(), but using only the aforementioned example as a base, you'd get an error if you tried:
apply!(s, sX, [2])

I now know that I should do apply!(s, sX(2)), but this perhaps is unintuitive to some people. My suggestion is to either make applying various gates with apply!() uniform in how to provide arguments, or at least provide some more explanations and simple examples early in the documentation so that it's more clear how to apply gates.

The test coverage is only 80%. While it covers most of the "math" algorithms, it does not cover well the various warnings and error messages raised on inconsistent input. Both can be improved (but improving the Experimental modules is not worth it yet).

CliffordOperator(s::Stabilizer) = CliffordOperator(tab(MixedDestabilizer(s))) would be wrong because the row operations involved in creating a MixedDestabilizer changes the order of the qubits. Create a version of canonicalize_gott! that knows how to undo row operations and use it for the creation of CliffordOperator. Check that there are no lingering bugs in logdot.

For each stabilizer code there is a encoding and decoding Clifford operator. Provide functionality that gives these circuits.

In the code block

initial_state = good_bell_state⊗good_bell_state
circuit = [
    SparseGate(CPHASE, [1,3]),
    SparseGate(CPHASE, [2,4]),
    BellMeasurement([X,X], [3,4]),
    VerifyOp(good_bell_state, [1,2])
]
displaycircuit(circuit)

The following error occurs:

MethodError: no method matching length(::SparseGate)

Stacktrace:
[1] _similar_for(::UnitRange{Int64}, ::Type{Any}, ::SparseGate, ::Base.HasLength) at .\array.jl:597
[2] _collect(::UnitRange{Int64}, ::SparseGate, ::Base.HasEltype, ::Base.HasLength) at .\array.jl:630
[3] collect(::SparseGate) at .\array.jl:624
[4] broadcastable(::SparseGate) at .\broadcast.jl:682
[5] broadcasted at .\broadcast.jl:1255 [inlined]
[6] circuitwidth(::SparseGate) at C:\Users\lcg75.julia\packages\Quantikz\Ni7Ey\src\Quantikz.jl:321
[7] _broadcast_getindex_evalf at .\broadcast.jl:648 [inlined]
[8] _broadcast_getindex at .\broadcast.jl:621 [inlined]
[9] getindex at .\broadcast.jl:575 [inlined]
[10] copy at .\broadcast.jl:876 [inlined]
[11] materialize at .\broadcast.jl:837 [inlined]
[12] circuitwidth(::Array{AbstractOperation,1}) at C:\Users\lcg75.julia\packages\Quantikz\Ni7Ey\src\Quantikz.jl:321
[13] displaycircuit(::Array{AbstractOperation,1}; scale::Int64, kw::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\lcg75.julia\packages\Quantikz\Ni7Ey\src\Quantikz.jl:449
[14] displaycircuit(::Array{AbstractOperation,1}) at C:\Users\lcg75.julia\packages\Quantikz\Ni7Ey\src\Quantikz.jl:449
[15] top-level scope at In[38]:8
[16] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091

When I tried apply!(state,tHadamard,[2,3]), I got an error due to using [2,3]. This was my mistake, however afterwards, the state had changed. Perhaps checking for this error before modifying the state would allow for greater robustness?

This feature is useful for simulating Clifford circuit. I am trying to implement one.

Hi Stefan:

I just update QuantumClifford to v"0.5.2". But now I face a problem when I try to generate all 2-qubit gates by collect(enumerate_cliffords(2)). The error I got is

julia> cliffords = collect(enumerate_cliffords(2))
ERROR: MethodError: no method matching var"#rowswap!#43"(::Bool, ::typeof(QuantumClifford.rowswap!)                                                             , ::Stabilizer{Vector{UInt8}, Matrix{UInt64}}, ::Int64, ::Int64)
Closest candidates are:
  var"#rowswap!#43"(::Val{B}, ::typeof(QuantumClifford.rowswap!), ::Stabilizer, ::Any, ::Any) where                                                              B at ~/.julia/packages/QuantumClifford/jpKJt/src/QuantumClifford.jl:895
Stacktrace:
 [1] symplecticGS(pauli::PauliOperator{Array{UInt8, 0}, Vector{UInt64}}; padded_n::Int64)
   @ QuantumClifford ~/.julia/packages/QuantumClifford/jpKJt/src/enumeration.jl:107
 [2] enumerate_cliffords_slow(n::Int64, i::BigInt; padded_n::Int64, onlycoset::Bool)
   @ QuantumClifford ~/.julia/packages/QuantumClifford/jpKJt/src/enumeration.jl:133
 [3] #enumerate_cliffords#148
   @ ~/.julia/packages/QuantumClifford/jpKJt/src/enumeration.jl:124 [inlined]
 [4] #152
   @ ./none:0 [inlined]
 [5] iterate
   @ ./generator.jl:47 [inlined]
 [6] collect(itr::Base.Generator{UnitRange{BigInt}, QuantumClifford.var"#152#154"{Int64, Int64}})
   @ Base ./array.jl:724
 [7] top-level scope
   @ REPL[6]:1

Is anything changed in the new version of QuantumClifford?

Thanks.

• https://www.cgranade.com/python-quaec/
• https://qecsim.github.io/
• https://github.com/qiskit-community/qiskit-qec
• https://panqec.readthedocs.io/en/latest/index.html
• https://github.com/PECOS-packages/PECOS

This can be based around package extensions and PythonCall interfaces, similar to how we are doing python-based decoders. E.g. see how we have implemented the pymatching decoder interface:

• In ECC/decoder_pipeline.jl

  QuantumClifford.jl/src/ecc/decoder_pipeline.jl

  Line 272 in 64f695c

  """A perfect matching decoder built around tools from the python package `pymatching` available from the julia package `PyQDecoders.jl`."""

• actually implemented in an extensions

  QuantumClifford.jl/ext/QuantumCliffordPyQDecodersExt/QuantumCliffordPyQDecodersExt.jl

  Line 92 in 64f695c

  function PyMatchingDecoder(c; weights=nothing)

• with the extension being dependent on this library that uses PythonCall to install what is needed https://github.com/QuantumSavory/PyQDecoders.jl

As a newbie, I have some problems understanding a few of the descriptions. Some of the methods could be documented based on their actions and the expected results.

• the MixedDestabilizer could be phrased as:
  "Keeping track of the r constraints over n qubits, resulting in r qubits that trigger check on the unique corresponding constraint in the matrix, and (n-r) qubits that do not trigger check in any of the constraints"

• the sMZ gate, make_noisy, and UnbiasedUncorrelatedNoise can also be defined more specifically

It seems that Quantiks does not compile because Tectronic does not compile on OSX M1 chip (arm64)

      Info Packages marked with ⌃ and ⌅ have new versions available, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m`
    Building GR → `~/.julia/scratchspaces/44cfe95a-1eb2-52ea-b672-e2afdf69b78f/cf0a9940f250dc3cb6cc6c6821b4bf8a4286cf9c/build.log`
Precompiling project...
  ✗ QuantumClifford
  ✗ QuantumCliffordPlots
  71 dependencies successfully precompiled in 56 seconds. 366 already precompiled. 2 skipped during auto due to previous errors.
  2 dependencies errored. To see a full report either run `import Pkg; Pkg.precompile()` or load the packages
[ Info: Precompiling QuantumClifford [0525e862-1e90-11e9-3e4d-1b39d7109de1]
ERROR: LoadError: Cannot locate artifact 'tectonic_bin' for aarch64-apple-darwin-libgfortran5-cxx11-julia_version+1.8.

some is already implemented, but document the use of different approaches:

• https://arxiv.org/abs/1406.2170
• https://arxiv.org/pdf/2003.09412.pdf

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

E.g. like here https://github.com/ali-ramadhan/DocumenterCitations.jl/blob/master/docs/example.bib

E.g. something like sMX, sMZ, and sMY which serve as ways to measure/project a qubit without employing the heavier project!.

Currently the code branches a lot on the size of the operator being generated (using ints for n<30, float for n<500, bigint ootherwise). ILog2 takes care of some of that branching. See how much you can simplify the code by using it. See #47 (comment)

random_destabilizer causes high TTFX

using SnoopCompileCore
invalidations = @snoopr begin
    using QuantumClifford
    #QuantumClifford._precompile()
end
using SnoopCompile
##
length(uinvalidated(invalidations))
##
trees = invalidation_trees(invalidations)
##
ftrees = filtermod(QuantumClifford, trees)
##
tinf = @snoopi_deep QuantumClifford._my_precompile_() # InferenceTimingNode: 6.178558/22.063535 on Core.Compiler.Timings.ROOT() with 54 direct children
##
staletrees = precompile_blockers(trees, tinf) 
##
@time using SnoopCompile
@time using QuantumClifford
@time using ProfileView
##
fg = flamegraph(tinf)
ProfileView.view(fg)

The large bunch at the start is due to random_destabilizer.

Fix all warnings and errors in the documentation and add a test so that the CI warns if there are new documentation issues.

Hi, I want to use this package to simulate Clifford circuits. And I want to know how to compute the inner product of two stablizer states with this package.

using QuantumClifford
using QuantumClifford: mixed_stab_looks_good, mixed_destab_looks_good
using Random
using LinearAlgebra
##
test_sizes = [10]
for N in test_sizes, lala in 1:100
    for R in [rand(N÷2:N*2÷3), N]
        if N<10
            @test_broken error("can not process empty stab")
            continue
        end
        s = random_stabilizer(R,N)
        newstate = random_stabilizer(rand(N÷4:N*2÷3))
        perm = randperm(N)[1:nqubits(newstate)]
        # Testing MixedDestabilizer
        md = MixedDestabilizer(s)
        mdr1 = reset_qubits!(copy(md), newstate,perm)
        if !mixed_destab_looks_good(mdr1)
            error("bug")
        end
    end
end

caught by the CI

Stim is great at:

• symbolic left multiplication
• symbolic right multiplication
• pauli frames #32

We have similar performance on some of these topics, but a complete comparison would be useful

• stim comparison #16

See comments in d3e84c1#commitcomment-67183201 for code references

The standard count_ones function should be used instead of our custom written xor_bits function.

Hi Stefan:

I am thinking about generating all 24 single-qubit Clifford gates in matrix form. Is there a function to do that? Thanks

Hi,

How to get a maximally mixed state? I think this is a stabilizer state, right? The code state = one(MixedDestabilier,n,n) gives me a pure state with all spins up, right? (I apologize if I made mistakes. I am still slowly learning this field). Also, which reference talks about the destabilizers of mixed state? I understand the pure state case but for the mixed state case, it is not quite clear to me. Thanks

Like sCNOT and sSWAP in https://github.com/Krastanov/QuantumClifford.jl/blob/master/src/symbolic_cliffords.jl