pure states are not pure

[ziofil]

I tried to boil it down as much as possible. For some reason, states that should be pure are identified as mixed (in both the tf and fock backends):

import strawberryfields as sf

circuit = sf.Program(2)
eng  = sf.Engine("fock", backend_options={'cutoff_dim': 3})

with circuit.context as q:
    sf.ops.Fock(1) | q[0]
    sf.ops.Fock(1) | q[1]
    
output = eng.run(circuit)
output.state.is_pure

# False

I haven't investigated the issue yet

[co9olguy]

Hi @ziofil, I think what is happening here is that you have two local state preparations. IIRC the way things rae coded is that local state preparations force the simulator to use a mixed state representation.

When you are preparing a state locally, you have two options for how you might implement: either do something like
|new_state><current_state|,
or you trace out the current subsystem and replace it with the new state (which will necessarily make it mixed).

Obviously the first solution only works if you have a factorized state (and furthermore know the local pure state), so it was decided to do the second solution in general.

There is the option we could consider to---in this case at least---check whether you are still in the vacuum. If so, then the update would be much easier to implement

[ziofil]

Oh I see, thank you.
I agree with your last point, maybe it would be a good idea to allow for a pure separable state preparation (even I managed to get confused!)

[josh146]

👍, this is something that trips me up somewhat often (I always forget that Coherent, Fock etc result in a mixed state)

[glassnotes]

@ziofil since there is not a bug I'm going to close this, but have transferred the details into a potential new feature / enhancement (#490). Please feel free to add any additional details and ideas over there 👍

[ziofil]

If anyone stumbles here, this is a possible workaround if your initial state is simple enough to be written in this way:

import strawberryfields as sf

modes = 2
cutoff = 3
circuit = sf.Program(modes)
eng  = sf.Engine("fock", backend_options={'cutoff_dim': cutoff})

ket = np.zeros([cutoff]*modes, dtype=np.complex128)
ket[1,1] = 1.0 + 0.0j

with circuit.context as q:
    sf.ops.Ket(ket) | q
    
output = eng.run(circuit)
output.state.is_pure

# True