Bug in mean photon number

[trevor-vincent]

Issue description

Here is a minimal example that reproduces the bug:

import strawberryfields as sf
from strawberryfields import ops

modes = 4
fock_cutoff = 4
prog = sf.Program(modes)

with prog.context as q:
    ops.Sgate(.54) | q[0]
    ops.Sgate(.593) | q[1]
    ops.Sgate(.23) | q[2]
    ops.Sgate(.89) | q[3]

eng = sf.Engine("tf", backend_options={"cutoff_dim": fock_cutoff})
result = eng.run(prog)    
state = result.state
for i in range(0,modes):
    print("mean photon of mode " + str(i) + " = " + str((state.mean_photon(i)[0].numpy())))

This produces the output:

mean photon of mode 0 = 0.17980294
mean photon of mode 1 = 0.20575231
mean photon of mode 2 = 0.20575231
mean photon of mode 3 = 0.20575231

The mean photon of mode 2 and 3 should be different from mode 1 because the squeezing parameters are different.

• Expected behavior: (What you expect to happen)

mean photon of mode 0 = 0.17980294
mean photon of mode 1 = 0.20575231
mean photon of mode 2 = some different number
mean photon of mode 3 = some different number

• Actual behavior: (What actually happens)

mean photon of mode 0 = 0.17980294
mean photon of mode 1 = 0.20575231
mean photon of mode 2 = 0.20575231
mean photon of mode 3 = 0.20575231

• Reproduces how often: (What percentage of the time does it reproduce?)

100%

• System information: (post the output of import strawberryfields as sf; sf.about())

Strawberry Fields: a Python library for continuous-variable quantum circuits.
Copyright 2018-2020 Xanadu Quantum Technologies Inc.

Python version:            3.8.2
Platform info:             Linux-5.6.8-arch1-1-x86_64-with-glibc2.2.5
Installation path:         /home/tvincent/.local/lib/python3.8/site-packages/strawberryfields
Strawberry Fields version: 0.15.1
Numpy version:             1.18.4
Scipy version:             1.4.1
SymPy version:             1.6
NetworkX version:          2.4
The Walrus version:        0.13.0
Blackbird version:         0.2.4
TensorFlow version:        2.3.1

[co9olguy]

Thanks for discovering and reporting this @trevor-vincent!

I've determined that the source of this bug is the function tfbackend/ops.py:reduced_density_matrix (used internally when computing mean_photon). The final for loop in this function fails to account for the fact that removing subsystems in reduced_state will cause the mode indices m to be "off by one" for each subsystem that has been previously removed during the loop.

This is a bug at a core part of the codebase, so is high priority to fix. Thankfully, the fix looks relatively simple, something along the lines of:

for m in range(num_modes):
    if m != mode:
        reduced_state = partial_trace(reduced_state, m - removed_cnt, False, batched)
        removed_cnt += 1