{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Phase-Only Compressive Sensing\n", "\n", "This notebook aims at testing the possibility to recovery a signal from phase-only observations, as described in the papers:\n", "\n", "* [1] L. Jacques, T. Feuillen, \"The importance of phase in complex compressive sensing\", 2020 ([arXiv](https://arxiv.org/abs/2001.02529))\n", "* [2] L. Jacques, T. Feuillen, \"Keep the phase! Signal recovery in phase-only compressive sensing\", 2020 ([arXiv](https://arxiv.org/abs/2011.06499))\n", "\n", "The code relies on the \"basis pursuit denoising\" (BPDN) solver implemented by [\"SPGL1: A solver for sparse least squares\"](https://friedlander.io/spgl1/), and in particular its [port to Python](https://github.com/drrelyea/spgl1)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As described in [1,2], the phase-only compressive sensing is defined from the following model:\n", "\n", "$$\n", "\\boldsymbol z = {\\rm sign}_{\\mathbb C}(\\boldsymbol \\Phi \\boldsymbol x)\n", "\\tag{PO-CS}\n", "$$\n", "\n", "with ${\\rm sign}_{\\mathbb C}(\\lambda) = \\lambda/|\\lambda|$ the complex sign operator, a complex Gaussian random matrix $\\boldsymbol \\Phi \\sim \\mathbb C \\mathcal N^{m \\times n}(0,2)$, and $\\boldsymbol x$ an $k$-sparse vector in $\\mathbb R^n$.\n", "\n", "The papers above develop a method to recover, up to a global amplitude, $\\boldsymbol x$ from $\\boldsymbol z$.\n", "\n", "**Remark:** All the codes below are based on a complex Gaussian matrix $\\boldsymbol \\Phi$. However, complex Bernoulli and random partial Fourier sensing are also implemented through the function generate_sensing_mtx(m,n,model='gaussian'), with model in {'gaussian','bernouilli','fourier'}. Feel free to play with them." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Requirements**:\n", "* Numpy, SciPy, MatplotLib\n", "* [SPGL1](https://github.com/drrelyea/spgl1)\n", "* [tqdm](https://github.com/tqdm/tqdm) (just for progress bars ;-)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Preliminaries**: Loading important modules. If you cannot run or install the progress bar of tqdm, just comment the lines where it's used. See below." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy import interpolate\n", "from scipy.sparse.linalg import LinearOperator\n", "\n", "import matplotlib.pyplot as plt\n", "import matplotlib.colors as mcolors\n", "import matplotlib.cm as cm\n", "\n", "#from scipy.sparse.linalg import LinearOperator\n", "#from scipy.sparse import spdiags\n", "import spgl1\n", "\n", "# for nice progress bars (comment it if it doesn't work)\n", "from tqdm.notebook import tqdm\n", "\n", "\n", "# Initialize random number generators\n", "np.random.seed(43273289)\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Noiseless PO-CS reconstruction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us define a few important functions, namely the complex sign operator ${\\rm sign}_{\\mathbb C}$, the constant $\\kappa = \\sqrt{\\pi/2}$, the signal-to-noise ratio metric, as well as a shortcut to call the BPDN solver from SPGL1." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# numerical error constant: 10*np.finfo(np.float32).eps\n", "epsilon = 10*np.finfo(np.float32).eps \n", "\n", "# complex sign operator\n", "signc = lambda x: x / (np.absolute(x) + (x == 0)*epsilon) \n", "\n", "# the kappa constant (see [1])\n", "kappa = np.sqrt(np.pi/2)\n", "\n", "# the SNR metric between two signals (in dB)\n", "snr = lambda a,b: 20*np.log10(np.linalg.norm(a)/np.linalg.norm(a-b)) \n", "\n", "# simplifying the call to BP from SPGL1 and storing parameters in the same time\n", "#bpalg = lambda A, b: spgl1.spg_bp(A, b, opt_tol=1e-7, bp_tol=1e-7)[0]\n", "#bpalg = lambda A, b: spgl1.spg_bp(A, b, opt_tol=1e-8, bp_tol=1e-8, iscomplex=False)[0]\n", "bpalg = lambda A, b: spgl1.spg_bp(A, b, opt_tol=1e-7, bp_tol=1e-7)[0]\n", "\n", "# simplifying the call to BPDN from SPGL1 and storing parameters in the same time\n", "bpdnalg = lambda A, b, sigma: spgl1.spg_bpdn(A, b, sigma)[0]\n", "\n", "# a function to generate the j^th canonical vector in d dimensions\n", "def generate_canonic_vec(d,j):\n", " e_j = np.zeros(d)\n", " e_j[j] = 1.0\n", " return e_j" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A function to generate real $k$-sparse signals in $\\mathbb R^n$." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def generate_sparse_sig(k,n):\n", " idx = np.random.permutation(n)\n", " x = np.zeros(n)\n", " x[idx[0:k]] = np.random.randn(k)\n", " return x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's define several functions generating $m \\times n$ random sensing matrix constructions, namely complex Gaussian, complex Bernoulli, and partial Fourier random matrices." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def generate_complex_Gaussian_mtx(m,n):\n", " Phi = np.random.randn(m, n) + np.random.randn(m, n) *1j;\n", " A = Phi/np.sqrt(m)\n", " return A\n", "\n", "def generate_complex_Bernoulli_mtx(m,n):\n", " Phi = np.sign(np.random.randn(m, n)) + np.sign(np.random.randn(m, n)) *1j;\n", " A = Phi/np.sqrt(m)\n", " return A\n", "\n", "class generate_partial_Fourier_mtx(LinearOperator):\n", " def __init__(self, m, n):\n", " self.m = m\n", " self.n = n\n", " self.shape = (m, n)\n", " self.dtype = np.complex128\n", " self.idx = np.random.permutation(n)[0:m]\n", " \n", " def _matvec(self, x):\n", " z = np.fft.fft(x) / np.sqrt(self.m) \n", " return z[self.idx]\n", " \n", " def _rmatvec(self, y):\n", " z = np.zeros(self.n, dtype=complex)\n", " z[self.idx] = y\n", " return np.fft.ifft(z) * self.n / np.sqrt(self.m)\n", "\n", "def generate_sensing_mtx(m,n,model='gaussian'):\n", " if model == 'gaussian':\n", " return generate_complex_Gaussian_mtx(m,n)\n", " elif model == 'bernoulli':\n", " return generate_complex_Bernoulli_mtx(m,n)\n", " elif model == 'fourier':\n", " return generate_partial_Fourier_mtx(m,n)\n", " else:\n", " print(\"unknown model\")\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A function to turn complex linear operators (applied to real vectors) into real linear operators with two times as many rows." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "class Op2Real(LinearOperator):\n", " def __init__(self, A):\n", " self.A = A\n", " m, n = A.shape\n", " self.shape = [2*m, n]\n", " self.dtype = np.float64\n", " \n", " def _matvec(self, x):\n", " return np.concatenate(((self.A @ x).real,\n", " (self.A @ x).imag))\n", " \n", " def _rmatvec(self, y):\n", " m = self.A.shape[0]\n", " return (self.A.T @ y[0:m]).real + (self.A.T @ y[m:]).imag\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Following [1], this function constructs the equivalent matrix $\\boldsymbol A_{\\boldsymbol z}$ from both the phase measurements $\\boldsymbol z$ and the sensing matrix $\\boldsymbol A$." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "class generate_Az(LinearOperator):\n", " def __init__(self, A, z):\n", " self.A = A\n", " self.ismat = (type(A) == np.ndarray)\n", " self.z = z\n", " self.D = np.diag(z)\n", " self.shape = (A.shape[0]+2, A.shape[1])\n", " self.dtype = np.float64 \n", " \n", " def _matvec(self, x):\n", " \n", " m = self.A.shape[0]\n", " n = self.A.shape[1]\n", " \n", " H_z_x = self.D.conj().T @ (self.A @ x)\n", " \n", " alpha_adj_x = np.array([np.ones(m) @ H_z_x / (kappa * np.sqrt(m))])\n", " \n", " return np.concatenate([alpha_adj_x.real, alpha_adj_x.imag, H_z_x.imag])\n", " \n", " def _rmatvec(self, y):\n", " # passed the adjoint test\n", " m = self.A.shape[0]\n", " n = self.A.shape[1]\n", " \n", " if self.ismat:\n", " alpha = (self.A.conj().T @ (self.D @ np.ones(m))) / (kappa * np.sqrt(m))\n", " return (alpha.real * y[0]) - (alpha.imag * y[1]) - (self.A.conj().T @ (self.D @ y[2:])).imag\n", "\n", " else:\n", " alpha = (self.A.H @ (self.D @ np.ones(m))) / (kappa * np.sqrt(m))\n", " return (alpha.real * y[0]) - (alpha.imag * y[1]) - (self.A.H @ (self.D @ y[2:])).imag\n", " \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.1 Single reconstruction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now define the sensing model." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Dimensions\n", "n = 128 # space dimension\n", "m = 50 # number of measurements\n", "k = 10 # sparsity level\n", "\n", "# Generating the complex Gaussian matrix\n", "#A = generate_complex_Gaussian_mtx(m,n)\n", "sensing_model = 'gaussian'\n", "\n", "A = generate_sensing_mtx(m, n, model = sensing_model)\n", "\n", "# Generating a k-sparse signal, random k-length support, and Gaussian entries\n", "x = generate_sparse_sig(k,n)\n", "\n", "# Linear measurements\n", "y = A @ x\n", "\n", "# Phase only measurements\n", "z = signc(y)\n", "\n", "# Definition of the matrix A_z for sparse signal recovery in PO-CS, such that A_z x0 = (1, 0, ..., 0) \n", "A_z = generate_Az(A,z)\n", "\n", "#A_z = build_Az(A,z)\n", "\n", "# Let x0 such that ||Ax0||_1 = kappa sqrt(m) (see [1])\n", "# This is the signal we hope to estimate, as the amplitude of x is lost (see [1])\n", "x0 = x * kappa * np.sqrt(m) / np.linalg.norm(y,1)\n", "\n", "# The measurement vector e_1 such that A_z x0 = e_1\n", "e_1 = generate_canonic_vec(m+2,0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can compute that $\\|\\boldsymbol A \\boldsymbol x_0 - \\boldsymbol e_1\\|$ is 0, up to numerical errors:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "||Ax0||_1 = 1.0000000000000002\n", "||Ax0 - e_1|| = 2.056073113926681e-16\n" ] } ], "source": [ "print(f\"||Ax0||_1 = {np.linalg.norm(A @ x0,1)/(kappa*np.sqrt(m))}\")\n", "print(f\"||Ax0 - e_1|| = {np.linalg.norm((A_z @ x0) - e_1)}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now solve the linear and phase-only basis pursuit (BP) problems, defined by, respectively\n", "\n", "$$\n", "\\boldsymbol x_{\\rm lin} \\in \\arg\\min_{\\boldsymbol u} \\|\\boldsymbol u\\|_1\\ {\\rm s.t.}\\ \\boldsymbol A \\boldsymbol u = \\boldsymbol y,\n", "\\tag{BP/Lin}\n", "$$\n", "\n", "$$\n", "\\boldsymbol x_{\\rm PO-CS} \\in \\arg\\min_{\\boldsymbol u} \\|\\boldsymbol u\\|_1\\ {\\rm s.t.}\\ \\boldsymbol A_z \\boldsymbol u = \\boldsymbol e_1,\n", "\\tag{BP/PO-CS}\n", "$$\n", "\n", "rembering that $\\boldsymbol x_{\\rm lin}$ should estimate $\\boldsymbol x$, and $\\boldsymbol x_{\\rm PO-CS}$ should estimate $\\boldsymbol x_0$." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The SNR between x_lin and x is 138.6 dB\n", "The SNR between x_pocs and x0 is 110.8458148 dB\n" ] } ], "source": [ "# Linear model (note: running SPGL1 on complex mtx seems wrong). \n", "A_equiv = Op2Real(A)\n", "y_equiv = A_equiv @ x\n", "x_lin = bpalg(A_equiv, y_equiv)\n", "\n", "# Phase-only model\n", "x_pocs = bpalg(A_z, e_1)\n", "\n", "# SNRs\n", "print(f\"The SNR between x_lin and x is {snr(x_lin, x):.4} dB\")\n", "print(f\"The SNR between x_pocs and x0 is {snr(x_pocs,x0):.10} dB\")" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAApw0lEQVR4nO3dfXQV9Z348fcnyQ0xgKiASA0k0AKK8mhwfS4VqqAVZNEfcLIVjnVTtXZxd6tLT7aurebUVVdba2ubWp/TqsUWo9VVGq26PhIUUXkQRB6CKCEoxQRIQj6/P2YmuQk3T8xNZm7m8zrnnjsz95uZb+4lHz738535jqgqxhhjer+0oDtgjDGmZ1jAN8aYiLCAb4wxEWEB3xhjIsICvjHGRERG0B1oy6BBgzQvLy/obhhjTEpZuXLlLlUdnOi10Ab8vLw8Kioqgu6GMcakFBHZ0tZrVtIxxpiIsIBvjDERYQHfGGMiIrQ1/ETq6+uprKxk//79QXcl8rKyssjJySEWiwXdFWNMJ6VUwK+srKR///7k5eUhIkF3J7JUlerqaiorKxkxYkTQ3THGdFJKlXT279/PwIEDLdgHTEQYOHCgfdMKWmkp5OVBWprzXFoadI9MyKVUhg9YsA8J+xwCVloKhYVQW+usb9nirAMUFATXLxNqKZXhG2NcRUVobS0P8W1qOcLZVlsLRUXB9suEmgX8FHPjjTdy++23H7J92bJlrFmzpsv727x5M7///e+b1h944AGuueYaX300PWDrVjYwioU8xJPMbrHdmLb07oAfUI2zoaGhR44Tr72A315/Wgd8kyKGD+fvHAnAXvq32G5MW3pvwPdqnFu2gGpzjdNn0L/pppsYM2YMZ511FgsWLGjKtqdOncq1115Lfn4+P//5zykvL2fSpEmMGzeOyy+/nAMHDgDOlBG7du0CoKKigqlTpwJO5n755ZczdepURo4cyV133dV0zOLiYkaPHs1ZZ53F+vXrD+nTa6+9RllZGddddx0TJ07ko48+OqQ/ixYtYunSpU0/069fPwCWLFnCK6+8wsSJE7nzzjsB+OSTT5gxYwajRo3i+uuv9/V+mW5SXExNn4EA1NDX2ZadDcXFAXbKhF3KDdp2WlFR84CWx6txHuag1ooVK3jiiSd49913qa+vZ/LkyZxyyilNr9fV1VFRUcH+/fsZNWoU5eXljB49mssuu4x77rmHa6+9tt39r1u3jhdffJG9e/cyZswYrrrqKlavXs2jjz7KqlWraGhoOOSYAGeccQazZs3iW9/6Fpdccskh/QFYtGhRwmPecsst3H777Tz99NOAU9JZtWoV77zzDn369GHMmDF8//vfZ9iwYYfxjpluU1BA7btfgduglr6Qm+sEexuwNe3ovRl+W7VMHzXOV199ldmzZ5OVlUX//v256KKLWrw+b948ANavX8+IESMYPXo0AAsXLuTll1/ucP8XXnghffr0YdCgQRx77LF89tlnvPLKK8yZM4fs7GyOPPJIZs2a1en+ev3pqmnTpjFgwACysrIYO3YsW7a0OReTCVDNqd9wnn94M2zebMHedCgpAV9EZojIehHZKCJL2mk3V0RURPKTcdx2tVXL7MYaZ9++fTtsk5GRQWNjI8Ah57H36dOnaTk9Pd33WEB8f+KP29jYSF1dXZs/l+x+mO7hfYFt/UXWmLb4Dvgikg78EpgJjAUWiMjYBO36A4uBN/0es1OKi52aZjyfNc4zzzyTp556iv379/Pll182lUFaGzNmDJs3b2bjxo0APPzww3z9618HnBr+ypUrAXjiiSc6POY555zDsmXL2LdvH3v37uWpp55K2K5///7s3bu3zf3EH7esrIz6+vpO/ZwJr5qals+hZheJhUIyMvxTgY2quklV64BHIf48sSY3Af8N9MzlmQUFUFLi1DZFnOeSEl9fe6dMmcKsWbMYP348M2fOZNy4cQwYMOCQdllZWdx///1ceumljBs3jrS0NK688koA/uu//ovFixeTn59Penp6h8ecPHky8+bNY8KECcycOZMpU6YkbDd//nxuu+02Jk2axEcffXTI6//8z//MSy+9xIQJE3j99debsv/x48eTnp7OhAkTmgZtTWpImQy/m06gMIdBVX09gEuAe+PWvw3c3arNZOAJd/lvQH4b+yoEKoCK4cOHa2tr1qw5ZFtP27t3r6qq1tTU6CmnnKIrV64MuEfBCcPnEWU//rEqqM6eHXRPOpCbqwp6Kz/QNZzgdBqc7SbpgAptI153+6CtiKQBdwD/3lFbVS1R1XxVzR88OOEdugJXWFjIxIkTmTx5MnPnzmXy5MlBd8lEVMpk+Fu3UkeM67mNP7CgxXbTs5JxWuZ2IP6cvRx3m6c/cDLwN3f+leOAMhGZpaopdw9Du0jJhEXK1PCHD6dui3PtiV0kFqxkZPgrgFEiMkJEMoH5QJn3oqruUdVBqpqnqnnAG0BKBntjwiRlMvziYuqOOAqAL3Eu+LOLxILhO8NX1QYRuQZ4DkgH7lPVD0TkJzi1pLL292CMORwpE/ALCqjfkwXfcwO+XSQWmKRcaauqzwDPtNp2QxttpybjmMZEXcqUdIC6i+bC92DvhQvg6QUd/4DpFr33SltjermUyfAB7zq/L78Mth9RZwG/Gz3wwAN88sknTetXXHHFYU1h3JrNcGkgxTJ8C/ihYAG/G7UO+Pfeey9jxx5yEXKXWcA30JzZ19VB2Ge/8AK+XdQdLAv4h+GRRx7h1FNPZeLEiXz3u9/l4MGDLFq0iJNPPplx48Zx5513snTpUioqKigoKGDixIns27ePqVOnNs1e2a9fP6677jpOOukkpk+fzltvvdU0NXJZmTPOvXnzZs4++2wmT57M5MmTee2114BDpzQ+ePAg1113HVOmTGH8+PH85je/Cey9MT0nPrMPe1nHncnDMvyApez0yNdeC6tWJXefEyfCz37Wfpu1a9fy2GOP8eqrrxKLxbj66qu5+eab2b59O++//z4AX3zxBUcddRR33303t99+O/n5h84VV1NTw7nnnsttt93GnDlz+M///E+WL1/OmjVrWLhwIbNmzeLYY49l+fLlZGVlsWHDBhYsWEBFRcUhUxqXlJQwYMAAVqxYwYEDBzjzzDM577zzGDFiRHLfIBMq8UG+thaOPDK4vnTESjrhkLIBPyjl5eWsXLmyaU6bffv2MWPGDDZt2sT3v/99LrzwQs4777wO95OZmcmMGTMAGDduHH369CEWizFu3Dg2b94MQH19Pddccw2rVq0iPT2dDz/8MOG+nn/+eVavXt10g5M9e/awYcMGC/i9XE0NHHMM7N4d/jp+fMBXdaa3Mj0vZQN+R5l4d1FVFi5cyE9/+tMW24uLi3nuuef49a9/zeOPP859993X7n5isRjulcekpaU1TUmclpbWNB3xnXfeyZAhQ3j33XdpbGwkKyurzT794he/4Pzzz/f765kUUlsLX/2qE/DDXtLxAn5jI+zbd+hEtqZnWA2/i6ZNm8bSpUvZuXMnALt372bLli00NjYyd+5cbr75Zt5++23A/9TDe/bsYejQoaSlpfHwww9z8ODBhPs9//zzueeee5qmPP7www+pCXvKZ3zxBmq9KafC/nHH337ByjrBSdkMPyhjx47l5ptv5rzzzqOxsZFYLMYdd9zBnDlzmm4w4mX/ixYt4sorr+SII47g9ddf7/Kxrr76aubOnctDDz3EjBkzEk5pvGjRIhYvXszmzZuZPHkyqsrgwYNZtmxZ0n5nEz5eRu8F/LBn+N6gLTgB/9hjg+tLlIkzm2b45Ofnq3dGi2ft2rWceOKJAfXItGafR3C2b4ecHGda+ZISKCuDVnfcDJXHHoP5853lVatgwoRAu9OrichKVU14V0Er6RiTglItw7eSTjhYwDcmBXk1e6vhm65IuRq+qjad3RKI6mrn+3RdHWRmwvHHw8CBwfUnIGEtBUZFqmX48TV8u9o2OCmV4WdlZVFdXR1csKmudu7H6aUrdXXOenV1MP0JiKpSXV3d5mmipvtZhm8OR0pl+Dk5OVRWVlJVVRVMByor0YMH2cMA+vMl6TinSVJd7YygRUhWVhY5Efudw8TL6I85puV6WFnAD4eUCvixWCzYq0dPOon1OoozWM/dfI/v8Stnu4hzRYkxPcTL6Pv1cy5iSqUM30o6wUlKSUdEZojIehHZKCJLErz+byKyRkRWi0i5iOQm47g9bvhwqnC+Q3vP3nZjepKX0WdnQ9++qZPhp6dbhh8k3wFfRNKBXwIzgbHAAhFpPQfwO0C+qo4HlgK3+j1uIIqLqcp0yhhNAd/uzWkC4GX0ffumRoZfXw8ZGc43Egv4wUlGhn8qsFFVN6lqHfAoMDu+gaq+qKpeDvIGkJrF34ICdn37WgB2Mci5N2dJid2b0/S4VMzwMzOhf38r6QQpGTX844FtceuVwD+00/47wLNJOG4gqkae5jx/Yx68MC/g3pioqqmBtDTo0yc1Mnwv4FuGH6weHbQVkX8C8oGvt/F6IVAIMDykdXHvBKGgThQyBpyMPjvbOV8glTJ8C/jBSkZJZzswLG49x93WgohMB4qAWap6INGOVLVEVfNVNX/w4MGJmgRu166Wz8YEoabGCfSQGhl+fb2VdMIgGQF/BTBKREaISCYwHyiLbyAik4Df4AT7nUk4ZmC8zH7XLudGDsYEwcvwIXUy/FjMMvyg+Q74qtoAXAM8B6wFHlfVD0TkJyIyy212G9AP+KOIrBKRsjZ2F3pewG9ogC++CLQrJsJqa1tm+KkQ8L0M3wJ+cJJSw1fVZ4BnWm27IW55ejKOEwa7djmZSn29s3z00UH3yERRTU3LDD/sJZ34Gr6VdIKTUnPphEFVFYwe3bxsTBBSNcO3kk6wLOB3QU2Ncz/Ose5lZRbwTVBaZ/i1teEeU4oftK2psZlIgmIBvwu8M3O8mzzZmTomKK0zfFXYvz/YPrUnftAWwl+C6q0s4HeBl9F7Ad8yfBOU1hm+ty2s4ks6YGWdoFjA7wIvo8/NhSOOsAzfBKd1hu9tC6v4s3TAAn5QLOB3gZfRDx7sPCzDTyGlpZCX58xHkJfnrKew+Azfew5zhu/V8L0M387UCUZKzYcfNC/ADxrkPCzgp4jSUigsbE6Bt2xx1iElJ75rbHTq9V6G7z2HPcOPxSzDD5pl+F3gnYM/YICT4VtJJ0UUFUFtLT9lCWVc5GyrrXW2p6D4mTLjn8Oc4VsNPxws4HdBVZWT2YtYSSelbN1KHTFu5EZ+zuIW21ORF/BTLcO3kk7wLOB3gRfwwUo6KWX4cN7nZOrow0pOoRFp2p6KvEw+FTN8K+kEywJ+F+za5WT24Dx7F2KZkCsupiLzDAD2cBQf8dWUvlNZKmb4rQdtLeAHwwJ+F1RVtQz4YHX8lFBQQMWZixGcyztXDJqZ0ncqS9UMP/7CKyvpBMMCfhe0Lul420z4VXwxim+cm0ZWFlR8+66UDfZw6KBtKmT4XkknFnPu0mUZfjAs4HdSQwN8/rll+Klo/3547z047TSYNAkqKoLukT/xNzCH8Gf4qs0BH2wCtSBZwO+k6mrnuXXAT7kMv5ddgNQZq1c7/2Hn5zuPt9+GgweD7tXha53hZ2Q4wTSsGX5Dg/McH/CtpBMMC/idFH/RVfxzSgV87wKkLVuctMu7AKmXB30vo/cCfk0NrF8fbJ/8aJ3hQ7hvc1hf7zx7Ad9ughIcC/id5JVuvMz+6KMhPT3FSjpFRWytHcg8HmUVE5xtKXwBUmdVVMCxx0JOjhPwvW2pqnWGD+G+zWFdnfMciznPVtIJTlICvojMEJH1IrJRRJYkeL2PiDzmvv6miOQl47gJeSULEee7bvyzNydCR9sSvFY1bR4AgwvOg9JS0v5QykCtoqr4N0nZf0+0//uW3XyLp3mceVzIX9jOV5z3bMuWlOj/4e6r4v7V5O99Afl9KWNW/p6+UkPFwl+E73frZPuaa/4DgL5Txjr/3ktLyf5sEzW/+0PX9++V9RL93SSj/3l51P1+KeBm+KWl9F/1Mnuffy05/WmvfRg+S690evXVzaXUrn42yaSqvh5AOvARMBLIBN4FxrZqczXwa3d5PvBYR/s95ZRTtMseeUQ1O1vVKVgk9fFLrlJQ3cEQ1VhMNTNTx/K+zuGJbjlesh/1pOtM/qLp1OvP+Bftx991Eit1L30D71t3Pmo4QtNo0Bu4selzO5uX9HReDbxvh/u4kRsUVBtIa/qdJvCOzmLZ4e3T3Ud39Xdb1tcUVH/7nddVs7N1Dk/oON5NXn+6uf+BPrKznbjWBUBFW3E1Iwn/Z5wKbFTVTQAi8igwG1gT12Y2cKO7vBS4W0TE7VzyFBVRVZtNLskvrNcTQ2hkINVQ74xCHctOlnEx2YS0eBqnkTQOkMVvKKSQ3zKKDVzEUwykmnRSeASzA42k0Ug6p7CyqZicTwV38m8p8bklUkcmR1BLOo1Q71xb0J+9PMVFXfqdruFubuU/mt6X03mNd71SXxL8O//DTdxA/X7n7yVz2eNQW8uR/J33GH9IX6/nVm7kx1BfTwPpTGK1c5FcG37IT/kRN0N9PXXEmMgHbCYvaf33K0Y9TzCX6ZQDsIoJTKOcfRzR4c9OYQUvMbW55Jqk04iTEfCPB7bFrVcC/9BWG1VtEJE9wECgRQVcRAqBQoDhh3PZ+9atHEE213B313+2E05gHTEamtZ/wg085U3GlQLGs5p/wvmKeAHP8hQX8TemBtupHtCPLzmf55rW/4W76MMBDpIeYK/8Gcd7LdZ/wg08y8xO//wfWEAF+S22vck/cBpvcBb/57t/pRTwphsG6nBGa2PVnwLwr9zJsexs0f6PXMoLnOsEfGAHQ3mfcVzAXziJDw7Z/x+5lOc5zwn4wGbyWMtYLqKME1jnu//JcCf/yguc2xTwX+YcdjOQxfyMTOra/dnhxM3zlMw5n/yUc9wE/RLg3rj1bwN3t2rzPpATt/4RMKi9/R5WSSc3N/ivYPawRwo8vkG5ns1LTesHEQXVH/OjpOx/Gsv1DP5PFXQ1JyuoLh303Tbb/xMPaR6bmtZf5XQF1Wc5P2H7hdyvx7Otaf1/OU9B9WXOCvy99R4nsEYv5k9N61fyKz2K3drY1X3l5nYpDNJOSScZg7bbgWFx6znutoRtRCQDGABUJ+HYLRUXtzx1obvEYs3nmKWiVO//4eqNv/dh/k4ZNFBPrGkf9bG+TduTIZtaanD2WZd1JACZl81v8+8zh0q2c7wzsV0sxrb0EU3bExnJJraTw376QCzGpvTRTdvD4kTWso4TmtbXcQInstabuq9zkjznUzIC/gpglIiMEJFMnEHZslZtyoCF7vIlwAvu/0TJVVDgzJGSm+usp6e3fB440Hl0tK2913Jz4f774b77Dj1OMvbf3e1Tvf+Hu6/Wv7dI+H63HvwsMzhIAxlN+2i457fu9oak9L+vuAE/N5e6/7gBgMzzpjb/fbZ6/4fJdurJZGfOKXD//VTO+zdnO9sS7n+EbAFgy1fOgPvv5+MZV9GHAwxlRzg+m4EDOeGILWzka9QP/ypcdRXr0k/iBNZ3fv+5ucmf8+lwyjitH8AFwIc4pZoid9tPgFnuchbwR2Aj8BYwsqN9HlZJxxjTKbNmqU6Y0Lz++edO9eCOO5Kz/+98R3XoUGf5Jbdy9MILbbd/8kmnzVtvOeuLF6v266fa2Ji4/f851SJ95hlnfe5c1TFjktP3ZHnoIaePa9c2v7+33tr9x6Wbz9JBVZ8Bnmm17Ya45f3Apck4ljHGv4yM5ikPoHk5IykRwbkQzLvyt/WFV4kMc4vClZUwZYrzPGyY80UgkRFOxYePP25+9raFxQluNWftWvjii5bbgpKkj9cYk0piscQBv72g3BXxV/56Ab+9oYacHOd527bmZ29bIscd58y6GR/wTzvNX5+TzQvu69bBnj3O8oknBtcfsIBvTCRlZDTPcQPNy8nM8BsanGDfmYA/aJATwOMD/sknt90+Lc3J6DdtcrLnzz8PX4bfvz8cf3xzhp+Z6Vw8GyQL+MZEUE+UdMAp67SePC0RESejr6x02n/6afsZPjgB/uOPm7P8kSP99zvZTjyxOcMfPTp57+/hssnTjImgWCxxhp/Mkg44Ab8zGT44Nftt2+CTT5wT0IcNa7/9yJFOhr/JPRMzbBk+OGWddetgzZrg6/dgAd+YSOrJDL8zg7bQnOFXVjavt2fECCdzfvttZz2sGf7evbBxY/D1e7CAb0wk9cSgLXQ9w9++3Zm81Vtvj5fRl5c705UPGHD4/e0u8Vm9ZfjGmEB096Bt/G0XuxLwGxpg5crm9fZ4Gf2KFeEs5wCcuOaJ5uUfXBj4zYYs4BsTQW2VdLojw+/MoC00l3Bef905w+XII9tv7wX5xsZwlnMoLeW46y/jSJxzMkfv+Fvgd5izgG9MBLU1aBtkDd/L6Feu7Di7B6eEc8wxznIoM/yiImRfLSewjuFsoS+1gd9hzk7LNCaCMjKcG7mrOqdEdtegbW1t1wZtwWnf0YAtAKWljPhyHLsZz8jfFcGEscmdd8Yvd1rjJdzCHgYcsj0IluEbE0Fe8D3o3vumu0/LTE9vnhOsLQMHQlaWs9xhhl9aCoWFjKhz7kY/YndF4OWSQ7j39JjDMhbx4CHbg2AB35gI8jJ5L9B392mZnZnB2bv4CjqR4RcVQW1t03TII9kUeLnkEImma0/ydMddZQHfmAjyArsX6Lt70LZTU/aXljJs22sADPvF9e1n625Z5EL+wgyeJY/NLbaHQvx07SLdM91xF1nANyaCvMDuZfjJHrRNT3fmxvEy/A7/I3FLNDkHNgKQs/vd9ks0blnkHF7hWS5ovvVogOWShAoKYPNm51SizZsDH2OwgG9MBLWV4SdzrhdviuROlXTcEs0w9/bYw9jWfokmhOWSVGAB35gI8jJuL9Ane9AWnPjb6YDvlmLO5FW+ysaOSzQhLJekAl8BX0SOEZHlIrLBfT46QZuJIvK6iHwgIqtFZJ6fYxpj/OvuQVtozvA7VcN3SzEX8CwbGUU2+1psTyhk5ZJU4DfDXwKUq+oooNxdb60WuExVTwJmAD8TkaN8HtcY40PrDD/Zg7bQxZKOlWh6hN+APxuaTjB9ELi4dQNV/VBVN7jLnwA7gcE+j2uM8aF1hp/sQVtovutVpwZtrUTTI/x+vENUdYe7/CkwpL3GInIqkIlzs/NErxcChQDDwzbabkwv0lODttXVztk6nTots6DAAnw36/DjFZG/AscleKnF8LmqqohoO/sZCjwMLFTVxkRtVLUEKAHIz89vc1/GGH96YtC2b19nzLVfv04GfNPtOgz4qjq9rddE5DMRGaqqO9yAvrONdkcCfwGKVPWNw+6tMSYpenrQtk+f5O3XHD6/NfwyYKG7vBB4snUDEckE/gw8pKpLfR7PGJMEPZXhd2VqBdP9/Ab8W4BvisgGYLq7jojki8i9bpv/B5wDLBKRVe5jos/jGmN86MkMv1ODtqZH+Pp4VbUamJZgewVwhbv8CPCIn+MYY5Ir0aCtSMczWnZF376wbx/s328ZfljYlbbGRFCikk4ys3toPq1+zx4L+GFhAd+YCEpU0kl2wPdmzPz8cwv4YWEB35gISpThJ7vO7gX8/futhh8WFvCNiaCezPDBMvywsIBvTAQlGrTtrgwfLOCHhQV8YyKoJwZtLeCHjwV8YyLISjrRZAHfmAjqyUHb+OOZYFnANyaCLMOPJgv4xkSQDdpGkwV8YyKoJ6+0BQv4YWEB35gISlTSSXaGn5XlzM8DFvDDwgK+MRHUExm+SHNZxwZtw8ECvjER5M2K2Z2DttAc8C3DDwcL+MZEUFqa8+jO0zLBAn7YWMA3JqJisZZn6ViG3/v5CvgicoyILBeRDe7z0e20PVJEKkXkbj/HNMYkR0ZG9w7agtXww8Zvhr8EKFfVUUC5u96Wm4CXfR7PGJMk8Rl+dwzagmX4YeM34M8GHnSXHwQuTtRIRE4BhgDP+zyeMSZJWmf4FvB7P78Bf4iq7nCXP8UJ6i2ISBrwP8APOtqZiBSKSIWIVFRVVfnsmjGmPa0zfBu07f06/D9dRP4KHJfgpaL4FVVVEdEE7a4GnlHVSvGuwmiDqpYAJQD5+fmJ9mWMSZKMjO4ftPWutrWAHw4dfsSqOr2t10TkMxEZqqo7RGQosDNBs9OBs0XkaqAfkCkiX6pqe/V+Y0w3s0Hb6PFb0ikDFrrLC4EnWzdQ1QJVHa6qeThlnYcs2BsTvG4ftC0tpe8DvwQg8/xvQGlpkg9guspvwL8F+KaIbACmu+uISL6I3Ou3c8aY7tOtGX5pKRQW0vfvnwCQ+cnHUFhoQT9gvgK+qlar6jRVHaWq01V1t7u9QlWvSND+AVW9xs8xjTHJ0a0ZflER1NbSlxoAMqmD2lpnuwmMXWlrTER162mZW7cCcA4vcxFlDKaqxXYTDAv4xkRU/Fk6ST8tc/hwACbzDmXMJkZDi+0mGBbwjYmobp1Lp7i45R1QwFkvLk7iQUxXWcA3JqK8kk5jI6gmOcMvKICSEsjNdSbGz8111gsKkngQ01XdcKmFMSYVxGJw4EBzHT/pp2UWFFiADxnL8I2JKC/D98o63XGlrQkXC/jGRJQ3aOtl+HY1bO9nAd+YiPIGbS3Djw4L+MZEVOuSjmX4vZ8FfGMiysvwu23Q1oSOBXxjIsoy/OixgG9MRLUetLUMv/ezgG9MRNmgbfRYwDcmoqykEz0W8I2JKBu0jR4L+MZElGX40eMr4IvIMSKyXEQ2uM9Ht9FuuIg8LyJrRWSNiOT5Oa4xxj8btI0evxn+EqBcVUcB5e56Ig8Bt6nqicCpJL7ZuTGmB1lJJ3r8BvzZwIPu8oPAxa0biMhYIENVlwOo6peqWuvzuMYYn7wAf+CA82wlnd7Pb8Afoqo73OVPgSEJ2owGvhCRP4nIOyJym4ikJ9qZiBSKSIWIVFRVVfnsmjGmPV6A37fPebYMv/fr8CMWkb8CxyV4qcXdiFVVRUTbOMbZwCRgK/AYsAj4XeuGqloClADk5+cn2pcxJkm8AO8FfMvwe78OA76qTm/rNRH5TESGquoOERlK4tp8JbBKVTe5P7MMOI0EAd8Y03Msw48evyWdMmChu7wQeDJBmxXAUSIy2F0/F1jj87jGGJ+8AL9/f8t103v5Dfi3AN8UkQ3AdHcdEckXkXsBVPUg8AOgXETeAwT4rc/jGmN8spJO9Pj6P11Vq4FpCbZXAFfErS8Hxvs5ljEmuaykEz12pa0xEWUZfvRYwDcmoizDjx4L+MZEVOtBW8vwez8L+MZEVOuSjmX4vZ8FfGMiyko60WMB35iIskHb6LGAb0xEWYYfPRbwjYmo+EFbEUhPOKWh6U0s4BsTUfElHcvuo8ECvjERFV/SsYAfDRbwjYmo+AzfBmyjwQK+MRFlGX70WMA3JqLiB20tw48GC/jGRJQN2kaPBXxjIsrL6g8csAw/KizgGxNR8Vm9ZfjR4Cvgi8gxIrJcRDa4z0e30e5WEflARNaKyF0iIn6Oa4zxLz6rt4AfDX4z/CVAuaqOAsrd9RZE5AzgTJw7Xp0MTAG+7vO4xhif4oO8lXSiwW/Anw086C4/CFycoI0CWUAm0AeIAZ/5PK4xxicr6USP34A/RFV3uMufAkNaN1DV14EXgR3u4zlVXZtoZyJSKCIVIlJRVVXls2vGmPbEZ/WW4UdDh/+vi8hfgeMSvFQUv6KqKiKa4Oe/BpwI5LiblovI2ar6Suu2qloClADk5+cfsi9jTPJYhh89HX7Mqjq9rddE5DMRGaqqO0RkKLAzQbM5wBuq+qX7M88CpwOHBHxjTM+xgB89fks6ZcBCd3kh8GSCNluBr4tIhojEcAZsE5Z0jDE9Jy3NeYCVdKLCb8C/BfimiGwAprvriEi+iNzrtlkKfAS8B7wLvKuqT/k8rjEmCbxAbxl+NPj6mFW1GpiWYHsFcIW7fBD4rp/jGGO6R0aGXWkbJXalrTER5mX2luFHgwV8YyLMy+wtw48GC/jGRJhl+NFiAd+YCLNB22ixgG9MhHmB3ko60WAB35gIs5JOtFjANybCbNA2WizgGxNhluFHiwV8YyLMBm2jxQK+MRFmg7bRYgHfmAizkk60WMA3JsJs0DZaLOAbE2GW4UeLBXxjIswy/GixgG9MhFmGHy0W8I2JMAv40eIr4IvIpSLygYg0ikh+O+1miMh6EdkoIkv8HNMYkzxW0okWvxn++8A/Ai+31UBE0oFfAjOBscACERnr87jGmCSwDD9a/N7icC2AiLTX7FRgo6pucts+CswG1vg5tjHGP8vwo6UnavjHA9vi1ivdbcaYgFmGHy0dfswi8lfguAQvFanqk8nsjIgUAoUAw4cPT+aujTEJ2Fw60dLhx6yq030eYzswLG49x92W6FglQAlAfn6++jyuMaYDNpdOtPRESWcFMEpERohIJjAfKOuB4xpjOmAlnWjxe1rmHBGpBE4H/iIiz7nbvyIizwCoagNwDfAcsBZ4XFU/8NdtY0wy2KBttPg9S+fPwJ8TbP8EuCBu/RngGT/HMsYkn2X40WJX2hoTYTZoGy0W8I2JMBu0jRYL+MZEmJV0osUCvjERZoO20WIB35gIsww/WizgGxNVpaXEbvoRALGZ06G0NOAOme5mAd+YKCothcJCMnbvBCBjx1YoLLSg38tZwDcmioqKoLaWGPUAZNAAtbXOdtNrWcA3Joq2bgVgJs9SxM3ksbnFdtM7WcA3Jorc2WiH8ik38yPS0BbbTe9kAd+YKCouhuzsltuys53tpteygG9MFBUUQEkJ5OaCiPNcUuJsN72WnX1rTFQVFFiAjxjL8I0xJiIs4BtjTERYwDfGmIiwgG+MMRFhAd8YYyJCVDXoPiQkIlXAFh+7GATsSlJ3gmD9D04q9x2s/0ELuv+5qjo40QuhDfh+iUiFquYH3Y/DZf0PTir3Haz/QQtz/62kY4wxEWEB3xhjIqI3B/ySoDvgk/U/OKncd7D+By20/e+1NXxjjDEt9eYM3xhjTBwL+MYYExG9LuCLyAwRWS8iG0VkSdD96YiIDBORF0VkjYh8ICKL3e3HiMhyEdngPh8ddF/bIyLpIvKOiDztro8QkTfdz+ExEckMuo9tEZGjRGSpiKwTkbUicnoqvf8i8q/uv533ReQPIpIV5vdfRO4TkZ0i8n7ctoTvtzjucn+P1SIyObieN/U1Uf9vc//9rBaRP4vIUXGv/dDt/3oROT+QTrt6VcAXkXTgl8BMYCywQETGBturDjUA/66qY4HTgO+5fV4ClKvqKKDcXQ+zxcDauPX/Bu5U1a8BnwPfCaRXnfNz4H9V9QRgAs7vkRLvv4gcD/wLkK+qJwPpwHzC/f4/AMxota2t93smMMp9FAL39FAf2/MAh/Z/OXCyqo4HPgR+COD+Lc8HTnJ/5ldunApErwr4wKnARlXdpKp1wKPA7ID71C5V3aGqb7vLe3GCzfE4/X7QbfYgcHEgHewEEckBLgTuddcFOBdY6jYJbf9FZABwDvA7AFWtU9UvSKH3H+e+FkeISAaQDewgxO+/qr4M7G61ua33ezbwkDreAI4SkaE90tE2JOq/qj6vqg3u6htAjrs8G3hUVQ+o6sfARpw4FYjeFvCPB7bFrVe621KCiOQBk4A3gSGqusN96VNgSFD96oSfAdcDje76QOCLuD+AMH8OI4Aq4H63JHWviPQlRd5/Vd0O3A5sxQn0e4CVpM7772nr/U7Fv+nLgWfd5VD1v7cF/JQlIv2AJ4BrVfXv8a+pc+5sKM+fFZFvATtVdWXQfTlMGcBk4B5VnQTU0Kp8E/L3/2icLHIE8BWgL4eWG1JKmN/vjohIEU6ZtjToviTS2wL+dmBY3HqOuy3URCSGE+xLVfVP7ubPvK+u7vPOoPrXgTOBWSKyGaeEdi5OTfwot8QA4f4cKoFKVX3TXV+K8x9Aqrz/04GPVbVKVeuBP+F8Jqny/nvaer9T5m9aRBYB3wIKtPkCp1D1v7cF/BXAKPcMhUycwZKygPvULrfe/TtgrareEfdSGbDQXV4IPNnTfesMVf2hquaoah7O+/2CqhYALwKXuM3C3P9PgW0iMsbdNA1YQ4q8/zilnNNEJNv9t+T1PyXe/zhtvd9lwGXu2TqnAXviSj+hISIzcMqas1S1Nu6lMmC+iPQRkRE4g89vBdFHAFS1Vz2AC3BGyT8CioLuTyf6exbO19fVwCr3cQFOHbwc2AD8FTgm6L524neZCjztLo/E+Ye9Efgj0Cfo/rXT74lAhfsZLAOOTqX3H/gxsA54H3gY6BPm9x/4A854Qz3ON6zvtPV+A4Jz5t1HwHs4ZyOFsf8bcWr13t/wr+PaF7n9Xw/MDLLvNrWCMcZERG8r6RhjjGmDBXxjjIkIC/jGGBMRFvCNMSYiLOAbY0xEWMA3xpiIsIBvjDER8f8BCRccTauvimMAAAAASUVORK5CYII=\n", "text/plain": [ "