Curve analysis with uncertainties package

qiskit_experiments/curve_analysis/__init__.py

@@ -101,6 +101,13 @@
     plot_curve_fit
     plot_errorbar
     plot_scatter
+
+Utilities
+*********
+.. autosummary::
+    :toctree: ../stubs/
+
+    check_if_nominal_significant
 """
 from .curve_analysis import CurveAnalysis
 from .curve_data import CurveData, SeriesDef, FitData, ParameterRepr, FitOptions
@@ -110,6 +117,7 @@
     process_curve_data,
     process_multi_curve_data,
 )
+from .utils import check_if_nominal_significant
 from .visualization import plot_curve_fit, plot_errorbar, plot_scatter, FitResultPlotters
 from . import guess
 from . import fit_function

qiskit_experiments/curve_analysis/curve_data.py

@@ -107,8 +107,8 @@ def fitval(self, key: str) -> UFloat:
             key: Name of parameters to extract.
 
         Returns:
-            UFloat object. This object can be operated as standard
-            Python float object with automatic error propagation.
+            A UFloat object which functions as a standard Python float object
+            but with automatic error propagation.
 
         Raises:
             ValueError: When specified parameter is not defined.

qiskit_experiments/curve_analysis/curve_fit.py

@@ -140,7 +140,10 @@ def fit_func(x, *params):
         # 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]
+        # Ignore correlations, add standard error if finite.
+        fit_params = [
+            ufloat(n, s if np.isfinite(s) else np.nan) for n, s in zip(popt, np.sqrt(np.diag(pcov)))
+        ]
 
     # Calculate the reduced chi-squared for fit
     yfits = fit_func(xdata, *popt)

qiskit_experiments/curve_analysis/fit_function.py

@@ -16,28 +16,27 @@
 # pylint: disable=invalid-name, no-member
 
 import functools
-from typing import Callable
+from typing import Callable, Union
 
 import numpy as np
-from uncertainties import unumpy as unp
+from uncertainties import unumpy as unp, UFloat
 
 
-def calc_uncertainties(fit_func: Callable) -> Callable:
-    """Decolator that typecast y values to float array if input parameters have no error.
+def typecast_returns(fit_func: Callable) -> Callable:
+    """A decorator to typecast y values to a float array if the input parameters have no error.
 
     Args:
-        fit_func: Fit function that may return ufloat array.
+        fit_func: Fit function that returns a ufloat array or an array of float.
 
     Returns:
         Fit function with typecast.
     """
 
     @functools.wraps(fit_func)
-    def _wrapper(x, *args, **kwargs) -> np.ndarray:
+    def _wrapper(x, *args, **kwargs) -> Union[float, UFloat, np.ndarray]:
         yvals = fit_func(x, *args, **kwargs)
         try:
             if isinstance(x, float):
-                # single value
                 return float(yvals)
             return yvals.astype(float)
         except TypeError:
@@ -46,7 +45,7 @@ def _wrapper(x, *args, **kwargs) -> np.ndarray:
     return _wrapper
 
 
-@calc_uncertainties
+@typecast_returns
 def cos(
     x: np.ndarray,
     amp: float = 1.0,
@@ -63,7 +62,7 @@ def cos(
     return amp * unp.cos(2 * np.pi * freq * x + phase) + baseline
 
 
-@calc_uncertainties
+@typecast_returns
 def sin(
     x: np.ndarray,
     amp: float = 1.0,
@@ -80,7 +79,7 @@ def sin(
     return amp * unp.sin(2 * np.pi * freq * x + phase) + baseline
 
 
-@calc_uncertainties
+@typecast_returns
 def exponential_decay(
     x: np.ndarray,
     amp: float = 1.0,
@@ -97,7 +96,7 @@ def exponential_decay(
     return amp * base ** (-lamb * x + x0) + baseline
 
 
-@calc_uncertainties
+@typecast_returns
 def gaussian(
     x: np.ndarray, amp: float = 1.0, sigma: float = 1.0, x0: float = 0.0, baseline: float = 0.0
 ) -> np.ndarray:
@@ -109,7 +108,7 @@ def gaussian(
     return amp * unp.exp(-((x - x0) ** 2) / (2 * sigma ** 2)) + baseline
 
 
-@calc_uncertainties
+@typecast_returns
 def cos_decay(
     x: np.ndarray,
     amp: float = 1.0,
@@ -127,7 +126,7 @@ def cos_decay(
     return exponential_decay(x, lamb=1 / tau) * cos(x, amp=amp, freq=freq, phase=phase) + baseline
 
 
-@calc_uncertainties
+@typecast_returns
 def sin_decay(
     x: np.ndarray,
     amp: float = 1.0,
@@ -145,7 +144,7 @@ def sin_decay(
     return exponential_decay(x, lamb=1 / tau) * sin(x, amp=amp, freq=freq, phase=phase) + baseline
 
 
-@calc_uncertainties
+@typecast_returns
 def bloch_oscillation_x(
     x: np.ndarray, px: float = 0.0, py: float = 0.0, pz: float = 0.0, baseline: float = 0.0
 ):
@@ -163,7 +162,7 @@ def bloch_oscillation_x(
     return (-pz * px + pz * px * unp.cos(w * x) + w * py * unp.sin(w * x)) / (w ** 2) + baseline
 
 
-@calc_uncertainties
+@typecast_returns
 def bloch_oscillation_y(
     x: np.ndarray, px: float = 0.0, py: float = 0.0, pz: float = 0.0, baseline: float = 0.0
 ):
@@ -181,7 +180,7 @@ def bloch_oscillation_y(
     return (pz * py - pz * py * unp.cos(w * x) - w * px * unp.sin(w * x)) / (w ** 2) + baseline
 
 
-@calc_uncertainties
+@typecast_returns
 def bloch_oscillation_z(
     x: np.ndarray, px: float = 0.0, py: float = 0.0, pz: float = 0.0, baseline: float = 0.0
 ):

qiskit_experiments/curve_analysis/standard_analysis/decay.py

@@ -105,7 +105,7 @@ def _evaluate_quality(self, fit_data: curve.FitData) -> Union[str, None]:
 
         criteria = [
             fit_data.reduced_chisq < 3,
-            tau.std_dev is None or tau.std_dev < tau.nominal_value,
+            curve.check_if_nominal_significant(tau),
         ]
 
         if all(criteria):

qiskit_experiments/curve_analysis/standard_analysis/oscillation.py

@@ -115,7 +115,7 @@ def _evaluate_quality(self, fit_data: curve.FitData) -> Union[str, None]:
         criteria = [
             fit_data.reduced_chisq < 3,
             1.0 / 4.0 < fit_freq.nominal_value < 10.0,
-            (np.isnan(fit_freq.std_dev) or (fit_freq.std_dev < fit_freq.nominal_value)),
+            curve.check_if_nominal_significant(fit_freq),
         ]
 
         if all(criteria):
@@ -263,8 +263,8 @@ def _evaluate_quality(self, fit_data: curve.FitData) -> Union[str, None]:
 
         criteria = [
             fit_data.reduced_chisq < 3,
-            np.isnan(tau.std_dev) or tau.std_dev < tau.nominal_value,
-            np.isnan(freq.std_dev) or freq.std_dev < freq.nominal_value,
+            curve.check_if_nominal_significant(tau),
+            curve.check_if_nominal_significant(freq),
         ]
 
         if all(criteria):

qiskit_experiments/curve_analysis/standard_analysis/resonance.py

@@ -134,7 +134,7 @@ def _evaluate_quality(self, fit_data: curve.FitData) -> Union[str, None]:
             1.5 * freq_increment < fit_sigma.nominal_value,
             fit_width_ratio < 0.25,
             fit_data.reduced_chisq < 3,
-            (np.isnan(fit_sigma.std_dev) or fit_sigma.std_dev < fit_sigma.nominal_value),
+            curve.check_if_nominal_significant(fit_sigma),
             snr > 2,
         ]
 

qiskit_experiments/curve_analysis/utils.py

@@ -0,0 +1,47 @@
+# 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.
+
+"""
+A set of utility functions.
+"""
+
+from typing import Union, Optional
+
+import numpy as np
+from uncertainties.core import UFloat
+
+
+def check_if_nominal_significant(
+    val: Union[float, UFloat],
+    fraction: float = 1.0,
+    absolute: Optional[float] = None,
+) -> bool:
+    """Check if the nominal part of the given value is larger than the standard error.
+
+    Args:
+        val: Input value to evaluate.
+        fraction: Valid fraction of the nominal part to its standard error.
+            This function returns ``False`` if the nominal part is
+            smaller than the error by this fraction.
+        absolute: Use this value as a threshold if given.
+
+    Returns:
+        ``True`` if the nominal part is significant.
+    """
+    if isinstance(val, float):
+        return True
+
+    threshold = absolute if absolute is not None else fraction * val.nominal_value
+    if np.isnan(val.std_dev) or val.std_dev < threshold:
+        return True
+
+    return False

qiskit_experiments/curve_analysis/visualization/curves.py

@@ -72,10 +72,19 @@ def plot_curve_fit(
     # Confidence interval of 1 sigma
     stdev_arr = unp.std_devs(ys_fit_with_error)
     if np.isfinite(stdev_arr).all():
+        # 1 sigma plot
         ax.fill_between(
             xs,
             y1=unp.nominal_values(ys_fit_with_error) - stdev_arr,
             y2=unp.nominal_values(ys_fit_with_error) + stdev_arr,
+            alpha=0.3,
+            color=plot_opts["color"],
+        )
+        # 3 sigma plot
+        ax.fill_between(
+            xs,
+            y1=unp.nominal_values(ys_fit_with_error) - 3 * stdev_arr,
+            y2=unp.nominal_values(ys_fit_with_error) + 3 * stdev_arr,
             alpha=0.1,
             color=plot_opts["color"],
         )

qiskit_experiments/database_service/db_fitval.py

@@ -23,7 +23,7 @@ class FitVal:
     This data is serializable with the Qiskit Experiment json serializer.
     """
 
-    # TODO deprecate this (replace with UFloat)
+    # TODO deprecate this (replace with UFloat), see qiskit-experiments#559 for details.
 
     value: float
     stderr: Optional[float] = None

qiskit_experiments/framework/base_analysis.py

@@ -14,7 +14,9 @@
 """
 
 from abc import ABC, abstractmethod
-from typing import List, Tuple
+import copy
+from collections import OrderedDict
+from typing import List, Tuple, Union, Dict
 from uncertainties.core import UFloat
 
 from qiskit_experiments.database_service.device_component import Qubit
@@ -193,7 +195,7 @@ def _format_analysis_result(self, data, experiment_id, experiment_components=Non
             device_components = experiment_components
 
         # Convert ufloat to FitVal so that database service can parse
-        # TODO completely deprecate FitVal. We can store UFloat in database.
+        # TODO completely deprecate FitVal. We can store UFloat in database (qiskit-experiments#559)
         if isinstance(data.value, UFloat):
             value = FitVal(
                 value=data.value.nominal_value,

qiskit_experiments/library/characterization/analysis/drag_analysis.py

@@ -186,7 +186,7 @@ def _evaluate_quality(self, fit_data: curve.FitData) -> Union[str, None]:
         criteria = [
             fit_data.reduced_chisq < 3,
             fit_beta.nominal_value < 1 / fit_freq0.nominal_value,
-            fit_beta.std_dev < abs(fit_beta.nominal_value),
+            curve.check_if_nominal_significant(fit_beta),
         ]
 
         if all(criteria):

qiskit_experiments/library/characterization/analysis/remsey_xy_analysis.py

@@ -162,7 +162,7 @@ def _evaluate_quality(self, fit_data: curve.FitData) -> Union[str, None]:
 
         criteria = [
             fit_data.reduced_chisq < 3,
-            fit_freq.std_dev < abs(fit_freq.nominal_value),
+            curve.check_if_nominal_significant(fit_freq),
         ]
 
         if all(criteria):

qiskit_experiments/library/characterization/analysis/t1_analysis.py

@@ -68,9 +68,9 @@ def _evaluate_quality(self, fit_data: curve.FitData) -> Union[str, None]:
             fit_data.reduced_chisq < 3,
             abs(amp.nominal_value - 1.0) < 0.1,
             abs(base.nominal_value) < 0.1,
-            amp.std_dev is None or amp.std_dev < 0.1,
-            tau.std_dev is None or tau.std_dev < tau.nominal_value,
-            base.std_dev is None or base.std_dev < 0.1,
+            curve.check_if_nominal_significant(amp, absolute=0.1),
+            curve.check_if_nominal_significant(tau),
+            curve.check_if_nominal_significant(base, absolute=0.1),
         ]
 
         if all(criteria):

qiskit_experiments/library/characterization/analysis/t2ramsey_analysis.py

@@ -72,9 +72,9 @@ def _evaluate_quality(self, fit_data: curve.FitData) -> Union[str, None]:
 
         criteria = [
             fit_data.reduced_chisq < 3,
-            amp.std_dev is None or amp.std_dev < 0.1 * amp.nominal_value,
-            tau.std_dev is None or tau.std_dev < 0.1 * tau.nominal_value,
-            freq.std_dev is None or freq.std_dev < 0.1 * freq.nominal_value,
+            curve.check_if_nominal_significant(amp, fraction=0.1),
+            curve.check_if_nominal_significant(tau, fraction=0.1),
+            curve.check_if_nominal_significant(freq, fraction=0.1),
         ]
 
         if all(criteria):

qiskit_experiments/library/quantum_volume/qv_analysis.py

@@ -20,7 +20,7 @@
 import numpy as np
 from uncertainties import ufloat
 
-from qiskit_experiments.framework import BaseAnalysis, AnalysisResultData
+from qiskit_experiments.framework import BaseAnalysis, AnalysisResultData, Options
 from qiskit_experiments.curve_analysis import plot_scatter, plot_errorbar
 
 

qiskit_experiments/library/randomized_benchmarking/rb_utils.py

@@ -16,10 +16,11 @@
 
 from typing import Tuple, Dict, Optional, List, Union, Sequence
 import numpy as np
-from uncertainties.core import UFloat
+from uncertainties.core import UFloat, ufloat
 from qiskit import QiskitError, QuantumCircuit
 from qiskit.providers.backend import Backend
 from qiskit_experiments.database_service.device_component import Qubit
+from qiskit_experiments.database_service.db_fitval import FitVal
 from qiskit_experiments.framework import DbAnalysisResultV1, AnalysisResultData
 
 
@@ -260,7 +261,12 @@ def calculate_2q_epg(
                     gate = result.name.replace("EPG_", "")
 
                     # This keeps variance of previous experiment to propagate
-                    epg_1_qubit_dict[qubit][gate] = result.value
+                    analysis_value = result.value
+                    if isinstance(analysis_value, FitVal):
+                        analysis_value = ufloat(
+                            nominal_value=analysis_value.value, std_dev=analysis_value.stderr
+                        )
+                    epg_1_qubit_dict[qubit][gate] = analysis_value
 
         for key in gate_error_ratio:
             qubits, gate = key

releasenotes/notes/upgrade-curve-fit-4dc01b1db55ee398.yaml

@@ -1,8 +1,9 @@
 ---
 others:
   - |
-    Curve analysis has been updated to consider the fit parameter covariance in the error propagation.
-    Now all result figures show the confidence interval of 1 sigma that is visualized
-    as a semitransparent filled-area around the bold fit line. The width of the area corresponds to the
-    standard deviation computed by the uncertainties of fit parameters with the covariance matrix.
-    All analysis outcomes derived from the fit parameters also have a standard error computed in the same way.
+    Curve analysis is updated to use the covariance between fit parameters in the error propagation.
+    Now, all result figures show a 1 sigma and 3 sigma confidence intervals with
+    semitransparent filled-area with different transparency around the bold fit line.
+    The width of the area corresponds to the standard deviation computed by the
+    uncertainties of fit parameters with the covariance matrix. All analysis outcomes derived from the 
+    fit parameters also have a standard error computed in the same way.