entropy.jl

using Random
using Yao
using Plots
using LinearAlgebra
using QuantumInformation

function run_circuit!(r::AbstractRegister, θs::Matrix, p::Real=0.0)
    n = nqubits(r)
    for j in 1:size(θs, 2) # each layer
        for i in 1:size(θs, 1)÷2 # each wire
            locs = (2i-1, mod1(2i, n))
            apply!(r, put(locs=>rot(kron(2, X, X), θs[i, j])))
        end

        if j != size(θs, 2)
            for i in 1:size(θs, 1)
                if rand() < p
                    measure!(r, i)
                end
            end
        end

        for i in 1:size(θs, 1)÷2 # each wire
            locs = (2i, mod1(2i+1, n))
            apply!(r, put(locs=>rot(kron(2, X, X), θs[i, j])))
        end

        if j != size(θs, 2)
            for i in 1:size(θs, 1)
                if rand() < p
                    measure!(r, i)
                end
            end
        end
    end
    return r
end

run_circuit(θs, p=0.0) = run_circuit!(zero_state(size(θs, 1)), θs, p)

function entropy_vn(r::ArrayReg{1}, A, B)
    ρ = @views r.state[:, 1] * r.state[:, 1]'
    n = nqubits(r)
    return vonneumann_entropy(ptrace(ρ, [1<<A, 1<<B], [1, ]))
end

function scaling(n, depth, p::Real=-1.0; seed=42)
    Random.seed!(seed)
    θs = randn(n, depth)
    r = run_circuit(θs, p)
    return [entropy_vn(r, A, n-A) for A in 1:n÷2]
end