OnWavelets.ipynb

{
 "cells": [
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "We demonstrate the use of PCA to compress the wavelets of light curves. Here we read in the wavelets from a file into\na numpy array called `Data` which has the shape (`Nsamps`, `Nfeats`), where `Nsamps` is the number of objects, and\n`Nfeats` is the number of features. We would like to represent this in terms of `ncomps << Nfeats`, and test how well\nwe can reconstruct `Data` from these compressed vectors."
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "import numpy as np\nfrom snmachine.snfeatures import WaveletFeatures\nimport snmachine as sm",
   "execution_count": 1,
   "outputs": [
    {
     "output_type": "stream",
     "text": "PyMultinest not found. If you would like to use, please install Mulitnest with 'sh install/multinest_install.sh'\n",
     "name": "stdout"
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "!which python",
   "execution_count": 35,
   "outputs": [
    {
     "output_type": "stream",
     "text": "/Users/rbiswas/soft/mypython3/envs/snmachineHG/bin/python\r\n",
     "name": "stdout"
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "%matplotlib inline\nimport matplotlib.pyplot as plt\n\n",
   "execution_count": 37,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "import seaborn as sns\nsns.set_style('whitegrid')\nsns.set_context('notebook')",
   "execution_count": 34,
   "outputs": [
    {
     "output_type": "error",
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'seaborn'",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-34-be7d6e56ef25>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpyplot\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0;32mimport\u001b[0m \u001b[0mseaborn\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0msns\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      5\u001b[0m \u001b[0msns\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_style\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'whitegrid'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0msns\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_context\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'notebook'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'seaborn'"
     ]
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "import os",
   "execution_count": 2,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "fname = os.path.join('/Users/rbiswas/doc/projects/DimRed/dimred/example_data/wavelets.npz')",
   "execution_count": 3,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "wavelets =  np.load(fname)",
   "execution_count": 4,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "Data, _ = wavelets['arr_0']",
   "execution_count": 5,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "Data.shape",
   "execution_count": 6,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 6,
     "data": {
      "text/plain": "(7848, 26400)"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Do the pca "
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "wf = WaveletFeatures()",
   "execution_count": 9,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "obj_names = np.array(list('objid_{}'.format(i) for i in range(7848)))",
   "execution_count": 10,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "# Normalized and use tolerance\nwfpcaNt = wf.extract_pca(obj_names, Data, normalize_variance=True, method='svd', tol=0.99)",
   "execution_count": 11,
   "outputs": [
    {
     "output_type": "stream",
     "text": "OUTPUT ROOT: None\n\nRunning PCA...\nTook 8.233584880828857 secs for normalization\nShape of reduced data matrix X (7848, 26400)\nTook 678.9785442352295 secs for svd\nU shape is  (7848, 7848)\nUsing number of components =  325\n shape of U is  (7848, 7848)\nTime for PCA 688.9222347736359\n",
     "name": "stdout"
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "wavsNt, valsNt, vecNt, MNt, sNt = wfpcaNt",
   "execution_count": 12,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "wavsNt.to_pandas().set_index('Object').values.shape",
   "execution_count": 19,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 19,
     "data": {
      "text/plain": "(7848, 325)"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "obj_names.shape",
   "execution_count": 13,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 13,
     "data": {
      "text/plain": "(7848,)"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "# Unnormalized & use Tolerance\nwfpcaUt = wf.extract_pca(obj_names, Data, normalize_variance=False, method='svd',\n                         tol=0.99)",
   "execution_count": 14,
   "outputs": [
    {
     "output_type": "stream",
     "text": "OUTPUT ROOT: None\n\nRunning PCA...\nTook 3.5742199420928955 secs for normalization\nShape of reduced data matrix X (7848, 26400)\nTook 3284.9300050735474 secs for svd\nU shape is  (7848, 7848)\nUsing number of components =  5\n shape of U is  (7848, 7848)\nTime for PCA 3288.6212899684906\n",
     "name": "stdout"
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "# Unnormalized & use same number of components\nwfpcaUn = wf.extract_pca(obj_names, Data, normalize_variance=False, method='svd',\n                         ncomp=325, tol=None)",
   "execution_count": 15,
   "outputs": [
    {
     "output_type": "stream",
     "text": "OUTPUT ROOT: None\n\nRunning PCA...\nTook 8.207170724868774 secs for normalization\nShape of reduced data matrix X (7848, 26400)\nTook 2461.414072036743 secs for svd\nU shape is  (7848, 7848)\nUsing number of components =  325\n shape of U is  (7848, 7848)\nTime for PCA 2472.3972289562225\n",
     "name": "stdout"
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "wavsUt, valsUt, vecUt, MUt, sUt = wfpcaUt",
   "execution_count": 26,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "wavsUn, valsUn, vecUn, MUn, sUn = wfpcaUn",
   "execution_count": 27,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "def delta_Data_len(Data, wav, vec, M, s):\n    Z = wav.to_pandas().set_index('Object').values\n    D = WaveletFeatures.reconstruct_datamatrix_lossy(Z, vec, M, s)\n    Delta = D - Data\n    return np.sum(Delta**2, axis=1)",
   "execution_count": 22,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "lNt = delta_Data_len(Data, wavsNt, vecNt, MNt, sNt)",
   "execution_count": 28,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "lUt = delta_Data_len(Data, wavsUt, vecUt, MUt, sUt)",
   "execution_count": 29,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "lUn = delta_Data_len(Data, wavsUn, vecUn, MUn, sUn)",
   "execution_count": 30,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "fig. ax = plt.subplots()",
   "execution_count": 47,
   "outputs": [
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAADYBJREFUeJzt3HGI33d9x/Hny8ROprWO5QRJou1YuhrKoO7oOoRZ0Y20fyT/FEmguEppwK0OZhE6HCr1rylDELJptolT0Fr9Qw+J5A9X6RAjudJZmpTALTpzROhZu/5TtGZ774/fT++4XHLf3v3uLt77+YDA7/v7fX6/e+fD3TO/fH/3+6WqkCRtf6/a6gEkSZvD4EtSEwZfkpow+JLUhMGXpCYMviQ1sWrwk3wuyXNJnrnC7Uny6SRzSZ5O8rbJjylJWq8hz/A/Dxy4yu13AfvGf44C/7T+sSRJk7Zq8KvqCeBnV1lyCPhCjZwC3pDkTZMaUJI0GTsn8Bi7gQtLjufH1/1k+cIkRxn9L4DXvva1f3TLLbdM4MtLUh9PPvnkT6tqai33nUTws8J1K35eQ1UdB44DTE9P1+zs7AS+vCT1keS/13rfSfyWzjywd8nxHuDiBB5XkjRBkwj+DPDe8W/r3AG8WFWXnc6RJG2tVU/pJPkycCewK8k88FHg1QBV9RngBHA3MAe8BLxvo4aVJK3dqsGvqiOr3F7AX01sIknShvCdtpLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDUxKPhJDiQ5l2QuycMr3P7mJI8neSrJ00nunvyokqT1WDX4SXYAx4C7gP3AkST7ly37O+CxqroNOAz846QHlSStz5Bn+LcDc1V1vqpeBh4FDi1bU8Drx5dvAC5ObkRJ0iQMCf5u4MKS4/nxdUt9DLg3yTxwAvjASg+U5GiS2SSzCwsLaxhXkrRWQ4KfFa6rZcdHgM9X1R7gbuCLSS577Ko6XlXTVTU9NTX1yqeVJK3ZkODPA3uXHO/h8lM29wOPAVTV94DXALsmMaAkaTKGBP80sC/JTUmuY/Si7MyyNT8G3gWQ5K2Mgu85G0m6hqwa/Kq6BDwInASeZfTbOGeSPJLk4HjZQ8ADSX4AfBm4r6qWn/aRJG2hnUMWVdUJRi/GLr3uI0sunwXePtnRJEmT5DttJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNDAp+kgNJziWZS/LwFda8J8nZJGeSfGmyY0qS1mvnaguS7ACOAX8GzAOnk8xU1dkla/YBfwu8vapeSPLGjRpYkrQ2Q57h3w7MVdX5qnoZeBQ4tGzNA8CxqnoBoKqem+yYkqT1GhL83cCFJcfz4+uWuhm4Ocl3k5xKcmClB0pyNMlsktmFhYW1TSxJWpMhwc8K19Wy453APuBO4AjwL0necNmdqo5X1XRVTU9NTb3SWSVJ6zAk+PPA3iXHe4CLK6z5RlX9sqp+CJxj9A+AJOkaMST4p4F9SW5Kch1wGJhZtubrwDsBkuxidIrn/CQHlSStz6rBr6pLwIPASeBZ4LGqOpPkkSQHx8tOAs8nOQs8Dnyoqp7fqKElSa9cqpafjt8c09PTNTs7uyVfW5J+UyV5sqqm13Jf32krSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSE4OCn+RAknNJ5pI8fJV19ySpJNOTG1GSNAmrBj/JDuAYcBewHziSZP8K664H/hr4/qSHlCSt35Bn+LcDc1V1vqpeBh4FDq2w7uPAJ4CfT3A+SdKEDAn+buDCkuP58XW/luQ2YG9VffNqD5TkaJLZJLMLCwuveFhJ0toNCX5WuK5+fWPyKuBTwEOrPVBVHa+q6aqanpqaGj6lJGndhgR/Hti75HgPcHHJ8fXArcB3kvwIuAOY8YVbSbq2DAn+aWBfkpuSXAccBmZ+dWNVvVhVu6rqxqq6ETgFHKyq2Q2ZWJK0JqsGv6ouAQ8CJ4Fngceq6kySR5Ic3OgBJUmTsXPIoqo6AZxYdt1HrrD2zvWPJUmaNN9pK0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqYlDwkxxIci7JXJKHV7j9g0nOJnk6ybeTvGXyo0qS1mPV4CfZARwD7gL2A0eS7F+27Clguqr+EPga8IlJDypJWp8hz/BvB+aq6nxVvQw8ChxauqCqHq+ql8aHp4A9kx1TkrReQ4K/G7iw5Hh+fN2V3A98a6UbkhxNMptkdmFhYfiUkqR1GxL8rHBdrbgwuReYBj650u1VdbyqpqtqempqaviUkqR12zlgzTywd8nxHuDi8kVJ3g18GHhHVf1iMuNJkiZlyDP808C+JDcluQ44DMwsXZDkNuCzwMGqem7yY0qS1mvV4FfVJeBB4CTwLPBYVZ1J8kiSg+NlnwReB3w1yX8mmbnCw0mStsiQUzpU1QngxLLrPrLk8rsnPJckacJ8p60kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNDAp+kgNJziWZS/LwCrf/VpKvjG//fpIbJz2oJGl9Vg1+kh3AMeAuYD9wJMn+ZcvuB16oqt8HPgX8/aQHlSStz5Bn+LcDc1V1vqpeBh4FDi1bcwj4t/HlrwHvSpLJjSlJWq+dA9bsBi4sOZ4H/vhKa6rqUpIXgd8Ffrp0UZKjwNHx4S+SPLOWobehXSzbq8bci0XuxSL3YtEfrPWOQ4K/0jP1WsMaquo4cBwgyWxVTQ/4+tuee7HIvVjkXixyLxYlmV3rfYec0pkH9i453gNcvNKaJDuBG4CfrXUoSdLkDQn+aWBfkpuSXAccBmaWrZkB/mJ8+R7g36vqsmf4kqSts+opnfE5+QeBk8AO4HNVdSbJI8BsVc0A/wp8Mckco2f2hwd87ePrmHu7cS8WuReL3ItF7sWiNe9FfCIuST34TltJasLgS1ITGx58P5Zh0YC9+GCSs0meTvLtJG/Zijk3w2p7sWTdPUkqybb9lbwhe5HkPePvjTNJvrTZM26WAT8jb07yeJKnxj8nd2/FnBstyeeSPHel9ypl5NPjfXo6ydsGPXBVbdgfRi/y/hfwe8B1wA+A/cvW/CXwmfHlw8BXNnKmrfozcC/eCfz2+PL7O+/FeN31wBPAKWB6q+fewu+LfcBTwO+Mj9+41XNv4V4cB94/vrwf+NFWz71Be/GnwNuAZ65w+93Atxi9B+oO4PtDHnejn+H7sQyLVt2Lqnq8ql4aH55i9J6H7WjI9wXAx4FPAD/fzOE22ZC9eAA4VlUvAFTVc5s842YZshcFvH58+QYuf0/QtlBVT3D19zIdAr5QI6eANyR502qPu9HBX+ljGXZfaU1VXQJ+9bEM282QvVjqfkb/gm9Hq+5FktuAvVX1zc0cbAsM+b64Gbg5yXeTnEpyYNOm21xD9uJjwL1J5oETwAc2Z7RrzivtCTDsoxXWY2Ify7ANDP57JrkXmAbesaETbZ2r7kWSVzH61NX7NmugLTTk+2Ino9M6dzL6X99/JLm1qv5ng2fbbEP24gjw+ar6hyR/wuj9P7dW1f9t/HjXlDV1c6Of4fuxDIuG7AVJ3g18GDhYVb/YpNk222p7cT1wK/CdJD9idI5yZpu+cDv0Z+QbVfXLqvohcI7RPwDbzZC9uB94DKCqvge8htEHq3UzqCfLbXTw/ViGRavuxfg0xmcZxX67nqeFVfaiql6sql1VdWNV3cjo9YyDVbXmD426hg35Gfk6oxf0SbKL0Sme85s65eYYshc/Bt4FkOStjIK/sKlTXhtmgPeOf1vnDuDFqvrJanfa0FM6tXEfy/AbZ+BefBJ4HfDV8evWP66qg1s29AYZuBctDNyLk8CfJzkL/C/woap6fuum3hgD9+Ih4J+T/A2jUxj3bccniEm+zOgU3q7x6xUfBV4NUFWfYfT6xd3AHPAS8L5Bj7sN90qStALfaStJTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ18f+GmWq6NWLIwgAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "fig, ax = plt.subplots(1, 3, sharey=True, squeeze=True)\nax[0].plot(lNt, label='normalized var, tol=0.99, 325 comp')\nax[0].set_title('normed, tol=0.99')\nax[1].plot(lUt,label='Un-normalized var, tol=0.99, 325 comp')\nax[1].set_title('un-normed, tol=0.99')\nax[2].plot(lUn, label='Un-normalized var, 325 comp')\nax[2].set_ylim(0., 1.e11)\nax[2].set_title('un-normed, 325 comp')\n\n",
   "execution_count": 57,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 57,
     "data": {
      "text/plain": "Text(0.5, 1.0, 'un-normed, 325 comp')"
     },
     "metadata": {}
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 3 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEICAYAAACwDehOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcHGWdx/HPj4QEIVwxQTEBgoIiortoAJXVxYtTYVdRYXEliqK76rIrHkEUXQQVXMUDREAwXsilQoRguCEQApkg5CIhk3uSkJncCUlmMslv/6iapqane7pnprqruur7fr3mNd1VTz/1VD/Vz6+epy5zd0RERAB2S7oAIiKSHgoKIiJSoKAgIiIFCgoiIlKgoCAiIgUKCiIiUqCgUIaZTTCzyxotb4mPmbmZHdZoeUt3+q77RkEhBvVq5M3sfWY2z8y2mtnDZnZIL2nfaWZPm9lmM5tpZv8UmWdmdrGZLTOzTWZ2i5ntU+vyZ1k9Gp6w3q4ws7Xh35VmZr2kLVvHZjbKzO4ys3Vm1mJmn69l2RtF+LtqC7+z58zsjMi808zscTPbYGYvmtkNZrZ3ZP4EM+swsy2Rv0HJrEn/NWRQMLPBSZeh3sxsBPBn4FvAcKAJuLVM2uHAROCHwH7AlcBfzWz/MMkngX8HjgdeA7wC+Hktyx+3PG4DwPnAvwD/ALwF+CDwuTJpK9Xx74HFwKuA04Dvmdl7alPs8lJYjxcAB7r7PgTf9+/N7MBw3r7AZQTf5xuB0QS/sagr3X1Y5G9nvQoeG3ev2x+wBPgKMBPYSNCo7RGZ/1mgGVhH0Ki9JjLPgS8AC4DFkWn/GU7bDHwXeB3wJLAJuA0YEsnjg8CzwAZgKvCWyLyjgWfCfG4FbgEuq2Kdzgd2AB3AFuCv4fQ3Ao+Ey5oDnB75zIRq8i6xnKmR93sB24AjSqT9IDCnaNoLwHnh6zuAr0bmvRPYDuxZIi8HDitVduAEoAW4EGgFVgGf6mUdek1P8KP7LdAGLAW+CewWzhsHPAFcFW4flxVN2wAsCtdlHLA8XMa5kfyHAv8HLANWA78EXhGZ/9WwTCuBTxevey/r9ViY9qVwG/h4ldtzxbyLljMVOD/y/jxgWpm0xXXsQDuwJzAsfP/jSL1sIdjRyG09lljPYwl+F8eWmf9hYFZ/f9fh9vE8QZszF3hrlW3HL4B7wzp7Ang18BNgPTAPODqSfglwUZj/euDXRNrckuXq6xc1kL+wgE8TRNrh4Rfy+XDee4E1wFvDSv858FjRRn1/+LlXRKZNBPYB3kSw0T8IvDbcMOd2bUxhvq3AccAg4NywPEOBIeHG+z/A7sCZBA19VRVcvDGEeTQD3wjzfm9Y8W8oTg8cHFZ+ub9/C9P9FLi2aLmzgY+UKM+HgLlF0xYAV4Wv/wR8LTLv+PC7/IcSeVUKCp3ApeE6nwpsBfbvpTEpm56gIbkL2BsYQ/dANi787JeAwQR7vl3TPhXW6WUEDcU1Yb2eGH7vw8I8fhJuL8PDZfwV+H4472SCBuYogoB7c/G6V9gGir+narbnw8LX43vbBiKf2QgcF3k/FthcpjzFdexddRyue3FQ2EXQkOa6HsM87iYIBg78jTCglUj3E+CWot/GuvBvBiV+m5G0HwVWAMcABhwGHEJ1bcca4G3AHsBDBD2+T0a+u4eL2tzZwEHh9/UEFdq1JILCJyLvZxLs7c4GbiToenXNG0bQMI+JbNQvhBvPmZFpx4cV1/Uj+nokjx8BPwlfXwt8t6g884F/Bt5NsFdhkXlTK315RRtDNCi8C3gxujEBfwS+Uyp9lcu4EfhB0bQngHEl0r4y/C7ODjeycwl+9NeF8z8TfpdjCILnxPC7fEeJvCoFhW3A4Mj8VuDtvTQmJdOHG3Q7cGRk3ueAR/zlxmRZUX7jgAWR928Oy/uqyLS1wD8S/PBeAl4XmfcOXu513hT9foHXF697hfop/p6q2Z772lPYSaRnCBwe5mMl0hbXcVdQeEc4fzUwnaBh+Ww4b37e6zHyud2BU4D/KTP/AwR73q+PTHsrwW9vMEGg3AwcX+bzk4ELSkyvpu24ITLvS8DzRd9ddEdiCeGOd/j+VGBhb+uexHjei5HXMwi6h4cQ9B6e6Zrh7lvMbC0wimDFAL4MfKwov9UE43p7Ar8J33fZRtC1IlzGuWb2pcj8IeFyHVjh4bcWWtrXFYt4DbDc3XcV5TdqAHluIegRRe1DsOF14+5rwwNk/0ewtzUZeICgyw/BD+cggi7qYILg+aHI/L5Y6+6dkfdbgWFmdjBBT62rTMN6Sw+M4OUeW5fi72x5ieUX1zfuXjxtGDCSYBuZETk2awSNGAR1NqNo2QNRzfbcV8XbwD7AlqLttktxHXfpquPHCBrZ5QTDPC/Rvf5zXY/uvgO418wuMLOF7j6xsDCztxP0QM509xcin3kmksUkM/sDwRDTEyUWcRCwsMT0atqO4u+l1PcUFf2+l4bLKCvpA81LCcbiIdhTP8TMXmdmfzOzvwMHEOzJdJlHsMfbjbs/SInGschy4HJ33y/yt6e7/5Gg2zyq6EyOg/uwHsU/ypXAQWYW/X4PJugudmNmBxedrVD8d06YdA5B17/rc3sRHD+ZU7JA7o+6+zHuPpzggOMbCIbucPdd7v5tdx/j7qPDPFaUKh/Bj33PyPtXl0hTavnLPHLArYqPrCHYkz4kMq34OyvV+FVrDcEP5k2R+t83UrZVBD/U6LIHYiWRdQnr65WU3ga+0ds2EEnabRsIX5er/+I6bifYm+9a/t7Are4+kuC43G6E20dRPnmvx8EEvzMAzOxogp71p8N2pzdOELBKWR7NN6LqtqMPir+Plb0lTjooRN1MMKZ4M0GP4DGC4aWLY8r/BuDzZnZceLreXuEpZnsTHJjuBP7LzAab2YcJDjIVhKccnlAm79UExzG6PEWw5/U1M9s9/NyHCA5ed1P8oyvx94cw6V+Ao8zsI2a2B3AJMNPd55UqkJkdHS57H4IeQ4u7Tw7nDQ+Dr5nZkcCPgUuL9k66PAv8m5kNMrOTCYbbYufBWRq3AZeb2d4WnG77ZYKzZOLIfxfBNnCVmR0AhdMyTwqT3AaMM7MjzWxP4NvRz5vZODNb0ssiireBm4FPmdk/mtlQ4HvAU+7eIw93/15v20Ak6W+BL4flfg3Bgd4JpQpToo53EvRcLKzHE4AhZjaEYChkKMF2MCCNXI9mdoSZnWJmrwh/O58gGFp+NJx/FMFQ9Zfc/a8lPn+mmQ0zs93M7ETgEwQBpJRfAV8xs7eFdXRY+F1V3Xb0wRfMbHR4VuI3KHPWYpfUBIUw6l5GcOBlJsE44+7Agb18rC/5NxGMnV5NMBbYHC4Dd+8g6OaNC+d9nOD0TwDMbDRB131WmexvBI604PzlO8P8TicYk1xDcLbAJ8s14FWWvw34CHB5WMbjgLMiZfylmf0y8pGvhcteTvAd/mtk3ghgEsHGdy9wk7tfX2bRFxBslBuAc4A7+7sOVfhSWKZFwOMEDetNMeb/dYJ6n2ZmmwiG1N4A4O73Ehw4fChM81DRZw+i9DBAl+8Avwm3gY+F2/O3CA74riLYKzyrl89X4zqCg6qzCI7D3RNOAyDsWbwrfFtcx9fx8kkN5xD0ML5AsC2dDqwJt7E4NGo9GkE9thIMqV1AcCZZ17DQhQTDVzdGenLRntoFBHv0GwiGtD/r7o+UWpC7307wW76ZYJTjTmB4LdqOcBn3EdTHIoJ2tiwrPRxZP2Y2Brjb3Y8K92rnu3vZQGBmE8L0dxRNPwH4irt/sAZl/ARBd/WiuPOWxmBm9xEcGHw+6bJI/+WtHsNe0Wfc/YFqP5OqC0fcfZOZLTazj7r77WZmBNcSPJdwuWLp+krjcvcTky6DDJzqsbKKw0dmdpOZtZrZ7DLzzcx+ZmbNFtxO4a3VLtzM/kgwnv8GCy61P4+ga3uemT1H0MU9I0x7jJm1EJzfe12022ZmU4DbgfeF+ZxUvCwREams4vCRmb2bYDz9t+5+VIn5pxKMIZ5KMM79U3c/rgZlFRGRGqvYU3D3xwiu0CvnDIKA4e4+DdjPXr5XiIiINJA4jimMovvFES3htFXFCc3sfIJ7+LDXXnu97Ygjjohh8TJQM2bMWBOer94vlep1/dYOWtZv6/G5vYcOZsyIvZi1YmNh2ptH7Vt4/+ZR+wJ0mx9N1+XFTdtp29xeeH/oiL1YvOYl9hwyiK0d5e9HFl1Wf+cvWL2F7Z07OfyAYeyx+6Be0wMM3s3o3FW6d/6m1+zDnJWbev18sZF7DwVgzeZ2jop8J1D7epVkDLReK4kjKJS6OKPkVh+e9ng9wNixY72pqSmGxctAmdmArt6tVK9/mtHChbf3PFfgvUccwE3jjmHM+HsK05p+cFrhfdMPTgPoNj+arssVf5vHtY+8fHHojZ85jnN+9RRjD9mfpqXry5Y7uqz+zj/pqseYv3ozd/z3uzji1fv0mh7ggL2H0hoJYFFT/vck3vTtyb1+vth/nvA6HLhxymKaLj+l27xa16skY6D1Wkkc1ym00P2KudFUuGJORETSKY6gMBH4ZHgW0tuBje7eY+hIRETSr+LwUXja6AnAiPCU0G8TXGmMu/+S4KrJUwmuHtxKcKsKERFpQBWDgrufXWG+E1wuLyIiDS419z4SEZHkKSiIiEiBgoJIDBK+r6RIbBQURAbAyj1CRaRBKShIZmnnXaTvFBQkc7p23mf0cjWziJSmoCDSwNQbkrgpKEjNadxdpHEoKIiISIGCgoiIFCgoiNSRhtIk7RQURGIQx8VrChiSBgoKkipTF65JuggiuaagIKnSuqn0U8lEpD4UFEREpEBBQbJHY/Mi/aagICIiBQoKDay5dTO7dmX3Rgetm7YzZvw9SRdDJFcUFBrU7BUbef+PH+P6KYuSLkrNzF+9uX8fzG6c7EHPcZC4KSg0qOXrtgLw7LINCZekMp1/L9I4FBREYuB56p5IpikoiAyAqRskGaOgICIiBQoKIiJSoKAgIiIFCgoidbRa93aSlFNQkOzRsV+RflNQEBGRAgUFkRjE8pAddXEkBRQUJLUaoZFMfwlF+kZBQaSB6UpqiZuCgtRcI+zxi0hAQUFERAoUFCS1NDQiUn8KCiIiUlBVUDCzk81svpk1m9n4EvMPNrOHzezvZjbTzE6Nv6giIlJrFYOCmQ0CrgFOAY4EzjazI4uSfRO4zd2PBs4CfhF3QUWqpQPbIv1XTU/hWKDZ3Re5ewdwC3BGURoH9glf7wusjK+IIiJSL9UEhVHA8sj7lnBa1HeAT5hZCzAJ+FKpjMzsfDNrMrOmtra2fhRX0ijOem20Z9Y0Wnn7Qr/XfKomKJTa7ItPCzkbmODuo4FTgd+ZWY+83f16dx/r7mNHjhzZ99JKD2k4Q0f1mk2q13yqJii0AAdF3o+m5/DQecBtAO7+JLAHMCKOAkppWd5D7aJjAyL1V01QmA4cbmaHmtkQggPJE4vSLAPeB2BmbyQICupviog0mIpBwd07gS8Ck4HnCc4ymmNml5rZ6WGyC4HPmtlzwB+Bce5x3DdSpO+qHVJrWb+1xiWpA/3KJGaDq0nk7pMIDiBHp10SeT0XOD7eomVLV4y0PIz7FOnvKtf6eMk/XfFwTfMXaUS6orlOvnDzMxx60aTKCUVEEqSgUCeTZr2YdBEk5XLYiZQUUlDIsEfmt7Jha0fSxai7JM5a0hE0yQoFhYzauG0H4349nfN+05R0UTJNe/eSNQoKGdW5cxcAi9e8lHBJykv7CWpTFuisaskfBQWRMv79xqeTLoJI3SkoSGLatrTn8piHSJopKEhiZq/YxNHfvb/s/Lzd5iLlo2mSEwoKkig1hCLpoqAgqZWGO8CK5I2CgkgMkgpgCpsSNwUFkQHI23EPyT4FhYxL+7UAIpIuCgoZlce7sXZJYtX/4/fP0Ny6uf4LFomZgoJIDFZs2MZFf56VdDFEBkxBQUREChQURESkQEGhwWX5OLLO7BGpPwWFhqUGU0Tip6AgqaUrmkXqT0FBZACip79meShP8kNBIePUTolIXygoZJSOOOSDrliXuCkoSMNTABSJj4KCiIgUKCiIiEiBgoJITDS6L1mQu6Dg7mzfsTPpYgCwcsM2Nm3fkXQxREQKchcUfjN1CUd862+s2rgt6aLwzh88xMlXPZZ0MVJLt7kQqb/cBYV7Zq0CYPm65IMCwMqN25MuQmo1whXNCluSNbkLCnmj09gbRyMEQck+BYWMyvGD10RkABQURESkQEFBMiepTpJuOSFZUFVQMLOTzWy+mTWb2fgyaT5mZnPNbI6Z3RxvMUWkFMUhidvgSgnMbBBwDfABoAWYbmYT3X1uJM3hwEXA8e6+3swOqFWBpfGYDnCINIxqegrHAs3uvsjdO4BbgDOK0nwWuMbd1wO4e2u8xZSetIsoIvGrJiiMApZH3reE06JeD7zezJ4ws2lmdnKpjMzsfDNrMrOmtra2/pV4gLLW3U7DTnga6lXip3rNp2qCQqlmp7hpHQwcDpwAnA38ysz26/Eh9+vdfay7jx05cmRfyxqrNDSmcUhDkEtTvdZdVjakEnJdrzlWTVBoAQ6KvB8NrCyR5i533+Hui4H5BEFCaqa6xqiRz4hptKI3WHFFSqomKEwHDjezQ81sCHAWMLEozZ3AewDMbATBcNKiOAsat0ZrcPpK9w0Skf6oGBTcvRP4IjAZeB64zd3nmNmlZnZ6mGwysNbM5gIPA19197W1KvRApLW3P6tlI2PG38PTi9clXRRJSNZ3VKQxVDwlFcDdJwGTiqZdEnntwJfDP+mHx5vXAPDgvNUce+jwhEsjInmlK5pFRKRAQUEyRxfLifSfgoKIiBQoKEjmzF+9uW7LivZJkjhQrGPTErfcBQWd4VF/tRjMuebh5rLzvnXn7BosUSQfchcUuuRl2LmRY2BvZf/h5Pl1K0e97Ni5K+kiiOQ3KGReToJellxy15ykiyCioCDplbe49tTiVF7vKTmT26CgYwsiIj3lLijk5ViC1J/2MyQLchcURNJKvVdJAwUFEREpUFAQGQANR0rW5C4oNGoX/bK753LiVY8mXQwRybjcBYUupfbwHpnfSmdKLyD61eOLeWH1lqSLISIZl9ugUOzRF9oY9+vpXN3L7ROkvhqhUxfteT63fMPA8hrg8kXioKAQatvcDsCydVsTLkn2aNxdpHEoKDS4LO8oKpaI1J+CQoPS3reI1IKCQtpkeddfRFIvt0EhbQfotOfff1vaO5MugkhmDE66APWWh8b3tunL2XPooKSLUTfNrdk4VTdtOyqST7kLCnn44X3tTzOTLkJueB42KMmV3A4fZa3HcN2jCxkz/h42b9+RdFFEpIHlNiiU1aA7fn94ahkA617q6D6jQddHRJKhoBCqZcdh+46dNcw9W+6Y0VJ4rXgmUn8KCnVw3PceTLoIDWPKgjVJF0Ek1xQU6mDjNo3zS224+lMSMwUFqTnr5+BctZ/K2kkDIklSUBAZCEUkyRgFBRERKchtUNA1RyIiPeU2KIiISE+5DQrlhoKz1oHI2vpkm2pLkpfboFBMxwtFRKoMCmZ2spnNN7NmMxvfS7ozzczNbGx8RZS8qna/ub+nvIpITxWDgpkNAq4BTgGOBM42syNLpNsb+C/gqbgLGRd3p2np+qSLISKSWtX0FI4Fmt19kbt3ALcAZ5RI913gSmB7jOWL1SMvtCVdhAHTfZSya82WjsqJRGqsmqAwClgeed8STisws6OBg9z97t4yMrPzzazJzJra2urfQLfv2FX3Zcbt7pmrur1Pw6m1SddrkrI8cJXnes2zaoJCqe2+0BSZ2W7AVcCFlTJy9+vdfay7jx05cmT1pZQe0tQY1ape07SOeaTfaz5VExRagIMi70cDKyPv9waOAh4xsyXA24GJOtgsUntp6ClKtlQTFKYDh5vZoWY2BDgLmNg10903uvsIdx/j7mOAacDp7t5UkxI3gO/f+zxnXf9k1ek7d+5i567s/rp1um/tTJi6hNWbUnsYTxpQxWc0u3unmX0RmAwMAm5y9zlmdinQ5O4Te88hncq1U3E8c/e6Rxf1Kf1hF9874GWWk4dnCOc96CxYvSXpIkiGVAwKAO4+CZhUNO2SMmlPGHixaq+4qRxow7Jzl/Pdu+fy2Xe/dmAZSUPJe0CS7KkqKEhlTUvWMWHqEua9uCnWfJ9cuJb75r4Ya54iIuUoKPRD2+Z2Onft4sB9X1GY1tXziPvQwNk3TKsq3bJ1W+NdcApkf+ArHnr6msRJQaEfjrn8AQCW/OC0hEsiIhKvXN0Qr97jvw/Pa+Wl9s76LlTqSocUJGvUU6ihT02YzmlvPjDpYkjG5eAEM6mjXPUUkrB4zUtJF0FEpGoKCpJaGpoRqb/cBoW8NDgaWRCRvshtUGhUebhCWfpGm4TEKbNBYc7KjazYsC3pYkgdmC4rFolNZs8+Ou1njwO6lkBEpC8y21OoRD3u9FMdidRfboNCo1ODKV20LUicchUU6jLyXOdfqIbTk9Wxs/Ef8SoSlYmgMG3RWsaMv4eFbcndV15tcz7NXhHvXXFFkpaJoHDXs8HTQactWptwSepny3bdU0kCOk1Z4pSJoNAfDzy/mglPLO4xvVF+Xl//00yAis9vUHshIn2RqVNS+9IAdj0yc9zxhwJgKRkAqnYVOsMHN2zr0Ji2iMQnEz0FHWwVEYlHJoJCvml8SPsEIvHJVFCIo3mcs7I+Z5MsXVv6ltrrXuqgvXNnXcogIlIsE0Ehzj3F5tb6nNbadaC42B0zWjhvQlNdylAv/a0fnVVTnXkvbk66CJIhmQgKBRlpRB5vXtOn9OpZiEhcMhEUqj3Q3Je7aW7t6KRtc3s/S1Q/7nQr5/OrtNcoIv2XiaBQCx/6+eMcc/kDSRejzz7/+xlJF0FEGli2rlPox2fOuPpxhu0xmI++7aBu0xe26dnKIpI/megpDOTCs+daNvJEc35uj5FFuk5FJD6ZCApxqHfDEtcxcSe7Tx6rZr3cnUfmt9WhNCL5kKmgkJGTjwpO+9mUpIuQqD9MW1oxzQurk7szrkgWZSIoVLuj/EQfT/VMWr0upEur++aurpjGdUW3SKwyERS6VLrYacLUJfUpSB3pAq/k6LuXLMpEUGjEEfU4DwNUyuqGxxbxlu9Mjm+BKbF9R7IX7V14+3OJLl+kFjJ1SqqUdvmk55MuQk3Mf3Eze+w+KLHl//mZFYktW6RWMtFT6NJbZ15d/eTU6uQo1ahI/KoKCmZ2spnNN7NmMxtfYv6XzWyumc00swfN7JD4i9pr+Sqm6XpkZ1rEeUqqiEhcKgYFMxsEXAOcAhwJnG1mRxYl+zsw1t3fAtwBXDmQQrVu3h77nv2iNdm9Qrm3mJjls3PU+xOJXzU9hWOBZndf5O4dwC3AGdEE7v6wu28N304DRve3QLNXbOTYyx/k1unL+5uFpMzkOZVPLe2vjF63J5KYaoLCKCDaQreE08o5D7i31AwzO9/Mmsysqa2t9FWoXc8zeHJROm89MWflxros5/J75pac/tU7uj+HodLO8vYdtX+Gc6V6vWfmqpqXQeJXze9VsqeaoFBqX6xkU2RmnwDGAj8sNd/dr3f3se4+duTIkdWXMgZx7FBuae/ktJ89PqA8nl68rqp0N0xZXFW6NAyhJFWvya95tiX5e5XkVBMUWoDoLURHAz2O2prZ+4GLgdPdfcAPIuhPW1fr9vHDv3hiwHl87LonmTznxRhK87KB3BCwLmpYvJSved107Kx9j1DyoZqgMB043MwONbMhwFnAxGgCMzsauI4gILQOpEBpHiOu5j471RzYXbF+WxzFyb0UdJJEMqdiUHD3TuCLwGTgeeA2d59jZpea2elhsh8Cw4DbzexZM5tYJrua6q2NqHWwSepOpU767+lUq2/mI9dO5dEXNNYtEqeqrmh290nApKJpl0Revz/mcvVpvLia9rgWe5VrtrQzYtjQ+DPug83bO3N9u4XfPln5TqoiUr1MXdE8EBf9eVafPzP2sv49rjNvox5pHhIUke4yFRR6OxOnUsO0tePlm6vNaqnPaad5UcsD4cvWba2cSESqltqg0JdTLeNudD509cBOO5Xu1FMQaRypCwqN/mjJ6UvWJ12E1In2wkQk3VIXFGolbefy5+24gog0htQ+T6HaRvMLNz/DtIXpvCVGOQvbtigqiEgqpa6nUGl/fvWm7bRu3l54f8/MVax9qQNI5mKmO//e9wet3PzUMp5eUt3tLuL05m9P5obHFtV9uSLSOFIXFCo57nsPcuzlDyZdjIL/vvXZpItQtc3tnXV/Clsa7s0kItVL1fCRu3PHjJaki9FDe+dOtnfo3jL9MfG5dD3cSER6l6qg8PTidVXftqC9cydDB3d/Pm+tHihzzg1P0bRUZxX1x+IMP9xIJItSNXx0w5Tqx7u/8efZfcp7IGe6KiD0X0enelgijSQ1QWHj1h088HzkBqsVdvqfXLiGpWu1F5p2v3hkYdJFEJE+SE1Q6NzVc49yxtL13PL0srKf+ecfPtLtvY5piogMTGqOKewq0aB/5NqpAJx17MF1Lo2ISD6lpqdQfJC40kHjUrfDUEdBRGRg0hMUYmzRt+leOyIi/ZKpoGDAlAVtvPGSvw08MxGRHEpPUCga/Jk0q+8Pt3dgapn7IKXrdngDc63O6BGRGklPUKiipzBm/D39zr+z1JHsOmhu2xJ7nlf8bV7seYqIQIqCwq6YDiqU6xH89MEFseTfV7p4S0QaSWqCQl9jQqkrlEvloRuyiYhUL1NBodT0Tds7+1cgEZEcSk1QiGP4yPHUPWFNRKSRpOaK5loN8lz8l1lMntP3M5lERPIoPUGhRmP/d89cVZN8RUSyKD1BoY/pSw0TJXVM+ffTlrLnkEGVE4qIpFx6gkIfW/Rl67aWnJ7Ek76+eWffnu0gIpJWqTnQHMde/sK2LWWDRX9tadfZSyKSH6kJCnFccLy1Pf4b4R317cmx5ykiklYpCgoDjwpPLip93yMREalOaoJCHMNHG7ftGHgmIiI5lp6goEfkiIgkLjVBoWX9tqSLICKSe6kICqs2buNzv5tRdv62jp18/Y6ZdSyRiEg+peI6hbVbOnqd/64rH2bI6/hBAAAFbklEQVTNlvY6lUbi0t6px6KKNJqqegpmdrKZzTezZjMbX2L+UDO7NZz/lJmN6Ush9ti992IoIDQm3bU827Z2dLJkzUvdpnV07mJXeH75kwvX8tC81dw9cyUvbtzOr59Y3O26n4fntfLQvNUsD68t6u/vfMPWDp5bvqHH9O07djKrZWOvn23v3Mn5v23i4fmtLF0brMt3Js7hpKseA+DTE6bz8eueZFHbFmYsXc/3732+8NmnF6/jV1MW8dSitYV1KLZg9Wb+vmx94f2UBW3c9eyKwvtJs1Yxe8VGduzcxZQFbYXpnTt30bkzmWexVOwpmNkg4BrgA0ALMN3MJrr73Eiy84D17n6YmZ0FXAF8vNpCDB2sW0Rk0aDddMfarIo+BfFP//FOPnLt1G7z33X4CKYsWNPjc//717kcNPwV3HTuMXxqwvQe889822jumNHCRaccwQ1TFhcCxdDBu9HeuYsvvucwrn64uZD+6IP34+/LXg4IF37g9Rywz1BOftOBHPf9B9i+o3TDetgBw3jrwftxW1MLAPfNXQ3AwcP3LFwAe/Sl97F+a3BG43t/9Gjhs9c9uohX77MHL27a3iPfV+41BAfWvdT76McFtzxbcvo5xx3MH55aVnh/+b8exTnHHdJrXnGrZvjoWKDZ3RcBmNktwBlANCicAXwnfH0HcLWZmVd574rBg9R4ZNGgcg+9kEwpDghAyYDQZfm6bXwg3BMvdseMoJH+/r3dHznbHj7BMBoQgG4BAeBH978AwNf/NKvXMje3bqG5teejcqN3ROgKCKWUCggAaysEg0qiAQHg4r/M5n1HvIpX77vHgPLti2qCwihgeeR9C3BcuTTu3mlmG4FXAt22DDM7Hzg/fLvFzOYX5TOi+DMZlPg62hU9Jg1oV0T1CqRgHVWvNZH4Oh4Yc71WUk1QKLW7V9wDqCYN7n49cH3ZBZk1ufvYKsrUsLK4jqrXbK6j6jUf61ismgPNLcBBkfejgeJbkRbSmNlgYF9gXRwFFBGR+qkmKEwHDjezQ81sCHAWMLEozUTg3PD1mcBD1R5PEBGR9Kg4fBQeI/giMBkYBNzk7nPM7FKgyd0nAjcCvzOzZoIewln9LE/ZrmqG5GEdi+VhnfOwjsXysM55WMduTDv0IiLSJRW3uRARkXRQUBARkYJUBIVKt9FIOzNbYmazzOxZM2sKpw03s/vNbEH4f/9wupnZz8J1nWlmb43kc26YfoGZnVtueY1C9VrIR/WaIqrXCtw90T+Cg9cLgdcCQ4DngCOTLlcf12EJMKJo2pXA+PD1eOCK8PWpwL0E13a8HXgqnD4cWBT+3z98vX/S66Z6Vb2qXrNfr9G/NPQUCrfRcPcOoOs2Go3uDOA34evfAP8Smf5bD0wD9jOzA4GTgPvdfZ27rwfuB06ud6FjpHpVvTaSvNdrQRqCQqnbaIxKqCz95cB9ZjYjvDUAwKvcfRVA+P+AcHq59c3C9xCVhfVRvfaUhfVRvfYiDc9TqOoWGSl3vLuvNLMDgPvNbF4vacutbxa+h6gsrI/qtacsrI/qtRdp6ClUcxuNVHP3leH/VuAvBF3s1WE3k/B/a5i83Po2/PdQpOHXR/VaUsOvj+q1d2kICtXcRiO1zGwvM9u76zVwIjCb7rf+OBe4K3w9EfhkeFbD24GNYXd1MnCime0fnvlwYjitUaleVa+po3qtQtJHuiNH+F8gOKvh4qTL08eyv5bgDIzngDld5Se4dfiDwILw//BwuhE8tGghMAsYG8nr00Bz+PeppNdN9ap6Vb3mp167/nSbCxERKUjD8JGIiKSEgoKIiBQoKIiISIGCgoiIFCgoiIhIgYKCiIgUKCiIiEjB/wNOd46h10+RDgAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "fig, ax = plt.subplots(1, 3, sharey=True, squeeze=True)\nax[0].plot(lNt, label='normalized var, tol=0.99, 325 comp')\nax[0].set_title('normed, tol=0.99')\nax[1].plot(lUt,label='Un-normalized var, tol=0.99, 325 comp')\nax[1].set_title('un-normed, tol=0.99')\nax[2].plot(lUn, label='Un-normalized var, 325 comp')\nax[2].set_ylim(0., 1.e9)\nax[2].set_title('un-normed, 325 comp')\n\n\n#plt.legend()",
   "execution_count": 58,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 58,
     "data": {
      "text/plain": "Text(0.5, 1.0, 'un-normed, 325 comp')"
     },
     "metadata": {}
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 3 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEICAYAAACwDehOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcHHWd//HXJwkhQBIgJKAkQBCiEFEXzHKIB7srGEGDu7oKXkRBdH+i7MrqBg/kh4CA+xM8WCUKIiqEgKsECUTuSCCQCeQgF5kcMJNrch/knMzn90fVNDWd7umemeqpb/e8n4/HPKa76ttVn6pvdX3q+62jzd0REREB6JV1ACIiEg4lBRERyVFSEBGRHCUFERHJUVIQEZEcJQUREclRUijCzO40s2urbdqSHjNzMzu+2qYtbWldd4ySQgq6aydvZv9kZgvNbLuZPWlmx7RT9j1m9oKZbTWzOWb23sQ4M7PvmNlrZrbFzCaY2cBKx1/LumPHE9fbjWa2Pv67ycysnbJF69jMhprZA2a2wcwazewrlYy9WsTfq7XxOpttZucnxp1nZs+Y2SYzW21mvzKzAYnxd5rZbjPblvjrnc2SdF5VJgUz65N1DN3NzAYD/wt8DxgE1AH3Fik7CJgE/Ag4BLgJeNDMDo2LfB74HHAmcCRwAPCzSsaftp64DQCXAh8D3gW8E/gI8OUiZUvV8e+BZcARwHnA9Wb2D5UJu7gA6/Fy4M3uPpBoff/ezN4cjzsYuJZofZ4IDCP6jiXd5O79E397uyvw1Lh7t/0By4H/BOYAm4l2av0S478E1AMbiHZqRybGOfBVYDGwLDHs/8TDtgI/AI4DngO2ABOBvolpfASYBWwCngXemRh3MvBiPJ17gQnAtWUs06XAHmA3sA14MB5+IvBUPK95wJjEZ+4sZ9oF5vNs4v1BwA7ghAJlPwLMyxv2CnBx/Pp+4JuJce8BdgIHFpiWA8cXih04C2gErgCagFXAF9pZhnbLE33p7gLWAq8C3wV6xePGAtOAm+Pt49q8YZuApfGyjAUa4nlclJj+/sB/A68Ba4BfAgckxn8zjmkl8MX8ZW9nuabGZV+Pt4FPlbk9l5x23nyeBS5NvL8YmF6kbH4dO7ALOBDoH7//caJethEdaPTYeiywnKcSfS9OLTL+X4C5nf1ex9vHAqJ9znzglDL3Hf8DPBzX2TTgTcAtwEZgIXByovxy4Mp4+huB35DY5xaMq6Mrqit/cYAvEGXaQfEK+Uo87h+BdcApcaX/DJiat1E/Gn/ugMSwScBA4O1EG/3jwFviDXN+68YUT7cJOA3oDVwUx7M/0DfeeP8D2A/4BNGOvqwKzt8Y4mnUA9+Op/2PccW/Lb88cHRc+cX+Ph2X+wnwi7z5vgx8vEA8HwXm5w1bDNwcv/4j8K3EuDPjdfmuAtMqlRSagWviZT4X2A4c2s7OpGh5oh3JA8AAYDhtE9nY+LNfA/oQHfm2DvtCXKfXEu0obo3r9Zx4vfePp3FLvL0MiufxIPDDeNxooh3MSUQJ9+78ZS+xDeSvp3K25+Pj1+Pa2wYSn9kMnJZ4PwrYWiSe/Dr21jqOlz0/KbQQ7Uh7dD3G0/gLUTJw4BHihFag3C3AhLzvxob4byYFvpuJsv8KrAD+HjDgeOAYytt3rAPeDfQDniBq8X0+se6ezNvnvgwcFa+vaZTYr2WRFD6beD+H6Gj3ZeB2oqZX67j+RDvm4YmN+sX4M08RNd0cODPxmZnAfyXe/z/glvj1L4Af5MWzCPgA8H6iowpLjHu21MrL2xiSSeF9wOrkxgTcA1xdqHyZ87gduCFv2DRgbIGyhxHtUC6MN7KLiL70t8XjLyH6og4nSp6T4nV5RoFplUoKO4A+ifFNwOnt7EwKlo836F3AyMS4LwNPxa/HAq/lTW8ssDjx/h1xvEckhq0H/o7oi/c6cFxi3Bm80eq8I7l+gbfmL3uJ+slfT+Vszx1tKewl0TIERsTTsQJl8+u4NSmcEY9fA8wg2rF8KR63qKfXY+Jz+wEfBv6jyPiziY6835oYdgrRd68PUaLcSmL/lPf5KcDlBYaXs+/4VWLc14AFeesueSCxnPjAO35/LrCkvWXP4pzC6sTrmUSZDqLWw6utI9x9G9GGMDRRfoK7v5PoCOWH8bA1ifE7CrzvH78+BrgiPkm0ycw2EWXPI+O/FR6vtdirdN6RQIO7t+RNb2iR8uXYRtQiShpItOG14e7rgfOBbxCtj9HAY0RNfoi+OPcQJdd5wJPx8EY6br27Nyfebwf6m9nRyRNupcoDg3mjxdYqf501FJh/fn3j7oW2gSFEXSczE/X/SDwc4jrLm3dXlLM9d1T+NjAQ2Ja33bbKr+NWrXU8leh8UwNRC/l12tZ/j65Hd9/j7g8DHzKzMclxZnY6UQvkE+7+SuIzL7r7endvdvfJwB+IupgKOQpYUmB4OfuOcvd5rfLXx5FFYgKyP9H8KlFfPERH6seY2XFm9oiZvQQcTnQk0+rZ+P+TRDu9jmgArnP3QxJ/B7r7PUTN5qF5V3Ic3YFp538pVwJHmVly/R5N1FxsI/9LV+DvM3HReURN/9bPHUR0/mRe/jQB3P1pd/97dx9EdMLxbURdd7h7i7t/392Hu/uweBorCsVH9GU/MPH+TUXXQtv5v+aJE25lfGQd0ZH0MYlh+eus0M6vXOuIvjBvT9T/wYnYVhF9UZPz7oqVJJYlrq/DKLwNfLu9bSBRtM02EL8uVv/5dbyL6Gi+df4DgHvdfQjReblexNtH3nR6ej32IfqeAWBmJxO1rL/o7o+X+KwTtWwKaUhON6HsfUcH5K+Ple0VzjopJN1N1Kd4N9ER7lSirqLvJMp8KP7/z0QbdUf8CviKmZ0WX653UHyJ2QCiE9PNwNfNrI+Z/QvRSaac+JLDs4pMew3ReYxWzxMdeX3LzPaLP/dRopPXbeR/6Qr8/SEu+ifgJDP7uJn1A64C5rj7wkIBmdnJ8bwHEp2Ua3T3KfG4QXHyNTMbCfwYuCbv6KTVLODTZtbbzEYTdbelzqOrNCYC15nZAIsut/0G0VUyaUy/hWgbuNnMDofcZZmt29REYKyZjTSzA4HvJz9vZmPNbHk7s8jfBu4GvmBmf2dm+wPXA8+7+z7TcPfr29sGEkXvAr4Rx30k0YneOwsFU6CO9xJ1v1pcj2cBfc2sL1FXyP5E20GXVHM9mtkJZvZhMzsg/u58lqhr+el4/ElErZKvufuDBT7/CTPrb2a9zOwc4LNECaSQXwP/aWbvjuvo+Hhdlb3v6ICvmtkwi65K/DZFrlpsFUxSiLPutUQnXuYQ9TPuB7w5UezUuAXxATqYOd29jqjv9OdEfYH18Txw991Ezbyx8bhPEV3+CYCZDSNqus8tMvnbgZFxc/bP8fTGEPVJriO6WuDzxXbgZca/Fvg4cF0c42nABYkYf2lmv0x85FvxvBuI1uE/J8YNBiYTbXwPA3e4+/gis76caKPcBHwG+HNnl6EMX4tjWgo8Q7RjvSPF6f8XUb1PN7MtRF1qbwOIuwpuIerOrOeNbs1WRxGdwynmauC38TbwyXh7/h7RCd9VREeFF7Tz+XLcRnRSdS7RebiH4mEAxC2L98Vv8+v4Nt64qOEzRC2MrxJtS2OAdfE2loZqrUcjqscmoiunLie6kuzFePwVRN1UtydacsmW2uVE+6VNRJeqfsndnyo0I3e/j+i7fDdRF/CfgUGV2HfE8/grUX0sJdrPFmWFuyO7j5kNB/7i7ifFR7WL3P3NJT7TH1gYN4srLj5ieLu7X9kd85PwmNlfiU4MLsg6Fum8nlaPcavoEnd/rNzPBHXjiLtvMbNlZvav7n6fmRnRvQSzLbp5a0PcfLySdI88SsWVStNXqpe7n5N1DNJ1qsfSSnYfmdkdZtZkZi8XGW9m9lMzq7focQqnlDtzM7uHqD//bRbdan8xUdP2YjObTdTEbT2hfBawyMxeIboL87py5yMiIuUp2X1kZu8n6k+/y91PKjD+XKI+xHOJ+rl/4u6nVSBWERGpsJItBXefSnSHXjHnEyUMd/fpwCH2xrNCRESkiqRxTmEobW+OaIyHrcovaGaXEj3Dh4MOOujdJ5xwAgA79+xlcdO2/OJSIe8YenCb9zNnzlwXX6/eKcXq1R1eXrm5C5FKR3RXvUq2ulqvpaSRFArdnFGwTyq+7HE8wKhRo7yurg6A+Su3cO5P/5ZCKFKOuhvOa/PezLp0926xet25Zy8nfO+RrkxaOqC76lWy1dV6LSWN+xQaaXvH3DBK3DEnIiJhSiMpTAI+H1+FdDqw2d336ToSEZHwlew+ii8bPQsYbGaNRLeN7wfg7r8kumvyXKK7B7cTPapCRESqUMmk4O4XlhjvRLfLi4hIlQvm2UciIpI9JQUREclRUhARkZwgkoIV+xkKERHpVkEkBRERCYOSgoiI5CgpiIgUsXD1FuY29qzndwX1IzsiIiEZfUv0TLblec+VqmVqKYiIJLS0OHtbsv2Z4iwpKYiIJHzujuc57tuTsw4jM0oKIiIJ0+rXZx1CppQUREQkJ4ikoJvXRETCEERSEBGRMCgpiIhIjpKCiIjkKCmIiEiOkoKIiOQoKYiISI6SgoiI5CgpiIhIThBJwdDdayIiIQgiKUht0p3qItVHSUFERHKUFEREJEdJQUREcpQURKTHcnd+8dQSNm3fnXUowVBSEJEe67ml67nxkYVc+b9zsw4lGEoKItJj7dkb/Rbztl3NGUcSjiCSgi5dFBEJQxBJQUQkdC8s28DVk+ZlHUbFKSmIiJThk7c9x53PLs86jIpTUhARkRwlBRERyVFSEBGRnLKSgpmNNrNFZlZvZuMKjD/azJ40s5fMbI6ZnZt+qCIi4Zg0eyXfun921mGkrmRSMLPewK3Ah4GRwIVmNjKv2HeBie5+MnAB8D9pByoiEpKv3/MSE+sasw4jdeW0FE4F6t19qbvvBiYA5+eVcWBg/PpgYGV6IYqISHcpJykMBRoS7xvjYUlXA581s0ZgMvC1QhMys0vNrM7M6tauXfvG8I5ELMEpVq9S3VSvPVM5SaHQPtvz3l8I3Onuw4Bzgd+Z2T7Tdvfx7j7K3UcNGTKk49FKkFSvtUn12jOVkxQagaMS74exb/fQxcBEAHd/DugHDE4jQBGRWvHvE17ioz97Jusw2lVOUpgBjDCzY82sL9GJ5El5ZV4D/gnAzE4kSgpqb4qIJPx51krmrticdRjtKpkU3L0ZuAyYAiwgusponpldY2Zj4mJXAF8ys9nAPcBYd8/vYhIRkcD1KaeQu08mOoGcHHZV4vV84Mx0QxMRke6mO5pFpMfSlY/7UlIQEZEcJQURkYw8/cpaFq/ZmnUYbZR1TqHS9MtrItITXXTHCwAsv+G8jCN5g1oKIiKSo6QgFWM6jSdSdZQURCQTi9dspWnLzqzDkDxBnFMQkZ7n7Jun0ruXseR6/fxKSNRSEJHM7G3Rgw9Co6QgIiI5SgoiIpITSFLQVSoitei+ugYu+W1d1mFIB+hEs4hUzDfvn5N1CNJBgbQURKTWvLBsQ9YhSCcoKYhIRXzytueyDiFzl094iTE/D/uX1vKp+0hEpEIemJX/y8XhU0tBRHo8/U7kG5QURKTH0hOa96WkICIiOUoKIiKSE0RSUBNOpLat2LSD6x6aT4uedRS8IJKCiNS2r9/zEr/62zLmrNicdShBqlu+geHjHmLt1l1Zh6KkICKV17y3pWLT3t3cwvpt2e9Mu+LXf1sGRMkha0oKIlLVLrv7Rd597WNdmsYz9et4bP6alCKqbkoKIlLV/prSzvySu/TgPlBSEJGM7alg15J0nJKCiGTq7d+fknUIHeI1fvuzkoKIZGp3c2VaCg/MWsHwcQ+xbVdzRabfamJdA+OnLqnoPLqTkoKIVFwWx9Y/f6IegJWbdlR0Pt+6fw7XT16YyrRCaIMEkRR075pIz6DvemEh3cAbRFIQEZEwKCmIiEiOkoJUTEhNYklfS4uzY/ferMOQlCkpiEinfO+BlznxqkfYG/hD7mr8CtLUlZUUzGy0mS0ys3ozG1ekzCfNbL6ZzTOzu9MNU0RCc++MBgBaAt3rltNStU6c+q7k4oawKkv+RrOZ9QZuBc4GGoEZZjbJ3ecnyowArgTOdPeNZnZ4pQIWkeoTws6uK4aPe6ii0w+pq7WclsKpQL27L3X33cAE4Py8Ml8CbnX3jQDu3pRumCJSC0La+Ulh5SSFoUBD4n1jPCzprcBbzWyamU03s9GFJmRml5pZnZnVrV27tnMRS3BUr7UpzXptPe+wXSemg1dOUiiU2/Mbg32AEcBZwIXAr83skH0+5D7e3Ue5+6ghQ4a8MQMdPlS1YvUq1S3Nep2/agsAP/7rK2mEJhVUTlJoBI5KvB8GrCxQ5gF33+Puy4BFRElCRCSn0s8hKsSDeHhE9SgnKcwARpjZsWbWF7gAmJRX5s/APwCY2WCi7qSlaQYqImEK9SRyZ64syloICaxkUnD3ZuAyYAqwAJjo7vPM7BozGxMXmwKsN7P5wJPAN919faWCFpHspdnr6+6sqPCD60IWUgIreUkqgLtPBibnDbsq8dqBb8R/IiIdcv/MRr55/xzu/8oZjBo+KOtw2pX9sXxl6Y5mEek2xW50m/nqRgAWN23rznCkACUFEek2C1dvzTqENkLoww+NkoKI1LRQT4SHKoikEM4pFhGpFbr9qXOCSAoiUr3UBZOeEFo1Sgoi0imVuIwyhJ1iJgJq1SgpiEhqdjV37tlG1dTV492QuT7z6+n82+9nVnw+hSgpiEgqJs9dxdu++wgLV2/JOpSqN61+PQ+/vDqTeSspiEin7N7b0ub9YwvWAPD0orCelNvegf3CVWFdIhsCJQURSdUPH16YdQhlu27ygqxDCI6Sgoh0SY89OVwBIaxKJQURkYyFdJ49iKRQTVceiEjl5N/zcMlvZ/CV32VzFU5PVdZTUkVEKqvwkeFjC7r+c+8h3Fy3u7ml3YPfnXvC+ZlSJQURScVDc1ZlHUIbXf2Z32n16woO70yKOeF7DzPooL5Fx6eR/NISRPeRiFS/Xc0tpQtVkaXrXk9tWi0O67btTm16laSkIBWjU0UiHdMdd0uXoqQgIiI5SgoiEowADpR7PCUFEclcJS9LV6LpGCUFEalJOqfVOUEkhUo8l11EpNVTi5rYsTucewFCFkRSEJGeY+WmHd06v/qmbYz9zQy+86e5qUyv1rujlBREpFu954YnunV+W3fuAWDZ+vTuOyjlc7c/zyMZ/R5CV+mOZhHpkjSPnGvlIPxvi9fxt8WF74gOnVoKIpK5Sp5VrJVE012UFESkwz7z6+lZh7CP9dt2tXlf6DLXJxau4dVu7EaqRkoKItJhr67f3uZ9JU4e3/Xc8g6Vf/e1j7G8xPOKvnhnHR/40VOdjqknUFIQkS57bsn61Kd51QPzOvyZho3bi46rlquGPnbrtEznH0RS0I/siEiltO5eZjVsajM8hIfPFZIfZ3cLIimISHVJ7k/T+BGbpq3R+YAp3XgZ5x3TlnfqcyH8aE8lKSmISIetyDuH0NXd5IJVWwB4psgP23TFbVOXFhz+wKwVqc+rFigpiEhNa9zYvXdQVzslBRGRAnrqqU4lBRGRDMxt3Jx1CAWVlRTMbLSZLTKzejMb1065T5iZm9mo9EIUkdDlX8nTsKH4paHVYsbyDRWd/kd//sw+w0K4IKpkUjCz3sCtwIeBkcCFZjayQLkBwNeB59MOUkTCVejR9++76ckMIknXA7NWZh1CJsppKZwK1Lv7UnffDUwAzi9Q7gfATcDOFOMTkcClcYlmOfcqLV/3eqqPqMjynMGu5nB/26GcpDAUaEi8b4yH5ZjZycBR7v6X9iZkZpeaWZ2Z1a1du7bDwUqYVK+1qSP1OrWLTwQt54e2zvrvp0o+oiLZ/VKpm2LT6OK5YuLsrk+kQspJCoVWbW61mFkv4GbgilITcvfx7j7K3UcNGTKk/CglaKrX2tSRen1wds/saumsxxc0lVXu9V3NFY5kX+UkhUbgqMT7YUByCxgAnAQ8ZWbLgdOBSTrZLCJpeXlFulfq/HFmY6rTq5R/+8OL3T7PcpLCDGCEmR1rZn2BC4BJrSPdfbO7D3b34e4+HJgOjHH3uopELCI157USVyt95Gf7XqnTFVfcl373zd4WZ29LupcPzazwFVCFlEwK7t4MXAZMARYAE919npldY2ZjKh2giEg1OPOGJ3jH1VO6NI0QnqtU1s9xuvtkYHLesKuKlD2r62FJLTA9/lZ6kNVbauPCS93RLCJd0t03XD02fw0/eWxxyXIlr2jSQUtBZbUURERCccld0enKyz84IvVpN7bzIz1ZyKIzSS0FEZHYe2+s/juxuyqIpKBWnIh01PBxD+V+h0HSo+4jEalaf3qp+38o51Pjp/OOoQMrMu0QHoinpCAiVWt83q+qJfepleqBmN2widld/B3lHXuq+9lHIiJFBXBw2ynV0GudRctBSUFEgnL2j59mTRVf8z/yqkeyDqFLlBREJCiLm7ZxX11D6YKB2r473K6hcigpiIhIjpKCiEggQrj6SElBRGpSNZxILiWLB+QFkRT04DQR6W7a7RQWRFIQkerlFejzeGH5xtSnKeVRUhCR4Ex9Rb/1nRUlBRERyVFSEBEJlO5oFhFJSxWeSQ7gilQlBRGpTUvXbss6hKqkpCAiNWnrzuasQ6hKQSSF6mvkiUirELo8JD1BJAUREdmXfqNZREQypaQgIhKIStwd3lFKCiIikqOkICI9UjVc4JJFjEoKUjHV8KWT2hJC90uadjW3dPs8lRREpEuqdT+8c0/373CrgZKCiPRI81dtyTqEIAWRFKrwESUiIjUpiKQgIiJh3B2upCAiXfLyis1Zh7CPuY3hxVQtlBREpEtCTAqvbdiedQidEkJPupKCiHRJCF0eteKPLzZmHUJ5ScHMRpvZIjOrN7NxBcZ/w8zmm9kcM3vczI5JP1QRkdo2femGrEMonRTMrDdwK/BhYCRwoZmNzCv2EjDK3d8J3A/clHagIiJSeeW0FE4F6t19qbvvBiYA5ycLuPuT7t7aiTcdGJZumCISqhD6wfO5OrU6rZykMBRoSLxvjIcVczHwcKERZnapmdWZWd3atWvLj1KCpnqtTeXW6+otO7sxqvYpFXRdOUmh0IFAwXVvZp8FRgE/KjTe3ce7+yh3HzVkyJDEDEI81pByFatXqW7l1utvpi3vvqDKpH1K5/Upo0wjcFTi/TBgZX4hM/sg8B3gA+6+K53wRESkO5XTUpgBjDCzY82sL3ABMClZwMxOBm4Dxrh7U/phiohIdyiZFNy9GbgMmAIsACa6+zwzu8bMxsTFfgT0B+4zs1lmNqnI5EREJGDldB/h7pOByXnDrkq8/mDKcYmIdJquPuo83dEsIiI5QSSFvdX6Kx0iEpRla19nV/NeXX3UBWV1H1Xapu27sw5BRGrANX+Zz9wVm/ngiUdkHUrVCqKloKwuIml5dsm6rEOoakEkBRGRNOlEc+cpKYhIt9u5Z2/Fpq1TlF0TRFLQbzSL9CyXT3gp6xCkiCCSgoj0LM8tWV+xaTdt3aXzlF2gpCAiIjlKCiIikhNEUtA5hdqkepWsNLe0ZB1C1QoiKYiIpOm6hxZkHULVCiIp6KSQiKSpaat+0qWzgkgKIiISBiUFEel2W3Y2Zx2CFKGkICIiOUoKIiKSo6QgIiI5QSQFXc8uIhKGIJKCiIiEQUlBRERygkgK6j0SEQlDGElBWUFEJAhBJAUREQmDkoKIiOQoKYiISI6SgoiI5CgpiIhIjpKCiIjkBJIUdE2qiEgIgkgKuk9BRCQMQSQFEREJg5KCiIjkKCmIiEhOEElBpxRERApz926dX1lJwcxGm9kiM6s3s3EFxu9vZvfG4583s+FpByoi0hPt2du9SaFPqQJm1hu4FTgbaARmmNkkd5+fKHYxsNHdjzezC4AbgU9VImCRfN8970S+cOaxbNvVzLyVmzn+8P689NomNr6+mwP69ubyCbM4etCBvG/EYL78/uOYWNfA7MZNrNmyk1fWbONNA/vxuTOO4dOnHs0TC5s447jDOKx/X7bubGa/3r1o2rKTXzy9hDmNmznrrUM4rP/+nPCmATw4ZyXvGnYIy9a9TuPGHTy2YA23fvoUTnvLIP704go+d8YxTJm3msaNO+jTy3jTwf0Y2G8/nl2yjg2v7+HqMSMZ0G8/tu1q5qTvTwHg2+eewIjDB7Bmy04eX9jEo/PXAPB/x7ydIwbuz7BDD+Tf753FsYMPyo2b+d0Pclj//TNb/1JZ67bt4shDDui2+VmppomZnQFc7e4fit9fCeDuP0yUmRKXec7M+gCrgSHezsRHjRrldXV1ACxes5Wzb57a1WWRMi2/4bw2781spruPSmPayXrNt313Mwf2feM4pGHDdgb2248B/fowu3ETJx99KO7Onr1O3z5tG7F7W5wde/bSf/8+tLQ4vXqp0/G19dvZf79eHDGwX8HxlazX4eMeSmOyUobfX3wa7x0xOPc+zXotpGRLARgKNCTeNwKnFSvj7s1mthk4DFiXLGRmlwKXxm+3mdmivOkMzv9MDcp8Ge3GfQYd06XpqV4hzGVUvXZd5sv4vpS/r6WUkxQKHZLltwDKKYO7jwfGF52RWV0lM2AIanEZVa+1uYyq156xjPnKOdHcCByVeD8MWFmsTNx9dDCwIY0ARUSk+5STFGYAI8zsWDPrC1wATMorMwm4KH79CeCJ9s4niIhImEp2H8XnCC4DpgC9gTvcfZ6ZXQPUufsk4Hbgd2ZWT9RCuKCT8RRtqtaQnrCM+XrCMveEZczXE5a5JyxjGyWvPhIRkZ4jiDuaRUQkDEoKIiKSE0RSKPUYjdCZ2XIzm2tms8ysLh42yMweNbPF8f9D4+FmZj+Nl3WOmZ2SmM5FcfnFZnZRsflVC9Vrbjqq14CoXktw90z/iE5eLwHeAvQFZgMjs46rg8uwHBicN+wmYFz8ehxwY/z6XOBhons7Tgeej4cPApbG/w+NXx+a9bKpXlWvqtfar9fkXwgthVOBendf6u67gQnA+RnHlIbzgd/Gr38LfCwx/C6PTAcOMbM3Ax8CHnX0rEJ4AAABWklEQVT3De6+EXgUGN3dQadI9ap6rSY9vV5zQkgKhR6jMTSjWDrLgb+a2cz40QAAR7j7KoD4/+Hx8GLLWwvrIakWlkf1uq9aWB7VazvKecxFpZX1iIzAnenuK83scOBRM1vYTtliy1sL6yGpFpZH9bqvWlge1Ws7QmgplPMYjaC5+8r4fxPwJ6Im9pq4mUn8vykuXmx5q3495Kn65VG9FlT1y6N6bV8ISaGcx2gEy8wOMrMBra+Bc4CXafvoj4uAB+LXk4DPx1c1nA5sjpurU4BzzOzQ+MqHc+Jh1Ur1qnoNjuq1DFmf6U6c4X+F6KqG72QdTwdjfwvRFRizgXmt8RM9OvxxYHH8f1A83Ih+tGgJMBcYlZjWF4H6+O8LWS+b6lX1qnrtOfXa+qfHXIiISE4I3UciIhIIJQUREclRUhARkRwlBRERyVFSEBGRHCUFERHJUVIQEZGc/w+GWtysm9sdywAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### Do this as a function of dimensions retained\n\nUn-normalized"
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "def get_components(data, normalize_variance=False):\n    X, M, s = WaveletFeatures.normalize_datamatrix(Data, normalize_variance=normalize_variance)\n    U, SDiag, VT = WaveletFeatures.get_svd(X)\n    vec = VT.T\n    Z = np.dot(X, vec[:, :7848])\n    return Z, vec, M, s, SDiag*SDiag",
   "execution_count": 98,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "Z, vec, M, s, vals = get_components(Data, normalize_variance=False)",
   "execution_count": 100,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "isinstance(Z, np.ndarray)",
   "execution_count": 68,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 68,
     "data": {
      "text/plain": "True"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "def delta_Data_len(Data, wav, vec, M, s):\n    if not isinstance(wav, np.ndarray):\n        Z = wav.to_pandas().set_index('Object').values\n        \n    else:\n        Z = wav\n    D = WaveletFeatures.reconstruct_datamatrix_lossy(Z, vec, M, s)\n    Delta = D - Data\n    return np.sqrt(np.sum(Delta**2, axis=1))",
   "execution_count": 101,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "ncomp = 1200\nl = delta_Data_len(Data, Z[:, :ncomp], vec[:, :ncomp], M, s)\nfig, ax = plt.subplots()\nax.plot(l)\nprint(np.median(l), l.mean(), l.min(), l.max())",
   "execution_count": 102,
   "outputs": [
    {
     "output_type": "stream",
     "text": "1.59764903133 8.34934252527 0.0031154524072 167.550477798\n",
     "name": "stdout"
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcFPWZ+PHPw+UBHiAjEhUBxSMnGtZojMYrnlmNazaruzGa40XMrdnfZjHXGk3UaDTGqCiKmnjgfXIoyCmgwMAAM9zDMDAXMwPDXMw98/390dXQDH3W0VVd87xfr3lNT3V11TPd1U9961vfQ4wxKKWUCq9+fgeglFLKW5rolVIq5DTRK6VUyGmiV0qpkNNEr5RSIaeJXimlQk4TvVJKhVzKRC8iz4hIjYgUxSx7RURWWz+lIrLaWj5aRFpjnnvCy+CVUkqlNiCNdZ4DHgX+GV1gjPmP6GMReRBoiFl/qzFmvFsBKqWUciZlojfGLBKR0fGeExEBvgVc7CSI4cOHm9Gj4+5CKaVUAitXrtxljMlLtV46JfpkzgeqjTFbYpaNEZECoBH4rTHmo3gvFJGJwESAUaNGkZ+f7zAUpZTqW0RkezrrOb0ZeyMwLebvKmCUMeZM4JfASyJyZLwXGmOmGGMmGGMm5OWlPCEppZSyyXaiF5EBwL8Br0SXGWPajTG7rccrga3AqU6DVEopZZ+TEv2lwEZjTHl0gYjkiUh/6/FYYBxQ4ixEpZRSTqTTvHIa8DFwmoiUi8j3radu4MBqG4ALgLUisgZ4HbjVGFPnZsBKKaUyk06rmxsTLL8lzrI3gDech6WUUsot2jNWKaVCThO9UkqFnCb6EFlbXk9heUPqFZVSfYrTDlMqQK55dAkApfdd7XMkSqkg0RK9UkqFnCZ6pZQKOU30SikVcprolVIq5DTRK6VUyGmiV0qpkNNEr5RSIaeJXimlQk4TvVJKhZwmeqWUCjlN9EopFXKa6JVSKuQ00eewlo4uHpqzmc7uHr9DUUoFmCb6HPa3uVt4ZO4WXssvT72yUqrP0kSfw9o6ugHo6Or2ORKlVJBpoldKqZBLmehF5BkRqRGRophld4pIhYistn6uinnuDhEpFpFNInK5V4ErpZRKTzol+ueAK+Is/6sxZrz1MxNARD4N3AB8xnrN4yLS361glbLrP578mHtnbvA7DKV8kTLRG2MWAXVpbu9a4GVjTLsxZhtQDJztID6lXLFsWx1PLirxOwylfOGkjv6nIrLWqtoZai07HiiLWafcWqaUUsondhP9ZOBkYDxQBTxoLZc465p4GxCRiSKSLyL5tbW1NsNQSimViq1Eb4ypNsZ0G2N6gKfYXz1TDpwYs+oJQGWCbUwxxkwwxkzIy8uzE4ZSSqk02Er0IjIy5s/rgGiLnHeBG0TkEBEZA4wDljsLUSmllBMDUq0gItOAC4HhIlIO/B9woYiMJ1ItUwr8EMAYs05EXgXWA13AT4wx2ptHKaV8lDLRG2NujLN4apL1/wT8yUlQSikVZL97u4gxwwfzva+M8TuUtKRM9EoppQ70/CfbAXIm0esQCEopFXKa6JUjnd09PDR7E83tXX6HopRKQBO9cuTtggoemVfMg7M3+R2KUioBTfQqqdGTZvB/7xQlfL7DmvSkrVMnP1EqqDTR55BH5m5hZmHVQcvjdj120T8+3u7xHpRSXtJWNznkoTmbASi972qfI1FK5RIt0YdAvAGGlFIqShN9CKzYvsfvEJRSAaaJPgRmrD243l4ppaI00SvlsZ4ew63Pr2Tl9nTn71HKXZrolfLYruZ23l+3k1tfWOV3KKqP0kSvlFIhp4leKaVCThO9UkqFnCZ6pTzmdc9lpVLRRK9UlmjHNuUXTfTKJVpuVSqoNNErR0TLqWnTU6HyiyZ65YjR9JWSngqV30KX6I0xvLumkrbObr9D6WM0nSWip0Llt5SJXkSeEZEaESmKWfaAiGwUkbUi8paIHG0tHy0irSKy2vp5wsvg41lSvJufTyvgz+9vzPauVcB8UrKb0ZNmUNPY5ncogJ4KlX/SKdE/B1zRa9kc4LPGmM8Dm4E7Yp7baowZb/3c6k6Y6Wto7QSgOiBfbuWf55aUArBSR/dUfVzKRG+MWQTU9Vo22xgTnQ36E+AED2JTSinlAjfq6L8HzIr5e4yIFIjIQhE5P9GLRGSiiOSLSH5tba0LYSillIrHUaIXkd8AXcCL1qIqYJQx5kzgl8BLInJkvNcaY6YYYyYYYybk5eU5CUMlUFzTzJ69HX6H0Sd8++lljLljRtznTEjvxnb3GEZPmsGzS7b5HYpKwXaiF5Gbga8D/2VM5FA2xrQbY3Zbj1cCW4FT3QhUZe7ShxZyxd8W+R2G77KRZxcX7wptQk+ks7sHgPtmacOHoLOV6EXkCuB/gWuMMS0xy/NEpL/1eCwwDihxI1BlT3Vju98h+EYC0swlKHGovmtAqhVEZBpwITBcRMqB/yPSyuYQYI5EjuJPrBY2FwB3iUgX0A3caozRaXWUUspHKRO9MebGOIunJlj3DeANp0Ep5Ya+VpWiVCKh6xmrVG9ac6L6Ok30KvT8LtjrlYXymyZ65ZLgZbOg3QQNWjy5rKG1k46uHr/DyBma6JUjOkyx8sMX/jCbW19Y6XcYOUMTvXJdfUsHRusrlMfmbayx/doeq7NXXxn8MHSJXsdHz67e73dZXQvj75rD1MXaW7I3PfcFR7f1YTy1qG908wlFon+roJxnNLH4LFKFs6Mu0n/OSWkr19U0tnHnu+vo6tY6ZBUMoUj0t7+yhrumrwe0zlj579dvFfLc0lIWbTlwsD69Gav8EopEr1SQdPdoHY0KlpxP9HrTTymlksv5RP9WQYXfIaiAKqxo8DsEpQIh5xN9UUWj3yGogCrf0wr439pFW4Ipv+V8ou/LRO/u5RRtKKD8ook+h+n9idzSu2RfVNHAj15Yqc0wA6aktpkP1u30OwxXaaJXoef3hU+ikvzPpxUwq2gn2+ta4j6v/HHxgwv54fOZDa/Q0tHFi8u2B7bwlXI8+qDz+0usVFBNXrCVgf2FH5w/1u9QQu/u6RuYtnwHJw49nAtODd4c2Dmf6JVKxe9Cll83Y6PjuGii997u5siUnS0d3T5HEp9W3WRZT49hdVm932F4IJiXrEGiN2OVX0KX6BOVnlaX1XPNo4tp6/T3jDt54Va+8dgS8kvDMZWuJq/E/L6SUCoqdIk+qncC+sN761hb3sC6Sn/b3a+viuy/qqHN1zi81peTnDZ77bsKduzxO4S40kr0IvKMiNSISFHMsmEiMkdEtli/h1rLRUQeEZFiEVkrImd5FTzofKAqeLxoebGruZ2/ztlMj46jE2hPLioJZHPZdEv0zwFX9Fo2CZhrjBkHzLX+BrgSGGf9TAQmOw8zc9obMTsSvc/ZKNSW1bVQFuCmiZm+B13dPWytbY773K9eX8vf5m5heY5V+fX0GKYs2kpze5ffoWRNEDNPWoneGLMI6H2EXQv8w3r8D+AbMcv/aSI+AY4WkZFuBJsOrTP2S/bf9/Pvn8/598/P+n698uf3N3LJgwvjnrz2Woky16rE5myo5p6ZG/nTjA1+h9KnOamjH2GMqQKwfh9rLT8eKItZr9xa5qtc+4Kovmd5aaR+d5fVVC9W9Pjtl2PlmGjjh9gS/dLiXZx33zxaA9oUMYy8uBkb71A8KM2KyEQRyReR/Nra2jgvUSp9RRUNPLlwq99heKbHyvT9ci3Tx/HHGRuoqG9NWE2l3Ock0VdHq2Ss39G548qBE2PWOwGo7P1iY8wUY8wEY8yEvDzve5JpQ4hw+/rfF3PvrNyc6Dmdq83oHKdO8nxXdw8dXcG7UeiHFda9ji4Pbm5v372Xbbv2ur5dJ5wk+neBm63HNwPvxCz/jtX65hygIVrFkysq61v1CxEift+Yd6PasGdf1Y39TH/No0s49beznAeT4zq6evjPp5bt+zteVZkTlz60iIv+ssDVbTqVbvPKacDHwGkiUi4i3wfuA74mIluAr1l/A8wESoBi4Cngx65HbUO6X7a2zm6+fN88/uf1Nd4GpPqcRDk6ndxt9pXo7Sf6aB+Ovq6nVzJoagt/i6C0xroxxtyY4KlL4qxrgJ84CSoTmR73qdbvsNrAzttQk3xFlXXtXd20tHczdPAgv0PJuh4XEr3qu0LbMzYRu5fR762pZPrag241qCz67rMrOPPuORm/zssmt+1d3YyeNCPuc72PNSdVONE+OHby/Nryel7NL0u9ogqtPjN6pdOC0M+mFQDw9c9/ytmGbH7Z443R09e62i/dutvvEA7S0NJ50LLen4sbH5OTqptrHl2S1nqtHd18sG4n3zjT99bQymV9pkRvpzT1an4Ze/Z2uB+MDb94ueCgZUGd5CBognYz9qlFJZwWwJuid01fx22vrGZZibcn1D5WPgmE0CX6VF/qdA+ypvYufvX6Wn4eJ8E6YvMgn78x6H0NDnzf9Rx0sOix96eZG2gPYKuundZAe3s73L85qYUSf+V8ok9UfeFWvWxtk7tNrzLx67cKWbxll2/7T0fv91kLa5rUVPDkfKIPs5eW7eDbU5elXlElZTfvfrx1N9c9voROm6MRRssgQZ11KJv62v2koAltop9RmFN9tEIjiGXZldvtjRH+P6+voWBH/b4qDbu+99yKlOtU1reyJpQzjx1ML3iyL7SJPgjG3jGDH72Q/mzyj80vZtPOJg8jyp54BbjVZfU8/3FptkNJu0T+82kFcT8vu4kp+rodaQylHJ3fNZX/fm0NW6rDcYyo7AldondaN+9maaPHwKyinWmt29ndwwMfbOK6x9NrChcUmbRo+cZjS/jdO+s8jMaZd9dUHvB52a1t2FwdGaxr2vL02q5ncsxtqGrkZ9MK2NvexRf+MJtFm4N+k14FQegSfa7L3TF2wlsHm/RkFuffrqhvBWDRlhRJ2MFbVlzTTENrJ3+Zvcn+Rnyi1fXZ12cSfbqFJicH4brKBqoaWu1vIIaO1e0/ryexKakN1giHKrz6TKLPhqsfWcy5985zvJ3O7h7O+P37LkSk3JC0asWFqr7ehYuOrh6+dM+HzF63vxpJb2AqJ3I+0adb5grq1WJPjzloeIN4Nw/97t3pl8cXFHOTwyamdt85x1UMNndc3dROdWM7d01f7zAAe/SkEj45n+iDaOX29CdwvvO9dZz+u/fp9mAChDC4//1NfORzpzG3Pplk542gFkSi3i+qojRgk2nkkhKfZ9PSRO+B6yd/zPbdCb4UMVlj8ZZd/PPj7ZHFKbKJ2/XFxhjunr6ejTvDP0a53Xcu6MnXK/GuZG59YRWXPrTQ9jbj9RbOxpXDNycvZcoif6eYnL1uJxc/uJCZPvbt6TOJPtvl5cbW5OOF/GxaQdq9Xnc2tO0bJz9dNY3JO/nsau5g6uJt/NdT2vM2FWMMr+aXMXrSjLijiHq335jHAai6c2vavWy2usnfvod7Zvo7xeRGq2/MBh8nfukziT4q02NsT0sH9S3+jmC5dGvmVRdTl2w7aFl1YxvrKhsA9rUO8vpL19DayehJM3hzVXnKdRvbDh7yNygetJox1vceljiD9y/tll/pb5KG1uC+Z7HiDYHwr48upiXJAGpfvncuv3x1tZdh9Rm5n+h7HT+pSj6ZlkmqG9sZf1fmk10klOJb7GW57bz75nH1I4uB9McotytaGi2zeoU+/dHBJ55Yq8vq+fyds5nlweWtF+/prMIqmtuzNwVdoqq77btT97oNgrq98QcH3N2cuBBV2dDGm6sqXNn/1MXbUl7lhlnuJ/o0bN+9l6r69Nq3p6o3PP/+eXyUqiNML5f9NXXdZjYuZ72Y8X6/yLbt/huFFZErjcXF6V293J2FFimJBuLatLOJH724iklvrM1se24EZfG/IiczS4r3j3HvR6ueu6ev58cvrsr+jgOiTyT6rz6wgJokww2vr2xkxtr0SpJlda38cfqGjPYf7RKfjkxa37gzIqCzbXjdqSiRqYuTXyG4qfcnEi3JV6QqPNh8a+J9rEGoo891fWES8ERCl+hfzd9fF3zpQwtp7zrw5lm8795Vj3zET17y92x/+u/2d5Da0uvEkKgEFNuS4YEPDrzh5FcCDhO338FEqdpOCTdon64xhmnLd7A3i9VZKn22E72InCYiq2N+GkXkNhG5U0QqYpZf5WbAqcQO8lRc07yvjjiX2Cl5PDbf3yZkUblQ7ly5vY66DKaI1A5EqX1SUscdbxby+wAPWteXr4psJ3pjzCZjzHhjzHjgi0AL8Jb19F+jzxljZroRqF29xyLv/VF39Wq2GMQBl7yMyattB/F9jLp+8sd868mPU68YoP8h6Ceb1s5I4STRTddYQT42wsqtqptLgK3GmO0ubc81qerHq2xMKtHj87du+tpKCssbfI0h1xXXZNJT0e6A9Af+mbRnbICyX9BPKrluxtoqbnN7LuoU3Er0NwDTYv7+qYisFZFnRGRovBeIyEQRyReR/Nra7I2p7cbXaUtGScK53l+8n75UwL8+ujjpa5zmjQdnb+L37xSlXM/Ly2Ev5l7NZJuxb2H0Zc8tLbW/70T7CU6Ot61ub+r2/EvTbFHllXQ++mzcY/jJS6t4e3Wl5/uJ5TjRi8gg4BrgNWvRZOBkYDxQBTwY73XGmCnGmAnGmAl5eXlOw0goaKWTht4dbjzi9P/++7zifcMzpCd5trITzmsrU3eyytR1jy/N+DWx7+UTC7fanlow088k9qRkd97aZIprmhk9aQYLe01eYufE8/9eWxN3eeyVSrKWb0Fx/eTMj49c4EaJ/kpglTGmGsAYU22M6TbG9ABPAWe7sI+EYluXOO2ePuaOGayr8LabcrrtxP3y3ppKRk+a4XcYQGRybrettjkva2zy67YSsFeFiHgtpjq73d9ZdPC9GWu9K116cVXmpY0eTuXZ0NrJLc8u92z7ybiR6G8kptpGREbGPHcdkPr63yVvFxzci6531UKygYWMgReXBe42g+uX9u8nmd7wvTXxv/RvFZSz2MEokn7WTuRYrlEh9fKKMhZs8mfqR0eJXkQOB74GvBmz+H4RKRSRtcBFwO1O9pFZPKnXeSpFV/wg5gRb7ayTvBe3xkyAnW4Cvv2VNWkPwuYWv6uuU90gTXW82b1/Ed1uEI9FN3jRx6OtszuHp+H03gAnLzbGtADH9Fp2k6OIXKalufBbk6I6xukVkR5C7rJzAnz6oxJOOXYIF552bNznT//d+3zqqENZesclTsPbp6qhlZFHHeba9vw8kHK+Z6zTL/FLy3e4E0iOamrrYpMH9ZLZPMEW7NiTeiUbnJY7tXeye/44YwO3PLsi6TqVNm+SJ+LGtKBBkfOJPpVPShLf0Cura2HygmD0KPVLa2c3lz+8aN/fTk+crqY2lzbm9KST7PWZJPN0S7Lx9pdrNzWDyO93MNM5JdwU+kSf7C56aaJZoEIgbGXJZEMWePUFdnwTPM3X+53DN1c3HzARuVfCOARBdY4MfRyqRB87oFk6bpp6cFOnsCXIdO1uTq+Nc+KSpftf4tjS8ll3uzgnQIYMJqvJOP7old5ZXVbPxOdXMt+nFiG57Ev3zN03mU+QhSrR9x7Xxo7enUf6il+/VZjWer0vP9OturCTqJxUdWQqXpv9RP/bz6cV7Ntv0hh7PZW8Csgd762pdFzK9LJ/QGF5A0Ue91WBSCucif/MP2CZ3eqv3c3trNqxh3dWVzD+rtkH9TMpqQ1+zYCjVje57P2i+O3pwzBm9YJNNRm/pq0zksC96JiTyNryek4cejhDBw/K2j7jmb1uJxOfX5nweWP8GaYg07zU1tnNz6YVMDZvMPP++0JPYnIq1dAdblm0uZbZ66td2dZ1jy9lR10Lhw3sT2sW5wx2U6hK9JmYvLDE7xAcS9TOO1XrhER2NrQxb2PmJ4l0xIv0mkeXJO1ynq1WK5Vpzj6WLXb/6+hge1X1zkr0bp7UvB6s7a731jPfRsEmEztSDHWeC3ce+myiVwdLdUCnI9ODvmRX8C97vZZOLnSzOqW4ponRk2awrKTO0XaiVUTtPnZUembJNr5rs2CTqVy+mayJXnkiWfLaleaNX79HdUy1f696xka1d3XzzOJtdPeYA7bltJS81LofMd3hROzpTJGpzUKDIecTfV9tJQPJv0R2ckG2vpQT/vghJbXZHeo5mVSJ06129OmKhrOnpZO7pq/nlRVlB+wnKMkzKHHksk9Kdqeee9gFOZ/oVXb1/m7bLbUWVbrX8sLPdJOsE0zvk4DdOH/9VqGr1QZ+5We/qz6CeFq6YconXPTAAs/3o4k+S5xMaOD0AG1oTW8MfGf7OTCpdXUbRk+awfMZjWmfbItOOPvPkiWo8+6bl/TvA7YTsEyjg4BlxqvPLxs9Zvtkom/tcLeJVDqTQjz84WZX95mJ377t3kjR6Z409nZETmyv5Je5tu9EvKpCiFbp1O3tiPslL9hx8GBqdi7Df/X6Wtq69h+TXp8Q/Kpy0bF//JPzid5OXfQZv38/5YiH6VpX2cC438xiToo2u9F26tkS+6VqTDM5p+OHSdqbB8GmnU1p3+xN101Tl7OnJfEQDOlKdKyu3L6HmYX7hyAo3xOs5p7LtzlrnRMYDs9viT4/Y4ztCW2yJecTvZ+qGlp53BoUzav252CvCsOr+tDyPek1wTy4Lj87Ln94EZc8uLDXUnslydhXZbMj2Z9mbjhomZ/VPk9/5F2fk+27nTfpzZZkn4GdaRLvifM5e0UTvQOXPbSIGWudNVELGzcvzu22Iky3einbah3Mmepmnu/qSb41r04qzXHuU6UzmY3dOXqDbsqi7HXa1ETvQFMWZoyH7JWGvSw19s7Z7VntSp78H/O7vX62zd3gbU/SWLFNV+2egM+5d65b4TgSsHvpGdFE75JpaUxgks2bYNm/8ZXZ//Y/r6/1KI7g8Ls5YSKp4orm5gWbarI+DK8xhqqGYN2jCANN9FliMFm7Auht4ebag0bcc4uXJ5RsnKwaWjrJL40/6mmQSvqxZQSvx4+JuuXZFVzj4iBk6RR0nl1Syrn3zvNk1rNs82rmMzs00ecAO19rOyWxoJZAvXTTM8t4d01l3OeClOhjpUqYbn6O1Y3tjt6HTK9io8MzbPdhUqAdLt8Yvu7xxAP2ZZvjYYpFpBRoArqBLmPMBBEZBrwCjAZKgW8ZY4JzevOBILbrwO287LWVmU3Ckr7sZb90E4yTGrGiiuBPGtHbxp1NzHQwTk3QOm7F40WIqbZ5wQPznW0gwNwq0V9kjBlvjJlg/T0JmGuMGQfMtf5WLvuoeJcPe83e0e53QnJadeRu1dOBb8aTSVpstHX28HZBRZpbyp5sVTllm9/HaTq8qrq5FviH9fgfwDc82k/O9LZranO/yV8QZ7bJgWM+J2WaTG57ZbU3gQSAnYl1XJGow1QOHPVuJHoDzBaRlSIy0Vo2whhTBWD9PtaF/eS0t1dX9qmxRXoX3ux8FdKuunH5i9be1c3W2mYKA1St4+pkIO5tynXp/J92J9ax0+qtLHaOhuDn84TcmErwPGNMpYgcC8wRkY3pvMg6KUwEGDVqlAthBF86Y+L4ysMDeUOV9/OEuuWONwp5M0nVhx/crB5wa1NubWfTzibWVTbwb2ed4NIWk8vkvTz//v319tkYfMwrjkv0xphK63cN8BZwNlAtIiMBrN8HXWsZY6YYYyYYYybk5eU5DSMnxOvaHlZdDoYM2LZrL8YYR6XYv324xfZr/bn3kZyr5+AUGzMmu61eLn94Eb98dU3W9tcXOUr0IjJYRI6IPgYuA4qAd4GbrdVuBt5xsp/kMXi1ZfcFebgEu+9jomqTTdX22kGv3L6Hi/6ygBeWpe6AlsxfY0YLzbQ0HK/Xbi4cZ27er/pqzBjpsT1aO7p6ko7+muyGayZVJ0EfPC/XOK26GQG8ZX24A4CXjDHvi8gK4FUR+T6wA/h3h/vp09wcfTIRY9wdj96u6Byya8rqGdjfn+za2HZwx7Z+InRneMYwGNd6lrrZqzrTexrbYub1vX7yUgorGii97+r42+4Vp9stbb715Meubi8bHprj3xDlUY4SvTGmBPhCnOW7gUucbFvtd9f09X6H4BN3koQbucbWCKIGulMMIJb2ttJaJ719pTpn9H6/Yq8U/L5BnWrI5JLa5qQjjZbu2ptySPFMJXs/27I6plNibtyMVX3Irubk47LbSYi3vVzAiCMPtRdQgOVCdY/b1pT7eyK42BqiespNX4z7/DefWHrQMfztp5fxu6+fYXufyRL9Xz7YZHu7btJEb1PdXucTUYTR1MXbMn7N26vjD0GQrupGdycaicfvlnVutrrx8gTU+3sRW5WTTjWOl+dGA9S3HFwNWlHfyq0vrPJkn40e9J+xQxO9TWfdPcezbb+9OlhN+zIx16UJWDJJbKmqR/zquZjNznyPzS9mjxY+9mkNSJVJUGiiD6DXPRunJjlnCdH9bNoXqz4SSXYz1hjDAxlUEWTrxNfa0c1THs5OlcwvXs5ez+DYm9VBpaNXqn2a2zO/zIwttb6XYBRIOzTJHyhZbn7+k+2ubctND3+4maKK4HSU8+oE9+j8Ym827CJN9AqIJFan9ZQ/m1bgOA47w0S43YoiniCfd1YkGE8/kUybau7MsInoXmvehd4zSmVz4p1Esl2A2J2i8UK25HyiD/IXMJcE4Dt4kHQ/24r65DMSheHqIJf+hTXl9QC8vKLM50j859Y9K6dyPtEr/3R19/DRllq/w0gp0UlsyqKt6W/DtVjsbcnN83Cq1i/phnjAgF8ueP6T7exweZsqQm/GKtsembuFWUU7Pdm2MfaqcTJxz8yN9EuzqGyn41O8Dkwrt9ucf8fNQc1cunxLdqN1vo2hhH/3dpGTcFQSWqJXtm1zeeq13rybJSs7egz88+MDb5S+47DPQDybM5xfddWOetdj6O32EI+Hn4tyP9GHoQI2APzu2t5bGD7W7h7DEwsPrB7KtIVMVLIyuN0B5BLx6r33e4Ypg8mZiYrclvuJPqD+9dHFfoeQETs9ffvmV8YfQWix0ltfTZq5SBO9ckWqMXD8FMAcGWrrK4PTdl5FaKJXttS3dFBUub+6Jxtt2fuyXDpX/XGXTh9EAAANVklEQVRG/Al2gnhV0ldoole2fOvJjz2dnDxZThg9aUZG20pWNZwrqSfXcmSiwcP8ZAx9tr5RE72yZXN1s98huCJXEqjbE6C7Id2mqVHTlgegA1Xw3sas0ESvlDpAOuPTbKluorLBndmzsqV8Tys9uXJmd5kmehVIbrbEC8N3u63T285jsdKZ+u6ZJaXeB+KBLpdm/Mo1muhV6K3aYbM3qkooDP0c+pKcT/R6vIWTm6XwLTXhuJ+glF05n+iVUtmnBazcYjvRi8iJIjJfRDaIyDoR+YW1/E4RqRCR1dbPVe6Fq/oKt0dGVO7Sqpvc4mT0yi7gv40xq0TkCGCliEQnUv2rMeYvzsNTfVWzNXmFCiYd/iC32E70xpgqoMp63CQiG4Dj3QosXVqyCCf9XJVyjyt19CIyGjgTWGYt+qmIrBWRZ0RkaILXTBSRfBHJr60N/uQVSqn99EScWxwnehEZArwB3GaMaQQmAycD44mU+B+M9zpjzBRjzARjzIS8vDzb+2/t7Lb9WhVcYWj7rlRQOEr0IjKQSJJ/0RjzJoAxptoY022M6QGeAs52HmZiTy5MPMuNyl2a54NNC/S5xUmrGwGmAhuMMQ/FLB8Zs9p1gM4PFjJ/fn+j3yEon/k9iYjKjJNWN+cBNwGFIhKdN+zXwI0iMp5IoawU+KGjCFXgTF6Q/qTadmkaUco9TlrdLCb+93Gm/XCUUkq5TXvGqkDSOnql3KOJXimVseeWlvodgsqAJnoVSBuqdN5RpdyiiV4ppUJOE71SSoWcJnqllAo5TfRKKRVymuiVUirkNNErpVTIaaJXSqmQ00SvlFIhl9OJ3uig5UoplVKOJ3q/I1BKqeDL7UTvdwBKKZUDcjvRa5FeKaVSyu1E73cASimVA3I70WumV0qplHI60Xf3aKZXSqlUcjrRf7Bup98hKKVU4OV0ou/o6vE7BKWUCjzPEr2IXCEim0SkWEQmebEPo7djlVIqJU8SvYj0Bx4DrgQ+DdwoIp92ez/le1rd3qRSSoWOVyX6s4FiY0yJMaYDeBm41u2dNLZ2ur1JpZQKHa8S/fFAWczf5dYyV13x2ZFub1IppbLq+18Z4/k+Bni0XYmz7IAKdRGZCEwEGDVqlK2dnHvyMfz4wpMprGigtqmdU44dwqadTWypaT5o3X8ZPZS8Iw5hZuFOvnjSUAorGjglbwgAhw7sR21zO8bsrw464tABjD5mMEWVDQzs14+O7gNv/PYT6DEw5JABNLd3AXDqiCHUNrVz2MD+HHvkoZycN4SSXc0U7Kjn+KMPo6K+le+eN5pnl5Qy4shDqG5sB+DKzx5HUWUDFXta6TEwatjhHH/0YVx35vEMGzyI5vYuCisamLp4Gxeffiz9BD7cUAPATeecREV9K2OHD2ZW0U4q6iPxP/wf47ntldWcnDeYrbV7ATjuyEO57DMj6Ow2vJpfxqVnHEttUzu793awt72LW796MjMKq7j90lMZNngQJbv2snhLLRecmsfa8gbK6lq49IwR1Ld28tSiEn5w/hiqGto4a9RQdjW3896aSk4cdjgjjzqUivpWTjvuCArLGyjd3cL/XnEan5TUUdvUzhuryve9jxMvGMuURSUA3P/Nz/Pu6krWlNXT1N7F9WedQGd3D1d9biSziqo47bgj2NXUwZwNOymra+XwQf1p6egm74hDOCVvCEWVDTS1dTF8yCAaWjsZNexwbjlvDDWNbTy5sISvnpbH+eOG00+EmsY2FmyuZcghAxibN5gRRxxKZUMbK0rrGNS/H7+64jSKa5pZvq2Owwb157CB/dnb0c0tXz6Jtwoq+Mopw/loyy7+ZfQwVpfVU1nfSnVTO+2d3dS3dDLqmMP50phhLNpcS0V9G7d/bRyTF2xl8KABHHvkIWza2UT/fkJXj+H6s06gdNdeTh0xhIKyeoYcMoDrzzqBwooGRKC5rYvtdS18UrKb804ezqeOPowNVY30GMO144+nor6FFaV7WL6tjm+fM4qW9m6Wl9bx4L9/gW5j+M+nlnFy3mB+c/UZvFVQyXtrKrn0jBHcdO5JrNhWR0NrJ9WNbZx0zOF88aSh/Or1tbR2dvObq85geWkdP71oHJuqG3n6o21cO/5TdPeACAzs348nF27lys8ex4cbavjKKcN5e3UF7V09XP6ZERx92CBOGn44g/r348jDBnLhqXn89u0iZq+v5vxxw7n6cyOZ9GYhXzllOLVN7Xzm+COpbWpn4gVjqW/pZGZhFbOKdvKtCSdw5qihNLR2ct+sjZww9DDK97TypTHDOOXYIVQ1tDFq2OEU7NjD7V87lZOOGcw1f1/MGZ86ks8dfxQAA/oL762u5DtfHk1LRzdPLNzKLy4Zx8D+wj0zN3Ldmccz+pjBrKts4NQRR3DUYQOZUVgFwNi8wRwyoD8Xn34shRUNPDJ3C9859yQq9rSye28HV3z2OE7JG8IP/pnP6ccdwdljhrGkeBefO/4oLjr9WN5bU0VtUxunH3ckm6qbOO+UYzjvlOHM31jD+ePy+PLJx9jKf5kQL4YREJFzgTuNMZdbf98BYIy5N976EyZMMPn5+a7HoZRSYSYiK40xE1Kt51XVzQpgnIiMEZFBwA3Aux7tSymlVBKeVN0YY7pE5KfAB0B/4BljzDov9qWUUio5r+roMcbMBGZ6tX2llFLpyemesUoppVLTRK+UUiGniV4ppUJOE71SSoWcJnqllAo5TzpMZRyESC2w3cEmhgO7XArHTUGNCzQ2u4IaW1DjAo3NrnRiO8kYk5dqQ4FI9E6JSH46vcOyLahxgcZmV1BjC2pcoLHZ5WZsWnWjlFIhp4leKaVCLiyJforfASQQ1LhAY7MrqLEFNS7Q2OxyLbZQ1NErpZRKLCwleqWUUgnkdKLPxgTkcfb5jIjUiEhRzLJhIjJHRLZYv4day0VEHrHiWysiZ8W85mZr/S0icrMLcZ0oIvNFZIOIrBORXwQotkNFZLmIrLFi+4O1fIyILLP284o1pDUicoj1d7H1/OiYbd1hLd8kIpc7jS1mu/1FpEBEpgcpNhEpFZFCEVktIvnWsiB8pkeLyOsistE65s4NSFynWe9V9KdRRG4LQmzWNm+3vgNFIjLN+m54f6wZY3Lyh8jwx1uBscAgYA3w6Szs9wLgLKAoZtn9wCTr8STgz9bjq4BZRGbcOgdYZi0fBpRYv4daj4c6jGskcJb1+AhgM5GJ2YMQmwBDrMcDgWXWPl8FbrCWPwH8yHr8Y+AJ6/ENwCvW409bn/MhwBjr8+/v0uf6S+AlYLr1dyBiA0qB4b2WBeEz/QfwA+vxIODoIMTVK8b+wE7gpCDERmQ61W3AYTHH2C3ZONZceUP9+AHOBT6I+fsO4I4s7Xs0Byb6TcBI6/FIYJP1+Engxt7rATcCT8YsP2A9l2J8B/ha0GIDDgdWAV8i0hlkQO/Pk8g8BudajwdY60nvzzh2PYcxnQDMBS4Gplv7CkpspRyc6H39TIEjiSQsCVJcceK8DFgSlNjYP5f2MOvYmQ5cno1jLZerbrIyAXmaRhhjqgCs38dayxPF6Gns1iXemURKzoGIzaoaWQ3UAHOIlELqjTFdcfazLwbr+QbgGK9iAx4GfgVEJwY+JkCxGWC2iKyUyDzL4P9nOhaoBZ61qrueFpHBAYirtxuAadZj32MzxlQAfwF2AFVEjp2VZOFYy+VEn3IC8gBIFKNnsYvIEOAN4DZjTGNQYjPGdBtjxhMpPZ8NnJFkP1mLTUS+DtQYY1bGLg5CbJbzjDFnAVcCPxGRC5Ksm63YBhCpvpxsjDkT2EukOsTvuPbvMFLPfQ3wWqpVE8TgxbE2FLiWSHXLp4DBRD7XRPtxLbZcTvTlwIkxf58AVPoUS7WIjASwftdYyxPF6EnsIjKQSJJ/0RjzZpBiizLG1AMLiNSHHi0i0VnOYvezLwbr+aOAOo9iOw+4RkRKgZeJVN88HJDYMMZUWr9rgLeInCT9/kzLgXJjzDLr79eJJH6/44p1JbDKGFNt/R2E2C4Fthljao0xncCbwJfJwrGWy4k+SBOQvwtE78rfTKR+PLr8O9ad/XOABuuy8QPgMhEZap3lL7OW2SYiAkwFNhhjHgpYbHkicrT1+DAiB/wGYD7wzQSxRWP+JjDPRCoj3wVusFojjAHGAcudxGaMucMYc4IxZjSRY2ieMea/ghCbiAwWkSOij4l8FkX4/JkaY3YCZSJymrXoEmC933H1ciP7q22iMfgd2w7gHBE53Pq+Rt837481t258+PFD5I75ZiL1vb/J0j6nEalf6yRyZv0+kXqzucAW6/cwa10BHrPiKwQmxGzne0Cx9fNdF+L6CpHLt7XAauvnqoDE9nmgwIqtCPi9tXysdYAWE7nEPsRafqj1d7H1/NiYbf3GinkTcKXLn+2F7G9143tsVgxrrJ910WM8IJ/peCDf+kzfJtIyxfe4rG0eDuwGjopZFpTY/gBstL4HzxNpOeP5saY9Y5VSKuRyuepGKaVUGjTRK6VUyGmiV0qpkNNEr5RSIaeJXimlQk4TvVJKhZwmeqWUCjlN9EopFXL/H6O8Ij3zS82tAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "ncomp = 100\nl = delta_Data_len(Data, Z[:, :ncomp], vec[:, :ncomp], M, s)\nfig, ax = plt.subplots()\nax.plot(l)\nprint(np.median(l), l.mean(), l.min(), l.max())",
   "execution_count": 93,
   "outputs": [
    {
     "output_type": "stream",
     "text": "1353.73984793 6573.77512973 143.372970994 159126.749138\n",
     "name": "stdout"
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAD8CAYAAACLrvgBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXmUVdWV/z9bEJwFBY0RDBiJBk2ihqiZXElMFJO02N2mF3a6JWk7ZDA9pX/donZao9KaxNbEOCQGiKBGpY0GFBSZFFAECpCZgqKYiqGqoKiiqKKGV7V/f9xT8ijeqzff4b39WatW3bvvuffsd4fzPfMRVcUwDMMw0uG4oB0wDMMwooOJhmEYhpE2JhqGYRhG2phoGIZhGGljomEYhmGkjYmGYRiGkTYmGoZhGEbamGgYhmEYaWOiYRiGYaRN76AdyDcDBgzQIUOGBO2GYRhGpFi+fPk+VR2YKlzRicaQIUMoKysL2g3DMIxIISLb0wln1VOGYRhG2phoGIZhGGljomEYhmGkTUrREJFJIlIjImu72f9JRMpFZJ2I/CLOfoeIVLhj18XZRzpbhYiMi7MPFZElIrJZRF4UkT7O3tftV7jjQ/Lxgw3DMIzsSaek8TQwMt4gIl8GRgGfVNWLgYecfTgwGrjYnfOEiPQSkV7A48D1wHDgZhcW4OfAI6o6DDgA3OrstwIHVPUC4BEXzjAMwwiQlKKhqguAum7mHwIPqmqrC1Pj7KOAF1S1VVW3AhXAFe6vQlUrVbUNeAEYJSICfAV4yZ0/Gbgx7lqT3fZLwDUuvGEYhhEQ2bZpfAz4oqs2eltEPuPs5wI748JVOVsy+5lAvarGutmPupY73uDCG4ZhGAGRrWj0BvoDVwH/AUx1pYBEJQHNwk6KY0chImNFpExEympra1P5bhiGEQgVNYd4r3J/0G7kRLaiUQW8rB5LgU5ggLMPjgs3CNjdg30f0E9EenezE3+OO346x1aTAaCqT6nqCFUdMXBgygGNhmEYgfDVh99m9FPvBe1GTmQrGn/Ga4tARD4G9METgOnAaNfzaSgwDFgKLAOGuZ5SffAay6erqgLzgZvcdccA09z2dLePOz7PhTcMwzACIuU0IiLyPPAlYICIVAF3A5OASa4bbhswxiXo60RkKrAeiAG3qWqHu86PgVlAL2CSqq5zUdwOvCAi9wMrgYnOPhF4RkQq8EoYo/Pwew3DMIwcSCkaqnpzkkN/lyT8eGB8AvtMYGYCeyVe76ru9hbgW6n8K0UOtrRzQu9e9OltYzMNw/AXS3UiyCfveZO/m7AkaDd85a3yGoaMm8H2/U1Bu2IYJY2JRkRZui1hn4Ci5ZWVuwBYuaM+YE+MUqKzUxk7pSzyPZ7yiYmGYRhGEhpbYry5vpqxU2y5hS5MNAzDMIy0MdEwDMMw0sZEwzAMw0gbEw3DMEoOVWXO+mpiHZ1BuxI5TDQMwyg55m2s4R+nlPH4/C1BuxI5TDQMwyg59h1qBWBXfXPAnkQPEw3DMELJi8t2UFHTGLQbRjdMNEJMrKOTK/9nDq+u2p06sGEUGbf/aQ3XPrIgaDeMbphohJjGlhjVB1v56bS1qQMbRhHSafNahw4TDcMwDCNtTDQMwzCMtDHRMAzDMNLGRMMwDMNIm5SiISKTRKTGrdLX/dj/ExEVkQFuX0TkURGpEJHVInJ5XNgxIrLZ/Y2Js39aRNa4cx4VEXH2M0Rktgs/W0T65+cnG6XCKyurONjSHrQbhlFUpFPSeBoY2d0oIoOBrwE74szX460LPgwYCzzpwp6Bt0zslXir9N0dJwJPurBd53XFNQ6Yq6rDgLlu3zDSYuPeg/zbi6v4f1NXBe2KEUL+r6wKgJh1z8qYlKKhqgvw1ujuziPAfwLxd30UMEU93gP6icg5wHXAbFWtU9UDwGxgpDt2mqoudmuMTwFujLvWZLc9Oc5uGCk53NYBQHVja8CeGGGkbPsBAJpaYwF7Ej2yatMQkRuAXaraPRt3LrAzbr/K2XqyVyWwA5ytqnsA3P+zsvHVKA7UMoSh488rdzFk3AwamoOpAox1dHLfa+vZf6hwGQPFXrzuZCwaInIScBfw34kOJ7BpFvZMfRorImUiUlZbW5vp6YZhZMGERZUA7KgLZv6mORuqmbhoK/e8ur7gcbmmVoPsShofBYYCq0RkGzAIWCEiH8IrKQyOCzsI2J3CPiiBHaDaVV/h/tckc0hVn1LVEao6YuDAgVn8JCPs2DdrdKdrVvOOzsJPb65W1P2AjEVDVdeo6lmqOkRVh+Al/Jer6l5gOnCL60V1FdDgqpZmAdeKSH/XAH4tMMsdaxSRq1yvqVuAaS6q6UBXL6sxcXbDZ9pinXxm/BzeWLs3aFeMEFLMVTiSsDKktEmny+3zwGLgQhGpEpFbewg+E6gEKoDfAz8CUNU64D5gmfu719kAfghMcOdsAV539geBr4nIZrxeWg9m9tOMfFHT2EJtYyv3vVb4agAjOliCWpr0ThVAVW9OcXxI3LYCtyUJNwmYlMBeBlySwL4fuCaVf4ZhGIZ/2IhwwwiYCQsrmbCwMmg3ssaq+0sLEw3DCJj7Z2zg/hkbgnYjY1J1TvjE3bP4uwlL/HEmS6yKLXNMNAzDKAiNrTEWVewL2o0eKeZG/EJhohEBrPhvhBl7PbPjjbV7+OidMz+YvSAqmGgYhpEVVrGTG7+cVU5Hp7KrPpjBkdliohEBbGCbYRhhwUTDKFl21jXz0KxyG+2bJvPLa/jiL+bRGju6OsXuX2lhomGULD94djmPza+gouZQ0K7kzHf+sJSP3jmzoHHcPW0dO+sOs7ehxTNYEbgkMdEwSpa2mDdn0b5DbQF7kjtvldfSYWtDGD5gomEUN2lUndz8+/d8cKR4MakqLUw0DMPICqucKk1MNIzixurdC0bQJQw/BubZ4L9jMdEwcuKdin0s3ZpoNeBooap0WptAVpSCLNsiTEcw0UhBbWMr+wq4nGTU+faEJfzN7xYXPJ5C9+r8wbPLOb/AvY/CSk1jC0PGzWD6qt2pAycgKKm1eaOCwUQjBZ8ZP4cR988J2g2jwMxaVx20C4GxudrrcvzC0h0ZnWdJdmliopElh9s62BnQ2silSJRqBxpb2nloVjmxjvwvQ9rS3sHnH5zH/PIjqx+/G/JJAYsBG8B4hHRW7pskIjUisjbO9ksR2Sgiq0XkFRHpF3fsDhGpEJFyEbkuzj7S2SpEZFycfaiILBGRzSLyooj0cfa+br/CHR+Srx+dD7779FK++Iv5aYePdXRy6b1v8vKKqgJ6ZYSBX84q57H5FUx7P7vqnp7YVX+YXfWHue/VI6sortnVkPd4MqGY01OrAjuWdEoaTwMju9lmA5eo6ieBTcAdACIyHBgNXOzOeUJEeolIL+Bx4HpgOHCzCwvwc+ARVR0GHAC6lpO9FTigqhcAj7hwoeG9yswaf5taO6hvbuee6esK5JERFlravWk2Yp35L2mEibCW/tpindz32noamtuDdqUoSSkaqroAqOtme1NVY273PWCQ2x4FvKCqraq6FW/d7yvcX4WqVqpqG/ACMEq8LglfAV5y508Gboy71mS3/RJwjZRwF4ZnFm9jyLgZkZtGOXCKORtsJGTa+7uYuGgrD76xMWhXipJ8tGn8A/C62z4X2Bl3rMrZktnPBOrjBKjLftS13PEGF74k+e3b3nKg+5usJ1c6lHD+IgDCJcydLqPQEfKSXrjuWvrkJBoichcQA57rMiUIplnYe7pWIj/GikiZiJTV1tb27LRREpRSw2WyX7pgUy0Hmgo3r1aiD7SjUz+onosC1maROVmLhoiMAb4JfFuPfKFVwOC4YIOA3T3Y9wH9RKR3N/tR13LHT6dbNVkXqvqUqo5Q1REDBw7M9icZxUiJljia22LcMmkp33l6ma/x/vPzK7nop2/4GmdUieqbmZVoiMhI4HbgBlWN73c6HRjtej4NBYYBS4FlwDDXU6oPXmP5dCc284Gb3PljgGlx1xrjtm8C5mkpZR8NIwdibnR7ZR6nfU82pUb8VzljzZ68xWeEk96pAojI88CXgAEiUgXcjddbqi8w29Udv6eqP1DVdSIyFViPV211m6p2uOv8GJgF9AImqWpXN6LbgRdE5H5gJTDR2ScCz4hIBV4JY3Qefq9hGHmiGNqNbG6pzEkpGqp6cwLzxAS2rvDjgfEJ7DOBY+ZpUNVKvN5V3e0twLdS+RdVWmMd9O3dK2g3jCLAr6Q7yvX/exoO09gS42Nnnxq0K5HHRoQHwNpdDVz4X28wa93eoF0xCsjtf1rD9v1NQbtRcKKQV//sA/O49pEFvsU3cdFWFm4uzk45JhoB8P7OegDe3lScL5VxhF/P3VyQ61bua+KVld7sAsFNGGgk477X1vP3E5cG7UZBMNEweG7J9uIdOFjEfSf+7cVVvsZn9f8GmGgYwBPztwDFNXCwq5G2KQsh7D7WoL65jbkbSncW3GQELSF+LsJ0sCVGXQHHvEQJE40ACPpjKwW6emdX1BzKeDbi/3hp1VFjDb7/zHJunVzG/jTXVYlyg3FPJPtdxflrj+Xy+2b3eHxJ5X6GjJvBnobDPnkUDCYaAVIqH1vQ7MhQNF5eseuo/W2uMbu9w+Q+TIRNnJ9d4q1HUgwrWfaEiYbPRDnZsbGV6WF1/0YxY6KRI7PXV7O6qj6zk1ziWwRjo3zD9Cr/rKlqoLaxeNqxDH8w0ciR700p44bH3snonCMzMkZLNYphBLBxhL94bBFfe+TttMNHrQRlGY3CYKLhM5bsZkcQelUKqyzW20JFRoaYaARAV3fOqgO2xniYmbBwa9AuZE0hctnRKxmnDjNrXTVDxs2gszPxDcvHb/7fN8t54q2KY+xRLQilnHvKyD/zNtYAML/cRoQbmZEocbNqmNzpUOW4Aonib+Z5gvGjL11QkOv7jZU0jKKkmNpfGg63U9/sDSx7cdnOFKE9CvnzCzX9zfRVuxkybgY79pdGCTyqb6iJhlGUFFP34E/97E0uvdcbWNbYEksRunB0NYT/97R1KUJmx/T3vfXXNu49mJfrFdErECpMNALAXubMWLnjAFUHinuUbS4EXajK5nVWjdaysMYRTDQCwDQjM/7yiXdZvv1A0G5kRxE97O6Nwrlo1R/e2cZFP32D6oMtuTnVA0GLabGSUjREZJKI1IjI2jjbGSIyW0Q2u//9nV1E5FERqRCR1SJyedw5Y1z4zW598S77p0VkjTvnUXGV0cniKAoilpBYySgz/OxlFOvojOTzeXW1VxVlJcjokU5J42lgZDfbOGCuqg4D5rp9gOvx1gUfBowFngRPAPCWib0Sb5W+u+NE4EkXtuu8kSniMHzmSw+9FbQLgZJTolxg/eiIomIUKcXUjtYTKUVDVRfgrdEdzyhgstueDNwYZ5+iHu8B/UTkHOA6YLaq1qnqAWA2MNIdO01VF6t3x6d0u1aiOIw4ahtbeWHpjoLG0ZGkD3upUerVHWEeEb5x70Eu/K/Xj5phtkTScN/Jtk3jbFXdA+D+n+Xs5wLxfQKrnK0ne1UCe09xRJ58fnw/fHY5415ek/H032GhvaPzg5UMw44lQuHlmcXbaY11MmdDTWCDEIupm3dP5LshPNFd0yzsmUUqMlZEykSkrLa2tAbM7XNrPMQiWhr45axybnz8nbx1s8yVIeNmHGMrkbQgY4J+45JlvvKZKQtz6SooshWNale1hPtf4+xVwOC4cIOA3SnsgxLYe4rjGFT1KVUdoaojBg4cmOVPMrqjqkmnV8gXa3c1ALD/UH5XRYvP9f39xCV5vXYitu5rornNG0NRyIQmlYA9t2Q7Cwo0+K6nRZhqDrYkFFw/SOqXiX1ByFY0pgNdPaDGANPi7Le4XlRXAQ2uamkWcK2I9HcN4NcCs9yxRhG5yvWauqXbtRLFEWmUwlRzzFlfTVNrfgd+/fbtSs6/c2Zer5kp2TYuxp+Xre4p0BbrpPpg6unDv/zQW/zj5LLsIsojd72yllsmLc36/E3VjazYkXn35nV7gispJi1xWCGhIKTT5fZ5YDFwoYhUicitwIPA10RkM/A1tw8wE6gEKoDfAz8CUNU64D5gmfu719kAfghMcOdsAV539mRxRJO4XE8h3uXxMzcw7uU1eb3mi8sK28CeDk++vSXQ+CcuSn/Swne37C+gJ/5w7SML+Ksn3g3ajbRIXvLJvIiRrcC0tHcc0564u75wY0/CQMoJC1X15iSHrkkQVoHbklxnEjApgb0MuCSBfX+iOIzkZLqsaRSYvzFprWTBEeBgi00d3oXV7x/L96aUsXDzPrY9+I0PbD9/Y2PCsM8s3uaPUwXGRoQn4Lkl248Zgfy7XHO8cd9bqfTnLgT/+uL7QbsQKPbqpMZPcVu4eV/aYX/abc6uqD5KE40E3PXKWv76yaOL6A+8vpHWWO5z5VjbnBFVwrqexqbqQzy35Niq1EwawoNtNA/nfU2GiYZhlBCxjs7M17TPgjfW7i14HF1U1BxKaI9OqSwyjgImGoEQrVfEyJSw5sgBfvlmOTc89g4b8tjbKdH7HORg0zDf/3ii4eWxmGgEQHRyQMGT7b0KzejckD3rdbs8segaFJoJ3dsK/L7DLy7b4bsYRUWA/MREIwBClo4YRuhpi3Vy+5/WcNNvw9sd2I8OLg2H232pXuwJE40AsLxLNAm68HLP9HW8WJbecq9+8nZ5Le9WJO9F9OjczTnH0VXKOdCUugt0MXcNvmXSUm547J1AfTDRKFJa2juoacxtkFGUP73QdGvuQWjerdjHqgwma3z63W0J7X4lksmqah6bX8HvFyYfBPnw7E1586GtozOrqrV8EfRiYJm8L4XCRCMA/PjEvzeljCvGz/UhpvSIdXRy0U9fZ2oIc8pB8bcTljDq8WBzjRC9NrZdKRZuymc7xIPdBup174p/VLxBF0V9wkQjCHr4SitqGlm/O/eeLZkMOvKDptYOWto7uf+19QWN5421e6msTdwF08iNYq72ScbzGaxVsz9BCWhNVUM+3QkFJhoh46sPL+Drjy4M2g0gmm0vP3h2OV/537cLcm2/cuT1zW3UN+d31t+eiEIGOQq9mPYmWO/8Lx5bVJC4dtcfpjGgKW5MNAKg9PJrpUWuOfJL753NpffOTivso3M38+KyYKr8wvoel0KJ6HMPzgusQTzlhIVGfinU1OhRYsOeg+yuP8w1Hz87aFciz+Pz8zcLcKpZeguf28/sw8jrYksR/Ca37msKJF4rafhF+EvXvnH9rxdyax7XnujsVGas3nNUj6lSaZQMK+Nnbggs7ihUZUUZEw0jI2IdneHpzuqYsngbt/1xBf9XVpUybBfNbTE2J5mzKF0amtuZvb467fD7D7WG7t5lQuGrfZIn9pYHCA8mGgEQlTrXbfuPnrJBVbngrtf52avp9YDyq6GuutHrtVKbQf/9u15Zm3O8t/1xBd+bUkZ1twbQRDndippDfPr+OUxZvD3neA2PfJYoTJTSJyfREJF/E5F1IrJWRJ4XkRNEZKiILBGRzSLyooj0cWH7uv0Kd3xI3HXucPZyEbkuzj7S2SpEZFwuvoaJqGY2u5ZNffrdbUxftbvHsMu31/GJe97MKCfexa76w0xYWAlk39iaKke/bX/u9cFd12iLdaYO6+qfC7V+d7ET1W+mJ6L6k7IWDRE5F/hnYISqXgL0AkYDPwceUdVhwAHgVnfKrcABVb0AeMSFQ0SGu/MuBkYCT4hILxHpBTwOXA8MB252YaNJgd6QoOru//n5lbR3JE8sV+7wRq4uzmIJ1O/+YSn3z9jA3obiXjbTSJ/yvY1Bu1BAolXMybV6qjdwooj0Bk4C9gBfAV5yxycDN7rtUW4fd/wa8VK8UcALqtqqqlvx1gq/wv1VqGqlqrYBL7iwkUYozlxTtnS6m7G/6ci4hMaW2FHHIk0R/IQugmxgvn9GYQeF5oUietY9kbVoqOou4CFgB55YNADLgXpVjblgVcC5bvtcYKc7N+bCnxlv73ZOMnvkyee7FeWGVYD3KusAuPPlNQF7kpoo1Xu3tqeuMsuVMN6OIH3K/kv0zmyLdTJk3IwPqmbDSi7VU/3xcv5DgQ8DJ+NVJXWn614mep6ahT2RL2NFpExEymprg68zfnj2JoaMm8EjSSZqC+PHFjSH23NfSjcKzN1Y40s8nxk/J6F9UQ+z0aYiCh04svXwrfJa5m3MvP2ti/2HWjMeN9E9HWhq9fLaj82vyNoPP8ileuqrwFZVrVXVduBl4HNAP1ddBTAI6GoxrQIGA7jjpwN18fZu5ySzH4OqPqWqI1R1xMCBA3P4SamZsLCSsm11PYbpmgr610mmhA7/p2dkQsQLe0VPOon596aU8Q9PZz92qCrFJIrFRC6isQO4SkROcm0T1wDrgfnATS7MGGCa257u9nHH56lXtzIdGO16Vw0FhgFLgWXAMNcbqw9eY/n0HPzNC/fP2MBNv12c0zWiXqVUKFSVjs783JtEHQTqm9sYMm4GL69IfzxHT3T5+urqnnuShZ1ifx2//NBbQbtQVOTSprEEr0F7BbDGXesp4HbgJyJSgddmMdGdMhE409l/Aoxz11kHTMUTnDeA21S1w7V7/BiYBWwAprqwkWL59joTiTT5j5dW89E7Zxbs+tvduJNk61JkSp1rvP/lrPLkgSJeF9nZqWyq9nouFbIhPNbRyZwN2VcPGf6R09xTqno3cHc3cyVez6fuYVuAbyW5znhgfAL7TKBwqYgP/PWTi7l31MWM+lRRtOEXlJeW56cEYOSPJ9/ewi9nlfPqj79Q0HgenRfuenzjCDYi3Ae25DhdRSmQj9JYfI6/EKW7YiowTohbaa+nXmFdK8Xtqj9c0IbwqgNHZh+I4rxhEXQ5a0w0AiDTxKeYEqt0KLXfGwRv+zwyfcLCSpZv77kDSRdWnRtubGp0n8n359B9fqioIiKmFkXM/TO8WW+3PfiNrM63NyM8WEnDL3IovpZC0ffoac0T28NGKTyXnijkk/njkvSXWTX8xUQjAKIwSCrfhO03b8rDXEbFvG7D/a+tP2aW4vjf+/kH57F8+4GCxX/nK6lnCAjTO1XM70J3rHoqAEKceQ4Mv6unmtpKYwR6tkxYtPWoJPnpd7eyo+5IVeiu+tIZzJYOYRKwQmMljZDyTsU+2nqYRTYqTHrH66VT6UblJtOFfFdD9dQDJ19RFXv1VCzu/euaIyzs+PVIWto7GDJuRk7XiKrMWEkjpHx7whK++clzAov/jbV783q9sVOWA9DYGvvA1hrzN7f/YtnO1IF6oLtIpCU+UU0ZIkCqKqFC3vp9GSz4lZpo5T6spOEXmnCzR7bU5jYBWi784NnlebxaYjbuOdKu4Eff/K7G1WIvIRhGIbGShs9462kUJg9ULJnafIxj8aPXVb7Fp6k1RixPc2+FmZ+/sfGDRbryTfzds7xBYbCSRgZYA7a/lFqJ4Mr/mcunfvZm0G4UnCff2hK0C6HiN/M2883fLMz6/Jb2Du57bT2H4qp+C4mVNIy0CWsPkaDEJVG0B1vaae/MrgODXx99ukRxOo+gyKbLbdcZ097PbZbkZ9/bzsRFW+nT+zhuH3lRTtdKBxONDMjXN5SvpLehuT11oAgTxZLdJ+9Jr6Qwc82eAnuSPlG8zz3xm/neOjZF9rOS0jVFf76WFUiFVU9FmC89ND9oF/JGshHhkUcS97T50XMrAnAmGuRawtlZF80xJMmS/Oa2mG+CkA4mGiGgNdbB6qrMGwYPRKqkEZ6XPl+kk7it2H6AEffPYfqq6C3U9Fa5P0vTdidhJ4ZupkOtMZrbeq7Oi2reo765/agJJYf/9yzG/Wl1gB4djYlGEHT7AO6Zvp4bHnuHHUUy+WC6HNXTJUkCnK3U5KM+fv2eg0c12mZTjdM1oeTSreEdHJfsVoVlMsxYgkGul9w9i/94KTwJaT6IfwxjJi096tj/hWitmZxEQ0T6ichLIrJRRDaIyGdF5AwRmS0im93//i6siMijIlIhIqtF5PK464xx4TeLyJg4+6dFZI0751Epkpa57mnPml1eKaPhcJRKDobhD9bbKlzkWtL4NfCGql4EfApvWdZxwFxVHQbMdfsA1+Ot/z0MGAs8CSAiZ+Ct/ncl3op/d3cJjQszNu68kTn6Gzg9ZVbD2jspPxyr937mAH77dv4Tnkz8j+qzLbTXU5elHqWf7TxXft7xqD7fbMhaNETkNOBq3BrgqtqmqvXAKGCyCzYZuNFtjwKmqMd7QD8ROQe4DpitqnWqegCYDYx0x05T1cXqVXJOibuWL4z81YK0w3amaqiKyHTfYaC5AJMJPvj6xrxf08idXKd2yTcVtspmSnIpaZwP1AJ/EJGVIjJBRE4GzlbVPQDu/1ku/LlA/BtS5Ww92asS2I9BRMaKSJmIlNXW5m9Fso0ZTJ99/p0zqT7YklbYTCWjWDUm2QfaVU2XTmXkhj0HE6694IswZ1DUiOrU2VHzuru/mZZSvvrw21nGm/udynZ8j9/kIhq9gcuBJ1X1MqCJI1VRiUh0VzUL+7FG1adUdYSqjhg4cGDPXheQ7QE3HIbpA08nzY6fajtbrv/1Qt/agvLVopaoYdfIP4K37kdUuGL83KBdSItcRKMKqFLVJW7/JTwRqXZVS7j/NXHhB8edPwjYncI+KIE9tKSbu8007Uk3sSp03rr7z+vp546fuaGwzkSM+Drvq39RPONrwsZ//Xlt0C4UPVmLhqruBXaKyIXOdA2wHpgOdPWAGgNMc9vTgVtcL6qrgAZXfTULuFZE+rsG8GuBWe5Yo4hc5XpN3RJ3rUBIpQnpJtohGqdzFGGrBgubP7lQffDIAL/dDelVYxoe2b4Hfr4+mZZCDzS10XA4XNPGpEuu04j8E/CciPQBKoHv4gnRVBG5FdgBfMuFnQl8HagAml1YVLVORO4Dlrlw96pqV6f2HwJPAycCr7u/0JLuy51pT4v65nZO7uPfjC+FmQOp599snQOORVX52avrg3bDSECur+tl983OjyMBkFNKpKrvAyMSHLomQVgFbktynUnApAT2MuCSXHw00qcrtzRx4dZgHckDfgzpKXQMdU1tPP3utgLHEn7y+ShbYx307d0rrbCHM+jFV0p5HhsRnkebteo7AAAUxUlEQVR6LEFE6KUq9j7n+VvuNUxdD4qHRPc1/p3cUptet9hET2fayvSbRdfsakh+7QAefVhmQTbRyCdpJEbeIkxJTi/utLpHLAE20uUnU1dlfW46GaJUVaWJer/58fp+LcvuwPnGRKMbK3ccCNoFI45CCKnp09GELbOSr/atQv2sB3wcKNoUV7rYE5IOFCYacayuqucvn3g3aDeKlCwWqckhcbeGdaNQLNh07ADiQr1u970Wvo4QJhpx1Bw8dt2DTKg/3M663cnrQbso1Au2amc9Syr353ydqI5e9hsrsfhDtt/Lpur0Z3QIK3VNbUG7cAwmGnnkR8+t4BuPLsr6/HxoSS71vfnib3//3jG2RAsRBcX+Q/n5EItFM4pV/FZXpc7A5YtC3cMwlpdNNAKgkFUnBb32MfuJ43p3S+6lnUKS9ayp3UfE58EXw4gaJho+01NCE0SG7421e7JORHMlVaKbi/5Zb6z88Ob66qBdKNqSUFQx0fALn1787unsM+9t54rxc5KGn7J4e2Ed6oH4WxKmNZCTYYmX4TezQyDa3fFvbgojEH5qE7gZRkaoWgahJ6ykYRhGqPG79/RfPLaIB3qYpblTld8vqPTRo3BhJY0ASPYNhKWCxnJZyYl1dCIi9DrObpKf5Es4ttQ2pRzAu273QdbtPsg1Hz876TW21Dblx6EIYiWNAOgsaO+pno83HG7n4TfLi74NoVC9yD5xz5tc9UA0FssxjuWpBZVFNYC3Js3VQvOJiUaJMX7Geh6dV8GsdXuDdiWSHG7voLYxPGNOihHr+ZY+fzthSepAecZEo8Q43O5NttYekSVHw1oesmlKCsfMNXuCdiEy7Drgf3d5a9MIgEKmN35Oa57r74hysjv0jplBu1C0vLJy11H7IrDCJhINDTmXNESkl4isFJHX3P5QEVkiIptF5EW3qh8i0tftV7jjQ+KucYezl4vIdXH2kc5WISLjcvU1aiTKzVbnOD9Wd8I+ejtbrIojOuw71Mq2fckblvcdagvtEsmlSD6qp/4FiO+f9nPgEVUdBhwAbnX2W4EDqnoB8IgLh4gMB0YDFwMjgSecEPUCHgeuB4YDN7uw0SSLl37j3mMnXPvJ1Pfz4IxHpou6+F0jo+pNwmgUlqBr2j734DxW7Ej+nJvbwrH4kN+ku9iU3+QkGiIyCPgGMMHtC/AV4CUXZDJwo9se5fZxx69x4UcBL6hqq6puxVtD/Ar3V6GqlaraBrzgwkYaIbdqmdb2/LVFdHQk9iTdPHpFTW4vtZUFDIC2WDTa1/zm7mnrgnYhIbmWNH4F/CfQ9dTPBOpVtStrUAWc67bPBXYCuOMNLvwH9m7nJLMfg4iMFZEyESmrrT12rvuwUdA2DR9zjd/8TfYz+hqGEU2yFg0R+SZQo6rL480JgmqKY5najzWqPqWqI1R1xMCBA3vwOjfylyAnvpDfPXKKfS1wIz2s+cfIhFx6T30euEFEvg6cAJyGV/LoJyK9XWliENC1knsVMBioEpHewOlAXZy9i/hzktlLlr0pBvMUkwyEcXCf4T8mauEi65KGqt6hqoNUdQheQ/Y8Vf02MB+4yQUbA0xz29PdPu74PPW+7OnAaNe7aigwDFgKLAOGud5YfVwc07P1N+oUYvryYk5Xt/bQG8cwokA6NQFB1BYUYpzG7cALInI/sBKY6OwTgWdEpAKvhDEaQFXXichUYD0QA25T1Q4AEfkxMAvoBUxS1XC2DGVAmNLpQk5nEjRhWMHQyA9ReE2DKg21tHfywOsbfY0zL6Khqm8Bb7ntSryeT93DtADfSnL+eGB8AvtMoDhGUcW9VFH4CPzAboMRZjKZn60Q37SEtH+hTSMSAIVMLDN5eZMFtTpkw4CbE6x1b5hoBII10nqYNhnpYJmYcGGiEQLW7joYSLyZa5eJnZEb2WSY6praC+BJfimEsIW1S7yJRsBU+jxVQPyU6GEv8YTcPSMLpr2fea/5DXuCyVRlwt4G/9e1CAoTjYBpbuvwJZ63ymtYvr0uUlM27CmhD7GQ/PbtLSzYFI6ZEor1md46uSzv1wxrQ7hNjR4AQSz3+p0/LDtqv7K2iavOP7OAMR7h+89k90H99ZPFs8JakDzoumRue/AbAXtiFANW0siAQtcxfq8AuZVk/Hru5h6no84ns9ZV+xKPkR1+VgOGtZ7eSB8TDcfbm2qZ9M5WX+JK9pHub2rzJf4ufregMqHdr7UoLPkoPI1pTH//9Lv+vPdGZoRVYE00HHPWVxd2QaK451/YBuj0rx31EeHrd4e/gTQKzNlQ41tcYa2nN9LH2jR8Jtf1NAxPdFdVNfi+OI8ld4afhFVgraTh8HMAUWNLOFYiC7qgke0tn/zuNm58/B0Wbt6XV38Mw0iNiYYPTF683cfYwpk7ySfl1d7Ylp11zQF7YmRKbWN+17gvZvYdCue9MtFwFH9SeyyZFjSCLpkETYn//LzgV2eTYmDj3sagXUiIiYajeHoMZTIzZ7DJYKwzu4GGzy/dkWdPDMNIFxMNvyjFokwKHp+/JWgXDKNoeO697ewuwGJt3THRKGFKvbrJMIqJprYONtcUfi67rEVDRAaLyHwR2SAi60TkX5z9DBGZLSKb3f/+zi4i8qiIVIjIahG5PO5aY1z4zSIyJs7+aRFZ4855VPyqQzIigWmeYfhPLiWNGPDvqvpx4CrgNhEZDowD5qrqMGCu2we4Hm/972HAWOBJ8EQGuBu4Em/Fv7u7hMaFGRt33sgc/O2RUpSjfA7u+8+XbHlVwygFshYNVd2jqivcdiOwATgXGAVMdsEmAze67VHAFPV4D+gnIucA1wGzVbVOVQ8As4GR7thpqrpYvRbbKXHXCoRc0tjrf7Ugf47kiUSiUbatjsfmVWR8rallVflwKSNmrN7je5yGEWb8yPvmZUS4iAwBLgOWAGer6h7whEVEznLBzgV2xp1W5Ww92asS2BPFPxavRMJ5552X3W8o8O3eHcIpobtrxvefWc46m5rDMIweyLkhXEROAf4E/Kuq9pTiJEqVNQv7sUbVp1R1hKqOGDhwYCqXDUf3m5lKMJp8nrbDMIzM8KOaPSfREJHj8QTjOVV92ZmrXdUS7n/XbGhVwOC40wcBu1PYByWwF4RiadPIqAotw+q2lvboLOBUDGzb78/U9YaRCbn0nhJgIrBBVR+OOzQd6OoBNQaYFme/xfWiugpocNVYs4BrRaS/awC/FpjljjWKyFUurlvirpV3gtKM/5m5IaCYoz/Lrd/4/Y68sXZv6kCG4TO5tGl8Hvh7YI2IvO9sdwIPAlNF5FZgB/Atd2wm8HWgAmgGvgugqnUich/QtbTcvapa57Z/CDwNnAi87v6KinxPx76/qY2KmvSmHzDJMIziwo+ZcbMWDVVdRPLM1zUJwitwW5JrTQImJbCXAZdk62Mm1B9uTxnmnYpozKo68lcL0woX9DQihmFEDxsR7nhpeeouo2OfWe6DJ7kT60xPDNIMZgSEabqRKaFvCDcyJ0zt7ZYmZcYin0uaExbZjLBG+DDRKGUsK5sRP5lqo94Nw0TDZ8KUTlv1VGnRkEa7nRFt/KjJMNHwmcbW8AyQsy63pcWnfvZm0C4YRYCJRgmTb814q7wmdSDDMAqHNYQbhSTf5Yzv/GEZi/M87sQwjHBholHCFGKcRl1TW96vaRhGeDDRKGE6CtASvnLHgbxf0zCM9PBjRLiJRglzqACN8k+/uy3v1zQMIzyYaJQwe0K4xodhGNljI8INwzCMUGGiYRiGYaSNiYaRV9KdLNEwjPxjI8INwzCMtBEfGjVCLxoiMlJEykWkQkTGBe2PYRhGWDnh+MIn6aEWDRHpBTwOXA8MB24WkeHBemUYhhFOTurTq+BxhFo0gCuAClWtVNU24AVgVMA+GYZhhJLje5V4SQM4F9gZt1/lbIZhGEY3Tj3h+ILHkfUa4T6RqFXnmO45IjIWGAtw3nnnZRXRotu/zBd+Pj+tsCccfxwt7Z1ZxVMqXHPRWXzl42cx4iNnMGdDNZcO7kdl7SHeqdhPXVMbS7fVMfbq83mnYh+/Hn0pY6cs59qLP4QINLfGmF9eywN/9Qm+/8xyLh3cjwPNbVzy4dP59Ef60+s44c5X1tAa857BWaf25crzz6ThcDu7DjTzhQsG8P7OegafcRJfuegsppbt5L3KOr7zuSF8+8rz+P6zy/nbK87jknNP58d/XMmXLhzIGSf34akFlQw4pS+xzk6+8YlzGDrgZFZVNfDqqt1c9KFTuey8/ohAa3snG/ce5OqPDeTPK3dx+onH09zWgaKc2vd4Bp9xIrEOpbElxoY9Bxl79flUHTjMJYNOB1UmvbONk/v24rGbL+edLfu4bHB//rh0O81tHQhCw+F2VlfVU9PYysBT+3LmyX3Y09BCw+F2BpzSl8vP68eb66u59QtD2X+oldVVDZw/8BRWV9Xzl5edS11TG2ed1pc1uw4CcNngfizcXMuu+sN874vns373QZZsreOGSz/Mk29tAeBvRgzi5L69WV3VwJqqBto6OvnisAFcfl5/BpzSh6s/NpAHZm5kR10zV39sII0t7Ty3ZAcA466/iE9/pD+PzavghOOPo6MThp9zKk1tHby8ooqzTzuBg4fbufL8M/nMkDNoao3x2po99DvxeK46/0xeX7uHwWecRJ9ex9Ea62DXgcPUH26n30l9WLWzni8OG8CaXQ3c9fWP09GpPLWgktFXDOZ/Zm78IP6m1hivrtrNtv3N3m8+rx8fOu0EPv2R/ry+di+xTuWKIf1pjXVyoLmdLTWHOO3E3lTWNjHs7FM4vtdxvFVeC8DFHz6NYWedwqjLzuXuaevYUdfM728ZwbT3d3HqCcdT19TKrHXV3PMXw6nc18TSrXUMPLUvHZ3KRR86jQ+d3pdn39vBjrpmPn/Bmew/1MbGvY1853NDWLe7gYf/5lJ+OauchsPtbNhzkOa2Dq67+EO8vamGf7/2Ql5YtpPVVfV8/qMDaGyNcemg01mweR83Xnouj8zZxOknHs8XLhhAU1uM0088nmVb67j5ivP48kVnccbJfQr+bUshJq3LFyLyWeAeVb3O7d8BoKoPJDtnxIgRWlZW5pOHhmEYxYGILFfVEanChb16ahkwTESGikgfYDQwPWCfDMMwSpZQV0+pakxEfgzMAnoBk1R1XcBuGYZhlCyhFg0AVZ0JzAzaD8MwDCP81VOGYRhGiDDRMAzDMNLGRMMwDMNIGxMNwzAMI21MNAzDMIy0CfXgvmwQkVpge5anDwD25dGdfGK+ZUdYfQurX2C+ZUtYfUvXr4+o6sBUgYpONHJBRMrSGREZBOZbdoTVt7D6BeZbtoTVt3z7ZdVThmEYRtqYaBiGYRhpY6JxNE8F7UAPmG/ZEVbfwuoXmG/ZElbf8uqXtWkYhmEYaWMlDcMwDCNtTDQcIjJSRMpFpEJExvkU5yQRqRGRtXG2M0Rktohsdv/7O7uIyKPOv9UicnncOWNc+M0iMiYPfg0WkfkiskFE1onIv4TItxNEZKmIrHK+/czZh4rIEhfPi24qfUSkr9uvcMeHxF3rDmcvF5HrcvXNXbOXiKwUkddC5tc2EVkjIu+LSJmzBf483TX7ichLIrLRvXOfDYNvInKhu19dfwdF5F/D4Ju75r+5b2CtiDzvvo3Cv2+qWvJ/eNOubwHOB/oAq4DhPsR7NXA5sDbO9gtgnNseB/zcbX8deB1vNcOrgCXOfgZQ6f73d9v9c/TrHOByt30qsAkYHhLfBDjFbR8PLHFxTgVGO/tvgR+67R8Bv3Xbo4EX3fZw95z7AkPd8++Vh2f6E+CPwGtuPyx+bQMGdLMF/jzddScD/+i2+wD9wuJbnI+9gL3AR8LgG96y11uBE+Pes+/48b7l5YZG/Q/4LDArbv8O4A6f4h7C0aJRDpzjts8Byt3274Cbu4cDbgZ+F2c/KlyefJwGfC1svgEnASuAK/EGL/Xu/jzx1mL5rNvu7cJJ92ccHy4HfwYBc4GvAK+5eAL3y11nG8eKRuDPEzgNL/GTsPnWzZ9rgXfC4hueaOzEE6Le7n27zo/3zaqnPLoeQBdVzhYEZ6vqHgD3/yxnT+ZjQX13xdjL8HL0ofDNVQG9D9QAs/FyR/WqGksQzwc+uOMNwJkF8u1XwH8CXQvInxkSvwAUeFNElovIWGcLw/M8H6gF/uCq9SaIyMkh8S2e0cDzbjtw31R1F/AQsAPYg/f+LMeH981Ew0MS2MLWrSyZjwXzXUROAf4E/KuqHgyLb6raoaqX4uXsrwA+3kM8vvgmIt8EalR1ebw5aL/i+LyqXg5cD9wmIlf3ENZP33rjVdE+qaqXAU14VT5h8M2L0GsXuAH4v1RBk/iQd99cO8oovCqlDwMn4z3bZPHkzTcTDY8qYHDc/iBgd0C+VIvIOQDuf42zJ/OxIL6LyPF4gvGcqr4cJt+6UNV64C28+uN+ItK1EmV8PB/44I6fDtQVwLfPAzeIyDbgBbwqql+FwC8AVHW3+18DvIIntmF4nlVAlaoucfsv4YlIGHzr4npghapWu/0w+PZVYKuq1qpqO/Ay8Dl8eN9MNDyWAcNcz4M+eEXR6QH5Mh3o6l0xBq89oct+i+uhcRXQ4IrGs4BrRaS/y31c62xZIyICTAQ2qOrDIfNtoIj0c9sn4n08G4D5wE1JfOvy+SZgnnqVt9OB0a5XyVBgGLA0W79U9Q5VHaSqQ/Den3mq+u2g/QIQkZNF5NSubbznsJYQPE9V3QvsFJELnekaYH0YfIvjZo5UTXX5ELRvO4CrROQk97123bfCv2/5aiiK+h9ez4dNePXjd/kU5/N49ZHteIp/K14941xgs/t/hgsrwOPOvzXAiLjr/ANQ4f6+mwe/voBXRF0NvO/+vh4S3z4JrHS+rQX+29nPdy97BV41Ql9nP8HtV7jj58dd6y7nczlwfR6f65c40nsqcL+cD6vc37qu9zsMz9Nd81KgzD3TP+P1MAqLbycB+4HT42xh8e1nwEb3HTyD1wOq4O+bjQg3DMMw0saqpwzDMIy0MdEwDMMw0sZEwzAMw0gbEw3DMAwjbUw0DMMwjLQx0TAMwzDSxkTDMAzDSBsTDcMwDCNt/j+BOs83W/qFcgAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "ncomp = 500\nl = delta_Data_len(Data, Z[:, :ncomp], vec[:, :ncomp], M, s)\nfig, ax = plt.subplots()\nax.plot(l)\nprint(np.median(l), l.mean(), l.min(), l.max())",
   "execution_count": 94,
   "outputs": [
    {
     "output_type": "stream",
     "text": "126.894620301 589.325757575 2.31503205509 9609.70162146\n",
     "name": "stdout"
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD8CAYAAACcjGjIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl0HNWd6PHvT4vlfcOCgA0YgrOQTBbHAbK8bCQsIRnyXsgMSc5AMuQwkzBvskwW85IJ2RggYYsTQkJYTcJiDMEOBhtjbLCNsS15lzfJtmTLkqx9X1u6748umUbdLXVX19r9+5yjo+7q6qpfd1fdX91bt26JMQallFIqVp7fASillAoeTQ5KKaXiaHJQSikVR5ODUkqpOJoclFJKxdHkoJRSKo4mB6WUUnE0OSillIozZnIQkQdFpF5E9sRMmykiq0Wk3Po/w5ouIrJIRCpEZJeIzI95z7XW/OUicm3M9A+IyG7rPYtERJz+kEoppdIjY10hLSIfAzqBxcaYd1vTfg00G2NuFZGFwAxjzI9E5LPA/wU+C1wI/NYYc6GIzARKgAWAAUqBDxhjWkRkC/Bt4HXgeWCRMeaFsQKfNWuWmTt3rq0PrZRSuai0tLTRGFOcyrwFY81gjHlVROaOmHwl8Anr8SPAOuBH1vTFJppxXheR6SJyujXvamNMM4CIrAYuE5F1wFRjzCZr+mLgC8CYyWHu3LmUlJSMNZtSSimLiFSlOq/dcw6nGWNqAaz/p1rTZwPHYuartqaNNr06wfSEROR6ESkRkZKGhgaboSullBqL0yekE50vMDamJ2SMuc8Ys8AYs6C4OKWakVJKKRvsJocTVnMR1v96a3o1cGbMfHOAmjGmz0kwXSmllI/sJoflwHCPo2uBZTHTr7F6LV0EtFnNTquAS0RkhtWz6RJglfVah4hcZPVSuiZmWUoppXwy5glpEXmc6AnlWSJSDdwE3AosEZHrgKPAl6zZnyfaU6kC6Aa+DmCMaRaRXwJbrfl+MXxyGvgm8DAwgeiJ6DFPRiullHLXmF1Zg2rBggVGeysppVTqRKTUGLMglXn1CmmllFJxNDkolYMGhwxLth5jcCicLQfKfZoclMpBf91cxQ+f3sXiTZV+h6ICSpODUjmouasfgJbuAZ8jCY+BwSHWHagfe8YsoclBKaVS8NuXyvnaQ1t57VCj36F4QpODUkql4EhTFwBNnf0+R+INTQ4Oau8dYO7CFbywu9bvUJRSKiOaHBx0pCF6ZHHvK4d8jkQppTKjyUEplVXaewcI68W9QaLJQSmVNaqaunjPz15k8aaUb1ugktDkoJTKGpVN3QC8tO+E7WW8WFbH64ebnAoptMYceE8ppXLJ9Y+WAlB56xU+R+IvrTkopZSKo8lBKaVUHE0OSiml4mhyyAIdvQO09+oYOUop5+gJ6SzwDz97EdATaEqV1bT5HULW0JqDUirUhoYMFfWdAPx65QGfo8kemhyUUqH2+7UVfPrOV9hf1+53KFlFk4NSuSwLhpkorWoBoLat1+dIsosmB6VykCB+hxAqfZFB2ntyq9OHnpBWKgcZwl9j8NI1D2xh85Fmv8PwlNYclMplojWIVORaYgBNDkoppRLQ5KBywpp9J5i7cAX1HXrSUqlUaHJQOWF4fP+yGu3uqFQqNDkopZSKo8lBKaVUHE0OSiml4mhyUEp5YmtlM7VtPX6HoVKkyUGFxlMlx/jdmnK/w1A2femPm/j0Ha/4GkNrdz93vniAwSG9CHAsmhxUaPxg6S7uWH3Q7zBUBrr6B31d/03Ly1j0cgVr99f7GkcYZJQcROS7IlImIntE5HERGS8i54jIZhEpF5EnRWScNW+R9bzCen1uzHJutKYfEJFLM/tISimVWI+VnCJacxiT7eQgIrOB/wQWGGPeDeQDVwO3AXcZY+YBLcB11luuA1qMMecBd1nzISLnW+97F3AZ8AcRybcbl1IqR2VBeT934Qr+a8lOv8MAMm9WKgAmiEgBMBGoBT4FLLVefwT4gvX4Sus51usXi4hY058wxvQZY44AFcAFGcallMoR2TY81NPbqv0OAcggORhjjgO3A0eJJoU2oBRoNcZErNmqgdnW49nAMeu9EWv+U2KnJ3jPm4jI9SJSIiIlDQ0NdkOnp3+Q6xeXcLxVe05ks8jgkN8hBF8W3M8hiCKDQyfvMxFWmTQrzSB61H8OcAYwCbg8wazDW1+i/G5GmR4/0Zj7jDELjDELiouL0w/a8uLeOl7ce4JbX9hvexkq2Cobuzjvxy/w7PbjfocSSNl0P4cg5rc7Vh/ki/e+xu7q8N7TOpNmpU8DR4wxDcaYAeAZ4MPAdKuZCWAOUGM9rgbOBLBenwY0x05P8B6lbNlf1wHA87trfY5EeSZA+W5/bXQMr4bO5AM9Dg0Z1u6vxwQxu5FZcjgKXCQiE61zBxcDe4G1wFXWPNcCy6zHy63nWK+/bKLfynLgaqs30znAPGBLBnEpZdux5m7mLlzBsh1a4wgzCcGJiIdfq+TrD2/luV3BPIDJ5JzDZqInlrcBu61l3Qf8CPieiFQQPafwgPWWB4BTrOnfAxZayykDlhBNLCuBG4wx/naGVjlrn3XE9/edwdlhe3y+NiCMgno0Hqu6JXrO80R7MIeRz+g2ocaYm4CbRkw+TILeRsaYXuBLSZZzM3BzJrEo/9zywj6mFBXwrU+cR01bD3NmTPQ7pKxxqKGTi+94hdu/9F6u+sAcx5artwlVY9F7SKuM/emVwwAMDsFdLx1k3fc/wdxZk3yOKjsctM6dvLT3hKPJ4aQANb80d/UDMHPSOJ8jUaDDZygHbTzUCEBdQKvJEG2iWbnHuyajspo2WqxCT41u/i9XM/+Xq/0OY0y5UufSmoMLQtDcmTOGDCxaU05Ld7SA/smze2ju6ueer8znhse2AfDQ1z7IJ99xqivrv2LRBuaeMpF1P/ikK8tXyi05XXMIw0krN3T2RXhp7wm/w/DEpkON3Ln6ILus/ubDTReryupOzvODpe4OV1DZ1G37vbm5hQZbcBri3JWTycHtbm4BasZN6PtLdvKNxSVUNnb5HYrrBrJkgLVk21TvwCDf/EupXu2vHJeTyWHYc7tq+fAta/wOw3OVTdGk0K1dJOOsPRAdynnLkSafI0nNwqd38cKeOv7n+X1+h6KyTE4nB4CatuCePFXeiG1dfP1wMwDtvZEkcwfHKwcbeHZH8AYTuPulg/z872V+h6EylPPJQeWmbGhsaurs8zuEOENDhrtfKuehjZXerzzBj1pa1cLHfr2Wrr7gJ/ug0eSglHJMjQ/3iB7tHN9tK/dztLmb3cfDOwCeXzQ5AEtKjvGvD29NaV5jDI+8Vklb90DKyz/U0MnOY612w1MuC3oHAqX8oMkB+OHSXbyc4j1ltx1t5ablZSx8ZlfKy7/4jle48p6Nb5rWFxnk+0/tDOy4KrkkR3s0p+3Vgw1c9/DWrOkC3h8Z4pt/KaWivtPxZWfDN6QXwaWpLxLt4TN8UZVdL+2tZ2lpNd39Ef7w1Q8AsPlwE0WF+bzvzOkZx6myg+vlcBor+MYjJfQPDtE/OERRQbDv5JtKd/VtR1t4YU8dTV39LPm3D3kQVbhozSFA/vm+1/nCiBqGckfYjn6dbvoKys1+fv9yOTct2+P6erz+uYPx7WZGk4NSyje3v3iQRzZVZbSMRAX/cPLPhkLaL5oclFLZQTOBozQ5qJw0fLFbtpu7cAW3vKBXT6v0aXJQOakxgBeQuWX4fhu5TG9ulD5NDmOoa+vlb9ur/Q7DFbrDZCdjDIvWlOfEwIpjCdI1LD39gyfH7goDTQ5j+Or9r/PdJ3fSmUWX34fh5uu5rLMvwv66dtvvr+/o487VB7nmwS1J59EDA+/9+NndfP2hrRw80THmvFcsWu/7SLuaHMZQ3x5tfhgKWdfH0YStG2euue7hrVx293qMMbYK8eGfd/ianFHpgYJnjlg1uY4UBnUsq2nnkdcqXY5odDmZHNzeHcJS9galr3u2uW3lftsngY0xbD4Sf7Jcf6sMubRPVjV1sbcm3VpeOAqInEwObtGDMP/UtPYwd+EKdgRgDKt71x3y5CSwX9tbWA5+wP2k+vHfrOOzi9anGEu4aHJwUJh2mmyzvrwBgMc2p39BVVh/Ns+3t6CXbokuhgvgrxuWckKTgwtGO6JLqR1YKZWylGpQAUhsYesIosnBJjvZ/9WDDbz9JysprcqNC7Bynd8nFFWMAB2tByiUUWlySFMmbZgbKhoBKKlscSocT60qq6M8hW54qajv6CUyOOTIsoJid3UbAzGf6ablWXarzLCUarxxlO7XifxEX1W46g2aHNI2sg3zgQ1HmLtwRdYVdIn826OlfOauVzNeTlvPABfcvIZfPLfXgaiC4XBDJ5///QZuXpH9Q1Wk0zpy//rDzF24gt6B3GxOTZSc9JxDlhveQW5fdQCA/hxIDk7p6I3eRW/NPjevFrW3B9o9umvuit7fw+nbUfpRkKw7UM8ehz7HH185BKTWtz8bxR5MhuyUg97sJ4iGhgx5eSHbkpS73NocEmSfrz0UvWVu5a1XuLRS94280DNIHbvCchGq1hwCqDtHq+DOsleajrXbdvdHGBwKx849mly5qM6No/Uf/223rfeF7TvX5OCgrZXRXki7qp1tWgibTA6MuvoiDAW48D3/p6v43pIdfoehfNSeYRNZcLfuN9Pk4JD2nojvSWFwyHDP2oqUBwkM2gVC3f0R3nXTKm5bud/vUEa1bEeN3yE4ZsCjROzmtuZUK81/Pr6d/3hsmzMLSyRcFYfMkoOITBeRpSKyX0T2iciHRGSmiKwWkXLr/wxrXhGRRSJSISK7RGR+zHKuteYvF5FrM/1QfrjoljW09Qy8aVpJZTNX3fuaZzGs3FPHb1Yd4JbnR+8x4/bFOHYX32kdkT2z/biD0YSXm03Tde3RET/vXXco7femF5eHJeIoq0ol5uU7a3huV61z8YRcpjWH3wIrjTHvAN4L7AMWAmuMMfOANdZzgMuBedbf9cC9ACIyE7gJuBC4ALhpOKGETeuI5PCjp3dRUuXdNQ3DV1939+s5i3SE5QShk1q6BsaeaaSQHflC+HoIBYnt5CAiU4GPAQ8AGGP6jTGtwJXAI9ZsjwBfsB5fCSw2Ua8D00XkdOBSYLUxptkY0wKsBi6zG5dSmfA7UaSz+iAWfAODuZdo0xWWY5FMag7nAg3AQyKyXUTuF5FJwGnGmFoA6/+p1vyzgWMx76+2piWbrlzidwGYqrKaNl6zripP5msPbeGDN7/kUUTeC8cv9Qa7PXmc8IOndvm27lQEMJePKpPkUADMB+41xrwf6OKNJqREEn03ZpTp8QsQuV5ESkSkpKGhId14Y5Zj+62jC0mhOyzoXeuuWLSBr9y/edR51h1ooKEjs/tBh2FAtOBHGPXaoabEL3iwa4x2X/Agdb4IUiyjySQ5VAPVxpjhvXcp0WRxwmouwvpfHzP/mTHvnwPUjDI9jjHmPmPMAmPMguLi4gxC91c4No3cEVuTCkOicIJXBVSOfJ0pCdt3YTs5GGPqgGMi8nZr0sXAXmA5MNzj6FpgmfV4OXCN1WvpIqDNanZaBVwiIjOsE9GXWNNUQNR39PodQuAFMeGHrCJrWyqFrhe15JTHVwvJ75Jpb6X/C/xVRHYB7wP+B7gV+IyIlAOfsZ4DPA8cBiqAPwPfAjDGNAO/BLZaf7+wpgVawh0vxUOD2LmCvgOvL2/ggpvX8GJZnd+hBNpwDzEv7zU+1qrqU2xuC+I5qMWbKjnW3G37/WPXAJ39zKvK6jjvxy+wvy75LUOD3ow7UkZjKxljdgALErx0cYJ5DXBDkuU8CDyYSSxOaOsZYNqEwlHnCdsPnKnhC/u2B+D2m0H2m1XRC/e2Hw3f97T7eBvvmTN97BkzKE/TeWt77wA/XVbGmTMPs/6Hn7K/Ug+t3nsCiO4v73jL1FE/b/BScWJ6hXSMT92+zvN1hq0d0guZnmDOhN2fo6vPm2tLymraHB8ePpLmVdJub7PG+nht3TauxUi2zLgi2ZsdL/ZgMmz7uiaHGE1d/Zz/05UcaexKOs+oJ/ICWD0PK23GinegroMrFm3gjtUHPVvnugP1LNl6bOwZx5BOudjQmfk5rlcPNljrDV6JHJZiQpPDCN39gywtHXtnCNtRwFiqmrpOVo2D4PpHSzNcQkj2wDHEng8Y7hiw2+ExvEYrrL720FZ++PSbrx9wu3D79J2Z31AqKPR+DiqUYjfcT96+jiHjzBj+YTkyCpN0uti61k01ZIWbn0a9n0NIDly05hAQK2IG/NpV3crchSvYV5u850M6Nh9u4l0/XXmyDXe4oHmqpPrkPE4MzhmsssOd+zlkJ/c/td/faxAOWNJt4vK7F5kmhxSVVrVw/eIShly6G+g9aytOPl65J9renu7Q1clO5P7u5Qq6+gfZdfzNPWkefq2SnhQH6UtlO01lU95Y0fimpKSSS/Z9bqhopKa1J+n7UimEvGqL9/uAwfOmnAAkIadockjRN/9Syot7T9DU5V1PmnUH0hsi5JYXRk8mwyfpYo1VxbWzc432nq/ev5nfrilPf6E+OHiiw+8QkvqvJTsdW1YQjqrt+PpDW2wNOe4GN3KQ31fra3LIIX9ef8TvEFyTtIBLo8ZzzYNb3tRT7ZK7vD8xmmo5bTCOHRWHNDew9kBD4G8MlUhYkrEmhwTsVrlD8ptnvS1Hmnl6W/pNV4lqVk6pakrePTrsHt54hAN1/tWynGrydFvYeitpckggLL0J7PL7RJcbYne8f/rTpje65Xq0Q471nX78N+tcXX8mHzPTzeFnf9/LpXePXsuysw6/m1XcEpa9T7uy5oCxkl11S/IxbLIwj2SFkb9LUH+mTMr3MB3EpBPp1x7awlumjufyd5/uWjxO0JpDAn5fVen12j9629ox58nSg7jQ6BvQW78mk+m2WdPaw7t+utKZYFJgDNS2BX+k45xMDn4X/rGe3x0/TMQrLrZ92+XWQdzu6jb21jhzPYefUmkCea2ikbkLV1A5yvAsyby498Sow7pkIkxH6G6obumhK8P7rqdSooStmUyblWxyc3+qbLI/VLHT3N6eP//7DYAzV2YH3TPbjwOwpbKZubMmpf3+8hOdToeUscMNqcUU1vN4Rxq7aO9xZgDAcKWGHK05ZKKlq9/vEAJry5HA34YjVFbvraOzL3LyeWyiNjh3gJLJYj51xyvOBOE5Q3kK17F88vZ1XHnPxozX1jswaKtFYPPhJj7+m7UpX6zqJK05pKh3IHpp9M/+vteV5Tt5XHW4oTOle1OoYPv3v2xzbdluHMn//O9lSV/zuym3tKqFj7/tjVsLLymp5uX99aO8IzMjv127tY+bn99HVVM3B0908N4zU7jnhoO05pCFdla38cV7X/M7DOWwthQLmFQK/thah1M1kIc2VjqzIBcsGnFVflmNsyPbJvPcrlp2hvRGWZocHJTJsZHTx1UV9cFrn1aZWbaj5uTjTLeXdM8lba0MZ5NhspPAXtVkntl23HazlN8dBTQ52BSyjgdv4naviSB1fumPuDRSos+c/IpTqWlkevQb1hPS6Ur6OUNYXmhy8MkNj7nXnpyOdAf3y0ZBSmZOSPuo2MXPn8kRepi6fo4M1enI/dhE9YR0Ana3yTCWMd/6a2pJKleO/OxKpwlgX207DR19FE8pcjGixF7ae4KWbu1x5zU7SVLE39P4WnNIwO6RZFCPQO3E1d47wIJfrWZbVbQ54U+vHHY4quDw+gD1oY2VfPDmlxK+lk6SsRP3NxaX8IOlu95UWKWyxjAdxcdKFnVIP46nNDmETGNnH0tL3b9ZzvajrTR29vPSvvTvK51qr5pc1xfJoO/6iBJ9xa5a5i5ckVXX4fh9QjYdY4VqNxntdPh+4enQ5JBAkI8qrl9cwvef2kltW/I7gfnt3/9SCkQLv5IQ9HIJURk0psMuDbGRqWz6jkfjVkOQH4lSk0PI1Fu3Ao0M2t9YvNrQfvXcPq7646aUrkRNx9+2V4fqqNJLyQ5s/uWBzXzlz697GwzBPtBy0qC1PWbTuTlNDsqWjt4BzrlxBS/tTd7stL8uOqBeS7ezzUzffXInz+447ugyM2WnTX44vxlj+N2aco41Zz6mVrIo1pc38tqhplHjGGu5fpXzx5q7GRpyttB1+thi+P7t24/a6/K753gbpVUtSV/345yPJoeASHdb9fvAuaK+E2OiJzhT1dhp//7bI2sKLV3+ntd4cuvRkzGVVrVkNGLq8dYe7lh9kK8/vDXl9yQ7QrWVpAJ8tFvV1MX/+vVa2/cd97pMfWLrMVtNvp/73YbAjWqgySFkwlpN7+6PsOBXiXvojKUvMphwaHM//ejp3Ty25SgAd64+kNGyhvNerwP3bEh18/A6ISRa2/Kdb1zxnexe0MP3Pdh0OHGtZ8z1JrsmzcX9qKsvQvOIjgFh3G1z8joHJzaMRBtdJjtcuiEF+WgPoLath62Vb1STu/rsF3yX3vVqoIYxH7axopGvXng2rQ43m6Ui2YlPO9t2sgJ0ZDOHnS1utHD+8/HtJx/fu+5Q/HvDeiQElDs8fI2ekFZjstMbIllbc0ZxjLHjfvEPzlWREyWGYKdG/6S6faQynxPNHE79Tg0dfXzy9nUcTbAt+HGVf2RwiAc3HEl5/jAmOk0ODvJ7WOIgqQnBbRCHjXVU9uimSp7cejRuup+/t8EkPOIXSb9W6UWizfSbOtLYxZHGLhZvqkxvvS79RE+WHOMXzyUfvj+lk/wBLy5ysllpLAH/zYA3Nj47Qw+H8SjGT/+9LHqfgn/+4Flvmp5LX2O2fdZMW2k6eyOjvu5Es6/f3bUzrjmISL6IbBeR56zn54jIZhEpF5EnRWScNb3Iel5hvT43Zhk3WtMPiMilmcaUqTA1WVyxaIMv6/W7rHByx7GbLPsCOOLr+vLGtN/jRSEUlH3KrySXaLXpfO1+fH9ONCt9G9gX8/w24C5jzDygBbjOmn4d0GKMOQ+4y5oPETkfuBp4F3AZ8AcRyXcgLkd09I1+hOC2ZMMhZNuRXBitHuUaDzvSLaO/t2Rn3LT9de1pN3e5WfBkspm6kbSGF+n2/uN0k2NnbySz4VZsyCg5iMgc4Argfuu5AJ8CllqzPAJ8wXp8pfUc6/WLrfmvBJ4wxvQZY44AFcAFmcSVqVR+1lQ2rkMNmQ9lsKEi8ZGg39c5ZBO/q++xnkpx3KykXTRTXE/Qe7uNpqmzL+XfLOnAe86FA4xdY7OTjGJrtNc8uIW3/2Qlz+2qGeUdzsr0nMPdwA+BKdbzU4BWY8zw4XY1MNt6PBs4BmCMiYhImzX/bCD2uv7Y96gRsqnGkEmhnI3nTUTgxmd2e7/eNOfxM5ceONHBB2xeLxPL7Y/Q48A1K4k8XVrN595zhivLHsl2zUFEPgfUG2NKYycnmNWM8dpo7xm5zutFpEREShoa9CY1uSxIR/tB4mbSDEI6PpRFt799cGPqXWH9kEmz0keAfxSRSuAJos1JdwPTRWS4RjIHGK4HVQNnAlivTwOaY6cneM+bGGPuM8YsMMYsKC4uziD04El7+AxXovBWNh79pyNI+e0P6yo41OBswRs7HpJTyTzdpXT1e9tOb0d3QGO0nRyMMTcaY+YYY+YSPaH8sjHmq8Ba4CprtmuBZdbj5dZzrNdfNtEtZjlwtdWb6RxgHrDFblyp+MO6CjcXb8u+2vaE092+/aBdOV6uO+bR16toyGDMqZHs/CydfRF+vfIAX7hnI0+VHHMslt7IYEYHANleCxp25+qDfoeQkBvXOfwIeEJEfgVsBx6wpj8APCoiFURrDFcDGGPKRGQJsBeIADcYY1xNpXuOJy6IM2X3JF9diC4Yc5PfPcO8dqSxi/9+do/fYbCxInoFfUdvhB8s3eVzNN7w+iJNp3oveVnbdiQ5GGPWAeusx4dJ0NvIGNMLfCnJ+28GbnYiFi9Vt3Q7cuJpYDD9/vKZVNOD2l7/+d/5c82GX/pt/O5jslF2PL4l/urvoArKpvtYiL4zu3T4jAx89La1rq8jKDvDSG4MHZFJ26vT1xz4Iai/NUBVhveaCPJns6NqlIEgE27HDu0uL++vd2ZBKdDkEFJ+n8zt6A3WfaJf2uf+ThNJcKT/6kH7veZG/oLHW1O/D0DJKDeGSW3d6W0/D22szGh9ueSm5WV+h+AITQ6JpFDwenUklCwUv5uGfvNiZvcwCKP71h+Om3bNg/b7Trj1C7p1gVu6Sz3c0JVWwssWJxKczwhjBw5NDokEuA7sd41hWH8AxhVaX+7ttS61rYlPYtr9SVxJ8AHadN24s1lde+503nD61qjp0lFZQ2LZjuMU5DmTy4OSYDI9mf8vD8QftfsxGF6AjyV8Fftb2LsJUTi/2ERR2/ko96dxvwg3aHJIxGbh6ea2/O0ndgBwzqxJ0XW5tyrPjLyVohOW7Uh/7JlUv8uuEHS19WK7sNNs5Vc5H4zDoHDSZqWAG3niMCgbe0AqH57acazV7xBS4tZNiML2k6/cU8vmI81+hwHAugPe9TJyitYckiitaqatJ1g9ctLV1j3A1Anp/cQ9/YPsq3PnIkEnOXkk6lehF7ban9dNdu29EVbsqrX9/n//yzYHo8nMcM0/TDQ5JPHFezc5vsxaG1dlJqvCp1I4vvcXL3Lz/353Wuv7wdKdPLerlt9c9Z603qeUGx59vcrR5XlR47WznweRNivZZGcja+5ybgydVP34b+kNz7DnePS2o2NdkDZW00WYmp3CdgTvhDDfzyETIT3H7YucSw4v7LZfTfVDskI4VIWvG3f08qFwC0O5EtYePip4ci45JLq14kgPuzzO+i0v7Lf93sON0bvLaRngvWTnoOxf55BBMMmW6fwiVY7KueSQivZeb7ssujIAm89iy0s3rqtwq0fOaJJ1vdVErbKRJgebXj/sXBe57wSoJ0PlKAOKDSurafMgEmVHqonKj+SqwiXnkoObbfV2l737ePLCVsT5duRMl3fForGH1nb7yv/sOKHq32c4muEoq+nS2lX45FxyyHWDDpXaYTohHsvPuF8/3OTJelJJnHuT3HkwpeUbeOd/r2Txpkrby1DBp8kh4Lr6IhkPzxzHgsMEAAAQNklEQVRr5Z46x5aVNTw6qq1u8WaEUrcvuDJEx8X66bLsGJpaJaYXwTnIjapz8ts22ltZX4b39VXOcaupZUtAhoxQ4ZZzNQc3i8Vsalcd67yE3+nF0e/apw/jxuYS1E0wjMcjjZ3eX7QaJDmXHLKH/b2t14H7XqcjDBdmVdR32H5vkAq+oH7XAQ1rVP/0R+eH0AkTTQ45KJUCJFubnhJ+KgPLd4bryvlkQlgGB9bwBae5KueSQ7YWekEWiu88jIe2Pglq7UQ5K+eSg5/CdOGR172agnDbUbvslpVhLWP/53n7w7+o8NDk4KGgXLiVylXQm8bqk59mbaCqafQq+tt+8kJay7M7PlXCWoxoc0w6Ht9y1O8QlAdyLjmE59h9LOEqzv714a1+hwAkaRIJ11c5uoB+luOt3lzjMZYwtHAGRc4lh2zxb4+Wet7rKBPZOLhgpnKp7f7Lf37d7xCA8Dbl+UGTg4MyGZIgXYcauthQ3ujZ+rJFkE6OazmlgkyTg7IlOEWsM4YyOKQMUL4JzHktv82aXOR3CKGXe8khQDtyppwcc8lNQW8+6eiLcM/aQ36H4YiAf9WeufDcmX6HEHq5lxyyyB9fyY4CzUvpHhu4eVtZLchVkGlyULak05SyYnctvQPhPCH9zb9u8zuEtGjCicqiBoI4R1Poiu6EnEsO2bzReCmdQuiJrcfcCyQAglQgB+n8h3LHKwfrPVmP7eQgImeKyFoR2SciZSLybWv6TBFZLSLl1v8Z1nQRkUUiUiEiu0RkfsyyrrXmLxeRazP/WEol5kbh2d0fnC7FQUpUyh1u32VxWCY1hwjwX8aYdwIXATeIyPnAQmCNMWYesMZ6DnA5MM/6ux64F6LJBLgJuBC4ALhpOKGo4NIj1DccsTlAm/YsUnZk0rMuHbaTgzGm1hizzXrcAewDZgNXAo9Ysz0CfMF6fCWw2ES9DkwXkdOBS4HVxphmY0wLsBq4zG5cYwlSP/ecEaAyMFBH1i7EsrIsHHf6q+/o9TuE0PKq5uDIneBEZC7wfmAzcJoxphaiCURETrVmmw3ENj5XW9OSTc86kcEglUy5KeLVnqWS+uVze3lgwxFX15HoVw56l+qgyfiEtIhMBp4GvmOMGe0S4USH7GaU6YnWdb2IlIhISUNDQ/rB+uxosze9DLyQVv1LK2sqhtuJAQhUbdVpXu1OGSUHESkkmhj+aox5xpp8wmouwvo/fGq9Gjgz5u1zgJpRpscxxtxnjFlgjFlQXFxsM2Zbb1PKcV+5f7PfIWStbD6f41UZlklvJQEeAPYZY+6MeWk5MNzj6FpgWcz0a6xeSxcBbVbz0yrgEhGZYZ2IvsSapgJMz90o5Q+v9rxMzjl8BPgXYLeI7LCm/T/gVmCJiFwHHAW+ZL32PPBZoALoBr4OYIxpFpFfAsNjOv/CGNOcQVxKqRyX6PRCQ0cfFfWd3gcTUraTgzFmA8mT2MUJ5jfADUmW9SDwoN1YlPf05J4Kmw/f+rJ2SEiDXiGtlMoJ2ZIYvGrSzbnkoJyh5xxUkGVzxfam5WWerEeTg0vWHYgf/ySbNti0UkMWfW4VDtncW8krOZccvDri3X601ZP1+EUrDirIsulAzC85lxy8kmjbrGvPniEDJJ26gyYSpUJHk4NLevojfofgLi3wVYD1DARnpNywyrnk4FWZ9uf1HgwRYMOBug7P19kfCeeNflR4rS9v9DuE0Mu55JDr/vTqYUeWoxUHpbKbJgellFJxNDkoW7S3klLZLeeSgxZqzkirt5JSKnRyLjkoZ5zIom65Sql4OZccGjv7/Q4hKxy2ed9kpVQ45FxyUEopNTZNDkoppeJoclBKKRVHk4NSSqk4mhyUUkrF0eSglFIqjiYHpZRScTQ5KKWUiqPJQSmlVBxNDkoppeJoclBKKRVHk4NSSqk4mhyUUkrF0eSglFIqjiYHpZRScTQ5KKWUiqPJQSmlVBxNDkoppeJoclBKKRUnMMlBRC4TkQMiUiEiC/2ORymlclkgkoOI5AP3AJcD5wNfFpHz/Y1KKaVyVyCSA3ABUGGMOWyM6QeeAK70OSallMpZQUkOs4FjMc+rrWlKKaVizJpc5Ml6CjxZy9gkwTQTN5PI9cD1AGeddZatFa3/4SfZcayVyeMLKMgThgy0dPVzrLmb9545neOtPVQ1dfO595yOCJRUtjBlfAFVTd3MmjyO2TMmUFHfyVumTaC1u5/jrT0U5ecxe8YETp06nq6+CH0DQ9S29dAXGWJ8YT7FU4po6uxnyvgCTp82nvXljRTkCf8wZxrbj7by1uJJvGv2NJo7+2nrGaCquZvJRfkUFeRz9ikTOdHey+AQtPUMUJgvnH/GVI42dSMCR5u7iQwZPvLWWRTm5zFlfAEF+cKsyUUs21HDjImFFBXk09DRS29kiHeePhUBJhUVIAKllS1MLMrnM+efxroDDfQNDDKuII9jzT1MGJfP+8+cTmFBHgV5QkNHH4X5eZTVtDF1QiHzz5rB9ImFzJkxkc6+CMdbeth9vI0Lz5nJuII8Djd0MXfWRKaOL6S5q5/iKUWUn+hk3mmTqWrqpqmzj02Hm/jM+adhDBRPKaIgT6hu7eFtp02htbuf/bUddPVH+Oh5s9hypJmJRQW0dvdT3dLD599zBuMLo8c37b0RCvKE2rZezjt1MgC9A4MUTymip3+Qv++q4ZNvP5WO3gh7jrdxTvEkpowvoG9giCnjC1hVdoLegUH+z/zZjC/Mp66tlzwRRODUqUVEBg31HX0U5AmRIcOpU4oQgcaOfnojg5w+bTwAMyaOo6qpm/beAfJEeOfpU2jtHmDK+AK6+geJDA7R1jPA1PGFTJ9YSGVTN1VNXcw7dQrjCvKo7+iloaOP+WfNYPuxVt75linMmDSOjt4IzV199A1E3/+xtxXT3juAMbDpUBPzz55Bc1c/hflCZ1+EBWfPpLy+g6njC+mLDDI4BIcaOvnwW0+hvL6T/DyhsaOPM2dOZEJhPtMmFNI9MMi0CYWcaO+lb2CIGZMKGYgYmrr6mFxUwLnFk2nq7KOivpN3nTGNcQV5TBiXz5YjzeTnCTMnjaOyqYsFZ89gXEEeR5u6KczPo2dgkAmF+fRFhmjq7GPm5HGML8inIF8or+9k1qQiBo2hqCCPM6ZNYMgY2noGmDgun+rWHrYcaeaic09h3qmT2XO8jcnjC9hY0cjH3lZMU2c/806bzOCQoaqpmwmF+UyfWEhBfh6dvRHy82DiuALK6zt5a/EkTrT30tTZz7nFk5gwroDTphSRJ8KemjaGDEybUMj4wjxOmzKeF/ee4C3TxnPG9PEMDUFVUxeTigpo7urnfWdNp7Kxi4aOPgaHDO94y1QmFeXT3NXPpKJo2TJ1QiGt3QPsOd7GWadM5HhLD8VTiugdGORtb5lCS1d0f8/PEzp6IxRPKaKoII/mrn62HGnmI+dF9+mO3gHeeupk2roHTu7bXhBj4spgz4nIh4CfGWMutZ7fCGCMuSXZexYsWGBKSko8ilAppcJPREqNMQtSmTcozUpbgXkico6IjAOuBpb7HJNSSuWsQDQrGWMiIvIfwCogH3jQGFPmc1hKKZWzApEcAIwxzwPP+x2HUkqp4DQrKaWUChBNDkoppeJoclBKKRVHk4NSSqk4mhyUUkrFCcRFcHaISANQZfPts4BGB8NxUlBjC2pcoLHZpbGlL6hxQWqxnW2MKU5lYaFNDpkQkZJUrxL0WlBjC2pcoLHZpbGlL6hxgfOxabOSUkqpOJoclFJKxcnV5HCf3wGMIqixBTUu0Njs0tjSF9S4wOHYcvKcg1JKqdHlas1BKaXUKHIqOYjIZSJyQEQqRGShR+t8UETqRWRPzLSZIrJaRMqt/zOs6SIii6z4donI/Jj3XGvNXy4i1zoU25kislZE9olImYh8Owjxich4EdkiIjutuH5uTT9HRDZb63jSGt4dESmynldYr8+NWdaN1vQDInJpJnGNiDFfRLaLyHNBik1EKkVkt4jsEJESa1pQtrfpIrJURPZb29yHghCbiLzd+r6G/9pF5DsBie271j6wR0Qet/YNb7Y1Y0xO/BEdCvwQcC4wDtgJnO/Bej8GzAf2xEz7NbDQerwQuM16/FngBaJ3xrsI2GxNnwkctv7PsB7PcCC204H51uMpwEHgfL/js5Y/2XpcCGy21rcEuNqa/kfgm9bjbwF/tB5fDTxpPT7f+p2LgHOs3z/fod/1e8BjwHPW80DEBlQCs0ZMC8r29gjwDevxOGB6UGKLiTEfqAPO9js2ordKPgJMiNnGvubVtubIFxqGP+BDwKqY5zcCN3q07rm8OTkcAE63Hp8OHLAe/wn48sj5gC8Df4qZ/qb5HIxzGfCZIMUHTAS2ARcSvcCnYOTvSfQ+IB+yHhdY88nI3zh2vgxjmgOsAT4FPGetKyixVRKfHHz/PYGpRAs6CVpsI+K5BNgYhNiIJodjRJNNgbWtXerVtpZLzUrDX/SwamuaH04zxtQCWP9PtaYni9H12K0q6PuJHqX7Hp/VbLMDqAdWEz3aaTXGRBKs4+T6rdfbgFPciMtyN/BDYMh6fkqAYjPAiyJSKtF7rkMAfk+iNfYG4CGrOe5+EZkUkNhiXQ08bj32NTZjzHHgduAoUEt02ynFo20tl5KDJJgWtK5ayWJ0NXYRmQw8DXzHGNM+2qxJ4nA8PmPMoDHmfUSP0i8A3jnKOjyLS0Q+B9QbY0pjJwchNstHjDHzgcuBG0TkY6PM62VsBUSbV+81xrwf6CLaVBOE2KIrjLbd/yPw1FizJonB0discxxXEm0KOgOYRPR3TbYOR+PKpeRQDZwZ83wOUONTLCdE5HQA63+9NT1ZjK7FLiKFRBPDX40xzwQtPmNMK7COaNvudBEZvnth7DpOrt96fRrQ7FJcHwH+UUQqgSeINi3dHZDYMMbUWP/rgb8RTaxB+D2rgWpjzGbr+VKiySIIsQ27HNhmjDlhPfc7tk8DR4wxDcaYAeAZ4MN4tK3lUnLYCsyzzvSPI1p9XO5TLMuB4Z4M1xJt6x+efo3VG+IioM2qzq4CLhGRGdbRxCXWtIyIiAAPAPuMMXcGJT4RKRaR6dbjCUR3kn3AWuCqJHENx3sV8LKJNq4uB662enGcA8wDttiNC8AYc6MxZo4xZi7RbehlY8xXgxCbiEwSkSnDj4n+DnsIwPZmjKkDjonI261JFwN7gxBbjC/zRpPScAx+xnYUuEhEJlr76vB35s225tSJnDD8Ee1lcJBo+/WPPVrn40TbCweIZvDriLYDrgHKrf8zrXkFuMeKbzewIGY5/wpUWH9fdyi2jxKtXu4Cdlh/n/U7PuA9wHYrrj3AT63p51obdQXRqn+RNX289bzCev3cmGX92Ir3AHC5w7/tJ3ijt5LvsVkx7LT+yoa3cb9/z5hlvg8osX7XZ4n26AlKbBOBJmBazDTfYwN+Duy39oNHifY48mRb0yuklVJKxcmlZiWllFIp0uSglFIqjiYHpZRScTQ5KKWUiqPJQSmlVBxNDkoppeJoclBKKRVHk4NSSqk4/x9UR0W9E3Q/vwAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "var_explained = vals[:ncomp].sum()/totvar",
   "execution_count": 107,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "ncomp",
   "execution_count": 108,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 108,
     "data": {
      "text/plain": "1055"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "WaveletFeatures.ncompsForTolerance(vals, tol=0.99)",
   "execution_count": 104,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 104,
     "data": {
      "text/plain": "5"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "totvar = vals.sum()",
   "execution_count": 106,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "",
   "execution_count": null,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "var_explained = []\nfor ncomp in np.arange(5, 3000, 50):\n    l = delta_Data_len(Data, Z[:, :ncomp], vec[:, :ncomp], M, s)\n    vars = vals[:ncomp].sum()/totvar\n    var_explained.append([ncomp, np.median(l), l.mean(), l.min(), l.max(), vars])",
   "execution_count": 114,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "import pandas as pd",
   "execution_count": 112,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df = pd.DataFrame(np.array(var_explained), columns=['ncomp', 'median', 'mean', 'min', 'max', 'explained_var'])",
   "execution_count": 115,
   "outputs": []
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df.head()",
   "execution_count": 120,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 120,
     "data": {
      "text/plain": "   ncomp       median          mean          min           max  explained_var\n0    5.0  5336.252777  28577.553909  2343.339213  7.372932e+06       0.992006\n1   55.0  2130.934527   9749.476386   318.607789  2.912553e+05       0.999826\n2  105.0  1300.875342   6337.828447   138.797217  1.554230e+05       0.999933\n3  155.0   878.173163   4474.958087    86.786669  9.810235e+04       0.999968\n4  205.0   644.139044   3303.645053    38.675692  6.670458e+04       0.999984",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>ncomp</th>\n      <th>median</th>\n      <th>mean</th>\n      <th>min</th>\n      <th>max</th>\n      <th>explained_var</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>5.0</td>\n      <td>5336.252777</td>\n      <td>28577.553909</td>\n      <td>2343.339213</td>\n      <td>7.372932e+06</td>\n      <td>0.992006</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>55.0</td>\n      <td>2130.934527</td>\n      <td>9749.476386</td>\n      <td>318.607789</td>\n      <td>2.912553e+05</td>\n      <td>0.999826</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>105.0</td>\n      <td>1300.875342</td>\n      <td>6337.828447</td>\n      <td>138.797217</td>\n      <td>1.554230e+05</td>\n      <td>0.999933</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>155.0</td>\n      <td>878.173163</td>\n      <td>4474.958087</td>\n      <td>86.786669</td>\n      <td>9.810235e+04</td>\n      <td>0.999968</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>205.0</td>\n      <td>644.139044</td>\n      <td>3303.645053</td>\n      <td>38.675692</td>\n      <td>6.670458e+04</td>\n      <td>0.999984</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df.columns",
   "execution_count": 126,
   "outputs": [
    {
     "output_type": "execute_result",
     "execution_count": 126,
     "data": {
      "text/plain": "Index(['ncomp', 'median', 'mean', 'min', 'max', 'explained_var'], dtype='object')"
     },
     "metadata": {}
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "fig, ax = plt.subplots(2, sharex=True)\nax[0].plot(df.ncomp, df['median'],'.')\nax[1].plot(df.ncomp, df.explained_var, '.')\nax[0].set_ylabel('median loss')\nax[1].set_ylabel('explained_var')\nax[0].set_ylim(0, 100)\nfig.savefig('loss_comp.pdf')",
   "execution_count": 139,
   "outputs": [
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "<Figure size 432x288 with 2 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAD8CAYAAACl69mTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuYXVWZ5/HvLzeuwYSkQCCQCx1hsBtCUkLZOs1NMOAFUFQYQAbFODYo6tgjSD8QmGFQ28s0SGNHDIQ2glckzQQBM+FiSylVEEIAIyFSEBIhhEhQwKSod/7Y68ChOHU5VWefW/0+z3Oes8/al/Uuzkm9rLXX3lsRgZmZWR5G1ToAMzNrXk4yZmaWGycZMzPLjZOMmZnlxknGzMxy4yRjZma5qeskI2mhpGckrSoq21XS7ZIeTe8TU7kkXS5pjaSVkmbXLnIzM4M6TzLAtcDcXmXnAcsiYiawLH0GOBaYmV7zgKuqFKOZmfWhrpNMRNwFPNer+HhgUVpeBJxQVH5dZNqBCZL2qE6kZmZWyphaBzAEu0fEBoCI2CBpt1S+F/Bk0XbrUtmG3geQNI+st8NOO+00Z//99883YjOzJtPZ2flsRLQMtF0jJpm+qERZyXvmRMQCYAFAa2trdHR05BmXmVnTkdQ1mO3qerisD08XhsHS+zOpfB2wd9F2U4D1VY7NzMyKNGKSWQKckZbPAG4qKv9ommXWBjxfGFYzM7PaqOvhMknXA4cDkyWtAy4Cvgz8UNLHgSeAD6XNlwLHAWuAF4Ezqx6wmZm9Tl0nmYg4pY9VR5XYNoCz843IzMzK0YjDZVYlnV2buXL5Gjq7Ntc6FDNrUHXdk7Ha6ezazKlXt7O1u4dxY0ax+Kw25kydWOuwzKzBuCdjJbWv3cTW7h56ArZ199C+dlOtQzKzBuQkYyW1zZjEuDGjGC0YO2YUbTMm1TokM2tAHi6zkuZMncjis9poX7uJthmTPFRmZkPiJGN9mjN1opOLmQ2Lh8vMzCw3TjJmZpYbJxkzM8uNk4yZmeXGScbMzHLjJGNmZrlxkjEzs9w4yZiZWW6cZMzMLDdOMmZmlpuqJBlJ+0raLi0fLukzkiYM43j7SVpR9Noi6bOS5kt6qqj8uMq1wszMylWtnsxPgFck/RXwXWA68P2hHiwiVkfErIiYBcwhe9zyjWn1NwvrImLpcAM3M7Ohq1aS6YmIbuBE4P9ExOeAPSp07KOAxyKiq0LHMzOzCqlWktkm6RTgDODmVDa2Qsc+Gbi+6PM5klZKWiip5C2EJc2T1CGpY+PGjRUKw8zMeqtWkjkTeDtwaUT8XtJ04HvDPaikccD7gR+loquAfYFZwAbg66X2i4gFEdEaEa0tLS3DDcPMzPpQlefJRMTDwGcAUu9ifER8uQKHPha4LyKeTvU8XVgh6Tu81msyM7MaqNbssjsk7SJpV+AB4BpJ36jAoU+haKhMUvF5nhOBVRWow8zMhqhaw2VviogtwAeAayJiDvCu4RxQ0o7A0cBPi4q/KulBSSuBI4DPDacOMzMbnmo9fnlM6mV8GLigEgeMiBeBSb3KTq/Esc3MrDKq1ZO5BLiVbKrxvZJmAI9WqW4zM6uRap34/xGvzQAjItYCH6xG3VZZnV2baV+7ibYZk5gzteQMcTOzV1UlyUiaAlwBvAMI4JfAuRGxrhr1W2V0dm3m1Kvb2drdw7gxo1h8VpsTjZn1q1rDZdcAS4A9gb2Af09l1kDa125ia3cPPQHbuntoX7up1iGZWZ2rVpJpiYhrIqI7va4FfBVkg2mbMYlxY0YxWjB2zCjaZkwaeCczG9GqNbvsWUmn8do1LacA/t/gBjNn6kQWn9XmczJmNmjVSjIfA74FfJPsnMyvUpk1mDlTJzq5mNmgVWt22RNk9xgzM7MRJNckI+kKsp5LSRHxmTzrNzOz2sq7J9OR8/HNzKyO5ZpkImJRnsc3M7P6Vq0pzGZmNgI5yZiZWW6cZMzMLDfVundZC/AJYFpxnRHha2XMzJpYtS7GvAm4G/gF8EqV6jQzsxqrVpLZMSK+WMkDSnoceIEsaXVHRGt6vPMPyHpMjwMfjojNlazXzMwGr1rnZG6WdFwOxz0iImZFRGv6fB6wLCJmAsvSZzMzq5FqJZlzyRLNS5K2SHpB0pYc6jkeKFybswg4IYc6rJfOrs1cuXwNnV3uNJrZ61Xr3mXj8zgscJukAP41IhYAu0fEhlTnBkm7ldpR0jxgHsA+++yTQ2gjhx9kZmb9qdY5GSRNBGYC2xfKIuKuYRzyHRGxPiWS2yX9drA7poS0AKC1tbXPe6vZwEo9yMxJxswKqjWF+SyyIbMpwAqgDbgHOHKox4yI9en9GUk3AocAT0vaI/Vi9gCeGXbw1q/Cg8y2dff4QWZm9gbVPCfzNqArIo4ADgY2DvVgknaSNL6wDBwDrCJ7xPMZabMzyKZOW44KDzL7/DH7eajMzN6gWsNlL0fEy5KQtF1E/FbSfsM43u7AjZIga8P3I+Lnku4Ffijp48ATwIeGH7oNxA8yM7O+VCvJrJM0AfgZ2fmTzcD6oR4sItYCB5Uo3wQcNeQozcysoqo1u+zEtDhf0nLgTcDPq1G3mZnVTt5PxtwlIrakK/ELHkzvOwPP5Vm/mZnVVt49me8D7wU6ya5rUdG6AGbkXL+ZmdVQ3k/GfG96n55nPWZmVp/yHi6b3d/6iLgvz/rNzKy28h4u+3p63x5oBR4gGzI7EPg18M6c6zczsxrK9WLMiDgiXXzZBcyOiNaImEN2MeaaPOs2M7Paq9YV//tHRGFWGRGxCphVpbqtRnx3ZjOr1sWYj0i6Gvge2ayy04BHqlS31YDvzmxmUL2ezJnAQ2T3MPss8HAqsyZV6u7MZjbyVOuK/5clfRtYGhGrq1Gn1ZbvzmxmUL1b/b8f+CdgHDBd0izgkoh4fzXqt+or3J25fe0m2mZM8lCZ2QhVrXMyF5E97+UOgIhYIWlaleq2GvHdmc2sWudkuiPi+SrVZWZmdaJaPZlVkv4LMFrSTOAzwK+qVLeZmdVItXoynwbeCvwFuB7YQjbLbEgk7S1puaRHJD0k6dxUPl/SU5JWpNdxFYnezMyGpFqzy14ELkivSugG/ntE3Jcew9wp6fa07psR8bUK1WMV1tm12ZMBzEaQas0uawW+BEwrrjMiDhzK8SJiA7AhLb8g6RFgr+FHannyBZpmI0+1hssWA9cCHwTeV/QatjRL7WCyG24CnCNppaSFkkr+BZM0T1KHpI6NGzdWIgwbBF+gaTbyVCvJbIyIJRHx+4joKryGe1BJOwM/AT4bEVuAq4B9ye6LtoHX7gL9OhGxIN2ss7WlpWW4YdggFS7QHC18gabZCFG162TSvcuWkZ38ByAifjrUA0oaS5ZgFheOExFPF63/DnDzkCO2ivMFmmYjT7WSzJnA/sBYoCeVBTCkJCNJwHeBRyLiG0Xle6TzNQAnAquGHLHlwhdomo0s1UoyB0XE31TweO8ATgcelLQilX0JOCXdsiaAx4FPVrBOy4lnnJk1r2olmXZJB0TEw5U4WET8kuwJm70trcTxrXo848ysuVXrxP87gRWSVqeZXw9KWlmluq2OecaZWXOrVk9mbpXqsQbjRwKYNbdqXfE/7OnK1pw848ysuVWrJ2PWp1IzzjwZwKw5OMlY3fFkALPmUa0T/2aD5skAZs3DScbqjm8/Y9Y8PFxmdaevyQA+T2PWeJxkrC71ngzg8zRmjcnDZdYQfJ7GrDE5yVhD6Os8TWfXZq5cvobOrs01jtDMSvFwmTWEUudp+hpC87kbs/rhJGMNo/d5mr6G0Hzuxqx+eLjMGlapIbT+zt14aM2s+tyTsYbV11TnUjfc9NCaWW04yVhD6z2E1lfiKWdozYnHrHKaLslImgv8MzAauDoivlzjkKzKSt1ws9QjBYabeJyMzAY2ev78+bWOoWIkjQZ+DrwbuAy4/OKLL75r/vz5G/vaZ8GCBfPnzZtXrRCtRvacsANv33cye++6I+ce9RbmTJ3I6FHiphVPERGMHTOKc496C+1rN/HLR5+lJyAi2HvXHRk9Spx6dTu/fPRZblrxFG/fdzIbnn/5DWV7TtiBzq7N3Hj/U4weJfacsMOr9Zcqr0bZSK273uJplLrLcfHFF2+YP3/+goG2a7aezCHAmohYCyDpBuB4oCKPfbbGNtihtcH2eMrpBfUuL7VtpctGat31Fk+j1J1Xb1wRkcuBa0HSScDciDgrfT4dODQizum13Tyg0H3ZD1g9xConA88Ocd9600xtgWG0R+N22GnUuB3H92x98YXY+tKfNW6HncZO3PMtgIDYtnn97wB6l40at+P40TtN3CsriXjlz39c/8qfNv1h9M6T3ty7HCDvsvqte8JeSHUUj+subEt5/26mRkTLQBs1W09GJcrekEUjYgEwYDdvwMqkjohoHe5x6kEztQXcnnrWTG2B5mpPHm1ptutk1gF7F32eAqyvUSxmZiNesyWZe4GZkqZLGgecDCypcUxmZiNWUw2XRUS3pHOAW8mmMC+MiIdyrHLYQ251pJnaAm5PPWumtkBztafibWmqE/9mZlZfmm24zMzM6oiTjJmZ5cZJxszMcuMkY2ZmuXGSMTOz3DjJmJlZbpxkzMwsN04yZmaWGycZMzPLjZOMmZnlxknGzMxy4yRjZma5cZIxM7PcOMmYmVlumup5MkMxefLkmDZtWq3DMDNrKJ2dnc9GRMtA2+WeZCQtBN4LPBMRf11ivYB/Bo4DXgT+a0Tcl9adAfxj2vR/RcSiVD4HuBbYAVgKnBsRIWlX4AfANOBx4MMRsbm/+KZNm0ZHR8cwW2lmNrJI6hrMdtUYLrsWmNvP+mOBmek1D7gKICWMi4BDgUOAiyRNTPtclbYt7Fc4/nnAsoiYCSxLn0eEzq7NXLl8DZ1dmwcsr0ZZLeuut3hct38DjVJ3HnLvyUTEXZKm9bPJ8cB1kT2is13SBEl7AIcDt0fEcwCSbgfmSroD2CUi7knl1wEnALekYx2ejrsIuAP4YmVbVL7Ors20r91E24xJzJk6seJlnV2bOfXqdrZ29zBuzCgWn9XWZzmQe1kt6663eFy3fwONUnfh70ul1cM5mb2AJ4s+r0tl/ZWvK1EOsHtEbACIiA2SditVoaR5ZD0h9tlnnwo0IdNXQsj7h9S+dhNbu3voCdjW3UP72k3MmTqxZDmQe1kt6663eFy3fwONUnczJxmVKIshlA9aRCwAFgC0traWtW9f+upNVOPH1TZjEuPGjGJbdw9jx4yibcYkgD7Lq1FWy7rrLR7X7d9Ao9SdB2WjVPlKw2U393Hi/1+BOyLi+vR5NdmQ1+HA4RHxyeLt0mt5ROyfyk8pbFfYN/Vi9kjH3a+/2FpbW6MSJ/6vXL6Gr9+2mp6A0YLPH7MfZx/xV68mn8KXWdwbqVRZoUvcuxcF+Q/V9VVWy7rrLR7X7d9Ao9RdDkmdEdE64HZ1kGTeA5xDNrvsUODyiDgknfjvBGanTe8D5kTEc5LuBT4N/JpsdtkVEbFU0j8BmyLiy5LOA3aNiP/RX2yVSjKlkkm1f0hmZtVSN0lG0vVkvZLJwNNkM8bGAkTEt9MU5m+RzRB7ETgzIjrSvh8DvpQOdWlEXJPKW3ltCvMtwKfTFOZJwA+BfYAngA8VJg70pVJJBvzH38xGjrpJMvWukknGzGykGGyS8W1lzMwsN04yZmaWGycZMzPLjZOMmZnlxknGzMxy4yRjZma5GVSSUWbvvIMxM7PmMqgkk+6Q/LOcYzEzsyZTznBZu6S35RaJmZk1nXLuwnwE8Mn0NLQ/k90NOSLiwFwiMzOzhldOkjk2tyjMzKwpDTrJREQXQHoQ2Pa5RWRmZk1j0OdkJL1f0qPA74E7gcfJ7oBsZmZWUjkn/v8n0Ab8LiKmA0cB/5FLVGZm1hTKSTLbImITMErSqIhYDszKKS4zM2sC5Zz4/6OknYG7gMWSngG68wnLzMyaQTk9mePJnlz5OeDnwGPA+/IIyszMmkM5PZl5wI8iYh2wKKd4zMysiZTTk9kFuFXS3ZLOlrR7XkGZmVlzGHSSiYiLI+KtwNnAnsCdkn4x0H6S5kpaLWmNpPNKrJ8qaZmklZLukDSlaN1XJK1Kr48Uld8taUV6rZf0s1R+uKTni9ZdONj2mZlZ5ZUzXFbwDPAHYBOwW38bShoNXAkcDawD7pW0JCIeLtrsa8B1EbFI0pHAZcDpkt4DzCabwbYdWVK7JSK2RMR/LqrjJ8BNRce7OyLeO4R2mZlZhZVzMeanJN0BLAMmA58YxH3LDgHWRMTaiNgK3EA2gaDYAemYAMuL1h8A3BkR3RHxZ+ABYG6vmMYDR+I7RJuZ1aVyzslMBT4bEW+NiIt69UaQNLHEPnsBTxZ9XpfKij0AfDAtnwiMlzQplR8raUdJk8lu0Nn7mTYnAssiYktR2dslPSDpFklvLaN9ZmZWYeWckzkvIlb0s8myEmUqdahen78AHCbpfuAw4CmgOyJuA5YCvwKuB+7hjdflnJLWFdwHTI2Ig4Ar6KOHI2mepA5JHRs3buynSWZmNhyVfPxyqYSyjtf3PqYA64s3iIj1EfGBiDgYuCCVPZ/eL42IWRFxdDr+o69WlvV2DgH+b9GxtkTEn9LyUmBs6gW9TkQsiIjWiGhtaWkZWmvNzGxAlUwyvXsoAPcCMyVNlzQOOBlYUryBpMmSCnGcDyxM5aNTIkHSgcCBwG1Fu34IuDkiXi461pslKS0fQta+TZVonJmZlW8os8sGLSK6JZ0D3AqMBhZGxEOSLgE6ImIJcDhwmaQgu2XN2Wn3scDdKWdsAU6LiOLhspOBL/eq8iTgU5K6gZeAk9Ojo83MrAZUqb/Bku5PQ14NpbW1NTo6OmodhplZQ5HUGRGtA203YE9G0q79rY+I59LiUYOMzczMRojBDJd1kp1vEbAPsDktTwCeAKbD65KNmZkZMIgT/xExPSJmkJ1XeV9ETI6IScB7gZ/mHaCZmTWucmaXvS1NCwYgIm4hu67FzMyspHJmlz0r6R+B75ENn52GpwebmVk/yunJnAK0ADemV0sqMzMzK2nQPZl0Yv9cSTsXrqo3MzPrTzl3Yf5bSQ8DD6fPB0n6l9wiMzOzhlfOcNk3gXeTzsNExAPA3+URlJmZNYey7l0WEU/2KnqlgrGYmVmTKWd22ZOS/haIdLPLzwCP5BOWmZk1g3J6Mv+N7OaVe5Hdwn8Wr93M0szM7A3KmV32LHBqjrGYmVmTGXSSkdQCfAKYVrxfRHys8mGZmVkzKOeczE3A3cAv8Al/MzMbhHKSzI4R8cXcIjEzs6ZTzon/myUdl1skZmbWdMpJMueSJZqXJG2R9IKkLXkFZmZmja+c2WXj8wzEzMyaz4A9GUn7p/fZpV6D2H+upNWS1kg6r8T6qZKWSVop6Q5JU4rWfUXSqvT6SFH5tZJ+L2lFes1K5ZJ0eapr5WDiMzOz/AymJ/N5YB7w9RLrAjiyrx0ljQauBI4mu4DzXklLIuLhos2+BlwXEYskHQlcBpwu6T3AbLKLPrcD7pR0S0QUhuj+ISJ+3KvKY4GZ6XUocFV6NzOzGhgwyUTEvPR+xBCOfwiwJiLWAki6ATiedCfn5ADgc2l5OfCzovI7I6Ib6Jb0ADAX+GE/9R1PlrACaJc0QdIeEbFhCLGbmdkwlXWDTEl/LenDkj5aeA2wy15A8U0116WyYg8AH0zLJwLjJU1K5cdK2lHSZOAIYO+i/S5NQ2LflLRdGfWZmVmVlPM8mYuAK9LrCOCrwPsH2q1EWfT6/AXgMEn3A4cBTwHdEXEbsBT4FXA9cA/QnfY5H9gfeBuwK1C4fmcw9SFpnqQOSR0bN24coAlmZjZU5fRkTgKOAv4QEWcCB5GdK+nPOl7f+5gCrC/eICLWR8QHIuJg4IJU9nx6vzQiZkXE0WQJ5NFUviEyfwGuIRuWG1R9af8FEdEaEa0tLS2DaLqZmQ1FOUnmpYjoITs/sgvwDDBjgH3uBWZKmp4eD3AysKR4A0mTJRXiOB9YmMpHp2EzJB0IHAjclj7vkd4FnACsSvsvAT6aZpm1Ac/7fIyZWe2Uc1uZDkkTgO8AncCfgN/0t0NEdEs6B7gVGA0sjIiHJF0CdETEEuBw4DJJAdzFa48PGAvcneURtgCnpUkAAIvTDTsFrCB7DAFkw2vHAWuAF4Ezy2ifmZlVmLKJWGXuJE0DdomIlZUOqNpaW1ujo6Oj1mGYmTUUSZ0R0TrQdgP2ZPq7oFHS7Ii4r9zgzMxsZBjMcFmpizAL+r0Y08zMRrbBXIw5lIswzczMynoy5vbA3wPvJOvB3A18OyJezik2MzNrcOXMLrsOeIHsYkyAU4B/Az5U6aDMzKw5lJNk9ouIg4o+L0/3EzMzMyupnIsx708XOAIg6VDgPyofkpmZNYtyejKHkl1N/0T6vA/wiKQHgYiIAysenZmZNbRykszc3KIwM7OmVE6SmRkRvygukHRGRCyqcExmZtYkyjknc6GkqyTtJGl3Sf8OvC+vwMzMrPGVk2QOAx4juyHlL4HvR8RJuURlZmZNoZwkM5Hs5P9jwF+AqelW+2ZmZiWVk2TagVsiYi7ZEyn3xFOYzcysH+UkmXcB2yRdGBEvAV8DzssnLDMzawblJJnzgTay28lAdouZ/u7QbGZmI1xZF2NGxGxJ9wNExOb0SGUzM7OSyunJbJM0muwOzKTHH/fkEpWZmTWFcpLM5cCNwG6SLiWbxvy/c4nKzMyawqCHyyJisaRO4ChAwAkR8UhukZmZWcMrpydDRPw2Iq6MiG8NNsFImitptaQ1kt4wG03SVEnLJK2UdIekKUXrviJpVXp9pKh8cTrmKkkLJY1N5YdLel7SivS6sJz2mZlZZZWVZMqVzuFcCRwLHACcIumAXpt9Dbgu3cX5EuCytO97gNnALLKLQP9B0i5pn8XA/sDfADsAZxUd7+6ImJVel+TTMjMzG4xckwxwCLAmItZGxFbgBuD4XtscACxLy8uL1h8A3BkR3RHxZ+AB0p2gI2JpJMBvgCmYmVndyTvJ7AU8WfR5XSor9gDwwbR8IjBe0qRUfqykHSVNBo4A9i7eMQ2TnQ78vKj47ZIekHSLpLeWCkrSPEkdkjo2btw41LaZmdkA8k4ype5tFr0+fwE4LF1/cxjwFNAdEbcBS4FfAdcD9wDdvfb9F+CuiLg7fb4PmJoeE30F8LNSQUXEgohojYjWlpaWITTLzMwGI+8ks47X9z6mAOuLN4iI9RHxgYg4GLgglT2f3i9N51aOJktYjxb2k3QR0AJ8vuhYWyLiT2l5KTA29YLMzKwG8k4y9wIzJU1Pdwc4GVhSvIGkyZIKcZwPLEzlo9OwGZIOBA4EbkufzwLeDZwSET1Fx3pz4c7Qkg4ha9+mPBrW2bWZK5evobNrcx6HNzNrCuXcVqZsEdEt6RzgVmA0sDAiHpJ0CdAREUuAw4HLJAVwF3B22n0scHfKGVuA0yKiMFz2baALuCet/2maSXYS8ClJ3cBLwMlpckBFdXZt5tSr29na3cO4MaNYfFYbc6ZOrHQ1ZmYNL9ckA68OWy3tVXZh0fKPgR+X2O9lshlmpY5ZMu6I+BbwreHEOxjtazextbuHnoBt3T20r93kJGNmVkLew2VNqW3GJMaNGcVowdgxo2ibManWIZmZ1aXcezLNaM7UiSw+q432tZtomzHJvRgzsz44yQzRnKkTnVzMzAbg4TIzM8uNcph81VAkbSSbqTYUk4FnKxhOLTVTW8DtqWfN1BZorvaU05apETHg1ewjPskMh6SOiGitdRyV0ExtAbennjVTW6C52pNHWzxcZmZmuXGSMTOz3DjJDM+CWgdQQc3UFnB76lkztQWaqz0Vb4vPyZiZWW7ckzEzs9w4yZiZWW6cZIZA0lxJqyWtkXRereMZLEmPS3pQ0gpJHalsV0m3S3o0vU9M5ZJ0eWrjSkmzaxs9SFoo6RlJq4rKyo5f0hlp+0clnVFHbZkv6an0/ayQdFzRuvNTW1ZLendRec1/i5L2lrRc0iOSHpJ0bipv1O+mr/Y06vezvaTfKHti8EOSLk7l0yX9Ov23/oGyx7Egabv0eU1aP63oWCXb2a+I8KuMF9kjCx4DZgDjyB4TfUCt4xpk7I8Dk3uVfRU4Ly2fB3wlLR8H3EL2sLg24Nd1EP/fAbOBVUONH9gVWJveJ6bliXXSlvnAF0pse0D6nW0HTE+/v9H18lsE9gBmp+XxwO9SzI363fTVnkb9fgTsnJbHAr9O/91/SPY4FMgen/KptPz3wLfT8snAD/pr50D1uydTvkOANRGxNiK2AjcAx9c4puE4HliUlhcBJxSVXxeZdmCCpD1qEWBBRNwFPNeruNz43w3cHhHPRcRm4HZgbv7Rv14fbenL8cANEfGXiPg9sIbsd1gXv8WI2BAR96XlF4BHgL1o3O+mr/b0pd6/n4j0xGCyJDMWCOBIXnvMSu/vp/C9/Rg4SpLou539cpIp317Ak0Wf19H/D7CeBHCbpE5J81LZ7hGxAbJ/XMBuqbxR2llu/PXernPSENLCwvASDdSWNLRyMNn/LTf8d9OrPdCg34+yJw2vAJ4hS96PAX+M1x4EWRzbq3Gn9c8Dkxhie5xkyqcSZY0yD/wdETEbOBY4W9Lf9bNtI7cT+o6/ntt1FbAvMAvYAHw9lTdEWyTtDPwE+GxEbOlv0xJljdCehv1+IuKViJgFTCHrffynUpul94q2x0mmfOuAvYs+TwHW1yiWskTE+vT+DHAj2Y/t6cIwWHp/Jm3eKO0sN/66bVdEPJ3+GPQA3+G1oYi6b4uksWR/kBdHxE9TccN+N6Xa08jfT0FE/BG4g+yczARJhce9FMf2atxp/ZvIhnaH1B4nmfLdC8xMMzPGkZ0YW1LjmAYkaSdJ4wvLwDHAKrLYC7N4zgBuSstLgI+mmUBtwPOFoY86U278twLHSJqYhjuOSWU11+uc14lk3w9kbTk5zfqZDswEfkOd/BbTeP13gUci4htFqxryu+mrPQ38/bRImpCWdwDeRXaeaTlwUtqs9/dT+N5OAv5fZGf++2pn/6o906EZXmSzY35HNq55Qa3jGWTMM8hmhjwAPFSIm2ywuBfBAAAAq0lEQVSsdRnwaHrfNZULuDK18UGgtQ7acD3ZMMU2sv+r+vhQ4gc+RnbScg1wZh215d9SrCvTP+g9ira/ILVlNXBsPf0WgXeSDZusBFak13EN/N301Z5G/X4OBO5Pca8CLkzlM8iSxBrgR8B2qXz79HlNWj9joHb29/JtZczMLDceLjMzs9w4yZiZWW6cZMzMLDdOMmZmlhsnGTMzy42TjJmZ5cZJxszMcvP/ARaF2SgN/neNAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ]
  },
  {
   "metadata": {
    "trusted": true
   },
   "cell_type": "code",
   "source": "df.to_csv('diff_var.csv', index=False)",
   "execution_count": 116,
   "outputs": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python",
   "version": "3.6.7",
   "mimetype": "text/x-python",
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "pygments_lexer": "ipython3",
   "nbconvert_exporter": "python",
   "file_extension": ".py"
  },
  "gist_id": "ff2e9b19f488199b919c1736258f8133"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}