Support phase with expression (#66)

* ready for testing gadget extraction

* support plotting circuit during extraction

* support extract ZX-diagrams with phase gadgets

* add full reduction

* fix function name typos

* update for Expr

* fix test

* fix conversion from ZXDiagram to QCircuit

* rework for `Phase`

* put type args into field

* update simplification

* update show

Author: ChenZhao44

src/ZXCalculus.jl

@@ -1,5 +1,6 @@
 module ZXCalculus
 
+include("phase.jl")
 include("abstract_zx_diagram.jl")
 include("zx_layout.jl")
 include("zx_diagram.jl")

src/circuit_extraction.jl

@@ -7,34 +7,42 @@ Extract circuit from a graph-like ZX-diagram.
 """
 function circuit_extraction(zxg::ZXGraph{T, P}) where {T, P}
     nzxg = copy(zxg)
-    nbits = nqubits(zxg)
+    nbits = nqubits(nzxg)
+    gads = Set{T}()
+    for v in spiders(nzxg)
+        if spider_type(nzxg, v) == SpiderType.Z && degree(nzxg, v) == 1
+            push!(gads, v, neighbors(zxg, v)[1])
+        end
+    end
 
-    cir = ZXDiagram(nbits)
+    # TODO: extract a QCircuit instead
+    # cir = QCircuit(nbits)
     Outs = get_outputs(nzxg)
     Ins = get_inputs(nzxg)
-    if length(Outs) != length(Ins)
-        return cir
-    end
     if nbits == 0
         nbits = length(Outs)
     end
+    cir = ZXDiagram(nbits)
+    if length(Outs) != length(Ins)
+        return cir
+    end
     for v1 in Ins
         @inbounds v2 = neighbors(nzxg, v1)[1]
         if !is_hadamard(nzxg, v1, v2)
             insert_spider!(nzxg, v1, v2)
         end
     end
     @inbounds frontier = [neighbors(nzxg, v)[1] for v in Outs]
-
-    extracted = copy(Outs)
+    qubit_map = Dict(zip(frontier, 1:nbits))
 
     for i = 1:nbits
         @inbounds w = neighbors(nzxg, Outs[i])[1]
         @inbounds if is_hadamard(nzxg, w, Outs[i])
             pushfirst_gate!(cir, Val{:H}(), i)
         end
-        pushfirst_gate!(cir, Val{:Z}(), i, phase(nzxg, w))
-        set_phase!(nzxg, w, zero(P))
+        if phase(nzxg, w) != 0
+            pushfirst_gate!(cir, Val{:Z}(), i, phase(nzxg, w))
+            set_phase!(nzxg, w, zero(P)) end
         @inbounds rem_edge!(nzxg, w, Outs[i])
     end
     for i = 1:nbits
@@ -47,11 +55,21 @@ function circuit_extraction(zxg::ZXGraph{T, P}) where {T, P}
             end
         end
     end
-    extracted = [extracted; frontier]
 
-    while !isempty(setdiff(spiders(nzxg), extracted))
-        frontier = update_frontier!(nzxg, frontier, cir)
-        extracted = union!(extracted, frontier)
+    old_frontier = copy(frontier)
+    max_iter = 1000
+    current_iter = 1
+    while !isempty(frontier)
+        update_frontier!(nzxg, gads, frontier, qubit_map, cir)
+        if frontier != old_frontier
+            old_frontier = copy(frontier)
+            current_iter = 1
+        else
+            current_iter += 1
+            if current_iter > max_iter
+                error("Circuit extraction failed!")
+            end
+        end
     end
 
     frontier = T[]
@@ -61,8 +79,8 @@ function circuit_extraction(zxg::ZXGraph{T, P}) where {T, P}
             push!(frontier, nb[])
         end
     end
-    sort!(frontier, by = (v->qubit_loc(nzxg, v)))
-    M = biadjancency(nzxg, frontier, Ins)
+    sort!(frontier, by = (v->qubit_map[v]))
+    M = biadjacency(nzxg, frontier, Ins)
     M, steps = gaussian_elimination(M)
     for step in steps
         if step.op == :addto
@@ -84,28 +102,65 @@ function circuit_extraction(zxg::ZXGraph{T, P}) where {T, P}
 end
 
 """
-    update_frontier!(zxg, frontier, cir)
+    update_frontier!(zxg, frontier, qubit_map, cir)
 
 Update frontier. This is a important step in the circuit extraction algorithm.
 For more detail, please check the paper [arXiv:1902.03178](https://arxiv.org/abs/1902.03178).
 """
-function update_frontier!(zxg::ZXGraph{T, P}, frontier::Vector{T}, cir::ZXDiagram{T, P}) where {T, P}
-    frontier = frontier[[spider_type(zxg, f) == SpiderType.Z && (degree(zxg, f)) > 0 for f in frontier]]
+function update_frontier!(zxg::ZXGraph{T, P}, gads::Set{T}, frontier::Vector{T}, qubit_map::Dict{T, Int}, cir::ZXDiagram{T, P}) where {T, P}
+    # TODO: use inplace methods
+    deleteat!(frontier, [spider_type(zxg, f) != SpiderType.Z || (degree(zxg, f)) == 0 for f in frontier])
+
+    for i = 1:length(frontier)
+        v = frontier[i]
+        nb_v = neighbors(zxg, v)
+        u = findfirst([u in gads for u in nb_v])
+        if u !== nothing
+            u = nb_v[u]
+            gad_u = zero(T)
+            for w in neighbors(zxg, u)
+                if w in gads
+                    gad_u = w
+                    break
+                end
+            end
+            rewrite!(Rule{:pivot}(), zxg, [u, gad_u, v])
+            pop!(gads, u)
+            pop!(gads, gad_u)
+            frontier[i] = u
+            qubit_map[u] = qubit_map[v]
+            pushfirst_gate!(cir, Val(:H), qubit_map[u])
+            delete!(qubit_map, v)
+            for j = 1:length(frontier)
+                for k = j+1:length(frontier)
+                    if is_hadamard(zxg, frontier[j], frontier[k])
+                        pushfirst_gate!(cir, Val(:CZ), qubit_map[frontier[j]], qubit_map[frontier[k]])
+                        rem_edge!(zxg, frontier[j], frontier[k])
+                    end
+                end
+            end
+
+            return frontier
+        end
+    end
+
     SetN = Set{T}()
     for f in frontier
         union!(SetN, neighbors(zxg, f))
     end
     N = collect(SetN)
-    sort!(N, by = v -> qubit_loc(zxg, v))
-    M = biadjancency(zxg, frontier, N)
+
+    # TODO: qubit_loc is not necessary here
+    # sort!(N, by = v -> qubit_loc(zxg, v))
+    M = biadjacency(zxg, frontier, N)
     M0, steps = gaussian_elimination(M)
     ws = T[]
     @inbounds for i = 1:length(frontier)
         if sum(M0[i,:]) == 1
             push!(ws, N[findfirst(isone, M0[i,:])])
         end
     end
-    M1 = biadjancency(zxg, frontier, ws)
+    # M1 = biadjacency(zxg, frontier, ws)
     @inbounds for e in findall(M .== 1)
         if has_edge(zxg, frontier[e[1]], N[e[2]])
             rem_edge!(zxg, frontier[e[1]], N[e[2]])
@@ -117,65 +172,60 @@ function update_frontier!(zxg::ZXGraph{T, P}, frontier::Vector{T}, cir::ZXDiagra
 
     @inbounds for step in steps
         if step.op == :addto
-            ctrl = qubit_loc(zxg, frontier[step.r2])
-            loc = qubit_loc(zxg, frontier[step.r1])
+            ctrl = qubit_map[frontier[step.r2]]
+            loc = qubit_map[frontier[step.r1]]
             pushfirst_gate!(cir, Val{:CNOT}(), loc, ctrl)
         else
-            q1 = qubit_loc(zxg, frontier[step.r1])
-            q2 = qubit_loc(zxg, frontier[step.r2])
+            q1 = qubit_map[frontier[step.r1]]
+            q2 = qubit_map[frontier[step.r2]]
 
             pushfirst_gate!(cir, Val{:CNOT}(), q2, q1)
             pushfirst_gate!(cir, Val{:CNOT}(), q1, q2)
             pushfirst_gate!(cir, Val{:CNOT}(), q2, q1)
         end
     end
-    old_frontier = copy(frontier)
-    @inbounds for w in ws
-        nb_w = neighbors(zxg, w)
-        v = intersect(nb_w, old_frontier)[1]
-        if (degree(zxg, v)) == 1
-            qubit_v = qubit_loc(zxg, v)
-            qubit_w = qubit_loc(zxg, w)
-            pushfirst_gate!(cir, Val{:H}(), qubit_v)
-            if spider_type(zxg, w) == SpiderType.Z
-                pushfirst_gate!(cir, Val{:Z}(), qubit_v, phase(zxg, w))
+
+    for i in 1:length(frontier)
+        v = frontier[i]
+        if degree(zxg, v) > 1
+            continue
+        end
+        w = neighbors(zxg, v)[1]
+        if is_hadamard(zxg, v, w)
+            pushfirst_gate!(cir, Val(:H), qubit_map[v])
+        end
+        if spider_type(zxg, w) == SpiderType.Z
+            qubit_map[w] = qubit_map[v]
+            if phase(zxg, w) != 0
+                pushfirst_gate!(cir, Val{:Z}(), qubit_map[w], phase(zxg, w))
                 set_phase!(zxg, w, zero(P))
             end
-            if qubit_v != qubit_w && spider_type(zxg, w) == SpiderType.Z
-                loc_v = column_loc(zxg, v)
-                loc_w = column_loc(zxg, w)
-                set_loc!(zxg.layout, w, qubit_v, loc_v)
-                set_column!(zxg.layout, v, loc_v+1//2)
-                # deleteat!(zxg.layout.spider_seq[qubit_w], loc_w)
-                # insert!(zxg.layout.spider_seq[qubit_v], loc_v, w)
-            end
             rem_edge!(zxg, v, w)
-            if spider_type(zxg, w) == SpiderType.In
-                add_edge!(zxg, w, v, EdgeType.SIM)
-            end
-            deleteat!(frontier, frontier .== v)
-            push!(frontier, w)
+        else
+            rem_edge!(zxg, v, w)
+            add_edge!(zxg, v, w, EdgeType.SIM)
         end
+        frontier[i] = w
     end
-    @inbounds for i1 = 1:length(ws)
-        for i2 = i1+1:length(ws)
-            if has_edge(zxg, ws[i1], ws[i2])
-                pushfirst_gate!(cir, Val{:CZ}(), qubit_loc(zxg, ws[i1]),
-                    qubit_loc(zxg, ws[i2]))
-                rem_edge!(zxg, ws[i1], ws[i2])
+
+    @inbounds for i1 = 1:length(frontier)
+        for i2 = i1+1:length(frontier)
+            if has_edge(zxg, frontier[i1], frontier[i2])
+                pushfirst_gate!(cir, Val{:CZ}(), qubit_map[frontier[i1]],
+                    qubit_map[frontier[i2]])
+                rem_edge!(zxg, frontier[i1], frontier[i2])
             end
         end
     end
-    sort!(frontier, by = v -> qubit_loc(zxg, v))
     return frontier
 end
 
 """
-    biadjancency(zxg, F, N)
+    biadjacency(zxg, F, N)
 
-Return the biadjancency matrix of `zxg` from vertices in `F` to vertices in `N`.
+Return the biadjacency matrix of `zxg` from vertices in `F` to vertices in `N`.
 """
-function biadjancency(zxg::ZXGraph{T, P}, F::Vector{T}, N::Vector{T}) where {T, P}
+function biadjacency(zxg::ZXGraph{T, P}, F::Vector{T}, N::Vector{T}) where {T, P}
     M = zeros(Int, length(F), length(N))
 
     for i = 1:length(F)

src/phase.jl

@@ -0,0 +1,108 @@
+import Base: +, -, *, /, ==, isless, rem, convert, zero, one, iseven
+import Base: show
+
+"""
+    Phase
+The type supports manipulating phases as expressions.
+"""
+struct Phase
+    ex
+    type
+end
+
+Phase(p::T) where {T} = Phase(p, T)
+
+Phase(p::Phase) = p
+
+function show(io::IO, p::Phase)
+    if p.ex isa Number
+        print(io, "Phase($(p.ex))")
+    else
+        print(io, "Phase($(p.ex)::$(p.type))")
+    end
+end
+
+function +(p1::Phase, p2::Phase)
+    T1 = p1.type
+    T2 = p2.type
+    if p1.ex isa Number && p2.ex isa Number
+        return Phase(p1.ex + p2.ex)
+    end
+
+    T = Base.promote_op(+, T1, T2)
+    return Phase(Expr(:call, :+, p1.ex, p2.ex), T)
+end
++(p1::Phase, p2::Number) = p1 + Phase(p2)
++(p1::Number, p2::Phase) = Phase(p1) + p2
+
+function -(p1::Phase, p2::Phase)
+    T1 = p1.type
+    T2 = p2.type
+    if p1.ex isa Number && p2.ex isa Number
+        return Phase(p1.ex - p2.ex)
+    end
+
+    T = Base.promote_op(-, T1, T2)
+    return Phase(Expr(:call, :-, p1.ex, p2.ex), T)
+end
+-(p1::Phase, p2::Number) = p1 - Phase(p2)
+-(p1::Number, p2::Phase) = Phase(p1) - p2
+
+function *(p1::Phase, p2::Phase)
+    T1 = p1.type
+    T2 = p2.type
+    if p1.ex isa Number && p2.ex isa Number
+        return Phase(p1.ex * p2.ex)
+    end
+
+    T = Base.promote_op(*, T1, T2)
+    return Phase(Expr(:call, :*, p1.ex, p2.ex), T)
+end
+*(p1::Phase, p2::Number) = p1 * Phase(p2)
+*(p1::Number, p2::Phase) = Phase(p1) * p2
+
+function /(p1::Phase, p2::Phase)
+    T1 = p1.type
+    T2 = p2.type
+    if p1.ex isa Number && p2.ex isa Number
+        return Phase(p1.ex / p2.ex)
+    end
+
+    T = Base.promote_op(/, T1, T2)
+    return Phase(Expr(:call, :/, p1.ex, p2.ex), T)
+end
+/(p1::Phase, p2::Number) = p1 / Phase(p2)
+/(p1::Number, p2::Phase) = Phase(p1) / p2
+
+
+function -(p::Phase)
+    T0 = p.type
+    if p.ex isa Number
+        return Phase(-p.ex)
+    end
+
+    T = Base.promote_op(-, T0)
+    return Phase(Expr(:call, :-, p.ex), T)
+end
+
+==(p1::Phase, p2::Phase) = p1.ex == p2.ex
+==(p1::Phase, p2::Number) = (p1 == Phase(p2))
+==(p1::Number, p2::Phase) = (Phase(p1) == p2)
+
+isless(p1::Phase, p2::Number) = (p1.ex isa Number) && p1.ex < p2
+function rem(p::Phase, d::Number)
+    if p.ex isa Number
+        return Phase(rem(p.ex, d))
+    end
+    return p
+end
+
+convert(::Type{Phase}, p) = Phase(p)
+convert(::Type{Phase}, p::Phase) = p
+
+zero(::Phase) = Phase(0//1)
+zero(::Type{Phase}) = Phase(0//1)
+one(::Phase) = Phase(1//1)
+one(::Type{Phase}) = Phase(1//1)
+
+iseven(p::Phase) = (p.ex isa Number) && (-1)^p.ex > 0