Repository: pasqal-io/Pulser (language: python)

Issue:
Various formatting issues
There are various formatting issues:

- the code in notebooks is not formatted by black
- imports are sorted in a very conservative way by black, which let them messy
- some files could be formatted systematically, like markdown ones

I propose adding the following requirements to the project:

- `black[jupyter]`
- `isort`

and format some other files.


Discussion / hints:

Fix patch:
diff --git a/.flake8 b/.flake8
index ce68d0eac..29fee19f4 100644
--- a/.flake8
+++ b/.flake8
@@ -1,6 +1,11 @@
 [flake8]
-exclude = ./build, ./docs
 docstring-convention = google
+exclude = ./build, ./docs
+extend-ignore =
+    # D105 Missing docstring in magic method
+    D105,
+    # E203 whitespace before ':' (for compliance with black)
+    E203,
 per-file-ignores =
   # D100 Missing docstring in public module
   # D103 Missing docstring in public function
@@ -8,8 +13,3 @@ per-file-ignores =
   pulser/tests/*: D100, D103
   __init__.py: F401
   setup.py: D100
-extend-ignore =
-    # D105 Missing docstring in magic method
-    D105,
-    # E203 whitespace before ':' (for compliance with black)
-    E203,
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index bd8bc8131..c886f4327 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -8,10 +8,10 @@ The steps to take will depend on what you want to do, but generally you'll want
 
 1. Do a quick search for keywords over the existing issues to ensure yours has not been added yet.
 2. If you can't find your issue already listed, create a new one. Please try to be as clear and detailed as possible in your description.
+
 - If you just want to give a suggestion or report a bug, that's already excellent and we thank you for it! Your issue will be listed and, hopefully, someone will take care of it at some point.
 - However, you may also want to be the one solving your issue, which would be even better! In these cases, you would proceed by preparing a [Pull Request](#making-a-pull-request).
 
-
 ## Making a Pull Request
 
 We're thrilled that you want to contribute to Pulser! For general contributions, we use a combination of two Git workflows: the [Forking workflow](https://www.atlassian.com/git/tutorials/comparing-workflows/forking-workflow) and the [Gitflow workflow](https://www.atlassian.com/git/tutorials/comparing-workflows/gitflow-workflow). If you don't know what any of this means, don't worry, you should still be able to make your contribution just by following the instructions detailed below. Nonetheless, in a nutshell, this workflow will have you making a fork from the main Pulser repository and working off a branch from `develop` (**not** `master`). Thus, you'll start your branch from `develop` and end with a pull request that merges your branch back to `develop`. The only exception to this rule is when making a `hotfix`, but in these cases the Pulser development team will take care of it for you.
@@ -19,12 +19,15 @@ We're thrilled that you want to contribute to Pulser! For general contributions,
 Here are the steps you should follow to make your contribution:
 
 0. Fork the Pulser repository and add the main Pulser repository as the `upstream`. You only have to do this once and you do so by clicking the "Fork" button at the upper right corner of the [repo page](https://github.com/pasqal-io/Pulser). This will create a new GitHub repo at `https://github.com/USERNAME/Pulser`, where `USERNAME` is your GitHub ID. Then, `cd` into the folder where you would like to place your new fork and clone it by doing:
+
     ```bash
     git clone https://github.com/USERNAME/Pulser.git
     ```
+
     **Note**: `USERNAME` should be replaced by your own GitHub ID.
 
    Then, you'll want to go into the directory of your brand new Pulser fork and add the main Pulser repository as the `upstream` by running:
+
    ```bash
    git remote add upstream https://github.com/pasqal-io/Pulser.git
    ```
@@ -32,10 +35,12 @@ Here are the steps you should follow to make your contribution:
 1. Have the related issue assigned to you. We suggest that you work only on issues that have been assigned to you; by doing this, you make sure to be the only one working on this and we prevent everyone from doing duplicate work. If a related issue does not exist yet, consult the [section above](#reporting-a-bug-or-suggesting-a-feature) to see how to proceed.
 
 2. You'll want to create a new branch where you will do your changes. The starting point will be `upstream/develop`, which is where you'll ultimately merge your changes. Inside your fork's root folder, run:
+
     ```bash
     git fetch upstream
     git checkout -b branch-name-here upstream/develop
     ```
+
     This will create and checkout the new branch, where you will do your changes.
 
     **Note**: `branch-name-here` should be replaced by the name you'll give your branch. Try to be descriptive, pick a name that identifies your new feature.
@@ -43,12 +48,15 @@ Here are the steps you should follow to make your contribution:
 3. Do your work and commit the changes to this new branch. Try to make the first line of your commit messages short but informative; in case you want to go into more detail, you have the option to do so in the next lines.
 
 4. At this point, your branch might have drifted out of sync with Pulser's `develop` branch (the `upstream`). By running
+
     ```shell
     git pull upstream develop
     ```
+
    you will fetch the latest changes in `upstream/develop` and merge them with your working branch, at which point you'll have to solve any merge conflicts that may    arise. This will keep your working branch in sync with `upstream/develop`.
 
 5. Finally, you push your code to your local branch:
+
     ```bash
     git push origin branch-name-here
     ```
@@ -66,29 +74,37 @@ pip install -r requirements.txt
 ```
 
 - **Tests**: We use [`pytest`](https://docs.pytest.org/en/latest/) to run unit tests on our code. If your changes break existing tests, you'll have to update these tests accordingly. Additionally, we aim for 100% coverage over our code. Try to cover all the new lines of code with simple tests, which should be placed in the `Pulser/pulser/tests` folder. To run all tests and check coverage, run:
+
     ```bash
     pytest --cov pulser
     ```
+
     All lines that are not meant to be tested must be tagged with `# pragma: no cover`. Use it sparingly,
     every decision to leave a line uncovered must be well justified.
 
 - **Style**: We use [`flake8`](https://flake8.pycqa.org/en/latest/) and the `flake8-docstrings` extension to enforce PEP8 style guidelines. To lint your code with `flake8`, simply run:
+
     ```bash
     flake8 .
     ```
+
     To help you keep your code compliant with PEP8 guidelines effortlessly, we suggest you look into installing a linter for your text editor of choice.
 
 - **Format**: We use the [`black`](https://black.readthedocs.io/en/stable/index.html) auto-formatter to enforce a consistent style throughout the entire code base. It will also ensure your code is compliant with the formatting enforced by `flake8` for you. To automatically format your code with black, just run:
+
     ```bash
     black .
     ```
+
     Note that some IDE's and text editors support plug-ins which auto-format your code with `black` upon saving, so you don't have to worry about code format at all.
 
 - **Type hints**: We use [mypy](http://mypy-lang.org/) to type check the code. Your code should have type
 annotations and pass the type checks from running:
+
     ```bash
     mypy
     ```
+
     In case `mypy` produces a false positive, you can ignore the respective line by adding the `# type: ignore` annotation.
 
     **Note**: Type hints for `numpy` have only been added in version 1.20. Make sure you have `numpy >= 1.20`
diff --git a/README.md b/README.md
index ebd2dc48e..eb6cceed2 100644
--- a/README.md
+++ b/README.md
@@ -2,9 +2,9 @@
 
 Pulser is a framework for composing, simulating and executing **pulse** sequences for neutral-atom quantum devices.
 
-**Documentation** for the [latest release](https://pypi.org/project/pulser/) of `pulser` is available at https://pulser.readthedocs.io (for the docs tracking the `develop` branch of this repository, visit https://pulser.readthedocs.io/en/latest instead).
+**Documentation** for the [latest release](https://pypi.org/project/pulser/) of `pulser` is available at <https://pulser.readthedocs.io> (for the docs tracking the `develop` branch of this repository, visit <https://pulser.readthedocs.io/en/latest> instead).
 
-The source code can be found at https://github.com/pasqal-io/Pulser.
+The source code can be found at <https://github.com/pasqal-io/Pulser>.
 
 ## Overview of Pulser
 
@@ -51,6 +51,7 @@ Then, you can do the following to run the test suite and report test coverage:
 ```bash
 pytest --cov pulser
 ```
+
 ## Contributing
 
 Want to contribute to Pulser? Great! See [How to Contribute][contributing] for information on how you can do so.
diff --git a/docs/source/intro_rydberg_blockade.ipynb b/docs/source/intro_rydberg_blockade.ipynb
index 2e91174e5..ba7eb8459 100644
--- a/docs/source/intro_rydberg_blockade.ipynb
+++ b/docs/source/intro_rydberg_blockade.ipynb
@@ -88,8 +88,10 @@
     "from pulser import Pulse\n",
     "from pulser.waveforms import RampWaveform, BlackmanWaveform\n",
     "\n",
-    "duration = 1000 # Typical: ~1 µsec\n",
-    "pulse = Pulse(BlackmanWaveform(duration, np.pi), RampWaveform(duration, -5., 10.), 0)\n",
+    "duration = 1000  # Typical: ~1 µsec\n",
+    "pulse = Pulse(\n",
+    "    BlackmanWaveform(duration, np.pi), RampWaveform(duration, -5.0, 10.0), 0\n",
+    ")\n",
     "pulse.draw()"
    ]
   },
@@ -131,18 +133,18 @@
    "source": [
     "from pulser import Sequence\n",
     "\n",
-    "reg = Register.rectangle(1, 2, spacing=8, prefix='atom')\n",
+    "reg = Register.rectangle(1, 2, spacing=8, prefix=\"atom\")\n",
     "reg.draw()\n",
     "\n",
-    "pi_pulse = Pulse.ConstantDetuning(BlackmanWaveform(duration, np.pi), 0., 0.)\n",
+    "pi_pulse = Pulse.ConstantDetuning(BlackmanWaveform(duration, np.pi), 0.0, 0.0)\n",
     "\n",
     "seq = Sequence(reg, Chadoq2)\n",
     "\n",
-    "seq.declare_channel('ryd','rydberg_local','atom0')\n",
+    "seq.declare_channel(\"ryd\", \"rydberg_local\", \"atom0\")\n",
     "\n",
-    "seq.add(pi_pulse,'ryd')\n",
-    "seq.target('atom1', 'ryd')\n",
-    "seq.add(pi_pulse,'ryd')\n",
+    "seq.add(pi_pulse, \"ryd\")\n",
+    "seq.target(\"atom1\", \"ryd\")\n",
+    "seq.add(pi_pulse, \"ryd\")\n",
     "\n",
     "seq.draw()"
    ]
@@ -199,25 +201,27 @@
     "data = []\n",
     "distances = np.linspace(6.5, 14, 7)\n",
     "\n",
-    "r = [1,0] # |r>\n",
-    "rr = np.kron(r,r) # |rr>\n",
+    "r = [1, 0]  # |r>\n",
+    "rr = np.kron(r, r)  # |rr>\n",
     "occup = [np.outer(rr, np.conj(rr))]  # |rr><rr|\n",
     "\n",
-    "for i,R in enumerate(distances):\n",
+    "for i, R in enumerate(distances):\n",
     "    # Atom Register and Device\n",
-    "    reg = Register.rectangle(1,2,spacing=R, prefix='atom')\n",
+    "    reg = Register.rectangle(1, 2, spacing=R, prefix=\"atom\")\n",
     "\n",
     "    # Pulse Sequence\n",
     "    seq = Sequence(reg, Chadoq2)\n",
-    "    seq.declare_channel('ryd','rydberg_local','atom0')\n",
-    "    seq.add(pi_pulse,'ryd')\n",
-    "    seq.target('atom1', 'ryd')\n",
-    "    seq.add(pi_pulse,'ryd')\n",
+    "    seq.declare_channel(\"ryd\", \"rydberg_local\", \"atom0\")\n",
+    "    seq.add(pi_pulse, \"ryd\")\n",
+    "    seq.target(\"atom1\", \"ryd\")\n",
+    "    seq.add(pi_pulse, \"ryd\")\n",
     "\n",
     "    sim = Simulation(seq)\n",
-    "    \n",
+    "\n",
     "    res = sim.run()  # Returns a SimulationResults instance\n",
-    "    data.append(res.expect(occup)[0]) # Get expectation value for the occupation operator\n",
+    "    data.append(\n",
+    "        res.expect(occup)[0]\n",
+    "    )  # Get expectation value for the occupation operator\n",
     "print(\"...Simulation Complete!\")"
    ]
   },
@@ -247,10 +251,10 @@
     }
    ],
    "source": [
-    "for i,R in enumerate(distances):\n",
-    "    plt.plot(data[i], label=f'R={R}')\n",
-    "    plt.xlabel('Time (ns)', fontsize=14)\n",
-    "    plt.ylabel(r'Occupation of $|rr\\rangle$', fontsize=14)\n",
+    "for i, R in enumerate(distances):\n",
+    "    plt.plot(data[i], label=f\"R={R}\")\n",
+    "    plt.xlabel(\"Time (ns)\", fontsize=14)\n",
+    "    plt.ylabel(r\"Occupation of $|rr\\rangle$\", fontsize=14)\n",
     "    plt.legend()"
    ]
   },
diff --git a/pulser/__init__.py b/pulser/__init__.py
index fbebcfedf..eb8ea3b5a 100644
--- a/pulser/__init__.py
+++ b/pulser/__init__.py
@@ -15,11 +15,7 @@
 """A pulse-level composer for neutral-atom quantum devices."""
 
 from pulser._version import __version__
-
 from pulser.pulse import Pulse
-
 from pulser.register import Register, Register3D
-
 from pulser.sequence import Sequence
-
 from pulser.simulation import Simulation
diff --git a/pulser/_seq_drawer.py b/pulser/_seq_drawer.py
index d65d52614..e6e603b97 100644
--- a/pulser/_seq_drawer.py
+++ b/pulser/_seq_drawer.py
@@ -15,18 +15,18 @@
 from __future__ import annotations
 
 from collections import defaultdict
-from typing import Any, cast, Optional, Union
+from itertools import combinations
+from typing import Any, Optional, Union, cast
 
 import matplotlib.pyplot as plt
-from matplotlib.figure import Figure
 import numpy as np
+from matplotlib.figure import Figure
 from scipy.interpolate import CubicSpline
-from itertools import combinations
 
 import pulser
-from pulser.waveforms import ConstantWaveform, InterpolatedWaveform
-from pulser.pulse import Pulse
 from pulser import Register, Register3D
+from pulser.pulse import Pulse
+from pulser.waveforms import ConstantWaveform, InterpolatedWaveform
 
 
 def gather_data(seq: pulser.sequence.Sequence) -> dict:
diff --git a/pulser/channels.py b/pulser/channels.py
index 89015e9f6..1dd230b62 100644
--- a/pulser/channels.py
+++ b/pulser/channels.py
@@ -15,13 +15,13 @@
 
 from __future__ import annotations
 
-from dataclasses import dataclass
-from typing import cast, ClassVar, Optional
 import warnings
+from dataclasses import dataclass
+from typing import ClassVar, Optional, cast
 
 import numpy as np
 from numpy.typing import ArrayLike
-from scipy.fft import fft, ifft, fftfreq
+from scipy.fft import fft, fftfreq, ifft
 
 # Warnings of adjusted waveform duration appear just once
 warnings.filterwarnings("once", "A duration of")
diff --git a/pulser/devices/__init__.py b/pulser/devices/__init__.py
index 0660fdcfc..bf6d17fca 100644
--- a/pulser/devices/__init__.py
+++ b/pulser/devices/__init__.py
@@ -13,12 +13,10 @@
 # limitations under the License.
 """Valid devices for Pulser Sequence execution."""
 
+from pulser.devices._device_datacls import Device
 from pulser.devices._devices import Chadoq2
-
 from pulser.devices._mock_device import MockDevice
 
-from pulser.devices._device_datacls import Device
-
 # Registers which devices can be used to avoid definition of custom devices
 _mock_devices = (MockDevice,)
 _valid_devices = (Chadoq2,)
diff --git a/pulser/devices/_device_datacls.py b/pulser/devices/_device_datacls.py
index a904daa02..5003eef21 100644
--- a/pulser/devices/_device_datacls.py
+++ b/pulser/devices/_device_datacls.py
@@ -21,10 +21,10 @@
 from scipy.spatial.distance import pdist, squareform
 
 from pulser import Pulse
-from pulser.register.base_register import BaseRegister
 from pulser.channels import Channel
-from pulser.json.utils import obj_to_dict
 from pulser.devices.interaction_coefficients import c6_dict
+from pulser.json.utils import obj_to_dict
+from pulser.register.base_register import BaseRegister
 
 
 @dataclass(frozen=True, repr=False)
diff --git a/pulser/devices/_devices.py b/pulser/devices/_devices.py
index 0fa9e50e0..9a757ec91 100644
--- a/pulser/devices/_devices.py
+++ b/pulser/devices/_devices.py
@@ -14,9 +14,8 @@
 """Definitions of real devices."""
 import numpy as np
 
-from pulser.devices._device_datacls import Device
 from pulser.channels import Raman, Rydberg
-
+from pulser.devices._device_datacls import Device
 
 Chadoq2 = Device(
     name="Chadoq2",
diff --git a/pulser/devices/_mock_device.py b/pulser/devices/_mock_device.py
index 0b6e02216..01bf385e5 100644
--- a/pulser/devices/_mock_device.py
+++ b/pulser/devices/_mock_device.py
@@ -12,9 +12,8 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
+from pulser.channels import Microwave, Raman, Rydberg
 from pulser.devices._device_datacls import Device
-from pulser.channels import Rydberg, Raman, Microwave
-
 
 MockDevice = Device(
     name="MockDevice",
diff --git a/pulser/json/coders.py b/pulser/json/coders.py
index d4a79a999..011140d01 100644
--- a/pulser/json/coders.py
+++ b/pulser/json/coders.py
@@ -17,7 +17,7 @@
 
 import importlib
 import inspect
-from json import JSONEncoder, JSONDecoder
+from json import JSONDecoder, JSONEncoder
 from typing import Any, cast
 
 import numpy as np
diff --git a/pulser/parametrized/__init__.py b/pulser/parametrized/__init__.py
index 5bb0d5cd5..6f0775cf3 100644
--- a/pulser/parametrized/__init__.py
+++ b/pulser/parametrized/__init__.py
@@ -14,7 +14,5 @@
 """Classes for parametrized pulse-sequence building."""
 
 from pulser.parametrized.paramabc import Parametrized
-
-from pulser.parametrized.variable import Variable
-
 from pulser.parametrized.paramobj import ParamObj
+from pulser.parametrized.variable import Variable
diff --git a/pulser/parametrized/decorators.py b/pulser/parametrized/decorators.py
index ae27c70a2..a77df039a 100644
--- a/pulser/parametrized/decorators.py
+++ b/pulser/parametrized/decorators.py
@@ -22,7 +22,6 @@
 
 from pulser.parametrized import Parametrized, ParamObj
 
-
 F = TypeVar("F", bound=Callable)
 
 
diff --git a/pulser/parametrized/paramabc.py b/pulser/parametrized/paramabc.py
index fff36e656..f050ab736 100644
--- a/pulser/parametrized/paramabc.py
+++ b/pulser/parametrized/paramabc.py
@@ -16,7 +16,7 @@
 from __future__ import annotations
 
 from abc import ABC, abstractmethod
-from typing import Any, TYPE_CHECKING
+from typing import TYPE_CHECKING, Any
 
 if TYPE_CHECKING:
     from pulser.parametrized import Variable  # pragma: no cover
diff --git a/pulser/parametrized/paramobj.py b/pulser/parametrized/paramobj.py
index 33d173f73..0bcf0ad4a 100644
--- a/pulser/parametrized/paramobj.py
+++ b/pulser/parametrized/paramobj.py
@@ -15,13 +15,13 @@
 
 from __future__ import annotations
 
-from collections.abc import Callable
-from functools import partialmethod
-from itertools import chain
 import inspect
 import operator
 import warnings
-from typing import Any, Union, TYPE_CHECKING
+from collections.abc import Callable
+from functools import partialmethod
+from itertools import chain
+from typing import TYPE_CHECKING, Any, Union
 
 from pulser.json.utils import obj_to_dict
 from pulser.parametrized import Parametrized
diff --git a/pulser/parametrized/variable.py b/pulser/parametrized/variable.py
index 43bcd2914..743c337bb 100644
--- a/pulser/parametrized/variable.py
+++ b/pulser/parametrized/variable.py
@@ -16,16 +16,16 @@
 from __future__ import annotations
 
 import collections.abc  # To use collections.abc.Sequence
-from collections.abc import Iterable
 import dataclasses
-from typing import Union, Any, cast
+from collections.abc import Iterable
+from typing import Any, Union, cast
 
 import numpy as np
 from numpy.typing import ArrayLike
 
+from pulser.json.utils import obj_to_dict
 from pulser.parametrized import Parametrized
 from pulser.parametrized.paramobj import OpSupport
-from pulser.json.utils import obj_to_dict
 
 
 @dataclasses.dataclass(frozen=True, eq=False)
diff --git a/pulser/pulse.py b/pulser/pulse.py
index c3a2a75e7..e800bc7a5 100644
--- a/pulser/pulse.py
+++ b/pulser/pulse.py
@@ -15,19 +15,19 @@
 
 from __future__ import annotations
 
-from dataclasses import dataclass, field
 import functools
 import itertools
-from typing import Any, cast, Union
+from dataclasses import dataclass, field
+from typing import Any, Union, cast
 
 import matplotlib.pyplot as plt
 import numpy as np
 
 from pulser.channels import Channel
+from pulser.json.utils import obj_to_dict
 from pulser.parametrized import Parametrized, ParamObj
 from pulser.parametrized.decorators import parametrize
-from pulser.waveforms import Waveform, ConstantWaveform
-from pulser.json.utils import obj_to_dict
+from pulser.waveforms import ConstantWaveform, Waveform
 
 
 @dataclass(init=False, repr=False, frozen=True)
diff --git a/pulser/register/__init__.py b/pulser/register/__init__.py
index a57034208..084859366 100644
--- a/pulser/register/__init__.py
+++ b/pulser/register/__init__.py
@@ -14,7 +14,5 @@
 """Classes for qubit register definition."""
 
 from pulser.register.base_register import QubitId
-
 from pulser.register.register import Register
-
 from pulser.register.register3d import Register3D
diff --git a/pulser/register/base_register.py b/pulser/register/base_register.py
index 94ad37dae..d9437f192 100644
--- a/pulser/register/base_register.py
+++ b/pulser/register/base_register.py
@@ -16,13 +16,13 @@
 from __future__ import annotations
 
 from abc import ABC, abstractmethod
-from collections.abc import Mapping, Iterable
+from collections.abc import Iterable, Mapping
 from collections.abc import Sequence as abcSequence
-from typing import Any, cast, Optional, Union, TypeVar, Type
+from typing import Any, Optional, Type, TypeVar, Union, cast
 
 import matplotlib.pyplot as plt
-from matplotlib import collections as mc
 import numpy as np
+from matplotlib import collections as mc
 from numpy.typing import ArrayLike
 from scipy.spatial import KDTree
 
diff --git a/pulser/register/register3d.py b/pulser/register/register3d.py
index 04b5807b6..cbc5d7763 100644
--- a/pulser/register/register3d.py
+++ b/pulser/register/register3d.py
@@ -16,8 +16,8 @@
 from __future__ import annotations
 
 from collections.abc import Mapping
-from typing import Any, Optional
 from itertools import combinations
+from typing import Any, Optional
 
 import matplotlib.pyplot as plt
 import numpy as np
diff --git a/pulser/sequence.py b/pulser/sequence.py
index cdd85f12a..97df719c5 100644
--- a/pulser/sequence.py
+++ b/pulser/sequence.py
@@ -15,31 +15,31 @@
 
 from __future__ import annotations
 
+import copy
+import json
+import os
+import warnings
 from collections import namedtuple
 from collections.abc import Callable, Generator, Iterable
-import copy
 from functools import wraps
 from itertools import chain
-import json
 from sys import version_info
-from typing import Any, cast, NamedTuple, Optional, Tuple, TypeVar, Union
-import warnings
-import os
+from typing import Any, NamedTuple, Optional, Tuple, TypeVar, Union, cast
 
 import matplotlib.pyplot as plt
 import numpy as np
 from numpy.typing import ArrayLike
 
 import pulser
+from pulser._seq_drawer import draw_sequence
 from pulser.channels import Channel
 from pulser.devices import MockDevice
 from pulser.devices._device_datacls import Device
-from pulser.json.coders import PulserEncoder, PulserDecoder
+from pulser.json.coders import PulserDecoder, PulserEncoder
 from pulser.json.utils import obj_to_dict
 from pulser.parametrized import Parametrized, Variable
 from pulser.pulse import Pulse
 from pulser.register.base_register import BaseRegister
-from pulser._seq_drawer import draw_sequence
 
 if version_info[:2] >= (3, 8):  # pragma: no cover
     from typing import Literal, get_args
diff --git a/pulser/simulation/__init__.py b/pulser/simulation/__init__.py
index c32b17714..3d2ae5021 100644
--- a/pulser/simulation/__init__.py
+++ b/pulser/simulation/__init__.py
@@ -13,5 +13,5 @@
 # limitations under the License.
 """Classes for classical emulation of a Sequence."""
 
-from pulser.simulation.simulation import Simulation
 from pulser.simulation.simconfig import SimConfig
+from pulser.simulation.simulation import Simulation
diff --git a/pulser/simulation/simconfig.py b/pulser/simulation/simconfig.py
index 0222fd9f1..8d175541f 100644
--- a/pulser/simulation/simconfig.py
+++ b/pulser/simulation/simconfig.py
@@ -15,9 +15,9 @@
 
 from __future__ import annotations
 
-from sys import version_info
 from dataclasses import dataclass, field
-from typing import Union, Any
+from sys import version_info
+from typing import Any, Union
 
 import numpy as np
 import qutip
diff --git a/pulser/simulation/simresults.py b/pulser/simulation/simresults.py
index fc30c0079..f8c689462 100644
--- a/pulser/simulation/simresults.py
+++ b/pulser/simulation/simresults.py
@@ -15,17 +15,17 @@
 
 from __future__ import annotations
 
-from collections import Counter
 import collections.abc
 from abc import ABC, abstractmethod
+from collections import Counter
 from functools import lru_cache
-from typing import Optional, Union, cast, Tuple, Mapping
+from typing import Mapping, Optional, Tuple, Union, cast
 
 import matplotlib.pyplot as plt
-import qutip
-from qutip.piqs import isdiagonal
 import numpy as np
+import qutip
 from numpy.typing import ArrayLike
+from qutip.piqs import isdiagonal
 
 
 class SimulationResults(ABC):
diff --git a/pulser/simulation/simulation.py b/pulser/simulation/simulation.py
index aca3ecb19..a480eb3a0 100644
--- a/pulser/simulation/simulation.py
+++ b/pulser/simulation/simulation.py
@@ -15,30 +15,29 @@
 
 from __future__ import annotations
 
-from typing import Optional, Union, cast, Any
-from collections.abc import Mapping
 import itertools
+import warnings
 from collections import Counter
+from collections.abc import Mapping
 from copy import deepcopy
 from dataclasses import asdict
-import warnings
+from typing import Any, Optional, Union, cast
 
-import qutip
+import matplotlib.pyplot as plt
 import numpy as np
+import qutip
 from numpy.typing import ArrayLike
-import matplotlib.pyplot as plt
 
 from pulser import Pulse, Sequence
+from pulser._seq_drawer import draw_sequence
 from pulser.register import QubitId
+from pulser.sequence import _TimeSlot
+from pulser.simulation.simconfig import SimConfig
 from pulser.simulation.simresults import (
-    SimulationResults,
     CoherentResults,
     NoisyResults,
+    SimulationResults,
 )
-from pulser.simulation.simconfig import SimConfig
-from pulser._seq_drawer import draw_sequence
-from pulser.sequence import _TimeSlot
-
 
 SUPPORTED_NOISE = {
     "ising": {"dephasing", "doppler", "amplitude", "SPAM"},
diff --git a/pulser/waveforms.py b/pulser/waveforms.py
index 7daddbf0a..5a2757594 100644
--- a/pulser/waveforms.py
+++ b/pulser/waveforms.py
@@ -15,26 +15,26 @@
 
 from __future__ import annotations
 
-from abc import ABC, abstractmethod
 import functools
 import inspect
 import itertools
 import sys
+import warnings
+from abc import ABC, abstractmethod
 from sys import version_info
 from types import FunctionType
-from typing import Any, cast, Optional, Tuple, Union
-import warnings
+from typing import Any, Optional, Tuple, Union, cast
 
-from matplotlib.axes import Axes
 import matplotlib.pyplot as plt
 import numpy as np
-from numpy.typing import ArrayLike
 import scipy.interpolate as interpolate
+from matplotlib.axes import Axes
+from numpy.typing import ArrayLike
 
 from pulser.channels import Channel
+from pulser.json.utils import obj_to_dict
 from pulser.parametrized import Parametrized, ParamObj
 from pulser.parametrized.decorators import parametrize
-from pulser.json.utils import obj_to_dict
 
 if version_info[:2] >= (3, 8):  # pragma: no cover
     from functools import cached_property
diff --git a/pyproject.toml b/pyproject.toml
index a8f43fefd..1e75c2fbd 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,2 +1,5 @@
 [tool.black]
 line-length = 79
+
+[tool.isort]
+profile = "black"
diff --git a/requirements.txt b/requirements.txt
index c870d413d..df9ab926a 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -4,18 +4,20 @@ qutip
 scipy
 
 # version specific
-typing-extensions; python_version == '3.7'
 backports.cached-property; python_version == '3.7'
+typing-extensions; python_version == '3.7'
 
 # tests
-pytest
-pytest-cov
+black
+black[jupyter]
 flake8
 flake8-docstrings
+isort
 mypy == 0.921
-black
+pytest
+pytest-cov
 
 # tutorials
 notebook
-scikit-optimize
 python-igraph
+scikit-optimize
diff --git a/setup.py b/setup.py
index 33b94be63..db2211c6c 100644
--- a/setup.py
+++ b/setup.py
@@ -12,7 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-from setuptools import setup, find_packages
+from setuptools import find_packages, setup
 
 __version__ = ""
 exec(open("pulser/_version.py").read())
diff --git a/tutorials/advanced_features/Composite Waveforms.ipynb b/tutorials/advanced_features/Composite Waveforms.ipynb
index 3cb38b815..fa26b46b3 100644
--- a/tutorials/advanced_features/Composite Waveforms.ipynb
+++ b/tutorials/advanced_features/Composite Waveforms.ipynb
@@ -16,7 +16,12 @@
     "import numpy as np\n",
     "\n",
     "from pulser import Pulse\n",
-    "from pulser.waveforms import BlackmanWaveform, RampWaveform, CompositeWaveform, ConstantWaveform"
+    "from pulser.waveforms import (\n",
+    "    BlackmanWaveform,\n",
+    "    RampWaveform,\n",
+    "    CompositeWaveform,\n",
+    "    ConstantWaveform,\n",
+    ")"
    ]
   },
   {
@@ -37,7 +42,7 @@
    "outputs": [],
    "source": [
     "# Defining simple waveforms\n",
-    "pi_pulse = BlackmanWaveform(1000, np.pi)   # Blackman pi-pulse of 1us\n",
+    "pi_pulse = BlackmanWaveform(1000, np.pi)  # Blackman pi-pulse of 1us\n",
     "up = RampWaveform(500, 0, 5)\n",
     "down = RampWaveform(500, 5, 0)\n",
     "\n",
diff --git a/tutorials/advanced_features/Interpolated Waveforms.ipynb b/tutorials/advanced_features/Interpolated Waveforms.ipynb
index fdf8ea3f6..c3987ceb2 100644
--- a/tutorials/advanced_features/Interpolated Waveforms.ipynb
+++ b/tutorials/advanced_features/Interpolated Waveforms.ipynb
@@ -64,7 +64,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ts = np.r_[np.linspace(0., 0.5, num=len(values) - 1), 1]\n",
+    "ts = np.r_[np.linspace(0.0, 0.5, num=len(values) - 1), 1]\n",
     "int_wf_t = InterpolatedWaveform(duration, values, times=ts)\n",
     "int_wf_t.draw()"
    ]
@@ -103,7 +103,9 @@
    },
    "outputs": [],
    "source": [
-    "int_wf3 = InterpolatedWaveform(duration, values, interpolator=\"interp1d\", kind=\"cubic\")\n",
+    "int_wf3 = InterpolatedWaveform(\n",
+    "    duration, values, interpolator=\"interp1d\", kind=\"cubic\"\n",
+    ")\n",
     "int_wf3.draw()"
    ]
   },
@@ -163,7 +165,9 @@
     "det_vals = param_seq.declare_variable(\"det_vals\", size=4, dtype=float)\n",
     "\n",
     "amp_wf = InterpolatedWaveform(1000, amp_vals)\n",
-    "det_wf = InterpolatedWaveform(1000, det_vals, interpolator=\"interp1d\", kind=\"cubic\")\n",
+    "det_wf = InterpolatedWaveform(\n",
+    "    1000, det_vals, interpolator=\"interp1d\", kind=\"cubic\"\n",
+    ")\n",
     "pls = Pulse(amp_wf, det_wf, 0)\n",
     "\n",
     "param_seq.add(pls, \"rydberg_global\")"
@@ -182,7 +186,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "seq1 = param_seq.build(amp_vals=[0, 2, 1, 2, 0], det_vals = [0, -5, 5, -5])\n",
+    "seq1 = param_seq.build(amp_vals=[0, 2, 1, 2, 0], det_vals=[0, -5, 5, -5])\n",
     "seq1.draw()"
    ]
   },
diff --git a/tutorials/advanced_features/Parametrized Sequences.ipynb b/tutorials/advanced_features/Parametrized Sequences.ipynb
index df69949b5..e0c2273b1 100644
--- a/tutorials/advanced_features/Parametrized Sequences.ipynb
+++ b/tutorials/advanced_features/Parametrized Sequences.ipynb
@@ -35,10 +35,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "reg = Register.square(2, prefix='q')\n",
+    "reg = Register.square(2, prefix=\"q\")\n",
     "seq = Sequence(reg, Chadoq2)\n",
-    "seq.declare_channel('rydberg', 'rydberg_global')\n",
-    "seq.declare_channel('raman', 'raman_local')"
+    "seq.declare_channel(\"rydberg\", \"rydberg_global\")\n",
+    "seq.declare_channel(\"raman\", \"raman_local\")"
    ]
   },
   {
@@ -61,9 +61,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "Omega_max = seq.declare_variable('Omega_max')\n",
-    "ts = seq.declare_variable('ts', size=2, dtype=int)\n",
-    "last_target = seq.declare_variable('last_target', dtype=str)"
+    "Omega_max = seq.declare_variable(\"Omega_max\")\n",
+    "ts = seq.declare_variable(\"ts\", size=2, dtype=int)\n",
+    "last_target = seq.declare_variable(\"last_target\", dtype=str)"
    ]
   },
   {
@@ -85,7 +85,7 @@
     "U = Omega_max / 2.3\n",
     "delta_0 = -6 * U\n",
     "delta_f = 2 * U\n",
-    "t_sweep = (delta_f - delta_0)/(2 * np.pi * 10) * 1000"
+    "t_sweep = (delta_f - delta_0) / (2 * np.pi * 10) * 1000"
    ]
   },
   {
@@ -177,7 +177,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "generic_pulse = Pulse.ConstantPulse(100, 2*np.pi, 2, 0.)\n",
+    "generic_pulse = Pulse.ConstantPulse(100, 2 * np.pi, 2, 0.0)\n",
     "seq.add(generic_pulse, \"rydberg\")\n",
     "seq.target(\"q0\", \"raman\")\n",
     "seq.add(generic_pulse, \"raman\")\n",
@@ -256,7 +256,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "built_seq = seq.build(Omega_max = 2.3 * 2*np.pi,  ts = [200, 500], last_target=\"q3\")\n",
+    "built_seq = seq.build(\n",
+    "    Omega_max=2.3 * 2 * np.pi, ts=[200, 500], last_target=\"q3\"\n",
+    ")\n",
     "built_seq.draw()"
    ]
   },
@@ -273,7 +275,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "alt_seq = seq.build(Omega_max = 2*np.pi,  ts = [400, 100], last_target=\"q2\")\n",
+    "alt_seq = seq.build(Omega_max=2 * np.pi, ts=[400, 100], last_target=\"q2\")\n",
     "alt_seq.draw()"
    ]
   }
diff --git a/tutorials/advanced_features/Phase Shifts and Virtual Z gates.ipynb b/tutorials/advanced_features/Phase Shifts and Virtual Z gates.ipynb
index ff724dcec..3422c18e6 100644
--- a/tutorials/advanced_features/Phase Shifts and Virtual Z gates.ipynb
+++ b/tutorials/advanced_features/Phase Shifts and Virtual Z gates.ipynb
@@ -124,7 +124,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "reg = Register({'q0': (0, 0)})\n",
+    "reg = Register({\"q0\": (0, 0)})\n",
     "device = MockDevice\n",
     "seq = Sequence(reg, device)\n",
     "seq.available_channels"
@@ -136,7 +136,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "seq.declare_channel('ch0', 'raman_local', initial_target = 'q0')"
+    "seq.declare_channel(\"ch0\", \"raman_local\", initial_target=\"q0\")"
    ]
   },
   {
@@ -147,7 +147,8 @@
    "source": [
     "# Defining the waveform for a pi/2 pulse\n",
     "from pulser.waveforms import BlackmanWaveform\n",
-    "pi2_wf = BlackmanWaveform(1000, np.pi/2)   # Duration: 1us, Area: pi/2\n",
+    "\n",
+    "pi2_wf = BlackmanWaveform(1000, np.pi / 2)  # Duration: 1us, Area: pi/2\n",
     "pi2_wf.draw()"
    ]
   },
@@ -160,7 +161,7 @@
    "outputs": [],
    "source": [
     "# 2. Create the pi/2 pulse\n",
-    "pi_2 = Pulse.ConstantDetuning(pi2_wf, detuning=0, phase=np.pi/2)\n",
+    "pi_2 = Pulse.ConstantDetuning(pi2_wf, detuning=0, phase=np.pi / 2)\n",
     "pi_2.draw()"
    ]
   },
@@ -172,11 +173,11 @@
    },
    "outputs": [],
    "source": [
-    "#3. Applying the H gate\n",
+    "# 3. Applying the H gate\n",
     "\n",
-    "seq.add(pi_2, 'ch0')   # The first pi/2-pulse\n",
+    "seq.add(pi_2, \"ch0\")  # The first pi/2-pulse\n",
     "# Now the phase shift of pi on 'q0', for the 'digital' basis, which is usually where phase shifts are useful\n",
-    "seq.phase_shift(np.pi, 'q0', basis='digital') \n",
+    "seq.phase_shift(np.pi, \"q0\", basis=\"digital\")\n",
     "\n",
     "seq.draw(draw_phase_shifts=True)"
    ]
@@ -198,9 +199,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "h = Pulse.ConstantDetuning(pi2_wf, detuning=0, phase=np.pi/2, post_phase_shift=np.pi)\n",
+    "h = Pulse.ConstantDetuning(\n",
+    "    pi2_wf, detuning=0, phase=np.pi / 2, post_phase_shift=np.pi\n",
+    ")\n",
     "\n",
-    "seq.add(h, 'ch0')\n",
+    "seq.add(h, \"ch0\")\n",
     "seq.draw(draw_phase_shifts=True)"
    ]
   },
@@ -233,7 +236,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "seq.add(h, 'ch0')\n",
+    "seq.add(h, \"ch0\")\n",
     "print(seq)"
    ]
   },
@@ -264,12 +267,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "reg = Register({'q0': (0, 0), 'q1': (5, 5)})\n",
+    "reg = Register({\"q0\": (0, 0), \"q1\": (5, 5)})\n",
     "device = MockDevice\n",
     "seq = Sequence(reg, device)\n",
-    "seq.declare_channel('raman', 'raman_local', initial_target = 'q0')\n",
-    "seq.declare_channel('ryd1', 'rydberg_local', initial_target = 'q0')\n",
-    "seq.declare_channel('ryd2', 'rydberg_local', initial_target = 'q0')\n",
+    "seq.declare_channel(\"raman\", \"raman_local\", initial_target=\"q0\")\n",
+    "seq.declare_channel(\"ryd1\", \"rydberg_local\", initial_target=\"q0\")\n",
+    "seq.declare_channel(\"ryd2\", \"rydberg_local\", initial_target=\"q0\")\n",
     "seq.declared_channels"
    ]
   },
@@ -286,8 +289,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "seq.add(h, 'raman')\n",
-    "seq.add(h, 'ryd1')\n",
+    "seq.add(h, \"raman\")\n",
+    "seq.add(h, \"ryd1\")\n",
     "seq.draw(draw_phase_shifts=True)"
    ]
   },
@@ -320,7 +323,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "seq.add(pi_2, 'ryd2')\n",
+    "seq.add(pi_2, \"ryd2\")\n",
     "print(seq)"
    ]
   },
@@ -339,10 +342,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "seq.target('q1', 'raman')\n",
-    "seq.add(h, 'raman')\n",
-    "seq.target('q1', 'ryd1')\n",
-    "seq.add(h, 'ryd1')\n",
+    "seq.target(\"q1\", \"raman\")\n",
+    "seq.add(h, \"raman\")\n",
+    "seq.target(\"q1\", \"ryd1\")\n",
+    "seq.add(h, \"ryd1\")\n",
     "print(seq)\n",
     "seq.draw(draw_phase_shifts=True)"
    ]
@@ -362,7 +365,7 @@
    },
    "outputs": [],
    "source": [
-    "seq.target('q1', 'ryd2')\n",
+    "seq.target(\"q1\", \"ryd2\")\n",
     "seq.draw(draw_phase_shifts=True)"
    ]
   },
diff --git a/tutorials/advanced_features/Serialization.ipynb b/tutorials/advanced_features/Serialization.ipynb
index 295ed5082..c9883f621 100644
--- a/tutorials/advanced_features/Serialization.ipynb
+++ b/tutorials/advanced_features/Serialization.ipynb
@@ -42,10 +42,12 @@
     "seq.declare_channel(\"digital\", \"raman_local\", initial_target=\"control\")\n",
     "seq.declare_channel(\"rydberg\", \"rydberg_local\", initial_target=\"control\")\n",
     "\n",
-    "half_pi_wf = BlackmanWaveform(200, np.pi/2)\n",
+    "half_pi_wf = BlackmanWaveform(200, np.pi / 2)\n",
     "\n",
-    "ry = Pulse.ConstantDetuning(amplitude=half_pi_wf,detuning=0,phase=-np.pi/2)\n",
-    "ry_dag = Pulse.ConstantDetuning(amplitude=half_pi_wf,detuning=0,phase=np.pi/2)\n",
+    "ry = Pulse.ConstantDetuning(amplitude=half_pi_wf, detuning=0, phase=-np.pi / 2)\n",
+    "ry_dag = Pulse.ConstantDetuning(\n",
+    "    amplitude=half_pi_wf, detuning=0, phase=np.pi / 2\n",
+    ")\n",
     "\n",
     "seq.add(ry, \"digital\")\n",
     "seq.target(\"target\", \"digital\")\n",
@@ -55,7 +57,7 @@
     "pi_pulse = Pulse.ConstantDetuning(pi_wf, 0, 0)\n",
     "\n",
     "max_val = Chadoq2.rabi_from_blockade(8)\n",
-    "two_pi_wf = BlackmanWaveform.from_max_val(max_val, 2*np.pi)\n",
+    "two_pi_wf = BlackmanWaveform.from_max_val(max_val, 2 * np.pi)\n",
     "two_pi_pulse = Pulse.ConstantDetuning(two_pi_wf, 0, 0)\n",
     "\n",
     "seq.align(\"digital\", \"rydberg\")\n",
diff --git a/tutorials/advanced_features/Simulating Sequences with Errors and Noise.ipynb b/tutorials/advanced_features/Simulating Sequences with Errors and Noise.ipynb
index fecbad9e7..155822b7f 100644
--- a/tutorials/advanced_features/Simulating Sequences with Errors and Noise.ipynb
+++ b/tutorials/advanced_features/Simulating Sequences with Errors and Noise.ipynb
@@ -69,7 +69,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "reg = Register.from_coordinates([(0,0)], prefix='q')"
+    "reg = Register.from_coordinates([(0, 0)], prefix=\"q\")"
    ]
   },
   {
@@ -99,10 +99,10 @@
    ],
    "source": [
     "seq = Sequence(reg, Chadoq2)\n",
-    "seq.declare_channel('ch0', 'rydberg_global')\n",
+    "seq.declare_channel(\"ch0\", \"rydberg_global\")\n",
     "duration = 2500\n",
-    "pulse = Pulse.ConstantPulse(duration, 2*np.pi, 0., 0.)\n",
-    "seq.add(pulse, 'ch0')\n",
+    "pulse = Pulse.ConstantPulse(duration, 2 * np.pi, 0.0, 0.0)\n",
+    "seq.add(pulse, \"ch0\")\n",
     "seq.draw()"
    ]
   },
@@ -136,7 +136,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "obs = qutip.basis(2,0).proj()"
+    "obs = qutip.basis(2, 0).proj()"
    ]
   },
   {
@@ -216,7 +216,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "config_spam = SimConfig(noise=('SPAM'), runs = 30, samples_per_run = 5)\n",
+    "config_spam = SimConfig(noise=(\"SPAM\"), runs=30, samples_per_run=5)\n",
     "sim.set_config(config_spam)"
    ]
   },
@@ -284,7 +284,12 @@
     }
    ],
    "source": [
-    "cfg2 = SimConfig(noise=('SPAM', 'dephasing', 'doppler'), eta=0.8, temperature=1000, runs=10000)\n",
+    "cfg2 = SimConfig(\n",
+    "    noise=(\"SPAM\", \"dephasing\", \"doppler\"),\n",
+    "    eta=0.8,\n",
+    "    temperature=1000,\n",
+    "    runs=10000,\n",
+    ")\n",
     "sim.add_config(cfg2)\n",
     "sim.show_config()"
    ]
@@ -316,7 +321,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sim.evaluation_times = .8"
+    "sim.evaluation_times = 0.8"
    ]
   },
   {
@@ -401,7 +406,7 @@
     }
    ],
    "source": [
-    "res.plot(obs, fmt='.')\n",
+    "res.plot(obs, fmt=\".\")\n",
     "plt.show()"
    ]
   },
@@ -424,7 +429,7 @@
     }
    ],
    "source": [
-    "res.plot(obs, error_bars=False, fmt='.')"
+    "res.plot(obs, error_bars=False, fmt=\".\")"
    ]
   },
   {
@@ -467,7 +472,7 @@
     "res_spam = sim.run()\n",
     "res_spam.plot(obs)\n",
     "sim.reset_config()\n",
-    "sim.eval_times = 'Full'\n",
+    "sim.eval_times = \"Full\"\n",
     "res_clean = sim.run()\n",
     "res_clean.plot(obs)\n",
     "plt.show()"
@@ -499,7 +504,7 @@
     }
    ],
    "source": [
-    "config_spam_mod = SimConfig(noise=('SPAM'), eta=0.4, runs = 100)\n",
+    "config_spam_mod = SimConfig(noise=(\"SPAM\"), eta=0.4, runs=100)\n",
     "sim.set_config(config_spam_mod)\n",
     "sim.evaluation_times = 0.5\n",
     "res_large_eta = sim.run()\n",
@@ -546,12 +551,14 @@
     }
    ],
    "source": [
-    "plt.figure(figsize=(10,5))\n",
+    "plt.figure(figsize=(10, 5))\n",
     "res_clean.plot(obs)\n",
-    "for eta in np.linspace(0,0.99,4):\n",
-    "    config_spam_eta = SimConfig(noise = 'SPAM', eta=eta, runs = 50, epsilon=0, epsilon_prime=0)\n",
+    "for eta in np.linspace(0, 0.99, 4):\n",
+    "    config_spam_eta = SimConfig(\n",
+    "        noise=\"SPAM\", eta=eta, runs=50, epsilon=0, epsilon_prime=0\n",
+    "    )\n",
     "    sim.set_config(config_spam_eta)\n",
-    "    sim.run().plot(obs, label=f'eta = {eta}')\n",
+    "    sim.run().plot(obs, label=f\"eta = {eta}\")\n",
     "plt.legend()\n",
     "plt.show()"
    ]
@@ -596,12 +603,14 @@
     }
    ],
    "source": [
-    "plt.figure(figsize=(10,5))\n",
+    "plt.figure(figsize=(10, 5))\n",
     "res_clean.plot(obs)\n",
-    "for eps in np.linspace(0,.99,4):\n",
-    "    config_spam_eps = SimConfig(noise = 'SPAM', eta=0, runs = 50, epsilon=eps, epsilon_prime=0)\n",
+    "for eps in np.linspace(0, 0.99, 4):\n",
+    "    config_spam_eps = SimConfig(\n",
+    "        noise=\"SPAM\", eta=0, runs=50, epsilon=eps, epsilon_prime=0\n",
+    "    )\n",
     "    sim.set_config(config_spam_eps)\n",
-    "    sim.run().plot(obs, label=f'epsilon = {eps}')\n",
+    "    sim.run().plot(obs, label=f\"epsilon = {eps}\")\n",
     "plt.legend()\n",
     "plt.show()"
    ]
@@ -646,12 +655,14 @@
     }
    ],
    "source": [
-    "plt.figure(figsize=(10,5))\n",
+    "plt.figure(figsize=(10, 5))\n",
     "res_clean.plot(obs)\n",
-    "for eps_p in np.linspace(0,.99,4):\n",
-    "    config_spam_eps_p = SimConfig(noise = 'SPAM', eta=0, runs = 50, epsilon=0, epsilon_prime=eps_p)\n",
+    "for eps_p in np.linspace(0, 0.99, 4):\n",
+    "    config_spam_eps_p = SimConfig(\n",
+    "        noise=\"SPAM\", eta=0, runs=50, epsilon=0, epsilon_prime=eps_p\n",
+    "    )\n",
     "    sim.set_config(config_spam_eps_p)\n",
-    "    sim.run().plot(obs, label=f'epsilon = {eps_p}')\n",
+    "    sim.run().plot(obs, label=f\"epsilon = {eps_p}\")\n",
     "plt.legend()\n",
     "plt.show()"
    ]
@@ -703,7 +714,9 @@
     }
    ],
    "source": [
-    "config_doppler = SimConfig(noise='doppler', runs=100, temperature = 5000, samples_per_run=1)\n",
+    "config_doppler = SimConfig(\n",
+    "    noise=\"doppler\", runs=100, temperature=5000, samples_per_run=1\n",
+    ")\n",
     "sim.set_config(config_doppler)\n",
     "sim.show_config()"
    ]
@@ -763,28 +776,34 @@
    "outputs": [],
    "source": [
     "# Parameters in rad/µs and ns\n",
-    "Omega_max = 2.3 * 2*np.pi \n",
+    "Omega_max = 2.3 * 2 * np.pi\n",
     "U = Omega_max / 2.3\n",
     "delta_0 = -6 * U\n",
     "delta_f = 2 * U\n",
     "t_rise = 252\n",
     "t_fall = 500\n",
-    "t_sweep = (delta_f - delta_0)/(2 * np.pi * 10) * 1000\n",
+    "t_sweep = (delta_f - delta_0) / (2 * np.pi * 10) * 1000\n",
     "R_interatomic = Chadoq2.rydberg_blockade_radius(U)\n",
     "\n",
     "N_side = 3\n",
-    "reg = Register.rectangle(N_side, N_side, R_interatomic, prefix='q')\n",
+    "reg = Register.rectangle(N_side, N_side, R_interatomic, prefix=\"q\")\n",
     "\n",
-    "rise = Pulse.ConstantDetuning(RampWaveform(t_rise, 0., Omega_max), delta_0, 0.)\n",
-    "sweep = Pulse.ConstantAmplitude(Omega_max, RampWaveform(t_sweep, delta_0, delta_f), 0.)\n",
-    "fall = Pulse.ConstantDetuning(RampWaveform(t_fall, Omega_max, 0.), delta_f, 0.)\n",
+    "rise = Pulse.ConstantDetuning(\n",
+    "    RampWaveform(t_rise, 0.0, Omega_max), delta_0, 0.0\n",
+    ")\n",
+    "sweep = Pulse.ConstantAmplitude(\n",
+    "    Omega_max, RampWaveform(t_sweep, delta_0, delta_f), 0.0\n",
+    ")\n",
+    "fall = Pulse.ConstantDetuning(\n",
+    "    RampWaveform(t_fall, Omega_max, 0.0), delta_f, 0.0\n",
+    ")\n",
     "\n",
     "seq = Sequence(reg, Chadoq2)\n",
-    "seq.declare_channel('ising', 'rydberg_global')\n",
+    "seq.declare_channel(\"ising\", \"rydberg_global\")\n",
     "\n",
-    "seq.add(rise, 'ising')\n",
-    "seq.add(sweep, 'ising')\n",
-    "seq.add(fall, 'ising')"
+    "seq.add(rise, \"ising\")\n",
+    "seq.add(sweep, \"ising\")\n",
+    "seq.add(fall, \"ising\")"
    ]
   },
   {
@@ -793,9 +812,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "config_all_noise = SimConfig(noise=('SPAM', 'doppler', 'amplitude'),\n",
-    "                             runs=100, samples_per_run=10)\n",
-    "simul = Simulation(seq, sampling_rate=0.05, evaluation_times=0.2, config=config_all_noise)\n",
+    "config_all_noise = SimConfig(\n",
+    "    noise=(\"SPAM\", \"doppler\", \"amplitude\"), runs=100, samples_per_run=10\n",
+    ")\n",
+    "simul = Simulation(\n",
+    "    seq, sampling_rate=0.05, evaluation_times=0.2, config=config_all_noise\n",
+    ")\n",
     "spam_results = simul.run()\n",
     "simul.reset_config()\n",
     "clean_results = simul.run()"
@@ -827,17 +849,19 @@
     }
    ],
    "source": [
-    "plt.figure(figsize=(20,5))\n",
+    "plt.figure(figsize=(20, 5))\n",
     "spam_count = spam_results.sample_final_state(N_samples=1e5)\n",
     "clean_count = clean_results.sample_final_state(N_samples=1e5)\n",
     "\n",
-    "clean_most_freq = {k:v for k,v in clean_count.items() if v>500}\n",
-    "spam_most_freq = {k:v for k,v in spam_count.items() if v>500}\n",
+    "clean_most_freq = {k: v for k, v in clean_count.items() if v > 500}\n",
+    "spam_most_freq = {k: v for k, v in spam_count.items() if v > 500}\n",
     "\n",
-    "plt.bar(list(clean_most_freq.keys()), list(clean_most_freq.values()), width=0.9)\n",
+    "plt.bar(\n",
+    "    list(clean_most_freq.keys()), list(clean_most_freq.values()), width=0.9\n",
+    ")\n",
     "plt.bar(list(spam_most_freq.keys()), list(spam_most_freq.values()), width=0.5)\n",
     "\n",
-    "plt.xticks(rotation='vertical')\n",
+    "plt.xticks(rotation=\"vertical\")\n",
     "plt.show()"
    ]
   },
diff --git a/tutorials/advanced_features/State Preparation with the SLM Mask.ipynb b/tutorials/advanced_features/State Preparation with the SLM Mask.ipynb
index cbe7cb6e8..6f32fd09a 100644
--- a/tutorials/advanced_features/State Preparation with the SLM Mask.ipynb
+++ b/tutorials/advanced_features/State Preparation with the SLM Mask.ipynb
@@ -36,7 +36,7 @@
     "from pulser.simulation import Simulation\n",
     "\n",
     "# Qubit register\n",
-    "qubits = {\"q0\": (-5,0), \"q1\": (0,0), \"q2\": (5,0)}\n",
+    "qubits = {\"q0\": (-5, 0), \"q1\": (0, 0), \"q2\": (5, 0)}\n",
     "reg = Register(qubits)\n",
     "reg.draw()"
    ]
@@ -58,9 +58,9 @@
     "seq = Sequence(reg, MockDevice)\n",
     "\n",
     "# Declare a global XY channel and add the pi pulse\n",
-    "seq.declare_channel('ch', 'mw_global')\n",
+    "seq.declare_channel(\"ch\", \"mw_global\")\n",
     "pulse = Pulse.ConstantDetuning(BlackmanWaveform(200, np.pi), 0, 0)\n",
-    "seq.add(pulse, 'ch')"
+    "seq.add(pulse, \"ch\")"
    ]
   },
   {
@@ -193,7 +193,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "seq.add(pulse, 'ch')"
+    "seq.add(pulse, \"ch\")"
    ]
   },
   {
diff --git a/tutorials/applications/Control-Z Gate Sequence.ipynb b/tutorials/applications/Control-Z Gate Sequence.ipynb
index 438a2812f..8e0aaad92 100644
--- a/tutorials/applications/Control-Z Gate Sequence.ipynb
+++ b/tutorials/applications/Control-Z Gate Sequence.ipynb
@@ -58,7 +58,7 @@
     "from pulser import Pulse, Sequence, Register\n",
     "from pulser.devices import Chadoq2\n",
     "from pulser.simulation import Simulation\n",
-    "from pulser.waveforms import BlackmanWaveform,ConstantWaveform"
+    "from pulser.waveforms import BlackmanWaveform, ConstantWaveform"
    ]
   },
   {
@@ -86,10 +86,12 @@
    "outputs": [],
    "source": [
     "Rabi = np.linspace(1, 10, 10)\n",
-    "R_blockade = [Chadoq2.rydberg_blockade_radius(2.*np.pi*rabi) for rabi in Rabi]\n",
+    "R_blockade = [\n",
+    "    Chadoq2.rydberg_blockade_radius(2.0 * np.pi * rabi) for rabi in Rabi\n",
+    "]\n",
     "\n",
     "plt.figure()\n",
-    "plt.plot(Rabi, R_blockade,'--o')\n",
+    "plt.plot(Rabi, R_blockade, \"--o\")\n",
     "plt.xlabel(r\"$\\Omega/(2\\pi)$ [MHz]\", fontsize=14)\n",
     "plt.ylabel(r\"$R_b$ [$\\mu\\.m$]\", fontsize=14)\n",
     "plt.show()"
@@ -109,9 +111,10 @@
    "outputs": [],
    "source": [
     "# Atom Register and Device\n",
-    "q_dict = {\"control\":np.array([-2,0.]),\n",
-    "          \"target\": np.array([2,0.]),\n",
-    "          }\n",
+    "q_dict = {\n",
+    "    \"control\": np.array([-2, 0.0]),\n",
+    "    \"target\": np.array([2, 0.0]),\n",
+    "}\n",
     "reg = Register(q_dict)\n",
     "reg.draw()"
    ]
@@ -144,22 +147,24 @@
    "outputs": [],
    "source": [
     "def build_state_from_id(s_id, basis_name):\n",
-    "    if len(s_id) not in {2,3}:\n",
+    "    if len(s_id) not in {2, 3}:\n",
     "        raise ValueError(\"Not a valid state ID string\")\n",
-    "        \n",
-    "    ids = {'digital': 'gh', 'ground-rydberg': 'rg', 'all': 'rgh'}\n",
+    "\n",
+    "    ids = {\"digital\": \"gh\", \"ground-rydberg\": \"rg\", \"all\": \"rgh\"}\n",
     "    if basis_name not in ids:\n",
-    "        raise ValueError('Not a valid basis')\n",
-    "        \n",
-    "    pool = {''.join(x) for x in product(ids[basis_name], repeat=len(s_id))}\n",
+    "        raise ValueError(\"Not a valid basis\")\n",
+    "\n",
+    "    pool = {\"\".join(x) for x in product(ids[basis_name], repeat=len(s_id))}\n",
     "    if s_id not in pool:\n",
-    "        raise ValueError('Not a valid state id for the given basis.')\n",
+    "        raise ValueError(\"Not a valid state id for the given basis.\")\n",
     "\n",
-    "    ket = {op: qutip.basis(len(ids[basis_name]), i) \n",
-    "           for i, op in enumerate(ids[basis_name])}\n",
+    "    ket = {\n",
+    "        op: qutip.basis(len(ids[basis_name]), i)\n",
+    "        for i, op in enumerate(ids[basis_name])\n",
+    "    }\n",
     "    if len(s_id) == 3:\n",
-    "        #Recall that s_id = 'C1'+'C2'+'T' while in the register reg_id = 'C1'+'T'+'C2'.\n",
-    "        reg_id = s_id[0]+s_id[2]+s_id[1] \n",
+    "        # Recall that s_id = 'C1'+'C2'+'T' while in the register reg_id = 'C1'+'T'+'C2'.\n",
+    "        reg_id = s_id[0] + s_id[2] + s_id[1]\n",
     "        return qutip.tensor([ket[x] for x in reg_id])\n",
     "    else:\n",
     "        return qutip.tensor([ket[x] for x in s_id])"
@@ -178,7 +183,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "build_state_from_id('hg','digital')"
+    "build_state_from_id(\"hg\", \"digital\")"
    ]
   },
   {
@@ -194,8 +199,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "duration = 300 \n",
-    "pi_Y = Pulse.ConstantDetuning(BlackmanWaveform(duration, np.pi), 0., -np.pi/2)\n",
+    "duration = 300\n",
+    "pi_Y = Pulse.ConstantDetuning(\n",
+    "    BlackmanWaveform(duration, np.pi), 0.0, -np.pi / 2\n",
+    ")\n",
     "pi_Y.draw()"
    ]
   },
@@ -214,31 +221,37 @@
    "source": [
     "def preparation_sequence(state_id, reg):\n",
     "    global seq\n",
-    "    \n",
-    "    if not set(state_id) <= {'g','h'} or len(state_id) != len(reg.qubits):\n",
-    "        raise ValueError('Not a valid state ID')\n",
+    "\n",
+    "    if not set(state_id) <= {\"g\", \"h\"} or len(state_id) != len(reg.qubits):\n",
+    "        raise ValueError(\"Not a valid state ID\")\n",
     "\n",
     "    if len(reg.qubits) == 2:\n",
-    "        seq_dict = {'1':'target', '0':'control'}\n",
+    "        seq_dict = {\"1\": \"target\", \"0\": \"control\"}\n",
     "    elif len(reg.qubits) == 3:\n",
-    "        seq_dict = {'2':'target', '1':'control2', '0':'control1'}\n",
+    "        seq_dict = {\"2\": \"target\", \"1\": \"control2\", \"0\": \"control1\"}\n",
     "\n",
     "    seq = Sequence(reg, Chadoq2)\n",
-    "    if set(state_id) == {'g'}:\n",
-    "        basis = 'ground-rydberg'\n",
-    "        print(f'Warning: {state_id} state does not require a preparation sequence.')\n",
+    "    if set(state_id) == {\"g\"}:\n",
+    "        basis = \"ground-rydberg\"\n",
+    "        print(\n",
+    "            f\"Warning: {state_id} state does not require a preparation sequence.\"\n",
+    "        )\n",
     "    else:\n",
-    "        basis = 'all'\n",
+    "        basis = \"all\"\n",
     "        for k in range(len(reg.qubits)):\n",
-    "            if state_id[k] == 'h':\n",
-    "                if 'raman' not in seq.declared_channels:\n",
-    "                    seq.declare_channel('raman','raman_local', seq_dict[str(k)])\n",
+    "            if state_id[k] == \"h\":\n",
+    "                if \"raman\" not in seq.declared_channels:\n",
+    "                    seq.declare_channel(\n",
+    "                        \"raman\", \"raman_local\", seq_dict[str(k)]\n",
+    "                    )\n",
     "                else:\n",
-    "                    seq.target(seq_dict[str(k)],'raman')\n",
-    "                seq.add(pi_Y,'raman')\n",
+    "                    seq.target(seq_dict[str(k)], \"raman\")\n",
+    "                seq.add(pi_Y, \"raman\")\n",
+    "\n",
+    "    prep_state = build_state_from_id(\n",
+    "        state_id, basis\n",
+    "    )  # Raises error if not a valid `state_id` for the register\n",
     "\n",
-    "    prep_state = build_state_from_id(state_id, basis)  # Raises error if not a valid `state_id` for the register\n",
-    "    \n",
     "    return prep_state"
    ]
   },
@@ -256,7 +269,7 @@
    "outputs": [],
    "source": [
     "# Define sequence and Set channels\n",
-    "prep_state = preparation_sequence('hh', reg)\n",
+    "prep_state = preparation_sequence(\"hh\", reg)\n",
     "seq.draw(draw_phase_area=True)"
    ]
   },
@@ -280,8 +293,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pi_pulse = Pulse.ConstantDetuning(BlackmanWaveform(duration, np.pi), 0., 0)\n",
-    "twopi_pulse = Pulse.ConstantDetuning(BlackmanWaveform(duration, 2*np.pi), 0., 0)"
+    "pi_pulse = Pulse.ConstantDetuning(BlackmanWaveform(duration, np.pi), 0.0, 0)\n",
+    "twopi_pulse = Pulse.ConstantDetuning(\n",
+    "    BlackmanWaveform(duration, 2 * np.pi), 0.0, 0\n",
+    ")"
    ]
   },
   {
@@ -291,21 +306,25 @@
    "outputs": [],
    "source": [
     "def CZ_sequence(initial_id):\n",
-    "    \n",
+    "\n",
     "    # Prepare State\n",
-    "    prep_state = preparation_sequence(initial_id, reg) \n",
-    "    prep_time = max((seq._last(ch).tf for ch in seq.declared_channels), default=0)\n",
-    "    \n",
+    "    prep_state = preparation_sequence(initial_id, reg)\n",
+    "    prep_time = max(\n",
+    "        (seq._last(ch).tf for ch in seq.declared_channels), default=0\n",
+    "    )\n",
+    "\n",
     "    # Declare Rydberg channel\n",
-    "    seq.declare_channel('ryd', 'rydberg_local', 'control')\n",
-    "    \n",
+    "    seq.declare_channel(\"ryd\", \"rydberg_local\", \"control\")\n",
+    "\n",
     "    # Write CZ sequence:\n",
-    "    seq.add(pi_pulse, 'ryd', 'wait-for-all')  # Wait for state preparation to finish.\n",
-    "    seq.target('target', 'ryd')    # Changes to target qubit\n",
-    "    seq.add(twopi_pulse, 'ryd')\n",
-    "    seq.target('control', 'ryd')     # Changes back to control qubit\n",
-    "    seq.add(pi_pulse, 'ryd') \n",
-    "    \n",
+    "    seq.add(\n",
+    "        pi_pulse, \"ryd\", \"wait-for-all\"\n",
+    "    )  # Wait for state preparation to finish.\n",
+    "    seq.target(\"target\", \"ryd\")  # Changes to target qubit\n",
+    "    seq.add(twopi_pulse, \"ryd\")\n",
+    "    seq.target(\"control\", \"ryd\")  # Changes back to control qubit\n",
+    "    seq.add(pi_pulse, \"ryd\")\n",
+    "\n",
     "    return prep_state, prep_time"
    ]
   },
@@ -315,10 +334,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "prep_state, prep_time = CZ_sequence('gh') # constructs seq, prep_state and prep_time\n",
+    "prep_state, prep_time = CZ_sequence(\n",
+    "    \"gh\"\n",
+    ")  # constructs seq, prep_state and prep_time\n",
     "seq.draw(draw_phase_area=True)\n",
-    "print(f'Prepared state: {prep_state}')\n",
-    "print(f'Preparation time: {prep_time}ns')"
+    "print(f\"Prepared state: {prep_state}\")\n",
+    "print(f\"Preparation time: {prep_time}ns\")"
    ]
   },
   {
@@ -335,25 +356,27 @@
    "outputs": [],
    "source": [
     "CZ = {}\n",
-    "for state_id in {'gg','hg','gh','hh'}:\n",
+    "for state_id in {\"gg\", \"hg\", \"gh\", \"hh\"}:\n",
     "    # Get CZ sequence\n",
-    "    prep_state, prep_time = CZ_sequence(state_id) # constructs seq, prep_state and prep_time\n",
-    "    \n",
+    "    prep_state, prep_time = CZ_sequence(\n",
+    "        state_id\n",
+    "    )  # constructs seq, prep_state and prep_time\n",
+    "\n",
     "    # Construct Simulation instance\n",
     "    simul = Simulation(seq)\n",
     "    res = simul.run()\n",
-    "    \n",
-    "    data=[st.overlap(prep_state) for st in res.states]\n",
-    "    \n",
+    "\n",
+    "    data = [st.overlap(prep_state) for st in res.states]\n",
+    "\n",
     "    final_st = res.states[-1]\n",
     "    CZ[state_id] = final_st.overlap(prep_state)\n",
-    "    \n",
+    "\n",
     "    plt.figure()\n",
     "    plt.plot(np.real(data))\n",
     "    plt.xlabel(r\"Time [ns]\")\n",
-    "    plt.ylabel(fr'$ \\langle\\,{state_id} |\\, \\psi(t)\\rangle$')\n",
-    "    plt.axvspan(0, prep_time, alpha=0.06, color='royalblue')\n",
-    "    plt.title(fr\"Action of gate on state $|${state_id}$\\rangle$\")"
+    "    plt.ylabel(rf\"$ \\langle\\,{state_id} |\\, \\psi(t)\\rangle$\")\n",
+    "    plt.axvspan(0, prep_time, alpha=0.06, color=\"royalblue\")\n",
+    "    plt.title(rf\"Action of gate on state $|${state_id}$\\rangle$\")"
    ]
   },
   {
@@ -386,9 +409,11 @@
    "outputs": [],
    "source": [
     "# Atom Register and Device\n",
-    "q_dict = {\"control1\":np.array([-2.0, 0.]),\n",
-    "              \"target\": np.array([0., 2*np.sqrt(3.001)]),\n",
-    "              \"control2\": np.array([2.0, 0.])}\n",
+    "q_dict = {\n",
+    "    \"control1\": np.array([-2.0, 0.0]),\n",
+    "    \"target\": np.array([0.0, 2 * np.sqrt(3.001)]),\n",
+    "    \"control2\": np.array([2.0, 0.0]),\n",
+    "}\n",
     "reg = Register(q_dict)\n",
     "reg.draw()"
    ]
@@ -399,7 +424,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "preparation_sequence('hhh', reg)\n",
+    "preparation_sequence(\"hhh\", reg)\n",
     "seq.draw(draw_phase_area=True)"
    ]
   },
@@ -412,22 +437,26 @@
     "def CCZ_sequence(initial_id):\n",
     "    # Prepare State\n",
     "    prep_state = preparation_sequence(initial_id, reg)\n",
-    "    prep_time = max((seq._last(ch).tf for ch in seq.declared_channels), default=0)\n",
-    "    \n",
+    "    prep_time = max(\n",
+    "        (seq._last(ch).tf for ch in seq.declared_channels), default=0\n",
+    "    )\n",
+    "\n",
     "    # Declare Rydberg channel\n",
-    "    seq.declare_channel('ryd', 'rydberg_local', 'control1')\n",
-    "    \n",
+    "    seq.declare_channel(\"ryd\", \"rydberg_local\", \"control1\")\n",
+    "\n",
     "    # Write CCZ sequence:\n",
-    "    seq.add(pi_pulse, 'ryd', protocol='wait-for-all')  # Wait for state preparation to finish.\n",
-    "    seq.target('control2', 'ryd')\n",
-    "    seq.add(pi_pulse, 'ryd')\n",
-    "    seq.target('target','ryd')\n",
-    "    seq.add(twopi_pulse, 'ryd')\n",
-    "    seq.target('control2','ryd')\n",
-    "    seq.add(pi_pulse, 'ryd')\n",
-    "    seq.target('control1','ryd')\n",
-    "    seq.add(pi_pulse,'ryd')\n",
-    "    \n",
+    "    seq.add(\n",
+    "        pi_pulse, \"ryd\", protocol=\"wait-for-all\"\n",
+    "    )  # Wait for state preparation to finish.\n",
+    "    seq.target(\"control2\", \"ryd\")\n",
+    "    seq.add(pi_pulse, \"ryd\")\n",
+    "    seq.target(\"target\", \"ryd\")\n",
+    "    seq.add(twopi_pulse, \"ryd\")\n",
+    "    seq.target(\"control2\", \"ryd\")\n",
+    "    seq.add(pi_pulse, \"ryd\")\n",
+    "    seq.target(\"control1\", \"ryd\")\n",
+    "    seq.add(pi_pulse, \"ryd\")\n",
+    "\n",
     "    return prep_state, prep_time"
    ]
   },
@@ -437,7 +466,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "CCZ_sequence('hhh')\n",
+    "CCZ_sequence(\"hhh\")\n",
     "seq.draw(draw_phase_area=True)"
    ]
   },
@@ -450,25 +479,25 @@
    "outputs": [],
    "source": [
     "CCZ = {}\n",
-    "for state_id in {''.join(x) for x in product('gh', repeat=3)}:\n",
+    "for state_id in {\"\".join(x) for x in product(\"gh\", repeat=3)}:\n",
     "    # Get CCZ sequence\n",
     "    prep_state, prep_time = CCZ_sequence(state_id)\n",
-    "    \n",
+    "\n",
     "    # Construct Simulation instance\n",
     "    simul = Simulation(seq)\n",
-    "    \n",
+    "\n",
     "    res = simul.run()\n",
-    "    \n",
-    "    data=[st.overlap(prep_state) for st in res.states]\n",
+    "\n",
+    "    data = [st.overlap(prep_state) for st in res.states]\n",
     "    final_st = res.states[-1]\n",
     "    CCZ[state_id] = final_st.overlap(prep_state)\n",
-    "    \n",
+    "\n",
     "    plt.figure()\n",
     "    plt.plot(np.real(data))\n",
     "    plt.xlabel(r\"Time [ns]\")\n",
-    "    plt.ylabel(fr'$ \\langle\\,{state_id} | \\psi(t)\\rangle$')\n",
-    "    plt.axvspan(0, prep_time, alpha=0.06, color='royalblue')\n",
-    "    plt.title(fr\"Action of gate on state $|${state_id}$\\rangle$\")"
+    "    plt.ylabel(rf\"$ \\langle\\,{state_id} | \\psi(t)\\rangle$\")\n",
+    "    plt.axvspan(0, prep_time, alpha=0.06, color=\"royalblue\")\n",
+    "    plt.title(rf\"Action of gate on state $|${state_id}$\\rangle$\")"
    ]
   },
   {
diff --git a/tutorials/applications/Quantum Evolution Kernel.ipynb b/tutorials/applications/Quantum Evolution Kernel.ipynb
index 6d10aa0de..bc89f817e 100644
--- a/tutorials/applications/Quantum Evolution Kernel.ipynb
+++ b/tutorials/applications/Quantum Evolution Kernel.ipynb
@@ -76,26 +76,27 @@
     "\n",
     "def JSdiv(p1, p2):\n",
     "    \"\"\"Compute the Jensen-Shannon divergence between two distributions.\"\"\"\n",
-    "    q1 = np.array(p1)/np.sum(p1)\n",
-    "    q2 = np.array(p2)/np.sum(p2)\n",
+    "    q1 = np.array(p1) / np.sum(p1)\n",
+    "    q2 = np.array(p2) / np.sum(p2)\n",
     "    # Alowing for distributions to have different sizes\n",
     "    delta = len(q1) - len(q2)\n",
     "    if delta < 0:\n",
     "        q1 = np.concatenate((q1, np.array([0 for i in range(-delta)])))\n",
     "    elif delta > 0:\n",
     "        q2 = np.concatenate((q2, np.array([0 for i in range(delta)])))\n",
-    "    pq = (q1 + q2)/2\n",
+    "    pq = (q1 + q2) / 2\n",
     "\n",
     "    def entropy(pl_unscaled):\n",
     "        # Making sure the probability distributions are similarly normalized\n",
-    "        pl = np.array(pl_unscaled)/np.sum(pl_unscaled)\n",
+    "        pl = np.array(pl_unscaled) / np.sum(pl_unscaled)\n",
     "        res = 0\n",
     "        for p in pl:\n",
     "            if p > 0:\n",
-    "                res += p*np.log(p)\n",
+    "                res += p * np.log(p)\n",
     "        return -res\n",
-    "    out = entropy(pq)-(entropy(q1)+entropy(q2))/2\n",
-    "    return out\n"
+    "\n",
+    "    out = entropy(pq) - (entropy(q1) + entropy(q2)) / 2\n",
+    "    return out"
    ]
   },
   {
@@ -168,32 +169,32 @@
     "import matplotlib.pyplot as plt\n",
     "from IPython.display import Latex\n",
     "import scipy.special\n",
+    "\n",
     "# Load graph package\n",
     "import networkx as nx\n",
     "\n",
     "\n",
-    "def pk(G, theta=np.pi/4):\n",
+    "def pk(G, theta=np.pi / 4):\n",
     "    cnt = nx.degree_histogram(G)\n",
     "    kappamax = len(cnt)\n",
     "\n",
-    "    c = np.cos(theta)**2\n",
-    "    s = 1-c\n",
-    "    t = np.tan(theta)**2\n",
+    "    c = np.cos(theta) ** 2\n",
+    "    s = 1 - c\n",
+    "    t = np.tan(theta) ** 2\n",
     "    sp = 2 * c * s\n",
     "\n",
     "    res0 = 0\n",
     "    for kappa, m in enumerate(cnt):\n",
-    "        res0 += m * (1-c**kappa)\n",
+    "        res0 += m * (1 - c**kappa)\n",
     "    res = [(sp * res0)]\n",
     "    for k in range(1, kappamax):\n",
     "        res0 = 0\n",
     "        for kappa in range(k, kappamax):\n",
     "            m_kappa = cnt[kappa]\n",
     "            binom = scipy.special.comb(kappa, k, exact=True)\n",
-    "            res0 += m_kappa * binom * (c**(kappa+1-k))\n",
-    "        res.append(((s**(1+k)) * res0))\n",
-    "    return res\n",
-    "    "
+    "            res0 += m_kappa * binom * (c ** (kappa + 1 - k))\n",
+    "        res.append(((s ** (1 + k)) * res0))\n",
+    "    return res"
    ]
   },
   {
@@ -217,10 +218,7 @@
     "n_graphs = 100\n",
     "\n",
     "\n",
-    "def create_random_graphs(N_max=100,\n",
-    "                         n_graphs=100,\n",
-    "                         rho_low=0.35,\n",
-    "                         rho_high=0.65):\n",
+    "def create_random_graphs(N_max=100, n_graphs=100, rho_low=0.35, rho_high=0.65):\n",
     "    # Dataset with graphs of two different Erdős–Rényi classes\n",
     "    graphs = []\n",
     "    # Classes of these graphs\n",
@@ -229,8 +227,8 @@
     "    probability_distributions = []\n",
     "    for _ in range(n_graphs):\n",
     "        # Number of nodes in the graph in [N_max/2,N_max]\n",
-    "        N = np.random.randint(N_max//2, N_max+1)\n",
-    "        if np.random.rand() < .5:\n",
+    "        N = np.random.randint(N_max // 2, N_max + 1)\n",
+    "        if np.random.rand() < 0.5:\n",
     "            rho = rho_low\n",
     "            classes.append(0)\n",
     "        else:\n",
@@ -239,9 +237,9 @@
     "        G = nx.erdos_renyi_graph(N, rho)\n",
     "        graphs.append(G)\n",
     "        pdist = pk(G)\n",
-    "        probability_distributions.append(pdist/np.sum(pdist))\n",
+    "        probability_distributions.append(pdist / np.sum(pdist))\n",
     "\n",
-    "    return graphs, classes, probability_distributions\n"
+    "    return graphs, classes, probability_distributions"
    ]
   },
   {
@@ -266,14 +264,15 @@
    "outputs": [],
    "source": [
     "def kernel_matrix(pdist1, pdist2, mu=1):\n",
-    "    Kmat = np.array([[np.exp(-mu * JSdiv(p1, p2)) for p1 in pdist1]\n",
-    "                     for p2 in pdist2])\n",
+    "    Kmat = np.array(\n",
+    "        [[np.exp(-mu * JSdiv(p1, p2)) for p1 in pdist1] for p2 in pdist2]\n",
+    "    )\n",
     "    return Kmat\n",
     "\n",
     "\n",
     "graphs, classes, proba_dists = create_random_graphs()\n",
     "\n",
-    "Kmat = kernel_matrix(proba_dists, proba_dists)\n"
+    "Kmat = kernel_matrix(proba_dists, proba_dists)"
    ]
   },
   {
@@ -298,13 +297,13 @@
    "source": [
     "def plot_kernel_matrix(Kmat):\n",
     "    fig, ax = plt.subplots(figsize=(8, 8))\n",
-    "    im = ax.imshow(Kmat, cmap='OrRd')\n",
-    "    ax.set_xlabel('Graph #', fontsize=18)\n",
-    "    ax.set_ylabel('Graph #', fontsize=18)\n",
-    "    cbar = plt.colorbar(im, extend='max')\n",
+    "    im = ax.imshow(Kmat, cmap=\"OrRd\")\n",
+    "    ax.set_xlabel(\"Graph #\", fontsize=18)\n",
+    "    ax.set_ylabel(\"Graph #\", fontsize=18)\n",
+    "    cbar = plt.colorbar(im, extend=\"max\")\n",
     "\n",
     "\n",
-    "plot_kernel_matrix(Kmat)\n"
+    "plot_kernel_matrix(Kmat)"
    ]
   },
   {
@@ -327,13 +326,15 @@
    "outputs": [],
    "source": [
     "from sklearn import svm\n",
-    "from sklearn.metrics import f1_score, accuracy_score, recall_score, precision_score\n",
+    "from sklearn.metrics import (\n",
+    "    f1_score,\n",
+    "    accuracy_score,\n",
+    "    recall_score,\n",
+    "    precision_score,\n",
+    ")\n",
     "\n",
     "\n",
-    "scores_types = ['Accuracy ',\n",
-    "                'f1       ',\n",
-    "                'Precision',\n",
-    "                'Recall   ']\n",
+    "scores_types = [\"Accuracy \", \"f1       \", \"Precision\", \"Recall   \"]\n",
     "\n",
     "\n",
     "def trained_classifier_from_Kmat(Kmat, classes_train):\n",
@@ -341,7 +342,7 @@
     "    Create and train a classifier from the Kernel matrix `Kmat`\n",
     "    obtained from graphs of classes `classes_train`\n",
     "    \"\"\"\n",
-    "    classifier = svm.SVC(kernel='precomputed')\n",
+    "    classifier = svm.SVC(kernel=\"precomputed\")\n",
     "    classifier.fit(Kmat, classes_train)\n",
     "\n",
     "    return classifier\n",
@@ -357,11 +358,10 @@
     "\n",
     "    return trained_classifier_from_Kmat(Kmat, classes_train)\n",
     "\n",
-    "def test_classifier(classifier,\n",
-    "                    p_dist_train,\n",
-    "                    p_dist_test,\n",
-    "                    classes_test,\n",
-    "                    verbose=False):\n",
+    "\n",
+    "def test_classifier(\n",
+    "    classifier, p_dist_train, p_dist_test, classes_test, verbose=False\n",
+    "):\n",
     "    \"\"\"\n",
     "    Test a trained classifier `classifier` from the probability\n",
     "    distributions of the train and test data sets `p_dist_train`\n",
@@ -372,40 +372,37 @@
     "\n",
     "    predicted_classes = classifier.predict(X)\n",
     "\n",
-    "    scores = [accuracy_score(classes_test,\n",
-    "                             predicted_classes),\n",
-    "              f1_score(classes_test,\n",
-    "                       predicted_classes,\n",
-    "                       average='weighted'),\n",
-    "              precision_score(classes_test,\n",
-    "                              predicted_classes,\n",
-    "                              average='weighted',\n",
-    "                              zero_division=0),\n",
-    "              recall_score(classes_test,\n",
-    "                           predicted_classes,\n",
-    "                           average='weighted')]\n",
+    "    scores = [\n",
+    "        accuracy_score(classes_test, predicted_classes),\n",
+    "        f1_score(classes_test, predicted_classes, average=\"weighted\"),\n",
+    "        precision_score(\n",
+    "            classes_test,\n",
+    "            predicted_classes,\n",
+    "            average=\"weighted\",\n",
+    "            zero_division=0,\n",
+    "        ),\n",
+    "        recall_score(classes_test, predicted_classes, average=\"weighted\"),\n",
+    "    ]\n",
     "\n",
     "    if verbose:\n",
     "        for st, s in zip(scores_types, scores):\n",
-    "            print(f'{st} : {s:6.3}')\n",
+    "            print(f\"{st} : {s:6.3}\")\n",
     "\n",
     "    return scores\n",
     "\n",
     "\n",
-    "def train_and_test_classifier(p_dist_train,\n",
-    "                              classes_train,\n",
-    "                              p_dist_test,\n",
-    "                              classes_test,\n",
-    "                              verbose=False):\n",
+    "def train_and_test_classifier(\n",
+    "    p_dist_train, classes_train, p_dist_test, classes_test, verbose=False\n",
+    "):\n",
     "    \"\"\"\n",
     "    Train and test a classifier from test and\n",
     "    train probability distributions and classes\n",
     "    \"\"\"\n",
     "    classifier = trained_classifier_pdist(p_dist_train, classes_train)\n",
     "\n",
-    "    return classifier, test_classifier(classifier, p_dist_train,\n",
-    "                                       p_dist_test, classes_test,\n",
-    "                                       verbose=verbose)\n"
+    "    return classifier, test_classifier(\n",
+    "        classifier, p_dist_train, p_dist_test, classes_test, verbose=verbose\n",
+    "    )"
    ]
   },
   {
@@ -440,12 +437,9 @@
     "graphs_test, classes_test, p_dist_test = create_random_graphs(n_graphs=50)\n",
     "\n",
     "# Compute the score of the classifier\n",
-    "classifier, scores = train_and_test_classifier(p_dist_train,\n",
-    "                                               classes_train,\n",
-    "                                               p_dist_test,\n",
-    "                                               classes_test,\n",
-    "                                               verbose=True\n",
-    "                                               )\n"
+    "classifier, scores = train_and_test_classifier(\n",
+    "    p_dist_train, classes_train, p_dist_test, classes_test, verbose=True\n",
+    ")"
    ]
   },
   {
@@ -467,51 +461,51 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "prefix = './Fingerprint/Fingerprint_'\n",
+    "prefix = \"./Fingerprint/Fingerprint_\"\n",
     "\n",
     "graphs = {}\n",
     "node_to_graph = {}\n",
-    "class_count  = {}\n",
+    "class_count = {}\n",
     "\n",
-    "label_file = prefix + 'graph_labels' + '.txt'\n",
+    "label_file = prefix + \"graph_labels\" + \".txt\"\n",
     "with open(label_file) as f:\n",
     "    lines = f.readlines()\n",
     "    for i, line in enumerate(lines):\n",
     "        labl = int(line)\n",
-    "        graphs[i+1] = nx.Graph(label=labl)\n",
+    "        graphs[i + 1] = nx.Graph(label=labl)\n",
     "        if labl in class_count.keys():\n",
     "            class_count[labl] += 1\n",
     "        else:\n",
     "            class_count[labl] = 1\n",
     "\n",
     "\n",
-    "node_to_graph_file = prefix + 'graph_indicator' + '.txt'\n",
+    "node_to_graph_file = prefix + \"graph_indicator\" + \".txt\"\n",
     "with open(node_to_graph_file) as f:\n",
     "    lines = f.readlines()\n",
     "    for i, line in enumerate(lines):\n",
     "        gi = int(line)\n",
-    "        node_to_graph[i+1] = gi\n",
-    "        graphs[gi].add_node(i+1)\n",
+    "        node_to_graph[i + 1] = gi\n",
+    "        graphs[gi].add_node(i + 1)\n",
     "\n",
-    "adjacency_file = prefix + 'A' + '.txt'\n",
+    "adjacency_file = prefix + \"A\" + \".txt\"\n",
     "with open(adjacency_file) as f:\n",
     "    lines = f.readlines()\n",
     "    for line in lines:\n",
-    "        Ind = line.split(',')\n",
+    "        Ind = line.split(\",\")\n",
     "        i = int(Ind[0])\n",
     "        j = int(Ind[1])\n",
     "        gi = node_to_graph[i]\n",
     "        graphs[gi].add_edge(i, j)\n",
     "\n",
-    "coordinates_file = prefix + 'node_attributes' + '.txt'\n",
+    "coordinates_file = prefix + \"node_attributes\" + \".txt\"\n",
     "with open(coordinates_file) as f:\n",
     "    lines = f.readlines()\n",
-    "    for i,line in enumerate(lines):\n",
-    "        Ind = line.split(',')\n",
+    "    for i, line in enumerate(lines):\n",
+    "        Ind = line.split(\",\")\n",
     "        x = float(Ind[0])\n",
     "        y = float(Ind[1])\n",
-    "        gi = node_to_graph[i+1]\n",
-    "        nx.set_node_attributes(graphs[gi], {i+1: (x,y)}, \"coords\")"
+    "        gi = node_to_graph[i + 1]\n",
+    "        nx.set_node_attributes(graphs[gi], {i + 1: (x, y)}, \"coords\")"
    ]
   },
   {
@@ -556,20 +550,20 @@
     "count = {clas: 0 for clas in class_count.keys()}\n",
     "for g in graphs.values():\n",
     "    if Nmin <= g.number_of_nodes() <= Nmax:\n",
-    "        count[g.graph['label']] += 1\n",
+    "        count[g.graph[\"label\"]] += 1\n",
     "\n",
     "# Number of graphs in the most represented class\n",
     "size_of_largest_class = max(count.values())\n",
     "# Include only classes with at least 10% of the size of the largest one\n",
     "include_classes = {clas: False for clas in class_count.keys()}\n",
     "for clas, prop in count.items():\n",
-    "    if prop > .1*size_of_largest_class:\n",
+    "    if prop > 0.1 * size_of_largest_class:\n",
     "        include_classes[clas] = True\n",
     "\n",
     "\n",
     "data_preprocessed = []\n",
     "for g in graphs.values():\n",
-    "    labl = g.graph['label']\n",
+    "    labl = g.graph[\"label\"]\n",
     "    if Nmin <= g.number_of_nodes() <= Nmax and include_classes[labl]:\n",
     "        mapping = {l: i for i, l in enumerate(g.nodes())}\n",
     "        g_shift = nx.relabel_nodes(g, mapping)\n",
@@ -584,10 +578,11 @@
     "    if include_classes[clas]:\n",
     "        included_classes[clas] = icount\n",
     "\n",
-    "print(f'After preprocessing, the dataset now contains {len(data_preprocessed)} \\n' +\n",
-    "      f'graphs of at least {Nmin} and at most {Nmax} nodes, distributed \\n' +\n",
-    "      f'across the different classes in the following way {included_classes}\\n')\n",
-    "\n"
+    "print(\n",
+    "    f\"After preprocessing, the dataset now contains {len(data_preprocessed)} \\n\"\n",
+    "    + f\"graphs of at least {Nmin} and at most {Nmax} nodes, distributed \\n\"\n",
+    "    + f\"across the different classes in the following way {included_classes}\\n\"\n",
+    ")"
    ]
   },
   {
@@ -623,8 +618,8 @@
     "\n",
     "data_reduced = []\n",
     "for g in data_preprocessed:\n",
-    "    if(len(data_reduced) < dataset_targetsize):\n",
-    "        cls = g.graph['label']\n",
+    "    if len(data_reduced) < dataset_targetsize:\n",
+    "        cls = g.graph[\"label\"]\n",
     "        if cls in kept_classes.keys():\n",
     "            kept_classes[cls] += 1\n",
     "            data_reduced.append(g)\n",
@@ -632,10 +627,12 @@
     "# size of the dataset\n",
     "n_graphs = len(data_reduced)\n",
     "\n",
-    "print(f'After preprocessing, the dataset now contains {len(data_reduced)} ' +\n",
-    "      f'graphs of at least {Nmin} and at most {Nmax} nodes, distributed ' +\n",
-    "      f'across {len(kept_classes)} different classes in the following way ' +\n",
-    "      f'{kept_classes}')\n"
+    "print(\n",
+    "    f\"After preprocessing, the dataset now contains {len(data_reduced)} \"\n",
+    "    + f\"graphs of at least {Nmin} and at most {Nmax} nodes, distributed \"\n",
+    "    + f\"across {len(kept_classes)} different classes in the following way \"\n",
+    "    + f\"{kept_classes}\"\n",
+    ")"
    ]
   },
   {
@@ -672,10 +669,10 @@
     "\n",
     "\n",
     "def correct_coordinates(g):\n",
-    "    '''\n",
+    "    \"\"\"\n",
     "    Corrects the coordinates of the nodes so that the\n",
     "    atoms fit the hardware constraints.\n",
-    "    '''\n",
+    "    \"\"\"\n",
     "    n = g.number_of_nodes()\n",
     "\n",
     "    # Coordinates given in the dataset\n",
@@ -683,47 +680,50 @@
     "    r_list += -np.mean(r_list, axis=0)\n",
     "    scale = np.max([np.sqrt(r.dot(r)) for r in r_list])\n",
     "\n",
-    "    x0 = r_list.reshape(2*n)\n",
+    "    x0 = r_list.reshape(2 * n)\n",
     "\n",
     "    # Ensures the atoms are within range of the device\n",
-    "    xmax = device.max_radial_distance/np.sqrt(2)\n",
-    "    bounds = [(-xmax, xmax)] * (2*n)\n",
-    "    x0 *= xmax/scale\n",
+    "    xmax = device.max_radial_distance / np.sqrt(2)\n",
+    "    bounds = [(-xmax, xmax)] * (2 * n)\n",
+    "    x0 *= xmax / scale\n",
     "\n",
-    "    #Encode the constraint of a minimal distance bewteen atoms\n",
+    "    # Encode the constraint of a minimal distance bewteen atoms\n",
     "    def min_dist(params):\n",
     "        return np.min(pdist(params.reshape(n, 2)))\n",
     "\n",
-    "    dmin = 1.1*device.min_atom_distance\n",
+    "    dmin = 1.1 * device.min_atom_distance\n",
     "    nlc = NonlinearConstraint(min_dist, dmin, np.inf)\n",
     "\n",
     "    def cost_function(params):\n",
     "        return 1\n",
     "\n",
-    "    res = minimize(cost_function,\n",
-    "                   x0=x0,\n",
-    "                   bounds=bounds,\n",
-    "                   constraints=nlc,\n",
-    "                   method='SLSQP')\n",
+    "    res = minimize(\n",
+    "        cost_function, x0=x0, bounds=bounds, constraints=nlc, method=\"SLSQP\"\n",
+    "    )\n",
     "    x = res.x\n",
     "    rmax = device.max_radial_distance\n",
-    "    scale_diameter = .95 * rmax/np.max(pdist(x.reshape(n, 2))) \n",
-    "    x *= max(scale_diameter,1.)\n",
+    "    scale_diameter = 0.95 * rmax / np.max(pdist(x.reshape(n, 2)))\n",
+    "    x *= max(scale_diameter, 1.0)\n",
     "    r_list = x.reshape(n, 2)\n",
     "    r_list += -np.mean(r_list, axis=0)\n",
     "\n",
-    "\n",
     "    for node, r in zip(g.nodes(), r_list):\n",
     "        g.nodes[node][\"coords\"] = r\n",
     "\n",
     "\n",
     "def max_edge_length(g):\n",
-    "    '''\n",
+    "    \"\"\"\n",
     "    Computes the maximal distance between nodes connected by an edge\n",
     "    of the graph\n",
-    "    '''\n",
+    "    \"\"\"\n",
     "    n = g.number_of_nodes()\n",
-    "    edges = np.array([1 if (i, j) in g.edges() else 0  for i in range(n) for j in range(i+1, n)])\n",
+    "    edges = np.array(\n",
+    "        [\n",
+    "            1 if (i, j) in g.edges() else 0\n",
+    "            for i in range(n)\n",
+    "            for j in range(i + 1, n)\n",
+    "        ]\n",
+    "    )\n",
     "    r_list = np.array([g.nodes[node][\"coords\"] for node in g.nodes()])\n",
     "\n",
     "    distances = pdist(r_list)\n",
@@ -746,9 +746,10 @@
     "r_max = device.max_radial_distance\n",
     "d_min = device.min_atom_distance\n",
     "\n",
-    "omega_max = device.channels['rydberg_global'].max_amp\n",
+    "omega_max = device.channels[\"rydberg_global\"].max_amp\n",
     "min_bond_length = device.rabi_from_blockade(omega_max)\n",
     "\n",
+    "\n",
     "def reg_from_data(data_reduced):\n",
     "    # The list of registers for each graph\n",
     "    reg_list = []\n",
@@ -764,10 +765,10 @@
     "    label_list = []\n",
     "    for g in data_reduced:\n",
     "\n",
-    "        label_list.append(g.graph['label'])\n",
+    "        label_list.append(g.graph[\"label\"])\n",
     "\n",
     "        correct_coordinates(g)\n",
-    "        graph_dict= {i:g.nodes[i][\"coords\"] for i in g.nodes()}\n",
+    "        graph_dict = {i: g.nodes[i][\"coords\"] for i in g.nodes()}\n",
     "\n",
     "        edges_list.append(g.edges)\n",
     "\n",
@@ -782,6 +783,7 @@
     "\n",
     "    return reg_list, rabi_list, edges_list, label_list\n",
     "\n",
+    "\n",
     "reg_list, rabi_list, edges_list, label_list = reg_from_data(data_reduced)"
    ]
   },
@@ -854,65 +856,56 @@
    "source": [
     "from pulser import Pulse, Sequence, Simulation\n",
     "\n",
-    "def pulse_seqence(reg,\n",
-    "                  t_1=100,\n",
-    "                  omega=omega_max, # amplitude of the initial and final pulses\n",
-    "                  omega_g=0, # amplitude in the \"free evolution\" parts\n",
-    "                  total_time=512): # total duration of the pulse\n",
+    "\n",
+    "def pulse_seqence(\n",
+    "    reg,\n",
+    "    t_1=100,\n",
+    "    omega=omega_max,  # amplitude of the initial and final pulses\n",
+    "    omega_g=0,  # amplitude in the \"free evolution\" parts\n",
+    "    total_time=512,\n",
+    "):  # total duration of the pulse\n",
     "    seq = Sequence(reg, device)\n",
-    "    seq.declare_channel('Channel 0','rydberg_global')\n",
+    "    seq.declare_channel(\"Channel 0\", \"rydberg_global\")\n",
     "\n",
-    "    # making sure that the value of omega does not exceed the \n",
-    "    # maximal value, and that it doesn't lead to a pulse \n",
+    "    # making sure that the value of omega does not exceed the\n",
+    "    # maximal value, and that it doesn't lead to a pulse\n",
     "    # duration that is too short\n",
-    "    omega = min([omega,1000*np.pi/2,omega_max])\n",
+    "    omega = min([omega, 1000 * np.pi / 2, omega_max])\n",
     "\n",
     "    # Set the initial and final pulse times to the optimal value\n",
     "    # be careful about the units : Omega(rad/μs) -> t (ns)\n",
-    "    t = 1000*np.pi/(2*omega)\n",
+    "    t = 1000 * np.pi / (2 * omega)\n",
     "    # Set the total_time\n",
-    "    tau = (total_time - 2*t - t_1)/2\n",
+    "    tau = (total_time - 2 * t - t_1) / 2\n",
     "    # No detuning needed here\n",
     "    delta = 0\n",
     "    # We want the pulse to be along sigma_y\n",
-    "    phi=np.pi/2\n",
-    "\n",
-    "    initial_pulse = Pulse.ConstantPulse(t,\n",
-    "                                        omega,\n",
-    "                                        delta,\n",
-    "                                        phase=phi)\n",
-    "    if total_time > t_1 + 2*t:\n",
-    "        Hg_pulse = Pulse.ConstantPulse(tau,\n",
-    "                                       omega_g,\n",
-    "                                       delta,\n",
-    "                                       phase=phi)\n",
+    "    phi = np.pi / 2\n",
+    "\n",
+    "    initial_pulse = Pulse.ConstantPulse(t, omega, delta, phase=phi)\n",
+    "    if total_time > t_1 + 2 * t:\n",
+    "        Hg_pulse = Pulse.ConstantPulse(tau, omega_g, delta, phase=phi)\n",
     "    if t_1 > 0:\n",
-    "        middle_pulse = Pulse.ConstantPulse(t_1,\n",
-    "                                           omega,\n",
-    "                                           delta,\n",
-    "                                           phase=phi)\n",
-    "    final_pulse = Pulse.ConstantPulse(t,\n",
-    "                                      omega,\n",
-    "                                      delta,\n",
-    "                                      phase=phi)\n",
-    "\n",
-    "    seq.add(initial_pulse, 'Channel 0')\n",
-    "    if total_time > t_1 + 2*t:\n",
-    "        seq.add(Hg_pulse, 'Channel 0')\n",
+    "        middle_pulse = Pulse.ConstantPulse(t_1, omega, delta, phase=phi)\n",
+    "    final_pulse = Pulse.ConstantPulse(t, omega, delta, phase=phi)\n",
+    "\n",
+    "    seq.add(initial_pulse, \"Channel 0\")\n",
+    "    if total_time > t_1 + 2 * t:\n",
+    "        seq.add(Hg_pulse, \"Channel 0\")\n",
     "    if t_1 > 0:\n",
-    "        seq.add(middle_pulse, 'Channel 0')\n",
-    "    if total_time > t_1 + 2*t:\n",
-    "        seq.add(Hg_pulse, 'Channel 0')\n",
-    "    seq.add(final_pulse, 'Channel 0')\n",
+    "        seq.add(middle_pulse, \"Channel 0\")\n",
+    "    if total_time > t_1 + 2 * t:\n",
+    "        seq.add(Hg_pulse, \"Channel 0\")\n",
+    "    seq.add(final_pulse, \"Channel 0\")\n",
     "\n",
-    "    seq.measure(basis='ground-rydberg')\n",
+    "    seq.measure(basis=\"ground-rydberg\")\n",
     "\n",
     "    return seq\n",
     "\n",
     "\n",
     "# Illustrate the pulse on a register containing a single atom\n",
     "reg = Register.from_coordinates([(0, 0)])\n",
-    "pulse_seqence(reg, t_1=160).draw()\n"
+    "pulse_seqence(reg, t_1=160).draw()"
    ]
   },
   {
@@ -933,30 +926,34 @@
     "from tqdm.auto import tqdm\n",
     "\n",
     "\n",
-    "def proba_distributions(t_1=100,  # duration of the central pulse\n",
-    "                        omega=omega_max,  # amplitude of the pulses\n",
-    "                        omega_g_factor=1,  # set to 1 if the Amplitude is non\n",
-    "                                           # zero during the \"free evolution\"\n",
-    "                        total_time=512,  # total duration of the pulse\n",
-    "                        Nsamples=1000,\n",
-    "                        indices=list(range(n_graphs))):  # graphs to be used\n",
-    "    '''\n",
+    "def proba_distributions(\n",
+    "    t_1=100,  # duration of the central pulse\n",
+    "    omega=omega_max,  # amplitude of the pulses\n",
+    "    omega_g_factor=1,  # set to 1 if the Amplitude is non\n",
+    "    # zero during the \"free evolution\"\n",
+    "    total_time=512,  # total duration of the pulse\n",
+    "    Nsamples=1000,\n",
+    "    indices=list(range(n_graphs)),\n",
+    "):  # graphs to be used\n",
+    "    \"\"\"\n",
     "    Compute the probability distributions for a given pulse\n",
     "    for all graphs in `indices`\n",
-    "    '''\n",
+    "    \"\"\"\n",
     "\n",
-    "    bins = np.linspace(0, Nmax*Nmax, Nmax*Nmax + 1)\n",
+    "    bins = np.linspace(0, Nmax * Nmax, Nmax * Nmax + 1)\n",
     "    histograms = []\n",
     "    for i in tqdm(indices):\n",
     "        reg, rabi, edges = reg_list[i], rabi_list[i], edges_list[i]\n",
-    "        seq = pulse_seqence(reg,\n",
-    "                            t_1=t_1,\n",
-    "                            omega=omega,\n",
-    "                            omega_g=omega_g_factor*rabi,\n",
-    "                            total_time=total_time)\n",
+    "        seq = pulse_seqence(\n",
+    "            reg,\n",
+    "            t_1=t_1,\n",
+    "            omega=omega,\n",
+    "            omega_g=omega_g_factor * rabi,\n",
+    "            total_time=total_time,\n",
+    "        )\n",
     "\n",
     "        # Simulate and sample\n",
-    "        simul = Simulation(seq, evaluation_times=.5, sampling_rate=.1)\n",
+    "        simul = Simulation(seq, evaluation_times=0.5, sampling_rate=0.1)\n",
     "        results = simul.run()\n",
     "        sampling = results.sample_final_state(N_samples=Nsamples)\n",
     "\n",
@@ -968,19 +965,18 @@
     "            ie_weights.append(num)\n",
     "\n",
     "        # Create histogram of the measurements and append to list\n",
-    "        ncount, b = np.histogram(ie_meas,\n",
-    "                                 bins=bins,\n",
-    "                                 density=True,\n",
-    "                                 weights=ie_weights)\n",
+    "        ncount, b = np.histogram(\n",
+    "            ie_meas, bins=bins, density=True, weights=ie_weights\n",
+    "        )\n",
     "        histograms.append(ncount)\n",
     "    return histograms\n",
     "\n",
     "\n",
     "def compute_ising_energy(outcome, edges):\n",
-    "    '''\n",
+    "    \"\"\"\n",
     "    Computes the Ising energy (i.e. the observable\n",
     "    used by the kernel) from a measure bitstgring/state\n",
-    "    '''\n",
+    "    \"\"\"\n",
     "    # split outcome string in a list\n",
     "    outcome_ls = [char for char in outcome]\n",
     "\n",
@@ -990,9 +986,9 @@
     "        i = int(edge[0])\n",
     "        j = int(edge[1])\n",
     "        if i < j:\n",
-    "            energy += int(outcome_ls[i])*int(outcome_ls[j])\n",
+    "            energy += int(outcome_ls[i]) * int(outcome_ls[j])\n",
     "\n",
-    "    return energy\n"
+    "    return energy"
    ]
   },
   {
@@ -1053,14 +1049,14 @@
    "source": [
     "n_graphs = len(data_reduced)\n",
     "\n",
-    "#sample 150 graphs and train on 100 of them\n",
+    "# sample 150 graphs and train on 100 of them\n",
     "n_train = 150\n",
     "n_test = 50\n",
     "\n",
     "# randomize graph order\n",
     "indices_all = list(range(n_graphs))\n",
     "indices_train = indices_all[0:n_train]\n",
-    "indices_test = indices_all[n_train:n_train+n_test]\n",
+    "indices_test = indices_all[n_train : n_train + n_test]\n",
     "\n",
     "# Labels of the sampled graphs\n",
     "train_classes = np.array([label_list[i] for i in indices_train])\n",
@@ -1068,17 +1064,17 @@
     "\n",
     "\n",
     "# Probability distributions obtained after the pulse\n",
-    "print('Training in progress...')\n",
+    "print(\"Training in progress...\")\n",
     "probas_train = proba_distributions(t_1=0, indices=indices_train)\n",
-    "print('Testing in progress...')\n",
+    "print(\"Testing in progress...\")\n",
     "probas_test = proba_distributions(t_1=0, indices=indices_test)\n",
     "\n",
     "# Resulting kernel matrix\n",
     "Kmat = kernel_matrix(probas_train, probas_train)\n",
     "\n",
     "fig, ax = plt.subplots(figsize=(8, 8))\n",
-    "im = ax.imshow(Kmat, cmap='OrRd')\n",
-    "cbar = plt.colorbar(im, extend='max')\n"
+    "im = ax.imshow(Kmat, cmap=\"OrRd\")\n",
+    "cbar = plt.colorbar(im, extend=\"max\")"
    ]
   },
   {
@@ -1099,11 +1095,9 @@
     }
    ],
    "source": [
-    "classifier, scores = train_and_test_classifier(probas_train,\n",
-    "                                               train_classes,\n",
-    "                                               probas_test,\n",
-    "                                               test_classes,\n",
-    "                                               verbose=True)\n"
+    "classifier, scores = train_and_test_classifier(\n",
+    "    probas_train, train_classes, probas_test, test_classes, verbose=True\n",
+    ")"
    ]
   },
   {
@@ -1138,15 +1132,17 @@
     "M = 1\n",
     "\n",
     "\n",
-    "def score_function(t_1=100,\n",
-    "                   total_time=512,\n",
-    "                   repetitions=M,\n",
-    "                   nblocks=N,\n",
-    "                   label_list=label_list,\n",
-    "                   indices=list(range(n_graphs))):  # list of graphs included\n",
-    "    '''\n",
+    "def score_function(\n",
+    "    t_1=100,\n",
+    "    total_time=512,\n",
+    "    repetitions=M,\n",
+    "    nblocks=N,\n",
+    "    label_list=label_list,\n",
+    "    indices=list(range(n_graphs)),\n",
+    "):  # list of graphs included\n",
+    "    \"\"\"\n",
     "    Computes the accuracy, f1, precision and recall\n",
-    "    '''\n",
+    "    \"\"\"\n",
     "\n",
     "    accuracy = []\n",
     "    f1 = []\n",
@@ -1155,65 +1151,74 @@
     "\n",
     "    n_g = len(indices)\n",
     "\n",
-    "    block_size = n_g//nblocks\n",
-    "\n",
+    "    block_size = n_g // nblocks\n",
     "\n",
-    "    # Compute the probability distributions of all \n",
+    "    # Compute the probability distributions of all\n",
     "    # graphs in the data set\n",
     "    start_time = time.time()\n",
-    "    probas_all = proba_distributions(t_1=t_1,\n",
-    "                                     total_time=total_time,\n",
-    "                                     Nsamples=1000,\n",
-    "                                     indices=indices)\n",
+    "    probas_all = proba_distributions(\n",
+    "        t_1=t_1, total_time=total_time, Nsamples=1000, indices=indices\n",
+    "    )\n",
     "\n",
-    "    print(f'  Probability lists were computed in {time.time() - start_time:4.1f} seconds')\n",
+    "    print(\n",
+    "        f\"  Probability lists were computed in {time.time() - start_time:4.1f} seconds\"\n",
+    "    )\n",
     "\n",
     "    classes = np.array([label_list[i] for i in indices])\n",
     "    start_time = time.time()\n",
     "\n",
     "    for r in range(repetitions):\n",
-    "        #divide data in training set and test set\n",
+    "        # divide data in training set and test set\n",
     "        indices_all = np.array(list(range(n_g)))\n",
     "        np.random.shuffle(indices_all)\n",
     "\n",
-    "        mean_scores = np.zeros((4, ))\n",
+    "        mean_scores = np.zeros((4,))\n",
     "        for iblock in range(nblocks):\n",
-    "            indices_test = [indices_all[(iblock * block_size + i) % n_g]\n",
-    "                            for i in range(block_size)]\n",
-    "            indices_train = [indices_all[((iblock + 1) * block_size + i) % n_g]\n",
-    "                             for i in range(n_g - block_size)]\n",
-    "\n",
-    "            train_classes = np.array([label_list[indices[i]]\n",
-    "                                      for i in indices_train])\n",
-    "            test_classes = np.array([label_list[indices[i]]\n",
-    "                                     for i in indices_test])\n",
-    "\n",
-    "            # create probability histogram for train and test data \n",
+    "            indices_test = [\n",
+    "                indices_all[(iblock * block_size + i) % n_g]\n",
+    "                for i in range(block_size)\n",
+    "            ]\n",
+    "            indices_train = [\n",
+    "                indices_all[((iblock + 1) * block_size + i) % n_g]\n",
+    "                for i in range(n_g - block_size)\n",
+    "            ]\n",
+    "\n",
+    "            train_classes = np.array(\n",
+    "                [label_list[indices[i]] for i in indices_train]\n",
+    "            )\n",
+    "            test_classes = np.array(\n",
+    "                [label_list[indices[i]] for i in indices_test]\n",
+    "            )\n",
+    "\n",
+    "            # create probability histogram for train and test data\n",
     "            probas_train = np.array([probas_all[i] for i in indices_train])\n",
     "            probas_test = np.array([probas_all[i] for i in indices_test])\n",
     "\n",
-    "            classifier, scores = train_and_test_classifier(probas_train,\n",
-    "                                                           train_classes,\n",
-    "                                                           probas_test,\n",
-    "                                                           test_classes,\n",
-    "                                                           verbose=False)\n",
+    "            classifier, scores = train_and_test_classifier(\n",
+    "                probas_train,\n",
+    "                train_classes,\n",
+    "                probas_test,\n",
+    "                test_classes,\n",
+    "                verbose=False,\n",
+    "            )\n",
     "            mean_scores += scores\n",
     "\n",
     "        # calculate score metrics\n",
-    "        accuracy.append(mean_scores[0]/nblocks)\n",
-    "        f1.append(mean_scores[1]/nblocks)\n",
-    "        precision.append(mean_scores[2]/nblocks)\n",
-    "        recall.append(mean_scores[3]/nblocks)\n",
+    "        accuracy.append(mean_scores[0] / nblocks)\n",
+    "        f1.append(mean_scores[1] / nblocks)\n",
+    "        precision.append(mean_scores[2] / nblocks)\n",
+    "        recall.append(mean_scores[3] / nblocks)\n",
     "\n",
     "    A = (np.mean(accuracy), np.std(accuracy))\n",
     "    B = (np.mean(f1), np.std(f1))\n",
     "    C = (np.mean(precision), np.std(precision))\n",
     "    D = (np.mean(recall), np.std(recall))\n",
     "\n",
-    "    print(f'  Kernel scores computed in {time.time() - start_time:4.1f} seconds')\n",
+    "    print(\n",
+    "        f\"  Kernel scores computed in {time.time() - start_time:4.1f} seconds\"\n",
+    "    )\n",
     "\n",
-    "    return A, B, C, D\n",
-    "    "
+    "    return A, B, C, D"
    ]
   },
   {
@@ -1235,49 +1240,58 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def scan_scores(M=2,\n",
-    "                N=4,\n",
-    "                indices=list(range(n_graphs)),\n",
-    "                durations=[512],\n",
-    "                ):\n",
+    "def scan_scores(\n",
+    "    M=2,\n",
+    "    N=4,\n",
+    "    indices=list(range(n_graphs)),\n",
+    "    durations=[512],\n",
+    "):\n",
     "\n",
     "    scores_dict = {}\n",
     "\n",
     "    for s in scores_types:\n",
     "        scores_dict[s] = []\n",
     "\n",
-    "    print(' ------------------------------------------------')\n",
-    "    print(f'|   Max. duration of the middle pulse: {durations[-1]:4d} ns   |')\n",
-    "    print(f'|   Total duration of the pulse:       {total_time:4d} ns   |')\n",
-    "    print(f'|   Using {N:2d} blocks of {len(indices)//N:4d} graphs each          |')\n",
-    "    print(' ------------------------------------------------')\n",
+    "    print(\" ------------------------------------------------\")\n",
+    "    print(f\"|   Max. duration of the middle pulse: {durations[-1]:4d} ns   |\")\n",
+    "    print(f\"|   Total duration of the pulse:       {total_time:4d} ns   |\")\n",
+    "    print(\n",
+    "        f\"|   Using {N:2d} blocks of {len(indices)//N:4d} graphs each          |\"\n",
+    "    )\n",
+    "    print(\" ------------------------------------------------\")\n",
     "\n",
     "    for t_1 in durations:\n",
-    "        print(f'  Duration of the middle pulse: {t_1:4d} ns')\n",
-    "        score_inst = score_function(t_1=t_1,\n",
-    "                                    total_time=total_time,\n",
-    "                                    repetitions=M,\n",
-    "                                    nblocks=N,\n",
-    "                                    indices=indices_in)\n",
+    "        print(f\"  Duration of the middle pulse: {t_1:4d} ns\")\n",
+    "        score_inst = score_function(\n",
+    "            t_1=t_1,\n",
+    "            total_time=total_time,\n",
+    "            repetitions=M,\n",
+    "            nblocks=N,\n",
+    "            indices=indices_in,\n",
+    "        )\n",
     "\n",
     "        for sc, st in zip(score_inst, scores_types):\n",
     "            scores_dict[st].append(sc)\n",
-    "            print(f'    > {st}: {sc[0]:6.3} +/- {sc[1]:6.3}')\n",
+    "            print(f\"    > {st}: {sc[0]:6.3} +/- {sc[1]:6.3}\")\n",
     "        print()\n",
     "    return scores_dict\n",
     "\n",
+    "\n",
     "def plot_scores(scores_dict):\n",
     "    fig, ax = plt.subplots(figsize=(9, 5))\n",
     "    for k in scores_dict.keys():\n",
-    "        ax.errorbar(list(durations), [s[0] for s in scores_dict[k]],\n",
-    "                    yerr=[s[1] for s in scores_dict[k]],\n",
-    "                    label=k)\n",
-    "    ax.set_title('Score vs duration $t_1$ of the central pulse', fontsize=16)\n",
-    "    ax.set_ylabel(r'Score', fontsize=16)\n",
-    "    ax.set_xlabel(r'$t_1$ (ns)', fontsize=16)\n",
+    "        ax.errorbar(\n",
+    "            list(durations),\n",
+    "            [s[0] for s in scores_dict[k]],\n",
+    "            yerr=[s[1] for s in scores_dict[k]],\n",
+    "            label=k,\n",
+    "        )\n",
+    "    ax.set_title(\"Score vs duration $t_1$ of the central pulse\", fontsize=16)\n",
+    "    ax.set_ylabel(r\"Score\", fontsize=16)\n",
+    "    ax.set_xlabel(r\"$t_1$ (ns)\", fontsize=16)\n",
     "    ax.legend()\n",
     "\n",
-    "    plt.show()\n"
+    "    plt.show()"
    ]
   },
   {
@@ -1447,13 +1461,13 @@
    ],
    "source": [
     "# Duration of the initial and final pulses\n",
-    "t_1 = 4*round(1000*np.pi/(4*2*omega_max))\n",
+    "t_1 = 4 * round(1000 * np.pi / (4 * 2 * omega_max))\n",
     "\n",
     "# Total duration of the pulse\n",
-    "total_time = 2*t_1 + 256\n",
+    "total_time = 2 * t_1 + 256\n",
     "\n",
     "# duration of the middle pulse\n",
-    "durations = range(0, total_time-2*round(t_1)-32, 32)\n",
+    "durations = range(0, total_time - 2 * round(t_1) - 32, 32)\n",
     "\n",
     "\n",
     "M = 4\n",
@@ -1465,11 +1479,7 @@
     "np.random.shuffle(indices_all)\n",
     "indices_in = indices_all[0:n_g]\n",
     "\n",
-    "scores_2layers = scan_scores(M=M,\n",
-    "                             N=N,\n",
-    "                             indices=indices_in,\n",
-    "                             durations=durations\n",
-    "                             )\n"
+    "scores_2layers = scan_scores(M=M, N=N, indices=indices_in, durations=durations)"
    ]
   },
   {
diff --git a/tutorials/applications/Using QAOA to solve a MIS problem.ipynb b/tutorials/applications/Using QAOA to solve a MIS problem.ipynb
index d2f72c265..e5cc7fe06 100644
--- a/tutorials/applications/Using QAOA to solve a MIS problem.ipynb
+++ b/tutorials/applications/Using QAOA to solve a MIS problem.ipynb
@@ -126,10 +126,14 @@
    "outputs": [],
    "source": [
     "def pos_to_graph(pos):\n",
-    "    rb = Chadoq2.rydberg_blockade_radius(1.)\n",
+    "    rb = Chadoq2.rydberg_blockade_radius(1.0)\n",
     "    g = igraph.Graph()\n",
     "    N = len(pos)\n",
-    "    edges = [[m,n] for m,n in combinations(range(N), r=2) if np.linalg.norm(pos[m] - pos[n]) < rb]           \n",
+    "    edges = [\n",
+    "        [m, n]\n",
+    "        for m, n in combinations(range(N), r=2)\n",
+    "        if np.linalg.norm(pos[m] - pos[n]) < rb\n",
+    "    ]\n",
     "    g.add_vertices(N)\n",
     "    g.add_edges(edges)\n",
     "    return g"
@@ -161,13 +165,17 @@
     }
    ],
    "source": [
-    "pos = np.array([[0., 0.], [-4, -7], [4,-7], [8,6], [-8,6]])\n",
+    "pos = np.array([[0.0, 0.0], [-4, -7], [4, -7], [8, 6], [-8, 6]])\n",
     "\n",
-    "G = pos_to_graph(pos) \n",
+    "G = pos_to_graph(pos)\n",
     "qubits = dict(enumerate(pos))\n",
     "\n",
     "reg = Register(qubits)\n",
-    "reg.draw(blockade_radius=Chadoq2.rydberg_blockade_radius(1.), draw_graph=True, draw_half_radius=True)"
+    "reg.draw(\n",
+    "    blockade_radius=Chadoq2.rydberg_blockade_radius(1.0),\n",
+    "    draw_graph=True,\n",
+    "    draw_half_radius=True,\n",
+    ")"
    ]
   },
   {
@@ -203,23 +211,23 @@
     "\n",
     "# Parametrized sequence\n",
     "seq = Sequence(reg, Chadoq2)\n",
-    "seq.declare_channel('ch0','rydberg_global')\n",
+    "seq.declare_channel(\"ch0\", \"rydberg_global\")\n",
     "\n",
-    "t_list = seq.declare_variable('t_list', size=LAYERS)\n",
-    "s_list = seq.declare_variable('s_list', size=LAYERS)\n",
+    "t_list = seq.declare_variable(\"t_list\", size=LAYERS)\n",
+    "s_list = seq.declare_variable(\"s_list\", size=LAYERS)\n",
     "\n",
     "if LAYERS == 1:\n",
     "    t_list = [t_list]\n",
     "    s_list = [s_list]\n",
-    "    \n",
-    "for t, s in zip(t_list, s_list): \n",
-    "    pulse_1 = Pulse.ConstantPulse(1000*t, 1., 0., 0) \n",
-    "    pulse_2 = Pulse.ConstantPulse(1000*s, 1., 1., 0)\n",
     "\n",
-    "    seq.add(pulse_1, 'ch0')\n",
-    "    seq.add(pulse_2, 'ch0')\n",
-    "    \n",
-    "seq.measure('ground-rydberg')"
+    "for t, s in zip(t_list, s_list):\n",
+    "    pulse_1 = Pulse.ConstantPulse(1000 * t, 1.0, 0.0, 0)\n",
+    "    pulse_2 = Pulse.ConstantPulse(1000 * s, 1.0, 1.0, 0)\n",
+    "\n",
+    "    seq.add(pulse_1, \"ch0\")\n",
+    "    seq.add(pulse_2, \"ch0\")\n",
+    "\n",
+    "seq.measure(\"ground-rydberg\")"
    ]
   },
   {
@@ -246,10 +254,10 @@
     "    params = np.array(parameters)\n",
     "    t_params, s_params = np.reshape(params.astype(int), (2, LAYERS))\n",
     "    assigned_seq = seq.build(t_list=t_params, s_list=s_params)\n",
-    "    simul = Simulation(assigned_seq, sampling_rate=.01)\n",
+    "    simul = Simulation(assigned_seq, sampling_rate=0.01)\n",
     "    results = simul.run()\n",
-    "    count_dict = results.sample_final_state() #sample from the state vector \n",
-    "    return count_dict "
+    "    count_dict = results.sample_final_state()  # sample from the state vector\n",
+    "    return count_dict"
    ]
   },
   {
@@ -258,8 +266,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "guess = {'t': np.random.uniform(8, 10, LAYERS),\n",
-    "         's': np.random.uniform(1, 3, LAYERS)}"
+    "guess = {\n",
+    "    \"t\": np.random.uniform(8, 10, LAYERS),\n",
+    "    \"s\": np.random.uniform(1, 3, LAYERS),\n",
+    "}"
    ]
   },
   {
@@ -268,7 +278,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "example_dict = quantum_loop(np.r_[guess['t'], guess['s']])"
+    "example_dict = quantum_loop(np.r_[guess[\"t\"], guess[\"s\"]])"
    ]
   },
   {
@@ -286,13 +296,13 @@
    "source": [
     "def plot_distribution(C):\n",
     "    C = dict(sorted(C.items(), key=lambda item: item[1], reverse=True))\n",
-    "    indexes = ['01011', '00111']  # MIS indexes\n",
-    "    color_dict = {key:'r' if key in indexes else 'g' for key in C}\n",
-    "    plt.figure(figsize=(12,6))\n",
+    "    indexes = [\"01011\", \"00111\"]  # MIS indexes\n",
+    "    color_dict = {key: \"r\" if key in indexes else \"g\" for key in C}\n",
+    "    plt.figure(figsize=(12, 6))\n",
     "    plt.xlabel(\"bitstrings\")\n",
     "    plt.ylabel(\"counts\")\n",
-    "    plt.bar(C.keys(), C.values(), width=0.5, color = color_dict.values())\n",
-    "    plt.xticks(rotation='vertical')\n",
+    "    plt.bar(C.keys(), C.values(), width=0.5, color=color_dict.values())\n",
+    "    plt.xticks(rotation=\"vertical\")\n",
     "    plt.show()"
    ]
   },
@@ -357,12 +367,13 @@
     "    z = np.array(list(bitstring), dtype=int)\n",
     "    A = np.array(G.get_adjacency().data)\n",
     "    # Add penalty and bias:\n",
-    "    cost = penalty*(z.T @ np.triu(A) @ z) - np.sum(z)\n",
-    "    return cost \n",
+    "    cost = penalty * (z.T @ np.triu(A) @ z) - np.sum(z)\n",
+    "    return cost\n",
+    "\n",
     "\n",
-    "def get_cost(counter,G):\n",
-    "    cost = sum(counter[key] * get_cost_colouring(key,G) for key in counter) \n",
-    "    return cost / sum(counter.values()) # Divide by total samples"
+    "def get_cost(counter, G):\n",
+    "    cost = sum(counter[key] * get_cost_colouring(key, G) for key in counter)\n",
+    "    return cost / sum(counter.values())  # Divide by total samples"
    ]
   },
   {
@@ -384,7 +395,7 @@
     }
    ],
    "source": [
-    "get_cost_colouring('00111', G)"
+    "get_cost_colouring(\"00111\", G)"
    ]
   },
   {
@@ -413,10 +424,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def func(param,*args):\n",
+    "def func(param, *args):\n",
     "    G = args[0]\n",
     "    C = quantum_loop(param)\n",
-    "    cost = get_cost(C,G)\n",
+    "    cost = get_cost(C, G)\n",
     "    return cost"
    ]
   },
@@ -440,13 +451,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "res = minimize(func, \n",
-    "               args=G,\n",
-    "               x0=np.r_[guess['t'], guess['s']],\n",
-    "               method='Nelder-Mead',\n",
-    "               tol=1e-5,\n",
-    "               options = {'maxiter': 100}\n",
-    "              )"
+    "res = minimize(\n",
+    "    func,\n",
+    "    args=G,\n",
+    "    x0=np.r_[guess[\"t\"], guess[\"s\"]],\n",
+    "    method=\"Nelder-Mead\",\n",
+    "    tol=1e-5,\n",
+    "    options={\"maxiter\": 100},\n",
+    ")"
    ]
   },
   {
diff --git a/tutorials/creating_sequences.ipynb b/tutorials/creating_sequences.ipynb
index 923b79958..206d8b048 100644
--- a/tutorials/creating_sequences.ipynb
+++ b/tutorials/creating_sequences.ipynb
@@ -66,7 +66,7 @@
     "reg1 = Register(qubits)  # Copy of 'reg' to keep the original intact\n",
     "print(\"The original array:\")\n",
     "reg1.draw()\n",
-    "reg1.rotate(45)   # Rotate by 45 degrees\n",
+    "reg1.rotate(45)  # Rotate by 45 degrees\n",
     "print(\"The rotated array:\")\n",
     "reg1.draw()"
    ]
@@ -84,7 +84,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "reg2 = Register.from_coordinates(square, prefix='q')   # All qubit IDs will start with 'q'\n",
+    "reg2 = Register.from_coordinates(\n",
+    "    square, prefix=\"q\"\n",
+    ")  # All qubit IDs will start with 'q'\n",
     "reg2.draw()"
    ]
   },
@@ -105,7 +107,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "reg3 = Register.square(4, spacing=5)   # 4x4 array with atoms 5 um apart\n",
+    "reg3 = Register.square(4, spacing=5)  # 4x4 array with atoms 5 um apart\n",
     "reg3.draw()"
    ]
   },
@@ -174,11 +176,11 @@
    },
    "outputs": [],
    "source": [
-    "seq.declare_channel('ch0', 'raman_local')\n",
+    "seq.declare_channel(\"ch0\", \"raman_local\")\n",
     "print(\"Available channels after declaring 'ch0':\")\n",
     "pprint(seq.available_channels)\n",
     "\n",
-    "seq.declare_channel('ch1', 'rydberg_local', initial_target=4)\n",
+    "seq.declare_channel(\"ch1\", \"rydberg_local\", initial_target=4)\n",
     "print(\"\\nAvailable channels after declaring 'ch1':\")\n",
     "pprint(seq.available_channels)"
    ]
@@ -221,7 +223,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "seq.target(1, 'ch0')"
+    "seq.target(1, \"ch0\")"
    ]
   },
   {
@@ -253,7 +255,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "seq.add(simple_pulse, 'ch0')"
+    "seq.add(simple_pulse, \"ch0\")"
    ]
   },
   {
@@ -269,7 +271,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "seq.delay(100, 'ch1')"
+    "seq.delay(100, \"ch1\")"
    ]
   },
   {
@@ -288,8 +290,10 @@
     "from pulser.waveforms import RampWaveform, BlackmanWaveform\n",
     "\n",
     "duration = 1000\n",
-    "amp_wf = BlackmanWaveform(duration, np.pi/2)   # Duration: 1000 ns, Area: pi/2\n",
-    "detuning_wf = RampWaveform(duration, -20, 20)  # Duration: 1000ns, linear sweep from -20 to 20 rad/µs"
+    "amp_wf = BlackmanWaveform(duration, np.pi / 2)  # Duration: 1000 ns, Area: pi/2\n",
+    "detuning_wf = RampWaveform(\n",
+    "    duration, -20, 20\n",
+    ")  # Duration: 1000ns, linear sweep from -20 to 20 rad/µs"
    ]
   },
   {
@@ -323,7 +327,7 @@
    },
    "outputs": [],
    "source": [
-    "amp_wf.integral   # dimensionless"
+    "amp_wf.integral  # dimensionless"
    ]
   },
   {
@@ -356,7 +360,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "seq.add(complex_pulse, 'ch1')"
+    "seq.add(complex_pulse, \"ch1\")"
    ]
   },
   {
@@ -406,8 +410,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "seq.target(4, 'ch0')\n",
-    "seq.add(complex_pulse, 'ch0')\n",
+    "seq.target(4, \"ch0\")\n",
+    "seq.add(complex_pulse, \"ch0\")\n",
     "\n",
     "print(\"Current Schedule:\")\n",
     "print(seq)\n",
@@ -436,9 +440,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "seq.target(0, 'ch1')\n",
-    "seq.add(simple_pulse, 'ch1', protocol='min-delay')\n",
-    "seq.add(simple_pulse, 'ch1', protocol='wait-for-all')\n",
+    "seq.target(0, \"ch1\")\n",
+    "seq.add(simple_pulse, \"ch1\", protocol=\"min-delay\")\n",
+    "seq.add(simple_pulse, \"ch1\", protocol=\"wait-for-all\")\n",
     "\n",
     "print(\"Current Schedule:\")\n",
     "print(seq)\n",
@@ -465,8 +469,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "seq.target(0, 'ch0')\n",
-    "seq.add(complex_pulse, 'ch0', protocol='no-delay')\n",
+    "seq.target(0, \"ch0\")\n",
+    "seq.add(complex_pulse, \"ch0\", protocol=\"no-delay\")\n",
     "\n",
     "print(\"Current Schedule:\")\n",
     "print(seq)\n",
@@ -500,7 +504,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "seq.measure(basis='ground-rydberg')"
+    "seq.measure(basis=\"ground-rydberg\")"
    ]
   },
   {
diff --git a/tutorials/quantum_simulation/Bayesian Optimisation for antiferromagnetic state preparation.ipynb b/tutorials/quantum_simulation/Bayesian Optimisation for antiferromagnetic state preparation.ipynb
index 85a792908..551633ded 100644
--- a/tutorials/quantum_simulation/Bayesian Optimisation for antiferromagnetic state preparation.ipynb
+++ b/tutorials/quantum_simulation/Bayesian Optimisation for antiferromagnetic state preparation.ipynb
@@ -87,26 +87,29 @@
    "source": [
     "# Parameters in rad/µs and ns\n",
     "\n",
-    "T = 1000    # duration\n",
+    "T = 1000  # duration\n",
     "U = 2 * np.pi * 5.0\n",
     "\n",
-    "Omega_max = 0.5 * U \n",
+    "Omega_max = 0.5 * U\n",
     "\n",
-    "delta_0 = -1.0 * U \n",
-    "delta_f = 1.0 * U \n",
+    "delta_0 = -1.0 * U\n",
+    "delta_f = 1.0 * U\n",
     "\n",
     "R_interatomic = Chadoq2.rydberg_blockade_radius(Omega_max) / 1.2\n",
-    "print(f'Interatomic Radius is: {R_interatomic}µm.')\n",
+    "print(f\"Interatomic Radius is: {R_interatomic}µm.\")\n",
     "\n",
     "N_side = 4\n",
     "coords = (\n",
-    "    [R_interatomic * np.r_[x,0] for x in range(N_side-1)] \n",
-    "    + [R_interatomic * np.r_[N_side-1,y] for y in range(N_side-1)]\n",
-    "    + [R_interatomic * np.r_[N_side-1-x,N_side-1] for x in range(N_side-1)] \n",
-    "    + [R_interatomic * np.r_[0,N_side-1-y] for y in range(N_side-1)]\n",
+    "    [R_interatomic * np.r_[x, 0] for x in range(N_side - 1)]\n",
+    "    + [R_interatomic * np.r_[N_side - 1, y] for y in range(N_side - 1)]\n",
+    "    + [\n",
+    "        R_interatomic * np.r_[N_side - 1 - x, N_side - 1]\n",
+    "        for x in range(N_side - 1)\n",
+    "    ]\n",
+    "    + [R_interatomic * np.r_[0, N_side - 1 - y] for y in range(N_side - 1)]\n",
     ")\n",
-    "reg = Register.from_coordinates(coords, prefix='q')\n",
-    "N=len(coords)\n",
+    "reg = Register.from_coordinates(coords, prefix=\"q\")\n",
+    "N = len(coords)\n",
     "N_samples = 1000\n",
     "reg.draw()"
    ]
@@ -133,15 +136,15 @@
    ],
    "source": [
     "seq = Sequence(reg, Chadoq2)\n",
-    "seq.declare_channel('ising', 'rydberg_global')\n",
+    "seq.declare_channel(\"ising\", \"rydberg_global\")\n",
     "\n",
     "tol = 1e-6\n",
-    "max_amp = seq.declared_channels['ising'].max_amp * (1-tol)\n",
-    "max_det = seq.declared_channels['ising'].max_abs_detuning * (1-tol)\n",
-    "Omega_max=min(max_amp, Omega_max)\n",
-    "delta_0=np.sign(delta_0)*min(max_det, abs(delta_0))\n",
-    "delta_f=np.sign(delta_f)*min(max_det, abs(delta_f))\n",
-    "print(Omega_max/U, np.round(delta_0/U,2), delta_f/U)"
+    "max_amp = seq.declared_channels[\"ising\"].max_amp * (1 - tol)\n",
+    "max_det = seq.declared_channels[\"ising\"].max_abs_detuning * (1 - tol)\n",
+    "Omega_max = min(max_amp, Omega_max)\n",
+    "delta_0 = np.sign(delta_0) * min(max_det, abs(delta_0))\n",
+    "delta_f = np.sign(delta_f) * min(max_det, abs(delta_f))\n",
+    "print(Omega_max / U, np.round(delta_0 / U, 2), delta_f / U)"
    ]
   },
   {
@@ -178,8 +181,8 @@
     "m = 3\n",
     "\n",
     "# Random instance of the parameter space\n",
-    "amp_params = np.random.uniform(0,Omega_max,m)\n",
-    "det_params = np.random.uniform(delta_0,delta_f,m)"
+    "amp_params = np.random.uniform(0, Omega_max, m)\n",
+    "det_params = np.random.uniform(delta_0, delta_f, m)"
    ]
   },
   {
@@ -197,7 +200,7 @@
    "source": [
     "def create_interp_pulse(amp_params, det_params):\n",
     "    return Pulse(\n",
-    "        InterpolatedWaveform(T, [1e-9, *amp_params, 1e-9]), \n",
+    "        InterpolatedWaveform(T, [1e-9, *amp_params, 1e-9]),\n",
     "        InterpolatedWaveform(T, [delta_0, *det_params, delta_f]),\n",
     "        0,\n",
     "    )"
@@ -228,8 +231,8 @@
    ],
    "source": [
     "seq = Sequence(reg, Chadoq2)\n",
-    "seq.declare_channel('ising', 'rydberg_global')\n",
-    "seq.add(create_interp_pulse(amp_params, det_params),'ising')\n",
+    "seq.declare_channel(\"ising\", \"rydberg_global\")\n",
+    "seq.add(create_interp_pulse(amp_params, det_params), \"ising\")\n",
     "seq.draw()"
    ]
   },
@@ -279,37 +282,43 @@
    "outputs": [],
    "source": [
     "def occupation(j, N):\n",
-    "    up = qutip.basis(2,0)\n",
+    "    up = qutip.basis(2, 0)\n",
     "    prod = [qutip.qeye(2) for _ in range(N)]\n",
     "    prod[j] = up * up.dag()\n",
     "    return qutip.tensor(prod)\n",
     "\n",
+    "\n",
     "def get_corr_pairs(k, N):\n",
-    "    corr_pairs = [[i,(i+k)%N] for i in range(N)]\n",
+    "    corr_pairs = [[i, (i + k) % N] for i in range(N)]\n",
     "    return corr_pairs\n",
     "\n",
+    "\n",
     "def get_corr_function(k, N, state):\n",
     "    corr_pairs = get_corr_pairs(k, N)\n",
     "    operators = [occupation(j, N) for j in range(N)]\n",
     "    covariance = 0\n",
     "    for qi, qj in corr_pairs:\n",
-    "        covariance += qutip.expect(operators[qi]*operators[qj], state)\n",
-    "        covariance -= qutip.expect(operators[qi], state)*qutip.expect(operators[qj], state)\n",
-    "    return covariance/len(corr_pairs)   \n",
+    "        covariance += qutip.expect(operators[qi] * operators[qj], state)\n",
+    "        covariance -= qutip.expect(operators[qi], state) * qutip.expect(\n",
+    "            operators[qj], state\n",
+    "        )\n",
+    "    return covariance / len(corr_pairs)\n",
+    "\n",
     "\n",
     "def get_full_corr_function(reg, state):\n",
     "    N = len(reg.qubits)\n",
     "    correlation_function = {}\n",
-    "    for k in range(-N//2, N//2+1):\n",
+    "    for k in range(-N // 2, N // 2 + 1):\n",
     "        correlation_function[k] = get_corr_function(k, N, state)\n",
     "    return correlation_function\n",
     "\n",
+    "\n",
     "def get_neel_structure_factor(reg, state):\n",
     "    N = len(reg.qubits)\n",
     "    st_fac = 0\n",
-    "    for k in range(-N//2, N//2+1):\n",
+    "    for k in range(-N // 2, N // 2 + 1):\n",
     "        kk = np.abs(k)\n",
-    "        st_fac += 4 * (-1)**kk * get_corr_function(k, N, state)\n",
+    "        st_fac += 4 * (-1) ** kk * get_corr_function(k, N, state)\n",
     "    return st_fac"
    ]
   },
@@ -328,7 +337,11 @@
    "source": [
     "def proba_from_state(results, min_p=0.1):\n",
     "    sampling = results.sample_final_state(N_samples=N_samples)\n",
-    "    return {k: f'{100*v/N_samples}%' for k, v in sampling.items() if v/N_samples > min_p}"
+    "    return {\n",
+    "        k: f\"{100*v/N_samples}%\"\n",
+    "        for k, v in sampling.items()\n",
+    "        if v / N_samples > min_p\n",
+    "    }"
    ]
   },
   {
@@ -358,11 +371,11 @@
     "AF2 = qutip.tensor([qutip.basis(2, (k + 1) % 2) for k in range(N)])\n",
     "AF_state = (AF1 + AF2).unit()\n",
     "\n",
-    "t1=time.process_time()\n",
+    "t1 = time.process_time()\n",
     "S_max = get_neel_structure_factor(reg, AF_state)\n",
-    "print('S_Neel(AF state) =', S_max)\n",
-    "t2=time.process_time()\n",
-    "print('computed in', (t2-t1),'sec')"
+    "print(\"S_Neel(AF state) =\", S_max)\n",
+    "t2 = time.process_time()\n",
+    "print(\"computed in\", (t2 - t1), \"sec\")"
    ]
   },
   {
@@ -387,16 +400,21 @@
    "source": [
     "def score(params):\n",
     "    seq = Sequence(reg, Chadoq2)\n",
-    "    seq.declare_channel('ising', 'rydberg_global')\n",
-    "    seq.add(create_interp_pulse(params[:m], params[m:]),'ising')\n",
-    "    \n",
+    "    seq.declare_channel(\"ising\", \"rydberg_global\")\n",
+    "    seq.add(create_interp_pulse(params[:m], params[m:]), \"ising\")\n",
+    "\n",
     "    simul = Simulation(seq, sampling_rate=0.5)\n",
     "    results = simul.run()\n",
     "\n",
     "    sampling = results.sample_final_state(N_samples=N_samples)\n",
-    "    sampled_state = sum([np.sqrt(sampling[k]/N_samples)*qutip.ket(k) for k in sampling.keys()])\n",
+    "    sampled_state = sum(\n",
+    "        [\n",
+    "            np.sqrt(sampling[k] / N_samples) * qutip.ket(k)\n",
+    "            for k in sampling.keys()\n",
+    "        ]\n",
+    "    )\n",
     "\n",
-    "    F = get_neel_structure_factor(reg, sampled_state)/S_max\n",
+    "    F = get_neel_structure_factor(reg, sampled_state) / S_max\n",
     "\n",
     "    return 1 - F"
    ]
@@ -472,7 +490,9 @@
     "n_r = 30\n",
     "n_c = 120\n",
     "\n",
-    "RESULT = gp_minimize(score, bounds, n_random_starts=n_r, n_calls=n_c, verbose=False)"
+    "RESULT = gp_minimize(\n",
+    "    score, bounds, n_random_starts=n_r, n_calls=n_c, verbose=False\n",
+    ")"
    ]
   },
   {
@@ -498,10 +518,10 @@
     "def sort_improv(RESULT):\n",
     "    score_vals = RESULT.func_vals\n",
     "    min = score_vals[0]\n",
-    "    score_list=[]\n",
+    "    score_list = []\n",
     "    for s in score_vals:\n",
-    "        if s<min:\n",
-    "            min=s\n",
+    "        if s < min:\n",
+    "            min = s\n",
     "            score_list.append(min)\n",
     "        else:\n",
     "            score_list.append(min)\n",
@@ -528,11 +548,11 @@
    ],
    "source": [
     "fig = plt.figure()\n",
-    "plt.semilogy(range(n_c), sort_improv(RESULT), 'r-')\n",
-    "plt.xlabel(r'$n_{calls}$', fontsize=14)\n",
-    "plt.ylabel('Score reached', fontsize=14)\n",
-    "plt.legend(['1 run'],fontsize=12)\n",
-    "plt.xlim(0,n_c)\n",
+    "plt.semilogy(range(n_c), sort_improv(RESULT), \"r-\")\n",
+    "plt.xlabel(r\"$n_{calls}$\", fontsize=14)\n",
+    "plt.ylabel(\"Score reached\", fontsize=14)\n",
+    "plt.legend([\"1 run\"], fontsize=12)\n",
+    "plt.xlim(0, n_c)\n",
     "plt.show()"
    ]
   },
@@ -579,9 +599,9 @@
    ],
    "source": [
     "seq = Sequence(reg, Chadoq2)\n",
-    "seq.declare_channel('ising', 'rydberg_global')\n",
+    "seq.declare_channel(\"ising\", \"rydberg_global\")\n",
     "P = create_interp_pulse(RESULT.x[:m], RESULT.x[m:])\n",
-    "seq.add(P,'ising')\n",
+    "seq.add(P, \"ising\")\n",
     "seq.draw()"
    ]
   },
@@ -613,21 +633,21 @@
    "source": [
     "omega = P.amplitude.samples / U\n",
     "delta = P.detuning.samples / U\n",
-    "        \n",
+    "\n",
     "fig, ax = plt.subplots()\n",
-    "ax.grid(True, which='both')\n",
+    "ax.grid(True, which=\"both\")\n",
     "\n",
     "ax.set_ylabel(r\"$\\hbar\\delta(t)/U$\", fontsize=14)\n",
     "ax.set_xlabel(r\"$\\hbar\\Omega(t)/U$\", fontsize=14)\n",
-    "ax.set_xlim(0, Omega_max/U*1.1)\n",
-    "ax.axhline(y=0, color='k')\n",
-    "ax.axvline(x=0, color='k')\n",
+    "ax.set_xlim(0, Omega_max / U * 1.1)\n",
+    "ax.axhline(y=0, color=\"k\")\n",
+    "ax.axvline(x=0, color=\"k\")\n",
     "\n",
     "y = np.arange(0.0, 3, 0.01)\n",
-    "x = 0.5 * (1 - 1.0 * (y - 1)**2)\n",
+    "x = 0.5 * (1 - 1.0 * (y - 1) ** 2)\n",
     "ax.fill_between(x, y, alpha=0.4)\n",
     "\n",
-    "ax.plot(omega,delta, 'red', lw=2)\n",
+    "ax.plot(omega, delta, \"red\", lw=2)\n",
     "plt.show()"
    ]
   },
@@ -648,13 +668,15 @@
    "source": [
     "simul = Simulation(seq)\n",
     "results = simul.run()\n",
-    "print('final =', proba_from_state(results, min_p=0.05))\n",
+    "print(\"final =\", proba_from_state(results, min_p=0.05))\n",
     "\n",
     "sampling = results.sample_final_state(N_samples=1000)\n",
-    "sampled_state = sum([np.sqrt(sampling[k]/1000) *qutip.ket(k)for k in sampling.keys()])\n",
+    "sampled_state = sum(\n",
+    "    [np.sqrt(sampling[k] / 1000) * qutip.ket(k) for k in sampling.keys()]\n",
+    ")\n",
     "\n",
-    "s_neel_sampled = np.round(get_neel_structure_factor(reg, sampled_state),3)\n",
-    "print(f'S_Neel (final_sampled) = {s_neel_sampled}')"
+    "s_neel_sampled = np.round(get_neel_structure_factor(reg, sampled_state), 3)\n",
+    "print(f\"S_Neel (final_sampled) = {s_neel_sampled}\")"
    ]
   },
   {
@@ -673,7 +695,10 @@
     "# Skip parameter. Only 1 out of 10 pts are kept to enhance the speed.\n",
     "a = 10\n",
     "time_domain_reduced = np.arange(T)[::a]\n",
-    "S_Neel_states=[get_neel_structure_factor(reg, results.states[t]) for t in time_domain_reduced]"
+    "S_Neel_states = [\n",
+    "    get_neel_structure_factor(reg, results.states[t])\n",
+    "    for t in time_domain_reduced\n",
+    "]"
    ]
   },
   {
@@ -696,9 +721,9 @@
    ],
    "source": [
     "plt.figure()\n",
-    "plt.plot(time_domain_reduced,S_Neel_states)\n",
-    "plt.xlabel('Time (ns)',fontsize=12)\n",
-    "plt.ylabel(r'$S_{Néel}$',fontsize=14)\n",
+    "plt.plot(time_domain_reduced, S_Neel_states)\n",
+    "plt.xlabel(\"Time (ns)\", fontsize=12)\n",
+    "plt.ylabel(r\"$S_{Néel}$\", fontsize=14)\n",
     "\n",
     "plt.show()"
    ]
@@ -733,14 +758,16 @@
     "final = results.get_final_state()\n",
     "correlation_function = get_full_corr_function(reg, final)\n",
     "\n",
-    "A = np.reshape(list(correlation_function.values()), (1, N+1))\n",
-    "A[0, N//2] = None\n",
+    "A = np.reshape(list(correlation_function.values()), (1, N + 1))\n",
+    "A[0, N // 2] = None\n",
     "\n",
-    "plt.figure(figsize=(10,5))\n",
-    "plt.imshow(A, cmap='coolwarm')\n",
-    "plt.xlabel('k',fontsize=12)\n",
-    "plt.xticks(range(A.shape[1]), ['{}'.format(i) for i in range(-N//2 , N//2+1)])\n",
-    "plt.title(r'$g^{(2)}(k)$ after simulation', fontsize=14)\n",
+    "plt.figure(figsize=(10, 5))\n",
+    "plt.imshow(A, cmap=\"coolwarm\")\n",
+    "plt.xlabel(\"k\", fontsize=12)\n",
+    "plt.xticks(\n",
+    "    range(A.shape[1]), [\"{}\".format(i) for i in range(-N // 2, N // 2 + 1)]\n",
+    ")\n",
+    "plt.title(r\"$g^{(2)}(k)$ after simulation\", fontsize=14)\n",
     "plt.yticks([])\n",
     "plt.colorbar()\n",
     "plt.show()"
diff --git a/tutorials/quantum_simulation/Building 1D Rydberg Crystals.ipynb b/tutorials/quantum_simulation/Building 1D Rydberg Crystals.ipynb
index c3cecfa19..87c5c4bc8 100644
--- a/tutorials/quantum_simulation/Building 1D Rydberg Crystals.ipynb
+++ b/tutorials/quantum_simulation/Building 1D Rydberg Crystals.ipynb
@@ -51,10 +51,10 @@
    "outputs": [],
    "source": [
     "def occupation(reg, j):\n",
-    "    r = qutip.basis(2,0)\n",
+    "    r = qutip.basis(2, 0)\n",
     "    N = len(reg.qubits)\n",
     "    prod = [qutip.qeye(2) for _ in range(N)]\n",
-    "    prod[j] = r*r.dag()\n",
+    "    prod[j] = r * r.dag()\n",
     "    return qutip.tensor(prod)"
    ]
   },
@@ -79,30 +79,39 @@
     }
    ],
    "source": [
-    "Omega_max = 2 * 2*np.pi\n",
+    "Omega_max = 2 * 2 * np.pi\n",
     "R_blockade = MockDevice.rydberg_blockade_radius(Omega_max)\n",
-    "print(f'Blockade Radius is: {R_blockade}µm.')\n",
+    "print(f\"Blockade Radius is: {R_blockade}µm.\")\n",
     "groups = 3\n",
     "\n",
+    "\n",
     "def blockade_cluster(N):\n",
-    "    \n",
-    "    atom_coords = [((R_blockade/N)*x+40*group,0) for group in range(groups) for x in range(1,N+1)]\n",
-    "    reg = Register.from_coordinates(atom_coords, prefix='q')\n",
-    "    reg.draw(blockade_radius=R_blockade, draw_half_radius=True, draw_graph=False)\n",
-    "    \n",
-    "    resonant_pulse = Pulse.ConstantPulse(1500, Omega_max, 0., 0.)\n",
-    "    \n",
+    "\n",
+    "    atom_coords = [\n",
+    "        ((R_blockade / N) * x + 40 * group, 0)\n",
+    "        for group in range(groups)\n",
+    "        for x in range(1, N + 1)\n",
+    "    ]\n",
+    "    reg = Register.from_coordinates(atom_coords, prefix=\"q\")\n",
+    "    reg.draw(\n",
+    "        blockade_radius=R_blockade, draw_half_radius=True, draw_graph=False\n",
+    "    )\n",
+    "\n",
+    "    resonant_pulse = Pulse.ConstantPulse(1500, Omega_max, 0.0, 0.0)\n",
+    "\n",
     "    seq = Sequence(reg, MockDevice)\n",
-    "    seq.declare_channel('ising', 'rydberg_global')\n",
-    "    seq.add(resonant_pulse, 'ising')\n",
-    "    \n",
-    "    simul = Simulation(seq, sampling_rate=.7)\n",
-    "    \n",
-    "    obs = [sum(occupation(reg, j) for j in range(i,i+N)) for i in range(0,groups*N,N)]    \n",
-    "    \n",
-    "    res = simul.run(progress_bar=True, method='bdf')\n",
-    "    return res.expect(obs)\n",
-    "    "
+    "    seq.declare_channel(\"ising\", \"rydberg_global\")\n",
+    "    seq.add(resonant_pulse, \"ising\")\n",
+    "\n",
+    "    simul = Simulation(seq, sampling_rate=0.7)\n",
+    "\n",
+    "    obs = [\n",
+    "        sum(occupation(reg, j) for j in range(i, i + N))\n",
+    "        for i in range(0, groups * N, N)\n",
+    "    ]\n",
+    "\n",
+    "    res = simul.run(progress_bar=True, method=\"bdf\")\n",
+    "    return res.expect(obs)"
    ]
   },
   {
@@ -203,7 +212,7 @@
     }
    ],
    "source": [
-    "data = [blockade_cluster(N) for N in [1,2,3]]  "
+    "data = [blockade_cluster(N) for N in [1, 2, 3]]"
    ]
   },
   {
@@ -257,11 +266,11 @@
    ],
    "source": [
     "for N, expectation in enumerate(data):\n",
-    "    plt.xlabel('Time (µs)', fontsize=14)\n",
-    "    plt.ylabel('Probability of Rydberg state', fontsize=14)\n",
-    "    plt.title(f'Atoms per cluster N = {N+1}', fontsize = 15)\n",
-    "    avg = sum(expectation)/groups\n",
-    "    plt.plot(np.arange(len(avg))/1000, avg)\n",
+    "    plt.xlabel(\"Time (µs)\", fontsize=14)\n",
+    "    plt.ylabel(\"Probability of Rydberg state\", fontsize=14)\n",
+    "    plt.title(f\"Atoms per cluster N = {N+1}\", fontsize=15)\n",
+    "    avg = sum(expectation) / groups\n",
+    "    plt.plot(np.arange(len(avg)) / 1000, avg)\n",
     "    plt.show()"
    ]
   },
@@ -293,9 +302,9 @@
    "outputs": [],
    "source": [
     "# Parameters in rad/µs and ns\n",
-    "delta_0 = -6 * 2*np.pi\n",
-    "delta_f = 10 * 2*np.pi\n",
-    "Omega_max = 2 * 2*np.pi  \n",
+    "delta_0 = -6 * 2 * np.pi\n",
+    "delta_f = 10 * 2 * np.pi\n",
+    "Omega_max = 2 * 2 * np.pi\n",
     "t_rise = 500\n",
     "t_stop = 4500\n",
     "t_end = 5000"
@@ -335,10 +344,10 @@
    ],
    "source": [
     "R_blockade = MockDevice.rydberg_blockade_radius(Omega_max)\n",
-    "a = 7.\n",
+    "a = 7.0\n",
     "\n",
-    "reg = Register.rectangle(1,11, spacing=a, prefix='q')\n",
-    "print(f'Blockade Radius is: {R_blockade}µm.')\n",
+    "reg = Register.rectangle(1, 11, spacing=a, prefix=\"q\")\n",
+    "print(f\"Blockade Radius is: {R_blockade}µm.\")\n",
     "reg.draw(blockade_radius=R_blockade, draw_half_radius=True)"
    ]
   },
@@ -357,8 +366,10 @@
    "source": [
     "hold = ConstantWaveform(t_rise, delta_0)\n",
     "excite = RampWaveform(t_stop - t_rise, delta_0, delta_f)\n",
-    "sweep = Pulse.ConstantAmplitude(Omega_max, CompositeWaveform(hold, excite), 0.)\n",
-    "stay = Pulse.ConstantPulse(t_end - t_stop, 0., delta_f, 0.)"
+    "sweep = Pulse.ConstantAmplitude(\n",
+    "    Omega_max, CompositeWaveform(hold, excite), 0.0\n",
+    ")\n",
+    "stay = Pulse.ConstantPulse(t_end - t_stop, 0.0, delta_f, 0.0)"
    ]
   },
   {
@@ -379,12 +390,12 @@
    ],
    "source": [
     "seq = Sequence(reg, MockDevice)\n",
-    "seq.declare_channel('ising', 'rydberg_global')\n",
+    "seq.declare_channel(\"ising\", \"rydberg_global\")\n",
     "\n",
-    "seq.add(sweep, 'ising')\n",
-    "seq.add(stay, 'ising')\n",
+    "seq.add(sweep, \"ising\")\n",
+    "seq.add(stay, \"ising\")\n",
     "\n",
-    "#print(seq)\n",
+    "# print(seq)\n",
     "seq.draw()"
    ]
   },
@@ -403,34 +414,33 @@
    "source": [
     "def phase_diagram(seq):\n",
     "    ratio = []\n",
-    "    for x in seq._schedule['ising']:\n",
-    "        if isinstance(x.type,Pulse):\n",
-    "            ratio += list(x.type.detuning.samples/Omega_max)\n",
+    "    for x in seq._schedule[\"ising\"]:\n",
+    "        if isinstance(x.type, Pulse):\n",
+    "            ratio += list(x.type.detuning.samples / Omega_max)\n",
     "\n",
     "    fig, ax = plt.subplots()\n",
-    "    ax.grid(True, which='both')\n",
+    "    ax.grid(True, which=\"both\")\n",
     "\n",
     "    ax.set_ylabel(r\"Interaction Range $R_b/a$\", fontsize=14)\n",
     "    ax.set_xlabel(r\"Detuning $\\delta/\\Omega$\", fontsize=14)\n",
-    "    ax.set_xlim(-4,6)\n",
-    "    ax.set_ylim(0,3.2)\n",
-    "    ax.axhline(y=0, color='k')\n",
-    "    ax.axvline(x=0, color='k')\n",
+    "    ax.set_xlim(-4, 6)\n",
+    "    ax.set_ylim(0, 3.2)\n",
+    "    ax.axhline(y=0, color=\"k\")\n",
+    "    ax.axvline(x=0, color=\"k\")\n",
     "\n",
     "    y = np.arange(0.0, 5, 0.01)\n",
-    "    x = 2*(0.6+8*(y-1.2)**2)\n",
+    "    x = 2 * (0.6 + 8 * (y - 1.2) ** 2)\n",
     "    ax.fill_between(x, y, alpha=0.4)\n",
     "\n",
     "    y = np.arange(0.0, 5, 0.01)\n",
-    "    x = 2*(0.8+50*(y-2.45)**2)\n",
+    "    x = 2 * (0.8 + 50 * (y - 2.45) ** 2)\n",
     "    ax.fill_between(x, y, alpha=0.4)\n",
     "\n",
     "    y = np.arange(0.0, 5, 0.01)\n",
-    "    x = 2*(1.+170*(y-3.06)**2)\n",
+    "    x = 2 * (1.0 + 170 * (y - 3.06) ** 2)\n",
     "    ax.fill_between(x, y, alpha=0.4)\n",
     "\n",
-    "\n",
-    "    ax.plot(np.array(ratio), np.full(len(ratio), R_blockade/a),'red',lw=2)\n",
+    "    ax.plot(np.array(ratio), np.full(len(ratio), R_blockade / a), \"red\", lw=2)\n",
     "    plt.show()"
    ]
   },
@@ -493,7 +503,7 @@
     }
    ],
    "source": [
-    "simul = Simulation(seq, sampling_rate=.1)\n",
+    "simul = Simulation(seq, sampling_rate=0.1)\n",
     "\n",
     "occup_list = [occupation(reg, j) for j in range(len(reg.qubits))]\n",
     "\n",
@@ -515,11 +525,11 @@
    "outputs": [],
    "source": [
     "def plot_evolution(results):\n",
-    "    plt.figure(figsize=(10,5))\n",
-    "    plt.xlabel('Time (µs)', fontsize=14)\n",
-    "    plt.ylabel('Rydberg Occupation Probabilty', fontsize=14)\n",
+    "    plt.figure(figsize=(10, 5))\n",
+    "    plt.xlabel(\"Time (µs)\", fontsize=14)\n",
+    "    plt.ylabel(\"Rydberg Occupation Probabilty\", fontsize=14)\n",
     "    for expv in results:\n",
-    "        plt.plot(np.arange(len(expv))/1000,expv)\n",
+    "        plt.plot(np.arange(len(expv)) / 1000, expv)\n",
     "\n",
     "    plt.show()"
    ]
@@ -559,29 +569,31 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def heat_detuning(data,start,end):\n",
+    "def heat_detuning(data, start, end):\n",
     "    N = len(reg.qubits)\n",
     "    time_window = []\n",
-    "    x =[]\n",
-    "    detunings = simul.samples['Global']['ground-rydberg']['det'][[int(1000*t) for t in simul.sampling_times]]\n",
+    "    x = []\n",
+    "    detunings = simul.samples[\"Global\"][\"ground-rydberg\"][\"det\"][\n",
+    "        [int(1000 * t) for t in simul.sampling_times]\n",
+    "    ]\n",
     "\n",
-    "    for t,d in enumerate(detunings):\n",
+    "    for t, d in enumerate(detunings):\n",
     "        if start <= d <= end:\n",
     "            time_window.append(t)\n",
-    "            x.append(d/(2*np.pi))\n",
+    "            x.append(d / (2 * np.pi))\n",
     "\n",
-    "    y = np.arange(1,N+1)\n",
+    "    y = np.arange(1, N + 1)\n",
     "\n",
     "    X, Y = np.meshgrid(x, y)\n",
-    "    Z = np.array(data)[:,time_window]\n",
+    "    Z = np.array(data)[:, time_window]\n",
     "\n",
-    "    plt.figure(figsize=(14,3))\n",
-    "    plt.pcolormesh(X,Y,Z, cmap='hot', shading='auto')\n",
-    "    plt.xlabel('Detuning/2π (MHz)', fontsize=14)\n",
-    "    plt.ylabel('Atom in array', fontsize=14)\n",
-    "    plt.yticks(range(1,N+1), [f'q{i}' for i in range(N)], va='center')\n",
+    "    plt.figure(figsize=(14, 3))\n",
+    "    plt.pcolormesh(X, Y, Z, cmap=\"hot\", shading=\"auto\")\n",
+    "    plt.xlabel(\"Detuning/2π (MHz)\", fontsize=14)\n",
+    "    plt.ylabel(\"Atom in array\", fontsize=14)\n",
+    "    plt.yticks(range(1, N + 1), [f\"q{i}\" for i in range(N)], va=\"center\")\n",
     "    plt.colorbar(fraction=0.047, pad=0.015)\n",
-    "    \n",
+    "\n",
     "    plt.show()"
    ]
   },
@@ -693,11 +705,11 @@
    ],
    "source": [
     "a = 3.5\n",
-    "reg = Register.rectangle(1, 10, spacing=a, prefix='q')\n",
+    "reg = Register.rectangle(1, 10, spacing=a, prefix=\"q\")\n",
     "\n",
-    "delta_0 = -4 * 2*np.pi\n",
-    "delta_f = 10 * 2*np.pi\n",
-    "Omega_max = 2.0 * 2*np.pi  # btw 1.8-2 * 2pi MHz\n",
+    "delta_0 = -4 * 2 * np.pi\n",
+    "delta_f = 10 * 2 * np.pi\n",
+    "Omega_max = 2.0 * 2 * np.pi  # btw 1.8-2 * 2pi MHz\n",
     "t_rise = 600\n",
     "t_stop = 2500\n",
     "t_end = 3000\n",
@@ -707,25 +719,27 @@
     "#\n",
     "hold = ConstantWaveform(t_rise, delta_0)\n",
     "excite = RampWaveform(t_stop - t_rise, delta_0, delta_f)\n",
-    "sweep = Pulse.ConstantAmplitude(Omega_max, CompositeWaveform(hold, excite), 0.)\n",
-    "stay = Pulse.ConstantPulse(t_end - t_stop, 0., delta_f, 0.)\n",
+    "sweep = Pulse.ConstantAmplitude(\n",
+    "    Omega_max, CompositeWaveform(hold, excite), 0.0\n",
+    ")\n",
+    "stay = Pulse.ConstantPulse(t_end - t_stop, 0.0, delta_f, 0.0)\n",
     "\n",
     "#\n",
     "seq = Sequence(reg, MockDevice)\n",
-    "seq.declare_channel('ising', 'rydberg_global')\n",
+    "seq.declare_channel(\"ising\", \"rydberg_global\")\n",
     "\n",
-    "seq.add(sweep, 'ising')\n",
-    "seq.add(stay, 'ising')\n",
+    "seq.add(sweep, \"ising\")\n",
+    "seq.add(stay, \"ising\")\n",
     "\n",
-    "#print(seq)\n",
-    "#seq.draw()\n",
+    "# print(seq)\n",
+    "# seq.draw()\n",
     "phase_diagram(seq)\n",
     "\n",
     "simul = Simulation(seq, sampling_rate=0.1)\n",
     "\n",
     "occup_list = [occupation(reg, j) for j in range(simul._size)]\n",
     "\n",
-    "res = simul.run(progress_bar=True, method='bdf')\n",
+    "res = simul.run(progress_bar=True, method=\"bdf\")\n",
     "\n",
     "occupations = res.expect(occup_list)\n",
     "\n",
@@ -819,13 +833,13 @@
    ],
    "source": [
     "a = 2.8\n",
-    "reg = Register.rectangle(1,9, spacing=a, prefix='q')\n",
+    "reg = Register.rectangle(1, 9, spacing=a, prefix=\"q\")\n",
     "\n",
     "# Parameters in rad/µs and ns\n",
     "\n",
-    "delta_0 = -4 * 2*np.pi\n",
-    "delta_f = 10 * 2*np.pi\n",
-    "Omega_max = 2.0 * 2*np.pi  # btw 1.8-2 2pi*MHz\n",
+    "delta_0 = -4 * 2 * np.pi\n",
+    "delta_f = 10 * 2 * np.pi\n",
+    "Omega_max = 2.0 * 2 * np.pi  # btw 1.8-2 2pi*MHz\n",
     "t_rise = 600\n",
     "t_stop = 2500\n",
     "t_end = 3000\n",
@@ -835,28 +849,30 @@
     "#\n",
     "hold = ConstantWaveform(t_rise, delta_0)\n",
     "excite = RampWaveform(t_stop - t_rise, delta_0, delta_f)\n",
-    "sweep = Pulse.ConstantAmplitude(Omega_max, CompositeWaveform(hold, excite), 0.)\n",
-    "stay = Pulse.ConstantPulse(t_end - t_stop, 0., delta_f, 0.)\n",
+    "sweep = Pulse.ConstantAmplitude(\n",
+    "    Omega_max, CompositeWaveform(hold, excite), 0.0\n",
+    ")\n",
+    "stay = Pulse.ConstantPulse(t_end - t_stop, 0.0, delta_f, 0.0)\n",
     "\n",
     "#\n",
     "seq = Sequence(reg, MockDevice)\n",
-    "seq.declare_channel('ising', 'rydberg_global')\n",
+    "seq.declare_channel(\"ising\", \"rydberg_global\")\n",
     "\n",
-    "seq.add(sweep, 'ising')\n",
-    "seq.add(stay, 'ising')\n",
+    "seq.add(sweep, \"ising\")\n",
+    "seq.add(stay, \"ising\")\n",
     "\n",
     "phase_diagram(seq)\n",
     "\n",
     "simul = Simulation(seq, sampling_rate=0.6)\n",
-    "    \n",
+    "\n",
     "occup_list = [occupation(reg, j) for j in range(simul._size)]\n",
     "\n",
     "#\n",
-    "res = simul.run(progress_bar=True, method='bdf')\n",
+    "res = simul.run(progress_bar=True, method=\"bdf\")\n",
     "occupations = res.expect(occup_list)\n",
     "\n",
-    "plot_evolution(occupations)        \n",
-    "heat_detuning(occupations,delta_0,delta_f)\n"
+    "plot_evolution(occupations)\n",
+    "heat_detuning(occupations, delta_0, delta_f)"
    ]
   }
  ],
diff --git a/tutorials/quantum_simulation/Microwave-engineering of programmable XXZ Hamiltonians in arrays of Rydberg atoms.ipynb b/tutorials/quantum_simulation/Microwave-engineering of programmable XXZ Hamiltonians in arrays of Rydberg atoms.ipynb
index 0fec7863c..45e0ae741 100644
--- a/tutorials/quantum_simulation/Microwave-engineering of programmable XXZ Hamiltonians in arrays of Rydberg atoms.ipynb
+++ b/tutorials/quantum_simulation/Microwave-engineering of programmable XXZ Hamiltonians in arrays of Rydberg atoms.ipynb
@@ -52,10 +52,16 @@
     "# Times are in ns\n",
     "t_pulse = 26\n",
     "\n",
-    "X_pulse = Pulse.ConstantDetuning(BlackmanWaveform(t_pulse, np.pi/2.), 0, 0)\n",
-    "Y_pulse = Pulse.ConstantDetuning(BlackmanWaveform(t_pulse, np.pi/2.), 0, -np.pi/2)\n",
-    "mX_pulse = Pulse.ConstantDetuning(BlackmanWaveform(t_pulse, np.pi/2.), 0, np.pi)\n",
-    "mY_pulse = Pulse.ConstantDetuning(BlackmanWaveform(t_pulse, np.pi/2.), 0, np.pi/2)"
+    "X_pulse = Pulse.ConstantDetuning(BlackmanWaveform(t_pulse, np.pi / 2.0), 0, 0)\n",
+    "Y_pulse = Pulse.ConstantDetuning(\n",
+    "    BlackmanWaveform(t_pulse, np.pi / 2.0), 0, -np.pi / 2\n",
+    ")\n",
+    "mX_pulse = Pulse.ConstantDetuning(\n",
+    "    BlackmanWaveform(t_pulse, np.pi / 2.0), 0, np.pi\n",
+    ")\n",
+    "mY_pulse = Pulse.ConstantDetuning(\n",
+    "    BlackmanWaveform(t_pulse, np.pi / 2.0), 0, np.pi / 2\n",
+    ")"
    ]
   },
   {
@@ -72,19 +78,19 @@
    "outputs": [],
    "source": [
     "def Floquet_XXZ_cycles(n_cycles, tau_1, tau_2, t_pulse):\n",
-    "    t_half = t_pulse/2.\n",
-    "    tau_3 = tau_2 \n",
-    "    tc = 4*tau_2 + 2*tau_1\n",
+    "    t_half = t_pulse / 2.0\n",
+    "    tau_3 = tau_2\n",
+    "    tc = 4 * tau_2 + 2 * tau_1\n",
     "    for _ in range(n_cycles):\n",
-    "        seq.delay(tau_1-t_half, 'MW')\n",
-    "        seq.add(X_pulse, 'MW')\n",
-    "        seq.delay(tau_2-2*t_half, 'MW')\n",
-    "        seq.add(mY_pulse, 'MW')\n",
-    "        seq.delay(2*tau_3-2*t_half, 'MW')\n",
-    "        seq.add(Y_pulse, 'MW')\n",
-    "        seq.delay(tau_2-2*t_half, 'MW')\n",
-    "        seq.add(mX_pulse, 'MW')\n",
-    "        seq.delay(tau_1-t_half, 'MW')"
+    "        seq.delay(tau_1 - t_half, \"MW\")\n",
+    "        seq.add(X_pulse, \"MW\")\n",
+    "        seq.delay(tau_2 - 2 * t_half, \"MW\")\n",
+    "        seq.add(mY_pulse, \"MW\")\n",
+    "        seq.delay(2 * tau_3 - 2 * t_half, \"MW\")\n",
+    "        seq.add(Y_pulse, \"MW\")\n",
+    "        seq.delay(tau_2 - 2 * t_half, \"MW\")\n",
+    "        seq.add(mX_pulse, \"MW\")\n",
+    "        seq.delay(tau_1 - t_half, \"MW\")"
    ]
   },
   {
@@ -106,18 +112,18 @@
     "reg = Register(qubits)\n",
     "\n",
     "seq = Sequence(reg, MockDevice)\n",
-    "seq.declare_channel('MW', 'mw_global')\n",
-    "seq.set_magnetic_field(0., 0., 1.)\n",
+    "seq.declare_channel(\"MW\", \"mw_global\")\n",
+    "seq.set_magnetic_field(0.0, 0.0, 1.0)\n",
     "\n",
     "tc = 300\n",
-    "seq.delay(3 * tc, 'MW')\n",
-    "Floquet_XXZ_cycles(4, tc/6., tc/6., t_pulse)\n",
-    "seq.delay(6 * tc, 'MW')  \n",
+    "seq.delay(3 * tc, \"MW\")\n",
+    "Floquet_XXZ_cycles(4, tc / 6.0, tc / 6.0, t_pulse)\n",
+    "seq.delay(6 * tc, \"MW\")\n",
     "\n",
     "# Here are our evaluation times\n",
-    "t_list= []\n",
+    "t_list = []\n",
     "for p in range(13):\n",
-    "    t_list.append(tc/1000.*p)"
+    "    t_list.append(tc / 1000.0 * p)"
    ]
   },
   {
@@ -154,7 +160,7 @@
    "outputs": [],
    "source": [
     "sim = Simulation(seq, sampling_rate=1.0, config=None, evaluation_times=t_list)\n",
-    "psi_y = (qutip.basis(2, 0)+1j*qutip.basis(2, 1)).unit()\n",
+    "psi_y = (qutip.basis(2, 0) + 1j * qutip.basis(2, 1)).unit()\n",
     "sim.initial_state = qutip.tensor(psi_y, psi_y)\n",
     "res = sim.run()"
    ]
@@ -167,7 +173,7 @@
    "source": [
     "sy = qutip.sigmay()\n",
     "Id = qutip.qeye(2)\n",
-    "Sigma_y = (qutip.tensor(sy, Id)+qutip.tensor(Id, sy))/2.\n",
+    "Sigma_y = (qutip.tensor(sy, Id) + qutip.tensor(Id, sy)) / 2.0\n",
     "Sigma_y_res = res.expect([Sigma_y])"
    ]
   },
@@ -194,13 +200,13 @@
     "\n",
     "# Showing the Hamiltonian engineering period.\n",
     "line1 = 0.9\n",
-    "line2 = 2.1 \n",
-    "plt.axvspan(line1, line2, alpha=.1, color='grey')\n",
-    "plt.text(1., 0.5, r\"$H_{XX} \\to H_{XXX}$\", fontsize=14)\n",
+    "line2 = 2.1\n",
+    "plt.axvspan(line1, line2, alpha=0.1, color=\"grey\")\n",
+    "plt.text(1.0, 0.5, r\"$H_{XX} \\to H_{XXX}$\", fontsize=14)\n",
     "\n",
-    "plt.plot(sim.evaluation_times, Sigma_y_res[0], 'o')\n",
+    "plt.plot(sim.evaluation_times, Sigma_y_res[0], \"o\")\n",
     "plt.xlabel(r\"Time [µs]\", fontsize=16)\n",
-    "plt.ylabel(fr'$ (\\langle \\sigma_1^y + \\sigma_2^y \\rangle)/2$', fontsize=16)\n",
+    "plt.ylabel(rf\"$ (\\langle \\sigma_1^y + \\sigma_2^y \\rangle)/2$\", fontsize=16)\n",
     "plt.show()"
    ]
   },
@@ -227,14 +233,14 @@
    "outputs": [],
    "source": [
     "def Floquet_XX2Z_cycles(n_cycles, t_pulse):\n",
-    "    t_half = t_pulse/2.\n",
-    "    tau_3 = tau_2 = tc/4.\n",
+    "    t_half = t_pulse / 2.0\n",
+    "    tau_3 = tau_2 = tc / 4.0\n",
     "    for _ in range(n_cycles):\n",
-    "        seq.delay(tau_2-t_half, 'MW')\n",
-    "        seq.add(mY_pulse, 'MW')\n",
-    "        seq.delay(2*tau_3-2*t_half, 'MW')\n",
-    "        seq.add(Y_pulse, 'MW')\n",
-    "        seq.delay(tau_2-t_half, 'MW')"
+    "        seq.delay(tau_2 - t_half, \"MW\")\n",
+    "        seq.add(mY_pulse, \"MW\")\n",
+    "        seq.delay(2 * tau_3 - 2 * t_half, \"MW\")\n",
+    "        seq.add(Y_pulse, \"MW\")\n",
+    "        seq.delay(tau_2 - t_half, \"MW\")"
    ]
   },
   {
@@ -244,10 +250,10 @@
    "outputs": [],
    "source": [
     "N_at = 10\n",
-    "# Number of Floquet cycles \n",
-    "N_cycles = 20 \n",
+    "# Number of Floquet cycles\n",
+    "N_cycles = 20\n",
     "# In the following, we will take 1000 projective measurements of the system at the final time.\n",
-    "N_samples = 1000 "
+    "N_samples = 1000"
    ]
   },
   {
@@ -265,8 +271,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "initial_pi_pulse = Pulse.ConstantDetuning(BlackmanWaveform(2*t_pulse, np.pi), 0, 0)\n",
-    "masked_indices = np.arange(N_at//2)"
+    "initial_pi_pulse = Pulse.ConstantDetuning(\n",
+    "    BlackmanWaveform(2 * t_pulse, np.pi), 0, 0\n",
+    ")\n",
+    "masked_indices = np.arange(N_at // 2)"
    ]
   },
   {
@@ -279,25 +287,32 @@
     "reg = Register.rectangle(1, N_at, 10)\n",
     "magnetizations_obc = np.zeros((N_at, N_cycles), dtype=float)\n",
     "correl_obc = np.zeros(N_cycles, dtype=float)\n",
-    "for m in range(N_cycles):     # Runtime close to 2 min!\n",
+    "for m in range(N_cycles):  # Runtime close to 2 min!\n",
     "    seq = Sequence(reg, MockDevice)\n",
-    "    seq.declare_channel('MW', 'mw_global')\n",
-    "    seq.set_magnetic_field(0., 0., 1.)\n",
+    "    seq.declare_channel(\"MW\", \"mw_global\")\n",
+    "    seq.set_magnetic_field(0.0, 0.0, 1.0)\n",
     "    # Configure the SLM mask, that will prevent the masked qubits from interacting with the first global \\pi pulse\n",
     "    seq.config_slm_mask(masked_indices)\n",
-    "    seq.add(initial_pi_pulse, 'MW')\n",
-    "    seq.add(X_pulse, 'MW')\n",
+    "    seq.add(initial_pi_pulse, \"MW\")\n",
+    "    seq.add(X_pulse, \"MW\")\n",
     "    Floquet_XX2Z_cycles(m, t_pulse)\n",
-    "    seq.add(mX_pulse, 'MW')\n",
+    "    seq.add(mX_pulse, \"MW\")\n",
     "    sim = Simulation(seq)\n",
     "    res = sim.run()\n",
     "    samples = res.sample_final_state(N_samples)\n",
-    "    correl = 0.\n",
+    "    correl = 0.0\n",
     "    for key, value in samples.items():\n",
     "        for j in range(N_at):\n",
-    "            correl -= (2*float(key[j])-1)*(2*float(key[(j+1)%N_at])-1)*value/N_samples \n",
-    "            magnetizations_obc[j][m] += (2*float(key[j])-1)*value/N_samples\n",
-    "    correl_obc[m] = N_at/2+correl/2"
+    "            correl -= (\n",
+    "                (2 * float(key[j]) - 1)\n",
+    "                * (2 * float(key[(j + 1) % N_at]) - 1)\n",
+    "                * value\n",
+    "                / N_samples\n",
+    "            )\n",
+    "            magnetizations_obc[j][m] += (\n",
+    "                (2 * float(key[j]) - 1) * value / N_samples\n",
+    "            )\n",
+    "    correl_obc[m] = N_at / 2 + correl / 2"
    ]
   },
   {
@@ -307,29 +322,49 @@
    "outputs": [],
    "source": [
     "# Circular geometry\n",
-    "coords = 10.*N_at/(2*np.pi)*np.array([(np.cos(theta*2*np.pi/N_at), np.sin(theta*2*np.pi/N_at)) for theta in range(N_at)])\n",
+    "coords = (\n",
+    "    10.0\n",
+    "    * N_at\n",
+    "    / (2 * np.pi)\n",
+    "    * np.array(\n",
+    "        [\n",
+    "            (\n",
+    "                np.cos(theta * 2 * np.pi / N_at),\n",
+    "                np.sin(theta * 2 * np.pi / N_at),\n",
+    "            )\n",
+    "            for theta in range(N_at)\n",
+    "        ]\n",
+    "    )\n",
+    ")\n",
     "reg = Register.from_coordinates(coords)\n",
     "\n",
     "magnetizations_pbc = np.zeros((N_at, N_cycles), dtype=float)\n",
     "correl_pbc = np.zeros(N_cycles, dtype=float)\n",
-    "for m in range(N_cycles):     # Runtime close to 2 min!\n",
+    "for m in range(N_cycles):  # Runtime close to 2 min!\n",
     "    seq = Sequence(reg, MockDevice)\n",
-    "    seq.declare_channel('MW', 'mw_global')\n",
-    "    seq.set_magnetic_field(0., 0., 1.)\n",
+    "    seq.declare_channel(\"MW\", \"mw_global\")\n",
+    "    seq.set_magnetic_field(0.0, 0.0, 1.0)\n",
     "    seq.config_slm_mask(masked_indices)\n",
-    "    seq.add(initial_pi_pulse, 'MW')\n",
-    "    seq.add(X_pulse, 'MW')\n",
+    "    seq.add(initial_pi_pulse, \"MW\")\n",
+    "    seq.add(X_pulse, \"MW\")\n",
     "    Floquet_XX2Z_cycles(m, t_pulse)\n",
-    "    seq.add(mX_pulse, 'MW')\n",
+    "    seq.add(mX_pulse, \"MW\")\n",
     "    sim = Simulation(seq)\n",
     "    res = sim.run()\n",
     "    samples = res.sample_final_state(N_samples)\n",
-    "    correl = 0.\n",
+    "    correl = 0.0\n",
     "    for key, value in samples.items():\n",
     "        for j in range(N_at):\n",
-    "            correl -= (2*float(key[j])-1)*(2*float(key[(j+1)%N_at])-1)*value/N_samples \n",
-    "            magnetizations_pbc[j][m] += (2*float(key[j])-1)*value/N_samples\n",
-    "    correl_pbc[m] = N_at/2+correl/2"
+    "            correl -= (\n",
+    "                (2 * float(key[j]) - 1)\n",
+    "                * (2 * float(key[(j + 1) % N_at]) - 1)\n",
+    "                * value\n",
+    "                / N_samples\n",
+    "            )\n",
+    "            magnetizations_pbc[j][m] += (\n",
+    "                (2 * float(key[j]) - 1) * value / N_samples\n",
+    "            )\n",
+    "    correl_pbc[m] = N_at / 2 + correl / 2"
    ]
   },
   {
@@ -370,22 +405,22 @@
     }
    ],
    "source": [
-    "fig, ax = plt.subplots(1,1)\n",
-    "img = ax.imshow(magnetizations_obc, cmap=plt.get_cmap('RdBu'))\n",
-    "plt.title('OBC',fontsize=16)\n",
-    "ax.set_xlabel('Cycle',fontsize=16)\n",
-    "ax.set_ylabel('Atom number',fontsize=16)\n",
+    "fig, ax = plt.subplots(1, 1)\n",
+    "img = ax.imshow(magnetizations_obc, cmap=plt.get_cmap(\"RdBu\"))\n",
+    "plt.title(\"OBC\", fontsize=16)\n",
+    "ax.set_xlabel(\"Cycle\", fontsize=16)\n",
+    "ax.set_ylabel(\"Atom number\", fontsize=16)\n",
     "cbar = fig.colorbar(img, shrink=0.7)\n",
-    "cbar.set_label(r'$\\langle \\sigma^z \\rangle$', fontsize=16)\n",
+    "cbar.set_label(r\"$\\langle \\sigma^z \\rangle$\", fontsize=16)\n",
     "\n",
     "\n",
-    "fig, ax = plt.subplots(1,1)\n",
-    "img = ax.imshow(magnetizations_pbc, cmap=plt.get_cmap('RdBu'))\n",
-    "plt.title('PBC',fontsize=16)\n",
-    "ax.set_xlabel('Cycle',fontsize=16)\n",
-    "ax.set_ylabel('Atom number',fontsize=16)\n",
+    "fig, ax = plt.subplots(1, 1)\n",
+    "img = ax.imshow(magnetizations_pbc, cmap=plt.get_cmap(\"RdBu\"))\n",
+    "plt.title(\"PBC\", fontsize=16)\n",
+    "ax.set_xlabel(\"Cycle\", fontsize=16)\n",
+    "ax.set_ylabel(\"Atom number\", fontsize=16)\n",
     "cbar = fig.colorbar(img, shrink=0.7)\n",
-    "cbar.set_label(r'$\\langle \\sigma^z \\rangle$', fontsize=16)"
+    "cbar.set_label(r\"$\\langle \\sigma^z \\rangle$\", fontsize=16)"
    ]
   },
   {
@@ -424,12 +459,15 @@
     }
    ],
    "source": [
-    "fig, ax = plt.subplots(1,1)\n",
-    "plt.title(r'Evolution of $\\langle N_{flip} \\rangle$ in time for OBC (red) and PBC (blue).', fontsize=16)\n",
-    "ax.set_xlabel('Cycle',fontsize=16)\n",
-    "ax.set_ylabel(r'$\\langle N_{flip} \\rangle$',fontsize=14)\n",
-    "ax.plot(correl_pbc,'--o',color='blue')\n",
-    "ax.plot(correl_obc,'--o',color='red')"
+    "fig, ax = plt.subplots(1, 1)\n",
+    "plt.title(\n",
+    "    r\"Evolution of $\\langle N_{flip} \\rangle$ in time for OBC (red) and PBC (blue).\",\n",
+    "    fontsize=16,\n",
+    ")\n",
+    "ax.set_xlabel(\"Cycle\", fontsize=16)\n",
+    "ax.set_ylabel(r\"$\\langle N_{flip} \\rangle$\", fontsize=14)\n",
+    "ax.plot(correl_pbc, \"--o\", color=\"blue\")\n",
+    "ax.plot(correl_obc, \"--o\", color=\"red\")"
    ]
   },
   {
@@ -445,27 +483,29 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "N_cycles = 26  \n",
+    "N_cycles = 26\n",
     "magnetizations_pbc = np.zeros((N_at, N_cycles), dtype=float)\n",
     "samples_evol = []\n",
     "masked_indices = [0, 1, 2]\n",
-    "for m in range(N_cycles): # Runtime close to 4 min!\n",
+    "for m in range(N_cycles):  # Runtime close to 4 min!\n",
     "    seq = Sequence(reg, MockDevice)\n",
-    "    seq.set_magnetic_field(0., 0., 1.)\n",
-    "    seq.declare_channel('MW', 'mw_global')\n",
+    "    seq.set_magnetic_field(0.0, 0.0, 1.0)\n",
+    "    seq.declare_channel(\"MW\", \"mw_global\")\n",
     "    seq.config_slm_mask(masked_indices)\n",
-    "    seq.add(initial_pi_pulse, 'MW')\n",
-    "    seq.add(X_pulse, 'MW')\n",
+    "    seq.add(initial_pi_pulse, \"MW\")\n",
+    "    seq.add(X_pulse, \"MW\")\n",
     "    Floquet_XX2Z_cycles(m, t_pulse)\n",
-    "    seq.add(mX_pulse, 'MW')\n",
+    "    seq.add(mX_pulse, \"MW\")\n",
     "    sim = Simulation(seq)\n",
     "    res = sim.run()\n",
     "    samples = res.sample_final_state(N_samples)\n",
     "    samples_evol.append(samples)\n",
-    "    correl = 0.\n",
+    "    correl = 0.0\n",
     "    for key, value in samples.items():\n",
     "        for j in range(N_at):\n",
-    "            magnetizations_pbc[j][m] += (2*float(key[j])-1)*value/N_samples"
+    "            magnetizations_pbc[j][m] += (\n",
+    "                (2 * float(key[j]) - 1) * value / N_samples\n",
+    "            )"
    ]
   },
   {
@@ -487,12 +527,12 @@
     }
    ],
    "source": [
-    "fig, ax = plt.subplots(1,1)\n",
-    "img = ax.imshow(magnetizations_pbc, cmap=plt.get_cmap('RdBu'))\n",
-    "ax.set_xlabel('Cycle', fontsize=16)\n",
-    "ax.set_ylabel('Atom number', fontsize=16)\n",
+    "fig, ax = plt.subplots(1, 1)\n",
+    "img = ax.imshow(magnetizations_pbc, cmap=plt.get_cmap(\"RdBu\"))\n",
+    "ax.set_xlabel(\"Cycle\", fontsize=16)\n",
+    "ax.set_ylabel(\"Atom number\", fontsize=16)\n",
     "cbar = fig.colorbar(img)\n",
-    "cbar.set_label(r'$\\langle \\sigma^z \\rangle$', fontsize=16)"
+    "cbar.set_label(r\"$\\langle \\sigma^z \\rangle$\", fontsize=16)"
    ]
   },
   {
@@ -665,14 +705,33 @@
     }
    ],
    "source": [
-    "dw_preserved = ['0001111111', '1000111111', '1100011111', '1110001111', \n",
-    "                 '1111000111', '1111100011', '1111110001', '1111111000', '0111111100', '0011111110']\n",
+    "dw_preserved = [\n",
+    "    \"0001111111\",\n",
+    "    \"1000111111\",\n",
+    "    \"1100011111\",\n",
+    "    \"1110001111\",\n",
+    "    \"1111000111\",\n",
+    "    \"1111100011\",\n",
+    "    \"1111110001\",\n",
+    "    \"1111111000\",\n",
+    "    \"0111111100\",\n",
+    "    \"0011111110\",\n",
+    "]\n",
     "\n",
-    "for n_cycle in [2*k for k in range(int(N_cycles/2))]:    # Runtime close to 2 min !\n",
-    "    color_dict = {key: 'red' if key in dw_preserved else 'black' for key in samples_evol[n_cycle]}\n",
+    "for n_cycle in [\n",
+    "    2 * k for k in range(int(N_cycles / 2))\n",
+    "]:  # Runtime close to 2 min !\n",
+    "    color_dict = {\n",
+    "        key: \"red\" if key in dw_preserved else \"black\"\n",
+    "        for key in samples_evol[n_cycle]\n",
+    "    }\n",
     "    plt.figure(figsize=(16, 5))\n",
-    "    plt.title(r'Cycle $= {}$'.format(n_cycle), fontsize=18)\n",
-    "    plt.bar(samples_evol[n_cycle].keys(), samples_evol[n_cycle].values(), color=color_dict.values())\n",
+    "    plt.title(r\"Cycle $= {}$\".format(n_cycle), fontsize=18)\n",
+    "    plt.bar(\n",
+    "        samples_evol[n_cycle].keys(),\n",
+    "        samples_evol[n_cycle].values(),\n",
+    "        color=color_dict.values(),\n",
+    "    )\n",
     "    plt.xlabel(\"bitstrings\", fontsize=16)\n",
     "    plt.ylabel(\"counts\", fontsize=16)\n",
     "    plt.xticks(rotation=90)\n",
diff --git a/tutorials/quantum_simulation/Preparing state with antiferromagnetic order in the Ising model.ipynb b/tutorials/quantum_simulation/Preparing state with antiferromagnetic order in the Ising model.ipynb
index 6700ab7b0..e1a2f0803 100644
--- a/tutorials/quantum_simulation/Preparing state with antiferromagnetic order in the Ising model.ipynb
+++ b/tutorials/quantum_simulation/Preparing state with antiferromagnetic order in the Ising model.ipynb
@@ -59,7 +59,7 @@
    "outputs": [],
    "source": [
     "# Parameters in rad/µs and ns\n",
-    "Omega_max = 2.3 * 2*np.pi \n",
+    "Omega_max = 2.3 * 2 * np.pi\n",
     "U = Omega_max / 2.3\n",
     "\n",
     "delta_0 = -6 * U\n",
@@ -67,13 +67,13 @@
     "\n",
     "t_rise = 252\n",
     "t_fall = 500\n",
-    "t_sweep = (delta_f - delta_0)/(2 * np.pi * 10) * 1000\n",
+    "t_sweep = (delta_f - delta_0) / (2 * np.pi * 10) * 1000\n",
     "\n",
     "R_interatomic = Chadoq2.rydberg_blockade_radius(U)\n",
     "\n",
     "N_side = 3\n",
-    "reg = Register.square(N_side, R_interatomic, prefix='q')\n",
-    "print(f'Interatomic Radius is: {R_interatomic}µm.')\n",
+    "reg = Register.square(N_side, R_interatomic, prefix=\"q\")\n",
+    "print(f\"Interatomic Radius is: {R_interatomic}µm.\")\n",
     "reg.draw()"
    ]
   },
@@ -97,9 +97,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "rise = Pulse.ConstantDetuning(RampWaveform(t_rise, 0., Omega_max), delta_0, 0.)\n",
-    "sweep = Pulse.ConstantAmplitude(Omega_max, RampWaveform(t_sweep, delta_0, delta_f), 0.)\n",
-    "fall = Pulse.ConstantDetuning(RampWaveform(t_fall, Omega_max, 0.), delta_f, 0.)"
+    "rise = Pulse.ConstantDetuning(\n",
+    "    RampWaveform(t_rise, 0.0, Omega_max), delta_0, 0.0\n",
+    ")\n",
+    "sweep = Pulse.ConstantAmplitude(\n",
+    "    Omega_max, RampWaveform(t_sweep, delta_0, delta_f), 0.0\n",
+    ")\n",
+    "fall = Pulse.ConstantDetuning(\n",
+    "    RampWaveform(t_fall, Omega_max, 0.0), delta_f, 0.0\n",
+    ")"
    ]
   },
   {
@@ -109,11 +115,11 @@
    "outputs": [],
    "source": [
     "seq = Sequence(reg, Chadoq2)\n",
-    "seq.declare_channel('ising', 'rydberg_global')\n",
+    "seq.declare_channel(\"ising\", \"rydberg_global\")\n",
     "\n",
-    "seq.add(rise, 'ising')\n",
-    "seq.add(sweep, 'ising')\n",
-    "seq.add(fall, 'ising')\n",
+    "seq.add(rise, \"ising\")\n",
+    "seq.add(sweep, \"ising\")\n",
+    "seq.add(fall, \"ising\")\n",
     "\n",
     "seq.draw()"
    ]
@@ -140,25 +146,25 @@
    "source": [
     "delta = []\n",
     "omega = []\n",
-    "for x in seq._schedule['ising']:\n",
-    "    if isinstance(x.type,Pulse):\n",
+    "for x in seq._schedule[\"ising\"]:\n",
+    "    if isinstance(x.type, Pulse):\n",
     "        omega += list(x.type.amplitude.samples / U)\n",
     "        delta += list(x.type.detuning.samples / U)\n",
-    "        \n",
+    "\n",
     "fig, ax = plt.subplots()\n",
-    "ax.grid(True, which='both')\n",
+    "ax.grid(True, which=\"both\")\n",
     "\n",
     "ax.set_ylabel(r\"$\\hbar\\delta(t)/U$\", fontsize=16)\n",
     "ax.set_xlabel(r\"$\\hbar\\Omega(t)/U$\", fontsize=16)\n",
     "ax.set_xlim(0, 3)\n",
-    "ax.axhline(y=0, color='k')\n",
-    "ax.axvline(x=0, color='k')\n",
+    "ax.axhline(y=0, color=\"k\")\n",
+    "ax.axvline(x=0, color=\"k\")\n",
     "\n",
     "y = np.arange(0.0, 6, 0.01)\n",
-    "x = 1.522 * (1 - 0.25 * (y - 2)**2)\n",
+    "x = 1.522 * (1 - 0.25 * (y - 2) ** 2)\n",
     "ax.fill_between(x, y, alpha=0.4)\n",
     "\n",
-    "ax.plot(omega,delta, 'red', lw=2)\n",
+    "ax.plot(omega, delta, \"red\", lw=2)\n",
     "plt.show()"
    ]
   },
@@ -208,9 +214,9 @@
    "source": [
     "count = results.sample_final_state()\n",
     "\n",
-    "most_freq = {k:v for k,v in count.items() if v>10}\n",
+    "most_freq = {k: v for k, v in count.items() if v > 10}\n",
     "plt.bar(list(most_freq.keys()), list(most_freq.values()))\n",
-    "plt.xticks(rotation='vertical')\n",
+    "plt.xticks(rotation=\"vertical\")\n",
     "plt.show()"
    ]
   },
@@ -227,8 +233,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def occupation(j,N):\n",
-    "    up = qutip.basis(2,0)\n",
+    "def occupation(j, N):\n",
+    "    up = qutip.basis(2, 0)\n",
     "    prod = [qutip.qeye(2) for _ in range(N)]\n",
     "    prod[j] = up * up.dag()\n",
     "    return qutip.tensor(prod)"
@@ -240,7 +246,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "occup_list = [occupation(j, N_side*N_side) for j in range(N_side*N_side)]"
+    "occup_list = [occupation(j, N_side * N_side) for j in range(N_side * N_side)]"
    ]
   },
   {
@@ -260,8 +266,8 @@
     "    corr_pairs = []\n",
     "    for i, qi in enumerate(register.qubits):\n",
     "        for j, qj in enumerate(register.qubits):\n",
-    "            r_ij = register.qubits[qi]-register.qubits[qj]\n",
-    "            distance = np.linalg.norm(r_ij - R_interatomic*np.array([k, l]))\n",
+    "            r_ij = register.qubits[qi] - register.qubits[qj]\n",
+    "            distance = np.linalg.norm(r_ij - R_interatomic * np.array([k, l]))\n",
     "            if distance < 1:\n",
     "                corr_pairs.append([i, j])\n",
     "    return corr_pairs"
@@ -283,22 +289,27 @@
     "def get_corr_function(k, l, reg, R_interatomic, state):\n",
     "    N_qubits = len(reg.qubits)\n",
     "    corr_pairs = get_corr_pairs(k, l, reg, R_interatomic)\n",
-    "    \n",
+    "\n",
     "    operators = [occupation(j, N_qubits) for j in range(N_qubits)]\n",
     "    covariance = 0\n",
     "    for qi, qj in corr_pairs:\n",
-    "        covariance += qutip.expect(operators[qi]*operators[qj], state)\n",
-    "        covariance -= qutip.expect(operators[qi], state)*qutip.expect(operators[qj], state)\n",
-    "    return covariance/len(corr_pairs)\n",
-    "        \n",
+    "        covariance += qutip.expect(operators[qi] * operators[qj], state)\n",
+    "        covariance -= qutip.expect(operators[qi], state) * qutip.expect(\n",
+    "            operators[qj], state\n",
+    "        )\n",
+    "    return covariance / len(corr_pairs)\n",
+    "\n",
+    "\n",
     "def get_full_corr_function(reg, state):\n",
     "    N_qubits = len(reg.qubits)\n",
-    "    \n",
+    "\n",
     "    correlation_function = {}\n",
     "    N_side = int(np.sqrt(N_qubits))\n",
-    "    for k in range(-N_side+1, N_side):\n",
-    "        for l in range(-N_side+1, N_side):\n",
-    "            correlation_function[(k, l)] = get_corr_function(k, l, reg, R_interatomic, state)\n",
+    "    for k in range(-N_side + 1, N_side):\n",
+    "        for l in range(-N_side + 1, N_side):\n",
+    "            correlation_function[(k, l)] = get_corr_function(\n",
+    "                k, l, reg, R_interatomic, state\n",
+    "            )\n",
     "    return correlation_function"
    ]
   },
@@ -326,12 +337,14 @@
    "outputs": [],
    "source": [
     "expected_corr_function = {}\n",
-    "xi = 1 # Estimated Correlation Length\n",
-    "for k in range(-N_side+1,N_side):\n",
-    "    for l in range(-N_side+1,N_side):\n",
+    "xi = 1  # Estimated Correlation Length\n",
+    "for k in range(-N_side + 1, N_side):\n",
+    "    for l in range(-N_side + 1, N_side):\n",
     "        kk = np.abs(k)\n",
     "        ll = np.abs(l)\n",
-    "        expected_corr_function[(k, l)] = (-1)**(kk + ll) * np.exp(-np.sqrt(k**2 + l**2)/xi)"
+    "        expected_corr_function[(k, l)] = (-1) ** (kk + ll) * np.exp(\n",
+    "            -np.sqrt(k**2 + l**2) / xi\n",
+    "        )"
    ]
   },
   {
@@ -342,22 +355,34 @@
    },
    "outputs": [],
    "source": [
-    "A = 4*np.reshape(list(correlation_function.values()), (2*N_side-1, 2*N_side-1))\n",
-    "A = A/np.max(A)\n",
-    "B = np.reshape(list(expected_corr_function.values()), (2*N_side-1, 2*N_side-1))\n",
-    "B = B*np.max(A)\n",
-    "\n",
-    "for i, M in enumerate([A.copy(),B.copy()]):\n",
-    "    M[N_side-1, N_side-1] = None\n",
-    "    plt.figure(figsize=(3.5,3.5))\n",
-    "    plt.imshow(M, cmap='coolwarm', vmin=-.6, vmax=.6)\n",
-    "    plt.xticks(range(len(M)), [f'{x}' for x in range(-N_side + 1, N_side)])\n",
-    "    plt.xlabel(r'$\\mathscr{k}$', fontsize=22)\n",
-    "    plt.yticks(range(len(M)), [f'{-y}' for y in range(-N_side + 1, N_side)])\n",
-    "    plt.ylabel(r'$\\mathscr{l}$', rotation=0, fontsize=22, labelpad=10)\n",
+    "A = 4 * np.reshape(\n",
+    "    list(correlation_function.values()), (2 * N_side - 1, 2 * N_side - 1)\n",
+    ")\n",
+    "A = A / np.max(A)\n",
+    "B = np.reshape(\n",
+    "    list(expected_corr_function.values()), (2 * N_side - 1, 2 * N_side - 1)\n",
+    ")\n",
+    "B = B * np.max(A)\n",
+    "\n",
+    "for i, M in enumerate([A.copy(), B.copy()]):\n",
+    "    M[N_side - 1, N_side - 1] = None\n",
+    "    plt.figure(figsize=(3.5, 3.5))\n",
+    "    plt.imshow(M, cmap=\"coolwarm\", vmin=-0.6, vmax=0.6)\n",
+    "    plt.xticks(range(len(M)), [f\"{x}\" for x in range(-N_side + 1, N_side)])\n",
+    "    plt.xlabel(r\"$\\mathscr{k}$\", fontsize=22)\n",
+    "    plt.yticks(range(len(M)), [f\"{-y}\" for y in range(-N_side + 1, N_side)])\n",
+    "    plt.ylabel(r\"$\\mathscr{l}$\", rotation=0, fontsize=22, labelpad=10)\n",
     "    plt.colorbar(fraction=0.047, pad=0.02)\n",
-    "    if i == 0 :plt.title(r'$4\\times\\.g^{(2)}(\\mathscr{k},\\mathscr{l})$ after simulation', fontsize=14)\n",
-    "    if i == 1 :plt.title(r'Exponential $g^{(2)}(\\mathscr{k},\\mathscr{l})$ expected', fontsize=14)\n",
+    "    if i == 0:\n",
+    "        plt.title(\n",
+    "            r\"$4\\times\\.g^{(2)}(\\mathscr{k},\\mathscr{l})$ after simulation\",\n",
+    "            fontsize=14,\n",
+    "        )\n",
+    "    if i == 1:\n",
+    "        plt.title(\n",
+    "            r\"Exponential $g^{(2)}(\\mathscr{k},\\mathscr{l})$ expected\",\n",
+    "            fontsize=14,\n",
+    "        )\n",
     "    plt.show()"
    ]
   },
@@ -411,13 +436,17 @@
     "    N_side = int(np.sqrt(N_qubits))\n",
     "\n",
     "    st_fac = 0\n",
-    "    for k in range(-N_side+1, N_side):\n",
-    "        for l in range(-N_side+1, N_side):\n",
+    "    for k in range(-N_side + 1, N_side):\n",
+    "        for l in range(-N_side + 1, N_side):\n",
     "            kk = np.abs(k)\n",
     "            ll = np.abs(l)\n",
     "            if not (k == 0 and l == 0):\n",
-    "                st_fac += 4 * (-1)**(kk + ll) * get_corr_function(k, l, reg, R_interatomic, state)\n",
-    "    return st_fac "
+    "                st_fac += (\n",
+    "                    4\n",
+    "                    * (-1) ** (kk + ll)\n",
+    "                    * get_corr_function(k, l, reg, R_interatomic, state)\n",
+    "                )\n",
+    "    return st_fac"
    ]
   },
   {
@@ -426,36 +455,44 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def calculate_neel(det, N, Omega_max = 2.3 * 2 * np.pi):\n",
-    "    #Setup:\n",
+    "def calculate_neel(det, N, Omega_max=2.3 * 2 * np.pi):\n",
+    "    # Setup:\n",
     "    U = Omega_max / 2.3\n",
     "    delta_0 = -6 * U\n",
-    "    delta_f = det * U              \n",
-    "    \n",
+    "    delta_f = det * U\n",
+    "\n",
     "    t_rise = 252\n",
     "    t_fall = 500\n",
-    "    t_sweep = int((delta_f - delta_0)/(2 * np.pi * 10) * 1000)\n",
-    "    t_sweep += 4 - t_sweep % 4 # To be a multiple of the clock period of Chadoq2 (4ns)\n",
-    "    \n",
-    "    R_interatomic = Chadoq2.rydberg_blockade_radius(U) \n",
+    "    t_sweep = int((delta_f - delta_0) / (2 * np.pi * 10) * 1000)\n",
+    "    t_sweep += (\n",
+    "        4 - t_sweep % 4\n",
+    "    )  # To be a multiple of the clock period of Chadoq2 (4ns)\n",
+    "\n",
+    "    R_interatomic = Chadoq2.rydberg_blockade_radius(U)\n",
     "    reg = Register.rectangle(N, N, R_interatomic)\n",
     "\n",
-    "    #Pulse Sequence\n",
-    "    rise = Pulse.ConstantDetuning(RampWaveform(t_rise, 0., Omega_max), delta_0, 0.)\n",
-    "    sweep = Pulse.ConstantAmplitude(Omega_max, RampWaveform(t_sweep, delta_0, delta_f), 0.)\n",
-    "    fall = Pulse.ConstantDetuning(RampWaveform(t_fall, Omega_max, 0.), delta_f, 0.)\n",
+    "    # Pulse Sequence\n",
+    "    rise = Pulse.ConstantDetuning(\n",
+    "        RampWaveform(t_rise, 0.0, Omega_max), delta_0, 0.0\n",
+    "    )\n",
+    "    sweep = Pulse.ConstantAmplitude(\n",
+    "        Omega_max, RampWaveform(t_sweep, delta_0, delta_f), 0.0\n",
+    "    )\n",
+    "    fall = Pulse.ConstantDetuning(\n",
+    "        RampWaveform(t_fall, Omega_max, 0.0), delta_f, 0.0\n",
+    "    )\n",
     "\n",
     "    seq = Sequence(reg, Chadoq2)\n",
-    "    seq.declare_channel('ising', 'rydberg_global')\n",
-    "    seq.add(rise, 'ising')\n",
-    "    seq.add(sweep, 'ising')\n",
-    "    seq.add(fall, 'ising')\n",
+    "    seq.declare_channel(\"ising\", \"rydberg_global\")\n",
+    "    seq.add(rise, \"ising\")\n",
+    "    seq.add(sweep, \"ising\")\n",
+    "    seq.add(fall, \"ising\")\n",
     "\n",
     "    simul = Simulation(seq, sampling_rate=0.2)\n",
     "    results = simul.run()\n",
-    "    \n",
+    "\n",
     "    final = results.states[-1]\n",
-    "    return get_neel_structure_factor(reg, R_interatomic, final)                                                                                                                                                                      "
+    "    return get_neel_structure_factor(reg, R_interatomic, final)"
    ]
   },
   {
@@ -467,19 +504,21 @@
    "outputs": [],
    "source": [
     "N_side = 3\n",
-    "occup_list = [occupation(j, N_side*N_side) for j in range(N_side*N_side)]\n",
+    "occup_list = [occupation(j, N_side * N_side) for j in range(N_side * N_side)]\n",
     "\n",
     "detunings = np.linspace(-1, 5, 20)\n",
-    "results=[]\n",
+    "results = []\n",
     "for det in detunings:\n",
-    "    print(f'Detuning = {np.round(det,3)} x 2π Mhz.')\n",
+    "    print(f\"Detuning = {np.round(det,3)} x 2π Mhz.\")\n",
     "    results.append(calculate_neel(det, N_side))\n",
-    "plt.xlabel(r'$\\hbar\\delta_{final}/U$')\n",
-    "plt.ylabel(r'Néel Structure Factor $S_{Neel}$')\n",
-    "plt.plot(detunings, results, 'o', ls='solid')\n",
+    "plt.xlabel(r\"$\\hbar\\delta_{final}/U$\")\n",
+    "plt.ylabel(r\"Néel Structure Factor $S_{Neel}$\")\n",
+    "plt.plot(detunings, results, \"o\", ls=\"solid\")\n",
     "plt.show()\n",
     "max_index = results.index(max(results))\n",
-    "print(f'Max S_Neel {np.round(max(results),2)} at detuning = {np.round(detunings[max_index],2)} x 2π Mhz.')"
+    "print(\n",
+    "    f\"Max S_Neel {np.round(max(results),2)} at detuning = {np.round(detunings[max_index],2)} x 2π Mhz.\"\n",
+    ")"
    ]
   }
  ],
diff --git a/tutorials/quantum_simulation/Shadow estimation for VQS.ipynb b/tutorials/quantum_simulation/Shadow estimation for VQS.ipynb
index d8ad3c4f2..b7bd17a06 100644
--- a/tutorials/quantum_simulation/Shadow estimation for VQS.ipynb
+++ b/tutorials/quantum_simulation/Shadow estimation for VQS.ipynb
@@ -115,13 +115,13 @@
     "num_qubits = 4\n",
     "zero_state = qutip.basis(2, 0).proj()\n",
     "one_state = qutip.basis(2, 1).proj()\n",
-    "hadamard = 1/np.sqrt(2) * qutip.Qobj([[1., 1.], [1., -1.]])\n",
-    "h_mul_phase = qutip.Qobj(np.array([[1., 1], [1.j, -1.j]])) / np.sqrt(2)\n",
+    "hadamard = 1 / np.sqrt(2) * qutip.Qobj([[1.0, 1.0], [1.0, -1.0]])\n",
+    "h_mul_phase = qutip.Qobj(np.array([[1.0, 1], [1.0j, -1.0j]])) / np.sqrt(2)\n",
     "unitary_ensemble = [hadamard, h_mul_phase, qutip.qeye(2)]\n",
     "\n",
-    "g = qutip.basis(2,1)\n",
-    "r = qutip.basis(2,0)\n",
-    "n = r*r.dag()\n",
+    "g = qutip.basis(2, 1)\n",
+    "r = qutip.basis(2, 0)\n",
+    "n = r * r.dag()\n",
     "\n",
     "sx = qutip.sigmax()\n",
     "sy = qutip.sigmay()\n",
@@ -175,11 +175,11 @@
    "source": [
     "def compute_shadow_size(delta, epsilon, observables):\n",
     "    \"\"\"Helper function.\n",
-    "    \n",
-    "    Computes both the number of shadows needed as well as the size of blocks needed \n",
+    "\n",
+    "    Computes both the number of shadows needed as well as the size of blocks needed\n",
     "    for the median_of_means method in order to approximate the expectation value of M\n",
     "    (linear) observables with additive error epsilon and fail probability delta.\n",
-    "    \n",
+    "\n",
     "    Args:\n",
     "        delta (float): Failure probability.\n",
     "        epsilon (float): Additive error on expectation values.\n",
@@ -190,8 +190,9 @@
     "    shadow_norm = (\n",
     "        lambda op: np.linalg.norm(\n",
     "            op - np.trace(op) / 2 ** int(np.log2(op.shape[0])), ord=np.inf\n",
-    "            ) ** 2\n",
     "        )\n",
+    "        ** 2\n",
+    "    )\n",
     "    # Theoretical number of shadows per cluster in the median of means procedure :\n",
     "    # N = 34 * max(shadow_norm(o) for o in observables) / epsilon ** 2\n",
     "    # We use N = 20 here to allow for quick simulation\n",
@@ -225,7 +226,9 @@
     "    unitary_ids = np.random.randint(0, 3, size=(shadow_size, num_qubits))\n",
     "    outcomes = []\n",
     "    for ns in range(shadow_size):\n",
-    "        unitmat = qutip.tensor([unitary_ensemble[unitary_ids[ns, i]] for i in range(num_qubits)])\n",
+    "        unitmat = qutip.tensor(\n",
+    "            [unitary_ensemble[unitary_ids[ns, i]] for i in range(num_qubits)]\n",
+    "        )\n",
     "        outcomes.append(measure_bitstring(unitmat.dag() * rho * unitmat))\n",
     "\n",
     "    # combine the computational basis outcomes and the sampled unitaries\n",
@@ -257,18 +260,18 @@
     "\n",
     "    Args:\n",
     "        outcome_ns: Bitstring at ns\n",
-    "        unitary_ids_ns: Rotation applied at ns. \n",
+    "        unitary_ids_ns: Rotation applied at ns.\n",
     "\n",
     "    Returns:\n",
     "        Reconstructed snapshot.\n",
     "    \"\"\"\n",
     "    state_list = []\n",
-    "    \n",
+    "\n",
     "    for k in range(num_qubits):\n",
     "        op = unitary_ensemble[unitary_ids_ns[k]]\n",
-    "        b = zero_state if outcome_ns[k] == '0' else one_state\n",
+    "        b = zero_state if outcome_ns[k] == \"0\" else one_state\n",
     "        state_list.append(3 * op * b * op.dag() - qutip.qeye(2))\n",
-    "    \n",
+    "\n",
     "    return qutip.tensor(state_list)"
    ]
   },
@@ -294,7 +297,7 @@
     "    # computing and saving mean per block\n",
     "    means = []\n",
     "    for block in range(K):\n",
-    "        states = [snap_list[i] for i in np.where(indic==block)[0]]\n",
+    "        states = [snap_list[i] for i in np.where(indic == block)[0]]\n",
     "        exp = qutip.expect(obs, states)\n",
     "        means.append(np.mean(exp))\n",
     "    return np.median(means)"
@@ -351,11 +354,20 @@
    "source": [
     "num_qubits = 2\n",
     "shadow_size = 10000\n",
-    "rho_1 = (qutip.tensor([qutip.basis(2,0), qutip.basis(2,0)]) + qutip.tensor([qutip.basis(2,0), qutip.basis(2,1)])).proj().unit()\n",
+    "rho_1 = (\n",
+    "    (\n",
+    "        qutip.tensor([qutip.basis(2, 0), qutip.basis(2, 0)])\n",
+    "        + qutip.tensor([qutip.basis(2, 0), qutip.basis(2, 1)])\n",
+    "    )\n",
+    "    .proj()\n",
+    "    .unit()\n",
+    ")\n",
     "print(\"Original density matrix :\")\n",
     "print(rho_1.full())\n",
     "outcomes, unitary_ids = calculate_classical_shadow(rho_1, shadow_size)\n",
-    "snapshots = [snapshot_state(outcomes[ns], unitary_ids[ns]) for ns in range(shadow_size)]\n",
+    "snapshots = [\n",
+    "    snapshot_state(outcomes[ns], unitary_ids[ns]) for ns in range(shadow_size)\n",
+    "]\n",
     "print(\"Shadow reconstruction :\")\n",
     "print(np.around(state_reconstruction(snapshots).full(), 2))\n",
     "\n",
@@ -363,7 +375,10 @@
     "shadow_sizes = [100, 1000, 2000, 5000, 10000]\n",
     "for i, shadow_size in enumerate(shadow_sizes):\n",
     "    outcomes, unitary_ids = calculate_classical_shadow(rho_1, shadow_size)\n",
-    "    snapshots = [snapshot_state(outcomes[ns], unitary_ids[ns]) for ns in range(shadow_size)]\n",
+    "    snapshots = [\n",
+    "        snapshot_state(outcomes[ns], unitary_ids[ns])\n",
+    "        for ns in range(shadow_size)\n",
+    "    ]\n",
     "    dist[i] = qutip.tracedist(state_reconstruction(snapshots), rho_1)\n",
     "num_qubits = 4"
    ]
@@ -447,7 +462,7 @@
     "    if end != -1:\n",
     "        o = o[:end]\n",
     "    for i, x in enumerate(o):\n",
-    "        if not(x == p[i] or x == \"1\"):\n",
+    "        if not (x == p[i] or x == \"1\"):\n",
     "            return False\n",
     "    return True"
    ]
@@ -486,7 +501,7 @@
     "def cond_conf(o, P_sharp):\n",
     "    \"\"\"Returns the (modified) conditionned expectation value of the cost function depending\n",
     "    on already chosen Paulis in P_sharp.\n",
-    "    \n",
+    "\n",
     "    Args:\n",
     "        o (list[str]): list of Pauli strings to be measured\n",
     "        P_sharp (list[str]): list of already chosen Paulis\n",
@@ -495,14 +510,20 @@
     "    eta = 0.9\n",
     "    nu = 1 - np.exp(-eta / 2)\n",
     "    L = len(o)\n",
-    "    m = len(P_sharp) - 1 # index of last chosen Pauli string\n",
-    "    k = len(P_sharp[-1]) - 1 # index of last chosen Pauli matrix in mth Pauli string\n",
+    "    m = len(P_sharp) - 1  # index of last chosen Pauli string\n",
+    "    k = (\n",
+    "        len(P_sharp[-1]) - 1\n",
+    "    )  # index of last chosen Pauli matrix in mth Pauli string\n",
     "    result = 0\n",
     "    for l in range(0, L):\n",
     "        v = 0\n",
-    "        for m_prime in range(0,m):\n",
+    "        for m_prime in range(0, m):\n",
     "            v += (eta / 2) * int(hits(P_sharp[m_prime], o[l]))\n",
-    "        v -= np.log(1 - (nu / 3**(weight(o[l], start=k+1))) * hits(P_sharp[m], o[l], end=k+1))\n",
+    "        v -= np.log(\n",
+    "            1\n",
+    "            - (nu / 3 ** (weight(o[l], start=k + 1)))\n",
+    "            * hits(P_sharp[m], o[l], end=k + 1)\n",
+    "        )\n",
     "        result += np.exp(-v)\n",
     "    return result"
    ]
@@ -523,7 +544,7 @@
     "def derandomization(M, o):\n",
     "    \"\"\"Derandomization algorithm returning best Pauli indices according to a greedy algorithm\n",
     "    that aims at minimizing the cost function above.\n",
-    "    \n",
+    "\n",
     "    Args:\n",
     "        M (int): number of measurements\n",
     "        n (int): number of qubits (size of Pauli strings)\n",
@@ -593,7 +614,7 @@
     "    assuming the state is already rotated in the needed eigenbases of all single-qubit Paulis.\n",
     "\n",
     "    NB : Faster than using qutip.measure due to not returning the eigenstates...\n",
-    "    \n",
+    "\n",
     "    Args:\n",
     "        x (str): input bitstring\n",
     "        sigma (str): input Pauli string to be measured on |x>\n",
@@ -615,9 +636,9 @@
    "source": [
     "def classical_shadow_derand(rho, measurements):\n",
     "    \"\"\"Returns the n-strings of ±1 corresponding to measurements in the input list on state rho.\n",
-    "    \n",
+    "\n",
     "    Args:\n",
-    "        rho (qutip.Qobj): input state as a density matrix \n",
+    "        rho (qutip.Qobj): input state as a density matrix\n",
     "        measurements (list[str]): derandomized measurement bases in which to measure state rho\n",
     "\n",
     "    Returns:\n",
@@ -629,7 +650,12 @@
     "    outcomes = []\n",
     "    for ns in range(shadow_size):\n",
     "        # multi-qubit change of basis\n",
-    "        unitmat = qutip.tensor([unitary_ensemble[_pauli_index(measurements[ns][i])]for i in range(num_qubits)])\n",
+    "        unitmat = qutip.tensor(\n",
+    "            [\n",
+    "                unitary_ensemble[_pauli_index(measurements[ns][i])]\n",
+    "                for i in range(num_qubits)\n",
+    "            ]\n",
+    "        )\n",
     "        x = measure_bitstring(unitmat.dag() * rho * unitmat)\n",
     "        outcomes.append(pauli_string_value(x, measurements[ns]))\n",
     "    # ±1 strings\n",
@@ -658,7 +684,8 @@
     "                break\n",
     "            if pauli != \"1\":\n",
     "                product *= single_measurement[i][1]\n",
-    "        if not_match: continue\n",
+    "        if not_match:\n",
+    "            continue\n",
     "\n",
     "        sum_product += product\n",
     "        cnt_match += 1\n",
@@ -719,7 +746,10 @@
    "outputs": [],
    "source": [
     "def pauli(positions=[], operators=[]):\n",
-    "    op_list = [operators[positions.index(j)] if j in positions else qutip.qeye(2) for j in range(num_qubits)]\n",
+    "    op_list = [\n",
+    "        operators[positions.index(j)] if j in positions else qutip.qeye(2)\n",
+    "        for j in range(num_qubits)\n",
+    "    ]\n",
     "    return qutip.tensor(op_list)"
    ]
   },
@@ -729,24 +759,30 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "coeff_fact = [0.81261,\n",
-    "              0.171201,\n",
-    "              0.2227965,\n",
-    "              0.16862325,\n",
-    "              0.174349,\n",
-    "              0.12054625,\n",
-    "              0.165868,\n",
-    "              0.04532175]\n",
+    "coeff_fact = [\n",
+    "    0.81261,\n",
+    "    0.171201,\n",
+    "    0.2227965,\n",
+    "    0.16862325,\n",
+    "    0.174349,\n",
+    "    0.12054625,\n",
+    "    0.165868,\n",
+    "    0.04532175,\n",
+    "]\n",
     "\n",
-    "paulis = [pauli(),\n",
-    "          pauli([0], [sz]) + pauli([1], [sz]),\n",
-    "          pauli([2], [sz]) + pauli([3], [sz]),\n",
-    "          pauli([1, 0], [sz, sz]),\n",
-    "          pauli([3, 2], [sz, sz]),\n",
-    "          pauli([2, 0], [sz, sz]) + pauli([3, 1], [sz, sz]),\n",
-    "          pauli([2, 1], [sz, sz]) + pauli([3, 0], [sz, sz]),\n",
-    "          pauli([3, 2, 1, 0], [sx, sx, sy, sy]) + pauli([3, 2, 1, 0], [sy, sy, sx, sx]),\n",
-    "          pauli([3, 2, 1, 0], [sx, sy, sy, sx]) + pauli([3, 2, 1, 0], [sy, sx, sx, sy])]"
+    "paulis = [\n",
+    "    pauli(),\n",
+    "    pauli([0], [sz]) + pauli([1], [sz]),\n",
+    "    pauli([2], [sz]) + pauli([3], [sz]),\n",
+    "    pauli([1, 0], [sz, sz]),\n",
+    "    pauli([3, 2], [sz, sz]),\n",
+    "    pauli([2, 0], [sz, sz]) + pauli([3, 1], [sz, sz]),\n",
+    "    pauli([2, 1], [sz, sz]) + pauli([3, 0], [sz, sz]),\n",
+    "    pauli([3, 2, 1, 0], [sx, sx, sy, sy])\n",
+    "    + pauli([3, 2, 1, 0], [sy, sy, sx, sx]),\n",
+    "    pauli([3, 2, 1, 0], [sx, sy, sy, sx])\n",
+    "    + pauli([3, 2, 1, 0], [sy, sx, sx, sy]),\n",
+    "]"
    ]
   },
   {
@@ -772,7 +808,14 @@
    "source": [
     "# H2 Molecule : 4 qubits in Jordan-Wigner mapping of the Hamiltonian\n",
     "a = 10\n",
-    "reg = Register.from_coordinates([[0, 0], [a, 0], [0.5*a, a*np.sqrt(3)/2], [0.5*a, -a*np.sqrt(3)/2]])\n",
+    "reg = Register.from_coordinates(\n",
+    "    [\n",
+    "        [0, 0],\n",
+    "        [a, 0],\n",
+    "        [0.5 * a, a * np.sqrt(3) / 2],\n",
+    "        [0.5 * a, -a * np.sqrt(3) / 2],\n",
+    "    ]\n",
+    ")\n",
     "reg.draw()"
    ]
   },
@@ -798,17 +841,20 @@
    ],
    "source": [
     "def cost_hamiltonian_JW():\n",
-    "    H = - coeff_fact[0] * paulis[0] \\\n",
-    "        + coeff_fact[1] * paulis[1] \\\n",
-    "        - coeff_fact[2] * paulis[2] \\\n",
-    "        + coeff_fact[3] * paulis[3] \\\n",
-    "        + coeff_fact[4] * paulis[4] \\\n",
-    "        + coeff_fact[5] * paulis[5] \\\n",
-    "        + coeff_fact[6] * paulis[6] \\\n",
-    "        - coeff_fact[7] * paulis[7] \\\n",
+    "    H = (\n",
+    "        -coeff_fact[0] * paulis[0]\n",
+    "        + coeff_fact[1] * paulis[1]\n",
+    "        - coeff_fact[2] * paulis[2]\n",
+    "        + coeff_fact[3] * paulis[3]\n",
+    "        + coeff_fact[4] * paulis[4]\n",
+    "        + coeff_fact[5] * paulis[5]\n",
+    "        + coeff_fact[6] * paulis[6]\n",
+    "        - coeff_fact[7] * paulis[7]\n",
     "        + coeff_fact[7] * paulis[8]\n",
+    "    )\n",
     "    return H\n",
     "\n",
+    "\n",
     "global H\n",
     "H = cost_hamiltonian_JW()\n",
     "exact_energy, ground_state = cost_hamiltonian_JW().groundstate()\n",
@@ -851,25 +897,31 @@
     "        in_state (qubit.Qobj): initial state.\n",
     "    \"\"\"\n",
     "    seq = Sequence(r, Chadoq2)\n",
-    "    seq.declare_channel('ch0','rydberg_global')\n",
-    "    middle = len(param)//2\n",
-    "    \n",
+    "    seq.declare_channel(\"ch0\", \"rydberg_global\")\n",
+    "    middle = len(param) // 2\n",
+    "\n",
     "    for tau, t in zip(param[middle:], param[:middle]):\n",
-    "        pulse_1 = Pulse.ConstantPulse(tau, 1., 0, 0) \n",
-    "        pulse_2 = Pulse.ConstantPulse(t, 1., 1., 0) \n",
-    "        seq.add(pulse_1, 'ch0')\n",
-    "        seq.add(pulse_2, 'ch0')\n",
-    "    \n",
-    "    seq.measure('ground-rydberg')\n",
-    "    simul = Simulation(seq, sampling_rate=.05)\n",
+    "        pulse_1 = Pulse.ConstantPulse(tau, 1.0, 0, 0)\n",
+    "        pulse_2 = Pulse.ConstantPulse(t, 1.0, 1.0, 0)\n",
+    "        seq.add(pulse_1, \"ch0\")\n",
+    "        seq.add(pulse_2, \"ch0\")\n",
+    "\n",
+    "    seq.measure(\"ground-rydberg\")\n",
+    "    simul = Simulation(seq, sampling_rate=0.05)\n",
     "    simul.initial_state = in_state\n",
     "    results = simul.run()\n",
     "    return results.expect([H])[-1][-1]\n",
     "\n",
+    "\n",
     "def loop_JW(param, in_state):\n",
-    "    res = minimize(quantum_loop, param, method='Nelder-Mead', args=in_state,\n",
-    "                   options={'return_all':True, 'maxiter':200, 'adaptive':True})\n",
-    "    return(res)"
+    "    res = minimize(\n",
+    "        quantum_loop,\n",
+    "        param,\n",
+    "        method=\"Nelder-Mead\",\n",
+    "        args=in_state,\n",
+    "        options={\"return_all\": True, \"maxiter\": 200, \"adaptive\": True},\n",
+    "    )\n",
+    "    return res"
    ]
   },
   {
@@ -887,7 +939,7 @@
    "source": [
     "# Setup for VQS\n",
     "layers = 5\n",
-    "param = [2000]*layers + [4000]*layers"
+    "param = [2000] * layers + [4000] * layers"
    ]
   },
   {
@@ -912,6 +964,7 @@
    ],
    "source": [
     "import warnings\n",
+    "\n",
     "# Ignore the warnings\n",
     "warnings.filterwarnings(\"ignore\", category=UserWarning)\n",
     "\n",
@@ -957,8 +1010,10 @@
     }
    ],
    "source": [
-    "plt.plot([quantum_loop(pars, gggg) for pars in loop_ising_results.allvecs], 'k')\n",
-    "plt.axhline(exact_energy, color='red')"
+    "plt.plot(\n",
+    "    [quantum_loop(pars, gggg) for pars in loop_ising_results.allvecs], \"k\"\n",
+    ")\n",
+    "plt.axhline(exact_energy, color=\"red\")"
    ]
   },
   {
@@ -997,15 +1052,17 @@
    "outputs": [],
    "source": [
     "def exp_value_JW(exp_values):\n",
-    "    return (- coeff_fact[0] * exp_values[0] \\\n",
-    "            + coeff_fact[1] * exp_values[1] \\\n",
-    "            - coeff_fact[2] * exp_values[2] \\\n",
-    "            + coeff_fact[3] * exp_values[3] \\\n",
-    "            + coeff_fact[4] * exp_values[4] \\\n",
-    "            + coeff_fact[5] * exp_values[5] \\\n",
-    "            + coeff_fact[6] * exp_values[6] \\\n",
-    "            - coeff_fact[7] * exp_values[7] \\\n",
-    "            + coeff_fact[7] * exp_values[8])"
+    "    return (\n",
+    "        -coeff_fact[0] * exp_values[0]\n",
+    "        + coeff_fact[1] * exp_values[1]\n",
+    "        - coeff_fact[2] * exp_values[2]\n",
+    "        + coeff_fact[3] * exp_values[3]\n",
+    "        + coeff_fact[4] * exp_values[4]\n",
+    "        + coeff_fact[5] * exp_values[5]\n",
+    "        + coeff_fact[6] * exp_values[6]\n",
+    "        - coeff_fact[7] * exp_values[7]\n",
+    "        + coeff_fact[7] * exp_values[8]\n",
+    "    )"
    ]
   },
   {
@@ -1021,19 +1078,19 @@
     "        in_state (qubit.Qobj): initial state.\n",
     "    \"\"\"\n",
     "    seq = Sequence(r, Chadoq2)\n",
-    "    seq.declare_channel('ch0','rydberg_global')\n",
-    "    middle = len(param)//2\n",
-    "    \n",
+    "    seq.declare_channel(\"ch0\", \"rydberg_global\")\n",
+    "    middle = len(param) // 2\n",
+    "\n",
     "    for tau, t in zip(param[middle:], param[:middle]):\n",
-    "        pulse_1 = Pulse.ConstantPulse(tau, 1., 0, 0) \n",
-    "        pulse_2 = Pulse.ConstantPulse(t, 1., 1., 0)\n",
-    "        seq.add(pulse_1, 'ch0')\n",
-    "        seq.add(pulse_2, 'ch0')\n",
-    "    \n",
-    "    seq.measure('ground-rydberg')\n",
-    "    simul = Simulation(seq, sampling_rate=.01)\n",
+    "        pulse_1 = Pulse.ConstantPulse(tau, 1.0, 0, 0)\n",
+    "        pulse_2 = Pulse.ConstantPulse(t, 1.0, 1.0, 0)\n",
+    "        seq.add(pulse_1, \"ch0\")\n",
+    "        seq.add(pulse_2, \"ch0\")\n",
+    "\n",
+    "    seq.measure(\"ground-rydberg\")\n",
+    "    simul = Simulation(seq, sampling_rate=0.01)\n",
     "    simul.initial_state = in_state\n",
-    "    \n",
+    "\n",
     "    # Classical shadow estimation\n",
     "    # Theoretical shadow size and number of clusters :\n",
     "    # shadow_size, K = compute_shadow_size(0.1, 0.5, paulis)\n",
@@ -1041,14 +1098,23 @@
     "    K = 4\n",
     "    rho = simul.run().get_final_state().proj()\n",
     "    outcomes, unitary_ids = calculate_classical_shadow(rho, shadow_size)\n",
-    "    snapshots = [snapshot_state(outcomes[ns], unitary_ids[ns]) for ns in range(shadow_size)]\n",
+    "    snapshots = [\n",
+    "        snapshot_state(outcomes[ns], unitary_ids[ns])\n",
+    "        for ns in range(shadow_size)\n",
+    "    ]\n",
     "    meds = [_median_of_means(obs, snapshots, K) for obs in paulis]\n",
     "    return exp_value_JW(meds)\n",
     "\n",
+    "\n",
     "def loop_JW_shadows(param, in_state, shadow_size=20):\n",
-    "    res = minimize(quantum_loop_shadows, param, method='Nelder-Mead', args=(in_state, shadow_size),\n",
-    "                   options={'return_all':True, 'maxiter':100, 'adaptive':True})\n",
-    "    return(res)"
+    "    res = minimize(\n",
+    "        quantum_loop_shadows,\n",
+    "        param,\n",
+    "        method=\"Nelder-Mead\",\n",
+    "        args=(in_state, shadow_size),\n",
+    "        options={\"return_all\": True, \"maxiter\": 100, \"adaptive\": True},\n",
+    "    )\n",
+    "    return res"
    ]
   },
   {
@@ -1057,10 +1123,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "shadow_sizes = [10,20,40,60,80,100]\n",
+    "shadow_sizes = [10, 20, 40, 60, 80, 100]\n",
     "energies = []\n",
     "for shadow_size in shadow_sizes:\n",
-    "    energies.append(abs(loop_JW_shadows(param, gggg, shadow_size=shadow_size).fun - exact_energy))"
+    "    energies.append(\n",
+    "        abs(\n",
+    "            loop_JW_shadows(param, gggg, shadow_size=shadow_size).fun\n",
+    "            - exact_energy\n",
+    "        )\n",
+    "    )"
    ]
   },
   {
@@ -1092,10 +1163,10 @@
     }
    ],
    "source": [
-    "plt.figure(figsize=(8,5))\n",
+    "plt.figure(figsize=(8, 5))\n",
     "plt.xlabel(\"Shadow size\", fontsize=15)\n",
     "plt.ylabel(r\"$|\\frac{E - E_{ground}}{E_{ground}}|$\", fontsize=20)\n",
-    "plt.plot(shadow_sizes, [-e/exact_energy for e in energies])"
+    "plt.plot(shadow_sizes, [-e / exact_energy for e in energies])"
    ]
   },
   {
@@ -1125,24 +1196,41 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "coeff_non_fact = [-0.81261,\n",
-    "                  0.171201,\n",
-    "                  0.171201,\n",
-    "                  -0.2227965,\n",
-    "                  -0.2227965,\n",
-    "                  0.16862325,\n",
-    "                  0.174349,\n",
-    "                  0.12054625,\n",
-    "                  0.12054625,\n",
-    "                  0.165868,\n",
-    "                  0.165868,\n",
-    "                  -0.04532175,\n",
-    "                  -0.04532175,\n",
-    "                  0.04532175,\n",
-    "                  0.04532175]\n",
+    "coeff_non_fact = [\n",
+    "    -0.81261,\n",
+    "    0.171201,\n",
+    "    0.171201,\n",
+    "    -0.2227965,\n",
+    "    -0.2227965,\n",
+    "    0.16862325,\n",
+    "    0.174349,\n",
+    "    0.12054625,\n",
+    "    0.12054625,\n",
+    "    0.165868,\n",
+    "    0.165868,\n",
+    "    -0.04532175,\n",
+    "    -0.04532175,\n",
+    "    0.04532175,\n",
+    "    0.04532175,\n",
+    "]\n",
     "\n",
-    "paulis_str = [\"1111\", \"Z111\", \"1Z11\", \"11Z1\", \"111Z\", \"ZZ11\", \"11ZZ\", \"Z1Z1\", \"1Z1Z\", \"1ZZ1\",\n",
-    "              \"Z11Z\", \"YYXX\", \"XXYY\", \"XYYX\", \"YXXY\"]"
+    "paulis_str = [\n",
+    "    \"1111\",\n",
+    "    \"Z111\",\n",
+    "    \"1Z11\",\n",
+    "    \"11Z1\",\n",
+    "    \"111Z\",\n",
+    "    \"ZZ11\",\n",
+    "    \"11ZZ\",\n",
+    "    \"Z1Z1\",\n",
+    "    \"1Z1Z\",\n",
+    "    \"1ZZ1\",\n",
+    "    \"Z11Z\",\n",
+    "    \"YYXX\",\n",
+    "    \"XXYY\",\n",
+    "    \"XYYX\",\n",
+    "    \"YXXY\",\n",
+    "]"
    ]
   },
   {
@@ -1152,7 +1240,12 @@
    "outputs": [],
    "source": [
     "def exp_value_JW_non_fact(outcomes):\n",
-    "    return sum([c*exp_value(sigma, outcomes) for c, sigma in zip(coeff_non_fact, paulis_str)])"
+    "    return sum(\n",
+    "        [\n",
+    "            c * exp_value(sigma, outcomes)\n",
+    "            for c, sigma in zip(coeff_non_fact, paulis_str)\n",
+    "        ]\n",
+    "    )"
    ]
   },
   {
@@ -1177,9 +1270,11 @@
    ],
    "source": [
     "measurements = derandomization(60, paulis_str)\n",
-    "print(f\"ZZZZ measurements : {measurements.count('ZZZZ')}, XXYY measurements : {measurements.count('XXYY')}, \" +\n",
-    "      f\"YXXY measurements : {measurements.count('YXXY')}, XYYX measurements : {measurements.count('XYYX')}, \" +\n",
-    "      f\"YYXX measurements : {measurements.count('YYXX')} : total = 60 measurements\")"
+    "print(\n",
+    "    f\"ZZZZ measurements : {measurements.count('ZZZZ')}, XXYY measurements : {measurements.count('XXYY')}, \"\n",
+    "    + f\"YXXY measurements : {measurements.count('YXXY')}, XYYX measurements : {measurements.count('XYYX')}, \"\n",
+    "    + f\"YYXX measurements : {measurements.count('YYXX')} : total = 60 measurements\"\n",
+    ")"
    ]
   },
   {
@@ -1202,28 +1297,34 @@
     "        in_state (qubit.Qobj): initial state.\n",
     "    \"\"\"\n",
     "    seq = Sequence(r, Chadoq2)\n",
-    "    seq.declare_channel('ch0','rydberg_global')\n",
-    "    middle = len(param)//2\n",
-    "    \n",
+    "    seq.declare_channel(\"ch0\", \"rydberg_global\")\n",
+    "    middle = len(param) // 2\n",
+    "\n",
     "    for tau, t in zip(param[middle:], param[:middle]):\n",
-    "        pulse_1 = Pulse.ConstantPulse(tau, 1., 0, 0) \n",
-    "        pulse_2 = Pulse.ConstantPulse(t, 1., 1., 0)\n",
-    "        seq.add(pulse_1, 'ch0')\n",
-    "        seq.add(pulse_2, 'ch0')\n",
-    "    \n",
-    "    seq.measure('ground-rydberg')\n",
-    "    simul = Simulation(seq, sampling_rate=.05)\n",
+    "        pulse_1 = Pulse.ConstantPulse(tau, 1.0, 0, 0)\n",
+    "        pulse_2 = Pulse.ConstantPulse(t, 1.0, 1.0, 0)\n",
+    "        seq.add(pulse_1, \"ch0\")\n",
+    "        seq.add(pulse_2, \"ch0\")\n",
+    "\n",
+    "    seq.measure(\"ground-rydberg\")\n",
+    "    simul = Simulation(seq, sampling_rate=0.05)\n",
     "    simul.initial_state = in_state\n",
-    "    \n",
+    "\n",
     "    # Classical shadow estimation\n",
     "    rho = simul.run().get_final_state().proj()\n",
     "    outcomes = classical_shadow_derand(rho, measurements)\n",
     "    return exp_value_JW_non_fact(outcomes)\n",
     "\n",
+    "\n",
     "def loop_JW_derand(param, in_state):\n",
-    "    res = minimize(quantum_loop_derand, param, method='Nelder-Mead', args=in_state,\n",
-    "                   options={'return_all':True, 'maxiter':150, 'adaptive':True})\n",
-    "    return(res)"
+    "    res = minimize(\n",
+    "        quantum_loop_derand,\n",
+    "        param,\n",
+    "        method=\"Nelder-Mead\",\n",
+    "        args=in_state,\n",
+    "        options={\"return_all\": True, \"maxiter\": 150, \"adaptive\": True},\n",
+    "    )\n",
+    "    return res"
    ]
   },
   {
@@ -1232,11 +1333,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "measurement_sizes = [20,30,40,60,80,100]\n",
+    "measurement_sizes = [20, 30, 40, 60, 80, 100]\n",
     "energies_derand = []\n",
     "for meas_size in measurement_sizes:\n",
-    "    measurements=derandomization(meas_size, paulis_str)\n",
-    "    energies_derand.append(abs(loop_JW_derand(param, gggg).fun - exact_energy) / abs(exact_energy))"
+    "    measurements = derandomization(meas_size, paulis_str)\n",
+    "    energies_derand.append(\n",
+    "        abs(loop_JW_derand(param, gggg).fun - exact_energy) / abs(exact_energy)\n",
+    "    )"
    ]
   },
   {
@@ -1268,7 +1371,7 @@
     }
    ],
    "source": [
-    "plt.figure(figsize=(8,5))\n",
+    "plt.figure(figsize=(8, 5))\n",
     "plt.xlabel(\"Measurement size\", fontsize=15)\n",
     "plt.ylabel(r\"$|\\frac{E - E_{ground}}{E_{ground}}|$\", fontsize=20)\n",
     "plt.plot(measurement_sizes, energies_derand)"
diff --git a/tutorials/quantum_simulation/Spin chain of 3 atoms in XY mode.ipynb b/tutorials/quantum_simulation/Spin chain of 3 atoms in XY mode.ipynb
index e805f0b50..4c23f0381 100644
--- a/tutorials/quantum_simulation/Spin chain of 3 atoms in XY mode.ipynb
+++ b/tutorials/quantum_simulation/Spin chain of 3 atoms in XY mode.ipynb
@@ -74,7 +74,7 @@
     "\n",
     "reg = Register(qubits)\n",
     "seq = Sequence(reg, MockDevice)\n",
-    "seq.declare_channel('MW', 'mw_global')"
+    "seq.declare_channel(\"MW\", \"mw_global\")"
    ]
   },
   {
@@ -90,9 +90,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "simple_pulse = Pulse.ConstantPulse(4000, 2*np.pi*4.6, 0, 0)\n",
-    "seq.add(simple_pulse, 'MW')\n",
-    "seq.measure(basis='XY')\n",
+    "simple_pulse = Pulse.ConstantPulse(4000, 2 * np.pi * 4.6, 0, 0)\n",
+    "seq.add(simple_pulse, \"MW\")\n",
+    "seq.measure(basis=\"XY\")\n",
     "\n",
     "sim = Simulation(seq)\n",
     "\n",
@@ -117,11 +117,12 @@
     "    prod[j] = (qutip.sigmaz() + qutip.qeye(2)) / 2\n",
     "    return qutip.tensor(prod)\n",
     "\n",
+    "\n",
     "magn = magnetization(0, 1)\n",
     "plt.figure(figsize=[16, 6])\n",
     "results.plot(magn)\n",
-    "plt.xlabel('Pulse duration (ns)', fontsize='x-large')\n",
-    "plt.ylabel('Excitation of the atom', fontsize='x-large')\n",
+    "plt.xlabel(\"Pulse duration (ns)\", fontsize=\"x-large\")\n",
+    "plt.ylabel(\"Excitation of the atom\", fontsize=\"x-large\")\n",
     "plt.show()"
    ]
   },
@@ -145,24 +146,24 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "coords = np.array([[-8., 0], [0, 0], [8., 0]])\n",
+    "coords = np.array([[-8.0, 0], [0, 0], [8.0, 0]])\n",
     "qubits = dict(enumerate(coords))\n",
     "\n",
     "reg = Register(qubits)\n",
     "seq = Sequence(reg, MockDevice)\n",
-    "seq.declare_channel('ch0', 'mw_global')\n",
+    "seq.declare_channel(\"ch0\", \"mw_global\")\n",
     "reg.draw()\n",
     "\n",
     "# State preparation using SLM mask\n",
     "masked_qubits = [1, 2]\n",
     "seq.config_slm_mask(masked_qubits)\n",
     "masked_pulse = Pulse.ConstantDetuning(BlackmanWaveform(200, np.pi), 0, 0)\n",
-    "seq.add(masked_pulse, 'ch0')\n",
+    "seq.add(masked_pulse, \"ch0\")\n",
     "\n",
     "# Simulation pulse\n",
     "simple_pulse = Pulse.ConstantPulse(7000, 0, 0, 0)\n",
-    "seq.add(simple_pulse, 'ch0')\n",
-    "seq.measure(basis='XY')\n",
+    "seq.add(simple_pulse, \"ch0\")\n",
+    "seq.measure(basis=\"XY\")\n",
     "\n",
     "sim = Simulation(seq, sampling_rate=1)\n",
     "results = sim.run(nsteps=5000)"
@@ -181,17 +182,17 @@
     "plt.figure(figsize=[16, 18])\n",
     "plt.subplot(311)\n",
     "plt.plot(expectations[0])\n",
-    "plt.ylabel('Excitation of atom 0', fontsize='x-large')\n",
-    "plt.xlabel('Time (ns)', fontsize='x-large')\n",
+    "plt.ylabel(\"Excitation of atom 0\", fontsize=\"x-large\")\n",
+    "plt.xlabel(\"Time (ns)\", fontsize=\"x-large\")\n",
     "plt.subplot(312)\n",
     "plt.plot(expectations[1])\n",
-    "plt.ylabel('Excitation of atom 1', fontsize='x-large')\n",
-    "plt.xlabel('Time (ns)', fontsize='x-large')\n",
+    "plt.ylabel(\"Excitation of atom 1\", fontsize=\"x-large\")\n",
+    "plt.xlabel(\"Time (ns)\", fontsize=\"x-large\")\n",
     "plt.ylim([0, 1])\n",
     "plt.subplot(313)\n",
     "plt.plot(expectations[2])\n",
-    "plt.ylabel('Excitation of atom 2', fontsize='x-large')\n",
-    "plt.xlabel('Time (ns)', fontsize='x-large')\n",
+    "plt.ylabel(\"Excitation of atom 2\", fontsize=\"x-large\")\n",
+    "plt.xlabel(\"Time (ns)\", fontsize=\"x-large\")\n",
     "plt.show()"
    ]
   },
@@ -234,13 +235,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "coords = np.array([[-1., 0], [0, 0], [np.sqrt(2/3), np.sqrt(1/3)]]) * 8.\n",
+    "coords = np.array([[-1.0, 0], [0, 0], [np.sqrt(2 / 3), np.sqrt(1 / 3)]]) * 8.0\n",
     "qubits = dict(enumerate(coords))\n",
     "\n",
     "reg = Register(qubits)\n",
     "seq = Sequence(reg, MockDevice)\n",
-    "seq.declare_channel('ch0', 'mw_global')\n",
-    "seq.set_magnetic_field(0., 1., 0)\n",
+    "seq.declare_channel(\"ch0\", \"mw_global\")\n",
+    "seq.set_magnetic_field(0.0, 1.0, 0)\n",
     "reg.draw()"
    ]
   },
@@ -261,12 +262,12 @@
     "masked_qubits = [1, 2]\n",
     "seq.config_slm_mask(masked_qubits)\n",
     "masked_pulse = Pulse.ConstantDetuning(BlackmanWaveform(200, np.pi), 0, 0)\n",
-    "seq.add(masked_pulse, 'ch0')\n",
+    "seq.add(masked_pulse, \"ch0\")\n",
     "\n",
     "# Simulation pulse\n",
     "simple_pulse = Pulse.ConstantPulse(7000, 0, 0, 0)\n",
-    "seq.add(simple_pulse, 'ch0')\n",
-    "seq.measure(basis='XY')\n",
+    "seq.add(simple_pulse, \"ch0\")\n",
+    "seq.measure(basis=\"XY\")\n",
     "\n",
     "sim = Simulation(seq, sampling_rate=1)\n",
     "results = sim.run(progress_bar=True, nsteps=5000)"
@@ -285,18 +286,18 @@
     "plt.figure(figsize=[16, 18])\n",
     "plt.subplot(311)\n",
     "plt.plot(expectations[0])\n",
-    "plt.ylabel('Excitation of atom 0', fontsize='x-large')\n",
-    "plt.xlabel('Time ($\\mu$s)', fontsize='x-large')\n",
+    "plt.ylabel(\"Excitation of atom 0\", fontsize=\"x-large\")\n",
+    "plt.xlabel(\"Time ($\\mu$s)\", fontsize=\"x-large\")\n",
     "plt.ylim([0, 1])\n",
     "plt.subplot(312)\n",
     "plt.plot(expectations[1])\n",
-    "plt.ylabel('Excitation of atom 1', fontsize='x-large')\n",
-    "plt.xlabel('Time ($\\mu$s)', fontsize='x-large')\n",
+    "plt.ylabel(\"Excitation of atom 1\", fontsize=\"x-large\")\n",
+    "plt.xlabel(\"Time ($\\mu$s)\", fontsize=\"x-large\")\n",
     "plt.ylim([0, 1])\n",
     "plt.subplot(313)\n",
     "plt.plot(expectations[2])\n",
-    "plt.ylabel('Excitation of atom 2', fontsize='x-large')\n",
-    "plt.xlabel('Time ($\\mu$s)', fontsize='x-large')\n",
+    "plt.ylabel(\"Excitation of atom 2\", fontsize=\"x-large\")\n",
+    "plt.xlabel(\"Time ($\\mu$s)\", fontsize=\"x-large\")\n",
     "plt.ylim([0, 1])"
    ]
   },
diff --git a/tutorials/simulating_sequences.ipynb b/tutorials/simulating_sequences.ipynb
index 367c62590..03fa994e7 100644
--- a/tutorials/simulating_sequences.ipynb
+++ b/tutorials/simulating_sequences.ipynb
@@ -39,7 +39,7 @@
     "# Setup\n",
     "L = 14\n",
     "\n",
-    "Omega_max = 2.3 * 2*np.pi \n",
+    "Omega_max = 2.3 * 2 * np.pi\n",
     "U = Omega_max / 2.3\n",
     "\n",
     "delta_0 = -3 * U\n",
@@ -47,15 +47,24 @@
     "\n",
     "t_rise = 2000\n",
     "t_fall = 2000\n",
-    "t_sweep = (delta_f - delta_0)/(2 * np.pi * 10) * 5000\n",
+    "t_sweep = (delta_f - delta_0) / (2 * np.pi * 10) * 5000\n",
     "\n",
     "# Define a ring of atoms distanced by a blockade radius distance:\n",
     "R_interatomic = MockDevice.rydberg_blockade_radius(U)\n",
-    "coords = R_interatomic/(2*np.tan(np.pi/L)) * np.array([(np.cos(theta*2*np.pi/L), np.sin(theta*2*np.pi/L)) for theta in range(L)])\n",
-    "                                  \n",
-    "reg = Register.from_coordinates(coords, prefix='atom')\n",
+    "coords = (\n",
+    "    R_interatomic\n",
+    "    / (2 * np.tan(np.pi / L))\n",
+    "    * np.array(\n",
+    "        [\n",
+    "            (np.cos(theta * 2 * np.pi / L), np.sin(theta * 2 * np.pi / L))\n",
+    "            for theta in range(L)\n",
+    "        ]\n",
+    "    )\n",
+    ")\n",
     "\n",
-    "reg.draw(blockade_radius=R_interatomic, draw_half_radius=True, draw_graph = True)"
+    "reg = Register.from_coordinates(coords, prefix=\"atom\")\n",
+    "\n",
+    "reg.draw(blockade_radius=R_interatomic, draw_half_radius=True, draw_graph=True)"
    ]
   },
   {
@@ -78,16 +87,22 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "rise = Pulse.ConstantDetuning(RampWaveform(t_rise, 0., Omega_max), delta_0, 0.)\n",
-    "sweep = Pulse.ConstantAmplitude(Omega_max, RampWaveform(t_sweep, delta_0, delta_f), 0.)\n",
-    "fall = Pulse.ConstantDetuning(RampWaveform(t_fall, Omega_max, 0.), delta_f, 0.)\n",
+    "rise = Pulse.ConstantDetuning(\n",
+    "    RampWaveform(t_rise, 0.0, Omega_max), delta_0, 0.0\n",
+    ")\n",
+    "sweep = Pulse.ConstantAmplitude(\n",
+    "    Omega_max, RampWaveform(t_sweep, delta_0, delta_f), 0.0\n",
+    ")\n",
+    "fall = Pulse.ConstantDetuning(\n",
+    "    RampWaveform(t_fall, Omega_max, 0.0), delta_f, 0.0\n",
+    ")\n",
     "\n",
     "seq = Sequence(reg, MockDevice)\n",
-    "seq.declare_channel('ising', 'rydberg_global')\n",
+    "seq.declare_channel(\"ising\", \"rydberg_global\")\n",
     "\n",
-    "seq.add(rise, 'ising')\n",
-    "seq.add(sweep, 'ising')\n",
-    "seq.add(fall, 'ising')\n",
+    "seq.add(rise, \"ising\")\n",
+    "seq.add(sweep, \"ising\")\n",
+    "seq.add(fall, \"ising\")\n",
     "\n",
     "seq.draw()"
    ]
@@ -158,7 +173,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "results.states[23] # Given as a `qutip.Qobj` object"
+    "results.states[23]  # Given as a `qutip.Qobj` object"
    ]
   },
   {
@@ -176,9 +191,9 @@
    "source": [
     "counts = results.sample_final_state(N_samples=1000)\n",
     "\n",
-    "large_counts = {k:v for k,v in counts.items() if v > 5}\n",
+    "large_counts = {k: v for k, v in counts.items() if v > 5}\n",
     "\n",
-    "plt.figure(figsize=(15,4))\n",
+    "plt.figure(figsize=(15, 4))\n",
     "plt.xticks(rotation=90, fontsize=14)\n",
     "plt.title(\"Most frequent observations\")\n",
     "plt.bar(large_counts.keys(), large_counts.values())"
@@ -209,6 +224,7 @@
     "    prod[j] = qutip.sigmaz()\n",
     "    return qutip.tensor(prod)\n",
     "\n",
+    "\n",
     "magn_list = [magnetization(j, L) for j in range(L)]"
    ]
   },


Test patch:
diff --git a/pulser/tests/test_channels.py b/pulser/tests/test_channels.py
index 4c25202df..532be943f 100644
--- a/pulser/tests/test_channels.py
+++ b/pulser/tests/test_channels.py
@@ -17,7 +17,7 @@
 
 import pulser
 from pulser.channels import Raman, Rydberg
-from pulser.waveforms import ConstantWaveform, BlackmanWaveform
+from pulser.waveforms import BlackmanWaveform, ConstantWaveform
 
 
 def test_device_channels():
diff --git a/pulser/tests/test_json.py b/pulser/tests/test_json.py
index 5eddc1dc7..482b87c00 100644
--- a/pulser/tests/test_json.py
+++ b/pulser/tests/test_json.py
@@ -17,9 +17,9 @@
 import numpy as np
 import pytest
 
-from pulser import Sequence, Register, Register3D
+from pulser import Register, Register3D, Sequence
 from pulser.devices import Chadoq2, MockDevice
-from pulser.json.coders import PulserEncoder, PulserDecoder
+from pulser.json.coders import PulserDecoder, PulserEncoder
 from pulser.parametrized.decorators import parametrize
 from pulser.waveforms import BlackmanWaveform
 
diff --git a/pulser/tests/test_parametrized.py b/pulser/tests/test_parametrized.py
index c0bee6d0f..233869a6f 100644
--- a/pulser/tests/test_parametrized.py
+++ b/pulser/tests/test_parametrized.py
@@ -21,7 +21,6 @@
 from pulser.parametrized import Variable
 from pulser.waveforms import BlackmanWaveform, CompositeWaveform
 
-
 a = Variable("a", float)
 b = Variable("b", int, size=2)
 b._assign([-1.5, 1.5])
diff --git a/pulser/tests/test_paramseq.py b/pulser/tests/test_paramseq.py
index f210cbedc..04bbf527b 100644
--- a/pulser/tests/test_paramseq.py
+++ b/pulser/tests/test_paramseq.py
@@ -17,7 +17,7 @@
 import numpy as np
 import pytest
 
-from pulser import Sequence, Register, Pulse
+from pulser import Pulse, Register, Sequence
 from pulser.devices import Chadoq2, MockDevice
 from pulser.parametrized import Variable
 from pulser.waveforms import BlackmanWaveform
diff --git a/pulser/tests/test_pulse.py b/pulser/tests/test_pulse.py
index b63793d77..595ce469e 100644
--- a/pulser/tests/test_pulse.py
+++ b/pulser/tests/test_pulse.py
@@ -18,7 +18,7 @@
 import pytest
 
 from pulser import Pulse
-from pulser.waveforms import ConstantWaveform, BlackmanWaveform, RampWaveform
+from pulser.waveforms import BlackmanWaveform, ConstantWaveform, RampWaveform
 
 cwf = ConstantWaveform(100, -10)
 bwf = BlackmanWaveform(200, 3)
diff --git a/pulser/tests/test_sequence.py b/pulser/tests/test_sequence.py
index 1899b1ef3..0fac5dfbc 100644
--- a/pulser/tests/test_sequence.py
+++ b/pulser/tests/test_sequence.py
@@ -18,16 +18,16 @@
 import pytest
 
 import pulser
-from pulser import Sequence, Pulse, Register, Register3D
-from pulser.channels import Rydberg, Raman
+from pulser import Pulse, Register, Register3D, Sequence
+from pulser.channels import Raman, Rydberg
 from pulser.devices import Chadoq2, MockDevice
 from pulser.devices._device_datacls import Device
 from pulser.sequence import _TimeSlot
 from pulser.waveforms import (
     BlackmanWaveform,
     CompositeWaveform,
-    RampWaveform,
     InterpolatedWaveform,
+    RampWaveform,
 )
 
 reg = Register.triangular_lattice(4, 7, spacing=5, prefix="q")
diff --git a/pulser/tests/test_simresults.py b/pulser/tests/test_simresults.py
index f8219c884..ec0aecadf 100644
--- a/pulser/tests/test_simresults.py
+++ b/pulser/tests/test_simresults.py
@@ -11,20 +11,19 @@
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
-from copy import deepcopy
-
 from collections import Counter
+from copy import deepcopy
 
 import numpy as np
 import pytest
 import qutip
 from qutip.piqs import isdiagonal
 
-from pulser import Sequence, Pulse, Register
+from pulser import Pulse, Register, Sequence
 from pulser.devices import Chadoq2, MockDevice
-from pulser.waveforms import BlackmanWaveform
-from pulser.simulation import Simulation, SimConfig
+from pulser.simulation import SimConfig, Simulation
 from pulser.simulation.simresults import CoherentResults, NoisyResults
+from pulser.waveforms import BlackmanWaveform
 
 np.random.seed(123)
 q_dict = {
diff --git a/pulser/tests/test_simulation.py b/pulser/tests/test_simulation.py
index a3707909d..9db3e9f08 100644
--- a/pulser/tests/test_simulation.py
+++ b/pulser/tests/test_simulation.py
@@ -12,18 +12,17 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-from unittest.mock import patch
-
 from collections import Counter
+from unittest.mock import patch
 
 import numpy as np
 import pytest
 import qutip
 
-from pulser import Sequence, Pulse, Register
+from pulser import Pulse, Register, Sequence
 from pulser.devices import Chadoq2, MockDevice
-from pulser.waveforms import BlackmanWaveform, RampWaveform, ConstantWaveform
 from pulser.simulation import SimConfig, Simulation
+from pulser.waveforms import BlackmanWaveform, ConstantWaveform, RampWaveform
 
 q_dict = {
     "control1": np.array([-4.0, 0.0]),
diff --git a/pulser/tests/test_waveforms.py b/pulser/tests/test_waveforms.py
index fa10d51b3..3c140ddb3 100644
--- a/pulser/tests/test_waveforms.py
+++ b/pulser/tests/test_waveforms.py
@@ -18,19 +18,19 @@
 
 import numpy as np
 import pytest
-from scipy.interpolate import interp1d, PchipInterpolator
+from scipy.interpolate import PchipInterpolator, interp1d
 
-from pulser.json.coders import PulserEncoder, PulserDecoder
 from pulser.channels import Rydberg
-from pulser.parametrized import Variable, ParamObj
+from pulser.json.coders import PulserDecoder, PulserEncoder
+from pulser.parametrized import ParamObj, Variable
 from pulser.waveforms import (
-    ConstantWaveform,
-    KaiserWaveform,
-    RampWaveform,
     BlackmanWaveform,
-    CustomWaveform,
     CompositeWaveform,
+    ConstantWaveform,
+    CustomWaveform,
     InterpolatedWaveform,
+    KaiserWaveform,
+    RampWaveform,
 )
 
 np.random.seed(20201105)