Frequency custom operator fails for some states

[stavros11]

Describe the bug
The frequency custom operator gives wrong frequencies for some state vectors / probability distributions, such as the ones in the example below. Thanks to @igres26 for finding this.

To Reproduce

from qibo import models, gates

c = models.Circuit(2)
c.add(gates.H(0))
c.add(gates.CNOT(0, 1))
c.add(gates.M(0, 1))
result = c(nshots=int(1e6))

print(result.state())
print(result.frequencies())

While the final state vector is correct the frequencies contain only the '00' bitstring while the '11' is also to be measured with equal probability. The problem does not appear for nshots < 1e5 where the old tf native approach is used.

A simpler code to reproduce:

import numpy as np
from qibo import K

probs = np.zeros(4)
probs[0], probs[3] = 0.5, 0.5
frequencies = K.sample_frequencies(probs, nshots=int(1e6))
print(frequencies)

[scarrazza]

Thanks for reporting this. We should include the relevant missing test that did not spot this issue on time.
This may require changes at the op and python level, right?

[stavros11]

I think the issue is related to Metropolis only. It fails in cases where the allowed bitstrings of non-zero probability have Hamming distance greater than 1. In the above example the allowed bitstrings are '00' and '11' so the Metropolis starts in '00' and gets stuck there forever because it cannot to '11' directly with single bitflips but cannot cross '01' or '10' either because these have zero probability.

An easy solution is to replace the mechanism that generates new bitstrings from the current:

// Generate random index to flip its bit
int flip_index = ((int) rand_r(&seed) % nqubits);
// Flip the corresponding bit
int current_value = ((int64) shot >> flip_index) % 2;
int64 new_shot = shot + ((int64)(1 - 2 * current_value)) * ((int64) 1 << flip_index);

which adds only powers of 2 (bitflips) to

int64 new_shot = (shot + ((int64) rand_r(&seed) % nstates)) % nstates;

which adds any random integer allowing to go from '00' to '11' directly. I tested this and it fixes the issue but I am not sure how it affects other distributions.

My main concern is that I don't know how to make sure that the algorithm will be robust and work properly for all cases of distributions that may come up here. We have tests for random distributions but these give some probability to every bitstring so they did not capture the issue with the "sparse" distributions found here. We can add tests for sparse distributions but we may still be missing a different case.

[scarrazza]

@stavros11 ok, could you please check with @igres26 if this approach works for the tests he is doing?
What happens for other channels with these changes?