Add GaussianSquare pulse shape

[Jacfomg]

This should add the flat top gaussian pulses or gaussian square. I would advise to test the devices output on the scope just to be sure the drivers behave properly for readout and flux pulses.

Checklist:

• [ ] Reviewers confirm new code works as expected.
• [ ] Tests are passing.
• [ ] Coverage does not decrease.
• [ ] Documentation is updated.

[codecov[bot]]

Codecov Report

Attention: 2 lines in your changes are missing coverage. Please review.

Comparison is base (05d0dc5) 62.18% compared to head (a735683) 62.43%.

Files	Patch %	Lines
src/qibolab/instruments/zhinst.py	0.00%	1 Missing :warning:
src/qibolab/pulses.py	97.61%	1 Missing :warning:

Additional details and impacted files

```diff @@ Coverage Diff @@ ## main #753 +/- ## ========================================== + Coverage 62.18% 62.43% +0.24% ========================================== Files 47 47 Lines 5837 5880 +43 ========================================== + Hits 3630 3671 +41 - Misses 2207 2209 +2 ``` | [Flag](https://app.codecov.io/gh/qiboteam/qibolab/pull/753/flags?src=pr&el=flags&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=qiboteam) | Coverage Δ | | |---|---|---| | [unittests](https://app.codecov.io/gh/qiboteam/qibolab/pull/753/flags?src=pr&el=flag&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=qiboteam) | `62.43% <95.34%> (+0.24%)` | :arrow_up: | Flags with carried forward coverage won't be shown. [Click here](https://docs.codecov.io/docs/carryforward-flags?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=qiboteam#carryforward-flags-in-the-pull-request-comment) to find out more.

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[alecandido]

To point it out now (later would be too late): I understand you're giving it the name used in Zurich, and it makes sense, but we have already another name around (in Qibosoq), i.e. FlatTop.

I'm not sure which is the most used name in the literature, but it'd be better to take a decision now, and stick to that name (possibly asking @rodolfocarobene to switch name even in Qibosoq, since it's still in the Qibo domain).

[rodolfocarobene]

I used the name "flattop" as it was used in qick, but I have no problem changing it in something else. I agree that we should all use the same one (if the definition is the same)

[Jacfomg]

To point it out now (later would be too late): I understand you're giving it the name used in Zurich, and it makes sense, but we have already another name around (in Qibosoq), i.e. FlatTop.

I'm not sure which is the most used name in the literature, but it'd be better to take a decision now, and stick to that name (possibly asking @rodolfocarobene to switch name even in Qibosoq, since it's still in the Qibo domain).

Yeah, I would not have any problem changing it either, I just went with the name I used when I implemented it a while ago, GaussianSquare, on the zurich driver. But, I would say that from FlatTop you could not infer that has a gaussian rise and fall as I suppuse you could also play with other functions there, as well as being the rectangular pulse also a flattop.

I think, the most used name is flattop gaussian in the literature, but it quite a long one.

[Jacfomg]

I used the name "flattop" as it was used in qick, but I have no problem changing it in something else. I agree that we should all use the same one (if the definition is the same)

For the definition I had two options, to use a width that tells the percentage of the pulse that is flat or just ask for the gaussian duration and flattop duration directly. I would prefer the first one as it could be simpler for duration sweeps. Which one did you use ?

[alecandido]

Btw: instead of glueing (i.e. concatenate) arrays, I'd suggest to use indexing, in the spirit of vectorized functions.

You could do something like:

def fvec(x):
    x = np.array(x)
    a = np.zeros_like(x)
    a[x > 0] = x[x > 0] + 100
    a[x < -1] = -x[x < -1]
    return a

The advantage is that this function could be called with an arbitrarily shaped array (including scalars), and give back a homogeneous result. Of course, we will call the function with an array of times, and at the moment you will have to supply also that (but in a later stage I'm planning to drop num_samples and just provide the callable as a function of time, the array of times can be generated in a single place for all the waveforms).

[alecandido]

@Jacfomg you could solve the vectorization and the size problem with a single (arguably simple) approach:

t = np.arange(num_samples)
length = num_samples
pulse = np.ones_like(t)
pulse[t < gaussian_samples] = ...
pulse[t > length - gaussian_samples] = ...

[Jacfomg]

@Jacfomg you could solve the vectorization and the size problem with a single (arguably simple) approach:

t = np.arange(num_samples)
length = num_samples
pulse = np.ones_like(t)
pulse[t < gaussian_samples] = ...
pulse[t > length - gaussian_samples] = ...

Yeah, I realized I will do this now, I started with the other comment as it was simpler before noticing this :)

[alecandido]

Ok, I'm definitely using this one as a complicate prototype of what I want to do to all the others, but I solved the problem of the elementwise application:;

def gaussian(t):
    return np.exp(...)

def fvec(t, gaussian_samples, length=None):
    if length is None:
        length = t.shape[0]

    pulse = 50 * np.ones_like(t)
    rise = t < gaussian_samples
    fall = t > length - gaussian_samples
    pulse[rise] = gaussian(t[rise])
    pulse[fall] = - gaussian(t[fall])
    return pulse

I.e.: if you want to do it, you will make my life slightly easier when I will attempt to replace Waveform and PulseShape with the vectorized callables. Otherwise, I will do it later anyhow.

[Jacfomg]

Ok, I'm definitely using this one as a complicate prototype of what I want to do to all the others, but I solved the problem of the elementwise application:;

def gaussian(t):
    return np.exp(...)

def fvec(t, gaussian_samples, length=None):
    if length is None:
        length = t.shape[0]

    pulse = 50 * np.ones_like(t)
    rise = t < gaussian_samples
    fall = t > length - gaussian_samples
    pulse[rise] = gaussian(t[rise])
    pulse[fall] = - gaussian(t[fall])
    return pulse

I.e.: if you want to do it, you will make my life slightly easier when I will attempt to replace Waveform and PulseShape with the vectorized callables. Otherwise, I will do it anyhow later.

I'm fine with it. that way I know a bit how things will work afterwards. I must say this seems to be closer to what zurich does for their waveforms as well, maybe you could have a look there.

[alecandido]

I'm fine with it. that way I know a bit how things will work afterwards. I must say this seems to be closer to what zurich does for their waveforms as well, maybe you could have a look there.

Yeah, I could guess it (and I will have a look). My feeling is that, despite all the changes in laboneq, Zurich has one of the most advanced interfaces (in a different direction from QM, in a certain sense). And Qibolab is really at the same level of laboneq, not a higher one. So, we're duplicating part of their work, for the sake of standardization with the other instruments.

[Jacfomg]

Ok, I'm definitely using this one as a complicate prototype of what I want to do to all the others, but I solved the problem of the elementwise application:;

def gaussian(t):
    return np.exp(...)

def fvec(t, gaussian_samples, length=None):
    if length is None:
        length = t.shape[0]

    pulse = 50 * np.ones_like(t)
    rise = t < gaussian_samples
    fall = t > length - gaussian_samples
    pulse[rise] = gaussian(t[rise])
    pulse[fall] = - gaussian(t[fall])
    return pulse

I.e.: if you want to do it, you will make my life slightly easier when I will attempt to replace Waveform and PulseShape with the vectorized callables. Otherwise, I will do it later anyhow.

That should be it, I will leave the functions inside the GaussianSquare for now then.

def gaussian(t, rel_sigma, gaussian_samples):
  mu = (2 * gaussian_samples - 1) / 2
  sigma = (2 * gaussian_samples) / rel_sigma
  return np.exp(-0.5 * ((t - mu) / sigma) ** 2)

def fvec(t, gaussian_samples, rel_sigma, length=None):
  if length is None:
      length = t.shape[0]

  pulse = np.ones_like(t, dtype=float)
  rise = t < gaussian_samples
  fall = t > length - gaussian_samples -1 
  pulse[rise] = gaussian(t[rise], rel_sigma, gaussian_samples)
  pulse[fall] = gaussian(t[rise], rel_sigma, gaussian_samples)[::-1]
  return pulse

num_samples = int(np.rint(self.pulse.duration * sampling_rate))
gaussian_samples = num_samples * (1 -self.width) // 2
t = np.arange(0, num_samples)

pulse = fvec(t, gaussian_samples, rel_sigma=self.rel_sigma)

[rodolfocarobene]

For the definition I had two options, to use a width that tells the percentage of the pulse that is flat or just ask for the gaussian duration and flattop duration directly. I would prefer the first one as it could be simpler for duration sweeps. Which one did you use ?

I'm using the QICK definition, so you give the duration of the pulse and the sigma of the gaussian... To be fair I don't fully understand how this works

[Jacfomg]

For the definition I had two options, to use a width that tells the percentage of the pulse that is flat or just ask for the gaussian duration and flattop duration directly. I would prefer the first one as it could be simpler for duration sweeps. Which one did you use ?

I'm using the QICK definition, so you give the duration of the pulse and the sigma of the gaussian... To be fair I don't fully understand how this works

Me neither you would need at least 3 parameters to describe the pulse to my understanding. Pulse duration, gaussian duration or width and the sigma.

[alecandido]

That should be it, I will leave the functions inside the GaussianSquare for now then.

Indeed, since the current interface requires the envelopes methods, and return Waveforms, this is definitely what you should do.

Me neither you would need at least 3 parameters to describe the pulse to my understanding. Pulse duration, gaussian duration or width and the sigma.

Agreed, it can't be fewer than that.

I'm using the QICK definition, so you give the duration of the pulse and the sigma of the gaussian... To be fair I don't fully understand how this works

To me, it seems that you're only making a Gaussian pulse, rather than an actual FlatTop... https://github.com/qiboteam/qibosoq/blob/0dcd37c08f173f28cf7a83e662bb3300fa26f7d8/src/qibosoq/programs/base.py#L135-L138 https://github.com/openquantumhardware/qick/blob/c1e3c0dba8bbf5e48c79ae1101f3f186c2a83eec/qick_lib/qick/asm_v2.py#L499-L501

It looks like QICK wants to generate the Gaussian, and later split it in half (not sure when or where, since they are filling an internal dictionary, not supposed to be touched outside https://github.com/openquantumhardware/qick/blob/c1e3c0dba8bbf5e48c79ae1101f3f186c2a83eec/qick_lib/qick/qick_asm.py#L785)

[rodolfocarobene]

To me, it seems that you're only making a Gaussian pulse, rather than an actual FlatTop...

It is actually correct to just load a Gaussian waveform, since the constant part does not require it. The flat top is then "created" from the style parameter.

I guess that they decide the ratio in duration between gaussian and constant part from the sigma. But I was not able to find exactly where ...

[rodolfocarobene]

Or maybe my "length=length/2" is actually enforcing that the gaussian part has an equal length to the constant one... I will have to test it and fix it, but if that was the case it should be very easy