reset_qubits! bug with MixedDestabilizer

Could be seen in comparing
```
mds = S"- YZ_
        - ZYZ
        + _ZY"
md = MixedDestabilizer(mds)
ns = S"+ XY
       - ZX"
nsp = S"+ XY_
        - ZX_"

tmd = traceout!(copy(md), [1,2])
project!(tmd, nsp[1])
project!(tmd, nsp[2])
QuantumClifford.mixed_destab_looks_good(nsp)

nmd = reset_qubits!(copy(md), ns,1:2)
QuantumClifford.mixed_destab_looks_good(nmd)
```

src/QuantumClifford.jl

@@ -1631,6 +1631,12 @@ Base.zero(::Type{<:CliffordOperator}, n) = CliffordOperator(zero(Stabilizer, 2n,
     p
 end
 
+@inline function zero!(s::Stabilizer,i)
+    s.xzs[:,i] .= zero(eltype(s.xzs))
+    s.phases[i] = 0x0
+    s
+end
+
 # TODO make faster by using fewer initializations, like in Base.zero above
 function Base.one(::Type{<:Stabilizer}, n; basis=:Z) # TODO support `basis` in all other `one(::[Mixed][De]Stabilizer)` functions
     if basis==:X

src/project_trace_reset.jl

@@ -452,9 +452,8 @@ Trace out a qubit.
 """ # TODO all of these should raise an error if length(qubits)>rank
 function traceout!(s::Stabilizer, qubits; phases=true, rank=false)
     _,i = canonicalize_rref!(s,qubits;phases=phases)
-    idpaulis = zero(PauliOperator,nqubits(s))
     for j in i+1:size(s,1)
-        s[j] = idpaulis # TODO - this can be done without creating/allocating an idpaulis object
+        zero!(s,j)
     end
     if rank return (s, i) else return s end
 end
@@ -481,18 +480,20 @@ end
 $TYPEDSIGNATURES
 
 Reset a given set of qubits to be in the state `newstate`.
+These qubits are traced out first, which could lead to "nonlocal" changes in
+the tableau.
 """
 function reset_qubits!(s::Stabilizer, newstate, qubits; phases=true)
-    # TODO raise error if sizes and length of qubits do not match
+    nqubits(newstate)==length(qubits) || throw(DimensionMismatch("`qubits` and `newstate` have to be of consistent size"))
+    length(qubits) <= nqubits(s) || throw(DimensionMismatch("the stabilizer is not big enough to contain the new state"))
     n = nqubits(s)
     s, x, z = canonicalize!(s,ranks=true) # TODO this is unnecessary, but it provides for nicely formatted tableaux; consider removing it for speed reasons
     _, rref_i = canonicalize_rref!((@view s[1:z]),qubits,phases=phases)
     for row in 1:length(newstate)
         s[row+rref_i] = _expand_pauli(newstate[row], qubits, n) # TODO do something that does not alocate temporary arrays
     end
-    idpaulis = zero(PauliOperator, n)
     for row in rref_i+length(newstate)+1:z
-        s[row] = idpaulis # TODO - this can be done without creating/allocating an idpaulis object
+        zero!(s,row)
     end
     s
 end
@@ -501,7 +502,8 @@ end
 $TYPEDSIGNATURES
 """
 function reset_qubits!(s::MixedStabilizer, newstate, qubits; phases=true) # TODO create the necessary interfaces so that Stabilizer and MixedStabilizer share this code
-    # TODO raise error if sizes and length of qubits do not match
+    nqubits(newstate)==length(qubits) || throw(DimensionMismatch("`qubits` and `newstate` have to be of consistent size"))
+    length(qubits) <= nqubits(s) || throw(DimensionMismatch("the stabilizer is not big enough to contain the new state"))
     n = nqubits(s)
     sv = stabilizerview(s)
     sv, rref_i = canonicalize_rref!(sv,qubits,phases=phases)
@@ -515,16 +517,17 @@ end
 """
 $TYPEDSIGNATURES
 """
-function reset_qubits!(s::MixedDestabilizer, newstate::Stabilizer, qubits; phases=true) # TODO this is really inefficient
-    # TODO raise error if sizes and length of qubits do not match
-    for pauli in newstate
+function reset_qubits!(s::MixedDestabilizer, newstate::AbstractStabilizer, qubits; phases=true) # TODO this is really inefficient
+    nqubits(newstate)==length(qubits) || throw(DimensionMismatch("`qubits` and `newstate` have to be of consistent size"))
+    length(qubits) <= nqubits(s) || throw(DimensionMismatch("the stabilizer is not big enough to contain the new state"))
+    newstatestab = stabilizerview(newstate)
+    traceout!(s,qubits)
+    for pauli in newstatestab
         expanded = _expand_pauli(pauli, qubits, nqubits(s)) # TODO, use a sparse project that does not require this expand
         _, anticomm, res = project!(s,expanded, phases=phases) # TODO make an `apply_measurement_phase!(project!(...), phase)`
         sv =  stabilizerview(s)
-        if anticomm!=0 # Does not commute with the stabilizer
+        if anticomm!=0 # Does not commute with the stabilizer or logical ops
             sv.phases[anticomm] = pauli.phase[]
-        elseif isnothing(res) # Is not in the stabilizer
-            sv.phases[s.rank] = pauli.phase[]
         else # Commutes with everyone
             if res!=0 && phases # TODO many of the checks below were already done by project!; find a way to not repeat them
                 destab = destabilizerview(s)

test/test_trace.jl

@@ -52,22 +52,27 @@ function test_trace()
         end
     end
     @testset "Qubit resets" begin
+        @test_throws DimensionMismatch reset_qubits!(S"XX YY",S"X",[1,2])
+        @test_throws DimensionMismatch reset_qubits!(S"X",S"XX YY",[1,2])
+        @test_throws DimensionMismatch reset_qubits!(MixedStabilizer(S"XX YY"),S"X",[1,2])
+        @test_throws DimensionMismatch reset_qubits!(MixedStabilizer(S"X"),S"XX YY",[1,2])
+        @test_throws DimensionMismatch reset_qubits!(MixedDestabilizer(S"XX YY"),S"X",[1,2])
+        @test_throws DimensionMismatch reset_qubits!(MixedDestabilizer(S"X"),S"XX YY",[1,2])
         for N in test_sizes
             for R in [rand(N÷2:N*2÷3), N]
-                if N<10
-                    @test_broken error("can not process empty stab")
-                    continue
-                end
+                R > 0 || continue
                 s = random_stabilizer(R,N)
-                newstate = random_stabilizer(rand(N÷4:N*2÷3))
+                nnew = rand(N÷4:N*2÷3)
+                nnew > 0 || continue
+                newstate = random_stabilizer(nnew)
                 perm = randperm(N)[1:nqubits(newstate)]
                 to_trace = setdiff(1:N,perm)
                 # Testing MixedDestabilizer
                 md = MixedDestabilizer(s)
                 mdr1 = reset_qubits!(copy(md), newstate,perm)
                 @test mixed_destab_looks_good(mdr1)
                 mdr2 = reset_qubits!(copy(mdr1),newstate,perm)
-                @test mdr1==mdr2
+                @test canonicalize!(copy(stabilizerview(mdr1)))==canonicalize!(copy(stabilizerview(mdr2)))
                 traceout!(mdr2,to_trace)
                 mdr2v = stabilizerview(mdr2)
                 @test canonicalize!(copy(mdr2v)[:,perm]) == canonicalize!(copy(newstate))
@@ -91,6 +96,8 @@ function test_trace()
                 c, x, z = canonicalize!(ssr2v, ranks=true)
                 @test canonicalize!(copy(ssr2v)[:,perm])[1:z] == canonicalize!(copy(newstate))
                 @test canonicalize!(msr2v) == c[1:z]
+                # Compare different datastractures
+                @test canonicalize!(copy(stabilizerview(mdr1)))==canonicalize!(copy(stabilizerview(msr1)))==canonicalize!(ssr1[1:mdr1.rank])
             end
         end
     end