Tutorial_Group_Convolution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "N2dSAwLACjSN"
   },
   "source": [
    "# Tutorial: (Relaxed) Group Convolution\n",
    "\n",
    "\n",
    "# [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/iaifi/summer-school-2023/blob/main/Tutorial_Symmetry.ipynb)\n",
    "\n",
    "This google colab tutorial is supporting material for the lectures by Robin Walters on *Incorporating Symmetry into Deep Dynamics Models* at [IAIFI Summer School](https://iaifi.org/phd-summer-school.html). We will 1) build regular group  convolution neural net layers; 2) define approximate equivariance and implement relaxed group convolution for translation group and C4 rotation group; 3) introduce basic ideas of steerable convolutions and [e2cnn](https://github.com/QUVA-Lab/e2cnn). 4) introduce several methods that can improve models' long-term dynmaics prediction accuracy.\n",
    "\n",
    "This tutorial is based on the following materials:\n",
    "\n",
    "* [Dr. Denis Boyda's tutorial](https://www.dropbox.com/sh/450ao3m5v4qt6o3/AABe1B2jEcyM5ZyoY_pmOStfa?dl=0) from IAIFI Summer School 2022.\n",
    "\n",
    "* [UvA Deep Learning Tutorials](https://uvadlc-notebooks.readthedocs.io/en/latest/index.html).\n",
    "\n",
    "* Prof. Erik Bekkers's course [Group Equivariant Deep Learning](https://www.youtube.com/@erikbekkers6398/videos).\n",
    "\n",
    "* Cohen et al. [Group Equivariant Convolutional Networks](https://arxiv.org/pdf/1602.07576.pdf).\n",
    "\n",
    "* Cohen et al. [Steerable CNNs](https://arxiv.org/abs/1612.08498).\n",
    "\n",
    "* [E2CNN](https://github.com/QUVA-Lab/e2cnn) & [ESCNN](https://github.com/QUVA-Lab/escnn) & [E3NN](https://e3nn.org/)\n",
    "\n",
    "* Wang et al. [Discovering Symmetry Breaking in Physical Systems with Relaxed Group Convolution](https://arxiv.org/abs/2310.02299).\n",
    "\n",
    "* Wang and Walters et al. [Approximately Equivariant Networks for Imperfectly Symmetric Dynamics](https://proceedings.mlr.press/v162/wang22aa.html).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7YNFbmARCvTf"
   },
   "source": [
    "# 1. Group Equivariant Convolutional Networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "p2xNgzMjGpeP"
   },
   "source": [
    "## 1.1 Group, Equivariance and Group Equivariant Convolutional Networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "DGxx9GLXY6_M"
   },
   "source": [
    "### 1.1.1 Definition of group\n",
    "\n",
    "A Group (or Abstract group) $G = (\\Omega, \\cdot)$ is a set with a binary operation, called group multiplication, that satisfies the following four axioms:\n",
    "\n",
    "1.   <ins>Closure</ins>: For any $g_1, g_2 \\in G$ the product $g_1 \\cdot g_2 = g_3 \\in G$.\n",
    "2.   <ins>Associativity</ins>: For any $g_1, g_2, g_3 \\in G$: ($g_1 \\cdot g_2) \\cdot g_3 = g_1 \\cdot (g_2 \\cdot g_3)$.\n",
    "3.   <ins>Identity</ins>: There exists an identity element $e$ s.t. $e \\cdot g = g \\cdot e = g$ for each $g \\in G$.\n",
    "4.   <ins>Inverse</ins>: For any $g \\in G$ there exists an inverse element $g^{-1} \\in G$, s.t. $g \\cdot g^{-1} = e$.\n",
    "\n",
    "E.g. Translation Group $(R^2, +)$ consists of all possible translations in $R^2$ and is equipped with group product and group inverse,\n",
    "\n",
    "$$g \\cdot g' = (x + x') ; \\;\\;\\ g^{-1} = -x $$\n",
    "\n",
    "with $g = (x), g' = (x')$ and $x, x' \\in \\mathbb{R}^2$.\n",
    "\n",
    "E.g. Roto-translation group $SE(2) = \\mathbb{R}^2 \\rtimes SO(2)$ consists of the coupled space of translation vectors and rotations and is equipped with group product and group inverse,\n",
    "$$g \\cdot g' = (x, R_{\\theta}) \\cdot (x', R_{\\theta'}) = (R_{\\theta}x' + x, R_{\\theta+\\theta'}) ; \\;\\;\\ g^{-1} = (-R_{\\theta}^{-1}x, -R_{\\theta}^{-1}) $$\n",
    "\n",
    "with $g = (x, R_{\\theta}), g' = (x', R_{\\theta'})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EdQvrsc4aPzR"
   },
   "source": [
    "### 1.1.2 Equivariance and Invariance\n",
    "\n",
    "Formally, a function $f \\colon X \\to Y$ may be described as respecting the symmetry coming from a group $G$ using the notion of equivariance.  Assume a group representation $\\rho_{\\text{in}}$ of $G$ acts on $X$ and $\\rho_{\\text{out}}$ acts on $Y$. We say a function $f$ is <ins>$G$-equivariant</ins> if\n",
    "\n",
    "$$ f( \\rho_{\\text{in}}(g)(x)) = \\rho_{\\text{out}}(g) f(x) $$\n",
    "for all $x \\in X$ and $g \\in G$. The function $f$ is <ins>$G$-invariant</ins> if $f( \\rho_{\\text{out}}(g)(x)) = f(x)$ for all $x \\in X$ and $g \\in G$.  This is a special case of equivariance for the case $\\rho_{\\mathrm{out}}(g) = 1$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5d7OEDyP1e6s"
   },
   "source": [
    "### 1.1.3 Recap on CNN\n",
    "\n",
    "<ins>CNN is nothing else but group convolution with the translation group. </ins>\n",
    "\n",
    "At each layer of a regular convnet convolves or correlates stack of feature maps $f : \\mathbb{Z}^2 \\rightarrow \\mathbb{R}^K$ with a filter (or a set of filters) $\\psi : \\mathbb{Z}^2 \\rightarrow \\mathbb{R}^K$.\n",
    "\n",
    "\n",
    "$$[f * \\psi] (\\mathbf{x}) = \\sum_{\\mathbf{y} \\in \\mathbb{Z}^2}f(\\mathbf{y})\\psi(\\mathbf{y}-\\mathbf{x}) = \\sum_{\\mathbf{y} \\in \\mathbb{Z}^2}f(\\mathbf{y})T_{\\mathbf{x}}\\psi(\\mathbf{y}) $$\n",
    "\n",
    "A left regular representation is a representation that transforms the function $f$ by transforming their domains via inverse group actions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "cFQhHnEo--_2"
   },
   "source": [
    "### 1.1.4 Group Equivariant Convolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9HVYP4R7_pJ7"
   },
   "source": [
    "<ins>G-CNN is a generalization of CNN for the case of an arbitrary group $G$.</ins>  \n",
    "\n",
    "By replacing the shift with a more general transformation from some group $G$, We can make models become equivariant to other groups.\n",
    "\n",
    "* The first layer of a G-CNN (Lifting Convolution): $$[f * \\psi] (g) =  \\sum_{\\mathbf{y} \\in \\mathbb{Z}^2}f(\\mathbf{y})\\psi(g^{-1}\\mathbf{y})$$\n",
    "\n",
    "* For all layers after the first, both the input and the filters are functions on $G$:\n",
    "$$[f * \\psi] (g) =  \\sum_{h \\in G}f(h)\\psi(g^{-1}h)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1uNB5aq6EVzH"
   },
   "source": [
    "<ins>E.g.</ins> $G=\\mathbb{Z}^2 \\rtimes C_4$\n",
    "\n",
    "* <ins>Lifting Convolution</ins> maps signals $f_{in}: \\mathbb{Z}^2 → R$ to signals on arbitrary group G: $f_{out}: G→R$.  Group elements $g \\in G = \\mathbb{Z}^2 \\rtimes C_4$ can be parametrized using $\\mathbf{x} \\in R^2$ and $r \\in C_4$ such that $g=(\\mathbf{x}, r)$. For such groups lifting convolution has form\n",
    "$$\n",
    "\\begin{aligned}\n",
    "f_{out}(\\mathbf{x}, r) = (f_{in} * \\psi) (\\mathbf{x})\n",
    "    &=\n",
    "\\sum_{\\mathbf{y} \\in \\mathbb{Z}^2} f_{in}(\\mathbf{y}) L_{r}  L_{\\mathbf{x}} \\psi(\\mathbf{y}),\\\\\n",
    "    &=\n",
    "    \\sum_{\\mathbf{y} \\in \\mathbb{Z}^2}  f_{in}(\\mathbf{y}) L_{r} \\psi(\\mathbf{y} -\\mathbf{x}),\\\\\n",
    "    &=\n",
    "    \\sum_{\\mathbf{y} \\in \\mathbb{Z}^2}  f_{in}(\\mathbf{y})\\psi_{r}(\\mathbf{y} -\\mathbf{x}),\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "<ins>So lifting convolution is the same as conventional convolutions with a stack of filter bank </ins>$\\{k_{r}(\\mathbf{x})\\}_{r \\in C_4}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "aEAK-TdRWn65"
   },
   "source": [
    "* <ins>Group Convolution:</ins> Note that the input feature map, $f_{in}$, now has one additional group dimension defined over $C_4$, besides the usual spatial dimensions defined over $\\mathbb{R}^2$,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "f_{out}(\\mathbf{x}, r) = (f_{in} * \\psi) (\\mathbf{x}, r)\n",
    "    &=\n",
    "\\sum_{\\mathbf{y} \\in \\mathbb{Z}^2} \\sum_{r' \\in C^4} f_{in}(\\mathbf{y} , r) L_{r} L_{\\mathbf{x}} \\psi(\\mathbf{y},r'),\\\\\n",
    "    &=\n",
    "    \\sum_{\\mathbf{y} \\in \\mathbb{Z}^2} \\sum_{r' \\in C^4} f_{in}(\\mathbf{y} , r')  L_{r} \\psi(\\mathbf{y-x},r'),\\\\\n",
    "    &=\n",
    "    \\sum_{\\mathbf{y} \\in \\mathbb{Z}^2} \\sum_{r' \\in C^4}  f_{in}(\\mathbf{y}, r') \\psi_r(\\mathbf{y-x},r')\\\\\n",
    "    &=\n",
    "    \\sum_{\\mathbf{y} \\in \\mathbb{Z}^2} \\sum_{r' \\in C^4}  f_{in}(\\mathbf{y}, r') \\psi(r^{-1}(\\mathbf{y-x}),r^{-1}r')\\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "id": "WFUJF8OYeJS1",
    "tags": []
   },
   "outputs": [],
   "source": [
    "# %%capture\n",
    "# !rm symm4ml_obfuscated.zip\n",
    "# !wget -O symm4ml_obfuscated.zip https://symm4ml.mit.edu/_static/symm4ml_s24/symm4ml/symm4ml_obfuscated.zip\n",
    "# !unzip -o symm4ml_obfuscated.zip\n",
    "# import sys\n",
    "# Update to proper folder if doing this locally\n",
    "# sys.path.append('/pscratch/sd/r/rwang2/SymmetryCourse/symm4ml_obfuscated/symm4ml')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "ux6sZ6fIeKmm",
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from PIL import Image\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import torch\n",
    "import torchvision\n",
    "from torchvision.datasets import MNIST\n",
    "from torchvision import transforms\n",
    "\n",
    "from torch import nn\n",
    "from torch import optim\n",
    "import torch.nn.functional as F\n",
    "\n",
    "# import symm4ml\n",
    "# import symm4ml.groups as groups\n",
    "# import symm4ml.vis as vis\n",
    "# import symm4ml.rep as rep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import group_conv\n",
    "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "V31DXglsOlXQ"
   },
   "source": [
    "## 1.2 C4 Cyclic Group"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oRpZpbuvImdI"
   },
   "source": [
    "**Let's first implement a function that can rotate images**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "fMrerIAaMr5N",
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from PIL import Image\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import math\n",
    "import random\n",
    "import torch\n",
    "import torchvision\n",
    "from torch import nn\n",
    "from torch import optim\n",
    "import torch.nn.functional as F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "id": "gogWJEQMNHQI",
    "tags": []
   },
   "outputs": [],
   "source": [
    "def rot_img(x, theta):\n",
    "    \"\"\" Rotate 2D images\n",
    "    Args:\n",
    "        x : input images with shape [N, C, H, W]\n",
    "        theta: angle\n",
    "    Returns:\n",
    "        rotated images\n",
    "    \"\"\"\n",
    "    # Rotation Matrix (2 x 3)\n",
    "    rot_mat = torch.FloatTensor([[np.cos(theta), -np.sin(theta), 0],\n",
    "                                 [np.sin(theta), np.cos(theta), 0]])\n",
    "\n",
    "    # The affine transformation matrices should have the shape of N x 2 x 3\n",
    "    rot_mat = rot_mat.repeat(x.shape[0],1,1)\n",
    "\n",
    "    # Obtain transformed grid\n",
    "    # grid is the coordinates of pixels for rotated image\n",
    "    # F.affine_grid assumes the origin is in the middle\n",
    "    # and it rotates the positions of the coordinates\n",
    "    # r(f(x)) = f(r^-1 x)\n",
    "    grid = F.affine_grid(rot_mat.to(x.device), x.size(), align_corners=False).float()\n",
    "    x = F.grid_sample(x, grid)\n",
    "    return x.float()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 362
    },
    "id": "aOHncymzNPo9",
    "outputId": "3b792031-825d-4461-a19f-7ebe649b1249",
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/global/homes/r/rwang2/.conda/envs/e3nn/lib/python3.9/site-packages/torch/nn/functional.py:4236: UserWarning: Default grid_sample and affine_grid behavior has changed to align_corners=False since 1.3.0. Please specify align_corners=True if the old behavior is desired. See the documentation of grid_sample for details.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABlgAAAGCCAYAAACM8s6KAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAv80lEQVR4nO3dfZBV9X0/8M/ysLtY2SVEYUEQ8QmMD4CouJgRHKlEqRM6TWpMGoxR0iTQarBJpdNGY2a6vzQaM01JjM0obQyjiUZtjTFFEBwjauRhqsbQQCigPxaSqLtAdHnY8/sjP9es7ALf3XP3PvB6zZyZ3Lvn3P18c8/1vfe+ufdWZVmWBQAAAAAAAIetX7EHAAAAAAAAKDcKFgAAAAAAgEQKFgAAAAAAgEQKFgAAAAAAgEQKFgAAAAAAgEQKFgAAAAAAgEQKFgAAAAAAgEQKFgAAAAAAgEQKFgAAAAAAgEQKFgAAAAAAgEQFK1hee+21+NjHPhZ1dXUxZMiQuOaaa2LXrl0HPWb69OlRVVXVafv0pz9dqBEBKBMyBYA8yBMA8iJTAIiIqMqyLCvEDV966aWxbdu2+Pa3vx179+6Nq6++Os4999xYsmRJt8dMnz49Tj311Ljllls6rjvqqKOirq6uECMCUCZkCgB5kCcA5EWmABARMaAQN/ryyy/HY489Fj/72c/inHPOiYiIb3zjG3HZZZfFrbfeGiNHjuz22KOOOioaGhoKMRYAZUimAJAHeQJAXmQKAG8rSMGyatWqGDJkSEfIRETMmDEj+vXrF88++2z86Z/+abfHfu9734t77rknGhoa4vLLL49/+Id/iKOOOqrb/dva2qKtra3jcnt7e7z22mvx3ve+N6qqqvJZEMARIsuy2LlzZ4wcOTL69SuNr+mSKQDlqdQyRZ4AlKdSy5OIvssUeQKQn0LlSUEKlubm5hg2bFjnXzRgQAwdOjSam5u7Pe6jH/1ojBkzJkaOHBn//d//HX/7t38b69evjx/+8IfdHtPU1BRf+tKXcpsdgIitW7fGqFGjij1GRMgUgHJXKpkiTwDKW6nkSUTfZYo8Achf3nmSVLDceOON8ZWvfOWg+7z88ss9HuZTn/pUx/8+88wzY8SIEXHxxRfHxo0b46STTurymIULF8aCBQs6Lre0tMTxxx8fI2+7MfoNqu3xLBTPyZ9dW+wR6KXvr19X7BHooZ272uOkyVtj8ODBBf9d5ZQpG1ePjsFHl8a/liPNy3v2FXsEeume1xqLPQI9tGf33vi3yx4ueKaUU55c9egHo/qPBvZ4Fornv1/r/uN+KA9Hf6F/sUegh/btb4uVv/rmEfkcpbs8mXbiZ2NA/5oez0Hx3P+jh4s9Ar105uNzij0CPdT+5lvxf2/4P7nnSVLBcsMNN8QnPvGJg+5z4oknRkNDQ+zYsaPT9fv27YvXXnst6XMmp0yZEhERGzZs6PbJS01NTdTUHBgq/QbVKljK1IAqTzrLXd1gL0SXu754u3k5Zcrgo/s5r8vU0Xvcb+Wueo+/C8pdoTOlnPKk+o8GRvXRzulyNKDNC5nlbkB/BUu5OxKfo3SXJwP61yhYypTnleXP683lL+88SSpYjj322Dj22GMPuV9jY2O88cYbsXr16pg8eXJERCxfvjza29s7wuNwrFu3LiIiRowYkTImAGVApgCQB3kCQF5kCgCpClKbnnbaafGBD3wg5s6dG88991z89Kc/jfnz58dHPvKRGDny92+tfvXVV2P8+PHx3HPPRUTExo0b48tf/nKsXr06/vd//zf+4z/+I+bMmRMXXnhhnHXWWYUYE4AyIFMAyIM8ASAvMgWAtxXsfWnf+973Yvz48XHxxRfHZZddFu9///vjzjvv7Pj53r17Y/369fG73/0uIiKqq6vj8ccfj0suuSTGjx8fN9xwQ/zZn/1Z/Od//mehRgSgTMgUAPIgTwDIi0wBICLxI8JSDB06NJYsWdLtz0844YTIsqzj8ujRo2PlypWFGgeAMiZTAMiDPAEgLzIFgIgCvoMFAAAAAACgUilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEilYAAAAAAAAEvVJwbJo0aI44YQTora2NqZMmRLPPffcQff/wQ9+EOPHj4/a2to488wz49FHH+2LMQEocfIEgDzIEwDyIlMAjmwFL1juu+++WLBgQdx0002xZs2amDBhQsycOTN27NjR5f5PP/10XHnllXHNNdfE2rVrY/bs2TF79ux48cUXCz0qACVMngCQB3kCQF5kCgBVWZZlhfwFU6ZMiXPPPTf+5V/+JSIi2tvbY/To0fFXf/VXceONNx6w/xVXXBG7d++ORx55pOO6888/PyZOnBh33HHHIX9fa2tr1NfXx6hv3hz9BtXmtxD6zKmffL7YI9BLj7y6utgj0EOtO9tj2LjN0dLSEnV1dcUep5O+zpOIdzJlx/oxUTfYp2qWo5f27Cv2CPTSXb99f7FHoIf27Nob/zrt/pLLlGLmydyVH4rqowfmsxD61Lrfjir2CPTS0X/dv9gj0EP79rfFsl/eXnJ5ElG817wuPuVzMaB/TX4Loc88+sT9xR6BXhr72LXFHoEean/zrXjlszfnnicFfbVoz549sXr16pgxY8Y7v7Bfv5gxY0asWrWqy2NWrVrVaf+IiJkzZ3a7f1tbW7S2tnbaAKgsfZEnETIFoNLJEwDy4jUvACIKXLD85je/if3798fw4cM7XT98+PBobm7u8pjm5uak/ZuamqK+vr5jGz16dD7DA1Ay+iJPImQKQKWTJwDkxWteAET00ZfcF9LChQujpaWlY9u6dWuxRwKgTMkUAPIgTwDIgzwBKH0DCnnjxxxzTPTv3z+2b9/e6frt27dHQ0NDl8c0NDQk7V9TUxM1NT53EqCS9UWeRMgUgEonTwDIi9e8AIgo8DtYqqurY/LkybFs2bKO69rb22PZsmXR2NjY5TGNjY2d9o+IWLp0abf7A1D55AkAeZAnAORFpgAQUeB3sERELFiwIK666qo455xz4rzzzouvf/3rsXv37rj66qsjImLOnDlx3HHHRVNTU0REXHfddTFt2rS47bbbYtasWXHvvffG888/H3feeWehRwWghMkTAPIgTwDIi0wBoOAFyxVXXBG//vWv44tf/GI0NzfHxIkT47HHHuv4Uq8tW7ZEv37vvJFm6tSpsWTJkvj7v//7+Lu/+7s45ZRT4qGHHoozzjij0KMCUMLkCQB5kCcA5EWmAFCVZVlW7CHy1NraGvX19THqmzdHv0G1xR6HHjj1k88XewR66ZFXVxd7BHqodWd7DBu3OVpaWqKurq7Y4xTd25myY/2YqBtc0E/VpEBe2rOv2CPQS3f99v3FHoEe2rNrb/zrtPtlSryTJ3NXfiiqjx5Y7HHogXW/HVXsEeilo/+6f7FHoIf27W+LZb+8XZ7EO3ly8SmfiwH9fTdLOXr0ifuLPQK9NPaxa4s9Aj3U/uZb8cpnb849T7xaBAAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkEjBAgAAAAAAkKhPCpZFixbFCSecELW1tTFlypR47rnnut138eLFUVVV1Wmrra3tizEBKHHyBIA8yBMA8iJTAI5sBS9Y7rvvvliwYEHcdNNNsWbNmpgwYULMnDkzduzY0e0xdXV1sW3bto5t8+bNhR4TgBInTwDIgzwBIC8yBYCCFyxf+9rXYu7cuXH11VfH+973vrjjjjviqKOOirvuuqvbY6qqqqKhoaFjGz58eKHHBKDEyRMA8iBPAMiLTAGgoAXLnj17YvXq1TFjxox3fmG/fjFjxoxYtWpVt8ft2rUrxowZE6NHj44PfvCD8dJLLxVyTABKnDwBIA/yBIC8yBQAIiIGFPLGf/Ob38T+/fsPaOOHDx8ev/jFL7o8Zty4cXHXXXfFWWedFS0tLXHrrbfG1KlT46WXXopRo0YdsH9bW1u0tbV1XG5tbY2IiBdm/HvUDe6Tr5ghZzOmf7LYI9BLL+15ttgj0EO79rQXe4Qu9UWeRHSfKX8+bmIMqBqY02roS/unn13sEeilPQtfL/YI9NC+3W2H3qmPFTtP/mLoqjjac5Sy9Dd/85lij0AvPfpE9+8ooLS17myP95xa7CkOVMzXvHb90/4Y8Ef7c1wNfWXGR73mVe5OXfF8sUegh/Zle+OVAtxuyf1139jYGHPmzImJEyfGtGnT4oc//GEce+yx8e1vf7vL/ZuamqK+vr5jGz16dB9PDEApSs2TCJkCwIHkCQB58ZoXQOUpaMFyzDHHRP/+/WP79u2drt++fXs0NDQc1m0MHDgwJk2aFBs2bOjy5wsXLoyWlpaObevWrb2eG4DS0hd5EiFTACqdPAEgL17zAiCiwAVLdXV1TJ48OZYtW9ZxXXt7eyxbtiwaGxsP6zb2798fL7zwQowYMaLLn9fU1ERdXV2nDYDK0hd5EiFTACqdPAEgL17zAiCiwN/BEhGxYMGCuOqqq+Kcc86J8847L77+9a/H7t274+qrr46IiDlz5sRxxx0XTU1NERFxyy23xPnnnx8nn3xyvPHGG/HVr341Nm/eHNdee22hRwWghMkTAPIgTwDIi0wBoOAFyxVXXBG//vWv44tf/GI0NzfHxIkT47HHHuv4ErAtW7ZEv37vvJHm9ddfj7lz50Zzc3O85z3vicmTJ8fTTz8d73vf+wo9KgAlTJ4AkAd5AkBeZAoAVVmWZcUeIk+tra1RX18fr//PiVE3uKCfgEaBzPjoJ4s9Ar10693fKvYI9NCune0x7cxXo6WlxdvP451MmR4fjAFVA4s9Dj2wf/rZxR6BXtqz8PVij0AP7dvdFs/+6T/LlHgnT1a+cFwc7TlKWfqbqz9T7BHopceX3FXsEeih1p3t8Z5TfyVP4p08mfLgX8eAP6op9jj0QHXTe4o9Ar3Uf8WaYo9AD+3L9saKeDj3PPHXPQAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQCIFCwAAAAAAQKKCFixPPvlkXH755TFy5MioqqqKhx566JDHrFixIs4+++yoqamJk08+ORYvXlzIEQEoA/IEgLzIFADyIE8AiChwwbJ79+6YMGFCLFq06LD237RpU8yaNSsuuuiiWLduXVx//fVx7bXXxk9+8pNCjglAiZMnAORFpgCQB3kCQETEgELe+KWXXhqXXnrpYe9/xx13xNixY+O2226LiIjTTjstnnrqqbj99ttj5syZhRoTgBInTwDIi0wBIA/yBICIEvsOllWrVsWMGTM6XTdz5sxYtWpVkSYCoBzJEwDyIlMAyIM8AahMBX0HS6rm5uYYPnx4p+uGDx8era2t8eabb8agQYMOOKatrS3a2to6Lre2thZ8TgBKW0/yJEKmAHAgz1EAyIM8AahMJfUOlp5oamqK+vr6jm306NHFHgmAMiVTAMiDPAEgD/IEoPSVVMHS0NAQ27dv73Td9u3bo66urtt/bbxw4cJoaWnp2LZu3doXowJQwnqSJxEyBYADeY4CQB7kCUBlKqmPCGtsbIxHH32003VLly6NxsbGbo+pqamJmpqaQo8GQBnpSZ5EyBQADuQ5CgB5kCcAlamg72DZtWtXrFu3LtatWxcREZs2bYp169bFli1bIuL3TfycOXM69v/0pz8dv/rVr+ILX/hC/OIXv4hvfvOb8f3vfz8+97nPFXJMAEqcPAEgLzIFgDzIEwAiClywPP/88zFp0qSYNGlSREQsWLAgJk2aFF/84hcjImLbtm0dwRMRMXbs2PjRj34US5cujQkTJsRtt90W3/nOd2LmzJmFHBOAEidPAMiLTAEgD/IEgIgCf0TY9OnTI8uybn++ePHiLo9Zu3ZtAacCoNzIEwDyIlMAyIM8ASCixL7kHgAAAAAAoBwoWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIpWAAAAAAAABIVtGB58skn4/LLL4+RI0dGVVVVPPTQQwfdf8WKFVFVVXXA1tzcXMgxAShx8gSAvMgUAPIgTwCIKHDBsnv37pgwYUIsWrQo6bj169fHtm3bOrZhw4YVaEIAyoE8ASAvMgWAPMgTACIiBhTyxi+99NK49NJLk48bNmxYDBkyJP+BAChL8gSAvMgUAPIgTwCIKHDB0lMTJ06Mtra2OOOMM+Lmm2+OCy64oNt929raoq2treNya2trRER8aNYHY0D/moLPSv72/PPrxR6BXrrrt+8v9gj00J5deyPi/mKPkZuUPInoPlM2fHNS9BtUW9BZKYxNH/hOsUegly676EPFHoEe2re/f7FHyFUez1FOqx4QddW+BrMc9V+xptgj0EtjH7u22CPQQ+1vvhURNxd7jNzkkSdnDf2/UX30wILPSv4+efe9xR6BXjq9uiRfTucwtO5sj2Hj8r/dkvrrfsSIEXHHHXfEAw88EA888ECMHj06pk+fHmvWdP/HbFNTU9TX13dso0eP7sOJAShFPcmTCJkCwIE8RwEgD/IEoDKVVOU2bty4GDfunRpp6tSpsXHjxrj99tvju9/9bpfHLFy4MBYsWNBxubW1VeAAHOF6kicRMgWAA3mOAkAe5AlAZSqpgqUr5513Xjz11FPd/rympiZqanwUGAAHd6g8iZApABwez1EAyIM8ASh/JfURYV1Zt25djBgxothjAFDm5AkAeZEpAORBngCUv4K+g2XXrl2xYcOGjsubNm2KdevWxdChQ+P444+PhQsXxquvvhr//u//HhERX//612Ps2LFx+umnx1tvvRXf+c53Yvny5fFf//VfhRwTgBInTwDIi0wBIA/yBICIAhcszz//fFx00UUdl9/+3MirrroqFi9eHNu2bYstW7Z0/HzPnj1xww03xKuvvhpHHXVUnHXWWfH44493ug0AjjzyBIC8yBQA8iBPAIiIqMqyLCv2EHlqbW2N+vr6uPiUz8WA/j6nshzt+uf9xR6BXpr43leKPQI9tGfX3vjXafdHS0tL1NXVFXucons7U0Z98+boN6i22OPQA5s+8J1ij0AvXXbRh4o9Aj20b39bLPvl7TIl3smTHevHRN3gkv+UZrrwJ8dNLvYI9NL/3HVOsUegh9rffCte+ezN8iTeyZO5Kz8U1UcPLPY49MAn33vw7wWl9J1eXfJfaU43Wne2x7Bxm3PPE3/dAwAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJFKwAAAAAAAAJCpowdLU1BTnnntuDB48OIYNGxazZ8+O9evXH/K4H/zgBzF+/Piora2NM888Mx599NFCjglAiZMnAORBngCQF5kCQESBC5aVK1fGvHnz4plnnomlS5fG3r1745JLLondu3d3e8zTTz8dV155ZVxzzTWxdu3amD17dsyePTtefPHFQo4KQAmTJwDkQZ4AkBeZAkBERFWWZVlf/bJf//rXMWzYsFi5cmVceOGFXe5zxRVXxO7du+ORRx7puO7888+PiRMnxh133HHI39Ha2hr19fVx8SmfiwH9a3Kbnb6z65/3F3sEemnie18p9gj00J5de+Nfp90fLS0tUVdXV+xxutUXeRLxTqaM+ubN0W9QbS6z07c2feA7xR6BXrrsog8VewR6aN/+tlj2y9tLOlP6Ok92rB8TdYN9SnM5+pPjJhd7BHrpf+46p9gj0EPtb74Vr3z25pLOk4i+fc1r7soPRfXRA3Obnb7zyfc+VewR6KXTqwcUewR6qHVnewwbtzn3POnTv+5bWloiImLo0KHd7rNq1aqYMWNGp+tmzpwZq1at6nL/tra2aG1t7bQBUNkKkScRMgXgSCNPAMiL17wAjkx9VrC0t7fH9ddfHxdccEGcccYZ3e7X3Nwcw4cP73Td8OHDo7m5ucv9m5qaor6+vmMbPXp0rnMDUFoKlScRMgXgSCJPAMiL17wAjlx9VrDMmzcvXnzxxbj33ntzvd2FCxdGS0tLx7Z169Zcbx+A0lKoPImQKQBHEnkCQF685gVw5OqTD42bP39+PPLII/Hkk0/GqFGjDrpvQ0NDbN++vdN127dvj4aGhi73r6mpiZoa37UCcCQoZJ5EyBSAI4U8ASAvXvMCOLIV9B0sWZbF/Pnz48EHH4zly5fH2LFjD3lMY2NjLFu2rNN1S5cujcbGxkKNCUCJkycA5EGeAJAXmQJARIHfwTJv3rxYsmRJPPzwwzF48OCOz5Ssr6+PQYMGRUTEnDlz4rjjjoumpqaIiLjuuuti2rRpcdttt8WsWbPi3nvvjeeffz7uvPPOQo4KQAmTJwDkQZ4AkBeZAkBEgd/B8q1vfStaWlpi+vTpMWLEiI7tvvvu69hny5YtsW3bto7LU6dOjSVLlsSdd94ZEyZMiPvvvz8eeuihg35JGACVTZ4AkAd5AkBeZAoAEQV+B0uWZYfcZ8WKFQdc9+EPfzg+/OEPF2AiAMqRPAEgD/IEgLzIFAAiCvwOFgAAAAAAgEqkYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEikYAEAAAAAAEhU0IKlqakpzj333Bg8eHAMGzYsZs+eHevXrz/oMYsXL46qqqpOW21tbSHHBKDEyRMA8iBPAMiLTAEgosAFy8qVK2PevHnxzDPPxNKlS2Pv3r1xySWXxO7duw96XF1dXWzbtq1j27x5cyHHBKDEyRMA8iBPAMiLTAEgImJAIW/8scce63R58eLFMWzYsFi9enVceOGF3R5XVVUVDQ0NhRwNgDIiTwDIgzwBIC8yBYCIAhcs79bS0hIREUOHDj3ofrt27YoxY8ZEe3t7nH322fGP//iPcfrpp3e5b1tbW7S1tR3wO/btb+tyf0rfvt37iz0CvbSnZm+xR6CH9uz+/X2XZVmRJzm4QuRJRPeZ0v7mWzlMTTG07mwv9gj0kr/pytfb910pZ0pf58nOXf6bVK72Zf6+LXf+nitfb993pZwnEX37mtfbz9soP7uq/S1Q7lrdh2Xr7b/Fc8+TrI/s378/mzVrVnbBBRccdL+nn346+7d/+7ds7dq12YoVK7I/+ZM/yerq6rKtW7d2uf9NN92URYTNZrPZctw2btxYiCjIRaHyJMtkis1msxViK9VMkSc2m81WXlup5kmWec3LZrPZymnLO0+qsqxv/gnAZz7zmfjxj38cTz31VIwaNeqwj9u7d2+cdtppceWVV8aXv/zlA37+7jb/jTfeiDFjxsSWLVuivr4+l9lLSWtra4wePTq2bt0adXV1xR6nICp9jdZX3ip9fS0tLXH88cfH66+/HkOGDCn2OF0qVJ5EyJRKY33lzfrKX6lnijzJT6Wfz9ZX/ip9jZW+vlLPkwiveeWl0s/liMpfo/WVt0pfX6HypE8+Imz+/PnxyCOPxJNPPpkUNBERAwcOjEmTJsWGDRu6/HlNTU3U1NQccH19fX1Fnghvq6urq+j1RVT+Gq2vvFX6+vr161fsEbpUyDyJkCmVyvrKm/WVv1LMFHlSGJV+Pltf+av0NVb6+koxTyK85lUIlX4uR1T+Gq2vvFX6+vLOk4KmU5ZlMX/+/HjwwQdj+fLlMXbs2OTb2L9/f7zwwgsxYsSIAkwIQDmQJwDkQZ4AkBeZAkBEgd/BMm/evFiyZEk8/PDDMXjw4Ghubo6I3zftgwYNioiIOXPmxHHHHRdNTU0REXHLLbfE+eefHyeffHK88cYb8dWvfjU2b94c1157bSFHBaCEyRMA8iBPAMiLTAEgosAFy7e+9a2IiJg+fXqn6+++++74xCc+ERERW7Zs6fS2nNdffz3mzp0bzc3N8Z73vCcmT54cTz/9dLzvfe87rN9ZU1MTN910U5dvoawElb6+iMpfo/WVN+srjmLkSUTp/v+RF+srb9ZX3ip9fRGluUZ5UhjWV94qfX0Rlb9G6ysOr3nlr9LXF1H5a7S+8mZ9PdNnX3IPAAAAAABQKUrzG8IAAAAAAABKmIIFAAAAAAAgkYIFAAAAAAAgkYIFAAAAAAAgUUUULK+99lp87GMfi7q6uhgyZEhcc801sWvXroMeM3369Kiqquq0ffrTn+6jiQ9u0aJFccIJJ0RtbW1MmTIlnnvuuYPu/4Mf/CDGjx8ftbW1ceaZZ8ajjz7aR5P2TMr6Fi9efMD9VFtb24fTpnnyySfj8ssvj5EjR0ZVVVU89NBDhzxmxYoVcfbZZ0dNTU2cfPLJsXjx4oLP2VOp61uxYsUB919VVVU0Nzf3zcCJmpqa4txzz43BgwfHsGHDYvbs2bF+/fpDHlcuj8GerK/cHoO9VWl5EiFT/lC5nc8ypbNyypRKz5MImXI4Ki1T5Mk7yu1cliedlVOeRFR+psiTQ5Mn5XEu/6FKzRR50pk8KS3FzJOKKFg+9rGPxUsvvRRLly6NRx55JJ588sn41Kc+dcjj5s6dG9u2bevY/umf/qkPpj24++67LxYsWBA33XRTrFmzJiZMmBAzZ86MHTt2dLn/008/HVdeeWVcc801sXbt2pg9e3bMnj07XnzxxT6e/PCkri8ioq6urtP9tHnz5j6cOM3u3btjwoQJsWjRosPaf9OmTTFr1qy46KKLYt26dXH99dfHtddeGz/5yU8KPGnPpK7vbevXr+90Hw4bNqxAE/bOypUrY968efHMM8/E0qVLY+/evXHJJZfE7t27uz2mnB6DPVlfRHk9BnurkvIkQqZ0pZzOZ5nStXLIlErPkwiZcjgqKVPkyYHK6VyWJ10rhzyJqPxMkSeHJk/K41x+WyVnijzpmjwpDUXNk6zM/fznP88iIvvZz37Wcd2Pf/zjrKqqKnv11Ve7PW7atGnZdddd1wcTpjnvvPOyefPmdVzev39/NnLkyKypqanL/f/8z/88mzVrVqfrpkyZkv3lX/5lQefsqdT13X333Vl9fX0fTZeviMgefPDBg+7zhS98ITv99NM7XXfFFVdkM2fOLOBk+Tic9T3xxBNZRGSvv/56n8yUtx07dmQRka1cubLbfcrtMfiHDmd95fwYTFVpeZJlMuXdyvl8linlnSmVnidZJlPerdIyRZ50Vs7nsjwp7zzJssrPFHnSmTwpv3P5SMkUeSJPSl1f5knZv4Nl1apVMWTIkDjnnHM6rpsxY0b069cvnn322YMe+73vfS+OOeaYOOOMM2LhwoXxu9/9rtDjHtSePXti9erVMWPGjI7r+vXrFzNmzIhVq1Z1ecyqVas67R8RMXPmzG73L6aerC8iYteuXTFmzJgYPXp0fPCDH4yXXnqpL8btE+V0//XGxIkTY8SIEfHHf/zH8dOf/rTY4xy2lpaWiIgYOnRot/uU8314OOuLqOzH4B+qpDyJkCndqeTzuZzuv94ox0yp9DyJkCnvVkmZIk+6Vsnncjndf71RjnkSUfmZIk86kyfldS7LlM7K7f7rKXlSmvoyT8q+YGlubj7grVcDBgyIoUOHHvQz7z760Y/GPffcE0888UQsXLgwvvvd78Zf/MVfFHrcg/rNb34T+/fvj+HDh3e6fvjw4d2upbm5OWn/YurJ+saNGxd33XVXPPzww3HPPfdEe3t7TJ06NV555ZW+GLngurv/Wltb48033yzSVPkZMWJE3HHHHfHAAw/EAw88EKNHj47p06fHmjVrij3aIbW3t8f1118fF1xwQZxxxhnd7ldOj8E/dLjrq/TH4B+qpDyJkCldqfTzWaaUpkrPkwiZ0pVKyhR5cqBKP5flSemq9EyRJweSJ+V1LsuUzuRJ6ZInv5fX429AbwculBtvvDG+8pWvHHSfl19+uce3/4efV3nmmWfGiBEj4uKLL46NGzfGSSed1OPbJV+NjY3R2NjYcXnq1Klx2mmnxbe//e348pe/XMTJOBzjxo2LcePGdVyeOnVqbNy4MW6//fb47ne/W8TJDm3evHnx4osvxlNPPVXsUQricNdXCY9BecLbKuF8PpKVa6ZUep5EyJR3kymVrxLO5SNZueZJROVnijzpTJ4cGSrhfD5SyZPS1dd5UrIFyw033BCf+MQnDrrPiSeeGA0NDQd8UdS+ffvitddei4aGhsP+fVOmTImIiA0bNhQtbI455pjo379/bN++vdP127dv73YtDQ0NSfsXU0/W924DBw6MSZMmxYYNGwoxYp/r7v6rq6uLQYMGFWmqwjrvvPNK/j/g8+fP7/jywFGjRh1033J6DL4tZX3vVo6PwSMxTyJkyuEox/P5YGRK6an0PImQKV2ptEyRJ4dWjufywciT0lTpmSJPDiRPyu9climdyZPSJE+619PHX8l+RNixxx4b48ePP+hWXV0djY2N8cYbb8Tq1as7jl2+fHm0t7d3BMjhWLduXUT8/u1dxVJdXR2TJ0+OZcuWdVzX3t4ey5Yt69Sm/aHGxsZO+0dELF26tNv9i6kn63u3/fv3xwsvvFDU+ylP5XT/5WXdunUle/9lWRbz58+PBx98MJYvXx5jx4495DHldB/2ZH3vVo6PwSMxTyJkyuEox/P5YMrp/stLqWZKpedJhEw5kjJFnhxaOZ7LB1NO919eSjVPIio/U+SJPKmUPImQKe9WbvdfHuRJ8RQ1T7r44vuy84EPfCCbNGlS9uyzz2ZPPfVUdsopp2RXXnllx89feeWVbNy4cdmzzz6bZVmWbdiwIbvllluy559/Ptu0aVP28MMPZyeeeGJ24YUXFmsJHe69996spqYmW7x4cfbzn/88+9SnPpUNGTIka25uzrIsyz7+8Y9nN954Y8f+P/3pT7MBAwZkt956a/byyy9nN910UzZw4MDshRdeKNYSDip1fV/60peyn/zkJ9nGjRuz1atXZx/5yEey2tra7KWXXirWEg5q586d2dq1a7O1a9dmEZF97Wtfy9auXZtt3rw5y7Isu/HGG7OPf/zjHfv/6le/yo466qjs85//fPbyyy9nixYtyvr375899thjxVrCQaWu7/bbb88eeuih7Je//GX2wgsvZNddd13Wr1+/7PHHHy/WEg7qM5/5TFZfX5+tWLEi27ZtW8f2u9/9rmOfcn4M9mR95fYY7K1KypMskynlfj7LlPLNlErPkyyTKYejkjJFnpT3uSxPyjdPsqzyM0WeHJo8KY9z+W2VnCnyRJ6U8mOwmHlSEQXLb3/72+zKK6/Mjj766Kyuri67+uqrs507d3b8fNOmTVlEZE888USWZVm2ZcuW7MILL8yGDh2a1dTUZCeffHL2+c9/PmtpaSnSCjr7xje+kR1//PFZdXV1dt5552XPPPNMx8+mTZuWXXXVVZ32//73v5+deuqpWXV1dXb66adnP/rRj/p44jQp67v++us79h0+fHh22WWXZWvWrCnC1IfniSeeyCLigO3tNV111VXZtGnTDjhm4sSJWXV1dXbiiSdmd999d5/PfbhS1/eVr3wlO+mkk7La2tps6NCh2fTp07Ply5cXZ/jD0NXaIqLTfVLOj8GerK/cHoO9VWl5kmUypZzPZ5lSvplS6XmSZTLlcFRapsiTqzoul9u5LE/KN0+yrPIzRZ4cmjwpj3P5D1VqpsgTeVLKj8Fi5knV/x8AAAAAAACAw1Sy38ECAAAAAABQqhQsAAAAAAAAiRQsAAAAAAAAiRQsAAAAAAAAiRQsAAAAAAAAiRQsAAAAAAAAiRQsAAAAAAAAiRQsAAAAAAAAiRQsAAAAAAAAiRQsAAAAAAAAiRQsAAAAAAAAiRQsAAAAAAAAif4f/4ppu/hcT5MAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 2000x500 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# C4 rotations: e, r, r^2, r^3\n",
    "x = torch.rand(1,1,3,3)\n",
    "fig=plt.figure(figsize=(20, 5))\n",
    "for i in range(1, 5):\n",
    "    fig.add_subplot(1, 4, i)\n",
    "    plt.imshow(rot_img(x, np.pi/2*i)[0,0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vDNYJxKIPn9b"
   },
   "source": [
    "## 1.3 Implementing the Lift Convolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sWpduUdeCEYR"
   },
   "source": [
    "$$f_{out}(\\mathbf{x}, r) = (f_{in} * \\psi) (\\mathbf{x})\n",
    "= \\sum_{\\mathbf{y} \\in \\mathbb{Z}^2}  f_{in}(\\mathbf{y})\\psi_{r}(\\mathbf{y} -\\mathbf{x})$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "d2IlLWUZMsyd"
   },
   "source": [
    "**In lift convolution, the input convolves with a stack of rotated filters.**\n",
    "\n",
    "*The reason why we call it lifting is that we add a new axis $r$ in the output feature maps.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Icc31YdfEcXn"
   },
   "source": [
    "![picture](https://drive.google.com/uc?export=view&id=17RprAkuSe43q_gjsi1al2RX18mZ0o00q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "id": "j9UENyEg0ucj",
    "tags": []
   },
   "outputs": [],
   "source": [
    "class LiftingConvolution(nn.Module):\n",
    "    \"\"\"Lifting Convolution Layer for finite rotation group\n",
    "\n",
    "    Attributes:\n",
    "        in_channels: number of input channels\n",
    "        out_channels: number of output channels\n",
    "        kernel_size: kernel size\n",
    "        group_order: the order of rotation groups\n",
    "        activation: whether to use relu.\n",
    "    \"\"\"\n",
    "    def __init__(self,\n",
    "                 in_channels,\n",
    "                 out_channels,\n",
    "                 kernel_size,\n",
    "                 group_order,\n",
    "                 activation = True\n",
    "                 ):\n",
    "        super(LiftingConvolution, self).__init__()\n",
    "\n",
    "        self.in_channels = in_channels\n",
    "        self.out_channels = out_channels\n",
    "        self.kernel_size = kernel_size\n",
    "        self.group_order = group_order\n",
    "        self.activation = activation\n",
    "\n",
    "        # Initialize an unconstrained kernel.\n",
    "        self.kernel = torch.nn.Parameter(torch.zeros(self.out_channels,\n",
    "                                                     self.in_channels,\n",
    "                                                     self.kernel_size,\n",
    "                                                     self.kernel_size))\n",
    "\n",
    "        # Initialize weights\n",
    "        torch.nn.init.kaiming_uniform_(self.kernel.data, a=math.sqrt(5))\n",
    "\n",
    "    def generate_filter_bank(self):\n",
    "        # Obtain a stack of rotated filters\n",
    "        # Rotate kernels by 0, 90, 180, and 270 degrees\n",
    "        # ==============================\n",
    "        filter_bank = torch.stack([rot_img(self.kernel, -np.pi*2/self.group_order*i)\n",
    "                                   for i in range(self.group_order)])\n",
    "        # ==============================\n",
    "\n",
    "        # [#out, group_order, #in, k, k]\n",
    "        filter_bank = filter_bank.transpose(0,1)\n",
    "        return filter_bank\n",
    "\n",
    "\n",
    "    def forward(self, x):\n",
    "        # input shape: [bz, #in, h, w]\n",
    "        # output shape: [bz, #out, group order, h, w]\n",
    "\n",
    "        # generate filter bank given input group order\n",
    "        filter_bank = self.generate_filter_bank()\n",
    "\n",
    "        # concatenate the first two dims before convolution.\n",
    "        # ==============================\n",
    "        x = F.conv2d(\n",
    "            input=x,\n",
    "            weight=filter_bank.reshape(\n",
    "                self.out_channels * self.group_order,\n",
    "                self.in_channels,\n",
    "                self.kernel_size,\n",
    "                self.kernel_size\n",
    "            ),\n",
    "            padding = (self.kernel_size-1)//2\n",
    "        )\n",
    "        # ==============================\n",
    "\n",
    "        # reshape output signal to shape [bz, #out, group order, h, w].\n",
    "        # ==============================\n",
    "        x = x.view(\n",
    "            x.shape[0],\n",
    "            self.out_channels,\n",
    "            self.group_order,\n",
    "            x.shape[-1],\n",
    "            x.shape[-2]\n",
    "        )\n",
    "        # ==============================\n",
    "        if self.activation:\n",
    "            return F.leaky_relu(x)\n",
    "        return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AQKnUvTTQhZx"
   },
   "source": [
    "## 1.4 Implementing Group Convolution Kernel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "B__anHOmzvHD"
   },
   "source": [
    "**After the first layer, both the filter and input are now functions on $G=\\mathbb{Z}^2 \\rtimes C_4$**\n",
    "\n",
    "*Be careful about this, rotation not only acts on the spatial dimension but also the new $r$ axis. A rotation means a planar rotation and periodic shft.*\n",
    "\n",
    "*A single channel from regular CNNs is expanded to four channels since we apply kernel with rotations four times.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Lj36BPj1CLaO"
   },
   "source": [
    "$$f_{out}(\\mathbf{x}, r) = (f_{in} * \\psi) (\\mathbf{x}, r)\n",
    "    = \\sum_{\\mathbf{y} \\in \\mathbb{Z}^2} \\sum_{r' \\in C^4}  f_{in}(\\mathbf{y}, r') \\psi(r^{-1}(\\mathbf{y-x}),r^{-1}r')\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "PZ_bF30ScLwX"
   },
   "source": [
    "![picture](https://drive.google.com/uc?export=view&id=1ADYVLN8veOHRgK_l50p7Vw3-Nvn2E0P9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "id": "dzTItCkechNT",
    "tags": []
   },
   "outputs": [],
   "source": [
    "class GroupConvolution(nn.Module):\n",
    "    \"\"\"Group Convolution Layer for finite rotation group\n",
    "\n",
    "    Attributes:\n",
    "        in_channels: number of input channels\n",
    "        out_channels: number of output channels\n",
    "        kernel_size: kernel size\n",
    "        group_order: the order of rotation groups\n",
    "        activation: whether to use relu.\n",
    "    \"\"\"\n",
    "    def __init__(self,\n",
    "                 in_channels,\n",
    "                 out_channels,\n",
    "                 kernel_size,\n",
    "                 group_order,\n",
    "                 activation = True,\n",
    "                 ):\n",
    "        super(GroupConvolution, self).__init__()\n",
    "        self.in_channels = in_channels\n",
    "        self.out_channels = out_channels\n",
    "        self.kernel_size = kernel_size\n",
    "        self.group_order = group_order\n",
    "        self.activation = activation\n",
    "\n",
    "        # Initialize an unconstrained kernel.\n",
    "        # the weights have an additional group order dimension.\n",
    "        self.kernel = torch.nn.Parameter(torch.zeros(self.out_channels,\n",
    "                                                     self.in_channels,\n",
    "                                                     self.group_order, # this is different from the lifting convolution\n",
    "                                                     self.kernel_size,\n",
    "                                                     self.kernel_size))\n",
    "\n",
    "        torch.nn.init.kaiming_uniform_(self.kernel.data, a=math.sqrt(5))\n",
    "\n",
    "\n",
    "    def generate_filter_bank(self):\n",
    "        # Obtain a stack of rotated and cyclic shifted filters\n",
    "        filter_bank = []\n",
    "        filter = self.kernel.reshape(self.out_channels*self.in_channels,\n",
    "                                     self.group_order,\n",
    "                                     self.kernel_size,\n",
    "                                     self.kernel_size)\n",
    "\n",
    "        for i in range(self.group_order):\n",
    "            # planar rotation\n",
    "            rotated_filter = rot_img(filter, -np.pi*2/self.group_order*i)\n",
    "\n",
    "            # cyclic shift\n",
    "            shifted_indices = torch.roll(torch.arange(0, self.group_order, 1), shifts = i)\n",
    "            shifted_rotated_filter = rotated_filter[:,shifted_indices]\n",
    "\n",
    "\n",
    "            filter_bank.append(shifted_rotated_filter.reshape(self.out_channels,\n",
    "                                                              self.in_channels,\n",
    "                                                              self.group_order,\n",
    "                                                              self.kernel_size,\n",
    "                                                              self.kernel_size))\n",
    "        # reshape output signal to shape [#out, g_order, #in, g_order, k, k].\n",
    "        filter_bank = torch.stack(filter_bank).transpose(0,1)\n",
    "        return filter_bank\n",
    "\n",
    "    def forward(self, x):\n",
    "        # input shape: [bz, in, group order, x, y]\n",
    "        # output shape: [bz, out, group order, x, y]\n",
    "\n",
    "        # Generate filter bank with shape [#out, g_order, #in, g_order, h, w]\n",
    "        filter_bank = self.generate_filter_bank()\n",
    "\n",
    "        # Reshape filter_bank to use F.conv2d\n",
    "        # [#out, g_order, #in, g_order, h, w] -> [#out*g_order, #in*g_order, h, w]\n",
    "        # ==============================\n",
    "        x = torch.nn.functional.conv2d(\n",
    "            input=x.reshape(\n",
    "                x.shape[0],\n",
    "                x.shape[1] * x.shape[2],\n",
    "                x.shape[3],\n",
    "                x.shape[4]\n",
    "                ),\n",
    "            weight=filter_bank.reshape(\n",
    "                self.out_channels * self.group_order,\n",
    "                self.in_channels * self.group_order,\n",
    "                self.kernel_size,\n",
    "                self.kernel_size\n",
    "            ),\n",
    "            padding = (self.kernel_size-1)//2\n",
    "        )\n",
    "\n",
    "        # Reshape signal back [bz, #out * g_order, h, w] -> [bz, out, g_order, h, w]\n",
    "        x = x.view(x.shape[0], self.out_channels, self.group_order, x.shape[-2], x.shape[-1])\n",
    "        # ========================\n",
    "\n",
    "        if self.activation:\n",
    "            return F.leaky_relu(x)\n",
    "        return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "hLppDfHeQx8_"
   },
   "source": [
    "## 1.5 Test Equivariance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GZA_6iEseHZA"
   },
   "source": [
    "$$ \\|f\\|_{EE} = sup_{x,g}\\|f(\\rho_X(g)x) - \\rho_Y(g)f(x)\\| = 0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UR6gszpjKT1Z"
   },
   "source": [
    "**Let's build a group equivariant CNNs.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "id": "skRK8RplcuTm",
    "tags": []
   },
   "outputs": [],
   "source": [
    "class GroupEquivariantCNN(torch.nn.Module):\n",
    "    def __init__(self,\n",
    "                 in_channels,\n",
    "                 out_channels,\n",
    "                 kernel_size,\n",
    "                 hidden_dim,\n",
    "                 group_order,\n",
    "                 num_gconvs, # number of group convolution layers.\n",
    "                 classifer = False,\n",
    "                 sigmoid = False\n",
    "                 ):\n",
    "        super().__init__()\n",
    "\n",
    "        self.gconvs = []\n",
    "        # First Layer\n",
    "        self.gconvs.append(LiftingConvolution(in_channels = in_channels,\n",
    "                                               out_channels = hidden_dim,\n",
    "                                               kernel_size = kernel_size,\n",
    "                                               group_order = group_order,\n",
    "                                               activation = True))\n",
    "        # Middle Layers\n",
    "        \n",
    "        for i in range(num_gconvs-2):\n",
    "            self.gconvs.append(GroupConvolution(in_channels = hidden_dim,\n",
    "                                                out_channels = hidden_dim,\n",
    "                                                kernel_size = kernel_size,\n",
    "                                                group_order = group_order,\n",
    "                                                activation = True))\n",
    "    \n",
    "        # Final Layer # To generate equivariant outputs\n",
    "        self.gconvs.append(GroupConvolution(in_channels = hidden_dim,\n",
    "                                            out_channels = out_channels,\n",
    "                                            kernel_size = kernel_size,\n",
    "                                            group_order = group_order,\n",
    "                                            activation = False))\n",
    "        self.gconvs = nn.Sequential(*self.gconvs)\n",
    "        self.classifer = classifer\n",
    "        self.sigmoid = sigmoid\n",
    "        \n",
    "    def forward(self, x):\n",
    "        # x [bz, c_in, h, w]\n",
    "        out = self.gconvs(x)\n",
    "        \n",
    "        # functions on (g,x,y) -> functions on (x,y)\n",
    "        # [bz, c_out, |G|, H, W] -> [bz, c_out, H, W]\n",
    "        out = torch.mean(out, dim = 2)\n",
    "        \n",
    "        # If we want to have a invariant classifer, we can average over the spatial dimensions.\n",
    "        if self.classifer:\n",
    "            out = torch.mean(out, dim = (2,3))\n",
    "        if self.sigmoid:\n",
    "            out = out.sigmoid()\n",
    "        return out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test whether the model is equivariant under C4 transformations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def test_equivariance(model, group_order, plot = False):\n",
    "    # Download mnist dataset.\n",
    "    train_ds = torchvision.datasets.MNIST(root=\"data\", train=True, download=False)\n",
    "    image = torch.from_numpy(train_ds.data[0].numpy()).unsqueeze(0).unsqueeze(0).float()\n",
    "    image = (image - image.min())/(image.max() - image.min())\n",
    "\n",
    "\n",
    "    ## |rf(x) - f(rx)|\n",
    "    # test whether the model is equivariant under C4 transformations.\n",
    "    for i in range(group_order):\n",
    "        theta = np.pi*2/group_order*i\n",
    "        x = image\n",
    "        f_x = model(x)\n",
    "        r_f_x = rot_img(f_x, theta).data.numpy()\n",
    "        r_x = rot_img(x, theta)\n",
    "        f_r_x = model(r_x).data.numpy()\n",
    "        equiv_error = np.mean(np.abs(r_f_x - f_r_x))\n",
    "        print(\"Equivariance Error with Rotation {:0.0f} Degrees: {:0.5f}\".format(theta/np.pi*2/group_order*360, equiv_error))\n",
    "        \n",
    "    if plot:\n",
    "        titles = [\"$x$\", \"$f(x)$\", \"$gx$\", \"$f(gx)$\", \"$g(f(x))$\"]\n",
    "        vis_imgs = [x[0,0], f_x[0,0].data.numpy(), r_x[0,0].data.numpy(), f_r_x[0,0], r_f_x[0,0]]\n",
    "        fig=plt.figure(figsize=(10, 5))\n",
    "        idx = 0\n",
    "        for i in range(6):\n",
    "            if i == 1:\n",
    "                continue\n",
    "            fig.add_subplot(2, 3, i+1)\n",
    "            plt.imshow(vis_imgs[idx])\n",
    "            plt.title(titles[idx], fontsize = 20)\n",
    "            plt.axis('off')\n",
    "            idx += 1\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "id": "Fp6RplwdcwNS",
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Equivariance Error with Rotation 0 Degrees: 0.00000\n",
      "Equivariance Error with Rotation 90 Degrees: 0.00000\n",
      "Equivariance Error with Rotation 180 Degrees: 0.00000\n",
      "Equivariance Error with Rotation 270 Degrees: 0.00000\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuYAAAG0CAYAAAB+N+IyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA89UlEQVR4nO3de3hU5bn//8+amUwOJCSEQwAxgBytgKhUVFpb259CtZ5QPOCpetWWotbuVmt/blu0P091W3c9VGu1Vve2flVorbpbrbT11HqCCgJWEDmJnCIEQshxMrN+f/glG8v9hKwwSR6Y9+u6uC78zKznWTNk7twss26CMAxDAQAAAOhWse4+AQAAAAA05gAAAIAXaMwBAAAAD9CYAwAAAB6gMQcAAAA8QGMOAAAAeIDGHAAAAPAAjTkAAADgARpzAAAAwAM05gAAAIAHaMwBAAD2A88++6xOPvlkDRw4UMlkUkEQKAgCLVmyRM3NzRoxYoSCINCcOXM6Zf/LLrtMQRDooosu6pT1c0EQhmHY3ScBAACAjvvhD3+oG2+8cbc8mUyqrq5O//mf/6nvf//7GjNmjBYtWqQgCLJ+Dh9++KFGjBihVCqlefPm6Ygjjsj6Hvs7GnMAAIB92MKFC3X44YcrDENNmjRJ3//+9zV48GDF43Hl5eVp4MCBGjp0qLZs2aLHH39cZ599dqedyze+8Q098MAD+spXvqI//vGPnbbP/orGHAAAYB/2rW99S7/4xS/Us2dPrV69Wr169frU47fddpuuueYaVVZWatWqVYrFOu8nmZctW6bRo0dLkubPn89V84j4GXMAAIB92AsvvCBJmjJlym5NeTqd1j333CNJOvfcczu1KZekUaNG6fDDD5ck3X333Z261/6IxhwAAGAftXnzZq1cuVKSNGnSpN0enzt3rtauXStJOu+887rknHbuM3v2bNXW1nbJnvsLGvMuVFdXp4qKCgVBoIMOOkipVMp8XkNDg4455hgFQaCCggK98sorXXymAADAZ+eff76CIFDfvn1bsyuvvLJ1EksQBHrsscf05JNPSpJGjBihsWPHOtfLZo9yxhlnSJLq6+v19NNP783LzDk05l2oR48euvbaayVJq1at0iOPPLLbczKZjM477zy9/vrrisVievTRR3Xsscd29akCAACPvffee3t8zmc+8xm9+OKLkqSjjjqqzedms0cZPHiw+vfvL0l67rnn9nie+F805l1sxowZqqyslCTddNNNu/2N9Dvf+Y6eeuopSdIdd9yhM888s8vPEQAA+O2xxx7T4sWLNXXqVElSRUWFFi9e/Klfffr00erVqyVJn/3sZ/e4ZjZ7lCOPPFKS9PLLL0d+bbmMxryL5efn60c/+pEkafXq1fr1r3/d+thPf/rT1hslrrrqKl155ZXdco4AAMBvo0aN0pgxY7Ru3TpJ0vjx4zVmzJhP/Xrttddan3/YYYftcc1s9ig7p7GsW7dOmzZtivbichiNeTf42te+ppEjR0qSbr75ZqVSKT355JO6+uqrJX1y1/Rtt93WnacIAAA8l8lktHjxYkmfNOb/6qOPPmr9fb9+/dq1ZrZ6lF3323lzKvaMxrwbxONx/fjHP5YkrVmzRjNnztSFF16oMAx13HHH6eGHH+6Uf5ELAADsP5YvX676+npJ0qGHHrrb4x9//HHr7/91jKJLtnqU8vLy1t9v3LixXXuDxrzbnHXWWa1/u33wwQfV1NSksWPH6qmnnlIymezekwMAAN5buHBh6++tK+bV1dWtv29vYy5lp0fZdb+6urp2753raMy7SRAEuvTSS1v/u3///nruuedUWlrajWcFAAD2FTsb88LCwtYfP9lVQUFB6+8bGhravW42epRd98vLy2v3cbmOxrybLF++XLNmzWr977q6OuXn53fjGQEAgH3JO++8I0kaO3as4vH4bo/vOuN816vne5KNHmXX/crKyiIdm8tozLtBVVWVpkyZos2bN6t3796SpNraWt1yyy3dfGYAAGBfsfOKufXz5dKnG/OtW7e2a81s9Si77rdzBCP2jMa8i9XV1emkk07SypUrVVxcrLlz5+q0006TJN17772tY48AAABcqqqqtGHDBkn2z5dL+tS/9Pn+++/vcc1s9ig798vPz9fw4cPbfVyuozHvQi0tLZo2bZrmz5+vRCKhOXPm6LDDDtP111+vIAjU2NjYeic0AACAy55u/JSkCRMmtP6c+bx589pcL9s9ys79DjvsMH7GPAIa8y40Y8aM1n+a9v7779fkyZMlffK/oE4//XRJ0kMPPaQVK1Z02zkCAAD/7fz58iAING7cOPM5yWRSEydOlCS99dZbba6XzR6lqalJixYtkiSdcMIJ7Xg12InGvItcf/31+tWvfiVJmjVrli655JLdHg+CQC0tLa3/6hYAAIBl5xXzYcOGqbi42Pm8U089VdInjXltba35nGz3KK+88opSqZQktTb1aB8a8y7wq1/9SjfccIMk6ZJLLtH111+/23PGjh2rM888U5L0+OOPt/5LXgAAAP9qZ2Pu+jGWnS688ELl5+ersbFRTz311G6Pd0aP8thjj0mSDjnkkD2eHz6NxryT/fGPf9SMGTMkSZMnT9b999/vfO6sWbMUi8WUyWR03XXXddUpAgCAfUhjY6OWLVsmac+Nee/evTV16lRJ/9sw79QZPUpjY6N+97vfSZJmzpy5x9eCTwvCMAy7+yQAAADQOd58800dddRRisfjWrFihQYPHtxpez366KO64IIL1Lt3b61evbrNH7PB7rhiDgAAsB+bOHGipk6dqnQ63an/Zkomk9HNN98sSbr66qtpyjuAK+YAAAD7uWXLlmnMmDGKxWJasWKFBg0alPU9nnjiCZ1zzjmqrKzU0qVLVVhYmPU99neJ7j4BAAAAdK5Ro0a1jjv88MMPO6UxT6fTmjVrlr70pS/RlHcQV8wBAAAAD/Az5gAAAIAHaMwBAAAAD9CYAwAAAB5o982fx8emdeZ5AJ1mbmZ2d58CACALpoz+QXefAtAhzy+9tV3P44o5AAAA4AEacwAAAMADNOYAAACAB2jMAQAAAA/QmAMAAAAeoDEHAAAAPEBjDgAAAHiAxhwAAADwAI05AAAA4AEacwAAAMADNOYAAACABxLdfQIAAABoQxC4H4s5HovZ117DhCNPulvChopC+5i4/fx4U2jmhWu3O/eoHt/L3ruPfb7FGzJmnlebdu4RS9vnpYydJ7c02us0tzj3UKqNx9qBK+YAAACAB2jMAQAAAA/QmAMAAAAeoDEHAAAAPEBjDgAAAHiAxhwAAADwAOMSAQAA9sQ1stA1rlBSpiDpyO32q6Ukz8wbe7nbtcAxAbC52D6vWMp+fjzlWEhSosEeTdhjzQ77gNBeq6VngXOP5A57j16/X2zmmbo6M4+PGu7cIyy0/zxcIyRjOxrsvYvt8ZGSFGNcIgAAALDvozEHAAAAPEBjDgAAAHiAxhwAAADwAI05AAAA4AGmsgAAAOwUt69Zhvn2xJR0kZ1LUs3wIjOPtdhTS5K19mSSgi2OUSqSklsa7T2a7ekgoWO6TGxbrXOP6s8faOau6TI1w+zXXbqi3rlH7QH2WsUF+Wa+6ZJDzbzko7RzD5etI+Nm3m+BPUUmcPz5SVJ+nf3n0V5cMQcAAAA8QGMOAAAAeIDGHAAAAPAAjTkAAADgARpzAAAAwANMZekEQcJ+W+N9+2Rtj2VXDTHzdJF9R/fgYVXOtYpm2ndob7wjaeZvT3jCzDen65x7TJz9PTMf/t03nMcAANDlYo6pLAk7z1u/1blUUYk9saXw9fftAw6oMOOth/V27hFP2d/3tw7uaea93tth5tVj7ckrkpTqYfcJ719kTy0pXmW/V9uHFjv3KF1pv44P7rHP68SRb5n5/ywd69yj5PVCM698drOZN1fY55vY4Z6Ss7e4Yg4AAAB4gMYcAAAA8ACNOQAAAOABGnMAAADAAzTmAAAAgAdybipL/OARZh7m23dOr/9CmXOthqPsKSTlpXb+6qH2NJOu8Fx9ifOxn9wzxczfHPuYma9KNZj5rZuOd+4x8NWwjbMDAMAPYWBPIGkuyzfzTIU9SUWSmkviZp53yFAzrxleZObbRju30MZj7f5FSXtySHh2o5l/vMl+fZJ0wMBqMz8wsL+3n//FN838lldOcu7ROCJt5oPKt5v5kAJ7ksqRQ1Y791j6N/uNbCmzp7UkN9aaeabYnkaTDVwxBwAAADxAYw4AAAB4gMYcAAAA8ACNOQAAAOABGnMAAADAAzTmAAAAgAf223GJ6S8ebuZ3PPxzMx+Zl+zM0+kyqdAeN/Sju7/mPCZRZ487Onr25WZesq7FzPM322MUJalovj06CQAAnwTN9pjBwuVVZh4WuUfn1X6+j5l/fIU9kvErI+eb+cod9jqSdErFO2beO77DzN9rHGjmY0Z+5NxjYGKrmX+ctkcxn1Rkj2S8NeW+HjzwBfux2gEDzPzBHnY+4O/1zj36ZOzHmnvZPWAm3x53mdhhf41kA1fMAQAAAA/QmAMAAAAeoDEHAAAAPEBjDgAAAHiAxhwAAADwwH47lSV/2Xoz/0fjgWY+Mm9TZ55Om7634SjnY647sR8eNsfMazL2hJWKu16LfmIR2TsDALDvCIvyzbxuZG8z7/H+FudaLYX29JXTRtuTVHol7Kkh0yvfcO4xa9WpZn7v8MfN/NonzzPz2Ch7iosktawqts9r8itmfvkLk8y87D17yokk9Vhn71+0wX4PE1vq7IW2bXfu0TC+0syTW5vtPbbZ0+bCWOdd1+aKOQAAAOABGnMAAADAAzTmAAAAgAdozAEAAAAP0JgDAAAAHgjCMGzXMI3jY9M6+1y6RPXFR5v59in23b3xRfadyJL0zsy7I+194+ZxZj7vC/bkFUlKb6sx8/DoQ8189bftdYaea98BngvmZmZ39ykAALJgyugfdPoeYTLPzGsOKTPz0qX292lJauzfw8zr+tt71J9qTxTJ+0upc4/p3/qTmf9i4bFmXlhkTyDpf7c9jUaSUj3sIX6ZpD0xJXBMiGvq6Z7KUrbUnsrSUpI089DeWi093HuEcfugwg2NZh6vt9+roLnFuYfSGTN+fumt7mN2wRVzAAAAwAM05gAAAIAHaMwBAAAAD9CYAwAAAB6gMQcAAAA8YN9mux8r//XrZt732d5mnt5S7VzrkDGXmPm7xz5k5s/88gtm3m/ba849XILX7SkrQ+2XBwAA2iFIp828cHPKzFt6FjjXirXY00kKq+09et5lT0ZZMd3eW5J+d9PxZj7/tjvM/OiHrjLzhr72NBFJKtrQZOYbjik088pfvmfmG68e7dyj17v2/s097Va1oNqemBLfYr+3khRvtKepxBzTV1wTVtS+gYYdwhVzAAAAwAM05gAAAIAHaMwBAAAAD9CYAwAAAB6gMQcAAAA8QGMOAAAAeCDnxiW6pDdviXxMansy0vMPOe+fZv7xfXH3QRn32B8AAJBljhF5edWNZh5rtkfwSVK6xB5/GG+0v7cnqrab+cE/c7drG28NzPyIZ/7NzI853u5Ftj7Qw7lH48j+Zl660n4da2YcbOZ9F7h7miDlGFO5yX7fE1vq7IU6cZRhV+CKOQAAAOABGnMAAADAAzTmAAAAgAdozAEAAAAP0JgDAAAAHmAqy144+Jr3zfzisV82818P/ouZf2HaZc49Sp54I/qJAQCArIo1pcw8TLRxjTNjTwhpGGBPa8kc0M/MCzbbe0tS7Fl7Qtw3v/1nM//Nw8fbe99V49yjZE6emffY0Gzm6Xx7UkxevT3xRpJ2DCs188KN9lSWMM9uYYOUe0rOvjCxhSvmAAAAgAdozAEAAAAP0JgDAAAAHqAxBwAAADxAYw4AAAB4gKkseyG9zb6Decu3DjbzD59pMPMf3Phfzj3+37NON/NwgX338oE3vW4vtA/ciQwAgLda0mYcOHJJijfZE0KKHM9v6ltg5nnbm5x79P3Ne2b+SD97+srdM+8380v/51LnHkdfudjMX3xlnJkPesl+3clqe4qLJDX2syfVNPWx84K0PeEl1mRPhJGkoMG9vy+4Yg4AAAB4gMYcAAAA8ACNOQAAAOABGnMAAADAAzTmAAAAgAeYytIJMu/Yd0ifc8PVZv6bWbc711p4lGNiy1F2fEiPy818xAMbnHu0rFztfAwAAHSQYyJarD5l5sltcTOvHlPi3KLPmh5mPuivdWb+zbJvmPnDZ9zn3GPG2+eb+RNn3mnmFx50sZlX3GdPnZGkHqu2m3njQPu1N5c7JtjU2u+tJCWa7Wkxckx46Q5cMQcAAAA8QGMOAAAAeIDGHAAAAPAAjTkAAADgARpzAAAAwAM05gAAAIAHgjB0zPL5F8fHpnX2ueSscNJ452M9b/3IzP/PQX+KtMfoF7/ufGzUDTVmnl6+MtIevpqbmd3dpwAAyIIpo3/Q3afQqcLCpJk39y5yHtPY25583fN9e/zgjoN62vkAe1SjJB198dtm/tdVI8x8zAB7RHMscLecVTcdZOZBxj4mWdNs5vUDCp17lLy/zX4g5RijmEXPL721Xc/jijkAAADgARpzAAAAwAM05gAAAIAHaMwBAAAAD9CYAwAAAB6wb+VFlwr+vtD5WP2Z/cz8s2dfYeZvXnOnmS897kHnHucNOcHMaz7nPAQAAGRZ0GxPB4k3uKeG5G+186qjysy8qCpj5skd7okpa04rN/Pmm/PMfPH6gWb+vXF/du5x6zmDzXzYQ/Z5xbc3mnlquHuCTabAbntjXTCVpb24Yg4AAAB4gMYcAAAA8ACNOQAAAOABGnMAAADAAzTmAAAAgAeYyuK59KYqM6+4y84bv2/fWVwUJJ17PDDkf8z8q6d/x17rqTedawEAgP8rCMw4zLPbr0yRPeWkudz9Pbyh3F6reH3azBN1dl7y/g7nHo2jB5h56Zv5Zj7s3PfNPBPa74ckTRi2xsxXjBpp5n0ck2oKttqvT5JiO5qcj/mCK+YAAACAB2jMAQAAAA/QmAMAAAAeoDEHAAAAPEBjDgAAAHiAqSweyHxuvPOxFdMKzHzM+NVm3tb0FZe7qw+z13p6fuS1AADYLzkmrCjuvsaZKbSnrLSU2NNMGvrZ38NjLaFzj97zN9t7lBWZearU3mPVGb2de5SuyJh585dr7Dxjt5czytY593hqg92LZJL2+95YUWjmacfzJSkssP88goZm5zFdjSvmAAAAgAdozAEAAAAP0JgDAAAAHqAxBwAAADxAYw4AAAB4gMYcAAAA8ADjEjtBMGGMmb//bXtE0QOTHnGudWxBdkb4NIUp52NvVA+1H8hsyMreAADsK1wj9dLF9ojDHQfaY/skKdXDHt2XV2ePPyzaZH/Pz9tc79yjbngvM68ebbd49QfYow8HjN7o3GPqWQvN/G/Vw8y8ssdWM/99XbFzjzV/P9DMe1fZ51v89kdmHpaVOPcIXSMvPcIVcwAAAMADNOYAAACAB2jMAQAAAA/QmAMAAAAeoDEHAAAAPMBUlj1IDB3sfGzFxQPN/PqzHzfzM4o3Z+Wc2nLtpglm/vKdRzmP6fXI6511OgAAdL64fZ0xTLrbnFQve5rKjgPsCWqxtL1O0Ub39LS8rQ32WtW1Zl5z5AFmvmVMuXOP2sH2hJeB4+zJav8+5CUzbwztaTSSdFCyysxjgT0x5c63vmzmz9Qe4dxjxP/sMPNNE+0pK/mH2D1Y2MbglcIPa9wPeoIr5gAAAIAHaMwBAAAAD9CYAwAAAB6gMQcAAAA8QGMOAAAAeCDnprIkhlSaec0RA8z87B8/71xrRtnvsnJObfneBnuayuv32tNXyh9+y8x7ZZi8AgDYP4XxuJkHLfbUEEkK0vY0k7x6Oy95camZ1xw/2rlHLGXvv/64XmbeUuRY57PbnHv826hXzbwsXm/ma1P2hJdM6L5We+M/TjLz+KoCMy9fa6/Td/525x7xLfakmvJl+WZe8JG9VkupPW1HkpR2fz34givmAAAAgAdozAEAAAAP0JgDAAAAHqAxBwAAADxAYw4AAAB4YJ+eypIY0N/5WPVDPcz8W0NfNvNzSzZl5Zzacvm6z5n52/eNdx7TZ84SMy+vZcoKAACSlCnKM/OWkqT7mDz72mRhVZOZp8YMNfP6CnsijCStP94+r7yeO8z8ynEvmnnvuP18SSqIpcy83HHMLe9OMfNwfqlzj2J7YIrK37Pfq+SWRjMP0mnnHlsmDTRz15ScvGJ7Wku8zj6nfQVXzAEAAAAP0JgDAAAAHqAxBwAAADxAYw4AAAB4gMYcAAAA8ACNOQAAAOABr8YlNk+eYOf/Vm3m1w7/o3OtEwrrsnJObdmUbjDzY5/5npmPvm6pmZdvc48+zEQ/LQAAckssMOP8je4xgy2lBWa++dAiM996WIuZX/f53zn3qE7bo5vP7/mOma9pKTTz29dNdu7hsvilEWaerLXfq/7z7BGHkpTcZL+PmSJ7HOX2kSVmHrZxObjoY/v9za+qN/Og2X6+WtwjGfcFXDEHAAAAPEBjDgAAAHiAxhwAAADwAI05AAAA4AEacwAAAMADXk1lWX2a/feE98fOztoeP982zMzvfPkEMw/S9t3LkjT6xlVmPmLTm2a+b98nDACAn4KU/R22pcyeciJJia32tI9Mnj1JxeXJ9fZEOUnKT9iTQx794Egzb17Qy8wT9qlKkpI1oZn33WLPdSvc1GTm8bqUc4/qw8rNPOZobAo222vlb7an2Un775SVqLhiDgAAAHiAxhwAAADwAI05AAAA4AEacwAAAMADNOYAAACAB7yayjLyW2+Z+Ve/dUTn7y1777bk1n3CAAD4KVVWYD9gDyyRJH10Rh8zL11hTzPpMde+llm1qNK5R9nyZjPv7Tivugr7gdKV7mkmjX2TZl78txVmXnvscDMPg3znHj022lNW8qobzTzW5JjwkmMTVjqCK+YAAACAB2jMAQAAAA/QmAMAAAAeoDEHAAAAPEBjDgAAAHjAq6ksAAAAUTmng+xwTzOp/GOh/YA9lEXxbTvMvCzV4twjbLTPK+hRZOZ528rs57cxXSZvuz3pZOsJI+zn19svsGi9+72K1dtTVoJmx/QVdBhXzAEAAAAP0JgDAAAAHqAxBwAAADxAYw4AAAB4gMYcAAAA8ACNOQAAAOABxiUCAIB9WpCyRwZmejpGIkpq7m0/Fgb289OVPcy8qTTu3CN0XP5MFdubxBvsuYjJOve8xESDPf6w5wd1Zs7oQ79xxRwAAADwAI05AAAA4AEacwAAAMADNOYAAACAB2jMAQAAAA8EYRi6b/UFAAAA0CW4Yg4AAAB4gMYcAAAA8ACNOQAAAOABGnMAAADAAzTmANAFnn32WZ188skaOHCgksmkgiBQEARasmSJmpubNWLECAVBoDlz5nT3qTpddtllCoJAF110UXefCpDzOlI32qpDHV0ziig1JGfrTQgA6FTXXXddKGm3X8lkMkylUuFtt90WSgrHjBkTZjKZ7j5dpzVr1oTJZDIMgiCcP39+d58OkNOi1o091aGOrBlVlBqSq/WGK+YA0IkWLlyom266SZI0adIkPf3001q4cKEWL16sRYsWqaGhQT/5yU8kSdddd52CIOjO021TZWWlLrroIoVhqB/+8IfdfTpAzqqtrY1UN/ZUhxKJROQ1OyJKDcnZetPNfzEAgP3ajBkzQklhz549w+rq6t0e/8lPfhJKCisrK8N0Ot0NZxjN0qVLW6+05dJVLMAnUevGnupQR9bsqCg1JBfrDVfMAaATvfDCC5KkKVOmqFevXp96LJ1O65577pEknXvuuYrF/C/Jo0aN0uGHHy5Juvvuu7v5bIDc05G60VYd6uiaHRWlhuRivfH/uwAA7KM2b96slStXSvrkfx//q7lz52rt2rWSpPPOO69Lz21v7DzX2bNnq7a2tpvPBsgtUevGnupQR9bcW1FqSK7VGxrzTpDJZPToo4/qhBNOUN++fdWjRw+NHz9ed911l9LptHbs2NF6J/Ttt9/eelxdXZ0qKioUBIEOOuggpVIpc/2GhgYdc8wxCoJABQUFeuWVV7rqpQFoh/PPP19BEKhv376t2ZVXXtn6uQ+CQI899piefPJJSdKIESM0duzYdq3d0foiZa/GnHHGGZKk+vp6Pf300+1+XwB8Ym8+x+2tG+2tQ1HW7I4aknP1prt/lmZ/s2nTpvDoo48273yWFH71q18N582b1/rfc+fO/dTxP/vZz1ofe+CBB3ZbP51Oh6effnooKYzFYuHs2bO76qUBaKfDDz/cWQN2/lqwYEE4ZMiQUFJ4wQUXtGvdva0vYZi9GtO/f/9QUjh9+vRobw6Q4/b2c9zeutHeOhRlzTDsnhqSS/WGxjyLamtrw4MPPjiUFAZBEE6fPj185plnwn/84x/h7NmzWz8kxx13XOsX9ccff/ypNRobG8PKyspQUjhkyJCwubn5U49fccUVrcf+7Gc/68qXB6Cdli5dGi5evDicOnVqKCmsqKgIFy9e/Klfa9eubf0s33XXXXtcMxv1JQyzV2NOOeWUUFJ4wAEHRH+DgBy1t5/jKHWjPXUolUpFrkXdUUNyqd7QmGfReeedF0oKE4lE+PTTT+/2+I4dO8JBgwa1fsEOHDjQXOfBBx9sfc7999/fmt9+++2t+VVXXdVprwNAdkycODGUFE6ePHm3x5544onWz/Orr766x7WyVV/CMDs15oYbbmh97saNG/d4/gD2/nMctW6EYdt1qKNrdnUNyaV6Q2OeJa+88krrF82NN97ofN5//Md/tD7vxBNPNJ/T0tISjhw5MpQUDh48OGxubg6feOKJMAiCUFJ47rnnev2PkAD45H/nFhUVhZLCa665ZrfHf/rTn7bWgmXLlrW5VjbrSxhmp8bcd999rXu99tprbT4XQHY+x1HqRhjuuQ51ZM0w7Poakkv1hps/s+THP/6xpE8G4l9zzTXO5x1yyCGtvz/00EPN58Tj8db11qxZo5kzZ+rCCy9UGIY67rjj9PDDD3v9j5AAkJYvX676+npJ9mf9448/bv29Nb5sV9msL1J2akx5eXnr7zdu3NjmcwFk53McpW5Ie65DHVlT6voakkv1hsY8C9avX68///nPkqSZM2cqkUg4n1taWtr6+/Hjxzufd9ZZZ7U+/uCDD6qpqUljx47VU089pWQymZXzBtB5Fi5c2Pp767NeXV3d+vu2vhl2Rn2R9r7G7HrOdXV1e3w+kMuy9Tlub93YaU91qCNr7tSVNSSX6g2NeRY8//zzrb8/8cQT23zurh+Atq5oBUGgSy+9tPW/+/fvr+eee+5TH1gA/tr5DbGwsFAjR47c7fGCgoLW3zc0NDjX6Yz6Iu19jdn1nPPy8tp1DJCrsvU5bm/d2GlPdagja+7UlTUkl+oNjXkW7PzCz8/P3+Ms4iVLlkiSioqKNGLECOfzli9frlmzZrX+d11dnfLz8/f+ZAF0iXfeeUeSNHbsWMXj8d0e33W28K7fiP9VZ9QXae9rzK7nXFZW1u7jgFyUrc9xe+vGTnuqQx1Zc6eurCG5VG9ozLNg58877frF7TJ37lxJn3xIXP/kbVVVlaZMmaLNmzerd+/ekqTa2lrdcsstWTpjAJ1t5zdi15XrXevF1q1bnetku75I2akxu55zZWVlu48DclG2PsftrRs77akOdWRNqetrSC7VGxrzLGhqapIkbd++vc3nLV26VC+++KIk98961dXV6aSTTtLKlStVXFysuXPn6rTTTpMk3XvvvVq3bl3WzhtA56iqqtKGDRskuT/ru141e//9951rZbO+SNmrMTvPOT8/X8OHD2/XMUCuytbnuL11Q2pfHYq6ptQ9NSSX6g2NeRb069dP0icfuLVr15rPyWQyuuKKKxSGoST7b68tLS2aNm2a5s+fr0QioTlz5uiwww7T9ddfryAI1NjY2HoXNAB/teeGqwkTJrT+bOe8efOca2WrvkjZrTE7z/mwww7b73/mE9hb2foct7duSO2rQ1HX7K4akkv1hsY8C4466qjW31v/GyedTuvyyy9vvSNbsj8kM2bM0HPPPSdJuv/++zV58mRJn3w4Tz/9dEnSQw89pBUrVmTz9AFk2c6f6wyCQOPGjTOfk0wmNXHiREnSW2+95VwrW/VFyl6NaWpq0qJFiyRJJ5xwQpvPBZC9z3F764bUvjoUdc3uqCE5V2+6c4j6/qKmpibs1atX6/D7r3/96+Ff//rX8K233goffvjhcMKECaGk1n/CNgiCsLa29lNrzJo1q/X4WbNm7bbHokWLWgf3T58+vYteGYCOmD59eigpHD58eJvPu+OOO0JJYUFBQbh9+3bzOdmoL2GY3RrzwgsvtK61YMGCNp8LIHuf4zBsX90Iw/bXofau2V01JNfqDY15lvz2t78NE4lE6xfPrr8SiUT4wx/+MLz88stDSeGYMWM+deyu/7TtJZdc4txj2rRpoaQwFouFixYt6uyXBKCDPvOZz4SSwjPPPLPN523evDnMz88PJYWPPPKI83l7U1/CMPs15mtf+1ooKTzkkEPafH0A/tfefo53am/daG8das+a3VlDcq3e0Jhn0WuvvRaeeOKJYVlZWVhQUBAOGzYs/OY3v9n6xTlu3LhQUjhz5szWY/7whz+0flAnT54cplIp5/pLliwJY7FYKCk85ZRTOv31AIiuoaEhjMfje/xnt3c699xzWz//belIfQnD7NeYhoaGsGfPnqGk8Oc///keXx+A/9XRz/G/2lPdiFqH2lqzO2tILtYbGvMusmDBgta/bf7lL3/p7tMB4Ik33ngjlBTG4/Fw9erVHVqjK+vLf//3f4eSwt69ezv/VzuA6KJ8jrNRN7piTUuUGpKL9YabP7tAJpPRd7/7XUnSqFGjdNxxx3XzGQHwxcSJEzV16lSl0+kO/VsFXVlfMpmMbr75ZknS1VdfreLi4k7bC8glUT/He1s3umrNfxWlhuRsvenuvxnsDz744APnYw0NDeGFF17Y+rfgZ555pgvPDMC+YOnSpWEikQiTyWS4du3aTz3mU315/PHHW29Qq6+v79S9gP1JZ3yO26obHdUZa+4qSg3J1XqT6La/EexHpk2bpvz8fJ1zzjkaP368SktLtXXrVr3xxhv65S9/qdWrV0uSvvvd7+rkk0/u3pMF4J1Ro0a1jhj78MMPNWjQoNbHfKov6XRas2bN0pe+9CUVFhZ26l7A/qQzPsdt1Y2O6ow1dxWlhuRqvQnC8P9OskeHtLS0qLi4uPVf9bIkEgldf/31uvbaaxUEQReeHYB9GfUF2PfxOUYUNOZ7qaWlRXPmzNGzzz6rt99+Wx9//LFqamrUs2dPDR8+XF/+8pc1Y8YMVVZWdvepAtjHUF+AfR+fY0RBYw4AAAB4gKksAAAAgAfaffPn8bFpnXkeQKeam5nd3acASVNG/6C7TwHosOeX3trdp5DzqCHYV7W3fnDFHAAAAPAAjTkAAADgARpzAAAAwAM05gAAAIAHaMwBAAAAD7R7KgsA7K0wmWfneXHnManyAjPP29YYba2M+59siNWnzDxIp+0D0hnnWgA6T9Qa4qofUvZqiKt+SNQQRMcVcwAAAMADNOYAAACAB2jMAQAAAA/QmAMAAAAeoDEHAAAAPMBUFgBdJsy3px3Etje4D+plT1XIFNjTGRLb7LWa+hc7t8ivazbzMGmXyKDBfj6AzhW5hjjqh5S9GuKqHxI1BNFxxRwAAADwAI05AAAA4AEacwAAAMADNOYAAACAB2jMAQAAAA/QmAMAAAAeYFwigC7T3LvQzJNtHLNhUr6ZD/ntZjNvKbP3yOS5r0OkyovMPIwHZp5k1BnQLaLWEFf9kLJXQ1z1Q6KGIDqumAMAAAAeoDEHAAAAPEBjDgAAAHiAxhwAAADwAI05AAAA4IFumcry/n1HmvmQ32fMPPmn+Z15OgC6SGgPKJDsj74kqXSF/eD6/6ePmcdSoZn3WdTg3KOhnz25IVVkn3Cv6jznWkFjyvnYPiVuX7cJ43H7+Qn3dZ6gvsleK2m/j0Fo/xlKkjKOL5Z0G19E2G9ErSGu+iFlr4a46ocUvYbkbP2QnDUk1+oHV8wBAAAAD9CYAwAAAB6gMQcAAAA8QGMOAAAAeIDGHAAAAPBA90xlOeU+M//i0GlmnvxTZ54NgK6SLrCvBcS37XAeE0sVm3nNmGjTC3ovdj9W8uJSM9825WAzTxe7pzAkmlrsB9qaFNDZEu5JCJl8e7JBqrzAzBv62M/Pq3NPNSheXGvmjZVlZp7cZk9hkKSgxd4naGi2D+jO9x1ZF7WGuOqHlL0a4qofUvQa4mX9kJw1JFv1Q3LXkFyrH1wxBwAAADxAYw4AAAB4gMYcAAAA8ACNOQAAAOABGnMAAADAA90ylcXl2uF/NPN7B3zReUzLho2ddDYAsq2p1L6zvyTlmEQgqb7CPmbksHVmftbA+WZ+z6Kpzj2KBoy28432nfo7Dix0rlW2rcF+oCXtPKazhQn3NZhYs/3ex5rt6QWlS2vM3DnVQFJYZE9oyOQ7Jmxs3u5cKzWwl5knmhwTNtJMZdmfRK0hrvohZa+GuOqHFL2G+Fg/JHcNyVb9kNw1JNfqB1fMAQAAAA/QmAMAAAAeoDEHAAAAPEBjDgAAAHiAxhwAAADwAI05AAAA4AGvxiWeUFhn5jc+1MN5TM+vdNbZ5I7EgP7Ox5pGDTTz+Etvd9LZYH8WOi4FhI2NzmPKltsjtPIT9piu6rRdL7aNc4zDktT3Efux+I4mM08dlO9cK1OYZ+axHfb4MIVZHMcVBGacKbDPSZIUcxyTtP+wGvvb72/1aHsMmSQlGuzXmL/dzoNR/Zxr1VfYr6XXYsfYtjr31xb2PVFriKt+SNmrIa76IUWvIZHrh5S9GuKoH1IbNSRL9UNy15Bcqx9cMQcAAAA8QGMOAAAAeIDGHAAAAPAAjTkAAADgARpzAAAAwAPdMpXl59uGmfllZSvM/FtDX3au9cSQY8y8ZfWH0U9sP1d98dFmfsU1s53HHFGw1syv/sqFZp5+b3n0E0POSBXbd/AHPYqcxwSOgQNrttp38N8/9B0zf6DnJOcemz5bbOaDfrvZzPPqSp1rtZTY0xbymtJmHjS7Jzq4hEl7skCmyM7rD3C/v4VV9tSIlqK4mTeU2/kJF73u3OP3Sw8186F323+4iXdXOddq/OpnzDxTYH87i9nDvrCPilpDXPVDyl4NcdUPKXoNiVo/pOg1JGr9kNw1JFv1Q3LXkFyrH1wxBwAAADxAYw4AAAB4gMYcAAAA8ACNOQAAAOABGnMAAADAA90yleXOl08w88tOvc/Mzy3Z5Fzrl0cMMPMeOTyVJd6nt5lvn2LfXtzW+yslzTTMd9+9DbjEG+y76FsqypzH1FXYX2vNC+yJCmvGFZr5leNedO5x76KTzbzmyAPMvGhTs3Othn72ZyZeb09OiKftaQth0l2e08X25IbmMse0hTx7koUkJaq2O44pM/MtF9t1pFei3rnHV0a+a+ZvDp9g5r23VTjXStZmzLylxH7tiS3OpbAPilpDXPVDyl4NcdUPKXoNiVo/pOg1JGr9kNw1JFv1Q3LXkFyrH1wxBwAAADxAYw4AAAB4gMYcAAAA8ACNOQAAAOABGnMAAADAA90ylSVIuycERHX2j5838+ffHm/mLavWZG1vX6W3VJt5fNEo+4DPRd9j/RfKzLz/wuhrIXck6+yJCoEdS5JKVzaYeX3/IjO/fd1kMz+z3z+ce8Q+u83M67aXmXnJshrnWk3l9t39TX0LzDwZt6+PpAvd5bm20p7cUL6k1szjPe0pDJKkPHufFdPjZl72l1Iz/9L4fzq3uPHDr5r5ttH282Npe7KUJBVW2dMpGnvb73vBGsf3m7CNLzp4K2oNcdUPKXs1xFU/pOg1JGr9kKLXkKj1Q2qjhmSpfkjuGpJr9YMr5gAAAIAHaMwBAAAAD9CYAwAAAB6gMQcAAAA8QGMOAAAAeIDGHAAAAPBAt4xLHH3jKjPfdKo91qgiXuhca0bZSjP/+cUnmvngH+1b4xLjZe7RQu/9ZKSZ5/VsNvN3j707K+ckSQ1H1dkP3Jm1LbAfSjRkzDxVYo/vkqRUT3vsVrIm2riqgpg9JkuSZox61cz/8wN7TFfPoT2da5W8Z48r3TKhj5nHmuz3pKmXuzz3+cMHZh7W2Z/LmvMOda6V+Y49Dqz//7Hf96n//iczn7XqVOceZwx828xv6X2gmRevc/9Zxbfb9c01Zk4xx7izNOMS90VRa4irfkjZqyGu+iFFryFR64cUvYZErR+Su4Zkq35I7hqSa/WDK+YAAACAB2jMAQAAAA/QmAMAAAAeoDEHAAAAPEBjDgAAAHigW6aypDdVmfmxz3zPzJedfm/kPa4/+3Ezf/iZk5zHhPOXRN4nW2KHHmzm5fdtcB7z/uBfdNbp7FF5qX33dpBwf0mFLS2ddTrYRxRusCcvbR/ew3lMrxeWm3nmC8PNfPFLI8y8/Lw/O/eoy+Sb+cBxG828umqAc638zUVmXrAtbeY1w+xpEv3e2Obco2X4QDP/6Ev2+3jROXOda93/ypfM/B+33WHmE16+zMxf+Lx76tPCJvt8lbSnSWwfXOBcq/c/6s28udgxPSHmuP6UtveG36LWEFf9kLJXQ1z1Q4peQ6LWDyl6DYlaPyR3DclW/ZDcNSTX6gdXzAEAAAAP0JgDAAAAHqAxBwAAADxAYw4AAAB4gMYcAAAA8EC3TGVxGX3dUjO//MjPOY+554C/mfkZxZvN/Npv23cvS9KIC9s4uSypvuRoM3981n+YeWWi0LnWxWu+bObv/uYzZn7KN1428+v6LHLu4fLqoU+Y+Wl93VNvWjbYd6cjd8TqU2aeV+++w732WHtyQuGmJjNP1tpTDWYuPM+5x0Uj37CPGfKSmf/7AWc710qV2jUmUWdPVSiqsqcB7Diop3OPDcfY11Tun3q/mf/bvd90rnXMGf8086MfusrMC8fWmPmJ/321c48LTnnRzPtW2Gv1es89VSEM7PcrZn9pKUzY71XgeD78FrWGuOqHlL0a4qofUvQaErV+SNFrSNT6IblrSLbqh+SuIblWP7hiDgAAAHiAxhwAAADwAI05AAAA4AEacwAAAMADNOYAAACAB7yaypLeZt9h+/Z99iQTSdKN9lQWlwcmPeJ87NbPXWDmsb8tjLRHOGm887HfzLrdzN9t7mfmF3zfPSqmZPY8M++Xec3M5/2mj5l/789HOff46QD33eZAVEGzfSt74foG5zF1B9oTEuJ19loV8xrNfEOi1LlHZoR9jaIxzDPzAaOrnGttrO5v5kN/u8XMdwwqN/PmYnt6gCQ9dMa9Zv61Z2eYedEk9ySErWf1MPPex9hTLgpfyjfzdTN3OPcYU/iRmf9qkz1xKz7c/a2p/NW19jGpMjMPk/ZaQUOzcw/4K2oNcdUPKXs1xFU/pOg1JGr9kKLXkKj1Q3LXkGzVD8ldQ3KtfnDFHAAAAPAAjTkAAADgARpzAAAAwAM05gAAAIAHaMwBAAAAD9CYAwAAAB7walyiS585S5yPfe8ye9Sfa8zfsQXuETeXTisw8xGOiYzxCnvEYdGt9mgfSRqasPc4ff7pZj7oieyNK3SNo1y5Y0DW9gA6IlZvjy2TpPxtLWa+7TMlZl6+oNrM82rtz54kPbDIHrv10NG/NvOpgxY613pyxQlm3tS/2Mz7/GWNmQ/+vf06JGnG2+eb+ZePWmzmS+4c61yrcaT9vhRtaDLz5lJ7/FvLKvcYtIETtpr5AQPt19jYo8K5VvXnDzTzvHp7PFtDRaGZ96ipd+6BfY+rhrjqh5S9GuKqH1L0GhK1fkjRa0jU+iG5a0i26ofkriG5Vj+4Yg4AAAB4gMYcAAAA8ACNOQAAAOABGnMAAADAAzTmAAAAgAf2iaksmdpa52Ov33u0/cD/F32ayZjxq83cvrdYWn/2cDOfd9DdkfcOF5RGPibyHkcfauYPD7u3jaPckywsy64a4nxs2Pc2RloLuSNodk9lyatuNPNYX/vu/kxR0szL33N9kqWGCvvr/O3xQ8z8b9XDnGvVTrXrVdNf7M94zYX9zXzFqt7OPX7z2V+Z+dlzrjTzIY4JCZLU0Md+H2sOipt5r+VpM58++RXnHh+n7ekXsSA08y1H2HtI0shHHF8Pjfb0jbrB7mkW2H+4aoirfkjZqyGu+iFFryFR64cUvYZErR+Su4Zkq35I7hqSa/WDK+YAAACAB2jMAQAAAA/QmAMAAAAeoDEHAAAAPEBjDgAAAHhgn5jK0pbyh98y82svn2DmN1fMd6712+F/MPP6j5rNvCCY51jJ/fedQ/7rcjMfetPrzmOyZfW37bwmY9/ZLElFgX13c15g33GdLspEPi+gLbEme9pCwWY73z7Mvou+5/vu6U7Fa+0pDHe+9WUzP3XcO861WkL7s5E8d5OZL14/0MzHDNjg3OPCty8280Ev2ZMFdhxgvz5Jqj3QrleDf/GevdYXRpj5f/19knOPG05718zXVb5p5nf/5TTnWjXDisy815IaM483uesb9n+u+iFlr4a46ocUvYZErR9S9BoStX5I7hqSrfohuWtIrtUPrpgDAAAAHqAxBwAAADxAYw4AAAB4gMYcAAAA8ACNOQAAAOCBfX4qizL21JCX7zzazJvamH6SH+SZeVHgvuPaMvrFr7sfe8C+S7ol7Pw7f4eea98FPlOfcx6z6dvHmPmPr3jYzAcPq4p8XkCbWuzPeP7mBjNv6Gt/joO0vY4k9Z2/3cy3HmJPZ/hiz6XOtX42YIeZ/2LbAWZ+Qp9/mvmftxzs3KPivgIzT9TZE6Qkdw1zTR344OrRZt7L8dLL3g2cexyU+aaZh3n2FKfCNkpu6Yp6M08X2QcVrrX/bJEjHPVDyl4NcdUPKXoNiVo/pOg1JHr9kFw1JFv1Q3LXkFyrH1wxBwAAADxAYw4AAAB4gMYcAAAA8ACNOQAAAOABGnMAAADAAzTmAAAAgAf2/XGJDr0esccinnPx6c5jjipfZeZXlC8w80tXf9XMR91Q49yjZeVq52M+StTZo5B+smKKmZde7h776B5aBUQXNLeYeeEWO986tsy5Vu+/rzfzEY/a1y5+sO0C51r3TVpr5qX59mi2+SsGm3nl43HnHoVVtWZeN7SnmYcx9yjDnsvselWy0t4/jDuu57SxR4+NhWa+/hS7KpQ69pakqsPt8XMl6+w/9+aD7L3L533s3AO5IVs1xFU/pOg1JGr9kKLXkKj1Q3LXkKzVD8lZQ3KtfnDFHAAAAPAAjTkAAADgARpzAAAAwAM05gAAAIAHaMwBAAAAD+y3U1lcUl/c4Hzs1VgPM//Tqd8x86Kn3nSstCXiWe17Xho728wn3jHdeUzfUzrrbJCTWuw78vOr6s081aPEuVTDyH5mXj0qaea9l2Sca1VVVZr55mZ7YtGwZY1m3lLgnnDUONB+LYHjtAqq7D0kKV2cb+YNFQVmntxuTy+IN9i5JNUNsKckDOpfZeYbT3ZPhxh+uT21ImxsMvOarx/qXAs5Lks1xFU/pOg1JGr9kKLXkKj1Q3LXkGzVD8ldQ3KtfnDFHAAAAPAAjTkAAADgARpzAAAAwAM05gAAAIAHaMwBAAAAD+TcVJY2Zew7tN3TV/Z/Jevsu6Q/amkw87cnPOFca7LGZ+OUgDYFzfbXbNE6e9KCJMV32Hfkl8u+u7/w3fXOtRKHDzLzgk32Z8a1d6ynPdVAkprL7EkIrj1i9SnnWkHKfr8K8uzrNmlHnthS59wj1cOeeDV10AIzX92nj3Otl6YfaeYVd79m5sUb2hg1ARii1hDXZ1iKXkOi1o+29nfVkKj1Q3LXkGzVD8ldQ3KtfnDFHAAAAPAAjTkAAADgARpzAAAAwAM05gAAAIAHaMwBAAAADzCVBW3K32zfpX3rpuPN/KYBf+7M0wH2rMWerhRrdE8maSktNPMwsJ8flpU410on7YOa+toTEhpH2hMH8upC5x6FVY4pDE32hATX5ARJUujYJ23nLWVx+/nbtju3GPD3YjN/65Shdr56iHutj+w/3/io4WaeV2s/H3CKWENc9UOKXkOi1g8peg2JWj+kNmpItuqH5KwhuVY/uGIOAAAAeIDGHAAAAPAAjTkAAADgARpzAAAAwAM05gAAAIAHaMwBAAAADzAuEW0K5y8x83/eMtHMJ37+UOdaw/VGVs4J6IigjXGJCcdjiW32tYsw7h75VfKBPfLLtX9Bcb6ZxxrbGHGYydh52pF3QNyxfxi3z7dhfKVzrbxa+7UvfXS0mZfEHTPmJEn2+LKwMGnmMcfYNiAq12fYVT+k6DUkav2QOlBDPKwfkruG5Fr94Io5AAAA4AEacwAAAMADNOYAAACAB2jMAQAAAA/QmAMAAAAeYCoLOqToqTfNfPhTXXwiQGdyTCkIsji9IFbbmLW1silW32zmhRscE2naGITQ3MueeNDnnXozT2xrcK714cl9zLzHWsd1pgxTWdCNcrSGRK4fkrOG5Fr94Io5AAAA4AEacwAAAMADNOYAAACAB2jMAQAAAA/QmAMAAAAeYCoLAGB3jqkRcce0hSCVdi6VybcnMSS22lMVmvuXONfqt8DeP1Zjr5VMM5UF6HIR64fkriG5Vj+4Yg4AAAB4gMYcAAAA8ACNOQAAAOABGnMAAADAAzTmAAAAgAeYygIA2F1oTyMImlvsp8ftyQmSlNiRMvNMcUGk50tSqiTpWKvQzGOO8wXQiSLWD8ldQ3KtfnDFHAAAAPAAjTkAAADgARpzAAAAwAM05gAAAIAHaMwBAAAAD9CYAwAAAB5gXCIAoP3SGTMOHLkkBc3u8WVR5dfUZ20tAF2srTrhqi05Vj+4Yg4AAAB4gMYcAAAA8ACNOQAAAOABGnMAAADAAzTmAAAAgAeCMAzD7j4JAAAAINdxxRwAAADwAI05AAAA4AEacwAAAMADNOYAAACAB2jMAQAAAA/QmAMAAAAeoDEHAAAAPEBjDgAAAHiAxhwAAADwwP8P0WZx6vOWbGsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1000x500 with 5 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "GCNN = GroupEquivariantCNN(in_channels = 1,\n",
    "                           out_channels = 1,\n",
    "                           kernel_size = 3,\n",
    "                           hidden_dim = 4,\n",
    "                           group_order = 4,\n",
    "                           num_gconvs = 3)\n",
    "test_equivariance(GCNN, 4, True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nTcOF87al9i8"
   },
   "source": [
    "## 1.6 Experiments on MNIST data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "torch.manual_seed(0)\n",
    "random.seed(0)\n",
    "np.random.seed(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 211
    },
    "id": "PHRCjKgdl89g",
    "outputId": "d32ffd09-a112-424a-b010-b09449629480",
    "tags": []
   },
   "outputs": [],
   "source": [
    "# We normalize the training data.\n",
    "train_transform = torchvision.transforms.Compose([torchvision.transforms.ToTensor(),\n",
    "                                                  torchvision.transforms.Normalize((0.1307,), (0.3081,))\n",
    "                                                  ])\n",
    "\n",
    "# To demonstrate the generalization capabilities our rotation equivariant layers bring, we apply a random\n",
    "# rotation between 0 and 360 deg to the test set.\n",
    "test_transform = torchvision.transforms.Compose([torchvision.transforms.ToTensor(),\n",
    "                                                 torchvision.transforms.RandomRotation(\n",
    "                                                     [0, 360],\n",
    "                                                     torchvision.transforms.InterpolationMode.BILINEAR,\n",
    "                                                     fill=0),\n",
    "                                                 torchvision.transforms.Normalize((0.1307,), (0.3081,))\n",
    "                                                 ])\n",
    "\n",
    "# We demonstrate our models on the MNIST dataset.\n",
    "train_ds = torchvision.datasets.MNIST(root=\"data/MNIST\", train=True, transform=train_transform, download=False)\n",
    "test_ds = torchvision.datasets.MNIST(root=\"data/MNIST\", train=False, transform=test_transform)\n",
    "train_loader = torch.utils.data.DataLoader(torch.utils.data.Subset(train_ds, list(range(50000))), batch_size=64, shuffle=True)\n",
    "valid_loader = torch.utils.data.DataLoader(torch.utils.data.Subset(train_ds, list(range(50000, 60000))), batch_size=64, shuffle=False)\n",
    "test_loader = torch.utils.data.DataLoader(test_ds, batch_size=64, shuffle=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def train_classifer(model, train_loader, valid_loader, test_loader, learning_rate = 0.001, num_epoch = 100):\n",
    "    print(\"number of paramters:\", sum(p.numel() for p in model.parameters() if p.requires_grad))\n",
    "    for epoch in range(num_epoch):\n",
    "        optimizer = torch.optim.Adam(model.parameters(), learning_rate)\n",
    "        scheduler = torch.optim.lr_scheduler.StepLR(optimizer, step_size=1, gamma=0.95)\n",
    "        loss_fun = nn.CrossEntropyLoss()\n",
    "\n",
    "        train_acc = []\n",
    "        for x, y in train_loader:\n",
    "            x = x.to(device)\n",
    "            y = y.to(device)\n",
    "            pred = model(x)\n",
    "            loss = loss_fun(pred, y)  \n",
    "            acc = (y == pred.argmax(dim=-1)).float().mean().cpu().data.numpy()\n",
    "            train_acc.append(acc)\n",
    "            optimizer.zero_grad()\n",
    "            loss.backward()\n",
    "            optimizer.step()\n",
    "\n",
    "        valid_acc = []\n",
    "        for x, y in valid_loader:\n",
    "            x = x.to(device)\n",
    "            y = y.to(device)\n",
    "            pred = model(x)\n",
    "            acc = (y == pred.argmax(dim=-1)).float().mean().cpu().data.numpy()\n",
    "            valid_acc.append(acc)\n",
    "            \n",
    "        test_acc = []\n",
    "        for x, y in test_loader:\n",
    "            x = x.to(device)\n",
    "            y = y.to(device)\n",
    "            pred = model(x)\n",
    "            acc = (y == pred.argmax(dim=-1)).float().mean().cpu().data.numpy()\n",
    "            test_acc.append(acc)\n",
    "\n",
    "        print(\"Epoch {} | Train Accuracy: {:0.5f} | Valid Accuracy: {:0.5f} | Test Accuracy: {:0.5f}\".format(epoch+1, np.mean(train_acc), np.mean(valid_acc), np.mean(test_acc)))\n",
    "            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's train the conventional CNN first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of paramters: 10842\n",
      "Epoch 1 | Train Accuracy: 0.62604 | Valid Accuracy: 0.86425 | Test Accuracy: 0.28165\n",
      "Epoch 2 | Train Accuracy: 0.88621 | Valid Accuracy: 0.92207 | Test Accuracy: 0.30314\n",
      "Epoch 3 | Train Accuracy: 0.91932 | Valid Accuracy: 0.93999 | Test Accuracy: 0.31559\n",
      "Epoch 4 | Train Accuracy: 0.93354 | Valid Accuracy: 0.94984 | Test Accuracy: 0.32604\n",
      "Epoch 5 | Train Accuracy: 0.94100 | Valid Accuracy: 0.95611 | Test Accuracy: 0.34047\n"
     ]
    }
   ],
   "source": [
    "class CNN(torch.nn.Module):\n",
    "    def __init__(self, in_channels, out_channels, kernel_size, hidden_dim):\n",
    "        super().__init__()\n",
    "\n",
    "        self.model = nn.Sequential(\n",
    "            nn.Conv2d(in_channels, hidden_dim, kernel_size, padding=(kernel_size-1)//2),\n",
    "            nn.ReLU(),\n",
    "            nn.Conv2d(hidden_dim, hidden_dim, kernel_size, padding=(kernel_size-1)//2),\n",
    "            nn.ReLU(),\n",
    "            nn.Conv2d(hidden_dim, out_channels, kernel_size, padding=(kernel_size-1)//2),\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        out = self.model(x)\n",
    "        # Apply average pooling over spatial dimensions.\n",
    "        return out.mean((2,3))\n",
    "    \n",
    "CNN = CNN(in_channels = 1, \n",
    "          out_channels = 10, \n",
    "          kernel_size = 5,\n",
    "          hidden_dim = 16).to(device)\n",
    "train_classifer(CNN, train_loader, valid_loader, test_loader, num_epoch = 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Next, we train the GCNN!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of paramters: 14600\n",
      "Epoch 1 | Train Accuracy: 0.26065 | Valid Accuracy: 0.51055 | Test Accuracy: 0.35460\n",
      "Epoch 2 | Train Accuracy: 0.57225 | Valid Accuracy: 0.63167 | Test Accuracy: 0.47193\n",
      "Epoch 3 | Train Accuracy: 0.66466 | Valid Accuracy: 0.73806 | Test Accuracy: 0.55195\n",
      "Epoch 4 | Train Accuracy: 0.70542 | Valid Accuracy: 0.70233 | Test Accuracy: 0.55693\n",
      "Epoch 5 | Train Accuracy: 0.72456 | Valid Accuracy: 0.76294 | Test Accuracy: 0.51612\n"
     ]
    }
   ],
   "source": [
    "GCNN = GroupEquivariantCNN(in_channels = 1,\n",
    "                           out_channels = 10,\n",
    "                           kernel_size = 5,\n",
    "                           hidden_dim = 8,\n",
    "                           group_order = 4,\n",
    "                           num_gconvs = 3,\n",
    "                           classifer = True).to(device)\n",
    "train_classifer(GCNN, train_loader, valid_loader, test_loader, learning_rate = 0.01, num_epoch = 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation with the exercise codes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "class Image_GConv_Layer(nn.Module):\n",
    "    def __init__(self,\n",
    "                 in_channels,\n",
    "                 out_channels,\n",
    "                 kernel_size,\n",
    "                 group_order,\n",
    "                 activation = True,\n",
    "                 lift_conv = False\n",
    "                 ):\n",
    "        super(Image_GConv_Layer, self).__init__()\n",
    "        \n",
    "        self.lift_conv = lift_conv\n",
    "        self.kernel_size = kernel_size\n",
    "        self.activation = activation\n",
    "        self.group_order = group_order\n",
    "        self.rep_2D = np.array([[[np.cos(theta), -np.sin(theta)],\n",
    "                                 [np.sin(theta), np.cos(theta)]] \n",
    "                                for theta in [np.pi*2/self.group_order*i\n",
    "                                              for i in range(group_order)]])\n",
    "\n",
    "        # Compute the multiplication tables given group elements\n",
    "        self.mul_table = group_conv.make_multiplication_table(self.rep_2D)\n",
    "        \n",
    "        # Compute regular representations of all group elements\n",
    "        self.rep_reg = torch.from_numpy(group_conv.regular_representation(self.mul_table)).float()\n",
    "        \n",
    "        # Find inverse of each element.\n",
    "        self.inverses_indices = group_conv.inverses(self.mul_table)\n",
    "        \n",
    "   \n",
    "        self.kernel = torch.nn.Parameter(torch.zeros(out_channels,\n",
    "                                                     in_channels,\n",
    "                                                     group_order, \n",
    "                                                     kernel_size,\n",
    "                                                     kernel_size))\n",
    "\n",
    "        torch.nn.init.kaiming_uniform_(self.kernel.data, a=math.sqrt(5))\n",
    "\n",
    "    def image2D_group_convolution_filter_bank(self, rep_2D, filter):\n",
    "        \"\"\"Creates rotated filter bank for 2D image convolution\n",
    "        Input:\n",
    "            rep_2D: torch.Tensor of shape [|G|, 2, 2] of the group representation as 2D rotations and mirrors\n",
    "            filter: torch.Tensor of shape [channel_out, channel_in, rep_reg_filter, kernel_height, kernel_width]\n",
    "        \"\"\"\n",
    "        (j, _, _) = rep_2D.shape\n",
    "        (d, c, i, kh, kw) = filter.shape\n",
    "        affine_mat = torch.zeros(j, 2, 3)\n",
    "        affine_mat[:, :2, :2] = rep_2D\n",
    "        grid = F.affine_grid(affine_mat, [j, ((d * c) * i), kh, kw], align_corners=False).float().to(filter.device)\n",
    "        filter_bank = F.grid_sample(filter.reshape(1, ((d * c) * i), kh, kw).repeat(j, 1, 1, 1), grid, align_corners=False)\n",
    "        return filter_bank.reshape(j, d, c, i, kh, kw)\n",
    "\n",
    "\n",
    "    def image2D_group_convolution(self, rep_2D, rep_reg, inverses, input, filter):\n",
    "        \"\"\"Performs group convolution of inputs and filters over the regular representation\n",
    "        Input:\n",
    "            rep_2D: torch.Tensor of shape [|G|, 2, 2] of the group representation as 2D rotations and mirrors\n",
    "            rep_reg: torch.Tensor of shape [|G|, |G|, |G|] of the left regular representation\n",
    "            inverse: torch.LongTensor of shape [|G|] with the indices of the inverses of the group elements\n",
    "            input: torch.Tensor of shape [batch, channel_in, rep_reg_in, height, width]\n",
    "            filter: torch.Tensor of shape [channel_out, channel_in, rep_reg_filter, kernel_height, kernel_width]\n",
    "        Output:\n",
    "            output: torch.Tensor of shape [batch, channel_out, rep_reg_out]\n",
    "        \"\"\"\n",
    "        (k, i, j) = rep_2D.shape\n",
    "        (z, c, i, h, w) = input.shape\n",
    "        (d, c, j, kh, kw) = filter.shape\n",
    "        filter_bank = self.image2D_group_convolution_filter_bank(torch.from_numpy(rep_2D).float(), filter)\n",
    "        filter_bank_contract = torch.einsum('kij,kdcjhw->dkcihw', rep_reg.to(filter.device), filter_bank[inverses])\n",
    "        conv = F.conv2d(input.reshape(z, (c * i), h, w), \n",
    "                        filter_bank_contract.reshape((d * k), (c * i), kh, kw), \n",
    "                        padding=(self.kernel_size-1)//2)\n",
    "        return conv.reshape(z, d, k, h, w)\n",
    "    \n",
    "    def forward(self, x):\n",
    "        if self.lift_conv:\n",
    "            x = x.unsqueeze(2).repeat(1, 1, self.group_order, 1, 1)\n",
    "        out = self.image2D_group_convolution(self.rep_2D, self.rep_reg, self.inverses_indices, x, self.kernel)\n",
    "        if self.activation:\n",
    "            return F.leaky_relu(out)\n",
    "        else:\n",
    "            return out\n",
    "        \n",
    "class Image_GConvNet(torch.nn.Module):\n",
    "    def __init__(self,\n",
    "                 in_channels,\n",
    "                 out_channels,\n",
    "                 kernel_size,\n",
    "                 hidden_dim,\n",
    "                 group_order,\n",
    "                 num_gconvs #number of group convolution layers.\n",
    "                 ):\n",
    "        super().__init__()\n",
    "\n",
    "        self.gconvs = []\n",
    "        # First Layer\n",
    "        self.gconvs.append(Image_GConv_Layer(in_channels, hidden_dim, kernel_size, group_order, lift_conv = True, activation = True))\n",
    "\n",
    "        # Middle Layers\n",
    "        \n",
    "        for i in range(num_gconvs-2):\n",
    "            self.gconvs.append(Image_GConv_Layer(hidden_dim, hidden_dim, kernel_size, group_order, lift_conv = False, activation = True))\n",
    "            \n",
    "        # Final Layer\n",
    "        self.gconvs.append(Image_GConv_Layer(hidden_dim, out_channels, kernel_size, group_order, lift_conv = False, activation = False))\n",
    "        self.gconvs = nn.Sequential(*self.gconvs)\n",
    "        \n",
    "    def forward(self, x):\n",
    "\n",
    "        out = self.gconvs(x)\n",
    "        \n",
    "        # functions on (g,x,y) -> functions on (x,y)\n",
    "        out = torch.mean(out, dim = 2)\n",
    "\n",
    "        # If we want to have a invariant classifer, we can average over the spatial dimensions.\n",
    "        out = torch.mean(out, dim = (2,3))\n",
    "        return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of paramters: 15200\n",
      "Epoch 1 | Train Accuracy: 0.17022 | Valid Accuracy: 0.19646 | Test Accuracy: 0.20332\n",
      "Epoch 2 | Train Accuracy: 0.22618 | Valid Accuracy: 0.34942 | Test Accuracy: 0.31628\n",
      "Epoch 3 | Train Accuracy: 0.35688 | Valid Accuracy: 0.58260 | Test Accuracy: 0.47383\n",
      "Epoch 4 | Train Accuracy: 0.58764 | Valid Accuracy: 0.75279 | Test Accuracy: 0.53643\n",
      "Epoch 5 | Train Accuracy: 0.71821 | Valid Accuracy: 0.83818 | Test Accuracy: 0.60689\n"
     ]
    }
   ],
   "source": [
    "GConvNet = Image_GConvNet(in_channels = 1,\n",
    "                          out_channels = 10,\n",
    "                          kernel_size = 5,\n",
    "                          hidden_dim = 8,\n",
    "                          group_order = 4,\n",
    "                          num_gconvs = 3).to(device)\n",
    "train_classifer(GConvNet, train_loader, valid_loader, test_loader, learning_rate = 0.01, num_epoch = 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "YQBueEWbDLlg"
   },
   "source": [
    "# 2.  Relaxed Group Convolution Networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "iUrJKGvtJUfU"
   },
   "source": [
    "## 2.1 Situations where perfectly equivariant models may not work"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lmDJuWxv3H7H"
   },
   "source": [
    "### 2.1.1 If the image of 6 is perfectly rotated version of 9, Can GCNN distinguish them?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "id": "WxRSGrUBJT7L",
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.5, 15.5, -0.5, 15.5)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAD1CAYAAADNj/Z6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAEmElEQVR4nO3dwU0CURhG0YFQBVXQhKECq7QCYxNWYRk+GzCKZPBl5p6zZvE2TG7+zXcYY4wFAMg6zn4AADCXGACAODEAAHFiAADixAAAxIkBAIgTAwAQJwYAIE4MAEDc6dYfPh2fH/kO4AZvny+zn/Bnvh0w32/fDpcBAIgTAwAQJwYAIE4MAECcGACAODEAAHFiAADixAAAxIkBAIgTAwAQJwYAIE4MAECcGACAODEAAHFiAADixAAAxIkBAIgTAwAQJwYAIE4MAECcGACAODEAAHFiAADixAAAxIkBAIgTAwAQJwYAIE4MAECcGACAODEAAHFiAADixAAAxIkBAIgTAwAQJwYAIE4MAECcGACAODEAAHFiAADixAAAxIkBAIgTAwAQJwYAIE4MAECcGACAODEAAHFiAADixAAAxIkBAIgTAwAQJwYAIE4MAECcGACAODEAAHFiAADiTrMfAFD1+vE++wm7cj1fZj9hs1wGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4qwWrsDy2Losj8HP/EdYm8sAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnNXCifa+PGbNEWAbXAYAIE4MAECcGACAODEAAHFiAADixAAAxIkBAIgTAwAQJwYAIE4MAECcGACAODEAAHFiAADirBYCbIxF0O/tfQn2kVwGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4qwWTmR5DLjH3tf5fBv/n8sAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnNXCFex9QexelscAtsFlAADixAAAxIkBAIgTAwAQJwYAIE4MAECcGACAODEAAHFiAADixAAAxIkBAIgTAwAQJwYAIM5qIcDGWARlbS4DABAnBgAgTgwAQJwYAIA4MQAAcWIAAOLEAADEiQEAiBMDABAnBgAgTgwAQJwYAIA4MQAAcVYLV2BBDLjH9XyZ/QRYlsVlAADyxAAAxIkBAIgTAwAQJwYAIE4MAECcGACAODEAAHFiAADixAAAxIkBAIgTAwAQJwYAIE4MAECcCeMVmCEFYMtcBgAgTgwAQJwYAIA4MQAAcWIAAOLEAADEiQEAiBMDABAnBgAgTgwAQJwYAIA4MQAAcWIAAOLEAADEiQEAiBMDABAnBgAgTgwAQJwYAIA4MQAAcWIAAOLEAADEiQEAiBMDABAnBgAgTgwAQJwYAIA4MQAAcWIAAOLEAADEiQEAiBMDABAnBgAgTgwAQJwYAIA4MQAAcWIAAOLEAADEiQEAiBMDABAnBgAgTgwAQJwYAIA4MQAAcWIAAOLEAADEiQEAiBMDABAnBgAgTgwAQJwYAIA4MQAAcWIAAOLEAADEiQEAiDuMMcbsRwAA87gMAECcGACAODEAAHFiAADixAAAxIkBAIgTAwAQJwYAIE4MAEDcFxyfKvViyvExAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_6, x_9 = np.zeros((16, 16)), np.zeros((16, 16))\n",
    "points_6 = [(5, 5),(5, 6),(5, 7),(5, 8),(5, 9),(5, 10),(6, 5),(6, 10),(7, 5),(7, 10),(8, 5),(8, 6),\n",
    "            (8, 7),(8, 8),(8, 9),(8, 10),(9, 5),(10, 5),(11, 5),(11, 6),(11, 7),(11, 8),(11, 9),(11, 10)]\n",
    "points_9 = [(4, 5),(4, 6),(4, 7),(4, 8),(4, 9),(4, 10),(5, 10),(6, 10),(7, 5),(7, 6),(7, 7),(7, 8),(7, 9),\n",
    "            (7, 10),(8, 5),(8, 10),(9, 5),(9, 10),(10, 5),(10, 6),(10, 7), (10, 8),(10, 9),(10, 10)]\n",
    "for x,y in points_6:\n",
    "    x_6[x,y] = 1.0\n",
    "for x,y in points_9:\n",
    "    x_9[x,y] = 1.0\n",
    "f, axarr = plt.subplots(1,2) \n",
    "axarr[0].imshow(x_9, origin=\"lower\")\n",
    "axarr[0].axis(\"off\")\n",
    "axarr[1].imshow(x_6, origin=\"lower\")\n",
    "axarr[1].axis(\"off\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def train_model(model, inp, tgt, loss_fun, lr = 1e-3, n_epochs = 200):\n",
    "    optimizer = torch.optim.Adam(model.parameters(), lr)\n",
    "    best_loss = 1e6\n",
    "    for epoch in range(n_epochs):\n",
    "        output = model(inp)\n",
    "        loss = loss_fun(output.reshape(tgt.shape), tgt)\n",
    "        optimizer.zero_grad()\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "         # print statistics\n",
    "        if epoch % max((n_epochs//10), 50) == 0:\n",
    "            print(\"Epoch: {} | Loss: {:0.5f} \".format(epoch,  loss.item()))\n",
    "    if output.device.type == 'cpu':\n",
    "        return output.data.numpy()\n",
    "    else:\n",
    "        return output.cpu().data.numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KFIBKyVxUJ09"
   },
   "source": [
    "**You'll see the loss doesn't converge because the C4-GCNN can't distinguish these two digits.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "id": "XD0WZaiSNiL1",
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch: 0 | Loss: 0.69318 \n",
      "Epoch: 50 | Loss: 0.69315 \n",
      "Epoch: 100 | Loss: 0.69315 \n",
      "Epoch: 150 | Loss: 0.69315 \n",
      "Epoch: 200 | Loss: 0.69315 \n",
      "Epoch: 250 | Loss: 0.69315 \n",
      "Epoch: 300 | Loss: 0.69315 \n",
      "Epoch: 350 | Loss: 0.69315 \n",
      "Epoch: 400 | Loss: 0.69315 \n",
      "Epoch: 450 | Loss: 0.69315 \n"
     ]
    }
   ],
   "source": [
    "GCNN = GroupEquivariantCNN(in_channels = 1,\n",
    "                           out_channels = 1,\n",
    "                           kernel_size = 3,\n",
    "                           hidden_dim = 4,\n",
    "                           group_order = 4,\n",
    "                           num_gconvs = 1, \n",
    "                           classifer = True, \n",
    "                           sigmoid = True)\n",
    "# input images 6 and 9\n",
    "inps = torch.stack([torch.from_numpy(x_6), torch.from_numpy(x_9)]).unsqueeze(1).float()\n",
    "# labels\n",
    "tgts = torch.FloatTensor([0,1])\n",
    "out = train_model(GCNN, inps, tgts, nn.BCELoss(), n_epochs = 500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1.2 When the outputs have lower symmetry than input.\n",
    "\n",
    "*Smidt et al. [Finding Symmetry Breaking Order Parameters with Euclidean Neural Networks](https://arxiv.org/abs/2007.02005)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAD1CAYAAADNj/Z6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAEqklEQVR4nO3YwWnDQBBAUSeoClWhJoIqSJWpwKQJV+Eyohx88yUJRFrQf+9s2FnWNp952bZtuwAAWa+jBwAAxhIDABAnBgAgTgwAQJwYAIA4MQAAcWIAAOLEAADEiQEAiJt++8G31/c95wB+4fPrY/QIf+a/A8b76b/DZgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgLhp9AAAPFzvt93PWOdl9zMul3PdpcBmAADixAAAxIkBAIgTAwAQJwYAIE4MAECcGACAODEAAHFiAADixAAAxIkBAIgTAwAQJwYAIE4MAECcGACAODEAAHFiAADixAAAxIkBAIgTAwAQJwYAIE4MAECcGACAODEAAHHT6AGOdr3fdj9jnZfdzzgTbwIPZ/qenukuBTYDABAnBgAgTgwAQJwYAIA4MQAAcWIAAOLEAADEiQEAiBMDABAnBgAgTgwAQJwYAIA4MQAAcWIAAOLEAADEiQEAiBMDABAnBgAgTgwAQJwYAIA4MQAAcWIAAOLEAADEiQEAiJtGD3C0dV5Gj8ATbwIwls0AAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAEDeNHuBo1/tt9zPWedn9jDPxJgBj2QwAQJwYAIA4MQAAcWIAAOLEAADEiQEAiBMDABAnBgAgTgwAQJwYAIA4MQAAcWIAAOLEAADEiQEAiBMDABAnBgAgTgwAQJwYAIA4MQAAcWIAAOLEAADEiQEAiBMDABAnBgAgbho9wNHWeRk9Ak+8CTxc77fdzzjq93amuxTYDABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBuGj0AAA/rvIwe4d+c6S4FNgMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMSJAQCIe9m2bRs9BAAwjs0AAMSJAQCIEwMAECcGACBODABAnBgAgDgxAABxYgAA4sQAAMR9A0l5LHfcJwJqAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# D4\n",
    "square_image = torch.zeros((9,9))\n",
    "square_image[2,2] = 1.0\n",
    "square_image[2,6] = 1.0\n",
    "square_image[6,2] = 1.0\n",
    "square_image[6,6] = 1.0\n",
    "square_image = square_image[None, None]\n",
    "\n",
    "# D2\n",
    "rect_image = torch.zeros((9,9))\n",
    "rect_image[1,3] = 1.0\n",
    "rect_image[1,5] = 1.0\n",
    "rect_image[7,3] = 1.0\n",
    "rect_image[7,5] = 1.0\n",
    "rect_image = rect_image[None, None]\n",
    "\n",
    "# Random\n",
    "rand_image = torch.zeros((9,9))\n",
    "rand_image[1,3] = 1.0\n",
    "rand_image[1,7] = 1.0\n",
    "rand_image[7,3] = 1.0\n",
    "rand_image[7,5] = 1.0\n",
    "rand_image = rand_image[None, None]\n",
    "\n",
    "f, axarr = plt.subplots(1,2) \n",
    "axarr[0].imshow(square_image[0,0], origin=\"lower\")\n",
    "axarr[0].axis(\"off\")\n",
    "axarr[1].imshow(rect_image[0,0], origin=\"lower\")\n",
    "axarr[1].axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch: 0 | Loss: 0.05063 \n",
      "Epoch: 50 | Loss: 0.03455 \n",
      "Epoch: 100 | Loss: 0.02485 \n",
      "Epoch: 150 | Loss: 0.02469 \n",
      "Epoch: 200 | Loss: 0.02469 \n",
      "Epoch: 250 | Loss: 0.02469 \n",
      "Epoch: 300 | Loss: 0.02469 \n",
      "Epoch: 350 | Loss: 0.02469 \n",
      "Epoch: 400 | Loss: 0.02469 \n",
      "Epoch: 450 | Loss: 0.02469 \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAGFCAYAAAASI+9IAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAFg0lEQVR4nO3XwY3aYBRGUWbkKlyFm4hcQapMBShNUAVlxNndLSyAB6Nz1l58+iXr6n0dx3GcAOB0On1PDwDgfYgCABEFACIKAEQUAIgoABBRACCiAECWez/89f37mTsAeLK///7c/MalAEBEAYCIAgARBQAiCgBEFACIKAAQUQAgogBARAGAiAIAEQUAIgoARBQAiCgAEFEAIKIAQEQBgIgCABEFACIKAEQUAIgoABBRACCiAEBEAYCIAgARBQAiCgBEFACIKAAQUQAgogBARAGAiAIAEQUAIgoARBQAiCgAEFEAIMv0AN7P+XqZnnDTvm7TE+7iLfk0LgUAIgoARBQAiCgAEFEAIKIAQEQBgIgCABEFACIKAEQUAIgoABBRACCiAEBEAYCIAgARBQAiCgBEFACIKAAQUQAgogBARAGAiAIAEQUAIgoARBQAiCgAEFEAIKIAQEQBgIgCABEFACIKAEQUAIgoABBRACCiAEBEAYCIAgBZpgfwfvZ1m57wY3hLPo1LAYCIAgARBQAiCgBEFACIKAAQUQAgogBARAGAiAIAEQUAIgoARBQAiCgAEFEAIKIAQEQBgIgCABEFACIKAEQUAIgoABBRACCiAEBEAYCIAgARBQAiCgBEFACIKAAQUQAgogBARAGAiAIAEQUAIgoARBQAiCgAEFEAIKIAQJbpAY90vl6mJ9y0r9v0BPg4/u3XcSkAEFEAIKIAQEQBgIgCABEFACIKAEQUAIgoABBRACCiAEBEAYCIAgARBQAiCgBEFACIKAAQUQAgogBARAGAiAIAEQUAIgoARBQAiCgAEFEAIKIAQEQBgIgCABEFACIKAEQUAIgoABBRACCiAEBEAYCIAgARBQAiCgBEFADIMj3gkfZ1m54APIF/+3VcCgBEFACIKAAQUQAgogBARAGAiAIAEQUAIgoARBQAiCgAEFEAIKIAQEQBgIgCABEFACIKAEQUAIgoABBRACCiAEBEAYCIAgARBQAiCgBEFACIKAAQUQAgogBARAGAiAIAEQUAIgoARBQAiCgAEFEAIKIAQEQBgIgCABEFALJMD3ik8/UyPeGmfd2mJ8DH8W+/jksBgIgCABEFACIKAEQUAIgoABBRACCiAEBEAYCIAgARBQAiCgBEFACIKAAQUQAgogBARAGAiAIAEQUAIgoARBQAiCgAEFEAIKIAQEQBgIgCABEFACIKAEQUAIgoABBRACCiAEBEAYCIAgARBQAiCgBEFACIKAAQUQAgogBARAGALNMDHmlft+kJwBP4t1/HpQBARAGAiAIAEQUAIgoARBQAiCgAEFEAIKIAQEQBgIgCABEFACIKAEQUAIgoABBRACCiAEBEAYCIAgARBQAiCgBEFACIKAAQUQAgogBARAGAiAIAEQUAIgoARBQAiCgAEFEAIKIAQEQBgIgCABEFACIKAEQUAIgoABBRACDL9ADez/l6mZ5w075u0xPu4i35NC4FACIKAEQUAIgoABBRACCiAEBEAYCIAgARBQAiCgBEFACIKAAQUQAgogBARAGAiAIAEQUAIgoARBQAiCgAEFEAIKIAQEQBgIgCABEFACIKAEQUAIgoABBRACCiAEBEAYCIAgARBQAiCgBEFACIKAAQUQAgogBARAGAiAIAWaYH8H72dZue8GN4Sz6NSwGAiAIAEQUAIgoARBQAiCgAEFEAIKIAQEQBgIgCABEFACIKAEQUAIgoABBRACCiAEBEAYCIAgARBQAiCgBEFACIKAAQUQAgogBARAGAiAIAEQUAIgoARBQAiCgAEFEAIKIAQEQBgIgCABEFACIKAEQUAIgoABBRACCiAEC+juM4pkcA8B5cCgBEFACIKAAQUQAgogBARAGAiAIAEQUAIgoA5D9LtCqbkN7iFAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "GCNN = GroupEquivariantCNN(in_channels = 1,\n",
    "                           out_channels = 1,\n",
    "                           kernel_size = 3,\n",
    "                           hidden_dim = 16,\n",
    "                           group_order = 4,\n",
    "                           num_gconvs = 3)\n",
    "out = train_model(GCNN, square_image, rect_image, nn.MSELoss(), n_epochs = 500)\n",
    "\n",
    "# visualize the prediction\n",
    "plt.imshow(out[0,0])\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1.3 When the physical systems do not have perfect symmetry or have lower symmetry than expected.\n",
    "\n",
    "* Wang et al. [Discovering Symmetry Breaking in Physical Systems with Relaxed Group Convolution](https://arxiv.org/abs/2310.02299).*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rvVCi9UmHoas"
   },
   "source": [
    "## 2.2 Approximate Equivariance\n",
    "Let $f \\colon X \\to Y$ be a function and $G$ be a group. Assume that $G$ acts on $X$ and $Y$ via representations $\\rho_{X}$ and $\\rho_{Y}$. We say $f$ is $\\epsilon$-approximately $G$-equivariant if for any $g \\in G$,\n",
    "\n",
    "$$ \\|f(\\rho_{X}(g)(x)) - \\rho_{Y}(g)f(x)\\| \\leq \\epsilon. $$\n",
    "\n",
    "Note that strictly equivariant functions are $\\epsilon = 0$ **approximately equivariant**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KNwaV2dyIG3U"
   },
   "source": [
    "## 2.3 Relaxed Group Convolution\n",
    "**We relax the weight-sharing constraints in Equivariant Networks by introducing group element dependent parameters.**\n",
    "\n",
    "Group Convolution: $$[f * \\psi] (g) =  \\sum_{h \\in G}f(h)\\psi(g^{-1} h)$$\n",
    "The $G$-equivariance of group convolution results from the shared kernel $\\psi(g^{-1}h)$.  To relax this and consequently relax the $G$-equivariance, we replace the single kernel $\\psi$ with a set of kernels $\\lbrace \\psi_l \\rbrace_{l=1}^L$. We define the new kernel $\\psi$ as a linear combination of $\\psi_l$ with coefficients that vary with $h$ and thus introduce symmetry-breaking dependence on the specific pair $(g,h)$,\n",
    "\n",
    "$$ \\psi(g,h) = \\sum_{l=1}^L w_l(h) \\psi_l(g^{-1}h). $$\n",
    "We define relaxed group convolution by multiplication with $\\psi$ as such\n",
    "\n",
    "$$[f {\\Large\\star} \\psi](g) = \\sum_{h\\in G}f(h)\\psi(g, h)\n",
    "  = \\sum_{h\\in G}\\sum_{l=1}^L f(h) w_l(h) \\psi_l(g^{-1}h) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "masbn5CZIc4w"
   },
   "source": [
    "By varying the number of kernels $L$, we can control the degree of relaxation of equivariance. The weights $w_i(h) \\in \\mathbb{R}$ and the kernels $\\psi_l(g^{-1}h) \\in \\mathbb{R}^{c_\\mathrm{out} \\times c_\\mathrm{in}}$ can be learnt from data.  Relaxed group convolution reduces to group convolution and is fully equivariant if and only if $g_1^{-1}h_1 = g_2^{-1} h_2$ implies $\\psi(g_1,h_1) = \\psi(g_2,h_2)$.  In particular, this occurs if $w_l(h_1) =  w_l(h_2)$ for all $h_1,h_2 \\in G$ and for all $l$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WPS02LvmLAbp"
   },
   "source": [
    "e.g. To relax translation equivariance, instead of using a single filter, we use a linear combination of multiple filters where the coefficients can be varying across the input location."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wHHYkI1tTX3B"
   },
   "source": [
    "![picture](https://drive.google.com/uc?export=view&id=1oNcmOFXuedlOkX-rQ8Pq4_fssEw4TOCm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oC2ZtqGCRzQ-"
   },
   "source": [
    "## 2.4 Implementing Relaxed Rotaion Equivariant Layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "##### 2D Relaxed Rotation Lifting Convolution Layer #####\n",
    "class RelaxedRotLiftConv2d(torch.nn.Module):\n",
    "    \"\"\"Relaxed lifting convolution Layer for 2D finite rotation group\"\"\"\n",
    "    def __init__(self,\n",
    "                 in_channels,\n",
    "                 out_channels,\n",
    "                 kernel_size,\n",
    "                 group_order, # the order of 2d finite rotation group\n",
    "                 num_filter_banks,\n",
    "                 activation = True # whether to apply relu in the end\n",
    "                 ):\n",
    "        super(RelaxedRotLiftConv2d, self).__init__()\n",
    "\n",
    "        self.num_filter_banks = num_filter_banks\n",
    "        self.in_channels = in_channels\n",
    "        self.out_channels = out_channels\n",
    "        self.kernel_size = kernel_size\n",
    "        self.group_order = group_order\n",
    "        self.activation = activation\n",
    "  \n",
    "        # The relaxed weights are initialized as equal\n",
    "        # they do not need to be equal across different filter bank\n",
    "        self.relaxed_weights = torch.nn.Parameter(torch.ones(num_filter_banks, group_order).float())\n",
    "\n",
    "        # Initialize an unconstrained kernel.\n",
    "        self.kernel = torch.nn.Parameter(torch.zeros(self.num_filter_banks, # Additional dimension\n",
    "                                                     self.out_channels,\n",
    "                                                     self.in_channels,\n",
    "                                                     self.kernel_size,\n",
    "                                                     self.kernel_size))\n",
    "        torch.nn.init.kaiming_uniform_(self.kernel.data, a=math.sqrt(5))\n",
    "        \n",
    "    def generate_filter_bank(self):\n",
    "        \"\"\" Obtain a stack of rotated filters\"\"\"\n",
    "        weights = self.kernel.reshape(self.num_filter_banks*self.out_channels,\n",
    "                                      self.in_channels,\n",
    "                                      self.kernel_size,\n",
    "                                      self.kernel_size)\n",
    "        filter_bank = torch.stack([rot_img(weights, -np.pi*2/self.group_order*i)\n",
    "                                   for i in range(self.group_order)])\n",
    "        filter_bank = filter_bank.transpose(0,1).reshape(self.num_filter_banks, # Additional dimension\n",
    "                                                         self.out_channels,\n",
    "                                                         self.group_order,\n",
    "                                                         self.in_channels,\n",
    "                                                         self.kernel_size,\n",
    "                                                         self.kernel_size)\n",
    "        return filter_bank\n",
    "\n",
    "\n",
    "    def forward(self, x):\n",
    "        # input shape: [bz, #in, h, w]\n",
    "        # output shape: [bz, #out, group order, h, w]\n",
    "\n",
    "        # generate filter bank given input group order\n",
    "        filter_bank = self.generate_filter_bank()\n",
    "\n",
    "        # for each rotation, we have a linear combination of multiple filters with different coefficients.\n",
    "        relaxed_conv_weights = torch.einsum(\"na, noa... -> oa...\", self.relaxed_weights, filter_bank)\n",
    "\n",
    "        # concatenate the first two dims before convolution.\n",
    "        # ==============================\n",
    "        x = F.conv2d(\n",
    "            input=x,\n",
    "            weight=relaxed_conv_weights.reshape(\n",
    "                self.out_channels * self.group_order,\n",
    "                self.in_channels,\n",
    "                self.kernel_size,\n",
    "                self.kernel_size\n",
    "            ),\n",
    "            padding = (self.kernel_size-1)//2\n",
    "        )\n",
    "        # ==============================\n",
    "\n",
    "        # reshape output signal to shape [bz, #out, group order, h, w].\n",
    "        # ==============================\n",
    "        x = x.view(\n",
    "            x.shape[0],\n",
    "            self.out_channels,\n",
    "            self.group_order,\n",
    "            x.shape[-1],\n",
    "            x.shape[-2]\n",
    "        )\n",
    "        # ==============================\n",
    "\n",
    "        if self.activation:\n",
    "            return F.leaky_relu(x)\n",
    "        return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "##### 2D Relaxed Rotation Group Convolution Layer #####\n",
    "class RelaxedRotGroupConv2d(torch.nn.Module):\n",
    "    \"\"\"Relaxed group convolution Layer for 2D finite rotation group\"\"\"\n",
    "    def __init__(self,\n",
    "                 in_channels,\n",
    "                 out_channels,\n",
    "                 kernel_size,\n",
    "                 group_order, # the order of 2d finite rotation group\n",
    "                 num_filter_banks,\n",
    "                 activation = True # whether to apply relu in the end\n",
    "                ):\n",
    "\n",
    "        super(RelaxedRotGroupConv2d, self).__init__()\n",
    "\n",
    "        self.num_filter_banks = num_filter_banks\n",
    "        self.in_channels = in_channels\n",
    "        self.out_channels = out_channels\n",
    "        self.kernel_size = kernel_size\n",
    "        self.group_order = group_order\n",
    "        self.activation = activation\n",
    "\n",
    "\n",
    "        # Initialize weights\n",
    "        # If relaxed_weights are equal values, then the model is still equivariant\n",
    "        # Relaxed weights do not need to be equal across different filter bank\n",
    "        self.relaxed_weights = torch.nn.Parameter(torch.ones(group_order, num_filter_banks).float())\n",
    "        self.kernel = torch.nn.Parameter(torch.randn(self.num_filter_banks, # additional dimension\n",
    "                                                     self.out_channels,\n",
    "                                                     self.in_channels,\n",
    "                                                     self.group_order,\n",
    "                                                     self.kernel_size,\n",
    "                                                     self.kernel_size))\n",
    "\n",
    "        torch.nn.init.kaiming_uniform_(self.kernel.data, a=math.sqrt(5))\n",
    "\n",
    "        \n",
    "    def generate_filter_bank(self):\n",
    "        \"\"\" Obtain a stack of rotated and cyclic shifted filters\"\"\"\n",
    "        filter_bank = []\n",
    "        weights = self.kernel.reshape(self.num_filter_banks*self.out_channels*self.in_channels,\n",
    "                                      self.group_order,\n",
    "                                      self.kernel_size,\n",
    "                                      self.kernel_size)\n",
    "\n",
    "        for i in range(self.group_order):\n",
    "            # planar rotation\n",
    "            rotated_filter = rot_img(weights, -np.pi*2/self.group_order*i)\n",
    "\n",
    "            # cyclic shift\n",
    "            shifted_indices = torch.roll(torch.arange(0, self.group_order, 1), shifts = i)\n",
    "            shifted_rotated_filter = rotated_filter[:,shifted_indices]\n",
    "\n",
    "\n",
    "            filter_bank.append(shifted_rotated_filter.reshape(self.num_filter_banks,\n",
    "                                                              self.out_channels,\n",
    "                                                              self.in_channels,\n",
    "                                                              self.group_order,\n",
    "                                                              self.kernel_size,\n",
    "                                                              self.kernel_size))\n",
    "        # stack\n",
    "        filter_bank = torch.stack(filter_bank).permute(1,2,0,3,4,5,6)\n",
    "        return filter_bank\n",
    "    \n",
    "    def forward(self, x):\n",
    "\n",
    "        filter_bank = self.generate_filter_bank()\n",
    "\n",
    "        relaxed_conv_weights = torch.einsum(\"na, aon... -> on...\", self.relaxed_weights, filter_bank)\n",
    "\n",
    "        x = torch.nn.functional.conv2d(\n",
    "            input=x.reshape(\n",
    "                x.shape[0],\n",
    "                x.shape[1] * x.shape[2],\n",
    "                x.shape[3],\n",
    "                x.shape[4]\n",
    "                ),\n",
    "            weight=relaxed_conv_weights.reshape(\n",
    "                self.out_channels * self.group_order,\n",
    "                self.in_channels * self.group_order,\n",
    "                self.kernel_size,\n",
    "                self.kernel_size\n",
    "            ),\n",
    "            padding = (self.kernel_size-1)//2\n",
    "        )\n",
    "\n",
    "        # Reshape signal back [bz, #out * g_order, h, w] -> [bz, out, g_order, h, w]\n",
    "        x = x.view(x.shape[0], self.out_channels, self.group_order, x.shape[-2], x.shape[-1])\n",
    "        # ========================\n",
    "        if self.activation:\n",
    "            return F.leaky_relu(x)\n",
    "        return x\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "class RelaxedRotCNN2d(torch.nn.Module):\n",
    "    \"\"\"A small relaxed rotation 2d CNN model\"\"\"\n",
    "    def __init__(self, \n",
    "                 in_channels, \n",
    "                 out_channels, \n",
    "                 kernel_size, \n",
    "                 hidden_dim, \n",
    "                 group_order, # the order of 2d finite rotation group\n",
    "                 num_gconvs, # number of group conv layers\n",
    "                 num_filter_banks, \n",
    "                 classifier = False,\n",
    "                 sigmoid = False,\n",
    "                ):\n",
    "        super(RelaxedRotCNN2d, self).__init__()\n",
    "        \n",
    "        self.gconvs = []\n",
    "        self.classifier = classifier\n",
    "        self.sigmoid = sigmoid\n",
    "        \n",
    "        self.gconvs = [RelaxedRotLiftConv2d(in_channels, hidden_dim, kernel_size, group_order, num_filter_banks, True)]\n",
    "\n",
    "        for i in range(num_gconvs-2):\n",
    "            self.gconvs.append(RelaxedRotGroupConv2d(hidden_dim, hidden_dim, kernel_size, group_order, num_filter_banks, True))\n",
    "\n",
    "        self.gconvs.append(RelaxedRotGroupConv2d(hidden_dim, out_channels, kernel_size, group_order, num_filter_banks, False))\n",
    "\n",
    "        self.gconvs = torch.nn.Sequential(*self.gconvs)\n",
    "\n",
    "    def forward(self, x): \n",
    "        # average over h axis or not         \n",
    "        out = self.gconvs(x).mean(2)\n",
    "        \n",
    "        if self.classifier:\n",
    "            out = out.mean((2,3))\n",
    "        if self.sigmoid:\n",
    "            out = out.sigmoid()\n",
    "        return out\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Equivariance Error with Rotation 0 Degrees: 0.00000\n",
      "Equivariance Error with Rotation 90 Degrees: 0.00000\n",
      "Equivariance Error with Rotation 180 Degrees: 0.00000\n",
      "Equivariance Error with Rotation 270 Degrees: 0.00000\n"
     ]
    }
   ],
   "source": [
    "Relaxed_GConvNet = RelaxedRotCNN2d(in_channels = 1, \n",
    "                                   out_channels = 1, \n",
    "                                   kernel_size = 3, \n",
    "                                   hidden_dim = 8, \n",
    "                                   group_order = 4,\n",
    "                                   num_gconvs = 3, \n",
    "                                   num_filter_banks=1)\n",
    "test_equivariance(Relaxed_GConvNet, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0a5H_QTdRRt3"
   },
   "source": [
    "#### Let's train a Relaxed GCNN classifier to distinguish image 9 and 6. Does the loss converge?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch: 0 | Loss: 0.69315 \n",
      "Epoch: 50 | Loss: 0.69229 \n",
      "Epoch: 100 | Loss: 0.64459 \n",
      "Epoch: 150 | Loss: 0.36013 \n",
      "Epoch: 200 | Loss: 0.08413 \n",
      "Epoch: 250 | Loss: 0.02209 \n",
      "Epoch: 300 | Loss: 0.00982 \n",
      "Epoch: 350 | Loss: 0.00564 \n",
      "Epoch: 400 | Loss: 0.00368 \n",
      "Epoch: 450 | Loss: 0.00260 \n"
     ]
    }
   ],
   "source": [
    "Relaxed_GConvNet = RelaxedRotCNN2d(in_channels = 1, \n",
    "                                   out_channels = 1, \n",
    "                                   kernel_size = 3, \n",
    "                                   hidden_dim = 8, \n",
    "                                   group_order = 4,\n",
    "                                   num_gconvs = 3, \n",
    "                                   num_filter_banks=1,\n",
    "                                   classifier = True,\n",
    "                                   sigmoid = True)\n",
    "out = train_model(Relaxed_GConvNet, inps, tgts, nn.BCELoss(), n_epochs = 500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.5 Discovering Symmetry Breaking Factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch: 0 | Loss: 0.04936 \n",
      "Epoch: 50 | Loss: 0.00428 \n",
      "Epoch: 100 | Loss: 0.00005 \n",
      "Epoch: 150 | Loss: 0.00001 \n",
      "Epoch: 0 | Loss: 0.04938 \n",
      "Epoch: 50 | Loss: 0.03696 \n",
      "Epoch: 100 | Loss: 0.03429 \n",
      "Epoch: 150 | Loss: 0.01700 \n",
      "Epoch: 0 | Loss: 0.04931 \n",
      "Epoch: 50 | Loss: 0.03713 \n",
      "Epoch: 100 | Loss: 0.03381 \n",
      "Epoch: 150 | Loss: 0.02903 \n"
     ]
    }
   ],
   "source": [
    "# Square to Square\n",
    "Relaxed_GConvNet_SS = RelaxedRotCNN2d(in_channels = 1, out_channels = 1, kernel_size = 3, hidden_dim = 4, group_order = 4, num_gconvs = 3, num_filter_banks=1).to(device)\n",
    "out_square = train_model(model = Relaxed_GConvNet_SS, inp = square_image.to(device), tgt = square_image.to(device), loss_fun = nn.MSELoss())\n",
    "\n",
    "# Square to Rectangle\n",
    "Relaxed_GConvNet_SR = RelaxedRotCNN2d(in_channels = 1, out_channels = 1, kernel_size = 3, hidden_dim = 8, group_order = 4, num_gconvs = 3, num_filter_banks=1).to(device)\n",
    "out_rectangle = train_model(model = Relaxed_GConvNet_SR, inp = square_image.to(device), tgt = rect_image.to(device), loss_fun = nn.MSELoss())\n",
    "\n",
    "# # Square to Non-symmetric shape\n",
    "Relaxed_GConvNet_SN = RelaxedRotCNN2d(in_channels = 1, out_channels = 1, kernel_size = 3, hidden_dim = 8, group_order = 4, num_gconvs = 3, num_filter_banks=1).to(device)\n",
    "out_random = train_model(model = Relaxed_GConvNet_SN, inp = square_image.to(device), tgt = rand_image.to(device), loss_fun = nn.MSELoss())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0UAAAGRCAYAAACqtxilAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABDI0lEQVR4nO3dd3RUdf7/8ddNhjRSIJRQJKCAgBgpIupaiIK4gAVFJKCiICgqiEZZ1oKJ6IourAVRFkTF7qrYkK+u4KoogiuLoDTBQpEeCAkBElI+vz/4zZWQQqaEKff5OOceh5l777znM84r933vnTuWMcYIAAAAABwqItAFAAAAAEAg0RQBAAAAcDSaIgAAAACORlMEAAAAwNFoigAAAAA4Gk0RAAAAAEejKQIAAADgaDRFAAAAAByNpggAAACAo9EUAQDgR9nZ2bIsS+np6YEuBYCHbrjhBlmWpRtuuCHQpeA4oykKI8YYvf3227riiivUsmVLxcbGKj4+Xq1bt9a5556rzMxMvffee8rPzw90qQD8xL0BfvQUHR2tZs2a6eKLL9asWbNUXFwc6FIlScuXL1d2draefPLJQJcCwAuhljlATbkCXQD8Y+/everfv7++/PJL+z6Xy6W4uDht2rRJv/76qxYtWqQnnnhCL774IntAgDCUkpJi3963b5+2bdumbdu26dNPP9WMGTP06aefqn79+gGs8HBT9OCDD6ply5a64447AloLAN+EQuYANcWRojAxdOhQffnll4qMjNRdd92ldevWqaioSLt379bBgwe1YsUKPfbYY+rUqVOgSwVQS7Zv325P+/fv18aNGzVy5EhJ0tKlS3X77bcHuEIA4YTMQTihKQoD69ev19y5cyVJDz/8sKZMmaK2bdsqIuLw2+tyuXTaaafpL3/5i5YvX65BgwYFslwAx0lqaqpmzpypCy+8UJL01ltvqaCgIMBVAQhXZA5CGU1RGFi+fLl9+/LLLz/m/LGxsRXuKy0t1dNPP62uXbuqbt26Sk5OVnp6ut555x1JUnp6uizLUnZ2drnlNmzYYJ9PvGHDhiqfs1WrVrIsS7Nnz67w2MqVK5Wdna0LL7xQrVu3VmxsrBITE9WlSxfdf//9ysnJqdF6CwoK9MADDygtLU0JCQmV1rRo0SJde+21atmypWJiYpSUlKTu3bvrscceI7gRtv785z9Lkg4dOqT169dXeHzfvn169NFHdfbZZys5OVnR0dFq0aKFMjIytHjx4mOu/9NPP1VGRob9Xcbk5GSddtppGjNmTLnlLcvSsGHDJEkbN26s8J2EI/PlwIEDeuONNzR06FB17txZjRo1sr+z0L9/f3388cdV1jN79mxZlqVWrVpJkv73v//p6quvVtOmTRUdHa2TTjpJmZmZys3NrfZ1LVy4UJdeeqkaNmyo2NhYtWvXTvfdd58KCgoqPIendu3apfvvv19dunRRUlKSYmJidNJJJ+nGG2/UqlWrvFonECyqy5zt27fr6aef1uWXX64OHTooKSlJsbGxatOmjUaMGFHt//9HXwThnXfeUXp6upKTkxUXF6fOnTvrqaeeUllZWbX1vfbaazrnnHOUkJCgpKQknXnmmZo5c6aMMTV6fe+++64uueQSpaSkKCoqSikpKbrkkkv03nvv1bj22bNn6+yzz1ZSUpLq16+vXr16aeHChfb8JSUlevrpp3X66acrMTFRSUlJ6tu3r5YtW1ajGuEFg5D31ltvGUlGkvn00089Xr6wsNBcfPHF9joiIiJMvXr1jGVZRpIZP3686dGjh5FksrKyyi3722+/2cv99ttvVT5Hy5YtjSTz4osvVvmYJBMTE2OSk5Pt55ZkmjdvbtauXVvteqdMmWJOPvlkI8lERUWZevXqlauptLTU3H777fY6JZn4+HgTGRlp/7tdu3Zmw4YNHo8fEEhZWVn2/8NVeeyxx+x5vvvuu3KPff/99+aEE06wH4+MjDQJCQn2vy3LMo888kil692/f78ZOHBguc9VQkKCSUpKsv/dqVMne/6UlBSTmJho50xKSkq5afLkyfa8L774YrkakpKSTFxcXLnnuuuuuyqty71sy5YtzWuvvWbq1KljJJmkpCQTERFhL9+xY0ezb9++StcxderUcjmUlJRkoqKijCTToUMH88QTT9jPUdV70qNHj0rXPX/+fDujJJk6deqYunXr2v+OiooyL730UqXLAoHma+Zcf/319mMul8skJycbl8tl3xcdHW3eeeedStfrXvb66683t912W7ltliOzYejQoZUuX1ZWZoYNG1YuW+rXr2/nQkZGRrnnOFpRUZEZNGhQue2lI5eXZAYPHmwOHTpUbe3u2y6Xq1zeulwuM3fuXFNYWGh69+5t58GR+RAXF2eWLl1azTsEb9EUhYHffvvN/uOdlpZmfvrpJ4+Wv/POO+1wePjhh01eXp4xxpgdO3aYW265xd4gqK2maOjQoWb27Nlm48aN9n1FRUVmwYIFpnv37kaS6dq1a7XrjY+PN02aNDHvvfeeHUabN282+/fvN8YYc//99xtJpnHjxuaZZ54xu3fvNsYYc+jQIfP555+bLl262M9TWlpa06EDAq4mGygXXnih/RnPycmx79+6datp3LixkWSuvPJKs3TpUvvzs2PHDjNhwgR7Y+W9996rsN6rr77a3jAYP3682bx5s/3Yrl27zGuvvWZGjRpVbpkjG5bqvP/+++buu+82X3/9tf05dtf84IMP2o3OBx98UGFZ93PExcWZ6OhoM2LECLNp0yZjzOFGbtq0afbyEyZMqLD8okWL7I2ciy66yM7U4uJi8/bbb5vk5GRTv359r5qiH374wcTGxhpJZuTIkWb16tWmpKTEGGPMxo0bza233mpvHB29MQkEA18yxxhjHnroITN58mTz448/muLiYmPM4R2XK1euNNdcc42RZOrWrWu2bNlSYb3uZqJ+/fomKirKPP744/Y2S05OjhkxYoRd22effVZh+aeeesp+fPTo0WbXrl3GGGP27t1rsrOzjWVZdoNVWVN011132a9rwoQJJjc31xhjzJ49e8y9995rr3v8+PFV1l6vXj0TGxtrZsyYYQ4cOGCMMWbt2rXm9NNPN5JMq1atzOjRo01ycrJ56623zKFDh0xZWZlZunSpad26tZFkzjnnnCrHHt6jKQoTI0eOLLfno0uXLubWW281zz//vPnxxx9NWVlZpctt2bLF3uipbOPAGGMGDx5sr7s2mqLq7Nu3z6SkpBhJ5quvvqpyvZGRkWbZsmWVruO3334zkZGRJjY21ixfvrzSefLz8+295ZVt/AHBqroNlI0bN5bLhssuu6zc48OHDzeSzJAhQ6pc/+OPP17hiI8xxixYsMBe77PPPlvjemvaFB3L5MmTjSTTs2fPKp+jqg0bY4zJzMw0kkybNm0qPNazZ08jyZxyyimmsLCwwuP/+c9/7PV72hS5NxbvueeeKl+b+6j25ZdfXuU8QKD4kjk10a9fPyPJPPTQQxUeO/IoU1XbE+7mYsSIEeXuP3jwoElOTjaSzHXXXVfpsn/961+rzI7ff//d3l6q6vPrzpU6deqYrVu3Vln7q6++WmHZn3/+udzRrsq2eT777DP78SN3QsE/aIrCRHFxsZkwYUK5Q6xHTo0bNzZ33nmn2b59e7nl3HtNYmNj7b0tR/vpp58C1hQZ88fe6EmTJlW53n79+lW5vDvAr7jiimqfZ/To0UZShT3bQDA7cgPlyFPRjj7VrH379ub333+3lzt48KCJiYkxksyKFSuqXH9OTo69jiPzY8iQIUaSOfXUUz2q119N0erVq+2jQe4jLUc/hySzfv36Spf/8ssv7XmOPBK1e/du+8h7dXl13nnnedwUufPS5XLZR6srs3TpUntv+dGvDQg0bzOnpp599lkjyVx88cUVHnM3Fi1atKhyZ++DDz5oJJnu3buXu/+DDz44Zi7s3bvXzsWjmyL39lJMTEyV20t79uwx0dHRRpKZOnVqpbWnpqZWWXubNm2MJHPeeedV+nhJSYm9/nnz5lU6D7zH7xSFCZfLpYkTJ+quu+7S3Llz9eWXX+q7777TmjVrdOjQIe3cuVNPPPGEXnnlFc2bN0/du3eXdPiSmZLUrVs3JSYmVrruk08+Wc2bN9eWLVtqrf6PPvpIr7zyir777jvt2LFDBw4cqDDP77//XuXy55xzTpWPLVq0SNLhL4M3adKkyvncF1rYuHFjTcsGgsqOHTsqvX/o0KGaMWOGYmJi7Pv+97//qbCwUJLUu3fvGq1/48aN9u+SfPPNN5KkSy65xJeSq7Vjxw49++yz+vTTT7Vu3Trl5eWptLS03DwHDhxQbm6uGjZsWGH55ORktWnTptJ1N2vWzL6dm5uruLg4SdL3339vf9m6R48eVdaWnp6ur776yqPX486isrIynXLKKVXO536N+/fv1+7du9W4cWOPngc4XjzJnCOtWLFCM2bM0Ndff60NGzaooKCgwkUOqvubf8YZZ8iyrEofc3+29+zZU+5+9/ZOixYtqsyFpKQknX766fZntbLlzzjjjCq3l+rXr69u3bpp0aJF9vxH69atW5W1p6Sk6Oeff9YZZ5xR6eORkZFq2LChtmzZcswLxcBzNEVhJikpSddee62uvfZaSVJhYaG+/vprTZ06VXPnzlVOTo4GDBig9evXKyYmRjt37pQkNW/evNr1nnDCCbXSFJWVlenaa6/VG2+8Yd/ncrlUv359RUVFSZLy8vJUWFio/fv3V7me6jYatm7dKunwBkZ163CrrCEDQoF7o8IYo+3bt+vDDz/UX//6V7388stKS0vT3Xffbc/r/lxIVW/YHO3Iz8b27dslSS1btvRH6RUsXrxYffv21d69e+374uPjFRcXJ8uyVFpaal+Zcv/+/ZU2RQkJCVWu3+X6489fcXGxfXvXrl327SMbp6MdKzMr4x7zsrIyr8YcCDaeZI7btGnTNHbsWPsKcZZlKSkpSdHR0ZKkgwcPKj8/v9q/1zX5bB/5uZbk0fZOZTxd3j3/0WpSuzevD77jktxhLiYmRr169dKHH36o66+/XtLhvS+ffPJJgCs77Pnnn9cbb7yhyMhIPfDAA1q/fr2Kioq0Z88e+wfhrrrqKkmq9lKZkZGRVT7m3us6fvx4mcOnjFY7ffHFF359jcDxZlmWmjZtqptvvlnvvfeeLMvSX/7yF/3nP/+x5znyiMvBgwdr9NlIT08v9xy1paSkRIMHD9bevXvVuXNn/d///Z/y8/O1b98+7dixQ9u3b9eSJUvs+avLBl/4+zW6xzwlJaVG422M8fqS38DxVJPMkaQ1a9bojjvuUFlZmQYOHKj//ve/KiwsVG5urv03//HHH5dUe59roCo0RQ5y00032bd/+uknSX8cYTnWUaCqHj9yb6v7VJzK5OXlVXr/m2++KUkaMWKEHnzwQbVp08b+0Vk39x5pb7lPmeO0ODhRenq6rrvuOhljNGbMGHvD/MhTSb35bNTm52rx4sXauHGjIiMj9dFHH6lPnz4V9pz6mgtVadSokX37yKNpR/PmyLl7zHJycmp01BoIRVVljnT4d4VKS0vVoUMHvfnmmzrjjDPss0Lcauuz7ev2jnv56k7rO/JxTnsNPTRFDhIfH2/fdh+m7tatm6TD58pW9eOl69evrzIE6tevb9/evHlzpfOsW7eu3CkwR3Iv06VLl0ofLygo0LffflvpYzXl/r7RggULqm3cgHD1wAMPKDIyUqtXr9ZLL70kSeU2RubOnevxOv/0pz95tax7p0d1e4HdudCoUaMqT1VZsGCBR89bU126dLGPEFV31NibI8ruLCotLa32x2eBUFdZ5kh/fLY7depUYQeoW219tt3bO5s3b9Yvv/xS6Tz5+fn63//+V+3yS5curXJH7969e8t99wihhaYoDPz2229at27dMec7Mpi6du0qSRowYIAiIyN18OBBTZkypdLlJk6cWOU669atq9atW0uS5syZU+k8f/vb36pcPikpSdLhL11W5qGHHtK+ffuqXL4mhg8fLpfLpZycHGVlZVU776FDh6psDoFQ1bp1aw0aNEjS4c9UcXGx6tatqyFDhkiSHnvsMW3atKnadRz9peUbb7xRkrRq1SpNnz69xrW4v6Bc1Y4S6Y9c2LFjR6Xfvfn99981derUGj+nJ5KTk3XBBRdIkv7xj3/o0KFDFeZZuHChxxdZkKS2bdvapyDed999VW5YuR095kCoqCxzpD8+2z/++GOlO0Y+/vjjWjuF/aKLLrJ35D700EOVzvP3v/9dBw8erPSxAQMGyOVyqbCwUI899lil8zzyyCMqKipSnTp1NGDAAP8UjuOGpigMrFq1Sh06dFC/fv308ssva8OGDfZjxcXF+v777zVs2DD7PN3u3bvr3HPPlXT4C4O33XabpMMhMWnSJLsJ2bVrl0aPHq1XX33VDrLKDB48WJL0wgsv6Nlnn7UDZfPmzRoxYoT+9a9/2Vd2Otqf//xnSdJzzz2nmTNn2hsg27dv15133qm///3vatCggbdDI+lwOE+YMEHS4cAbOnSoVq5caT9eUlKi5cuXa+LEiWrTpo2WL1/u0/MBweiee+6RZVnasGGDnn/+eUmH/4A3a9ZMOTk5Ovvss/XKK6+U2wmxa9cuzZkzR1dccYX9OXe74IILlJGRIUkaPXq07rnnnnJHlHNycjRr1iy7eXI79dRTJR3eI/vWW29VWuu5556runXryhijq6++2t7pU1paqn//+99KT0+v1e80Pfjgg7IsSytXrtRll12m9evXSzqcFe+++64GDBhQ7ii5J55++mnFx8dr3bp1Ouuss/TBBx+UO4K9ZcsWvfLKK+rZs6fGjx/vl9cDBEJlmeP+m79q1SrddtttduO/f/9+zZgxQ1dddZXPf/OrEhsba28LvPTSS7rjjju0e/duSYfz6KGHHtIjjzyievXqVbp88+bNNXbsWEnSo48+qqysLHvnzt69ezVhwgRNnjxZkpSZmammTZvWyutALaqta33j+Pnkk08q/C5RVFSUSU5Otn9vwz117dq1wq9EHzx40PTq1cueJzIy0tSvX99edvz48aZHjx6V/k6RMYd/YPWUU06xl4+IiLB/EbpOnTrmjTfeqPJ3inJzc0379u0rLOt+7ptvvtm+tn9lP8JY098/KisrMxMmTCg3HrGxsaZBgwYmMjKy3Bh9/fXXHr4DQODU5Nfl3S6//HIjyZxwwgn2j5KuXr3anHzyyeU+g8nJyRV+86xXr14V1rd//35z5ZVXlpsvMTHRJCUl2f8++kdfjfnjx1ElmYSEBNOyZUvTsmVL88QTT9jzTJ8+vdx64+Pj7d8Padiwofnwww+r/I20mvwW0rF+Y+2JJ54o9/z16tWzfx/k1FNPtR9v165dhWWr+/FWY4z5+uuvTZMmTcplboMGDUxsbGy55zz6xyeBYOBr5mRkZFT4bLn/Dp9++unm6aefrvLzW932gFt1n//S0lJz3XXXlcu7+vXr28+fkZFR7XMUFRXZv5145PIRERH2fYMHDzaHDh3yqvbqtrXcfPndR1SPI0Vh4OKLL9b69ev11FNPaeDAgerQoYOio6O1d+9excXFqW3btrr66qv15ptv6rvvvqtwmdmYmBh9/PHHeuqpp9S5c2dFRUXJGKPzzjtPb731lh599NFqnz8+Pl5ff/21MjMzdeKJJ8rlctmHjhcvXmzvTa5MvXr19M033+iOO+5Qq1atFBkZKZfLpfT0dL3xxhv65z//6ZcxsixLEydO1A8//KBbb71VHTp0UGRkpPLy8lS/fn396U9/0rhx4/TNN99U+5tHQCi77777JB0+/WzGjBmSpA4dOuiHH37QjBkz1Lt3bzVs2FD5+fkyxqhNmzYaOHCgZs6cWelRnbi4OM2ZM0cfffSRrrjiCjVr1kyFhYVyuVw67bTTdPvtt2vmzJkVlnvnnXd055136uSTT1ZxcbE2btyojRs3ljulbtSoUZo3b57S09MVHx+vkpISNW/eXGPGjNGKFSuUlpZWO4P0/91xxx364osv1LdvX9WvX1+FhYVq1aqV7r//fi1ZssQ+9aeqvcrVOeecc7Ru3TpNmTJF559/vurVq6e9e/cqMjJSHTp00LXXXqvXXntNTz75pH9fFHCcVZY57v+3TzvtNEVHR6u0tFRpaWmaNGmSFi1aVO77z/4WERGhl19+WS+//LLOOussxcbGqqSkRF27dtU///lPvf7669UuHxUVpX/9619655131KdPHzVo0ED79u1TgwYN1KdPH7377rt6/fXXVadOnVp7Dag9ljFc8xDHlp6eri+//FJZWVnKzs4OdDkAEFDXXHONXn/9dQ0fPtw+NQgAELo4UgQAgAfWrVund999V9If35EAAIQ2miIAAI7ywAMPaNq0adq0aZPKysokHf4y+L/+9S9dcMEFKiwsVPv27dW/f//AFgoA8AvXsWcBAMBZfvjhB33wwQcaM2aM6tSpo4SEBO3du9dukJo3b663336b7w4AQJigKQIA4Ch33nmnmjVrpm+++Ubbtm3Tnj17lJCQoJNPPlmXXHKJRo8ereTk5ECXCQDwEy60AAAAAMDR+E4RAAAAAEejKQIAAADgaDRFAAAAAByNpggAAACAo9EUAQAAAHA0miIAAAAAjkZTBAAAAMDRaIoAAAAAOJpjmqIvvvhClmVpypQpgS6l1syYMUPXXHON2rdvr8jISFmW5fW6wn28tmzZokmTJqlHjx5q2rSp6tatq44dO2rcuHHavXu3x+sL9/HauXOnhg0bptNOO03JycmKiYlRmzZtdOONN+rnn38OdHnHRbi/xxIZ4gkyxDNkyGHh/j5L5IgnyBHP1HaOuPxQI4LEpEmTtHv3bnXp0kX79+/X77//HuiSgtbcuXOVnZ2tfv36ady4cUpISNB///tfPfnkk3rzzTf13XffqUmTJoEuM2jk5uZq3bp16t27t1q2bKnY2FitX79eL7zwgt5++20tWbJEp5xySqDLhI/IkJojQzxDhjgHOVJz5IhnajtHaIpCRHFxsUpLSxUTE1PlPF988YVSU1MVERGhSy65xNFBdKzxOu+887Rx48ZyYTNy5EideeaZGjlypKZMmRK2e1oqc6zxateunRYtWlTh/quuukrdu3fXtGnT9Oyzz9Z2mfABGeIZMsQzZIgzkCOeIUc8E+gccczpczW1b98+3X///TrzzDPVsGFDRUdHq02bNvrrX/+qAwcO2PN9//33sixL9913X6Xr6devnxITE7V//377vm3btumWW25RamqqoqKi1KxZM910003auXNnuWWzs7NlWZZWrVqlzMxMnXDCCYqJidGSJUuqrb1Vq1aKiDi+b2mojlfHjh0r3fsyaNAgSdLKlSs9GoeaCtXxqkrLli0lHd57g8NC+T0mQ8gQMiQ4hPL7TI6QI6GaIxwpOsqWLVs0a9YsDRgwQEOGDJHL5dKXX36pv//97/r+++/173//W5LUpUsXnX766XrppZc0ceJERUZGllvHv//9bw0fPlx169aVJG3atElnn322Dh06pBtvvFGtW7fWzz//rOnTp+vzzz/X0qVLlZSUVK6Wa665RrGxsbrrrrtkWZaaNm16/AaihsJtvNx7tFJSUrwdkmqF+ngVFxcrLy9PxcXF+vnnn5WdnS1J6tu3r59GKPSF+nt8vIXbeJEh1SNDaibU3+fjLdzGixypXq3liHGIzz//3EgykydPrna+oqIic+jQoQr333///UaS+fbbb+37ZsyYYSSZefPmlZv34YcfrjDvZZddZho1amQ2b95cbt7vvvvOREZGmqysLPu+rKwsI8n06NHDFBcXe/Iybf369TO+vL1OGy+3gQMHGknms88+82g5p4zX3LlzjSR7SklJMf/4xz88Wkeocsp77EaGeIcMqZ6TM8QY57zPbuSId8iR6tVWjnD63FGioqJUp04dSVJJSYlyc3OVk5OjXr16SZK+/fZbe94hQ4YoPj5ezz//vH2fMUYvvPCC0tLS1L17d0lSXl6ePvroI1122WWKiYlRTk6OPbVq1Upt2rTRp59+WqGWO+64Qy5XcB/MC6fx+sc//qG3335bN910ky688EKv11OdUB+vs846S/Pnz9eHH36oRx99VE2bNlVubq5KSko8HotwFerv8fEWTuNFhhwbGVIzof4+H2/hNF7kyLHVWo743FaFiJp2z8YY88wzz5i0tDQTERFRrhOVZB588MFy844cOdLUqVPH7Ny50xhjzH/+8x8jyTz55JP2PN9++22F9Rw9nXTSSfb87u551apVXr/e47V3xpjwGK/nnnvOWJZl+vXrV+nek2Nx2ni5bdmyxTRo0MDcdNNNPq8r2DntPSZDPEOGeMdJGWKM895ncsQz5Ih3/JUjwd36B8Djjz+uu+66S71799btt9+uZs2aKSoqSlu2bNENN9ygsrKycvPfdNNNeu655/Tyyy/rrrvu0vPPP6/o6Ghdd9119jzGGEnStddeq+uvv77S542Nja1wX1xcnB9fWe0Ih/F64YUXdNNNN6l3796aM2eOvfekNoTDeB2pWbNm6tWrl55//nlNnTpV0dHRPq8z1IXbe1zbwmG8yBDvkSGVC7f3ubaFw3iRI97zV47QFB3llVdeUatWrfTxxx+Xu3rKJ598Uun83bp1U5cuXfT888/rxhtv1Jw5c9S/f38lJyfb87Rp00aWZenQoUP2oclwEerj9cILL2jEiBHq1auX3n///Vr/gxzq41WZgwcPqrS0VPn5+WrUqNFxf/5gE47vcW0K9fEiQ3xHhlQUju9zbQr18SJHfOePHOE7RUdx//qyu+OVDp9v+eijj1a5zMiRI7VmzRqNGTNGhYWFGjFiRLnHGzRooL59++rdd9+t9FKDxhjt2rXLfy/iOArl8Zo9e7ZGjhypCy+8UB988EG1v7vgL6E6Xjt27Kj0/tWrV+uzzz5T69at2Zj5/0L1PQ6UUB4vMqTmyBDPhOr7HCihPF7kSM3Vdo447kjRZ599psLCwgr3N2zYUKNGjdJVV12le+65R3369NGVV16p/Px8vf7669Uexrzmmms0btw4vfrqqzrxxBPVs2fPCvNMnz5d5557rs4//3wNHTpUXbp0UVlZmX799Vd98MEHGjp0qH1JQW/NnTtXK1askCT9/PPPkqSHH35YklSvXj2NHj3a43WG63h9+OGHuvHGG5WYmKhBgwZpzpw55R6Pj49X//79PV5vuI7XpEmTNH/+fPXr10+tWrWSMUYrV67UK6+8ouLiYj3zzDNerzvUhOt7LJEhniBDPEOGlBeu77NEjniCHPFMreeIT99ICiHuL59VNbVr184YY0xJSYl55JFHTOvWrU1UVJRJTU0148aNM6tXrzaSyl1O8EjDhw83kszEiROrrGHXrl3m7rvvNm3btjXR0dEmKSnJnHrqqeb2228v90Uz95fPfvvtN49e4/XXX1/l62vZsqVH6wr38XIvw3jVzPz5882AAQNMy5YtTWxsrImKijInnniiueGGG8zKlStrvJ5QFu7vsTFkCBlChtS2cH+fjSFHyJHQzRHHNEW17ZZbbjGRkZEVrs2OyjFenmG8wh/vsWcYL88wXs7A++wZxssz4T5eljFHnFAIr+Tl5alFixbq0aOH5s6dG+hygh7j5RnGK/zxHnuG8fIM4+UMvM+eYbw844Txctx3ivxp5cqV+v777/XSSy+poKBA9957b6BLCmqMl2cYr/DHe+wZxsszjJcz8D57hvHyjKPGK9CHqkKZ+3zI5s2bm+nTpwe6nKDHeHmG8Qp/vMeeYbw8w3g5A++zZxgvzzhpvDh9DgAAAICj8TtFAAAAABwtJL5TVFZWpq1btyohIUGWZQW6HMDRjDHat2+fmjVrVu6XsIMZGQIEj1DMEIkcAYJJbeRISDRFW7duVYsWLQJdBoAjbN68WSeccEKgy6gRMgQIPqGUIRI5AgQjf+ZISDRFCQkJkqSff9ushMTEAFcT/Dqd2l7bt21Tk6ZNtWLl2kCXE/QYL8+knXKydu7YYX8uQ4G71qhTrpcVGRXgaoJf0ZrXpJIDkitO0R2vC3Q5Qa9o1StSyQFFJzXQ+VlvBbqcoPdl9lU6lJ8bUhki/ZEj//3xF8WHWO2BkH5mJ+3csU3x9Rtp5MyPA11O0Hvupj4qyN3FtkgN1ca2SEg0Re7D1AmJiUqkKTom92HEiIgIxqsGGC/PuMcrlE4fcddqRUbRFNWE+721LFmR0YGtJRS4//+yIlQnNj7AxQQ/ywq9DJH+qDc+IYEdtDVg/62IiFB0HJ+LY7HYFvFIbWyLhM7JvAAAAABQC2iKAAAAADgaTREAAAAAR6MpAgAAAOBoNEUAAAAAHI2mCAAAAICj0RQBAAAAcDSaIgAAAACORlMEAAAAwNFoigAAAAA4Gk0RAAAAAEejKQIAAADgaDRFAAAAAByNpggAAACAo9EUAQAAAHA0miIAAAAAjkZTBAAAAMDRaIoAAAAAOBpNEQAAAABHoykCAAAA4Gg0RQAAAAAczeem6ODBg3K5XLIsS4888og/agLgMOQIAF+QIQB85XNT9P3336u0tFSS1K1bN58LAuA85AgAX5AhAHzlc1O0dOlS+zZBBMAb5AgAX5AhAHzlc1N0++23yxgjY4ySk5P9URMAhyFHAPiCDAHgKy60AAAAAMDRaIoAAAAAOJpPTdHatWtlWZYsy9Kbb77pr5oAOAg5AsAXZAgAf/CpKVqxYoV9u3Pnzr7WAsCByBEAviBDAPiDX5qi2NhYtW3b1i8FAXAWcgSAL8gQAP7gl6YoLS1NkZGRfikIgLOQIwB8QYYA8AefmqLly5dLkjp16uSPWgA4EDkCwBdkCAB/8LopysnJ0datWyVxDi8A75AjAHxBhgDwF6+boiO/2MjeGQDeIEcA+IIMAeAvPjdFlmXptNNO81tBAJyDHAHgCzIEgL943RS5z+E96aSTlJCQ4K96ADgIOQLAF2QIAH/x+UgR5/AC8BY5AsAXZAgAf/GqKTp06JDWrFkjiXN4AXiHHAHgCzIEgD951RStWbNGxcXFktg7A8A75AgAX5AhAPzJq6bIfQ6vxN4ZAN4hRwD4ggwB4E9eNUXuc3jr16+v1NRUvxYEwBnIEQC+IEMA+JNPTRF7ZgB4ixwB4AsyBIA/+dQUcQ4vAG+RIwB8QYYA8CeXNwvl5OT4uw4ADkOOAPAFGQLAn7z+nSIAAAAACAc0RQAAAAAcjaYIAAAAgKPRFAEAAABwNJoiAAAAAI5GUwQAAADA0WiKAAAAADgaTREAAAAAR6MpAgAAAOBoNEUAAAAAHI2mCAAAAICj0RQBAAAAcDSaIgAAAACORlMEAAAAwNFoigAAAAA4Gk0RAAAAAEejKQIAAADgaDRFAAAAAByNpggAAACAo9EUAQAAAHA0miIAAAAAjkZTBAAAAMDRaIoAAAAAOBpNEQAAAABHoykCAAAA4Gg0RQAAAAAcjaYIAAAAgKPRFAEAAABwNMsYYwJdxNGKiopUVFRk/zs/P18tWrRQk6ZNFRFBH3cs27dtU1lZmSIiItSkadNAlxP0GC/PbNu6VcYY5eXlKTExMdDlVKqqDJErTrKsAFYWIooPSDKSLKlOXKCrCX7u8bIiFJPUINDVBL3CvTmSgjtDpKpzpHEK2yI1sXPH4b+tVkSE4us3CnQ5Qa8gd5cM2yI1VivbIiYIZWVlGR3+i8zExBSkU15eXqCjokpkCBNT8E/BnCHGkCNMTKEw+TNHOFIUhjjy4RnGyzOhfKQopQkZUhM7tvOZ8IQ7QyRLiooPdDnB79A+SQrqDJGqzpH6DVPIkRrIzdnxx+eCI87HduQR+ph6AS4mBBTmSvJvjrj8shY/i46OVnR0dIX7V6xcG9QBGixatzpBW7dsUZOmTfXLht8DXU7QY7w8c2JqM23fti3QZVSrqgz5ZumPSiBDjqlT+xO1bSufiZpyZ4ii4hVz5p2BLifoFS55XCouCHQZx1RVjjz74VeKi08IQEWhZXivrtq9c5tUJ04xHW8IdDlBr3DVbKl4vxRTTzF9Jwe6nKBXOO9uqWivX9fJrg4AAAAAjkZTBAAAAMDRaIoAAAAAOBpNEQAAAABHoykCAAAA4Gg0RQAAAAAcjaYIAAAAgKPRFAEAAABwNJoiAAAAAI5GUwQAAADA0WiKAAAAADgaTREAAAAAR6MpAgAAAOBoNEUAAAAAHI2mCAAAAICj0RQBAAAAcDSaIgAAAACORlMEAAAAwNFoigAAAAA4Gk0RAAAAAEejKQIAAADgaDRFAAAAAByNpggAAACAo9EUAQAAAHA0miIAAAAAjkZTBAAAAMDRaIoAAAAAOJpfmqKysjK9+uqr6t27txo1aqS6deuqc+fOmjp1qkpLS1VQUCDLsmRZlqZMmeKPpwQQZsgRAL4gQwD4wuXrCnbu3Kn+/ftr8eLF5e5fsWKFxo4dq/nz5ysrK8u+v3Pnzr4+JYAwQ44A8AUZAsBXPjVFBQUFSk9P15o1a2RZlgYPHqyMjAw1b95cv/76qyZNmqSPPvpI+/fvt5chiAAciRwB4AsyBIA/+NQUjRo1SmvWrJHL5dKcOXN02WWX2Y917dpVffr0Ufv27fX5559Lkpo1a6aGDRv6VjGAsEKOAPAFGQLAH7z+TtFXX32l1157TZKUnZ1dLoTc6tatq7Fjx9r/Zs8MgCORIwB8QYYA8Bevm6KJEydKklJTUzV+/Pgq5+vYsaN9u1OnTt4+HYAwRI4A8AUZAsBfvGqKtm7dqgULFkiSbr31VrlcVZ+Fl5SUZN9m7wwAN3IEgC/IEAD+5FVT9Mknn9i3+/btW+28e/bssW+zdwaAGzkCwBdkCAB/8qopWr58uSQpOjpaaWlp1c67cuVKSVJcXJzatm3rzdMBCEPkCABfkCEA/Mmrpmj79u2SpEaNGh1z3vnz50uS0tLSFBHhl9+KBRAGyBEAviBDAPiTV8lQVFQkScrPz692vrVr19qXwOQcXgBHIkcA+IIMAeBPXjVFjRs3lnQ4iDZv3lzpPGVlZRozZoyMMZI4hxdAeeQIAF+QIQD8yaum6KyzzrJvT5o0qcLjpaWlGj16tH1VGIm9MwDKI0cA+IIMAeBPVV+/shoDBw7UuHHjlJubq+nTp6u4uFhDhgxRfHy8Vq9erWnTpmnp0qVKTU3Vpk2bZFnWMb8ECcBZyBEAviBDAPiTV01RYmKiZs2apUGDBqmkpESzZs3SrFmz/lipy6UJEyYoNzdX06ZNU8eOHRUfH++3ogGEPnIEgC/IEAD+5PUlWK688kotXLhQffv2Vb169RQTE6PWrVvr5ptv1rJlyzRx4kQtXLhQknT++ef7rWAA4YMcAeALMgSAv3h1pMjt7LPP1rx58yp9bPny5frhhx8kSQMGDPDlaQCEMXIEgC/IEAD+UCsX6y8rK1NmZqYkqV27drrgggtq42kAhDFyBIAvyBAAnvCqKfrll1+qfKywsFDDhg2zfxNg8uTJsizLu+oAhC1yBIAvyBAA/uT11eeio6OVkZGhzp07KykpSbm5uVqyZIlmzpypDRs2SJIyMzN16aWX+rNeAGGCHAHgCzIEgD953BSVlJRo9erVKioq0pIlSypfqcul7Oxs3XvvvT4XCCD8kCMAfEGGAPA3r44UzZ49W3PnztWyZcu0a9cu5eXlKTExUW3atFHPnj01atQopaam+rtWAGGEHAHgCzIEgD953BS5XC5lZGQoIyOjNuoB4ADkCABfkCEA/K1Wrj4HAAAAAKGCpggAAACAo9EUAQAAAHA0miIAAAAAjkZTBAAAAMDRaIoAAAAAOBpNEQAAAABHoykCAAAA4Gg0RQAAAAAcjaYIAAAAgKPRFAEAAABwNJoiAAAAAI5GUwQAAADA0WiKAAAAADgaTREAAAAAR6MpAgAAAOBoNEUAAAAAHI2mCAAAAICj0RQBAAAAcDSaIgAAAACORlMEAAAAwNFoigAAAAA4mmWMMYEu4mhFRUUqKiqy/52fn68WLVqoSdOmioigjzuW7du2qaysTBEREWrStGmgywl6jJdntm3dKmOM8vLylJiYGOhyKlVVhqQ0IUNqYsd2PhOecGeIZElR8YEuJ/gd2idJQZ0hUtU5Ur9hCjlSA7k5O/74XNSJC3Q5wa/4gCQjyZJi6gW4mBBQmCvJzzliglBWVpbR4f8zmJiYgnTKy8sLdFRUiQxhYgr+KZgzxBhyhIkpFCZ/5ghHisIQRz48w3h5JpSPFMkVJ1lWACsLEUfusWQP77G5x8uKUExSg0BXE/QK9+ZICu4MkarJkZh6ssiRYzIH90oyioiIUEoT/rYei/sIvawIxdZrGOhygt7B3F3yd464/LIWP4uOjlZ0dHSF+1esXBvUARosWrc6QVu3bFGTpk31y4bfA11O0GO8PHNiajNt37Yt0GVUq6oMie5wjazIqABUFFoKV82WivdLdeIUk3ZjoMsJeoU/Pi8V71dMUgNdOOmjQJcT9D77a18V5e0OdBnHVFWO1O0zSVad2ABUFFr2z71T5mCuUpo01Y/rNgS6nKCXdnIrbdu6RbH1GqrvlP8LdDlBb95df/7/O1j8h8MuAAAAAByNpggAAACAo9EUAQAAAHA0miIAAAAAjkZTBAAAAMDRaIoAAAAAOBpNEQAAAABHoykCAAAA4Gg0RQAAAAAcjaYIAAAAgKPRFAEAAABwNJoiAAAAAI5GUwQAAADA0WiKAAAAADgaTREAAAAAR6MpAgAAAOBoNEUAAAAAHC2gTdGWLVv05JNPqnfv3kpNTVVUVJSaNGmiAQMG6Ntvvw1kaQBCBDkCwBdkCABJcgXyyZ9++mk99thjat26tXr37q1GjRpp/fr1ev/99/X+++/r9ddf16BBgwJZIoAgV5Mc6dOnT6DLBBCk2BYBIAW4Kerevbu++OIL9ejRo9z9X331lXr27KlbbrlF/fv3D0xxAEJCTXLkp59+ClB1AIId2yIApACfPnfllVdWCCFJOu+883TBBRcoNzdXP/74YwAqAxAqapIjq1atCkBlAEIB2yIApCC+0EKdOnUkSS5XQA9mAQhh5AgAX5AhgHMEZVO0adMmLViwQE2bNlVaWlqgywEQgo7MkY4dOwa6HAAhhm0RwFmCrikqLi7Wddddp6KiIj322GOKjIwMdEkAQgw5AsAXZAjgPEHVFJWVlemGG27QwoULNXLkSF133XWBLglAiCFHAPiCDAGcKWiaorKyMg0fPlyvv/66rr32Wv3zn/8MdEkAQgw5AsAXZAjgXH5pisrKyvTqq6/a1/evW7euOnfurKlTp6q0tFQFBQWyLEuWZWnKlCmVLj9s2DC99NJLGjx4sGbPnq2IiKDp1wAcB+QIAF+QIQB84fPlVHbu3Kn+/ftr8eLF5e5fsWKFxo4dq/nz5ysrK8u+v3PnzuXmc4fQyy+/rEGDBumVV17h3F3AYcgRAL4gQwD4yqemqKCgQOnp6VqzZo0sy9LgwYOVkZGh5s2b69dff9WkSZP00Ucfaf/+/fYyRwaR+zD1yy+/rIEDB+rVV18lhACHIUcA+IIMAeAPPjVFo0aN0po1a+RyuTRnzhxddtll9mNdu3ZVnz591L59e33++eeSpGbNmqlhw4b2PBMnTtRLL72k+Ph4nXzyyXr44YcrPEf//v110kkn+VImgCB2PHKkV69etf9CAAQE2yIA/MHrpuirr77Sa6+9JknKzs4uF0JudevW1dixYzVu3DhJFQ9Xb9iwQdLhvTx/+9vfKn2eVq1aEURAmDpeOZKSkuK/ogEEDbZFAPiL198gnDhxoiQpNTVV48ePr3K+I380sVOnTuUemz17towx1U433HCDtyUCCHLHK0euueaa2nkBAAKKbREA/uJVU7R161YtWLBAknTrrbfK5ar6gFNSUpJ9++i9MwCcixwB4AsyBIA/edUUffLJJ/btvn37Vjvvnj177NtH750B4FzkCABfkCEA/Mmrpmj58uWSpOjoaKWlpVU778qVKyVJcXFxatu2rTdPByAMkSMAfEGGAPAnr5qi7du3S5IaNWp0zHnnz58vSUpLS+NH0ADYyBEAviBDAPiTV8lQVFQkScrPz692vrVr19qXwOQcXgBHIkcA+IIMAeBPXjVFjRs3lnQ4iDZv3lzpPGVlZRozZoyMMZI4hxdAeeQIAF+QIQD8yaum6KyzzrJvT5o0qcLjpaWlGj16tH1VGIm9MwDKI0cA+IIMAeBPXv1468CBAzVu3Djl5uZq+vTpKi4u1pAhQxQfH6/Vq1dr2rRpWrp0qVJTU7Vp0yZZlnXML0ECcBZyBIAvyBAA/uRVU5SYmKhZs2Zp0KBBKikp0axZszRr1qw/VupyacKECcrNzdW0adPUsWNHxcfH+61oAKGPHAHgCzIEgD95fQmWK6+8UgsXLlTfvn1Vr149xcTEqHXr1rr55pu1bNkyTZw4UQsXLpQknX/++X4rGED4IEcA+IIMAeAvXh0pcjv77LM1b968Sh9bvny5fvjhB0nSgAEDfHkaAGGMHAHgCzIEgD/UysX6y8rKlJmZKUlq166dLrjggtp4GgBhjBwB4AsyBIAnvGqKfvnllyofKyws1LBhw+zfBJg8ebIsy/KuOgBhixwB4AsyBIA/eX31uejoaGVkZKhz585KSkpSbm6ulixZopkzZ2rDhg2SpMzMTF166aX+rBdAmCBHAPiCDAHgTx43RSUlJVq9erWKioq0ZMmSylfqcik7O1v33nuvzwUCCD/kCABfkCEA/M2rI0WzZ8/W3LlztWzZMu3atUt5eXlKTExUmzZt1LNnT40aNUqpqan+rhVAGCFHAPiCDAHgTx43RS6XSxkZGcrIyKiNegA4ADkCwBdkCAB/q5WrzwEAAABAqKApAgAAAOBoNEUAAAAAHI2mCAAAAICj0RQBAAAAcDSaIgAAAACORlMEAAAAwNFoigAAAAA4Gk0RAAAAAEejKQIAAADgaDRFAAAAAByNpggAAACAo9EUAQAAAHA0miIAAAAAjkZTBAAAAMDRaIoAAAAAOBpNEQAAAABHoykCAAAA4Gg0RQAAAAAcjaYIAAAAgKPRFAEAAABwNJoiAAAAAI5GUwQAAADA0SxjjAl0EUcrKipSUVGR/e/8/Hy1aNFCTZo2VUQEfdyxbN+2TWVlZYqIiFCTpk0DXU7QY7w8s23rVhljlJeXp8TExECXU6mqMkSuOMmyAlhZiCg+IMlIsqQ6cYGuJvi5x8uKUExSg0BXE/QK9+ZICu4MkarJkZh6ssiRYzIH90oyioiIUEoT/rYey47th7dFZEUotl7DQJcT9A7m7pLfc8QEoaysLKPDf5GZmJiCdMrLywt0VFSJDGFiCv4pmDPEGHKEiSkUJn/mCEeKwhBHPjzDeHkmlI8UpTQhQ2rCvceSz0TNuDNEsqSo+ECXE/wO7ZOkoM4QqZojRVEJHHGuiaJ9koysiAglJDcKdDVBb9+eXTLkSM3VQo64/LIWP4uOjlZ0dHSF+1esXBvUARosWrc6QVu3bFGTpk31y4bfA11O0GO8PHNiajNt37Yt0GVUq6oM+Wbpj0ogQ46pU/sTtW0rn4macmeIouIVc+adgS4n6BUueVwqLgh0GcdUVY5En3mHLFdMACoKLYXfPCYV5SshuZHufv3rQJcT9KYMOVf5OTsO58jptwW6nKBXuHSa33OEXaYAAAAAHI2mCAAAAICj0RQBAAAAcDSaIgAAAACORlMEAAAAwNFoigAAAAA4Gk0RAAAAAEejKQIAAADgaDRFAAAAAByNpggAAACAo9EUAQAAAHA0miIAAAAAjkZTBAAAAMDRaIoAAAAAOBpNEQAAAABHoykCAAAA4Gg0RQAAAAAcLaBNUWFhoTIzM3X++eerWbNmiomJUZMmTXTOOefoxRdfVHFxcSDLAxACyBEAviBDAEgBbooKCgo0ffp0WZalfv36KTMzU1dccYW2bNmi4cOH65JLLlFZWVkgSwQQ5MgRAL4gQwBIkiuQT56cnKy8vDxFRUWVu7+kpEQXXXSRPv30U3388cc677zzAlQhgGBXkxyZP39+gKoDEOzYFgEgBfhIUURERIUQkiSXy6UrrrhCkvTzzz8f77IAhJCa5Mivv/56vMsCECLYFgEgBemFFsrKyvTJJ59Ikk499dQAVwMgFB2ZIx06dAhwNQBCDdsigLME9PQ5t0OHDumRRx6RMUa7d+/WZ599prVr12rYsGHq2bOn8vPzA10igCBXXY6kp6cHujwAQY5tEcDZgqYpevDBB+1/W5alu+++W5MmTQpgVQBCSXU5cuDAgQBWBiAUsC0COFtQnD4XHx8vY4xKS0u1efNmPfPMM5o1a5bS09PZMwOgRsgRAL4gQwBnC4qmyC0iIkInnHCCbrnlFs2cOVOLFi3S3/72t0CXBSCEVJYjU6ZMCXRZAEIE2yKAM/mlKSorK9Orr76q3r17q1GjRqpbt646d+6sqVOnqrS0VAUFBbIsS5Zl1XjjpHfv3pKkL774wh8lAghytZkjX3/9dW2WDiAIsC0CwBc+f6do586d6t+/vxYvXlzu/hUrVmjs2LGaP3++srKy7Ps7d+5co/Vu3bpVklSnTh1fSwQQ5MgRAL4gQwD4yqemqKCgQOnp6VqzZo0sy9LgwYOVkZGh5s2b69dff9WkSZP00Ucfaf/+/fYyRwbR6tWr1apVK8XFxZVb74EDB5SZmSlJ6tu3ry8lAghyxyNHLrroIi1ZsuS4vB4AxxfbIgD8waemaNSoUVqzZo1cLpfmzJmjyy67zH6sa9eu6tOnj9q3b6/PP/9cktSsWTM1bNjQnuett97S448/rnPPPVetWrVSYmKitmzZoo8//li7d+/WeeedpzvvvFPFxcW+lAkgiB2PHLntttv00EMPHffXBqD2sS0CwB+8boq++uorvfbaa5Kk7OzsciHkVrduXY0dO1bjxo2TVPFw9SWXXKKtW7fqm2++0eLFi1VQUKCkpCSddtppysjI0PDhw+VyuQgiIEwdrxzhktxAeGJbBIC/eN0UTZw4UZKUmpqq8ePHVzlfx44d7dudOnUq91i3bt3UrVs3b0sAEOLIEQC+IEMA+ItXV5/bunWrFixYIEm69dZb5XJV3VslJSXZt2v6xUYA4Y8cAeALMgSAP3nVFH3yySf27WN9+XDPnj327aP3zgBwLnIEgC/IEAD+5FVTtHz5cklSdHS00tLSqp135cqVkqS4uDi1bdvWm6cDEIbIEQC+IEMA+JNXTdH27dslSY0aNTrmvPPnz5ckpaWlKSLCL78VCyAMkCMAfEGGAPAnr5KhqKhIkpSfn1/tfGvXrrUvgck5vACORI4A8AUZAsCfvGqKGjduLOlwEG3evLnSecrKyjRmzBgZYyRxDi+A8sgRAL4gQwD4k1dN0VlnnWXfnjRpUoXHS0tLNXr0aPuqMBJ7ZwCUR44A8AUZAsCfvPqdooEDB2rcuHHKzc3V9OnTVVxcrCFDhig+Pl6rV6/WtGnTtHTpUqWmpmrTpk2yLOuYX4IE4CzkCABfkCEA/MmrpigxMVGzZs3SoEGDVFJSolmzZmnWrFl/rNTl0oQJE5Sbm6tp06apY8eOio+P91vRAEIfOQLAF2QIAH/y+hIsV155pRYuXKi+ffuqXr16iomJUevWrXXzzTdr2bJlmjhxohYuXChJOv/88/1WMIDwQY4A8AUZAsBfvDpS5Hb22Wdr3rx5lT62fPly/fDDD5KkAQMG+PI0AMIYOQLAF2QIAH+olYv1l5WVKTMzU5LUrl07XXDBBbXxNADCGDkCwBdkCABPeNUU/fLLL1U+VlhYqGHDhtm/CTB58mRZluVddQDCFjkCwBdkCAB/8vrqc9HR0crIyFDnzp2VlJSk3NxcLVmyRDNnztSGDRskSZmZmbr00kv9WS+AMEGOAPAFGQLAnzxuikpKSrR69WoVFRVpyZIlla/U5VJ2drbuvfdenwsEEH7IEQC+IEMA+JtXR4pmz56tuXPnatmyZdq1a5fy8vKUmJioNm3aqGfPnho1apRSU1P9XSuAMEKOAPAFGQLAnzxuilwulzIyMpSRkVEb9QBwAHIEgC/IEAD+VitXnwMAAACAUEFTBAAAAMDRaIoAAAAAOBpNEQAAAABHoykCAAAA4Gg0RQAAAAAcjaYIAAAAgKPRFAEAAABwNJoiAAAAAI5GUwQAAADA0WiKAAAAADgaTREAAAAAR6MpAgAAAOBoNEUAAAAAHI2mCAAAAICj0RQBAAAAcDSaIgAAAACORlMEAAAAwNFoigAAAAA4Gk0RAAAAAEejKQIAAADgaDRFAAAAAByNpggAAACAo1nGGBPoIo5WVFSkoqIi+995eXlKTU1V45QURUTQxx3Lju3bZYyRZVlKadIk0OUEPcbLM9u3bZMk7d27V0lJSQGupnJVZUijxmRITezcwWfCE+4MkSTViQ9sMaGguEBScGeIVHWOqE68ZFkBrCxEHNp3+L+WpYT6DQNbSwjYl5sjkSM1Vxs5YoJQVlaWkcTExBTE0y+//BLoqKgSGcLEFPxTMGeIMeQIE1MoTP7MkZA4UlRWVqY9e/aoQYMGstg7c0z5+flq0aKFNm/erMTExECXE/QYL8+495bm5uaqXr16gS6nUmSIb/hMeIbx8kwoZIhEjviKz4VnGC/P1EaOuPyyFj+Ljo5WdHR0ufuCOTiDVWJiIh8sDzBengnm09DIEP/gM+EZxsszwZwhEjniL3wuPMN4ecafORLciQQAAAAAtYymCAAAAICj0RSFoejoaGVlZVU47I/KMV6eYbzCH++xZxgvzzBezsD77BnGyzO1MV5BeaEFAAAAADheOFIEAAAAwNFoigAAAAA4Gk0RAAAAAEejKQIAAADgaDRFAAAAAByNpggAAACAo9EUAQAAAHA0miIAAAAAjkZTBAAAAMDRaIoAAAAAOBpNEQAAAABHoykCAAAA4Gg0RQAAAAAcjaYIAAAAgKPRFAEAAABwNJoiAAAAAI5GUwQAAADA0WiKAAAAADgaTREAAAAAR6MpAgAAAOBoNEUAAAAAHI2mCAAAAICj0RQBAAAAcDSaIgAAAACORlMEAAAAwNFoigAAAAA4Gk0RAAAAAEf7fzuJjYUK/JXwAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1000x500 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "relaxed_weights_ss = torch.stack([Relaxed_GConvNet_SS.gconvs[i].relaxed_weights.data.reshape(-1) for i in range(3)]).T.cpu().data.numpy()\n",
    "relaxed_weights_sr = torch.stack([Relaxed_GConvNet_SR.gconvs[i].relaxed_weights.data.reshape(-1) for i in range(3)]).T.cpu().data.numpy()\n",
    "relaxed_weights_sn = torch.stack([Relaxed_GConvNet_SN.gconvs[i].relaxed_weights.data.reshape(-1) for i in range(3)]).T.cpu().data.numpy()\n",
    "relaxed_weights = [relaxed_weights_ss, relaxed_weights_sr, relaxed_weights_sn]\n",
    "titles = [\"Square\", \"Rectangle\", \"Random\"]\n",
    "fig=plt.figure(figsize=(10, 5))\n",
    "for i in range(1, 4):\n",
    "    fig.add_subplot(1, 3, i)\n",
    "    plt.imshow(relaxed_weights[i-1], cmap = 'Blues')\n",
    "    plt.xticks(np.arange(0.5,3,1), color='w')\n",
    "    plt.yticks(np.arange(0.5,4,1), color='w')\n",
    "    plt.text(-0.9, 0.1, \"$i$\", fontsize=20)\n",
    "    plt.text(-0.9, 1.1, \"$g$\", fontsize=20)\n",
    "    plt.text(-0.9, 2.1, \"$g^2$\", fontsize=20)\n",
    "    plt.text(-0.9, 3.1, \"$g^3$\", fontsize=20)\n",
    "    plt.text(-0.3, -0.6, \"Layer 1\", fontsize=13)\n",
    "    plt.text(0.7, -0.6, \"Layer 2\", fontsize=13)\n",
    "    plt.text(1.7, -0.6, \"Layer 3\", fontsize=13)\n",
    "    plt.grid(color='black', linestyle='-', linewidth=2)\n",
    "    plt.title(titles[i-1], size = 18, pad=25)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.6 Projecting the relaxed weights onto the irreps of the group"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Character Table\n",
    "c4_ft = np.array([[1, 1, 1, 1],\n",
    "                  [1, -1, 1, -1],\n",
    "                  [1, 0+1j, -1, 0-1j],\n",
    "                  [1, 0-1j, -1, 0+1j]])\n",
    "\n",
    "fourier_components_ss = np.round(relaxed_weights_ss[:,0] @ c4_ft.T, decimals = 3)\n",
    "fourier_components_sr = np.round(relaxed_weights_sr[:,0] @ c4_ft.T, decimals = 3)\n",
    "fourier_components_sn = np.round(relaxed_weights_sn[:,0] @ c4_ft.T, decimals = 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Only trivial irrep A is non zero when the output is an square"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.238+0.j, 0.   +0.j, 0.   +0.j, 0.   +0.j])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fourier_components_ss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Both A and B are non zero when the output is a rectangle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.192+0.j, 0.785+0.j, 0.   +0.j, 0.   +0.j])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fourier_components_sr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "e3nn",
   "language": "python",
   "name": "e3nn"
  },
  "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.9.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}