diff --git a/.gitignore b/.gitignore
index ebcf4e4b5..27086e2de 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,4 +3,5 @@ Manifest.toml
 LocalPreferences.toml
 */.*swp
 scratch/
-*.cov
\ No newline at end of file
+*.cov
+*.DS_Store
diff --git a/Project.toml b/Project.toml
index 1293bad0f..79a82305a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -12,12 +12,14 @@ HostCPUFeatures = "3e5b6fbb-0976-4d2c-9146-d79de83f2fb0"
 ILog2 = "2cd5bd5f-40a1-5050-9e10-fc8cdb6109f5"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LDPCDecoders = "3c486d74-64b9-4c60-8b1a-13a564e77efb"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 QuantumInterface = "5717a53b-5d69-4fa3-b976-0bf2f97ca1e5"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SIMD = "fdea26ae-647d-5447-a871-4b548cad5224"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SumTypes = "8e1ec7a9-0e02-4297-b0fe-6433085c89f2"
 
 [weakdeps]
diff --git a/src/ecc/ECC.jl b/src/ecc/ECC.jl
index d6b6de969..ff0837672 100644
--- a/src/ecc/ECC.jl
+++ b/src/ecc/ECC.jl
@@ -6,6 +6,8 @@ using QuantumClifford: AbstractOperation, AbstractStabilizer
 import QuantumClifford: Stabilizer, MixedDestabilizer
 using DocStringExtensions
 using Combinatorics: combinations
+using LDPCDecoders
+using SparseArrays
 
 abstract type AbstractECC end
 
@@ -291,6 +293,7 @@ function isdegenerate(H::Stabilizer, d::Int=1)
 end
 
 include("circuits.jl")
+include("pipelines.jl")
 
 include("codes/bitflipcode.jl")
 include("codes/fivequbit.jl")
diff --git a/src/ecc/pipelines.jl b/src/ecc/pipelines.jl
new file mode 100644
index 000000000..22e6d1967
--- /dev/null
+++ b/src/ecc/pipelines.jl
@@ -0,0 +1,512 @@
+"""Generate a lookup table for decoding single qubit errors."""
+function create_lookup_table(code::Stabilizer)
+    lookup_table = Dict()
+    constraints, qubits = size(code)
+    # In the case of no errors
+    lookup_table[ zeros(UInt8, constraints) ] = zero(PauliOperator, qubits)
+    # In the case of single bit errors
+    for bit_to_be_flipped in 1:qubits
+        for error_type in [single_x, single_y, single_z]
+            # Generate e⃗
+            error = error_type(qubits, bit_to_be_flipped)
+            # Calculate s⃗
+            # (check which stabilizer rows do not commute with the Pauli error)
+            syndrome = comm(error, code)
+            # Store s⃗ → e⃗
+            lookup_table[syndrome] = error
+        end
+    end
+    lookup_table
+end
+
+"""Given a syndrome circuit, returns the fault tolerant encoding circuit. Basically just a copy of the syndrome circuit that throws away the measurement results."""
+function fault_tolerant_encoding(scirc)
+    ecirc = Vector{QuantumClifford.AbstractOperation}()
+    for gate in scirc
+        if isa(gate,sMRZ)
+            push!(ecirc, sMZ(gate.qubit))
+        elseif isa(gate,sMRX)
+            push!(ecirc, sMX(gate.qubit))
+        elseif isa(gate,sMRY)
+            push!(ecirc, sMY(gate.qubit))
+        else
+            push!(ecirc, gate)
+        end
+    end
+    return ecirc
+end
+
+""" Most simple simulation function
+* Uses a lookup table decoder over single qubit errors ['create_lookup_table'](@ref)
+* Uses ['naive_syndrome_circuit'](@ref)
+Returns a logical X and Z error for the provided p_error - physical error rate per qubit
+"""
+function naive_error_correction_pipeline(checks::Stabilizer, p_error; nframes=10_000, ecirc=nothing, encoding_locs=nothing, scirc=nothing)
+    if isnothing(scirc)
+        scirc , _= naive_syndrome_circuit(checks)
+    end
+    lookup_table = create_lookup_table(checks)
+    O = faults_matrix(checks)
+    circuit_Z = Base.copy(scirc)
+    circuit_X = Base.copy(scirc)
+
+    s, n = size(checks)
+    k = n-s
+
+    if isnothing(encoding_locs)
+        pre_X = [sHadamard(i) for i in n-k+1:n]
+    else
+        pre_X = [sHadamard(i) for i in encoding_locs]
+    end
+    
+    md = MixedDestabilizer(checks)
+    logview_Z = [ logicalzview(md);]
+    logcirc_Z, _ = naive_syndrome_circuit(logview_Z) # numLogBits shoudl equal k
+
+    logview_X = [ logicalxview(md);]
+    logcirc_X, _ = naive_syndrome_circuit(logview_X)
+
+    # Z logic circuit
+    for gate in logcirc_Z
+        type = typeof(gate)
+        if type == sMRZ
+            push!(circuit_Z, sMRZ(gate.qubit+s, gate.bit+s))
+        else
+            push!(circuit_Z, type(gate.q1, gate.q2+s))
+        end
+    end
+
+   # X logic circuit
+    for gate in logcirc_X
+        type = typeof(gate)
+        if type == sMRZ
+            push!(circuit_X, sMRZ(gate.qubit+s, gate.bit+s))
+        else
+            push!(circuit_X, type(gate.q1, gate.q2+s))
+        end
+    end
+
+    # Z simulation
+    errors = [PauliError(i,p_error) for i in 1:n]
+    if isnothing(ecirc)
+        ecirc_z = fault_tolerant_encoding(circuit_Z) # double syndrome encoding
+    else
+        ecirc_z = ecirc
+    end
+    fullcircuit_Z = vcat(ecirc_z, errors, circuit_Z)
+
+    frames = PauliFrame(nframes, n+s+k, s+k)
+    pftrajectories(frames, fullcircuit_Z)
+    syndromes = pfmeasurements(frames)[:, 1:s]
+    logicalSyndromes = pfmeasurements(frames)[:, s+1: s+k]
+
+    decoded = 0
+    for i in 1:nframes
+        row = syndromes[i,:]
+        guess = get(lookup_table,row,nothing)
+        if isnothing(guess)
+            continue
+        else
+            result_Z = (O * stab_to_gf2(guess))[k+1:2k]
+            if result_Z == logicalSyndromes[i,:]
+                decoded += 1
+            end
+        end
+    end
+    z_error = 1 - decoded / nframes
+
+    # X simulation
+    if isnothing(ecirc)
+        ecirc_x = fault_tolerant_encoding(circuit_X)  # double syndrome encoding
+    else
+        ecirc_x = ecirc
+    end
+    fullcircuit_X = vcat(pre_X, ecirc_x, errors, circuit_X)
+    frames = PauliFrame(nframes, n+s+k, s+k)
+    pftrajectories(frames, fullcircuit_X)
+    syndromes = pfmeasurements(frames)[:, 1:s]
+    logicalSyndromes = pfmeasurements(frames)[:, s+1: s+k]
+
+    decoded = 0
+    for i in 1:nframes
+        row = syndromes[i,:]
+        guess = get(lookup_table,row,nothing)
+        if isnothing(guess)
+            continue
+        else
+            result_X = (O * stab_to_gf2(guess))[1:k]
+            if result_X == logicalSyndromes[i, :]
+                decoded += 1
+            end
+        end
+    end
+    x_error = 1 - decoded / nframes
+
+    return x_error, z_error
+end
+
+""" Similiar to ['naive_error_correction_pipeline'](@ref) but now uses shor style fault tolerant syndrome measurement
+* Uses a lookup table decoder over single qubit errors ['create_lookup_table'](@ref)
+* Uses ['shor_syndrome_circuit'](@ref)
+Returns a logical X and Z error for the provided p_error - physical error rate per qubit
+"""
+function shor_error_correction_pipeline(checks::Stabilizer, p_init; nframes=10_000, cat=nothing, scirc=nothing, ecirc=nothing, encoding_locs=nothing)
+    if isnothing(scirc) || isnothing(cat)
+        cat, scirc, _ = shor_syndrome_circuit(checks)
+    end
+
+    lookup_table = create_lookup_table(checks)
+    O = faults_matrix(checks)
+    circuit_Z = Base.copy(scirc)
+    circuit_X = Base.copy(scirc)
+
+    s, n = size(checks)
+    k = n-s
+
+    if isnothing(encoding_locs)
+        pre_X = [sHadamard(i) for i in n-k+1:n]
+    else
+        pre_X = [sHadamard(i) for i in encoding_locs]
+    end
+
+    anc_qubits = 0
+    for pauli in checks
+        anc_qubits += mapreduce(count_ones,+, xview(pauli) .| zview(pauli))
+    end
+    regbits = anc_qubits + s
+
+    md = MixedDestabilizer(checks)
+    logview_Z = logicalzview(md)
+    logcirc_Z, _ = naive_syndrome_circuit(logview_Z)
+
+    logview_X = logicalxview(md)
+    logcirc_X, _ = naive_syndrome_circuit(logview_X)
+
+    # Z logic circuit
+    for gate in logcirc_Z
+        type = typeof(gate)
+        if type == sMRZ
+            push!(circuit_Z, sMRZ(gate.qubit+anc_qubits, gate.bit+regbits))
+        else
+            push!(circuit_Z, type(gate.q1, gate.q2+anc_qubits))
+        end
+    end
+
+   # X logic circuit
+    for gate in logcirc_X
+        type = typeof(gate)
+        if type == sMRZ
+            push!(circuit_X, sMRZ(gate.qubit+anc_qubits, gate.bit+regbits))
+        else
+            push!(circuit_X, type(gate.q1, gate.q2+anc_qubits))
+        end
+    end
+
+    # Z simulation
+    errors = [PauliError(i,p_init) for i in 1:n]
+    if isnothing(ecirc)
+        ecirc_z = fault_tolerant_encoding(vcat(cat,circuit_Z)) # double syndrome encoding
+        fullcircuit_Z = vcat(ecirc_z, errors, circuit_Z)
+    else
+        ecirc_z = ecirc
+        fullcircuit_Z = vcat(ecirc_z, errors, cat, circuit_Z)
+    end
+    
+    frames = PauliFrame(nframes, n+anc_qubits+k, regbits+k)
+    pftrajectories(frames, fullcircuit_Z)
+    syndromes = pfmeasurements(frames)[:, anc_qubits+1:regbits]
+    logicalSyndromes = pfmeasurements(frames)[:, regbits+1:regbits+k]
+
+    decoded = 0
+    for i in 1:nframes
+        row = syndromes[i,:]
+        guess = get(lookup_table,row,nothing)
+        if isnothing(guess)
+            continue
+        else
+            result_Z = (O * stab_to_gf2(guess))[k+1:2k]
+            if result_Z == logicalSyndromes[i,:]
+                decoded += 1
+            end
+        end
+    end
+    z_error = 1 - decoded / nframes
+
+    # X simulation
+    if isnothing(ecirc)
+        ecirc_x = fault_tolerant_encoding(vcat(cat,circuit_X)) # double syndrome encoding
+        fullcircuit_X = vcat(pre_X, ecirc_x, errors, circuit_X)
+    else
+        ecirc_x = ecirc
+        fullcircuit_X = vcat(pre_X, ecirc_x, errors, cat, circuit_X)
+    end
+    
+    frames = PauliFrame(nframes, n+anc_qubits+k, regbits+k)
+    pftrajectories(frames, fullcircuit_X)
+    syndromes = pfmeasurements(frames)[:, anc_qubits+1:regbits]
+    logicalSyndromes = pfmeasurements(frames)[:, regbits+1:regbits+k]
+
+    decoded = 0
+    for i in 1:nframes
+        row = syndromes[i,:]
+        guess = get(lookup_table,row,nothing)
+        if isnothing(guess)
+            continue
+        else
+            result_X = (O * stab_to_gf2(guess))[1:k]
+            if result_X == logicalSyndromes[i, :]
+                decoded += 1
+            end
+        end
+    end
+    x_error = 1 - decoded / nframes
+
+    return x_error, z_error
+end
+
+struct CSS; tableau; cx; cz; end
+
+"""Naive syndrome measurement on a CSS ECC with a Cx and Cz matrix
+- Only wroks with fault tolerant encoding. 
+"""
+function CSS_naive_error_correction_pipeline(code::CSS, p_init; nframes=1_000, scirc=nothing, max_iters = 25)
+    if isnothing(scirc)
+        scirc , _= naive_syndrome_circuit(code.tableau)
+    end
+
+    O = faults_matrix(code.tableau)
+    circuit_Z = Base.copy(scirc)
+    circuit_X = Base.copy(scirc)
+
+    @assert size(code.cx, 2) == size(code.cz, 2) == nqubits(code.tableau)
+    @assert size(code.cx, 1) + size(code.cz, 1) == length(code.tableau)
+
+    s, n = size(code.tableau)
+    _, _, r = canonicalize!(Base.copy(code.tableau), ranks=true)
+    k = n - r
+
+    # Krishna decoder
+    log_probabs = zeros(n)
+    channel_probs = fill(p_init, n)
+
+    numchecks_X = size(code.cx)[1]
+    b2c_X = zeros(numchecks_X, n)
+    c2b_X = zeros(numchecks_X, n)
+
+    numchecks_Z = size(code.cz)[1]
+    b2c_Z = zeros(numchecks_Z, n)
+    c2b_Z = zeros(numchecks_Z, n)
+    err = zeros(n)
+
+    pre_X = [sHadamard(i) for i in n-k+1:n]
+
+    md = MixedDestabilizer(code.tableau)
+    logview_Z = logicalzview(md)
+    logcirc_Z, numLogBits_Z, _ = naive_syndrome_circuit(logview_Z)
+    @assert numLogBits_Z == k
+
+    logview_X = logicalxview(md)
+    logcirc_X, numLogBits_X, _ = naive_syndrome_circuit(logview_X)
+    @assert numLogBits_X == k
+
+    # Z logic circuit
+    for gate in logcirc_Z
+        type = typeof(gate)
+        if type == sMRZ
+            push!(circuit_Z, sMRZ(gate.qubit+s, gate.bit+s))
+        else
+            push!(circuit_Z, type(gate.q1, gate.q2+s))
+        end
+    end
+
+   # X logic circuit
+    for gate in logcirc_X
+        type = typeof(gate)
+        if type == sMRZ
+            push!(circuit_X, sMRZ(gate.qubit+s, gate.bit+s))
+        else
+            push!(circuit_X, type(gate.q1, gate.q2+s))
+        end
+    end
+
+    errors = [PauliError(i,p_init) for i in 1:n]
+
+    # Z simulation
+    ecirc = fault_tolerant_encoding(circuit_Z)
+    fullcircuit_Z = vcat(ecirc, errors, circuit_Z)
+
+    frames = PauliFrame(nframes, n+s+k, s+k)
+    pftrajectories(frames, fullcircuit_Z)
+    syndromes = pfmeasurements(frames)[:, 1:s]
+    logicalSyndromes = pfmeasurements(frames)[:, s+1: s+k]
+
+    decoded = 0
+    for i in 1:nframes
+        row = syndromes[i,:]
+        row_x = row[1:numchecks_X]
+        row_z = row[numchecks_X+1:numchecks_X+numchecks_Z]
+
+        KguessX, success = syndrome_decode(sparse(code.cx), sparse(code.cx'), row_x, max_iters, channel_probs, b2c_X, c2b_X, log_probabs, Base.copy(err))
+        KguessZ, success = syndrome_decode(sparse(code.cz), sparse(code.cz'), row_z, max_iters, channel_probs, b2c_Z, c2b_Z, log_probabs, Base.copy(err))
+        guess = vcat(KguessZ, KguessX)
+        
+        result_Z = (O * (guess))[k+1:2k]
+        if result_Z == logicalSyndromes[i,:]
+            decoded += 1
+        end
+    end
+    z_error = 1 - decoded / nframes
+
+    # X simulation
+    ecirc = fault_tolerant_encoding(circuit_X)
+    fullcircuit_X = vcat(pre_X, ecirc, errors, circuit_X)
+
+    frames = PauliFrame(nframes, n+s+k, s+k)
+    pftrajectories(frames, fullcircuit_X)
+    syndromes = pfmeasurements(frames)[:, 1:s]
+    logicalSyndromes = pfmeasurements(frames)[:, s+1: s+k]
+
+    decoded = 0
+    for i in 1:nframes
+        row = syndromes[i,:]
+        row_x = row[1:numchecks_X]
+        row_z = row[numchecks_X+1:numchecks_X+numchecks_Z]
+
+        KguessX, success = syndrome_decode(sparse(code.cx), sparse(code.cx'), row_x, max_iters, channel_probs, b2c_X, c2b_X, log_probabs, Base.copy(err))
+        KguessZ, success = syndrome_decode(sparse(code.cz), sparse(code.cz'), row_z, max_iters, channel_probs, b2c_Z, c2b_Z, log_probabs, Base.copy(err))
+        guess = vcat(KguessZ, KguessX)
+        
+        result_X = (O * (guess))[1:k]
+        if result_X == logicalSyndromes[i, :]
+            decoded += 1
+        end
+    end
+    x_error = 1 - decoded / nframes
+
+    return x_error, z_error
+end
+
+"""Shor syndrome measurement on a CSS ECC with a Cx and Cz matrix
+- Only wroks with fault tolerant encoding. 
+"""
+function CSS_shor_error_correction_pipeline(code::CSS, p_init;  nframes=10_000, cat=nothing, scirc=nothing, max_iters = 25)
+    if isnothing(scirc) || isnothing(cat)
+        cat, scirc, _ = shor_syndrome_circuit(code.tableau)
+    end
+
+    O = faults_matrix(code.tableau)
+    circuit_Z = Base.copy(scirc)
+    circuit_X = Base.copy(scirc)
+
+    @assert size(code.cx, 2) == size(code.cz, 2) == nqubits(code.tableau)
+    @assert size(code.cx, 1) + size(code.cz, 1) == length(code.tableau)
+
+    s, n = size(code.tableau)
+    _, _, r = canonicalize!(Base.copy(code.tableau), ranks=true)
+    k = n - r
+
+    # Krishna decoder
+    log_probabs = zeros(n)
+    channel_probs = fill(p_init, n)
+
+    numchecks_X = size(code.cx)[1]
+    b2c_X = zeros(numchecks_X, n)
+    c2b_X = zeros(numchecks_X, n)
+
+    numchecks_Z = size(code.cz)[1]
+    b2c_Z = zeros(numchecks_Z, n)
+    c2b_Z = zeros(numchecks_Z, n)
+    err = zeros(n)
+
+    pre_X = [sHadamard(i) for i in n-k+1:n]
+
+    anc_qubits = 0
+    for pauli in code.tableau
+        anc_qubits += mapreduce(count_ones,+, xview(pauli) .| zview(pauli))
+    end
+
+    regbits = anc_qubits + s
+
+    md = MixedDestabilizer(code.tableau)
+    logview_Z = logicalzview(md)
+    logcirc_Z, _ = naive_syndrome_circuit(logview_Z)
+
+    logview_X = logicalxview(md)
+    logcirc_X, _ = naive_syndrome_circuit(logview_X)
+
+    # Z logic circuit
+    for gate in logcirc_Z
+        type = typeof(gate)
+        if type == sMRZ
+            push!(circuit_Z, sMRZ(gate.qubit+anc_qubits, gate.bit+regbits))
+        else
+            push!(circuit_Z, type(gate.q1, gate.q2+anc_qubits))
+        end
+    end
+
+   # X logic circuit
+    for gate in logcirc_X
+        type = typeof(gate)
+        if type == sMRZ
+            push!(circuit_X, sMRZ(gate.qubit+anc_qubits, gate.bit+regbits))
+        else
+            push!(circuit_X, type(gate.q1, gate.q2+anc_qubits))
+        end
+    end
+
+    errors = [PauliError(i,p_init) for i in 1:n]
+
+    # Z simulation
+    ecirc = fault_tolerant_encoding(vcat(cat, circuit_Z))# notice that the ecirc now contains the cat state
+    fullcircuit_Z = vcat(ecirc, errors, circuit_Z)
+
+    frames = PauliFrame(nframes, n+anc_qubits+k, regbits+k)
+    pftrajectories(frames, fullcircuit_Z)
+    syndromes = pfmeasurements(frames)[:, anc_qubits+1:regbits]
+    logicalSyndromes = pfmeasurements(frames)[:, regbits+1:regbits+k]
+
+    decoded = 0
+    for i in 1:nframes
+        row = syndromes[i,:]
+        row_x = row[1:numchecks_X]
+        row_z = row[numchecks_X+1:numchecks_X+numchecks_Z]
+
+        KguessX, success = syndrome_decode(sparse(code.cx), sparse(code.cx'), row_x, max_iters, channel_probs, b2c_X, c2b_X, log_probabs, Base.copy(err))
+        KguessZ, success = syndrome_decode(sparse(code.cz), sparse(code.cz'), row_z, max_iters, channel_probs, b2c_Z, c2b_Z, log_probabs, Base.copy(err))
+        guess = vcat(KguessZ, KguessX)
+        
+        result_Z = (O * (guess))[k+1:2k]
+        if result_Z == logicalSyndromes[i,:]
+            decoded += 1
+        end
+    end
+    z_error = 1 - decoded / nframes
+
+    # X simulation
+    ecirc = fault_tolerant_encoding(vcat(cat, circuit_X))
+    fullcircuit_X = vcat(pre_X, ecirc, errors, circuit_X) # notice that the ecirc now contains the cat state
+
+    frames = PauliFrame(nframes, n+anc_qubits+k, regbits+k)
+    pftrajectories(frames, fullcircuit_X)
+    syndromes = pfmeasurements(frames)[:, anc_qubits+1:regbits]
+    logicalSyndromes = pfmeasurements(frames)[:, regbits+1:regbits+k]
+
+    decoded = 0
+    for i in 1:nframes
+        row = syndromes[i,:]
+        row_x = row[1:numchecks_X]
+        row_z = row[numchecks_X+1:numchecks_X+numchecks_Z]
+
+        KguessX, success = syndrome_decode(sparse(code.cx), sparse(code.cx'), row_x, max_iters, channel_probs, b2c_X, c2b_X, log_probabs, Base.copy(err))
+        KguessZ, success = syndrome_decode(sparse(code.cz), sparse(code.cz'), row_z, max_iters, channel_probs, b2c_Z, c2b_Z, log_probabs, Base.copy(err))
+        guess = vcat(KguessZ, KguessX)
+        
+        result_X = (O * (guess))[1:k]
+        if result_X == logicalSyndromes[i, :]
+            decoded += 1
+        end
+    end
+    x_error = 1 - decoded / nframes
+
+    return x_error, z_error
+end
\ No newline at end of file