diff --git a/src/nonclifford.jl b/src/nonclifford.jl
index 594cb2b05..2f69dcbd8 100644
--- a/src/nonclifford.jl
+++ b/src/nonclifford.jl
@@ -206,11 +206,11 @@ with ϕᵢⱼ | Pᵢ | Pⱼ:
 
 """
 function project!(sm::GeneralizedStabilizer, p::PauliOperator)
-    newstab, anticom_idx, res = project!(sm.stab, p)
+    newstab, anticom, res = project!(sm.stab, p)
     χ′ = expect(p, sm)
     n = nqubits(newstab)
     newsm = GeneralizedStabilizer(newstab, DefaultDict(0.0im, (falses(n),falses(n))=>χ′))
-    return newsm, anticom_idx, res
+    return newsm, anticom, res
 end
 
 """
@@ -225,21 +225,24 @@ function projectrand!(sm::GeneralizedStabilizer, p::PauliOperator)
 end
 
 function _proj₊(sm::GeneralizedStabilizer, p::PauliOperator)
-    newstab, anticom_idx, res = project!(sm.stab, p)
-    sm.stab = newstab
+    newstab, res = projectrand!(sm.stab, p)
+    χ′ = expect(p, sm)
     n = nqubits(newstab)
-    newsm = GeneralizedStabilizer(newstab, DefaultDict(0.0im, (falses(n),falses(n))=>expect(p, sm)))
-    return newsm, anticom_idx, res
+    newsm = GeneralizedStabilizer(newstab, DefaultDict(0.0im, (falses(n),falses(n))=>χ′))
+    return newsm, res
 end
 
 function _proj₋(sm::GeneralizedStabilizer, p::PauliOperator)
-    newstab, anticom_idx, res = project!(sm.stab, -p)
-    sm.stab = newstab
+    newstab, res = projectrand!(sm.stab, -p)
+    χ′ = expect(p, sm)
     n = nqubits(newstab)
-    newsm = GeneralizedStabilizer(newstab, DefaultDict(0.0im, (falses(n),falses(n))=>expect(p, sm)))
-    return newsm, anticom_idx, res
+    newsm = GeneralizedStabilizer(newstab, DefaultDict(0.0im, (falses(n),falses(n))=>χ′))
+    return newsm, res
 end
 
+Base.copy(sm::GeneralizedStabilizer) = GeneralizedStabilizer(copy(sm.stab),copy(sm.destabweights))
+Base.:(==)(sm₁::GeneralizedStabilizer, sm₂::GeneralizedStabilizer) = sm₁.stab==sm₂.stab && sm₁.destabweights==sm₂.destabweights
+
 abstract type AbstractPauliChannel <: AbstractOperation end
 
 """A Pauli channel datastructure, mainly for use with [`GeneralizedStabilizer`](@ref)
diff --git a/test/test_nonclifford_quantumoptics.jl b/test/test_nonclifford_quantumoptics.jl
index 445379f81..c89092411 100644
--- a/test/test_nonclifford_quantumoptics.jl
+++ b/test/test_nonclifford_quantumoptics.jl
@@ -52,15 +52,45 @@ qo_tgate.data[2,2] = exp(im*pi/4)
 end
 
 @testset "project!" begin
-    for s in [S"X", S"Y", S"Z", S"-X", S"-Y", S"-Z"]
-        for p in [P"X", P"Y", P"Z", P"-X", P"-Y", P"-Z"]
+    for n in 1:4
+        for repetition in 1:5
+            s = random_stabilizer(n)
+            p = random_pauli(n)
             gs = GeneralizedStabilizer(s)
-            apply!(gs, pcT)
-            ρ = Operator(gs)
-            gs_proj, _, _ = project!(gs, p)
-            ρ_proj = Operator(gs_proj)
-            @test isapprox(Operator(gs_proj), ρ_proj; atol=1e-5)
-            @test isapprox(expect(p, gs_proj), expect(Operator(p), ρ_proj); atol=1e-5)
+            for i in 1:rand(1:3)
+                apply!(copy(gs), embed(n, i, pcT))
+            end
+            qo_state = Operator(gs)
+            project!(copy(gs), copy(p))[1]
+            qo_state_after_proj = Operator(gs)
+            qo_pauli = Operator(gs)
+            qo_proj1 = (identityoperator(qo_pauli) - qo_pauli)/2
+            qo_proj2 = (identityoperator(qo_pauli) + qo_pauli)/2
+            result1  = qo_proj1*qo_state*qo_proj1'
+            result2  = qo_proj2*qo_state*qo_proj2'
+            @test qo_state_after_proj ≈ result2 || qo_state_after_proj ≈ result1
+        end
+    end
+end
+
+@testset "projectrand!" begin
+    for n in 1:4
+        for repetition in 1:5
+            s = random_stabilizer(n)
+            p = random_pauli(n)
+            gs = GeneralizedStabilizer(s)
+            for i in 1:rand(1:3)
+                apply!(copy(gs), embed(n, i, pcT))
+            end
+            qo_state = Operator(gs)
+            projectrand!(copy(gs), copy(p))[1]
+            qo_state_after_proj = Operator(gs)
+            qo_pauli = Operator(gs)
+            qo_proj1 = (identityoperator(qo_pauli) - qo_pauli)/2
+            qo_proj2 = (identityoperator(qo_pauli) + qo_pauli)/2
+            result1  = qo_proj1*qo_state*qo_proj1'
+            result2  = qo_proj2*qo_state*qo_proj2'
+            @test qo_state_after_proj ≈ result2 || qo_state_after_proj ≈ result1
         end
     end
 end