Curve analysis with uncertainties package

qiskit_experiments/curve_analysis/curve_analysis.py

@@ -44,7 +44,6 @@
     BaseAnalysis,
     ExperimentData,
     AnalysisResultData,
-    FitVal,
     Options,
 )
 
@@ -72,8 +71,6 @@ class CurveAnalysis(BaseAnalysis, ABC):
             - ``name``: Name of the curve. This is arbitrary data field, but should be unique.
             - ``plot_color``: String color representation of this series in the plot.
             - ``plot_symbol``: String formatter of the scatter of this series in the plot.
-            - ``plot_fit_uncertainty``: A Boolean signaling whether to plot fit uncertainty
-              for this series in the plot.
 
         - ``__fixed_parameters__``: A list of parameter names fixed during the fitting.
             These parameters should be provided in some way. For example, you can provide
@@ -578,22 +575,17 @@ def _is_target_series(datum, **filters):
 
         x_key = self._get_option("x_key")
         try:
-            x_values = np.asarray([datum["metadata"][x_key] for datum in data], dtype=float)
+            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
 
         if isinstance(data_processor, DataProcessor):
-            y_data = data_processor(data)
-
-            y_nominals = unp.nominal_values(y_data)
-            y_stderrs = unp.std_devs(y_data)
+            ydata = data_processor(data)
         else:
             y_nominals, y_stderrs = zip(*map(data_processor, data))
-
-            y_nominals = np.asarray(y_nominals, dtype=float)
-            y_stderrs = np.asarray(y_stderrs, dtype=float)
+            ydata = unp.uarray(y_nominals, y_stderrs)
 
         # Store metadata
         metadata = np.asarray([datum["metadata"] for datum in data], dtype=object)
@@ -602,7 +594,7 @@ def _is_target_series(datum, **filters):
         shots = np.asarray([datum.get("shots", np.nan) for datum in data])
 
         # Find series (invalid data is labeled as -1)
-        data_index = np.full(x_values.size, -1, dtype=int)
+        data_index = np.full(xdata.size, -1, dtype=int)
         for idx, series_def in enumerate(self.__series__):
             data_matched = np.asarray(
                 [_is_target_series(datum, **series_def.filter_kwargs) for datum in data], dtype=bool
@@ -612,9 +604,9 @@ def _is_target_series(datum, **filters):
         # Store raw data
         raw_data = CurveData(
             label="raw_data",
-            x=x_values,
-            y=y_nominals,
-            y_err=y_stderrs,
+            x=xdata,
+            y=unp.nominal_values(ydata),
+            y_err=unp.std_devs(ydata),
             shots=shots,
             data_index=data_index,
             metadata=metadata,
@@ -983,7 +975,7 @@ def _run_analysis(
             analysis_results.append(
                 AnalysisResultData(
                     name=PARAMS_ENTRY_PREFIX + self.__class__.__name__,
-                    value=FitVal(fit_result.popt, fit_result.popt_err),
+                    value=[p.nominal_value for p in fit_result.popt],
                     chisq=fit_result.reduced_chisq,
                     quality=quality,
                     extra={
@@ -1010,7 +1002,8 @@ def _run_analysis(
                         unit = None
                     result_entry = AnalysisResultData(
                         name=p_repr,
-                        value=fit_result.fitval(p_name, unit),
+                        value=fit_result.fitval(p_name),
+                        unit=unit,
                         chisq=fit_result.reduced_chisq,
                         quality=quality,
                         extra=self._get_option("extra"),

qiskit_experiments/curve_analysis/curve_data.py

@@ -18,9 +18,9 @@
 from typing import Any, Dict, Callable, Union, List, Tuple, Optional, Iterable
 
 import numpy as np
+from uncertainties.core import UFloat
 
 from qiskit_experiments.exceptions import AnalysisError
-from qiskit_experiments.framework import FitVal
 
 
 @dataclasses.dataclass(frozen=True)
@@ -42,9 +42,6 @@ class SeriesDef:
     # Symbol to represent data points of this line.
     plot_symbol: str = "o"
 
-    # Whether to plot fit uncertainty for this line.
-    plot_fit_uncertainty: bool = False
-
     # Latex description of this fit model
     model_description: Optional[str] = None
 
@@ -83,14 +80,11 @@ class FitData:
     """Set of data generated by the fit function."""
 
     # Order sensitive fit parameter values
-    popt: np.ndarray
+    popt: List[UFloat]
 
     # Order sensitive parameter name list
     popt_keys: List[str]
 
-    # Order sensitive fit parameter uncertainty
-    popt_err: np.ndarray
-
     # Covariance matrix
     pcov: np.ndarray
 
@@ -106,26 +100,22 @@ class FitData:
     # Y data range
     y_range: Tuple[float, float]
 
-    def fitval(self, key: str, unit: Optional[str] = None) -> FitVal:
+    def fitval(self, key: str) -> UFloat:
         """A helper method to get fit value object from parameter key name.
 
         Args:
             key: Name of parameters to extract.
-            unit: Optional. Unit of this value.
 
         Returns:
-            FitVal object.
+            UFloat object. This object can be operated as standard
+            Python float object with automatic error propagation.
 
         Raises:
             ValueError: When specified parameter is not defined.
         """
         try:
             index = self.popt_keys.index(key)
-            return FitVal(
-                value=self.popt[index],
-                stderr=self.popt_err[index],
-                unit=unit,
-            )
+            return self.popt[index]
         except ValueError as ex:
             raise ValueError(f"Parameter {key} is not defined.") from ex
 

qiskit_experiments/curve_analysis/curve_fit.py

@@ -15,6 +15,7 @@
 # pylint: disable = invalid-name
 
 from typing import List, Dict, Tuple, Callable, Optional, Union
+from uncertainties import correlated_values, ufloat
 
 import numpy as np
 import scipy.optimize as opt
@@ -62,8 +63,9 @@ def curve_fit(
         ``xrange`` the range of xdata values used for fit.
 
     Raises:
-        AnalysisError: if the number of degrees of freedom of the fit is
-                       less than 1, or the curve fitting fails.
+        AnalysisError:
+            When the number of degrees of freedom of the fit is
+            less than 1, or the curve fitting fails.
 
     .. note::
         ``sigma`` is assumed to be specified in the same units as ``ydata``
@@ -92,7 +94,7 @@ def fit_func(x, *params):
             return func(x, **dict(zip(param_keys, params)))
 
     else:
-        param_keys = None
+        param_keys = [f"p{i}" for i in range(len(p0))]
         param_p0 = p0
         if bounds:
             param_bounds = bounds
@@ -134,7 +136,11 @@ def fit_func(x, *params):
             "scipy.optimize.curve_fit failed with error: {}".format(str(ex))
         ) from ex
 
-    popt_err = np.sqrt(np.diag(pcov))
+    if np.isfinite(pcov).all():
+        # Keep parameter correlations in following analysis steps
+        fit_params = correlated_values(nom_values=popt, covariance_mat=pcov, tags=param_keys)
+    else:
+        fit_params = [ufloat(nom, np.nan) for nom in popt]
 
     # Calculate the reduced chi-squared for fit
     yfits = fit_func(xdata, *popt)
@@ -148,9 +154,8 @@ def fit_func(x, *params):
     ydata_range = np.min(ydata), np.max(ydata)
 
     return FitData(
-        popt=popt,
-        popt_keys=param_keys,
-        popt_err=popt_err,
+        popt=list(fit_params),
+        popt_keys=list(param_keys),
         pcov=pcov,
         reduced_chisq=reduced_chisq,
         dof=dof,

qiskit_experiments/curve_analysis/fit_function.py

@@ -13,11 +13,40 @@
 """
 A library of fit functions.
 """
-# pylint: disable=invalid-name
+# pylint: disable=invalid-name, no-member
+
+import functools
+from typing import Callable
 
 import numpy as np
+from uncertainties import unumpy as unp
+
+
+def calc_uncertainties(fit_func: Callable) -> Callable:
+    """Decolator that typecast y values to float array if input parameters have no error.
+
+    Args:
+        fit_func: Fit function that may return ufloat array.
+
+    Returns:
+        Fit function with typecast.
+    """
+
+    @functools.wraps(fit_func)
+    def _wrapper(x, *args, **kwargs) -> np.ndarray:
+        yvals = fit_func(x, *args, **kwargs)
+        try:
+            if isinstance(x, float):
+                # single value
+                return float(yvals)
+            return yvals.astype(float)
+        except TypeError:
+            return yvals
+
+    return _wrapper
 
 
+@calc_uncertainties
 def cos(
     x: np.ndarray,
     amp: float = 1.0,
@@ -31,9 +60,10 @@ def cos(
         y = {\rm amp} \cdot \cos\left(2 \pi {\rm freq} \cdot x
             + {\rm phase}\right) + {\rm baseline}
     """
-    return amp * np.cos(2 * np.pi * freq * x + phase) + baseline
+    return amp * unp.cos(2 * np.pi * freq * x + phase) + baseline
 
 
+@calc_uncertainties
 def sin(
     x: np.ndarray,
     amp: float = 1.0,
@@ -47,9 +77,10 @@ def sin(
         y = {\rm amp} \cdot \sin\left(2 \pi {\rm freq} \cdot x
             + {\rm phase}\right) + {\rm baseline}
     """
-    return amp * np.sin(2 * np.pi * freq * x + phase) + baseline
+    return amp * unp.sin(2 * np.pi * freq * x + phase) + baseline
 
 
+@calc_uncertainties
 def exponential_decay(
     x: np.ndarray,
     amp: float = 1.0,
@@ -66,6 +97,7 @@ def exponential_decay(
     return amp * base ** (-lamb * x + x0) + baseline
 
 
+@calc_uncertainties
 def gaussian(
     x: np.ndarray, amp: float = 1.0, sigma: float = 1.0, x0: float = 0.0, baseline: float = 0.0
 ) -> np.ndarray:
@@ -74,9 +106,10 @@ def gaussian(
     .. math::
         y = {\rm amp} \cdot \exp \left( - (x - x0)^2 / 2 \sigma^2 \right) + {\rm baseline}
     """
-    return amp * np.exp(-((x - x0) ** 2) / (2 * sigma ** 2)) + baseline
+    return amp * unp.exp(-((x - x0) ** 2) / (2 * sigma ** 2)) + baseline
 
 
+@calc_uncertainties
 def cos_decay(
     x: np.ndarray,
     amp: float = 1.0,
@@ -94,6 +127,7 @@ def cos_decay(
     return exponential_decay(x, lamb=1 / tau) * cos(x, amp=amp, freq=freq, phase=phase) + baseline
 
 
+@calc_uncertainties
 def sin_decay(
     x: np.ndarray,
     amp: float = 1.0,
@@ -111,6 +145,7 @@ def sin_decay(
     return exponential_decay(x, lamb=1 / tau) * sin(x, amp=amp, freq=freq, phase=phase) + baseline
 
 
+@calc_uncertainties
 def bloch_oscillation_x(
     x: np.ndarray, px: float = 0.0, py: float = 0.0, pz: float = 0.0, baseline: float = 0.0
 ):
@@ -123,11 +158,12 @@ def bloch_oscillation_x(
     where :math:`\omega = \sqrt{p_x^2 + p_y^2 + p_z^2}`. The `p_i` stands for the
     measured probability in :math:`i \in \left\{ X, Y, Z \right\}` basis.
     """
-    w = np.sqrt(px ** 2 + py ** 2 + pz ** 2)
+    w = unp.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
+    return (-pz * px + pz * px * unp.cos(w * x) + w * py * unp.sin(w * x)) / (w ** 2) + baseline
 
 
+@calc_uncertainties
 def bloch_oscillation_y(
     x: np.ndarray, px: float = 0.0, py: float = 0.0, pz: float = 0.0, baseline: float = 0.0
 ):
@@ -140,11 +176,12 @@ def bloch_oscillation_y(
     where :math:`\omega = \sqrt{p_x^2 + p_y^2 + p_z^2}`. The `p_i` stands for the
     measured probability in :math:`i \in \left\{ X, Y, Z \right\}` basis.
     """
-    w = np.sqrt(px ** 2 + py ** 2 + pz ** 2)
+    w = unp.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
+    return (pz * py - pz * py * unp.cos(w * x) - w * px * unp.sin(w * x)) / (w ** 2) + baseline
 
 
+@calc_uncertainties
 def bloch_oscillation_z(
     x: np.ndarray, px: float = 0.0, py: float = 0.0, pz: float = 0.0, baseline: float = 0.0
 ):
@@ -157,6 +194,6 @@ def bloch_oscillation_z(
     where :math:`\omega = \sqrt{p_x^2 + p_y^2 + p_z^2}`. The `p_i` stands for the
     measured probability in :math:`i \in \left\{ X, Y, Z \right\}` basis.
     """
-    w = np.sqrt(px ** 2 + py ** 2 + pz ** 2)
+    w = unp.sqrt(px ** 2 + py ** 2 + pz ** 2)
 
-    return (pz ** 2 + (px ** 2 + py ** 2) * np.cos(w * x)) / (w ** 2) + baseline
+    return (pz ** 2 + (px ** 2 + py ** 2) * unp.cos(w * x)) / (w ** 2) + baseline

qiskit_experiments/curve_analysis/standard_analysis/decay.py

@@ -57,7 +57,6 @@ class DecayAnalysis(curve.CurveAnalysis):
             ),
             plot_color="blue",
             model_description=r"amp \exp(-x/tau) + base",
-            plot_fit_uncertainty=True,
         )
     ]
 
@@ -106,7 +105,7 @@ def _evaluate_quality(self, fit_data: curve.FitData) -> Union[str, None]:
 
         criteria = [
             fit_data.reduced_chisq < 3,
-            tau.stderr is None or tau.stderr < tau.value,
+            tau.std_dev is None or tau.std_dev < tau.nominal_value,
         ]
 
         if all(criteria):

qiskit_experiments/curve_analysis/standard_analysis/error_amplification_analysis.py

@@ -191,12 +191,12 @@ def _evaluate_quality(self, fit_data: curve.FitData) -> Union[str, None]:
             - a measured angle error that is smaller than the allowed maximum good angle error.
               This quantity is set in the analysis options.
         """
-        fit_d_theta = fit_data.fitval("d_theta").value
+        fit_d_theta = fit_data.fitval("d_theta")
         max_good_angle_error = self._get_option("max_good_angle_error")
 
         criteria = [
             fit_data.reduced_chisq < 3,
-            abs(fit_d_theta) < abs(max_good_angle_error),
+            abs(fit_d_theta.nominal_value) < abs(max_good_angle_error),
         ]
 
         if all(criteria):