diff --git a/exercises/cv_submissions/cv_tutorial_pijus.ipynb b/exercises/cv_submissions/cv_tutorial_pijus.ipynb new file mode 100644 index 0000000..ccee954 --- /dev/null +++ b/exercises/cv_submissions/cv_tutorial_pijus.ipynb @@ -0,0 +1,698 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "5c61d218", + "metadata": {}, + "source": [ + "## Phase space and the Wigner function\n", + "\n", + "Classically, the state of an oscillator is a point $(x,p)$ in phase space, and an ensemble\n", + "with some uncertainty is a probability distribution over that plane. A quantum state does not have the same kind of picture, $x$ and $p$\n", + "have no simultaneous sharp values, so there is nothing for a joint distribution $p(x,p)$ to\n", + "be a distribution *over*. \n", + "\n", + "The standard workaround is the Wigner function, a *quasi*probability distribution:\n", + "\n", + "$$W_\\rho(x,p) = \\frac{1}{2\\pi}\\int_{-\\infty}^{\\infty}\\Big\\langle x + \\tfrac{q}{2}\\Big|\\,\\hat\\rho\\,\\Big|x - \\tfrac{q}{2}\\Big\\rangle e^{ipq}\\,dq$$\n", + "\n", + "### It behaves like a probability distribution:\n", + "\n", + "- It is real — no imaginary parts to explain away.\n", + "- It is normalised: $\\int W\\,dx\\,dp = 1$. Every figure below prints this as a numerical sanity check.\n", + "\n", + "### And in one crucial way it does not\n", + "\n", + "$W$ can go negative. No probability ever can. When it happens it is not a numerical artefact — it is the signature of a state with no classical description at all.\n", + "\n", + "States whose Wigner function is Gaussian (vacuum, coherent, squeezed, thermal) are called Gaussian states. They are the easy ones to make and to calculate with. Everything else — Fock states, cat states — is hard to produce in a lab, and is exactly what is valuable for quantum technology." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Imports OK ✓\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.gridspec as gridspec\n", + "from matplotlib import cm, colors\n", + "from matplotlib.ticker import MaxNLocator\n", + "from scipy.special import factorial\n", + "import warnings\n", + "warnings.filterwarnings('ignore')\n", + "\n", + "try:\n", + " import ipywidgets as w\n", + "except ImportError:\n", + " %pip install -q ipywidgets\n", + " import ipywidgets as w\n", + "import io, pickle, multiprocessing\n", + "import concurrent.futures as cf\n", + "from pathlib import Path\n", + "from IPython.display import display\n", + "\n", + "%matplotlib inline\n", + "plt.rcParams.update({\n", + " 'figure.dpi': 120,\n", + " 'font.family': 'DejaVu Sans',\n", + " 'axes.titlesize': 11,\n", + " 'axes.labelsize': 10,\n", + "})\n", + "print('Imports OK ✓')" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "23a50aac", + "metadata": {}, + "outputs": [], + "source": [ + "GRID = 120 # grid points per axis (increase for higher resolution)\n", + "XLIM = 5.0 # phase-space extent ±XLIM\n", + "\n", + "xvec = np.linspace(-XLIM, XLIM, GRID)\n", + "pvec = np.linspace(-XLIM, XLIM, GRID)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1 Core Quantum-Optics Utilities" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "595c7101", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.special import hermite as H_phys\n", + "\n", + "# ─────────────────────────────────────────────────────────────────────────────\n", + "# Position-space wavefunction ψ_n(x) of the quantum harmonic oscillator\n", + "# ψ_n(x) = (2^n n! √π)^{-1/2} · H_n(x) · exp(-x²/2)\n", + "# where H_n is the physicists' Hermite polynomial.\n", + "# ─────────────────────────────────────────────────────────────────────────────\n", + "def qho_wavefunction(n, x):\n", + " \"\"\"Harmonic-oscillator eigenfunction ψ_n evaluated on array x.\"\"\"\n", + " norm = 1.0 / np.sqrt(2.0**n * factorial(n, exact=False) * np.sqrt(np.pi))\n", + " return norm * H_phys(n)(x) * np.exp(-0.5 * x**2)\n", + "\n", + "\n", + "# ─────────────────────────────────────────────────────────────────────────────\n", + "# Wigner function via direct integration (Wigner 1932 definition)\n", + "#\n", + "# W(x, p) = (1/π) ∫_{-∞}^{∞} ⟨x+y|ρ|x-y⟩ e^{2ipy} dy\n", + "#\n", + "# In the Fock basis:\n", + "# ⟨x+y|ρ|x-y⟩ = Σ_{mn} ρ_{mn} ψ_m(x+y) ψ_n*(x-y)\n", + "#\n", + "# Algorithm:\n", + "# 1. Build y-integration grid.\n", + "# 2. Precompute ψ_n(x+y) and ψ_n(x-y) for all n, x, y → shape (dim, Nx, Ny)\n", + "# 3. Contract with ρ to form the kernel K(x,y) = ⟨x+y|ρ|x-y⟩ → (Nx, Ny)\n", + "# 4. For each p: W(x,p) = (1/π) Re[∫ K(x,y) e^{2ipy} dy] via np.trapz\n", + "# ─────────────────────────────────────────────────────────────────────────────\n", + "def wigner_function(rho, xvec, pvec, Ny=400, y_extent=None):\n", + " \"\"\"\n", + " Compute W(x, p) by numerically integrating the Wigner definition.\n", + "\n", + " Parameters\n", + " ----------\n", + " rho : (N, N) complex array – density matrix in Fock basis\n", + " xvec : 1-D real array – position quadrature grid\n", + " pvec : 1-D real array – momentum quadrature grid\n", + " Ny : int – number of points on the y integration grid (default 400)\n", + " y_extent : float – half-width of y integration domain.\n", + " Defaults to max(|xvec|, |pvec|) + 3 to capture the tails.\n", + "\n", + " Returns\n", + " -------\n", + " W : (len(pvec), len(xvec)) real array\n", + " \"\"\"\n", + " dim = rho.shape[0]\n", + " Nx = len(xvec)\n", + " Np = len(pvec)\n", + "\n", + " if y_extent is None:\n", + " y_extent = float(max(np.abs(xvec).max(), np.abs(pvec).max()) + 3.0)\n", + "\n", + " yvec = np.linspace(-y_extent, y_extent, Ny) # integration variable\n", + "\n", + " # ── Step 1: grids x±y with shape (Nx, Ny) ────────────────────────────\n", + " xp = xvec[:, None] + yvec[None, :] # x + y\n", + " xm = xvec[:, None] - yvec[None, :] # x - y\n", + "\n", + " # ── Step 2: precompute ψ_n on both grids → (dim, Nx, Ny) ─────────────\n", + " psi_plus = np.stack([qho_wavefunction(n, xp) for n in range(dim)]) # ψ_n(x+y)\n", + " psi_minus = np.stack([qho_wavefunction(n, xm) for n in range(dim)]) # ψ_n(x-y)\n", + "\n", + " # ── Step 3: kernel K(x,y) = Σ_{mn} ρ_{mn} ψ_m(x+y) ψ_n*(x-y) ────────\n", + " # einsum: 'mn, mxy, nxy -> xy' where x=position index, y=integration index\n", + " # K is genuinely complex whenever ρ has imaginary coherences (e.g. |0>+i|1>);\n", + " # keep it complex — Hermiticity K(x,-y) = K(x,y)* makes the y-integral real.\n", + " K = np.einsum('mn, mxy, nxy -> xy', rho, psi_plus, psi_minus.conj(),\n", + " optimize='optimal') # shape (Nx, Ny), complex\n", + "\n", + " # ── Step 4: integrate W(x,p) = (1/π) ∫ K(x,y) e^{2ipy} dy over y ───\n", + " # exp_factor shape: (Np, Ny); K shape: (Nx, Ny)\n", + " # integrand: K[x,y] * exp(2i p y) → broadcast to (Nx, Np, Ny)\n", + " exp_factor = np.exp(2j * pvec[:, None] * yvec[None, :]) # (Np, Ny)\n", + " integrand = K[:, None, :] * exp_factor[None, :, :] # (Nx, Np, Ny)\n", + " W = np.real(np.trapezoid(integrand, yvec, axis=-1)) / np.pi # (Nx, Np)\n", + "\n", + " return W.T # → (Np, Nx) to match meshgrid convention W[p_idx, x_idx]" + ] + }, + { + "cell_type": "markdown", + "id": "22fb755b", + "metadata": {}, + "source": [ + "### First: what is $\\hat\\rho$?\n", + "\n", + "The code below builds **density matrices**, written $\\hat\\rho$. This is the general way to write down a quantum state.\n", + "\n", + "- If you know the state exactly — a **pure** state $|\\psi\\rangle$ — then $\\hat\\rho = |\\psi\\rangle\\langle\\psi|$.\n", + "- If you only know a *classical* probabilistic mixture — \"half the time it's $|0\\rangle$, half the time it's $|1\\rangle$\" — then $\\hat\\rho = \\tfrac12|0\\rangle\\langle0| + \\tfrac12|1\\rangle\\langle1|$. This is a **mixed** state and it cannot be written as any single $|\\psi\\rangle$.\n", + "\n", + "The test to see if the state is pure or mixed: $\\mathrm{Tr}(\\hat\\rho^2)$: 1 for a pure state, less than 1 for a mixed one. Every figure below prints it.\n", + "\n", + "All matrices here are written in the **Fock basis**: the row and column indices *are* photon numbers, so $\\rho_{mn} = \\langle m|\\hat\\rho|n\\rangle$. The `dim` argument is where we chop off the infinite ladder; it must sit comfortably above the largest photon number the state actually contains, or the state gets silently mangled." + ] + }, + { + "cell_type": "markdown", + "id": "6484391e", + "metadata": {}, + "source": [ + "### Fock states $|n\\rangle$\n", + "\n", + "A Fock state (or *number state*) is a single rung of the energy ladder: exactly $n$ photons, with no uncertainty in energy whatsoever. It is built by applying the creation operator $n$ times to the vacuum:\n", + "\n", + "$$|n\\rangle = \\frac{(\\hat a^\\dagger)^n}{\\sqrt{n!}}\\,|0\\rangle$$\n", + "\n", + "Fock state Wigner function is a perfectly **rotationally symmetric** set of rings,\n", + "\n", + "$$W_{|n\\rangle}(x,p) = \\frac{1}{\\pi}e^{-x^2-p^2}(-1)^n L_n(2x^2+2p^2)$$\n", + "\n", + "with $L_n$ the $n$-th Laguerre polynomial.\n", + "\n", + "\n", + "### Coherent states $|\\alpha\\rangle$ \n", + "\n", + "A coherent state is as close as quantum mechanics gets to a classical wave, and it is what comes out of a laser. Three equivalent ways to think about it:\n", + "\n", + "1. Removing a photon changes nothing. Formally, $|\\alpha\\rangle$ is an eigenstate of the annihilation operator:\n", + "\n", + "$$\\hat a |\\alpha\\rangle = \\alpha|\\alpha\\rangle$$\n", + "\n", + "This is why a laser beam does not change character as it is absorbed along the way — take a photon out and you have the same state back.\n", + "\n", + "2. It is a displaced vacuum. $|\\alpha\\rangle = \\hat D(\\alpha)|0\\rangle$. The complex number $\\alpha$ *is* the location in phase space:\n", + "\n", + "$$\\langle \\hat x\\rangle = \\sqrt2\\,\\mathrm{Re}\\,\\alpha, \\qquad \\langle \\hat p\\rangle = \\sqrt2\\,\\mathrm{Im}\\,\\alpha$$\n", + "\n", + "$\\alpha = 0$ gives back the vacuum exactly.\n", + "\n", + "3. Poissonian distribution of Fock states. In the photon-number basis:\n", + "\n", + "$$|\\alpha\\rangle = e^{-|\\alpha|^2/2}\\sum_{n=0}^{\\infty}\\frac{\\alpha^n}{\\sqrt{n!}}\\,|n\\rangle$$\n", + "\n", + "so the photon number is not definite. Counting photons gives a Poisson distribution with mean $\\langle\\hat n\\rangle = |\\alpha|^2$. " + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "def coherent_state_rho(alpha_c, dim=30):\n", + " \"\"\"Density matrix of coherent state |alpha_c>.\"\"\"\n", + " n = np.arange(dim)\n", + " psi = np.exp(-0.5 * np.abs(alpha_c)**2) * (alpha_c**n) / np.sqrt(factorial(n))\n", + " psi /= np.linalg.norm(psi)\n", + " return np.outer(psi, psi.conj())\n", + "\n", + "\n", + "def fock_state_rho(n_fock, dim=30):\n", + " \"\"\"Density matrix of Fock (number) state |n>.\"\"\"\n", + " rho = np.zeros((dim, dim), dtype=complex)\n", + " rho[n_fock, n_fock] = 1.0\n", + " return rho\n", + "\n", + "\n", + "def cat_state_rho(alpha_c, dim=40):\n", + " \"\"\"Even cat state ∝ |alpha> + |-alpha>.\"\"\"\n", + " n = np.arange(dim)\n", + " def coh(a):\n", + " psi = np.exp(-0.5*np.abs(a)**2) * (a**n) / np.sqrt(factorial(n))\n", + " return psi\n", + " psi = coh(alpha_c) + coh(-alpha_c)\n", + " psi /= np.linalg.norm(psi)\n", + " return np.outer(psi, psi.conj())\n", + "\n", + "\n", + "def thermal_state_rho(n_bar, dim=30):\n", + " \"\"\"Thermal state with mean photon number n_bar.\"\"\"\n", + " n = np.arange(dim)\n", + " diag = (n_bar**n) / ((1 + n_bar)**(n + 1))\n", + " diag /= diag.sum()\n", + " return np.diag(diag.astype(complex))\n", + "\n", + "\n", + "def squeezed_vacuum_rho(r, dim=40):\n", + " \"\"\"Squeezed vacuum state with squeezing parameter r.\"\"\"\n", + " n = np.arange(dim)\n", + " # Only even Fock states contribute\n", + " psi = np.zeros(dim, dtype=complex)\n", + " tanh_r = np.tanh(r)\n", + " sech_r = 1.0 / np.cosh(r)\n", + " for k in range(0, dim, 2):\n", + " psi[k] = (np.sqrt(factorial(k)) / (2**(k/2) * factorial(k//2))) \\\n", + " * ((-tanh_r)**(k//2)) * np.sqrt(sech_r)\n", + " psi /= np.linalg.norm(psi)\n", + " return np.outer(psi, psi.conj())\n", + "\n", + "\n", + "def visualize_state(state:str, density=None):\n", + " if state == 'coherent':\n", + " rho = coherent_state_rho(alpha_c=2.0, dim=30)\n", + " state_label = r'Coherent State $|\\alpha{=}2\\rangle$'\n", + "\n", + " elif state == 'fock':\n", + " rho = fock_state_rho(n_fock=3, dim=20)\n", + " state_label = r'Fock State $|n{=}3\\rangle$'\n", + "\n", + " elif state == 'cat':\n", + " rho = cat_state_rho(alpha_c=2.0, dim=40)\n", + " state_label = r'Even Cat State $|\\alpha{=}2\\rangle + |{-}\\alpha\\rangle$'\n", + "\n", + " elif state == 'thermal':\n", + " rho = thermal_state_rho(n_bar=2.0, dim=30)\n", + " state_label = r'Thermal State $\\bar{n}{=}2$'\n", + "\n", + " elif state == 'squeezed':\n", + " rho = squeezed_vacuum_rho(r=1.0, dim=40)\n", + " state_label = r'Squeezed Vacuum $r{=}1$'\n", + " elif state == 'custom':\n", + " rho = density\n", + " state_label = r'custom density matrix'\n", + "\n", + " else:\n", + " raise ValueError(f'Unknown STATE={state!r}')\n", + "\n", + " dim = verify_rho(rho, state_label)\n", + " W, dx = compute_wigner(rho)\n", + " visualize(rho, state_label, W, dim, dx)\n", + " \n", + "def verify_rho(rho, state_label):\n", + " assert np.allclose(rho, rho.conj().T, atol=1e-10), \"rho must be Hermitian!\"\n", + " assert np.isclose(np.trace(rho).real, 1.0, atol=1e-8), \"Tr(rho) must equal 1!\"\n", + " eigvals = np.linalg.eigvalsh(rho)\n", + " assert np.all(eigvals > -1e-8), \"rho must be positive semi-definite!\"\n", + "\n", + " dim = rho.shape[0]\n", + " print(f'State: {state_label}')\n", + " print(f'Hilbert space dimension: {dim}')\n", + " print(f'Tr(rho) = {np.trace(rho).real:.6f}')\n", + " print(f'Tr(rho²) = {np.trace(rho @ rho).real:.6f} (1 = pure, <1 = mixed)')\n", + " return dim\n", + "\n", + "\n", + "def compute_wigner(rho):\n", + " print(f'Computing Wigner function on {GRID}×{GRID} grid … ', end='', flush=True)\n", + " W = wigner_function(rho, xvec, pvec)\n", + "\n", + " # Normalisation check: ∫W dx dp = 1 (Δx·Δp ≈ dx²)\n", + " dx = xvec[1] - xvec[0]\n", + " norm_W = np.sum(W) * dx**2\n", + " print(f'∫W dx dp = {norm_W:.4f} (should be ≈ 1.0)')\n", + " print(f'W_min = {W.min():.4f}, W_max = {W.max():.4f}')\n", + " if W.min() < 0:\n", + " print(' → Negative values detected: state is non-classical')\n", + " return W, dx" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "b696ab2e", + "metadata": {}, + "outputs": [], + "source": [ + "def visualize(rho, state_label, W, dim, dx, grid=None):\n", + " # grid=(x, p) draws on a window other than the notebook-wide xvec/pvec; the\n", + " # interactive explorer needs a per-state window and renders off-thread.\n", + " xvec, pvec = grid if grid is not None else (globals()['xvec'], globals()['pvec'])\n", + " W_ABS_MAX = np.max(np.abs(W)) * 1.05\n", + " CMAP_W = 'RdBu_r'\n", + " CMAP_RHO = 'seismic'\n", + "\n", + " X_grid, P_grid = np.meshgrid(xvec, pvec)\n", + "\n", + " # ─────────────────────────────────────────────────────────────────────────────\n", + " # FIGURE LAYOUT\n", + " # Row 0 : [3-D Wigner] [2-D Wigner] [Fock distribution]\n", + " # Row 1 : [Re(rho)] [Im(rho)] [Marginals]\n", + " # ─────────────────────────────────────────────────────────────────────────────\n", + " fig = plt.figure(figsize=(18, 11))\n", + " fig.suptitle(f'Quantum State Visualization — {state_label}',\n", + " fontsize=14, fontweight='bold', y=0.98)\n", + "\n", + " gs = gridspec.GridSpec(2, 3, figure=fig,\n", + " hspace=0.45, wspace=0.38,\n", + " left=0.06, right=0.97,\n", + " top=0.93, bottom=0.06)\n", + "\n", + " # ── 4-A 3-D Wigner surface ───────────────────────────────────────────────────\n", + " ax3d = fig.add_subplot(gs[0, 0], projection='3d')\n", + "\n", + " norm_surf = colors.Normalize(vmin=-W_ABS_MAX, vmax=W_ABS_MAX)\n", + " \n", + " ax3d.contourf(X_grid, P_grid, W, zdir='z',\n", + " offset=W.min() - 0.05 * (W.max() - W.min()),\n", + " levels=30, cmap=CMAP_W, alpha=0.5)\n", + " ax3d.set_xlabel('x (Position)', labelpad=4, fontsize=9)\n", + " ax3d.set_ylabel('p (Momentum)', labelpad=4, fontsize=9)\n", + " ax3d.set_zlabel('W(x, p)', labelpad=4, fontsize=9)\n", + " ax3d.set_title('3-D Wigner Function', pad=6)\n", + " ax3d.view_init(elev=28, azim=-55)\n", + " ax3d.tick_params(labelsize=7)\n", + " m = cm.ScalarMappable(cmap=CMAP_W, norm=norm_surf)\n", + " m.set_array([])\n", + " fig.colorbar(m, ax=ax3d, shrink=0.5, pad=0.12, label='W')\n", + "\n", + " # ── 4-B 2-D Wigner contour map ───────────────────────────────────────────────\n", + " ax2d = fig.add_subplot(gs[0, 1])\n", + " im_w = ax2d.pcolormesh(X_grid, P_grid, W,\n", + " cmap=CMAP_W,\n", + " vmin=-W_ABS_MAX, vmax=W_ABS_MAX,\n", + " shading='auto')\n", + " ax2d.set_xlabel('x (Position)')\n", + " ax2d.set_ylabel('p (Momentum)')\n", + " ax2d.set_title('Wigner Function W(x, p)')\n", + " ax2d.set_aspect('equal')\n", + " ax2d.axhline(0, color='k', lw=0.5, ls='--', alpha=0.4)\n", + " ax2d.axvline(0, color='k', lw=0.5, ls='--', alpha=0.4)\n", + " cb_w = fig.colorbar(im_w, ax=ax2d, fraction=0.046, pad=0.04)\n", + " cb_w.set_label('W(x, p)')\n", + "\n", + " # ── 4-C Fock-state distribution ─────────────────────────────────────────────\n", + " ax_fock = fig.add_subplot(gs[0, 2])\n", + " fock_probs = np.real(np.diag(rho))\n", + " n_vals = np.arange(dim)\n", + " ax_fock.set_xlabel('Fock number n')\n", + " ax_fock.set_ylabel('Probability P(n)')\n", + " ax_fock.set_title('Fock State Distribution')\n", + " ax_fock.set_xlim(-0.5, min(dim - 0.5, dim - 0.5))\n", + " ax_fock.xaxis.set_major_locator(MaxNLocator(integer=True))\n", + " ax_fock.yaxis.grid(True, ls='--', alpha=0.4)\n", + " mean_n = np.sum(n_vals * fock_probs)\n", + " ax_fock.axvline(mean_n, color='crimson', lw=1.5, ls='--', label=f'⟨n⟩ = {mean_n:.2f}')\n", + " ax_fock.legend(fontsize=8)\n", + "\n", + " # Trim trailing zeros for display\n", + " last_nz = np.where(fock_probs > 1e-4)[0][-1]\n", + " ax_fock.set_xlim(-0.5, last_nz + 1.5)\n", + "\n", + " # ── 4-D Re(rho) heatmap ──────────────────────────────────────────────────────\n", + " ax_re = fig.add_subplot(gs[1, 0])\n", + " rho_max = np.max(np.abs(rho)) * 1.05\n", + " re_im = ax_re.imshow(np.real(rho),\n", + " cmap=CMAP_RHO,\n", + " vmin=-rho_max, vmax=rho_max,\n", + " interpolation='nearest', aspect='auto')\n", + " ax_re.set_xlabel('n')\n", + " ax_re.set_ylabel('m')\n", + " ax_re.set_title(r'Density Matrix — $\\mathrm{Re}(\\rho_{mn})$')\n", + " cb_re = fig.colorbar(re_im, ax=ax_re, fraction=0.046, pad=0.04)\n", + " cb_re.set_label(r'$\\mathrm{Re}(\\rho)$')\n", + " # Tick only every few steps for large matrices\n", + " tick_step = max(1, dim // 10)\n", + " ticks = np.arange(0, dim, tick_step)\n", + " ax_re.set_xticks(ticks); ax_re.set_yticks(ticks)\n", + "\n", + " # ── 4-E Im(rho) heatmap ──────────────────────────────────────────────────────\n", + " ax_im = fig.add_subplot(gs[1, 1])\n", + " im_im = ax_im.imshow(np.imag(rho),\n", + " cmap=CMAP_RHO,\n", + " vmin=-rho_max, vmax=rho_max,\n", + " interpolation='nearest', aspect='auto')\n", + " ax_im.set_xlabel('n')\n", + " ax_im.set_ylabel('m')\n", + " ax_im.set_title(r'Density Matrix — $\\mathrm{Im}(\\rho_{mn})$')\n", + " cb_im = fig.colorbar(im_im, ax=ax_im, fraction=0.046, pad=0.04)\n", + " cb_im.set_label(r'$\\mathrm{Im}(\\rho)$')\n", + " ax_im.set_xticks(ticks); ax_im.set_yticks(ticks)\n", + "\n", + " # ── 4-F Marginal distributions ───────────────────────────────────────────────\n", + " ax_marg = fig.add_subplot(gs[1, 2])\n", + "\n", + " # Marginals from Wigner function: integrate over one axis\n", + " # Renormalise so each integrates to exactly 1 (finite grid truncation\n", + " # causes small deviations from unity; shape is preserved exactly).\n", + " marg_x = np.sum(W, axis=0) * dx # ∫W dp → P(x)\n", + " marg_p = np.sum(W, axis=1) * dx # ∫W dx → P(p)\n", + " marg_x /= np.trapezoid(marg_x, xvec) # normalise: ∫P(x) dx = 1\n", + " marg_p /= np.trapezoid(marg_p, pvec) # normalise: ∫P(p) dp = 1\n", + "\n", + " ax_marg.plot(xvec, marg_x, color='steelblue',\n", + " lw=2, label=r'$P(x) = \\int W\\,dp$')\n", + " ax_marg.plot(pvec, marg_p, color='darkorange',\n", + " lw=2, ls='--', label=r'$P(p) = \\int W\\,dx$')\n", + " ax_marg.fill_between(xvec, marg_x, alpha=0.15, color='steelblue')\n", + " ax_marg.fill_between(pvec, marg_p, alpha=0.15, color='darkorange')\n", + " ax_marg.axhline(0, color='k', lw=0.5, ls='--', alpha=0.4)\n", + " ax_marg.set_xlabel('Quadrature value')\n", + " ax_marg.set_ylabel('Probability density')\n", + " ax_marg.set_title('Marginal Distributions')\n", + " ax_marg.legend(fontsize=8)\n", + " ax_marg.yaxis.grid(True, ls='--', alpha=0.4)\n", + "\n", + " # Annotate variances\n", + " var_x = np.trapezoid(xvec**2 * marg_x, xvec) - np.trapezoid(xvec * marg_x, xvec)**2\n", + " var_p = np.trapezoid(pvec**2 * marg_p, pvec) - np.trapezoid(pvec * marg_p, pvec)**2\n", + " heisenberg = np.sqrt(var_x * var_p)\n", + " ax_marg.text(0.97, 0.97,\n", + " f'$\\\\Delta x = {np.sqrt(var_x):.3f}$\\n'\n", + " f'$\\\\Delta p = {np.sqrt(var_p):.3f}$\\n'\n", + " f'$\\\\Delta x \\\\cdot \\\\Delta p = {heisenberg:.3f}$',\n", + " transform=ax_marg.transAxes, fontsize=8,\n", + " va='top', ha='right',\n", + " bbox=dict(boxstyle='round,pad=0.3', fc='white', alpha=0.8))\n", + "\n", + " plt.savefig('quantum_state_visualization.png',\n", + " dpi=150, bbox_inches='tight')\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "9dc2458f", + "metadata": {}, + "source": [ + "## 2 Interactive explorer\n", + "\n", + "The first run renders all 178 frames and pickles them to `wigner_frames_pijus.pkl` (~1 min on 8 cores, ~3 min elsewhere, 44 MB). Every run after that loads the file and skips the build. Once it is loaded the dropdown and slider are pure lookups, so you can drag continuously. Delete the file to re-render from scratch." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "77dfbbb0", + "metadata": {}, + "outputs": [], + "source": [ + "def qubit_rho(c0, c1, dim=10):\n", + " \"\"\"Pure state living in the lowest two Fock levels.\"\"\"\n", + " psi = np.zeros(dim, complex)\n", + " psi[0], psi[1] = c0, c1\n", + " return np.outer(psi, psi.conj()) / np.vdot(psi, psi).real\n", + " \n", + "STATES = {\n", + " 'coherent |α⟩': (lambda a: coherent_state_rho(a, 30), ('α', 0.0, 3.0, 0.1, 2.0)),\n", + " 'Fock |n⟩': (lambda n: fock_state_rho(int(n), 20), ('n', 0.0, 8.0, 1.0, 3.0)),\n", + " 'even cat |α⟩+|−α⟩': (lambda a: cat_state_rho(a, 40), ('α', 0.5, 3.0, 0.1, 2.0)),\n", + " 'thermal': (lambda nb: thermal_state_rho(nb, 30), ('n̄', 0.0, 5.0, 0.1, 2.0)),\n", + " 'squeezed vacuum': (lambda r: squeezed_vacuum_rho(r, 40), ('r', 0.0, 1.5, 0.05, 1.0)),\n", + " 'cat mixture ½(|α⟩⟨α|+|−α⟩⟨−α|)': (lambda a: 0.5 * coherent_state_rho(a, 40)\n", + " + 0.5 * coherent_state_rho(-a, 40),\n", + " ('α', 0.5, 3.0, 0.1, 2.0)),\n", + " '|+⟩ = (|0⟩+|1⟩)/√2': (lambda _: qubit_rho(r2, r2), None),\n", + " '|−⟩ = (|0⟩−|1⟩)/√2': (lambda _: qubit_rho(r2, -r2), None),\n", + " '|+i⟩ = (|0⟩+i|1⟩)/√2': (lambda _: qubit_rho(r2, 1j * r2), None),\n", + " 'mixed ½(|0⟩⟨0|+|1⟩⟨1|)': (lambda _: np.diag([0.5, 0.5] + [0.0] * 8).astype(complex), None),\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "1e837267", + "metadata": {}, + "outputs": [], + "source": [ + "r2 = 1 / np.sqrt(2)\n", + "\n", + "CACHE = Path('./wigner_frames_pijus.pkl')\n", + "DPI = 90 # 18×11 in at 90 dpi ≈ 250 kB per frame ≈ 44 MB for the full set\n", + "\n", + "_cache = pickle.loads(CACHE.read_bytes()) if CACHE.exists() else {}\n", + "\n", + "def _render(label, val):\n", + " build, spec = STATES[label]\n", + " rho = build(val)\n", + " title = label + (f' {spec[0]} = {val:g}' if spec else '')\n", + "\n", + " n_bar = float(np.real(np.diag(rho) @ np.arange(len(rho))))\n", + " lim = max(XLIM, 2.5 * np.sqrt(2 * (n_bar + 1)))\n", + " xs = np.linspace(-lim, lim, GRID)\n", + " dx = xs[1] - xs[0]\n", + " W = wigner_function(rho, xs, xs)\n", + "\n", + " old_show, old_save = plt.show, plt.savefig\n", + " plt.show = plt.savefig = lambda *a, **k: None\n", + " try:\n", + " visualize(rho, title, W, rho.shape[0], dx, grid=(xs, xs))\n", + " fig = plt.gcf()\n", + " buf = io.BytesIO()\n", + " fig.savefig(buf, format='png', dpi=DPI, bbox_inches='tight')\n", + " plt.close(fig)\n", + " finally:\n", + " plt.show, plt.savefig = old_show, old_save\n", + "\n", + " stats = (f'dim {rho.shape[0]} · Tr(ρ²) = {np.trace(rho @ rho).real:.3f} · '\n", + " f'∫W dx dp = {W.sum() * dx**2:.3f} · '\n", + " f'W ∈ [{W.min():.3f}, {W.max():.3f}]'\n", + " + (' · negative ⇒ non-classical' if W.min() < -1e-6 else ''))\n", + " return buf.getvalue(), stats\n", + "\n", + "def _render_key(key):\n", + " return _render(*key)\n", + "\n", + "def _combos():\n", + " \"\"\"Every (state, slider value) pair the widget can be asked to show.\"\"\"\n", + " out = []\n", + " for label, (_, spec) in STATES.items():\n", + " vals = np.arange(spec[1], spec[2] + spec[3] / 2, spec[3]) if spec else [0.0]\n", + " out += [(label, round(float(v), 6)) for v in vals]\n", + " return out\n", + "\n", + "def build_cache(force=False):\n", + " \"\"\"Render whatever is missing and save the lot. ~1 min on 8 cores, then never again.\"\"\"\n", + " todo = [k for k in _combos() if force or k not in _cache]\n", + " if not todo:\n", + " return\n", + " bar = w.IntProgress(min=0, max=len(todo), description=f'0/{len(todo)}',\n", + " style={'description_width': 'initial'})\n", + " display(bar)\n", + " \n", + " pool = (cf.ProcessPoolExecutor() if multiprocessing.get_start_method() == 'fork'\n", + " else None)\n", + " try:\n", + " rendered = pool.map(_render_key, todo, chunksize=2) if pool else map(_render_key, todo)\n", + " for key, frame in zip(todo, rendered):\n", + " _cache[key] = frame\n", + " bar.value += 1\n", + " bar.description = f'{bar.value}/{len(todo)}'\n", + " finally:\n", + " if pool:\n", + " pool.shutdown()\n", + " CACHE.write_bytes(pickle.dumps(_cache))\n", + " bar.bar_style = 'success'\n", + " bar.description = f'{len(_cache)} frames · {CACHE.stat().st_size / 1e6:.0f} MB'\n", + "\n", + "build_cache()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "f1b7f516", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "3a28bbad8e4a4b4b8944ff5c96bd6d7e", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "VBox(children=(HBox(children=(Dropdown(description='state:', layout=Layout(width='360px'), options=('coherent …" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# ── widgets ───────────────────────────────────────────────────────────────────\n", + "dd = w.Dropdown(options=list(STATES), description='state:', layout=w.Layout(width='360px'))\n", + "sl = w.FloatSlider(description='α', value=2.0, min=0.0, max=3.0, step=0.1,\n", + " continuous_update=True, readout_format='.2f')\n", + "img = w.Image(format='png', layout=w.Layout(width='100%'))\n", + "info = w.HTML()\n", + "\n", + "def _key():\n", + " return (dd.value, round(float(sl.value), 6) if STATES[dd.value][1] else 0.0)\n", + "\n", + "def draw(*_):\n", + " key = _key()\n", + " if key not in _cache: # slider landed off the prebuilt grid\n", + " _cache[key] = _render(*key)\n", + " img.value, info.value = _cache[key]\n", + "\n", + "def on_state(*_):\n", + " spec = STATES[dd.value][1]\n", + " sl.layout.display = 'none' if spec is None else ''\n", + " if spec:\n", + " with sl.hold_trait_notifications(): # set all bounds before revalidating\n", + " sl.description, sl.min, sl.max, sl.step, sl.value = spec\n", + " draw()\n", + "\n", + "dd.observe(on_state, names='value')\n", + "sl.observe(draw, names='value')\n", + "\n", + "display(w.VBox([w.HBox([dd, sl]), info, img]))\n", + "on_state()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/exercises/cv_submissions/wigner_frames_pijus.pkl b/exercises/cv_submissions/wigner_frames_pijus.pkl new file mode 100644 index 0000000..602dbfa Binary files /dev/null and b/exercises/cv_submissions/wigner_frames_pijus.pkl differ