From e10ab72eb711a4dc5d47d6f64a155a41d45f6177 Mon Sep 17 00:00:00 2001 From: Naoki Kanazawa Date: Mon, 22 Aug 2022 22:42:57 +0900 Subject: [PATCH 1/2] update standard analysis classes --- qiskit_experiments/curve_analysis/guess.py | 152 +++++++++++++++--- .../curve_analysis/standard_analysis/decay.py | 18 +-- .../standard_analysis/oscillation.py | 142 +++++++--------- test/curve_analysis/test_guess.py | 45 ++++-- 4 files changed, 224 insertions(+), 133 deletions(-) diff --git a/qiskit_experiments/curve_analysis/guess.py b/qiskit_experiments/curve_analysis/guess.py index fbfc7f981e..24551e2ee7 100644 --- a/qiskit_experiments/curve_analysis/guess.py +++ b/qiskit_experiments/curve_analysis/guess.py @@ -22,6 +22,7 @@ from scipy import signal from qiskit_experiments.exceptions import AnalysisError +from qiskit_experiments.warnings import deprecated_function def frequency( @@ -59,27 +60,26 @@ def frequency( Frequency estimation of oscillation signal. """ # to run FFT x interval should be identical - sampling_interval = np.unique(np.round(np.diff(x), decimals=20)) - - if len(sampling_interval) != 1: - # resampling with minimum xdata interval - sampling_interval = np.min(sampling_interval) - x_ = np.arange(x[0], x[-1], sampling_interval) - y_ = np.interp(x_, xp=x, fp=y) - else: - sampling_interval = sampling_interval[0] + intervals = np.unique(np.diff(x).round(decimals=4)) + try: + dt = float(intervals) x_ = x y_ = y + except TypeError: + # Interpolate when data is nonuniform series. + dt = np.min(intervals) + x_ = np.arange(min(x), max(x), dt) + y_ = np.interp(x_, xp=x, fp=y) fft_data = np.fft.fft(y_ - np.average(y_)) - freqs = np.fft.fftfreq(len(x_), sampling_interval) + freqs = np.fft.fftfreq(len(x_), dt) positive_freqs = freqs[freqs >= 0] positive_fft_data = fft_data[freqs >= 0] freq_guess = positive_freqs[np.argmax(np.abs(positive_fft_data))] - if freq_guess < 1.5 / (sampling_interval * len(x_)): + if freq_guess < 1.5 / (dt * len(x_)): # low frequency fit, use this mode when the estimate is near the resolution y_smooth = signal.savgol_filter(y_, window_length=filter_window, polyorder=filter_dim) @@ -90,7 +90,7 @@ def frequency( # no oscillation signal return 0.0 - freq_guess = max(np.abs(np.diff(y_smooth) / sampling_interval)) / (y_amp * 2 * np.pi) + freq_guess = max(np.abs(np.diff(y_smooth) / dt)) / (y_amp * 2 * np.pi) return freq_guess @@ -162,38 +162,61 @@ def get_height( def exp_decay(x: np.ndarray, y: np.ndarray) -> float: - r"""Get exponential decay parameter from monotonically increasing (decreasing) curve. + r"""Get exponential decay parameter. This assumes following function form. .. math:: - y(x) = e^{\alpha x} + y(x) = e^{\alpha x} F(x) - We can calculate :math:`\alpha` as + Where :math:`F(x)` might be a sinusoidal function or constant value. + + .. note:: + + This function assumes the data is roughly evenly spaced and that + the y data goes through a few periods so that the peak to peak + value early in the data can be compared to the peak to peak later + in the data to estimate the decay constant. + + This function splits the data at the middle of x data, and compare the + 10-90 percentile peak to peak value on the + left-hand :math:`y_{p-p}^R` and right-hand :math:`y_{p-p}^L` side. Namely, .. math:: - \alpha = \log(y(x)) / x + y_{p-p}^R = \exp(\alpha x_R) + y_{p-p}^L = \exp(\alpha x_L) - To find this number, the numpy polynomial fit with ``deg=1`` is used. + and the exponent :math:`\alpha` can be solved by + + .. math:: + + \alpha = \frac{\log y_{p-p}^R / y_{p-p}^L}{x_R - x_L} Args: x: Array of x values. y: Array of y values. Returns: - Decay rate of signal. + Exponent of curve. """ - inds = y > 0 - if np.count_nonzero(inds) < 2: - return 0 + x_median = np.median(x) + i_l = x < x_median + i_r = x > x_median + + y_l_min, y_l_max = np.percentile(y[i_l], [10, 90]) + y_r_min, y_r_max = np.percentile(y[i_r], [10, 90]) + dy_l = y_l_max - y_l_min + dy_r = y_r_max - y_r_min - coeffs = np.polyfit(x[inds], np.log(y[inds]), deg=1) + x_l = np.average(x[i_l]) + x_r = np.average(x[i_r]) - return float(coeffs[0]) + return np.log(dy_r / dy_l) / (x_r - x_l) +@deprecated_function(last_version="0.6", msg="Use exp_decay instead.") def oscillation_exp_decay( x: np.ndarray, y: np.ndarray, @@ -407,3 +430,86 @@ def rb_decay( dx = np.diff(x) return np.average(ry ** (1 / dx)) + + +def sinusoidal_freq_offset_amp( + x: np.ndarray, + y: np.ndarray, + delay: int, +) -> Tuple[float, float]: + r"""Get frequency, offset and amplitude of sinusoidal signal. + + This function simultaneously estimates the frequency and offset + by using two delayed data sets :math:`y(t-\lambda)` and :math:`y(t-2\lambda)` + generated from the input :math:`y(t)` values, where :math:`\lambda` is a + delay parameter specified by the function argument. + See the following paper for the details of the protocol. + + Khac, T.; Vlasov, S. and Iureva, R. (2021). + Estimating the Frequency of the Sinusoidal Signal using + the Parameterization based on the Delay Operators. + DOI: 10.5220/0010536506560660 + + Once offset :math:`y_0` is estimated, the amplitude of the sinusoidal + is obtained by :math:`A = \max |y - y_0|`. + + .. note:: + + This algorithm poorly works for y values containing only less than a half oscillation cycle. + When the parameter cannot be estimated, this switches to the standard FFT analysis + for the frequency and estimates the offset by simply averaging over y values. + + Args: + x: X values. + y: Y values. + delay: Integer parameter to specify the delay samples. Using smaller value + increases (decreases) sensitivity for high (low) frequency signal. + + Returns: + A tuple of frequency, offset and amplitude. + """ + intervals = np.unique(np.diff(x).round(decimals=4)) + try: + dt = float(intervals) + y_ = y + except TypeError: + # Interpolate when data is nonuniform series. + dt = np.min(intervals) + x_ = np.arange(min(x), max(x), dt) + y_ = np.interp(x_, xp=x, fp=y) + + y0 = y_[2 * delay:] + y1 = y_[delay:-delay] + y2 = y_[:-2 * delay] + + # Liner regression form + # psi = xi theta + xi = np.vstack((y1, np.ones_like(y1))).T + psi = (y0 + y2).reshape((y0.size, 1)) + + # Solve linear regression + # theta = Inv(xi^T xi) xi^T psi + xtx = np.dot(xi.T, xi) + theta = np.dot(np.dot(np.linalg.inv(xtx), xi.T), psi) + + theta1 = float(theta[0]) + theta2 = float(theta[1]) + + # Offset + if np.isclose(theta1, 2.0): + offset = np.average(y) + else: + offset = theta2 / (2 - theta1) + + # Frequency + a1 = theta1 / 2.0 + if np.abs(a1) >= 1.0: + # out of arc cosine domain, use FFT + freq = frequency(x, y - offset) + else: + freq = np.arccos(a1) / (delay * 2 * np.pi) / dt + + # Amplitude + amp = max_height(y - offset, absolute=True)[0] + + return freq, offset, amp diff --git a/qiskit_experiments/curve_analysis/standard_analysis/decay.py b/qiskit_experiments/curve_analysis/standard_analysis/decay.py index 276ac26296..e5ab271cb1 100644 --- a/qiskit_experiments/curve_analysis/standard_analysis/decay.py +++ b/qiskit_experiments/curve_analysis/standard_analysis/decay.py @@ -79,19 +79,11 @@ def _generate_fit_guesses( user_opt.p0.set_if_empty(base=curve.guess.min_height(curve_data.y)[0]) alpha = curve.guess.exp_decay(curve_data.x, curve_data.y) - - if alpha != 0.0: - user_opt.p0.set_if_empty( - tau=-1 / alpha, - amp=curve.guess.max_height(curve_data.y)[0] - user_opt.p0["base"], - ) - else: - # Likely there is no slope. Cannot fit constant line with this model. - # Set some large enough number against to the scan range. - user_opt.p0.set_if_empty( - tau=100 * np.max(curve_data.x), - amp=curve.guess.max_height(curve_data.y)[0] - user_opt.p0["base"], - ) + user_opt.p0.set_if_empty( + # alpha might be zero or positive value due to noisy outcome. + tau=-1 / min(alpha, -1 / (100 * max(curve_data.x))), + amp=curve.guess.max_height(curve_data.y)[0] - user_opt.p0["base"], + ) return user_opt def _evaluate_quality(self, fit_data: curve.CurveFitResult) -> Union[str, None]: diff --git a/qiskit_experiments/curve_analysis/standard_analysis/oscillation.py b/qiskit_experiments/curve_analysis/standard_analysis/oscillation.py index 393b16b6f8..c911698c88 100644 --- a/qiskit_experiments/curve_analysis/standard_analysis/oscillation.py +++ b/qiskit_experiments/curve_analysis/standard_analysis/oscillation.py @@ -35,24 +35,25 @@ class OscillationAnalysis(curve.CurveAnalysis): # section: fit_parameters defpar \rm amp: desc: Amplitude of the oscillation. - init_guess: Calculated by :func:`~qiskit_experiments.curve_analysis.guess.max_height`. - bounds: [-2, 2] scaled to the maximum signal value. + init_guess: Calculated by :func:`~qiskit_experiments.curve_analysis.\ + guess.sinusoidal_freq_offset_amp`. + bounds: [-2 * maximum Y, 2 * maximum Y]. defpar \rm base: desc: Base line. init_guess: Calculated by :func:`~qiskit_experiments.curve_analysis.\ - guess.constant_sinusoidal_offset`. - bounds: [-1, 1] scaled to the maximum signal value. + guess.sinusoidal_freq_offset_amp`. + bounds: [-maximum Y, maximum Y]. defpar \rm freq: desc: Frequency of the oscillation. This is the fit parameter of interest. init_guess: Calculated by :func:`~qiskit_experiments.curve_analysis.\ - guess.frequency`. + guess.sinusoidal_freq_offset_amp`. bounds: [0, inf]. defpar \rm phase: desc: Phase of the oscillation. - init_guess: Zero. + init_guess: Multiple points between [-pi, pi]. bounds: [-pi, pi]. """ @@ -92,19 +93,18 @@ def _generate_fit_guesses( phase=(-np.pi, np.pi), base=(-max_abs_y, max_abs_y), ) - user_opt.p0.set_if_empty( - freq=curve.guess.frequency(curve_data.x, curve_data.y), - base=curve.guess.constant_sinusoidal_offset(curve_data.y), - ) - user_opt.p0.set_if_empty( - amp=curve.guess.max_height(curve_data.y - user_opt.p0["base"], absolute=True)[0], - ) options = [] - for phase_guess in np.linspace(0, np.pi, 5): - new_opt = user_opt.copy() - new_opt.p0.set_if_empty(phase=phase_guess) - options.append(new_opt) + for delay_ratio in (0.1, 0.2, 0.3): + for phase_guess in np.linspace(-np.pi, np.pi, 5): + tmp_opt = user_opt.copy() + freq, base, amp = curve.guess.sinusoidal_freq_offset_amp( + x=curve_data.x, + y=curve_data.y, + delay=int(delay_ratio * len(curve_data.x)), + ) + tmp_opt.p0.set_if_empty(freq=freq, base=base, amp=amp, phase=phase_guess) + options.append(tmp_opt) return options @@ -146,31 +146,33 @@ class DampedOscillationAnalysis(curve.CurveAnalysis): # section: fit_parameters defpar \rm amp: - desc: Amplitude. Height of the decay curve. - init_guess: 0.5 - bounds: [0, 1.5], + desc: Amplitude of the oscillation. + init_guess: Calculated by :func:`~qiskit_experiments.curve_analysis.\ + guess.sinusoidal_freq_offset_amp`. + bounds: [-2 * maximum Y, 2 * maximum Y]. defpar \rm base: - desc: Offset. Base line of the decay curve. - init_guess: Determined by :py:func:`~qiskit_experiments.curve_analysis.\ - guess.constant_sinusoidal_offset` - bounds: [0, 1.5] + desc: Base line. + init_guess: Calculated by :func:`~qiskit_experiments.curve_analysis.\ + guess.sinusoidal_freq_offset_amp`. + bounds: [-maximum Y, maximum Y]. + + defpar \rm freq: + desc: Frequency of the oscillation. This is the fit parameter of interest. + init_guess: Calculated by :func:`~qiskit_experiments.curve_analysis.\ + guess.sinusoidal_freq_offset_amp`. + bounds: [0, inf]. defpar \tau: desc: Represents the rate of decay. init_guess: Determined by :py:func:`~qiskit_experiments.curve_analysis.\ - guess.oscillation_exp_decay` - bounds: [0, None] + guess.exp_decay` + bounds: [0, inf] - defpar \rm freq: - desc: Oscillation frequency. - init_guess: Determined by :py:func:`~qiskit_experiments.curve_analysis.guess.frequency` - bounds: [0, 10 freq] - - defpar \phi: - desc: Phase. Relative shift of the sinusoidal function from the origin. - init_guess: Set multiple guesses within [-pi, pi] - bounds: [-pi, pi] + defpar \rm phase: + desc: Phase of the oscillation. + init_guess: Multiple points between [-pi, pi]. + bounds: [-pi, pi]. """ def __init__( @@ -201,59 +203,35 @@ def _generate_fit_guesses( Returns: List of fit options that are passed to the fitter function. """ - user_opt.p0.set_if_empty( - amp=0.5, - base=curve.guess.constant_sinusoidal_offset(curve_data.y), - ) - - # frequency resolution of this scan - df = 1 / ((curve_data.x[1] - curve_data.x[0]) * len(curve_data.x)) - - if user_opt.p0["freq"] is not None: - # If freq guess is provided - freq_guess = user_opt.p0["freq"] - - freqs = [freq_guess] - else: - freq_guess = curve.guess.frequency(curve_data.x, curve_data.y - user_opt.p0["base"]) - - # The FFT might be up to 1/2 bin off - if freq_guess > df: - freqs = [freq_guess - df, freq_guess, freq_guess + df] - else: - freqs = [0.0, freq_guess] - - # Set guess for decay parameter based on estimated frequency - if freq_guess > df: - alpha = curve.guess.oscillation_exp_decay( - curve_data.x, curve_data.y - user_opt.p0["base"], freq_guess=freq_guess - ) - else: - # Very low frequency. Assume standard exponential decay - alpha = curve.guess.exp_decay(curve_data.x, curve_data.y) - - if alpha != 0.0: - user_opt.p0.set_if_empty(tau=-1 / alpha) - else: - # Likely there is no slope. Cannot fit constant line with this model. - # Set some large enough number against to the scan range. - user_opt.p0.set_if_empty(tau=100 * np.max(curve_data.x)) + max_abs_y, _ = curve.guess.max_height(curve_data.y, absolute=True) user_opt.bounds.set_if_empty( - amp=[0, 1.5], - base=[0, 1.5], + amp=(-2 * max_abs_y, 2 * max_abs_y), + freq=(0, np.inf), tau=(0, np.inf), - freq=(0, 10 * freq_guess), - phi=(-np.pi, np.pi), + phase=(-np.pi, np.pi), + base=(-max_abs_y, max_abs_y), ) + alpha = curve.guess.exp_decay(curve_data.x, curve_data.y) + tau_guess = -1 / min(alpha, -1 / (100 * max(curve_data.x))) - # more robust estimation options = [] - for freq in freqs: - for phi in np.linspace(-np.pi, np.pi, 5)[:-1]: - new_opt = user_opt.copy() - new_opt.p0.set_if_empty(freq=freq, phi=phi) - options.append(new_opt) + for delay_ratio in (0.1, 0.2, 0.3): + for phase_guess in np.linspace(0, np.pi, 5): + tmp_opt = user_opt.copy() + freq, base, amp = curve.guess.sinusoidal_freq_offset_amp( + x=curve_data.x, + y=curve_data.y, + delay=int(delay_ratio * len(curve_data.x)), + ) + tmp_opt.p0.set_if_empty( + freq=freq, + base=base, + amp=amp, + phase=phase_guess, + tau=tau_guess, + ) + options.append(tmp_opt) return options diff --git a/test/curve_analysis/test_guess.py b/test/curve_analysis/test_guess.py index bff8870628..e5255d9288 100644 --- a/test/curve_analysis/test_guess.py +++ b/test/curve_analysis/test_guess.py @@ -110,24 +110,21 @@ def test_exp_decay(self, alpha: float): self.assertAlmostEqualAbsolute(alpha_guess, alpha) - def test_exp_decay_with_invalid_y(self): - """Test when invalid y data is input to exp curve init guess.""" - x = np.array([9.0e-06, 1.9e-05, 2.9e-05, 3.9e-05]) - y = np.array([0.16455749, 0.07045296, 0.02702439, -0.00135192]) - - # The last point is excluded. This point might be some artifact due to filtering. - alpha_guess = guess.exp_decay(x, y) - - np.testing.assert_almost_equal(alpha_guess, -90326, decimal=0) - - @data([1.2, 1.4], [-0.6, 2.5], [0.1, 2.3], [3.5, 1.1], [-4.1, 6.5], [3.0, 1.2]) + @data( + [-1.2, 1.4, 0.0], + [-1.6, 2.5, 0.2], + [-1.3, 2.3, -0.3], + [-3.5, 1.1, 0.5], + [-4.1, 6.5, 2.4], + [-3.0, 1.2, -0.4], + ) @unpack - def test_exp_osci_decay(self, alpha, freq): + def test_exp_decay_with_oscillation(self, alpha, freq, base): """Test of exponential decay guess with oscillation.""" x = np.linspace(0, 1, 100) - y = np.exp(alpha * x) * np.cos(2 * np.pi * freq * x) + y = np.exp(alpha * x) * np.cos(2 * np.pi * freq * x) + base - alpha_guess = guess.oscillation_exp_decay(x, y) + alpha_guess = guess.exp_decay(x, y) self.assertAlmostEqualAbsolute(alpha_guess, alpha) @@ -206,4 +203,22 @@ def test_rb_decay(self, a, b, alpha): alpha_guess = guess.rb_decay(x, y, b=b) - self.assertAlmostEqual(alpha, alpha_guess, delta=alpha * 0.1) + self.assertAlmostEqualAbsolute(alpha, alpha_guess) + + @data( + [0.5, 0.0, 2.0, 0.0, 0.5], + [0.2, -0.1, 1.3, 0.3, 0.2], + [0.6, 0.0, 0.7, 1.5, -0.4], + [0.3, -0.9, 0.6, 2.3, 0.9], + [0.4, -0.1, 1.6, -0.3, -1.3], + ) + @unpack + def test_sinusoidal_freq_offset_amp(self, amp, alpha, freq, phase, base): + """Test simulataneous freq, offset and amp guess.""" + x = np.linspace(0, 1, 100) + y = amp * np.exp(alpha * x) * np.cos(2 * np.pi * freq * x + phase) + base + + freq_guess, base_guess, amp_guess = guess.sinusoidal_freq_offset(x, y, 15) + self.assertAlmostEqualAbsolute(freq, freq_guess) + self.assertAlmostEqualAbsolute(base, base_guess) + self.assertAlmostEqualAbsolute(amp, amp_guess) From 917f4ba1fec8d4ce2625bc64c3d11dffcd8014ef Mon Sep 17 00:00:00 2001 From: Naoki Kanazawa Date: Mon, 22 Aug 2022 23:14:43 +0900 Subject: [PATCH 2/2] Update experiment classes with one stop algorithm --- qiskit_experiments/curve_analysis/__init__.py | 4 + .../curve_analysis/curve_analysis.py | 1 - qiskit_experiments/curve_analysis/guess.py | 363 +++++++++++++++--- .../standard_analysis/bloch_trajectory.py | 2 +- .../curve_analysis/standard_analysis/decay.py | 1 - .../standard_analysis/oscillation.py | 89 +++-- .../analysis/drag_analysis.py | 62 +-- .../analysis/ramsey_xy_analysis.py | 121 +++--- .../analysis/t2ramsey_analysis.py | 19 +- .../library/characterization/ramsey_xy.py | 12 - .../notes/0_5_module_curve_analysis.yaml | 31 ++ test/curve_analysis/test_guess.py | 114 +++++- .../library/characterization/test_t2ramsey.py | 31 +- 13 files changed, 611 insertions(+), 239 deletions(-) create mode 100644 releasenotes/notes/0_5_module_curve_analysis.yaml diff --git a/qiskit_experiments/curve_analysis/__init__.py b/qiskit_experiments/curve_analysis/__init__.py index 1eb7333dd6..7544c1720f 100644 --- a/qiskit_experiments/curve_analysis/__init__.py +++ b/qiskit_experiments/curve_analysis/__init__.py @@ -535,6 +535,10 @@ def _create_analysis_results(self, fit_data, quality, **metadata): guess.max_height guess.min_height guess.oscillation_exp_decay + guess.frequency_lorentz_fit + guess.low_frequency_limit + guess.sinusoidal_freq_offset + guess.composite_sinusoidal_estimate Utilities ********* diff --git a/qiskit_experiments/curve_analysis/curve_analysis.py b/qiskit_experiments/curve_analysis/curve_analysis.py index 949ccfdd15..46c5294380 100644 --- a/qiskit_experiments/curve_analysis/curve_analysis.py +++ b/qiskit_experiments/curve_analysis/curve_analysis.py @@ -437,7 +437,6 @@ def _objective(_params): max=bounds[1], vary=name not in fixed_parameters, ) - try: new = lmfit.minimize( fcn=_objective, diff --git a/qiskit_experiments/curve_analysis/guess.py b/qiskit_experiments/curve_analysis/guess.py index 24551e2ee7..433be3a8a3 100644 --- a/qiskit_experiments/curve_analysis/guess.py +++ b/qiskit_experiments/curve_analysis/guess.py @@ -16,22 +16,51 @@ # pylint: disable=invalid-name import functools -from typing import Optional, Tuple, Callable +from typing import Optional, Tuple, Callable, Iterator import numpy as np -from scipy import signal +from scipy import signal, optimize from qiskit_experiments.exceptions import AnalysisError from qiskit_experiments.warnings import deprecated_function +def _require_uniform_data(func) -> Callable: + """A decorator to convert X data into uniformly spaced series. + + Y data is resampled with new X data with interpolation. + The first and second argument of the decorated function must be x and y. + + Args: + func: Function to decorate. + + Returns: + Decorated function with uniform X, Y input. + """ + + @functools.wraps(func) + def wrap_guess(x, y, *args, **kwargs): + intervals = np.unique(np.diff(x).round(decimals=16)) + if len(intervals) == 1: + x_int = x + y_int = y + else: + x_int = np.arange(min(x), max(x), np.min(intervals)) + y_int = np.interp(x_int, xp=x, fp=y) + return func(x_int, y_int, *args, **kwargs) + + return wrap_guess + + +@deprecated_function(last_version="0.6", msg="Use frequency_lorentz_fit instead.") +@_require_uniform_data def frequency( x: np.ndarray, y: np.ndarray, filter_window: int = 5, filter_dim: int = 2, ) -> float: - r"""Get frequency of oscillating signal. + r"""Deprecated. Get frequency of oscillating signal. First this tries FFT. If the true value is likely below or near the frequency resolution, the function tries low frequency fit with @@ -59,29 +88,18 @@ def frequency( Returns: Frequency estimation of oscillation signal. """ - # to run FFT x interval should be identical - intervals = np.unique(np.diff(x).round(decimals=4)) - try: - dt = float(intervals) - x_ = x - y_ = y - except TypeError: - # Interpolate when data is nonuniform series. - dt = np.min(intervals) - x_ = np.arange(min(x), max(x), dt) - y_ = np.interp(x_, xp=x, fp=y) - - fft_data = np.fft.fft(y_ - np.average(y_)) - freqs = np.fft.fftfreq(len(x_), dt) + dt = x[1] - x[0] + fft_data = np.fft.fft(y - np.average(y)) + freqs = np.fft.fftfreq(len(x), dt) positive_freqs = freqs[freqs >= 0] positive_fft_data = fft_data[freqs >= 0] freq_guess = positive_freqs[np.argmax(np.abs(positive_fft_data))] - if freq_guess < 1.5 / (dt * len(x_)): + if freq_guess < 1.5 / (dt * len(x)): # low frequency fit, use this mode when the estimate is near the resolution - y_smooth = signal.savgol_filter(y_, window_length=filter_window, polyorder=filter_dim) + y_smooth = signal.savgol_filter(y, window_length=filter_window, polyorder=filter_dim) # no offset is assumed y_amp = max(np.abs(y_smooth)) @@ -95,6 +113,122 @@ def frequency( return freq_guess +@_require_uniform_data +def frequency_lorentz_fit( + x: np.ndarray, + y: np.ndarray, + fit_range: int = 5, +) -> float: + """Get oscilaltion frequency of sinusoidal. + + This function estimates frequency with FFT. + The FFT peak location is fine tuned by the Lorentzian fitting to give precise estimate. + When the estimated frequency is smaller than 150 % of a FFT frequency bin, + this returns the estimate by :func:`low_frequency_limit` instead. + + .. note:: + + The offset of y data must be subtracted otherwise this may return zero frequency. + + Args: + x: Array of x values. + y: Array of y values. + fit_range: Data points used for Lorentzian fitting. + + Returns: + Frequency estimate. + """ + dt = x[1] - x[0] + + fft_y_re = np.fft.fft(np.real(y)) + fft_y_im = np.fft.fft(np.imag(y)) + + fft_x = np.fft.fftshift(np.fft.fftfreq(len(x), dt)) + fft_y = np.fft.fftshift(fft_y_re + 1j * fft_y_im) + + # Use only positive side. + pos_inds = fft_x > 0 + fft_x = fft_x[pos_inds] + fft_y = fft_y[pos_inds] + + peak_ind = np.argmax(np.abs(fft_y)) + + # Fit a small Lorentzian to the FFT to fine tune the frequency + inds = list(range(max(0, peak_ind - fit_range), min(fft_x.size, peak_ind + fit_range))) + fit_x = fft_x[inds] + fit_y = fft_y[inds] + + def objective(p): + y_complex = (p[0] + 1j * p[1]) / (1 + 1j * (fit_x - p[2]) / p[3]) + y_out_concat = np.concatenate((y_complex.real, y_complex.imag)) + y_ref_concat = np.concatenate((fit_y.real, fit_y.imag)) + return y_out_concat - y_ref_concat + + guess = [ + np.real(fft_y[peak_ind]), + np.imag(fft_y[peak_ind]), + fft_x[peak_ind], + (fft_x[1] - fft_x[0]) / 3, + ] + fit_out, _ = optimize.leastsq( + func=objective, + x0=guess, + ) + # In principle we can estimate phase from arctan of p1 and p0. + # However the accuracy is not quite good because of poor resolution of imaginary part + # due to limitation of the experimental resource, e.g. limited job circuits. + freq = float(fit_out[2]) + + if freq < 1.5 / (dt * len(x)): + return low_frequency_limit(x, y) + return freq + + +@_require_uniform_data +def low_frequency_limit( + x: np.ndarray, + y: np.ndarray, +) -> float: + r"""Get low frequency estimate of sinusoidal. + + When the measured data contains only less than 1 period of the sinusoidal signal, + usually FFT doesn't report precise estimate. Instead, this function estimate + such a low frequency as follows. + + .. math:: + + f_{\rm est} = \frac{1}{2\pi {\rm max}\left| y \right|} + {\rm max} \left| \frac{dy}{dx} \right| + + given :math:`y = A \cos (2\pi f x + phi)`. In this mode, y data points are + smoothed by the Savitzky-Golay filter to protect against outlier points. + + .. note:: + + The offset of y data must be subtracted otherwise this may return poor estimate. + + Args: + x: Array of x values. + y: Array of y values. + + Returns: + Frequency estimate. + """ + dt = x[1] - x[0] + + # Assume x is less than 1 period. A signal of the quarter period can be well approximated by + # the 2nd order polynominal, and thus the windows size can be 1/4 len(x). + y_ = signal.savgol_filter(y, window_length=int(x.size / 4), polyorder=2) + + y_amp = max(np.abs(y_)) + if np.isclose(y_amp, 0.0): + # no oscillation signal + return 0.0 + freq = max(np.abs(np.diff(y_) / dt)) / (y_amp * 2 * np.pi) + + return freq + + def max_height( y: np.ndarray, percentile: Optional[float] = None, @@ -201,15 +335,26 @@ def exp_decay(x: np.ndarray, y: np.ndarray) -> float: Returns: Exponent of curve. """ - x_median = np.median(x) - i_l = x < x_median - i_r = x > x_median + if x.size < 5: + return 0 + + if np.unique(np.diff(x).round(decimals=16)).size == 1: + x0 = np.median(x) + else: + x0 = np.min(x) + np.ptp(x) / 2 + + i_l = x < x0 + i_r = x > x0 y_l_min, y_l_max = np.percentile(y[i_l], [10, 90]) y_r_min, y_r_max = np.percentile(y[i_r], [10, 90]) dy_l = y_l_max - y_l_min dy_r = y_r_max - y_r_min + if np.isclose(dy_l, dy_r): + # Avoid ZeroDiv. When y is flat, dy_l ~ dy_r ~ 0. + return 0 + x_l = np.average(x[i_l]) x_r = np.average(x[i_r]) @@ -224,7 +369,7 @@ def oscillation_exp_decay( filter_dim: int = 2, freq_guess: Optional[float] = None, ) -> float: - r"""Get exponential decay parameter from oscillating signal. + r"""Deprecated. Get exponential decay parameter from oscillating signal. This assumes following function form. @@ -269,7 +414,13 @@ def oscillation_exp_decay( x_peaks = x[peak_pos] y_peaks = y_smoothed[peak_pos] - return exp_decay(x_peaks, y_peaks) + inds = y_peaks > 0 + if np.count_nonzero(inds) < 2: + return 0 + + coeffs = np.polyfit(x_peaks[inds], np.log(y_peaks[inds]), deg=1) + + return float(coeffs[0]) def full_width_half_max( @@ -432,12 +583,13 @@ def rb_decay( return np.average(ry ** (1 / dx)) -def sinusoidal_freq_offset_amp( +@_require_uniform_data +def sinusoidal_freq_offset( x: np.ndarray, y: np.ndarray, delay: int, ) -> Tuple[float, float]: - r"""Get frequency, offset and amplitude of sinusoidal signal. + r"""Get frequency and offset of sinusoidal signal. This function simultaneously estimates the frequency and offset by using two delayed data sets :math:`y(t-\lambda)` and :math:`y(t-2\lambda)` @@ -450,14 +602,9 @@ def sinusoidal_freq_offset_amp( the Parameterization based on the Delay Operators. DOI: 10.5220/0010536506560660 - Once offset :math:`y_0` is estimated, the amplitude of the sinusoidal - is obtained by :math:`A = \max |y - y_0|`. - .. note:: This algorithm poorly works for y values containing only less than a half oscillation cycle. - When the parameter cannot be estimated, this switches to the standard FFT analysis - for the frequency and estimates the offset by simply averaging over y values. Args: x: X values. @@ -466,21 +613,16 @@ def sinusoidal_freq_offset_amp( increases (decreases) sensitivity for high (low) frequency signal. Returns: - A tuple of frequency, offset and amplitude. + A tuple of frequency and offset estimate. + + Raises: + AnalysisError: When parameters cannot be estimated with current delay parameter. """ - intervals = np.unique(np.diff(x).round(decimals=4)) - try: - dt = float(intervals) - y_ = y - except TypeError: - # Interpolate when data is nonuniform series. - dt = np.min(intervals) - x_ = np.arange(min(x), max(x), dt) - y_ = np.interp(x_, xp=x, fp=y) + dt = x[1] - x[0] - y0 = y_[2 * delay:] - y1 = y_[delay:-delay] - y2 = y_[:-2 * delay] + y0 = y[2 * delay :] + y1 = y[delay:-delay] + y2 = y[: -2 * delay] # Liner regression form # psi = xi theta @@ -495,21 +637,124 @@ def sinusoidal_freq_offset_amp( theta1 = float(theta[0]) theta2 = float(theta[1]) - # Offset - if np.isclose(theta1, 2.0): - offset = np.average(y) - else: - offset = theta2 / (2 - theta1) + if np.abs(theta1) >= 2.0: + raise AnalysisError("Invalid estimate. Try another delay parameter.") - # Frequency - a1 = theta1 / 2.0 - if np.abs(a1) >= 1.0: - # out of arc cosine domain, use FFT - freq = frequency(x, y - offset) - else: - freq = np.arccos(a1) / (delay * 2 * np.pi) / dt + offset = theta2 / (2 - theta1) + freq = np.arccos(theta1 / 2.0) / (delay * 2 * np.pi) / dt + + return freq, offset + + +def composite_sinusoidal_estimate( + x: np.ndarray, + y: np.ndarray, + amp: Optional[float] = None, + freq: Optional[float] = None, + base: Optional[float] = None, + phase: Optional[float] = None, +) -> Iterator[Tuple[float, ...]]: + r"""Composite estimate function of full sinusoidal parameters. + + This function switches the guess algorithm based on the situation of pre estimates. + This function often generates multiple guesses, and thus the output is an iterator. + The behavior of this function is roughly classified into following patterns. + + 1. When amp, freq, base, phase are all provided. + + Return pre estimates as-is. + + 2. When freq and base are not knwon. - # Amplitude - amp = max_height(y - offset, absolute=True)[0] + In this situation the function uses :func:`sinusoidal_freq_offset` that + simulataneously estimates the freq and base offset. This generates multiple guesses with + different delay parameters (from 10 t0 30 percent of full X range) so that + it can cover wider frequency range. - return freq, offset, amp + 3. When freq is not known but base is provided. + + In this situation we can precisely eliminate base offset that may harm FFT analysis. + The function uses the :func:`frequency_lorentz_fit` to + get the freq estimate with unbiased Y values computed with the provided base. + + 4. When freq is provided but base is not known. + + In this situation we can estimate the period of the signal. + When X value range is larger than the signal period, + the middle value of the 5-95 percentile of Y values is used as the base. + Otherwise, multiple base candidates are generated by adding delta to the mid Y value. + This delta ranges from -30 to 30 percent of the maximum absolute Y value. + + When the amp is not knwon, it is estimated by :math:`\max |y - y_0|` where :math:`y_0` + is the estimated base offset. When the phase is not known, further multiple guesses + are generated with multiple phase candidates ranging from -pi to pi. + + Args: + x: X values. + y: Y values. + amp: Pre estimate of amplitude. + freq: Pre estimate of frequency. + base: Pre estimate of base offset. + phase: Pre estimate of phase offset. + + Yields: + amplitude, frequency, base offset, phase offset of sinusoidal. + """ + + def yield_multiple_phase(_amp, _freq, _base): + for phase_est in np.linspace(-np.pi, np.pi, 5): + yield _amp, _freq, _base, phase_est + + if freq is None: + if base is None: + # When offset value is not provided, FFT is not reasonable approach. + # Especially when the signal is low frequency, usually the peak position of FFT is + # affected by the offset which is hardly estimated precisely. + # Now we switch to a technique to simulataneously estimate offset and frequency. + for r in [0.1, 0.2, 0.3]: + # Use delay parameter of 10-30% of full range. + try: + freq_est, base_est = sinusoidal_freq_offset(x, y, int(x.size * r)) + except AnalysisError: + continue + amp_est = amp or max_height(y - base_est, absolute=True)[0] + if phase is not None: + yield amp_est, freq_est, base_est, phase + else: + yield from yield_multiple_phase(amp_est, freq_est, base_est) + else: + # Use FFT approach since we can propertly eliminate offset. + freq_est = frequency_lorentz_fit(x, y - base) + amp_est = amp or max_height(y - base, absolute=True)[0] + if phase is not None: + yield amp_est, freq_est, base, phase + else: + yield from yield_multiple_phase(amp_est, freq_est, base) + else: + if base is None: + min_y, max_y = np.percentile(y, [5, 95]) + mid_y = 0.5 * (max_y - min_y) + if np.ptp(x) > 1 / freq: + # It contains more than 1 cycle and thus this is half of peak-to-peak. + base_est = mid_y + amp_est = amp or max_height(y - base_est, absolute=True)[0] + if phase is not None: + yield amp_est, freq, base_est, phase + else: + yield from yield_multiple_phase(amp_est, freq, base_est) + else: + # This is likely low frequency. Provide multiple offset guess for safety. + abs_y = np.max(np.abs(y)) + for delta_y in 0.3 * np.linspace(-abs_y, abs_y, 5): + base_est = mid_y + delta_y + amp_est = amp or max_height(y - base_est, absolute=True)[0] + if phase is not None: + yield amp_est, freq, base_est, phase + else: + yield from yield_multiple_phase(amp_est, freq, base_est) + else: + amp_est = amp or max_height(y - base, absolute=True)[0] + if phase is not None: + yield amp_est, freq, base, phase + else: + yield from yield_multiple_phase(amp_est, freq, base) diff --git a/qiskit_experiments/curve_analysis/standard_analysis/bloch_trajectory.py b/qiskit_experiments/curve_analysis/standard_analysis/bloch_trajectory.py index 098f7e9a70..cba1474834 100644 --- a/qiskit_experiments/curve_analysis/standard_analysis/bloch_trajectory.py +++ b/qiskit_experiments/curve_analysis/standard_analysis/bloch_trajectory.py @@ -177,7 +177,7 @@ def _generate_fit_guesses( # oscillation amplitude might be almost zero, # then exclude from average because of lower SNR continue - fft_freq = curve.guess.frequency(data.x, data.y) + fft_freq = curve.guess.frequency_lorentz_fit(data.x, data.y) omega_xyz.append(fft_freq) if omega_xyz: omega = 2 * np.pi * np.average(omega_xyz) diff --git a/qiskit_experiments/curve_analysis/standard_analysis/decay.py b/qiskit_experiments/curve_analysis/standard_analysis/decay.py index e5ab271cb1..fb96479f60 100644 --- a/qiskit_experiments/curve_analysis/standard_analysis/decay.py +++ b/qiskit_experiments/curve_analysis/standard_analysis/decay.py @@ -14,7 +14,6 @@ from typing import List, Union, Optional import lmfit -import numpy as np import qiskit_experiments.curve_analysis as curve diff --git a/qiskit_experiments/curve_analysis/standard_analysis/oscillation.py b/qiskit_experiments/curve_analysis/standard_analysis/oscillation.py index c911698c88..28ff68ebdb 100644 --- a/qiskit_experiments/curve_analysis/standard_analysis/oscillation.py +++ b/qiskit_experiments/curve_analysis/standard_analysis/oscillation.py @@ -36,20 +36,21 @@ class OscillationAnalysis(curve.CurveAnalysis): defpar \rm amp: desc: Amplitude of the oscillation. init_guess: Calculated by :func:`~qiskit_experiments.curve_analysis.\ - guess.sinusoidal_freq_offset_amp`. - bounds: [-2 * maximum Y, 2 * maximum Y]. + guess.sinusoidal_freq_offset_amp`. + bounds: [0, 2 * maximum Y]. defpar \rm base: desc: Base line. init_guess: Calculated by :func:`~qiskit_experiments.curve_analysis.\ - guess.sinusoidal_freq_offset_amp`. + guess.sinusoidal_freq_offset_amp`. bounds: [-maximum Y, maximum Y]. defpar \rm freq: desc: Frequency of the oscillation. This is the fit parameter of interest. init_guess: Calculated by :func:`~qiskit_experiments.curve_analysis.\ - guess.sinusoidal_freq_offset_amp`. - bounds: [0, inf]. + guess.sinusoidal_freq_offset_amp`. + bounds: [0, Nyquist frequency] where Nyquist frequency is + a half of maximum sampling frequency of the X values. defpar \rm phase: desc: Phase of the oscillation. @@ -86,26 +87,27 @@ def _generate_fit_guesses( List of fit options that are passed to the fitter function. """ max_abs_y, _ = curve.guess.max_height(curve_data.y, absolute=True) + nyquist_freq = 1 / np.min(np.diff(curve_data.x)) / 2 user_opt.bounds.set_if_empty( - amp=(-2 * max_abs_y, 2 * max_abs_y), - freq=(0, np.inf), + amp=(0, 2 * max_abs_y), + freq=(0, nyquist_freq), phase=(-np.pi, np.pi), base=(-max_abs_y, max_abs_y), ) options = [] - for delay_ratio in (0.1, 0.2, 0.3): - for phase_guess in np.linspace(-np.pi, np.pi, 5): - tmp_opt = user_opt.copy() - freq, base, amp = curve.guess.sinusoidal_freq_offset_amp( - x=curve_data.x, - y=curve_data.y, - delay=int(delay_ratio * len(curve_data.x)), - ) - tmp_opt.p0.set_if_empty(freq=freq, base=base, amp=amp, phase=phase_guess) - options.append(tmp_opt) - + for amp, freq, base, phase in curve.guess.composite_sinusoidal_estimate( + x=curve_data.x, + y=curve_data.y, + **user_opt.p0, + ): + tmp_opt = user_opt.copy() + tmp_opt.p0.set_if_empty(amp=amp, freq=freq, base=base, phase=phase) + options.append(tmp_opt) + + if len(options) == 0: + return user_opt return options def _evaluate_quality(self, fit_data: curve.CurveFitResult) -> Union[str, None]: @@ -141,32 +143,33 @@ class DampedOscillationAnalysis(curve.CurveAnalysis): .. math:: F(x) = {\rm amp} \cdot e^{-x/\tau} - \cos(2\pi \cdot {\rm freq} \cdot t + \phi) + {\rm base} + \cos(2\pi \cdot {\rm freq} \cdot t + {\rm phase}) + {\rm base} # section: fit_parameters defpar \rm amp: desc: Amplitude of the oscillation. init_guess: Calculated by :func:`~qiskit_experiments.curve_analysis.\ - guess.sinusoidal_freq_offset_amp`. - bounds: [-2 * maximum Y, 2 * maximum Y]. + guess.sinusoidal_freq_offset_amp`. + bounds: [0, 2 * maximum Y]. defpar \rm base: desc: Base line. init_guess: Calculated by :func:`~qiskit_experiments.curve_analysis.\ - guess.sinusoidal_freq_offset_amp`. + guess.sinusoidal_freq_offset_amp`. bounds: [-maximum Y, maximum Y]. defpar \rm freq: desc: Frequency of the oscillation. This is the fit parameter of interest. init_guess: Calculated by :func:`~qiskit_experiments.curve_analysis.\ - guess.sinusoidal_freq_offset_amp`. - bounds: [0, inf]. + guess.sinusoidal_freq_offset_amp`. + bounds: [0, Nyquist frequency] where Nyquist frequency is + a half of maximum sampling frequency of the X values. defpar \tau: desc: Represents the rate of decay. init_guess: Determined by :py:func:`~qiskit_experiments.curve_analysis.\ - guess.exp_decay` + guess.exp_decay`. bounds: [0, inf] defpar \rm phase: @@ -182,7 +185,7 @@ def __init__( super().__init__( models=[ lmfit.models.ExpressionModel( - expr="amp * exp(-x / tau) * cos(2 * pi * freq * x + phi) + base", + expr="amp * exp(-x / tau) * cos(2 * pi * freq * x + phase) + base", name="cos_decay", ) ], @@ -204,10 +207,11 @@ def _generate_fit_guesses( List of fit options that are passed to the fitter function. """ max_abs_y, _ = curve.guess.max_height(curve_data.y, absolute=True) + nyquist_freq = 1 / np.min(np.diff(curve_data.x)) / 2 user_opt.bounds.set_if_empty( - amp=(-2 * max_abs_y, 2 * max_abs_y), - freq=(0, np.inf), + amp=(0, 2 * max_abs_y), + freq=(0, nyquist_freq), tau=(0, np.inf), phase=(-np.pi, np.pi), base=(-max_abs_y, max_abs_y), @@ -215,24 +219,19 @@ def _generate_fit_guesses( alpha = curve.guess.exp_decay(curve_data.x, curve_data.y) tau_guess = -1 / min(alpha, -1 / (100 * max(curve_data.x))) - options = [] - for delay_ratio in (0.1, 0.2, 0.3): - for phase_guess in np.linspace(0, np.pi, 5): - tmp_opt = user_opt.copy() - freq, base, amp = curve.guess.sinusoidal_freq_offset_amp( - x=curve_data.x, - y=curve_data.y, - delay=int(delay_ratio * len(curve_data.x)), - ) - tmp_opt.p0.set_if_empty( - freq=freq, - base=base, - amp=amp, - phase=phase_guess, - tau=tau_guess, - ) - options.append(tmp_opt) + pre_estimate = user_opt.p0.copy() + del pre_estimate["tau"] + options = [] + for amp, freq, base, phase in curve.guess.composite_sinusoidal_estimate( + x=curve_data.x, y=curve_data.y, **pre_estimate + ): + tmp_opt = user_opt.copy() + tmp_opt.p0.set_if_empty(amp=amp, freq=freq, base=base, phase=phase, tau=tau_guess) + options.append(tmp_opt) + + if len(options) == 0: + return user_opt return options def _evaluate_quality(self, fit_data: curve.CurveFitResult) -> Union[str, None]: diff --git a/qiskit_experiments/library/characterization/analysis/drag_analysis.py b/qiskit_experiments/library/characterization/analysis/drag_analysis.py index 35ad86537b..63f93840fb 100644 --- a/qiskit_experiments/library/characterization/analysis/drag_analysis.py +++ b/qiskit_experiments/library/characterization/analysis/drag_analysis.py @@ -111,40 +111,48 @@ def _generate_fit_guesses( # Use the highest-frequency curve to estimate the oscillation frequency. max_rep_model = self._models[-1] max_rep = max_rep_model.opts["data_sort_key"]["nrep"] - curve_data = curve_data.get_subset_of(max_rep_model._name) + max_rep_curve_data = curve_data.get_subset_of(max_rep_model._name) - x_data = curve_data.x - min_beta, max_beta = min(x_data), max(x_data) - - freqs_guess = curve.guess.frequency(curve_data.x, curve_data.y) / max_rep - user_opt.p0.set_if_empty(freq=freqs_guess) - - avg_x = (max(x_data) + min(x_data)) / 2 - span_x = max(x_data) - min(x_data) - beta_bound = max(5 / user_opt.p0["freq"], span_x) - - ptp_y = np.ptp(curve_data.y) + ptp_x = np.ptp(max_rep_curve_data.x) + ptp_y = np.ptp(max_rep_curve_data.y) + mid_x = np.min(max_rep_curve_data.x) + ptp_x / 2 + min_beta, max_beta = min(max_rep_curve_data.x), max(max_rep_curve_data.x) + nyquist_freq = 1 / np.min(np.diff(max_rep_curve_data.x)) / 2 user_opt.bounds.set_if_empty( amp=(-2 * ptp_y, 0), - freq=(0, np.inf), - beta=(avg_x - beta_bound, avg_x + beta_bound), - base=(min(curve_data.y) - ptp_y, max(curve_data.y) + ptp_y), + freq=(0, nyquist_freq), + base=(min(max_rep_curve_data.y) - ptp_y, max(max_rep_curve_data.y) + ptp_y), ) - base_guess = (max(curve_data.y) - min(curve_data.y)) / 2 - user_opt.p0.set_if_empty(base=(user_opt.p0["amp"] or base_guess)) - - # Drag curves can sometimes be very flat, i.e. averages of y-data - # and min-max do not always make good initial guesses. We therefore add - # 0.5 to the initial guesses. Note that we also set amp=-0.5 because the cosine function - # becomes +1 at zero phase, i.e. optimal beta, in which y data should become zero - # in discriminated measurement level. + + try: + user_freq = max_rep * user_opt.p0["freq"] + except TypeError: + # Not provided. + user_freq = None + + pre_estimate = { + "amp": user_opt.p0["amp"], + "freq": user_freq, + "base": user_opt.p0["base"], + "phase": 0.0, + } + options = [] - for amp_factor in (-1, -0.5, -0.25): + for amp, freq, base, _ in curve.guess.composite_sinusoidal_estimate( + x=max_rep_curve_data.x, y=max_rep_curve_data.y, **pre_estimate + ): for beta_guess in np.linspace(min_beta, max_beta, 20): - new_opt = user_opt.copy() - new_opt.p0.set_if_empty(amp=ptp_y * amp_factor, beta=beta_guess) - options.append(new_opt) + tmp_opt = user_opt.copy() + # amp (beta) is negative (positive) value to get P1 = 0 at optimal beta; x==beta. + tmp_opt.p0.set_if_empty( + amp=-1 * amp, freq=freq / max_rep, base=base, beta=beta_guess + ) + beta_bound = max(5 / tmp_opt.p0["freq"], ptp_x) + tmp_opt.bounds.set_if_empty(beta=(mid_x - beta_bound, mid_x + beta_bound)) + options.append(tmp_opt) + if len(options) == 0: + return user_opt return options def _run_curve_fit( diff --git a/qiskit_experiments/library/characterization/analysis/ramsey_xy_analysis.py b/qiskit_experiments/library/characterization/analysis/ramsey_xy_analysis.py index 5765b53f7c..327051e3ca 100644 --- a/qiskit_experiments/library/characterization/analysis/ramsey_xy_analysis.py +++ b/qiskit_experiments/library/characterization/analysis/ramsey_xy_analysis.py @@ -17,7 +17,10 @@ import lmfit import numpy as np +from qiskit.qobj.utils import MeasLevel + import qiskit_experiments.curve_analysis as curve +from qiskit_experiments.framework import Options, ExperimentData class RamseyXYAnalysis(curve.CurveAnalysis): @@ -36,27 +39,31 @@ class RamseyXYAnalysis(curve.CurveAnalysis): # section: fit_parameters defpar \rm amp: desc: Amplitude of both series. - init_guess: Half of the maximum y value less the minimum y value. When the - oscillation frequency is low, it uses an averaged difference of - Ramsey X data - Ramsey Y data. + init_guess: Calculated by :func:`~qiskit_experiments.curve_analysis.\ + guess.sinusoidal_freq_offset_amp`. When frequency is too low, + the difference of averages of Ramsey X and Y outcomes is set. bounds: [0, 2 * average y peak-to-peak] defpar \tau: desc: The exponential decay of the curve. - init_guess: The initial guess is obtained by fitting an exponential to the - square root of (X data)**2 + (Y data)**2. + init_guess: Determined by :py:func:`~qiskit_experiments.curve_analysis.\ + guess.exp_decay`. The square root of (X data)**2 + (Y data)**2 is used to + estimate the decay. When frequency is too low, large value is set. bounds: [0, inf] defpar \rm base: desc: Base line of both series. - init_guess: Roughly the average of the data. When the oscillation frequency is low, - it uses an averaged data of Ramsey Y experiment. + init_guess: Calculated by :func:`~qiskit_experiments.curve_analysis.\ + guess.sinusoidal_freq_offset_amp`. When frequency is too low, + the average of Ramsey Y outcomes is set. bounds: [min y - average y peak-to-peak, max y + average y peak-to-peak] defpar \rm freq: desc: Frequency of both series. This is the parameter of interest. - init_guess: The frequency with the highest power spectral density. - bounds: [-inf, inf] + init_guess: Calculated by :func:`~qiskit_experiments.curve_analysis.\ + guess.sinusoidal_freq_offset_amp`. When frequency is too low, zero is set. + bounds: [-Nyquist frequency, Nyquist frequency] where Nyquist frequency is + a half of maximum sampling frequency of the X values. defpar \rm phase: desc: Common phase offset. @@ -81,7 +88,7 @@ def __init__(self): ) @classmethod - def _default_options(cls): + def _default_options(cls) -> Options: """Return the default analysis options. See :meth:`~qiskit_experiment.curve_analysis.CurveAnalysis._default_options` for @@ -121,18 +128,23 @@ def _generate_fit_guesses( avg_y_ptp = 0.5 * (np.ptp(ramx_data.y) + np.ptp(ramy_data.y)) max_y = np.max(curve_data.y) min_y = np.min(curve_data.y) + nyquist_freq = 1 / np.min(np.diff(ramx_data.x)) / 2 user_opt.bounds.set_if_empty( amp=(0, full_y_ptp * 2), tau=(0, np.inf), base=(min_y - avg_y_ptp, max_y + avg_y_ptp), phase=(-np.pi, np.pi), + freq=(-nyquist_freq, nyquist_freq), ) + user_opt.p0.set_if_empty(phase=0.0) if avg_y_ptp < 0.5 * full_y_ptp: # When X and Y curve don't oscillate, X (Y) usually stays at P(1) = 1.0 (0.5). # So peak-to-peak of full data is something around P(1) = 0.75, while # single curve peak-to-peak is almost zero. + # Tau guess is not applicable because we cannot distinguish the decay curve + # from the sinusoidal oscillation. avg_x = np.average(ramx_data.y) avg_y = np.average(ramy_data.y) @@ -140,51 +152,48 @@ def _generate_fit_guesses( amp=np.abs(avg_x - avg_y), tau=100 * np.max(curve_data.x), base=avg_y, - phase=0.0, freq=0.0, ) return user_opt - base_guess_x = curve.guess.constant_sinusoidal_offset(ramx_data.y) - base_guess_y = curve.guess.constant_sinusoidal_offset(ramy_data.y) - base_guess = 0.5 * (base_guess_x + base_guess_y) - user_opt.p0.set_if_empty( - amp=0.5 * full_y_ptp, - base=base_guess, - phase=0.0, + alpha = curve.guess.exp_decay(ramx_data.x, np.sqrt(ramx_data.y**2 + ramy_data.y**2)) + tau_guess = -1 / min(alpha, -1 / (100 * max(curve_data.x))) + + # Generate guess iterators for Ramsey X and Y and average the estimate + pre_estimate = user_opt.p0.copy() + del pre_estimate["tau"] + + ramx_guess_iter = curve.guess.composite_sinusoidal_estimate( + x=ramx_data.x, + y=ramx_data.y, + **pre_estimate, + ) + ramy_guess_iter = curve.guess.composite_sinusoidal_estimate( + x=ramy_data.x, + y=ramy_data.y, + **pre_estimate, ) - # Guess the exponential decay by combining both curves - ramx_unbiased = ramx_data.y - user_opt.p0["base"] - ramy_unbiased = ramy_data.y - user_opt.p0["base"] - decay_data = ramx_unbiased**2 + ramy_unbiased**2 - if np.ptp(decay_data) < 0.95 * 0.5 * full_y_ptp: - # When decay is less than 95 % of peak-to-peak value, ignore decay and - # set large enough tau value compared with the measured x range. - user_opt.p0.set_if_empty(tau=1000 * np.max(curve_data.x)) - else: - user_opt.p0.set_if_empty(tau=-1 / curve.guess.exp_decay(ramx_data.x, decay_data)) - - # Guess the oscillation frequency, remove offset to eliminate DC peak - freq_guess_x = curve.guess.frequency(ramx_data.x, ramx_unbiased) - freq_guess_y = curve.guess.frequency(ramy_data.x, ramy_unbiased) - freq_val = 0.5 * (freq_guess_x + freq_guess_y) - - # FFT might be up to 1/2 bin off - df = 2 * np.pi / (np.min(np.diff(ramx_data.x)) * ramx_data.x.size) - freq_guesses = [freq_val - df, freq_val + df, freq_val] - - # Ramsey XY is frequency sign sensitive. - # Since experimental data is noisy, correct sign is hardly estimated with phase velocity. - # Try both positive and negative frequency to find the best fit. - opts = [] - for sign in (1, -1): - for freq_guess in freq_guesses: - opt = user_opt.copy() - opt.p0.set_if_empty(freq=sign * freq_guess) - opts.append(opt) - - return opts + options = [] + for ramx_guess, ramy_guess in zip(ramx_guess_iter, ramy_guess_iter): + amp_x, freq_x, base_x, _ = ramx_guess + amp_y, freq_y, base_y, _ = ramy_guess + for sign in (-1, 1): + # Ramsey XY is frequency sign sensitive. + # Since experimental data is noisy, correct sign is hardly estimated with + # phase velocity. Try both positive and negative frequency to find the best fit. + tmp_opt = user_opt.copy() + tmp_opt.p0.set_if_empty( + amp=0.5 * (amp_x + amp_y), + freq=sign * 0.5 * (freq_x + freq_y), + base=0.5 * (base_x + base_y), + tau=tau_guess, + ) + options.append(tmp_opt) + + if len(options) == 0: + return user_opt + return options def _evaluate_quality(self, fit_data: curve.CurveFitResult) -> Union[str, None]: """Algorithmic criteria for whether the fit is good or bad. @@ -204,3 +213,17 @@ def _evaluate_quality(self, fit_data: curve.CurveFitResult) -> Union[str, None]: return "good" return "bad" + + def _initialize( + self, + experiment_data: ExperimentData, + ): + super()._initialize(experiment_data) + + if experiment_data.metadata.get("meas_level", MeasLevel.CLASSIFIED) == MeasLevel.CLASSIFIED: + init_guess = self.options.get("p0", {}) + bounds = self.options.get("bounds", {}) + + init_guess["base"] = init_guess.get("base", 0.5) + bounds["base"] = bounds.get("base", (0.0, 1.0)) + self.set_options(p0=init_guess, bounds=bounds) diff --git a/qiskit_experiments/library/characterization/analysis/t2ramsey_analysis.py b/qiskit_experiments/library/characterization/analysis/t2ramsey_analysis.py index 3073e1793c..686636dcdc 100644 --- a/qiskit_experiments/library/characterization/analysis/t2ramsey_analysis.py +++ b/qiskit_experiments/library/characterization/analysis/t2ramsey_analysis.py @@ -13,8 +13,11 @@ T2Ramsey Experiment class. """ from typing import Union + +from qiskit.qobj.utils import MeasLevel + import qiskit_experiments.curve_analysis as curve -from qiskit_experiments.framework import Options +from qiskit_experiments.framework import Options, ExperimentData class T2RamseyAnalysis(curve.DampedOscillationAnalysis): @@ -66,3 +69,17 @@ def _evaluate_quality(self, fit_data: curve.CurveFitResult) -> Union[str, None]: return "good" return "bad" + + def _initialize( + self, + experiment_data: ExperimentData, + ): + super()._initialize(experiment_data) + + if experiment_data.metadata.get("meas_level", MeasLevel.CLASSIFIED) == MeasLevel.CLASSIFIED: + init_guess = self.options.get("p0", {}) + bounds = self.options.get("bounds", {}) + + init_guess["base"] = init_guess.get("base", 0.5) + bounds["base"] = bounds.get("base", (0.0, 1.0)) + self.set_options(p0=init_guess, bounds=bounds) diff --git a/qiskit_experiments/library/characterization/ramsey_xy.py b/qiskit_experiments/library/characterization/ramsey_xy.py index 0ede2f6fca..cadd32a4ad 100644 --- a/qiskit_experiments/library/characterization/ramsey_xy.py +++ b/qiskit_experiments/library/characterization/ramsey_xy.py @@ -18,7 +18,6 @@ from qiskit import QuantumCircuit from qiskit.circuit import Parameter from qiskit.providers.backend import Backend -from qiskit.qobj.utils import MeasLevel from qiskit_experiments.framework import BaseExperiment from qiskit_experiments.framework.restless_mixin import RestlessMixin @@ -227,17 +226,6 @@ def circuits(self) -> List[QuantumCircuit]: return circs - def _finalize(self): - # Set initial guess for sinusoidal offset when meas level is 2. - # This returns probability P1 thus offset=0.5 is obvious. - # This guarantees reasonable fit especially when data contains only less than half cycle. - meas_level = self.run_options.get("meas_level", MeasLevel.CLASSIFIED) - if meas_level == MeasLevel.CLASSIFIED: - init_guess = self.analysis.options.get("p0", {}) - if "base" not in init_guess: - init_guess["base"] = 0.5 - self.analysis.set_options(p0=init_guess) - def _metadata(self): metadata = super()._metadata() # Store measurement level and meas return if they have been diff --git a/releasenotes/notes/0_5_module_curve_analysis.yaml b/releasenotes/notes/0_5_module_curve_analysis.yaml new file mode 100644 index 0000000000..d238a7bae8 --- /dev/null +++ b/releasenotes/notes/0_5_module_curve_analysis.yaml @@ -0,0 +1,31 @@ +--- +features: + - | + Several new features of initial guess algorithms for the curve fitting have been introduced. + + * The one-stop initial guess algorithm for sinusoidal signal has been introduced as + :func:`~qiskit_experiments.curve_analysis.guess.composite_sinusoidal_estimate`. + This algorithm assumes a model characterized by the fit parameters (amp, freq, base, phase), + and internally assembles the best fitting algorithm depending on which + parameters are pre-determined. This function returns an iterator of multiple fit parameters. + * :func:`~qiskit_experiments.curve_analysis.guess.frequency_lorentz_fit` has been introduced. + This is a drop-in replacement of :func:`~qiskit_experiments.curve_analysis.guess.frequency`. + New function fine tunes the frequency estimate by Lorentzian fitting around the main FFT peak. + * :func:`~qiskit_experiments.curve_analysis.guess.low_frequency_limit` has been introduced. + This function is dedicated to the frequency estimate at the lower frequency limit. + Usually this function gives better frequency estimate around zero frequency. + * :func:`~qiskit_experiments.curve_analysis.guess.sinusoidal_freq_offset` has been introduced. + This algorithm uses two delayed data set to simulataneously estimate frequency and offset + of a sinusoidal signal. Because this function doesn't use FFT, it is more robust to + the error of base offset estimation error. + * :func:`~qiskit_experiments.curve_analysis.guess.exp_decay` has been upgraded to + new algorithm that allows analysis of a damped oscillation signal with finite base offset. + +deprecations: + - | + Several initial guess functions are merged and deprecated. + + * :func:`~qiskit_experiments.curve_analysis.guess.frequency` has been deprecated and + replaced with :func:`~qiskit_experiments.curve_analysis.guess.frequency_lorentz_fit`. + * :func:`~qiskit_experiments.curve_analysis.guess.oscillation_exp_decay` has been deprecated + and merged into :func:`~qiskit_experiments.curve_analysis.guess.exp_decay`. diff --git a/test/curve_analysis/test_guess.py b/test/curve_analysis/test_guess.py index e5255d9288..071fc88d38 100644 --- a/test/curve_analysis/test_guess.py +++ b/test/curve_analysis/test_guess.py @@ -36,7 +36,18 @@ def test_frequency(self, freq: float): x = np.linspace(-1, 1, 101) y = 0.3 * np.cos(2 * np.pi * freq * x + 0.5) + 1.2 - freq_guess = guess.frequency(x, y) + with self.assertWarns(DeprecationWarning): + freq_guess = guess.frequency(x, y) + + self.assertAlmostEqualAbsolute(freq_guess, np.abs(freq)) + + @data(1.1, 2.0, 1.6, -1.4, 4.5) + def test_lfit_frequency(self, freq: float): + """Test for frequency guess with Lorentzian fit.""" + x = np.linspace(-1, 1, 101) + y = 0.3 * np.cos(2 * np.pi * freq * x + 0.5) + + freq_guess = guess.frequency_lorentz_fit(x, y) self.assertAlmostEqualAbsolute(freq_guess, np.abs(freq)) @@ -46,7 +57,28 @@ def test_frequency_with_non_uniform_sampling(self, freq: float): x = np.concatenate((np.linspace(-1, 0, 15), np.linspace(0.1, 1, 30))) y = 0.3 * np.cos(2 * np.pi * freq * x + 0.5) + 1.2 - freq_guess = guess.frequency(x, y) + with self.assertWarns(DeprecationWarning): + freq_guess = guess.frequency(x, y) + + self.assertAlmostEqualAbsolute(freq_guess, np.abs(freq)) + + @data(1.1, 2.0, 1.6, -1.4, 4.5) + def test_lfit_frequency_with_non_uniform_sampling(self, freq: float): + """Test for frequency guess with Lorentzian fit with non uniform x value.""" + x = np.concatenate((np.linspace(-1, 0, 15), np.linspace(0.1, 1, 30))) + y = 0.3 * np.cos(2 * np.pi * freq * x + 0.5) + 1.2 + + freq_guess = guess.frequency_lorentz_fit(x, y) + + self.assertAlmostEqualAbsolute(freq_guess, np.abs(freq)) + + @data(0.2, 0.1, -0.2, 0.08, 0.6) + def test_low_frequency(self, freq: float): + """Test for frequency guess at low frequency limit.""" + x = np.linspace(-1, 1, 101) + y = 0.3 * np.cos(2 * np.pi * freq * x + 0.5) + + freq_guess = guess.low_frequency_limit(x, y) self.assertAlmostEqualAbsolute(freq_guess, np.abs(freq)) @@ -110,17 +142,33 @@ def test_exp_decay(self, alpha: float): self.assertAlmostEqualAbsolute(alpha_guess, alpha) + @data([1.2, 1.4], [-0.6, 2.5], [0.1, 2.3], [3.5, 1.1], [-4.1, 6.5], [3.0, 1.2]) + @unpack + def test_exp_decay_with_oscillation(self, alpha, freq): + """Test of exponential decay guess with oscillation.""" + x = np.linspace(0, 1, 100) + y = np.exp(alpha * x) * np.cos(2 * np.pi * freq * x) + + with self.assertWarns(DeprecationWarning): + alpha_guess = guess.oscillation_exp_decay(x, y) + + self.assertAlmostEqualAbsolute(alpha_guess, alpha) + @data( - [-1.2, 1.4, 0.0], - [-1.6, 2.5, 0.2], - [-1.3, 2.3, -0.3], - [-3.5, 1.1, 0.5], + [1.2, 1.4, 0.0], + [-0.6, 2.5, 0.2], + [0.1, 2.3, -0.3], + [3.5, 1.1, 0.5], [-4.1, 6.5, 2.4], - [-3.0, 1.2, -0.4], + [3.0, 1.2, -0.4], ) @unpack - def test_exp_decay_with_oscillation(self, alpha, freq, base): - """Test of exponential decay guess with oscillation.""" + def test_exp_decay_with_oscillation_unified(self, alpha, freq, base): + """Test of exponential decay guess with oscillation. + + Damped oscillation can be fit with exp_dacay with new implementation. + It is also capable of fitting a biased decay curve. + """ x = np.linspace(0, 1, 100) y = np.exp(alpha * x) * np.cos(2 * np.pi * freq * x) + base @@ -213,12 +261,52 @@ def test_rb_decay(self, a, b, alpha): [0.4, -0.1, 1.6, -0.3, -1.3], ) @unpack - def test_sinusoidal_freq_offset_amp(self, amp, alpha, freq, phase, base): - """Test simulataneous freq, offset and amp guess.""" + def test_sinusoidal_freq_offset(self, amp, alpha, freq, phase, base): + """Test simulataneous freq and offset guess.""" x = np.linspace(0, 1, 100) y = amp * np.exp(alpha * x) * np.cos(2 * np.pi * freq * x + phase) + base - freq_guess, base_guess, amp_guess = guess.sinusoidal_freq_offset(x, y, 15) + freq_guess, base_guess = guess.sinusoidal_freq_offset(x, y, 15) self.assertAlmostEqualAbsolute(freq, freq_guess) self.assertAlmostEqualAbsolute(base, base_guess) - self.assertAlmostEqualAbsolute(amp, amp_guess) + + def test_composite_sinusoidal(self): + """Test composite guess function of sinusoidal curve. + + Because this just internally calls other guess functions depending on the situation, + the generated guess values are not explicitly validated. + """ + x = np.linspace(0, 1, 100) + y = 0.1 * np.sin(2 * np.pi * 0.8 * x + 0.3) + 0.3 + + # test for case when all values are user-provided. + # this returns a single guess. + init_params = list( + guess.composite_sinusoidal_estimate(x, y, amp=0.1, freq=0.8, base=0.3, phase=0.3) + ) + self.assertEqual(len(init_params), 1) + + # test for case when frequency and base are not available. + # this generates three guesses for frequency with different delay parameters. + init_params = list(guess.composite_sinusoidal_estimate(x, y, amp=0.1, phase=0.3)) + self.assertEqual(len(init_params), 3) + + # test for case when base is provided. + # this generates a single guesses with FFT. + init_params = list(guess.composite_sinusoidal_estimate(x, y, amp=0.1, phase=0.3, base=0.3)) + self.assertEqual(len(init_params), 1) + + # test for case when frequency is provided and estimate period is longer than x range. + # this generates multiple guesses with different base. + init_params = list(guess.composite_sinusoidal_estimate(x, y, amp=0.1, freq=0.8, phase=0.3)) + self.assertEqual(len(init_params), 5) + + # test for case when frequency is provided and estimate period is shorter than x range. + # this generates a single guess for base by taking mid value. + init_params = list(guess.composite_sinusoidal_estimate(x, y, amp=0.1, freq=1.1, phase=0.3)) + self.assertEqual(len(init_params), 1) + + # test for case when phase is not known. + # this generates multiple guesses for phase. + init_params = list(guess.composite_sinusoidal_estimate(x, y, amp=0.1, freq=0.8, base=0.3)) + self.assertEqual(len(init_params), 5) diff --git a/test/library/characterization/test_t2ramsey.py b/test/library/characterization/test_t2ramsey.py index a1cdbacaef..eda2735d01 100644 --- a/test/library/characterization/test_t2ramsey.py +++ b/test/library/characterization/test_t2ramsey.py @@ -53,7 +53,7 @@ def test_t2ramsey_run_end2end(self): "amp": 0.5, "tau": estimated_t2ramsey[0], "freq": estimated_freq, - "phi": 0, + "phase": 0, "base": 0.5, } @@ -101,26 +101,6 @@ def test_t2ramsey_parallel(self): ) par_exp = ParallelExperiment([exp0, exp2]) - - exp0_p0 = { - "A": 0.5, - "T2star": t2ramsey[0], - "f": estimated_freq[0], - "phi": 0, - "B": 0.5, - } - - exp2_p0 = { - "A": 0.5, - "T2star": t2ramsey[1], - "f": estimated_freq[1], - "phi": 0, - "B": 0.5, - } - - exp0.analysis.set_options(p0=exp0_p0) - exp2.analysis.set_options(p0=exp2_p0) - expdata = par_exp.run(backend=backend, shots=2000, seed_simulator=1).block_for_results() self.assertExperimentDone(expdata) @@ -157,14 +137,6 @@ def test_t2ramsey_concat_2_experiments(self): backend = NoisyDelayAerBackend(t1=[2 * estimated_t2ramsey], t2=[estimated_t2ramsey]) exp0 = T2Ramsey(qubit, delays0, osc_freq=osc_freq) - default_p0 = { - "A": 0.5, - "T2star": estimated_t2ramsey, - "f": estimated_freq, - "phi": 0, - "B": 0.5, - } - exp0.analysis.set_options(p0=default_p0) # run circuits expdata0 = exp0.run(backend=backend, shots=1000, seed_simulator=1) @@ -174,7 +146,6 @@ def test_t2ramsey_concat_2_experiments(self): # second experiment delays1 = list(range(2, 65, 2)) exp1 = T2Ramsey(qubit, delays1, osc_freq=osc_freq) - exp1.analysis.set_options(p0=default_p0) expdata1 = exp1.run(backend=backend, analysis=None, shots=1000, seed_simulator=1) self.assertExperimentDone(expdata1)