Update computation of probability

qiskit_experiments/curve_analysis/standard_analysis/error_amplification_analysis.py

@@ -17,7 +17,6 @@
 import numpy as np
 
 import qiskit_experiments.curve_analysis as curve
-from qiskit_experiments.exceptions import CalibrationError
 
 
 class ErrorAmplificationAnalysis(curve.CurveAnalysis):
@@ -32,7 +31,29 @@ class ErrorAmplificationAnalysis(curve.CurveAnalysis):
 
         .. math::
             y = \frac{{\rm amp}}{2}\cos\left(x[{\rm d}\theta + {\rm apg} ] \
-            +{\rm phase\_offset}\right)+{\rm base}
+            -{\rm phase\_offset}\right)+{\rm base}
+
+        To understand how the error is measured we can transformed the function above into
+
+        .. math::
+            y = \frac{{\rm amp}}{2} \left(\
+            \cos\right({\rm d}\theta \cdot x\left)\
+            \cos\right({\rm apg} \cdot x - {\rm phase\_offset}\left) -\
+            \sin\right({\rm d}\theta \cdot x\left)\
+            \sin\right({\rm apg} \cdot x - {\rm phase\_offset}\left)
+            \right) + {\rm base}
+
+        When :math:`{\rm apg} \cdot x - {\rm phase\_offset} = (2n + 1) \pi/2` is satisfied the
+        fit model above simplifies to
+
+        .. math::
+            y = \mp \frac{{\rm amp}}{2} \sin\left({\rm d}\theta \cdot x\right) + {\rm base}
+
+        In the limit :math:`{\rm d}\theta \ll 1`, the error can be estimated from the curve data
+
+        .. math::
+            {\rm d}\theta \simeq \mp \frac{2(y - {\rm base})}{x \cdot {\rm amp}}
+
 
     # section: fit_parameters
         defpar \rm amp:
@@ -49,6 +70,8 @@ class ErrorAmplificationAnalysis(curve.CurveAnalysis):
             desc: The angle offset in the gate that we wish to measure.
             init_guess: Multiple initial guesses are tried ranging from -a to a
                 where a is given by :code:`max(abs(angle_per_gate), np.pi / 2)`.
+                Extra guesses are added based on curve data when either :math:`\rm amp` or
+                :math:`\rm base` is :math:`\pi/2`. See fit model for details.
             bounds: [-pi, pi].
 
     # section: note
@@ -63,7 +86,7 @@ class ErrorAmplificationAnalysis(curve.CurveAnalysis):
                 x,
                 amp=0.5 * amp,
                 freq=(d_theta + angle_per_gate) / (2 * np.pi),
-                phase=phase_offset,
+                phase=-phase_offset,
                 baseline=base,
             ),
             plot_color="blue",
@@ -86,7 +109,6 @@ def _default_options(cls):
                 :math:`\pi/2` if the square-root of X gate is added before the repeated gates.
                 This is decided for the user in :meth:`set_schedule` depending on whether the
                 sx gate is included in the experiment.
-            number_of_guesses (int): The number of initial guesses to try.
             max_good_angle_error (float): The maximum angle error for which the fit is
                 considered as good. Defaults to :math:`\pi/2`.
         """
@@ -96,7 +118,6 @@ def _default_options(cls):
         default_options.ylabel = "Population"
         default_options.angle_per_gate = None
         default_options.phase_offset = 0.0
-        default_options.number_guesses = 21
         default_options.max_good_angle_error = np.pi / 2
 
         return default_options
@@ -115,8 +136,6 @@ def _generate_fit_guesses(
         Raises:
             CalibrationError: When ``angle_per_gate`` is missing.
         """
-        n_guesses = self._get_option("number_guesses")
-
         curve_data = self._data()
         max_abs_y, _ = curve.guess.max_height(curve_data.y, absolute=True)
         max_y, min_y = np.max(curve_data.y), np.min(curve_data.y)
@@ -128,15 +147,35 @@ def _generate_fit_guesses(
             user_opt.p0.set_if_empty(amp=max_y - min_y)
             user_opt.bounds.set_if_empty(amp=(-2 * max_abs_y, 2 * max_abs_y))
 
-        # Base the initial guess on the intended angle_per_gate.
-        angle_per_gate = self._get_option("angle_per_gate")
-
-        if angle_per_gate is None:
-            raise CalibrationError("The angle_per_gate was not specified in the analysis options.")
+        # Base the initial guess on the intended angle_per_gate and phase offset.
+        apg = self._get_option("angle_per_gate")
+        phi = self._get_option("phase_offset")
+
+        # Prepare logical guess for specific condition (often satisfied)
+        d_theta_guesses = []
+
+        offsets = apg * curve_data.x + phi
+        amp = user_opt.p0.get("amp", self._get_option("amp"))
+        for i in range(curve_data.x.size):
+            xi = curve_data.x[i]
+            yi = curve_data.y[i]
+            if np.isclose(offsets[i] % np.pi, np.pi / 2) and xi > 0:
+                # Condition satisfied: i.e. cos(apg x - phi) = 0
+                err = -np.sign(np.sin(offsets[i])) * (yi - user_opt.p0["base"]) / (0.5 * amp)
+                # Validate estimate. This is just the first order term of Maclaurin expansion.
+                if np.abs(err) < 0.5:
+                    d_theta_guesses.append(err / xi)
+                else:
+                    # Terminate guess generation because larger d_theta x will start to
+                    # reduce net y value and underestimate the rotation.
+                    break
+
+        # Add naive guess for more coverage
+        guess_range = max(abs(apg), np.pi / 2)
+        d_theta_guesses.extend(np.linspace(-guess_range, guess_range, 11))
 
-        guess_range = max(abs(angle_per_gate), np.pi / 2)
         options = []
-        for d_theta_guess in np.linspace(-guess_range, guess_range, n_guesses):
+        for d_theta_guess in d_theta_guesses:
             new_opt = user_opt.copy()
             new_opt.p0.set_if_empty(d_theta=d_theta_guess)
             options.append(new_opt)

qiskit_experiments/data_processing/__init__.py

@@ -37,6 +37,7 @@
     :toctree: ../stubs/
 
     Probability
+    DirichletProbability
     ToImag
     ToReal
     SVD
@@ -47,6 +48,7 @@
 from .data_action import DataAction, TrainableDataAction
 from .nodes import (
     Probability,
+    DirichletProbability,
     ToImag,
     ToReal,
     SVD,

qiskit_experiments/data_processing/nodes.py

@@ -472,7 +472,7 @@ def _process(
 
             return np.array(populations), np.array(errors)
 
-    def _population_error(self, counts_dict) -> Tuple[float, float]:
+    def _population_error(self, counts_dict: Dict[str, int]) -> Tuple[float, float]:
         """Helper method"""
         shots = sum(counts_dict.values())
         p_mean = counts_dict.get(self._outcome, 0.0) / shots
@@ -481,6 +481,94 @@ def _population_error(self, counts_dict) -> Tuple[float, float]:
         return p_mean, np.sqrt(p_var)
 
 
+class DirichletProbability(Probability):
+    r"""Compute probabilities and variances from a count dictionary.
+
+    This node is a subtype of :py:class:`~qiskit_experiments.data_processing.nodes.Probability`.
+    In the Probability node the variance is computed based on a binomial distribution which can
+    result in a zero variance when the probability is either zero or one.
+    The DirichletProbability node avoids this singularity by assuming a Dirichlet distribution
+    with Bayes update taking a prior distribution. Namely, the mean value is replaced by
+    a mode that represents the most likely value to be sampled
+
+    .. math::
+
+        p = \frac{\alpha_i - 1}{\alpha_0 + K}
+
+    where :math:`\alpha_i = f_i + \theta_i`, :math:`\alpha_0 = \sum_{i=1}^K \alpha_i`,
+    :math:`\theta_i` is the prior distribution and :math:`f_i` is
+    the count number of the :math:`i`-th state out of :math:`K` measurable states.
+    For example, if one provides a 2 bit outcome label, i.e. two-qubit measurement,
+    the possible measurable states are (00, 01, 10, 11) and thus :math:`K` is 4.
+
+    One can provide an arbitrary prior as a dictionary keyed on the outcome label, or
+    a flat prior with an arbitrary float value. It defaults to the Jefferys prior.
+
+    The variance is computed by
+
+    .. math::
+
+        v = \frac{E[x] (1 - E[x])}{\alpha_0 + 1}
+
+    where :math:`E[x] = \alpha_i / \alpha_0` is the mean value of the outcome of interest.
+    With a finite prior, this node always returns a finite variance which prevents
+    unexpected zero divisions.
+    """
+
+    def __init__(
+        self,
+        outcome: str = "1",
+        prior: Union[Dict[str, float], float] = 0.5,
+        validate: bool = True,
+    ):
+        """Initialize a counts to probability data conversion.
+
+        Args:
+            outcome: The bitstring for which to compute the probability which defaults to "1".
+            prior: Prior distribution. This can be a float or a dictionary with keys
+                corresponding to the outcome bitstrings. If n-bit label is provided in
+                ``outcome``, the dimension of the prior distribution, i.e. dictionary length,
+                should be :math:`2^n`. If a float value is applied, this applies a flat prior
+                with the provided value. By default, this assumes a flat prior of 0.5
+                corresponding to the Jefferys prior.
+            validate: If set to False the DataAction will not validate its input.
+
+        Raises:
+            DataProcessorError: When the dimension of the prior and expected parameter vector
+                do not match.
+        """
+        self._dim = 2 ** len(outcome)
+        self._prior = prior
+        super().__init__(outcome=outcome, validate=validate)
+
+        if isinstance(prior, dict) and self._dim != len(prior):
+            raise DataProcessorError(
+                "Dimension of probability density function and prior distribution do not match."
+            )
+
+    def _population_error(self, counts_dict: Dict[str, int]) -> Tuple[float, float]:
+        """Helper method"""
+        shots = sum(counts_dict.values())
+        freq = counts_dict.get(self._outcome, 0.0)
+
+        if isinstance(self._prior, dict):
+            alpha_i = freq + self._prior[self._outcome]
+            alpha_0 = sum([v + self._prior[k] for k, v in counts_dict.items()])
+        else:
+            alpha_i = freq + self._prior
+            alpha_0 = shots + self._prior * self._dim
+
+        p_mean = alpha_i / alpha_0
+        p_var = p_mean * (1 - p_mean) / (alpha_0 + 1)
+        mode = (alpha_i - 1) / (alpha_0 - self._dim)
+
+        # If outcome count is zero or full, mode becomes < 0 or > 1 with a finite prior
+        mode = max(0.0, mode)
+        mode = min(1.0, mode)
+
+        return mode, np.sqrt(p_var)
+
+
 class BasisExpectationValue(DataAction):
     """Compute expectation value of measured basis from probability.
 

qiskit_experiments/data_processing/processor_library.py

@@ -16,7 +16,7 @@
 
 from qiskit_experiments.data_processing.exceptions import DataProcessorError
 from qiskit_experiments.data_processing.data_processor import DataProcessor
-from qiskit_experiments.data_processing.nodes import AverageData, Probability, SVD, MinMaxNormalize
+from qiskit_experiments.data_processing import nodes
 
 
 def get_processor(
@@ -38,16 +38,16 @@ def get_processor(
         DataProcessorError: if the measurement level is not supported.
     """
     if meas_level == MeasLevel.CLASSIFIED:
-        return DataProcessor("counts", [Probability("1")])
+        return DataProcessor("counts", [nodes.DirichletProbability("1")])
 
     if meas_level == MeasLevel.KERNELED:
         if meas_return == "single":
-            processor = DataProcessor("memory", [AverageData(axis=1), SVD()])
+            processor = DataProcessor("memory", [nodes.AverageData(axis=1), nodes.SVD()])
         else:
-            processor = DataProcessor("memory", [SVD()])
+            processor = DataProcessor("memory", [nodes.SVD()])
 
         if normalize:
-            processor.append(MinMaxNormalize())
+            processor.append(nodes.MinMaxNormalize())
 
         return processor
 

qiskit_experiments/library/calibration/analysis/drag_analysis.py

@@ -138,16 +138,19 @@ def _generate_fit_guesses(
             beta=(-freq_bound, freq_bound),
             base=(-max_abs_y, max_abs_y),
         )
-        user_opt.p0.set_if_empty(amp=0.5, base=0.5)
+        user_opt.p0.set_if_empty(base=0.5)
 
         # Drag curves can sometimes be very flat, i.e. averages of y-data
         # and min-max do not always make good initial guesses. We therefore add
-        # 0.5 to the initial guesses.
+        # 0.5 to the initial guesses. Note that we also set amp=-0.5 because the cosine function
+        # becomes +1 at zero phase, i.e. optimal beta, in which y data should become zero
+        # in discriminated measurement level.
         options = []
-        for beta_guess in np.linspace(min_beta, max_beta, 20):
-            new_opt = user_opt.copy()
-            new_opt.p0.set_if_empty(beta=beta_guess)
-            options.append(new_opt)
+        for amp_guess in (0.5, -0.5):
+            for beta_guess in np.linspace(min_beta, max_beta, 20):
+                new_opt = user_opt.copy()
+                new_opt.p0.set_if_empty(amp=amp_guess, beta=beta_guess)
+                options.append(new_opt)
 
         return options
 

qiskit_experiments/library/calibration/analysis/fine_amplitude_analysis.py

@@ -28,4 +28,6 @@ class FineAmplitudeAnalysis(ErrorAmplificationAnalysis):
     """
 
     # The intended angle per gat of the gate being calibrated, e.g. pi for a pi-pulse.
-    __fixed_parameters__ = ["angle_per_gate", "phase_offset"]
+
+    # TODO remove amp from fixed parameter.
+    __fixed_parameters__ = ["angle_per_gate", "phase_offset", "amp"]

qiskit_experiments/library/calibration/fine_amplitude.py

@@ -126,6 +126,16 @@ def _default_experiment_options(cls) -> Options:
 
         return options
 
+    @classmethod
+    def _default_analysis_options(cls) -> Options:
+        """Default analysis options."""
+        options = super()._default_analysis_options()
+
+        # TODO remove amp from fixed parameter.
+        options.amp = 1.0
+
+        return options
+
     def __init__(self, qubit: int):
         """Setup a fine amplitude experiment on the given qubit.
 

qiskit_experiments/library/calibration/fine_drag.py

@@ -160,7 +160,7 @@ def _default_analysis_options(cls) -> Options:
         options = super()._default_analysis_options()
         options.normalization = True
         options.angle_per_gate = 0.0
-        options.phase_offset = -np.pi / 2
+        options.phase_offset = np.pi / 2
         options.amp = 1.0
 
         return options

qiskit_experiments/library/characterization/cr_hamiltonian_analysis.py

@@ -199,7 +199,8 @@ def _default_options(cls):
         """Return the default analysis options."""
         default_options = super()._default_options()
         default_options.data_processor = dp.DataProcessor(
-            input_key="counts", data_actions=[dp.Probability("1"), dp.BasisExpectationValue()]
+            input_key="counts",
+            data_actions=[dp.DirichletProbability("1"), dp.BasisExpectationValue()],
         )
         default_options.curve_plotter = "mpl_multiv_canvas"
         default_options.xlabel = "Flat top width"
@@ -266,21 +267,23 @@ def _generate_fit_guesses(
 
         guesses = defaultdict(list)
         for control in (0, 1):
-            # start from Z oscillation
             x_data = self._data(series_name=f"x|c={control}")
             y_data = self._data(series_name=f"y|c={control}")
             z_data = self._data(series_name=f"z|c={control}")
 
             omega_xyz = []
             for data in (x_data, y_data, z_data):
+                ymin, ymax = np.percentile(data.y, [10, 90])
+                if ymax - ymin < 0.2:
+                    # oscillation amplitude might be almost zero,
+                    # then exclude from average because of lower SNR
+                    continue
                 fft_freq = curve.guess.frequency(data.x, data.y)
-                # oscillation amplitude might be almost zero, then exclude from average
-                if fft_freq > 0:
-                    omega_xyz.append(fft_freq)
+                omega_xyz.append(fft_freq)
             if omega_xyz:
                 omega = 2 * np.pi * np.average(omega_xyz)
             else:
-                omega = 0.0
+                omega = 1e-3
 
             zmin, zmax = np.percentile(z_data.y, [10, 90])
             theta = np.arccos(np.sqrt((zmax - zmin) / 2))

qiskit_experiments/library/randomized_benchmarking/rb_experiment.py

@@ -22,8 +22,8 @@
 from qiskit.quantum_info import Clifford
 from qiskit.circuit import Gate
 
+import qiskit_experiments.data_processing as dp
 from qiskit_experiments.framework import BaseExperiment, ParallelExperiment, Options
-from qiskit_experiments.curve_analysis.data_processing import probability
 from .rb_analysis import RBAnalysis
 from .clifford_utils import CliffordUtils
 from .rb_utils import RBUtils
@@ -88,7 +88,12 @@ def __init__(
 
         # Set configurable options
         self.set_experiment_options(lengths=list(lengths), num_samples=num_samples)
-        self.set_analysis_options(data_processor=probability(outcome="0" * self.num_qubits))
+        self.set_analysis_options(
+            data_processor=dp.DataProcessor(
+                input_key="counts",
+                data_actions=[dp.DirichletProbability(outcome="0" * self.num_qubits)],
+            )
+        )
 
         # Set fixed options
         self._full_sampling = full_sampling