This PR implements the following changes:

change the default measurement basis to the computational basis
change the default random circuit to layers of randU gates (no single-qubit gates)
the default random circuit is built by feeding random matrices generated at definition (not compilation) using a seed
added layered argument to randomcircuit to create a layered structure.
added runcircuit for layered circuits
added observer for circuits

Example 1: random circuits with seed

N = 10
depth = 10

ψ1 = runcircuit(randomcircuit(N, depth; seed = 1234))
ψ2 = runcircuit(randomcircuit(N, depth; seed = 1234))
@show fidelity(ψ1,ψ2)
# fidelity(ψ1, ψ2) = 1.0
# 1.0

Example 2: run layered circuit with observer. By default, the observer measures at every layer the bond dimension at every bond, as well as the maximum bond dimension at any bond (I know this is somewhat redundant, but I think it's better to have both). I also changed the observer to a dictionary, since there was no real reason for it to have its own type. At least, no obvious reason I could think of.

N = 10
depth = 10
circuit = randomcircuit(N, depth; layered = true, seed = 1234)

obs = Dict()
runcircuit(circuit; observer! = obs)
@show obs
# Dict{Any,Any} with 2 entries:
#   "χmax" => Any[2, 4, 8, 16, 32, 32, 32, 32, 32, 32]
#   "χ"    => Any[Any[2, 1, 2, 1, 2, 1, 2, 1, 2], Any[2, 4, 2, 4, 2, 4, 2, 4, 2], Any[2, 4, 8, 4, 8, 4, 8, 4, 2], ...]

For a given pre-defined function, we can pass an observable to the runcircuit function in the form of a NamedTuple (or array thereof) where f is the function, name the function name. In this example, we measure the norm of the MPS:

N = 10
depth = 10
circuit = randomcircuit(N, depth; layered = true, seed = 1234)

obs = Dict()
runcircuit(circuit; observer! = obs)

observable = (f = norm, name = "norm")
runcircuit(circuit; observer! = obs, observables = observable)
@show obs["norm"]
# 10-element Array{Any,1}:
#                1.0 - 1.3484329742722466e-17im
#  0.9999999999999998 - 8.634013806807024e-18im
#  0.9999999999999987 + 2.1350855145183682e-17im
#  0.9999999999999976 + 1.8858646465528207e-17im
# ...

One can also measure a function with a set of arguments, which are simply passed as an array args, where the order corresponds to the order of the function arguments. Here is an example where we define a function for measuring Pauli operators for a given site, and use it to measure X on site 2:

N = 10
depth = 10
circuit = randomcircuit(N, depth; layered = true, seed = 1234)

# measure a Pauli operator at some site
function measurePauli(ψ::MPS, site::Int, pauli::String)
  orthogonalize!(ψ, site)
  ϕ = ψ[site]
  obs_op = gate(pauli, firstsiteind(ψ, site))
  T = noprime(ϕ * obs_op)
  return real((dag(T) * ϕ)[])
end

# measure Pauli X at site 2
observables = [(f = norm, name = "norm"), (f = measurePauli, name = "pauliX", args = [2, "X"])]
runcircuit(circuit; observer! = obs, observables = observables)
@show obs["pauliX"]
#   "pauliX" => Any[0.0567889, 0.0601699, -0.0532313, -0.0621552,...]

To run a function on a set of qubits, we use the argument sites, which should be a range or an array. In this example we measure the pauli X operator at every odd site:

N = 10
depth = 10
circuit = randomcircuit(N, depth; layered = true, seed = 1234)

# measure Pauli X at each odd site
observables = [(f = norm, name = "norm"), (f = measurePauli, name = "pauliX", args = ["X"], sites = 1:2:N)]
runcircuit(circuit; observer! = obs, observables = observables)
@show obs["pauliX"]
# obs["pauliX"] = Any[Any[-0.29293573358703523, -0.28194365860837894, -0.3210432622035978, 
# -0.22128416400408593, 0.4248564211985235], Any[-0.29293573358703445, -0.2564886958676067, 
# 0.191272827386181, ...

Codecov Report

Merging #172 (1e2542e) into master (d0eb542) will decrease coverage by 1.76%. The diff coverage is 81.95%.

@@            Coverage Diff             @@
##           master     #172      +/-   ##
==========================================
- Coverage   87.24%   85.48%   -1.77%     
==========================================
  Files          14       18       +4     
  Lines        1647     1853     +206     
==========================================
+ Hits         1437     1584     +147     
- Misses        210      269      +59

Impacted Files	Coverage Δ
src/PastaQ.jl	`100.00% <ø> (ø)`
src/measurements.jl	`30.37% <30.37%> (ø)`
src/inputoutput.jl	`50.00% <50.00%> (ø)`
src/lpdo.jl	`78.43% <50.00%> (-8.67%)`	:arrow_down:
src/distances.jl	`86.66% <83.87%> (+3.33%)`	:arrow_up:
src/circuits/gates.jl	`64.49% <85.00%> (+3.47%)`	:arrow_up:
src/tomography/tomographyutils.jl	`85.41% <85.41%> (ø)`
src/circuits/runcircuit.jl	`93.60% <90.27%> (-2.90%)`	:arrow_down:
src/observer.jl	`94.73% <94.73%> (+35.47%)`	:arrow_up:
src/circuits/circuits.jl	`98.03% <98.03%> (+56.78%)`	:arrow_up:
... and 12 more

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update d0eb542...1e2542e. Read the comment docs.

These are all nice additions to have. Good idea to just pass a seed as an argument to randomcircuit. I've got a few comments:

Can we just use layered by default? I don't see much of a disadvantage, most users won't even really notice (besides how it gets printed, so we should make sure the printing is nice).
When implementing the nonlayered version of randomcircuit, we could just make it layered from the beginning, and then at the end check if !layered and use vcat(gates...). This would make the code logic much simpler. We can also point this out to people as a simple way to convert to a nonlayered gate list.
In the layered version of runcircuit you should make sure to set move_sites_back = false in the intermediate layers to avoid extraneous swaps.
Setting Random.seed! in randomcircuit affects other random number generation as well (it changes the global Julia seed). To make a seed that is local to randomcircuit what you could do is if !isnothing(seed), create a local RNG rng = MersenneTwister(seed), and then pass rng to all calls of rand and randn inside randomcircuit. If isnothing(seed), then you can use the default RNG from Julia, i.e. rng = Random.GLOBAL_RNG.
For the observer functionality, I think the fact that the current design requires passing two objects (observer! and observables) indicates that we should just make this into a single object. Here would be my proposal:
```
# We can use a different name besides RuncircuitObserver
struct RuncircuitObserver <: AbstractObserver
functions::Dict{String, Function}
results::Dict{String, Any}
end
```

RuncircuitObserver(functions::Dict{String, Function}) = RuncircuitObserver(functions, Dict{String, Any}())

RuncircuitObserver(functions::Vector{Pair{String, Function}) = RuncircuitObserver(Dict(functions))

function measure_pauli(ψ::MPS, site::Int, pauli::String)

Note that it wouldn't be good to call `orthogonalize!` directly

on ψ since that would modify the MPS in runcircuit

Copying an MPS is cheap since it doesn't copy all of the ITensor data

ψ = orthogonalize!(copy(ψ), site) ϕ = ψ[site] obs_op = gate(pauli, firstsiteind(ψ, site)) T = noprime(ϕ obs_op) return real((dag(T) ϕ)[]) end

pauliX2(ψ::MPS) = measure_pauli(ψ, 2, "X") pauliYs(ψ::MPS) = [measure_pauli(ψ, n, "Y") for n in 1:length(ψ)]

obs = RuncircuitObserver(["χs" => linkdims, "χmax" => maxlinkdim, "norm" => norm, "pauliX2" => pauliX2, "pauliYs" => pauliYs]) runcircuit(circuit; observer! = obs)

obs.results stores the results of the measurements, for example:

obs.results["χs"] # All of the link dimensions after each layer obs.results["norm"] # The norms of each layer obs.results["pauliYs"] # A vector of results for measuring "Y" on every site for each layer

Yes sounds good. But what do you mean with nice printing?
OK
I call runcircuit to simulate each layer, so that should already be taken into account right?
OK
Looks good. Can we just call that CircuitObserver?

1. Yes sounds good. But what do you mean with nice printing?
Just that it doesn't look bad when someone prints the circuit.
2. OK

3. I call `runcircuit` to simulate each layer, so that should already be taken into account right?
Not if you call it multiple times, since at the end of each call, if you don't set move_sites_back = false, then it might do a bunch of swaps to return the sites to their original locations. So in the loop:
for layer in gates
M = runcircuit(M, layer; kwargs...)
measure!(observer!, M, obs)
end
after every layer, it will put all of the sites back unless you add move_sites_back = false. You'll then need to make sure the last call uses move_sites_back = true/false depending on what the user chose. It looks like runcircuit doesn't currently accept the options move_sites_back so we'll have to add that.
4. OK

5. Looks good. Can we just call that `CircuitObserver`?
Sure.

# add gate-observables with their nametag as a dictionary.
# here we measure X on each qubit
CircuitObserver("Xn" => "X").results
# "Xn" => Any[]

# measure also Y on qubit 1
CircuitObserver(["Xn" => "X", "Y1" => ("Y",1)]).results
#  "Xn" => Any[]
#  "Y1" => Any[]

# automatic name assignment
CircuitObserver(["X", ("Y",1), ("Z",1:10)]).results
#  "Z1:10" => Any[]
#  "X(n)"  => Any[]
#  "Y1"    => Any[]

# same for functions
CircuitObserver("χ" => maxlinkdim).results
#  "χ" => Any[]

# functions without name
CircuitObserver([norm, maxlinkdim]).results
#  "maxlinkdim" => Any[]
#  "norm"       => Any[]

# also a mix of functions and observables,
# as long as they are input as two separate arguments
CircuitObserver([("X",1)],[norm, maxlinkdim]).results
#  "maxlinkdim" => Any[]
#  "norm"       => Any[]
#  "X1"       => Any[]

Full example:

circuit = randomcircuit(4,4)
obs = CircuitObserver([("Z","Z"),("X")],[norm, maxlinkdim])
runcircuit(circuit; observer! = obs)
# --------------
#  WARNING
#
# A custom function is being measured during the gate evolution,
# but the MPS sites at depth D are now being restored to their 
# location at depth D-1. If any custom measurement is dependent 
# on the specific site ordering of the MPS, it should take the 
# re-ordering into account.
# --------------

# max bond dimension at each depth
measurements(obs,"maxlinkdim")
#4-element Array{Any,1}:
# 2
# 4
# 4
# 4

# X on each qubit at increasing depths
measurements(obs,"X(n)")
#4-element Array{Any,1}:
# Any[-0.543057283464483, 0.9037207770535344, 0.532220011216392, -0.24511903671723423]
# Any[-0.5430572834644839, -0.20053884847886705, 0.22651395521938691, -0.24511903671723506]
# Any[-0.4723532255342678, 0.1561121098988508, 0.2327918571328619, 0.061708358367673355]
# Any[-0.47235322553426795, -0.017036919895715104, -0.15400504719470903, 0.06170835836767424]

# ZZ correlation functions at the end of the circuit
measurements(obs,"ZZ")[end]
#4×4 Array{Float64,2}:
# 1.0        0.354003   0.141141   0.125133
# 0.354003   1.0        0.169881  -0.152048
# 0.141141   0.169881   1.0       -0.12335
# 0.125133  -0.152048  -0.12335    1.0

Great! What an epic PR.

GTorlai / PastaQ.jl

Observer for layered circuits #172

Codecov Report

Note that it wouldn't be good to call `orthogonalize!` directly

on ψ since that would modify the MPS in runcircuit

Copying an MPS is cheap since it doesn't copy all of the ITensor data

obs.results stores the results of the measurements, for example:

GTorlai / PastaQ.jl

Observer for layered circuits #172

Codecov Report

Note that it wouldn't be good to call orthogonalize! directly

on ψ since that would modify the MPS in runcircuit

Copying an MPS is cheap since it doesn't copy all of the ITensor data

obs.results stores the results of the measurements, for example:

Note that it wouldn't be good to call `orthogonalize!` directly