Stark P1 spectroscopy experiment (#1010)

### Summary

This PR adds `StarkP1Spectroscopy` experiment and a dedicated analysis.
User can set amplitudes to scan to run this experiment, and the analysis
class can convert the amplitudes into frequencies if Stark coefficients
are provided. These coefficients can be calibrated with
`StarkRamseyXYAmpScan` experiment introduced in #1009 .

### Details and comments

The test experiment result is available in
[6799e7ae-414f-4816-b887-b4c3ba498624](https://quantum-computing.ibm.com/experiments/6799e7ae-414f-4816-b887-b4c3ba498624).
This experiment result was obtained with the following experiment code:

```python
from qiskit_experiments.library import StarkP1Spectroscopy

exp = StarkP1Spectroscopy((0, ), backend)
exp_data = exp.run().block_for_results()
```

The Stark coefficients can be directly set in the analysis options. By
the default setting, the analysis class searches for the coefficients in
the experiment service. If analysis results for all coefficients are
found (i.e. previously saved in the service), it automatically converts
the amplitudes into frequencies for visualization.

A public class method
`StarkP1SpectAnalysis.retrieve_coefficients_from_service` is also
offered so that a user can set a coefficients dictionary in advance.
This is convenient when the experiment instance is run repeatedly,
because retrieving analysis data from the service causes a communication
overhead. For example,

```python
from qiskit_experiments.library.characterization.analysis import StarkP1SpectAnalysis

# overhead happens only once
coeffs = StarkP1SpectAnalysis.retrieve_coefficients_from_service(
    service=service,
    qubit=0,
    backend=backend.name,
)
exp.analysis.set_options(stark_coefficients=coeffs)

for _ in range(10):
    exp_data = exp.run().block_for_result()
    exp_data.save()
```

User can make a subclass of this analysis class `StarkP1SpectAnalysis`
to perform custom analysis. This built-in class doesn't perform any
analysis except for visualization. Instead, this class provides a
convenient hook `._run_spect_analysis` that takes (x, y, y_err) data and
returns a list of `AnalysisResultData`.

---------

Co-authored-by: Yael Ben-Haim <[email protected]>

qiskit_experiments/library/__init__.py

@@ -160,6 +160,7 @@ class instance to manage parameters and pulse schedules.
 )
 from .characterization import (
     T1,
+    StarkP1Spectroscopy,
     T2Hahn,
     T2Ramsey,
     Tphi,

qiskit_experiments/library/characterization/__init__.py

@@ -24,6 +24,7 @@
     :template: autosummary/experiment.rst
 
     T1
+    StarkP1Spectroscopy
     T2Ramsey
     T2Hahn
     Tphi
@@ -62,6 +63,7 @@
 
     T1Analysis
     T1KerneledAnalysis
+    StarkP1SpectAnalysis
     T2RamseyAnalysis
     T2HahnAnalysis
     TphiAnalysis
@@ -84,6 +86,7 @@
     FineAmplitudeAnalysis,
     RamseyXYAnalysis,
     StarkRamseyXYAmpScanAnalysis,
+    StarkP1SpectAnalysis,
     T2RamseyAnalysis,
     T1Analysis,
     T1KerneledAnalysis,
@@ -98,7 +101,7 @@
     MultiStateDiscriminationAnalysis,
 )
 
-from .t1 import T1
+from .t1 import T1, StarkP1Spectroscopy
 from .qubit_spectroscopy import QubitSpectroscopy
 from .ef_spectroscopy import EFSpectroscopy
 from .t2ramsey import T2Ramsey

qiskit_experiments/library/characterization/analysis/__init__.py

@@ -17,8 +17,7 @@
 from .ramsey_xy_analysis import RamseyXYAnalysis, StarkRamseyXYAmpScanAnalysis
 from .t2ramsey_analysis import T2RamseyAnalysis
 from .t2hahn_analysis import T2HahnAnalysis
-from .t1_analysis import T1Analysis
-from .t1_analysis import T1KerneledAnalysis
+from .t1_analysis import T1Analysis, T1KerneledAnalysis, StarkP1SpectAnalysis
 from .tphi_analysis import TphiAnalysis
 from .cr_hamiltonian_analysis import CrossResonanceHamiltonianAnalysis
 from .readout_angle_analysis import ReadoutAngleAnalysis

qiskit_experiments/library/characterization/analysis/t1_analysis.py

@@ -12,14 +12,20 @@
 """
 T1 Analysis class.
 """
-from typing import Union
+from typing import Union, Tuple, List, Dict
 
 import numpy as np
+from qiskit_ibm_experiment import IBMExperimentService
+from qiskit_ibm_experiment.exceptions import IBMApiError
 from uncertainties import unumpy as unp
 
 import qiskit_experiments.curve_analysis as curve
-from qiskit_experiments.framework import Options
+import qiskit_experiments.data_processing as dp
+import qiskit_experiments.visualization as vis
 from qiskit_experiments.curve_analysis.curve_data import CurveData
+from qiskit_experiments.data_processing.exceptions import DataProcessorError
+from qiskit_experiments.database_service.device_component import Qubit
+from qiskit_experiments.framework import BaseAnalysis, ExperimentData, AnalysisResultData, Options
 
 
 class T1Analysis(curve.DecayAnalysis):
@@ -142,3 +148,212 @@ def _format_data(
 
             return super()._format_data(new_curve_data)
         return super()._format_data(curve_data)
+
+
+class StarkP1SpectAnalysis(BaseAnalysis):
+    """Analysis class for StarkP1Spectroscopy.
+
+    # section: overview
+
+        The P1 landscape is hardly predictable because of the random appearance of
+        lossy TLS notches, and hence this analysis doesn't provide any
+        generic mathematical model to fit the measurement data.
+        A developer may subclass this to conduct own analysis.
+
+        This analysis just visualizes the measured P1 values against Stark tone amplitudes.
+        The tone amplitudes can be converted into the amount of Stark shift
+        when the calibrated coefficients are provided in the analysis option,
+        or the calibration experiment results are available in the result database.
+
+    # section: see_also
+        :class:`qiskit_experiments.library.characterization.ramsey_xy.StarkRamseyXYAmpScan`
+
+    """
+
+    stark_coefficients_names = [
+        "stark_pos_coef_o1",
+        "stark_pos_coef_o2",
+        "stark_pos_coef_o3",
+        "stark_neg_coef_o1",
+        "stark_neg_coef_o2",
+        "stark_neg_coef_o3",
+        "stark_ferr",
+    ]
+
+    @property
+    def plotter(self) -> vis.CurvePlotter:
+        """Curve plotter instance."""
+        return self.options.plotter
+
+    @classmethod
+    def _default_options(cls) -> Options:
+        """Default analysis options.
+
+        Analysis Options:
+            plotter (Plotter): Plotter to visualize P1 landscape.
+            data_processor (DataProcessor): Data processor to compute P1 value.
+            stark_coefficients (Union[Dict, str]): Dictionary of Stark shift coefficients to
+                convert tone amplitudes into amount of Stark shift. This dictionary must include
+                all keys defined in :attr:`.StarkP1SpectAnalysis.stark_coefficients_names`,
+                which are calibrated with :class:`.StarkRamseyXYAmpScan`.
+                Alternatively, it searches for these coefficients in the result database
+                when "latest" is set. This requires having the experiment service set in
+                the experiment data to analyze.
+            x_key (str): Key of the circuit metadata to represent x value.
+        """
+        options = super()._default_options()
+
+        p1spect_plotter = vis.CurvePlotter(vis.MplDrawer())
+        p1spect_plotter.set_figure_options(
+            xlabel="Stark amplitude",
+            ylabel="P(1)",
+            xscale="quadratic",
+        )
+
+        options.update_options(
+            plotter=p1spect_plotter,
+            data_processor=dp.DataProcessor("counts", [dp.Probability("1")]),
+            stark_coefficients="latest",
+            x_key="xval",
+        )
+        return options
+
+    # pylint: disable=unused-argument
+    def _run_spect_analysis(
+        self,
+        xdata: np.ndarray,
+        ydata: np.ndarray,
+        ydata_err: np.ndarray,
+    ) -> List[AnalysisResultData]:
+        """Run further analysis on the spectroscopy data.
+
+        .. note::
+            A subclass can overwrite this method to conduct analysis.
+
+        Args:
+            xdata: X values. This is either amplitudes or frequencies.
+            ydata: Y values. This is P1 values measured at different Stark tones.
+            ydata_err: Sampling error of the Y values.
+
+        Returns:
+            A list of analysis results.
+        """
+        return []
+
+    @classmethod
+    def retrieve_coefficients_from_service(
+        cls,
+        service: IBMExperimentService,
+        qubit: int,
+        backend: str,
+    ) -> Dict:
+        """Retrieve stark coefficient dictionary from the experiment service.
+
+        Args:
+            service: A valid experiment service instance.
+            qubit: Qubit index.
+            backend: Name of the backend.
+
+        Returns:
+            A dictionary of Stark coefficients to convert amplitude to frequency.
+            None value is returned when the dictionary is incomplete.
+        """
+        out = {}
+        try:
+            for name in cls.stark_coefficients_names:
+                results = service.analysis_results(
+                    device_components=[str(Qubit(qubit))],
+                    result_type=name,
+                    backend_name=backend,
+                    sort_by=["creation_datetime:desc"],
+                )
+                if len(results) == 0:
+                    return None
+                result_data = getattr(results[0], "result_data")
+                out[name] = result_data["value"]
+        except (IBMApiError, ValueError, KeyError, AttributeError):
+            return None
+        return out
+
+    def _convert_axis(
+        self,
+        xdata: np.ndarray,
+        coefficients: Dict[str, float],
+    ) -> np.ndarray:
+        """A helper method to convert x-axis.
+
+        Args:
+            xdata: An array of Stark tone amplitude.
+            coefficients: Stark coefficients to convert amplitudes into frequencies.
+
+        Returns:
+            An array of amount of Stark shift.
+        """
+        names = self.stark_coefficients_names  # alias
+        positive = np.poly1d([coefficients[names[idx]] for idx in [2, 1, 0, 6]])
+        negative = np.poly1d([coefficients[names[idx]] for idx in [5, 4, 3, 6]])
+
+        new_xdata = np.where(xdata > 0, positive(xdata), negative(xdata))
+        self.plotter.set_figure_options(
+            xlabel="Stark shift",
+            xval_unit="Hz",
+            xscale="linear",
+        )
+        return new_xdata
+
+    def _run_analysis(
+        self,
+        experiment_data: ExperimentData,
+    ) -> Tuple[List[AnalysisResultData], List["matplotlib.figure.Figure"]]:
+
+        x_key = self.options.x_key
+
+        # Get calibrated Stark tone coefficients
+        if self.options.stark_coefficients == "latest" and experiment_data.service is not None:
+            # Get value from service
+            stark_coeffs = self.retrieve_coefficients_from_service(
+                service=experiment_data.service,
+                qubit=experiment_data.metadata["physical_qubits"][0],
+                backend=experiment_data.backend_name,
+            )
+        elif isinstance(self.options.stark_coefficients, dict):
+            # Get value from experiment options
+            missing = set(self.stark_coefficients_names) - self.options.stark_coefficients.keys()
+            if any(missing):
+                raise KeyError(
+                    "Following coefficient data is missing in the "
+                    f"'stark_coefficients' dictionary: {missing}."
+                )
+            stark_coeffs = self.options.stark_coefficients
+        else:
+            # No calibration is available
+            stark_coeffs = None
+
+        # Compute P1 value and sampling error
+        data = experiment_data.data()
+        try:
+            xdata = np.asarray([datum["metadata"][x_key] for datum in data], dtype=float)
+        except KeyError as ex:
+            raise DataProcessorError(
+                f"X value key {x_key} is not defined in circuit metadata."
+            ) from ex
+        ydata_ufloat = self.options.data_processor(data)
+        ydata = unp.nominal_values(ydata_ufloat)
+        ydata_err = unp.std_devs(ydata_ufloat)
+
+        # Convert x-axis of amplitudes into Stark shift by consuming calibrated parameters.
+        if stark_coeffs:
+            xdata = self._convert_axis(xdata, stark_coeffs)
+
+        # Draw figures and create analysis results.
+        self.plotter.set_series_data(
+            series_name="stark_p1",
+            x_formatted=xdata,
+            y_formatted=ydata,
+            y_formatted_err=ydata_err,
+            x_interp=xdata,
+            y_interp=ydata,
+        )
+        analysis_results = self._run_spect_analysis(xdata, ydata, ydata_err)
+
+        return analysis_results, [self.plotter.figure()]

qiskit_experiments/library/characterization/ramsey_xy.py

@@ -420,7 +420,7 @@ def parameters(self) -> np.ndarray:
             return np.arange(0, max_period, interval)
         return opt.delays
 
-    def parameterized_circuits(self) -> Tuple[QuantumCircuit, QuantumCircuit]:
+    def parameterized_circuits(self) -> Tuple[QuantumCircuit, ...]:
         """Create circuits with parameters for Ramsey XY experiment with Stark tone.
 
         Returns:
@@ -538,10 +538,11 @@ def circuits(self) -> List[QuantumCircuit]:
 
     def _metadata(self) -> Dict[str, any]:
         """Return experiment metadata for ExperimentData."""
-        return {
-            "stark_amp": self.experiment_options.stark_amp,
-            "stark_freq_offset": self.experiment_options.stark_freq_offset,
-        }
+        metadata = super()._metadata()
+        metadata["stark_amp"] = self.experiment_options.stark_amp
+        metadata["stark_freq_offset"] = self.experiment_options.stark_freq_offset
+
+        return metadata
 
 
 class StarkRamseyXYAmpScan(BaseExperiment):
@@ -691,7 +692,7 @@ def parameters(self) -> np.ndarray:
 
         return params
 
-    def parameterized_circuits(self) -> Tuple[QuantumCircuit, QuantumCircuit]:
+    def parameterized_circuits(self) -> Tuple[QuantumCircuit, ...]:
         """Create circuits with parameters for Ramsey XY experiment with Stark tone.
 
         Returns:
@@ -816,7 +817,10 @@ def circuits(self) -> List[QuantumCircuit]:
 
     def _metadata(self) -> Dict[str, any]:
         """Return experiment metadata for ExperimentData."""
-        return {
-            "stark_length": self._timing.pulse_time(time=self.experiment_options.stark_length),
-            "stark_freq_offset": self.experiment_options.stark_freq_offset,
-        }
+        metadata = super()._metadata()
+        metadata["stark_length"] = self._timing.pulse_time(
+            time=self.experiment_options.stark_length
+        )
+        metadata["stark_freq_offset"] = self.experiment_options.stark_freq_offset
+
+        return metadata