Merge pull request #262 from nonhermitian/integrate-bitstrings

Integrate new balanced calibrations into M3

docs/advanced.rst

@@ -41,13 +41,15 @@ When calibrating the mitigator, it is possible to vary the number of shots per c
 ``method``
 ~~~~~~~~~~
 
-There are three ways to do the calibration.  The default ``balanced`` method executes :math:`2N` circuits
+There are three ways to do the calibration.  The default ``balanced`` method a set of circuits
 with varying bitstring patterns such that the :math:`|0\rangle` and :math:`|1\rangle` states are each
-prepared :math:`N` times and averaged over.  For example, the balanced bit-strings over four qubits are
+prepared an even number of times, as is their pair-wise correlations.  For example, the 
+balanced bit-strings over four qubits are:
 
 .. jupyter-execute::
 
-    mthree.circuits.balanced_cal_strings(4)
+    gen = mthree.generators.HadamardGenerator(4)
+    list(gen)
 
 The ``independent`` method also sends :math:`2N` circuits but measures only a single qubit at a time.
 As such, this is a truely uncorrelated calibration process.

docs/balanced.rst

@@ -4,41 +4,41 @@
 Balanced calibrations
 #####################
 
-The default calibration method for M3 is what is called "balanced" calibration.  We came up
-with this method as an intermediary step between truly "independent" calibrations that run
-two circuits for each qubit to get the error rates for :math:`|0\rangle` :math:`|1\rangle`,
-and "marginal" calibrations that execute only two circuits :math:`|0\rangle^{\otimes N}`
-and :math:`|1\rangle^{\otimes N}`.  These two methods are expensive, or can lead to inaccurate
-results when state prep errors are present, respectively.
+The default calibration method for M3 is what is called "balanced" calibration.  Balanced calibrations
+sample all independent and pair-wise readout error rates evenly; hence the name "balanced".  In M3 v3
+and higher, this routine has been updated to use a method from Bravyi et al, Phys. Rev. A 103, 042605 (2021).
 
-Balanced calibrations run :math:`2N` circuits for :math:`N` measured qubits, but the calibration
-circuits are chosen in such a way as to sample each error rate :math:`N` times.  For example,
-consider the balanced calibration circuits for 5 qubits:
+To see the bit-patterns used to generate the calibration circuits one can call the generator explicitly.  E.g. to
+see the strings for 5-qubits one would do 
 
 .. jupyter-execute::
 
-    from qiskit_ibm_runtime.fake_provider import FakeAthensV2
     import mthree
 
-    mthree.circuits.balanced_cal_strings(5)
+    gen = mthree.generators.HadamardGenerator(5)
+    list(gen)
 
 
-For every position in the bit-string you will see that a `0` or `1` appears `N` times.
+For every position in the bit-string you will see that `0` or `1` appear an equal number of times,
+and the same is true for all pair-wise combinations of `0` and `1`.
 If there is a `0`, then that circuit samples the :math:`|0\rangle` state for that qubit,
-similarly for the `1` element.  So when we execute the `2N` balanced calibration circuits
-using `shots` number of samples, then each error rate in the calibration data is actually
-being sampled `N*shots` times.  Thus, when you pass the `shots` value to M3, in the balanced
-calibration mode internally it divides by the number of measured qubits so that the precision
-matches the precison of the other methods.  That is to say that the following:
+similarly for the `1` element.  So when we execute the balanced calibration circuits
+using `shots` number of samples, each error rate in the calibration data is actually
+being sampled more times than requested.  Thus, when you pass the `shots` value to M3, in the balanced
+calibration mode internally it divides by the number of measured qubits so that the precision of each
+error rate matches the precison of the other methods.  That is to say that the following:
 
  .. jupyter-execute::
 
+    from qiskit_ibm_runtime.fake_provider import FakeAthensV2
+
     backend = FakeAthensV2()
     mit = mthree.M3Mitigation(backend)
     mit.cals_from_system(method='balanced')
 
-Will sample each qubit error rate `10000` times regardless of which method is used.  Moreover,
-this also yields a calibration process whose overhead is independent of the number of qubits
-used.  Note that, when using a simulator or "fake" device, M3 defaults to `independent`
-calibration mode for efficiency.  As such, to enable `balanced` calibration on a simulator
-one must explicitly set the `method`` as done above.
+Will sample each indepdendent qubit error rate `10000` times (or the max allowed by the target system if less) 
+regardless of which method is used. All pair-wise correlations are measured half this number of times. This
+yields a calibration process whose  overhead is independent of the number of qubits used; there is no additional
+cost to compute the calibration over a full device.  Note that, when using a simulator or "fake" device, 
+M3 defaults to `independent` calibration mode for efficiency.  As such, to enable `balanced` calibration on a 
+simulator one must explicitly set the `method`` as done above.

mthree/circuits.py

@@ -48,33 +48,11 @@ def _marg_meas_states(qubits, num_system_qubits, initial_reset=False):
     return [qc0, qc1]
 
 
-def balanced_cal_strings(num_qubits):
-    """Compute the 2*num_qubits strings for balanced calibration.
-
-    Parameters:
-        num_qubits (int): Number of qubits to be measured.
-
-    Returns:
-        list: List of strings for balanced calibration circuits.
-    """
-    strings = []
-    for rep in range(1, num_qubits + 1):
-        str1 = ""
-        str2 = ""
-        for jj in range(int(np.ceil(num_qubits / rep))):
-            str1 += str(jj % 2) * rep
-            str2 += str((jj + 1) % 2) * rep
-
-        strings.append(str1[:num_qubits])
-        strings.append(str2[:num_qubits])
-    return strings
-
-
-def balanced_cal_circuits(cal_strings, layout, system_qubits, initial_reset=False):
+def balanced_cal_circuits(generator, layout, system_qubits, initial_reset=False):
     """Build balanced calibration circuits.
 
     Parameters:
-        cal_strings (list): List of strings for balanced cal circuits.
+        generator (HadamardGenerator): Generator of balanced cal bit-strings for circuits
         layout (list): Logical to physical qubit layout
         initial_reset (bool): Use resets at beginning of circuit.
         system_qubits (int): Number of qubits in system
@@ -83,16 +61,16 @@ def balanced_cal_circuits(cal_strings, layout, system_qubits, initial_reset=Fals
         list: List of balanced cal circuits.
     """
     circs = []
-    num_active_qubits = len(cal_strings[0])
-    for string in cal_strings:
+    num_active_qubits = generator.num_qubits
+    for bit_array in generator:
         qc = QuantumCircuit(system_qubits, num_active_qubits)
         if initial_reset:
             qc.barrier()
             qc.reset(range(system_qubits))
             qc.reset(range(system_qubits))
             qc.reset(range(system_qubits))
-        for idx, bit in enumerate(string[::-1]):
-            if bit == "1":
+        for idx, bit in enumerate(bit_array[::-1]):
+            if bit == 1:
                 qc.x(layout[idx])
         qc.measure(layout, range(num_active_qubits))
         circs.append(qc)

mthree/generators/__init__.py

@@ -0,0 +1,16 @@
+# This code is part of Mthree.
+#
+# (C) Copyright IBM 2024.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Generators
+"""
+
+from .hadamard import HadamardGenerator

mthree/generators/hadamard.pyx

@@ -0,0 +1,93 @@
+# This code is part of Mthree.
+#
+# (C) Copyright IBM 2023.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+# pylint: disable=no-name-in-module
+"""mthree Hadamard array generator"""
+cimport cython
+
+from libc.stdlib cimport malloc, free
+from libc.math cimport log2, floor
+
+import numpy as np
+cimport numpy as np
+
+from mthree.exceptions import M3Error
+
+
+cdef class HadamardGenerator():
+    """Hadamard calibration generator"""
+    cdef public str name
+    cdef unsigned int p
+    cdef unsigned char * integer_bits
+    cdef unsigned char * out_bits
+    cdef public unsigned int num_qubits
+    cdef public unsigned int length
+    cdef unsigned int _iter_index
+
+    @cython.boundscheck(False)
+    def __cinit__(self, unsigned int num_qubits):
+        """Hadamard calibration generator
+        
+        Generates a set of bit-arrays that evenly
+        sample all independent and pair-wise correlated
+        measurement errors.
+
+        References:
+            Bravyi et al, Phys. Rev. A 103, 042605 (2021)
+        """
+        self.name = 'hadamard'
+        self.num_qubits = num_qubits
+        self.p = <unsigned int>floor(log2(num_qubits)+1)
+        self.length = <unsigned int>(2**self.p)
+        self.integer_bits = <unsigned char *>malloc(self.p*sizeof(unsigned char))
+        # output set of bitstrings
+        self.out_bits = <unsigned char *>malloc(num_qubits*sizeof(unsigned char))
+        self._iter_index = 0
+
+    def __dealloc__(self):
+        if self.integer_bits is not NULL:
+            free(self.integer_bits)
+        if self.out_bits is not NULL:
+            free(self.out_bits)
+
+    def __iter__(self):
+        self._iter_index = 0
+        return self
+
+    def __next__(self):
+        if self._iter_index < self.length:
+            self._iter_index += 1
+            return self._generate_array(self._iter_index-1)
+        else:
+            raise StopIteration
+
+    @cython.boundscheck(False)
+    def _generate_array(self, unsigned int index):
+        if index > self.length-1:
+            raise M3Error('Index must within generator length {}'.format(self.length))
+        cdef size_t kk, jj
+        cdef unsigned int tot
+        cdef list out = []
+
+        # Set the bitstrings for the integer_bits
+        for kk in range(self.p):
+            self.integer_bits[self.p-kk-1] = (index >> kk) & 1
+
+        for kk in range(self.num_qubits):
+            tot = 0
+            for jj in range(self.p):
+                tot += self.integer_bits[self.p-jj-1] and ((kk+1) >> jj) & 1
+            self.out_bits[kk] = tot % 2
+
+        # Need to return copy since the underlying memory will be reused
+        # It turns out that it is faster to copy the NumPy array then it is
+        # to copy the underlying MemoryView
+        return np.asarray((<np.uint8_t[:self.num_qubits]> self.out_bits)).copy()

mthree/mitigation.py

@@ -25,10 +25,10 @@
 import runningman as rm
 from runningman.utils import is_ibm_backend
 
+from mthree.generators import HadamardGenerator
 from mthree.circuits import (
     _tensor_meas_states,
     _marg_meas_states,
-    balanced_cal_strings,
     balanced_cal_circuits,
 )
 from mthree.direct import direct_solver as direct_solve
@@ -320,10 +320,9 @@ def _grab_additional_cals(
             raise M3Error("System is not set.  Use 'cals_from_file'.")
         if self.single_qubit_cals is None:
             self.single_qubit_cals = [None] * self.num_qubits
-        if self.cal_shots is None:
-            if shots is None:
-                shots = min(self.system_info["max_shots"], 10000)
-            self.cal_shots = shots
+        if shots is None:
+            shots = min(self.system_info["max_shots"], 10000)
+        self.cal_shots = shots
         if self.rep_delay is None:
             self.rep_delay = rep_delay
 
@@ -361,7 +360,7 @@ def _grab_additional_cals(
             )
 
         num_cal_qubits = len(qubits)
-        cal_strings = []
+        generator = None
         # shots is needed here because balanced cals will use a value
         # different from cal_shots
         shots = self.cal_shots
@@ -371,14 +370,14 @@ def _grab_additional_cals(
                 qubits, self.num_qubits, initial_reset=initial_reset
             )
         elif method == "balanced":
-            cal_strings = balanced_cal_strings(num_cal_qubits)
+            generator = HadamardGenerator(num_cal_qubits)
             trans_qcs = balanced_cal_circuits(
-                cal_strings, qubits, self.num_qubits, initial_reset=initial_reset
+                generator, qubits, self.num_qubits, initial_reset=initial_reset
             )
-            shots = self.cal_shots // num_cal_qubits
-            if self.cal_shots / num_cal_qubits != shots:
+            shots = 2 * self.cal_shots // generator.length
+            if 2 * self.cal_shots / generator.length != shots:
                 shots += 1
-            self._balanced_shots = shots * num_cal_qubits
+            self._balanced_shots = shots * generator.length
         # Independent
         else:
             trans_qcs = []
@@ -421,12 +420,12 @@ def _grab_additional_cals(
         if async_cal:
             thread = threading.Thread(
                 target=_job_thread,
-                args=(jobs, self, qubits, num_cal_qubits, cal_strings),
+                args=(jobs, self, qubits, num_cal_qubits, generator),
             )
             self._thread = thread
             self._thread.start()
         else:
-            _job_thread(jobs, self, qubits, num_cal_qubits, cal_strings)
+            _job_thread(jobs, self, qubits, num_cal_qubits, generator)
 
         return jobs
 
@@ -733,15 +732,15 @@ def _thread_check(self):
             raise self._job_error  # pylint: disable=raising-bad-type
 
 
-def _job_thread(jobs, mit, qubits, num_cal_qubits, cal_strings):
+def _job_thread(jobs, mit, qubits, num_cal_qubits, generator):
     """Run the calibration job in a different thread and post-process
 
     Parameters:
         jobs (list): A list of job instances
         mit (M3Mitigator): The mitigator instance
         qubits (list): List of qubits used
         num_cal_qubits (int): Number of calibration qubits
-        cal_strings (list): List of cal strings for balanced cals
+        generator (None or HadamardGenerator): Generator for bit-arrays for balenced cals
     """
     counts = []
     for job in jobs:
@@ -818,20 +817,20 @@ def _job_thread(jobs, mit, qubits, num_cal_qubits, cal_strings):
     else:
         cals = [np.zeros((2, 2), dtype=np.float32) for kk in range(num_cal_qubits)]
 
-        for idx, count in enumerate(counts):
-            target = cal_strings[idx][::-1]
+        for idx, target in enumerate(generator):
+            count = counts[idx]
             good_prep = np.zeros(num_cal_qubits, dtype=np.float32)
             # divide by 2 since total shots is double
-            denom = mit._balanced_shots
-
+            denom = mit._balanced_shots / 2
+            target = target[::-1]
             for key, val in count.items():
                 key = key[::-1]
                 for kk in range(num_cal_qubits):
-                    if key[kk] == target[kk]:
+                    if int(key[kk]) == target[kk]:
                         good_prep[kk] += val
 
             for kk, cal in enumerate(cals):
-                if target[kk] == "0":
+                if target[kk] == 0:
                     cal[0, 0] += good_prep[kk] / denom
                 else:
                     cal[1, 1] += good_prep[kk] / denom