- {
- “cells”: [
- {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“# Constructing Equivariant Models”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Previously we showed examples of finding equivariant bases for different groups and representations, now we’ll show how these bases can be assembled into equivariant neural networks such as EMLP. n”, “n”, “We will give examples at a high level showing how the specific EMLP model can be applied to different groups and input-output types, and later in a lower level showing how models like EMLP can be constructed with equivariant layers and making use of the equivariant bases.”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## Using EMLP with different groups and representations (high level)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“A basic EMLP is a sequence of EMLP layers (containing G-equivariant linear layers, bilinear layers incorporated with a shortcut connection, and gated nonlinearities. While our numerical equivariance solver can work with any finite dimensional linear representation, for EMLP we restrict ourselves to _tensor_ representations.”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“By tensor representations, we mean all representations which can be formed by arbitrary combinations of $\oplus$,$\otimes$,$^*$ (+,`*`,`.T`) of a base representation $\rho$. This is useful because it simplifies the construction of our bilinear layer, which is a crucial ingredient for expressiveness and universality in EMLP.n”, “n”, “Following the $T_{(p,q)}=V^{\otimes p}\otimes (V^*)^{\otimes p}$ notation in the paper, we provide the convenience function for constructing higher rank tensors.”
]
}, {
“cell_type”: “code”, “execution_count”: 1, “metadata”: {}, “outputs”: [
- {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“V⊗V⊗V*⊗V*⊗V*n”, “V²⊗V*³n”
]
}
], “source”: [
“from emlp.reps import Vn”, “from emlp.groups import SO13n”, “n”, “def T(p,q=0):n”, ” return (V**p*V.T**q)n”, “n”, “print(T(2,3))n”, “print(T(2,3)(SO13()))”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Lets get started with a toy dataset: learning how an inertia matrix depends on the positions and masses of 5 point masses distributed in different ways. The data consists of mappings (positions, masses) –> (inertia matrix) pairs, and has an $G=O(3)$ symmetry (3D rotation and reflections). If we rotate all the positions, the resulting inertia matrix should be correspondingly rotated.”
]
}, {
“cell_type”: “code”, “execution_count”: 2, “metadata”: {}, “outputs”: [
- {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“Input type: 5V⁰+5V, output type: V²n”
]
}
], “source”: [
“from emlp.datasets import Inertian”, “from emlp.groups import SO,O,S,Zn”, “n”, “trainset = Inertia(1000) # Initialize dataset with 1000 examplesn”, “testset = Inertia(2000)n”, “G = SO(3)n”, “print(f"Input type: {trainset.rep_in(G)}, output type: {trainset.rep_out(G)}")”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“For convenience, we store in the dataset the types for the input and the output. 5V⁰ are the $5$ mass values and 5V are the position vectors of those masses, V² is the matrix type for the output, equivalent to $T_2$. To initialize the [EMLP](https://emlp.readthedocs.io/en/latest/package/emlp.nn.html#emlp.nn.EMLP), we just need these input and output representations, the symmetry group, and the size of the network as parametrized by number of layers and number of channels (the dimension of the feature representation).”
]
}, {
“cell_type”: “code”, “execution_count”: 3, “metadata”: {}, “outputs”: [], “source”: [
“import emlpn”, “model = emlp.nn.EMLP(trainset.rep_in,trainset.rep_out,group=G,num_layers=3,ch=384)n”, “# uncomment the following line to instead try the MLP baselinen”, “#model = emlp.nn.MLP(trainset.rep_in,trainset.rep_out,group=G,num_layers=3,ch=384)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## Example Objax Training Loop”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“We build our EMLP model with [objax](https://objax.readthedocs.io/en/latest/) because we feel the object oriented design makes building complicated layers easier. Below is a minimal training loop that you could use to train EMLP.”
]
}, {
“cell_type”: “code”, “execution_count”: 4, “metadata”: {}, “outputs”: [], “source”: [
“BS=500n”, “lr=3e-3n”, “NUM_EPOCHS=500n”, “n”, “import objaxn”, “import jax.numpy as jnpn”, “import numpy as npn”, “from tqdm.auto import tqdmn”, “from torch.utils.data import DataLoadern”, “n”, “n”, “opt = objax.optimizer.Adam(model.vars())n”, “n”, “@objax.Jitn”, “@objax.Function.with_vars(model.vars())n”, “def loss(x, y):n”, ” yhat = model(x)n”, ” return ((yhat-y)**2).mean()n”, “n”, “grad_and_val = objax.GradValues(loss, model.vars())n”, “n”, “@objax.Jitn”, “@objax.Function.with_vars(model.vars()+opt.vars())n”, “def train_op(x, y, lr):n”, ” g, v = grad_and_val(x, y)n”, ” opt(lr=lr, grads=g)n”, ” return vn”, “n”, “trainloader = DataLoader(trainset,batch_size=BS,shuffle=True)n”, “testloader = DataLoader(testset,batch_size=BS,shuffle=True)”
]
}, {
“cell_type”: “code”, “execution_count”: 5, “metadata”: {}, “outputs”: [
- {
- “data”: {
- “application/vnd.jupyter.widget-view+json”: {
“model_id”: “d8674d2ad786477f88ae916672683edb”, “version_major”: 2, “version_minor”: 0
}, “text/plain”: [
” 0%| | 0/500 [00:00<?, ?it/s]”
]
}, “metadata”: {}, “output_type”: “display_data”
}
], “source”: [
“test_losses = []n”, “train_losses = []n”, “for epoch in tqdm(range(NUM_EPOCHS)):n”, ” train_losses.append(np.mean([train_op(jnp.array(x),jnp.array(y),lr) for (x,y) in trainloader]))n”, ” if not epoch%10:n”, ” test_losses.append(np.mean([loss(jnp.array(x),jnp.array(y)) for (x,y) in testloader]))”
]
}, {
“cell_type”: “code”, “execution_count”: 6, “metadata”: {}, “outputs”: [
- {
- “data”: {
“image/png”: “iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6a0lEQVR4nO3dd3hUVfrA8e+ZSe8VAkkggdADCZClCoIFUATWgosFQUVWfyqurn3X7lZdC+qq2HAtiCKgCIqCIiAIEnpJKIFAqCmk98z5/XEnIUACIZlkJjPv53nuMzN35t77nijvnDn3FKW1RgghhPMz2TsAIYQQLUMSvhBCuAhJ+EII4SIk4QshhIuQhC+EEC7Czd4BnEtYWJiOiYmxdxhCCNFqJCcnZ2mtw+t6z6ETfkxMDBs2bLB3GEII0WoopdLre0+adIQQwkVIwhdCCBchCV8IIVyEQ7bhK6XGAePi4uLsHYoQohlUVFSQkZFBaWmpvUNptby8vIiKisLd3b3BxyhHnksnKSlJy01bIZzP/v378ff3JzQ0FKWUvcNpdbTWZGdnU1BQQGxs7GnvKaWStdZJdR0nTTpCiBZXWloqyb4JlFKEhoZe8C8kSfhCCLuQZN80jfn7OWXCL1j5X/T2+fYOQwghHIpTJvwTP89i/YLXWLbzuL1DEUI4oOzsbBITE0lMTCQiIoLIyMia1+Xl5ec8dsOGDcyYMeOCrhcTE0NWVlZTQrYJh+yl0xQWi8atXW86Hf6FIZ8k8829w+gW4W/vsIQQDiQ0NJTNmzcD8PTTT+Pn58eDDz5Y835lZSVubnWnx6SkJJKS6rwn6vCcroZvMik69kgiXGcTYipizvqD9g5JCNEKTJ06lQceeICRI0fyyCOPsH79eoYMGULfvn0ZMmQIqampAKxYsYKrrroKML4sbrvtNkaMGEGnTp2YOXPmea/z0ksvER8fT3x8PK+88goARUVFjB07loSEBOLj45k7dy4Ajz76KD179qRPnz6nfSE1lkPW8JvcD79tLwDGtjnJb+knbReYEMLmnlm0g51H8m16zp7tA3hqXK8LPm737t0sW7YMs9lMfn4+K1euxM3NjWXLlvH444/z5ZdfnnVMSkoKP/30EwUFBXTr1o277rqr3r7xycnJfPDBB6xbtw6tNQMHDuTiiy8mLS2N9u3bs3jxYgDy8vLIyclhwYIFpKSkoJQiNzf3gstzJoes4WutF2mtpwcGBjbuBG2M/9AXBRxjx5E8issrbRidEMJZTZw4EbPZDBhJd+LEicTHx3P//fezY8eOOo8ZO3Ysnp6ehIWF0aZNG44fr//e4erVq7n66qvx9fXFz8+Pa665hlWrVtG7d2+WLVvGI488wqpVqwgMDCQgIAAvLy+mTZvG/Pnz8fHxaXL5HLKG32T+EeAdQqwlHYvuz4n8MmLCnLOoQrR2jamJNxdfX9+a50888QQjR45kwYIFHDhwgBEjRtR5jKenZ81zs9lMZWX9Fcz6Brp27dqV5ORklixZwmOPPcaoUaN48sknWb9+PcuXL+ezzz7j9ddf58cff2xcwawcsobfZEpB216EFOwBIL+0ws4BCSFam7y8PCIjIwGYPXu2Tc45fPhwFi5cSHFxMUVFRSxYsIBhw4Zx5MgRfHx8uPnmm3nwwQfZuHEjhYWF5OXlceWVV/LKK6/U3GRuCuet9rbthV/G/1BYyC+RJh0hxIV5+OGHmTJlCi+99BKXXHKJTc7Zr18/pk6dyoABAwCYNm0affv2ZenSpTz00EOYTCbc3d158803KSgoYMKECZSWlqK15uWXX27y9Z13Lp3kD2HRDIaXvcwjN1zB2D7tbBucEKLRdu3aRY8ePewdRqtX19/RNefSsfbU6a4OSpOOEELgzAk/vDsaRXd1iLwSSfhCCOG8Cd/TD4Jj6GE+RL4kfCGEcOKED6i2vehuOiRNOkIIgZMnfNr2ogPHKC4qsnckQghhdw6Z8JVS45RSs/Ly8pp2ora9MGPBv2CvbQITQohWzCETfpOnVqhmnWIhrGiPDaISQjiLpkyPDMYEamvWrKnzvdmzZ3PPPffYOmSbcN6BVwAhsZQrTyJK0+wdiRDCgZxveuTzWbFiBX5+fgwZMqSZImweDlnDtxmTmRNesXSo3G/vSIQQDi45OZmLL76Y/v37M3r0aI4ePQrAzJkza6YonjRpEgcOHOCtt97i5ZdfJjExkVWrVtV7zvT0dC699FL69OnDpZdeysGDxnTtX3zxBfHx8SQkJDB8+HAAduzYwYABA0hMTKRPnz7s2WP7lgnnruED2X5xdC5ehdZa1tAUwhF9+ygc22bbc0b0hiv+2eCPa6259957+eqrrwgPD2fu3Ln85S9/4f333+ef//wn+/fvx9PTk9zcXIKCgrjzzjsb9Kvgnnvu4ZZbbmHKlCm8//77zJgxg4ULF/Lss8+ydOlSIiMja6Y9fuutt7jvvvu46aabKC8vp6qqqil/gTo5dw0fKAjoRpjKoyz3mL1DEUI4qLKyMrZv387ll19OYmIizz//PBkZGQD06dOHm266iY8//rjeVbDqs3btWm688UYAJk+ezOrVqwEYOnQoU6dO5Z133qlJ7IMHD+bvf/87//rXv0hPT8fb29uGJTQ4fQ2/NKQb7IOSjK14Bct8OkI4nAuoiTcXrTW9evVi7dq1Z723ePFiVq5cyddff81zzz1X77z4DVHdyvDWW2+xbt06Fi9eTGJiIps3b+bGG29k4MCBLF68mNGjR/Puu+/abNK2ak5fw68MN3rqVB7dbudIhBCOytPTk8zMzJqEX1FRwY4dO7BYLBw6dIiRI0fy73//m9zcXAoLC/H396egoOC85x0yZAifffYZAJ988gkXXXQRAPv27WPgwIE8++yzhIWFcejQIdLS0ujUqRMzZsxg/PjxbN261ebldPqE7x3UlkwdiDqx096hCCEclMlkYt68eTzyyCMkJCSQmJjImjVrqKqq4uabb6Z379707duX+++/n6CgIMaNG8eCBQvOe9N25syZfPDBB/Tp04ePPvqIV199FYCHHnqI3r17Ex8fz/Dhw0lISGDu3LnEx8eTmJhISkoKt9xyi83L6bzTI1ttOniSgnfGkRimCbjvFxtFJoRoCpke2TZkeuQzBHq7k6qj8c3bAxbb3/UWQojWwukTfoA14ZstZZAjA7CEEK7L+RO+lzu7LB2MF8flxq0QjsKRm5Nbg8b8/Zw+4Xu4mTjs1gELJjguN26FcAReXl5kZ2dL0m8krTXZ2dl4eXld0HFO3w8fwNPbh0wVRdvjje8/K4SwnaioKDIyMsjMzLR3KK2Wl5cXUVFRF3RMiyV8pVQn4C9AoNb6upa6LhjNOhlVsbQ9IQlfCEfg7u5ObGysvcNwOU1q0lFKva+UOqGU2n7G/jFKqVSl1F6l1KMAWus0rfXtTbleYwV6u7PPFAMnD0DZ+QdLCCGEM2pqG/5sYEztHUopM/AGcAXQE7hBKdWziddpkuqeOgCcSLFnKEIIYTdNSvha65VAzhm7BwB7rTX6cuAzYEJDz6mUmq6U2qCU2mCr9r1Ab3e2VVrbuqSnjhDCRTVHL51I4FCt1xlApFIqVCn1FtBXKfVYfQdrrWdprZO01knh4eE2CSjMz4Pfcv2ocvMFmWJBCOGimiPh1zXpvNZaZ2ut79Rad9Za/6MZrluvGwd2RGNiY3kkZbu+g9L8lry8EEI4hOZI+BlAdK3XUcCRCzmBzRYxt4oN8+XTaQP5T+X1mPMzyPl4KlgsNjm3EEK0Fs2R8H8DuiilYpVSHsAk4OsLOYHNFjGvZUhcGDf+4SaerZxMSMZyWPF3m51bCCFag6Z2y5wDrAW6KaUylFK3a60rgXuApcAu4HOttUN0gB+f0J7SxNuYZxkJK1+AHQvtHZIQQrSYJg280lrfUM/+JcCSxp5XKTUOGBcXF9fYU9Tr/lHduGHvXcRXHaP7wrsgNA4i4m1+HSGEcDQOOZdOczTpVGsX6M2U4V2ZXHgP5e7+8NkNUJRt8+sIIYSjcciE39yuT4pG+7XlTzyILjgOX0yBqgp7hyWEEM3KJRO+r6cbL12fyE8FHXgn6D44sAq+nAbpa6Cq0t7hCSFEs3DI2TKbsw2/2vCu4Tw4uhvPfVNFcNhNXJvyOaadC8E7GLqMhm5XQNyl4OnfbDEIIURLcvo1bc9Fa83H6w7yxo97KczP4fcBqdwQuJ0eBb9iKj0JZg+I7A8BkeAfAf7tTj2GdIKAds0WmxBCNMa51rR16YRfrbSiird+3sec9Qc5nl9GW183/tDuML3y19ClMpUQSw7+5ZnGMok1lPErYPA90HEIqLoGGAshRMuShN9AFotmc0Yury3fw8GcYtzNJk4Wl3M8vwzQBFBMvH8RCUEl9K7awcV5i/CtyiPVFEdOwnT6X3ErHh4eLRavEEKcqdUl/Fpt+Hfs2bPH3uFQWlHFzqP5/LY/h5RjBew9UUh2YRnFxQWMtfzM7eYldDIdI8vcBs+L/4T/sP+TGr8Qwi5aXcKv1tI1/AtVVFbJ/qwiYkO9WbHoI9pse5vfmVLREz9E9fq9vcMTQrigcyV8l+yWaSu+nm7ERwbi6+XB2Im3kzJmDnst7Sn+/m8yOZsQwuFIwreh6wfE8IHb9fjm7UbvXGjvcIQQ4jSS8G3I081Mr1FT2GOJJP+758FSZe+QhBCihkMmfFvPh9+SbhgYy1dBNxNYuI9l897mcG6JvUMSQgjAQRN+c06e1tyUUky48W7SzR2I2f4aN8/6BYvFcW+MCyFch0Mm/NauS0QgHa95ljjTEXrn/siG9JP2DkkIISThN5seE7CE9+R+jwW8tizF3tEIIYQk/GZjMmEa+SixHCF0/yL2ZRbaOyIhhIuThN+cuo+jPKwXM9zm8+OOw/aORgjh4hwy4bfmXjqnMZnwuPRxOpmOobbPs3c0QggX55AJvzX30jlL97Ec8uzC6OyPZHEVIYRdOWTCdypKsS1uOtH6KMW7lto7GiGEC5OE3wJ8eo6hVLuTs+17e4cihHBhkvBbwKCukWwx9UCnrbB3KEIIFyYJvwV4uZspiRpGdMUB8rOkt44Qwj4k4beQ8IRRAOxdt8TOkQghXJUk/BbSLWEo+fiSvXUpjrzojBDCeTlkwneafvi1uLm7c7LNIHqUbmLP8QJ7hyOEcEEOmfCdqh9+LT7dLyFKZbF71zZ7hyKEcEEOmfCdVVhvox2/fM+Pdo5ECOGKJOG3IBXWhRxzOKEn1to7FCGEC5KE35KU4njoQHpXbCGvuMze0QghXIwk/BZmjhtBiCpk71ap5QshWpYk/BYW2W8MAEUpy+0ciRDC1UjCb2G+YdGkm6IJOrbG3qEIIVyMJHw7OBw8gC6l29CV0o4vhGg5kvDtoCTqIrwpI3eP1PKFEC1HEr4deHe9mCqtKNy5zN6hCCFciEMmfGecWqG22Kj2bNOd8Dy02t6hCCFciEMmfGedWqFaRIAXv5n6EJq7DUrz7R2OEMJFOGTCd3ZKKXLaDMZMFaRLO74QomVIwreTgK5DKdXulO2WeXWEEC1DEr6d9OrQlt8s3aja95O9QxFCuAhJ+HbSLcKfNZZ4fHJ3Q+EJe4cjhHABkvDtpI2/J1s8Eo0Xu7+zayxCCNcgCd9OlFKo9omkmWJg7Rtgsdg7JCGEk5OEb0eX9IhgZumVkJkCe763dzhCCCcnCd+OLunehm8sgyjyage/vGLvcIQQTk4Svh11CPHBZPZgTdtJcHAtHFxn75CEEE5MEr4dmU2K6BBvFpkuBa8gWDPT3iEJIZyYJHw7iw3zJTVHw4DpkLIYMnfbOyQhhJOShG9nsWG+HMguoiJpGrh5Si1fCNFsJOHbWd8OwZRVWtie6wF9b4atcyH/qL3DEkI4IUn4dva7mBAAfjuQA4PvBkslrHvTzlEJIZxRiyV8pZSvUupDpdQ7SqmbWuq6ji7c35PYMF/W7z8JIZ2g5wTY8AGUOudaAEII+2lSwldKva+UOqGU2n7G/jFKqVSl1F6l1KPW3dcA87TWdwDjm3JdZ/O7mGB+O5CDxaJh6H1Qlg/Js+0dlhDCyTS1hj8bGFN7h1LKDLwBXAH0BG5QSvUEooBD1o9VNfG6TmVAbCh5JRVsSD8J7ftC7MWw9r8gi5wLIWyoSQlfa70SyDlj9wBgr9Y6TWtdDnwGTAAyMJL+Oa+rlJqulNqglNqQmZnZlPBajSviI2gb4MmL36caO4beB4XHYPmzMseOEMJmmqMNP5JTNXkwEn0kMB+4Vin1JrCovoO11rO01kla66Tw8PBmCM/x+Hq6MXVILOv355CeXQSdL4F+U2Dt6zD3JlkGUQhhE82R8FUd+7TWukhrfavW+i6t9SfNcN1W7fd926MULNh0GJSCca/CFS/A7qXw7mWQtdfeIQohWrnmSPgZQHSt11HAkQs5gVJqnFJqVl6e6/RUaRfozZDOoSzYdBittZH0B06HW76C4ix45xLY84O9wxRCtGLNkfB/A7oopWKVUh7AJODrCzmB1nqR1np6YGBgM4TnuK7uG0V6djEbD548tTN2GExfAcEd4JOJsOol0NpuMQohWq+mdsucA6wFuimlMpRSt2utK4F7gKXALuBzrfWOpofq/MbER+DlbmL+xsOnvxHUAW77HnpdDcufgc9vkX76QogL1tReOjdordtprd211lFa6/es+5dorbtqrTtrrf92oed1xSYdAD9PN8b0iuCbrUcpqzyj56qHD1z3Plz+nDHJ2qwRcHSrXeIUQrRODjm1gqs26QBc1z+avJIKPlt/6Ow3lYKhM2DqYqgoMW7mJn8oTTxCiAZxyITvyobGhXJRXBgvfp9KVmE9A686DoY/rjIeF82ABXdCeVHLBiqEaHUk4TsYpRRPj+9FaUUV//w2pf4P+oXDzfNhxGPGDJvvXCJz6QshzskhE76rtuFXi2vjx+0XdWJecgY/phyv/4MmM4x4FCYvgKIseO9yyEhuuUCFEK2KQyZ8V27Dr/any7rQPcKfB7/YyrG80nN/uPNImLYMvALhf+Mh7eeWCVII0ao4ZMIX4OVu5vUb+1FSXsWf5m6iynKeG7MhsXDbUgiMNvrrpyxumUCFEK2GJHwHFtfGj2cn9OLXtBymvL+ekvLzTDIa0A5uXQIR8TB3Mmz5rGUCFUK0Cg6Z8F29Db+26/pH8dDobqzem8WTX21nXVo2x/JK66/x+4QY0zHEDIUFf4R1s1o2YCGEw1LagftwJyUl6Q0bNtg7DIfw4BdbmJecUfO6d2Qgn9wxkAAv97oPqCiFebdB6mJjuuURj4G7dwtFK4SwF6VUstY6qc73JOG3HtmFZWzJyGXviUL+viSFh0Z34+6RcfUfUFUJix+AjR9CYAcY/Tz0GG8M4BJCOKVzJXyHbNIRdQv18+SS7m2ZPrwzCdFBLN1x7NwHmN1g/ExjZK6nvzEHz/8mwIldLROwEMKhSMJvpX6f2J6tGXm8/fO+83845iL440pjfv2jm+HNofDdY1CSe/5jtTa6eX50DVUvdGX6v95jf5aM6m1Ji7ceJfVYgb3DEE7AIRO+3LQ9v1sGxzC2dzv+8W0K320/T00fjNr+wOlw7yboNxl+fRP+0x3m3ACbPjYGbtVWVQnb5xuTtP1vPBzbRmml5t/FT/DNd0uapUzibNmFZdz96UZu/WC9vUMRTsDN3gHURWu9CFiUlJR0h71jcVRmk+KVSYmkZRXx8LwtBPm40zsyEF/P8/wn9Q01VtNKug02fgSpS4xNmSB6IHS7Ety84Nc34OQBCOlsfL7PJH7dsI1u303itrT74UgXaJ/YEkV1aUt3GCOtvdzNdo5EOAOHrOGLhnE3m3jnlv60CfBi0qxf6fXUUtalZTfs4HYJMPZFuH8HTP8Zhj8EZYXwwxPw7UPgEwZ/+Bju+Q36TwV3Lwq92zOp/AmK8DbuBRzZ3JzFE0B6jtF81j5IeliJppOE38pFBfvw5Z1DGNYlDIBP1x+8sBMoZdTURz4Od62G+7Ya7f3TlkGPccZ8PVYl5VVk6HBmeD4HngGS9FtAeaUFgILSCjtHIpyBJHwnEOjjzke3D+SmgR34avMR3vipCQueB3c0av91dN0sso703VESDFO/kaTfAspqEn6lnSMRzkASvhN5clxPxie054Wlqdw/d/P559+5ACcKSim0Jp2C0kosgR1OT/qHfrPZtcQp1TX8fEn4wgYcMuFLL53G8XQz85/rE5hxSRwLNh3m1eV7bHLe/NIKBvxtOS8vOzXf/snicuPXwNRvwDsYPhwHe36wyfXEKdUJv7BMmnRE0zlkwpfpkRvP3WzigVHd+H1ie2Yu38M/luw6/6Rr55FTWH7WvmP51imbgzvC7d9DWBeYM0kmbLOx6rWNSyssVBSdhPnTIXm2fYMSrZZDJnzRdC9MTGDyoI68vTKNHk9+x9dbjlzQ8eWVlppkk1tydu3yeH6tOfr92hijeTsOMSZsW/Nak2IXp1TX8CPJRH0wxljdbMU/wdK0L3HhmiThOyl3s4lnJ/Tivku7ADBjziYG/2M5R3JLGnT8pS+tIPEZo4nmZPGpGr6Xu/G/zG2zN/DUV9tPHeAVADfNg56/h+//Ct8/IYur20B5lYXeKo0Fnk9hyj8CA++EgqNwYJW9QxOtkCR8J6aU4v7LuzJrcn8AjuaVsqhWTX/74Tzun7u5ziafQzkllFRYa/i1En5MqG/N8yXbj3Ha5HtunnDd+/C7O2DNTFh4F5TJlABN0btwDXM9nqMcN/aOmw+XPQ0e/rD1C3uHJlohSfguYFSvCL66eyiRQd7MWpnGp+sOorXmvdX7WbDpMK//VP/N3e2H8zhZdKpJJ8DLnV8fu5QnrupJZkEZR89cftFkhitfgJF/gS1z4KWeRm0/73BzFc95rXubh3KfI40ori57luNeMcYU1z3Hw66voaJhv9aEqCYJ30UkRAfx1s398fNy4/EF2xj54goWbDKS8LzkDCy1unBmF5bVPL/qtdUcrtUM5O6miAj0IqljMACbD+WefTGl4OKH4Y6fIO4yWPsGvNoHvrwDjm5pngI6m++fgG8fZp3b73jE/+9kEkROUbnR1bb3RCjLh93f2TtK0cpIwnchvaMC+fHPI7isR1s83U6NoD2eX8abP+/j3jmbKCqrpP/zy0477pe9pyZW8/Uw5urpFuGP2aTYdTS//gtG9oOJH8CMTTBgujFnz9vDYfZVsHlOw2brdEUpS4wmsf638qTnIwQFBQHwwtJUBv1jOfv9+4NfhDTriAvmkJOnKaXGAePi4s6xuIdoFLNJ8e6UJKosmn8vTWFMrwgeX7CdF5amAtAxxOesY1KOFdApzJfnr44nOth438vdTFSwN6/9uJdO4b5c3Teq/osGd4Qx/4CLHzEWY1n/Diy8E0zu0Hkk9JxgTNrmE9IsZW5VSvOMRWva9IIr/k3JjtW08fcCIOOk8Uvrvz/v54Xe18G6t6E4R/5uosEcsoYv/fCbn9mkeOyKHvTtEMzjV3av2f96PdMy9I4KZEjnMKJrfSH4WGv7//l+d53HnMU7yFhu8b6tMG05DLoTMlPgq7vhxS7w0TVG88+RTcb0zK7oh6eg8Djr+zxDzF9/4EheyVkzZWYVlhnNOpYK2LnQPnGKVskha/iiZQ3tHMbvE9tTXF7F9zuP1+yfOiSGI7klfL/zeE2bfW0vTuzD2JmrKa2woLVGNXTpRJMJopKM7fLnjAS/8yvYtQiWPm58xsMfOgw0+vZ3HAqRScac/s7swC+Q/AEMvoe/rPcAytEaPN1O1csig7yNaRba/Q7CuhnNOkm32S9m0ao4+b8g0RAmk+KVSX0pLq/k+rfXkldSwaNjenBl7wgyC8rw93JnfGLkWcf1ah/IM+N78dTXO9ifVUSncL8Lv7hSRlt/ZD+4/BnIPwLpa05ty581PuffDvrdAv2mQODZsbR6FaWwaAYEdYSRj7NvxYqatzxqJfye7QNIzy4y/m59JsKPz0PuQQjqYIegRWvjkE06wj58PNz45t5hrHxoJGP7tEMpRZsAL/5zfQKB3u51HjO6VwQebib+u6IBSy02REB76H0dXPUS3P0rPJQGE2dD23j4+d/wSjzMuRH2LAOLxTbXdAQ//wuy9xqLzXj4UnveOw+zic7hxviHIG938kuszV29JxqP2+TmrWgYqeGLszS4aQaICPRi8qCOfPDLfu4a0ZnOjanln4tvKPS62thOHjDmkdn4EaQuNmrDPcZBZH/jF0JQx7Onda6qhOPb4dB6OLTOGCfQb4rRVHQB5WxWR7fCL69C4s3QeeTpg9kwmnQWzxhGpUXz8g+7OVlcTnmlBY/gGIgeZDTrXPSA45RHOCxJ+KLJ7hrRmTnrD/LKsj28dkPf5rtQcIwx0nTE45CyyEj+62dBlXUksE8otLc2D1mqjAR/OBkqio33/dtDRZExH03beBhwh1FL9vCt54ItoKoSvr7HiH3Uc8DZUyF7uJlqbtz6e7lRVmmh61+/5cA/xxrNOov/bHypRfRu8fBF6yIJXzRZmJ8ntw6N4Y2f9nFtv0hGdGvTvBd084D4a42tshxO7IDDG+HIRuNx33JAQbs+Rrt/9ABjvd7AKCgvNppA1r8Di+6DH56EvpMh8UYI6wrmupuums2v/zUGo02cDT4h7D1RwGUvrTztI7Xb8AO8zoiv59Xw7SPGl5gkfHEekvCFTUwf1pmP1qZz6+zfmHfnYPp3bKG+4W4e0L6vsXG7sa+8CFDgcfaYAjx8oP8U44vg4K/GL4R1b8Ha141xASGdILyr0QMmvBu0SzReN4eUJcZN125jjUnnMMY8nKn2IDl/rzP+yfqGQtzlsO1LuOyZ05akFOJMkvCFTQT6uPPVPRcx+pWVzFqZxtuT7TgYqCFNNEpBx8HGln8U9v8MmamQtRtOpBjJWFsnlYu/Fi59yhhAZiu/vQdLHjS+UMa/VtP+7mE+ux9F9QylAKUVdUyL3Gci7P4W0n+B2OG2i1E4HUn4wmZiw3y58+LOzFy+h3Vp2QzsFGrvkBomoB0kTDp9X2U55KTBjvnwy0zY9Q0MuguGPQBeTRgQqLVRq1/1InQZbUw9UesLyr2OhB/q51nzvGOt2Uorqyy4mU3Q9Qpj3ELybEn44pwcslumLHHYev1xeCc6hvpw1ycb2X28FU+N7OYBbbrDyMfh3mSjlv/LKzCzr9H+Xz0SuLIccg8Za/ru/NroQXRiV91rAVRVwML/M5J9vykw6dOzfo1UVJ3d1TTMz6Pm+fCu4UweZPzSKLUujoKHDwyYBtu/NKZbEKIeDlnD11ovAhYlJSXdYe9YxIXx9XRj9q0DuPbNNYyduYrBncPIKijjf7cPIKxWTbVVCYyEq9+EgX80FndZ8iCsfMHoCVScVc8x0RB3qdG+HjvcaLKZOxnSfjKmjh7+0GndKMsrLVRaLJRV1pXwT/+7dW1rdH0tKa/Cz9P6T/iSJ4wmqe8eNa7d/UrblF04FYdM+KJ1iw3zZcmMYfz5i82s3J0JwMzle3h2QrydI2ui9okwZZEx6+fWz43F2/3bgX/EqUd3bziwGvYuM26kJs8Gk5vR7bIoCya8AX1vPuvUN7zzK8npJ3lxYsJZ74X4epz2urqL5mnt+SYzXPsuzB4L826DWxcb4xOEqEUSvmgWEYFefHjrABZvO8rCTYf539p0Lunepvm7bDY3paD7WGOrT1gXSLrVaO45tM5I/se2wqC7octldR6SnH4SqPum7Jnt+t4edSR8MJqHbvwc3r0UPp0E05bZ9kazaPUcsg1fOAc3s4kJiZH8eVQ3TAqmfvAbzyzaQV7x2YuiOyU3D4gdZswRNHlBvcm+trqadM7kZe2mWVJXjx2/NsbawlVl8MlEKDl5wWEL5yUJXzS7+MhA1j1+GVMGd+SDXw4w4O/LWJF6wt5hOaSislOjbAd3CiUh6uweQdU1/LrWIgaM8QOTPjV6Gc2dDJVldX9OuBxp0hEtItzfk2cmxHNd/2ge+Hwzd3+ykX9c24fxCe3tHZpDeemHU2sLzJk+qM7PVLfhf/DLgfq7vsZcBL//L8y/A94bZYwi9gkx7jt4Wx+9AsDNy7jv4OZlfe4FngHG2gXC6UjCFy2qd1Qgb0/uzwOfb+G+zzYxZ91Bru4byYmCUo7mlfK3q43pAS5ofn0ndOfFnet9z9ua8L/bcYyMk8VEBdcxohigz/XG2rcb/2fcSyg5abxuCJ9QCO0CYXHGl0VoF2MuI7M7KJNxL0OZjM3saTQlufB/r9ZCEr5ocZ3C/fhk2kBm/riHhZsO8/CXW2veaxfoxY4j+RzLL+XVP/TF28NMuH/ju3NWVlnQ1D2gyZE9ekX3et+rbtKBczTrVPvdNGOrVlVhLKNYnGMk/4oSo8mnssSYk7+yxPhiyN4LWXth91LY9PH5A/aLMEYtd7BubXvJNA8OSBK+sAtfTzceu6IHMy7pwqfrDvLp+oPszyrixVrLJQ5/4Sf8Pd3Y9szoRl/nypmrOJBdzO7nr7BF2A6h9lQLXyRn8NDobg3/QjO7g2+YsTVUyUkj+ecdtK5BoEFbTm1lhXB4A6SvhR0LjGM8AyDqd0YvIb+2xi8Av7bG5htu9Cgye4Cbp/Eovw5ahCR8YVe+nm7cMbwTdwzvxDdbj/DYl9sY2CmEZbuMm7oFZZXMWrmPSoumolJz32VdLuj8u48XNkfYNnXrB+sv6PPetda4nbUyjUBvd+4eGWfrsGpdMBiif2ds55N70Ej8B9dCxgY4uhmKs89/nMndSP6+YcbqXYEdICjaGEQWFG3cd3DztG5e1i8K62v5smgwSfjCYVzVpz1X9TFu4n7wy36eWbQTgL8vSan5zIUmfEdXVFbJT6mZF3RM9eLx1Y7kltgypKYJ6mBsCX84ta+qAooyofA4FJ4wHitKjHUMKstOPVaWGe/lHTLGLhQeO//13H2N6wV3tF7b+hgYZQyE8w1v/JTXliqj2cv77PWcWytJ+MIhTR0Sw8DYUK6cueq0/SXlVae1YTdUY49rbhsPXng/eQ83Ewv+bwhX/3cNAG4mB6/hmt2NpSsDLrBHVmUZ5GUYXwCl+dYvhlLrVmZ8aRRnw8l0yE03FoEvr2P+Jp/Q05uVfMOtzVrhpzaTG2TvgczdkJVqPGbvNcYzBERa111OgqgkY4ZTzwau7FZVYQy6O7wRSnKNBXgqSowpvCtKjK37WOh704X9bRpJEr5wSEoperYPYPfzV/CnuZtYss2o7Y165WemD+9cM4FYQ2UXlRFV1/z4dpZytHETzHWPCKh5bnL0hN9Ybp4Q2tnYGkJr435D7kHIP2z9NXHC+KVQ/cvi4FooyjYSb52U8WshrBvEXQI+YXBsm3GPYtci60dMxvshscYvicDoU4++ocbkeYfWGRPqHdlk3AivZnIHdx9jwjt3b+N5WctNMigJXzg0DzcTb9zYj02Hcrnmv2s4lFPCEwu3Nyjh114bNqeovP7ui3ZU52jZBqh94/a8PXVchVLGWAOfEGPeo3MpLzKamYqyjC+DqjIIjTM2d++6jynKNpbMPLzBWKXsZLoxb1JdXV1N7tAuwZhiI3qAcQPbr23Lr6h2hhZL+EqpTsBfgECt9XUtdV3R+iml6BsdxB+Sopm74RAAV766ilcnJdKlrX+9x9WepuDF73fz1s39+DUtm77RwQSfMSFZZkEZryzbzdje7RgSdwE9WJqorPLsZD11SMx5j6s9RuGz3w5x98g4okMc7wvNYXn4GltwTMOP8Q2FrqOMrbbSPGvTU4bxKyKsq9Hs4+5ly4htokF9uZRS7yulTiiltp+xf4xSKlUptVcp9ei5zqG1TtNa396UYIXrUkrxr+v6sOWpUSREB7HzaD6Xv7ySn3fXf8OzoNZi4Ct3Z9LzyaXcNnsDMz7bdNZnV6Se4JN1B/nT3M3NEX69SivOnjvn6fG9Lvg8d/xvgy3CEY3hFWiMO+g62lg6s8Mgh0z20PC5dGYDY2rvUEqZgTeAK4CewA1KqZ5Kqd5KqW/O2Fr5FInCUQR6u/PlnYOZMtho0pny/npmzNnEwezisz5bUFr3JG37TpzeVTOrsIxlu44Dp89l0xLqquE3VOfwU4unnCwut0U4wsk1KOFrrVcCOWfsHgDstdbcy4HPgAla621a66vO2GSmLGEzbmYTz0yIZ9vTo7h7ZGe+3nKEkf9ZwZfJGaQcy+eO/23gZFE506y13ofHdDvteE/303vr3DDrV5busCb88ioWbMqgpLyKmEcX84W1Cam5lNVRw2+o5X8eQdsAYxRy2wDHrFEKx9KU8eaRQO1/DRnWfXVSSoUqpd4C+iqlHjvH56YrpTYopTZkZl5Y/2ThWvy93HlodHe+vGswYX4e/PmLLYx5ZRU/7DxO3+d+IC3T6InRo13Aacd5up3+v/2eM2r898/dwoFs49hXlu1pxhLUWqawkT7/o1H2nCKp4Yvza0rCr6svWB0LeVrf0Dpba32n1rqz1vof5/jcLK11ktY6KTw8vAnhCVfRv2MIyx64mHtGxmFSEOrrcdo6sNFn9M4pKjeabbTWfLjmQJ3n3HXU6HlRexCnrmud2iYqO6OXzl/H9rig4zuG+nJt/yhO5Jc1S3zCuTSll04GEF3rdRRwpGnhGJRS44BxcXHNOFxcOBV/L3ceHN2NO0d0xtvdjElBcXkVeSUVtA/y5rs/DWNjei5PfrWdQzklTH5vHcXlVTUrTZ1p2+E8ADJOljBn/UHMJsXD87aS/NfLCLXh2ry1exI9NLob04Z1uuBzRAR4UV5l4VBOCR1CpaeOqF9Tavi/AV2UUrFKKQ9gEvC1LYLSWi/SWk8PDDx78QchzsXP0w2zSaGUwtfTjfZBRp/q7hEB3DiwA/df3hWAVXuyzkr2PdoF0KWNMYJyuzXhAzw2fxufrjsIQMox2wySyS+t4B9LdpFZcGpxksbW0IfGheHjYeaRWrOOClGXhnbLnAOsBboppTKUUrdrrSuBe4ClwC7gc631juYLVYimmzokpqaHT20J0UF8Nn0Q8+4aAnDWl4G/l/FjuLptv6l+Ts3k7ZVp7Dx6atCOpZEtMl3b+nP3yDjWpmVzKOfs3kpCVGtoL50btNbttNbuWusorfV71v1LtNZdre3yf2veUIVoOl9PN54e34v/3tQPgCt7RxDg5caz43sR6O1OgDWxn5l8V+3JAmD5Ltt0OKtrdGxVYzM+MLZ3OwBWnGNcghAOObWCtOGL5qSU4sre7djy1Ch8Pcy41ZpLXinFPSPjeP2nvXUe+2PKCeYlZ3Bd/6gmxVBXv3lLE266dgz1IdDbnVeX7eaiuDBiw3zPf5BwOQ65DJC04YuWEOjtflqyr/bg6G4c+OdYdj9/Bd/cexER1j7uvh5mYkJ9eOn7VE7klzbp2rklZw8Ka0oNXymFp5uJrMLyFh8tLFoPh0z4QjgCDzcT8ZGB/PTgCABuHtSRVyb1JbekghvfXUdhE0bl5tZRw69qYrfKsX2MZh0nnTtT2IBDJnyl1Dil1Ky8vLzzf1iIZubtYWbLU6N4eEx3EqODeOeWJPaeKOSx+dsobeRslyeLTtXw3c1Gik7qGNKkOB+7ogeJ0UFUVDVtMJdwXg7Zhq+1XgQsSkpKusPesQgBRvNPtaFxYUwZ3JEP16aTX1LBiG7h3DiwA55uZrZm5HIwp5h2gd6E+HoQE+pz2syW1XJLTtXwe0cGMuuWJMKa2L/fw81E9wj/muUhhTiTQyZ8IRzdMxPiKa2wMHfDIX7enUlaZhGPX9mD8a//ctrnXv5DAlf3Pf0Gb0l5FbnFFUSHeHMop4Ti8qomJ/tqbQK8yCos4w9vr2XuHwfb5JzCeThkk44QrcHfr+nN/P8bwjV9I/no13R6PPndWZ95b/V+ymuNps0qLKP300tJOVZAt7bGHD9NuVl7purlDtftzzntukKAJHwhGs1sUvTrEMx/rk/gvSlJ/OmyLrx0fcJpzT/bD+fz1NenlpHYmH6SSoume4Q/L/0hgT9e3Ik3b+5ns5gGdw6tef784p02O69wDsoRJ1yq1Q//jj17mne2QiGaS05ROfd9tolVe7IY26cdiVFBHMsv5cM1B9j+zGi83JtnUfUTBaUM+NtyAFY/MtIhl3YUzUcplay1TqrrPYes4Us/fOEMQnw9eOeWJKYOiWHZzuP8bcku3lu9nz5Rgc2W7AHa+HvxxZ1G+/2LS1NlFk1RQ27aCtGMvNzNPD2+F0+P78WiLUdIzy7imn5NG6XbEEkdg+kQ4sPCzUe4tEdbxiW0b/ZrCscnCV+IFtKSSVcpxeIZF3Hpf37mmUU7iQz2pl+H4Ba7vnBMDtmkIwOvhGg6fy933ry5P94eJia/u46Pf01v0hq6ovVzyIQvbfhC2Eb/jsHMu3MIkcHe/HXhdq57cy359SzuLpyfQyZ8IYTttA3w4tv7hvPqpES2Hc4j6bllPDB3MxsOSF99VyNt+EK4ALNJMSExknA/T97/ZT/zNx1m/qbDtA/0IikmhI6hPozuFUGv9gEopaisstQ5k6ho3RyyH361pKQkvWHDBnuHIYTT2Xwol1eX7Wb7kfzTllkc1CmEk0UV7M0s5JNpAxnUKfQcZxGO6Fz98CXhC+Hidh7Jp6Siim+3HWXVnizKqyzszyoiIsCLKUNiuHVoDFUWTWWVJq+kgn2ZhYzs3sbeYYt6SMIXQlyQDQdy+Pd3qaw/kINScGaamHFJHOVVmolJUXQO90NrTVF5FX6e0kpsb60u4cvUCkI4hl/Tslm9JwtfTzeO5JZwLL+UH3YeB0ApMClFpzBfsgrLKKmoYt6dQ4iPlN51ddFas/FgLmF+HnQMbb4lKFtdwq8mNXwhHE9BaQUWDXnFFby3Oo1Ve7NAQ1pWESYF/ToEM7J7GzILyrj9oliiQ2QuH4CtGbmMf/0Xgnzc2fzkqGa7zrkSvvz+EkJcEH8vYzbQQG93npkQX7P/UE4xc9Yf5MeUE7ywNBWA2WsOEB3izfThnRmf0P6sReNdyfF84+Z4bnEFpRVVPLNoJw9c3pVwf9ushdAQkvCFEDYRHeLDw2O689DobhzILmbnkXwem7+VQzklPLFwO09+tR2TUozsFo6vpxtpmUX0jgrk+QnxKAX/+X43A2JDGN413N5FaRb5tRauX7rjGHPWH6SiysKLExNaLAZJ+EIIm1JKERvmS2yYL5f3bAsYzRk/7DrOziP57DpaQFF5JW4mE9sO57HpYC6FZRUcyimBn+CTaQPp3zEYd7MJk6LOJSJbo7xaCb+swhjwVlRW2aIxSMIXQjQbDzej+SYpJoSkmNMXadda8/Gv6Xy95QjuZkVcuB97Mwu56d11uJsVZpMi2MeDfh2CSYgOZFTPCGLC6r7ZmVNUToivR7OXpzH2ZRZyPL+Un3dn1ux7Y8VewPgS+GHn8ZovxuYmN22FEA7jZFE532w9QurxAqoskHosn8O5JTXt39Eh3lydGEm4vycHc4oZ1SuCNXuzeXnZbr659yLiIwPZfbyA3OIKBsSGnOdq51dYVokCfBvR3XRrRi4/7DzOpoO5rN6bdc7PLv/zxXQO92tklKdrdb10pFumEKKaxaL5NS2bD9ceID27mNTjBWeNCwDwdjfz+NgePLHQWFIy5bkxTV5opvdTSzGbVaN61fR44jtKKqrwdjdTUnHuWUrn3DHotOUpm6LV9dLRWi8CFiUlJd1h71iEEPZlMimGxIUxJC4MMJpBcovL8XAzseNwPkfySli9J4uVezJrkj3ANf9dQ1wbPwZ1CmVoXCjeHmY2HDiJSSlW7cnk1qGxxLU5u1adXViGu5sJf083Cqxt7FrrmnsJ1765hq5t/bhtaCxd2vpjsWgsWp/V+6g6yZ8v2YOxuH1LcMiEL4QQ9Qn0dq9ZKL5doDcAtwyOoayyit3HCjmcW8z6/SdJOZbPr2nZfL3lSJ3nmfvbIV6d1Je+HYJo4+9JUXkVj8/fxuJtR4kJ9eHaWiuTPfvNToZ3CSchOojk9JMkp59kzvpD7P/Hldw7ZxPr9mez4a+XN7pMLZXwHbJJp5q04QshmsJi0aTnFLN813HySyoY0b0NS7cfY15yBkXllZRWXNj00BP7R/FFckbN66lDYpi95gAAW54cRaCPe817nR9fQpXl3Pn12/uGcdVrq7nz4k4Ul1fh7+lGuyBvxiW0b/Q0Fa2uSUcIIWzBZDK6iE4b1qlmX78OwTx2ZQ8qqyysTcvmYE4xy3Yep6C0kv4dgxncOZRVe7L4KeUE8ZGB7DlRyK6j+QCnJXugJtkDbMnI5aK4MP69NJVhXcJqFo8P9nHnZLHRJfOZ8b146usdNcf0aBdAqK8HP6VkstN6DTC6a9aO2VYk4QshXJKb2cSwLsYgr5sGdjztvRHd2vDEVT1rXq/fn8PyXcfJKixnbJ8Ibpu9gd/FBPP7vpFc1qMtA/++nFveX8+QzqGs2ZfN5xsOUV25nzasEy8sTeXJq3pyw4AONQnfZB1ecKKgjBMFpzfppB4raJYyS5OOEEI0gdaavy3exbur9wPU3Oz18TCz4qERtPH3qqntK6XYn1WEWSl8Pc2E+nnyn+9Tee1Ho1/+b3+5jPvnbqagtIKv7rmoUfFIk44QQjQTpRR/vaon91/elaU7jtG3QzAf/LKfa/tF0cbfq+Yz1WLPGDz251HdmDyoIynHCgj396RLWz8+W38Ii0VjMtl2lLEkfCGEsAFfTzeusfbsebbWpHIN0SbAizYBxpfDoE6h5BZXUFReWTNRna1IwhdCCAcyulcEo3tFNMu5XXOeUiGEcEEOmfCVUuOUUrPy8vLsHYoQQjgNh0z4WutFWuvpgYGyVJoQQtiKQyZ8IYQQticJXwghXIQkfCGEcBGS8IUQwkVIwhdCCBfh0HPpKKUygfRGHh4GnHtdMecjZXYNUmbX0Ngyd9Rah9f1hkMn/KZQSm2obwIhZyVldg1SZtfQHGWWJh0hhHARkvCFEMJFOHPCn2XvAOxAyuwapMyuweZldto2fCGEEKdz5hq+EEKIWiThCyGEi3C6hK+UGqOUSlVK7VVKPWrveGxFKfW+UuqEUmp7rX0hSqkflFJ7rI/Btd57zPo3SFVKjbZP1E2jlIpWSv2klNqllNqhlLrPut9py62U8lJKrVdKbbGW+RnrfqctczWllFkptUkp9Y31tVOXWSl1QCm1TSm1WSm1wbqvecustXaaDTAD+4BOgAewBehp77hsVLbhQD9ge619/wYetT5/FPiX9XlPa9k9gVjr38Rs7zI0osztgH7W5/7AbmvZnLbcgAL8rM/dgXXAIGcuc62yPwB8Cnxjfe3UZQYOAGFn7GvWMjtbDX8AsFdrnaa1Lgc+AybYOSab0FqvBHLO2D0B+ND6/EPg97X2f6a1LtNa7wf2YvxtWhWt9VGt9Ubr8wJgFxCJE5dbGwqtL92tm8aJywyglIoCxgLv1trt1GWuR7OW2dkSfiRwqNbrDOs+Z9VWa30UjOQItLHud7q/g1IqBuiLUeN16nJbmzY2AyeAH7TWTl9m4BXgYcBSa5+zl1kD3yulkpVS0637mrXMzraIuapjnyv2O3Wqv4NSyg/4EviT1jpfqbqKZ3y0jn2trtxa6yogUSkVBCxQSsWf4+OtvsxKqauAE1rrZKXUiIYcUse+VlVmq6Fa6yNKqTbAD0qplHN81iZldrYafgYQXet1FHDETrG0hONKqXYA1scT1v1O83dQSrljJPtPtNbzrbudvtwAWutcYAUwBucu81BgvFLqAEYz7CVKqY9x7jKjtT5ifTwBLMBoomnWMjtbwv8N6KKUilVKeQCTgK/tHFNz+hqYYn0+Bfiq1v5JSilPpVQs0AVYb4f4mkQZVfn3gF1a65dqveW05VZKhVtr9iilvIHLgBScuMxa68e01lFa6xiMf7M/aq1vxonLrJTyVUr5Vz8HRgHbae4y2/tOdTPc+b4SozfHPuAv9o7HhuWaAxwFKjC+7W8HQoHlwB7rY0itz//F+jdIBa6wd/yNLPNFGD9btwKbrduVzlxuoA+wyVrm7cCT1v1OW+Yzyj+CU710nLbMGD0Jt1i3HdW5qrnLLFMrCCGEi3C2Jh0hhBD1kIQvhBAuQhK+EEK4CEn4QgjhIiThCyGEi5CEL4QQLkISvhBCuIj/Bxk5SdAnc6urAAAAAElFTkSuQmCCn”, “text/plain”: [
“<Figure size 432x288 with 1 Axes>”
]
}, “metadata”: {
“needs_background”: “light”
}, “output_type”: “display_data”
}
], “source”: [
“import matplotlib.pyplot as pltn”, “plt.plot(np.arange(NUM_EPOCHS),train_losses,label=’Train loss’)n”, “plt.plot(np.arange(0,NUM_EPOCHS,10),test_losses,label=’Test loss’)n”, “plt.legend()n”, “plt.yscale(‘log’)”
]
}, {
“cell_type”: “code”, “execution_count”: 7, “metadata”: {}, “outputs”: [], “source”: [
“from jax import vmapn”, “def rel_err(a,b):n”, ” return jnp.sqrt(((a-b)**2).mean())/(jnp.sqrt((a**2).mean())+jnp.sqrt((b**2).mean()))#n”, “n”, “rin,rout = trainset.rep_in(G),trainset.rep_out(G)n”, “n”, “def equivariance_err(mb):n”, ” x,y = mbn”, ” x,y= jnp.array(x),jnp.array(y)n”, ” gs = G.samples(x.shape[0])n”, ” rho_gin = vmap(rin.rho_dense)(gs)n”, ” rho_gout = vmap(rout.rho_dense)(gs)n”, ” y1 = model((rho_gin@x[…,None])[…,0],training=False)n”, ” y2 = (rho_gout@model(x,training=False)[…,None])[…,0]n”, ” return rel_err(y1,y2)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“As expected, the network continues to be equivariant as it is trained.”
]
}, {
“cell_type”: “code”, “execution_count”: 8, “metadata”: {}, “outputs”: [
- {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“Average test equivariance error 4.34e-07n”
]
}
], “source”: [
“print(f"Average test equivariance error {np.mean([equivariance_err(mb) for mb in testloader]):.2e}")”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## Breaking EMLP down into equivariant layers (mid level)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Internally for EMLP, we use representations that [uniformly allocate dimensions](https://emlp.readthedocs.io/en/latest/package/emlp.models.mlp.html#emlp.models.uniform_rep) between different tensor representations.”
]
}, {
“cell_type”: “code”, “execution_count”: 9, “metadata”: {}, “outputs”: [
- {
“name”: “stdout”, “output_type”: “stream”, “text”: [
“122V⁰+40V+12V²+3V³+V⁴n”
]
}
], “source”: [
“from emlp.nn import uniform_repn”, “r = uniform_rep(512,G)n”, “print(r)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Below is a trimmed down version of EMLP, so you can see how it is built from the component layers Linear, BiLinear, and GatedNonlinearities. These layers can be constructed like ordinary objax modules, using the input and output representations.”
]
}, {
“cell_type”: “code”, “execution_count”: 12, “metadata”: {}, “outputs”: [], “source”: [
“from objax.module import Modulen”, “n”, “class EMLPBlock(Module):n”, ” """ Basic building block of EMLP consisting of G-Linear, biLinear,n”, ” and gated nonlinearity. """n”, ” def __init__(self,rep_in,rep_out):n”, ” super().__init__()n”, ” rep_out_wgates = emlp.nn.gated(rep_out)n”, ” self.linear = emlp.nn.Linear(rep_in,rep_out_wgates)n”, ” self.bilinear = emlp.nn.BiLinear(rep_out_wgates,rep_out_wgates)n”, ” self.nonlinearity = emlp.nn.GatedNonlinearity(rep_out)n”, ” def __call__(self,x):n”, ” lin = self.linear(x)n”, ” preact =self.bilinear(lin)+linn”, ” return self.nonlinearity(preact)n”, “n”, “class EMLP(Module):n”, ” def __init__(self,rep_in,rep_out,group,ch=384,num_layers=3):n”, ” super().__init__()n”, ” reps = [rep_in(group)]+num_layers*[uniform_rep(ch,group)]n”, ” self.network = emlp.nn.Sequential(n”, ” *[EMLPBlock(rin,rout) for rin,rout in zip(reps,reps[1:])],n”, ” emlp.nn.Linear(reps[-1],rep_out(group))n”, ” )n”, ” def __call__(self,x,training=True):n”, ” return self.network(x)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“The representations of the hidden layers (taking place of the number of channels in a standard MLP) is by default given by this uniform_rep shown above. Unlike this pedagogical implementation you can specify the representation of the hidden layers directly in the [full EMLP](https://emlp.readthedocs.io/en/latest/package/emlp.nn.html#emlp.nn.EMLP) by feeding in a representation to the ch argument, or even a list of representations to specify each hidden layer.”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Note that since we are using the GatedNonlinearity, additional scalar gate channels need to be added to the output representation for the layer directly before the nonlinearity (in this case the Linear layer) which can be achieved with the gated function.”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## The equivariant linear layers (low level)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“At a lower level, the implementation of the Linear is fairly straightforward. An unconstrained bias b and weight matrix w are initialized. The projection matrices $P_b$ and $P_w$ are computed which are used project onto the symmetric subspace for each. Finally, during the forward pass, the unconstrained parameters are reshaped to vectors, projected via the matrices, and reshaped back to the original sizes. Then these projected parameters are applied to the input like a standard linear layer.”
]
}, {
“cell_type”: “code”, “execution_count”: 13, “metadata”: {}, “outputs”: [], “source”: [
“from objax.variable import TrainVarn”, “from objax.nn.init import orthogonaln”, “n”, “class Linear(Module):n”, ” """ Basic equivariant Linear layer from repin to repout."""n”, ” def __init__(self, repin, repout):n”, ” nin,nout = repin.size(),repout.size()n”, ” self.b = TrainVar(objax.random.uniform((nout,))/jnp.sqrt(nout))n”, ” self.w = TrainVar(orthogonal((nout, nin)))n”, ” self.rep_W = rep_W = repout*repin.Tn”, ” n”, ” self.Pb = repout.equivariant_projector() # the bias vector has representation repoutn”, ” self.Pw = rep_W.equivariant_projector()n”, ” n”, ” def __call__(self, x):n”, ” W = (self.Pw@self.w.value.reshape(-1)).reshape(*self.w.value.shape)n”, ” b = self.Pb@self.b.valuen”, ” return x@W.T+b”
]
}
], “metadata”: {
- “kernelspec”: {
“display_name”: “Python 3”, “language”: “python”, “name”: “python3”
}, “language_info”: {
- “codemirror_mode”: {
“name”: “ipython”, “version”: 3
}, “file_extension”: “.py”, “mimetype”: “text/x-python”, “name”: “python”, “nbconvert_exporter”: “python”, “pygments_lexer”: “ipython3”, “version”: “3.8.5”
}
}, “nbformat”: 4, “nbformat_minor”: 2
}