noncliff: implement project!

src/nonclifford.jl

@@ -119,15 +119,75 @@ function _dictvaltype(dict)
     return eltype(dict).parameters[2] # TODO there must be a cleaner way to do this
 end
 
+"""
+Projects the state of a [`GeneralizedStabilizer`](@ref) onto the eigenspaces of a
+given Pauli operator `p`. The projection is determined by the expectation value of
+`p` with respect to the [`GeneralizedStabilizer`](@ref).
+
+If the expectation value `⟨p⟩` is close to `1`, the state is projected onto the `+1`
+ eigenspace of `p`. Conversely, if `⟨p⟩` is near `−1`, the state is projected onto
+the `−1` eigenspace of `p`. For intermediate values of `⟨p⟩`, the projection is
+performed probabilistically, with the likelihood of projecting onto the `+1` or
+`−1` eigenspace proportional to the expectation value.
+
+```jldoctest
+julia> sm = GeneralizedStabilizer(S"-X")
+A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
+𝒟ℯ𝓈𝓉𝒶𝒷
++ Z
+𝒮𝓉𝒶𝒷
+- X
+with ϕᵢⱼ | Pᵢ | Pⱼ:
+ 1.0+0.0im | + _ | + _
+
+julia> apply!(sm, 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> project!(sm, P"-Y")[1]
+A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
+𝒟ℯ𝓈𝓉𝒶𝒷
+- X
+𝒮𝓉𝒶𝒷
++ Y
+with ϕᵢⱼ | Pᵢ | Pⱼ:
+ 0.0+0.353553im | + _ | - X
+ 0.0-0.353553im | - X | + _
+ 0.853553+0.0im | + _ | + _
+ 0.146447+0.0im | - X | - X
+```
+
+"""
 function project!(sm::GeneralizedStabilizer, p::PauliOperator)
     eval = expect(p, sm)
-    prob₁ = (real(eval)+1)/2
-    error("This functionality is not implemented yet")
+    prob = (real(eval)+1)/2
+    if prob ≈ 0
+        return _proj₋(sm, p)
+    elseif prob ≈ 1
+        return _proj₊(sm, p)
+    else
+        return rand() < prob ? _proj₊(sm, p) : _proj₋(sm, p)
+    end
 end
 
-function _proj₋(sm::GeneralizedStabilizer, p::PauliOperator)
-end
 function _proj₊(sm::GeneralizedStabilizer, p::PauliOperator)
+    newstab, anticom_idx, res = project!(sm.stab, p)
+    sm.stab = newstab
+    return sm, anticom_idx, res
+end
+
+function _proj₋(sm::GeneralizedStabilizer, p::PauliOperator)
+    newstab, anticom_idx, res = project!(sm.stab, -p)
+    sm.stab = newstab
+    return sm, anticom_idx, res
 end
 
 abstract type AbstractPauliChannel <: AbstractOperation end
@@ -172,6 +232,7 @@ function embed(n::Int,idx,pc::PauliChannel)
 end
 
 nqubits(pc::PauliChannel) = nqubits(pc.paulis[1][1])
+nqubits(sm::GeneralizedStabilizer) = nqubits(sm.stab)
 
 function apply!(state::GeneralizedStabilizer, gate::PauliChannel)
     dict = state.destabweights

test/test_nonclifford_quantumoptics.jl

@@ -49,4 +49,18 @@ qo_tgate.data[2,2] = exp(im*pi/4)
            end
         end
     end
-end
+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"]
+            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)
+        end
+    end
+end