diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index 93758d30b..48d059964 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -15,14 +15,8 @@ steps:
           codecov: true
     command:
       - echo "Julia depot path $${JULIA_DEPOT_PATH}"
-      - mkdir -p "$${JULIA_DEPOT_PATH}/conda/3/x86_64"
       - julia --project='~' -e '
         using Pkg;
-        ENV["PYTHON"] = "";
-        Pkg.add("Conda");
-        Pkg.add("PyCall");
-        Pkg.build("PyCall");
-        println("and now we dev the current package");
         pkg"dev .";'
   - label: "Downstream Breakage CI on Buildkite"
     plugins:
@@ -31,18 +25,11 @@ steps:
       - QuantumSavory/julia-xvfb#v1:
     command:
       - echo "Julia depot path $${JULIA_DEPOT_PATH}"
-      - mkdir -p "$${JULIA_DEPOT_PATH}/conda/3/x86_64"
       - julia --project='~' -e '
         using Pkg;
-        ENV["PYTHON"] = "";
-        Pkg.add("Conda");
-        Pkg.add("PyCall");
-        Pkg.build("PyCall");
-        println("and now we dev the current package");
         pkg"dev .";'
       - julia --project='~' -e '
         using Pkg;
-        ENV["PYTHON"] = "";
         Pkg.add("QuantumSavory");
         Pkg.build("QuantumSavory");
         Pkg.test("QuantumSavory");
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 400ec1e2a..6243fe140 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -47,9 +47,6 @@ jobs:
       - uses: julia-actions/setup-julia@v1
         with:
           version: '1'
-      - name: "Manually add Conda.jl folder due to bug in its build step"
-        run: |
-          mkdir -p "/home/runner/.julia/conda/3/x86_64"
       - uses: julia-actions/cache@v1
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-docdeploy@v1
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 04ad40e5c..8bb53e67e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,10 @@
 
 # News
 
+## v0.8.22 - 2024-03-15
+
+- Bump the `PyQDecoders` dependency, switching to `PythonCall` behind the scenes for reliability.
+
 ## v0.8.21 - 2024-03-05
 
 - Bump the `LDPCDecoders` dependency.
diff --git a/Project.toml b/Project.toml
index 0e10d045a..0d53428ed 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "QuantumClifford"
 uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 authors = ["Stefan Krastanov <stefan@krastanov.org>"]
-version = "0.8.21"
+version = "0.8.22"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -56,7 +56,7 @@ Makie = "0.20"
 Nemo = "0.38, 0.39, 0.40, 0.41, 0.42, 0.43"
 Plots = "1.38.0"
 PrecompileTools = "1.2"
-PyQDecoders = "0.1.1"
+PyQDecoders = "0.2.0"
 Quantikz = "1.3.1"
 QuantumInterface = "0.3.3"
 QuantumOpticsBase = "0.4.18"
diff --git a/ext/QuantumCliffordPyQDecodersExt/QuantumCliffordPyQDecodersExt.jl b/ext/QuantumCliffordPyQDecodersExt/QuantumCliffordPyQDecodersExt.jl
index b617dfcf5..a666c36c9 100644
--- a/ext/QuantumCliffordPyQDecodersExt/QuantumCliffordPyQDecodersExt.jl
+++ b/ext/QuantumCliffordPyQDecodersExt/QuantumCliffordPyQDecodersExt.jl
@@ -1,13 +1,14 @@
 module QuantumCliffordPyQDecodersExt
 
-using PyQDecoders: np, sps, ldpc, pm
+using PyQDecoders: np, sps, ldpc, pm, PythonCall
 using SparseArrays
 using QuantumClifford
 using QuantumClifford.ECC
 import QuantumClifford.ECC: AbstractSyndromeDecoder, decode, parity_checks
 
+abstract type PyBP <: AbstractSyndromeDecoder end
 
-struct PyBeliefPropDecoder <: AbstractSyndromeDecoder
+struct PyBeliefPropDecoder <: PyBP
     code
     H
     Hx
@@ -19,23 +20,47 @@ struct PyBeliefPropDecoder <: AbstractSyndromeDecoder
     pyz
 end
 
-function PyBeliefPropDecoder(c)
-    Hx = parity_checks_x(c) |> collect
-    Hz = parity_checks_z(c) |> collect
+struct PyBeliefPropOSDecoder <: PyBP
+    code
+    H
+    Hx
+    Hz
+    nx
+    nz
+    faults_matrix
+    pyx
+    pyz
+end
+
+function PyBeliefPropDecoder(c; maxiter=nothing)
+    Hx = parity_checks_x(c) |> collect # TODO should be sparse
+    Hz = parity_checks_z(c) |> collect # TODO should be sparse
     H = parity_checks(c)
     fm = faults_matrix(c)
-    pyx = ldpc.bp_decoder(Hx)
-    pyz = ldpc.bp_decoder(Hz)
+    max_iter=isnothing(maxiter) ? 0 : maxiter
+    pyx = ldpc.bp_decoder(np.array(Hx); max_iter) # TODO should be sparse
+    pyz = ldpc.bp_decoder(np.array(Hz); max_iter) # TODO should be sparse
     return PyBeliefPropDecoder(c, H, Hx, Hz, size(Hx, 1), size(Hz, 1), fm, pyx, pyz)
 end
 
-parity_checks(d::PyBeliefPropDecoder) = d.H
+function PyBeliefPropOSDecoder(c; maxiter=nothing)
+    Hx = parity_checks_x(c) |> collect # TODO should be sparse
+    Hz = parity_checks_z(c) |> collect # TODO should be sparse
+    H = parity_checks(c)
+    fm = faults_matrix(c)
+    max_iter=isnothing(maxiter) ? 0 : maxiter
+    pyx = ldpc.bposd_decoder(np.array(Hx); max_iter) # TODO should be sparse
+    pyz = ldpc.bposd_decoder(np.array(Hz); max_iter) # TODO should be sparse
+    return PyBeliefPropOSDecoder(c, H, Hx, Hz, size(Hx, 1), size(Hz, 1), fm, pyx, pyz)
+end
+
+parity_checks(d::PyBP) = d.H
 
-function decode(d::PyBeliefPropDecoder, syndrome_sample)
-    row_x = syndrome_sample[1:d.nx] # TODO This copy is costly!
+function decode(d::PyBP, syndrome_sample)
+    row_x = syndrome_sample[1:d.nx] # TODO These copies and indirections might be costly!
     row_z = syndrome_sample[d.nx+1:end]
-    guess_x = d.pyx.decode(row_x)
-    guess_z = d.pyz.decode(row_z)
+    guess_x = PythonCall.PyArray(d.pyx.decode(np.array(row_x)))
+    guess_z = PythonCall.PyArray(d.pyz.decode(np.array(row_z)))
     vcat(guess_z, guess_x)
 end
 
@@ -72,8 +97,8 @@ parity_checks(d::PyMatchingDecoder) = d.H
 function decode(d::PyMatchingDecoder, syndrome_sample)
     row_x = syndrome_sample[1:d.nx] # TODO This copy is costly!
     row_z = syndrome_sample[d.nx+1:end]
-    guess_z_errors = d.pyx.decode(row_x)
-    guess_x_errors = d.pyz.decode(row_z)
+    guess_z_errors = PythonCall.PyArray(d.pyx.decode(row_x))
+    guess_x_errors = PythonCall.PyArray(d.pyz.decode(row_z))
     vcat(guess_x_errors, guess_z_errors)
 end
 
diff --git a/src/ecc/ECC.jl b/src/ecc/ECC.jl
index 745eea428..060c65c4a 100644
--- a/src/ecc/ECC.jl
+++ b/src/ecc/ECC.jl
@@ -24,7 +24,7 @@ export parity_checks, parity_checks_x, parity_checks_z,
     CommutationCheckECCSetup, NaiveSyndromeECCSetup, ShorSyndromeECCSetup,
     TableDecoder,
     BeliefPropDecoder,
-    PyBeliefPropDecoder, PyMatchingDecoder
+    PyBeliefPropDecoder, PyBeliefPropOSDecoder, PyMatchingDecoder
 
 """Parity check tableau of a code.
 
diff --git a/src/ecc/decoder_pipeline.jl b/src/ecc/decoder_pipeline.jl
index 1602cecf0..53bda0971 100644
--- a/src/ecc/decoder_pipeline.jl
+++ b/src/ecc/decoder_pipeline.jl
@@ -212,21 +212,30 @@ end
 # From extensions:
 
 """A simple Belief Propagation decoder built around tools from `LDPCDecoders.jl`."""
-function BeliefPropDecoder(args...)
+function BeliefPropDecoder(args...; kwargs...)
     ext = Base.get_extension(QuantumClifford, :QuantumCliffordLDPCDecodersExt)
     if isnothing(ext)
         throw("The `BeliefPropDecoder` depends on the package `LDPCDecoders` but you have not installed or imported `LDPCDecoders` yet. Immediately after you import `LDPCDecoders`, the `BeliefPropDecoder` will be available.")
     end
-    return ext.BeliefPropDecoder(args...)
+    return ext.BeliefPropDecoder(args...; kwargs...)
 end
 
 """A Belief Propagation decoder built around tools from the python package `ldpc` available from the julia package `PyQDecoders.jl`."""
-function PyBeliefPropDecoder(args...)
+function PyBeliefPropDecoder(args...; kwargs...)
     ext = Base.get_extension(QuantumClifford, :QuantumCliffordPyQDecodersExt)
     if isnothing(ext)
         throw("The `PyBeliefPropDecoder` depends on the package `PyQDecoders` but you have not installed or imported `PyQDecoders` yet. Immediately after you import `PyQDecoders`, the `PyBeliefPropDecoder` will be available.")
     end
-    return ext.PyBeliefPropDecoder(args...)
+    return ext.PyBeliefPropDecoder(args...; kwargs...)
+end
+
+"""A Belief Propagation decoder with ordered statistics decoding, built around tools from the python package `ldpc` available from the julia package `PyQDecoders.jl`."""
+function PyBeliefPropOSDecoder(args...; kwargs...)
+    ext = Base.get_extension(QuantumClifford, :QuantumCliffordPyQDecodersExt)
+    if isnothing(ext)
+        throw("The `PyBeliefPropOSDecoder` depends on the package `PyQDecoders` but you have not installed or imported `PyQDecoders` yet. Immediately after you import `PyQDecoders`, the `PyBeliefPropOSDecoder` will be available.")
+    end
+    return ext.PyBeliefPropOSDecoder(args...; kwargs...)
 end
 
 """A perfect matching decoder built around tools from the python package `pymatching` available from the julia package `PyQDecoders.jl`."""
diff --git a/src/symbolic_cliffords.jl b/src/symbolic_cliffords.jl
index 4ccee98e7..0d69879a9 100644
--- a/src/symbolic_cliffords.jl
+++ b/src/symbolic_cliffords.jl
@@ -313,10 +313,10 @@ function CliffordOperator(op::AbstractTwoQubitOperator, n; compact=false)
         for (i,q) in ((1,op.q1),(2,op.q2))
             for (ii,qq) in ((1,op.q1),(2,op.q2))
                 c[q,qq] = _c[i,ii] # TODO define an `embed` helper function
-                c[n+q,qq] = _c[n+i,ii]
+                c[n+q,qq] = _c[2+i,ii]
             end
             c.tab.phases[q] = _c.tab.phases[i] # TODO define a `phasesview` or `phases` helper function
-            c.tab.phases[n+q] = _c.tab.phases[n+i]
+            c.tab.phases[n+q] = _c.tab.phases[2+i]
         end
         return c
     end
diff --git a/test/Project.toml b/test/Project.toml
index b320f66b1..aad3cdc34 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -4,7 +4,6 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
-Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
diff --git a/test/runtests.jl b/test/runtests.jl
index 33b3222e7..c26eaa33c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -24,7 +24,7 @@ macro doset(descr)
 end
 
 println("Starting tests with $(Threads.nthreads()) threads out of `Sys.CPU_THREADS = $(Sys.CPU_THREADS)`...")
-println("ENV[\"PYTHON\"] = \"$(get(ENV,"PYTHON",nothing))\"")
+
 
 # in order to run the gpu tests automatically set GPU_TESTS to true in the .env file
 if get(ENV, "GPU_TESTS", "") == "true"
diff --git a/test/test_ecc_decoder_all_setups.jl b/test/test_ecc_decoder_all_setups.jl
index c8d06253a..84d0d6bc2 100644
--- a/test/test_ecc_decoder_all_setups.jl
+++ b/test/test_ecc_decoder_all_setups.jl
@@ -2,16 +2,10 @@ using Test
 using QuantumClifford
 using QuantumClifford.ECC
 
-@testset "is Conda ok" begin
-    # Trigger Python install if required. Required for Buildkite CI!
-    import Conda
-    Conda.list()
-end
-
-# Run this only after checking Conda works
 import PyQDecoders
+import LDPCDecoders
 
-@testset "table decoder" begin
+@testset "table decoder, good for small codes" begin
     codes = [
         Steane7(),
         Shor9(),
@@ -29,18 +23,47 @@ import PyQDecoders
 
     for c in codes
         for s in setups
-            e = evaluate_decoder(TableDecoder(c), s, 100000)
-            #@show c
-            #@show s
-            #@show e
-            @assert max(e...) < noise/4
+            for d in [TableDecoder]
+                e = evaluate_decoder(d(c), s, 100000)
+                #@show c
+                #@show s
+                #@show e
+                @assert max(e...) < noise/4
+            end
+        end
+    end
+end
+
+##
+
+@testset "belief prop decoders, good for small codes" begin
+    codes = [
+    ]
+
+    noise = 0.001
+
+    setups = [
+        CommutationCheckECCSetup(noise/2),
+        NaiveSyndromeECCSetup(noise, 0),
+        ShorSyndromeECCSetup(noise, 0),
+    ]
+
+    for c in codes
+        for s in setups
+            for d in [c->PyBeliefPropOSDecoder(c, maxiter=10)]
+                e = evaluate_decoder(d(c), s, 100000)
+                @show c
+                @show s
+                @show e
+                @assert max(e...) < noise/4
+            end
         end
     end
 end
 
 ##
 
-@testset "matching decoder" begin
+@testset "matching decoder, good as long as column weight of the code is limited" begin
     codes = [
         Toric(8,8),
         Toric(9,9)
diff --git a/test/test_symcliff.jl b/test/test_symcliff.jl
index 8d6d190bf..9cfddb976 100644
--- a/test/test_symcliff.jl
+++ b/test/test_symcliff.jl
@@ -59,4 +59,10 @@ end
         s = random_destabilizer(2)
         @test sop*s == cop*s == tcop*s == ccop*s
     end
+    for op in subtypes(QuantumClifford.AbstractTwoQubitOperator)
+        sop = op(1,10)
+        op1 = CliffordOperator(sop,20)
+        op2 = sSWAP(2,10)*(CliffordOperator(sop,2; compact=true)⊗tensor_pow(tId1, 18))*CliffordOperator(sSWAP(2,10),20)
+        @test op1 == op2
+    end
 end