A bug on fusion gate

[FumiKobayashi]

I may find a bug case in TensorNetworkCodes.jl.
Under fusion gate of below, the type of pure_errors::Vector{Vector{Int}} is estimated as Union{nothing, Vector{Vector{int}}} instead of Vector{Vector{int}} and causes TypeError on _remove_qubits!(). The reason for it seems that _make_ready() is not working well.

# using Pkg
# Pkg.add("TensorNetworkCodes")
using TensorNetworkCodes

function init_state(num_seedtensor::Int)
    stabilizers::Vector{Vector{Int}} = [[3, 1, 0], [1, 3, 1], [0, 1, 3]]
    logicals::Vector{Vector{Int}} = []
    code = SimpleCode("3-body_graph_state", stabilizers, logicals)
    seedstate = SimpleCode("3-body_graph_state", stabilizers, logicals)
    for _ in 1:num_seedtensor-1
        code = combine(code, seedstate)
    end
    return code
end

"""
    hadamard(code::SimpleCode, index::Int) -> SimpleCode

Given a simple code with an index to apply a Hadamard gate, return a new simple code rotated a Hadamard gate.
"""
function hadamard(code::SimpleCode, index::Int)
    g = deepcopy(code.stabilizers)
    l = deepcopy(code.logicals)
    e = deepcopy(code.pure_errors)
    n = num_qubits(code)
    # precondition
    (index in 1:n) || error("logical qubit index out of bounds!")
    for element in g
        element[index] = _hadamard(element[index])
    end
    for element in l
        element[index] = _hadamard(element[index])
    end
    for element in e
        element[index] = _hadamard(element[index])
    end
    name = "$(index)th-Hadamard $(code.name)"
    return SimpleCode(name, g, l, e)
end
function _hadamard(p::Int)
    # precondition
    (p in 0:3) || error("logical qubit index out of bounds!")
    if p == 1
        return 3
    elseif p == 3
        return 1
    else
        return p
    end
end

# run
code = init_state(5)
new_code = hadamard(code, 2)
fused_new_code = fusion(new_code, [2, 3])

Then they return those errors.

MethodError: no method matching _remove_qubits!(::Vector{Union{Nothing, Vector{Int64}}}, ::Vector{Int64})

Closest candidates are:
  _remove_qubits!(!Matched::AbstractVector{<:AbstractVector{Int64}}, ::AbstractVector{Int64})
   @ TensorNetworkCodes ~/Documents/Tensor-Network-Codes/TensorNetworkCodes/TensorNetworkCodes/src/contraction_functions.jl:453


Stacktrace:
 [1] fusion(code::SimpleCode, qubit_pair::Vector{Int64})
   @ TensorNetworkCodes ~/Documents/Tensor-Network-Codes/TensorNetworkCodes/TensorNetworkCodes/src/contraction_functions.jl:219
 [2] top-level scope
   @ ~/Documents/Tensor-Network-Codes/TensorNetworkCodes/TensorNetworkCodes/nbs/bug_case.ipynb:3

I temporarlly solved this problem by commenting out codes related to pure_errors because SimpleCode can calculate correct pure_errors automatically. But it is better to modify it properly.

function fusion(code::SimpleCode, qubit_pair::AbstractVector{Int})
    stabilizers = deepcopy(code.stabilizers)
-   pure_errors = deepcopy(code.pure_errors)
+   # pure_errors = deepcopy(code.pure_errors)
    logicals = deepcopy(code.logicals)

    if length(useful_inds) != 0
        sort!(useful_inds)
        useful_ops = [stabilizers[ind] for ind in useful_inds]
        deleteat!(stabilizers, sort!(useful_inds))
-       deleteat!(pure_errors, sort!(useful_inds))
+       # deleteat!(pure_errors, sort!(useful_inds))

        stabilizers = _make_ready(stabilizers, useful_ops, qubit_pair)
        logicals = _make_ready(logicals, useful_ops, qubit_pair)
-       pure_errors = _make_ready(pure_errors, useful_ops, qubit_pair)
+       # pure_errors = _make_ready(pure_errors, useful_ops, qubit_pair)
    end

    _remove_qubits!(stabilizers, qubit_pair)
-   _remove_qubits!(pure_errors, qubit_pair)
+   # _remove_qubits!(pure_errors, qubit_pair)
    if length(logicals) != 0
        _remove_qubits!(logicals, qubit_pair)
    end

    if length(useful_inds) != 2
        new_stabilizers = Array{Int64,1}[]
        new_pure_errors = Array{Int64,1}[]
        for α in 1:length(stabilizers)
            if pauli_are_independent(vcat(new_stabilizers, [stabilizers[α]]))
                push!(new_stabilizers, stabilizers[α])
-               push!(new_pure_errors, pure_errors[α])
+               # push!(new_pure_errors, pure_errors[α])
            end
        end
-       pure_errors = new_pure_errors
+       # pure_errors = new_pure_errors
        stabilizers = new_stabilizers
    end

-   return SimpleCode(" ", stabilizers, logicals, pure_errors)
+   return SimpleCode(" ", stabilizers, logicals)
end