Pauli iter all #598

fdmalone · 2023-07-21T17:53:54Z

Adds stim.PauliString.iter_all method.

Following the discussion in #397 the algorithm is a combination of finding the next lexicographically ordered permutation of w bits followed by iterating over all 3^w PauliStrings given this permutation of the qubit labels. For the first part I modified the bit twiddle algorithm to account for multiple words, for the second part I just iterate over 3^w integers and map this to a Pauli using the ternary representation of the integer. This was a bit trickier to get right than I expected, plenty of edge cases.

The modifications for the bit twiddle algorithm are quite specific and I was a bit torn between adding more general operators to simd_bits (like left / right shift + subtraction) which would be cleaner, and what I did, which is a bit clunky.

There are also a few optimizations that could be made which I haven't:

for small w I could avoid the repeated loop over 3^w and store the Pauli strings in the first pass, rather than repeating this work for each of the num_qubits C w permutations (combinations))
some bitwise algorithms are suboptimal (counting trailing zeros) and/or may already exist in stim but I couldn't find them. Happy to change these.
The decoded bit labels / locations could be stored rather than decode each time.

Draft for the moment until:

Fix pybind issue.
Profile against naive python implmentation.
Check docstrings / algorithm description

Strilanc

Looking good so far. A few ideas commented.

doc/python_api_reference_vDev.md

src/stim/py/numpy.pybind.h

Strilanc · 2023-07-21T22:30:35Z

src/stim/stabilizers/pauli_string.pybind.cc

+        pybind11::arg("num_qubits"),
+        pybind11::kw_only(),
+        pybind11::arg("min_weight") = pybind11::none(),
+        pybind11::arg("max_weight") = pybind11::none(),


A useful argument here would be allowed_paulis: str, so the user could restrict to X errors (allowed_paulis="X") or not-Y-errors (allowed_paulis="XZ").

Good point.

I still need to add this, and also the random selection of signs/phases.

src/stim/stabilizers/pauli_string.pybind.cc

src/stim/stabilizers/pauli_string_iter.inl

src/stim/stabilizers/pauli_string_iter.perf.cc

Strilanc · 2023-07-21T22:55:13Z

doc/python_api_reference_vDev.md

+) -> None:
+    """Seed the iterator with a given qubit pattern.


Not clear to me what this means.

I'm actually going to remove this function, after closer thought it's not super helpful. The idea was to be able to "seed" the iterator at a specific bit pattern which may be difficult to reach if w was large, but this can be tested on the C++ side more easily, and not sure if it's actually useful on the python side.

Actually scratch that, I needed it for testing random starting points for long strings with higher weight. I will change the name to something more descriptive

I renamed this to set_current_permutation.

fdmalone · 2023-07-21T23:01:24Z

Sorry, I forced pushed when it was still a draft to clean up the many doc related fixes. Should be ready for review now sans profiling.

Strilanc · 2023-07-21T23:17:02Z

I was a bit torn between adding more general operators to simd_bits (like left / right shift + subtraction) which would be cleaner, and what I did, which is a bit clunky

I'd be fine with those additions to simd_bits.

Strilanc · 2023-07-22T00:49:14Z

Here's bit twiddle python code that generates all the requested pauli strings of a given weight, assuming length<=32.

from typing import Iterable


def bits_to_paulis(x: int, n: int) -> str:
    bits = bin(x)[2:].rjust(n*2, '0')[-n*2:]
    return ''.join('_XYZ'[int(bits[k:k+2], 2)] for k in range(0, len(bits), 2))


def count_trailing_zeros(x: int) -> int:
    t = 0
    while not (x & 1):
        x >>= 1
        t += 1
    return t


def masked_increment(x: int, mask: int) -> int:
    return ((x | ~mask) + 1) & mask


def next_bitstring_of_same_hamming_weight(x: int) -> int:
    c1 = x | (x - 1)
    c2 = c1 + 1
    c3 = (~c1 & -~c1) - 1
    c4 = c3 >> (count_trailing_zeros(x) + 1)
    return c2 | c4


def iter_pauli_strings(weight: int, length: int) -> Iterable[str]:
    hamming_mask = (1 << weight) - 1
    while hamming_mask < 2**length:
        # Spread out bits into pairs.
        h = hamming_mask
        mh = 0b0000000000000000000000000000000011111111111111110000000000000000
        ml = 0b0000000000000000000000000000000000000000000000001111111111111111
        h = (h & ml) | ((h & mh) << 16)
        mh = 0b0000000000000000111111110000000000000000000000001111111100000000
        ml = 0b0000000000000000000000001111111100000000000000000000000011111111
        h = (h & ml) | ((h & mh) << 8)
        mh = 0b0000000011110000000000001111000000000000111100000000000011110000
        ml = 0b0000000000001111000000000000111100000000000011110000000000001111
        h = (h & ml) | ((h & mh) << 4)
        mh = 0b0000110000001100000011000000110000001100000011000000110000001100
        ml = 0b0000001100000011000000110000001100000011000000110000001100000011
        h = (h & ml) | ((h & mh) << 2)
        mh = 0b0010001000100010001000100010001000100010001000100010001000100010
        ml = 0b0001000100010001000100010001000100010001000100010001000100010001
        h = (h & ml) | ((h & mh) << 1)
        h |= h << 1

        # Iterate over non-00 values of masked bit pairs, with 00 elsewhere
        xz_mask = 0
        for _ in range(3**weight):
            xz_mask = masked_increment(xz_mask, h)
            xz_mask |= ~(xz_mask | (xz_mask >> 1)) & 0b0101010101010101010101010101010101010101010101010101010101010101
            xz_mask &= h
            yield bits_to_paulis(xz_mask, length)

        # Next mask
        hamming_mask = next_bitstring_of_same_hamming_weight(hamming_mask)


t = 0
for e in iter_pauli_strings(weight=3, length=20):
    print(e)
    t += 1
print(t)

Strilanc · 2023-07-22T02:21:05Z

A better one (in particular see pair_sat_increment):

import math
from typing import Iterable, Tuple


def bits_to_paulis(x: int, z: int, n: int) -> str:
    s = ''
    for k in range(n):
        s += '_XYZ'[(x & 1) + (z & 1) * 2]
        x >>= 1
        z >>= 1
    return s[::-1]


def count_trailing_zeros(x: int) -> int:
    t = 0
    while not (x & 1):
        x >>= 1
        t += 1
    return t


def masked_increment(x: int, mask: int) -> int:
    return ((x | ~mask) + 1) & mask


def next_bitstring_of_same_hamming_weight(x: int) -> int:
    c1 = x | (x - 1)
    c3 = ((c1 + 1) & ~c1) - 1
    c4 = c3 >> (count_trailing_zeros(x) + 1)
    return (c1 + 1) | c4


def pair_sat_increment(x: int, z: int, m: int) -> Tuple[int, int]:
    """Finds the next (x, z) such that x | z == m."""
    inc = x & z
    up = ~inc
    inc |= ~m
    inc += 1
    inc &= m
    up &= inc
    z &= inc | ~x
    z ^= x & up
    x ^= up
    return x, z


def iter_pauli_strings(weight: int, length: int) -> Iterable[str]:
    h = (1 << weight) - 1
    while h < 2**length:
        x, z = h, 0
        for _ in range(3**weight):
            yield bits_to_paulis(x, z, length)
            x, z = pair_sat_increment(x, z, h)
        h = next_bitstring_of_same_hamming_weight(h)


t = 0
seen = set()
for e in iter_pauli_strings(weight=4, length=10):
    assert e not in seen, e
    seen.add(e)
    print(e)
    t += 1
print(t, 3**4 * math.factorial(10) // math.factorial(6) // math.factorial(4))

fdmalone · 2023-07-24T18:35:10Z

Ok, I can replace my ternary iteration with the bit twiddle above.

fdmalone · 2023-07-25T22:29:53Z

In particular, I will add additional ops to simd_bits (left/right shift and an adder) and replace my ternary iteration with the bit twiddle. I may separate out the new operations in a separate PR.

Strilanc · 2023-07-25T22:58:27Z

SGTM

From #598 added +=, >>= and <<= to simd_bits. It wasn't obvious to me that these could use word level parallelism without using more memory? For example, the shifts could store the relevant carry masks and or these at the end but this would require a temporary of the same size as the simd_bits instance.

Rename some functions. WIP. Better tests. Update python test.

fdmalone · 2023-09-05T17:27:51Z

The windows failures seem to be due to a memory error

          pauli_it = list(
              stim.PauliString.iter_all(
  >               num_qubits, min_weight=min_weight, max_weight=max_weight
              )
          )
  E       MemoryError

The test_iter_all_random_permutation tests also seems to trigger test failures on windows periodically on win32 platforms. Not sure what's going on there.

Caught this when trying to address #598.

Strilanc

I think the major tasks left are debugging the windows crash and adding pytest unit tests using the python side of the API.

Strilanc · 2023-09-16T21:06:20Z

src/stim/stabilizers/pauli_string_iter.inl

+        // x ^= up
+        result.xs ^= up;
+        cur_k++;
+        return true;


Why is this a while loop if it always exits?

Strilanc · 2023-09-16T21:07:01Z

src/stim/stabilizers/pauli_string_iter.inl

+
+template <size_t W>
+bool PauliStringIterator<W>::pair_sat_increment() {
+    // This will overflow for large cur_w.


Overflow is bad or good here?

Bad. I think I should assert on max_weight <= 40 on the pybind side? Or catch the wrapping and exit? Or just use the too large value since it would be pretty hard to iterate through 10^19 values.

Strilanc · 2023-09-16T21:11:57Z

src/stim/stabilizers/pauli_string_iter.inl

+template <size_t W>
+bool PauliStringIterator<W>::pair_sat_increment() {
+    // This will overflow for large cur_w.
+    size_t num_terms = static_cast<size_t>(pow(3, cur_w));


Are you sure this actually returns the right answer for all relevant values? Doubles have 53 bits of precision; not enough for a 64 bit integer. It'd be safer to have a method like

uint64_t pow3(uint64_t p) { assert(p < 41); uint64_t r = 1; if (p & 1) r *= 3 if (p & 2) r *= 9; if (p & 4) r *= 81; if (p & 8) r *= 6561; if (p & 16) r *= 43046721; if (p & 32) r *= 1853020188851841; return r; }

Strilanc · 2023-09-16T21:14:17Z

src/stim/stabilizers/pauli_string_iter.inl

+        simd_bits<W> one(cur_perm.num_bits_padded());
+        one.u64[0] = uint64_t{1};
+        // inc += 1
+        inc += one;


Maybe add a ++x method to simd_bits?

Yeah, it probably makes sense to add -- / -= too. For example, there are a few awkward parts where I'm doing (1 << sum_number_greather_than_64) - 1 in a clunky way.

fdmalone · 2023-10-26T14:35:05Z

I'm going to close this as I won't have time to come back to it for another month or so.

Move numpy array size utilites to numpy.pybind.

c4e91c4

fdmalone force-pushed the pauli_iter_all branch from 3c29370 to eb46021 Compare July 21, 2023 20:15

Add PauliString iter_all.

c8ef141

fdmalone force-pushed the pauli_iter_all branch from eb46021 to c8ef141 Compare July 21, 2023 22:54

Strilanc reviewed Jul 21, 2023

View reviewed changes

fdmalone marked this pull request as ready for review July 21, 2023 23:00

fdmalone added 2 commits July 24, 2023 17:32

Address review comments.

6ddc2ad

Rename seed_iterator, make internal and regenerate docs.

1c5b6dc

Fix infinite loop caused by overflow of num_terms.

2f1a9b9

fdmalone mentioned this pull request Aug 10, 2023

Add left/right shift and addition to simd_bits. #603

Merged

fdmalone and others added 6 commits September 4, 2023 21:28

Merge branch 'main' into pauli_iter_all

b89611b

Fix bug with simd_bits +=.

2ec29e2

Incorporate pair_sat_increment from review comment.

801356c

Rename some functions. WIP. Better tests. Update python test.

Fix doctest.

6cbcdfe

Update api reference.

19c3632

More doc issues.

0d762a1

fdmalone mentioned this pull request Sep 5, 2023

Fix bug with simd_bits +=. #633

Merged

Strilanc pushed a commit that referenced this pull request Sep 11, 2023

Fix bug with simd_bits +=. (#633)

0fdddef

Caught this when trying to address #598.

Strilanc requested changes Sep 16, 2023

View reviewed changes

fdmalone closed this Oct 26, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Pauli iter all #598

Pauli iter all #598

fdmalone commented Jul 21, 2023 •

edited

Loading

Strilanc left a comment

Strilanc Jul 21, 2023

fdmalone Jul 24, 2023

fdmalone Sep 5, 2023

Strilanc Jul 21, 2023

fdmalone Jul 24, 2023

fdmalone Jul 24, 2023

fdmalone Sep 5, 2023

fdmalone commented Jul 21, 2023

Strilanc commented Jul 21, 2023

Strilanc commented Jul 22, 2023 •

edited

Loading

Strilanc commented Jul 22, 2023

fdmalone commented Jul 24, 2023

fdmalone commented Jul 25, 2023 •

edited

Loading

Strilanc commented Jul 25, 2023

fdmalone commented Sep 5, 2023

Strilanc left a comment

Strilanc Sep 16, 2023

Strilanc Sep 16, 2023

fdmalone Sep 21, 2023 •

edited

Loading

Strilanc Sep 16, 2023

Strilanc Sep 16, 2023

fdmalone Sep 21, 2023 •

edited

Loading

fdmalone commented Oct 26, 2023

Pauli iter all #598

Pauli iter all #598

Conversation

fdmalone commented Jul 21, 2023 • edited Loading

Strilanc left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fdmalone commented Jul 21, 2023

Strilanc commented Jul 21, 2023

Strilanc commented Jul 22, 2023 • edited Loading

Strilanc commented Jul 22, 2023

fdmalone commented Jul 24, 2023

fdmalone commented Jul 25, 2023 • edited Loading

Strilanc commented Jul 25, 2023

fdmalone commented Sep 5, 2023

Strilanc left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fdmalone Sep 21, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fdmalone Sep 21, 2023 • edited Loading

Choose a reason for hiding this comment

fdmalone commented Oct 26, 2023

fdmalone commented Jul 21, 2023 •

edited

Loading

Strilanc commented Jul 22, 2023 •

edited

Loading

fdmalone commented Jul 25, 2023 •

edited

Loading

fdmalone Sep 21, 2023 •

edited

Loading

fdmalone Sep 21, 2023 •

edited

Loading