Resonator spectroscopy (qiskit-community#583)

* This PR introduces the resonator spectroscopy experiment which scans the frequency 
  of a pulse applied to the readout resonator. It also switches the default resonance
  analysis to a Lorentzian fit instead of a Gaussian fit. Smaller changes allow developers
  to set resonators as target devices.

Co-authored-by: Naoki Kanazawa <[email protected]>
Co-authored-by: Will Shanks <[email protected]>

qiskit_experiments/curve_analysis/__init__.py

@@ -48,6 +48,7 @@
     DumpedOscillationAnalysis
     OscillationAnalysis
     ResonanceAnalysis
+    GaussianAnalysis
     ErrorAmplificationAnalysis
 
 Functions
@@ -73,6 +74,7 @@
     fit_function.cos_decay
     fit_function.exponential_decay
     fit_function.gaussian
+    fit_function.sqrt_lorentzian
     fit_function.sin
     fit_function.sin_decay
     fit_function.bloch_oscillation_x
@@ -120,5 +122,6 @@
     DumpedOscillationAnalysis,
     OscillationAnalysis,
     ResonanceAnalysis,
+    GaussianAnalysis,
     ErrorAmplificationAnalysis,
 )

qiskit_experiments/curve_analysis/fit_function.py

@@ -77,6 +77,17 @@ def gaussian(
     return amp * np.exp(-((x - x0) ** 2) / (2 * sigma**2)) + baseline
 
 
+def sqrt_lorentzian(
+    x: np.ndarray, amp: float = 1.0, kappa: float = 1.0, x0: float = 0.0, baseline: float = 0.0
+) -> np.ndarray:
+    r"""Square-root Lorentzian function for spectroscopy.
+
+    .. math::
+        y = {\rm amp}{\rm abs}\left(\frac{1}{1 + 2i(x - x0)/\kappa}\right) + {\rm baseline}
+    """
+    return amp * np.abs(1 / (1 + 2.0j * (x - x0) / kappa)) + baseline
+
+
 def cos_decay(
     x: np.ndarray,
     amp: float = 1.0,

qiskit_experiments/curve_analysis/standard_analysis/__init__.py

@@ -14,5 +14,6 @@
 
 from .oscillation import OscillationAnalysis, DumpedOscillationAnalysis
 from .resonance import ResonanceAnalysis
+from .gaussian import GaussianAnalysis
 from .error_amplification_analysis import ErrorAmplificationAnalysis
 from .decay import DecayAnalysis

qiskit_experiments/curve_analysis/standard_analysis/gaussian.py

@@ -0,0 +1,156 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021.
+#
+# 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.
+
+"""Resonance analysis class based on a Gaussian fit."""
+
+from typing import List, Union
+
+import numpy as np
+
+import qiskit_experiments.curve_analysis as curve
+from qiskit_experiments.framework import Options
+
+
+class GaussianAnalysis(curve.CurveAnalysis):
+    r"""A class to analyze a resonance, typically seen as a peak.
+
+    Overview
+        This analysis takes only single series. This series is fit by the Gaussian function.
+
+    Fit Model
+        The fit is based on the following Gaussian function.
+
+        .. math::
+
+            F(x) = a \exp(-(x-f)^2/(2\sigma^2)) + b
+
+    Fit Parameters
+        - :math:`a`: Peak height.
+        - :math:`b`: Base line.
+        - :math:`f`: Center frequency. This is the fit parameter of main interest.
+        - :math:`\sigma`: Standard deviation of Gaussian function.
+
+    Initial Guesses
+        - :math:`a`: Calculated by :func:`~qiskit_experiments.curve_analysis.guess.max_height`.
+        - :math:`b`: Calculated by :func:`~qiskit_experiments.curve_analysis.guess.\
+          constant_spectral_offset`.
+        - :math:`f`: Frequency at max height position calculated by
+          :func:`~qiskit_experiments.curve_analysis.guess.max_height`.
+        - :math:`\sigma`: Calculated from FWHM of peak :math:`w`
+          such that :math:`w / \sqrt{8} \ln{2}`, where FWHM is calculated by
+          :func:`~qiskit_experiments.curve_analysis.guess.full_width_half_max`.
+
+    Bounds
+        - :math:`a`: [-2, 2] scaled with maximum signal value.
+        - :math:`b`: [-1, 1] scaled with maximum signal value.
+        - :math:`f`: [min(x), max(x)] of frequency scan range.
+        - :math:`\sigma`: [0, :math:`\Delta x`] where :math:`\Delta x`
+          represents frequency scan range.
+
+    """
+
+    __series__ = [
+        curve.SeriesDef(
+            fit_func=lambda x, a, sigma, freq, b: curve.fit_function.gaussian(
+                x, amp=a, sigma=sigma, x0=freq, baseline=b
+            ),
+            plot_color="blue",
+            model_description=r"a \exp(-(x-f)^2/(2\sigma^2)) + b",
+        )
+    ]
+
+    @classmethod
+    def _default_options(cls) -> Options:
+        options = super()._default_options()
+        options.result_parameters = [curve.ParameterRepr("freq", "f01", "Hz")]
+        options.normalization = True
+        options.xlabel = "Frequency"
+        options.ylabel = "Signal (arb. units)"
+        options.xval_unit = "Hz"
+        return options
+
+    def _generate_fit_guesses(
+        self, user_opt: curve.FitOptions
+    ) -> Union[curve.FitOptions, List[curve.FitOptions]]:
+        """Compute the initial guesses.
+
+        Args:
+            user_opt: Fit options filled with user provided guess and bounds.
+
+        Returns:
+            List of fit options that are passed to the fitter function.
+        """
+        curve_data = self._data()
+        max_abs_y, _ = curve.guess.max_height(curve_data.y, absolute=True)
+
+        user_opt.bounds.set_if_empty(
+            a=(-2 * max_abs_y, 2 * max_abs_y),
+            sigma=(0, np.ptp(curve_data.x)),
+            freq=(min(curve_data.x), max(curve_data.x)),
+            b=(-max_abs_y, max_abs_y),
+        )
+        user_opt.p0.set_if_empty(b=curve.guess.constant_spectral_offset(curve_data.y))
+
+        y_ = curve_data.y - user_opt.p0["b"]
+
+        _, peak_idx = curve.guess.max_height(y_, absolute=True)
+        fwhm = curve.guess.full_width_half_max(curve_data.x, y_, peak_idx)
+
+        user_opt.p0.set_if_empty(
+            a=curve_data.y[peak_idx] - user_opt.p0["b"],
+            freq=curve_data.x[peak_idx],
+            sigma=fwhm / np.sqrt(8 * np.log(2)),
+        )
+
+        return user_opt
+
+    def _evaluate_quality(self, fit_data: curve.FitData) -> Union[str, None]:
+        """Algorithmic criteria for whether the fit is good or bad.
+
+        A good fit has:
+            - a reduced chi-squared less than 3,
+            - a peak within the scanned frequency range,
+            - a standard deviation that is not larger than the scanned frequency range,
+            - a standard deviation that is wider than the smallest frequency increment,
+            - a signal-to-noise ratio, defined as the amplitude of the peak divided by the
+              square root of the median y-value less the fit offset, greater than a
+              threshold of two, and
+            - a standard error on the sigma of the Gaussian that is smaller than the sigma.
+        """
+        curve_data = self._data()
+
+        max_freq = np.max(curve_data.x)
+        min_freq = np.min(curve_data.x)
+        freq_increment = np.mean(np.diff(curve_data.x))
+
+        fit_a = fit_data.fitval("a").value
+        fit_b = fit_data.fitval("b").value
+        fit_freq = fit_data.fitval("freq").value
+        fit_sigma = fit_data.fitval("sigma").value
+        fit_sigma_err = fit_data.fitval("sigma").stderr
+
+        snr = abs(fit_a) / np.sqrt(abs(np.median(curve_data.y) - fit_b))
+        fit_width_ratio = fit_sigma / (max_freq - min_freq)
+
+        criteria = [
+            min_freq <= fit_freq <= max_freq,
+            1.5 * freq_increment < fit_sigma,
+            fit_width_ratio < 0.25,
+            fit_data.reduced_chisq < 3,
+            (fit_sigma_err is None or fit_sigma_err < fit_sigma),
+            snr > 2,
+        ]
+
+        if all(criteria):
+            return "good"
+
+        return "bad"

qiskit_experiments/curve_analysis/standard_analysis/resonance.py

@@ -21,50 +21,50 @@
 
 
 class ResonanceAnalysis(curve.CurveAnalysis):
-    r"""A class to analyze a resonance, typically seen as a peak.
+    r"""A class to analyze a resonance peak with a square rooted Lorentzian function.
 
     Overview
-        This analysis takes only single series. This series is fit by the Gaussian function.
+        This analysis takes only single series. This series is fit to the square root of
+        a Lorentzian function.
 
     Fit Model
-        The fit is based on the following Gaussian function.
+        The fit is based on the following Lorentzian function.
 
         .. math::
 
-            F(x) = a \exp(-(x-f)^2/(2\sigma^2)) + b
+            F(x) = a{\rm abs}\left(\frac{1}{1 + 2i(x - x0)/\kappa}\right) + b
 
     Fit Parameters
         - :math:`a`: Peak height.
         - :math:`b`: Base line.
-        - :math:`f`: Center frequency. This is the fit parameter of main interest.
-        - :math:`\sigma`: Standard deviation of Gaussian function.
+        - :math:`x0`: Center value. This is typically the fit parameter of interest.
+        - :math:`\kappa`: Linewidth.
 
     Initial Guesses
         - :math:`a`: Calculated by :func:`~qiskit_experiments.curve_analysis.guess.max_height`.
         - :math:`b`: Calculated by :func:`~qiskit_experiments.curve_analysis.guess.\
           constant_spectral_offset`.
-        - :math:`f`: Frequency at max height position calculated by
+        - :math:`x0`: The max height position is calculated by the function
           :func:`~qiskit_experiments.curve_analysis.guess.max_height`.
-        - :math:`\sigma`: Calculated from FWHM of peak :math:`w`
-          such that :math:`w / \sqrt{8} \ln{2}`, where FWHM is calculated by
+        - :math:`\kappa`: Calculated from FWHM of the peak using
           :func:`~qiskit_experiments.curve_analysis.guess.full_width_half_max`.
 
     Bounds
         - :math:`a`: [-2, 2] scaled with maximum signal value.
         - :math:`b`: [-1, 1] scaled with maximum signal value.
-        - :math:`f`: [min(x), max(x)] of frequency scan range.
-        - :math:`\sigma`: [0, :math:`\Delta x`] where :math:`\Delta x`
-          represents frequency scan range.
+        - :math:`f`: [min(x), max(x)] of x-value scan range.
+        - :math:`\kappa`: [0, :math:`\Delta x`] where :math:`\Delta x`
+          represents the x-value scan range.
 
     """
 
     __series__ = [
         curve.SeriesDef(
-            fit_func=lambda x, a, sigma, freq, b: curve.fit_function.gaussian(
-                x, amp=a, sigma=sigma, x0=freq, baseline=b
+            fit_func=lambda x, a, kappa, freq, b: curve.fit_function.sqrt_lorentzian(
+                x, amp=a, kappa=kappa, x0=freq, baseline=b
             ),
             plot_color="blue",
-            model_description=r"a \exp(-(x-f)^2/(2\sigma^2)) + b",
+            model_description=r"a abs(1 / (1 + 2i * (x - x_0) / \kappa)) + b",
         )
     ]
 
@@ -94,7 +94,7 @@ def _generate_fit_guesses(
 
         user_opt.bounds.set_if_empty(
             a=(-2 * max_abs_y, 2 * max_abs_y),
-            sigma=(0, np.ptp(curve_data.x)),
+            kappa=(0, np.ptp(curve_data.x)),
             freq=(min(curve_data.x), max(curve_data.x)),
             b=(-max_abs_y, max_abs_y),
         )
@@ -106,9 +106,9 @@ def _generate_fit_guesses(
         fwhm = curve.guess.full_width_half_max(curve_data.x, y_, peak_idx)
 
         user_opt.p0.set_if_empty(
-            a=curve_data.y[peak_idx] - user_opt.p0["b"],
+            a=(curve_data.y[peak_idx] - user_opt.p0["b"]),
             freq=curve_data.x[peak_idx],
-            sigma=fwhm / np.sqrt(8 * np.log(2)),
+            kappa=fwhm,
         )
 
         return user_opt
@@ -124,7 +124,7 @@ def _evaluate_quality(self, fit_data: curve.FitData) -> Union[str, None]:
             - a signal-to-noise ratio, defined as the amplitude of the peak divided by the
               square root of the median y-value less the fit offset, greater than a
               threshold of two, and
-            - a standard error on the sigma of the Gaussian that is smaller than the sigma.
+            - a standard error on the kappa of the Lorentzian that is smaller than the kappa.
         """
         curve_data = self._data()
 
@@ -135,18 +135,18 @@ def _evaluate_quality(self, fit_data: curve.FitData) -> Union[str, None]:
         fit_a = fit_data.fitval("a").value
         fit_b = fit_data.fitval("b").value
         fit_freq = fit_data.fitval("freq").value
-        fit_sigma = fit_data.fitval("sigma").value
-        fit_sigma_err = fit_data.fitval("sigma").stderr
+        fit_kappa = fit_data.fitval("kappa").value
+        fit_kappa_err = fit_data.fitval("kappa").stderr
 
         snr = abs(fit_a) / np.sqrt(abs(np.median(curve_data.y) - fit_b))
-        fit_width_ratio = fit_sigma / (max_freq - min_freq)
+        fit_width_ratio = fit_kappa / (max_freq - min_freq)
 
         criteria = [
             min_freq <= fit_freq <= max_freq,
-            1.5 * freq_increment < fit_sigma,
+            1.5 * freq_increment < fit_kappa,
             fit_width_ratio < 0.25,
             fit_data.reduced_chisq < 3,
-            (fit_sigma_err is None or fit_sigma_err < fit_sigma),
+            (fit_kappa_err is None or fit_kappa_err < fit_kappa),
             snr > 2,
         ]
 

qiskit_experiments/data_processing/nodes.py

@@ -13,6 +13,7 @@
 """Different data analysis steps."""
 
 from abc import abstractmethod
+from enum import Enum
 from numbers import Number
 from typing import List, Union, Sequence
 
@@ -374,6 +375,22 @@ def _process(self, data: np.ndarray) -> np.ndarray:
         return data[..., 1] * self.scale
 
 
+class ToAbs(IQPart):
+    """IQ data post-processing. Take the absolute value of the IQ point."""
+
+    def _process(self, data: np.array) -> np.array:
+        """Take the absolute value of the IQ data.
+
+        Args:
+            data: An N-dimensional array of complex IQ point as [real, imaginary].
+
+        Returns:
+            A N-1 dimensional array, each entry is the absolute value of the given IQ data.
+        """
+        # pylint: disable=no-member
+        return unp.sqrt(data[..., 0] ** 2 + data[..., 1] ** 2) * self.scale
+
+
 class Probability(DataAction):
     r"""Compute the mean probability of a single measurement outcome from counts.
 
@@ -555,3 +572,12 @@ def _process(self, data: np.ndarray) -> np.ndarray:
             The data that has been processed.
         """
         return 2 * (0.5 - data)
+
+
+class ProjectorType(Enum):
+    """Types of projectors for data dimensionality reduction."""
+
+    SVD = SVD
+    ABS = ToAbs
+    REAL = ToReal
+    IMAG = ToImag