Fixing fft guess #915

andrea-pasquale · 2024-07-03T08:30:40Z

Closes #897.

Checklist:

Reviewers confirm new code works as expected.
Tests are passing.
Coverage does not decrease.
Documentation is updated.
Compatibility with Qibo modules (Please edit this section if the current pull request is not compatible with the following branches).
- Qibo: master
- Qibolab: main
- Qibolab_platforms_qrc: main

codecov · 2024-07-03T08:43:39Z

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 97.46%. Comparing base (b0dabce) to head (a5141fd).
Report is 3 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #915      +/-   ##
==========================================
+ Coverage   97.41%   97.46%   +0.04%     
==========================================
  Files         117      117              
  Lines        8905     8874      -31     
==========================================
- Hits         8675     8649      -26     
+ Misses        230      225       -5

Flag	Coverage Δ
unittests	`97.46% <100.00%> (+0.04%)`	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files	Coverage Δ
src/qibocal/protocols/drag.py	`96.33% <100.00%> (-0.17%)`	⬇️
src/qibocal/protocols/flipping.py	`97.89% <100.00%> (-0.17%)`	⬇️
src/qibocal/protocols/flipping_signal.py	`98.50% <100.00%> (-0.09%)`	⬇️
src/qibocal/protocols/rabi/amplitude.py	`96.87% <100.00%> (ø)`
src/qibocal/protocols/rabi/amplitude_frequency.py	`97.91% <100.00%> (-0.03%)`	⬇️
...bocal/protocols/rabi/amplitude_frequency_signal.py	`98.36% <100.00%> (-0.02%)`	⬇️
src/qibocal/protocols/rabi/amplitude_signal.py	`97.46% <100.00%> (+0.03%)`	⬆️
src/qibocal/protocols/rabi/length.py	`100.00% <100.00%> (+2.70%)`	⬆️
src/qibocal/protocols/rabi/length_frequency.py	`97.93% <100.00%> (ø)`
.../qibocal/protocols/rabi/length_frequency_signal.py	`98.37% <100.00%> (ø)`
... and 7 more

... and 1 file with indirect coverage changes

rodolfocarobene · 2024-07-03T12:53:16Z

Sorry, this PR closes which issue? #651 seems just an old PR.
Anyway, it looks good to me, but I don't understand exactly what it is doing, was the fft not working before?

andrea-pasquale · 2024-07-03T12:55:38Z

Sorry, this PR closes which issue? #651 seems just an old PR. Anyway, it looks good to me, but I don't understand exactly what it is doing, was the fft not working before?

Thanks @rodolfocarobene I linked the wrong issue.
It was working before, now it seems that during some refactor something went wrong...
Now it should work properly.

rodolfocarobene

Apart to the line I don't really understand, it looks good to me

rodolfocarobene · 2024-07-03T13:21:09Z

src/qibocal/protocols/rabi/utils.py

+    fft = np.fft.rfft(y)
+    fft_freqs = np.fft.rfftfreq(len(y), d=(x[1] - x[0]))
+    mags = abs(fft)
+    mags[0] = 0


Could you explain this?

Yes, when you run fft you expect to see a peak corresponding to the frequency of your signal.
However, you will always find a peak in the first element of the array which corresponds to the dc component, which can be linked to the sample average.
main was currently broken because we were not removing this first point (in this case I'm just setting in to zero).

I see, thanks for the explanation!

Thanks for reviewing!
I requested you for a review since you recently edited the Rabi protocols #857. I just wanted someone to double check that I didn't break anything :)

Why not fft[1:], instead of setting it to 0?

Ok, it's stupid, but it's just more readable (for me):

mags = abs(fft[1:])

(obviously cutting also fft_freqs)

It's not truly changing anything, so it definitely personal preference...

Edoardo-Pedicillo

Thank @andrea-pasquale ! You should propagate it also in other protocols like Ramsey

Edoardo-Pedicillo · 2024-07-04T06:18:21Z

src/qibocal/protocols/rabi/utils.py

    # 0.5 hardcoded guess for less than one oscillation
-    return x[index] / (x[1] - x[0]) if index is not None else 0.5
+    return 1 / dominant_freq if dominant_freq != 0 else 0.5


I think it is better to return 1 in stead of 0.5 when you acquire less than a period

0.5 is the Nyquist frequency, so it's the maximum value returned by np.fft.rfftfreq

First of all, it is returning a period so it should be 2. I am not fully familiar with this concepts, but the numpy function is returning that maximum value because it is expecting you chose the parameters wisely to completely reconstruct the function but in general this is not the case. Of course, if the frequency is zero and we are expecting an oscillation we are out of the Nyquist theorem, so we can return whatever value (or raise an error).

To be fair: this is fully redundant, because we could deal with this in the above conditional, instead of reflecting it here.
I.e. instead of

qibocal/src/qibocal/protocols/rabi/utils.py

Lines 281 to 286 in 26a43b5

if len(local_maxima) > 0:

dominant_freq = fft_freqs[np.argmax(mags)]

else:

dominant_freq = 0 # Default if no peaks are found

# 0.5 hardcoded guess for less than one oscillation

return 1 / dominant_freq if dominant_freq != 0 else 0.5

use:

if len(local_maxima) > 0: return 1 / fft_freqs[np.argmax(mags)] return 2.

Even returning that value when no maximum is found is quite arbitrary. In principle, we could just fail. But in this way you leave your optimizer a chance of finding on its own (in principle, we could even return 4., instead, as it would be at half the frequency range [0, 0.5], but it's just a random guess)

First of all, it is returning a period so it should be 2.

Correct.

I am not fully familiar with this concepts, but the numpy function is returning that maximum value

The frequency is half the sampling frequency, because you want to make sure to have more than one point per period (approximate explanation). If you take only one point per period, exactly spaced by one period, it would appear a constant. If your frequency is slightly less, you might see as a much wider oscillation (because if you move by 2pi - delta, it would appear as if you were moving by just -delta, think about when you look at the fan or helicopter blades, and they seem to move slowly in the opposite direction).
The Nyquist-Shannon theorem proves formally that this "stops" when you take at least two points per period, because you realize there has to be an oscillations in that time range.

So, yes, the period saturates to 2.

it is expecting you chose the parameters wisely to completely reconstruct the function but in general this is not the case.

The problem is that, above the Nyquist frequency, you simply don't know the frequency, if you're sampling.

You can be use a response function to modulate your signal, and move to a different Nyquist zone. But this is inductive bias about you knowing already something about your frequency, so you can distinguish aliases (@rodolfocarobene has also a nice plot to show how this work when you used during synthesis)

Even returning that value when no maximum is found is quite arbitrary. In principle, we could just fail. But in this way you leave your optimizer a chance of finding on its own (in principle, we could even return 4., instead, as it would be at half the frequency range [0, 0.5], but it's just a random guess)

Ok, but I don't think we can trust the parameters of the fit (otherwise this is against the Nyquist theorem).

The problem is that, above the Nyquist frequency, you simply don't know the frequency, if you're sampling.

Yes, so if the dominant frequency is zero, this means that your real frequency is above the Nyquist one, but you cannot easily infer which one, so 2 is as good as 3,4,5. For this reason I will suggest to raise an error.

I was not taking into account the noisy case, usually (especially with probabilities) it is not so common, but also the case above the Nyquist frequency is not, so it is better to take it into account and return a guess. Maybe I suggest to have a warning just to be aware of the situation.

As usual, I'm not a huge fun of warnings: they are ignored in programs (since side-effects), and then ignored also by users.

In case, much better to return None, and leave the choice to the caller. If this information should be called, it can be further propagated, and then made part of the Results.

(i.e. always prefer returning over warnings - otherwise use exceptions)

Actually by putting None I believe scipy will be put 1 instead.
Shall we stick to 2 or 4 as an initial guess?

Well, there are two considerations here:

if you return 2 or 4, it is indistinguishable from an actual outcome, so you'll never know (out of this function) whether the maximum was found or not - moreover, this function's sake is to give the best estimate of the sine frequency, not to provide an initial value for a parameter in a fit (a function should never know how it is used, it's definitely anti-modular)

if you return None, you will have to provide a fallback in every place the function is used - it's not terrible, since the "optimal" fallback value may depend on the usage, but if it's happening a lot could be pretty redundant

A modular design could be something like:

fallback_period(guess_period())

with the outer function to be something like lambda x: x if x is not None else 4. If the function is not going to be used many times, you can just inline it.

Thanks I implemented the fallback function in b3891dc

tests/test_rabi.py

src/qibocal/protocols/flipping_signal.py

src/qibocal/protocols/ramsey/ramsey.py

alecandido

I didn't test, but it should be mostly fine.

This PR is already doing more than one thing, since it started by fixing the FFT guess, and then deduplicated this guess in many places. Which is great, but it could be another PR, before or after this.
I believe there are more things related to sinusoidal fits to be deduplicated, but let's keep them for another PR.

src/qibocal/protocols/utils.py

src/qibocal/protocols/rabi/amplitude_frequency.py

Co-authored-by: Alessandro Candido <[email protected]>

src/qibocal/protocols/utils.py

Co-authored-by: Edoardo Pedicillo <[email protected]>

andrea-pasquale added 2 commits July 3, 2024 12:29

fix: Fix fft guess

e80f7ec

fix: Add missing file

cccb07a

test: Add fft test

26a43b5

andrea-pasquale marked this pull request as ready for review July 3, 2024 08:51

andrea-pasquale requested review from rodolfocarobene and Edoardo-Pedicillo July 3, 2024 09:23

andrea-pasquale self-assigned this Jul 3, 2024

andrea-pasquale added the bug Something isn't working label Jul 3, 2024

rodolfocarobene approved these changes Jul 3, 2024

View reviewed changes

Edoardo-Pedicillo reviewed Jul 4, 2024

View reviewed changes

andrea-pasquale added 7 commits July 4, 2024 14:17

feat: Progaate fft to other protocols

c606f52

feat: Fix flipping fit

8ee66b2

test: Fix tests

a6a7ec1

fix: Fix flipping in tests

9d7b893

feat: Change Ramsey to excited state and change initiall guess

0e71fbb

test: Test period > 1

d8f4a7c

feat: Change default to None

ea7bb67

andrea-pasquale commented Jul 4, 2024

View reviewed changes

src/qibocal/protocols/flipping_signal.py Outdated Show resolved Hide resolved

src/qibocal/protocols/ramsey/ramsey.py Show resolved Hide resolved

andrea-pasquale added 4 commits July 5, 2024 11:10

feat: Show opposite sign for Ramsey

9b07c05

refactor: Remove useless 0.5 shift

9a4112a

feat: Implement fallback_period function

b3891dc

Merge branch 'main' into fix_fft

d54b696

andrea-pasquale requested review from Edoardo-Pedicillo and alecandido July 9, 2024 11:16

alecandido approved these changes Jul 9, 2024

View reviewed changes

src/qibocal/protocols/utils.py Outdated Show resolved Hide resolved

src/qibocal/protocols/rabi/amplitude_frequency.py Outdated Show resolved Hide resolved

andrea-pasquale and others added 2 commits July 9, 2024 19:22

refactor: Make pylint happy

b095feb

Co-authored-by: Alessandro Candido <[email protected]>

refactor: f -> period

b0def1b

Edoardo-Pedicillo approved these changes Jul 10, 2024

View reviewed changes

src/qibocal/protocols/utils.py Outdated Show resolved Hide resolved

refactor: Remove redundant variable

a5141fd

Co-authored-by: Edoardo Pedicillo <[email protected]>

andrea-pasquale added this pull request to the merge queue Jul 10, 2024

Merged via the queue into main with commit d90d6dd Jul 10, 2024
21 checks passed

andrea-pasquale deleted the fix_fft branch July 10, 2024 08:07

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fixing fft guess #915

Fixing fft guess #915

andrea-pasquale commented Jul 3, 2024 •

edited

Loading

codecov bot commented Jul 3, 2024 •

edited

Loading

rodolfocarobene commented Jul 3, 2024

andrea-pasquale commented Jul 3, 2024

rodolfocarobene left a comment

rodolfocarobene Jul 3, 2024

andrea-pasquale Jul 3, 2024

rodolfocarobene Jul 3, 2024

andrea-pasquale Jul 3, 2024

alecandido Jul 3, 2024 •

edited

Loading

Edoardo-Pedicillo left a comment

Edoardo-Pedicillo Jul 4, 2024

alecandido Jul 4, 2024 •

edited

Loading

Edoardo-Pedicillo Jul 4, 2024 •

edited

Loading

alecandido Jul 4, 2024

Edoardo-Pedicillo Jul 4, 2024 •

edited

Loading

Edoardo-Pedicillo Jul 4, 2024 •

edited

Loading

alecandido Jul 4, 2024

andrea-pasquale Jul 4, 2024

alecandido Jul 5, 2024

andrea-pasquale Jul 9, 2024

alecandido left a comment

	if len(local_maxima) > 0:
	dominant_freq = fft_freqs[np.argmax(mags)]
	else:
	dominant_freq = 0 # Default if no peaks are found
	# 0.5 hardcoded guess for less than one oscillation
	return 1 / dominant_freq if dominant_freq != 0 else 0.5

Fixing fft guess #915

Fixing fft guess #915

Conversation

andrea-pasquale commented Jul 3, 2024 • edited Loading

codecov bot commented Jul 3, 2024 • edited Loading

Codecov Report

rodolfocarobene commented Jul 3, 2024

andrea-pasquale commented Jul 3, 2024

rodolfocarobene left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

alecandido Jul 3, 2024 • edited Loading

Choose a reason for hiding this comment

Edoardo-Pedicillo left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

alecandido Jul 4, 2024 • edited Loading

Choose a reason for hiding this comment

Edoardo-Pedicillo Jul 4, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Edoardo-Pedicillo Jul 4, 2024 • edited Loading

Choose a reason for hiding this comment

Edoardo-Pedicillo Jul 4, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

alecandido left a comment

Choose a reason for hiding this comment

andrea-pasquale commented Jul 3, 2024 •

edited

Loading

codecov bot commented Jul 3, 2024 •

edited

Loading

alecandido Jul 3, 2024 •

edited

Loading

alecandido Jul 4, 2024 •

edited

Loading

Edoardo-Pedicillo Jul 4, 2024 •

edited

Loading

Edoardo-Pedicillo Jul 4, 2024 •

edited

Loading

Edoardo-Pedicillo Jul 4, 2024 •

edited

Loading