Example Gallery

The scripts in examples/gallery are complete, runnable starting points for common open-system workflows. Each script accepts --fast for a small smoke run, --output-dir for generated figures/data, and --json for a compact machine-readable summary.

Run an example from a source checkout with:

python examples/gallery/qubit_decay.py --fast

The default output location is runs/example_gallery/<example-name>.

Deterministic Qubit Decay

Solves spontaneous emission with mesolve, saves the result as HDF5, plots the excited-state population, and compares against exp(-gamma * t).

examples/gallery/qubit_decay.py

"""Deterministic spontaneous emission of a two-level system."""

from __future__ import annotations

import argparse
from pathlib import Path

import numpy as np

import openquantumsim as oqs

try:
    from ._common import (
        Summary,
        add_common_arguments,
        configure_matplotlib,
        emit_summary,
        example_output_dir,
    )
except ImportError:  # pragma: no cover - used when executed as a script
    from _common import (  # type: ignore[no-redef]
        Summary,
        add_common_arguments,
        configure_matplotlib,
        emit_summary,
        example_output_dir,
    )


def run_example(
    *,
    output_dir: Path | None = None,
    fast: bool = False,
    save_plots: bool = True,
) -> Summary:
    """Run spontaneous emission and compare against the analytic curve."""
    atom = oqs.SpinSpace(0.5, label="atom")
    excited = oqs.basis(atom, "up")

    gamma = 0.25
    t_final = 0.4 if fast else 8.0
    time_points = 5 if fast else 201
    times = np.linspace(0.0, t_final, time_points)

    H = 0.0 * oqs.sigmaz(atom)
    rho0 = oqs.ket2dm(excited)
    collapse = np.sqrt(gamma) * oqs.sigmam(atom)
    projector = oqs.Operator(oqs.ket2dm(excited), atom, "P_excited")

    result = oqs.mesolve(
        H,
        rho0,
        times,
        c_ops=[collapse],
        e_ops=[projector],
        options=oqs.Options(rtol=1e-8, atol=1e-10),
    )

    expected = np.exp(-gamma * times)
    population = result.expect[0].real
    max_error = float(np.max(np.abs(population - expected)))

    target = example_output_dir("qubit_decay", output_dir)
    result.save_hdf5(target / "qubit_decay.h5")
    if save_plots:
        plt = configure_matplotlib()
        fig, ax = plt.subplots(figsize=(6.0, 3.6))
        ax.plot(times, population, label="OpenQuantumSim")
        ax.plot(times, expected, "--", label="analytic")
        ax.set_xlabel("time")
        ax.set_ylabel("excited population")
        ax.legend()
        fig.tight_layout()
        fig.savefig(target / "qubit_decay.png", dpi=160)
        plt.close(fig)

    return {
        "example": "qubit_decay",
        "time_points": len(times),
        "max_abs_error": max_error,
        "final_population": float(population[-1]),
        "output_dir": str(target),
    }


def main() -> None:
    """Run the example from the command line."""
    parser = argparse.ArgumentParser(description=__doc__)
    add_common_arguments(parser)
    args = parser.parse_args()
    summary = run_example(
        output_dir=args.output_dir,
        fast=args.fast,
        save_plots=not args.no_plots,
    )
    emit_summary(summary, as_json=args.json)


if __name__ == "__main__":
    main()

Driven Qubit

Builds a time-dependent Hamiltonian from an interpolated drive envelope and tracks the excited population and Bloch components of a dissipative qubit.

examples/gallery/driven_qubit.py

"""Time-dependent resonant drive of a dissipative two-level system."""

from __future__ import annotations

import argparse
from pathlib import Path

import numpy as np

import openquantumsim as oqs

try:
    from ._common import (
        Summary,
        add_common_arguments,
        configure_matplotlib,
        emit_summary,
        example_output_dir,
    )
except ImportError:  # pragma: no cover - used when executed as a script
    from _common import (  # type: ignore[no-redef]
        Summary,
        add_common_arguments,
        configure_matplotlib,
        emit_summary,
        example_output_dir,
    )


def run_example(
    *,
    output_dir: Path | None = None,
    fast: bool = False,
    save_plots: bool = True,
) -> Summary:
    """Run a qubit with a smooth pulsed drive and weak spontaneous emission."""
    atom = oqs.SpinSpace(0.5, label="atom")
    ground = oqs.basis(atom, "down")
    excited = oqs.basis(atom, "up")

    gamma = 0.04
    omega = 1.2
    t_final = 1.0 if fast else 10.0
    time_points = 7 if fast else 251
    times = np.linspace(0.0, t_final, time_points)

    envelope = oqs.InterpolatedCoefficient(
        [0.0, 0.25 * t_final, 0.75 * t_final, t_final],
        [0.0, 1.0, 1.0, 0.0],
    )
    H = oqs.time_dependent_hamiltonian(
        0.0 * oqs.sigmaz(atom),
        [(0.5 * omega * oqs.sigmax(atom), envelope)],
    )
    rho0 = oqs.ket2dm(ground)
    collapse = np.sqrt(gamma) * oqs.sigmam(atom)
    excited_projector = oqs.Operator(oqs.ket2dm(excited), atom, "P_excited")

    result = oqs.mesolve(
        H,
        rho0,
        times,
        c_ops=[collapse],
        e_ops=[excited_projector, oqs.sigmax(atom), oqs.sigmaz(atom)],
        options=oqs.Options(rtol=1e-8, atol=1e-10),
    )

    target = example_output_dir("driven_qubit", output_dir)
    result.save_hdf5(target / "driven_qubit.h5")
    if save_plots:
        plt = configure_matplotlib()
        fig, ax = plt.subplots(figsize=(6.0, 3.6))
        ax.plot(times, result.expect[0].real, label="P(excited)")
        ax.plot(times, result.expect[1].real, label="<sigma_x>")
        ax.plot(times, result.expect[2].real, label="<sigma_z>")
        ax.set_xlabel("time")
        ax.set_ylabel("expectation")
        ax.legend()
        fig.tight_layout()
        fig.savefig(target / "driven_qubit.png", dpi=160)
        plt.close(fig)

    return {
        "example": "driven_qubit",
        "time_points": len(times),
        "final_excited_population": float(result.expect[0].real[-1]),
        "peak_excited_population": float(np.max(result.expect[0].real)),
        "output_dir": str(target),
    }


def main() -> None:
    """Run the example from the command line."""
    parser = argparse.ArgumentParser(description=__doc__)
    add_common_arguments(parser)
    args = parser.parse_args()
    summary = run_example(
        output_dir=args.output_dir,
        fast=args.fast,
        save_plots=not args.no_plots,
    )
    emit_summary(summary, as_json=args.json)


if __name__ == "__main__":
    main()

Jaynes-Cummings Dynamics

Uses the built-in Jaynes-Cummings system helper to simulate energy exchange between a truncated cavity and a two-level atom.

examples/gallery/jaynes_cummings.py

"""Damped Jaynes-Cummings dynamics in a truncated cavity."""

from __future__ import annotations

import argparse
from pathlib import Path

import numpy as np

import openquantumsim as oqs

try:
    from ._common import (
        Summary,
        add_common_arguments,
        configure_matplotlib,
        emit_summary,
        example_output_dir,
    )
except ImportError:  # pragma: no cover - used when executed as a script
    from _common import (  # type: ignore[no-redef]
        Summary,
        add_common_arguments,
        configure_matplotlib,
        emit_summary,
        example_output_dir,
    )


def run_example(
    *,
    output_dir: Path | None = None,
    fast: bool = False,
    save_plots: bool = True,
) -> Summary:
    """Run a small damped Jaynes-Cummings model and plot energy exchange."""
    cavity_dim = 3 if fast else 8
    t_final = 0.8 if fast else 30.0
    time_points = 7 if fast else 301
    times = np.linspace(0.0, t_final, time_points)

    system = oqs.jaynes_cummings_system(
        cavity_dim,
        cavity_frequency=1.0,
        atom_frequency=1.0,
        coupling=0.08,
        cavity_decay=0.02,
        atom_decay=0.01,
        initial_photon=0,
        atom_state="up",
    )
    result = oqs.mesolve(
        system.H,
        system.rho0,
        times,
        c_ops=system.c_ops,
        e_ops=system.e_ops,
        options=oqs.Options(rtol=1e-8, atol=1e-10),
    )

    target = example_output_dir("jaynes_cummings", output_dir)
    result.save_hdf5(target / "jaynes_cummings.h5")
    if save_plots:
        plt = configure_matplotlib()
        fig, ax = plt.subplots(figsize=(6.0, 3.6))
        ax.plot(times, result.expect[0].real, label="<n_cavity>")
        ax.plot(times, result.expect[1].real, label="P(atom excited)")
        ax.set_xlabel("time")
        ax.set_ylabel("expectation")
        ax.legend()
        fig.tight_layout()
        fig.savefig(target / "jaynes_cummings.png", dpi=160)
        plt.close(fig)

    return {
        "example": "jaynes_cummings",
        "time_points": len(times),
        "cavity_dim": cavity_dim,
        "max_photon_number": float(np.max(result.expect[0].real)),
        "final_atom_excited_population": float(result.expect[1].real[-1]),
        "output_dir": str(target),
    }


def main() -> None:
    """Run the example from the command line."""
    parser = argparse.ArgumentParser(description=__doc__)
    add_common_arguments(parser)
    args = parser.parse_args()
    summary = run_example(
        output_dir=args.output_dir,
        fast=args.fast,
        save_plots=not args.no_plots,
    )
    emit_summary(summary, as_json=args.json)


if __name__ == "__main__":
    main()

Quantum Trajectories

Runs a finite Monte Carlo wave-function ensemble and plots the sample mean with standard-error bands.

examples/gallery/quantum_trajectory.py

"""Monte Carlo wave-function simulation of finite trajectory noise."""

from __future__ import annotations

import argparse
from pathlib import Path

import numpy as np

import openquantumsim as oqs

try:
    from ._common import (
        Summary,
        add_common_arguments,
        configure_matplotlib,
        emit_summary,
        example_output_dir,
    )
except ImportError:  # pragma: no cover - used when executed as a script
    from _common import (  # type: ignore[no-redef]
        Summary,
        add_common_arguments,
        configure_matplotlib,
        emit_summary,
        example_output_dir,
    )


def run_example(
    *,
    output_dir: Path | None = None,
    fast: bool = False,
    save_plots: bool = True,
) -> Summary:
    """Estimate spontaneous emission from a finite trajectory ensemble."""
    atom = oqs.SpinSpace(0.5, label="atom")
    excited = oqs.basis(atom, "up")

    gamma = 0.3
    t_final = 0.4 if fast else 6.0
    time_points = 5 if fast else 151
    n_traj = 12 if fast else 500
    times = np.linspace(0.0, t_final, time_points)

    H = 0.0 * oqs.sigmaz(atom)
    collapse = np.sqrt(gamma) * oqs.sigmam(atom)
    projector = oqs.Operator(oqs.ket2dm(excited), atom, "P_excited")
    result = oqs.mcsolve(
        H,
        excited,
        times,
        c_ops=[collapse],
        e_ops=[projector],
        n_traj=n_traj,
        options=oqs.Options(seed=2026, max_step=0.02, n_jobs=1),
    )

    population = result.expect[0].real
    stderr = result.expect_stderr[0]
    expected = np.exp(-gamma * times)
    mean_abs_error = float(np.mean(np.abs(population - expected)))

    target = example_output_dir("quantum_trajectory", output_dir)
    result.save_hdf5(target / "quantum_trajectory.h5")
    if save_plots:
        plt = configure_matplotlib()
        fig, ax = plt.subplots(figsize=(6.0, 3.6))
        ax.plot(times, population, label=f"{n_traj} trajectories")
        ax.fill_between(
            times,
            population - 2.0 * stderr,
            population + 2.0 * stderr,
            alpha=0.25,
            label="2 standard errors",
        )
        ax.plot(times, expected, "--", label="analytic")
        ax.set_xlabel("time")
        ax.set_ylabel("excited population")
        ax.legend()
        fig.tight_layout()
        fig.savefig(target / "quantum_trajectory.png", dpi=160)
        plt.close(fig)

    return {
        "example": "quantum_trajectory",
        "time_points": len(times),
        "n_traj": n_traj,
        "mean_abs_error": mean_abs_error,
        "final_standard_error": float(stderr[-1]),
        "output_dir": str(target),
    }


def main() -> None:
    """Run the example from the command line."""
    parser = argparse.ArgumentParser(description=__doc__)
    add_common_arguments(parser)
    args = parser.parse_args()
    summary = run_example(
        output_dir=args.output_dir,
        fast=args.fast,
        save_plots=not args.no_plots,
    )
    emit_summary(summary, as_json=args.json)


if __name__ == "__main__":
    main()

Phase-Space Distributions

Computes Wigner and Husimi-Q functions for a finite Fock-space cat state. This example is pure Python and does not require loading the Julia backend.

examples/gallery/phase_space.py

"""Wigner and Husimi-Q functions for a finite Fock-space state."""

from __future__ import annotations

import argparse
from pathlib import Path

import numpy as np

import openquantumsim as oqs

try:
    from ._common import (
        Summary,
        add_common_arguments,
        configure_matplotlib,
        emit_summary,
        example_output_dir,
    )
except ImportError:  # pragma: no cover - used when executed as a script
    from _common import (  # type: ignore[no-redef]
        Summary,
        add_common_arguments,
        configure_matplotlib,
        emit_summary,
        example_output_dir,
    )


def run_example(
    *,
    output_dir: Path | None = None,
    fast: bool = False,
    save_plots: bool = True,
) -> Summary:
    """Evaluate Wigner and Q functions for a coherent-state superposition."""
    space = oqs.FockSpace(12 if fast else 30, label="cavity")
    alpha = 1.2 + 0.4j
    cat = oqs.coherent(space, alpha) + oqs.coherent(space, -alpha)
    cat = cat / np.linalg.norm(cat)
    rho = oqs.ket2dm(cat)

    grid_points = 31 if fast else 151
    x, p = oqs.phase_space_grid(xlim=(-4.0, 4.0), points=grid_points)
    wigner = oqs.wigner(rho, x, p)
    q_values = oqs.q_function(rho, x, p)
    negativity = float(np.sum(np.abs(wigner) - wigner) * (x[1] - x[0]) * (p[1] - p[0]))

    target = example_output_dir("phase_space", output_dir)
    np.savez(target / "phase_space.npz", x=x, p=p, wigner=wigner, q_function=q_values)
    if save_plots:
        plt = configure_matplotlib()
        fig, axes = plt.subplots(1, 2, figsize=(8.0, 3.6), constrained_layout=True)
        vmax = float(np.max(np.abs(wigner)))
        axes[0].imshow(
            wigner,
            origin="lower",
            extent=(float(x[0]), float(x[-1]), float(p[0]), float(p[-1])),
            cmap="RdBu_r",
            vmin=-vmax,
            vmax=vmax,
            aspect="auto",
        )
        axes[0].set_title("Wigner")
        axes[0].set_xlabel("x")
        axes[0].set_ylabel("p")
        axes[1].imshow(
            q_values,
            origin="lower",
            extent=(float(x[0]), float(x[-1]), float(p[0]), float(p[-1])),
            cmap="viridis",
            aspect="auto",
        )
        axes[1].set_title("Husimi Q")
        axes[1].set_xlabel("x")
        fig.savefig(target / "phase_space.png", dpi=160)
        plt.close(fig)

    return {
        "example": "phase_space",
        "fock_dim": space.dim,
        "grid_points": grid_points,
        "wigner_min": float(np.min(wigner)),
        "q_max": float(np.max(q_values)),
        "wigner_negativity": negativity,
        "output_dir": str(target),
    }


def main() -> None:
    """Run the example from the command line."""
    parser = argparse.ArgumentParser(description=__doc__)
    add_common_arguments(parser)
    args = parser.parse_args()
    summary = run_example(
        output_dir=args.output_dir,
        fast=args.fast,
        save_plots=not args.no_plots,
    )
    emit_summary(summary, as_json=args.json)


if __name__ == "__main__":
    main()

Parameter Sweep

Runs a restartable sweep over a qubit decay rate, writing a manifest, one HDF5 result per point, summary.csv, and summary.h5.

examples/gallery/parameter_sweep.py

"""Restartable sweep over a qubit decay rate."""

from __future__ import annotations

import argparse
from pathlib import Path

import numpy as np

import openquantumsim as oqs

try:
    from ._common import (
        Summary,
        add_common_arguments,
        configure_matplotlib,
        emit_summary,
        example_output_dir,
    )
except ImportError:  # pragma: no cover - used when executed as a script
    from _common import (  # type: ignore[no-redef]
        Summary,
        add_common_arguments,
        configure_matplotlib,
        emit_summary,
        example_output_dir,
    )


def run_example(
    *,
    output_dir: Path | None = None,
    fast: bool = False,
    save_plots: bool = True,
) -> Summary:
    """Sweep spontaneous-emission rates and write restartable outputs."""
    target = example_output_dir("parameter_sweep", output_dir)
    times = np.linspace(0.0, 0.4 if fast else 4.0, 5 if fast else 101)
    gamma_values = [0.1, 0.3] if fast else [0.05, 0.1, 0.2, 0.4]

    sweep = oqs.ParameterSweep(
        base_system={"model": "qubit_decay", "time_points": len(times)},
        params={"gamma": gamma_values},
    )
    run = sweep.run(
        lambda point: _run_decay_point(point, times),
        output_dir=target,
        summarize=_summarize_decay_point,
        restart=True,
    )

    final_values = [float(row["final_population"]) for row in run.summary]
    if save_plots:
        plt = configure_matplotlib()
        fig, ax = plt.subplots(figsize=(5.4, 3.4))
        ax.plot(gamma_values, final_values, marker="o")
        ax.set_xlabel("decay rate gamma")
        ax.set_ylabel("final excited population")
        fig.tight_layout()
        fig.savefig(target / "parameter_sweep.png", dpi=160)
        plt.close(fig)

    return {
        "example": "parameter_sweep",
        "points": len(run.summary),
        "min_final_population": float(min(final_values)),
        "max_final_population": float(max(final_values)),
        "manifest": str(run.manifest_path),
        "summary_csv": str(run.summary_csv_path),
        "output_dir": str(target),
    }


def _run_decay_point(point: oqs.SweepPoint, times: np.ndarray) -> oqs.Result:
    gamma = float(point.params["gamma"])
    atom = oqs.SpinSpace(0.5, label="atom")
    excited = oqs.basis(atom, "up")
    H = 0.0 * oqs.sigmaz(atom)
    collapse = np.sqrt(gamma) * oqs.sigmam(atom)
    projector = oqs.Operator(oqs.ket2dm(excited), atom, "P_excited")
    return oqs.mesolve(
        H,
        oqs.ket2dm(excited),
        times,
        c_ops=[collapse],
        e_ops=[projector],
        options=oqs.Options(rtol=1e-8, atol=1e-10),
    )


def _summarize_decay_point(
    point: oqs.SweepPoint,
    result: object,
) -> dict[str, float | int | str]:
    if not isinstance(result, oqs.Result):
        msg = "parameter-sweep example expected an OpenQuantumSim Result."
        raise TypeError(msg)
    population = result.expect[0].real
    expected = np.exp(-float(point.params["gamma"]) * result.times)
    return {
        "id": point.id,
        "gamma": float(point.params["gamma"]),
        "final_population": float(population[-1]),
        "max_abs_error": float(np.max(np.abs(population - expected))),
    }


def main() -> None:
    """Run the example from the command line."""
    parser = argparse.ArgumentParser(description=__doc__)
    add_common_arguments(parser)
    args = parser.parse_args()
    summary = run_example(
        output_dir=args.output_dir,
        fast=args.fast,
        save_plots=not args.no_plots,
    )
    emit_summary(summary, as_json=args.json)


if __name__ == "__main__":
    main()