Run failure on DensityMatrixSimulator

[glanzz]

Description of the issue
Running a circuit imported from QASM is failing when run on DensityMatrixSimulator
How to reproduce the issue

from cirq.contrib.qasm_import import circuit_from_qasm
from cirq import DensityMatrixSimulator, measure_each, X
import cirq as c


circuit = None
with open("mycircuit.qasm", "r") as file:
    circuit = circuit_from_qasm(file.read())


s = DensityMatrixSimulator()
REPEATS = 1100
circuit.append(measure_each(*circuit.all_qubits()))
result = s.run(circuit,repetitions=REPEATS)

mycircuit.qasm

OPENQASM 2.0;
include "qelib1.inc";
qreg q[2];
x q[1];
cz q[1],q[0];
h q[0];
sdg q[1];
cz q[1],q[0];
y q[1];
h q[0];
cy q[0],q[1];
x q[0];
s q[1];
z q[0];
u2(3.2476939604661648,2.3022570151458246) q[1];
cx q[1],q[0];
cx q[0],q[1];
cx q[1],q[0];
swap q[1],q[0];
u2(-5.443849668735618,-0.10610130687637165) q[1];
z q[0];
sdg q[1];
x q[0];
cy q[0],q[1];
h q[0];
y q[1];
cz q[1],q[0];
s q[1];
h q[0];
cz q[1],q[0];

Traceback (most recent call last):
  File "/Users/xxx/Projects/myproject/cirq-sample.py", line 13, in 
    result = s.run(circuit,repetitions=REPEATS)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/xxx/Projects/myproject/cirqenv/lib/python3.11/site-packages/cirq/work/sampler.py", line 63, in run
    return self.run_sweep(program, param_resolver, repetitions)[0]
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/xxx/Projects/myproject/cirqenv/lib/python3.11/site-packages/cirq/sim/simulator.py", line 72, in run_sweep
    return list(self.run_sweep_iter(program, params, repetitions))
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/xxx/Projects/myproject/cirqenv/lib/python3.11/site-packages/cirq/sim/simulator.py", line 103, in run_sweep_iter
    records = self._run(
              ^^^^^^^^^^
  File "/Users/xxx/Projects/myproject/cirqenv/lib/python3.11/site-packages/cirq/sim/simulator_base.py", line 256, in _run
    return step_result.sample_measurement_ops(
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/xxx/Projects/myproject/cirqenv/lib/python3.11/site-packages/cirq/sim/simulator.py", line 809, in sample_measurement_ops
    indexed_sample = self.sample(measured_qubits, repetitions, seed=seed)
                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/xxx/Projects/myproject/cirqenv/lib/python3.11/site-packages/cirq/sim/simulator_base.py", line 400, in sample
    return self._sim_state.sample(qubits, repetitions, seed)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/xxx/Projects/myproject/cirqenv/lib/python3.11/site-packages/cirq/sim/simulation_product_state.py", line 165, in sample
    column = v.sample(qs, repetitions, seed)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/xxx/Projects/myproject/cirqenv/lib/python3.11/site-packages/cirq/sim/simulation_state.py", line 147, in sample
    return self._state.sample(self.get_axes(qubits), repetitions, seed)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/xxx/Projects/myproject/cirqenv/lib/python3.11/site-packages/cirq/sim/density_matrix_simulation_state.py", line 222, in sample
    return sim.sample_density_matrix(
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/xxx/Projects/myproject/cirqenv/lib/python3.11/site-packages/cirq/sim/density_matrix_utils.py", line 86, in sample_density_matrix
    result = prng.choice(len(probs), size=repetitions, p=probs)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "numpy/random/mtrand.pyx", line 973, in numpy.random.mtrand.RandomState.choice
ValueError: probabilities are not non-negative

Cirq version
You can get the cirq version by printing cirq.__version__. From the command line:

Cirq Version: 1.4.1

[myrimur]

It seems that the resulting probabilities can become slightly less than zero due to imprecision in floating-point calculations (but this is just my guess). cirq/sim/density_matrix_utils.py, line 86, in sample_density_matrix. For example, I have encountered probs = [ 7.5000006e-01 -8.4293717e-08 -2.1073429e-08 2.5000003e-01]

[mhucka]

Discussed during bit.ly/cirq-cynq today (2024-10-16): although the original example came from a QASM circuit definition, the problem itself is not related to QASM.