Two-Bit-Err-vs-Noise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Small-scale synthetic experiments with 2-bit sparse vectors\n",
    "### Example: estimation error as a function of noise level\n",
    "\n",
    "This notebook presents a small-scale example that demonstrates the advantage of an exhaustive method for solving our QUBO formualtion to the sparse coding problem. For demonstration purposes the code below runs on a classic computer. It is important to emphasize that our QUBO problem can be minimized on quantum computers or other quantum-inspired solvers, see our paper for larger-scale experiments.\n",
    "\n",
    "We start by generating a matrix $A$ with $L_2$-normalized columns such that its mutual coherence is minimized.\n",
    "Then, given $A$, we generate a vector $\\bar{b}=Ax$ by drawing an $x$ with $k$ randomly chosen non-zero  **2-bit** entries that are equal to one, and then form its noisy version $b = \\bar{b} + v$, where the noise component $v$ is sampled from a normal distribution with zero mean  standard deviation (STD) $\\sigma$. Now, given $A$ and the noisy $b$, we attempt to reconstruct $x$ via (1) lasso, (2) OMP, (3) an exhaustive algorithm for solving the original sparse coding problem, and (4) an exhaustive QUBO solver for the original sparse coding problem. \n",
    "\n",
    "Due to the high computational complexity of the exhaustive solvers, we run this experiment only for a matrix $A$ with a small the number of columns, which we set to $N=8$. The hyper-parameter of each method is tuned by an oracle that has access to the true $x$, optimizing each of the two error metrics presented in the graphs.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Type A:\n",
      "  low_coherence\n",
      "Type x0:\n",
      "  fixed_point\n",
      "N = 8 ,  K = 3 ,  M =  4\n",
      "x_vals =\n",
      "  [0 1 2 3]\n",
      "N = 8 ,  K = 3 ,  M =  4\n",
      "x0_min = 0 ,  x0_max = 3 ,  dx0 = 1\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from tqdm import tqdm\n",
    "import os\n",
    "import sys\n",
    "\n",
    "base_code_path = \"./\"\n",
    "\n",
    "sys.path.append(base_code_path)\n",
    "sys.path.append(base_code_path + \"/solvers\")\n",
    "sys.path.append(base_code_path + \"/utils\")\n",
    "sys.path.append(base_code_path + \"/solvers_experiments\")\n",
    "\n",
    "pd.set_option('display.max_columns', None)\n",
    "\n",
    "from utils import utils\n",
    "\n",
    "#########################\n",
    "# Experiment parameters #\n",
    "#########################\n",
    "\n",
    "#Type of A\n",
    "type_A = \"low_coherence\"\n",
    "\n",
    "# Number of (features) variables\n",
    "N = 8\n",
    "# Number of measurements\n",
    "M = 4\n",
    "# Number of non-zero variables (cardinality)\n",
    "K = 3\n",
    "\n",
    "# Set whether and how much Gaussian noise N(0, 1) to add to the measurement\n",
    "# vector\n",
    "noise_vec = [0.0, 0.05, 0.1, 0.15, 0.2, 0.3]\n",
    "\n",
    "# Set the threshold to the absolute value of and element of a solution\n",
    "# vector, under which this element is considered as zero for cardinality\n",
    "L0_th = 1e-1\n",
    "\n",
    "# number of trials per experiemnt\n",
    "num_exps = 20\n",
    "\n",
    "# binary x0 (0) or P-bits (1)?\n",
    "type_x0_id = 1\n",
    "TYPE_X_LIST =  [\"binary\",\n",
    "                \"fixed_point\"]\n",
    "type_x0 = TYPE_X_LIST[type_x0_id]\n",
    "\n",
    "print(\"Type A:\\n  \" + type_A)\n",
    "print(\"Type x0:\\n  \" + type_x0)\n",
    "sys.stdout.flush()\n",
    "\n",
    "# internal normalization of A (lasso)\n",
    "lasso_normalize = True\n",
    "\n",
    "# Set the parameters for uniform discrete\n",
    "# The values of the solution must include 0 and their number must be\n",
    "# a positive power of 2\n",
    "if type_x0 == \"fixed_point\":\n",
    "    x0_min = 0\n",
    "    d_x0 = 1\n",
    "    P = 2\n",
    "    w_fp = 10\n",
    "\n",
    "    x_ip_0_spins = utils.check_fixed_point_parameters(x0_min, d_x0, P)\n",
    "\n",
    "    # Set the parameters to generate a valid solution to the linear system\n",
    "    params_x0 = (x0_min, d_x0, P)\n",
    "\n",
    "    # %% Display info\n",
    "    print(\"N =\", N, \",  K =\", K, \",  M = \", M)\n",
    "    x_vals = x0_min + d_x0 * np.arange(2 ** P)\n",
    "    print(\"x_vals =\")\n",
    "    print(\" \", x_vals)\n",
    "    x0_max = max(x_vals)\n",
    "    num_bits = P\n",
    "\n",
    "elif type_x0 == \"binary\":\n",
    "    x0_min = 0\n",
    "    x0_max = 1\n",
    "    num_bits = 1\n",
    "    d_x0 = 1\n",
    "    w_fp = 50\n",
    "else:\n",
    "    1/0\n",
    "\n",
    "params_x0 = (x0_min, x0_max + 1e-10, d_x0)\n",
    "\n",
    "print(\"N =\", N, \",  K =\", K, \",  M = \", M)\n",
    "print(\"x0_min =\", x0_min, \",  x0_max =\", x0_max, ',  dx0 =', d_x0)\n",
    "\n",
    "\n",
    "# arrrange params for the solver\n",
    "params = dict()\n",
    "params[\"x0_min\"]   = x0_min\n",
    "params[\"x0_max\"]   = x0_max\n",
    "params[\"num_bits\"] = num_bits\n",
    "params[\"d_x0\"]     = d_x0\n",
    "params[\"w_fp\"]     = w_fp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "noise_id: 1/6 exp_id: 1/20\n",
      "Init loss = 0.20562297105789185 , final loss = 0.1249999850988388\n",
      "noise_id: 1/6 exp_id: 2/20\n",
      "Init loss = 0.181496262550354 , final loss = 0.125\n",
      "noise_id: 1/6 exp_id: 3/20\n",
      "Init loss = 0.2060803771018982 , final loss = 0.125\n",
      "noise_id: 1/6 exp_id: 4/20\n",
      "Init loss = 0.33006545901298523 , final loss = 0.12499997764825821\n",
      "noise_id: 1/6 exp_id: 5/20\n",
      "Init loss = 0.19870223104953766 , final loss = 0.125\n",
      "noise_id: 1/6 exp_id: 6/20\n",
      "Init loss = 0.2521541714668274 , final loss = 0.125\n",
      "noise_id: 1/6 exp_id: 7/20\n",
      "Init loss = 0.1775641143321991 , final loss = 0.1249999925494194\n",
      "noise_id: 1/6 exp_id: 8/20\n",
      "Init loss = 0.16283631324768066 , final loss = 0.125\n",
      "noise_id: 1/6 exp_id: 9/20\n",
      "Init loss = 0.22313155233860016 , final loss = 0.1249999925494194\n",
      "noise_id: 1/6 exp_id: 10/20\n",
      "Init loss = 0.1852712631225586 , final loss = 0.125\n",
      "noise_id: 1/6 exp_id: 11/20\n",
      "Init loss = 0.2721059322357178 , final loss = 0.125\n",
      "noise_id: 1/6 exp_id: 12/20\n",
      "Init loss = 0.20138363540172577 , final loss = 0.125\n",
      "noise_id: 1/6 exp_id: 13/20\n",
      "Init loss = 0.29642254114151 , final loss = 0.125\n",
      "noise_id: 1/6 exp_id: 14/20\n",
      "Init loss = 0.19850453734397888 , final loss = 0.125\n",
      "noise_id: 1/6 exp_id: 15/20\n",
      "Init loss = 0.2513588070869446 , final loss = 0.125\n",
      "noise_id: 1/6 exp_id: 16/20\n",
      "Init loss = 0.15429526567459106 , final loss = 0.1249999925494194\n",
      "noise_id: 1/6 exp_id: 17/20\n",
      "Init loss = 0.24222373962402344 , final loss = 0.125\n",
      "noise_id: 1/6 exp_id: 18/20\n",
      "Init loss = 0.24154523015022278 , final loss = 0.1249999925494194\n",
      "noise_id: 1/6 exp_id: 19/20\n",
      "Init loss = 0.17316149175167084 , final loss = 0.1249999925494194\n",
      "noise_id: 1/6 exp_id: 20/20\n",
      "Init loss = 0.19232887029647827 , final loss = 0.1249999850988388\n",
      "noise_id: 2/6 exp_id: 1/20\n",
      "Init loss = 0.20562297105789185 , final loss = 0.1249999850988388\n",
      "noise_id: 2/6 exp_id: 2/20\n",
      "Init loss = 0.181496262550354 , final loss = 0.125\n",
      "noise_id: 2/6 exp_id: 3/20\n",
      "Init loss = 0.2060803771018982 , final loss = 0.125\n",
      "noise_id: 2/6 exp_id: 4/20\n",
      "Init loss = 0.33006545901298523 , final loss = 0.12499997764825821\n",
      "noise_id: 2/6 exp_id: 5/20\n",
      "Init loss = 0.19870223104953766 , final loss = 0.125\n",
      "noise_id: 2/6 exp_id: 6/20\n",
      "Init loss = 0.2521541714668274 , final loss = 0.125\n",
      "noise_id: 2/6 exp_id: 7/20\n",
      "Init loss = 0.1775641143321991 , final loss = 0.1249999925494194\n",
      "noise_id: 2/6 exp_id: 8/20\n",
      "Init loss = 0.16283631324768066 , final loss = 0.125\n",
      "noise_id: 2/6 exp_id: 9/20\n",
      "Init loss = 0.22313155233860016 , final loss = 0.1249999925494194\n",
      "noise_id: 2/6 exp_id: 10/20\n",
      "Init loss = 0.1852712631225586 , final loss = 0.125\n",
      "noise_id: 2/6 exp_id: 11/20\n",
      "Init loss = 0.2721059322357178 , final loss = 0.125\n",
      "noise_id: 2/6 exp_id: 12/20\n",
      "Init loss = 0.20138363540172577 , final loss = 0.125\n",
      "noise_id: 2/6 exp_id: 13/20\n",
      "Init loss = 0.29642254114151 , final loss = 0.125\n",
      "noise_id: 2/6 exp_id: 14/20\n",
      "Init loss = 0.19850453734397888 , final loss = 0.125\n",
      "noise_id: 2/6 exp_id: 15/20\n",
      "Init loss = 0.2513588070869446 , final loss = 0.125\n",
      "noise_id: 2/6 exp_id: 16/20\n",
      "Init loss = 0.15429526567459106 , final loss = 0.1249999925494194\n",
      "noise_id: 2/6 exp_id: 17/20\n",
      "Init loss = 0.24222373962402344 , final loss = 0.125\n",
      "noise_id: 2/6 exp_id: 18/20\n",
      "Init loss = 0.24154523015022278 , final loss = 0.1249999925494194\n",
      "noise_id: 2/6 exp_id: 19/20\n",
      "Init loss = 0.17316149175167084 , final loss = 0.1249999925494194\n",
      "noise_id: 2/6 exp_id: 20/20\n",
      "Init loss = 0.19232887029647827 , final loss = 0.1249999850988388\n",
      "noise_id: 3/6 exp_id: 1/20\n",
      "Init loss = 0.20562297105789185 , final loss = 0.1249999850988388\n",
      "noise_id: 3/6 exp_id: 2/20\n",
      "Init loss = 0.181496262550354 , final loss = 0.125\n",
      "noise_id: 3/6 exp_id: 3/20\n",
      "Init loss = 0.2060803771018982 , final loss = 0.125\n",
      "noise_id: 3/6 exp_id: 4/20\n",
      "Init loss = 0.33006545901298523 , final loss = 0.12499997764825821\n",
      "noise_id: 3/6 exp_id: 5/20\n",
      "Init loss = 0.19870223104953766 , final loss = 0.125\n",
      "noise_id: 3/6 exp_id: 6/20\n",
      "Init loss = 0.2521541714668274 , final loss = 0.125\n",
      "noise_id: 3/6 exp_id: 7/20\n",
      "Init loss = 0.1775641143321991 , final loss = 0.1249999925494194\n",
      "noise_id: 3/6 exp_id: 8/20\n",
      "Init loss = 0.16283631324768066 , final loss = 0.125\n",
      "noise_id: 3/6 exp_id: 9/20\n",
      "Init loss = 0.22313155233860016 , final loss = 0.1249999925494194\n",
      "noise_id: 3/6 exp_id: 10/20\n",
      "Init loss = 0.1852712631225586 , final loss = 0.125\n",
      "noise_id: 3/6 exp_id: 11/20\n",
      "Init loss = 0.2721059322357178 , final loss = 0.125\n",
      "noise_id: 3/6 exp_id: 12/20\n",
      "Init loss = 0.20138363540172577 , final loss = 0.125\n",
      "noise_id: 3/6 exp_id: 13/20\n",
      "Init loss = 0.29642254114151 , final loss = 0.125\n",
      "noise_id: 3/6 exp_id: 14/20\n",
      "Init loss = 0.19850453734397888 , final loss = 0.125\n",
      "noise_id: 3/6 exp_id: 15/20\n",
      "Init loss = 0.2513588070869446 , final loss = 0.125\n",
      "noise_id: 3/6 exp_id: 16/20\n",
      "Init loss = 0.15429526567459106 , final loss = 0.1249999925494194\n",
      "noise_id: 3/6 exp_id: 17/20\n",
      "Init loss = 0.24222373962402344 , final loss = 0.125\n",
      "noise_id: 3/6 exp_id: 18/20\n",
      "Init loss = 0.24154523015022278 , final loss = 0.1249999925494194\n",
      "noise_id: 3/6 exp_id: 19/20\n",
      "Init loss = 0.17316149175167084 , final loss = 0.1249999925494194\n",
      "noise_id: 3/6 exp_id: 20/20\n",
      "Init loss = 0.19232887029647827 , final loss = 0.1249999850988388\n",
      "noise_id: 4/6 exp_id: 1/20\n",
      "Init loss = 0.20562297105789185 , final loss = 0.1249999850988388\n",
      "noise_id: 4/6 exp_id: 2/20\n",
      "Init loss = 0.181496262550354 , final loss = 0.125\n",
      "noise_id: 4/6 exp_id: 3/20\n",
      "Init loss = 0.2060803771018982 , final loss = 0.125\n",
      "noise_id: 4/6 exp_id: 4/20\n",
      "Init loss = 0.33006545901298523 , final loss = 0.12499997764825821\n",
      "noise_id: 4/6 exp_id: 5/20\n",
      "Init loss = 0.19870223104953766 , final loss = 0.125\n",
      "noise_id: 4/6 exp_id: 6/20\n",
      "Init loss = 0.2521541714668274 , final loss = 0.125\n",
      "noise_id: 4/6 exp_id: 7/20\n",
      "Init loss = 0.1775641143321991 , final loss = 0.1249999925494194\n",
      "noise_id: 4/6 exp_id: 8/20\n",
      "Init loss = 0.16283631324768066 , final loss = 0.125\n",
      "noise_id: 4/6 exp_id: 9/20\n",
      "Init loss = 0.22313155233860016 , final loss = 0.1249999925494194\n",
      "noise_id: 4/6 exp_id: 10/20\n",
      "Init loss = 0.1852712631225586 , final loss = 0.125\n",
      "noise_id: 4/6 exp_id: 11/20\n",
      "Init loss = 0.2721059322357178 , final loss = 0.125\n",
      "noise_id: 4/6 exp_id: 12/20\n",
      "Init loss = 0.20138363540172577 , final loss = 0.125\n",
      "noise_id: 4/6 exp_id: 13/20\n",
      "Init loss = 0.29642254114151 , final loss = 0.125\n",
      "noise_id: 4/6 exp_id: 14/20\n",
      "Init loss = 0.19850453734397888 , final loss = 0.125\n",
      "noise_id: 4/6 exp_id: 15/20\n",
      "Init loss = 0.2513588070869446 , final loss = 0.125\n",
      "noise_id: 4/6 exp_id: 16/20\n",
      "Init loss = 0.15429526567459106 , final loss = 0.1249999925494194\n",
      "noise_id: 4/6 exp_id: 17/20\n",
      "Init loss = 0.24222373962402344 , final loss = 0.125\n",
      "noise_id: 4/6 exp_id: 18/20\n",
      "Init loss = 0.24154523015022278 , final loss = 0.1249999925494194\n",
      "noise_id: 4/6 exp_id: 19/20\n",
      "Init loss = 0.17316149175167084 , final loss = 0.1249999925494194\n",
      "noise_id: 4/6 exp_id: 20/20\n",
      "Init loss = 0.19232887029647827 , final loss = 0.1249999850988388\n",
      "noise_id: 5/6 exp_id: 1/20\n",
      "Init loss = 0.20562297105789185 , final loss = 0.1249999850988388\n",
      "noise_id: 5/6 exp_id: 2/20\n",
      "Init loss = 0.181496262550354 , final loss = 0.125\n",
      "noise_id: 5/6 exp_id: 3/20\n",
      "Init loss = 0.2060803771018982 , final loss = 0.125\n",
      "noise_id: 5/6 exp_id: 4/20\n",
      "Init loss = 0.33006545901298523 , final loss = 0.12499997764825821\n",
      "noise_id: 5/6 exp_id: 5/20\n",
      "Init loss = 0.19870223104953766 , final loss = 0.125\n",
      "noise_id: 5/6 exp_id: 6/20\n",
      "Init loss = 0.2521541714668274 , final loss = 0.125\n",
      "noise_id: 5/6 exp_id: 7/20\n",
      "Init loss = 0.1775641143321991 , final loss = 0.1249999925494194\n",
      "noise_id: 5/6 exp_id: 8/20\n",
      "Init loss = 0.16283631324768066 , final loss = 0.125\n",
      "noise_id: 5/6 exp_id: 9/20\n",
      "Init loss = 0.22313155233860016 , final loss = 0.1249999925494194\n",
      "noise_id: 5/6 exp_id: 10/20\n",
      "Init loss = 0.1852712631225586 , final loss = 0.125\n",
      "noise_id: 5/6 exp_id: 11/20\n",
      "Init loss = 0.2721059322357178 , final loss = 0.125\n",
      "noise_id: 5/6 exp_id: 12/20\n",
      "Init loss = 0.20138363540172577 , final loss = 0.125\n",
      "noise_id: 5/6 exp_id: 13/20\n",
      "Init loss = 0.29642254114151 , final loss = 0.125\n",
      "noise_id: 5/6 exp_id: 14/20\n",
      "Init loss = 0.19850453734397888 , final loss = 0.125\n",
      "noise_id: 5/6 exp_id: 15/20\n",
      "Init loss = 0.2513588070869446 , final loss = 0.125\n",
      "noise_id: 5/6 exp_id: 16/20\n",
      "Init loss = 0.15429526567459106 , final loss = 0.1249999925494194\n",
      "noise_id: 5/6 exp_id: 17/20\n",
      "Init loss = 0.24222373962402344 , final loss = 0.125\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "noise_id: 5/6 exp_id: 18/20\n",
      "Init loss = 0.24154523015022278 , final loss = 0.1249999925494194\n",
      "noise_id: 5/6 exp_id: 19/20\n",
      "Init loss = 0.17316149175167084 , final loss = 0.1249999925494194\n",
      "noise_id: 5/6 exp_id: 20/20\n",
      "Init loss = 0.19232887029647827 , final loss = 0.1249999850988388\n",
      "noise_id: 6/6 exp_id: 1/20\n",
      "Init loss = 0.20562297105789185 , final loss = 0.1249999850988388\n",
      "noise_id: 6/6 exp_id: 2/20\n",
      "Init loss = 0.181496262550354 , final loss = 0.125\n",
      "noise_id: 6/6 exp_id: 3/20\n",
      "Init loss = 0.2060803771018982 , final loss = 0.125\n",
      "noise_id: 6/6 exp_id: 4/20\n",
      "Init loss = 0.33006545901298523 , final loss = 0.12499997764825821\n",
      "noise_id: 6/6 exp_id: 5/20\n",
      "Init loss = 0.19870223104953766 , final loss = 0.125\n",
      "noise_id: 6/6 exp_id: 6/20\n",
      "Init loss = 0.2521541714668274 , final loss = 0.125\n",
      "noise_id: 6/6 exp_id: 7/20\n",
      "Init loss = 0.1775641143321991 , final loss = 0.1249999925494194\n",
      "noise_id: 6/6 exp_id: 8/20\n",
      "Init loss = 0.16283631324768066 , final loss = 0.125\n",
      "noise_id: 6/6 exp_id: 9/20\n",
      "Init loss = 0.22313155233860016 , final loss = 0.1249999925494194\n",
      "noise_id: 6/6 exp_id: 10/20\n",
      "Init loss = 0.1852712631225586 , final loss = 0.125\n",
      "noise_id: 6/6 exp_id: 11/20\n",
      "Init loss = 0.2721059322357178 , final loss = 0.125\n",
      "noise_id: 6/6 exp_id: 12/20\n",
      "Init loss = 0.20138363540172577 , final loss = 0.125\n",
      "noise_id: 6/6 exp_id: 13/20\n",
      "Init loss = 0.29642254114151 , final loss = 0.125\n",
      "noise_id: 6/6 exp_id: 14/20\n",
      "Init loss = 0.19850453734397888 , final loss = 0.125\n",
      "noise_id: 6/6 exp_id: 15/20\n",
      "Init loss = 0.2513588070869446 , final loss = 0.125\n",
      "noise_id: 6/6 exp_id: 16/20\n",
      "Init loss = 0.15429526567459106 , final loss = 0.1249999925494194\n",
      "noise_id: 6/6 exp_id: 17/20\n",
      "Init loss = 0.24222373962402344 , final loss = 0.125\n",
      "noise_id: 6/6 exp_id: 18/20\n",
      "Init loss = 0.24154523015022278 , final loss = 0.1249999925494194\n",
      "noise_id: 6/6 exp_id: 19/20\n",
      "Init loss = 0.17316149175167084 , final loss = 0.1249999925494194\n",
      "noise_id: 6/6 exp_id: 20/20\n",
      "Init loss = 0.19232887029647827 , final loss = 0.1249999850988388\n"
     ]
    }
   ],
   "source": [
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(action='ignore')\n",
    "\n",
    "###################\n",
    "# Core execution #\n",
    "###################\n",
    "\n",
    "from solvers_experiments import random_linear_system as RLS\n",
    "from solvers import annealer\n",
    "from solvers import others\n",
    "\n",
    "def get_stats(s_L, r_L, r_L_test):\n",
    "    # Find the best solution, minimizing the L2 distance between x_hat and x_0\n",
    "    tmp = s_L[:,2]\n",
    "    ind_best_solution_x = np.argmin(tmp)\n",
    "    rec_err_x = tmp[ind_best_solution_x]\n",
    "    rmse = r_L[ind_best_solution_x]\n",
    "\n",
    "    # Find the best solution, minimizing the L0 distance between x_hat and x_0\n",
    "    tmp = s_L[:,0]\n",
    "    ind_best_solution_x = np.argmin(tmp)\n",
    "    cardinality_err_x = s_L[ind_best_solution_x,0]\n",
    "\n",
    "\n",
    "    # Find the best CV solution, minimizing the test system error\n",
    "    ind_best_solution_cv_sys = np.argmin(r_L_test)\n",
    "    cv_rec_err_x = s_L[ind_best_solution_cv_sys,2]\n",
    "    cv_cardinality_err_x = s_L[ind_best_solution_cv_sys,0]\n",
    "    cv_rmse = r_L_test[ind_best_solution_cv_sys]\n",
    "\n",
    "    return rec_err_x, cardinality_err_x, rmse, cv_rec_err_x, cv_cardinality_err_x, cv_rmse\n",
    "\n",
    "\n",
    "# Initialize table of results\n",
    "results = pd.DataFrame()\n",
    "\n",
    "\n",
    "methods_to_run = [\"exhaustive_L0\", \"exhaustive_qubo_L0\", \"lasso_L1\", \"omp_L0\"]\n",
    "\n",
    "for noise_id, b_noise_coef in enumerate(noise_vec):\n",
    "    for exp_id in range(num_exps):\n",
    "        print(\"noise_id: \" + str(noise_id+1) + \"/\" + str(len(noise_vec)) + \" exp_id: \" + str(exp_id+1) + \"/\" + str(num_exps))\n",
    "        sys.stdout.flush()\n",
    "\n",
    "        # %% Form the system of equations\n",
    "        [A, b, x0] = \\\n",
    "            RLS.generate_random_linear_system(N=N, M=M, K=K,\n",
    "                                              type_A=type_A,\n",
    "                                              type_x0=type_x0, params_x0=params_x0,\n",
    "                                              b_noise_coef=b_noise_coef, seed=exp_id)\n",
    "\n",
    "        A_test=A\n",
    "        b_test=b\n",
    "\n",
    "\n",
    "        # %% Apply L0 regularization (classical)\n",
    "        if \"exhaustive_L0\" in methods_to_run:\n",
    "            x_L0, r_L0, r_L0_test, k_L0, s_L0 = others.exhaustive_L0(A, b,\n",
    "                                                                     A_test, b_test,\n",
    "                                                                     L0_th=L0_th, x0=x0)\n",
    "\n",
    "\n",
    "            rec_err_L0, cardinality_err_L0, rmse_L0, \\\n",
    "                cv_rec_err_L0, cv_cardinality_err_L0, cv_rmse_L0 = \\\n",
    "                    get_stats(s_L0, r_L0, r_L0_test)\n",
    "\n",
    "            # Store results\n",
    "            results = results.append({'Experiment':exp_id,\n",
    "                                      'Cardinality': K,\n",
    "                                      'Noise STD': b_noise_coef,\n",
    "                                      'Method':'Exhaustive L0',\n",
    "                                      'Infeasible X-Recovery L2 Err':rec_err_L0,\n",
    "                                      'Infeasible X-Recovery L0 Err':cardinality_err_L0,\n",
    "                                      'Infeasible System RMSE':rmse_L0,\n",
    "                                      'CV X-Recovery L2 Err':cv_rec_err_L0,\n",
    "                                      'CV X-Recovery L0 Err':cv_cardinality_err_L0,\n",
    "                                      'CV System RMSE':cv_rmse_L0,\n",
    "                                      'Type A':type_A,\n",
    "                                      'Normalize A':lasso_normalize,\n",
    "                                      'L0 Thresh':L0_th,\n",
    "                                      'N':N,\n",
    "                                      'M':M,\n",
    "                                      'x0_min':x0_min,\n",
    "                                      'x0_max':x0_max,\n",
    "                                      'num_bits':num_bits,\n",
    "                                      'd_x0':d_x0,\n",
    "                                      'w_fp':w_fp,\n",
    "                                      'type_x0':type_x0},\n",
    "                                      ignore_index=True)\n",
    "\n",
    "\n",
    "\n",
    "        # %% Apply L0 regularization (QUBO)\n",
    "        if \"exhaustive_qubo_L0\" in methods_to_run:\n",
    "            lambda_max = sum(np.abs(A.T@b))\n",
    "            eps = 1e-8\n",
    "            n_grid = 20\n",
    "            lam_L0 = np.logspace(np.log10(lambda_max * eps), np.log10(lambda_max), num=n_grid)[::-1]\n",
    "            np.append(lam_L0, 1e-10)\n",
    "            np.append(lam_L0, 1e-15)\n",
    "\n",
    "\n",
    "            x_qubo, r_qubo, r_qubo_test, k_qubo, s_qubo = \\\n",
    "                annealer.exhaustive_qubo_Lp(A, b,\n",
    "                                            A_test, b_test,\n",
    "                                            norm=\"L0\",\n",
    "                                            params=params,\n",
    "                                            lam_Lp=lam_L0,\n",
    "                                            L0_th=L0_th, x0=x0)\n",
    "\n",
    "\n",
    "\n",
    "            rec_err_qubo, cardinality_err_qubo, rmse_qubo, \\\n",
    "                cv_rec_err_qubo, cv_cardinality_err_qubo, cv_rmse_qubo = \\\n",
    "                    get_stats(s_qubo, r_qubo, r_qubo_test)\n",
    "\n",
    "            # Store results\n",
    "            results = results.append({'Experiment':exp_id,\n",
    "                                      'Cardinality': K,\n",
    "                                      'Noise STD': b_noise_coef,\n",
    "                                      'Method':'QUBO L0',\n",
    "                                      'Infeasible X-Recovery L2 Err':rec_err_qubo,\n",
    "                                      'Infeasible X-Recovery L0 Err':cardinality_err_qubo,\n",
    "                                      'Infeasible System RMSE':rmse_qubo,\n",
    "                                      'CV X-Recovery L2 Err':cv_rec_err_qubo,\n",
    "                                      'CV X-Recovery L0 Err':cv_cardinality_err_qubo,\n",
    "                                      'CV System RMSE':cv_rmse_qubo,\n",
    "                                      'Type A':type_A,\n",
    "                                      'Normalize A':lasso_normalize,\n",
    "                                      'L0 Thresh':L0_th,\n",
    "                                      'N':N,\n",
    "                                      'M':M,\n",
    "                                      'x0_min':x0_min,\n",
    "                                      'x0_max':x0_max,\n",
    "                                      'num_bits':num_bits,\n",
    "                                      'd_x0':d_x0,\n",
    "                                      'w_fp':w_fp,\n",
    "                                      'type_x0':type_x0},\n",
    "                                      ignore_index=True)\n",
    "\n",
    "\n",
    "        # Apply L1 regularization\n",
    "        if \"lasso_L1\" in methods_to_run:\n",
    "            x_L1, r_L1, r_L1_test, k_L1, s_L1 = others.lasso_L1(A, b,\n",
    "                                                                A_test, b_test,\n",
    "                                                                lasso_normalize,\n",
    "                                                                L0_th=L0_th, x0=x0)\n",
    "\n",
    "            rec_err_L1, cardinality_err_L1, rmse_L1, \\\n",
    "                cv_rec_err_L1, cv_cardinality_err_L1, cv_rmse_L1 = \\\n",
    "                    get_stats(s_L1, r_L1, r_L1_test)\n",
    "\n",
    "            # Store results\n",
    "            results = results.append({'Experiment':exp_id,\n",
    "                                      'Cardinality': K,\n",
    "                                      'Noise STD': b_noise_coef,\n",
    "                                      'Method':'Lasso',\n",
    "                                      'Infeasible X-Recovery L2 Err':rec_err_L1,\n",
    "                                      'Infeasible X-Recovery L0 Err':cardinality_err_L1,\n",
    "                                      'Infeasible System RMSE':rmse_L1,\n",
    "                                      'CV X-Recovery L2 Err':cv_rec_err_L1,\n",
    "                                      'CV X-Recovery L0 Err':cv_cardinality_err_L1,\n",
    "                                      'CV System RMSE':cv_rmse_L1,\n",
    "                                      'Type A':type_A,\n",
    "                                      'Normalize A':lasso_normalize,\n",
    "                                      'L0 Thresh':L0_th,\n",
    "                                      'N':N,\n",
    "                                      'M':M,\n",
    "                                      'x0_min':x0_min,\n",
    "                                      'x0_max':x0_max,\n",
    "                                      'num_bits':num_bits,\n",
    "                                      'd_x0':d_x0,\n",
    "                                      'w_fp':w_fp,\n",
    "                                      'type_x0':type_x0},\n",
    "                                      ignore_index=True)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "        # %%\n",
    "        # Apply L0 regularization using OMP\n",
    "        if \"omp_L0\" in methods_to_run:\n",
    "            x_omp, r_omp, r_omp_test, k_omp, s_omp = others.omp_L0(A, b,\n",
    "                                                                   A_test, b_test,\n",
    "                                                                   L0_th=L0_th, x0=x0)\n",
    "\n",
    "            rec_err_omp, cardinality_err_omp, rmse_omp, \\\n",
    "                cv_rec_err_omp, cv_cardinality_err_omp, cv_rmse_omp = \\\n",
    "                    get_stats(s_omp, r_omp, r_omp_test)\n",
    "\n",
    "            # Store results\n",
    "            results = results.append({'Experiment':exp_id,\n",
    "                                      'Cardinality': K,\n",
    "                                      'Noise STD': b_noise_coef,\n",
    "                                      'Method':'OMP',\n",
    "                                      'Infeasible X-Recovery L2 Err':rec_err_omp,\n",
    "                                      'Infeasible X-Recovery L0 Err':cardinality_err_omp,\n",
    "                                      'Infeasible System RMSE':rmse_omp,\n",
    "                                      'CV X-Recovery L2 Err':cv_rec_err_omp,\n",
    "                                      'CV X-Recovery L0 Err':cv_cardinality_err_omp,\n",
    "                                      'CV System RMSE':cv_rmse_omp,\n",
    "                                      'Type A':type_A,\n",
    "                                      'Normalize A':lasso_normalize,\n",
    "                                      'L0 Thresh':L0_th,\n",
    "                                      'N':N,\n",
    "                                      'M':M,\n",
    "                                      'x0_min':x0_min,\n",
    "                                      'x0_max':x0_max,\n",
    "                                      'num_bits':num_bits,\n",
    "                                      'd_x0':d_x0,\n",
    "                                      'w_fp':w_fp,\n",
    "                                      'type_x0':type_x0},\n",
    "                                      ignore_index=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEICAYAAABF82P+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABBg0lEQVR4nO3dd3wUZf7A8c83u6n0AEoJvam00FEUQZGmB8op6Hme7c7efjZsp3ieHtj1bIechfMEFMVyIiInyCki0iwI0qSEnoSetuX5/TGbZJNskt3N7mZ3832/XpidmWdmnsdNnu/M88w8jxhjUEoppcpLqO0MKKWUik4aIJRSSvmkAUIppZRPGiCUUkr5pAFCKaWUT/bazkCoNGvWzLRv3762s6GUUjFl1apV2caY5r62xU2AaN++PStXrqztbCilVEwRke2VbdMmJqWUUj5pgFBKKeWTBgillFI+xU0fhFKq9jgcDrKysigoKKjtrKhKpKSkkJGRQWJiot/7aIBQStVYVlYWDRo0oH379ohIbWdHlWOMIScnh6ysLDp06OD3ftrEpJSqsYKCApo2barBIUqJCE2bNg34Dk8DhFIqJDQ4RLdgvh8NEEopv+197DH2PvZYbWdDeXHs2YNjz56wHFsDhFLKb4XrN1C4fkNtZ6NaIsJll11Wsux0OmnevDnnnXdelfutXbuW+fPnlyxPmTKFJ598Muh81HR/f7jzC3Dnh+fhAA0Q6FWRCi/9/Yq8evXq8dNPP5Gfnw/A559/TuvWravdr3yAqOs0QBA7V0XViaeKKJ7KEi+/X7FmzJgxfPLJJwDMmjWLSy65pGTb8ePHueqqqxgwYAB9+vThww8/pKioiAcffJA5c+aQmZnJnDlzAPj5558ZNmwYHTt25Pnnny85xtNPP02PHj3o0aMHzz77bMn6Rx99lG7dujFixAh++eWXyBQ2TDRAxJF4qojiqSyqdlx88cXMnj2bgoICfvjhBwYNGlSy7dFHH+Wss87iu+++Y/Hixdx11104HA7+8pe/MGnSJNauXcukSZMA2LBhA5999hkrVqzg4YcfxuFwsGrVKl5//XW+/fZbli9fzquvvsqaNWtYtWoVs2fPZs2aNbz//vt89913tVX8kND3IJRScalXr15s27aNWbNmMXbs2DLbFi5cyEcffVTSP1BQUMCOHTt8Hufcc88lOTmZ5ORkTjjhBPbt28dXX33FBRdcQL169QCYMGEC//vf/3C73VxwwQWkpaUBMG7cuDCWMPw0QCil/OLMzaX5HbcDUJS1i4S0VOzp6bWcq6qNGzeOO++8kyVLlpCTk1Oy3hjDe++9R7du3cqk//bbbyscIzk5ueSzzWbD6XRijKn0nJF83Nc4ndhbnAiAu6gISUhA7KGr1rWJSakwKq5Um99xO0VZu3Dm5tZ2loLmzstn+8WXsP3iS9gyYgTuvPzazlK1rrrqKh588EF69uxZZv2oUaP4+9//XlLRr1mzBoAGDRpw9OjRao87dOhQPvjgA/Ly8jh+/Djz5s3jjDPOYOjQocybN4/8/HyOHj3Kxx9/HPpCeTFuN0Vbt1K0dSuFGzdi3O6QHr/O30HE4lWRih3FlWqxTosWgf56RUxGRga33nprhfV//vOfue222+jVqxfGGNq3b89//vMfhg8fztSpU8nMzOTee++t9Lh9+/bliiuuYODAgQD88Y9/pE+fPgBMmjSJzMxM2rVrxxlnnBGegkWIVHWrFEv69+9vgpkwqChrF1tGjChZ7rRoEUkZ1T8OF22cubkUedpQ7c2ax3Sgi6eyhOP3yxgDLpd1teh0lv50uTAul7XN66dxOsHtxjhd4HJiXG7PT690The4rZ/GVZo+pfspxSfFeeAAO//4p5J8tJnxKvZmzQDYWlTESV264Dp2DAx4/lMu45Wsr2pbAOtNFdtCfv6q8lBZneprfRDntZ94IiR4Gn+czpK/FYCkDh0Qr8H4yjc5rV+/npNPPrnM8URklTGmv89z+cyBijnxdKUaTWUxLhemoAB3fj7uggLceXnWcl4+7oJ8TH6+50Wl0vX1zxwKNptVqWZnlzle4caNOLMPlFQW7uPH2T9tmlUxu13gLFe5e1f4xcHA5YpY+dvOfJMdf7jc5zbvYOF88QVMp0449+2LVNbqLhGKtm71uano11/LLCd37UpNekTqZIBw5uaWtJ86du4ss82RtRPjKLIWRJCkJHC7kcRExG4v/ZeYCImJOv5MLTMOh6fizscU5FsVuffn4so7v7iSz8fkeSr7StcXV/z5mMLCgPOUNmhgpZVq1g03lFluO/NNCjdtDqrsSoVbnQwQ7rz8Mrf93nZccWWZ5aquoACw2coEDxLtnuXEkkBS7bZEO9h9bCteb08sk7b4c3L37iVXB84DB8pky7FzB6bQ6/X7hASKtm2z7lqN27oSdRswbqs5ofxnt7tCuoA+u90YU/rZ2lbJZ2NoNG68VV7A6fW0CUDhpo1W+Yx1DtfRo+y+627cBQXgcPjzlcefhATEZrN+/2w26/ckIQHsNsRW/NnuSZNgrSuT3oYk2BC7DapIn+B5XNMfkpBQ0uQEgu9L18rWV9wmlaz361iVXbgFuL7KC8CAjuV/ucuvr7DJZqs8TyFWJwNESBU3BdTCqasKXjuuvKpC2qzrb/CZNho0HDOG7V5vunorn++2M9/E7ceTJqEiqakkpKSQkJpqffYsS1oqCSmpnvUpJKSmkZCagq1hQ7+PbW/enA4ffojYbWUrcJsdsSVYy/byFbgtYneuztxcq4kPcO7Zw3av8Y3a/etf2Fu2BGDzwVwkMZHEFi0ikq+6zDidJHftan12OMo0K/nqg6gJDRDVEJuNxFatMA4Hxuks/edwgNNZ29lTCQlWBe2rsk5JISEtFfGsT0hN8VT2qWXXp6UiKaUVvHcgkOTkgP/I/K1UgajvgLenp1fa/2Nv2bKkw12OHgnp8/eqcmK3l9xVlH+oVRITSUhKCtm56uQ3mpCWGtAfcOcv/uvzOMYYq+PQO3g4ioOHoySQlKx3OkoCS8k2R3HQ8drm8D6mo/QcjrIBKiE11f8yJydT78yhiCRYT0AkiH+fExJAyn9OQBLE5+dq90d87mNr1NjvsiSe2IKu3y5HUlOt5rYo6wfyt1JVKtrVyQARqj9gEbE6qgOY4zWUAr1SbfuPf0Q8j/4KtCy2Ro0inkcV3Ww2W5kX4i6++GLuueeeStPXr1+fY8eOhS0/zz77LNdcc03JsBtjx47l7bffpnHjxjU6bvv27Vm5ciXNSvp7rIvVO6dOZeHy5aTVq8cbb7xB3759a3QeqKMBIl7E05VqPJUlXiWkpdJu9iyg9P2UYH2wZhdPfPYLuw/l06pxKneN6sb5fWr2HaemprJ27doaHSOUnn32WX7/+9+XBIhwDCMuCQkkdezI/AUL2JqdzcYNG1ixahXXX3+9z2FDAqVDbSgVRsWVarvZs+i0aFGNKtXaZk9P58BTT3PgqadJymgddN/JB2t2ce/7P7LrUD4G2HUon3vf/5EP1uwKbYaBw4cP061bt5Jhty+55BJeffXVku33338/vXv3ZvDgwezzvMPx8ccfM2jQIPr06cOIESNK1pef/KdHjx5s27aN48ePc+6559K7d2969OjBnDlzeP7559m9ezfDhw9n+PDhgHXln52dzeTJk3nppZdKjjNlyhSeeuopAJ544gkGDBhAr169eOihh6otn9jtOPfu48N353L5FVeQkJjI4MGDOXToEHtCMMtcnb+DCOVVkVLl2dPT2XXrbQC0+9fM2s1MhDz88Tp+3n2k0u1rdhyiyFW2ezXf4eLuuT8wa4XvEVVPadWQh37Tvcrz5ufnk5mZWbJ87733MmnSJF544QWuuOIKbr31Vg4ePMif/mS94Hf8+HEGDx7Mo48+yt13382rr77KAw88wOmnn87y5csREWbMmMHjjz9eUoH7smDBAlq1alUy98Thw4dp1KgRTz/9NIsXLy7TFARW09dtt93GDZ53Yt555x0WLFjAwoUL2bRpEytWrMAYw7hx41i6dClDhw6tstwJqSnsyT5AmzZtStZlZGSwa9cuWno1zQajzgeIePkDjqdAF09lURWVDw7VrfdXZU1M55xzDu+++y433ngj33//fcn6pKSkkilI+/Xrx+effw5AVlYWkyZNYs+ePRQVFdGhQ4cqz9uzZ0/uvPNOJk+ezHnnnVft+Et9+vRh//797N69mwMHDtCkSRPatm3L888/z8KFC0vGdDp27BibNm2qNkAktmwJPp5cCsXDG3U+QMSLeAl0EF9lAUg++aTazkLI+FOW6q70h0z9gl2HKo4E27pxKnOuPTXovFXG7Xazfv16UlNTyc3NJSMjA4BEryfgiofxBrj55pu5/fbbGTduHEuWLGHKlCkA2O123F6jpRYUWC+idu3alVWrVjF//nzuvfdeRo4cyYMPPlhlni688ELmzp3L3r17ufjiiwGro/nee+/l2muvDbiMGRkZ7PQaFSIrK4tWrVoFfJzytA9CqTBrcd99tLjvvtrORkiEoix3jepGamLZt4FTE23cNapbJXvUzDPPPMPJJ5/MrFmzuOqqq3BU8/b94cOHS+avfvPNN0vWt2/fntWrVwOwevVqfvW8oLZ7927S0tL4/e9/z5133lmSpqqhw4tnu5s7dy4XXnghYA1B/tprr5U8WbVr1y7279/vVxnHjRvHzJkzMcawfPlyGjVqVOPmJdA7CBWl4umqW5VV/LRSqJ9iKt8HMXr0aK666ipmzJjBihUraNCgAUOHDuWvf/0rDz/8cKXHmTJlChdddBGtW7dm8ODBJYHgt7/9LTNnziQzM5MBAwbQ1fM2848//shdd91FQkICiYmJvPzyywBcc801jBkzhpYtW7J48eIy5+jevTtHjx6ldevWJRX5yJEjWb9+Paeeat1F1a9fn7feeosTTjihQh579epFgucFzokTJ/LUU08xf/58OnfuTFpaGq+//nqQ/xfLqvPDfQPsfewxgJi/ytt+2R+A+GiWUbHF1zDSKvrocN9BiPXAUEyvupVSoVQrAUJERgPPATZghjFmarntjYC3gLZYeXzSGBOae6Y4Fi+BTikVHSLeSS0iNuBFYAxwCnCJiJxSLtmNwM/GmN7AMOApEQndCFRKKaWqVRtPMQ0ENhtjthpjioDZwPhyaQzQQKxn0OoDuYAOnaqUUhFUGwGiNeA9jVuWZ523F4CTgd3Aj8CtxpgKb9GIyDUislJEVh4oN2GOUkqpmqmNAOHr9b7yj1KNAtYCrYBM4AURqTALizFmujGmvzGmf/PmzUOdT6WUqtNqI0BkAW28ljOw7hS8XQm8byybgV8BfURHKVWprKwsxo8fT5cuXejYsSM33XQThZ45xd944w1uuummMumHDRtG8aPx7du3p2fPnmRmZtKzZ08+/PDDknTr1q3jrLPOomvXrnTp0oVHHnkEX68HLFmypGToDm+//vorgwYNokuXLkyaNImioqJQFjusaiNAfAd0EZEOno7ni4GPyqXZAZwNICInAt2ArRHNpVIqfH5dCi8NhqP7yn4OkjGGCRMmcP7557Np0yY2bdpEfn4+d999t9/HWLx4MWvXrmXu3LnccsstgPXy3bhx47jnnnvYuHEj33//PcuWLSszGmt1Jk+ezP/93/+xadMmmjRpwj//+c+Ay1dbIh4gjDFO4CbgM2A98I4xZp2IXCci13mSPQKcJiI/Av8FJhtjsiOdV6VUGPy6FN6eCAc2wXt/LP385bSgD/nFF1+QkpLClVdeCVhjKz3zzDPMnDkz4EmBjhw5QpMmTQB4++23GTJkCCNHjgQgLS2NF154galTp1Z1iBLGGL744ouS4TQuv/xyPvjgg4DyU5tq5T0IY8x8YH65da94fd4NjIx0vpRSIfDpPbD3x8q3714NDs9gfdu/guLnT1a/CQd+8b1Pi54wpvJKed26dfTr16/MuoYNG9K+fXs2b97sV7aHDx+OMYatW7fyzjvvVHrcTp06cezYMY4cOULDhhW6RsvIycmhcePG2D3zdRcPwx0r9E1qpVRkndADcn6BwqOlwUESIL1z0Ic0xvgc3rq4r6Cyoa+91xfP3bBlyxbOPvtshg0bVulxqzqmr/MHul+00AChlAqtKq70gdImJu8n1xPs0H4InPd0UKfs3r077733Xpl1R44cYd++fXTr1o2srCwOHjxYZntubm6FyXzAukM48cQT+fnnn+nevTtLly4ts33r1q3Ur1+fBg0aVJuvZs2acejQIZxOJ3a7PWTDcEeKDvetlIqsT+8Gp2fIbXsq2JLAVQQ/fxD0Ic8++2zy8vKYOdMaqNLlcnHHHXdw0003kZqayoABA/j666/Zu3cvACtXrqSwsLDMLGzF9u/fz6+//kq7du249NJL+eqrr1i0aBFgdVrfcsstfnd+iwjDhw9n7ty5gDV8+Pjx5d8Ljl4aIJRSkXXZh9DvckhrChe8An0usz5f9EbQhxQR5s2bx9y5c+nSpQtNmzYlISGB+++/H4ATTzyR5557jrFjx5KZmcltt93GrFmzSobMBqsPIjMzk+HDhzN16lROPPFEUlNT+fDDD/nrX/9Kt27d6NmzJwMGDKjwyGyx//73v2RkZJT8++abb5g2bRpPP/00nTt3Jicnh6uvvjrockaaDvetlKqxaBvue9myZVxyySW8//77FTqZ67KIDfctIs2NMQc8n08yxmwI9lhKKRVKp512Gtu3b6/tbMS8oJqYROSfwF4R+UBE7MDC0GZLKaVUbQu2D2IC1pvOh4A38D2+klJKqRgWbIA4AnwFXA00xxpUTymlVBwJNkA8AGQYY1zARcCUkOVIKaVUVKi2k1pE/lDJpqEiMtTzebuI/MEYMzN0WVNKKVWb/HmK6Uo/j2UADRBKqVpRv379gAfmU1WrNkAYY4ZHIiNKKaWii999ECKS6GNd1UMZKqWUD8PmDKPnmz1L/g2bMyws5/n4448ZNGgQffr0YcSIEezbZ8058eWXX5KZmUlmZiZ9+vTh6NGj7Nmzh6FDh5KZmUmPHj343//+B8CsWbPo2bMnPXr0YPLkyWHJZ7QK5EW51SJyuTFmNYCIjAGmU3Z2OKVUHTdtxTQ25Fb93mxOQU6F5SsXVN6afVL6SUweGHjlfPrpp7N8+XJEhBkzZvD444/z1FNP8eSTT/Liiy8yZMgQjh07RkpKCtOnT2fUqFHcf//9uFwu8vLy2L17N5MnT2bVqlU0adKEkSNH8sEHH3D++ecHnJdYFEiAaAosF5FpQEvgKkDfnlZKRa2srCwmTZrEnj17KCoqokOHDgAMGTKE22+/nUsvvZQJEyaQkZHBgAEDuOqqq3A4HJx//vlkZmbyxRdfMGzYMIrnvL/00ktZunSpBggfTgJmAPd5ll8E7gx5jpRSMc2fK/2eb/assO710a+HPC8333wzt99+O+PGjWPJkiVMmTIFgHvuuYdzzz2X+fPnM3jwYBYtWsTQoUNZunQpn3zyCZdddhl33XVXtRMCxbtA3oO4GDgHOA44gQuA0eHIlFIqvjVNaVrlcqgcPnyY1q1bA9ZQ28W2bNlCz549mTx5Mv3792fDhg1s376dE044gT/96U9cffXVrF69mkGDBvHll1+SnZ2Ny+Vi1qxZnHnmmWHJazQK5A7iFWAR8EcgHWuIjfcBW+izpZSKZ0smLQn5MfPy8sjIyChZvv3225kyZQoXXXQRrVu3ZvDgwfz6668APPvssyxevBibzcYpp5zCmDFjmD17Nk888QSJiYnUr1+fmTNn0rJlS/72t7+VTEc6duzYmJrPoab8Hu5bRK43xrzstWwH/myMeShcmQuEDvetVO2JtuG+lW9hG+7bOzh4lp1AVAQHpZRSoVeToTYq0KE2lFIqfuhQG0oppXzyJ0AsBeYDK0y8zE+qlFKqWv485toU+DewX0TeEpFLRSQ8z6QppZSKGtUGCGPMTcaYzsAQYAXwe6zhvZeLyEMiMjDcmVRKKRV5fr8oZ4zZaIx53hgzBmgGPIx1d/FWuDKnlFKByMrKYvz48XTp0oVOnTpx6623UlRUxJIlSxAR/vnPf5akXbNmDSLCk08+CcAVV1xBhw4dyMzMpG/fvnzzzTe1VYyoEdSMcsaYAmClMeYWY0xXETkpxPlSSsW5jaefwfqTTmbj6WeE5HjGGCZMmMD555/Ppk2b2LhxI8eOHeP+++8HoGfPnsyZM6ck/ezZs+ndu3eZYzzxxBOsXbuWqVOncu2114YkX7EsqAAhIv8E9orIB54X5haGNltKqXjnys4u87OmvvjiC1JSUrjySuvBS5vNxjPPPMNrr71GXl4ebdu2paCggH379mGMYcGCBYwZM8bnsYYOHcrmzZtDkq9YFshQG94mAGcDV2ANuSEhyo9SKsbtfewxCtdXPtBz3po14HSCCBgDIqw/6WSw20nr08fnPsknn0SL++7zua3YunXr6NevX5l1DRs2pG3btiWV/YUXXsi7775Lnz596Nu3L8nJyT6P9fHHH9OzZ8UBBeuaoO4ggCPAV8DVQHOgVchypJSKb06n9bP4qfnin8Xrg2SMQaTitar3+okTJ/Luu+8ya9YsLrnkkgpp77rrLjIzM5k+fXqZ/oq6Ktg7iAeADGPMNhG5CLg1hHlSSsWw6q70N55+Bq7sbCQpCVNUVPLT1qwZ7f4V/Lu23bt357333iuz7siRI+zcuZNOnTpZeWvRgsTERD7//HOee+45li1bVib9E088wYUXXhh0HuJNsJ3U/zLGbPN8PmKMeSSkuVJKxa2uX/2PkzesxxQVAWCKijh5w3q6fvW/Gh337LPPJi8vj5kzrSDjcrm44447uOKKK0hLSytJ95e//IVp06Zhs+lA1NXRsZiUUrXC1qwZruxsbM2aheR4IsK8efO44YYbeOSRR3C73YwdO5bHHnuszCOrp512WkjOVxdUO9y3iCz281jGGHOWXycVGQ08hzWXxAxjzFQfaYYBzwKJQLYxpspZOnS4b6Vqjw73HRtCPty3MWZ4iPJWnBkb1nSl5wBZwHci8pEx5mevNI2Bl4DRxpgdInJCKPOglFKqetX2QRSPuyQiiT62BTNh60BgszFmqzGmCJgNlJ+i6XfA+8aYHQDGmP1BnEcppVQNVBkgRORKYIFncbWI9PXaNgZYF8Q5WwM7vZazPOu8dQWaiMgSEVlVWT+IiFwjIitFZOWBAweCyIpSKlR0sOfoFsz3U10T04tAD8/npsByEZkGtASuAip/G6Zyvl6qK59zO9AP62W8VOAbEVlujNlYZidjpgPTweqDCCIvSqkQSElJIScnh6ZNm/p8F0HVLmMMOTk5pKSkBLRfdQFiMlYT0EDgJGAGUPyQ84vAnQHmE6w7hjZeyxnAbh9pso0xx4HjIrIU6A1sRCkVdTIyMsjKykLv5KNXSkoKGRkZAe1TZYAwxvxdROZ6Fi/G6lg+DiQDFwCLgA8DzOd3QBcR6QDs8hz3d+XSfAi84BnnKQkYBDwT4HmUUhGSmJhIhw4dajsbKsT8mQ9ij+fjK1iVew+sCjsbeD/QExpjnMBNwGfAeuAdY8w6EblORK7zpFmP1ffxA9YcFDOMMT8Fei6llIp301ZMY9qKaWE5diBDbdxojHnZ83mHiPQH/hzMSY0x87GmMfVe90q55SeAJ4I5vlJK1RUbcoPpCvaP3wHCKzgULzuBh0KeI6WUUlEh2NFclVJKxbmAA4SI/EFEGoQjM0qp6BbO9m4VfYIZ7vt1YDlwNMR5UapEcSU0eeDkWs6J8hbO9m4VfYIJEPoWjAq7eKqINNipWBXshEFKKT/FU7BTdYt2UscRbR9WSoWS3kEQP00AeqWqlAqlYALELsAR6ozUJq1YlVKqooADhDGmTfWplFJKxTrtg1BKKeWTBgillFI++R0gQjjlqFJKqRgQSB/EahG53BizGkqmHJ1O2cl/lFJKhdnBgoPkOfMAuK3fbQDsOrarZHuaPY0mKU1qfJ5AAkSophxVSilVA3nOPEa/N7rS7Qt+u4AmRDZAhGrKUaWUqvOMMTjdTgpcBRS6Cq1/zkIKXAUUuYqs9c7Ckm0l650FDGk9JCJ5DCRAhGrKUaVUDIlUc0Ztchs3BU6vitmPCruk4nYWlF1XrlIvdBZS6C57rOLtbuMOKr+9mvcK8f8B3wIJEK8A/wWuBtKBN7CmHLWFPltKqWgRqeYMqPqqOtCr7Koq7fJpHO64evc3ZGplylGl4l08XHW73C6/rnLzHfnM2zSv8ivqyip1H1fVBc4CDCZCJax99gQ7ybZkkm3JpNhSSLIlkWJPKVmXbEsmxe5Zbytd3yApMlPyVBsgROQPno/HvT4X2+K9zhgzM5SZUypWhfKq223cFa6Cq6p8i6+Oi9w+rqSrSl+uecXpdgLw2qjXqszfwcKDPLjsQf//50SpFFsKyfZkkhOSrZ/lKunqKu1ku1clX3wsW+XHSbIlYU8Ibjg874uNcPInd1f6eSwDaIBQQYv2q+7iK2rvirWyZoxOjTtVeayDBQd5auVTVTZ/FJ+nrjV/2BPslVa0lV1l+5Omqko9MSEREZ3qprxqA4QxZngkMqKCE+2VaiD8uepubBrjdDvLVKaVVdLlK/Piq+ky68r9rOpKvfiK2h/VXXXnO/P5fPvnfh+vNghCij2FBKn6fdpUeyrjO40P6qq6QkVvS8aWoN2a1Umzp7HgtwsAyM7PBqBZarMy20NBh/uOcZHsQCzPu+mjuk5FX52D5Zf/cEr5Fsyy9hzbw9j3xwb95EcsK26SKF/5+rpqripddVfT3sezJ9gRkWqbM5qkNOGvp/81Qv8nFFj/z4v/rh/46gEAXh/9esjPU2cDRDxdeVelyFXEij0rKnQQBlOBlw8AoW76mNhtYpXbDabWg4P31XCSLclnBZxkSyLVnlrlcRonN+aJM5+oWEEnlLvqtqeQlJCkzR+qVtTZABGqK29jDE7jxOFylGnWKHQVUuQuKv3sKvvZ+/E67/3KpylyFVHk9rHO8/m54c9Vmb/s/GyuXnh1wP9/opld7D4r5ZJKu7iiraTyLn+1Xf7qurIr70Daqau76k5LTGN0+8p//6JJpJozVPTxK0CI9VeRYYzZGeb8RI3c/FweWf5ISUVcPgB4V/61eVVb248EBtK+XN1yekp6ledqWa8lay5bE/STHyo4kWrOUNHHr780Y4wRkQ+AfuHNTvQocBXw9a6vazsbNZaYkEj/E/uXNFsEUmFX2CfMTR/VXXWLSMwEB73qVvEgkL+25SIywBjzXdhyE6PsYifJllTyz7tJwrtpo/z2JFsSSQlJJZ+905VZl1DuuLbS4x4tOlpl3pqnNdervVqgV90qHgQSIIYD14rIdqzxmATr5iIyg4JEWHpKOi+d/VK1lXZNXnYJhWOOY7V27lDTq26lArct53jYjh1IH8R1wPaw5STKpNhTOCPjjNrORrXiqVLVq26lApdX5ArbsQPpg3jGGFNn+iBihVaqSqlwqbN9EPF05a2UUuEQaB/EdSKyjTjog9ArbxUp4WwjjrR4KouqXiABYkyoTioio4HnsOaSmGGMmVpJugHAcmCSMWZuqM6vVCSFs4040uKpLKp6VY/CBYjI3QDGmO3AQGPM9uJ/wLWBnlBEbFjTlY4BTgEuEZFTKkk3Dfgs0HMopZSquWoDBNZUo8XuLbctmLECBgKbjTFbjTFFwGxgvI90NwPvAfuDOIdSStUJw1qPYlzHC8JybH8ChFTy2deyP1oD3kN2ZHnWlR5UpDXWnNevVJkxkWtEZKWIrDxw4EAQWVFKqdh0ON/Bztw8WsiZHMnuxctLtpB1MI/D+aEbRNOfPghTyWdfy/7wFVTKH+dZYLIxxlXVUA7GmOnAdID+/fvXnXkKlVJ1Ws7xQh7/9BfeWbUT41XzPf7ZBib2a8PdY7rRtF5yjc/jT4DoLSJHsCr2VM9nPMspQZwzC2jjtZwB7C6Xpj8w2xMcmgFjRcRpjPkgiPMppVTcOJzv4PFPf2HOyopjpxpDyfr7zj2ZRqmJNTqXPzPKhXp6p++ALiLSAdiF1cfxu3Ln7FD8WUTeAP6jwUGp2jes9SgaJus7QqFmjKHI5Sa/yEWe55/12Umew1Wyvsjpon+7dN5ZVfXA2u+s2snNZ3cOf4AINWOMU0Ruwno6yQa8ZoxZJyLXebZX2e+glIq8w/kOjuQ7aCFnsje7kJeXbOE3vVvSICWxxpVQrDDGUOBwW5V2kYt8R3Fl7iypwH1V6tb68vt4rStykedw4XJX30o+vNsJ7D5cUKZZyXde4T8/7OG6M6ueG706tTLKnDFmPjC/3DqfgcEYc0Uk8qRUuKS421SfKIpFqr07FNxuU6biLn817quC9q6kK67zqvwdrmor5nBrnJbIgaOFfqXdmZNX4/PFxuD6qs45Kf2k2s5CyIxpdS31k2PzTy0c7d1Ol7vcFbaz7NW2w+vKunxl7vCq7MvsY60rcMTufOWJNiE10UZakp20JBupSTbPTztpidbnzDaNOZhX5Nfx2jSteVNgjX5rRaSFMWZvjXOhVDmTB06u7SzUWHGzTLI9gZ0H82OqWcblNuQVOTmc5/CrvfvK09vzzOe/cKSgXHNLuco+v8hFkSt2K/EkewJpSTbSEosrcHtJRV7P63Nqko20xLIVfUlln2TzBIKy+yfa/HnrALIO5vHsfzdVeTcjAuf1alnj8tb0suafwLk1zoVScSZSzTLGWE0qxwutyvd4kZO8IifHC0uvzo8XucgrdFZYPl7ktU+h9bN4ucDhZni3E+jeuqFf7d2f/LCHbTl5LPml9t9HKq58fV2Bl1bWnorZq6KvUJl7KvjS9XZsCaGbQTFYDVISmdivjc+7umIT+7WhQUrNL0JqFCCMMRocokg8NcvEMn+aZQyGa4d2Yu+RgtKKu9Drp6O00s4rtJpWiiv1krb1QqszNFzt4oG0dx84Wkjj1CS/0iYIZa6cfV1Ne1fQFa7AE8tW2t5pUuw2EqKgEg+nRqmJ3D2mG0CFCxARKzhMHtMtJHepfgcIEZlmjJlc3bpYFC8Vazw0y8SDowXVN8u8uyqLiwe05ao3vqPQGZ1NLgUOFxlNUv1K27pJKqe0bMiQzk19X417KvbUJBvJ9oSQzmVeFzWtl8x9557MzWd3Zu6qLPYdLqBds3qc1yu0TZiB3EGcA5Svgcb4WBdztGJVofTx97v9apZZ/Mt+Tu3UtMbNMimJCSWVcr0kO2nJpVfk9ZJspCVbP1PLLZfsk1yc1rqCr5dceiWedTCPFxZvrra9+4I+rcloou9HRFKjVCsQ5Be5aNUktcaPtPpSbYAQkeuBG4BOIvKD16YGwLKQ50ipGJVX5OTbrTls8/PxwgNHC+nfrgk2Ea9KvLSSL6mwvZaLr8jrJdtLgkA428Uj2d6tgrN25yEAbj6rS8iP7c8dxNvAp8DfgHu81h81xuSGPEdKxZi8Iif/Xr6DfyzdQs/WjeneuqFf+7VvVi8sV32hFMn2bhV9/Blq4zBwWETeB3KNMUdF5AGgr4g8YoxZE/ZcKhWF8otcvLV8O/9YuoXsY9az6cu2ZHPzWZ150Y9mmVA8hhgJkWrvVtEnkD6IPxtj3hWR04FRwJNYw3EPCkvOlIpS+UUu/v3tdl75sjQwFGtaLwm3MXHXLBOJ9m4VfQIJEMVzDZ4LvGyM+VBEpoQ+S0pFp6oCQ8tGKdw4vDMX9c8g2W6jQ/N6QPw1y4SzvVsF55RW/jVpBiOQALFLRP4BjACmiUgy/k04pFRMKw0MW8k+Vva9gJaNUrhheGcmegJDMW2WUZHy0G+6h+3YgQSIiVhTjD5pjDkkIi2Bu8KTLaVqXzCBwZs2y6hY53eAMMbkAe97Le8B9oQjU0rVpqoCQ4uGKdw4vBMTB7SpNDCUp80yKlYF8iZ1MvBboL33fsaYv4Q+W0pFXoHDxb+/3cErX26pMMREMIFBqVgXSBPTh8BhYBXg3wAtSsUADQxK+RZIgMgwxowOW06U8vLwx+uA8HbAFThcvP3tDl6uJDDcMLwTE/u3ISVRA0OxcD4xo6JPIAFimYj0NMb8GLbcKOXx8+4jYTt2VYHhxIbJ3Di8swaGSoQzYKvoE0iAOB24QkR+xWpiEsAYY3qFJWdKhViBw8WsFTt4eckW9vsIDDcM68ykARoYlCoWSIAYE7ZcKBVGGhiUCk4gj7luF5HewBmeVf8zxnwfnmwpVXMFDhezV+zgJR+B4YQGydwwrBMXD2yrgUGpSgTymOutwJ8ofRfiLRGZboz5e1hyplSQigPDy19uYd+R2g8M2rGrYlUgTUxXA4OMMcfBmk0O+AbQAKGiQoHDxZzvdvLSks0+A8P1wzpxSS3cMWjHropVgQQIoXTAPjyfdd5AVeuiNTAoFesCCRCvA9+KyDzP8vnAayHPkVJ+qiowNPc0JWlgUCp4gXRSPy0iS7AedxXgSp0sSIVLn7aNqZ/s+9ezwOHinZU7eWnxFvYeKSizrXmDZK4/sxO/G6SBQamaCqST+k3gVmPMas9yExF5zRhzVdhyp+qcw/kOjuQ7SLYnsPNgPi8v2cJveltDZCfbEzQwKBVBgTQx9TLGHCpeMMYcFJE+oc+Sqqtyjhfy+Ke/VJhk5/HPNvDbvq25oE8Gj36ynkKnu2Rbs/pWH8OlGhiUCrlAAkSCiDQxxhwEEJH0APdXYRaJ8YvC5XC+g8c//cXnNJ3GwNxVu3C54bYRXZi24JeSwPC7gW1JTdLAoFQ4BFLBP4U1HtNcwGBNIPRoWHKlghLO8YvC7WiBg3dWVT6HM8AHa3cx7/ohPPQbGxcP0MCgVLgF0kk9U0RWAmdhdVJPMMb8HLacqTrjwNEC5q3eVaZZyRdjYPmvOTozm1IRElATkScgaFBQNXLgaCHLtmTzzZYcvt6STf926STZ/ZvefGdOXphzp5QqFshTTAJcCnQ0xvxFRNoCLYwxK8KWOxUXjhQ4WLE1l6+3ZLNscw6/7DtaZnvn5g66t/ZvOIo2TdPCkUWllA+B3EG8BLixmpj+AhwF3gMGBHpSERkNPAfYgBnGmKnltl8KTPYsHgOu14EBY0eBw8Xq7Qf5eks2X2/O4cddh3G5K28/WrPjIHeM7MqLizdX2cwkAuf1ahmGHCulfAkkQAwyxvQVkTVQ8phrUqAnFBEb8CJwDpAFfCciH5Xrz/gVONNzjjHAdGBQoOdSkeF0uflx12GWbclh2ZZsVm47WOZR1PISbUKfNk04rXNTTuvUjMw2jcl3uJjYr43Pp5iKTezXhgYpieEoglLKh0AChMNTuRsAEWmOdUcRqIHAZmPMVs9xZgPj8erbMMYs80q/HMgI4jwqTIwxbNp/jK83W3cI327N4Wihs9L0InBKy4YM6dyM0zo1ZWCHdNKSyv7qJdkTuHtMN4AK70GIWMFh8phuNErVAKFUpAQSIJ4H5gEniMijwIXAA0GcszXgfZmYRdV3B1cDnwZxHhVCO3PzWOZpMlq2JYfsY4VVpu/YrB6ndW7KkE7NGNyxKU3qVX+z2bReMvedezI3n92Zuauy2He4gHbN6nFeL+tNag0OSvnw6T3WzzFTq04XhEAec/23iKwCzvasugjoGcQ5fY0A67PlWUSGYwWI0yvZfg1wDUDbtm2DyIqqTPaxQqvJaHM2y7bksCO36qeHWjRMKWkyGtK5KS0bpQZ13kapViDIL3LRqkmqPtKqVHX2/hi2Q1cbIESkIXAj1pX/R1id1TcBHwNrgX8HeM4soI3Xcgaw28d5ewEzgDHGmBxfBzLGTMfqn6B///7VPEUf/6oa4K46RwscfLs1t6QfYcPeo1Wmb5SayKkdmzKkc1NO69yMjs3qYT3oFhprdx4C4OazuoTsmEqpwPhTm/wLOIg1OdAfgbuAJGC8MWZtEOf8DugiIh2AXcDFwO+8E3geoX0fuMwYszGIc9QpVQ1wV1mzTIHDxeodB1m22XoX4Yesqp80Sk20MaBDOkM6NWVI52ac3LIhtgSdDkSpeOZPgOhojOkJICIzgGygrTGm6kvMShhjnCJyE/AZ1mOurxlj1onIdZ7trwAPAk2BlzxXpU5jTP9gzhfvqhrgbmK/Ntw9phtN6yXjcht+3HWYrzdn+/WkkT1B6NO2Mad1sjqW+7Rt4vfLbEqp+OBPgHAUfzDGuETk12CDg9dx5gPzy617xevzH7HuVlQVqhvgbs7KnbiMm6FdmnP/Bz9xtKDyJ42g+Ekjq8loYPt06gXZXKWUig/+1AC9RaR4FDgBUj3LAhhjjM7IXkv8GeDuvdW7+N3AdhT5uFvo0Kwep3majAZ3bEq6H08aKaXqjmoDhDFGh8yMUh9/v8evAe4W/7KfUzs1Zf2eIwzp1IxTPUGhVePgnjRSStUN2oYQw3bkHvcrXfaxQqZO6MWJDZND+qSRUqoWHc+BlEZl17mcUHAIUhqDrebVu/Y6xrDWft4BtGtajxaNUmIqOJzSqiGntNLWS6V8OrwL/v1b2PgZnHoj1GsOebmwcQFMH2atz8ut8Wn0DiIGHS908viCDYzPbI3Ixrgc4C4WZ8VTKiKO58CcS2H3GpjzO+gyEobeCR9cbwUIsNYPug6G3QupjYM+lQaIGLNyWy53vPs923PyaNEohfMzWzNvza5K0+sAd0rFmZRGcMZdVhAA2LTQ+uctvSOcdkuNggNoE1PMKHS6+Nun67noH9+w3TNpzrOLNvGHU9sxsX8G5VuPRGBSfx3gTqm4k38QGre17hx8Se8Il/8HGrWu8an0DiIGrNt9mNvnfF9mop0GyXamjOtOZpvGdGxen1vO7qID3CkVj4yBPd9bdwkbP4Ndq6BRBlzwSsU7B4AJ00MSHEADRFRzuty88uUWnl20CafXMBind27G4xf2KnlMVQe4UyrOFB6FrUusgLDpczi2t+z2YZPhq2d87/v+NXoHEe+2HDjGHe98XzJoHUBKYgL3jT2Z3w9qR4KPcZB0gDulYljOFk9A+Ay2LwNXke90XUZCgh02L/K9PXcrvHleSIKEBogo43Yb3vxmG9MWbKDAUfr2c5+2jXnqot50bF6/FnOnlAoZZxHsWAYbF1pBIWdz5WlTm0Dnc6DrKGh/Bsw4u3Rbekc49yn45A4rOID18+ObYcIMSEsPOosaIKLIrkP53PXu9yzbUjq6eaJNuG1EV64d2hG7TZ8pUCqmHd3neeroM9iyBIqqGNbuxJ7QdSR0GQUZ/SHBM6iFywmjp1lPMaV3hIvegB/fs+4Y3jzPCg7pHeE3f69RcAANEFHBGMPcVVn85eOfy0zdeVKLBjw9MVNfGFMqVrnd1vsKmz6zmo/2rK08bWIadDjTExRGWh3Rvtjs0O5U6z2H026BRVPgyG6rOeny/1h3Dr/5u/ZBxIMDRwu59/0fWbR+X8m6BIFrz+zEbSO6kGzXobCUiikFh2HLF1bT0ebP4fiBytM2bmc1G3UZBe1Ph8QU/86Rll76EtwRr/nWGrWucbOSNw0QtejTH/dw/wc/kXu8tDOqfdM0nprYm37tQvMFK6XCzBjI3lR6l7DjG3BXMrR+gh3anmrdIXQdBc26UuElJn9V9hJciIIDaICoFYfzHDz00U98sLbsTKt/OLUd94w5ibQk/VqUimqOAtj+VWkH88FtladNa+YJCCOh01kVB9iLYloTRdjSjQe4e+4P7D1SULKuZaMUHr+wF2d0aV6jY2tfhVJhdGS35zHUhdY7Co68ytO27G01G3UdBa36QkJsPmCiASJC8oqcPDZ/PW8t31Fm/YS+rXnoN91D8sazDnCnVAi5XdZby8XvJuz9sfK0SfWh4zBPf8JIaNAiYtkMJw0QEeA9wF6xpvWSePSCnozuER+/SErFhfyDsPm/1l3C5kWQl1N52vROpQGh3WlgT45cPiNEA0QYFTpdPPP5JqYv3YLXSBmMPOVEHpvQk2b14+8XSqmYYgzsX+/pYF4IO78F4/KdNiER2g8pbTpqGiVD2rToGbZDa4AIk3W7D3PHO9+zYW/FAfYm9G0dU5P3KBVXHPnw69LScY4O76g8bf0Tocs50HW01YSU3CBi2fTbmKlhO7QGiBArHmDvuf9uwuEqvW0Y0rkpj1/Y2+9Z4JRSIXRoZ+ldwq9LwZlfSUKB1n09dwkjoUXvmO1gDgUNECG09cAxbvcxwN69Y07mssG+B9hTSoWBywlZK0qfOtr/c+Vpkxtaj592HWWNd1S/Zk8TxhMNECHgdhtmfrONqTrAnlK153iO1bG86TOro7ngUOVpm3X1vJswGtoOBpvOm+KLBoga0gH2lKolxsC+n0rvErK+A+P2ndaWbA1lUfzUUXqHyOY1RmmACJIxhvdW7+Lhj9bpAHtKRUrRcdj6pXWXsOlzOFL5fOw0aFU6GmrHMyGpXuTyGSc0QAQh+5g1wN7nP+sAe0qFXe6vpdNtbvsKXIW+00kCZAwoHefoxB7Bj3OkAA0QAVvw0x7um6cD7CkVNi6HNeBdcdNR9sbK06Y0gs4jrL6ETmdDvaaRy2cdoAHCT4fzHTz80TreX1P2llYH2FMqBI4dsIbG3viZNVR24ZHK055wSuldQsZAa34EFRb6f9YP/9tkDbC353DpAHstGqbwxEU1H2BPqTrJ7Ya935eOhrprNWB8p7WnlJ1Ip3HbiGa1LtMAUYW8Iid/m7+Bfy3fXmb9hD6teWhcaAbYU6rOKDwKWxaXdjAf21d52kZtSzuYO5wBifqCaW3QAFGJVdtzueOd79nmNcBeer0kHrugB6N7tKzFnCkVQ3K2lI6Guu1rcDt8pxMbtBlkNRt1HQXNT9IO5ihQpwOE61g2ua5U5q7Zy47c47RNr8eFfVqQ5j7KjC+3lwkO55xyIo9d0JPmDaJwgL3jOVZnnXdbrMtpvSiU0ji22mi1LNHJ37I4i2D716VPHeVuqfyYqenWOEddRkLnsyG1SThLoIJQZ9/ich3Kwvb2hdi3LKRT/SKWbsymU/0i7Js/o94bI3i+337+ObETrRql8ORFvZl+Wb/oDA6Hd8G/f2v9MeblWuvycmHjApg+rOz6aKdliU7VluVT2PM9vPcneLwD/Ot8WP6S7+DQoieccSdc/TnctRkmTIeeF2pwiFJiTCUdQ+E8qcho4DnABswwxkwtt10828cCecAVxpjVVR2zf//+ZuXKlX6d33UsG9vbF8LuNQA4Oo/iyJAHaPjVX0nc8llJOvfA63CccTfJDaL00bnjOdYfrqccdB0DI6bAooesP95ig64rneA8WmlZolOlZXnQChjFBv7JGu109qVl909Mg47DSzuYG7aKVM6Vn0RklTGmv89tkQ4QImIDNgLnAFnAd8AlxpifvdKMBW7GChCDgOeMMYOqOm4gASL78HFsWxbS5KMrKk+U3hH3xLdIWP8R5GX7ddyIa9bNqlze/1PladI7wqR/w88fRm85QMsSrfwty4Tp8P41kLsVmrQvHQ213emQmBKp3KogVBUgaqMRdCCw2RizFUBEZgPjAe/hFscDM40VvZaLSGMRaWmM2ROKDLy7Zi+d65/EsM6jSNz8WcUEnl/4hHd+b/3CR7PxL1rtuJs+r7it+A93zqXRXw7QskQrf8ry7avQ/yorMDTroh3McaI2+iBaAzu9lrM86wJNg4hcIyIrRWTlgQMH/M7AjtzjTFm0hyNDHvCdYMw0WPK32PjjXTIVhtzme1sslQO0LNGqqrKMf9F6ce230+G0m6F5Vw0OcaQ27iB8/faUb+fyJw3GmOnAdLCamPzNQNv0ejw0IpGGXz3kO8Gnkzl8wVtsbHgOA1pE+ZMmGf1h8d98b/t0MlwyB3quqfrN1GihZYlOVZXlwxvh8v/oQHhxqjZqvyygjddyBrA7iDRBu6hPC2xbFpbpkC4jdyuN5v2ePn/4GBpnhOq0oedylj5j7kvuVpg1yfoDblThBiy6aFmikz9lefO82CiLClhtNDF9B3QRkQ4ikgRcDHxULs1HwB/EMhg4HKr+B4AmCXk0+fLPpSvSO3Lwd59a7anFcrdi/88t0f0oYsEhWDC5dDm9I/xxUYVy8PHN0V0O0LJEq3gqiwpYxAOEMcYJ3AR8BqwH3jHGrBOR60TkOk+y+cBWYDPwKnBDKPNgS2uCe7Tnydr0juy5YC63fWVnzwVzS3/x0zvCb/4OaVE8QmtKYxg9zfqc3tG6issYYP2MpXKAliVaxVNZVOCMMXHxr1+/fiYgx3OMe/7dxnFwp3n72+3mvvd+MG9/u904Du405l8XGHMoK7Dj1ZbjOcbMv7tifg9lxVY5jNGyRKt4KouqAFhpKqlXa+VFuXAI5D2IEvmHfL+klJcbW1dD8VIO0LJEq3gqiyqjqvcg6uxQG0Dlb7DG2i98vJQDtCzRKp7KovxWtwOEUkqpSmmAUEop5VPc9EGIyAFge7UJK9cMiOJBcfwWL+UALUu0ipeyxEs5oGZlaWeM8Tk1ZtwEiJoSkZWVddTEkngpB2hZolW8lCVeygHhK4s2MSmllPJJA4RSSimfNECUml7bGQiReCkHaFmiVbyUJV7KAWEqi/ZBKKWU8knvIJRSSvmkAUIppZRPcR8gRGS0iPwiIptF5B4f20VEnvds/0FE+vq7b6TVsCzbRORHEVkrIgEOWhV6fpTlJBH5RkQKReTOQPaNpBqWI9a+k0s9v1c/iMgyEent776RVsOyRM334kc5xnvKsNYzu+bp/u7rl8pG8YuHf4AN2AJ0BJKA74FTyqUZC3yKNYvdYOBbf/eNlbJ4tm0DmtX2dxJAWU4ABgCPAncGsm8slCNGv5PTgCaez2Ni/G/FZ1mi6Xvxsxz1Ke1L7gVsCOV3Eu93EAOBzcaYrcaYImA2ML5cmvHATGNZDjQWkZZ+7htJNSlLtKm2LMaY/caY7wBHoPtGUE3KEW38KcsyY8xBz+JyrJke/do3wmpSlmjiTzmOGU9EAOpROjVzSL6TeA8QrYGdXstZnnX+pPFn30iqSVnA+sVZKCKrROSasOXSPzX5fxtN30tN8xLL38nVWHerwewbbjUpC0TP9+JXOUTkAhHZAHwCXBXIvtWpjTmpI0l8rCv/XG9lafzZN5JqUhaAIcaY3SJyAvC5iGwwxiwNaQ79V5P/t9H0vdQ0LzH5nYjIcKxKtbi9O5q+E6hZWSB6vhe/ymGMmQfME5GhwCPACH/3rU6830FkAW28ljOA3X6m8WffSKpJWTDGFP/cD8zDugWtLTX5fxtN30uN8hKL34mI9AJmAOONMTmB7BtBNSlLNH0vAf1/9QSxTiLSLNB9K1XbHTHh/Id1h7QV6EBpR033cmnOpWzH7gp/942hstQDGnh9XgaMjuayeKWdQtlO6qj5XmpYjpj7ToC2WPPEnxbs/4cYKEvUfC9+lqMzpZ3UfYFdnr//kHwntfIFRvh/8lhgI1aP/v2eddcB13k+C/CiZ/uPQP+q9o3FsmA9yfC959+6GClLC6yroCPAIc/nhtH2vQRbjhj9TmYAB4G1nn8rq9o3FssSbd+LH+WY7MnnWuAb4PRQfic61IZSSimf4r0PQimlVJA0QCillPJJA4RSSimfNEAopZTySQOEUkopnzRAKKWU8kkDhFJKKZ80QKg6RUSMiDzltXyniEypZp9lITjv/SKyzmvs/kEi0tTzea2I7BWRXV7LSSLi8nxeJyLfi8jtIqJ/sypi4n2wPqXKKwQmiMjfjDHZ/uxgjDmtJicUkVOB84C+xphCz1g5ScYa/yfTk2YKcMwY86TXfvnGmOLtJwBvA42Ah2qSH6X8pVcjqq5xAtOB/yu/wXOF/pPn321e6495ftYTkU88V/M/icgkz/rfi8gKz9X+P0TEVu7QLYFsY0whgDEm23gGhPOXsQaOuwa4SUQqjNQpIr1FZKmI/Cwibs+d0sOBnEOp8jRAqLroReBSEWlUvEJE+gFXAoOwBjr8k4j0KbffaGC3Maa3MaYHsEBETgYmYQ0RnQm4gEvL7bcQaCMiG0XkJRE5M5hMG2O2Yv3NnuC9XkRSgDlYgwGegjV73ZNYAwQqFTQNEKrOMcYcAWYCt3itPh2YZ4w5bow5BrwPnFFu1x+BESIyTUTOMMYcBs4G+gHfichaz3LHcuc75klzDXAAmCMiVwSZfV/j/I8AVhtjVniWfwDSjQ60pmpI+yBUXfUssBp43bPsq+Itwxiz0XOnMRb4m4gsxBoR9E1jzL3V7OsClgBLRORH4HLgjUAyLCIdse5Q9pfb1AMreBXri1U2pWpE7yBUnWSMyQXewZpNDGApcL6IpIlIPeAC4H/e+4hIKyDPGPMWVhNOX+C/wIWeTmREJF1E2pXbr5uIdPFalQlsDyS/ItIceAV4wcedQQ7WhPWISFdgAtYcxErViN5BqLrsKeAmAGPMahF5AyhupplhjFlTLn1P4AkRcQMO4HpjzM8i8gDWHMYJnvU3UjYA1Af+LiKNsTrJN2M1N1Un1dNslejZ71/A0z7SzQLGichPQDZwifGaIU2pYOl8EEoppXzSJiallFI+aYBQSinlkwYIpZRSPmmAUEop5ZMGCKWUUj5pgFBKKeWTBgillFI+/T+sLH6kIFhBywAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEICAYAAACqMQjAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABDgklEQVR4nO3dd3hUVfrA8e/JpJMQUqgJoVcBQ7ciWFAUaRZ0seJa1v6zoWJB14JldS2rK4uK6IpgAQQUywqioiK9hi4QehIIpE85vz/uZNJmkplkJlPyfp4nT+aee+6952Qm9517zr3nKK01QgghhDNh/i6AEEKIwCVBQgghhEsSJIQQQrgkQUIIIYRLEiSEEEK4FO7vAnhTSkqKbt++vb+LIYQQQWXVqlXZWuvmztaFVJBo3749K1eu9HcxhBAiqCil9rhaJ81NQgghXJIgIYQQwiUJEkIIIVwKqT4JIYT/mM1msrKyKC4u9ndRhAvR0dGkpaURERHh9jYSJIQQXpGVlUV8fDzt27dHKeXv4ogqtNbk5OSQlZVFhw4d3N5OmpuEEF5RXFxMcnKyBIgApZQiOTnZ4ys9CRJCCK+RABHY6vL+SJAQogEceu45Dj33nL+L4RWhVJdQYT54EPPBgz7ZtwQJIRpAyZZMSrZk+rsYXhEsdVFKce211zqWLRYLzZs3Z+TIkTVut3btWr766ivH8pQpU3j55ZfrXI76bu8OW1ExtiLf3DAgQUIIEZKaNGnCxo0bKSoqAuC7774jNTW11u2qBonGToKEECJkjRgxgkWLFgEwa9Ysrr76ase6goICJk6cyMCBA+nbty/z58+ntLSUJ554gtmzZ5ORkcHs2bMB2Lx5M0OHDqVjx468/vrrjn288sor9OrVi169evHPf/7Tkf7ss8/SrVs3zj//fLZu3dowlfURCRJCiJB11VVX8cknn1BcXMz69esZPHiwY92zzz7Lueeeyx9//MGSJUt48MEHMZvNPP3004wfP561a9cyfvx4ADIzM/nmm29YsWIFTz31FGazmVWrVvH+++/z+++/89tvv/Gf//yHNWvWsGrVKj755BPWrFnDF198wR9//OGv6nuFPCdhV9YR1+rRR/1cElFG3hNRX3369OHPP/9k1qxZXHzxxZXWffvtt3z55ZeO/oLi4mL27t3rdD+XXHIJUVFRREVF0aJFCw4fPszPP//M2LFjadKkCQDjxo3jp59+wmazMXbsWGJjYwEYNWqUD2voexIk7IKhI66xkfdEeMOoUaN44IEHWLp0KTk5OY50rTWff/453bp1q5T/999/r7aPqKgox2uTyYTFYkFr7fKYoXQrsDQ3CSFC2sSJE3niiSfo3bt3pfQLL7yQN954w3GyX7NmDQDx8fGcPHmy1v0OGTKEefPmUVhYSEFBAXPnzuXss89myJAhzJ07l6KiIk6ePMmCBQu8X6kGJEFCCBHS0tLSuOeee6qlP/7445jNZvr06UOvXr14/PHHARg2bBibN2+u1HHtTL9+/bjhhhsYNGgQgwcP5q9//St9+/alX79+jB8/noyMDC677DLOPvtsn9WtIUhzU4iRdnwhDPn5+dXShg4dytChQwGIiYnhnXfeqZYnKSmpxs7mjRs3Ol7fd9993HfffdXyTJ48mcmTJ9eh1IFHggRgyc2l+f3GG12atZ+w2BjCk5L8XKq6kXb8wBNKn69Qqkuo0BYL4a1aAmArLUWFhaHCvXdqlyAB2AqL2HNV+f3Tnb7/HuRzL7wklD5foVSXUKFtNkp37XIsR3Xtije7zf3SJ6GUaquUWqKU2qKU2qSUqtZgqAyvK6V2KKXWK6X6+aOsQgjRmPnrSsIC3K+1Xq2UigdWKaW+01pvrpBnBNDF/jMYeNv+2zsFyM3FVmg8rl/655+V1pmz9oHVAmEmUEouqYVHtM2GJTsbW2EhWK2YDx2utL50925s+Scr5S/45ZeGLqbbmpx5JirM+D5pOZpdaV3FumizGVtREVYnfQHCu0xxcY7X2mKptE6bzdgqLNe3+ckvQUJrfRA4aH99Uim1BUgFKgaJ0cBMbdyf9ptSqplSqrV923qzFRax8/zzna7be8ONlZbb/fcjsu68i/CUFMJTkjElJxOebLwOT0nBlJJCeHIyYTEx3iiaCEDaZsN67BiW7BysOdlYcnKwHM3GkpONNTvHWM7JwZJ9FGvuMdLff4+9113vdF/7br650nL6zA84+o9XGqIadRJz6qlu1cXyrzfRnTphOXzYaV7hPWGxsZTu3u10XdX0+jY/+b1PQinVHugLVH2CJRXYV2E5y55WKUgopW4BbgFIT0/3SRm11UrR6tW15guLjTUChj1omFKqBJPkZMe6MPvTmN4USp2KDVEXbbViPX4cS3Y2luxsrDk5WLJz7Cf+bPtr+3JOLthste9UiBDj1yChlIoDPgfu1VqfqLraySbVHnHUWk8DpgEMGDDA9SOQDcBWWIht717MLh7tr0jFxjoCRqWrk+Ypla9UkpMJsz/2X/vxQ6dTsa510Var8Y0/J8c48Vc42Vtzsu3f/u3LuT488XvwxK0pLo7kv97km3J4QcWmjdqosDDCU1J8WJqaRaem0qtHD8fylaNH89Bdd7nMn9ipE8d27vRZeV6fNo2/XnNN+RAdEyYw8623aJaQUK/9durZk58++oiUxERHmtaaB6ZO5ZuffiI2Opp3nnmGvj171us44McgoZSKwAgQ/9Vaf+EkSxbQtsJyGnDAW8cPi40xTjyA5eBB9lQYd77dhx8S3ro1oMGmIdxE+gcfGE0Jlb5t5jhORpacHDCb3T6+LizE7ElASU6ucnVSoenLHmwIC+FnI21WijMzy5t7HCf/yu+D9dgxn534wxIS3HofTMnJ2AoK3Ph82fcbG0OLBx7wSZm9wZKb61ZddhzLRUVEENGqlVv7nbdmPy99s5UDx4to0yyGBy/sxpi+tQ/lXZOYmBjWVXiOoVZKuV3eunjjvfe4/vbbibAHzq//9z/v7Fgpojp1IiolBW02U7p7N9/89BM79uxh6+bNrFizhv+7/35+/flnR39SXfklSChjYJN3gS1aa1eNsV8CdyqlPsHosM7zVn8EYDRduPhmGt66NZFplT+skW3a1Lg/rTW2EyecfoN1nMjswcSanY32NKAUFmLet6/autTXX0dbLJgPH8Z6PK/SupLMTCyHDzmWrfn5HP77M24ft6G1fPwxx7dW67HjldYVZ27F1CwBFRlBeOvWqJhYjrxU/4lcwhISqlzR2V9XuaIzJScTFhnp/n6jojz6fAUyd/9X1MkTbneQzluzn0e+2ECR2QrA/uNFPPLFBoB6B4qq8vLyGDRoEF9++SXdunXj6quv5txzz+Vme3/K5MmTWbhwITExMcyfP5+WLVuyYMECnnnmGUpLS0lOTua///0vLVu2ZMqUKcTFxfGAPaj36tWLhQsX0rx5c6688kqysrKwWq08/vjjHD58mAMHDjBs2DBSUlJYsmQJ7du3Z+XKlbz00ku0a9eO22+/HTAmJoqPj+f+++/npZdeYs6cOZSUlDB27Fieeuopp/UKi4wkLDLS0Um9cMkSJowaRVhkJGecfTbH8/I4nJND6wpfSOrCX1cSZwLXAhuUUmvtaY8C6QBa638DXwEXAzuAQuDG6rsJHEopTAkJmBISiOrYsca8WmtsJ08aQaTq1UnZ6+zyTlFdWupyX6ZmCS47FbPuvLPScvrMDzBnZXleuQYSFhPDngnXOF23/+67Ky2nz/zA5X5MCQm19A3ZXycloTw48Qv3PbVgE5sPVG1BLrdm73FKrZWv+IrMVh76bD2zVji/uu7ZpilPXnpKjcctKioiIyPDsfzII48wfvx43nzzTW644Qbuuecejh075ggQBQUFnHbaaTz77LM89NBD/Oc//+Gxxx7jrLPO4rfffkMpxfTp03nxxRf5xz/+4fK4ixcvpk2bNo65K/Ly8khISOCVV15hyZIlpFRpgrvqqqu49957HUFizpw5LF68mG+//Zbt27ezYsUKtNaMGjWKZcuWMWTIkBrrDXDgyBHSKlwVpaWlsX///uAMElrrn3He51AxjwbuaJgSNSylFKamTTE1bUpUxw415tVaY8vPL786qdLcFRYd3UClDixhTZrQ/L77Kn/rT0khPDFRTvxBoGqAqC3dXTExMaxdu7Za+gUXXMCnn37KHXfcwbp16xzpkZGRjulM+/fvz3fffQdAVlYW48eP5+DBg5SWltKhQ83/p7179+aBBx5g0qRJjBw5stbxmvr27cuRI0c4cOAAR48eJTExkfT0dF5//XW+/fZb+vbtCxhDi2zfvt2tIOFsVFpvjEbr97ubAkFYbAztPpkFQHhKc8JiA+dWVqUUpvh4TPHx4OSDWpq13+19RbRsSafvvvVm8bzLg74EU0IzUm65ufaMASCQP1+ecrcutX3jP3PqD+w/XlQtPbVZDLNvPb3+Ba3CZrOxZcsWYmJiyM3NJS0tDYCIiAjHibRsCHCAu+66i/vuu49Ro0axdOlSpkyZAkB4eDi2Cp/T4mJjXumuXbuyatUqvvrqKx555BGGDx/OE088UWOZLr/8cj777DMOHTrEVVddBRgn+kceeYRbb73V7bqpsDAiO3YkvUsXDinl6IPIysqiTS3N5O6QIIHR5rr/nnsBaPfhTP8WxkPudcCX5w3kW2Ld7SAFgupEG8yfr6q8VZcHL+xWqU8CICbCxIMXdqthq7p79dVX6dGjB8899xwTJ07k119/JSIiwmX+vLw8x3zYH3xQ3rTZvn17Fi5cCMDq1avZbX8m4cCBAyQlJXHNNdcQFxfHjBkzgPJhx6s2N4HR5HTzzTeTnZ3Njz/+CBjDlz/++ONMmDCBuLg49u/fT0REBC1atHBZVhUejmXvPi4+4wze+fhj/nLNNfz2228kJCTUu6kJJEgEPU874ANZKNVF1Kysc9rbdzdV7ZO46KKLmDhxItOnT2fFihXEx8czZMgQnnnmGZcdwmB0JF9xxRWkpqZy2mmnOYLBZZddxsyZM8nIyGDgwIF07doVgA0bNvDggw8SFhZGREQEb7/9NgC33HILI0aMoHXr1ixZsqTSMU455RROnjxJamqq42Q+fPhwtmzZwumnG1dTcXFxfPTRR06DRJ8+fQizXzVcduGFvPDoZL5bvZrOnTsTGxvL+++/X8e/YmWqptmVgs2AAQP0ypUr67TtnmuvA4L7m15p1v5KT5F3+v77oD2xhlJdIDQ+X2Vc1WXLli30qPCMgmg45oPGjZ8Rblw5OHuflFKrtNYDnOWXKwkhhAhy7gSHupIgEUIaYwepEMK3JEiEEOkgFUJ4WwiP4yCEEKK+5ErCLqpHd38XQYSwUPp8hVJdRO0kSNi1evRRfxdBhLBQ+nyFUl1E7aS5SQgRMrKyshg9ejRdunShY8eO3HnnnZSUlAAwY8YM7qwyntnQoUMpu22+ffv29O7dm4yMDHr37s38+fMd+TZt2sS5555L165d6dKlC3//+9+dDoOxdOlSxzAfFe3evZvBgwfTpUsXxo8fT2kN47EFGgkSQgj/2L0M3joNTh6u/LqOtNaMGzeOMWPGsH37drZv305RUREPPfSQ2/tYsmQJa9eu5bPPPuNu+6CSRUVFjBo1iocffpht27axbt06li9fzltvveX2fidNmsT//d//sX37dhITE3n33Xc9rp+/SJAIMVE9uodMm3Eo1UVUsXsZfHwlHN0On/+1/PWPL9R5lz/88APR0dHceKMxYLTJZOLVV19l5syZ5Hs47/aJEydItE/o8/HHH3PmmWcyfPhwAGJjY3nzzTeZOnWqW/vSWvPDDz9w+eWXA3D99dczb948j8rjT9InEWJCqb04lOrS6Hz9MBza4Hr9gdVgtg/wt+dn0PZB81Z/AEe3Ot+mVW8Y4frEvGnTJvr3718prWnTprRv354dO3a4Vexhw4ahtWbXrl3MmTPH5X47depEfn4+J06coGnTpjXuMycnh2bNmhFun2ujbAjvYCFBQgjR8Fr0gpytUHKyPECoMEjqXOddaq2dDo1d1nfgatjsiullcz/s3LmT8847j6FDh7rcb037dHZ8T7cLFBIkhBDeV8M3fqC8uUlXGB4+LBzanwkjXU1WWbNTTjmFzz//vFLaiRMnOHz4MN26dSMrK4tjx45VWp+bm+t0hNZOnTrRsmVLNm/ezCmnnMKyZcsqrd+1axdxcXHEx8fXWq6UlBSOHz+OxWIhPDzca0N4NxTpkxBCNLyvHwKLfQrf8BgwRYK1FDbPq/MuzzvvPAoLC5k503hC32q1cv/993PnnXcSExPDwIED+eWXXzh0yJjSd+XKlZSUlNC2bdtq+zpy5Ai7d++mXbt2TJgwgZ9//pnv7cPYFxUVcffdd7vdIa6UYtiwYXz22WeAMfT46NGj61zPhiZBQgjR8K6dD/2vh9hkGPtv6Hut8fqKGXXepVKKuXPn8tlnn9GlSxeSk5MJCwtj8uTJALRs2ZLXXnuNiy++mIyMDO69915mzZrlGG4bjD6JjIwMhg0bxtSpU2nZsqVj7utnnnmGbt260bt3bwYOHFjtdtoy//vf/0hLS3P8/Prrr7zwwgu88sordO7cmZycHG666aY617OhyVDhQgivCLShwpcvX87VV1/NF198Ua3juTGTocKFEAI444wz2LNnj7+LEfSkuUkIIYRLEiSEEEK4JEFCCCGESxIkhBBCuCRBQgghhEsSJIQQISMuLs7fRQg5EiSEEEK4JEFCCOEXQ2cPpfcHvR0/Q2cP9clxFixYwODBg+nbty/nn38+hw8bc1b8+OOPZGRkkJGRQd++fTl58iQHDx5kyJAhZGRk0KtXL3766ScAZs2aRe/evenVqxeTJk3ySTkDlTxMJ4TwuhdWvEBmbmaNeXKKc6ot37j4Rpf5uyd1Z9Igz0/QZ511Fr/99htKKaZPn86LL77IP/7xD15++WX+9a9/ceaZZ5Kfn090dDTTpk3jwgsvZPLkyVitVgoLCzlw4ACTJk1i1apVJCYmMnz4cObNm8eYMWM8LkswkiAhhAhpWVlZjB8/noMHD1JaWkqHDh0AOPPMM7nvvvuYMGEC48aNIy0tjYEDBzJx4kTMZjNjxowhIyODH374gaFDh9K8eXMAJkyYwLJlyyRICCFEXbnzjb/3B72rpb1/0fteL8tdd93Ffffdx6hRo1i6dClTpkwB4OGHH+aSSy7hq6++4rTTTuP7779nyJAhLFu2jEWLFnHttdfy4IMP1jqpUCA4cNyYwKlNsxiv71v6JIQQfpEcnVzjsrfk5eWRmpoKGMN0l9m5cye9e/dm0qRJDBgwgMzMTPbs2UOLFi24+eabuemmm1i9ejWDBw/mxx9/JDs7G6vVyqxZszjnnHN8Uta6KjJbKTJbfbLvOl9JKKXCgfFAR2AXMEdrbfZWwYQQoW3p+KVe32dhYSFpaWmO5fvuu48pU6ZwxRVXkJqaymmnncbu3bsB+Oc//8mSJUswmUz07NmTESNG8Mknn/DSSy8RERFBXFwcM2fOpHXr1jz//POOqU0vvvjioJoPor7qPFS4UupTYACQBbQFVmqtL/di2TwmQ4UL4T+BNlR4Y7LzaD4AnZrX/pxIQw4VfgbQVmttU0qZgH312JcQIe2pBZsAePLSU/xckvoLpbqI2tUnSGwGflFK7QXS7ctCCCc2Hzjh7yJ4TSjVRdSuPh3XlwNzgVz7b782NQkhhPC+Ol9JaK3zgBfrur1S6j1gJHBEa93LyfqhwHxgtz3pC63103U9nhBCCM+5FSSUUte5u0Ot9Uw3s84A3gRqyv+T1nqku8cWQgjhXe5eSbh+Vr4yTc0n/fKMWi9TSrV3c7+iEXphxQuAew9mBbpD4bPtr073azm8IZTqEirM5NpfeX8UXLf6JLTWw9z8OdfL5TtdKbVOKfW1UsrprRRKqVuUUiuVUiuPHj3q5cMLf8rMzax1/J9gURy2j+Kw0LgBMNDrkpWVxejRo+nSpQudOnXinnvuobS0lKVLl6KU4t1333XkXbNmDUopXn75ZQBuuOEGOnToQEZGBv369ePXX3/1VzU8olUpWpX6ZN+1BgmlVLL99+lKqdeUUhE+KUl1q4F2WutTgTeAec4yaa2naa0HaK0HlI2tIoQIDtvOOpst3Xuw7ayzvbI/rTXjxo1jzJgxbN++nW3btpGfn8/kyZMB6N27N7Nnz3bk/+STTzj11FMr7eOll15i7dq1TJ06lVtvvdUr5QpmNQYJpdSNwGL74kRgK2BTSrWyf8P/VCkV7YuCaa1PaK3z7a+/AiKUUim+OJYQwj+s2dmVftfXDz/8QHR0NDfeaLSQm0wmXn31Vd577z0KCwtJT0+nuLiYw4cPo7Vm8eLFjBgxwum+hgwZwo4dO7xSrmBWW5/Ev4CyO4+sWuu3AJRST2F8sy8BngK83mislGoFHNZaa6XUIIyAllPLZkKIAHDoueco2eK6qbBwzRqwWEAp0BqUYkv3HhAeTmzfvk63ierRnVaPPlrjcTdt2kT//v0rpTVt2pT09HTHCf/yyy/n008/pW/fvvTr14+oqCin+1qwYAG9e1cfhLCxqS1ITAI+AQYBB5VSp2MMwzFAa90fQCk1qy4Htm83FEhRSmUBTwIRAFrrf2M8d/E3pZQFKAKu0nUdQ6QRCaXOXhHCLBbjd9m/dNnvsvQ60lqjlKox/corr2T8+PFkZmZy9dVXs3z58kp5H3zwQZ555hmaN29eqf+isaoxSGit31BKfWZffBHjobnOGE1PZQrrcmCt9dW1rH8T4xZZ4YFQ6egVwa22b/zbzjoba3Y2KjISXVrq+G1KSaHdh+7eRV/dKaecwueff14p7cSJE+zbt49OnToZZWvVioiICL777jtee+21akHipZde4vLL5dngMrV2XGutD9p/F2mtL9Jad9ZaLwNQSsUC1cO2EELUoOvPP9Ejcwu61LgjR5eW0iNzC11//qle+z3vvPMoLCxk5kwj0FitVu6//35uuOEGYmNjHfmefvppXnjhBUwmU72O1xjUeVgO+6B+fwEe9l5xhBCNiSklpdLv+lJKMXfuXD799FO6dOlC165diY6O5rnnnquU74wzzmg0M8vVV32G5bAC071YFr+StnwhGl59rxycadu2LQsWLKiWPnToUIYOHVotvWymOoAZM2Z4vTwNIS6yCSblmznkZPpSO2nLF0IEG4vNhs2miaAZFovmyMlimsVEEBamCA/zTtCoU5Cwj+U0V2t90iulEEII4RGL1cahvGJyCys/aX0or5ik2EhaJUQTbqp/oKjrHt4HWtf76EKIkCJ3qTcMi815gCiTW1jKobxiLDZbpfS6vD91bW6SO5qE8MCZaf1oGhXv72J4hau6REdHk5OTQ3JystNnFUTNtNbYNNi0xmar8Lrsx1a2DPHR4S4DRJncwlJaNI1yXAporcnJySE62rNBMqRPQggfyisyc6LITCqjOXLUyttLd3Lpqa2Jj44gIaahhkHzjtrqkpaWRlZWFqE+0KbWoNHGb63RgK3stX1d5WXj5K4rpNmosL7CPtwRHRFGpCmME8W1P3hYfDSc+Ojyz1l0dDRpaWke1VeChBA+klNQwotfb2XOqn1UvMp/8ZtMruzflodGdCO5ifMhIQKNu3Xp0KGD/wpZhdWmKSy1UFhqpbDUSkGJhSKz8bsszbG+xEJB1bRSCwUlVopKrRSUWhy/i8222g/uQ2P7phIZDrP/OFhr3gmD0nl2XI96HU+ChBA+kFdk5sWvtzJ7ZfUhtbXGkf7oJT0C/orC13WxWG0Umu0n4woncMeJ2dlJ3X4Cd3pStweCEot/T+Z1ERtpsv+EO143iQonJsL+O9LEqWkJHDpR7Nb+2ibH1p6pFnUNEvsBc72PLkSIOllsZs6qmudcmLNqH3cM60xRqbWBSlU3pRabW3W5eUgHFqw7wIliM4UuTuCFZuNbuyOt1Epp0J7MK5/Iq57gy0/uFfOWraueFhNhIizMvb6crGOF/PP77dTUD60UjOxT//uL6hQktNZt631kIULYgnUHa/wHBuNb+Ker9rFhfx5LtwZmO/6wbi04JbWpW3WZv/ZAQNVFKYiNMBEbFV79pF72usLJukmk8U29SZSJmIjwSifyiuuiw90/mftKfHQEV/Zv6/TqrsyV/dtW6o+oK2luEsIH9uYWuJXv6MkSmsVE+rg0ddcsNoKjJ0vcylvXuihF+Um44sk8Ktx+kq9wUq9wci/fpvwEXjEtOiIsZO+ySoiJ4KER3QCq9RMpZQSISSO6eaUp0+0goZQ6AxgPPKC1lqYm4RPHio9RaDEGFr63/70A7M/f71gfGx5LYnSiP4rmkfSkJrxzfTea1TDl8PF82JRlZV9uIS2bBmYHtkLTIj7KrbpkHrCRGBtJy4QoYqt8E3ec4Cud8I0mmajw0D2Z+1KYqZA7hydxxwVJFJRasNo04aYwYiNMKAXKVAjU/3PlyZXEjcAa7DPTAd8A24Brtdbu9aIIUYtCSyEXfX6Ry/WLL1tMIoEdJIrNVoZ0SSHfdpS//jDWZb7p587lygFtSUusf+eiL2UdK2R//oFa63JZv7SAr0soKbQUcvFc57Pqgff+VzwJEg06M51oXMw2M8eKj2G21nyRerLkJG9uD9xpRkosNhauO0D7lFjuPu2yGvNGRpQwe8d0osJ9MzCbt5RYbJzb9rwa8wRLXULJ+ennN8hxPAkSXp2ZToQ+s81MblEuOcU55BTlOH5nF2WTU5xTad2xkmMAvHfhezXu86T5JO+sf6chil93EZCdBxZcf8sDsFDIjM3/aaBC1c+QtqfXuD6Y6hIqBrce3CDH8SRIeHVmOhGcyk782cXZxom/wsm/0uviHI6XHPd3cYUQ9eR2kNBaFwGVGotlZrrAUN/OXrPVbJzcnZ3si3LKA0JxDnkleT6pg0LRLKoZ4WE1fyTjI+O5PeN2n5Shrg7lFTNvzX6KzMbzDkophvdsSXxkzWM1BWJdXAmluoSK2t4Tb6nzLbAyM13gqK2zd97oeSzdt9TRzFO16edE6QmflEuhSIxOJCk6iZSYFJJjkkmOTnb8rpiWGJ1IeFh4peDmTHxkPH879W8+KW9d/LjtKG98vorC0q4AxESYeOuafgzr1iLo6lKTUKpLqKjtPfEWmZkuiJRaS6s16WQXZXNW6lk1bpdbnMsTy5/wShnKTvyVTvRlJ/4qJ393rgyC2bw1+3ng03VY7COzNYuN4L0bBtIvPbDvvhLCE279B9snGXKL1npm3YvT+JRYS5y252cXZVdLO1nqfI6nvi361qsMYSqMxKhEp9/wK/2OSSYxKhFTmO8mj48Nj2XxZYsByC7KBiAlJqXS+kAw/addPLNoi2O5TUI0M28aTOcW5Q8TVKzLwfwjALSOa1FpfbAIpbqEioZ6T9z9mnejm/k0EDRBwlcPbhVbip2271ds7sktziWnKIeTZt9P7hdliuLSjpc6Tv5J0UmVAkKzqGY+PfF7IjE60XFv92M/PwbA+xe9788iVaK1ZurXmbyzbJcjrWvLOD6YOIjWCTGV8lasy7hP7wbg9xs/b7jCelEo1SVUNNR74laQ0FoP8/qRA4AnD24VWYqqfbN3dgWQU5RDvjnfJ+U1KVOlE3zZ76aRTWvcLjkmmefOfs4nZWpMzFYbkz5fzxery79IDGiXyLvXDyQhNrBHchWirkK3wdgLsguzueXbW8gpzqHA7N5YPJ4yKZPjhJ8Uk+R4nRKdUqmdv6yNP0xVf1ipoTqwGrPCUgt3/Hc1SyoMXnd+j5a8+Ze+REcExlWYEL4gQaIGpbZS9p7c6/F24Sq80gnf1R09yTHJJEQlOD3xi8BxrKCUG2f8wdp9xx1p4we05dmxvbwy0bwQgcyTAf4UkKa1rnlg+RAVHhZefitnlZN/1bSmUU0b9MQfLJ29wWj/8SKue/d3dh4tv5K8c1hn7h/eVQalE42CJw/TaaXUPKC/74oTWJJjkpk/ej7JMUa7f6CeFAK9szdYbT10kuvfW+GYBUwpeHJkT244M3Cm6BTC1zxtbvpNKTVQa/2HT0oTYKJMUaTGpfq7GMIP/vgzl5tm/OGYbD7CpHh1fAYj+7Txc8mEaFieBolhwK1KqT1AAcaQHFpr3cfrJRPCT77bfJg7P17tmCO5SaSJadcN4MzOKbVsKUTo8bRP4jZgj++K07CkLV9UNfuPvTzyxQbsD1GTEhfJjBsH0Ss1wb8FE8JPPO2TeLVsiPBQIG35oozWmreW7uSlb7Y60tKTYpk5cRDtU5r4sWRC1C420ne3YUufhGj0bDbN0ws3M2P5n460nq2bMmPiQFrER/uvYEK4qX2y777ISJ+EaNRKLFbun7OOhesPOtLO6JTMO9f2Jz5anqIWwtMgUfNUW0IEkfwSC7d+uJJfduQ40i7p3ZpXxp9KVLg8RS0EeBgktNYh02ktGrejJ0u4ccYKNu4vn0vjutPb8eSlp2AK8/7zML5sM25ooVQXUTu3HgtWSj1U4fUVVdbJyHEiqOzNKeTyfy+vFCDuv6ArT43yTYAAo83Yl+3GDSmU6iJq5+7YEVdVeP1IlXWuh1GtgVLqPaXUEaXURhfrlVLqdaXUDqXUeqVUv7ocR4iKNu7PY9zby9mTYwwRH6bg+XG9ueu8LgH7RL0Q/uRukFAuXjtbdtcMag4wI4Au9p9bgLfreBwhAFi+M5urpv1Gdn4JAFHhYbx9TX+uHpTu55IJEbjc7ZPQLl47W3Zvh1ovU0q1ryHLaGCm1lpj3HrbTCnVWmt9sIZtRAj5M8d7w7MvWn+Q/5u9llKr8RR1fHQ4714/kEEdkrx2DCFCkbtB4lSl1AmMq4YY+2vsy766kTwVqDjibJY9rVKQUErdgnGlQXq6fCPsntTd30XwGluxd8ZJ+vDXP3niy01o+9eZlk2j+GDiILq3qnmyJm8KpfcllOoSKnz5nrg7M50/bmdw1oxV7apFaz0NmAYwYMCAOl3VhJJJgyb5uwhe08oyvl7ba6159bttvP7DDkdax+ZNmDlxEGmJDTvkSii9L6FUl1Dhy/ckkCcdygLaVlhOAw74qSwiyFisNh6fv5FZK8ovRk9t24z3bxhIUpNIP5ZMiOASyNNqfQlcZ7/L6TQgT/ojhDuKzVZu/+/qSgHinK7NmXXzYAkQQnjIb1cSSqlZwFAgRSmVBTwJRABorf8NfAVcDOwACoEb/VNSEUzyiszc/MFKVvyZ60gb2zeVFy/vQ4RMNSqExzwKEkqpF7TWk2pLc4fW+upa1mvgDk/3KxqvwyeKuf69FWQeOulIu/nsDjwyogdhPnpITohQ5+lXqwucpMl4TsLvdh7NZ9xbyysFiEcv7s7kS3pKgBCiHty6klBK/Q24HeiklFpfYVU88IsvCiaEu9buO86N76/gWKEZgPAwxYuX92FcvzQ/l0yI4Oduc9PHwDfAdCr3DZzUWuc630QI31u69Qh/+2g1RWYrADERJt66ph/DurXwc8mECA3uPieRB+QppZqF6kiw8oBQ8Jm7JosHP12PxT7XaGJsBO/dMJC+6Yl+LpkQocPTu5t+DdWZ6eQBoeAy/addPLNoi2M5tVkMH0wcROcWcX4slRChpy4z092mlPoTmZlO+IHNpnlhcSbvLNvlSOvWMp4PJg6iVYJMNSqEt8nMdCJg9U1vRlxU+UfUbLUx6fP1fLF6vyNtYPtEpl83kIRYmWpUCF/w6BZYe39EM+BS+0/I9lEI/8krMrMvt5Co8DD2HSvi7aU72ZdbyKL1B1lUYS7q83u05MObBkuAEMKHPH2Y7h7gZuALe9JHSqlpWus3vF4y0SjlFJTw4tdbmbNqn2PUVoAXv8lkTEYq067tzy0frmJs31SeGdOLcHmKWgif8rS56SZgsNa6AIynrYFfAQkSot7yisy8+PVWZq/cV22d1jB3jdHMNP36AZzVOUVmkhOiAXj6NUwB1grLVuo+M50QlZwsNjNnVfUAUdG8tfvpkNJEAoQQDcTTIPE+8LtSaopS6ingN+Bd7xdLNEYL1h2s1MTkjNawcL0MBixEQ/GouUlr/YpSailwlj3pRq31Gq+XSjRKe3Pdm650X06hj0sihCjjacd1NMbw3mcDNsCklNqitS72QdlEI2K22mjt5nMObZMbdlY5IRozT5ubZgKnAK8DbwI9gA+9XSjRuKzac4yx//qFszo3p7auBqVgZJ/WDVMwIYTHdzd101qfWmF5iVJqnTcLJBqPE8VmXlq8lY9+34PW8PvuHMZkpDruYnLmyv5tiY+W5yKEaCieBok1SqnTtNa/ASilBiNDhYs6+GbTIZ6Yv5HDJ0ocae8s28WcW08n0hRW7TkJpYwAMWlENxJiJEgI0VCUru12koqZldoCdAP22pPSgS0Y/RN+H8NpwIABeuXKlf4sgqjFobxinvxyI99sOlwpfVi35vx9TC/SEmPJKzJzstjMZ6uyOJxXTLuUJozs05r46AgJEEL4gFJqldZ6gLN1nl5JXOSF8ohGyGbT/HfFXl78OpOTJRZHekpcFFNG9eSS3q0dzz4kxBjBoKjUSpvEGG47p5O/ii1Eo+fpLbAyTpPw2LbDJ3nkiw2s2nOsUvpVA9vyyIgeLsdeWrvvOAB3ndvF10UUQrhQl1tgb8d4TkIDPwNvyy2wwplis5W3luzg7R93YraWN2t2TGnCc+N6c1rHZD+WTgjhDk+bm2YCJykfq+lqjFtgr/BmoUTw+31XDo/M3cCuo+UPyEWYFH87pxO3D+tMdITJj6UTQrhLboEVXpVXaOb5r7fwyR+Vx2Dql96MqZf1oWvLeD+VTAhRF3ILrPAKrTWLNhxkypebyc4vv601LiqcSSO6M2FQOmFhMiifEMHG0yAxGLhOKVXpFlil1AYC4BZY4R/7jxfxxLyN/C/zSKX04T1b8vToXjKtqBBBTG6BDTFPLdgEwJOXnuLzY1ltmpm//snL32yloLR8BPmWTaN4alQvLurVyudlEEL4lqdB4npniVrrp71QFuEFmw+caJDjbDl4goe/2MA6+22qZa45LZ2HLupOUxk6Q4iQ4GmQqDiWczQwEuOJa9FIFJutvPa/7fxn2S4stvLbWru0iOP5cb0Z0D7Jj6UTQnibpw/T/aPislLqZeBLr5ZIBKxfdmTz6NwN7Kkwn0OkKYw7z+3Mbed0IjJc5psWItR4eiVRVSzQ0RsFEYHrWEEpzyzawuersyqlD+qQxPPjetOpeZyfSiaE8DVPn7jegPGkNYAJaA5If0SI0lozf+0Bnl64mdyCUkd60+hwHr24B1cOaCu3tQoR4jy9khhZ4bUFOKy1trjKLILXvtxCJs/byLJtRyulX9KnNU9e2pMW8b6/rbVnm6Y+P4YQomZuBQml1EBgX9kAf0qp64DLgD1KqSla61wfllE0IIvVxnu/7OaV77ZRbLY50tskRPP3Mb04r0fLBitLQ9zGK4SombtXEu8A5wMopYYAU4G7gAxgGnC5LwonGtaGrDwe/mI9myrcRqsU3HBGex4Y3o0mUfXtwhJCBBt3/+tNFa4WxgPTtNafA58rpdb6pGSiwRSWWnjl222898tuKtzVSvdW8Uy9rA8ZbZv5rWxCCDd8/bDxe8RUr+/a7SChlAq39z+cB9xSh32IALR06xEmz93I/uNFjrSo8DDuPb8rfz27AxEmua1ViIB3aIPPdu3uCX4W8KNSKhsoAn4CUEp1BvJ8VDbhQ9n5Jfx94Wbmrz1QKf3Mzsk8O6Y37VOa+KlkQohA4laQ0Fo/q5T6H9Aa+FaXT4wdhtE34TGl1EXAaxi30k7XWk+tsn4oMB/YbU/6Qob/qD+tNZ+uyuLZRVvIKzI70pvFRvD4JT0Z1y/VMY2oEEK43VRUNjx4lbRtdTmoUsoE/Au4AMgC/lBKfam13lwl609a65HVdiDqZHd2AY9+sYFfd+VUSh/bN5XHLulBclyUn0omhAhU/upPGATs0FrvAlBKfQKMBqoGCeGhvunNiKtyF5LZamPasl289r/tlFrKb2tNS4zh2bG9Oadr84YuphAiSPgrSKQCFacuy8KYq6Kq0+0z3x0AHtBab6qaQSl1C/aO9PT0dB8UNTjkFZk5UWQmKjyMfceKeHvpTi49tTUlZhv3zl7Dhv3lt7WawhR/PasD95zfhdhIue9ACOGav84Qzhq9dZXl1UA7rXW+UupiYB7QpdpGWk/DeFaDAQMGVN1Ho5BTUMKLX29lzqp96Ap/gRe/yWRMRioPDO/GLR+uosRio1dqU6aO60Ov1AT/FVgIETT8FSSygLYVltMwrhYctNYnKrz+Sin1llIqRWud3UBlDAp5RWZe/Hors1fuq7ZOa5i7Zj8A9w/vSphS3HBGe8LltlYhhJv8dbb4A+iilOqglIoErqLKkONKqVbKfpuNUmoQRllzqu2pkTtZbGbOquoBoqJ5a/dz4Smt+OvZHSVACCE84pczhv2hvDuBbzAmLZqjtd6klLpNKXWbPdvlwEZ7n8TrwFUVbr0VdgvWHaS2v4rW8PXGQw1TICFESPFbr6XW+ivgqypp/67w+k3gzYYuV7DZm1tQeyZgX4WJgoQQwl3S9hDkWifEuJWvbXKsj0sihAhFEiSCVKnFxpPzN3JW5xRqe0BaKRjZp3XDFEwIEVIkSAShnPwSrnn3dz74dQ+/785hTEZqjfmv7N+W+OiIBiqdECKUSJAIMpsO5DHqzV9YsdsYuf2f32/n+jPaM35A22pXFErB+AFtmTSiGwkxEiSECDkFOWCtMjmo1QIF2dXT60getw0ii9Yf5IFP11FktgJGELjn/C6cmpZAh5Qe3HVeZz5blcXhvGLapTRhZJ/WxEdHSIAQIhTl7YfZE+DsB+H0O2D9HCjMhT3LYfHDcNEL0O50iE2q12EkSAQBm03zz++38foPOxxp8VHhvHZ1Bud2N6YTTYgxgkFRqZU2iTHcdk4nfxVXCOFrBTlGgDiwBmb/BboMh7MfgHl/g22LjTyz/wKDb4Ohj0BMszofSoJEgMsvsfB/s9fy3ebDjrT2ybFMv34AnVvEV8u/dt9xAO46t9oIJkKIUBGdYFxBzP6Lsbz9W+OnoqSOcMbd9QoQIH0SAW1PTgHj3vqlUoA4u0sK8+84y2mAEEI0AicPw5qZYAqHLhc4z5PUEa5fCAk139TiDrmSCFA/b8/mjo9XV5oY6OazOzDpou4ytIYQjU3ubshcCFsWwr7fAQ0JbWHsv2H7d9Xzj5vmlQABEiQCjtaa93/5k2e/2oLVZoy3ERkexvNje3NZ/zQ/l04I0SC0hsMbjaCQudB4XdXQSfDzq863/+IWuZIIRSUWK4/N3cinq7IcaS3io5h23QAy2jbzX8GEEL5ns0HWCtiywAgMx/50nk+FwRl3QUQs7PjeeZ7cXfDBSK8ECgkSAeLIiWJu/WgVa/Yed6RltG3GO9f2p2XTaP8VTAjhO5ZS+HOZPTB8BQVHnOczRUGnYdB9JHQbYaRNG1q+PqkjXPIPWHS/ESDA+L3gLhg3vV63wUqQCADr9h3n1g9XcehEsSPtsn5pPDu2F9ERJj+WTAjhdSX5xhVA5kLY9i2U5DnPFxkPXYcbgaHLBRBV4WYVq8V4DmL2X4wAccUM2PC5ceXwwUgjQCR1hEvfkOckgt28NfuZ9Pl6SuxzT4cpePTiHtx0VgdUbYMyCSGCQ2EubP3aCAw7fwBLsfN8sSnQ/RLocSl0GALhUc7zmcKNB+UG32bc5vr9FDhxwGhaun6hcQVx6RvSJxHMrDbNi4szeWfZLkda0+hw3vxLP4Z0be7HkgkhvCJvP2Qugi1fGk9Ba6vzfM3Soful0GMktB0MYW62HsQmlT8od6LCxJ4JqfVuYqpIgoQf5BWZueeTNSzdetSR1rlFHP+5bgAdUpr4sWRCiHo5ug0yFxh3JR1Y7Tpfi55GM1KPS6FVb2odytkVVw/KeSlAgASJBrfzaD43f7CSXdnlkwWd36MFr47PkJFahQg2WhtDY5Q9w5C91XXetEHG1UL3kZAcPMPmSJBoQEsyj3D3rDWcLCkfnfGOYZ24/4JuhIV5p/+hZ5umXtmPEMIFqwX2Lrc/w7AITmQ5zxcWDu3PNgJDt0ugaXDO6SJBogForXln2S5eWJzpmI86OiKMly4/lUtPbePVYz156Sle3Z8QAjAXw64lRmDY+hUU5TrPFx4DXc43+hi6DoeYxIYtpw9IkPCxYrOVSZ+vZ/7a8o6lNgnRTLtuAL1SE/xYMiFEjYrzjCEvtnwJ278Hs4v55KObGc8udB8Jnc6FyNCaKliChA8dzCvi1g9XsT6r/D7oge0Tefua/qTEubi1TQjhP/lHjCakzIWw60ewmZ3ni29dfqtquzPBFLr9iRIkfGTVnlxu/XA12fkljrSrB6Xz1KhTiAyXAfqECBjH/iwfI2nvb4B2ni+pkxEUelwKbfpBWOP4P5Yg4QNz/tjHY/M2Umo1HpALD1M8eWlPrjmtnTwgJ4S/aQ1HNhtDYWxZCIc3uM7b+tTyZxiad6/7rapBTIKEF1msNp5ZtIUZy/90pCXGRvDWhP6c3inZfwUTorGz2SDrj/JnGI7tdp5PhUH66Ub/QvdLILFdw5YzAEmQ8JJjBaXc8fFqlu/McaR1bxXPf64bQNuk0OrIEiIoWErhz5+MZqTMRZB/2Hk+UyR0HGa/VfViaJLSsOUMcBIkvGDroZPcPHMle3MLHWkjerXi5StOpUmU/ImFaDClBcbgeVsWwrZvahg8L86YF7rHSOh8AUTL80WuyBmsnr7ddIj/m72WgtLycVnuu6Ardw7r7LUH5IQQNSjMhW2LjcCw8381D57XbQT0GAUdz3E9eF4watXbZ7uWIFFHWmve+GEHr3y3zZHWJNLEK+MzuPCUVn4smRCNwIkD5YPn/fmL68HzEtLLh8JIP839wfOCzYipPtu1BIk6KCy18MCn6/hqwyFHWtukGKZfN5BureJr2FIIUWfZ28tnbdu/ynW+5j2MwNDjUmjVp1HekeRNEiQ8tC+3kJtnriTz0ElH2hmdkvnXX/qR2CTSjyUTIsRoDQfXlj/DcDTTdd60geWjqgbR4HnBQIKEB37blcPt/11NbkGpI+2GM9oz+ZIeRJgax4M1QviUzQp7f7VfMSyCvH3O84WFQ/uzym9VberdMdBEOQkSbvrwtz089eUmLDbjacwIk+KZMb0YPzDdzyUTIsiZi2HXUuMZhq1fQ2GO83zhMdD5PONqoeuFITF4XjCQIFGLUouNKQs28fHvex1pKXFRvHNtP/q3897EHkI0KsUnYPu3RjPS9u+gNN95vugE6DrC6GPodF7IDZ4XDCRI1CA7v4TbP1rNij/LhwXunZrAO9f2p02zGD+WTIgglH8Uti4y+hh2/wjWUuf54lrZB88baczHEMKD5wUDCRIubDqQxy0zV7H/eJEjbXRGG164rA/RESF6G50Q3nZsT/msbXt/xfXgeR3tHc+jILV/oxk8Lxg0+iBhzc8m1xrDZ2sOsTe3gPSkJozu04I2EQX0btOE/ceLUAomXdSdW4d0DNwB+gpyjEtzU4W31GqB4uPGePemIHqrpS6ByZ26aA1Httg7nhfAoRoGz2vVx+hf6D4SWvSQW1UDVKMO19bjWZg+vpzwnd/SKa6UZduy6RRXSszu70j8aDhv9D/CjKs6M/PGQdx2TqfADRB5++G/lxnDEBTam8bKnkKdNrRyeqCTugSm2uqSuRD2LId3hsDbp8PS55wECAXpZ8CFz8E96+C2n+Cch6BlTwkQAUxp7eLyz9cHVuoi4DXABEzXWk+tsl7Z118MFAI3aK1X17TPAQMG6JUrV7p1fGt+NqaPLzcmMQfMnS/kxJmP0fTnZ4jY+Y0jnx58G2roIxDTzN2qNayCHOOf114Puo6A86fA908a/8BlBt8GgVwPkLoEKpd1ecIIGmUG3Qwdh8InE8rTTJFGWnf74HlxzRuw4MJdSqlVWusBTtf5I0gopUzANuACIAv4A7haa725Qp6LgbswgsRg4DWt9eCa9utJkMjOK8C081sSv7zBdaakjtiu/IiwLV9CYbZb+21wKd2ME8wXN7vOk9QRxv8XNs8P3HqA1CVQuVuXcdPgi1uM2d26XGAEhi7DZfC8IFBTkPBXg+ggYIfWeheAUuoTYDSwuUKe0cBMbUSx35RSzZRSrbXWB71RgE/XHKJzXHeGdr6QiB3fVM9g/9CHzbkGcnd545C+M/pfxj/l9u+qryv75509IfDrAVKXQFVrXaYbfRAXTYUO50BEdMOXUfiEv/okUoGKj1Jm2dM8zYNS6hal1Eql1MqjR4+6XYC9uQVM+f4gJ858zHmGES/A0ueD4x946VQ4817n64KpHiB1CVQ11WX0W9A6w2iC6nqhBIgQ468rCWe9VFXbvdzJg9Z6GjANjOYmdwuQntSEJ8+PoOnPTzrP8PUk8sZ+xLamFzCwVYDfgZI2AJY873zd15Pg6tnQew2UnGjYctWF1CUw1VSX+bfD9Qshodp3OBEC/HX2ywLaVlhOAw7UIU+dXdG3Faad31bqpK4kdxcJc6+h73ULoFmatw7rfVaL0Xm43XU9mDU+OP6JpS6ByZ26fDAyOOoiPOav5qY/gC5KqQ5KqUjgKuDLKnm+BK5ThtOAPG/1RwAkhhWS+OPj5QlJHTn2l6+N9tUyubsIX3h3YN+mWHwcFk8qX07qCH/9vlo9WHBXYNcDpC6BKpTqIjzmlyChtbYAdwLfAFuAOVrrTUqp25RSt9mzfQXsAnYA/wFu92YZTLGJ2C6y33Wb1JGDYz/j3p/DOTj2s/IPf1JHuPQNiA3gMZqim8FFLxivkzoa3+bSBhq/g6keIHUJVKFUF+E5rXXI/PTv3197pCBH2756SJuP7dMf/75HP/r5ev3x73u0+dg+rT8cq/XxLM/25y8FOVp/9VD18h7PCq56aC11CVShVBdRDbBSuziv+u1hOl/w5DkJh6Ljzh9kKswNrm9FoVIPkLoEqlCqi6ikpuckGvWwHIDrJ12D7UMfKvUAqUugCqW6CLdJkBBCCOGSBAkhhBAuhVSfhFLqKLCnHrtIAQJ4EB23hUo9QOoSqEKlLqFSD6hfXdpprZ2OvhhSQaK+lFIrXXXeBJNQqQdIXQJVqNQlVOoBvquLNDcJIYRwSYKEEEIIlyRIVDbN3wXwklCpB0hdAlWo1CVU6gE+qov0SQghhHBJriSEEEK4JEFCCCGES40iSCilLlJKbVVK7VBKPexkvVJKvW5fv14p1c/dbRtaPevyp1Jqg1JqrVLKw0GuvM+NunRXSv2qlCpRSj3gybYNqZ71CLb3ZIL9c7VeKbVcKXWqu9s2tHrWJWDeFzfqMdpeh7X2WTrPcndbt7ga+S9UfgATsBPoCEQC64CeVfJcDHyNMRveacDv7m4bLHWxr/sTSPH3e+JBXVoAA4FngQc82TYY6hGk78kZQKL99Ygg/19xWpdAel/crEcc5f3LfYBMb74njeFKYhCwQ2u9S2tdCnwCjK6SZzQwUxt+A5oppVq7uW1Dqk9dAk2tddFaH9Fa/wGYPd22AdWnHoHGnbos11ofsy/+hjFjpFvbNrD61CWQuFOPfG2PCkATyqd59sp70hiCRCqwr8Jylj3NnTzubNuQ6lMXMD483yqlVimlbvFZKd1Tn79tIL0v9S1LML8nN2FctdZlW1+rT10gcN4Xt+qhlBqrlMoEFgETPdm2Nv6a47ohKSdpVe/7dZXHnW0bUn3qAnCm1vqAUqoF8J1SKlNrvcyrJXRfff62gfS+1LcsQfmeKKWGYZxYy9q/A+k9gfrVBQLnfXGrHlrrucBcpdQQ4O/A+e5uW5vGcCWRBbStsJwGHHAzjzvbNqT61AWtddnvI8BcjMtRf6nP3zaQ3pd6lSUY3xOlVB9gOjBaa53jybYNqD51CaT3xaO/qz2QdVJKpXi6rUv+7pjx9Q/G1dIuoAPlnTenVMlzCZU7e1e4u20Q1aUJEF/h9XLgokCuS4W8U6jccR0w70s96xF07wmQjjHv/Bl1/TsEQV0C5n1xsx6dKe+47gfst///e+U98csb6Ic/9MXANoye/sn2tNuA2+yvFfAv+/oNwICatg3GumDc4bDO/rMpSOrSCuPb0AnguP1100B7X+pajyB9T6YDx4C19p+VNW0bjHUJtPfFjXpMspdzLfArcJY33xMZlkMIIYRLjaFPQgghRB1JkBBCCOGSBAkhhBAuSZAQQgjhkgQJIYQQLkmQEEII4ZIECSGEEC5JkBCNjlJKK6X+UWH5AaXUlFq2We6F405WSm2qMPb/YKVUsv31WqXUIaXU/grLkUopq/31JqXUOqXUfUop+b8VDaYxDPAnRFUlwDil1PNa62x3NtBan1GfAyqlTgdGAv201iX2sXUitTFeUIY9zxQgX2v9coXtirTWZetbAB8DCcCT9SmPEO6SbySiMbIA04D/q7rC/k19o/3n3grp+fbfTZRSi+zf6jcqpcbb069RSq2wf+t/RyllqrLr1kC21roEQGudre2DyLlLG4PN3QLcqZSqNsKnUupUpdQypdRmpZTNfsX0lCfHEKIqCRKisfoXMEEplVCWoJTqD9wIDMYYHPFmpVTfKttdBBzQWp+qte4FLFZK9QDGYwwvnQFYgQlVtvsWaKuU2qaUekspdU5dCq213oXxf9uiYrpSKhqYjTGAYE+MWfBexhhUUIg6kyAhGiWt9QlgJnB3heSzgLla6wKtdT7wBXB2lU03AOcrpV5QSp2ttc4DzgP6A38opdbalztWOV6+Pc8twFFgtlLqhjoW39k8AecDq7XWK+zL64EkLYOziXqSPgnRmP0TWA28b192dvKtRGu9zX7FcTHwvFLqW4yRRD/QWj9Sy7ZWYCmwVCm1AbgemOFJgZVSHTGuVI5UWdULI4CV6YdRNyHqRa4kRKOltc4F5mDMSgawDBijlIpVSjUBxgI/VdxGKdUGKNRaf4TRnNMP+B9wub1jGaVUklKqXZXtuimlulRIygD2eFJepVRz4N/Am06uEHKAPvZ8XYFxGHMaC1EvciUhGrt/AHcCaK1XK6VmAGVNNtO11muq5O8NvKSUsgFm4G9a681Kqccw5kQOs6ffQeUgEAe8oZRqhtFxvgOj6ak2MfYmrAj7dh8CrzjJNwsYpZTaCGQDV+sKM60JUVcyn4QQQgiXpLlJCCGESxIkhBBCuCRBQgghhEsSJIQQQrgkQUIIIYRLEiSEEEK4JEFCCCGES/8PoJsMYygdBOEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "\n",
    "# plot recovery error as a function of noise std\n",
    "df = results\n",
    "df = df[df[\"Cardinality\"]==K]\n",
    "df = df[df[\"M\"]==M]\n",
    "df = df[df[\"Noise STD\"]<=0.3]\n",
    "df = df[df[\"Method\"].isin(['QUBO L0', 'Exhaustive L0', 'Lasso', 'OMP'])]\n",
    "\n",
    "\n",
    "sns.lineplot(x = \"Noise STD\", y = \"Infeasible X-Recovery L2 Err\", ci=68, err_style='bars', style=\"Method\",markers=True, dashes=False,markersize=11,linewidth=3,\n",
    "             data = df[[\"Experiment\", \"Method\",\"Noise STD\",\"Infeasible X-Recovery L2 Err\"]],hue='Method')\n",
    "plt.ylabel(r'Reconst. Err $\\mathbf{\\|x-\\hat{x}\\|_2/\\|x\\|_2}$')\n",
    "plt.xlabel(r'Noise STD $\\sigma$')\n",
    "plt.show()\n",
    "\n",
    "# plot recovery error as a function of noise std\n",
    "df = results\n",
    "df = df[df[\"Cardinality\"]==K]\n",
    "df = df[df[\"M\"]==M]\n",
    "df = df[df[\"Noise STD\"]<=0.3]\n",
    "df = df[df[\"Method\"].isin(['QUBO L0', 'Exhaustive L0', 'Lasso', 'OMP'])]\n",
    "\n",
    "\n",
    "sns.lineplot(x = \"Noise STD\", y = \"Infeasible X-Recovery L0 Err\", ci=68, err_style='bars', style=\"Method\",markers=True, dashes=False,markersize=11,linewidth=3,\n",
    "             data = df[[\"Experiment\", \"Method\",\"Noise STD\",\"Infeasible X-Recovery L0 Err\"]],hue='Method')\n",
    "plt.ylabel(r'Support Err $\\mathbf{\\|\\mathcal{S} - \\hat{\\mathcal{S}} \\|_0}$')\n",
    "plt.xlabel(r'Noise STD $\\sigma$')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}