### Install and set up pre-commit

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/community/contributing.md

Installs pre-commit and sets it up to automatically run cleaning tasks (ruff + codespell) before each commit.

```shell
pip install pre-commit
pre-commit install
```

--------------------------------

### Install dependencies with uv

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/community/contributing.md

Installs the library with all dependencies and developer dependencies using uv.

```shell
uv sync --extra dev
```

--------------------------------

### Install Dynamiqs from source (GitHub)

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/getting_started/installation.md

Install the latest development version directly from the GitHub repository.

```shell
pip install git+https://github.com/dynamiqs/dynamiqs.git
```

--------------------------------

### Install Dynamiqs with pip

Source: https://github.com/dynamiqs/dynamiqs/blob/main/README_PYPI.md

Install the Dynamiqs package using pip. For GPU support, refer to the JAX installation documentation for device-specific instructions.

```shell
pip install dynamiqs
```

--------------------------------

### Preview documentation with live reload

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/community/contributing.md

Starts a live preview server that rebuilds the documentation website upon modifications.

```shell
task docserve
```

--------------------------------

### Install dependencies with pip in editable mode

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/community/contributing.md

Installs the library in editable mode with developer dependencies using pip. A virtual environment is strongly recommended.

```shell
pip install -e ".[dev]"
```

--------------------------------

### Import Dynamiqs and JAX NumPy

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/performance.md

Imports the Dynamiqs library and JAX NumPy, required for all examples on this page.

```python
import dynamiqs as dq
import jax.numpy as jnp
```

--------------------------------

### Output of dq.floquet simulation

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/floquet-integration.md

Example output from the previous simulation. Shows the progress bar, the quasienergies [0.972, -0.972], and the final Floquet modes as a QArray with shape (2, 2, 1) and complex64 dtype.

```text
|██████████| 100.0% ◆ elapsed 1.99ms ◆ remaining 0.00ms
[ 0.972 -0.972]
QArray: shape=(2, 2, 1), dims=(2,), dtype=complex64, layout=dense
[[[ 0.998+0.j]
  [-0.055-0.j]]

 [[ 0.055-0.j]
  [ 0.998+0.j]]]
```

--------------------------------

### Time a simulation with jax.block_until_ready

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/performance.md

Times a JAX-compiled simulation, accounting for compilation on the first call and asynchronous dispatch. Use jax.block_until_ready on the result before stopping the clock to get accurate timing.

```python
import jax

simulate = jax.jit(lambda: dq.sesolve(H, psi0, tsave, progress_meter=False))
_ = jax.block_until_ready(simulate())  # compile
result = jax.block_until_ready(simulate())  # time this one
```

--------------------------------

### Simulate Rabi oscillations for a Kerr oscillator in the transmon regime

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/advanced_examples/kerr-oscillator.md

Simulates a driven-dissipative Kerr oscillator in the transmon regime (K >> |epsilon| >> kappa) using dq.mesolve, starting from vacuum. It computes the populations of Fock states |0> and |1> over time to observe damped Rabi oscillations and leakage to higher Fock states. The plot is saved with renderfig('rabi-oscillations-kerr-oscillator').

```python
# simulation parameters
n = 8
K = 200.0
epsilon = 40.0
kappa = 1.0
T = 10 * jnp.pi / epsilon
ntsave = 401

# operators
a = dq.destroy(n)
H = -K * a.dag() @ a.dag() @ a @ a + epsilon * (a + a.dag())
jump_ops = [jnp.sqrt(kappa) * a]

# initial state
psi0 = dq.basis(n, 0)

# save times
tsave = jnp.linspace(0.0, T, ntsave)

# expectation operator
exp_ops = [dq.proj(dq.basis(n, 0)), dq.proj(dq.basis(n, 1))]

# run simulation and extract observables
result = dq.mesolve(H, jump_ops, psi0, tsave, exp_ops=exp_ops)
pop_0 = result.expects[0].real  # population of |0>
pop_1 = result.expects[1].real  # population of |1>

# plot Rabi oscillations
plt.plot(tsave * epsilon / jnp.pi, pop_0, label=r'$|0\rangle$')
plt.plot(tsave * epsilon / jnp.pi, pop_1, label=r'$|1\rangle$')
plt.plot(
    tsave * epsilon / jnp.pi,
    1 - (pop_0 + pop_1),
    color='black',
    label=r'$|2\rangle$, $|3\rangle$, $\ldots$',
)
plt.xlabel(r'Time, $t\epsilon / \pi$')
plt.ylabel('Population')
plt.ylim(0, 1)
plt.xlim(0, T * epsilon / jnp.pi)
plt.legend(frameon=True)
renderfig('rabi-oscillations-kerr-oscillator')
```

--------------------------------

### Simulate a lossy quantum harmonic oscillator with dq.mesolve

Source: https://github.com/dynamiqs/dynamiqs/blob/main/README_PYPI.md

Simulates a lossy harmonic oscillator with Hamiltonian H = omega * a^dagger a and jump operator sqrt(kappa) * a, from time 0 to T, starting from a coherent state. Uses dq.mesolve and prints the result, which includes method, steps, and state shape.

```python
import dynamiqs as dq
import jax.numpy as jnp

# parameters
n = 16  # Hilbert space dimension
omega = 1.0  # frequency
kappa = 0.1  # decay rate
alpha0 = 1.0  # initial coherent state amplitude
T = 2 * jnp.pi  # total evolution time (one full revolution)

# initialize operators, initial state and saving times
a = dq.destroy(n)
H = omega * dq.dag(a) @ a
jump_ops = [jnp.sqrt(kappa) * a]
psi0 = dq.coherent(n, alpha0)
tsave = jnp.linspace(0, T, 101)

# run simulation
result = dq.mesolve(H, jump_ops, psi0, tsave)
print(result)
```

```text
|██████████| 100.0% ◆ elapsed 6.30ms ◆ remaining 0.00ms
==== MESolveResult ====
Method : Tsit5
Infos  : 40 steps (40 accepted, 0 rejected)
States : QArray complex64 (101, 16, 16) | 202.0 Kb
```

--------------------------------

### Apply TEST_INSTANT priority marker

Source: https://github.com/dynamiqs/dynamiqs/blob/main/CLAUDE.md

Example of applying a priority marker to a test, using the TEST_INSTANT tier from tests/order.py. Every test must carry a tier marker, either per test or via a module-level pytestmark.

```python
@pytest.mark.run(order=TEST_INSTANT)
```

--------------------------------

### Plot PWC pulse with dq.plot.pwc_pulse()

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/time-dependent-operators.md

Plot the PWC pulse values using dq.plot.pwc_pulse() with H.times and H.values. The example also sets the x-axis label and renders the figure with renderfig().

```python
dq.plot.pwc_pulse(H.times, H.values)
plt.gca().set(xlabel='Time $t$')
renderfig('tqarray_plot_pwc_pulse')
```

--------------------------------

### Create PWC operator with dq.pwc(times, values, qarray)

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/time-dependent-operators.md

Construct a piecewise constant (PWC) operator with dq.pwc() using times, values, and a qarray. The example defines H(t) = 3*sigma_z for t in [0,1) and -2*sigma_z for t in [1,2). The returned object is a PWCTimeQArray.

```pycon
>>> times = [0.0, 1.0, 2.0]
>>> values = [3.0, -2.0]
>>> qarray = dq.sigmaz()
>>> H = dq.pwc(times, values, qarray)
>>> H
PWCTimeQArray: shape=(2, 2), dims=(2,), dtype=complex64, layout=dia, ndiags=1
```

--------------------------------

### Build the static documentation website

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/community/contributing.md

Builds the static documentation website locally, creating a site/ directory.

```shell
task docbuild
```

--------------------------------

### Simulate a two-level system with sesolve

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/closed-systems.md

Demonstrates creating an initial state and Hamiltonian to run a simulation, followed by the resulting output.

```python
import jax.numpy as jnp
import dynamiqs as dq

jnp.set_printoptions(precision=3, suppress=True)  # set custom array print style

psi0 = dq.ground()  # initial state
H = dq.sigmaz()  # Hamiltonian
tsave = jnp.linspace(0, 1.0, 11)  # saving times
res = dq.sesolve(H, psi0, tsave)  # run the simulation
print(res.states[-1])  # print the final state
```

```text
|██████████| 100.0% ◆ elapsed 2.52ms ◆ remaining 0.00ms
QArray: shape=(2, 1), dims=(2,), dtype=complex64, layout=dense
[[0.  +0.j   ]
 [0.54+0.841j]]
```

--------------------------------

### Run dynamiqs benchmarks with uv run task bench

Source: https://github.com/dynamiqs/dynamiqs/blob/main/benchmarks/README.md

Run the benchmark suite with various options: select tiers, filter by case key, use quick mode, set repeats, or write results to a JSON file. Default tier is both, default repeats is 5.

```shell
uv run task bench                     # both tiers (~2 minutes on a laptop CPU)
uv run task bench --tier physics      # workload regressions only (~1 minute)
uv run task bench --tier features     # feature comparisons only (~1 minute)
uv run task bench --quick             # tiny sizes, a few seconds (sanity check)
uv run task bench --filter cavity     # only cases whose key contains "cavity"
uv run task bench --out results.json  # also write results to a JSON file
```

--------------------------------

### Set matplotlib style with dq.plot.mplstyle

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/advanced_examples/kerr-oscillator.md

Sets the plotting style for Dynamiqs, with higher resolution, for the figures in this example.

```python
import dynamiqs as dq
import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt

dq.plot.mplstyle(dpi=150)  # set custom matplotlib style
```

--------------------------------

### Clone the repository and navigate into it

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/community/contributing.md

Clones the Dynamiqs repository and changes into the project directory.

```shell
git clone git@github.com:dynamiqs/dynamiqs.git
cd dynamiqs
```

--------------------------------

### List available tasks

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/community/contributing.md

Displays the list of available tasks with descriptions.

```shell
> task --list
lint         lint the code (ruff)
format       auto-format the code (ruff)
codespell    check for misspellings (codespell)
type         check types (ty)
clean        clean the code (ruff + codespell + ty)
test         run the unit tests suite (pytest)
doctest-code check code docstrings examples (doctest)
doctest-docs check documentation examples (doctest)
doctest      check all examples (doctest)
docbuild     build the documentation website
docserve     preview documentation website with hot-reloading
all          run all tasks before a commit (ruff + codespell + ty + pytest + doctest)
ci           run all the CI checks
```

--------------------------------

### Shift a timeqarray with H.shift(0.5)

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/time-dependent-operators.md

Shifts the time dependence of a timeqarray using the .shift() method. The example shows that H_shift(0.0) equals H(0.5), and demonstrates the resulting qarray values.

```pycon
>>> f = lambda t: jnp.cos(2.0 * jnp.pi * t)
>>> H = dq.modulated(f, dq.sigmax())
>>> H_shift = H.shift(0.5)
>>> H_shift(0.0)  # == H(0.5)
QArray: shape=(2, 2), dims=(2,), dtype=complex64, layout=dia, ndiags=2
[[   ⋅    -1.+0.j]
 [-1.+0.j    ⋅   ]]
>>> H(0.5)
QArray: shape=(2, 2), dims=(2,), dtype=complex64, layout=dia, ndiags=2
[[   ⋅    -1.+0.j]
 [-1.+0.j    ⋅   ]]
```

--------------------------------

### Import Dynamiqs and JAX

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/batching-simulations.md

Imports the Dynamiqs library, JAX NumPy, and the timeit module for timing simulations.

```python
import dynamiqs as dq
import jax.numpy as jnp
import timeit
```

--------------------------------

### Run dq.mesolve() on a two-level system

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/open-systems.md

Sets up a two-level system with a Hamiltonian and a single jump operator, then runs a master equation simulation with dq.mesolve() and prints the final state. The output shows the progress bar and the resulting density matrix.

```python
import jax.numpy as jnp
import dynamiqs as dq

jnp.set_printoptions(precision=3, suppress=True)  # set custom array print style

psi0 = dq.excited()  # initial state
H = dq.sigmaz()  # Hamiltonian
jump_ops = [dq.sigmam()]  # list of jump operators
tsave = jnp.linspace(0, 1.0, 11)  # saving times
res = dq.mesolve(H, jump_ops, psi0, tsave)  # run the simulation
print(res.states[-1])  # print the final state
```

```text
|██████████| 100.0% ◆ elapsed 1.75ms ◆ remaining 0.00ms
QArray: shape=(2, 2), dims=(2,), dtype=complex64, layout=dense
[[0.368+0.j 0.   +0.j]
 [0.   +0.j 0.632+0.j]]
```

--------------------------------

### Run simulation with dq.sesolve and print result

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/workflow.md

Runs the simulation with `dq.sesolve()` and prints the result type, attributes, and a summary. The output shows the `SESolveResult` object with method `Dopri5`, 56 steps (48 accepted, 8 rejected), and the shapes of `States` and `Expects` arrays.

```python
# run simulation
result = dq.sesolve(H, psi0, tsave, exp_ops=exp_ops, method=method)

# print some information
print(f'`result` is of type {type(result)}.')
print(f'`result` has the following attributes:')
print(f'{[attr for attr in dir(result) if not attr.startswith("__")]}\n')
print(result)
```

```text
`result` is of type <class 'dynamiqs.result.SESolveResult'>.
`result` has the following attributes:
['_abc_impl', '_saved', '_str_parts', 'block_until_ready', 'expects', 'extra', 'final_state', 'gradient', 'infos', 'method', 'options', 'states', 'to_numpy', 'to_qutip', 'tsave']

|██████████| 100.0% ◆ elapsed 4.57ms ◆ remaining 0.00ms
==== SESolveResult ====
Method  : Dopri5
Infos   : 56 steps (48 accepted, 8 rejected)
States  : QArray complex64 (101, 2, 1) | 1.6 Kb
Expects : Array complex64 (1, 101) | 0.8 Kb
```

--------------------------------

### Clip a timeqarray with H.clip(0, 1)

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/time-dependent-operators.md

Clips a time-dependent operator to a given time interval, outside which the returned qarray is null. The example creates a modulated operator H, clips it to 0 <= t < 1, and shows that evaluating the clipped operator at t=2.0 returns a zero qarray.

```pycon
>>> f = lambda t: jnp.cos(2.0 * jnp.pi * t)
>>> H = dq.modulated(f, dq.sigmax())
>>> H(2.0)
QArray: shape=(2, 2), dims=(2,), dtype=complex64, layout=dia, ndiags=2
[[  ⋅    1.+0.j]
 [1.+0.j   ⋅   ]]
>>> H_clipped = H.clip(0, 1)  # clip to 0 <= t < 1
>>> H_clipped(2.0)
QArray: shape=(2, 2), dims=(2,), dtype=complex64, layout=dia, ndiags=2
[[  ⋅      ⋅   ]
 [  ⋅      ⋅   ]]
```

--------------------------------

### Simulate diffusive SSE trajectories for a qubit with dq.dssesolve()

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/advanced_examples/continuous-diffusive-measurement.md

Sets up a qubit with no Hamiltonian, a sigma_z jump operator, and a pure initial state. Uses EulerMaruyama method and dq.dssesolve() to simulate 5 trajectories. Prints the result object.

```python
# define Hamiltonian, jump operators, initial state
H = dq.zeros(2)
jump_ops = [dq.sigmaz()]
psi0 = (dq.ground() + dq.excited()).unit()

# define save times
tsave = np.linspace(0, 1.0, 101)
delta_t = tsave[1] - tsave[0]

# define a certain number of PRNG key, one for each trajectory
key = jax.random.PRNGKey(20)
ntrajs = 5
keys = jax.random.split(key, ntrajs)

# simulate trajectories
method = dq.method.EulerMaruyama(dt=1e-3)
result = dq.dssesolve(H, jump_ops, psi0, tsave, keys, method=method)
print(result)
```

```text
==== DSSESolveResult ====
Method       : EulerMaruyama
Infos        : 1000 steps | infos shape (5,)
States       : QArray complex64 (5, 101, 2, 1) | 7.9 Kb
Measurements : Array float32 (5, 1, 100) | 2.0 Kb
```

```pycon
>>> result.states.shape  # (ntrajs, ntsave, n, 1)
(5, 101, 2, 1)
>>> result.measurements.shape  # (ntrajs, nLm, ntsave-1)
(5, 1, 100)
```

--------------------------------

### Benchmark and compare two runs

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/performance.md

Runs the benchmark suite twice with the physics tier, saving results to before.json and after.json, then compares the two runs with the benchmarks compare command. Useful for measuring performance changes before and after a modification.

```shell
uv run task bench --tier physics --out before.json
uv run task bench --tier physics --out after.json
python -m benchmarks compare before.json after.json
```

--------------------------------

### Run narrow pytest targets with uv run

Source: https://github.com/dynamiqs/dynamiqs/blob/main/CLAUDE.md

Run the narrowest relevant test target instead of the full suite. Each command targets a single file, a single concern via -k, or a single module's docstring examples. Prefix with uv run as required by the environment.

```shell
uv run pytest tests/sesolve/test_adaptive.py -q     # one file
uv run pytest tests/mesolve -q -k gradient          # one concern
uv run pytest dynamiqs/utils/general.py -q          # one module's docstring examples
```

--------------------------------

### Simulate monitored oscillator with dq.dsmesolve

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/advanced_examples/continuous-diffusive-measurement.md

Sets up and runs the diffusive SME simulation for 1000 trajectories of a damped harmonic oscillator. Defines the Hamiltonian, jump operators, efficiencies, initial coherent state, and save times, then calls dq.dsmesolve with the Euler-Maruyama method and prints the result.

```python
# define Hamiltonian, jump operators, efficiencies, initial state
n = 16
a = dq.destroy(n)
kappa = 1.0
omega = 10.0
alpha0 = 2.0
H = omega * a.dag() @ a
jump_ops = [jnp.sqrt(kappa / 2) * a, jnp.sqrt(kappa / 2) * (-1j * a)]
etas = [1.0, 1.0]
psi0 = dq.coherent(n, alpha0)

# define save times
tsave = np.linspace(0, 1 / kappa, 101)
delta_t = tsave[1] - tsave[0]

# define a certain number of PRNG key, one for each trajectory
key = jax.random.PRNGKey(42)
ntrajs = 1000
keys = jax.random.split(key, ntrajs)

# simulate trajectories
method = dq.method.EulerMaruyama(dt=1e-3)
result = dq.dsmesolve(
    H, jump_ops, etas, psi0, tsave, keys, method=method, save_states=False
)
print(result)
```

--------------------------------

### Benchmark batched vs unbatched simulation with sesolve

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/batching-simulations.md

Defines 3000 Hamiltonians as a batched array and compares running sesolve in a nested for loop versus a single batched call. The output shows the batched simulation is about 60 times faster.

```ipython
# define 3000 Hamiltonians
omega = jnp.linspace(0.0, 10.0, 100)[:, None, None]
epsilon = jnp.linspace(0.0, 1.0, 30)[:, None, None, None]
H = omega * dq.sigmaz() + epsilon * dq.sigmax()  # (100, 30, 2, 2)

# other simulation parameters
psi0 = dq.basis(2, 0)
tsave = jnp.linspace(0.0, 1.0, 50)

# running the simulations successively
def run_unbatched():
    results = []
    for i in range(len(omega)):
        for j in range(len(epsilon)):
            result = dq.sesolve(H[i, j], psi0, tsave, progress_meter=False)
            results.append(result)
    return results

# running the simulations simultaneously
def run_batched():
    return dq.sesolve(H, psi0, tsave, progress_meter=False)

# exclude JIT time from benchmarking by running each function once first
%timeit -n1 -r1 -q run_unbatched()
%timeit -n1 -r1 -q run_batched()

# time functions
%timeit run_unbatched()
%timeit run_batched()
```

```text
2.59 s ± 52.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
44.1 ms ± 2.66 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
```

--------------------------------

### Raise ValueError with argument validation message

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/community/contributing.md

Use this pattern for exceptions arising from incorrect argument input. Messages should start with a capital letter, end with a period, and follow the format "Argument ... must ..., but ...". Use backticks for variable/function names and single quotes for strings.

```python
raise ValueError(
    f'Argument `H` must have shape (n, n), but has shape H.shape={H.shape}.'
)
raise ValueError(
    "Argument `matmul_precision` should be a string 'low', 'high', or"
    f" 'highest', but is '{matmul_precision}'."
)
```

--------------------------------

### Import Dynamiqs

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/getting_started/sharp-bits.md

Imports the Dynamiqs library, typically the first step in any script using Dynamiqs.

```python
import dynamiqs as dq
```

--------------------------------

### Simulate Kerr oscillator with dq.mesolve

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/advanced_examples/kerr-oscillator.md

Simulates the time evolution of a driven-dissipative Kerr oscillator with constant driving and coherent initial state, using dq.mesolve.

```python
# simulation parameters
n = 32  # Hilbert space size
K = 1.0  # Kerr non-linearity
epsilon = 3.0  # driving field
kappa = 1.5  # dissipation rate
alpha0 = 2.0  # initial coherent state amplitude
T = 5.0  # simulation time
ntsave = 201  # number of saved states

# operators
a = dq.destroy(n)
H = -K * a.dag() @ a.dag() @ a @ a + epsilon * (a + a.dag())
jump_ops = [jnp.sqrt(kappa) * a]

# initial state
psi0 = dq.coherent(n, alpha0)

# save times
tsave = jnp.linspace(0.0, T, ntsave)

# run simulation
result = dq.mesolve(H, jump_ops, psi0, tsave)
```

--------------------------------

### Repository layout overview

Source: https://github.com/dynamiqs/dynamiqs/blob/main/CLAUDE.md

Shows the top-level layout of the dynamiqs repository, including the main package directories, docs, tests, and benchmarks. Use this to understand where different components live, such as qarrays, integrators, and utilities.

```text
dynamiqs/
  qarrays/          # QArray: the quantum array type (dense + sparse DIA layouts)
  time_qarray.py    # TimeQArray: constant / pwc / modulated / timecallable
  integrators/
    apis/           # public solver entry points (sesolve, mesolve, ...)
    core/           # integrator implementations (diffrax, expm, rouchon, ...)
  utils/  plot/  random/
  method.py  gradient.py  options.py  result.py
  conftest.py       # sybil setup: docstring examples are executed as tests
docs/               # mkdocs sources; docs/conftest.py executes the markdown examples
tests/  benchmarks/
```

--------------------------------

### Batch sesolve over Hamiltonians and initial states

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/batching-simulations.md

Batches a Schrödinger equation simulation over three Hamiltonians and four initial states. The result.states shape (3, 4, 11, 2, 1) shows the batched dimensions: 3 Hamiltonians, 4 states, 11 saved times, and (2, 1) per state. All result attributes, such as expects when exp_ops is given, are batched accordingly.

```python
# define three Hamiltonians
H = [dq.sigmax(), dq.sigmay(), dq.sigmaz()]  # (3, 2, 2)

# define four initial states
g = dq.basis(2, 0)
e = dq.basis(2, 1)
plus = (g + e).unit()
minus = (g - e).unit()
psi = [g, e, plus, minus]  # (4, 2, 1)

# run the simulation
tsave = jnp.linspace(0.0, 1.0, 11)  # (11,)
result = dq.sesolve(H, psi, tsave)
print(f'Shape of result.states: {result.states.shape}')
```

```text
Shape of result.states: (3, 4, 11, 2, 1)
```

--------------------------------

### Simulate larger systems with TF_FORCE_UNIFIED_MEMORY and XLA_PYTHON_CLIENT_MEM_FRACTION

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/performance.md

Set TF_FORCE_UNIFIED_MEMORY and XLA_PYTHON_CLIENT_MEM_FRACTION to allow XLA to spill GPU allocations into system memory, enabling simulation of systems larger than GPU memory. XLA_PYTHON_CLIENT_MEM_FRACTION values greater than 1 let allocations overflow into host memory, trading speed for capacity. Works on NVIDIA GPUs, especially with a fast interconnect like Grace Hopper.

```shell
TF_FORCE_UNIFIED_MEMORY=1 XLA_PYTHON_CLIENT_MEM_FRACTION=4 python simulation.py
```

--------------------------------

### Compute the propagator with sepropagator

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/closed-systems.md

Shows how to compute the propagator for a system, continuing from the previous Hamiltonian and time definitions, followed by the output.

```python
res = dq.sepropagator(H, tsave)
print(res.propagators[-1])  # print the final propagator
```

```text
QArray: shape=(2, 2), dims=(2,), dtype=complex64, layout=dense
[[0.54-0.841j 0.  +0.j   ]
 [0.  +0.j    0.54+0.841j]]
```

--------------------------------

### Simulate Rabi model with dq.floquet and print quasienergies

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/floquet-integration.md

Simulates the Rabi model of a driven two-level system using Dynamiqs' floquet function. Sets a custom array print style, defines the drive period and frequency, constructs the drift and drive Hamiltonians, and saves times over one period. Prints the quasienergies and the final Floquet modes.

```python
import jax.numpy as jnp
import dynamiqs as dq

jnp.set_printoptions(precision=3, suppress=True)  # set custom array print style

T = 1.0  # drive period
omega_d = 2 * jnp.pi / T  # drive frequency
H0 = 0.1 * dq.sigmaz()  # drift Hamiltonian
f = lambda t: jnp.cos(omega_d * t)  # time-dependent drive
H1 = dq.modulated(f, dq.sigmax())  # drive Hamiltonian
tsave = jnp.linspace(0, T, 11)  # saving times
res = dq.floquet(H0 + H1, T, tsave)  # run the simulation
print(res.quasienergies)  # print the quasienergies
print(res.modes[-1])  # print the final floquet modes
```

--------------------------------

### Simulate diffusive SME trajectories with imperfect detection using dq.dsmesolve()

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/advanced_examples/continuous-diffusive-measurement.md

Sets up the same qubit but with an efficiency of 0.2 and uses dq.dsmesolve() to simulate mixed-state trajectories. Prints the result object.

```python
# define efficiencies
etas = [0.2]

# simulate trajectories
result = dq.dsmesolve(H, jump_ops, etas, psi0, tsave, keys, method=method)
print(result)
```

```text
==== DSMESolveResult ====
Method       : EulerMaruyama
Infos        : 1000 steps | infos shape (5,)
States       : QArray complex64 (5, 101, 2, 2) | 7.9 Kb
Measurements : Array float32 (5, 1, 100) | 2.0 Kb
```

```pycon
>>> result.states.shape  # (ntrajs, ntsave, n, n)
(5, 101, 2, 2)
>>> result.measurements.shape  # (ntrajs, nLm, ntsave-1)
(5, 1, 100)
```

--------------------------------

### Compare two benchmark runs (A/B)

Source: https://github.com/dynamiqs/dynamiqs/blob/main/benchmarks/README.md

Run benchmarks on two git branches and compare the results. The comparison table aligns cases by name and parameters and reports the change in median run time (positive = slower). The JSON files record device, precision, package versions, and git SHA; compare warns if device or precision differ.

```shell
git checkout main
uv run task bench --tier physics --out before.json
git checkout my-branch
uv run task bench --tier physics --out after.json
python -m benchmarks compare before.json after.json
```

--------------------------------

### Progress Meters

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/python_api/index.md

Progress meters for monitoring solver progress.

```APIDOC
## Progress Meters

### NoProgressMeter
No progress meter (silent).

### TextProgressMeter
Text-based progress meter.

### TqdmProgressMeter
Tqdm-based progress meter.
```

--------------------------------

### JSSESolveResult output

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/advanced_examples/continuous-jump-measurement.md

Output of the JSSESolveResult for the jump SSE simulation. Shows method, number of steps, states shape (3, 101, 2, 1) and clicktimes shape (3, 2, 10000).

```text
==== JSSESolveResult ====
Method     : EulerJump
Infos      : 1000 steps | infos shape (3,)
States     : QArray complex64 (3, 101, 2, 1) | 7.9 Kb
Clicktimes : Array float32 (3, 2, 10000) | 390.6 Kb
```

--------------------------------

### Simulate jump SSE trajectories with dq.jssesolve()

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/advanced_examples/continuous-jump-measurement.md

Simulates jump SSE trajectories for a qubit with perfect detection using dq.jssesolve(). The Hamiltonian is zero, jump operators are sigma_minus and sigma_plus, and the initial state is the excited state. Uses EulerJump method with dt=1e-3 and 3 trajectories.

```python
# define Hamiltonian, jump operators, initial state
H = dq.zeros(2)
jump_ops = [dq.sigmam(), dq.sigmap()]
psi0 = dq.excited()

# define save times
tsave = np.linspace(0, 2.0, 101)

# define a certain number of PRNG key, one for each trajectory
key = jax.random.PRNGKey(42)
ntrajs = 3
keys = jax.random.split(key, ntrajs)

# simulate trajectories
method = dq.method.EulerJump(dt=1e-3)
result = dq.jssesolve(H, jump_ops, psi0, tsave, keys, method=method)
print(result)
```

--------------------------------

### Compute propagator with dq.mepropagator()

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/basics/open-systems.md

Computes the propagator for the same Hamiltonian and jump operators using dq.mepropagator() and prints the final propagator. The output shows a 4x4 complex matrix with vectorized layout.

```python
res = dq.mepropagator(H, jump_ops, tsave)
print(res.propagators[-1])  # print the final propagator
```

```text
QArray: shape=(4, 4), dims=(2,), dtype=complex64, layout=dense, vectorized=True
[[ 0.368+0.j     0.   +0.j     0.   +0.j     0.   +0.j   ]
 [ 0.   +0.j    -0.252+0.552j  0.   +0.j     0.   +0.j   ]
 [ 0.   +0.j     0.   +0.j    -0.252-0.552j  0.   +0.j   ]
 [ 0.632+0.j     0.   +0.j     0.   +0.j     1.   +0.j   ]]
```

--------------------------------

### Sweep parameters and compute pi-pulse fidelity

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/advanced_examples/kerr-oscillator.md

Sets up the simulation parameters, defines the Hamiltonian with a modulated drive, and uses jax.vmap to compute the fidelity for a range of gate durations. The fidelity is estimated as the maximum overlap with the |1> state.

```python
from functools import partial

# simulation parameters
n = 8
K = 200.0
kappa = 1.0
ntsave = 401

# parameters to sweep
Ts = jnp.linspace(0.05, 0.5, 24)
sigmas = jnp.linspace(0.05, 0.2, 14)

# operators, initial state, and expectation operator
a = dq.destroy(n)
H0 = -K * a.dag() @ a.dag() @ a @ a
jump_ops = [jnp.sqrt(kappa) * a]
psi0 = dq.basis(n, 0)
exp_ops = [dq.proj(dq.basis(n, 0)), dq.proj(dq.basis(n, 1))]


@jax.vmap
def compute_fidelity(T):
    """Compute the fidelity of a pi-pulse for a given gate duration."""
    # time-dependent Hamiltonian, defined with functools.partial and broadcasting
    # `f` has signature (t: float) -> Array of shape (len(sigmas),)
    f = partial(pulse, T=T, sigma=sigmas)
    H = H0 + dq.modulated(f, a + a.dag())

    # save times
    tsave = jnp.linspace(0.0, T, ntsave)

    # run simulation
    result = dq.mesolve(H, jump_ops, psi0, tsave, exp_ops=exp_ops)

    # fidelity is roughly estimated as the largest overlap with |1>
    # in a proper study, we would need to compute a full process tomography
    return jnp.max(result.expects[:, 1, :].real, axis=-1)


# run simulation
fidelities = compute_fidelity(Ts)

# plot results
plt.pcolormesh(sigmas, Ts, jnp.log10(1 - fidelities), cmap='Blues_r')
plt.xlabel('Pulse width')
plt.ylabel('Gate duration')
plt.colorbar(label='log10(1 - Fidelity)')
plt.grid(False)
renderfig('pipulse-kerr-oscillator')
```

--------------------------------

### Simulate jump SME trajectories with dq.jsmesolve()

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/advanced_examples/continuous-jump-measurement.md

Simulates jump SME trajectories for a qubit with imperfect detection using dq.jsmesolve(). Dark count rates are set to 0.0 and efficiencies to 0.5 for both jump operators. Uses the same EulerJump method and keys as the perfect detection case.

```python
# define dark count rates
thetas = [0.0, 0.0]
# define efficiencies
etas = [0.5, 0.5]

# simulate trajectories
result = dq.jsmesolve(H, jump_ops, thetas, etas, psi0, tsave, keys, method=method)
print(result)
```

--------------------------------

### JSMESolveResult output

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/advanced_examples/continuous-jump-measurement.md

Output of the JSMESolveResult for the jump SME simulation. Shows method, number of steps, states shape (3, 101, 2, 2) and clicktimes shape (3, 2, 10000).

```text
==== JSMESolveResult ====
Method     : EulerJump
Infos      : 1000 steps | infos shape (3,)
States     : QArray complex64 (3, 101, 2, 2) | 9.5 Kb
Clicktimes : Array float32 (3, 2, 10000) | 234.4 Kb
```

--------------------------------

### Simulate quantum harmonic oscillator with continuous jump measurement

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/advanced_examples/continuous-jump-measurement.md

Defines the Hamiltonian, jump operator, initial Fock state, and simulation parameters for a quantum harmonic oscillator with a single loss channel. Uses the EulerJump method to simulate 100 trajectories and prints the result summary.

```python
# define Hamiltonian, jump operators, initial state
n = 16
a = dq.destroy(n)
kappa = 1.0
n0 = 10
H = dq.zeros(n)
jump_ops = [jnp.sqrt(kappa) * a]
psi0 = dq.fock(n, n0)

# define save times
tsave = np.linspace(0, 2.0, 101)

# define a certain number of PRNG key, one for each trajectory
key = jax.random.PRNGKey(42)
ntrajs = 100
keys = jax.random.split(key, ntrajs)

# simulate trajectories
method = dq.method.EulerJump(dt=1e-3)
exp_ops = [a.dag() @ a]
result = dq.jssesolve(
    H, jump_ops, psi0, tsave, keys, method=method, save_states=False, exp_ops=exp_ops
)
print(result)
```

```text
==== JSMESolveResult ====
Method     : EulerJump
Infos      : 1000 steps | infos shape (1000,)
States     : QArray complex64 (100, 16, 1) | 2.0 Mb
Expects    : Array complex64 (100, 1, 101) | 789.1 Kb
Clicktimes : Array float32 (100, 1, 10000) | 38.1 Mb
```

--------------------------------

### Compute gradients with jax.grad and dq.mesolve

Source: https://github.com/dynamiqs/dynamiqs/blob/main/README.md

Computes gradients of the final photon number with respect to omega, kappa, and alpha0 using jax.grad on a function that runs dq.mesolve. The output shows the gradient values; the note explains the analytical verification with the coherent state amplitude formula.

```python
import dynamiqs as dq
import jax.numpy as jnp
import jax

# parameters
n = 16  # Hilbert space dimension
omega = 1.0  # frequency
kappa = 0.1  # decay rate
alpha0 = 1.0  # initial coherent state amplitude
T = 2 * jnp.pi  # total evolution time (one full revolution)


def population(omega, kappa, alpha0):
    """Return the oscillator population after time evolution."""
    # initialize operators, initial state and saving times
    a = dq.destroy(n)
    H = omega * a.dag() @ a
    jump_ops = [jnp.sqrt(kappa) * a]
    psi0 = dq.coherent(n, alpha0)
    tsave = jnp.linspace(0, T, 101)

    # run simulation
    result = dq.mesolve(H, jump_ops, psi0, tsave)

    return dq.expect(dq.number(n), result.states[-1]).real


# compute gradient with respect to omega, kappa and alpha
grad_population = jax.grad(population, argnums=(0, 1, 2))
grads = grad_population(omega, kappa, alpha0)
print(f'Gradient w.r.t. omega : {grads[0]:.4f}')
print(f'Gradient w.r.t. kappa : {grads[1]:.4f}')
print(f'Gradient w.r.t. alpha0: {grads[2]:.4f}')
```

```text
|██████████| 100.0% ◆ elapsed 5.94ms ◆ remaining 0.00ms
Gradient w.r.t. omega : 0.0000
Gradient w.r.t. kappa : -3.3520
Gradient w.r.t. alpha0: 1.0670
```

--------------------------------

### Create Wigner function GIF with dq.plot.wigner_gif

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/advanced_examples/kerr-oscillator.md

Generates a GIF of the Wigner function evolution and renders it to a file. The 'rendergif' function is used for display in the documentation.

```python
gif = dq.plot.wigner_gif(result.states, ymax=3.0, gif_duration=10.0)
rendergif(gif, 'wigner-kerr-oscillator')
```

--------------------------------

### Simulate and plot Kerr oscillator revival

Source: https://github.com/dynamiqs/dynamiqs/blob/main/docs/documentation/advanced_examples/kerr-oscillator.md

Simulates the time evolution of a Kerr oscillator with photon loss and plots the absolute value of the cavity field over time. Uses Dynamiqs' mesolve with a jump operator and an expectation operator, and renders the figure with renderfig.

```python
# simulation parameters
n = 32
K = 1.0
kappa = 0.02
alpha0 = 2.0
ntsave = 201

# operators
a = dq.destroy(n)
H = -K * a.dag() @ a.dag() @ a @ a
jump_ops = [jnp.sqrt(kappa) * a]

# initial state
psi0 = dq.coherent(n, alpha0)

# save times
T = 5 * jnp.pi / K
tsave = jnp.linspace(0.0, T, ntsave)

# expectation operator
exp_ops = [a]

# run simulation
result = dq.mesolve(H, jump_ops, psi0, tsave, exp_ops=exp_ops)

# plot cavity field
plt.plot(tsave * K / jnp.pi, jnp.abs(result.expects[0]))
plt.xlabel(r'Time, $tK / \pi$')
plt.ylabel(r'$|\langle a(t) \rangle|$')
renderfig('photon-number-kerr-oscillator')
```