Merge remote-tracking branch 'upstream/main' into T2_Hahn_Echo

docs/_ext/custom_styles/formatter.py

@@ -79,6 +79,18 @@ def format_note(self, lines: List[str]) -> List[str]:
 
         return format_lines
 
+    @_check_no_indent
+    def format_see_also(self, lines: List[str]) -> List[str]:
+        """Format see also section."""
+        text = ".. seealso:: Module(s) "
+
+        modules = []
+        for line in lines:
+            modules.append(f":py:mod:`~{line.lstrip()}`")
+        text += ", ".join(modules)
+
+        return [text, ""]
+
     @_check_no_indent
     def format_tutorial(self, lines: List[str]) -> List[str]:
         """Format tutorial section."""

docs/_ext/custom_styles/styles.py

@@ -163,6 +163,7 @@ class ExperimentDocstring(QiskitExperimentDocstring):
         "run_opts": None,
         "example": load_standard_section,
         "note": load_standard_section,
+        "see_also": load_standard_section,
     }
 
     __formatter__ = ExperimentSectionFormatter
@@ -293,6 +294,7 @@ class AnalysisDocstring(QiskitExperimentDocstring):
         "analysis_opts": None,
         "example": load_standard_section,
         "note": load_standard_section,
+        "see_also": load_standard_section,
     }
 
     __formatter__ = AnalysisSectionFormatter

docs/_ext/custom_styles/utils.py

@@ -140,7 +140,12 @@ def _format_default_options(defaults: Dict[str, Any], indent: str = "") -> List[
         docstring_lines.append("")
         for par, value in defaults.items():
             if callable(value):
-                value_repr = f"Callable {value.__name__}"
+                if value.__class__.__name__ == "function":
+                    # callback function
+                    value_repr = f"Callable {value.__name__}"
+                else:
+                    # class instance with call method
+                    value_repr = repr(value)
             else:
                 value_repr = repr(value)
             docstring_lines.append(indent * 2 + f"{par:<25} := {value_repr}")

qiskit_experiments/curve_analysis/curve_analysis.py

@@ -309,7 +309,8 @@ def _default_options(cls) -> Options:
             axis (AxesSubplot): Optional. A matplotlib axis object to draw.
             xlabel (str): X label of fit result figure.
             ylabel (str): Y label of fit result figure.
-            ylim (Tuple[float, float]): Min and max height limit of fit plot.
+            xlim (Tuple[float, float]): Min and max value of horizontal axis of the fit plot.
+            ylim (Tuple[float, float]): Min and max value of vertical axis of the fit plot.
             xval_unit (str): SI unit of x values. No prefix is needed here.
                 For example, when the x values represent time, this option will be just "s"
                 rather than "ms". In the fit result plot, the prefix is automatically selected
@@ -347,6 +348,7 @@ def _default_options(cls) -> Options:
         options.axis = None
         options.xlabel = None
         options.ylabel = None
+        options.xlim = None
         options.ylim = None
         options.xval_unit = None
         options.yval_unit = None
@@ -730,6 +732,16 @@ def _transpile_options(self, index: int = -1) -> Dict[str, Any]:
             # Ignore experiment metadata or job metadata is not set or key is not found
             return None
 
+    def _extra_metadata(self) -> Dict[str, Any]:
+        """Returns extra metadata.
+
+        Returns:
+            Extra metadata explicitly added by the experiment subclass.
+        """
+        exclude = ["experiment_type", "num_qubits", "physical_qubits", "job_metadata"]
+
+        return {k: v for k, v in self.__experiment_metadata.items() if k not in exclude}
+
     def _data(
         self,
         series_name: Optional[str] = None,
@@ -929,6 +941,7 @@ def _run_analysis(
             fit_options_candidates = [
                 self._format_fit_options(**fit_options) for fit_options in fit_candidates
             ]
+
             fit_results = []
             for fit_options in fit_options_candidates:
                 fit_result = curve_fitter(
@@ -1028,15 +1041,16 @@ def _run_analysis(
         #
         if self._get_option("plot"):
             fit_figure = FitResultPlotters[self._get_option("curve_plotter")].value.draw(
-                curves=[
-                    (ser, self._data(ser.name, "raw_data"), self._data(ser.name, "fit_ready"))
-                    for ser in self.__series__
-                ],
+                series_defs=self.__series__,
+                raw_samples=[self._data(ser.name, "raw_data") for ser in self.__series__],
+                fit_samples=[self._data(ser.name, "fit_ready") for ser in self.__series__],
                 tick_labels={
                     "xval_unit": self._get_option("xval_unit"),
                     "yval_unit": self._get_option("yval_unit"),
                     "xlabel": self._get_option("xlabel"),
                     "ylabel": self._get_option("ylabel"),
+                    "xlim": self._get_option("xlim"),
+                    "ylim": self._get_option("ylim"),
                 },
                 fit_data=fit_result,
                 result_entries=analysis_results,

qiskit_experiments/curve_analysis/curve_data.py

@@ -46,6 +46,9 @@ class SeriesDef:
     # Latex description of this fit model
     model_description: Optional[str] = None
 
+    # Index of canvas if the result figure is multi-panel
+    canvas: Optional[int] = None
+
 
 @dataclasses.dataclass(frozen=True)
 class CurveData:

qiskit_experiments/curve_analysis/fit_function.py

@@ -13,7 +13,7 @@
 """
 A library of fit functions.
 """
-# pylint: disable=invalid-name
+# pylint: disable=invalid-name, line-too-long
 
 import numpy as np
 
@@ -103,3 +103,51 @@ def sin_decay(
         y = {\rm amp} e^{-x/\tau} \sin\left(2 \pi {\rm freq} x + {\rm phase}\right) + {\rm baseline}
     """
     return exponential_decay(x, lamb=1 / tau) * sin(x, amp=amp, freq=freq, phase=phase) + baseline
+
+
+def bloch_oscillation_x(
+    x: np.ndarray, px: float = 0.0, py: float = 0.0, pz: float = 0.0, baseline: float = 0.0
+):
+    r"""Bloch oscillation in x basis.
+
+    .. math::
+        y = \frac{\left( - p_z p_x + p_z p_x \cos (\omega x) + \omega p_y \sin (\omega x) \right)}{\omega^2} + {\rm baseline},
+
+    where :math:`\omega = \sqrt{p_x^2 + p_y^2 + p_z^2}`. The `p_i` stands for the
+    measured probability in `i \in \left\{ X, Y, Z \right\}` basis.
+    """
+    w = np.sqrt(px ** 2 + py ** 2 + pz ** 2)
+
+    return (-pz * px + pz * px * np.cos(w * x) + w * py * np.sin(w * x)) / (w ** 2) + baseline
+
+
+def bloch_oscillation_y(
+    x: np.ndarray, px: float = 0.0, py: float = 0.0, pz: float = 0.0, baseline: float = 0.0
+):
+    r"""Bloch oscillation in y basis.
+
+    .. math::
+        y = \frac{\left( p_z p_y - p_z p_y \cos (\omega x) - \omega p_x \sin (\omega x) \right)}{\omega^2} + {\rm baseline},
+
+    where :math:`\omega = \sqrt{p_x^2 + p_y^2 + p_z^2}`. The `p_i` stands for the
+    measured probability in `i \in \left\{ X, Y, Z \right\}` basis.
+    """
+    w = np.sqrt(px ** 2 + py ** 2 + pz ** 2)
+
+    return (pz * py - pz * py * np.cos(w * x) - w * px * np.sin(w * x)) / (w ** 2) + baseline
+
+
+def bloch_oscillation_z(
+    x: np.ndarray, px: float = 0.0, py: float = 0.0, pz: float = 0.0, baseline: float = 0.0
+):
+    r"""Bloch oscillation in z basis.
+
+    .. math::
+        y = \frac{\left( p_z^2 + (p_x^2 + p_y^2) \cos (\omega x) \right)}{\omega^2} + {\rm baseline},
+
+    where :math:`\omega = \sqrt{p_x^2 + p_y^2 + p_z^2}`. The `p_i` stands for the
+    measured probability in `i \in \left\{ X, Y, Z \right\}` basis.
+    """
+    w = np.sqrt(px ** 2 + py ** 2 + pz ** 2)
+
+    return (pz ** 2 + (px ** 2 + py ** 2) * np.cos(w * x)) / (w ** 2) + baseline

qiskit_experiments/curve_analysis/guess.py

@@ -22,41 +22,85 @@
 from scipy import signal
 
 from qiskit_experiments.exceptions import AnalysisError
+from .data_processing import mean_xy_data
 
 
-def frequency(x: np.ndarray, y: np.ndarray) -> float:
-    """Get frequency of oscillating signal.
+def frequency(
+    x: np.ndarray,
+    y: np.ndarray,
+    filter_window: int = 5,
+    filter_dim: int = 2,
+) -> float:
+    r"""Get frequency of oscillating signal.
+
+    First this tries FFT. If the true value is likely below or near the frequency resolution,
+    the function tries low frequency fit with
+
+    .. math::
+
+        f_{\rm est} = \frac{1}{2\pi {\rm max}\left| y \right|}
+            {\rm max} \left| \frac{dy}{dx} \right|
+
+    given :math:`y = A \cos (2\pi f x + phi)`. In this mode, y data points are
+    smoothed by a Savitzky-Golay filter to protect against outlier points.
 
     .. note::
 
         This function returns always positive frequency.
+        This function is sensitive to the DC offset.
 
     Args:
         x: Array of x values.
         y: Array of y values.
+        filter_window: Window size of Savitzky-Golay filter. This should be odd number.
+        filter_dim: Dimension of Savitzky-Golay filter.
 
     Returns:
         Frequency estimation of oscillation signal.
     """
-    sampling_rates = np.unique(np.diff(x))
-    sampling_rates = sampling_rates[sampling_rates > 0]
+    # TODO averaging and sort should be done by data processor
+
+    # average the same data points
+    x, y, _ = mean_xy_data(x, y)
+
+    # sorted by x to be monotonic increase
+    x, y = np.asarray(sorted(zip(x, y), key=lambda xy: xy[0])).T
+
+    # to run FFT x interval should be identical
+    sampling_interval = np.unique(np.round(np.diff(x), decimals=20))
 
-    if len(sampling_rates) != 1:
-        sampling_rate = np.min(sampling_rates)
-        x_ = np.arange(x[0], x[-1], sampling_rate)
+    if len(sampling_interval) != 1:
+        # resampling with minimum xdata interval
+        sampling_interval = np.min(sampling_interval)
+        x_ = np.arange(x[0], x[-1], sampling_interval)
         y_ = np.interp(x_, xp=x, fp=y)
     else:
-        sampling_rate = sampling_rates[0]
+        sampling_interval = sampling_interval[0]
         x_ = x
         y_ = y
 
     fft_data = np.fft.fft(y_ - np.average(y_))
-    freqs = np.fft.fftfreq(len(x_), sampling_rate)
+    freqs = np.fft.fftfreq(len(x_), sampling_interval)
 
     positive_freqs = freqs[freqs >= 0]
     positive_fft_data = fft_data[freqs >= 0]
 
-    return positive_freqs[np.argmax(np.abs(positive_fft_data))]
+    freq_guess = positive_freqs[np.argmax(np.abs(positive_fft_data))]
+
+    if freq_guess < 1.5 / (sampling_interval * len(x_)):
+        # low frequency fit, use this mode when the estimate is near the resolution
+        y_smooth = signal.savgol_filter(y_, window_length=filter_window, polyorder=filter_dim)
+
+        # no offset is assumed
+        y_amp = max(np.abs(y_smooth))
+
+        if np.isclose(y_amp, 0.0):
+            # no oscillation signal
+            return 0.0
+
+        freq_guess = max(np.abs(np.diff(y_smooth) / sampling_interval)) / (y_amp * 2 * np.pi)
+
+    return freq_guess
 
 
 def max_height(

qiskit_experiments/curve_analysis/visualization/__init__.py

@@ -25,3 +25,4 @@ class FitResultPlotters(Enum):
     """Map the plotter name to the plotters."""
 
     mpl_single_canvas = fit_result_plotters.MplDrawSingleCanvas
+    mpl_multiv_canvas = fit_result_plotters.MplDrawMultiCanvasVstack