3. Utility Example

using PauliPropagation
using Base.Threads
using Plots
IBM_mitigated_vals = [1.01688859, 1.00387483, 0.95615886, 0.95966435, 0.83946763,
    0.81185907, 0.54640995, 0.45518584, 0.19469377, 0.01301832,0.01016334]
IBM_angles = [0.    , 0.1   , 0.2   , 0.3   , 0.4   , 0.5   , 0.6   , 0.7   , 0.8   , 1.    , 1.5707]

tn_vals = [9.99999254e-01,  9.99593653e-01,  9.95720077e-01,  9.88301532e-01,
        9.78553511e-01,  9.58023054e-01,  9.21986059e-01,  8.81726079e-01,
        8.49816779e-01,  8.24900527e-01,  7.91257641e-01,  7.37435202e-01,
        6.68573798e-01,  5.88096040e-01,  4.81874079e-01,  3.50316579e-01,
        2.26709331e-01,  1.39724659e-01,  7.86639143e-02,  4.24124371e-02,
        1.90595136e-02,  6.18879050e-03, -8.27168956e-04, -4.63372099e-03,
       -7.05202121e-03, -7.68387421e-03, -6.33121142e-03, -4.32594440e-03,
        6.52050191e-04,  1.72598340e-04,  5.64696020e-05, -7.70582375e-07]
tn_angles = LinRange(0, π/2, length(tn_vals));


google_vals = [1, 0.9957248494988796,0.9785691808859784,0.9227771228638144,0.8549387301071374,
    0.7790027681081058,0.6093954499261214,0.4257824480390532,0.2090085348069826,0.0116245678393245, -4.8872551798147e-09]
google_angles = [0.000, 0.100, 0.200, 0.300, 0.400,0.500,0.600,0.700,0.800,1.000,1.571]
function kickedisingcircuit(nq, nl; topology=nothing)
    
    # just in case a topology is not provided
    if isnothing(topology)
        topology = bricklayertopology(nq)
    end
    
    # define a layer of parametrized Rx-gates
    xlayer(circuit) = append!(circuit, (PauliRotation([:X], [qind]) for qind in 1:nq))
    
    # define a layer of fixed RZZ(π/2)-gates
    zzlayer(circuit) = append!(circuit, (PauliRotation([:Z, :Z], pair, -π/2) for pair in topology))
    
    # define the empty circuit
    circuit = Gate[]
    # append to the circuit
    for _ in 1:nl
        xlayer(circuit)
        zzlayer(circuit)
    end

    return circuit
end
kickedisingcircuit (generic function with 1 method)
# number of qubits
nq = 127
# the IBM Eagle topology
topology = ibmeagletopology
# 20 layers for the hardest simulation
nl = 20

# create the circuit
circuit = kickedisingcircuit(nq, nl; topology);
# count the number of parameters. Here only the RX gates, the RZZ are frozen to π/2
nparams = countparameters(circuit)

# define the observable (note the index is 62 if you start with 0)
pstr = PauliString(nq, :Z, 63)
PauliString(nqubits: 127, 1.0 * IIIIIIIIIIIIIIIIIIII...)
# perform one simulation
# the reason we make this into a function is so that the Julia garbage collector gets rid of `psum` after we leave the function
# otherwise it floats around in the global name space and takes up RAM
function onesimulation(circuit, pstr, thetas; min_abs_coeff=0, max_weight=Inf)
    # propagate
    psum = propagate(circuit, pstr, thetas; min_abs_coeff, max_weight)
    # overlap with the zero-state
    return overlapwithzero(psum)
end
onesimulation (generic function with 1 method)
# set the truncations
min_abs_coeff = 1e-4
max_weight = 8

# prepare everything and run
x_angles = LinRange(0, π/2, 20)
expectations = zeros(length(x_angles))
# multi-thread over the different x_angles
t = @timed @threads for ii in eachindex(x_angles)
    angle = x_angles[ii]
    thetas = ones(nparams) * angle
    @time expectations[ii] = onesimulation(circuit, pstr, thetas; min_abs_coeff, max_weight)
end

print("Total simulation time: ", t.time, " seconds")
  3.502918 seconds (2.19 M allocations: 105.969 MiB, 1.33% gc time, 99.89% compilation time)
  0.023115 seconds (21.07 k allocations: 839.367 KiB)


  0.088549 seconds (21.20 k allocations: 1022.062 KiB)


  0.171156 seconds (21.23 k allocations: 1.639 MiB)


  0.511400 seconds (21.25 k allocations: 2.100 MiB)


  1.087727 seconds (21.30 k allocations: 3.263 MiB)


  1.948378 seconds (21.36 k allocations: 6.151 MiB)


  3.184862 seconds (21.43 k allocations: 5.514 MiB)


  7.520723 seconds (21.71 k allocations: 18.348 MiB)


  6.096936 seconds (21.62 k allocations: 13.216 MiB)


  7.174851 seconds (21.68 k allocations: 18.347 MiB)


  8.178048 seconds (21.75 k allocations: 18.350 MiB)


  5.595569 seconds (21.59 k allocations: 18.344 MiB)


  1.531474 seconds (21.33 k allocations: 5.510 MiB)


  1.200856 seconds (21.32 k allocations: 5.509 MiB)


  0.820974 seconds (21.28 k allocations: 2.942 MiB)


  0.364006 seconds (21.25 k allocations: 1.980 MiB)


  0.093983 seconds (21.21 k allocations: 1.138 MiB)
  0.010648 seconds (21.08 k allocations: 860.289 KiB)
  0.002969 seconds (21.05 k allocations: 819.789 KiB)
Total simulation time: 52.065341391

 seconds
pl = plot(xlabel="θ_x Angle", ylabel="⟨Z62⟩")
plot!(tn_angles, tn_vals, label="Flatiron Tensor Network", color="black", linewidth=3)
plot!(google_angles, google_vals, label="Google Lightcone Statevector", color="grey", linewidth=3)
plot!(x_angles, expectations, marker=:circle, label="Pauli Propagation", linewidth=3, color="royalblue")
scatter!(IBM_angles, IBM_mitigated_vals, label="IBM Quantum Computer", color="Red", ms=4)

svg