Source code for qiskit_experiments.curve_analysis.standard_analysis.gaussian

# 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, Optional

import lmfit
import numpy as np

import qiskit_experiments.curve_analysis as curve
from qiskit_experiments.framework import Options


[docs] 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. """ def __init__( self, name: Optional[str] = None, ): super().__init__( models=[ lmfit.models.ExpressionModel( expr="a * exp(-(x-freq)**2 / (2*sigma**2)) + b", name="gaussian", ) ], name=name, ) @classmethod def _default_options(cls) -> Options: options = super()._default_options() options.plotter.set_figure_options( xlabel="Frequency", ylabel="Signal (arb. units)", xval_unit="Hz", ) options.result_parameters = [curve.ParameterRepr("freq", "f01", "Hz")] options.normalization = True return options def _generate_fit_guesses( self, user_opt: curve.FitOptions, curve_data: curve.ScatterTable, ) -> Union[curve.FitOptions, List[curve.FitOptions]]: """Create algorithmic initial fit guess from analysis options and curve data. Args: user_opt: Fit options filled with user provided guess and bounds. curve_data: Formatted data collection to fit. Returns: List of fit options that are passed to the fitter function. """ 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.CurveFitResult) -> Union[str, None]: """Algorithmic criteria for whether the fit is good or bad. A good fit has: - a reduced chi-squared less than 3 and greater than zero, - 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. """ freq_increment = np.mean(np.diff(fit_data.x_data)) fit_a = fit_data.ufloat_params["a"] fit_b = fit_data.ufloat_params["b"] fit_freq = fit_data.ufloat_params["freq"] fit_sigma = fit_data.ufloat_params["sigma"] snr = abs(fit_a.n) / np.sqrt(abs(np.median(fit_data.y_data) - fit_b.n)) fit_width_ratio = fit_sigma.n / np.ptp(fit_data.x_data) criteria = [ fit_data.x_range[0] <= fit_freq.n <= fit_data.x_range[1], 1.5 * freq_increment < fit_sigma.n, fit_width_ratio < 0.25, 0 < fit_data.reduced_chisq < 3, curve.utils.is_error_not_significant(fit_sigma), snr > 2, ] if all(criteria): return "good" return "bad"