{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "varied-fossil",
   "metadata": {},
   "source": [
    "# Circular Convolutional Layers in the Fourier Domain\n",
    "--------------------------------------------------------\n",
    "This is a tutorial accompanying the ICLR 2021 paper *\"Orthogonalizing Convolutional Layers with the Cayley Transform\"* by Asher Trockman and Zico Kolter.\n",
    "\n",
    "This Jupyter notebook is best viewed with [nbviewer (click here)](https://nbviewer.jupyter.org/github/locuslab/orthogonal-convolutions/blob/main/FFT%20Convolutions.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "military-interstate",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [3, 3]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "institutional-graduation",
   "metadata": {},
   "source": [
    "### Quick start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "balanced-narrative",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try changing these parameters and running the notebook.\n",
    "cin = 3      # input channels\n",
    "cout = 3     # output channels\n",
    "n = 5        # spatial (input) size\n",
    "k = 3        # conv. kernel size\n",
    "batches = 1  # batches"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "loose-aside",
   "metadata": {},
   "source": [
    "We will implement the following circular convolution using FFT functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "worst-boulder",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = torch.randn(batches, cin, n, n)\n",
    "conv = nn.Conv2d(cin, cout, k, bias=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "positive-finance",
   "metadata": {},
   "source": [
    "Notice that we make the convolution \"circular\" using circular padding (below). The left/right and top/bottom paddings differ to account for even kernel sizes (they will all be the same for odd kernel sizes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "located-thing",
   "metadata": {},
   "outputs": [],
   "source": [
    "y1 = conv(F.pad(x, ((k - 1) // 2, k // 2, (k - 1) // 2, k // 2), mode=\"circular\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "favorite-palestine",
   "metadata": {},
   "source": [
    "First, we compute the 2D FFT of the input (applied only to the last two dimensions):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "lyric-nicaragua",
   "metadata": {},
   "outputs": [],
   "source": [
    "xfft = torch.fft.fft2(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "still-return",
   "metadata": {},
   "source": [
    "Then we compute the 2D FFT of the weights, also applied to only the last two dimensions.\n",
    "We have to do a couple things here, however: first, we have to pad the weights to have the\n",
    "same spatial size as the inputs.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "checked-membrane",
   "metadata": {},
   "outputs": [],
   "source": [
    "wpad = F.pad(conv.weight, (0, n - k, 0, n - k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "physical-killer",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAL4AAADSCAYAAAD0Qnq8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAKu0lEQVR4nO3dbYxU5RnG8f/FCgFBCyJRYRFMtLSERGwRTbCt8SVFq9X4xZdK/GDSpNEGW5PGpo3V1vZTozbR1phqqfGtVq01RmtpizUmviGiEdEEjVSQCgIbAYm83f0wZ+3pdoWDO+fMLPf1SyY5M8+c57mZvebhnJnd5ygiMMtmRKcLMOsEB99ScvAtJQffUnLwLSUH31Jy8NtA0nWS7t5L+zuSzviMfX/mfe3TpQ5+EartkrZKel/SIknjOl1XVUW9NzQwzqmS1tQ9TpNSB79wbkSMA74EzAF+3OF6rAEOfiEi1gJPALMkTZD0mKQNkjYX2739z5V0jKR/StoiaTFweLkvSQskrZa0UdKPBrSNkHSNpLeK9gckHVZl372RNF1SSLpM0r8kfVDevzgce1DSH4q6l0k6vtQeko4t3V8k6QZJY4vXZXLxP+NWSZOr1tWtHPyCpKnA2cDLtF6X3wHTgKOB7cAtpaffC7xEK/A/Ay4r9TMT+A2wAJgMTAR6S/t+Fzgf+FrRvhm4teK+VZwCzABOB66V9MVS23nAH4HDin/DI5JG7q2ziNgGnAW8FxHjitt7+1lT94mItDfgHWAr0AesBn4NjBnkebOBzcX20cAuYGyp/V7g7mL7WuD+UttYYAdwRnF/JXB6qf0oYCdw0L72HaSuRcANxfZ0IIDeUvsLwEXF9nXAc6W2EcA64CvF/QCO/ZS+TwXWdPrn1c7bQUN4zxwozo+Iv5UfkHQwcBMwH5hQPHyIpB6KWTpaM2G/1cDUYnsy8G5/Q0Rsk7Sx9NxpwJ8k7Sk9ths4osK+Vfy7tP0RUD5ZL/e9pzhhHfaHLZ+FD3UGdzWtw4WTIuJQ4KvF46I1S04ojn37HV3aXsd/3wT9b6KJpfZ3gbMiYnzpNjpa5xj72neoyn2PoHUY1X/Y8hFwcOm5R5a2D7hf4XXwB3cIreP6vuLE8yf9DRGxGlgKXC9plKRTgHNL+z4InCPpFEmjgJ/yv6/zbcDPJU0DkDRJ0nkV9x2qL0u6QNJBwFXAx8BzRdty4BJJPZLm0zoH6fc+MFHS59pYS0c5+IO7GRgDfEArGH8Z0H4JcBKwidab4q7+hohYAVxB67h/Ha2T1/Jn4L8CHgX+KmlL0f9JFfcdqj8DFxb9LgAuiIidRdtCWm/gPuBbwCOlf9MbwH3A25L6DoRPdVScvNgBTtJ1tE5eL+10Ld3AM76l5OBbSj7UsZQ841tKDr6lVMs3t6PGj4kxRx5aR9eV7d4wqqPjH9f7fkfHB1i58YhOl9BxOzdvYve2bRr4eC3BH3Pkocy7/cI6uq5s023TOjr+E7+8qaPjA5x41/c6XULHrbll8J+DD3UsJQffUnLwLSUH31Jy8C0lB99ScvAtJQffUnLwLSUH31Jy8C2lSsGXNF/Sm5JWSbqm7qLM6rbP4BdrydxKazWtmcDFxYpfZsNWlRl/LrAqIt6OiB3A/bSWojMbtqoEfwqlFbhoLXcxZeCTJH1b0lJJS3f0bW9XfWa1aNvJbUTcHhFzImLOqPFj2tWtWS2qBH8tpaXnaC07t7aecsyaUSX4LwLHFWvCjwIuorUSmNmwtc8/PYyIXZKuBJ4EeoA7i6XuzIatSn9zGxGPA4/XXItZY/zNraXk4FtKDr6l5OBbSg6+peTgW0oOvqXk4FtKDr6l5OBbSrUsE74nxPZdI+vourIPp3f2PX3iXd/v6Pi2d57xLSUH31Jy8C0lB99ScvAtJQffUnLwLSUH31Jy8C0lB99ScvAtJQffUqqyTPidktZLeq2JgsyaUGXGXwTMr7kOs0btM/gR8TSwqYFazBrjY3xLqW3BL18YYmffR+3q1qwWtVwYYuT4g9vVrVktfKhjKVX5OPM+4FlghqQ1ki6vvyyzelW5MMTFTRRi1iQf6lhKDr6l5OBbSg6+peTgW0oOvqXk4FtKDr6l5OBbSg6+peTgW0q1XBhixug+/j7z0Tq6ruzrZ8zu6PgjZn2ho+MDvHXJhE6X0LU841tKDr6l5OBbSg6+peTgW0oOvqXk4FtKDr6l5OBbSg6+peTgW0oOvqVUZSW1qZKWSHpd0gpJC5sozKxOVX47cxdwdUQsk3QI8JKkxRHxes21mdWmyoUh1kXEsmJ7C7ASmFJ3YWZ12q9jfEnTgROA5wdp+2R9/A0bd7epPLN6VA6+pHHAQ8BVEfHhwPby+viTJva0s0aztqsUfEkjaYX+noh4uN6SzOpX5VMdAXcAKyPixvpLMqtflRl/HrAAOE3S8uJ2ds11mdWqyoUhngHUQC1mjfE3t5aSg28pOfiWkoNvKTn4lpKDbyk5+JaSg28pOfiWkoNvKdWyPv5rGyfx+UXfqaPr6n7R2eGtu3nGt5QcfEvJwbeUHHxLycG3lBx8S8nBt5QcfEvJwbeUHHxLycG3lBx8S6nKSmqjJb0g6ZViffzrmyjMrE5VfjvzY+C0iNharKH5jKQnIuK5mmszq02VldQC2FrcHVncos6izOpWdbXkHknLgfXA4oj4v/XxzYaTSsGPiN0RMRvoBeZKmjXwOeULQ+zetq3NZZq11359qhMRfcASYP4gbZ9cGKJn7Ng2lWdWjyqf6kySNL7YHgOcCbxRc11mtaryqc5RwO8l9dB6ozwQEY/VW5ZZvap8qvMqrQu+mR0w/M2tpeTgW0oOvqXk4FtKDr6l5OBbSg6+peTgW0oOvqXk4FtKDr6l5OBbSg6+peTgW0oOvqXk4FtKDr6l5OBbSg6+peTgW0oOvqXk4FtKDr6l5OBbSpWDX6yY/LIkr6Jmw97+zPgLgZV1FWLWpKrr4/cC3wB+W285Zs2oOuPfDPwA2PNpT/D6+DacVFkm/BxgfUS8tLfneX18G06qzPjzgG9Kege4HzhN0t21VmVWs30GPyJ+GBG9ETEduAj4R0RcWntlZjXy5/iWUpUronwiIp4CnqqlErMGeca3lBx8S8nBt5QcfEvJwbeUHHxLycG3lBx8S8nBt5QcfEvJwbeUFBHt71TaAKweQheHAx+0qZzhOL5raN/40yJi0sAHawn+UElaGhFzso7vGuof34c6lpKDbyl1a/BvTz4+uIZax+/KY3yzunXrjG9Wq64KvqT5kt6UtErSNR0Y/05J6yW91vTYpRqmSloi6XVJKyQtbHj80ZJekPRKMf71TY4/oJbalq3smuBL6gFuBc4CZgIXS5rZcBmLgPkNjznQLuDqiJgJnAxc0fDr8DFwWkQcD8wG5ks6ucHxy2pbtrJrgg/MBVZFxNsRsYPWGj7nNVlARDwNbGpyzEFqWBcRy4rtLbR+8FMaHD8iYmtxd2Rxa/xEsO5lK7sp+FOAd0v319DgD7wbSZoOnAA83/C4PZKWA+uBxRHR6PiFm9nHspVD0U3BtxJJ44CHgKsi4sMmx46I3RExG+gF5kqa1eT4VZetHIpuCv5aYGrpfm/xWDqSRtIK/T0R8XCn6oiIPmAJzZ/31L5sZTcF/0XgOEnHSBpFa7nCRztcU+MkCbgDWBkRN3Zg/EmSxhfbY4AzgTearKGJZSu7JvgRsQu4EniS1gndAxGxoskaJN0HPAvMkLRG0uVNjl+YByygNcstL25nNzj+UcASSa/SmowWR8QBdxUcf3NrKXXNjG/WJAffUnLwLSUH31Jy8C0lB99ScvAtJQffUvoPP3aC0AYsZjoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 216x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title('Padded Input')\n",
    "plt.imshow(wpad.detach()[0, 0]); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cooperative-cooper",
   "metadata": {},
   "source": [
    "Next, in order to agree with the circular convolution implemented in PyTorch, we have to shift the kernel to the \"center\" of the input field, which is the top left corner. We could implement this using `torch.roll` operations, but we found that it is easier and more efficient to center the kernel\n",
    "in the Fourier domain using the [shift theorem](https://ccrma.stanford.edu/~jos/st/Shift_Theorem.html).\n",
    "\n",
    "We construct a \"shift matrix\" and then use it in an elementwise product in the Fourier domain. The *shift amount* ensures that the center of the kernel is in the top left of the input field; the remainder of the kernel wraps around as if the input were \"circular\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "permanent-possession",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fft_shift_matrix(n, shift_amount):\n",
    "    shift = torch.arange(0, n).repeat((n, 1))\n",
    "    shift = shift + shift.T\n",
    "    return torch.exp(1j * 2 * np.pi * shift_amount * shift / n)\n",
    "    \n",
    "shift_amount = (k - 1) // 2\n",
    "shift_matrix = fft_shift_matrix(n, -shift_amount)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "elementary-irrigation",
   "metadata": {},
   "source": [
    "Then, we also must take the complex conjugate. This is equivalent to flipping the kernel horizontally and vertically--\n",
    "convolution as implemented in PyTorch is actually \"cross-correlation\" mathematically.\n",
    "The difference is just that the kernel is flipped.\n",
    "Conjugation in the Fourier domain is equivalent to flipping the signal in the spatial domain.\n",
    "Refer to the [flip theorems](https://ccrma.stanford.edu/~jos/sasp/Flip_Theorems.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "thousand-steam",
   "metadata": {},
   "outputs": [],
   "source": [
    "wfft = shift_matrix * torch.fft.fft2(wpad).conj()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dying-there",
   "metadata": {},
   "source": [
    "The result of our Fourier-domain shifting and flipping operations can be seen below. Compare to the padded input above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "textile-security",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAN0AAADSCAYAAADOksXPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAPWUlEQVR4nO3de7BdZX3G8e+T5ITcuIMISSDItZERsMhlRGUi0EBR6EynBhHBSztaUmGgtIiVi4CjbbnYAYsImFpuImKlKEIsCRrljoEKwZlIExMIhJCkkMMtCb/+8b4nruyeyz6clXefvfN8ZvbM2vtd633fvdZ61u3ss5YiAjMrZ0SrO2C2uXHozApz6MwKc+jMCnPozApz6MwKKxY6SRdIuqGf8kWSjnybdb/tadudpJC0Z6vblnS1pC+3oA+nSppXut2hGDB0eYV+TdIaSS9ImiVpQonObSqSPi7pkfydlkm6S9LhNdTb74al1ZRcKuml/LqtiWnmSno9z6ue12GN40XE5yLiok3T87dH0pS8YRhVoK25kj7bzLjN7uk+EhETgPcCBwH/8HY712qSzgSuAL4K7ATsCnwTOL6F3QKgwMpxNPAJYH9gF+BbTU43MyImVF73b7IebgYGdXgZEc8CdwH7SdpW0p2SXpS0Kg9P6hlX0u6S7pP0iqTZwA7VuiSdLGlx3uJ+qaFshKRzJP0ul98qabtmpu2PpK2BrwCnRcTtEdEdEWsj4j8j4uyB2q5sOU+R9HtJK3ralzQdOBf4WN4bPN7TpqTr8h71WUkXSxqZy06V9EtJl0t6CbhA0haS/jnX/0I+bBtb+Q5n57qek/TpZr97thZ4DXg+It6IiNmDnL5P+Qjo4jx8hKSlks7N82iRpJMaxr1a0uy8ftwnabdK+b65bKWk30r6i0rZ9pLukPSypIeAPQbZx6sk/Ti3+6CkPSrlIekLkp7J/f4nSSNy2UZHMdW9qKRLgA8AV+Zlf2V//RhU6CRNBo4Ffp2n/Q6wG2lv8RpQbewm4FFS2C4CTqnUMxX4V+Bk0hZ3e2BSZdq/AU4APpTLVwFXNTltfw4DxgA/7GecPtuuOBzYB/gwcJ6kP4qIn5L2nt/Le4P987izgHXAnsCBpL1N9TDkEOAZ0l73EuBrwN7AAXmaicB5+btPB/4WOArYCxjseezTwHbAtT0r0yb0TtKyn0ha9tdI2qdSfhJpvdgBmA/cCCBpPDCbtP68A5gBfDMvd0jL4nVgZ+DT+TUYM4ALgW2BhaR5XvVnpKO595KOfgasPyK+BPyCPxwRzBxogn5fwCJgDbAaWEw6FBvby3gHAKvy8K6kFW18pfwm4IY8fB5wS6VsPPAmcGR+vwD4cKV8Z9JWetRA0w7wXU4ibeX7G6e/tqcAAUyqlD8EzMjDF/R8x/x+J+CN6vwCTgTm5OFTgd9XygR0A3tUPjsM+J88fD3wtUrZ3rk/ezbx3buA/yYdXv4o1zUil80jnUL0Nt1c4NW8/FcDj1XKNrRN2rhcnIeP6GX53wp8uTJudRlOANYDk4GPAb9o6MO3gPOBkXlZ7Fsp+yowr4++9yyvUZV2r62UHws83fB9plfe/zXwX30s28a65wKfHWg5RATNnkOcEBE/q34gaRxwOTCdtNUA2DIfOu1CCmB3ZZLFeaaSy5f0FEREdz686rEb8ENJb1U+W09aiQeatj8vATtIGhUR6/oYp7+2ezxfGX6VtNL0VVcXsExSz2cjqv1vGN4RGAc8WhlfpJUN0nd/tDL+4j7a7c00YHRE3KB0AeUu0h7vDGBfUvD68oWIuHYQbUHvy3+XyvvqMlwjaWUu3w04RNLqyrijgH8nzZ9RbDzPBjMPYOBl11j3LtRsKIcYZ5EOsQ6JiK2AD+bPBSwDts2HCj12rQwv4w8B7Anw9pXyJcAxEbFN5TUm0jnlQNP2537SnueEfsbpr+2BNP7LxpLc3g6VuraKiHf3Mc0K0mH6uyvjbx3pIhY0fHc2nqcDGUXaABARrwMfBd4DPEza66waRF3N6G35P1d5X12GE0iHvc+R5tl9DfN/QkR8HniRtAd9u/OgGY119/S5m7RB7PHOhuma/nedoYRuS9IKsjpfaDh/Q+sRi4FHgAsljVa6HP+RyrS3AcdJOlzSaNLFjWpfrgYu6Tm5lrSjpOObmTafxPc6AyLif0mHp1dJOkHSOEldko6R9I9NtD2QF4ApPedLEbEMuAe4VNJWShdp9pD0oT769xbwbeBySe/I7U+U9Cd5lFuBUyVNzRub86vT5wszi/ro2zxgjKSv5AszI4A5pEPUV5v8foPVs/w/ABwHfL9SdmxlGV4EPBARS4A7gb2VLpZ15df78nnzeuB20gWncfk875TGRofobKWLhJOB04Hv5c/nAx+UtKvSBbkvNkz3AvCuZhoYSuiuAMaSts4PAD9tKP846SLBStLK8d2egoh4EjiNdJ63jHSxYmll2m8AdwD3SHol139Ik9NOBn7VV6cj4lLgTNKfPV4kbVlnAv8xUNtN6FmpXpL0WB7+JDAaeCr39TbSeWJf/p50gv+ApJeBn5GOKIiIu0jz/d48zr0N004GftlbpXmDczRwKGnr/TvSEcLBwKck/WWT37FZz5O+73OkiySfi4inK+U3kdaLlcAfk841iYhXcj9n5GmfB74ObJGnm0k6JHyedI72nZr7/SPSIfx84MfAdblfs0kBfCKX39kw3TeAP1e6kv8v/TWg6LB/YpV0LfD9iLi71X0pTdI9wOkRsaDF/TiCdNGh16vKkmYBSyNiWP29Nx8h7RURCzdlO5v8L/WlRURTvwroRBFxdKv7YAPzD57NCuu4w0uz4c57OrPCHDqzwlpyIWXUuPHRtfV2A4+4Kfuw1dqWtr+2u6ul7VuydtVK1nd3a+Ax69OS0HVtvR27f+rMVjS9wU5HLh14pE1oyUMTW9q+JUuvvLx4mz68NCvMoTMrzKEzK8yhMyvMoTMrzKEzK8yhMyvMoTMrzKEzK8yhMyvMoTMrrJbQSZqe78S7UNI5ddRp1qmGHLp8n8urgGOAqcCJlbvxmlmDOvZ0BwMLI+KZiHgTuIVh8DAOs+GqjtBNZOO74i7Nn21E0l8pPZ7qkXWvdjcWm202il1IiYhrIuKgiDho1LjxA09g1qHqCN2zbHwr6kn5MzPrRR2hexjYS+l5dKNJd+a9o4Z6zTrSkG/XEBHrJM0E7iY9Xeb6fOtzM+tFLfdIiYifAD+poy6zTudfpJgV5tCZFebQmRXm0JkV5tCZFebQmRXm0JkV5tCZFebQmRXm0JkV1pJHZY18E7Za9FYrmt5g7KjWPp/Okoc/eVlL25928/LibXpPZ1aYQ2dWmENnVphDZ1aYQ2dWmENnVphDZ1aYQ2dWmENnVphDZ1aYQ2dWmENnVlhdz6e7XtJySb+poz6zTlbXnm4WML2musw6Wi2hi4ifAyvrqMus0/mczqywYqGrPhRy7et+KKRtvlryUMiuMX4opG2+fHhpVlhdfzK4Gbgf2EfSUkmfqaNes05U1/PpTqyjHrPNgQ8vzQpz6MwKc+jMCnPozApz6MwKc+jMCnPozApz6MwKc+jMCnPozApz6MwKa8lDIdeNhRXvUSua3mDFr6a0tH1L3vfdM1va/tKXLi/epvd0ZoU5dGaFOXRmhTl0ZoU5dGaFOXRmhTl0ZoU5dGaFOXRmhTl0ZoU5dGaFOXRmhQ05dJImS5oj6SlJT0o6vY6OmXWqOv7LYB1wVkQ8JmlL4FFJsyPiqRrqNus4Q97TRcSyiHgsD78CLAAmDrVes05V6zmdpCnAgcCDvZRteD7d+m4/n842X7WFTtIE4AfAGRHxcmN59fl0I8f7+XS2+arrUVldpMDdGBG311GnWaeq4+qlgOuABRFx2dC7ZNbZ6tjTvR84GZgmaX5+HVtDvWYdach/MoiIeUBr7zJk1kb8ixSzwhw6s8IcOrPCHDqzwhw6s8IcOrPCHDqzwhw6s8IcOrPCHDqzwhw6s8IcOrPCHDqzwhw6s8IcOrPCHDqzwhw6s8IcOrPCHDqzwhw6s8IcOrPCHDqzwhw6s8Lquq36GEkPSXo8P6PuwjrqNetEdTyfDuANYFpErMnPNZgn6a6IeKCm+s06Ri2hi4gA1uS3XfkVddRt1mnqfFTWSEnzgeXA7Ij4f8+oM7MaQxcR6yPiAGAScLCk/arlfiikWVL71cuIWA3MAaY3fO6HQppR39XLHSVtk4fHAkcBT9dRt1mnqevq5c7Av0kaSQryrRFxZ011m3WUuq5ePgEcWEddZp3Ov0gxK8yhMyvMoTMrzKEzK8yhMyvMoTMrzKEzK8yhMyvMoTMrzKEzK8yhMyusrh88D8roZ7vZ/dz7W9H0Bnc/N7+l7e896/MtbX+42OOmVS1t/8WV64u36T2dWWEOnVlhDp1ZYQ6dWWEOnVlhDp1ZYQ6dWWEOnVlhDp1ZYQ6dWWEOnVlhDp1ZYXU/tefXknxnZ7N+1LmnOx1YUGN9Zh2prgeITAL+FLi2jvrMOllde7orgL8D3uprhOrz6dbyRk3NmrWfIYdO0nHA8oh4tL/xqs+n62KLoTZr1rbq2NO9H/iopEXALcA0STfUUK9ZRxpy6CLiixExKSKmADOAeyPiE0PumVmH8t/pzAqr9cZEETEXmFtnnWadxns6s8IcOrPCHDqzwhw6s8IcOrPCHDqzwhw6s8IcOrPCHDqzwhw6s8IcOrPCFBHlG5VeBBYPoYodgBU1dcd9aN/26+jDbhGxY12daUZLQjdUkh6JiIPch9b2odXtD5c+DJYPL80Kc+jMCmvX0F3T6g7gPgyH9mF49GFQ2vKczqydteuezqxttV3oJE2X9FtJCyWd04L2r5e0XNJvSred258saY6kpyQ9Ken0FvRhjKSHJD2e+3Bh6T7kfrTlrfzbKnSSRgJXAccAU4ETJU0t3I1ZwPTCbVatA86KiKnAocBpLZgHbwDTImJ/4ABguqRDC/cB2vRW/m0VOuBgYGFEPBMRb5Lus3l8yQ5ExM+BlSXbbGh/WUQ8lodfIa10Ewv3ISJiTX7blV9FLw6086382y10E4EllfdLKbzCDSeSpgAHAg+2oO2RkuYDy4HZEVG6D1cwwK38h6t2C51lkiYAPwDOiIiXS7cfEesj4gBgEnCwpP1Ktd3srfyHq3YL3bPA5Mr7SfmzzYqkLlLgboyI21vZl4hYDcyh7HluW9/Kv91C9zCwl6TdJY0m3cb9jhb3qShJAq4DFkTEZS3qw46StsnDY4GjgKdLtd/ut/Jvq9BFxDpgJnA36QLCrRHxZMk+SLoZuB/YR9JSSZ8p2T5pK38yaes+P7+OLdyHnYE5kp4gbQhnR0RbXbZvJf8ixaywttrTmXUCh86sMIfOrDCHzqwwh86sMIfOrDCHzqwwh86ssP8Dy/K9LKsBS0gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 216x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title('Padded, Centered, & Flipped Input')\n",
    "plt.imshow(torch.fft.ifft2(wfft)[0, 0].real.detach()); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "composed-equity",
   "metadata": {},
   "source": [
    "Then, we implement the product described by this figure, using `einsum`.<br><br>\n",
    "<center>\n",
    "<img width=\"200\" src=\"https://github.com/locuslab/orthogonal-convolutions/blob/main/img/AnimationKeyVerySmall.gif?raw=true\" alt=\"Multi-Channel 2D Convolution Theorem\" /><br>\n",
    "    <p style=\"width: 50%; min-width: 400px;\">(<b>Colored slices:</b> pixels; <b>arrows:</b> dot products; <b>left:</b> Fourier-domain weights; <b>middle:</b> input tensor; <b>right:</b> output tensor; <b>bottom-right</b>: the matrix-vector product for each pixel.)</p>\n",
    "</center><br>\n",
    "That is: Each channel of each pixel of the output is the dot product of the corresponding input pixel and weight pixel.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "induced-inspector",
   "metadata": {},
   "outputs": [],
   "source": [
    "yfft = torch.einsum('dchw, bchw -> bdhw', wfft, xfft)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "pursuant-services",
   "metadata": {},
   "source": [
    "Now we apply the inverse Fourier transform, and we're done."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "specialized-walker",
   "metadata": {},
   "outputs": [],
   "source": [
    "y2 = torch.fft.ifft2(yfft)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "banned-compiler",
   "metadata": {},
   "source": [
    "The result is the same as the circular convolution we implemented with `nn.Conv2d`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "catholic-maldives",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.5329876532632625e-06"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(y1 - y2).norm().item()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "heated-employee",
   "metadata": {},
   "source": [
    "------------------------------------------------------------------------\n",
    "### Matrix Block-Diagonalization Perspective\n",
    "And how to orthogonalize FFT convolutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "tribal-toner",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [5, 5]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "liable-joseph",
   "metadata": {},
   "source": [
    "We can make a block diagonal Fourier-domain matrix out of the $n^2$ pixels of `wfft`.\n",
    "Then, we can implement circular convolution in the Fourier domain by multiplying the Fourier-domain\n",
    "inputs with the block diagonal matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "minute-fifteen",
   "metadata": {},
   "outputs": [],
   "source": [
    "D = torch.block_diag(*[wfft.reshape(cin, cout, n**2)[..., i] for i in range(n**2)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "unusual-specific",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATEAAAFBCAYAAAAfR5ipAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAcOUlEQVR4nO3de5TcZZ3n8fcn1Z2k0yC5QUhIJLhkwzAjBCYgHJ0ZLuKAorgecWDcNeuyy+roDh51Bced3Zkd9yyeXa+zo57MoGZcV7moAzKjM2yE9SCeaNSIQIghmEhCksZcOx1y6e7v/vH7pap+RV+q0tVV/XR/XufUqef5Pb+q37cu/e3neep3UURgZpaqae0OwMxsLJzEzCxpTmJmljQnMTNLmpOYmSXNSczMkuYkZmZJcxIzs6Q5iU0ikrZKem27n6Pm+Z6UdEWznq/Zmv16x2Kiv1cTlZNYC+R/KCFpUNKhvH6PpFc1+Bwv5o/fJ+nvJS0Zz7gbiKlX0n5Jj0l6l6Ty9yoifjMiHmljmGOSv8ZjkubXLP9p/pkurfM5Rk2Uqb9X7eIk1lp/D9wDHAVuBB6VdGMDj39jRJwCLAR2A3/Z/BAb9saIOBU4G7gTuB24q70hNd0vgZtPVCS9EpjVrCeX1NGs55qKnMRa666I+DfAbwJfAzqAz0tq6A8iIo4A9wHnD7eOpN+Q9EjeQ3pS0puq2pZI+oakFyTtkfS/RniOX0q6eaj2mpgORMQDwB8AqyT9Vv4c5V6IpDskbcl7bk9J+hc127s47+H0SrpX0t2SPlrH69kq6YOSHpd0IH/czKr2Ebdbhy8D76iqrwL+tib2Ibch6cvAy4Fv5b3oD1XFfLukx4E+SR0n3itJ/0zSXkkX5+suyj+rKxqMe2qICN/G+QZsBQJ4c9Wy8/JlAVwDfBb47CjP8dq8PAtYA/ztUOsAncAzwJ8A04GrgF5gOVACfgZ8EugGZgKvGeI5LgZ+BVxfT0w1y38FvHuIuG8EFpH98/wDoA9YmLdNB7YBt+XxvwU4Bnx0pNdTtY0f5s89F9gIvKsqnpG2O+RrGOL92AT8Rv7+bSfreQaw9GS2kS/bACwBuoZ4r/4d8FT+Wf8j8D/b/T2eqDf3xNpnW1X5jIj4o4j4o1Ee83eS9gMHyBLf/xhmvcuAU4A7I+JYRHwXeJBsSHQp2R/bf4yIvog4EhGP1jz+d4AHgHdExIMNvarM82TJpCAi7o2I5yNiMCLuBjbn8ZyIuQP4TEQcj4hvkCWm0V7PCZ/Jn3sv8C1gRZ3brdeJ3tg1ZElyRwOvbTifiYjnIuLF2oaI+GuyxL2ObPrgIw3GO2U4ibXP2VXlnjof8+aImE3We3ov8P8knTnEeouA5yJisGrZNuAssv/82yKif4TtvAt4LE5+kvksYG/tQknvkLQhHxLuB34LODFhvgjYEXk3JPdcVdtwr+eEXVXlw2RJr57t1uvLwB8C/5qaoeQYtvHcKO1/nT/PX0bE0QbjnTKcxNogn8j9L3l1L/D9Rh4fEQN5T2UAeM0QqzwPLKn+lZBsXmYH2R/Oy0eZTH5Xvs4nG4kLQNIlZMnl0ZrlZ5P9Ub4XmJcn4ycA5avsBM6SpKqHnfj1daTXM1o8o223LhGxjWyC//XANxrcxnAn7Rv2ZH6STgE+RfYjyZ9JeknP1jJOYq11i6QvAE8CNwH9ZHM3hxt5EmVuAOaQDW1qrSPrjXxIUmc+IfxGsh8TfkiWMO6U1C1ppqRX1zy+F7gW+F1Jd9YZ08skXZ9v439HxM9rVukm+6N9IV//nWS9jBN+QJaU35tPct9AZTg20usZzWjbbcQtwFUR0dfgNnYDr2hwW58G1kfEvyX7VfvzJxXxFOAk1lpvIJv0nUG2q8WrI+JeAEmflzTaF/Vbkg4BB4H/BqyKiCdrV4qIY2R/5NcBvyb70eAdEfF0RAzkbeeSTcBvz2OqfY79ZPM/10n6i1Fi6iXr4X0E+ATwziGe7yng42TJajfwSqp6oHnMbyFLFPuBf0k273V0pNczQlx1bbcREbElItafxDb+O/Cf8qHmB0fbTp7ArwXenS96P3CxpLefTNyTnYpTEGYTh6R1wOcj4ovtjsUmLvfEbMKQ9HuSzsyHk6uAC4DvtDsum9i8p7BNJMvJhtndwLPAWyNiZ3tDsonOw0kzS9qYhpOSrpW0SdIzku5oVlBmZvU66Z6YpBLwC7JfsLYDPwJuzn+pMTNribHMiV0KPBMRzwJI+hpwA9nxXkMqndodHafPLtd1pNIR7DxUTKbHTmtoX0Qzm+SObd/+64g4vXb5WJLYWRQPm9gOjHh+rI7TZ7PoL95Trs/4RVe5fOYPikdV/Oq6zjGEZmaTzS/f/8FtQy0f910sJN0qab2k9QMHa3d0NjMbm7EksR1Ujm0DWMwQx7JFxOqIWBkRK0sv6x7D5szMXmosw8kfAcsknUOWvG4iO8p/eANisK9qmFg1DbZ/2fSalb3rh5mN7qSTWET0S3ov2QnbSsAXhjqOz8xsPI1pj/2I+AfgH5oUi5lZw1p62NG0zkFOPbO3XH/xwOxy+ejc4vBRHk2aWR18ALiZJc1JzMyS5iRmZklr6ZzY4NESh7ecVrX1ysRX7RxY58Fifh2sWrerp3hI0qGlg5jZ1OSemJklzUnMzJI2Yc/semRh8bKIXdsroR4628NHM8u4J2ZmSXMSM7OkOYmZWdIm7JxY595SoX7kzIFyOboGCm2l/RP2ZZjZOHNPzMyS5iRmZklzEjOzpE3YyaTBGcXjkKYdrTrU6Ggx7KUPHiuXt76heIbYaYsPF593x6wmRWhmE4F7YmaWNCcxM0vahB1ONmLLjZWX0bWzmJf/+SW7C/UndpzTkpjMrDXcEzOzpDmJmVnSnMTMLGmTYk5s2Vcqu1hseffxQtuO3tNqVzezScQ9MTNLmpOYmSXNSczMkjYp5sSefcvMclnF3cLYt3tGoe5DlMwml1F7YpK+IKlH0hNVy+ZKekjS5vx+zviGaWY2tHqGk18Crq1ZdgewNiKWAWvzuplZy406nIyI70laWrP4BuCKvLwGeAS4vZmBjRcfomQ2uZzsxP6CiNiZl3cBC5oUj5lZQ8b862REBBDDtUu6VdJ6SesH+vrGujkzs4KTTWK7JS0EyO97hlsxIlZHxMqIWFnq7j7JzZmZDe1kd7F4AFgF3Jnf39+0iMaZD1Eym1zq2cXiq8APgOWStku6hSx5XSNpM/DavG5m1nL1/Dp58zBNVzc5FjOzhvmwIzNL2qQ47KgRjRyi1HmwkuMHO4o/wHb1VK6+dGjpYBMjNLNGuCdmZklzEjOzpE254WQjjizsL5e7thffqkNnewhpNhG4J2ZmSXMSM7OkOYmZWdI8JzaCzr2lcvnImQOFtuiq1Ev7/TaatYt7YmaWNCcxM0uax0EjGJxR2Ut/2lEVG49W3rqBU4pDzVlbOwv1UuXEGRxe5F0zzJrJPTEzS5qTmJklzUnMzJLmObEm0JFSod7xqn2Fet/m2S2MxmxqcU/MzJLmJGZmSXMSM7OkeU6sCTrOeLFQ77qveNWk3ssr+4ZNO1Kzv5mZjYl7YmaWNCcxM0uah5NNMPh8V6H+64uK7dVDyDnL9xbb7plXqM/aXbmg73Ov88djNhr3xMwsaU5iZpY0JzEzS5onXVrs6HfnF+q9l/cX6jNeqD6NT/GCvWb2UqP2xCQtkfSwpKckPSnptnz5XEkPSdqc388Z/3DNzIrqGU72Ax+IiPOBy4D3SDofuANYGxHLgLV53cyspUZNYhGxMyJ+kpd7gY3AWcANwJp8tTXAm8cpRjOzYTU0JyZpKXARsA5YEBE786ZdwILmhjY5HTqneCrrjgPF0/gMepbSrCF1/zop6RTg68D7IuJgdVtEBMPMQku6VdJ6SesH+vrGFKyZWa26kpikTrIE9pWI+Ea+eLekhXn7QqBnqMdGxOqIWBkRK0vd3c2I2cysbNTBiyQBdwEbI+ITVU0PAKuAO/P7+8clwklm2rHiWSyis9iBra5poLhu166qes3JMHwVJZuq6pmBeTXwr4CfS9qQL/sTsuR1j6RbgG3A28YlQjOzEYyaxCLiUV7yf7/s6uaGY2bWGB92ZGZJ8w/6E1jpxWIH+PjLKuUjZxYPVyr1+f+RTU3+5ptZ0pzEzCxpHk5OYN07irtf9F55uFyetrOrdnWzKck9MTNLmpOYmSXNSczMkuY5sQnswPKaY+qrrqpUu/fxmY8V191xTaU+f/H+Qtu+TXObEZ7ZhOCemJklzUnMzJLmJGZmSfOc2CTRc+OLhXpp+6xy+Wh/qXZ1s0nDPTEzS5qTmJklzcPJSWL2t4un/t6zorKLxaGDxUOU/J/LJhN/n80saU5iZpY0JzEzS5rnxCaJvRcUDztS1cWP9ML0Qtu8DcWDlnquOF4uz19QuKSoD1GyCc89MTNLmpOYmSXNSczMkuY5sSno+I17C/UZRypzZgd6fdprS4t7YmaWNCcxM0uah5NTUO9Txd0m5jxdKe993YuYpcQ9MTNL2qhJTNJMST+U9DNJT0r683z5OZLWSXpG0t2Spo/2XGZmzVZPT+wocFVEXAisAK6VdBnwMeCTEXEusA+4ZdyiNDMbxqhzYhERwKG82pnfArgK+MN8+Rrgz4DPNT9Ea7boLB6itPeVVZWdMwtt6q8cotTVUzxc6dRrdhXqu588ozkBmjWgrjkxSSVJG4Ae4CFgC7A/IvrzVbYDZ41LhGZmI6griUXEQESsABYDlwLn1bsBSbdKWi9p/UBf38lFaWY2jIZ2sYiI/ZIeBi4HZkvqyHtji4EdwzxmNbAaYMaSJTHUOjZxTauc4ILHP/DZQtv1v7iuUN/dioDMatTz6+Tpkmbn5S7gGmAj8DDw1ny1VcD94xSjmdmw6umJLQTWSCqRJb17IuJBSU8BX5P0UeCnwF3jGKeZ2ZDq+XXyceCiIZY/SzY/ZmbWNj7syEZ07PT+cvn3F60otG3+0qJC3Yd/WDv4e2dmSXMSM7OkOYmZWdI8J2YjKh0qlctbPn5ZoW3anuK6c5ZXnTH26/MKbX2Li4csHZs9iFkzuCdmZklzEjOzpHk4aU1z5Hvzy+U3/vGjhbavrisORUt9/v9pzeFvkpklzUnMzJLmJGZmSfOcmDXNofOOlcvf/9PiHNiMS4pftf5TvIuFNYd7YmaWNCcxM0uak5iZJc1zYtY0pX2Vr9OOK2tbi3NgF1yypVx+9t5lhbYZ+ytnMd+zwmc0t5G5J2ZmSXMSM7OkeThpbfHEY+eWy/3nFoeap22q/t/q4aSNzD0xM0uak5iZJc1JzMyS5jkxa4vjZ1QuLa7DpUJb6ajnwax+7omZWdKcxMwsaR5OWltU791fa//5leFkLDhaaOvYNrNQj6qR6GCnh6FTkXtiZpa0upOYpJKkn0p6MK+fI2mdpGck3S1p+viFaWY2tEZ6YrcBG6vqHwM+GRHnAvuAW5oZmJlZPepKYpIWA28A/iavC7gKuC9fZQ3w5nGIz6a46J9WuB2bN1C4dR5Q+WZTU709sU8BH6JyPpV5wP6I6M/r24GzmhuamdnoRk1ikq4HeiLixyezAUm3Slovaf1AX9/JPIWZ2bDq2cXi1cCbJL0emAm8DPg0MFtSR94bWwzsGOrBEbEaWA0wY8kS/wZuZk01ak8sIj4cEYsjYilwE/DdiHg78DDw1ny1VcD94xalTVlxbFrhNn1PqXDrOEz5ZlPTWPYTux14v6RnyObI7mpOSGZm9Wtoj/2IeAR4JC8/C1za/JDMzOrnw45sQiv1Fs9wMTCzOK3a9/JKvfNgcWDR31Vc9xWXPFcub9mwuFkhWpv5sCMzS5qTmJklzUnMzJLmOTGbNI6cfaxQv/2ybxfqf7Xp91oZjrWIe2JmljQnMTNLmpOYmSXNc2I2eRwr/k++/6bfKdT7/v2p5bL/e08e/izNLGlOYmaWNA8nbdIo9RX/J29eNbtQn3akUlZ/8UywXT2V+qnX7Cq07X7yjOYEaOPCPTEzS5qTmJklzUnMzJLmOTGbkqYdL9Yf/8Bny+Xrf3FdoW13KwKyk+aemJklzUnMzJLm4aRNScdO7y/Uf3/RinJ585cWFdr8n35i8+djZklzEjOzpDmJmVnSPCdmU1LpUPEqSls+flm5PG1Pcd2RrqJUfQUl8FWU2sE9MTNLmpOYmSXNSczMkuY5MbNRjHQVJV9Bqf3qSmKStgK9wADQHxErJc0F7gaWAluBt0XEvvEJ08xsaI0MJ6+MiBURsTKv3wGsjYhlwNq8bmbWUmMZTt4AXJGX1wCPALePMR6ziWeEC5BUX3wEPMncDvW+5wH8k6QfS7o1X7YgInbm5V3AgqZHZ2Y2inp7Yq+JiB2SzgAekvR0dWNEhKQY6oF50rsVoDRnzpiCNTOrVVdPLCJ25Pc9wDeBS4HdkhYC5Pc9wzx2dUSsjIiVpe7u5kRtZpYbtScmqRuYFhG9efl1wH8FHgBWAXfm9/ePZ6Bm7TLSVZSqr6AEEAuOFuod22ZW2opHOjHYOeTgxRpUz3ByAfBNSSfW/z8R8R1JPwLukXQLsA142/iFaWY2tFGTWEQ8C1w4xPI9wNXjEZSZWb38i7CZJc2HHZk1UfQX+wXH5g2Uy7O2Ff/cjs73nFgzuCdmZklzEjOzpHk4adZEUXOI0vQ9lf0qOg4X1y3ujGEnyz0xM0uak5iZJc1JzMyS5jkxsyYq9RaPLRqYWdmNou/lxV0qLrhkS7n87L3LCm0z9hfX3bPCu2MMxz0xM0uak5iZJc3DSbM2eeKxc8vl/nMHC22nbartX3g4ORz3xMwsaU5iZpY0JzEzS5rnxMza5PgZx8tlHS7umlE66jmwerknZmZJcxIzs6Q5iZlZ0jwnZtYmpX3D//ntP784JzZn+d5K5evzCm19i1WoH5td3OdssnNPzMyS5iRmZknzcNIsAUe+N79cfuMfP1po++q6ywr12ov9TnZT69Wa2aTjJGZmSXMSM7OkeU7MLAGHzjtWLn//T4tzYDMuKf4Z95/iXSxeQtJsSfdJelrSRkmXS5or6SFJm/P7OeMdrJlZrXqHk58GvhMR5wEXAhuBO4C1EbEMWJvXzcxaatQkJuk04HeBuwAi4lhE7AduANbkq60B3jw+IZqZDa+eObFzgBeAL0q6EPgxcBuwICJ25uvsAhaMT4hmVn2I0o4ra1uLc2ADpwyUy7O2dhafpzK1xuFFk2PurJ7hZAdwMfC5iLgI6KNm6BgRwTAnAZd0q6T1ktYP9PWNNV4zs4J6kth2YHtErMvr95Eltd2SFgLk9z1DPTgiVkfEyohYWerubkbMZmZlow4nI2KXpOckLY+ITcDVwFP5bRVwZ35//7hGamZ10ZHKWWI7XrWv0Na3eXaLoxl/9e4n9h+Ar0iaDjwLvJOsF3ePpFuAbcDbxidEM7Ph1ZXEImIDsHKIpqubGo2ZWYN82JGZJc2HHZlNMh1nvFgud913WqGt9/LKbhXTjhTPCJsq98TMLGlOYmaWNA8nzSaZwee7yuVfX1Rsqx5CzttQHE72XHG8UJ+/4GC5vG/T3CZG2FzuiZlZ0pzEzCxpTmJmljTPiZlNUcdv3FuozzgyvVA/0NtFCtwTM7OkOYmZWdKcxMwsaZ4TM5uiep8q7vs15+li+97XvUgK3BMzs6Q5iZlZ0jycNJuiorN4WYy9r6xZYefMcvHMx4rr7rimWJ+/eH+53OpDlNwTM7OkOYmZWdKcxMwsaZ4TM7NR9dxY3N2itH1WoX60v0S7uCdmZklzEjOzpDmJmVnSPCdmZqOa/e3uQn3PiuJ+YocOVk7b0+qekXtiZpY0JzEzS5qHk2Y2qr0XFIePGiy264XKWWE1ULyKUteuqnrN9XoPL6p5opPgnpiZJW3UJCZpuaQNVbeDkt4naa6khyRtzu/ntCJgM7NqoyaxiNgUESsiYgXw28Bh4JvAHcDaiFgGrM3rZmYt1eic2NXAlojYJukG4Ip8+RrgEeD25oVmZikqvVic+Dr+skr5yJn9xXX7xj6j1egz3AR8NS8viIideXkXsGCoB0i6VdJ6SesH+vpOMkwzs6HVncQkTQfeBNxb2xYRAcRLHpS1rY6IlRGxstTdPdQqZmYnrZHh5HXATyJid17fLWlhROyUtBDoaX54Zpaa7h3F/kzvlYfL5Wk7m39B3kaGkzdTGUoCPACsysurgPubFZSZWb3qSmKSuoFrgG9ULb4TuEbSZuC1ed3MrKXqGk5GRB8wr2bZHrJfK83M2saHHZlZUx1YXvMb3/OVebCao46Ys3xvoT7tnkpfadbu44W2Xw6zPR92ZGZJcxIzs6Q5iZlZ0jwnZmZtc/S78wv13ssrhyXNeKGzuPJ3hn4O98TMLGlOYmaWNA8nzaxtDp0zUKh3HKhchHewzuzknpiZJc1JzMyS5iRmZklTdiqwFm1MegHYBswHft2yDY/O8YxsosUDEy8mxzOyZsRzdkScXruwpUmsvFFpfUSsbPmGh+F4RjbR4oGJF5PjGdl4xuPhpJklzUnMzJLWriS2uk3bHY7jGdlEiwcmXkyOZ2TjFk9b5sTMzJrFw0kzS1pLk5ikayVtkvSMpLZcMVzSFyT1SHqiatlcSQ9J2pzfz2lhPEskPSzpKUlPSrqtnTFJminph5J+lsfz5/nycyStyz+7u/NL+LWMpJKkn0p6sN3xSNoq6eeSNkhany9r23co3/5sSfdJelrSRkmXt/E7tDx/b07cDkp633jF07IkJqkE/BXZpd/OB26WdH6rtl/lS8C1NcvuANZGxDJgbV5vlX7gAxFxPnAZ8J78fWlXTEeBqyLiQmAFcK2ky4CPAZ+MiHOBfcAtLYrnhNuAjVX1dsdzZUSsqNptoJ3fIYBPA9+JiPOAC8neq7bEFBGb8vdmBfDbwGHgm+MWT0S05AZcDvxjVf3DwIdbtf2aWJYCT1TVNwEL8/JCYFM74sq3fz/ZlaXaHhMwC/gJ8CqyHRU7hvosWxDH4vxLfxXwINmp2tsZz1Zgfs2ytn1ewGlkp6DXRImpKobXAd8fz3haOZw8C3iuqr49XzYRLIiInXl5F7CgHUFIWgpcBKxrZ0z50G0D2QWRHwK2APsj4sQZ61r92X0K+BAwmNfntTmeAP5J0o8l3Zova+d36BzgBeCL+ZD7b/LLLE6E7/VNVK5XOy7xeGK/RmT/Jlr+k62kU4CvA++LiIPtjCkiBiIbCiwGLgXOa9W2a0m6HuiJiB+3K4YhvCYiLiabGnmPpN+tbmzDd6gDuBj4XERcBPRRM1Rrx/c6n6d8E3BvbVsz42llEtsBLKmqL86XTQS7JS0EyO97WrlxSZ1kCewrEXHiAsVtjQkgIvYDD5MN12ZLOnGGp1Z+dq8G3iRpK/A1siHlp9sYDxGxI7/vIZvruZT2fl7bge0RsS6v30eW1Nr9HboO+ElE7M7r4xJPK5PYj4Bl+a9K08m6mQ+0cPsjeQBYlZdXkc1LtYQkAXcBGyPiE+2OSdLpkmbn5S6y+bmNZMnsra2OJyI+HBGLI2Ip2XfmuxHx9nbFI6lb0qknymRzPk/Qxu9QROwCnpO0PF90NfBUO2PK3UxlKMm4xdPiSb7XA78gm2P5SKsnGfMYvgrsBI6T/Qe7hWyOZS2wGfi/wNwWxvMasm7148CG/Pb6dsUEXAD8NI/nCeA/58tfAfwQeIZseDCjDZ/dFcCD7Ywn3+7P8tuTJ77H7fwO5dtfAazPP7e/A+a0+XvdDewBTqtaNi7xeI99M0uaJ/bNLGlOYmaWNCcxM0uak5iZJc1JzMyS5iRmZklzEjOzpDmJmVnS/j9PKnbhLQJi9QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(D.real.detach())\n",
    "plt.title(r'$\\bf{D}$: Block Diagonal Matrix')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fuzzy-monte",
   "metadata": {},
   "source": [
    "We have to reshape the input into a batch of vectors, where we lay out the channels first. That is,\n",
    "the first `cin` elements of the vector belong to the first pixel, the next `cin` to the second pixel,\n",
    "and so on.\n",
    "\n",
    "Then, we multiply by the diagonal matrix `D` and convert back to the original input shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "refined-bedroom",
   "metadata": {},
   "outputs": [],
   "source": [
    "yfft_matrix = (D @ xfft.reshape(batches, cin, n**2).permute(2, 1, 0).reshape(-1, batches)) \\\n",
    "                .reshape(n, n, cin, batches).permute(3, 2, 0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fallen-evening",
   "metadata": {},
   "source": [
    "Again, this is the same as doing circular convolution using `nn.Conv2d`, as well as our previous implementation of Fourier convolutions. We apply the inverse 2D FFT and compare:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "verified-carter",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.5019836610008497e-06"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y3 = torch.fft.ifft2(yfft_matrix)\n",
    "(y1 - y3).norm().item()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "exciting-start",
   "metadata": {},
   "source": [
    "-----------------------------------------\n",
    "We will now illustrate the \"block diagonalization\" part of this figure by constructing DFT matrices that convert the previous block diagonal matrix `D` into a block matrix of convolution matrices `C`:\n",
    "<center>\n",
    "    <img width=\"300\" src=\"https://github.com/locuslab/orthogonal-convolutions/blob/main/img/ICLR-Thumbnail-Fourier.001.png?raw=true\" alt=\"Block-Diagonalization in the Fourier Domain\" />\n",
    "</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "usual-cartoon",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Somewhat easier to build these matrices with np.block\n",
    "import numpy as np \n",
    "\n",
    "# Perfect shuffle matrices, implemented as described in the appendix of our paper\n",
    "def shuffle_matrix(p, q):\n",
    "    r = p * q\n",
    "    Ir = np.eye(r)\n",
    "    S = []\n",
    "    for i in range(q):\n",
    "        for x in Ir[i:r:q]:\n",
    "            S.append([x])\n",
    "    return np.block(S)\n",
    "\n",
    "# The \"script F\" matrix from Corollary A.1.1 of our paper\n",
    "def script_F(c):\n",
    "    global n\n",
    "    S = shuffle_matrix(c, n**2)\n",
    "    F = np.fft.fft(np.eye(n), norm=\"ortho\")\n",
    "    FF = np.kron(F, F)\n",
    "    return S @ np.kron(np.eye(c), FF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "colonial-action",
   "metadata": {},
   "outputs": [],
   "source": [
    "Fs_cout = torch.tensor(script_F(cout), dtype=torch.complex64)\n",
    "Fs_cin = torch.tensor(script_F(cin), dtype=torch.complex64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "polar-worship",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUEAAAFBCAYAAADpDh0xAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAArG0lEQVR4nO2de9wcZXn3vxdPSMjJhADNmQQEgggSMaUg1FIOCrwI9PO2VsCKbSylRYVGDEGsxVZbtFjFA2r6aov1AIggiHiACNZDBQIG5GAIp0CAJBIJxHAwCdf7x8yzz8zm2dl795nZnd35fT+f5/PMPTM7c+3MPdfude19/W5zd4QQoqrs0G0DhBCim8gJCiEqjZygEKLSyAkKISqNnKAQotLICQohKo2coBCi0sgJCiEqTSmdoJk9amZHl+U4eWNm95rZETkfc56ZrTCzTWb2njyP3YINDa93kfcij2PnbV8R9zhx7FL265FS5DXLom0naGaHmdm3zWyDmb1oZg+Z2afNbHTAax81sxfM7Ldm9oyZfcfMZrdry0iJ7fmdme1at/4XZuZmNreF4zTtnO7+ane/pT1rG7IYuNndJ7r7pxrY9kLsJDea2c/M7EwzK+UHYZKy9ZeEXaea2fLYrqfM7LtmdjgUdo9zpVl/7ZPnoiltPQBm9lbgR8AJwOPAfwMPA2cC4wIP82Z3nwBMB9YBn27Hlhx5BDhlsGFmBxD+XoIws1F5Hq+OOcC9TfZ5s7tPjPe9CDgP+GKBNuVJqfqLmS0CPgn8CzAV2B24FDipyeuK7ANF0OvPRXPcvaU/oguwAXAi57dDYtsrgdEBx3gUODrRPh54IGP7q4BbgI1ED/qJiW2zgauBX8d2fWa448THeAQ4pYE9HwBuT6y7GLggfp9zE+uXAA8Bm4D7gD+J1/838DLwAvBbYHHi2OcBdwMvAaMG7Yqv12+Ag+J9Z8Tv44gG123Y6wD8ENgGvBife59m1zxed3Bs8/6B19qBvRLt/wI+XHeO8+Pr8gzwn8BOw9yL9wHfrLPlU8Al3eovoX0l3j4pvs5/FtLHG/SBrH4bcp0Hjz1sf0zsd2583meBKxL3Y9j+WobnIuC5ngF8M972CPCexLbzgCdiO1YCRzX1R6HOL3GSY+IL4MC8BvtcClwa2EHGAZcBX26wfUfgQeD9wGjgyPgNzgMGgLuATwDjgZ2Aw+uPAxwEPAackGVPfNFeFR93DdE3pvqb/WfxTdgB+HNgMzA9w9E8CqyIb+rYYd7fX8edZhzwfeDiBjY2vA7x9luAd4Zc87r1jwF/G3iOkIfznvi9TgF+Ori97j1Pj6/b5Lg9ClgPvK5b/SW0r8T7HgtsBUa14ARrfSDLjjacYLP+eFu8fQpwP3Bmsz7R7eeiyX3aAbgD+GB8f/ckikLfFN/jx4EZ8b5zgVcW4QRPY8gJ7tTq6xNv9rdEn9RbgCeBAxrc5D8E1pL+xvl14ELgUKJPg2E7Y3ycD8U37oiAm/0B4F+JOvmNRA9n6mYP89oVwElNbvZfZTkk4Drgl0SfimManKfhdYiXb6E9J/hz4ILAc4Q8nMmH7HjgoQbv+bvAX8fLJwD3dbO/hPaVxDOwNqCPJ53gXyW2Neu3wU4woD++LbHtY8DnQ47TzeeiyX36A+CxunXnE0UdexF9mB4N7Jh1f5J/7eQE1yeW57Tx+kFOdvfJRF7+XcCPzGzaMPvNAB5395cT61YDM4k+RVa7+9aM85wJ/MzDEq7/DZwKvAP48nA7mNnb419hN5rZRmB/YNfh9k3weJPt/xEf59Pu/lKDfbKuw0iYSRSS53WO5HtdHR9zOC4D3hYvv43o2mNmp8U/NPzWzL6b2L8T/WW7vtLAng3Ari3mspLXJaTfBhHQH9cmlp8HJrRxmk4/F1nXZw4wY/A88bneD0x19weBc4g+8Nab2eVm1qj/1WjHCf6MKN8D8IHkr4tmNsfMdmzlYO6+zd2vJsppHT7MLk8Cs+t+xdydKO5/HNi9SWc8M97nEwG2rCbKMRxPlI9IYWZziBzWu4Bd4ofyHsAGD9Ho0I3OaWYTiBLsXwQuNLMpDXbNug5tYWa/T+QcfhJ4judJJ8WHc0LJX213j485HN8CXmNm+xN9E/wqgLt/1d0nxH/H1b+o4P6yXV9pYM//EuWxTm5wnOFI9oFmdoRc55D+2IpNjXfq/HORdX0eBx5x98mJv4nufnxs69fc/XCGQvaPNnt/LTtBd98MvJso2fk24E4zW2pm1wMPEMXwwVjEScDORDmLem4l6hSLzWzHeBzRm4HLifIdTwEXmdl4M9vJzA6re/0moq/xbzCziwJMWggcGb/PesYTXdhfx7b/JdEn3iDriHIUrXAJsNzd3wl8B/h8g/2yrkNLmNkrzOyE+LVfcfdfBp5jBXCqmQ2Y2bHAHw1z+LPMbFbszC8gSsZvh7u/CFwFfA24zd0fC7S9yP4S1Ffc/VminNRnzexkMxsXn+s4M/tYwNtoZscKml9naN4fm9FKf+3kc5F1fW4DNpnZeWY2Nr5G+5vZ71s0VvZIMxtD9CPhC0R+KpvQuHmY2PwNRA/tb4g+FR8CPkOUrPw8idxDg1zD4K9Fm4g+NU5rlKsAXk00JOdZtv8FbHeibxUbgKeBTw13HKLE8F3APzfKfQyzfrvcB/CR+D0/Dfx7bNc7420nESXVNwLnNsmHHB3v/wQwJV4/gSipf1qD65Z1HW6heU7whfh6P0v0beYsYKCFcywg+rV1E1GI9HUa/zq8kSjkHZdxHQ6Pr+9fNulrhfeX0L5SZ9dpwHKiHwHWEj0Prx/meMO996x+G3KdB4+d1R/rr8uFRB96g+3t+msZnouA6zMjviZriaLSnxM9T68hdpKxLdcT/0iS9WfxQYXoOGa2O/ArYJq7P9dte0Q1KX21gOhP4pzdIuByOUDRTXpt9LroA8xsPFGeaDVRDk6IrqFwWAhRaUYUDpvZsWa20sweNLMleRklhBCdou1vgmY2QDQk5hiiUfa3E9Va3pefeUIIUSwjyQkeDDzo7g8DmNnlRD+FN3SCA+PH+46Th8YCewvfQ0eNHRo8PvBEeizoy6OHDrRlQvY40TmTf11bXr1xt3ADOsTA74aWZ/zehtS2xzfu0vB1O47dkmrbuqFbu8OW9FCpF3cZaHiciRNeSLW3JW7S8y/VqaRtbXytd0ibw8stDaEvnh3qahEGEnU6W5qMdB0zbugmbXkufU1GP5N+4y/u1viNT530bG153aZXpLbNesUzqfaajY3G0MOYNemhey/NammobuHs+Nuh5YHn0v3rxRljG75uzMZ0v9302yefdvfcH9qROMGZpMte1hDV9TVkx8lTmH3W39fa28aGfwvdbb8h5zXpAzultj0/e+imP/mH2Z71syf/R235b77118Hn7xQTHxmy/x/PTlconXvd2+p3rzFj/3Wp9qiPDznMMevSD8mqt09qeJw/OuyeVPvZLUPX+s6H0lWSO2xo/ICPXZu+Dy9Maz5mtZOMeTpt3+SHhuxbl9mLYY/5Q0U6Ty5LyxrO+cbaVPuBM3+v4XH+/rjv1JYv/p90ccxFR6bHwJ/37VMbHueV7/15qv3QokMa7tsNpv9k6Dmf8P1fpratWvSahq/b47p0BekPb7lgdb6WRRQ+RMbMzrBIeHL5ts3DDTYXQojuMRIn+ATpOtFZDFPH6u5L3X2Buy8YGF+ur+lCCDGSH0ZGEf0wchSR87sdONXdG6obj99nuu//6dNr7d/cMRQqtBIafyER0gL821veWltOhsbQPDwuG8mwdtO3p6e2JcPjrNAY0mHtIx/cN7WtlfA4GdbOO/GB1LZWwuOykQxpAZ68aejzPBkaQ3Z4nAxpAS79yptT7WR4nBUa14e0D3z+4FT74kR4nBUal5FkWPvilHQOdbvw+J8bh8ePLDr3DndfkK91I8gJuvtWM3sXkRDoAPClLAcohBBlZEQVI+5+A3BDTrYIIUTH6WjFyJjZs33monNq7eQvvsnQGNoPj5OhMfR2eFz/i28yPG7ll+P6X3yzwuPQ0Biyw+NeCo2h7hffm9K/+OYVHrfyy3FWeHxxC78cl436X3yzwuP60LiocLh3PIIQQhSAnKAQotLICQohKk1Xc4JJkvlBaD9HmDV8BlqrLikbocNnoP0cYbvDZyCdI+zX4TPQWnVJMkfY7vAZSOcIs4bPQP/kCOuHz/xg85eVExRCiLyRExRCVJrShMP1ZIXHU5dvS23LCmuzwuMHT50YZEsnufjEr9SWWxFMqA+Pp/14SIUkK6SF7OqSR08d6h/NQtqs6pIVt+2V+dpOM/XWdDsrrM0Kj+dc1Z5gAqTD4396x1dS21oRTKgPjwc2l+u7zd7/cHdtOasiBLKrS3529fsUDgshRN7ICQohKo2coBCi0nQ0Jzhu6mzf+88X1dqtqKIkc2Drbp+W2laEOCuE5xrLKM6aFA1tpeyrPv+17aNTh46ZkzgrhJfYlV2cNSkYCq2poiTzX08ellZYblecFdICra0Mn9lOnPXj5RJnVdmcEEIUgJygEKLSyAkKISpNR3OCE/aZ5vM/+/Zau11pqKwxhHkpVENvl9jlJQ1VhEI19FGJ3QikoZLkpVAN4SV2vVReB8oJCiFEIcgJCiEqTVfL5vJSRSlCoRr6ZwKnvFRR8lKohvAJnHopNIbWVFHaDY+rqlCtcFgIIQpATlAIUWnkBIUQlaY0UlrNpKFCS+zyUqiG8Fnseik/CPlNPC6F6uaETjzeTGIqVKEawnOEvaZQrZygEEIUgJygEKLSyAkKISpNR3OCex0wzj/2rXm1drvy8Zv2CJdTqqpMf1I+vtmMaFkldrvfkL9MP4TLZ5Vdpj8pHQ/h0lmQzhFOvCfdT9uV6Ye0VH8/yfR3LSdoZl8ys/Vmdk9i3RQzu9HMVsX/d87bMCGE6AQhrv6/gGPr1i0Blrn73sCyuC2EED1HUDhsZnOB6919/7i9EjjC3Z8ys+nALe4+L+sY0Npsc1lMfCTtu/tFoRrKp1JdhEI1hJfY5aVQDekwu18VqiGtUt1PCtVlGyIz1d2fipfXAlOzdhZCiLIy4synR18lG36FMrMzzGy5mS3ftnlzo92EEKIrtOsE18VhMPH/9Y12dPel7r7A3RcMjB/faDchhOgK7eYE/w3Y4O4XmdkSYIq7L252nLxyglnDZ/JSqIbwHKEUqiNUYtc+eSlUQ3iJXV4K1dCZErtuDpH5OvC/wDwzW2NmC4GLgGPMbBVwdNwWQoieY1SzHdz9lAabjsrZFiGE6Di9FZsJIUTOlEZKayQUIdMPxcxiJ5n+iNBZ7PKS6Yf+yRGWTaYfOjOLXdnGCQohRF8gJyiEqDR9EQ4nyUuhGoqZxS5r+Az0T3hctuEz0D+z2GUNn4Hw8DgvhWrozCx2CoeFEKIA5ASFEJVGTlAIUWl6Mid48YlfSbX7RaEayqdSXYRCNYTLZ+WlUA3pHGG/KlRDWqW6EwrVkF1il5dCtXKCQghRAHKCQohK05PhcF5kKVRDeJidl0I1pFWq+1WhGsKrS/JSqIZ0mN2vCtWQVqnuJ4Xqm/wqhcNCCJE3coJCiEojJyiEqDSVzgnmVWKXl0I1hJfY9VJ5HfRWiV0vq89AeIldXgrVEF5iNxKFag2REUKIApATFEJUGjlBIUSlqXROsJ5OzGInheqI0Fns8lKohvAcoRSqt6dZiV3oLHYjUahWTlAIIQpATlAIUWkUDmcQOoGTFKpHRieGz0B4iZ0UqiO6PYFT/fCZt+x9p8JhIYTIGzlBIUSlkRMUQlSaSucEy6ZQDWmV6n5VqIZw+ay8FKohnSPsV4VqSKtUd0KhGrJL7PJSqH7sbxZ3JydoZrPN7GYzu8/M7jWzs+P1U8zsRjNbFf/fOW/jhBCiaELC4a3Ae919P+AQ4Cwz2w9YAixz972BZXFbCCF6iqZO0N2fcvc74+VNwP3ATOAk4LJ4t8uAkwuyUQghCqOlnKCZzQX+B9gfeMzdJ8frDXhmsN2IV4yf4Yfsd0at3Yo0VDIHVjbp+Lxk+iEt1T8Smf6ZE56tLbdS9lWf/yqbfHy7Mv2QzjXmJdP/wxX7pba1Kx//0McPyTxnNyibTP97XnVzd8cJmtkE4JvAOe7+XHKbR5502CfWzM4ws+VmtnzL1udHZKwQQuRNkBM0sx2JHOBX3f3qePU6M5seb58OrB/ute6+1N0XuPuCHUeNy8NmIYTIjabhcBzqXgb8xt3PSaz/N2CDu19kZkuAKe6+OOtY+xywk3/2urm1dr+oouSlUA35ldj1rSpKyRSqR6KKUnaKKLEbiUJ1USoyowL2OQz4C+CXZrYiXvd+4CLgSjNbCKwG3pK3cUIIUTRNnaC7/wSwBpuPytccIYToLL0TcwohRAF0tWyuX6WhyqZQ3dfSUIEK1RCeIxyJQnUr0lC9lCNsVmIXOovdSBSqpSwthBAFICcohKg0pVGRaaaK0ssTj4cqVEN4eJzX8Bnon4nHs4bPQHh1SV4K1RA+8XgvhcbQmQmc6ofPfP+BjykcFkKIvJETFEJUGjlBIUSl6WhOcMLOs3z+H59da7ernFw21eS8FKohXKW6mUL1llm/qy03y+tlldiVTTm5XYVqSOcI81KoPuszV6a2harGQDpHOLC5fN9HQhWqIbvELi+F6l99ZJFygkIIkTdygkKISiMnKISoNJ0dJ7jHLJ924btq7V1+OpQ72PiqJnZMG8o5TN3l2dSmP5r6YG35yh++PvMw+3zovtryA/+4X8ae3WHMhqHPpWlHrkltW33XjIav2zZpa6p9/dGfri2/457TU9ueWTml4XFu+NOPp9rvPPvva8uPn5DOVw4821h/Y+8LVqTaqz4yv+G+3WCv85an2iu/ML+2PLAxW1dk9g8S1/rcdG5266XTUu2nDm+kPQInHXFbbfman/5+atse16bv5+rjG+d19/5y+nloNnax02wbN9Rv9rgm3Ycee1Pja/1y3abV71HZnBBC5I6coBCi0nQ2HN5zps/8yFm19rbfDfngZGgM2eFxMqQFePbKXWvLydAYmofHZSMZ1o57OH1NkuFxVmgM6bB2i6c/61oJj5Nh7eZj08MkWgmPy0YqpAXGLFtRW06GxpAdHidDWoA7f5Mu1UuGx1mhcX1Iu3LhK1LtZHicFRqXkWRY66PTfaaV8FgqMkIIUQBygkKISiMnKISoNN2V0koMe0nmB6H9HGEyPwi9nSOsH/aSzBG2MnymfthLVo4wND8I2TnCXsoPQjpHmMwPQn45wlaGz2TlCFsZPlM26oe9ZOUI6/ODygkKIUQByAkKISqNnKAQotKURl4/mR+E7Bxhu2MIobUSu7IROoYQ2s8RtjuGENI5wn4dQwitldglc4TtjiGEdI4wawwh9E+OsH4M4Y++v0Q5QSGEyBs5QSFEpeloODxut9m+758MqZJkKsdkhMenzL89tS0rrM0Kj9fdm61y2w3mHPhkbbkV1Zj68Hjc65+uLWeFtJBdYnfCTe+uLTcLabNK7Mo2Q+C2yenrlxXWZoXHm9/8utS2UNUYSIfHO1yUTtu0ohpTHx7vsCU7tO40u39/6PpllcVBdondY2cuVjgshBB509QJmtlOZnabmd1lZvea2Yfi9XuY2a1m9qCZXWFmo5sdSwghykbIN8GXgCPd/UBgPnCsmR0CfBT4hLvvBTwDLCzMSiGEKIiWcoJmNg74CfC3wHeAae6+1cwOBS509zdlvX7s9Nk+d+GiWrsVaahUDmwgbXMRCtUQnmsso0J1Ujm5lbKv+vzXZ76Qv0I1hJfYlV2hOqmaDO0rJ8/Z76nUtnYVqiGtUt1PCtVdLZszswEzWwGsB24EHgI2uvvgFV4DzMzbOCGEKJogJ+ju29x9PjALOBjYN/sVQ5jZGWa23MyWb31+c3tWCiFEQbQ8RMbMPgi8AJxHi+HwmLmzfNoHh4ZctKuKkjV8Ji+Faujt6pK8VFGKUKiG/qkuGYkqSpK8FKohvLqklypLoIvhsJntZmaT4+WxwDHA/cDNwJ/Gu50OXJu3cUIIUTQhH7nTgcvMbIDIaV7p7teb2X3A5Wb2YeAXwBcLtFMIIQqhqRN097uB1w6z/mGi/KAQQvQsXVWRyUsVpQiFauifWezyUkXJS6Eawmex66X8ILSmitJujrCqCtVSlhZCiAKQExRCVBo5QSFEpSmNsnQzaajgErucFKohfBa7XsoPQnaOsN0xhCCF6uHIkoYKHUMI4QrVEJ4j7DWFauUEhRCiAOQEhRCVpqPh8MR503zB506rtdtVTh67Pm2zFKq3J6mc3GwyoKwSuw3f2rO2nJdCNYQrx5RdoTqpmgzhqjGQDo9n/SAdwrarUA1plep+UqhWOCyEEAUgJyiEqDRygkKISlOaITItHWdD2ncXoVAN4bnGvBSqoXwq1UUoVEN4iV1eCtWQzjX2q0I1pFWq+0mhWjlBIYQoADlBIUSlkRMUQlSanswJZo0hzEumH8JL7CTTH6ESu/bJS6Yfwkvs8pLph86U2CknKIQQBSAnKISoND0ZDtdThEI1FDOLnRSqI0JnsctLoRr6Jzwum0I1dGYWO4XDQghRAHKCQohKIycohKg0fZETTJKbQjUEz2KXl0I19E+OsGzDZ6B/ZrHLGj4D4TnCvBSqoTOz2CknKIQQBSAnKISoND0ZDs858MlUu18UqqF8KtVFKFRDuHJMXgrVkA6P+1WhGtIq1Z1QqIbs6pK8FKoVDgshRAEEO0EzGzCzX5jZ9XF7DzO71cweNLMrzGx0s2MIIUTZaOWb4NnA/Yn2R4FPuPtewDPAwjwNE0KIThCUEzSzWcBlwEeARcCbgV8D09x9q5kdClzo7m/KOk4nhsi0QpZCNbSQa8xJoRrSKtX9qlAN4SV2eSlUQzrX2K8K1ZBWqe4nheofrPjnruYEPwksBgav6C7ARncffEdrgJn5miaEEMXT1Ama2QnAene/o50TmNkZZrbczJZv27y5nUMIIURhhAybPww40cyOB3YCXgFcAkw2s1Hxt8FZwBPDvdjdlwJLIQqHc7FaCCFyoqVxgmZ2BHCuu59gZt8Avunul5vZ54G73f3SrNeXLSeYW4ldTgrVEF5i10vlddBbJXa9LMEF4SV2eSlUQ3iJ3UgUqss4TvA8YJGZPUiUI/xiPiYJIUTnaOkjzt1vAW6Jlx8GDs7fJCGE6Bw9WTZXFJ2YwEkK1RGhEzjlpVAN4eGxFKq3p1mJXegETiNRqC5jOCyEED2PnKAQotLICQohKo1yghmEzmKXl0I1hOcIpVC9PVnDZyC8xE4K1RHdnsWufvjMLTedr5ygEELkjZygEKLSyAkKISpNpXOCZZPph7RUf7/K9EO4fFZeMv2QzhH2q0w/pKX6OyHTD9kldnnJ9K9+9/uUExRCiLyRExRCVJqOhsMT9pnmB3xmaJhCK6ooyfCvbKrJuSlUQ0qleiQK1TO+MRR2tFL2VR/6lU05uV2FakiH2XkpVB977aLUtnaVk1e9fVLmObtB2RSqP3nQlQqHhRAib+QEhRCVRk5QCFFpOpoTnDRqVz90wkm1dr9IQ+WmUA3BJXbNFKr7VRqqdArVI5CGKjtFlNiNRKFaUlpCCFEAcoJCiErT1YqRflVFKZtCdT+rooQqVEN4eDwiheoWVFF6KTxuVl0SOoHTSBSqFQ4LIUQByAkKISqNnKAQotKURkWmmSpKL088HqpQDeE5wpEoVPfrxONZw2cgvMQuL4VqCJ94vJfyg9CZWezqh8/87GqpyAghRO7ICQohKo2coBCi0nQ0Jzj91Tv76V87stZuVzm5bKrJeSlUQ1qleiQK1WwbGmPVLK+XVWJXNuXkdhWqIZ0jzEuh+j0L35XaFiqdBekc4Q5bspWbu0GoQjVkl9jlpVD9P0d/vJCcYFDW28weBTYB24Ct7r7AzKYAVwBzgUeBt7j7M3kbKIQQRdLKx/wfu/v8hCdeAixz972BZXFbCCF6iqBwOP4muMDdn06sWwkc4e5Pmdl04BZ3n5d1nJ1mzPa57xxS4v3gX3y9tvyB6/8824jpL9YW571vfWqTf2Vo+YE7d888zOLjr6stf+yGE7PP2QVeHjN0P6wuRLKX6/ceYtv49MZX/cMjQ42X0qFzljL3nO9sSbVf/bG7a8ufmnF7atteXz+z4XHOPva7qfYl3zuu4b7d4Nzjvp1qX/yLNw411o7JfG0yVXP/v++T2jb7+vT3iqx0wodPuKK2vGxj+p6cvMsdqfbZ176j4XEmPpI+x6Y9MjpKF0imgOZcne7Ta44aaPi6V74v3d9u2nZFV4fIOPADM7vDzM6I101190Hd7LXA1LyNE0KIogkdCXu4uz9hZr8H3Ghmv0pudHc3s2G/UsZO8wyAUZN2HpGxQgiRN0HfBN39ifj/euAa4GBgXRwGE/9f3+C1S919gbsvGBg3Ph+rhRAiJ5rmBM1sPLCDu2+Kl28E/gk4Ctjg7heZ2RJgirsvzjrWmD1n+sx/+btae/R942rLyfwgZOcIk3k9gGuPnl9bTuYHoXmOsGwkc3u2tS4nmMgRZuUHIZ3bG/PLx9IbW8gRJnN7q15ID01qJUdYNuqHYJ3484dqy6n8IGTmCJN5PYDzf/x/U+1kjjArP1if11vwF3el2skcYVZ+sIwkc3sPLH1talsrOcKipLRCwuGpwDVmNrj/19z9e2Z2O3ClmS0EVgNvyds4IYQomqZO0N0fBg4cZv0Gom+DQgjRs5SrHEAIITpMd6W0EmP/kvlBaD9HmMwPQm/nCOvH/iVzhK2MIawf+5eVIwzND0J2jrCX8oOQzhEm84OQX46wlTGEWTnCVsYQlo36sX9ZOcL6/KDk9YUQogDkBIUQlaY0ytLJ0Biyw+N2h89AayV2ZSN0+AyMIDxuc/gMpMPjfh0+A62V2CXD43aHz0A6PM4aPgP9Ex7XD5/58XfOUzgshBB5IycohKg0coJCiErT0ZzgngeM9w9f/epaO1M+KyNHOPeoR1PbsnJ7WTnClYvKlxNMTmzWinRWfY5w3/ffX1vOyutBdondDXfdWFtultfLKrG74Uevy3xtx6lT5s7K7WXlCL/2/v+T2hYqnQXpHOFnjkiP5WpFOqs+R3jzTw5o+NpuMGvZttpyVlkcZJfYPbZwiXKCQgiRN3KCQohK09khMnNm+/Tzzh46eQuqKMnwb+yadMlzEQrVEB5ml1GhOqmc3ErFQ33ox5jEvjkpVEN4dUnZFarrJ85qVzm5PvRrV6Ea0irV/aRQrYoRIYQoADlBIUSlkRMUQlSaDucEZ/m0CxI5wTZVUbKGz+SlUA29XWKXlypKEQrV0D8ldiNRRUmSl0I1hJfY9VJ5HSgnKIQQhSAnKISoNHKCQohK01UprbykoYpQqIb+mcUuL2movBSqIXwWu17KD0Jr0lDt5girqlCtnKAQQhSAnKAQotKURlm6mSpKcIldTgrVED6BUy+FxpDfxONSqG5O6MTjzdRVQhWqITw87jWFaoXDQghRAHKCQohKIycohKg0nc0J7j7bZ7z3nKGTt6mc/C/HXJnaJoXqYUgqJzeZES2rxO66Q15ZW85LoRrC5bPKrlCdVE2GcOksSOcIRyKdVZ8jTKpU95NCdVdzgmY22cyuMrNfmdn9ZnaomU0xsxvNbFX8f+e8jRNCiKIJDYcvAb7n7vsCBwL3A0uAZe6+N7AsbgshRE/R1Ama2STgDcAXAdz9d+6+ETgJuCze7TLg5GJMFEKI4miaEzSz+cBS4D6ib4F3AGcDT7j75HgfA54ZbDcia5xgK7w8Jm1zETL9EJ5rzEumH8on1V+ITD8El9jlJdMP6Vxjv8r0Q365xrLJ9HczJzgKOAj4nLu/FthMXejrkScd1pua2RlmttzMlm/bvHmk9gohRK6EOME1wBp3vzVuX0XkFNeZ2XSA+P/64V7s7kvdfYG7LxgYPz4Pm4UQIjeChsiY2Y+Bd7r7SjO7EBj0Zhvc/SIzWwJMcffFWcfJKxzOGj6Tl0I1hJfYSaE6RiV2bZOXQjWEl9jlpVANnSmxKyocHtV8FwDeDXzVzEYDDwN/SfQt8kozWwisBt6St3FCCFE0QU7Q3VcAw3ngo3K1RgghOozK5oQQlaY0UlojoQiFaihmFjspVEeEzmKXl0I19E+OsGwK1dCZWewkpSWEEAUgJyiEqDR9EQ4nyU2hGoIncMpLoRr6Jzwu2/AZ6J8JnLKGz0B4eJyXQjV0ZgInhcNCCFEAcoJCiEojJyiEqDQ9mRP0OtfdLwrVUEKV6gIUqiFcOSYvhWpI5wj7VaEawkvj8lKohuwSu7wUqpUTFEKIApATFEJUGjlBIUSl6cmcYF5kKVRDeK4xL4VqSKtU96tCNYSX2OWlUA3pXGO/KlRDWqW6nxSq7/70IuUEhRAib+QEhRCVptLhcG4ldjkpVEN4iV0vlddBb5XY9bL6DISX2OWlUA3hJXYjUajWEBkhhCgAOUEhRKWRExRCVJpK5wTr6cQsdlKojgiexS4nhWoIzxFKoXp7mpXYhZbqjUShWjlBIYQoADlBIUSlkRMUQlQa5QQzCJ3FLi+ZfgjPEUqmf3syxxBCcImdZPojuj2LXf0YwhNfeY9ygkIIkTdygkKISlPpcLhsCtWQVqnuW4VqCFaOyUuhGtLhcb8qVENapboTCtWQXWKXl0L1lw6+TOGwEELkTVMnaGbzzGxF4u85MzvHzKaY2Y1mtir+v3MnDBZCiDxp6gTdfaW7z3f3+cDrgOeBa4AlwDJ33xtYFreFEKKnaCknaGZvBP7R3Q8zs5XAEe7+lJlNB25x93lZr580eqq/fupba+1WpKGSObCyqSbnpVANaZXqkShUH7/PvbXlVsq+6vNfZVNOblehGuqG4uSkUL3nv/4qta1d5eRNezQbZ9V5yqZQfdo+txeSExzVfJcUbwUGn8yp7v5UvLwWmDrcC8zsDOAMgJ0GJrZjoxBCFEbwDyNmNho4EfhG/TaPvk4O+5XS3Ze6+wJ3XzB6h7FtGyqEEEUQHA6b2UnAWe7+xrjdcjg8Z/+Jfv43D6q1+0UVJTeFagiuLmmmUN2vqihlU6geiSpK2SmiumQkCtVlUJE5haFQGOA64PR4+XTg2ryMEkKIThHkBM1sPHAMcHVi9UXAMWa2Cjg6bgshRE8R9MOIu28GdqlbtwE4qgijhBCiU3S1bK5fVVHKplDdz6oowQrVEJwjHIlCdSuqKL2UI2xWYhc6i91IFKrLkBMUQoi+Q05QCFFp5ASFEJWmNFJazaShWimxKxuhCtUQniMciUJ1Vo6wX8cQQniJXV4K1ZAtDdUvYwihmFns6scQ/vSaxcoJCiFE3sgJCiEqTUfD4UnzpvqhS4dUZNpVTi6banJeCtWQVqkeiUL1fX93aW25WUibVWJXOuXkNhWqIR0e56VQfdb//ji1LVQ1BtLh8c0/OaDh67pFqEI1ZJfY5aVQ/djCJQqHhRAib+QEhRCVRk5QCFFpOpoTNLNfA6uBXYGnO3bi5siebMpmD5TPJtmTTR72zHH33fIwJklHnWDtpGbLi0hwtovsyaZs9kD5bJI92ZTNniQKh4UQlUZOUAhRabrlBJd26byNkD3ZlM0eKJ9NsiebstlToys5QSGEKAsKh4UQlaajTtDMjjWzlWb2oJkt6eS5EzZ8yczWm9k9iXVTzOxGM1sV/9+5g/bMNrObzew+M7vXzM7upk1mtpOZ3WZmd8X2fChev4eZ3RrfuyviKVg7hpkNmNkvzOz6bttjZo+a2S/NbIWZLY/Xda0PxeefbGZXmdmvzOx+Mzu0i31oXnxtBv+eM7Nzun2NGtExJ2hmA8BngeOA/YBTzCxbh6gY/gs4tm7dEmCZu+8NLIvbnWIr8F533w84BDgrvi7dsukl4Eh3PxCYDxxrZocAHwU+4e57Ac8ACztkzyBnA/cn2t2254/dfX5i2Ec3+xDAJcD33H1f4ECia9UVm9x9ZXxt5gOvA54HrumWPU1x9478AYcC30+0zwfO79T562yZC9yTaK8EpsfL04GV3bArPv+1RDP7dd0mYBxwJ/AHRANdRw13Lztgxyyih+ZI4HrAumzPo8Cudeu6dr+AScAjxDn+MtiUsOGNwE/LYs9wf50Mh2cCjyfaa+J1ZWCquz8VL68FpnbDCDObC7wWuLWbNsWh5wpgPXAj8BCw0d23xrt0+t59ElgMDErw7NJlexz4gZndYWZnxOu62Yf2AH4N/GecMvh/8TS5ZejXb2VovvIy2LMd+mGkDo8+pjr+k7mZTQC+CZzj7s910yZ33+ZRKDMLOBjYt1PnrsfMTgDWu/sdTXfuHIe7+0FEqZ2zzOwNyY1d6EOjgIOAz7n7a4HN1IWa3ejXcZ72ROAb9du69ZwNRyed4BPA7ER7VryuDKwzs+kA8f/1nTy5me1I5AC/6u6DE9x31SYAd98I3EwUbk42s8F5qjt57w4DTjSzR4HLiULiS7poD+7+RPx/PVGu62C6e7/WAGvc/da4fRWRU+x2HzoOuNPd18XtbtszLJ10grcDe8e/6o0m+pp8XZPXdIrrgNPj5dOJ8nIdwcwM+CJwv7v/e7dtMrPdzGxyvDyWKD95P5Ez/NNO2+Pu57v7LHefS9Rnfujup3XLHjMbb2YTB5eJcl730MU+5O5rgcfNbF686ijgvm7aFHMKQ6EwJbBneDqcJD0eeIAox3RBN5KgRDflKWAL0SfoQqIc0zJgFXATMKWD9hxOFBbcDayI/47vlk3Aa4BfxPbcA3wwXr8ncBvwIFF4M6YL9+4I4Ppu2hOf9674797BftzNPhSffz6wPL5v3wJ27nK/Hg9sACYl1nX1GjX6U8WIEKLS6IcRIUSlkRMUQlQaOUEhRKWRExRCVBo5QSFEpZETFEJUGjlBIUSlkRMUQlSa/w+6fDw7+qEA/AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "C = Fs_cout.conj().T @ D @ Fs_cin\n",
    "plt.imshow(C.real.detach())\n",
    "plt.title(r'$\\bf{C}$: Block Matrix of Doubly-Block-Circulant Matrices')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "superior-arrow",
   "metadata": {},
   "source": [
    "This matrix `C` represents the circular convolution in the spatial domain.\n",
    "We reshape the input vector such that the $n^2$ pixels of the first channel come first,\n",
    "then the second channel, and so on:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "finnish-victory",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_matrix = (C.real @ x.reshape(batches, cin, n**2).permute(1, 2, 0).reshape(-1, batches)) \\\n",
    "                .reshape(cin, n, n, batches).permute(3, 0, 1, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "helpful-standard",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.3579970072896685e-06"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(y_matrix - y3).norm().item()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fallen-travel",
   "metadata": {},
   "source": [
    "We will now make the claim in the figure above explicit: applying the Cayley transform to `C` is equivalent to applying the Cayley transform to `D` and then applying the DFT matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "permanent-cricket",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert cin == cout, \"This section only applies for cin == cout, since we need square matrices\"\n",
    "cayley = lambda B: (torch.eye(cin * n**2) - B + B.conj().T) @ torch.inverse(torch.eye(cin * n**2) + B - B.conj().T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "lucky-martial",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.071995368169155e-06"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(cayley(C) - Fs_cin.conj().T @ cayley(D) @ Fs_cin).norm().item()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "mathematical-shame",
   "metadata": {},
   "source": [
    "What does an orthogonal convolution look like, anyways?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "artificial-amendment",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATEAAAExCAYAAAAUZZVoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4s0lEQVR4nO2deZQc1Xnof59mNKOZ0TpaRlsbCZAl8IKwZQwmzsFgYmwTvMaBOH48h3N4iR0/sPGCE8eJcRY7TmL7OV4OxzhxEj/whgOPOE6IAifeAgizGEvISEICCQktaKTROprRfX9Ma+p+d6Zvd/XcrqH7fL9zdKaqq7rqTnXNVX1ff/d3xTmHYRhGszJlshtgGIYxEawTMwyjqbFOzDCMpsY6McMwmhrrxAzDaGqsEzMMo6mZUCcmIpeJyEYR2SQiN6ZqlGEYRq1IvXViItIG/BK4FNgOPABc5Zxbn655hmEYcdon8N7zgE3OuS0AInIb8CagYifW2zvFlUrZKX+5u6/iwWVIrzuvpSeDVrdPy3Y+eUhvnHIiOM6U8ZdHTlp53ymDlfetdhwZ9vZtC84R7Ntx8OTo8uBMfeC24977gmswdb9u4Ik5Hdn5Z+mLeXJAv9m/1m3H9X9qrj1r4HCnPqf/e4XrQ7P0caYckYr7htfEv37hOcZc60h7xsQZXpOGp+v2Te2v/HuGn5E6RXB/Dc84qdbbB7JGhPetvy7Bs0T70WBf7xqJPgXDs/QLcjhykfzzhPd7+DnkeL5RbQp/F++emr7ksNr23KHpar3tSLYcXsvBJ5/Z65ybH557Ip3YEuBpb3078MrYG0qldv7t+/NG11/z5Q+MLk8JOq2Ofn0lBmdnV/zoAv3LzVn13Ojysf+ap7Z17dbHGerKlk/06E/x5FTdhuFp2fL0p/RxTmZ9xNjjdKhVOvqz5cHZelv4B1K6e2B0+elLZ6htM7dlv/exOfpGXXL7VrX+zFuWjS5Pe+Ozatuhe/R/Hp37s99t5pO6Mxycnd0i/Wfqu9x/H0DHwWz92dfr48x4aJpabz+c7TsUXL8hb9eOA2oTw1163b9+nQeCz3qaPu6UoWx7/wXH1bZF/y/78MPfM9Zxdu/S59z/mmNqvfc/s1/m6ALdnuO92XvDDnjeI/q4x+Zk7516SO/bf7nuGDrWZR1D+HeF/6cT/F6Ds/T6mP8UKh0HaPM+7rBjn705a8SFN/232nbbTy5Q670PZ43qv0hfyyd/6w+3jdeUhif2ReRaEVknIuv27TtZ/Q2GYRg5mEgntgMoeetLy68pnHM3O+fWOOfWzJ1rX4YahpGWiST22xlJ7F/CSOf1APBbzrlfVHpP16KSO/3qLIRc+77PjC6/5ssfUvvGwsswdj94RrbNDy0hHl4OBaFJLLwc1tGQCi+PLNTva9ORigov/dASxoaXy76+NWvfkrlqmx9e+qEljA0vezdkjdi/Sid5YuFlGCI+9+rsOIvv1PF2LLz0Q8vx2nfSa5IfWoIOL4eC694ZpBmGurN9w9A8Fl4OztT7HluZhS5+aAnx8PJEkPubtVEf1w8v/dASdHh54lwdI7otPWrdDy/90BLGhpcHz8i2t+tIU/9dhYFRJId4YgZxvGO1Bfnjo6uzBN+S2/S1jYWXfmgJ8NDNNzzonFsTnrrunJhzbkhEfh/4N6AN+FqsAzMMw2gEE0ns45z7PvD9RG0xDMPITd3hZD10zy+5lW97/+i6/42jH1pCPLyMfePoh5YQDy8791X+xhF0eBn75jIMH8Owxt8e++YSdKgcfuPoh5exby4BhrqyRoTfOMbCy/Cby+lPZ8cNv3GMhZdhWHqsV18U/1vHk0FJQ+yby7D0xv/G0Q8tIR5eHllQuW7CDy0hHl6GpRH+N46gw8vYN5fhN3rhN45+eBn75hJg2A+bg28c/fAy+s0l+l4N2xcNL4PjdHt/r+E3jrHwMvzmctt7PjRuOGmZdsMwmhrrxAzDaGqsEzMMo6kpNCcWllj4ZRN+fgziObIwlo+VTcRyZAMP6vKLXFX5XigfVpXHqvJj5RcAM56qrSo/Vn4B0OblMML8VCxHFpZfDP7zgtHlsGwiliMLyxLC3ztWlR8rvwhzUH7ZhJ8fg3iOLMxxxqryYzmywRn6jbGq/Fj5xawf6/KLWFV+rPwC4HCfP15OH8fPkUXLL4AT3dlyWDaRJ0d2bGFW+h+WTcRyZGH5xadX3245McMwWg/rxAzDaGoKDyeXXTP+oO/YgG/Q4eUbH3m32uaXTcTKL0CHl537dB8eq8qPll9UKZvww8tY+QVog0KsKj9WfgFw8BOZDiA24Bt0eBmWX5zwJAOxAd+gw8uw/OJor77Wsar8WPnFyWC0hl82ERvwDTq8nPGUHt3sh7+x0BJ0eClT9DliVfmx8ouOw/GyCT+8jJVfAHTt9uPmoPHe7xYrv4Dg7yEy4Bt0eBmGlm9+649Gl2MDvkGHl2H5xY/u/LCFk4ZhtB7WiRmG0dRYJ2YYRlNTaE5s5vQl7pUv/l+j63msDH4O6O/++9tqW6z8ol65ItSea0slVwyPNRG5Yr1WhpBOLz+VSq4YtqEIuSLoXFsqueKM7Tq3lsfK4Od/Zge5tHrliqBlhoXIFSGaa0slV3zssx+wnJhhGK2HdWKGYTQ11okZhtHUFJoTm9XR517Vd+Xoeh61jJ8jC/MZjTDEQu0TkKQyxIK2xE7EEOvnePKoZcL82ML70htiIT5EqRGGWNA5slSG2EPn62mJ8qhl/BzZgRX6nPUaYkFbYgsxxEJ0ApJUhtj1f2k5McMwWhDrxAzDaGqKNbv2ldyZV2bDjvJYGfzw8kRgJ2iEIRZqn4AklSEWgglyJ2CIPT6HisTCyzC069pbuWyiXkMsxIcoNcIQCzq8TGWIDcljZfDDy+6n9XWv1xALujShCEMsxCcgSWWI3finFk4ahtGCWCdmGEZTY52YYRhNTaE5sZ65JfeiN14/up5HLePnyNqCuLoRhliofRalVIZY0DmeiRhi/a+186hlwvKLg8sr56fqNcRCfGhRIwyxoHNkqQyxh5frMTV51DJ+juzuL1yottVriAVtiS3CEAvxWZRSGWI33mQ5McMwWhDrxAzDaGqsEzMMo6kpNCc2bbHWU/u1XzH1Dugc2fbX6sC5EZprqH0WpVSa63D7RDTXM7dluZo86p2whqxvfpZoSqW5Dt9bhOYadI4sleb68t/+kdoWG1oUqyHb/ps6OVSv5hq06roIzTXEZ1FKpbmuOycmIl8Tkd0i8pj3Wq+I3C0iT5R/RkorDcMwGkct4eTfA5cFr90IrHXOrQDWltcNwzAKp6ZwUkSWAXc5515cXt8IXOSc2ykii4B7nXMrqx2na2HJnf4/snCyXitD6e4Bta0RhlioPUxNZYiF2ocoVTPE1mtlCEOn/guyODWVIRbiQ5QaYYgFHaamMsQe15VAY4lYGfzQafZmXZdQryEWtCW2CEMsxMPUVIbYX/5x2hKLPufczvLyLqAvtrNhGEajmPC3k27kUa7i45yIXCsi60Rk3fDRw5V2MwzDqIt6O7Fny2Ek5Z+7K+3onLvZObfGObemraun0m6GYRh10V59l3G5E7ga+FT55x21vMlN0TNc+2UBYQ4sliNr37FPbSvdnS1XM8TueOuy0eVLvqBzYPdEhiiFpRp+jmzqgN62n161PudX944uH42UcYyQrYflFz6HXqC3hTmydu844bCZ47P0eqy8YNrGLIey89dDhY9uU8wQO+dxXQOy38tATI8YYg+ik31jDbHjn388+v41O9ax4Pt0Pw8WKnwgUD95uZqpOj07dtiMdx2Gw6FhHjuu1CUWP/74+Wr9ypt+Orp8G0EZx706aXfAy07PuUdve+7iyobY8Pf0DbF7iQ9R8nNks+/S+x48I1vOY4gNS2QqUUuJxa3AT4GVIrJdRK5hpPO6VESeAF5bXjcMwyicqk9izrmrKmy6JHFbDMMwcmPDjgzDaGoKHXbUtUgPO8qjlvFzZOEsRY3QXEPtsyil0lyDVl1PRHM9LUvD5VLLhDVkRxZUzsvVq7mG+BClRmiuQdeRpdJcnwximXr1y93BPVSv5hq06roIzTXEZ1FKpbl+7HOm4jEMowWxTswwjKam2NmOFpTcindk4WQeK4MfWs14qvLQolSG2JH21TaLUipDLGhL7EQMsX74lsfKEH4O/nFTGWIhPrSoEYZY0OFlKkPs0fn6guWxMvgcW5jGEAvaEluEIRbisyilMsQ+9jcWThqG0YJYJ2YYRlNjnZhhGE3NpJZY+GUTMfUO6BzGcJCjaIQhdqR9tc2ilMoQC9oSOxFDrDpmRL0D8dmvZzyV3hALunSiCEMs6BxZKkNs2J6YuTRWfvHmt6YxxIK2xBZhiIX4LEqpDLHr/9JyYoZhtCDWiRmG0dQUG06GZlfvKTgWOkGV8oIGGGKh9jA1lSEWtCV2IoZYv0o6VnEO8fKCYe+4qQyxEK/ub4QhFnSYmsoQ2/PzoMo9Zi6NhE6iKyzGUqMhFrQltghDLMTD1FSG2E3v+LiFk4ZhtB7WiRmG0dRYJ2YYRlNT+LCjF/7G+0fX81gZ/BzZkYXBUA8v3RKzX4DOkS37+la1LTZEKTaL0nCQj8ozRCnMkfmGjtgQpVgZB2hDRx4rQ1he0NlfufwiliMLyzgGZ+p9Y0OUYobYhfcFhtjIEKVYjmysIbZyji6WI+tfqRuYx8qgHiGCP8Na7RcwNkd2dPXR0eXYEKVYGQfAgRXZcrUhSrFZlHxDbMx+ATpHFpZxrPuHGywnZhhG62GdmGEYTY11YoZhNDXF5sTml9yqN2c5sTxqGT9HFmpUGmGIBZ0jK8IQC7UPUapmiD3p/W551DJhDZl47WtmQyzoHFkqQ2x4DfKoZfzcVnif1muIBW2JLcIQCzpH1ihD7CNfspyYYRgtiHVihmE0NYWGk9OWlNwLfjcbdpTHyuCHl2F41AhDLNQ+AUkqQyzUPgFJNUPs4NosvMxjZQjLL8RbTWWIhdoNGKkMsaDDy2SG2KNoclgZ/M/6RLfeVq8hFrQltghDLMQnIElliL3v1g9aOGkYRuthnZhhGE2NdWKGYTQ1xZtdfyfLifllEzH1DgRDlCJlE6kMsVD7LEqpDLGgLbETMcTO+qU3lCii3oH4xLF+qUYqQyzEhxb5pDLEgs6RpTLEDpSCOp3I0KJY+UU4zK5eQyxoS2wRhliIz6KUyhC7+UYrsTAMowWp2omJSElE7hGR9SLyCxG5rvx6r4jcLSJPlH/OaXxzDcMwNLU8iQ0BNzjnzgbOB94rImcDNwJrnXMrgLXldcMwjELJnRMTkTuAvy3/u8g5t1NEFgH3OudWxt7bWSq5Je+/fnQ9j1rGzwHFashSaa6h9lxbKs01aNX1RDTXs7ZUnqUoppYJ8z8zHsouaCrNNcSHKDWT5nr7JfqmyaOW8fM/Het0kq5ezTVUUV03QHMN8VxbKs31lg8lyImJyDLgXOA+oM85t7O8aRfQV+E914rIOhFZN3z48Hi7GIZh1E3NnZiITAe+C1zvnDvob3Mjj3PjPtI55252zq1xzq1p6+kZbxfDMIy6aa++C4jIVEY6sG84524vv/ysiCzywsnd1Y7Tdhxmbs4eF/fTO7o851f3qn2PRsPLwGIamTj20Av0Nj+8HO7V+4ZDlPwQMmaIbd+xT20r3a33jRlid7x1mVq/5AtZCHlPZIhSGD6G4eUzV2SxweI7dXvCsMun7191eHTM+7omDB/D8BKy7WGI0xm0zy/BOD5LHzdmiJ22UYcfO3/dH6KkzxkLL+c8HhhivUBiehVD7EGyaySn6+hiL7VbGWbfle178Azd9jyG2DB09/88xgxR8q7BcJA+CdlxZXYP/fjj56ttV970U7V+G94sSvfqz+iAl2Sac4/e9tzFsfCy8t+1Ty3fTgpwC7DBOfc33qY7gavLy1cDd9R0RsMwjITU8iR2IfAu4Oci8nD5tT8APgV8S0SuAbYB72hICw3DMCJU7cSccz+i8nPdJWmbYxiGkY9Chx31zCu5s67IzK551DJ+CUbnvhxlE3UaYsPtRRhiQVtiJ2KI9c2ledQyYTnBsd6sDakMsSPtrW0WpVSGWNC/ZypDbPfOoAQkh1rGz5ENh7ND1WmIBX2vFmGIhfgsSqkMsT+7xYYdGYbRglgnZhhGU1O4xeL0qzOLRR4rgx9eDjyoyy8aYYiF2icgSWWIBW2JnYghdpoXcuexMoRlCf7vncoQO3KsylX5jTDEgg4vUxlij9+pDbF5rAx+eHm4L/hF6jTEgrbEFmGIhfgEJKkMsVvf+YcWThqG0XpYJ2YYRlNjnZhhGE1NsSUWKxa5VZ//ndF1v2wiZq0AnSPr3Kf73kYYYqH2WZRSGWJBW2InYoi9+ItZjixmrYD4xLFHe7O2pzLEQnxoUSMMsaBzZKkMsV17w8+6srk0Vn7RtTtM/gWNr9EQC4HJpQBDLMRnUUpliH3oq1ZiYRhGC2KdmGEYTY11YoZhNDWF5sRmtc1z50+/YnR9y1eXjy733KNraFx7oGfpz4L7tuO6zYcXZXmRwSq1L0PTs/dO2xvobELDqJePORmMMvXzSsfm6uOENWXtR706tvbKdWsAA6dnOZ+FP9H7Hl5U+f+cqQO67etu+vLo8pqP/57aNtxZeYhLWOd06AXZdZ/3kN62f1X4GWXrL37rBrVt81dW6Td7bw3zn4eXZMs92/W2E9ODnJNXjxbmV6btCZRNM7PtR1cfUduW/lP2QTz7Cv2hTNXiXo4sya7Jyr/aqredU1Lrx2dn9+bul+v2Lb0nS3oef58eZnfwh1r/I15uq+cZ/XsNBLopP+8V3hdOPAWSTscy3KH39f8+Di3XybWZm/S96H+G4T3t596Onauvu2zVH/68R72hTmfqc2z85AcsJ2YYRuthnZhhGE1NsSUWc0vuRW+8fnR9zh2/GF32Q0uIh5eHl+g2z/9Z9qjrh5YQDy+7nw2G9cyqHF6GX/X74aUEX2OHJRd+eOmHliPHCUPRbPtzL9UH9sPLWGgJerKIh2/8ktoWCy/Dr+/PeWMWFt7/Uz0PTCy89ENLGBvWtB/1G6uP44cmfmgJMPVQZTvG1GCoUyy8PLgiuO8XZSUCfmgJ8fCy93Fdp9Dz0NNq3Q8v/dASdHjZcSCY6OVCbc7ww8vwfgvDy/4V2XLHQOXPwQ8tYWx46d9DYRlHLLwM0wP+30OYaomFl35oCXDfrR+0cNIwjNbDOjHDMJoa68QMw2hqCs2JdfeV3IrfzFQ8ftmEnx+DeI4sLGnwyyb8/BjEc2QqL8PY4Th+jixWfhFOVjrGJhvJCYQ5Mv+4YR7Jz5FVK7/w8xJh2UQsRxaWX/hDWMKyiViOLCy/CFU8fulE+DnEyi/C4V9+2URoj43lyMIcj1824efHIJ4jOzFdn8MvmwCdI4uVXwx36F8sLJvwc2Sx8gsIJpMOyib8HFms/AJgyEtLh/d0LEcWll/4Q7PCsolYjiwsv9j8URt2ZBhGC2KdmGEYTY11YoZhNDXF5sQWlNwL357NduTXfvn5MYjnyNoenV7xHEPTK9eQwdgcmU9saFGshiyc7Sg2tChWQwZEJz32cxixGjKAo3Oz/5/C2q9YjiysIRs4LVsOa79iObKwhiwcLhQbWhSrIRs4rXLtl58fg3iO7MjCsH4qW1b5MYjmyAZnVq79Ap0ji9WQ7Vmtk06xoUWxGjLQSqnY0KJYDRmAa8u2h7nJWI4srCFrO5IdJ6z9iuXIwhqyLVd9zHJihmG0HtaJGYbR1BQaTk47Y4krfep3R9f9somYtQJ0eLnz3S9R22JDi2Lh5fFZQVlCZGhRrPwiLLGIDS2KlV+AtsLmsVaE4WXvo9l7Y9YK0OFlWH5x/oeyzytmrQAdXoblF3MfDb6+jwwtipVfDHXrdb9sImatAB1ehsOF/LKJmLUCUOFl6R/0BxobWhQrv9j67tN1WyNDi2LlF6BnYIoNLYqVXwAcWVy5bCIWXo4pX3ll5bKJWHgZplrW/6VZLAzDaEGqdmIiMk1E7heRR0TkFyLyifLry0XkPhHZJCLfFJGOascyDMNITS1PYseBi51z5wCrgctE5Hzg08BnnXNnAvuBaxrWSsMwjArkyomJSDfwI+D3gH8BFjrnhkTkAuBPnHOvi72/8wUlt+jD142u51HL+Dmg+d98TG1rhCEWas+1pTLEgh5SNRFD7OHzsjxEHrVMmP/ZdNVXRpdTGWIhPkSpEYZY0Lm2VIZYgj+fPGoZP/+z/JYtalu9hljQltgiDLEQz7WlMsQ+/mcTyImJSJuIPAzsBu4GNgP9zrlTV287sKTC2w3DMBpGTZ2Yc27YObcaWAqcB6yKvyNDRK4VkXUism740KHqbzAMw8hB7hILEfk4cBT4CDnDye75JbfqLVnFfh4rg0/XnsrlF6kMsSNtqG0CklSGWNCP+xMxxA6cXp+VISwv8CfwTWWIhXh1fyMMsaDDy1SG2GPz4pPExKwMfnjZcVDX6dRriAVtiS3CEAvxCUhSGWLrLrEQkfkiMru83AVcCmwA7gHeXt7tauCOascyDMNITXv1XVgEfF1E2hjp9L7lnLtLRNYDt4nInwIPAbc0sJ2GYRjjUrUTc849Cpw7zutbGMmPGYZhTBqFDjvqWlhyZ74zM7vmsTL4ObIpesRIQwyxUPssSqkMsaDzBxMxxJ7oyY6bx8oQll/41zqVIRbiQ5QaYYgFnWNJZYid/YS+tnmsDH6OLPys6zXEgrbEFmGIhfgsSqkMsY9+wYYdGYbRglgnZhhGU2OdmGEYTc2k5sR8qqll/BzZsdk65m6EIRbyzaLkU68hFrQldiKGWDUEKIdaJqwhO7Dcs5gmMsRCfGhRQwyxoHJkqQyx4e+VRy3j58jm36mThvUaYkFbYoswxEKVWZQSGWK3XmczgBuG0YJYJ2YYRlNTaDg5fW7Jvfh114+u12sundqn44RWMcRCMERjAobYeT/PaiNi1gqITxw7+z+z5/9UhliIDy1qhCEWdHiZyhB7dE3lsgnQ4WWs/ML/vKB+QyxoS2wRhliIT0CSyhD72OesxMIwjBbEOjHDMJoa68QMw2hqii2xWFRyy/9nVmKRRy3jx8oDp2ttSSMMsVB7ri2VIRa0BXMihthBL4+TRy0T5n+Ork5viIX4EKVGGGJB59pSGWKHuvQ58qhl/Hs6HEpXryEWtCW2CEMsxHNtqQyxmz52g+XEDMNoPawTMwyjqbFOzDCMpqbYnFifHnaURy3j58jCHFQjNNdQ+yxKqTTXoFXXE9FcO289j1omrJE6uMJbT6S5hvgQpUZorkHnyFJprtuO6Z3zqGX841YbRlar5hq06roIzTXEZ1FKpbneeJPViRmG0YJYJ2YYRlNTaDjZvaDkVvxmFk7msTL44WUYSjXCEAu1z6KUyhALOqyZiCF2ymB9Voaw/MIPj1IZYiE+RKkRhljQ4WUqQ+xgcG3zWBn88LJrbxpDbNiGIgyxEJ9FKZUh9tEvWjhpGEYLYp2YYRhNjXVihmE0NYXnxFa+LZsBPI9aRuXIghyFTypDLNQ+i9Lz0RB7ZFF9apmw/OLIQr8sQbenXkMsxIcWNcIQCzpHlsoQG+bE8qhl/BxZ2xF9nHoNsaAtsUUYYiE+i1IqQ+wjX7JhR4ZhtCDWiRmG0dRMasW+XzYRs1aADi9nbgtDsvSGWKh9ApJUhljQltiJGGKnHqpcNhELL8Pyi97H0xtiIV6V3whDLOjwMpUhdu/L9bWMVeXHyi+OvDKNIRa0JbYIQyzEJyBJZYhd/2krsTAMowWpuRMTkTYReUhE7iqvLxeR+0Rkk4h8U0Q6GtdMwzCM8cnzJHYd4D+ffxr4rHPuTGA/cE3KhhmGYdRCTTkxEVkKfB34M+ADwK8De4CFzrkhEbkA+BPn3Otix+laWHJnvCvLieWxMvg5oHCYTyMMsVB7ri2VIRa0oWMihtg5672VHFaGMP/jlyKkMsRCfIhSIwyxoHNtqQyxP7/jrKDttVsZ/PyPb/SF+g2xoId4FWGIhXiuLZUhduMnJ5YT+xzwYeDUKecC/c65U63dDiwZ532GYRgNpWonJiKXA7udcw/WcwIRuVZE1onIuuEjh6u/wTAMIwft1XfhQuAKEXkDMA2YCXwemC0i7eWnsaXAjvHe7Jy7GbgZRsLJJK02DMMok6tOTEQuAj7onLtcRL4NfNc5d5uIfAV41Dn3pdj7u+eX3Kq3ZMOO8qhl/BxZGEc3whALtc+ilMoQC9oSOxFDbJeXn8qjlglrpKYeqlxDVrchFqJDlBphiAWdI0tliD0+u361jJ8jE52erdsQC/G/h0YYYiE+i1IqQ+y6r6cfdvQR4AMisomRHNktEziWYRhGXdQSTo7inLsXuLe8vAU4L32TDMMwaqfYYUcLS+6M3/ZKLHJYGfzwMnz0boQhFmqfgCSVIRZ0+chEDLEnO7P35rEyhKGJ87alMsRCfIhSIwyxoMPLVIbYtsHAPJvDyuCHl0P6NqjbEAvaEluEIRbiE5CkMsRu+AsbdmQYRgtinZhhGE2NdWKGYTQ1xZpd+0ruzKu82Y5yqGX8HJkLZKiNMMRC7bMopTLEQtwSm8cQq4Yv5VDLhNdk4LT0hliIDy1qhCEWdI4slSHWVydBPrWMKr9oS2OIBW2JLcIQC/FZlFIZYh//M8uJGYbRglgnZhhGU2OdmGEYTU2xObFg2JFfpxVT74DOkY0ZotEimmvQquuJaK6HuiqrW2I5srCGbKg7W06luYb40KJGaK5B58hSaa5nbg7yoZGhRbEasiOL47VftWquQauui9BcQ3wWpVSa6/sv+7TlxAzDaD2sEzMMo6kpNJyc3ltyL7n0utH1PFYGP3wKLZiNMMRC7WFqKkMsaEvsRAyxe1dny3msDGHotPkrq7KVRIZYiA9RaoQhFnSYmsoQ2/X6+q0MfujUNhico05DLOi/jyIMsRAPU1MZYn/4Lx+xcNIwjNbDOjHDMJoa68QMw2hqCs2JTVtackvfl5VY5FHL+HQ/W7n8IpUhFmqfRSmVIRa0JXYihtizrqpPLROWF/i5mVSGWIgPUWqEIRZ0jiyVIdYFH2cetYyfI+tfobfVa4gFXYJUhCEW4rMopTLEPvnB9GZXwzCMScc6McMwmppiw8nFJbfsmsxikcfK4IeXYVjTCEMs1D4BSSpDLOivtSdiiB3qqc/KEJZfiPerpDLEQry6vxGGWNDhZSpDbMdA5bIJiFsZ/PByzL1XpyEWtCW2CEMsxCcgSWWIve/WD1o4aRhG62GdmGEYTY11YoZhNDWTmhPzqWZlCHNkPo0wxEKOWZQSGWJBW2InYojll1liJI+VISy/8IcLpTLEQnxoUSMMsaBzZKkMscfm6fbksTL4ObKOfn3Oeg2xoC2xRRhiIT6LUipD7H/9qw07MgyjBbFOzDCMpsY6McMwmppJne2oXnNpOMNMIwyxI8eqbRalVIZY0JbYiRhi91/xomw5ot6B+OzXcx+tnF+p1xAL8aFFDTHEgsqRpTLEvuhv36O2xYYWxWrIjt+5QG2r1xAL2hJbhCEW4rMopTLE/uCZvx03J9YevjAeIrIVGACGgSHn3BoR6QW+CSwDtgLvcM7tr+V4hmEYqcgTTr7GObfa6wlvBNY651YAa8vrhmEYhVLTk1gF3gRcVF7+OnAv8JHYG9qOOnrXZ8+aRxdkz7PPnVX5a3aAmRsPjC4/8a45atvpt2fPwfteouOWg2fq4yy/8Sejy9s/+irdPv2NPMPTvOUOvW2q9+h9YLn+vyAMP/xSjtBw2vO0frw+diQLIec+pmOygWVZg4am6eNs/POz1Xrn3qxNs162R23b94T+/v6x288aXV4ahGuHFme3yIL/7tdtXajD3WdfkV2k5f+k6y+euXypWu/IPk4OLVObmL41C1WOz9G/58HT9C3bvcsr1QiubddO/bmcOJyFROF9sfiH2e89sFSfIwyNV96SlVi84vL1atsD95yl1rv2ZjfylCH9Wc+4J8unHLg4CKv+Q8d2+1dl9/XRvri9uHNf9nv37NI5kz3nZNu6dwV1MEFmqe2J7JxDPTpODu0dvhUkTB089brs955yIiwPqVyf1PPtYGzYhePvV+uTmAP+XUQeFJFry6/1Oed2lpd3AX3jv9UwDKNx1Pok9ivOuR0isgC4W0Qe9zc655yIjPsNQbnTuxagc9rsibTVMAxjDDU9iTnndpR/7ga+B5wHPCsiiwDKP3dXeO/Nzrk1zrk1U6f2jLeLYRhG3VQtsRCRHmCKc26gvHw3cBNwCbDPOfcpEbkR6HXOfTh2rFndi935L7xmdP34gqxT8/NjMDZH5rPsj36q1jd/5oLR5dNv10MewhzZQDaPJ93P6HOcCPrYzv2ebqdH7+vnyEJLaJij8HNk4VCnkIEzshxG22H9f8yyu7LkjJ8fg7E5suk7s0Y9u0Zf2xnnVc6R9Tytz1n6bpbbevZSXVMx//5+te7nyPz8GIwtZ+nwyiHC/MqRxdl1n7413KbXu3dmy2F5Q5gj83VAg7OCoTreJVr8Q/2BxnJkszbrhNmCz+gG+zmyRT/Wx/XLDQ4u0xchvE/8HJmfH4OxObKpB71zBMOXejdkH4SfH4OxObJj87z2BJ9RmCOb/7NseeAFeuezL984uvzgA1phG+bI2j0DcN86nROrNNtRLeFkH/A9GalbaQf+r3PuByLyAPAtEbkG2Aa8o4ZjGYZhJKVqJ+ac2wKcM87r+xh5GjMMw5g0bNiRYRhNTbEqniUlV3pPNtvR6d/OCvz9/BjEc2QuGC6x/MYsR+bnxyCeIwtnAJ8aDLvwc2R+fgx0jqya5trPkcVqyEa2e9qecwI1ipcj8/NjMDZHdqQv23fMUJ1IjiysIfPrzcLar1iOLKwhC9U3fv1XRzBcyM+/+PmxkfYEecyZ2XY/PwbxHFlXMNTJn20ovL9iObJQTz1jmy429HNkYQ2ZnyPb+1J90jAH5d8nsRoyCHTpwZ+3nyPz82MwNkc29XB2nFBFFcuR+fkxgM7+7DwL/mCL2hbLkbUHM2T98hMfMBWPYRith3VihmE0NZMaTvr4oSXEw8vB4Ktzv2zCDy0hHl7uepU+R2xoUaz8Ynha5XIQ0OFlrPwCYNr+7LE8LJvww8tY+QXAlrdn4WVYNhELL8Pyi57/kykdwrKJWHgZll9svmq2WvdLJ8KhRbHyizBU98sm/NASqoeXPn7ZRDiRbTS8DMO1wNLgh5ex8ouZm/RxwrIJ/zrEyi8Atr0hCy9jQ4ti5RcAwx3Ze8OyiVh4GZZfTD2YbQzLJmLhZVh+seVDNnmuYRgtiHVihmE0NdaJGYbR1BSaE+uZW3Ivfv31o+uxoUWxHFn7Yf2Vt1824efHIJ4jm7lZ7xsbWhQrv6hWNuETK78AWPzDLGcXG1oUK78AmPOSvaPLYdlELEcWll8s++fsc4ipd0DnyMLyi6Gu0D5aeWhRrPxiyonAAuvlR8fOMl45RzbzqcplE6F6J5YjC8trYkOLYuUXu284TW2LDS2KlV9AjqFFkfILgMFZlcsmYjmysH3DXdmJwrKJWI4sLL/Yer3NAG4YRgtinZhhGE3NRMyuuRnqgr0vzR4nY1aGPedpe6sfPi27Qx+3w6tyl+Br2a2f1CUWK/4xC49ihlioPUwNDbGx6v4xJRZBmPrMq7NzxsLUGZt1DBGGqf0/zxQEpTEV5/oaHZ2bNepku2771jdn1ygsvzgehKlPvy2zt4ZlHL4hFuCMW/tHl8eEqfOyMHXxXXFD7BTvNGF1fyxMPRjc+n75xdH54YTL+jh+mHpihj7njosC44UXPnX2BxP2nsw+wwV/vU1t2zamur/yZBthmLr94uweOu37tRtix5hGDniTfwzqjSd6glA9Eqb65zkxU6dBwt9l/V3ZRDUdQRlHJexJzDCMpsY6McMwmhrrxAzDaGoKLbHofEHJLfrIdaPreawMfnnB/Psrl1+kMsSCHqJUhCEWtCV2IobYY3Ozc+axMoTlBV370htiIT5EqRGGWNA5slSG2LPevFFty2Nl8MsLOg4EpQZ1GmJBW2KLMMRCfBalVIbYDX9hFgvDMFoQ68QMw2hqrBMzDKOpKTQn1j2/5Fa9JVPx5FHL+Dmy/St1zN0IQyzUPotSKkMs6JzPRAyx057LzpNHLRPO7OPng1IZYiE+RKkRhljQObJUhtjeDfWrZfwc2exN+iT1GmJBW2KLMMRCfKhTKkPs/d+wYUeGYbQg1okZhtHUFG6xeNEbrx9dz2Nl8MPLp17XVfEcqQyxUPsEJKkMsRC3xOYxxLYNeuaAHFaGsPzi6Uuz+COVIRZ06UQRhljQoUsqQ+yRvsplExC3Mvjh5ZJ742UTtRpiQVtiizDEQnwCklSG2CdvsHDSMIwWxDoxwzCaGuvEDMNoagrNic2cscS94tz3jK77ZROxmX1A58jmParb3AhDLNQ+i1IqQyzUPrSomiF2/qNZbiYsm4jlyMLyi9kNMMRCfGhRIwyxoHNkqQyxA8v1vrGhRbHyi/Yj+n31GmJBW2ILMcRCdBalVIbYJz5uw44Mw2hBaurERGS2iHxHRB4XkQ0icoGI9IrI3SLyRPnnnOpHMgzDSEutT2KfB37gnFsFnANsAG4E1jrnVgBry+uGYRiFUjUnJiKzgIeB0523s4hsBC5yzu0UkUXAvc65lRUOA0DXopJb9jsfGF3Po5bxc0C+ZgbimutYrm3ZHTrvkEfj4+faxmquKw9RqjYbk6+6jg1RitWigc755FHLhPmfk955YgofiM+iPWVInzM2RMnPtVXTXMeGKMVybaHm2s+/5NFcD+vbLZdaxs//DM7W74sNUao2i3bHlCyXFRuilEdzXW2IUq25tpjCB+Ka60e/WH9ObDmwB/g7EXlIRL4qIj1An3PuVIpzF9A33ptF5FoRWSci64aOHB5vF8MwjLqppRNrB14GfNk5dy5wmCB0LD+hjftI55y72Tm3xjm3pr27Z7xdDMMw6qaW2Y62A9udc/eV17/DSCf2rIgs8sLJ3VWP5LS59MT0bLmalcGf1WY4eGTfenk2DKmaIbb3kazf7twdPhnqTraXLLwMQ8stv5GFkGd8qHZDLOjH8Cc/pfftfiZbDksEdHgZt8kOeyOzwmEzUw9VLi+YelAfd3BWtu8zr9a3Syy8nP2EPmfXvmCIEvNHl+eGhlj88FKH+HkMsX0P6DY8+dtZCBkzxIam3kPLwvAyWw+H9bQHw+X2vMwvL9D3rR9evvwV+oKFQ5SGvL+VsHyFP9c5Ct8S+4rPbFDbHqCyIXbgNP234l+j7a/VN1gYXvr3dXhNOvqz5efO0n/Y8x+p3RBbiapPYs65XcDTInIq33UJsB64E7i6/NrVwB3jvN0wDKOh1Drv5PuAb4hIB7AFeDcjHeC3ROQaYBvwjsY00TAMozI1dWLOuYeBMd8KMPJUZhiGMWkUOuyoq6/kznxnVmKRRy3j58jCnE4jDLFQ+yxKqQyxoC2xEzHE1quWCcsvuvakN8RCfIhSIwyxoEswUhliO/qDIWY51DJ++UVnv84N1WuIBW2JLcIQC/FZlFIZYjd97AYbdmQYRuthnZhhGE3NpIaTPtWsDH54OW2/fixvhCEWap+AJJUhFrQldiKGWD+8zGNlCKv7fVIZYiFeld8IQyzo8DKVIbZrb+WyCYhbGfzwcupBvbFeQyxoS2wRhliIT0CSyhD7029/yMJJwzBaD+vEDMNoaqwTMwyjqSl28ty+kjvzyiwnVq+5dPEPdR6pVQyxoC2xEzHE+pPnxqwVEJ84duZTlcsm6jXEgi6dKMIQCzpHlsoQO7CsctkE6BxZrPxiuEvfi/UaYkFbYoswxEJ8FqVUhtiNN5nZ1TCMFsQ6McMwmhrrxAzDaGpqHQCeBDdF63f82qYw/zOmTszLAT3zah2fxwyxMzbrF/xc296X6n1jQ5T2nKftrdoQq4/TEQyLEq+OZ+sndQ5sxT/qHI9viT39dp34ymOI3f37mSF2ii45ippLw/zPQe8WCWvIjs4PckUHsuWwFq3/5/PUekkNUdLHOTo3O+fJdn2crW/Wn4NfR3Y8GKL09Nu0vdXPtYWG2DNu7R9dHmOInadzbYvvynJtT75Ln6NtUNc2nejx8o2R/E+YYzoxU19sP5e1/i4tUO4Icm2+JXbHRUEe08u1dfbrfOPgSf23suCvt40ubxszREnfVLFcm2+IPe37tRtiQ11SJexJzDCMpsY6McMwmpriSyyuykos8lgZ/PKCKfpbYxV6VjPE+gycoZ9XY0OUYhOQzL+/cvkF6CFKsTIOgGV/VNsQpVgZB8DUA/VZGcLygiOLs+WY/QLiE9D6hliID1HySznCMo6ufXrf2BCl2AQkeQyxsQlIOp/T7ctjZfDDy2M62o4OUao2Ae3Zl28cXY4NUYqVcYA2xFYbohSbgOTgsqx9MfsF6PAyDLE3/IWVWBiG0YJYJ2YYRlNjnZhhGE1NsSqecPLcOtUy4VevjTDEgs6RFWGIBZ0rmogh1s+p5FHLhOUXnXuz96YyxEJ8iFIjDLGgc2SpDLH9Zwa+pBxqGT9HNvVwGkMsaEtsEYZY0DmyRhlif/Y1M7sahtGCWCdmGEZTY52YYRhNTbE5sYUld8a7spxYHrWMnyMbnla51iqV5hpqn0UpleYatOp6IprrrW/K8gl51DJhzsK/nqk016Drv4rQXIPOkaXSXG95+2y1LY9axs+RDXcEivM6NdegVddFaK4hPotSKs312p/8keXEDMNoPawTMwyjqSl22NGCknvh298/uu6XTcRm9gEdXuYZWlSvIRZqn0UplSEWtCV2IobY6U9nbYhZKyA+ceyUE+kNsRAfWtQIQyzo8DKVITZMSUTNpZHyi8FZaQyxoC2xRRhiIT6LUipDrJldDcNoSap2YiKyUkQe9v4dFJHrRaRXRO4WkSfKP+dUO5ZhGEZqqnZizrmNzrnVzrnVwMuBI8D3gBuBtc65FcDa8rphGEah5MqJicivAX/snLtQRDYCFznndorIIuBe59zK2Ps7SyW35PosJ5ZHLePngLZ/9FVqW8wQm2cW7Xpzbcfm6msYG6JUbTamZXdkuYdqQ5R8wlzbnldmD8Z51DJh/ueZyzNzaUzhA3FDbGd/YLuNDFHyc21hGcfJ4POMDVGK5dqmDOlz+uUXMYUP6Fzb7pfpzyiPWqbWGYJA59qqzaJ9tC9bjg1RipVxgDbEVhuiFMu1dUzJfpeYwgfihti1P/pYkpzYlcCt5eU+59yp9O0uoG+8N4jItSKyTkTWDR8+PN4uhmEYdVNzJyYiHcAVwLfDbW7kcW7cRzrn3M3OuTXOuTVtPT3j7WIYhlE3NYeTIvIm4L3OuV8rr+cOJ2dKr3ulXDK6nsfK4IeX3c8Ekyo0wBAL8QlIVNsSGWJBW2InYog985Ydo8t5rAxheYFvDElliIV4dX8jDLGgw8tUhtgx5Q45rAx+eDn1oH5fvYZY0JbYIgyxoMPLRhli7//GByccTl5FFkoC3AlcXV6+GrhjzDsMwzAaTE2dmIj0AJcCt3svfwq4VESeAF5bXjcMwyiUmuaddM4dBuYGr+0DLhn/HYZhGMVQ6LCjaUtLbun/zkos8lgZ/BzZsbk6Hm+EIRZqn0UplSEWtCV2IobYtuPZeh4rQ1h+cWB5dp5UhliID1FqhCEWdI4slSG2a5cOZvJYGfwc2VB3vGyiVkMsaEtsEYZYiM+ilMoQ+8N/+YgNOzIMo/WwTswwjKbGOjHDMJqaYlU880tu5duynFgetYyfI9v1Kp23aYQhFmqfRSmVIRbiltg8hthpXh1UHrVMWEO2+arZo8upDLEQH1rUCEMs6BxZKkOs+46eujvP0CI/R7btDbqGrF5DLGhLbBGGWIjPopTKELvpD2y2I8MwWhDrxAzDaGoKDSdnTl/izjvn90bX/bKJmLUCdHg5c7Pet1UMsaAtsRMxxHbuydoUs1ZAfOLYoa70hliIDy1qhCEWdHiZyhBb+vfKZROgw8tY+UW1solaDbGgLbFFGGIhPgFJKkPsrRd81cJJwzBaD+vEDMNoaqwTMwyjqSk0JyYie4BtwDxgb5Xdi8TaE+f51h54/rXJ2hMnRXtOc87ND18stBMbPanIuvESdJOFtSfO86098Pxrk7UnTiPbY+GkYRhNjXVihmE0NZPVid08SeethLUnzvOtPfD8a5O1J07D2jMpOTHDMIxUWDhpGEZTU2gnJiKXichGEdkkIpMyY7iIfE1EdovIY95rvSJyt4g8Uf45J3aMxO0picg9IrJeRH4hItdNZptEZJqI3C8ij5Tb84ny68tF5L7yZ/fN8hR+hSEibSLykIjcNdntEZGtIvJzEXlYRNaVX5u0e6h8/tki8h0ReVxENojIBZN4D60sX5tT/w6KyPWNak9hnZiItAFfBF4PnA1cJSJnF3V+j78HLgteuxFY65xbAawtrxfFEHCDc+5s4HzgveXrMlltOg5c7Jw7B1gNXCYi5wOfBj7rnDsT2A9cU1B7TnEdsMFbn+z2vMY5t9orG5jMewjg88APnHOrgHMYuVaT0ibn3MbytVkNvBw4AnyvYe1xzhXyD7gA+Ddv/aPAR4s6f9CWZcBj3vpGYFF5eRGwcTLaVT7/HYzMLDXpbQK6gZ8Br2SkULF9vM+ygHYsLd/0FwN3ATLJ7dkKzAtem7TPC5gFPEk5x/18aJPXhl8DftzI9hQZTi4BnvbWt5dfez7Q55w75UvYBfRNRiNEZBlwLnDfZLapHLo9DOwG7gY2A/3OuVMKgqI/u88BHwZO6RnmTnJ7HPDvIvKgiFxbfm0y76HlwB7g78oh91fL0yw+H+7rK8nmq21IeyyxH+BG/pso/CtbEZkOfBe43jmn5oMuuk3OuWE3EgosBc4DVhV17hARuRzY7Zx7cLLaMA6/4px7GSOpkfeKyK/6GyfhHmoHXgZ82Tl3LnCYIFSbjPu6nKe8Avh2uC1le4rsxHYAJW99afm15wPPisgigPLP3UWeXESmMtKBfcM5d2qC4kltE4Bzrh+4h5FwbbaInJJ3FfnZXQhcISJbgdsYCSk/P4ntwTm3o/xzNyO5nvOY3M9rO7DdOXdfef07jHRqk30PvR74mXPu2fJ6Q9pTZCf2ALCi/K1SByOPmXcWeP4YdwJXl5evZiQvVQgiIsAtwAbn3N9MdptEZL6IzC4vdzGSn9vASGf29qLb45z7qHNuqXNuGSP3zH865945We0RkR4RmXFqmZGcz2NM4j3knNsFPC0iK8svXQKsn8w2lbmKLJSkYe0pOMn3BuCXjORY/rDoJGO5DbcCO4ETjPwPdg0jOZa1wBPAfwC9BbbnVxh5rH4UeLj87w2T1SbgpcBD5fY8Bny8/PrpwP3AJkbCg85J+OwuAu6azPaUz/tI+d8vTt3Hk3kPlc+/GlhX/tz+GZgzyfd1D7APmOW91pD2WMW+YRhNjSX2DcNoaqwTMwyjqbFOzDCMpsY6McMwmhrrxAzDaGqsEzMMo6mxTswwjKbGOjHDMJqa/w+VowHkEsGUJQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "C_ortho = cayley(C)\n",
    "plt.imshow(C_ortho.detach().numpy().real)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "resistant-ceiling",
   "metadata": {},
   "source": [
    "As an aside, note that the imaginary part of `C_ortho` (and `C`) is zero:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "heard-legislation",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.275218548675184e-06, 1.0707211686167284e-06)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C_ortho.imag.norm().item(), C.imag.norm().item()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dominican-member",
   "metadata": {},
   "source": [
    "Let's illustrate that `C_ortho` is indeed orthogonal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "complicated-entity",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATEAAAFDCAYAAABSjzmiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUEUlEQVR4nO3de6xlZXnH8e+vwwCCFxilkykz6dA4QqmRQadUI22nIBSsgdpYK7ZmbGmmaazRtKZAmzTatA2mF0vShmTqpTS1gqAIpd5wCrU2BhllVGCEQQpl6DCDFgLahAI+/WOvwcP2XPa57Mt7zveTnOy11r6sJ3M2P5733evdJ1WFJLXqh8ZdgCQthiEmqWmGmKSmGWKSmmaISWqaISapaYaYJkaStye5Ocl9SW7vtn933HVpssXrxDRpkrwf+LuqumXctWjy2YlpEp0M3DnuItQGQ0yT6HlV9fi4i1AbDLEVKMnzkvxZknuSPJ7kP5P8TZLjxn3uJBuAB4Zdx6CSHJukkvzouGvR9AyxFSbJMcC/AycB51bV84CfBlYDQ/0PdcBz/wRwxzDrmKfNwCNVdf+4C9H0Dht3ARq59wH/A7yhqr4HUFX7gN+akHP/BJM1H7YZ2D3mGjQLQ2wF6YZqbwF++lCITNq5q+ovh1zHDcDpM9z9hap6Xd+xUzHEJpohtrK8Bni4qr64kCcn2Qi8pKo+2+1vBV5XVe8a9rmXyjQhNZfNwJ8PoRQtEefEVpa1wH8t4vk/Bpw9pnOTZGOSs7vtrUn+YjGvN8D5jgB+HDuxiWaIrSz/BRyfZNrfe5JVSf4xyb8l+Zfuk7mtSf45ybXAVcCvdFfSr+me9tIk1yb5apKXTvcag5x7Lt3zNrLwED30Op9K8p0Zfj7V9/CXAk8DexZzTg2XIbay3NDdXprk+QBJXpLksiQvAl4P7KuqnwWuBN7ePf4FwC8BvwxcVVVbq+p/uvtWV9XrgYuB35jlNWY99wAB+lbgt+lCFFjD4AH6jKo6t6qeO8PPuX0PPxW4vaqeWtC/tkbCEFtBquox4AzgJcDeJI8C1wLfrapvAS8Gbu0efiuwqdveVTOvT9vd3T4AHDvTawxw7lkDtKo+CFxOF6L0PuUcNEAXajMOJSeeE/srTFXdDfziDHffA5wGfAz4SWBvd/zQp4lPAqv6X3LKdmZ5jbnO3R9+ZwOfZ3EBuqihZ1X9zmKer9GwE9NUnwA2JPk8cAHwN333fx14RZKruwtXF/IaMzkUfjB9gMIPhuhMAdr/GlrG/BYLTYQkhwH/AKwHvgP8GvAyplzC0c2l3QAcAD4MnF5V70ryUuBdwG/2v8aUuTstU4aYpKY5nJTUtEWFWJJzktzVfSPBxUtVlCQNasHDySSrgLuBs4B99D4NuqCqJmnxrqRlbjGXWJwG3FNV9wIkuRI4n1m+geBFa1bVxg2rn9m/+2tHLeL0klaSx3nkW1X1A995t5gQO55nf3ndPuCnZnvCxg2r+dJnNjyz//M/snkRp5e0knyurpn2O92GPrGfZHuSXUl2Pfztp4d9OkkrzGJC7EFgw5T99d2xZ6mqHVW1paq2HPfC/ou9JWlxFjOcvBXYlOQEeuH1JuDNsz3h7q8d9awh5Gf+e/cz2w4tJS3EgkOsqp5K8jvAZ+gtBflgVU3Sd6NLWgEWtQC8qj4JfHKJapGkeRvrt1jMNLTsv0+SZuKyI0lNM8QkNc0Qk9S0iflm1/45MOfIJA3CTkxS0wwxSU2bmOFkv9mGlw4tJR1iJyapaYaYpKYZYpKaNrFzYv1coiRpOnZikppmiElqmiEmqWnNzIlN5TVkkg6xE5PUNENMUtOaHE728/ILaeWyE5PUNENMUtMMMUlNWxZzYlP5DbHSymInJqlphpikphlikpq27ObE+rlESVre5uzEknwwycEkt085tibJjUn2drfHDrdMSZreIMPJvwfO6Tt2MbCzqjYBO7t9SRq5OYeTVfX5JBv7Dp8PbO22rwBuBi5aysKGxSVK0vKy0In9tVW1v9t+CFi7RPVI0rws+tPJqiqgZro/yfYku5LsepInFns6SXqWhYbYgSTrALrbgzM9sKp2VNWWqtqymiMWeDpJmt5CL7G4HtgGXNrdXrdkFY2QS5Sk9g1yicVHgC8CJybZl+RCeuF1VpK9wGu6fUkauUE+nbxghrvOXOJaJGneXHYkqWnLftnRfLhESWqPnZikphlikprmcHIWLlGSJp+dmKSmGWKSmmaISWqac2ID8vILaTLZiUlqmiEmqWkOJxfIyy+kyWAnJqlphpikphlikprmnNgS8BtipfGxE5PUNENMUtMMMUlNc05sCFyiJI2OnZikphlikprmcHIEXKIkDY+dmKSmGWKSmmaISWqac2Ij5hIlaWnZiUlq2pwhlmRDkpuS3JnkjiTv6I6vSXJjkr3d7bHDL1eSnm2QTuwp4Peq6mTglcDbkpwMXAzsrKpNwM5uX5JGas45saraD+zvth9Psgc4Hjgf2No97ArgZuCioVS5jLlESVqceU3sJ9kInArcAqztAg7gIWDtDM/ZDmwHOJKjFlyoJE1n4In9JM8FPga8s6oem3pfVRVQ0z2vqnZU1Zaq2rKaIxZVrCT1G6gTS7KaXoB9uKo+3h0+kGRdVe1Psg44OKwiVxKXKEnzM8inkwE+AOypqr+actf1wLZuextw3dKXJ0mzG6QTezXwFuDrSXZ3x/4AuBT4aJILgfuBNw6lQkmaxSCfTn4ByAx3n7m05UjS/LjsaIJ5+YU0N5cdSWqaISapaQ4nG+LlF9IPshOT1DRDTFLTDDFJTXNOrFF+Q6zUYycmqWmGmKSmGWKSmuac2DLhEiWtVHZikppmiElqmsPJZcolSlop7MQkNc0Qk9Q0Q0xS05wTWwFcoqTlzE5MUtMMMUlNM8QkNc05sRXIJUpaTuzEJDXNEJPUNIeTcomSmmYnJqlpc4ZYkiOTfCnJV5PckeQ93fETktyS5J4kVyU5fPjlStKzDdKJPQGcUVWnAJuBc5K8Engv8L6qejHwCHDh0KqUpBnMOSdWVQV8p9td3f0UcAbw5u74FcC7gcuXvkSNkpdfqDUDzYklWZVkN3AQuBH4JvBoVT3VPWQfcPxQKpSkWQwUYlX1dFVtBtYDpwEnDXqCJNuT7Eqy60meWFiVkjSDeV1iUVWPJrkJeBVwTJLDum5sPfDgDM/ZAewAeH7W1CLr1Yh5+YUm3SCfTh6X5Jhu+znAWcAe4CbgDd3DtgHXDalGSZrRIJ3YOuCKJKvohd5Hq+qGJHcCVyb5E+A24ANDrFOSpjXIp5NfA06d5vi99ObHJGlsXHakgfkNsZpELjuS1DRDTFLTDDFJTXNOTAvmEiVNAjsxSU0zxCQ1zeGkloxLlDQOdmKSmmaISWqaISapac6JaShcoqRRsROT1DRDTFLTDDFJTXNOTCPhEiUNi52YpKYZYpKa5nBSY+ESJS0VOzFJTTPEJDXNEJPUNOfENHZefqHFsBOT1DRDTFLTHE5q4nj5hebDTkxS0wYOsSSrktyW5IZu/4QktyS5J8lVSQ4fXpmSNL35dGLvAPZM2X8v8L6qejHwCHDhUhYmSYMYaE4syXrgF4A/BX43SYAzgDd3D7kCeDdw+RBq1ArmN8RqLoN2Yn8N/D7wvW7/hcCjVfVUt78POH5pS5Okuc0ZYkleBxysqi8v5ARJtifZlWTXkzyxkJeQpBkNMpx8NXBektcCRwLPBy4DjklyWNeNrQcenO7JVbUD2AHw/KypJalakjpzhlhVXQJcApBkK/CuqvrVJFcDbwCuBLYB1w2vTKnHJUrqt5jrxC6iN8l/D705sg8sTUmSNLh5XbFfVTcDN3fb9wKnLX1JkjQ4lx2paS5RksuOJDXNEJPUNENMUtOcE9Oy4RKllclOTFLTDDFJTTPEJDXNOTEtWy5RWhnsxCQ1zRCT1DSHk1oxXKK0PNmJSWqaISapaYaYpKY5J6YVycsvlg87MUlNM8QkNc3hpISXX7TMTkxS0wwxSU0zxCQ1zTkxqY/fENsWOzFJTTPEJDXNEJPUNOfEpDm4RGmyDRRiSe4DHgeeBp6qqi1J1gBXARuB+4A3VtUjwylTkqY3n+Hkz1XV5qra0u1fDOysqk3Azm5fkkZqMcPJ84Gt3fYVwM3ARYusR5p4LlGaLIN2YgV8NsmXk2zvjq2tqv3d9kPA2iWvTpLmMGgndnpVPZjkh4Ebk3xj6p1VVUlquid2obcd4EiOWlSxktRvoE6sqh7sbg8C1wKnAQeSrAPobg/O8NwdVbWlqras5oilqVqSOnN2YkmOBn6oqh7vts8G/hi4HtgGXNrdXjfMQqVJ5BKl8RtkOLkWuDbJocf/U1V9OsmtwEeTXAjcD7xxeGVK0vTmDLGquhc4ZZrj3wbOHEZRkjQolx1JaprLjqQl5BKl0bMTk9Q0Q0xS0xxOSkPkEqXhsxOT1DRDTFLTDDFJTXNOTBoRL78YDjsxSU0zxCQ1zeGkNCZefrE07MQkNc0Qk9Q0Q0xS05wTkyaA3xC7cHZikppmiElqmiEmqWnOiUkTyCVKg7MTk9Q0Q0xS0xxOSg1widLM7MQkNc0Qk9Q0Q0xS05wTkxrjEqVnsxOT1LSBQizJMUmuSfKNJHuSvCrJmiQ3Jtnb3R477GIlqd+gndhlwKer6iTgFGAPcDGws6o2ATu7fUkaqTnnxJK8APgZ4K0AVfV/wP8lOR/Y2j3sCuBm4KJhFClpZit9idIgndgJwMPAh5LcluT9SY4G1lbV/u4xDwFrp3tyku1JdiXZ9SRPLE3VktQZJMQOA14OXF5VpwLfpW/oWFUF1HRPrqodVbWlqras5ojF1itJzzLIJRb7gH1VdUu3fw29EDuQZF1V7U+yDjg4rCIlDW6lLVGasxOrqoeAB5Kc2B06E7gTuB7Y1h3bBlw3lAolaRaDXuz6duDDSQ4H7gV+nV4AfjTJhcD9wBuHU6IkzWygEKuq3cCWae46c0mrkaR5ctmRtIythMsvXHYkqWmGmKSmOZyUVpDlePmFnZikphlikppmiElqmnNi0gq1XL4h1k5MUtMMMUlNM8QkNc05MUlAu0uU7MQkNc0Qk9Q0h5OSptXKEiU7MUlNM8QkNc0Qk9Q058QkzWmSlyjZiUlqmiEmqWmGmKSmOScmad4maYmSnZikphlikprmcFLSoo1ziZKdmKSmzRliSU5MsnvKz2NJ3plkTZIbk+ztbo8dRcGSNNWcIVZVd1XV5qraDLwC+F/gWuBiYGdVbQJ2dvuSNFLznRM7E/hmVd2f5Hxga3f8CuBm4KKlK01Si0Z9+cV8Q+xNwEe67bVVtb/bfghYO90TkmwHtgMcyVELqVGSZjTwxH6Sw4HzgKv776uqAmq651XVjqraUlVbVnPEgguVpOnMpxM7F/hKVR3o9g8kWVdV+5OsAw4ufXmSWjfsyy/mc4nFBXx/KAlwPbCt294GXLfoaiRpngYKsSRHA2cBH59y+FLgrCR7gdd0+5I0UgMNJ6vqu8AL+459m96nlZI0Ni47kjQyw/iGWJcdSWqaISapaYaYpKY5JyZpbJZiiZKdmKSmGWKSmuZwUtLEmG2J0qp10z/HTkxS0wwxSU0zxCQ1Lb2vAhvRyZKHgfuBFwHfGtmJ52Y9s5u0emDyarKe2S1FPT9aVcf1HxxpiD1z0mRXVW0Z+YlnYD2zm7R6YPJqsp7ZDbMeh5OSmmaISWrauEJsx5jOOxPrmd2k1QOTV5P1zG5o9YxlTkySlorDSUlNG2mIJTknyV1J7kkylr8YnuSDSQ4muX3KsTVJbkyyt7s9doT1bEhyU5I7k9yR5B3jrCnJkUm+lOSrXT3v6Y6fkOSW7nd3Vfcn/EYmyaoktyW5Ydz1JLkvydeT7E6yqzs2tvdQd/5jklyT5BtJ9iR51RjfQyd2/zaHfh5L8s5h1TOyEEuyCvhben/67WTggiQnj+r8U/w9cE7fsYuBnVW1CdjZ7Y/KU8DvVdXJwCuBt3X/LuOq6QngjKo6BdgMnJPklcB7gfdV1YuBR4ALR1TPIe8A9kzZH3c9P1dVm6dcNjDO9xDAZcCnq+ok4BR6/1Zjqamq7ur+bTYDrwD+F7h2aPVU1Uh+gFcBn5myfwlwyajO31fLRuD2Kft3Aeu67XXAXeOoqzv/dfT+stTYawKOAr4C/BS9CxUPm+53OYI61ndv+jOAG4CMuZ77gBf1HRvb7wt4AfCfdHPck1DTlBrOBv5jmPWMcjh5PPDAlP193bFJsLaq9nfbDwFrx1FEko3AqcAt46ypG7rtpvcHkW8Evgk8WlVPdQ8Z9e/ur4HfB77X7b9wzPUU8NkkX06yvTs2zvfQCcDDwIe6Iff7uz+zOAnv6zfx/b9XO5R6nNjvU73/TYz8I9skzwU+Bryzqh4bZ01V9XT1hgLrgdOAk0Z17n5JXgccrKovj6uGaZxeVS+nNzXytiQ/M/XOMbyHDgNeDlxeVacC36VvqDaO93U3T3kecHX/fUtZzyhD7EFgw5T99d2xSXAgyTqA7vbgKE+eZDW9APtwVR36A8VjrQmgqh4FbqI3XDsmyaHvnxvl7+7VwHlJ7gOupDekvGyM9VBVD3a3B+nN9ZzGeH9f+4B9VXVLt38NvVAb93voXOArVXWg2x9KPaMMsVuBTd2nSofTazOvH+H5Z3M9sK3b3kZvXmokkgT4ALCnqv5q3DUlOS7JMd32c+jNz+2hF2ZvGHU9VXVJVa2vqo303jP/WlW/Oq56khyd5HmHtunN+dzOGN9DVfUQ8ECSE7tDZwJ3jrOmzgV8fyjJ0OoZ8STfa4G76c2x/OGoJxm7Gj4C7AeepPd/sAvpzbHsBPYCnwPWjLCe0+m11V8Ddnc/rx1XTcDLgNu6em4H/qg7/mPAl4B76A0PjhjD724rcMM46+nO+9Xu545D7+Nxvoe6828GdnW/t08Ax475fX008G3gBVOODaUer9iX1DQn9iU1zRCT1DRDTFLTDDFJTTPEJDXNEJPUNENMUtMMMUlN+39fmirXKsoBIgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow((C_ortho @ C_ortho.T).detach().numpy().real)\n",
    "plt.title(r\"$C_\\mathsf{ortho} C_\\mathsf{ortho}^T = I$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "worst-entity",
   "metadata": {},
   "source": [
    "--------------------------------------------------------\n",
    "## Orthogonalizing FFT Convolutions\n",
    "\n",
    "Putting the ideas illustrated above together, we want to orthogonalize\n",
    "a convolution in the Fourier domain without explicitly constructing a large matrix.\n",
    "Essentially, we just want to operate on the blocks of `D`.\n",
    "Recall that we extracted those blocks from the convolution weight tensor earlier:\n",
    "*we will just rearrange the tensor to make working with those blocks easier.*\n",
    "\n",
    "In particular, we will use the fact that each pixel of the output is given by a matrix vector product:\n",
    "we will arrange the input and weight tensors so that we can compute the convolution\n",
    "with a batch of such matrix multiplications."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "minor-jason",
   "metadata": {},
   "outputs": [],
   "source": [
    "xfft = torch.fft.fft2(x)\n",
    "xfft = xfft.permute(2, 3, 1, 0) # batches, cin, n, n -> n, n, cin, batches\n",
    "xfft = xfft.reshape(n**2, cin, batches) # n**2 input pixels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "infinite-liability",
   "metadata": {},
   "outputs": [],
   "source": [
    "wfft = shift_matrix * torch.fft.fft2(wpad).conj() # cout, cin, n, n\n",
    "wfft = wfft.reshape(cout, cin, n**2) # cout, cin, n**2 pixels\n",
    "wfft = wfft.permute(2, 0, 1) # n**2, cout, cin"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "specific-breakdown",
   "metadata": {},
   "source": [
    "Now, similarly to the block diagonal matrix `D` earlier,\n",
    "we have `xfft` of shape ($n^2$, `cin`, `batches`) and\n",
    "`wfft` of shape ($n^2$, `cout`, `cin`).\n",
    "The last two dimensions for each of these are compatible for matrix multiplication with output\n",
    "shape ($n^2$, `cout`, `batches`).\n",
    "\n",
    "We can do batch matrix multiplication, and reshape again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "insured-botswana",
   "metadata": {},
   "outputs": [],
   "source": [
    "yfft = wfft @ xfft # n**2, cout, batches\n",
    "yfft = yfft.reshape(n, n, cout, batches)\n",
    "yfft = yfft.permute(3, 2, 0, 1) # batches, cout, n, n\n",
    "y4 = torch.fft.ifft2(yfft)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "therapeutic-vietnamese",
   "metadata": {},
   "source": [
    "Again, this is equivalent to the original circular convolution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "short-malaysia",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.5329876532632625e-06"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(y1 - y4).norm().item()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "quality-retreat",
   "metadata": {
    "tags": []
   },
   "source": [
    "But this time, it's very easy to orthogonalize the convolution by orthogonalizing\n",
    "all of the \"blocks\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "invisible-quantity",
   "metadata": {},
   "outputs": [],
   "source": [
    "# \"Batched\" Cayley transform\n",
    "wfft_skew_sym = wfft - wfft.conj().transpose(1, 2)\n",
    "I = torch.eye(cin, dtype=wfft.dtype)[None, :, :]\n",
    "wfft_ortho = (I - wfft_skew_sym) @ torch.inverse(I + wfft_skew_sym)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "heated-ranking",
   "metadata": {},
   "outputs": [],
   "source": [
    "yfft_ortho = wfft_ortho @ xfft # n**2, cout, batches\n",
    "yfft_ortho = yfft_ortho.reshape(n, n, cout, batches)\n",
    "yfft_ortho = yfft_ortho.permute(3, 2, 0, 1) # batches, cout, n, n\n",
    "y_ortho = torch.fft.ifft2(yfft_ortho)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "lyric-pocket",
   "metadata": {},
   "source": [
    "We can see that this implements an orthogonal convolution because\n",
    "it doesn't expand the norm of input tensors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "completed-welsh",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y_ortho   norm:  8.151487350463867\n",
      "input (x) norm:  8.151487350463867\n"
     ]
    }
   ],
   "source": [
    "print('y_ortho   norm: ', y_ortho.norm().item())\n",
    "print('input (x) norm: ', x.norm().item())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "understood-casting",
   "metadata": {},
   "source": [
    "---------------------------------------------------------------------\n",
    "### Optimizations and More"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "prime-efficiency",
   "metadata": {},
   "source": [
    "Using the symmetry of the Fourier transform, we can implement FFT-based convolutions\n",
    "even more efficiently. We can also handle strided FFT-based convolutions using the *aliasing theorem*, which can be found [here](https://ccrma.stanford.edu/~jos/st/Downsampling_Theorem_Aliasing_Theorem.html).\n",
    "\n",
    "We provide a more fully-featured FFT-based convolution layer in `extras/fftconv.py`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "outside-declaration",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0860642305488e-06"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from extras.fftconv import FFTConv\n",
    "\n",
    "fft_conv = FFTConv(cin, cout, k, bias=False)\n",
    "fft_conv.weight.data = conv.weight.data\n",
    "\n",
    "(fft_conv(x) - y1).norm().item()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "conceptual-calcium",
   "metadata": {},
   "source": [
    "### Orthogonal Convolutions when `cin` != `cout`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "norwegian-dictionary",
   "metadata": {},
   "source": [
    "In the appendix of our paper, we describe a method to efficiently \"orthogonalize\" convolutions with different numbers of input and output channels. This is actually called *semi-orthogonalization*.\n",
    "More details can be found in `Orthogonal Convolutions.ipynb`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adaptive-bristol",
   "metadata": {},
   "source": [
    "-----------------------------------------------------------------------\n",
    "You can cite our work as follows:\n",
    "\n",
    "<pre>\n",
    "@inproceedings{trockman2021ortho,\n",
    "    title={Orthogonalizing Convolutional Layers with the Cayley Transform},\n",
    "    author={Asher Trockman and J. Zico Kolter},\n",
    "    booktitle={International Conference on Learning Representations},\n",
    "    year={2021},\n",
    "    url={https://openreview.net/forum?id=Pbj8H_jEHYv}\n",
    "}\n",
    "</pre>"
   ]
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}