{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a1b2c3d4",
   "metadata": {},
   "source": [
    "# Synthesize diagonal unitaries with Walsh coefficients\n",
    "\n",
    "**Download Notebook** - {nb-download}`diagonal_walsh.ipynb`\n",
    "\n",
    "Turn a list of diagonal phases into commuting Pauli-Z rotations:\n",
    "\n",
    "$$\n",
    "D=\\operatorname{diag}(e^{i\\theta_0},\\ldots,e^{i\\theta_{2^n-1}})\n",
    "=\\prod_{k=0}^{2^n-1}e^{ia_kW_k},\n",
    "\\qquad W_k=\\bigotimes_{r=0}^{n-1}Z_r^{k_r}.\n",
    "$$\n",
    "\n",
    "The normalized Walsh coefficients are\n",
    "\n",
    "$$\n",
    "a_k=2^{-n}\\sum_{j=0}^{2^n-1}(-1)^{\\sum_r j_rk_r}\\theta_j.\n",
    "$$\n",
    "\n",
    "`fast_walsh_hadamard_transform` computes the unnormalized transform; divide by $2^n$ to obtain $a_k$. `diagonal_unitary_walsh` builds the rotations and omits the $k=0$ global phase. Each remaining term uses a parity ladder, a Z rotation, and the inverse ladder.\n",
    "\n",
    "Run this notebook from a source checkout with the development dependencies installed. The repository default is little endian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b2c3d4e5",
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from guppyalgos.algorithms.state_preparation.diagonal import diagonal_unitary_walsh, fast_walsh_hadamard_transform\n",
    "from guppyalgos.tests.helpers import Endianness, assert_allclose_ignorephase, get_unitary"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3d4e5f6",
   "metadata": {},
   "source": [
    "## 1. One qubit: the S gate\n",
    "\n",
    "$$\n",
    "S=\\begin{pmatrix}1&0\\\\0&i\\end{pmatrix}\n",
    "=e^{i\\pi/4}e^{-i\\pi Z/4}.\n",
    "$$\n",
    "\n",
    "The normalized coefficients are $(\\pi/4,-\\pi/4)$. After removing the global phase, one Z rotation remains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d4e5f6a7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Phases:          [0.         1.57079633]\n",
      "Walsh coeffs:    [ 1.57079633 -1.57079633]\n",
      "Normalized a_k:  [ 0.78539816 -0.78539816]\n"
     ]
    }
   ],
   "source": [
    "diagonal_1q = np.array([1.0, 1j])\n",
    "phases_1q = np.angle(diagonal_1q)\n",
    "walsh_coeffs_1q = fast_walsh_hadamard_transform(phases_1q)\n",
    "print(\"Phases:         \", phases_1q)\n",
    "print(\"Walsh coeffs:   \", walsh_coeffs_1q)\n",
    "print(\"Normalized a_k: \", walsh_coeffs_1q / len(diagonal_1q))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e5f6a7b8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Synthesized unitary:\n",
      " [[7.07106781e-01-7.07106781e-01j 4.32978028e-17+4.32978028e-17j]\n",
      " [0.00000000e+00+0.00000000e+00j 7.07106781e-01+7.07106781e-01j]]\n",
      "Matches S gate up to global phase.\n"
     ]
    }
   ],
   "source": [
    "circ_1q = diagonal_unitary_walsh(diagonal_1q)\n",
    "unitary_1q = get_unitary(circ_1q, 1, endianness=Endianness.LITTLE)\n",
    "print(\"Synthesized unitary:\\n\", unitary_1q)\n",
    "assert_allclose_ignorephase(np.diag(diagonal_1q), unitary_1q)\n",
    "print(\"Matches S gate up to global phase.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6a7b8c9",
   "metadata": {},
   "source": [
    "## 2. Two qubits: the CZ gate\n",
    "\n",
    "$$\n",
    "\\mathrm{CZ}=\\operatorname{diag}(1,1,1,-1)\n",
    "=e^{i\\pi/4}e^{-i\\pi Z_0/4}e^{-i\\pi Z_1/4}e^{i\\pi Z_0Z_1/4}.\n",
    "$$\n",
    "\n",
    "The three non-global Walsh terms produce two single-qubit rotations and one two-qubit parity rotation. The check compares the synthesized matrix with CZ up to global phase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a7b8c9d0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Walsh coeffs:  [ 3.14159265 -3.14159265 -3.14159265  3.14159265]\n",
      "Non-trivial terms (excluding global phase): 3\n"
     ]
    }
   ],
   "source": [
    "diagonal_2q = np.array([1.0, 1.0, 1.0, -1.0])\n",
    "phases_2q = np.angle(diagonal_2q)\n",
    "walsh_coeffs_2q = fast_walsh_hadamard_transform(phases_2q)\n",
    "n_terms = np.sum(np.abs(walsh_coeffs_2q / len(diagonal_2q)) > 1e-10)\n",
    "print(\"Walsh coeffs: \", walsh_coeffs_2q)\n",
    "print(f\"Non-trivial terms (excluding global phase): {n_terms - 1}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b8c9d0e1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Synthesized unitary:\n",
      " [[-0.707107+0.707107j -0.      -0.j        0.      -0.j\n",
      "  -0.      -0.j      ]\n",
      " [ 0.      +0.j       -0.707107+0.707107j  0.      +0.j\n",
      "   0.      +0.j      ]\n",
      " [ 0.      +0.j       -0.      -0.j       -0.707107+0.707107j\n",
      "  -0.      -0.j      ]\n",
      " [ 0.      +0.j        0.      -0.j        0.      +0.j\n",
      "   0.707107-0.707107j]]\n",
      "Matches CZ gate up to global phase.\n"
     ]
    }
   ],
   "source": [
    "circ_2q = diagonal_unitary_walsh(diagonal_2q)\n",
    "unitary_2q = get_unitary(circ_2q, 2, endianness=Endianness.LITTLE)\n",
    "print(\"Synthesized unitary:\\n\", np.round(unitary_2q, 6))\n",
    "assert_allclose_ignorephase(np.diag(diagonal_2q), unitary_2q)\n",
    "print(\"Matches CZ gate up to global phase.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9d0e1f2",
   "metadata": {},
   "source": [
    "## 3. Trade accuracy for fewer rotations\n",
    "\n",
    "`truncation_threshold` drops small Walsh coefficients. Because the terms commute, the omitted coefficients bound the phase error:\n",
    "\n",
    "$$\n",
    "\\|D-\\widetilde D\\|_2\\leq\\sum_{k\\ \\mathrm{omitted}}|a_k|,\n",
    "$$\n",
    "\n",
    "using the same global-phase convention. The example below compares exact and truncated synthesis after aligning their global phases, and reports the Frobenius matrix error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d0e1f2a3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Frobenius error from truncation: 0.5445\n"
     ]
    }
   ],
   "source": [
    "rng = np.random.default_rng(42)\n",
    "diagonal_generic = np.exp(1j * rng.uniform(-np.pi, np.pi, 4))\n",
    "\n",
    "circ_exact = diagonal_unitary_walsh(diagonal_generic, truncation_threshold=0.0)\n",
    "circ_approx = diagonal_unitary_walsh(diagonal_generic, truncation_threshold=0.3)\n",
    "\n",
    "u_exact = get_unitary(circ_exact, 2, endianness=Endianness.LITTLE)\n",
    "u_approx = get_unitary(circ_approx, 2, endianness=Endianness.LITTLE)\n",
    "\n",
    "phase = np.angle(np.vdot(u_exact, u_approx))\n",
    "u_approx_aligned = u_approx * np.exp(-1j * phase)\n",
    "error = np.linalg.norm(u_exact - u_approx_aligned)\n",
    "print(f\"Frobenius error from truncation: {error:.4f}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "guppyalgos (3.13.5)",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}