QSimSimulator.simulate_expectation_values for identity operators is wrong when the coefficient is other than one

[Hosseinberg]

The QSimSimulator.simulate_expectation_values simply just returns one for the identity operators, which is correct as long as there is no coefficient. When the observable has a coefficient the expectation value should be equal to the coefficient and not one. The following example calculates the expectation value of an observable in three different ways.

• Using the qsimcirq simulator
• Using cirq simulator
• Calculating the expectation value from the Hamiltonian

import cirq
import qsimcirq

pauli = 'IIIIIII'
num_qubits = 7
moments = 9
density = 0.95
qubits = cirq.LineQubit.range(num_qubits)
cirq_circuit = cirq.testing.random_circuit(qubits=qubits,
                                      n_moments=moments,
                                      op_density=density)

coeff = 0.25
hamiltonian = float(coeff) * cirq.PauliString(
                cirq.I(cirq.LineQubit(i)) for i, p in enumerate(pauli))

qubits = cirq.LineQubit.range(num_qubits)
energy_ham = hamiltonian.expectation_from_state_vector(
    state_vector= cirq.final_state_vector(cirq_circuit), 
    qubit_map={q: i for i, q in enumerate(qubits)},
    atol=0.01)

qsimSim = qsimcirq.QSimSimulator()
qsimcirq_result = qsimSim.simulate_expectation_values(cirq_circuit, hamiltonian) 

cirqSim = cirq.Simulator()
cirq_result = cirqSim.simulate_expectation_values(cirq_circuit, hamiltonian)  

#print(cirq_circuit)
print(f'energy_ham = {energy_ham}')
print(f'qsimcirq_result = {qsimcirq_result}')
print(f'cirq_result = {cirq_result}')

which returns:

energy_ham = (0.24999994039535522+0j)
qsimcirq_result = [(1+0j)]
cirq_result = [(0.24999994039535522+0j)]

As can be seen qsimcirq returns 1 regardless of the coefficient which I think is related to #565 and #566. I'm using 0.14.1.dev20220804.