beginning of biased pauli noise for PF sim (QuantumSavory#295)

Co-authored-by: Stefan Krastanov <[email protected]>
Co-authored-by: Stefan Krastanov <[email protected]>

CHANGELOG.md

@@ -5,9 +5,14 @@
 
 # News
 
+## v0.9.8 - 2024-07-24
+
+- `PauliError` can now encode biased noise during Pauli frame simulation, i.e. one can simulate only X errors, or only Y errors, or only Z errors, or some weighted combination of these.
+
 ## v0.9.7 - 2024-07-23
 
 - **(fix `#320`)** Fix a serious correctness bug in the SIMD implementation of Pauli string multiplication (affects the correctness of canonicalization and traceout for tableaux bigger than ~500 qubits; does not affect symbolic gates or Pauli frame simulations of any scale)
+
 ## v0.9.6 - 2024-07-12
 
 - `inv` implementation for single-qubit "symbolic" Clifford operators (subtypes of `AbstractSingleQubitOperator`).

src/noise.jl

@@ -28,6 +28,19 @@ struct UnbiasedUncorrelatedNoise{T} <: AbstractNoise
 end
 UnbiasedUncorrelatedNoise(p::Integer) = UnbiasedUncorrelatedNoise(float(p))
 
+"""Pauli noise model with probabilities `px`, `py`, and `pz` respectively for the three types of Pauli errors."""
+struct PauliNoise{T} <: AbstractNoise
+    px::T
+    py::T
+    pz::T
+end
+function PauliNoise(px::Real, py::Real, pz::Real)
+    px, py, pz = float.((px, py, pz))
+    px, py, pz = promote(px, py, pz)
+    T = typeof(px)
+    return PauliNoise{T}(px, py, pz)
+end
+
 """A convenient constructor for various types of Pauli noise models.
 Returns more specific types when necessary."""
 function PauliNoise end
@@ -38,7 +51,6 @@ function PauliNoise(p)
 end
 
 function applynoise!(s::AbstractStabilizer,noise::UnbiasedUncorrelatedNoise,i::Int)
-    n = nqubits(s)
     infid = noise.p/3
     r = rand()
     if r<infid
@@ -51,6 +63,18 @@ function applynoise!(s::AbstractStabilizer,noise::UnbiasedUncorrelatedNoise,i::I
     s
 end
 
+function applynoise!(s::AbstractStabilizer,noise::PauliNoise,i::Int)
+    r = rand()
+    if r<noise.px
+        apply_single_x!(s,i)
+    elseif r<noise.px+noise.pz
+        apply_single_z!(s,i)
+    elseif r<noise.px+noise.pz+noise.py
+        apply_single_y!(s,i)
+    end
+    s
+end
+
 """An operator that applies the given `noise` model to the qubits at the selected `indices`."""
 struct NoiseOp{N, Q} <: AbstractNoiseOp where {N, Q}
     noise::N #<:AbstractNoise
@@ -74,6 +98,22 @@ function PauliError(qubits,p)
     NoiseOp(PauliNoise(p), qubits)
 end
 
+""""Construct a gate operation that applies a biased Pauli error on  qubit `q` with independent probabilities `px`, `py`, `pz`.
+Note that the probability of any error occurring is `px+py+pz`. Because of this, `PauliError(1, p)` is equivalent to `PauliError(1,p/3,p/3,p/3)`.
+Similarly, if one wanted to exclude Z errors from `PauliError(1,p/3,p/3,p/3)` while mainting the same rate of X errors, one could write
+`PauliError(1, p*2/3, 0, 0)` (in the sense that Y errors can be interpreted as an X and a Z happening at the same time)."""
+function PauliError(q::Int, px, py, pz)
+    NoiseOp(PauliNoise(px,py,pz), (q,))
+end
+
+""""Construct a gate operation that applies a biased Pauli error on all `qubits` independently, each with  probabilities `px`, `py`, `pz`.
+Note that the probability of any error occurring is `px+py+pz`. Because of this, `PauliError(1, p)` is equivalent to `PauliError(1,p/3,p/3,p/3)`.
+Similarly, if one wanted to exclude Z errors from `PauliError(1,p/3,p/3,p/3)` while mainting the same rate of X errors, one could write
+`PauliError(1, p*2/3, 0, 0)` (in the sense that Y errors can be interpreted as an X and a Z happening at the same time)."""
+function PauliError(qubits, px, py, pz)
+    NoiseOp(PauliNoise(px,py,pz), qubits)
+end
+
 """An operator that applies the given `noise` model to all qubits."""
 struct NoiseOpAll <: AbstractNoiseOp
     noise::AbstractNoise

src/pauli_frames.jl

@@ -144,6 +144,28 @@ function applynoise!(frame::PauliFrame,noise::UnbiasedUncorrelatedNoise,i::Int)
     return frame
 end
 
+function applynoise!(frame::PauliFrame,noise::PauliNoise,i::Int)
+    T = eltype(frame.frame.tab.xzs)
+
+    lowbit = T(1)
+    ibig = _div(T,i-1)+1
+    ismall = _mod(T,i-1)
+    ismallm = lowbit<<(ismall)
+
+    @inbounds @simd for f in eachindex(frame)
+        r = rand()
+        if  r < noise.px # X error
+            frame.frame.tab.xzs[ibig,f] ⊻= ismallm
+        elseif r < noise.px+noise.pz # Z error
+            frame.frame.tab.xzs[end÷2+ibig,f] ⊻= ismallm
+        elseif r < noise.px+noise.pz+noise.py # Y error
+            frame.frame.tab.xzs[ibig,f] ⊻= ismallm
+            frame.frame.tab.xzs[end÷2+ibig,f] ⊻= ismallm
+        end
+    end
+    return frame
+end
+
 """
 Perform a "Pauli frame" style simulation of a quantum circuit.
 """

src/sumtypes.jl

@@ -217,7 +217,10 @@ function concrete_typeparams(t::Type{ClassicalXOR})
 end
 
 function concrete_typeparams(t::Type{NoiseOp})
-    return [(UnbiasedUncorrelatedNoise{Float64}, i) for i in 1:8]
+    return [
+        [(UnbiasedUncorrelatedNoise{Float64}, i) for i in 1:8];
+        [(PauliNoise{Float64}, i) for i in 1:8];
+    ]
 end