From 811ffca704543908cdda3e66ec77843f41fa4b0e Mon Sep 17 00:00:00 2001
From: Theo Guyard <theo.guyard@insa-rennes.fr>
Date: Fri, 21 Mar 2025 10:50:21 -0400
Subject: [PATCH] Move repo
---
.gitignore | 15 --
LICENSE | 21 --
MANIFEST.in | 13 --
README.md | 87 +-------
exp/ICASSP/convergence.py | 86 --------
exp/ICASSP/performance_profiles.py | 118 -----------
exp/utils.py | 22 --
notebooks/examples.ipynb | 260 ------------------------
pyproject.toml | 3 -
requirements.txt | 4 -
setup.cfg | 2 -
setup.py | 32 ---
srenet/__init__.py | 1 -
srenet/algebra.py | 95 ---------
srenet/data.py | 88 --------
srenet/enet.py | 19 --
srenet/projgrad.py | 316 -----------------------------
srenet/utils.py | 130 ------------
tests/convergence.py | 64 ------
tests/identification.py | 48 -----
20 files changed, 1 insertion(+), 1423 deletions(-)
delete mode 100644 .gitignore
delete mode 100644 LICENSE
delete mode 100644 MANIFEST.in
delete mode 100644 exp/ICASSP/convergence.py
delete mode 100644 exp/ICASSP/performance_profiles.py
delete mode 100644 exp/utils.py
delete mode 100644 notebooks/examples.ipynb
delete mode 100644 pyproject.toml
delete mode 100644 requirements.txt
delete mode 100644 setup.cfg
delete mode 100644 setup.py
delete mode 100644 srenet/__init__.py
delete mode 100644 srenet/algebra.py
delete mode 100644 srenet/data.py
delete mode 100644 srenet/enet.py
delete mode 100644 srenet/projgrad.py
delete mode 100644 srenet/utils.py
delete mode 100644 tests/convergence.py
delete mode 100644 tests/identification.py
diff --git a/.gitignore b/.gitignore
deleted file mode 100644
index fb38690..0000000
--- a/.gitignore
+++ /dev/null
@@ -1,15 +0,0 @@
-.DS_Store
-.vscode
-.ipynb_checkpoints/
-
-# Python .gitignore
-build/
-dist/
-*.egg-info/
-*.egg
-*.py[cod]
-__pycache__/
-*.so
-*~
-.tox
-.cache
\ No newline at end of file
diff --git a/LICENSE b/LICENSE
deleted file mode 100644
index 0b6a1b6..0000000
--- a/LICENSE
+++ /dev/null
@@ -1,21 +0,0 @@
-MIT License
-
-Copyright (c) 2021 Guyard Theo
-
-Permission is hereby granted, free of charge, to any person obtaining a copy
-of this software and associated documentation files (the "Software"), to deal
-in the Software without restriction, including without limitation the rights
-to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-copies of the Software, and to permit persons to whom the Software is
-furnished to do so, subject to the following conditions:
-
-The above copyright notice and this permission notice shall be included in all
-copies or substantial portions of the Software.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-SOFTWARE.
diff --git a/MANIFEST.in b/MANIFEST.in
deleted file mode 100644
index 4ec2135..0000000
--- a/MANIFEST.in
+++ /dev/null
@@ -1,13 +0,0 @@
-include pyproject.toml
-
-# Include the README
-include *.md
-
-# Include the license file
-include LICENSE
-
-# Include setup.py
-include setup.py
-
-# Include the data files
-recursive-include data *
\ No newline at end of file
diff --git a/README.md b/README.md
index ca7f897..7a1b81e 100644
--- a/README.md
+++ b/README.md
@@ -1,86 +1 @@
-# Screen & Relax : A fast solving procedure for the Elastic-Net problem by safe identification of the solution support.
-
-This repository implements the Screen & Relax method solving the Elastic-Net problem. We introduced this method in the paper
-
-> Guyard, T., Herzet, C., & Elvira, C. (2022, May). Screen & relax: accelerating the resolution of Elastic-Net by safe identification of the solution support. In *ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)* (pp. 5443-5447). IEEE.
-
-available [here](https://arxiv.org/pdf/2110.07281.pdf). This paper contains theoretical and numerical results that can be reproduced with this toolbox. If you encounter a bug or something unexpected, please let me know by raising an issue on the project page or by contacting me by [mail](mailto:theo.guyard@insa-rennes.fr).
-
-## Requirements
-
-The `srenet` (**S**creen & **R**elax method for the **E**lastic-**NET** problem) package works with Python 3.9+. Please refer to the [Python download page](https://www.python.org/downloads/) for install instructions.
-
-**Dependencies** : [NumPy](http://www.numpy.org), [SciPy](https://www.scipy.org), [Matplotlib](http://matplotlib.org), [Pandas](https://pandas.pydata.org)
-
-## Install instructions
-
-1. Clone the repository
-```bash
-git clone https://gitlab.insa-rennes.fr/Theo.Guyard/screen-and-relax.git
-```
-
-2. Enter the project folder
-```bash
-cd screen-and-relax
-```
-
-3. Ensure that Python 3.9+ is installed
-```bash
-python -V
-> Python 3.9.1 (or higher)
-```
-
-4. Install external dependencies
-```bash
-pip install --upgrade pip
-pip install -r requirements.txt
-```
-
-5. Install the `srenet` package
-```bash
-pip install .
-```
-
-## ICASSP experiments
-
-Experiments presented in the Screen & Relax paper submitted to [ICASSP 2022](https://2022.ieeeicassp.org) are located in the folder `screen-and-relax/exp/ICASSP/`. The two available experiments are
-* `convergence.py` : convergence of the four compared methods on a single problem instance
-* `performance_profiles.py` : figure 1 of the "Screen & Relax" paper
-
-To run experiments, make sure to be located at the project root and run experiments as follows.
-```bash
-cd screen-and-relax
-python exp/ICASSP/convergence.py
-python exp/ICASSP/performance_profiles.py
-```
-
-You can also edit the experimental setup directly in the script, just follow the comments documenting the code.
-
-## Running demo examples
-
-We also provide a notebook explaining how to use the `srenet` package. You can use it as following :
-```bash
-cd notebooks
-jupyter notebook
-```
-You may have to install Jupyter. Please refer to the [Jupyter website](https://jupyter.org) for install instructions.
-
-## Licence
-
-This software is distributed under the [MIT Licence](https://mit-license.org).
-
-
-## Cite this work
-
-If you use this package for your own work, please consider citing it as :
-
-```bibtex
-@inproceedings{guyard2022screen,
- title={Screen \& relax: accelerating the resolution of Elastic-Net by safe identification of the solution support},
- author={Guyard, Th{\'e}o and Herzet, C{\'e}dric and Elvira, Cl{\'e}ment},
- booktitle={ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)},
- pages={5443--5447},
- year={2022},
- organization={IEEE}
-}
-```
+This repository is deprecated. Please refer to the new repository [here](https://github.com/TheoGuyard/code_2022_ICASSP_ScreenAndRelax).
diff --git a/exp/ICASSP/convergence.py b/exp/ICASSP/convergence.py
deleted file mode 100644
index f8f4cf1..0000000
--- a/exp/ICASSP/convergence.py
+++ /dev/null
@@ -1,86 +0,0 @@
-# ICASSP experient - Convergence
-
-import os, sys
-sys.path.append(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
-
-import warnings
-import numpy as np
-import pandas as pd
-import matplotlib.pyplot as plt
-from srenet.data import sample_data
-from srenet.projgrad import projgrad
-from utils import get_params_from_string, get_label_from_string, get_color_from_string
-
-
-################################ Configuration #################################
-
-# This is the only part of the code you may have to edit
-
-m, n = 100, 300 # Dictionnary dimension
-dic_gen = 'gaussian' # Dictionnary generation method : `gaussian`, `uniform`, `dct` or `toeplitz`
-obs_gen = 'gaussian' # Observation generation method : `gaussian` or `posgaussian`
-l1ratio = 0.2 # Ratio lambda/lambda_max where lambda_max = max(A.T * y)
-l2ratio = 0.5 # Ration epsilon/lambda_max where lambda_max = max(A.T * y)
-
-
-################################## Experiment ##################################
-
-print("-------------------------------")
-print("ICASSP experiment - Convergence\n")
-
-methods = ['aPG','aPGs','aPGr','S&R']
-res = {method: {} for method in methods}
-
-x0 = np.zeros(n)
-A, y = sample_data(m, n, dic_gen, obs_gen)
-lmax = np.max(A.T @ y)
-l1 = l1ratio * lmax
-l2 = l2ratio * lmax
-
-for method in methods:
- print(f"Method {method} is running ...")
- params = get_params_from_string(method)
- params.maxit = 10_000
- params.tol = 1e-15
- params.verbosity = False
- x, mnt = projgrad(x0,A,y,l1,l2,params)
- res[method]['x'] = x
- res[method]['mnt'] = mnt
-
-
-############################## Plot results ################################
-
-plt.rcParams["font.size"] = 12
-plt.rcParams["axes.labelsize"] = 14
-plt.rcParams["axes.titlesize"] = 14
-plt.rcParams["lines.linewidth"] = 2
-plt.rcParams["lines.markersize"] = 5
-plt.rcParams["xtick.labelsize"] = 12
-plt.rcParams["ytick.labelsize"] = 12
-plt.rcParams["text.usetex"] = True
-plt.rcParams["font.family"] = "serif"
-plt.rcParams["font.serif"] = ["Helvetica"]
-
-fig, ax = plt.subplots(1, 1, figsize=(6,4))
-
-for method in methods:
- ax.plot(
- res[method]['mnt'].flops,
- res[method]['mnt'].obj - res[method]['mnt'].obj[-1],
- label=get_label_from_string(method),
- color=get_color_from_string(method),
- marker= '.',
- )
-ax.set_xlabel(r"FLOPs")
-ax.set_ylabel(r"Objective error")
-ax.set_xscale('log')
-ax.set_yscale('log')
-
-plt.legend()
-plt.grid(b=True, which='major', axis='y')
-plt.grid(b=True, which='minor', axis='y', alpha=0.2)
-plt.minorticks_on()
-plt.tight_layout()
-plt.show()
-
-print("-------------------------------")
\ No newline at end of file
diff --git a/exp/ICASSP/performance_profiles.py b/exp/ICASSP/performance_profiles.py
deleted file mode 100644
index f4dcaa8..0000000
--- a/exp/ICASSP/performance_profiles.py
+++ /dev/null
@@ -1,118 +0,0 @@
-# ICASSP experient - Performance profiles
-
-import os, sys
-sys.path.append(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
-
-import warnings
-import numpy as np
-import pandas as pd
-import matplotlib.pyplot as plt
-from srenet.data import sample_data
-from srenet.projgrad import projgrad
-from utils import get_params_from_string, get_label_from_string, get_color_from_string
-
-
-################################ Configuration #################################
-
-# This is the only part of the code you may have to edit
-
-m, n = 100, 300 # Dictionnary dimension
-dic_gen = 'gaussian' # Dictionnary generation method : `gaussian`, `uniform`, `dct` or `toeplitz`
-obs_gen = 'gaussian' # Observation generation method : `gaussian` or `posgaussian`
-l1ratio = 0.2 # Ratio lambda/lambda_max where lambda_max = max(A.T * y)
-l2ratio = 0.5 # Ration epsilon/lambda_max where lambda_max = max(A.T * y)
-flops = 2_000_000 # FLOPs budget
-repeats = 100 # Number of instances tested
-
-
-################################## Experiment ##################################
-
-print("---------------------------------------")
-print("ICASSP experiment - Performance profiles\n")
-
-tolmax = 1.
-tolmin = 1e-15
-tolnb = 30
-tols = np.logspace(np.log10(tolmax), np.log10(tolmin), tolnb)
-methods = ['aPG','aPGs','aPGr','S&R']
-solved = np.full((len(methods), tolnb, repeats), 0)
-
-for rep in range(repeats):
- print(f"Repeat {rep+1}/{repeats}")
-
- x0 = np.zeros(n)
- A, y = sample_data(m, n, dic_gen, obs_gen)
- lmax = np.max(A.T @ y)
- l1 = l1ratio * lmax
- l2 = l2ratio * lmax
-
- for (mi,method) in enumerate(methods):
- try:
- print(f" Method {method} is running ...")
-
- params = get_params_from_string(method)
- params.maxit = 10_000
- params.tol = tolmin
- params.verbosity = False
-
- x, mnt = projgrad(x0,A,y,l1,l2,params)
-
- if not mnt.stopcrit.has_converged:
- warnings.warn("Method did not converge.")
- continue
-
- for ti, tol in enumerate(tols):
- its_solved = np.where(mnt.gap < tol)[0]
- if len(its_solved) > 0:
- it_solved = its_solved[0]
- if mnt.flops[it_solved] < flops:
- solved[mi,ti,rep] = 1
- else:
- solved[mi,ti,rep] = 0
-
- except Exception as e:
- warnings.warn("Repeat failed with the following error :")
- print(e)
-
-solved = pd.DataFrame(
- np.nanmean(solved, axis=-1),
- columns = ['{:.2e}'.format(tol) for tol in tols],
- index = [method for method in methods],
- )
-
-############################## Plot results ################################
-
-plt.rcParams["font.size"] = 12
-plt.rcParams["axes.labelsize"] = 18
-plt.rcParams["axes.titlesize"] = 18
-plt.rcParams["lines.linewidth"] = 2
-plt.rcParams["lines.markersize"] = 5
-plt.rcParams["xtick.labelsize"] = 12
-plt.rcParams["ytick.labelsize"] = 12
-plt.rcParams["text.usetex"] = True
-plt.rcParams["font.family"] = "serif"
-plt.rcParams["font.serif"] = ["Helvetica"]
-
-fig, ax = plt.subplots(1, 1, figsize=(6,4))
-
-for method in methods:
- ax.plot(
- tols,
- solved.loc[method,:],
- label=get_label_from_string(method),
- color=get_color_from_string(method),
- marker= '.',
- )
-ax.set_xlabel(r"$\tau$")
-ax.set_ylabel(r"$\rho_{\%}$")
-ax.set_xscale('log')
-ax.yaxis.set_major_formatter(lambda x,pos: r"{0:.0f}\%".format(x * 100))
-
-plt.legend()
-plt.grid(b=True, which='major', axis='y')
-plt.grid(b=True, which='minor', axis='y', alpha=0.2)
-plt.minorticks_on()
-plt.tight_layout()
-plt.show()
-
-print("---------------------------------------")
\ No newline at end of file
diff --git a/exp/utils.py b/exp/utils.py
deleted file mode 100644
index 01fda60..0000000
--- a/exp/utils.py
+++ /dev/null
@@ -1,22 +0,0 @@
-from srenet.utils import ProjgradParams
-
-def get_params_from_string(name):
- if name == 'aPG' : return ProjgradParams(screen=False, relax=False)
- elif name == 'aPGs': return ProjgradParams(screen=True , relax=False)
- elif name == 'aPGr': return ProjgradParams(screen=False, relax=True)
- elif name == 'S&R' : return ProjgradParams(screen=True , relax=True)
- else: raise ValueError("Method name not understood.")
-
-def get_label_from_string(name):
- if name == 'aPG' : return r"aPG"
- elif name == 'aPGs': return r"aPGs"
- elif name == 'aPGr': return r"aPGr"
- elif name == 'S&R' : return r"S\&R"
- else: raise ValueError("Method name not understood.")
-
-def get_color_from_string(name):
- if name == 'aPG' : return "royalblue"
- elif name == 'aPGs': return "navy"
- elif name == 'aPGr': return "darkorange"
- elif name == 'S&R' : return "firebrick"
- else: raise ValueError("Method name not understood.")
\ No newline at end of file
diff --git a/notebooks/examples.ipynb b/notebooks/examples.ipynb
deleted file mode 100644
index 5eb6b07..0000000
--- a/notebooks/examples.ipynb
+++ /dev/null
@@ -1,260 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "source": [
- "# Screen & Relax : A fast solving procedure for the Elastic-Net problem by safe identification of the solution support.\n",
- "## Example notebook\n",
- "---"
- ],
- "metadata": {}
- },
- {
- "cell_type": "markdown",
- "source": [
- "First of all, import dependencies"
- ],
- "metadata": {}
- },
- {
- "cell_type": "code",
- "execution_count": 4,
- "source": [
- "import numpy as np\n",
- "import matplotlib.pyplot as plt\n",
- "from srenet.data import sample_data\n",
- "from srenet.projgrad import projgrad\n",
- "from srenet.utils import ProjgradParams"
- ],
- "outputs": [],
- "metadata": {}
- },
- {
- "cell_type": "markdown",
- "source": [
- "The problem addressed is the non-negative Elastic-Net given by"
- ],
- "metadata": {}
- },
- {
- "cell_type": "markdown",
- "source": [
- "$$min_{\\mathbf{x}\\geq\\mathbf{0}} \\ \\Big\\{ \\ \\tfrac{1}{2}\\|\\mathbf{y}-\\mathbf{A}\\mathbf{x}\\|_2^2 + \\lambda \\|\\mathbf{x}\\|_1 + \\tfrac{\\epsilon}{2}\\|\\mathbf{x}\\|_2^2 \\ \\Big\\}$$"
- ],
- "metadata": {}
- },
- {
- "cell_type": "markdown",
- "source": [
- "Parameters of the problem are defined as"
- ],
- "metadata": {}
- },
- {
- "cell_type": "code",
- "execution_count": 6,
- "source": [
- "m = 100 # y size\n",
- "n = 300 # x size\n",
- "gen_dic = 'gaussian' # Generation method for A\n",
- "gen_obs = 'gaussian' # Generation method for y\n",
- "l1ratio = 0.5 # Ratio lambda / np.max(A.T * y)\n",
- "l2ratio = 0.2 # Ratio epsilon / np.max(A.T * y)"
- ],
- "outputs": [],
- "metadata": {}
- },
- {
- "cell_type": "markdown",
- "source": [
- "Data of the problem can be sampled with the following lines of code"
- ],
- "metadata": {}
- },
- {
- "cell_type": "code",
- "execution_count": 7,
- "source": [
- "x0 = np.zeros(n) # Initial point for the projected gradient\n",
- "A, y = sample_data(m, n, gen_dic, gen_obs) # Dictionary and observation\n",
- "l1 = l1ratio * np.max(A.T @ y) # Value of \\lambda\n",
- "l2 = l2ratio * np.max(A.T @ y) # Value of \\epsilon "
- ],
- "outputs": [],
- "metadata": {}
- },
- {
- "cell_type": "markdown",
- "source": [
- "Now, run a Projected Gradient method to solve the problem :\n",
- "* Without screening and relaxing tests\n",
- "* With screning tests only\n",
- "* With relaxing tests only\n",
- "* With both screening and relaxing tests (Screen & Relax method)"
- ],
- "metadata": {}
- },
- {
- "cell_type": "code",
- "execution_count": 8,
- "source": [
- "x, mnt = projgrad(x0, A, y, l1, l2, ProjgradParams(verbosity=False,screen=False,relax=False))\n",
- "x_s, mnt_s = projgrad(x0, A, y, l1, l2, ProjgradParams(verbosity=False,screen=True,relax=False))\n",
- "x_r, mnt_r = projgrad(x0, A, y, l1, l2, ProjgradParams(verbosity=False,screen=False,relax=True))\n",
- "x_sr, mnt_sr = projgrad(x0, A, y, l1, l2, ProjgradParams(verbosity=False,screen=True,relax=True))\n",
- "res = {\n",
- " 'aPG' : {'x': x, 'mnt': mnt},\n",
- " 'aPGs': {'x': x_s, 'mnt': mnt_s},\n",
- " 'aPGr': {'x': x_r, 'mnt': mnt_r},\n",
- " 'S&R' : {'x': x_sr, 'mnt': mnt_sr},\n",
- "}\n",
- "# x : optimizer of the problem\n",
- "# mnt : monitoring values"
- ],
- "outputs": [],
- "metadata": {}
- },
- {
- "cell_type": "markdown",
- "source": [
- "You can display some monitoring values that are outputted by the algorithm"
- ],
- "metadata": {}
- },
- {
- "cell_type": "code",
- "execution_count": 9,
- "source": [
- "print(\"Monitoring values\")\n",
- "print(\"-----------------\")\n",
- "for k, v in res.items():\n",
- " print(k)\n",
- " print(\" Stop crit : {}\".format(v['mnt'].stopcrit))\n",
- " print(\" Last Gap : {:.2e}\".format(v['mnt'].gap[-1]))\n",
- " print(\" Last Obj : {:.2e}\".format(v['mnt'].obj[-1]))\n",
- " print(\" Iters : {}\".format(len(v['mnt'].gap)))\n",
- " print(\" FLOPs : {:.2e}\".format(v['mnt'].flops[-1]))"
- ],
- "outputs": [
- {
- "output_type": "stream",
- "name": "stdout",
- "text": [
- "Monitoring values\n",
- "-----------------\n",
- "aPG\n",
- " Stop crit : StopCrit.CONVERGED\n",
- " Last Gap : 2.78e-16\n",
- " Last Obj : 4.68e-01\n",
- " Iters : 199\n",
- " FLOPs : 2.48e+07\n",
- "aPGs\n",
- " Stop crit : StopCrit.CONVERGED\n",
- " Last Gap : 2.78e-16\n",
- " Last Obj : 4.68e-01\n",
- " Iters : 199\n",
- " FLOPs : 2.25e+06\n",
- "aPGr\n",
- " Stop crit : StopCrit.CONVERGED\n",
- " Last Gap : 5.55e-17\n",
- " Last Obj : 4.68e-01\n",
- " Iters : 34\n",
- " FLOPs : 4.37e+06\n",
- "S&R\n",
- " Stop crit : StopCrit.CONVERGED\n",
- " Last Gap : 1.11e-16\n",
- " Last Obj : 4.68e-01\n",
- " Iters : 34\n",
- " FLOPs : 1.18e+06\n"
- ]
- }
- ],
- "metadata": {}
- },
- {
- "cell_type": "markdown",
- "source": [
- "You can also construct various graphics with monitoring values :"
- ],
- "metadata": {}
- },
- {
- "cell_type": "code",
- "execution_count": 10,
- "source": [
- "fig, axs = plt.subplots(1,3,figsize=(12,3))\n",
- "\n",
- "for k,v in res.items():\n",
- " axs[0].plot(v['mnt'].flops, v['mnt'].obj-v['mnt'].obj[-1], label=k, marker='.')\n",
- " axs[1].plot(v['mnt'].flops, v['mnt'].gap, label=k, marker='.')\n",
- "\n",
- "I_pos = np.sum(np.abs(res['S&R']['x']) > 1e-8)\n",
- "I_null = np.sum(np.abs(res['S&R']['x']) <= 1e-8)\n",
- "axs[2].plot(res['S&R']['mnt'].gap, res['S&R']['mnt'].sc_size/I_null, label=\"% Screened\", marker='.')\n",
- "axs[2].plot(res['S&R']['mnt'].gap, res['S&R']['mnt'].rc_size/I_pos, label=\"% Relaxed\", marker='.')\n",
- "\n",
- " \n",
- "axs[0].set_xlabel(\"FLOPs\")\n",
- "axs[0].set_ylabel(\"Objective error\")\n",
- "axs[0].set_xscale('log')\n",
- "axs[0].set_yscale('log')\n",
- "axs[0].grid(b=True, which='major', axis='y')\n",
- "axs[0].legend()\n",
- "axs[1].set_xlabel(\"FLOPs\")\n",
- "axs[1].set_ylabel(\"Duality Gap\")\n",
- "axs[1].set_xscale('log')\n",
- "axs[1].set_yscale('log')\n",
- "axs[1].grid(b=True, which='major', axis='y')\n",
- "axs[1].legend() \n",
- "axs[2].yaxis.set_major_formatter(lambda x,pos: r\"{0:.0f}%\".format(x * 100))\n",
- "axs[2].invert_xaxis()\n",
- "axs[2].set_xlabel(\"GAP\")\n",
- "axs[2].set_ylabel(\"Indices identified in S&R method\")\n",
- "axs[2].set_xscale('log')\n",
- "axs[2].grid(b=True, which='major', axis='y')\n",
- "axs[2].legend()\n",
- "\n",
- "plt.tight_layout()\n",
- "plt.show()"
- ],
- "outputs": [
- {
- "output_type": "display_data",
- "data": {
- "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAADSCAYAAABXXGGLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAACRFElEQVR4nOydd3iUZdaH7zOTRighEBJ6RzoEEooNoyugrgXRFcG1Y/tWV92GLva+uutaV0HBsivYsaAooiAoNfTeWygJJUBCSJs53x/vTOqkkExNnvu6xpn3edtJJM+85znn/I6oKgaDwWAwGAwGg8FgqD22QBtgMBgMBoPBYDAYDHUF42AZDAaDwWAwGAwGg5cwDpbBYDAYDAaDwWAweAnjYBkMBoPBYDAYDAaDlzAOlsFgMBgMBoPBYDB4CeNgGQwGg8FgMBgMBoOXMA6WwWAwGAwGg8FgMHiJsEAb4A/i4uK0Y8eOgTbDYKiXLF++/LCqtgi0Hf7EzDkGQ2A4nflGRJpVtl9Vj3rHKt9i5huDIXBUNOfUCwerY8eOpKamBtoMg6FeIiK7A22DvzFzjsEQGE5zvlkOKCBAeyDT9bkpsAfo5G37fIGZbwyGwFHRnGNSBA0Gg8FgMNQ7VLWTqnYG5gCXqWqcqjYHLgVmB9Y6g8EQyhgHy2AwGCpBRBqKyHsi8paIXBdoewwGg9cZqqrfujdUdRZwVgDtMRgMIY5xsAwGQ71DRKaKSIaIrCszfpGIbBaRbSLygGt4NPCpqt4GXO53Yw0Gg6/ZLyIPiUhH12sisD/QRhkMhtDFOFgGg6E+8i5wUckBEbEDrwMXA72AsSLSC2gL7HUd5vCjjQaDwT+MBVoAM1yveNdYpXhaqBGRZiLyg4hsdb3HusZFRF5xLd6sEZGBrvHuIrLcNXamayxMROaISLQPflaDweAH6oXIRW5OFptXzqvyuA0/TCd3ywaanTWckbc84nvDDAZDQFDV+SLSsczwYGCbqu4AEJEPgSuANCwnaxWVLEqJyO3A7QAt45vz+f9eq9SG4+vnE7ZvD87uiXQ4+9oa/iQGg6G2uNQC7xWRxtamZlfz1HeB14D3S4w9APyoqs+5ouAPABOwFm66uV5DgDdc73cA9wK7gJeBq4C7gP+pak7tfjIDwLQle/ho2R4iw2x0S2hM79YxZObkM7Rzc5I6xLJ8dyaLdxypcBvwOFYS9/7Y6IhKr12WTcvmkLV2Fk27DaVb4jCf/y7qMoqiaqnWONX92fXu+oy69mGNldq2jmDP2oXk7lhITL+L6TX04hrbUy8cLHbswTn2rioP6+F611+ms/SV6ZyMhrwoIT/SRmGUncKocJwNIqFhQ2yNGhMe05yoZi1o1KItzVp3oVXHnsQ0b4U9rPJf6/dTnyBzyU/EDrnAOHIGQ/DQhuJIFViO1RDgFeA1Efkt8HVFJ6vqZGAyQJ+oBtrzqderdVNd+jPZn/zMyYaQGyXkNbBTGBWGIzoSbRSNrXEMYU2bEx3XiphWHYlr14O2XfvSoGGTmv6cBoOhBCLSF8tJaubaPgzcqKrrKjuvgoWaK4AU1+f3gHlYDtYVwPuqqsBiEWkqIq2AAiDa9SoQkabAZZSJsBuKKem0OJxO5m46RM9WjWnbLJrs3EKycgvJyi0gO6+QhdsP89OmQ0XnLt2VWfTZbhOS2jdl+Z5jOJzqcfs3PeIB+HFTRqmxFo0ji65zKCuvaH9l1y57XuNDK/jzvvsJFwfseRt+9OVvre4jrheAvRbXae56z9s3jU32D+kx6MIaXadeOFg5MRGsvqRzpcc02LSbrjsc2LC832MxkNM4jPA8J9EnHUQdcdAgL5/ovJOA59YY6UCaHXIiITcK8iKF/CgbBZHFzpmcyqPnmpO0dULhL9P5Ji+H3971nLd/ZIPB4CVU9SRw8+mcU9WcU9F8E3nKQUxmIQ1yC2l4Kpcw53HgQKlzC7GWunMiIScKTkUJeQ2EgqgwChuE42zYABo1wh4TS1SzeHL270H37yNm2Eguvu2J0/zpDYZ6wSTgT6o6F0BEUrAWS2oidJGgqu4/2oNAguuzpwWcNlhpye8DkVjRrIeBZ1TVWdlNSkbMExISmDdvXg1MDW62ZTrYdNRBj2Z2nCipBx2cKlR+3e+ghC9TYxxOZfnuTBxa8fbPm9Ndn0uPRZZ4gs9zFO+v7Nplz3tWpxFut7LOHSrM1YEslv6WkyDFzkLRu5TZLjNujUn5Yz1dq8y4e6z4vOKdJR2XCu2Qqu2s/H54vJ+U+E+p8TIXspXerOJapcdKjiccX8mgguXYRbHjYMcvn3HwZM1cpXrhYDVr141rX/ym0mO+n/oEBS9OJ8wBhXbQG8cyykN0Ke9UDul7NpGRtoXjB/eQc3g/eZmHcWQdQ7OzkZxT2E/lE5ZbSHiekwYnHcSWcs6KiXBAx5e/ZOHULzkWI+TERFDQvAlhLdvQtGtfuiRfSIczBlYZETMYDF5hH9CuxHZb19hpU9WcU535xlFYyKH9Ozmwcx2Z+7aTnbGXvKMZOE5kQtZJbDmnCDtVQMSpQqJOKU0z84k+lU90/kngMJYbVowu+4T5b3/Ckfhw8uIaY2vdhmY9B9Jn2Chatu+BwVCPaeh2rgBUdZ6INKztRVVVRaRSV0BV9+CKeIlIV6x5Z6OI/BeIAB5W1S0eziuKmCcnJ2tKSkptzQ0alu/O5LMVaXySupcCh2IF+SrHJnDlgDaMHdyexlHhNIoKo3FUGF+v2s/EL8oHIm0CEWE2Hrm0N0/MXE9BoZNwD9sfjB8KwHVvLy41VjLdb/nuTK57ezH5BU6clVy71Hl7luB8NxWHQ1CEAsJoe9nfeaiG0RKDd9i0bA75M8cSroUUEEbnc66ix6CUGl0rJJ/cRWQU8FugCTBFVWvdr2LkLY/wPVSZuhfZIJr23QfSvvvA075H3qkcZv1nAp2mzsHuBKfA5l6R2BxKw+MFtN2VR8z6Q8AhYBV5/JfVEXA0BrJjwshrGo3Gx6EF+ciJLGLPu8SkGBoM3mMZ0E1EOmE5VtcC43xxo+rMN/awMFq270bL9t1O69rZx4+yb8caDu3ZwoH/vUGvtblFkTKnDZofKqD5tqPY9CiwlkzeY2dDyIy1cbJZJI6EOBp07EbbxBT6nPVb5k3/p0lpNtR1dojIw8B/Xdu/B3bU8FrpItJKVQ+4UgAzXOPVWcB5GngI+CPwNtYqyTNAnW8PsXx3Jp+vSCPjRC5zNx+isIIQlWCl36laaXeI4HBYDsy4IR3K1TldN7QDIlJpDVb3lo1L1UmV3Qb4YPzQCmupkjrEFu0vW4Pl6Vpk7oIPx2Fr2p4tSY9wZNsyYntdUONUNIP36DHoQjYxncwNP9X6/4lY6cD+Q0SmYjXxy1DVPiXGL8Iq8LQDb6tqlXlzLnWef6rqrZUdl5ycrMHU5byyGqzDB3ayaekPHNqYSt6+XdgOH6XBsVwaH3cQe9yKerlRYGdbIatrSxr1TSLp0ltp1cGsRBuCCxFZrqrJgbajJCIyHWvVOA4ru/dRVZ0iIpcAL2HNQ1NV9emaXD9Y5pzvpz5BQolIWfqfxjLylkc4kZnBugVfcmDtIvL37iDi0DEaHc2neabSMLf4/AIb2J3WQ43DBpuuTmT0I/81UXVD0FKT+cb1LPE4cI5raAHwmKpmVnxW0bkdgZnu5xkReQE4UkLkopmq/s1Vw3k3cAmu2k5VHVziOucBo1T1fhH5N/A5loP1iqpeWZkNwTLf1JTluzMZO3kx+WXz7LDmnnB7aUfqkUt7FzkxQKUiEkFH7nGYMgKyDsD4HyHu9BbQDMFHRXNOIBysYUA2VrGne0KyA1uA4Vh5ycuwJFLtwLNlLnGLqma4zvsX8IGqrqjsnqE++bgpyM/js/Hn0XfpcWyAE6sGo5HrgcgJZDSDoy0jKejQmriBw8jLPkbWqsVm9dkQMILRwfI1wTTnnI6ojqOwkF0bl7J18SyOb11Lk9QtdNivpfLmT0TDgTZh5HZpQ+tzLuXMy28jPCKywmsaDP7En/ONp4Ua4AvgY6A9sBu4RlWPiohgKQ5eBOQAN6tqqus6AswGxriO7Ql8gJVldJeq/lqZHcE039SE1+du44XvN5cacztWv0tux+iBbYEQc6Q84SiEadfAzp/h959D5/MCbZHBCwSNg+UypiOlV3zOxFotGunafhBAVcs6V+7zBXgO+EFV51RwTMkC0KQPP/zQ2z9GQNiz6BP6ffBT0Yr0musuIKbTADLXzcWetoMm6Sdome6kiUvc1f1/12GD1IvPoNNl9wfMdkP95PzzzzcOVohSMgLmsMHG5BgiTpwiIS2fZlnWMdlRsL+NnZxOLbE1bIxmHiF26IVmQccQEGoYwToD+AvQkRKlE6p6gXet8w2hPt8s353JNZMWFanwhdlgzKD2jB7YNnSdKU988xdY9hZc9gok3RhoawxeoqI5J1jyPCqSR66Ie4ALgRgR6aqqb5Y9oM4WgKak8H1C8Yr0DUUPMeOLDnEUFrJ20Tds/+fD9NhcgA0Ic8LQb7aQsfAuDnZuTKOzUki5YSK/fPKyqa8wGAweKVsrdm2JOWLlzzPY+t3/sG3ZTvy+PM7YbpWTKOBcMJ3pa5dw9T8+N9EtQyjwCfAmVt2TaSbuZxxOZ5FzZRd44oq+jBvSPsBWeZklkyzn6qx7jHNVTwgWB+u0UNVXsHrT1EtG3vIIVOIM2cPCSDz3CtK3ri5SKnPYYNOARjQ8lEPPVVmEL/+azW9+TZtCaKtQ+Ot0vndf22AwGFxUNN8MOO9KBpxXXBry4Q1Di9KXbQqJs3aw9JdE0no0IeHSsZx71d2mdssQrBSq6huBNqI+krrrKHdPW1lqLDMnP0DW+IitP8B3D0D338KFjwfaGoOfCJZvO6/JIxuKKbv6PMb1kJSxbzu/vPcMjb9fSJt0V65zIdjenc5XeTmMvPUJIiMiAmq7wWAILWJTLqFghWtBxw4bhjSj4f7j9Fp+grBlk5j/8iTSe7UgrF17Cvenmai5IeCISDPXx69F5P+AGUCee7+qem56afAK0xbvZuKX63BXqtgEwsNsReIVIcXepbBrATRoDqeOQMdzod1gWP0RfHU3xHaCq94CW21a4BpCiWCpwQrDErn4DZZjtQwYp6rrvXG/UM9P9hUl6ysACm0Q6YBDMbC7bzPC23eGfXvNg5ChVhiRi/qDJ0GNnRuWsXTK4zRZtYOO+6zvG8WqIT1w/zWMHG9WdA3e43TmGxHZifXPUTzsVlWtuFt4EBFq843V52ov05fsLaoTF+CcbnHcd+EZoVd3tXcpvHMxOAuLx8QGzbvCYVcLs7AouPFry+ky1CmCpgarpOqOiKRRLI98N/A9xfLIXnGuDBVTNsI1+Lc3MffVv9FgyVqSfzkKHEWBgl+m8606uORW8yBkMBgqxlM6Yadeg+j0r5kAfHT9UPoss9IIwx3Q5I2P+WDXRi6bMIkmjUPsocoQ8qhqJwARiVLV3JL7RCQqMFbVbdxNeXMLSkuy220Sms4VwK8vlXauANQJWQeLtx0FVoTLOFj1Br87WKo6toLxb4Fv/WxOvafsA9Hopyy1xY+vH0LvZSeKHoT0vx/zcMMw/nzFn2jaoNYN7g0GQz2k6fmXULDSipo7xWp8PPDTtayZfRY7hrShYbe+FKxfaaLmBn+zEBhYjTFDLVm84wh5ZZyrMJvwxBV9QtO52rMENn9nRawUwGl9tkfC8Cet2itHPtgjrLRBQ70hWGqwDEFGzPm/LXoQAuiQriQ8O433vvqQhg0ak3ZGF3pefCdj+pkJw2AwVI+yUfMLr3+A2S/dh33WPAb9sA/9YZ8VNTeiOwY/ICItsVSMG4jIAIpTBZsA0QEzrA4ztHNz67esEFGiz1VIOlcnDsDH10PT9nDJP+Hg6vI1WAm9rMiVe9tQbzAOlsEjZR+EhvU9iy1/vY8RKxwoxylYvILHC+9i3q5b+NfFfyA63EgxGwyGqikbNb/4r/+Bv8Jn1yTRc00ONiCiEHbOnUHWdX+lcWSDwBlrqOuMBG7CEtZ6scT4CeDvgTCorhMdYUcVLugRzx/O7xqajhVAYZ7lXOVlw/VfWI5UtwvLH9dusHGs6inGwTJUSNkHoTZnDiBjxjIEIdwBN/zoYMWet7lp9f/IaN+fJuGxXNfnShPVMhgMp02ji64gf8N0wgutBe5Bq3P5z52D6RvVEhl0Np82a8CIzmeZ+cXgNVT1PeA9EblKVT8LtD31gU9S0wi3C//8XX+aNQxRtWJV+ObPkLYMrnnfcq4MhjIYB8tQbaJTLoavluFwAgrd98MZ+50UhuXw2NhlbG0rPLniJ+BV8xBkMBhOi5JR84SO3Wn+9QIuX1SIkkb+go/IHGvnycyPgJfN/GLwNr+KyBSgtapeLCK9gDNVdUqgDatLLNlxhGlLd5PcITZ0nSuA1Cmw8r9w7p+h1xWBtsYQpNgCbYAhdIgeOY6OLz5C9vBu7OrXCEURBLsDeu9xi606eG/txwG102AwhCYjb3mEayfN4/wHJ9F62EBAi/r09d6jIAV8t31+oM001D3ewVIxbu3a3gLcFzBr6iDLd2fy+ylLyC1wsnx3Jst3ZwbapJqxeyHMmgDdRsD5EwNtjSGIMQ6W4bSIHjmOoa98xQW33o8IgCI22NDOjqpVH7wnbxELd28MqJ0GgyG0iT7/EsQOimIDVBURWHn0R+btXBNo8wx1izhV/RhwAqhqIeAIrEl1i8U7jlDgsBZiHU5l8Y4jAbaoBhzfBx/fAE07wGjTNNhQOcbBMtSI6JHjiLukPyC0Gn8Z11z+OkOaXsd1nSYi2Pm/OXez91gITqAGgyEoiB45jvYvPkL85QNxNIUxC508PqcRZ+w/xt0/38BVH/2FW754lo/WLAi0qYbQ56SINMcltC0iQ4HjgTWpbjG0U7Oiz+FhNktNMFjYuxQW/Mt697QNsHMBvH2hJWpx7TRo0DQgphpCB1ODZagxsX98nMPfXEnWD3O5rNcAxoyaAHuXknh4CH/L+plRn91Av2Znc1GXc71SM7F8dyaLdxxhaOfmoas8ZDAYqk30yHFEjxxHxFtPsO9f0+i5LJPHV8A/rmnJ8o7fo6dg6fKPgVdMXZahNvwJ+AroIiK/Ai2AqwNrUt0i1lVzdWGvBO46r0vwfIfvXQrvXmr1qrKFQb9rYM3HVuNgWxicc7913IJ/Ws2D7RGQdyKwNhtCAuNgGWpM/qYVgHJyRzY5f3qC9hO3w6+vMehQOH/uEsG/euxh2fE9LFv+CbV9AFqw9RC3Tv8IabCL1xZ25X+/HxM8E7TBYPAp+Tu2uj4JOJSR+yG1A4iAUshbq983DpahxqjqChE5D+iOJWK5WVULAmxWneLXbYcBeOiSnnSMaxhga0qw/Sdw5FmfnQWw6oPifc4CmP986eOdDquvlZFeN1SBSRE01JicebNcnwR1wJEPPmb3j804tKYRyV9F0C3NqplACpm06j2cTmdll6uQ3AIHf/rudSLaTyK8xffYW0/mi42/eu3nMBgMwU10ysWIHVwZXLSNawMa7qr7FNKdiznnvXFc/9mTJmXQcNqIiB24BPgNMAK4R0T+FFir6ha/bDtMm6YN6NA8yPo3F9VR2SAsEi542HoXu/X+u/esV1iUNWaPsJoGGwxVYCJYhhoTnXIx8uUy1OU3nUrLx/LZBZtT6bMbtrSxhC8O6RKGvDOazo37MrrHiHKrzR+tWcDsHQvL9blZuvMIf/3ye07FfOkS1bBWrMOid/rhJzQYQER6AvcCccCPqvpGgE2qd0SPHEfCjUs5OPU7AOI+X84bo5N4d1A/UtoP5r/rP2W/Yz4rs9aycsWn7D42kQ5NW3mcUwwGD3wN5AJrcQldGLxHocPJwu1HuKRPK8T9RR4sHFwHUbFw1t3QaZgVmeo0zIpSdTy3OFLVpHX5MYOhEoyDZagx1kPPYg6+MxsAR5574lRs4eEMGvFn9oWlc0GnwXy8YSbbdTbrc7azfsVXOPRlxvVPASzn6snlfwRxsGR5cZ+b5buO8vtpUwlv9RHitIPNgaoSYQvnih7mgclQNSIyFbgUyFDVPiXGLwJeBuzA26r6XEXXUNWNwJ0iYgPeB4yDFQAcmW7RHAFV4j5bzj9jWhA76nzm7VnK/mOCiAJO/rvzSVd0S1my/EPgFQDjcBkqoq2q9gu0EXWVtfuOk5VbyDnd4gJtSmnyc2DrbOg/Fob9pXi83eDyTpSnMYOhEoyDZagVjmPuXhalV6USxgyhx9U3cZlr+6edS9meV/wA9Fzqw6zPuIXFBxeSkb8OCSsErOjU7B0LGdPvXF6Y/xURbd63LqBhdI68kB15P3Bei1tIjE/0x49nCH3eBV7DcoyAonSg14HhQBqwTES+wnK2ni1z/i2qmiEilwN3Af/1h9GG8kSnXAxfLANVrPlGOTj1OyL7T2NE57NYsvwjlEJQO5HONuTZdxelKD+Zeh/YCnA7XIv2jierINs4WwY3s0RkhKrODrQhdZFftlr1V2d1CSLlQLCcq4Ic6D0q0JYY6iDGwTLUiuiUi2HGsnLjjvBWpbZLPwDZUMnhq/0vAqB2QG1YmRnKGc06sfHgEdbmTUXC3YXsDvq2bM2OXTbWpe/x/Q9mqBOo6nwR6VhmeDCwTVV3AIjIh8AVqvosVrTL03W+Ar4SkW+AaZ6OEZHbgdsBEhISmDdvnld+BoOLyNYkXDIAvlmJVYtlOVnp/3qSbn+4lzGN/o+VWVsZ0LgbAB9l/8eabxDE2RC1ZRY5XHMOvQnAkuXT+XLl5eQ68hjQuBvntuhadLttmQ42HXXQo5mdrrGm300dZzEwwxWlLsD1j0tVmwTWrLrBL9sO07t1E5o3igy0KaXZ8CVEx0H7swJtiaEOYhwsQ62IHjmOmMFTOL50H8VRLCF6+OhSx1mrxC8Xpeh8tfVHVmd95qqrEjpGnE+EPZItJ3/kf1tf5YMtk7FFZGKXMJzqJMIeztW9LuTXfUvZn72KnPxCoiPMP19DjWgD7C2xnQYMqehgEUkBRgORwLcVHaeqk4HJAMnJyZqSklJ7Sw2lSUlh7/ahZG86VjSUu0fp9NnTPPz0Z9BufNH4GWu6F803AE8uv9e1wAOIW4DHwRqdATbYkh1GVoPHaGY/AxSmrpuHRG3nm4NGtbQe8CJwJrBWVdUbFxSR+4HxWP/i1gI3A62AD4HmwHLgelXNF5F7gDuAPcAo19g5wFWqer837AkUOfmFrNiTyS1ndwq0KaUpOAVbvrdk2e3mWcLgfcy/KkOtaXrdrRxf9oRb4AtQmPMYxE0slbM8pl/pflirl3/leuAJ4/o+VzGm37n8edabfJ/+uqvMwsb1Z9xHTKMCkhOSSYxP5Jw2Z/PF7rf5y+e/cOuZ/cs99JheWQZvo6rzgHkBNsPgovldfyT7vidKpwoua0jkj58SfVPF8417gadZg6bM2j+pnLOlFDJz/4s4c9sgYSeJaLfDOk1/5Ok5TZh44cVmTqm77AXWedG5agP8EeilqqdE5GPgWiylwn+r6oci8iZwK1ZN53VAP+DvwEgRmQk8DIz1hj2BZMnOoxQ4NPjqr7b+AAUnTXqgwWcYB8tQa6JHjqPpkHc4tngv7ihWzrJUok9dAjd/67EwtGxEy/0gdCzvOO6HJlA2pB9kypkTis7r12wIX+x+m9k75zNnbQ6/S27H6IFtAfho2R4+X7EPpyoRYTY+GD/UPBAZPLEPaFdiu61rzBACRI8cR7OzJnH013TXiCV6kbN+F5UJQJd0uJLX9CrvbAG2iKPYI4/ifsx2O17rs2dx/WcruWnghTTULmYBp+6xA5gnIrOAPPegqr5Yi2uGAQ1EpACIBg4AFwDjXPvfAx7DcrAECHcdVwD8Hpilqkdrcf+g4Neth4kIszGoY7NAm1KaDV9AdHPocE6gLTHUUYyDZfAKMWNv5viyJ1CH9WRij3BYTfpWT69Qeaf8CnPZWq2wovQeN+mHm+EsbEBEswXk5sfxwRLlw6V7UMBZYu2xoNDJ4h1HzEOQwRPLgG4i0gnLsbqW4oceQwhg7z4Mfv3EtWX94Tu2LYa9S6ul9OXJ2cqXI6zM/K7oenax41SrLjQ8ZhUA7+/8gcLjyby2MNmkDdYtdrpeEa5XrVDVfSLyT6yUv1PAbKyUwGOqWug6LA0rXRksIZ7FwHrgV+BLYGRl9wiVms/vVuXQNUZY/Gvg+tO12v89rQ78gEPCyWnYjuyGHei6/UuONkti/YJfAmaXoW5jHCyDV4geOY6Ei9/l4ExLgCJ9ZQyRTQuJ5vQyLiqKbLlpEZeO7MlF7KeIbj+ZnD234zjVodx1wsNsDO0cZIpFBr8jItOBFCBORNKAR1V1iojcDXyPpRw4VVXXB9BMw2kSPXw0vPc5OBy4I95HN0XTuEyaYHVwO1urMlZx6/c/UeAsIMIezgODJ7Ah/SAfrV6MrfGaIrGdsKZLIGYFT/zQmEeGX2KcrDqAqj7uzeuJSCxwBdAJOAZ8AlxUyf3/i0uhVEQeweorcLGI3ICVvvhnVXWWOSfoaz4zsnJJ++5H/nbOGaSkdK36BF+Q+i5s+U/RZuyJDUWfWxxbSUqXaCO/bvAJIelguYrOn8Ra7fnQVSNhCDCO+CFYC3aCOiDnUAOi+59+YMBTZMtNlmzGJoKiIA4i479CT/ZBT3WlMMfK+nIo/Ot35euzDPUPVfVYw6Cq31KJYIUhuIkeMICWjzzCwccecwWcXGmCO45XmiZYGYnxiUwZ+Tap6alFNZ90hx5NE3l21X04NB9wpw0WsCXvM6597zjX9D2X0QPbmvnGUJILgZ2qeghARD4HzgaaikiYK4pVLjVZRFoDg1X1CRH5GSul8CHgN8AP/vwBvMHCbVbvunO7tgicERu/rHif02E1DzYOlsEH2Px9QxGZKiIZIrKuzPhFIrJZRLaJyANVXEaBbCAKK8xuCALsHfpQsh+Wo/sYr09cyQnJRNojsIsdQQiL3kd4i+9p0GEyyclz+evlUQAcOJ7r1fsaDIbgInbMGJqd1wF3vSaA3XG4VtdMjE9kfN/xpfrsjel3Lu9eNIVrzriGcFs4qOVk2RtuJbLd20xfs4Dr3l7M8t2ZFV/YUN/YAwwVkWgRESwHaQMwF7jadcyNWKmAJXkSeMT1uQHWP2wn1HjdIKB8uWofUeE2cgsdgTOi64Wex8UG9gjoaPrgGXyD3x0srMafpULlJRp/Xgz0AsaKSC8R6SsiM8u84oEFqnoxMAHwamjfUHMcx46X2j762XfkrFzp1Xskxify1oi3uHvA3Vx9xtWIy6FzaAEbsr/nza1/JqHTV3y9aZFX72swGIIPewN3uYxrHsgp8Ml9EuMTefjMh3nnoqn0ik1CXU4WUkBk3A8UhO1k8Y4jPrm3ITCISMOanquqS4BPgRVYEu02rHS+CcCfRGQbllT7lBL3G+A6d4VraJrr3LOB72pqS6BYvusoczcfIrfAyfVTlgRuAaJFD+u9eTfocDYk3wKXvgwXPAQ3fmWiVwaf4fcUQW81/nSRidWbphyhUgBalwgPDyNWBHHJJ6vDwYY33iB7nPf1A7rSFXuenTAJo0CLH6oc6iAnciHbdRk3vn2Sc5t3MU1CDYY6ir37OTBrC0URrFadfXq/xPhEHjzzPm75fjwFzjxXJGsbDaJ3sTCtOUN3G3XBUMMlqd4KWOPqPxUP3AfcBLSu6XVV9VHg0TLDO7CedzwdvxJLtt29/RLwUk3vH2jmbMwo+hxQ0al9ywGB236EqBj/399QbwmWGqzTbfw5GkthpymW+k45QqEAtM6RksLeX34ie2EqIAhKmyZhtPLR7z6FFAZkDODr7V8zY9sMCp2FrtosgEIWZW5hye52Rq7dYKijOKQpxWnJSu6Cr2DsOJ+uSifGJzJ15Nu8sfoNFu5fhIiiFLLixIf8/r+Z/O96oy4YKojIfcBEYBsQKSL/Af4BvA8kBdC0kCeusRVdtkmARafSUiHuDONcGfxOpSmCImJ3dSMPKlT1c1W9Q1XHGIGL4KL5MPeCnyI2iOnTyKf3c6fuTB05ld+d8TurRsJFYU7nopUzg8FQ94gePAjsxXWfx7dHkPPjpz6/b2J8Inf1v4sou5VA4Y5k2dtM5uk5s0w9VuhwO9BdVc8ERmEt2I5Q1ftV9UBALQtx3L3k7krpErhFTlXYlwptk/1/b0O9p1IHS1Ud+KeTuGn8WUeI/s3VNIgrAJuSkJRFdJ9efrlvSUerR0wioIQ1XoetwW72HztlHngMhjpI9IABNB7kqrFAUIWcHccrPcdbuOtB+8QOKlGTVci6oyuM6EXokOtu5quqe4DNqro8wDbVCbamZ9O8YQR/HdkjcBHdzF2QcwTamGCkwf9UJ0XwVxF5DfgIOOkeLFGI6Q1M4886Qs7hcE4djQSnk/QVjYh8969E7/4VrnrLL/dPjE/k7qTx3P3T3UQ2W4DGLuajNfDZik6lVtGW785k8Y4jxEZHsG7/cQSM1LLBEII0HdKBrMUbsaLmSjSr/HbvxPhEJgz9o1WT5cgDhMKczmiBk5fmbOG+C88AYPGOI0SF21iTdpwIu40B7WPJzMlnaGdTsxVg2orIKyW2W5XcVtU/BsCmOsGWjCy6Jfg2g6VK9rl8ZRPBMgSA6jhYia73J0qMKVZ/htPGNP6s2+QsXQZOt8gF5GREEL32Y2jSGob7R/Bx67Gt1gdXLZYtegcFRzsUFdlOW7KHh79ch8NZugnyR6l7GZPczjhaBkMI0TCpHzCL6Ph8WvQ7QbTtILx/Jdwwwy/3d9dkPTj/MfaeSMN5qg0AC7Ye5pethxFxTYkl+GS51V0kzCY8cUUfxg1p7xdbDeX4a5ltE73yAqrKtvRsrhzYJrCGpKVCWAOI7x1YOwz1kiodLFU935s3NI0/6zZWTYQNHE7EBtHxedaONR/7zcFKTkgm3BZOgbMA1EZhTmcEIetUAX/9ZBWfLt+Hejiv0KFMW7KHz1akGWEMgyFEkIzV2CKUyJgCouNciqI7foK9S/0mwZwYn8iEIfdxz0/30KdrOuu2WQ+WSnEtiicKncpDX6wFME5WAFDV9yraJyJR/rSlLnHgeC5ZeYV0S2gcWEP2pULrRLAHi56boT5RZR8sEYkRkRdFJNX1+peIGDkWg0eiBwwg7rx4AFoOyix+4GnoPwUhd21Eg7AGdGjYB+epDjhUeXP+Dj6pwLlyo0BegZPPVpj+1QZDaCDYI5w48st8na2e5lcrzm59Nk0imtCqzUbCbFL1CS6cCo98uc7UbAUBIrJMRO4XkVbAj4G2J1TZkp4FwBnxAUwRLMyDA2tM/ZUhYFSn0fBUIAu4xvU6AbzjS6MMoU10XD4A4dHO4sFDW6wVZT+RlJDEZZ0vY3/uJmy2fI/HhNmEO4d1ZtyQ9ozolVAkRqbAp8vTzAOPwRAK9B+LOiFrXyR7F8SSc9itJFp9J8cbhNvDGdFxBCuOzGf42WuIbLqUqIQZNGg1g4uTcnnmyr6MG9Ke64a0585hnSnpgxU6lZfmbDFzTuC5GIgBdmPVnRtqwNb0bIDARrAOrgNHnqm/MgSM6sRNu6jqVSW2HxeRVT6yx1AHsPe/Aj5+p/SKsqMAdi3wa9f0izpdxMdbPiayyRbyjvXBidWTI8wm/M5DrdXEGWv5YMkewGqM+NmKNJMmaDAEOZkLd1KYY32VZe+LIntfFB0uPEx0y/5+t+WM2DPId+Tzy+EPiGhVPL7oVCqtnKO55uzLSYxPBKB984Y88uU6Cl0FWr9sPcyyXUdNerIfEZF3gMdUdbdrKAb4HfA80C9ghoU4W9KziGsUQbOGEYEzYl+q9d7GOFiGwFCdCNYpETnHvSEiZwOnfGeSIdSxn/l7APYvbkbaInc2qRMa+LfR4MD4gcRExNChy69cdNZOLjpnFcMGrWD0hWtI7LWJlSc+Y1XGqqLjRw9sS4QrjKXAp6kmimUwBDtZs38osWX9/R7fGQ2n/N//Ljs/G/EQOXOog0+2fML42eOL5pxxQ9rz0R1nMqRTM8CkJweIgW7nSkSSsOrAb1fVhzAOVo3ZkpFNt/gA11+lpUKjBIhpG1g7DPWW6kSw7gTeL1F3lQnc6DuTDKFO+gv/BEAdkLU7mjSg7VlZfn/gWXt4LdkF2RzPP86+nEnFO7Lgmz0gCJH2SN4a8RaJ8YkkdYjld8ntiqJYhU5nkfKgwb+IyEDgHKznzl+93BbCUIdoPGI4J3/9tcyoQu4Jv9syqOUgIu2R5DvyceIstz/Pkcfrq17nD4l/KJpz/nZRD655cyEOLU5PvqpEdN3dUsJIuvsEFZFhQHvgaeASVV0vIhFAgD2E0MRSEMzi6qQAOzb7Uq3olfg3VdhgcFNpBEtE7MD1qtofazWnn6oOUNU1frHOEJKcXLSo1Hb2viiwhUHHc/1qR2p6KlqJhJei5DnyeGP1G0WryqMHtiUq3PqzcCrszzRNiv2NiDwCvAc0B+KAd0TkocBaZQhWYseMofFll7q2rL/3qNgC+PUl+OFRv9riFti5Z+A9PDL0Ea454xouaHcBEbbiVKnFBxZz2+zbiuacpA6xXFXiYdThsBZ2AD5YvJtrJi3ihe83c82bC7n9/VQzH3mXO4BHsPpufgXcLyI3ADNc24bTZP/xXE7mOwJbf5VzFI7ugLZG4MIQOCqNYKmqw50eqKr+Xw40hCQNzz2HrK9nFm2rw0ZOuhI951G48HG/1WElJyQTYY8oWk0WBEWL3sFyshbtX8SK9BW8NeItkjok8sH4oTw3ayPLdmUyfdkePltZLNu+fHcmn61I43BWHi0aR5qeWb7hOqC/quYCiMhzwCrgqUAaZQheorp2w9ItE0Bx5NutHb++BLGdIPkmv9mSGJ9YVGflZlXGKt5Y/QYL9y8ErEjW19u/LjpuzKD2fL5in6sey2opcdf/ljNr3cGiazgUZm9I58dNGTxpemd5BVVdAlzo3haRy4GRWA7WlEDZFcoUKQgG0sFyNxg29VeGAFKdFMGVIvIV8Alw0j2oqp/7zCpDSNP2hRfYvnUb+Zs2AQKq5GREEr17IUy9CG75zi9Olns1OTU9lZiIGI7nHy/1/t2u71h6cGlRJMv9wJPUIZZzusaxbFcmTi2ui1i//xiPfrmhlMz7h0v38uQo87DjZfYDUUCuazsS2Bc4cwzBjr1p6c4h9ghH8cbK9/3qYHkiMT6Ru/rfxbKDyyhwFqAoM7bN4LIulxXNOY9c1pNHvtxQ1FKiIhxO5ZEv19G9ZWOzuON9ZgF7gX1aWfqDoUK2FjlYAZRoT0sFBFoPCJwNhnpPdRysKOAIcEGJMQWMg2WokGZjr+Xgo4/hTtkpeuBRh9Wfxo8NQMuuJrvpFtuNm767CYc6UJQvtn1R9MBzTrcWvPLj1qK6iI+W7sHh4evWoeZhxwccB9aLyA9Yv/7hwFIReQVAVf/oT2NExAY8CTQBUitrTmoIDI5jx61aC9czcVEEC6BxqwrO8i+J8YmM6jqKT7Z8AkChs5DU9NSi+Skr1+GKv1WNw6lG5dQLiMibwKuuuqsYYBHgAJqJyF9UdXpgLQw9tqRn06JxJE2jA6wg2KIHRDUJnA2Gek91arCOqOrNZV63+Mk+Q4jiOHbc9ckqMC31wOPn/jQVkRifyOhuo4u23Q88YNVFjBrQpmifJ+eq+Dw1yl/eZQbwd2AuMA+YCHwJLHe9qo2ITBWRDBFZV2b8IhHZLCLbROSBKi5zBdAWKADM/+ggJHrwICTMtV4oZSJYzbsGxigPXN7lcqLsUYCVnrzl6JaiWqyhnZtjL9Og2C5w3ZD2PHNlX4b3SiiaORX4aNleprkEeQw15lxVXe/6fDOwRVX7AknA3wJnVuiyNT0rsNErVStF0NRfGQJMpQ6WqjqAs/1ki6EOET14kNV0CkXsSnR8XvHOAPSnqYjLu1xOhN1aaRMRkhOKc7bHDelAmK28MyhAz5aNSzUKNZLu3kNV36vsdZqXexe4qOSAa+Hodaymor2AsSLSS0T6isjMMq94oDuwUFX/BNxV+5/Q4G2iBwwg9hbXup9C+sqY4obDi17za5PzynCnLV/W+TIAZu2aVSR4kdQhlieu6EOYzRJ6D7MJT47qy9OuBsVv3ZBcKhXZnSpo5p1aUbIL/XDgCwBVPejxaEOlOJ3K1kBLtB/dAacyTf2VIeBUJ0VwlanBMpwu0QMGEN01gZyt+0kYcJzouALXHglIf5qKSIxPZMqIKTy68FH2Z+8vKkJ310WMGVQs2w6Wz/jUKOuBp2Rj4nyHaUzsLUSkG/AslvMT5R5X1c6ney1VnS8iHcsMDwa2qeoO1/0+BK5Q1WeBS8sci4ikUfwg5ii7v8RxtwO3AyQkJDBv3rzTNddQCxru3Im1bi6oE6vuM64AdRay86f32dMhJ8AWluB48cc8Rx4fL/qYYzHHaA08MCiSTUcd9Ghmp/WpHcybV1yP1QkHNrEUTsGKnr/29RJu7BOFoUYcE5FLsWo8zwZuBRCRMKBBIA0LRfYdO0VOviOwAhdprgbDbY2DZQgspgbL4BNyVq4kZ3sGqJC+MobIpoUuJ0v93nC4KhLjE7mux3U8ueRJ3lz9Ju+se6eoN9bogW35bEUa+QVObDbhiRLqXaMHtuWT1L3ku/IHP00t3b/GUGPeAR4F/g2cj5W6U52m6NWlDVYhu5s0YEglx38OvCoi5wLzKzpIVScDkwGSk5M1JSWl9pYaqs3JyEj2zJ5tbdhsRVFzATrnbaBzEP3/aJrRlNmzZ5PnyEMQouKjaNqlKYnxiaRUcl4KEN16D3+fsbZo7NeDyvBBnVm3/zjb0rPIK3QyZlB7I7xTPe4AXgFaAveViFz9BvgmYFaFKFszgkDgYl8qhEdDi56Bs8FgoBoOlqre7A9DDHWLnKXLXMusgjqlaDU52CJYbo7nW0vKilLgLCgqPk/qEMsH44d6bPRZtjFxgYlieYsGqvqjiIiq7gYeE5HlWP1q/I6q5uBa2TYELxJZHMWRss1F9y+Hz26Dq97ys1WeSYxP5O0RbzPxl4nsydrDZ1s+4+vtXxct7FTGuCHtWbbzCDNW7Qcgv9BZyuECWJ22lj1HTvLAJeYhszJUdQtlUohd498D3/vfotBmS3o2QGBTBNNSLfVAe3XiBwaD76hyVVhEzhCRH91F4iLSzzT9NFRF9OBBYLeELSQsrEQNlsKSSZD6bsBs88SgloMIE2tCDrOFlarFSuoQyx/O7+rRcRo9sC0RduthToFPl1u1WMt3Z3LT1CVc8M+5pjno6ZPnUu7bKiJ3i8iVgDeXRPcB7Upst8XIwIc8OcuWFX1Wh5OcjMjSB2ye5WeLKicxPpGhrYYC4MRZtLBTHQZ1qjoLYPKCHWbeCSAi0lREPhWRTSKyUUTOFJFmIvKDiGx1vce6jr1KRNaLyAIRae4a6yIiHwX2pzg9tqRnEd84kpjo8MAYUJALB9dCGyNwYQg81Um7eQt4EEtBC1VdA1zrS6MMoU/0gAHE/eH/AGjUv13pndnpMPPeoHKyEuMTeebcZwDo0axHtc9zR7Hc5Bc6uWfacq56YyHzthxmx+EcZm9I53dvLjSKX9XnXiAa+COWmtf1wI1evP4yoJuIdBKRCKz57CsvXt8QAKIHD7Kk2gGxU1pYB6BhiwBYVTmXdiku+Qu3hZda2KmMzJz8KrVYnYpRNw0sLwPfqWoPoD+wEXgA+FFVuwE/urYB7gEGAZOAca6xp4CQWszemp4d2Pqrg2vBWQBtBwXOBoPBRXUcrGhVLSvBVOgLYwx1C1uDaACylu9iz9zmxapebla+HwCrKqZVw1YIwupDq4uUvapDySgWwP7jeeWOcSpG8auaqOoyVc1W1TRXW4jRqrq4JtcSkelYvW26i0iaiNyqqoXA3VgpQBuBj0tINRtClOgBAwjv0J6ITp1o/8cLSgjruDi+N2jUBN0MiB9AckIyTSKaVCs90M3Qzs0Jt5d2sWwCHZpFlxr7aKmRcg8Erp5aw4ApAKqar6rHsFo+uJVQ3wNGuT47sRqqRwMFrnrPg6q61Y9m1wqnU9mWkU23QNdfgRG4MAQF1UlSPSwiXXD1PxSRq4EDPrXKUCfI37nT+qBapg7LRZA0AHVTMj0n35FfqgloZZStxaoI0xy0ckTkHKCzqr7v2v4UaOba/ZSq/nS611TVsRWMfwt8W1NbDcFJeFwLECH69tfhhzhYMhkKXeqBzgK/NjmvLsM7DCc1PZW4BnHVPiepQyzTbz+Tz1akcTgrjxaNIxntEtgpqW5qGqFXDxGJBK4COlLiuUhVn6jhJTsBh4B3RKQ/Vv++e4EEVXU/Px0EElyfnwXmAPuB32OpNleaKRRsqqUZOU5OFTjQY/uZN+9QQGzoueFbmkY0Z9GKLcCWgNhgMLipjoP1ByxlrB4isg/YCVznU6sMdYLo5GSOffQRiCBh9vIpO0HUABQgOSGZcFs4+c58bDZbtdN1oLyiIFjqZc0bR3AkKx/FWqH4OHWvURqsmMexUmXcdAduAhpiNR4+bQfLUL9wFhZSuH8/OStXEj38ccg7AalTiw9Iq16Nkz8Z3NJy+P6V+i9u7H1jtaNYSR1iK6wL/WjZXgpdWu7uhR1PzpihiC+xxPOXA+VTEE6fMGAgcI+qLhGRlylOBwRAVVVE1PX5B+AHABG5AWvx5wwR+QuQCdzrEtspeX5QqZb+uDEd5qdy2bAkkjo0q/oEX7D6XuhyFoH+XRgMUI0UQVXdoaoXAi2AHqp6jkvZK2CISHsR+UJEporIA1WfYQgE0YOtPOhGF1xA+/f/S3RSmbzoIGoACq6eWCOn0CSiCX2a96n2gw4UryiPG9KeEb0SuG5Iez696yxSJw5neK+EouMKHWrqIiqmiapuKLG9VVWXq+p8IICJ/YZQIGflSnLXrKEwI4M9N99CzsqV5ZuaH1wLPzwaGAMrILvAUl6bs2fOaaUmV0RSh1jGn9OpaFuBaUv2MG3JHmZvSOeDJXsYM3mRSVcuTVtVHaOqz6vqv9yvWlwvDUhT1SWu7U+xHK50EWkF4HrPKHmSiERjLSq9jrXgdCPwCyGwqO1WEOwaKAXBk4chc5dpMGwIGqrdW0ZVT6pqVm1v6HKKMtyqhCXGLxKRzSKyrRpOU1/gU1W9BRhQW5sMvsHe2JpoNc+1IHjh41b1uRunA3YtCIBlFZMYn8hlXS5j/ZH1vLHqjdN62EnqEMszV/Zl8g3JPH1l36IV4haNS6uZHc7yxgJpnaRpyQ1VHV1iMwGDoRKs1hBOADQ/39reNrv8ganv+NmyyimZmpzryOXr7V/X+pqNG4RXKoJR6FDe/Hl7re9Th1goIn29dTFXP629ItLdNfQbYAOWmI5bsOdGrMhZSf4KvKKqBViNjhWrPiuaIGdrehYtm0QR0yBACoL7llvvpv7KECR4s3lndXmXMn0nRMSOtWJzMdALGCsivUSkr4jMLPOKBxYDt4rIT8B3frbfUE1yN20G4OSvv1oryofDIblkWzWFjcHXy7Fto7YUOAt4Y/UbXllRHj2wLWG24sednzZncPv7qYx67Rcj416aTSLy27KDInIpsDkA9hhCCHvTmOINp9PazjpY/sC840EVOXenJrv5fNvntZ5zhnZujt1Wuc7gjxvTzbxTzDnActci7xoRWSsia2p5zXuAD1zXSQSeAZ4DhovIVuBC1zYAItIaGKyqX7iGXsVSPL0TmFZLW3zOqr2ZRIXbAvdvKi0VxAatEgNzf4OhDH53sFzpPkfLDA8GtrnSEfOBD4ErVHWtql5a5pUB3Aw8qqoXAOUeyAzBQVFfGtXiFWWno/RB7gagQcTJgpOA1XQ4z5FX6xXlpA6xjBlULOVe6FBmb0hnVdpxI+NemvuBF0XkHRG5x/V6F3jRtc9gqBDHseNFMu0AuRs2woAbPB+8OnieVxPjExnVdRTiijkVOgv5z6r/1MrJSuoQyxNX9MFeiY9lZNxLcTHQDRgBXAZc6nqvMaq6SlWTVbWfqo5S1UxVPaKqv1HVbqp6oaoeLXH8flX9bYntT1S1t6qeraqBUY2oJqm7jrLjcA67j+Rw3duLA+Nk7UuF+F4QGUAVQ4OhBFWKXLhygv8MtFfV20SkG9BdVWd60Y42wN4S22nAkEqO/w54TETGAbs8HRBsCjv1kQYZ6cWFM04n2zPSsecspWRptQKFG2bya/N5/javQiLyIrBhw4kTRfl8y+e0PdGWTpGdqj65AjrhwC5QQgOjFE6Fh79YS87+rXSNtXs+qI6jqttEpB9WvUFv1/B84E5VzQ2cZYZQoKi5eaHVReT4jBnEjHqX6Etfhh8esSJXRVTVRcq/XN7lcr7e/jW5Duuf+eIDi1mZsfK0pNvLMm5Ie7q3bFxK3GJrehZLdxU//Nb3dGURaaKqJ4Balz/UZ37aZJWSKVBQ6GTxjiP+FVFxOq0UwV6j/HdPg6EKqqMi+A6Wss6Zru19WBKi3nSwTgtVXQdcXcUxQaWwUx85vHkLRctuNhtd4hOItUVZWk0uBAiP6xJUqj8ppLB30V4+2fIJACqKo7WDlL4ptbgm7GRtpVLuToW8ph1ISQkudUV/oqp5wNQqDzQYyhA9YABNLrqIEzOtryZ1OMhZuozoO263Dph5b/HBR3f538BKSIxP5K0Rb/H6qtdZfGAxilLgLKh2q4iKKKs0uHx3JmMmL6LQtdLz48YMnvt2I9sPnyTjRC5jBrVn3JD2tf1xQolpWNGq5Vj+QUnPW4HOgTAq1OgU1xCwerGFh9kY2rm5fw04uh1yj5v6K0NQUZ0UwS6q+jxQAOCSCvX28t8+oF2J7bauMUMIEz14ENisf2ISEWFte0rZ6fobP1tWNZd3ubyoLsImpyfZXhFlGxIDxEQXr3EokHWqTHNUg8FQbZpc7srqEkHCw4uUTDm4qvSBO34KutTkxPhE/pD4B2xizZnhtnCvzDslSeoQy5jk4q9ahypvzt/BDxvSWZ12nL/PWMtz32706j2DGVW91PXeSVU7u97dL+NcVZM4l5DT2MHt+WD8UP+3AHC3XzAKgoYgojoOVr6IuNVscDUd9nZewTKgm4h0EpEIrAZ7X3n5HgY/Ez1gANGDrLSdhAcfIHrAAEi+CXpcWvrAIJNrh+IV5Qh7BPHR8V65Zkkp9+uGtOezu87i9nO7lDrm7V92msJzg6GGNBw61Ho/80zavzPVmnMAsj2UsKz9OCjnnau7WckZ/075d62iVxUxemDbSmuzJs3fYeYgw2lxyJVqeud5XQLTX21fKkQ0ghbdqz7WYPAT1XGwHsOqeWonIh8APwJ/q+kNRWQ6sAjoLiJpInKrqhYCdwPfAxuBj1V1fU3vYQgOclauJGf5cnA4SH/2OasvDcDZ94KtRHaqsxDmBFdvGgC72HE4HezL3sf42eNrrewFxVLubhn3oZ2bl1IYdDgt+eQrXvul3qoLishlIhIIhVNDiGOLiIDwcKL69Cl2rgAaVbBIEkRiF24u6mSJ7Dpx+uT6SR1i+U3PirseKBgJd8Np4Xaw4hpFVnGkj0hLhdYDwFY/65cNwUl1Gg3PBkZjNb+bDiSr6rya3lBVx6pqK1UNV9W2qjrFNf6tqp6hql1U9emaXt8QPOQsXQYOSzVQCwqsbYB2g+HMu0sfvHth0DUATU1PxanWQ06Bo6BUvxpv4akpqDtdZ/aGdK6ZtLA+OlljgK0i8ryI9Ai0MYbQQiIiOLlkSfGCDkD/saV78LnxFNkKML2b90YQ3lv/nlcWdTxxx3ldsFfy7W8k3A2nw6GsPBpHhtEgIgAOTsEpSF9n6q8MQUeVDpaIfI0lXTpPVWeq6mHfm2WoCxSpegESFlZcDwGQ5iE1Z83HfrKseiQnJBNhjwC8V4flicaVNGZ0OOvfarKq/h6rgfh24F0RWSQit4tI4ypONdRzclauRE+eJHf1aqv3ntvJajcYbvkOGgV/v+otmVsAWHpwqVf68HkiqUMsH99xFsN7JdC/bQx3Dutc6mHAqfVj3hGRZpW9Am1fqHAo21KpDAirP7SyYCJjqj7WYPAj1UnD+SdwLrBBRD4VkatFJMrHdhnqANEDBtB8/HgAWj//j9IpO5m7yp8QHlz/rBLjE3l7xNs0Cm9Ecstkn9RDAOXSBMtSH1eTXdLJn2L1xGsFXAmsEJF7AmqYIagpipJTJmoOlpPVdlDpE/Ys9pNl1Sc1PRW1Sp6LlAR9QVKHWN66IZkv7z6HBy7pyYW9SjufP2xIZ9BTP9T1VOXlQKrr/RCwBdjq+rw8gHaFFIey8oqELvzK3qXw7V+tzz8/F3Q1lYb6TXVSBH9W1f/DkiudBFwDZPjaMEPdIKqHVXQa0alMD6m+15Q/+OiuoJsgE+MTSU5I5lCO71KJ3E1BK/pjrC+ryW5E5AoRmQHMA8KBwap6MdAfqyefweCR6MGDipoNl1IRrIicw/D+lX6wrPokJyRjd6Uz+kJJsCLuOK9LuTnoUHZ+nU5VLqEWOAe4TFXjVLU5lnT77MBaFzq4+6z5nV0LwOlS3nUUWNsGQ5BQrUJyl4rgVcCdwCDgPV8aZag7SJQVldLcMn1ihz/uwclyBmXR+RnNzmDXiV3kOXzXlHPckPZ8cldxyk6vVqWz4X7alFEnH3AqYDTwb1Xtq6ovqGoGFLWIuDWwphmCmegBAwhr1YqIbt1Kqwi68SR2sWu+f4yrJonxiZzX9jwahDWoVaPh0yWpQ2y5KJabepCqPFRVv3VvqOos4KwA2hNSHMrKo0UgBC46ngtuPSR7hLVtMAQJ1anB+hhL2e8C4DWsvlgmTcdQLWxRDQDI/PiT0kXnAFe9Bd0vKT12aIufLKs+3WO741AH/0z9p8+KzqF0ys6To/piL6MuOOHT1UycsbY+OFoHVbXUU6+I/ANAVX8MjEmGUMHWoAGRnTuXd67AErsoi7Mw6CLnCQ0TiLBH+M25clOZ+EUdT1XeLyIPiUhH12sisD/QRoUCp/IdZOUVBiaC1W4wNGkDLXrAjV9Z2wZDkFCdCNYULKfqTlWdq6q+0Y411Eny0/YCcPyzz0oXnbuJO6P09u5fIfVd/xhXTdxKgh9t+shnRedlSeoQy23nlE6r3HboJB8s2cPv3lzItCV7fG5DABnuYexiv1thCE0EUPW8r91guPTl8uNBFjkPs4VR6Cz0+31Lil/ENYootc+p8NysOtuEeCzQApgBfO767MEbN5TlcLaV2REQB0sVco5C5/ONc2UIOip0sETkAtfHhsAVIjK65Ms/5hlCnbzNm60PquWLzgE2zyp/0uL/+N6w02BvluUkKurTovOyVKQu6FR46Iu6F8kSkbtEZC3QQ0TWlHjtBNYE0K5eIvKxiLwhIlcHyg5D9RCppIsuWM3OW/YtPbbjZ5/ZUxPCJAyH0xGQe7sj6akPDWdEmZTBZbsySXlhbp2be1T1qKreC5yjqgNV9T5VPRpou0KBjKwAOli5x6DgJMS08f+9DYYqqCyCdZ7r/TIPr0t9bJehjhDV1/UgY7N5Ljr39DBUmFt+LIAMamnZLIhfi86Hdm6OvYJnxToqfDENa375ktLzTZJLuv20EZGpIpIhIuvKjF8kIptFZJuIPFDFZS4GXlXVu4AbamKHwc9UFMFyc+p46e2jO4Iqch5mC6NQ/R/BKssd53Wh7BS060gOV72xkOe+rTvRLBE5S0Q2YJVDICL9RSS4VvqCFHeT4YDUYB3fZ703MQ6WIfio0MFSVXfX1ydU9eaSL+BJ/5hnCHWielh9YptcfLHnovMhd3k4qYkfLKs+ifGJRNgjGBA/wO9F50+O6ktFCu51sCZCVXUX8Acgq8SLWvSkeRe4qOSAiNiB17Ecp17AWFeUqq+IzCzzigf+C1wrIi8AzWtoh8FvCFCFgxURXX5swb98Yk1NsNvsFDoL0aocRR+T1CGWO4Z19rjvzfk76lKq8r+BkcARAFVdDQwLqEUhwiFXimB8ICJYJ1wOVkxb/9/bYKiCsGoc8xkwsMzYp0CS980x1FUaDx/uueg8+SZY+xHsXlg8dnC9VXQeRDnVdrHTr0U/vxedjxvSnu4tG/PZijS2pWexdFexQ+WOYr11Q53pYD8NKzq+HOsJuaRrqVitIk4LVZ0vIh3LDA8GtqnqDgAR+RC4QlWfpeLo/B9cjtnnFd1LRG4HbgdISEhg3rx5p2uuwQs0y8nBkXGIbZX8/lvFXsAZhzYBxe5YYVY6vwbJ/7O9x6y05Lnz5mKTaon9+oyh0fBLMxvrjpYvv37t+7W0PrUjAFZ5H1XdWya9NDA5miHGoaw8RKBZw4iqD/Y2x9OsdxPBMgQhFTpYItID6A3ElKm5agIEV0dYQ2hz4eMwZSTg/gJ3wq8vwbXBVXgeKJI6xJLUIRaA299PZfaG9KJ9P2xI55KX5/PkqL5Fx4Qqqnqp671TVcfWkjbA3hLbacCQig52OWh/x6pHfaGi41R1MjAZIDk5WVNSUrxgquF02fHvlwiPi2NApb//FHj5O8i0nAMBwp15pHSJDoqFnW1rt8EKOHvY2UTaAxAZKENKCtz34Uq+WFVaWG9/Duxv0JlxQ9oHxjDvsVdEzgJURMKBe3GlCxoq51BWHs0bRhJWkfykLzmxH8QOjVv6/94GQxVU9hfRHWs1tyml6yEGArf53DJD/aHd4PIT5KZvgk46ORjw1Ax0w4Esrn4j9JUFRWRgZa9A2aWqu1T1dlW9TlV/CZQdhmoiUnUNFsDoSeXHfn3J6+bUhDCx1j4DJXThiZeuHcAzV/YtN15HWkfciZWa3AbYByS6tg1VcChQTYbBShFs3Aps9sDc32CohAojWKr6JfCliJypqov8aJOhPpKfVX7MRLHK4W4GWjKKBVaK00NfrKV7y8ahHMmqrAhGsXrxeYN9QLsS221dY4a6QHUdrHaDoVFLyD5YPLZrYcXH+5Ewm/XVXOAsCLAlpRk3pD3zNmeUmn8UuGHKEib+tlfIRrJU9TBwXaDtCEUOZQfQwTqeZhQEDUFLdWK6d4pIU/eGiMSKyFTfmWSol5zhoc3R4W3+tyMEuOO8Lh6FL0JdWVBVz6/k5S3nCmAZ0E1EOolIBHAt8JUXr28IJFWotFdKbmZQqAnaXSvyDg2eCJYbT8qCJ/Md/H3G2pBTFhSRv7neXxWRV8q+Am1fKHA4Ky8wCoJgRbBM/ZUhSKmOyEU/VT3m3lDVTBHxoFYQWhQUFJCWlkZubnBJgvuKqKgo2rZtS3i4595KAeeqt2D7T5BzuHisYVzg7AlikjrE8smdZ/HQjLVsPFg68udWFgzhKBYAItIHS+GvqN5TVd+vwXWmAylAnIikAY+q6hQRuRv4HrADU1V1vVcMrwQz5/iR6qrvNWhaOoIFVh++5Ju8bdFpYReXgxVEKYJu3MqCb84vL27x5vwdDO/dMpTmnw2ud580N3QJ46QC+1T1UhHpBHyIpUa6HLheVfNF5B7gDmAPMMo1dg5wlare7wvbvIGqBi5FUNWqwerxW//f22CoBtVxsGwiEquqmVAkl1yd84KatLQ0GjduTMeOHatuTBniqCpHjhwhLS2NTp18rR9QC9oPhU0zA21FSJDUIZZZ9w1j2pI9/H3G2qLxuqAsKCKPYjlFvYBvseTUfwFO28FS1bEVjH/rurbfMHOOfxCqmSIIVpuImfeWHss97vlYPxJus5zSQmfge2F54oFLenLwRG450QuA52Zt5JM7zwqAVTViDDATaKqqL/vg+m6xDHfvkX8A/1bVD0XkTeBW4A2s9MR+WGI6I0VkJvAw4HH+ChZOnCok3+EMjIOVc8TqmdnESLQbgpPqpAj+C1gkIk+KyJPAQuB535rle3Jzc2nevHmdf9ABEBGaN28eeivnu381QhdVMG5IewZ3LL1a/MOGdO77cGWALPIKVwO/AQ66+u71B2ICa1LtMXOO326OVtUHy03yTdCgTLTFke91k04Xd4pgMDQbroiXrh3AnR56ZC3blRlKgjtJItIauMVV/tCs5Ks2FxaRtsBvgbdd24JVR/qp65D3gFHuw4FwIBooAH4PzFLVo7Wxwdccyrb+vgPiYLkl2k0NliFIqTISparvi0gqxQXmo1V1Q2XnhAr14UHHTUj8rI3iy4/NuAP+GNLOgs/pltC4VH8sgC9W7efoyXzev7VC9fFg5pSqOkWkUESaABmUFqUIWULi79BLBOxnra7IhZsOZ5eOnJ86atVhBTBN0K0iGKwRLDcPXNIToFy64MQZISO48ybwI1aPveV4ofdeCV4C/gY0dm03B46pFnnNaViqhQCvAYuB9cCvwJdYjY8rJBj67m08YqWw7tu2kXmZW/x67+aHl9AXWL71IFnp8/x6b4OhOlQ31a8ZcFJV3xGRFiLSSVV3+tIwQz2k/1hILaOfcnRH0DUdDjZGD2zL9KV7cJZ5ppy/9TDPfbux6CEohEh1Ceu8hfXQkw0YJVND9RChugEsAM6+t3xq8oJ/BdTBKhK5CMIarLI8cElPVuzJLLXIo8Ct7y5jyk2DgtrJUtVXgFdE5A1Vvctb1xWRS4EMVV0uIinVsOO/wH9d5z4CvAJcLCI3YPXs+7OqOsucE/C+e8dX7YNlqxh+7hC6xjfy782XbIF1kHT+5dA4wb/3NhiqQZUpgq56iAnAg66hcOB/vjSqvtOxY0f69u1Lv379GDFiBAcPWkXY2dnZ3HXXXXTp0oWBAweSlJTEW2+9FWBrvUi7wRDrYcEwSHrTBCtJHWJ5alT5/jRgrSyHWo8aVf0/VT2mqm8Cw4EbXamCBh9Q5+ab041gtRsMEY1Ljx3fE9D05KI+WEGoIuiJCRf3LKcseOxUAVe/sTCo5x9XhBxgYtn0wFqmCJ4NXC4iu7BELS4AXgaaioh7YbtcewhXuuJgVf0C+DNWjdgxrJTpoONQVh4QoBTBE2lgC4eGLfx/b4OhGlSnButK4HLgJICq7qc45O1zRKSziEwRkU8rG6trzJ07lzVr1pCcnMwzzzwDwPjx44mNjWXr1q2sWLGC7777jqNHgzpF+/Tx1PzTyLVXybgh7fnsrrNo2qC8YluoSbeLyDD3C2iP9VAyLNB21WXq7XzjJjyq/FgAF3bcfbCCPUXQjVtZsCwKPFRChCcIcTdaXI6l9re8xKvGyoKq+qCqtlXVjlhtIH5S1euAuVg1pgA3YqUCluRJ4BHX5wZYv0InVm1W0HEoO4+IMBtNogKge3Z8HzRpDbbqPMYaDP6nOv8y81VVcSVdiEjD6l5cRKaKSIaIrCszfpGIbBaRbSLyQGXXUNUdqnprVWP+YPnuTF6fu82rK3KjRo0iKSmJ3r17M3ny5HL7hw0bxrZt29i+fTtLly7lqaeewuaaUFq0aMGECRO8ZktQ0G4wNCvzRR0WERhbQoykDrFMuWlQufE5G9KDehXZA38t8XoY+Bp4LJAGBQpvzzn1Yr4RTi+CBZDoocfsgcA5BqEgclGWBy7pybBu5VtrbDyYFbSiF6p6qeu9k6p2dr27X7Wpv6qICcCfRGQbVk3WFPcOd/sbVV3hGpoGrMWKhn3nA1tqzSFXD6yA1Fue2AcxRkHQELxUZ9nhYxGZhLWKfBtwC1ZtRHV4F6t4s0he2dUX4nWs1J80YJmIfIXVj+bZMuffoqoZ1bxXjXn86/Vs2H+i0mOycgvYdDALp4JNoEfLxjSOqri/S6/WTXj0st5V3nvq1Kk0a9aMU6dOMWjQIK666qpS+2fOnEnfvn1Zv349/fv3L3rYqdPYyvxej6UFxo4QJKlDLCN6JTB7Q3rRmHsVedZ9oREEUtXLSm6LSDusgvE6Q6DmnPow35yWTLub4Y/DsimQX6KvXAAFSUItguXm/VuHcMOUJczferjU+EfL9jBuSPsAWVU1IvKjqv6mqrGaoKrzgHmuzzsAjwXFqroSS7bdvf0SQT7vBawHFlgRrPZDA3Nvg6EaVEdF8J8iMhw4AXQHHlHVH6pzcVWdLyIdywwPBra5JhpE5EPgClV9Frj0dIyvjKoUdmJiYsjKsr5MC/ILcDgqz3U/npNfJCLgVGs7Orzih4+C/IKi61fGCy+8wMyZVoH13r17WbVqFarKeeedh91up3fv3kyYMIGFCxdSWFhYdM0XXniBL774gkOHDrFlS/XUe3Jzc/2uNGTft484YP369eRFVW8i7k1T4nAtRAPkZrJnys3s7HKj7wytAofDwd69e5mXPS9gNlSXwY0d/EDpOv+NB7MY8Mg3/GFAFF1j7YEyraakASGn1FFbTuQWlppzTuQWVupgVYdXXnmFGTNmANZ8s3XrVgDOP/987HY7/fr146mnnmL+/Pmlznv66af55JNPyMjIYP/+8r2PggopmjlOj1b9rNYQbrLSKz7WxwRzo+GqeP/WIaS8MJddR3KKxtamHeeaNxcy4eKeQSV6ISJRWOl3cSISS7GKYBOKFf4MFXAoK492zQKQveh0QNZ+I9FuCGqqlTjrcqiq5VRVgzZYqjhu0oAKtaRFpDnwNDBARB5U1Wc9jXmwuVKFnY0bN9K4sVVK9tRViVUavXx3Jte9vZiCQifhYTZeGZdU6y+KefPmsWDBApYsWUJ0dDQpKSnY7XZEhJ9//pm4uOJ0i6SkJB588EEaNmyIzWbjiSee4IknnqBRo0ZFP0dVREVFMWDAgFrZfLrkbtnCTqB37940qa7KUZdomDIcKP6267D/Wzrc+o4vTKwW9g/stGvXjpTklIDZUF1SgOjWpRsQA2TmwzNLc/nkzrOC6iGnLCLyKsVPyDYgEVhR4QkhSHWi22XnnJevHVCr/2/z5s1jzpw5LFq0qGi+cfepmjt3bqn5plevXqxevRqn04nNZmPixIlMnDiRRo38rBRWE0TQ041gAbToXtrBcuTCv/vC/f5PFQzVCJabs7vGsetIcVqgE1i6K5Pfvbkw2OafO4D7gNaUlmk/gZV9Y6iEw9l5DAzE/8uTh8BZCE2Mg2UIXioMwYjIL673LBE54eG1U0T+z9cGquoRVb1TVbu4HSlPY74mqUMsH4wfyp9GdOeD8UO98gVx/PhxYmNjiY6OZtOmTSxevLjCY7t27UpycjIPPfRQUbQtNze3Zg8SwU67wWAvU3flyDVNh0+DcUPaMyqxdblxp8JnK4I+5bJksfkiYIKq/j6wJvkfb8859Wa+OV2Zdjf9x5YfC5CaYFEfrBCqwSrJ6IGea2OCbf5R1ZdVtRPwlzI1WP1V1ThYlVDocHLkZD4tGgWiybBLfNHUYBmCmAojWKp6juvdY3jEFUVaCPznNO+5j9JNQ8tJlQYrSR1ivbrydtFFF/Hmm2/Ss2dPunfvztChlecTv/322/z1r3+la9euNG/enAYNGvD88897zZ6gotcoWPtx6bHV00w/rNPgpWsHsOvwSValHS81vi296tTVQKKq74lIC9fnQ4G2J5B4c86pN/NNTWun2g2GmPaWU1WSAMw7odQHyxNJHWIZldiaL1aVTycNxvlHVV8VkbOAjpR4LlLV9ys8qZ5z5GQ+qgGUaAcTwTIENdVKERSRgcA5WOuCv6jqSlU9Up0Geh5YBnQTkU5YjtW1wLgaXCfkiYyMZNasWeXGd+3a5fH4Jk2aMGmSBxnzushVb8GuX6w8azfZ9fpZu0Z8cfc5XPLyfDYcKH6oSd2dyfLdmcGUpgOAWFJUjwJ3Y0XXRUQKgVdV9YmAGlcHqFfzTU0jbfevhadbQ8HJ4rEdP3vHptOgKEUwRCNYYC3wtGwSxZRfdlJQogv60l2ZQVePJSL/BboAqwC3V6uUEOgylCagPbBMBMsQAlSn0fAjwHtYkqJxwLsi8hCAqh6o4tzpWCk+3UUkTURuVdVCrAeo74GNwMequr52P4ahTtJmYOntfXWqDMdvPDmqb6k/dKfCDVOWBKN08v1YksSDVLWZqsZi1WeeLSL3B9Y0Q8hQE5n2kvQfU3r76A744dFamXS6FDUaDtEIlpsHLunJh3ecia1MUHHprkyumRRUTYiTgbNdTc7vcb3+GGijgpnANhneB2ENoEFwOOgGgyeqo8F7HdYDz6Oq+igwFLi+OhdX1bGq2kpVw11N96a4xr9V1TNcNVRP19x8Q53m2O7S21n74bPbAmNLCJPUIZbuLUtn+p7Md/D3GWuDzcm6HhirqjvdAy610d8DNwTMKkNIUSOZ9pIc3VV+bM3H5cd8iDuCVeAs8Ot9fUFSh1jO7lq+P5bDGVT1WOuAloE2IpQocrACUoOVZikIBrCVgsFQFdVxsPYDJdvcRxIiNVOGECfnSPmxzeVTnAxVExHm+U996i87/GxJpYSr6uGyg646rNrpkxvqD1JLB+vAyvJjDZvX/Ho1oKgGS0M7guVm79Ecj+OHXQ/pQUAcsEFEvheRr9yvQBsVzBzKDnAEy9RfGYKcylQEXxWRV4DjwHoReVdE3sFa6TnmJ/sM9Zm+15Qfa9jC/3a4CAkFtQoYM8hzk8/th04GU5pOfg33GQzFiKA1khF00XV4+bHIJjW/Xg0I5T5Ynriot+fg0I7DQTP/PAaMAp4B/lXiZaiAQ1l5NI4KIyo8AH0Vj+6E3GNGWdgQ1FQWwXJLJc8A/g7MxepGPhH40ueWGQzDH4fYzqXHmrQKjC0hzrgh7Xnmyr7ENSotf68EVZpO/wpaQmQBfQNtnCFEqKlMu5ur3io/7+z+1a8Pc6HeB6ssD1zSkzuHdaZBeOlHjm0Z2Vz95sKApyqr6s/ALqwo+s9YYlym6LcSDmXlBSZ6tXsR5ByGA2vgvcuNk2UIWip0sFT1PVV9D/iI4p40H5UYNxh8T5eU0tt+ftCpS4wb0p7Uh4YzuGPpwuAZK/bx3LcbA2RVMapqV9UmHl6NVdVvKYIi0llEpojIpyXGRonIWyLykYiM8JcthhrgjbqMs+8tP/bNn2p/3WoS6n2wPPHAJT353/ih2Ms8dajCQ1+sDWgkS0RuAz4F3LKZbYAvAmZQCHAoKy8w9Vcb3ZmbCo582LXA/zYYDNWgshTBMBF5HkjDUhF8H9grIs+LiKmH8CEdO3akb9++9OvXjxEjRnDw4EEAsrOzueuuu+jSpQsDBw4kKSmJt956K8DW+hhPzT/9+KBTF+mWUFrw4lSBgzfn7wgKJ6u2iMhUEckQkXVlxi8Skc0isk1EHqjsGqq6Q1VvLTP2hareBtwJjPF8ZmhSJ+eb2qbznvJQ/5m5q3bXPA3qWgTLTVKHWI9NiJ0Ki3d4+J37jz9gKZieAFDVrUB8IA0Kdg5lByiC5frbQOxgj4CO5/rfBoOhGlSWIvgC0AzopKpJqjoQq09EU+CffrCtXjN37lzWrFlDcnIyzzzzDADjx48nNjaWrVu3smLFCr777juOHj0aYEt9TLvB1iRakoNrTRSrFvRuHeNx/OPUvX62xCe8C1xUckBE7MDrwMVAL2CsiPQSkb4iMrPMq6qHqodc16pT1Kn5prYy7QANPIhaRPlPErquiVyUpKDQ6XE8NjrC47ifyFPVojpPEQmjdommdZ6ApQgeT4Po5nD+RLjxK783ATcYqktljYYvBc7QEpX9qnpCRO4CNgEecijqOHuXWuHojud67Y961KhR7N27l9zcXO69915uv/32UvuHDRvGK6+8wvbt21m6dCnTpk3DZrP84hYtWjBhwgQADhw4wJgxYzhx4gSFhYW88cYbnHtuXVnZ8ZDys3qamVhrSGaOZ72IzJyCoGxAfDqo6nwR6VhmeDCwzSX5joh8CFyhqs9izXNV4mqC/BwwS1U91maIyO3A7QAJCQnMmzev1P6YmBiysrI8nFkxtv3LCdu7iMJ2Z+JsnXRa53pi7Nix7Nu3j9zcXO666y5uvvlmVJXs7GwiIyNJTk7mzTffZPXq1SxZsoRJkyZx8qTVdDcqKor/+7//Iysri4MHD3LTTTeRlZVFYWEh//73vznrrLPK3S83N7fc78HXNM08huTn1+q+7Xen0onimUcBju9hy7QHOdB6ZO2NrAJ3auDW7VuZd3iez+/nTxZt8awo+MRXa/ll5Uau6RGAh3b4WUT+DjQQkeHA/wFfB8KQUCAnv5DsvEL/O1iqVplA5/Nh2J/9e2+D4TSpzMHSks5ViUGHiNStlZ1ZD1hRkcrIOwHp60CdIDZI6FO5slTLvnDxc1XeeurUqTRr1oxTp04xaNAgrrrqqlL7Z86cSd++fVm/fj39+/cvcq7KMm3aNEaOHMnEiRNxOBzk5Hj+EgtJel0Ba8v0oTm0JTC21AGGdm6OXcBR5q9YgSe+Xs8jl/UOaSfLA22AkuG5NKwGxh4RkebA08AAEXnQ5YjdA1wIxIhIV1V9s+x5qjoZmAyQnJysKSkppfZv3LiRxo1d6ZmnOedEemnOef/990vNN9dddx0iQqNGjWjcuDE//fQTAwYMYPfu3SQmJhIT4znaOXnyZC655JJS803Rz1aCqKgoBgwYUPnP6WX2vP9fnCdP0r/M7/+02BsNU6eBK4LkdrS6b3mD7ude6fPFHac64X1o36E9KYkpPr2Xv7kyZyNvzi/fHiLXAd/uKqR9+/Y8cElPf5v1AHArsBa4A/gWeNvfRoQKh7OsRTq/12Ad2Q7Z6dDxbP/e12CoAZWlCG4QkXLNPUXk91gRrPpF7nHLuQLrPfe4Vy77yiuv0L9/f4YOHcrevXvZunUrAOeffz6JiYmcOHGCBx98sNx5Tz/9NImJibRu3RqAQYMG8c477/DYY4+xdu1ajw87IUsQqHrVJZI6xPLxnWcxqGMsEfbS0cHVace5ZtLCYJFODgiqekRV73Q1Qn/WNfaKK1X6Tk/OlU/wwZxTL+ab2sq0g+VA3fId2KPK7FAreu5jbGLDJrY6JXLhxq0o2LKJ54fzQKQqq6pTVd9S1d+p6tWuz3VrIdmLHMrOBQLQA2v3L9Z7h3P8e1+DoQZUFsH6A/C5iNyCpSAIkAw0AK70tWF+pRqRJvYutSRBHflWTdBVb9d6FXPevHnMmTOHRYsWER0dTUpKCrm51sQ1d+5c4uLiio7t1asXq1evxul0YrPZmDhxIhMnTqRRo0aAlUo4f/58vvnmG2666Sb+9Kc/ccMN5fzj0OXse2FmmazUOY/Bzd8GxJxQJ6lDLJ/ceRbTluzh7zNKR1IcTvjHrI18fGf5lK8QZR/QrsR2WwLdLD0Ac069mW9qK9Pupt1gaN4FMtaXHs8+5IWLV02YhNU5kQs3D1zSk+G9W3L1GwvL/a86dsp/qcoispZK/rWoaj+fGxGCHMoKUJPhXb9Cw3iI6+bf+xoMNaAymfZ9qjoEeAKrP8Qu4AlVHayqgX04CQTtBlsFlRd4r7Dy+PHjxMbGEh0dzaZNm1i8eHGFx3bt2pXk5GQeeughHA4rbSU3N7eo+e3u3btJSEjgtttuY/z48axYUcdaeHhS9Tps0gRrS0X1WNsysv1siU9ZBnQTkU4iEgFcC3xVxTmBx8tzTr2Zb7wh0+4mvGwECzi223vXrwS7zV5nGg17YvGOIx49G6fChM/W+CuKfilwGfCd63Wd6zULK03Q4IGAOFju+quOZ3v3b9xg8BGVRbAAUNWfgJ/8YEvw026wV3PvL7roIt5880169uxJ9+7dGTp0aKXHv/322/z1r3+la9euNG/enAYNGvD8888D1ur0Cy+8QHh4OI0aNeL999/3mp1BQe6J8mONW/rfjjrG0M7NrQX/Mk86x3MLQ1LwQkSmAylAnIikAY+q6hQRuRv4HrADU1V1fSWXCR68OOfUq/nGW9ldzbrAvuWlxw6uhR8etRqh+5AwCauTKYJuhnZujk0sh6os2zKyuXbyIj68/UyfzkGquhtARIarasliwQkisgKrNstQhjVpxxBg95Ec4ht7WITwBZm74MQ+6GDqrwyhQZUOlsF3REZGMmvWrHLju3bt8nh8kyZNmDRpksd9N954IzfeeKM3zQsuTDNBn5DUIZbWMQ3Yd+xUqXGHU/nDB8v542/OYNyQ9gGy7vRRVQ+N00BVv6Wer0jXm/nGGzLtbnIOex5f+7HvHSxb3U0RBGvu6dsmhtVpnmsLCxzK4h1H/LXIIyJytqr+6to4i8pr1Osty3dn8vnK/Shw/ZQlfDB+qO/+H+1dClu+h7aDYOd8ayzKs/COwRBsmAnEEBo0blV+zL2SbKgVvVt7VqY7eCKPv89Yy7Qle/xskcFQCzyFZGtKzys8j0d76JPlZew2e53sg1WSMYMqX7zxY2+sW4H/iMguEdkN/Ae4paYXE5F2IjJXRDaIyHoRudc13kxEfhCRra73WNf4Va7jFrhUTBGRLiLykRd+Nq+yeMcRHK6wY0Gh03cNovcuhXcvgQX/hOljYLGr/eBXfzQCV4aQwESwDKHB2ffC5m+LVdXc+GElua5zx3ld+GlzBoVlddtdzFp3IKSiWIb6jeBFByv5Jlj7EexeWHq8Mrl8L1HXI1hA0bzy0bI9HiNZk+dvL3Wcr1DV5UB/EYlxbddWsrMQ+LOqrhCRxsByEfkBuAn4UVWfE5EHsFIQJ2C1gRgEjAbGAa8CT2E1Ng8qhnZubgWJgfAwG0M7+2ixYdcCcHj49+/It/bV0z6YBQUFpKWlFQkUGfxHVFQUbdu2JTw8vFrHGwfLEBq0Gwy//Xd5JcEwP+V/12GSOsTy0e1n8tmKNI/Rqt6tfP8waTB4DW/ItJek75jyDtbuXyH1XcsB8xF2qdsiF27GDWlP95aNueqNheX27TqSU6Ry6gsnS0R+r6r/E5E/lRkHQFVfrMl1VfUAcMD1OUtENmL147sCq0YU4D1gHpaD5QQigWigQETOBQ6q6taa3N+X9GjZGAXO7RrHfcPP8F16YMdzi6PR9ghAwFlofe54rm/uGQKkpaXRuHFjOnbsWPTv1OB7VJUjR46QlpZGp06dqnWOcbAMoUNCLyhaO3ORuctKF6inq1newv0l6cnBOpFXt1fRDXUMb8m0uzm4yvP44v/41MGqDxEsN1WlmX20bI+volgNXe8+a+QmIh2BAcASIMHlfAEcBBJcn58F5gD7gd8Dn2CpnVZ23duB2wESEhKYN2+et033yM7jltPfv1EWWTtXM2+nj26kTs6RME426sj2rrcC0PTYOo417cOJ7TmwfZ6PbhzcxMTE0Lx5c7Kz65TSb0gQERHBsWPHqv23ZhwsQ+iwawHlnpzUaTX+NA5WranoIefT1L3k5BXy0rUDPO43GIIKb6/qVtT3ykvN5iuirqsIlqSqWqvIMN+Ui6vqJNe7T/LMRaQR8Blwn6qeKBlxUFUVEXV9/gH4wXXODViCPGeIyF+ATOBeVc0pY/tkYDJAcnKypqSk+OJHKMfRFWmwaDWjLhhC13gfNhjP3AU/5xOT8gcGJt3ku/uEGBs3bqRJE5NVEiiioqIYMKB6z0JG5MIQOnQ8FyuCVYbl75uiVy9Q0UNOvkP5YtV+7vtwpZ8tMhhqiLdqsAAatfA8HuZbAQa7zV5vIlgV9eNzk1fo9FdfLK8hIuFYztUHqvq5azhdRFq59rcCMsqcE41Vp/U68DhwI/ALVm+uoGBbRjZhNqFD84ZVH1wb0l2dNBL6+PY+BoOPCHoHS0Q6i8gUEfm0zHhDEUkVkUsDZZuv6NixI3379qVfv36MGDGCgwcPApCdnc1dd91Fly5dGDhwIElJSbz11lsBttaPtBsMZ1xUflwdRsbdC1T1kDNnY7qfLDH4kzo333hTRRCg/zgQe/nx4/t8urATZgur8yqCboZ2bk5UuM3T8hkAq9OOM2byopBxssQKVU0BNpap4/oKy2nC9f5lmVP/CryiqgVAA6yUDSdWbVZQsDUjm45xDQm3+/jx8eA6QKBFD9/ex3BaHDp0iHPOOYc+ffrwxRdfFI1fccUV7N+/3+M5mzdvJiUlhcTERHr27Mntt9/uJ2trxk033cSnn35a9YFV4NO/EBGZKiIZIrKuzPhFIrJZRLa5lHQqRFV3qOqtHnZNAD72pr3BxNy5c1mzZg3Jyck888wzAIwfP57Y2Fi2bt3KihUr+O677zh69Gi5cwsL6/Cq57l/8jzewPeyyXUd90NORZNCTr4jZB5wDKdHnZpvvNkHC1wCOx60DtQBv77svfuUIUzC6oXIBVg1oB+MH8rYIe2xV+BlFTqUST9v969hNeds4HrgAhFZ5XpdAjwHDBeRrcCFrm0ARKQ1MFhVv3ANvQosA+4EpvnT+MrYnpFN1xaNfH+j9HXQrBNE+uFehmozffp07rzzTpYuXcpLL70EwNdff82AAQNo3bq1x3P++Mc/cv/997Nq1So2btzIPffcU+37ORyhOwf6OoL1LlAq5CAidqzw98VAL2CsiPQSkb4iMrPMK97TRUVkOLCBMuF1X7MqYxVvr32bVRmrvHbNUaNGkZSURO/evZk8eXK5/cOGDWPbtm1s376dpUuX8tRTT2GzWf/bWrRowYQJEwCYN28e5557Lpdffjm9evXymn1BR7vB0M1DFGvbbP/bUsdwP+T8eWR3zutWPi3KqfDE1+uNk+VHvD3n1If5RrwdwQJLzKKhh6+jtFTv3qcE9SlFEKz555kr+3JxXw89D12s2OubuUdE7hWRJmIxRURWiMiIml5PVX9RVVHVfqqa6Hp9q6pHVPU3qtpNVS9U1aMlztmvqr8tsf2JqvZW1bNVtYJCQP+SV+hg99EcuiX4w8Fab9IDvcTy3Zm8PnebV767w8PDycnJIS8vD7vdTmFhIS+99BJ/+9vfKjznwIEDtG3btmi7b9++gOU8/eUvf6FPnz7069ePV199FbCyKiZMmMDAgQP55JNPmD17NmeeeSYDBw7kd7/7XZHAx/LlyznvvPNISkpi5MiRHDhg6cekpKQwYcIEBg8ezBlnnMGCBQuK7vfXv/6VQYMG0a9fPyZNmgRYCoF333033bt358ILLyQjwzuuhU9FLlR1vktBpySDgW2qugNARD4ErlDVZ4HqpvulYKn/9AJOici3qmUbJFWffyz9B5uObqr0mOz8bDZnbkZRBKF7bHcaRVQ8yfRo1oMJgydUee+pU6fSrFkzTp06xaBBg7jqqqtK7Z85cyZ9+/Zl/fr19O/fv+hhxxMrVqxg3bp11ZaQDFnOuR+2fld6bNM3PpdNrg8kdYh1vQ7z89by3+nuVJ2Pbj/Td/K89YBAzTn1Y74pozTqDfYuhZOHy49np/tMxTTMVn9ELkrSv20MM9cc8LjvSFY+y3dn+mLuuUVVXxaRkUAsVvTpv4BZuSvBrsM5OJxK13gfO1g7foaj26HD2b69T4jz+Nfr2bD/RKXHZOUWsOlgFk4Fm1gy+42jKu7j1Kt1Ex69rHeF+8eNG8e4ceOYPHky//jHP/jPf/7D9ddfT3R0xVms999/PxdccAFnnXUWI0aM4Oabb6Zp06ZMnjyZXbt2sWrVKsLCwkplSDRv3pwVK1Zw+PBhRo8ezZw5c2jYsCH/+Mc/ePHFF3nwwQe55557+PLLL2nRogUfffQREydOZOrUqYCVWbF06VK+/fZbHn/8cebMmcOUKVOIiYlh2bJl5OXlcfbZZzNixAhWrlzJ5s2b2bBhA+np6fTq1Ytbbqlxn/EiAqEi2AbYW2I7DRhS0cGuruZPAwNE5EFVfVZVJ7r23QQc9uRcVSVhGhMTQ1ZWFgD5+flVhiGP5x0v6q2iKMfzjtPA3qDC4/Pz84uuXxkvvPACM2fOBGDv3r2sWrUKVeW8887DbrfTu3dvJkyYwMKFCyksLCy65gsvvMAXX3zBoUOH2LJlCzk5OSQlJREXF1fhfXNzc/0m5erGvm8fccD69evJi4r0yjVb7f+e7mXGFDj50wukZnf0yj3K4nA42Ju2l3kn5/nk+sHGvD0FFe4rdCivfb2EG/uYHmS+JKsgq9Sck1WQVamDVR1eeeUVZsyYAVjzzdatVpud888/H7vdTr9+/XjqqaeYP39+qfOefvppPvnkEzIyMory7AcPHhyEzhXeVxEEV41nBWt4Pmp6ahc7+Y7K6yLrIjsOn6x0/+IdR3zhYLn/0VwC/FdV14tpMlSObRlW5KCLL1ME9y6FD35nfV7zIQy83qgE14ITuYU4XetNTrW2K3OwqiImJoZvvvkGgMzMTJ577jlmzJjBbbfdRmZmJn/+858588wzS51z8803M3LkSL777ju+/PJLJk2axOrVq5kzZw533nknYWGWK9KsWbOic8aMGQPA4sWL2bBhA2efbTnb+fn5nHnmmWzevJl169YxfPhwwHpGa9WqOPo9evRoAJKSkti1axcAs2fPZs2aNUX1VcePH2fr1q3Mnz+fsWPHYrfbad26NRdccEGNfz8lCXqZdlU9gpWD7Gnfu5WcV6mE6caNG2nc2JIYffich6u0Y1XGKm6bfRsFzgLCbeE8f97zJMYnVu+HqIB58+axYMEClixZQnR0NCkpKdjtdkSEn3/+mbi4uKJjk5KSePDBB2nYsCE2m40nnniCJ554gkaNGtG4cWOio6Np0qRJ0c/kidORl/QWuVu2sBPo3bs3TbwlI/vf8nUPAjRq2BBfSdXaP7DTrm07Ugb55vrBxpQpSwAPK/YuWrdpQ0pKX/8ZVMeoTnS77Jzz3LnP1WrOmTdvHnPmzGHRokVF801ubi5g1WCVnG969erF6tWrcTqd2Gw2Jk6cyMSJE2nUqPjBqmFDH6uI1QL1dopgx3PBZodyNVHis6anYbb6U4NVElslfk24XRja2Sf1tstFZDbQCXhQRBpToUddf9makYWIjx2sXQvAvbDgdPhsAaMuUFmkyc3y3Zlc9/ZiCgqdhIfZePnaAV5boHjyySeZOHEi06dP55xzzuHqq69m9OjRfP/99+WObd26Nbfccgu33HILffr0Yd26dR6uWIz7+0VVGT58ONOnTy+1f+3atfTu3ZtFixZ5PD8y0lrQd6cxuq/16quvMnLkyFLHfvvtt9X7gU+TQKgI7gPaldhu6xoLahLjE3lrxFvcPeBu3hrxVq2dK7C859jYWKKjo9m0aROLFy+u8NiuXbuSnJzMQw89VBRty83N9f6DRCjQsp/n8TMu9q8ddZjerSrvs9G7dYyfLKm/eHvOqTfzjbcbDYP1gHeJB6ELH8Y4wqT+qAiWZPTAtkRU1PfKd0GlW4EHgEGuflMRwM2+ulmosi0jm7axDWgQ4UFV01t0PLf4/7M9wmcLGPUFd231n0Z054PxQ73mXG3dupW0tDRSUlLIycnBZrMhIpw6darcsd999x0FBVZWzMGDBzly5Aht2rRh+PDhTJo0qcgB8iSiNHToUH799Ve2bdsGwMmTJ9myZQvdu3fn0KFDRQ5WQUEB69evr9TmkSNH8sYbbxTZsmXLFk6ePMmwYcP46KOPcDgcHDhwgLlz59b8F1OCQDhYy4BuItJJRCKwupV/FQA7TpvE+ETG9x3vFecK4KKLLqKwsJCePXvywAMPMHTo0EqPf/vttzly5EjRw8/w4cN5/vnnvWJLSBFVwcN/mm97YanXn9qCl8YNKk8heOXHLUxbssdP1tRfvDnn1Jv5xhciF2DVd6Y8WH7cRy0iwmxhFDgrTtWtqyR1iGX6bUO5dlDbcvscDmeFDdFriWLVdP/Rtd0QMDnQZdiWkU03XzYXBmje1fr77TQMbvzKRK+8QFKHWP5wflevptZOnDiRp59+GoCxY8fyxhtvMGjQIO69995yx86ePZs+ffrQv39/Ro4cyQsvvEDLli0ZP3487du3p1+/fvTv359p08qLZbZo0YJ3332XsWPH0q9fP84880w2bdpEREQEn376KRMmTKB///4kJiaycOHCSm0eP348vXr1YuDAgfTp04c77riDwsJCrrzySrp160avXr244YYbyqU41hSfpgiKyHQsQYo4EUkDHlXVKSJyN/A9YAemqmrlbmcdJTIyklmzZpUbd+eLlqVJkyZFqidlSUlJ8Vl6XNDR8VywhUFZha3dC+GHR2H4416/pfhyqToIcUu25xU4PbqVB0/k8fcZawEYN6S9f40z1Ih6M994W6a9JH1/B/OeLd62R/pshd1us9fLCBZYD4QD2zflw2VpRWMChIfZfJUi+B+slMALgCeALKwmwYN8cbNQxOFUdhw+ybAzKmi87S22zQEULnwM2iT59l6GGvPxx8VdkuLj4yt1bl588UVefLF8BkBYWJjHfWW/ky644AKWLVtW7vzExMRy9cJAKa2BuLi4ouvZbDaeeeaZolYkJXnttdcqtL+m+DSCpapjVbWVqoaraltVneIa/1ZVz1DVLqr6tC9tMNRB2g2GgTd43remzrZG8ysl+9JU5lp+tMxEsQzBhU9k2t2kl1gLFBtc9JzPVtjtUr9k2suyYs+xUts9WjX2aopTGYao6h+AXABVzcRKEzS42Hs0h/xCp+97YG35zmqJ0Mq/9eIGg7cJRIqgwVB7+o/1PB5tGg57C3dfmu4tK04JSWhismgMwYYPHaz9K4o/q8Ipn6SrAfVX5MJN2VTABuF2X7aGKHD16FQAEWmBEbkohVtBsKsve2A5Cq0IVrcRUEmLCIMhFDD/gg2hSUWrxu2S/WtHPcDh9PywKsAd53XxrzF+QEQ6u5qNflpiLEVEFojImyKSEjjrDFXiS3XtiJKLDQq5lfegqQ1htrB6HcGKjS4dQEo/kevLJuevADOAeBF5GvgFKJ9HVI/Z6nawfNkDa+8SyD0OZ4ys+liDIcgxDpYhNEl91/N4pFG38ybTluwp+mIti80WfHVpIjJVRDJEZF2Z8YtEZLOIbBORByq7hqruUNVbyw4D2ViF72nlzzIEEz4TpCnr8Cx6zerb4wPCpH42GnaTmVO6B9i+Y7lc9/ZinzhZqvoB8DfgWeAAMEpVP/H6jUKYbRnZxDeOpEkteihVyZbvwBYOXc733T0MBj9hHCxDaLLxS8/jB9f41446zqx1Byrc53SqrxS9asO7wEUlB1ypP68DF2MphY0VkV4i0ldEZpZ5xVdw3QWqejEwAfC+iorBe/hCpt1Nl/MtgR036vSZimB9FrmAYqGdkss4BYW+UREUkaHAPlV9XVVfA/aJyBCv3yiE2ZaRRTdfpgemvgtLJ1sqwRkbfXcfg8FPBH2jYYPBIz2vgO0/eR43eI2L+7RiwVbPDYftNp81/awxqjpfRDqWGR4MbFPVHQAi8iFwhao+C1xazeu66zEygUhPx4jI7cDtAAkJCaWUjABiYmLIysqq3g9SR8jNzS33e/A1TTIyCM/J8dl9W3W9nTO2/AcAp4Sx+mhDTvjgXgeOHiAvP8/vv79g4i8DI/h1XwFz0yxH0y4QeWw38+Z5PYj8BjCwxHa2h7F6i6qy/dBJrhrYxjc3SH0XZrrkvQtz4Z1L4OZvjUS7IaQxDlaAefrpp5k2bRp2ux2bzcakSZMYPHgwd9xxB7/88gt2u53JkyeX0uXv2LEjjRs3RkSIjY3l/fffp0OHDgH8KQJA8k2w7QfYNLN4TGyQ0CtgJtVF3BLsk+dvZ9eRnFL7fFnq4mXaAHtLbKcBFa5Oi0hz4GlggIg8qKrPishoYCTQFPCo56qqk4HJAMnJyVpWxnzjxo00buzjHjJV4O/5JioqigED/KsGtu/bbzm1f78PZeRT4Llp0LwL9oueY6CPHgJXpK5g8abFwSuH7wdSgPHAiH//TF6BkxfHJPpK6EK0RBdtVXWKiHk+cnHwRC7ZeYW+q79aXab/kbPAigwbByvoOHToEFdeeSXHjh3jqaeeYtSoUQBcccUVvPHGG7Ru3brcOY899hhvvfUWLVq0ID8/n4cffpixYysQKnPRqFEjsrM9lyd4g5SUFP75z3+SnOy7un2TIhhAFi1axMyZM1mxYgVr1qxhzpw5tGvXjl9++YWtW7eyfv16li5dSufOncudO3fuXNasWUNKSgpPPfVUAKwPAs4u39DOV+k69ZlxQ9rzr2sSCbeX9qgcwZkiWGtU9Yiq3ulqI/Gsa+xzVb1DVceo6rwAm1gj6st841OZdjfhDSC+l08fAO22+i3TXpK4RpHEN4n0pYrgDhH5o4iEu173Ajt8dbNQY2u6W+DCRwtEzc8ovW0L91l/OUPtmD59OnfeeSdLly7lpZdeAuDrr79mwIABHp0rN/fffz+rVq3iyy+/5I477qCgoO43UTcO1mmQs3IlhydNJmflSq9c78CBA8TFxREZaWUcxcXF0bp1ayIiIkhPT6egoIAGDRqQkJBQ4TXOPPNM9u3b5xV7QpMSD/22MDMp+4ikDrH8ZUT3UmPhdp81/fQ2+4B2JbbbusaCHm/OOfVnvvGDg+V0wIE1PhO4AJeKoBaivv5ZQoDcAgd7jub4UkXwTuAsrHnBHeG+3Vc3CzW2+VpBMDoWxA6tB0KPS016oLfZuxQW/Msr81V4eDg5OTnk5eVht9spLCzkpZde4m9/+1u1zu/WrRvR0dFkZlp/yy+88AKDBg2iX79+PProo+WOz87O5je/+Q0DBw6kb9++fPmlVX+/bNky+vXrR25uLidPnqR3796sW7eOkydPcssttzB48GAGDBhQdPypU6e49tpr6dmzJ1deeSWnTp2q9e+iKkwIHDj4zDPkbdxU6TGO7GzyNm2yvrhFiOzRA3ujiiebyJ49aPn3v1d6zREjRvDEE09wxhlncOGFFzJmzBjOO+88EhISyMrK4qabbuKDDz6wVmQr4LvvvisK0dY7di2gVDV7z8vNpOxDIsNKr8fcNqyTL1eUvckyoJuIdMJ6gLoWGBdIgwIx59Sb+cbXuat7l8LJQ3AyA967HG78yifzjl3sADjVWfS5PrJ8dyar9h7DqXDd24t90mxYVTOw5gWDB7YdyiamQThxjXzUezl9HST0htvn+ub6dZVZD8DBtZUfk3fC+v2q01VG0Qcim1R8fMu+cPFzFe4eN24c48aNY/LkyfzjH//gP//5D9dffz3R0dHVMnnFihV069aN+Ph4Zs+ezdatW1m6dCmqyuWXX878+fMZNmxY0fFRUVHMmDGDJk2acPjwYYYOHcrll1/OoEGDuPzyy3nooYc4deoUv//97+nTpw9///vfueCCC5g6dSrHjh1j8ODBXHjhhUyaNIno6Gg2btzImjVrGDjQ9+WVJoJVTZwnThSviqpa27WkUaNGLF++nMmTJ9OiRQvGjBnDu+++y9VXX838+fOJjo7m/vvvB+APf/gDM2cW1xudf/75tGnThlmzZlWZy1pn6XgupSJYPX4bMFPqAxsPlv43fzwn+NKXRGQ6sAjoLiJpInKrqhYCdwPfAxuBj1V1fSDtrA7ennPq03zjM5l2KL2w48j3WVpymEutsD5LtYPVcNjdis/bKoIi8jfX+6si8krZl9duFOJsS8+mW3yjShdfasXBddaDvcH75B63nCuw3nOP1+pyMTExfPPNN6SmpjJw4EC+/vprrr76am677TauvvpqFi1a5PG8f//73/Tu3ZshQ4YwceJEAGbPns3s2bMZMGAAAwcOZNOmTWzdurXUearK3//+d/r168eFF17Ivn37SE9PB+CRRx7hhx9+IDU1tSiCNnv2bJ577jkSExNJSUkhNzeXPXv2MH/+fH7/+98D0K9fP/r161er30N1MBEsqDLSBFaqzp6bb0ELCpDwcFr/8wWivVC8bbfbSUlJISUlhb59+zJlyhQOHz5Mp06dmDRpEldddRWPP/44y5Yt4/nnny86b+7cuTRt2pTrrruORx99lBdffLHWtoQmJR6kioTeDL4gwl56Fb1hRPCtqquqx6d/Vf0W+NbP5lRIoOacejHf+FKmHayFHXuk1RPLHuGztOQwl8aCw+mA4PtT8xtDOzcnMsxGgcNJeJjX05LdeuCp3rxoXWPjwRN0at6Q5bszvZ+1kJ1hRYMT+nj3uvWBSiJNRexdakXaHfnWfHXV216LuD/55JNMnDiR6dOnc84553D11VczevRovv/++3LH3n///fzlL3/hq6++4tZbb2X79u2oKg8++CB33HFHhff44IMPOHToEMuXLyc8PJyOHTuSm5sLwJEjR8jOzqagoIDc3FwaNmyIqvLZZ5/RvXv3Cq/pL0wEq5pEDxhA+3em0uKPf6T9O1O94lxt3ry5lLe+atUqOnfujKoyd+7cIkWvl19+mYEDB9KwYcNS54eFhfHSSy/x/vvvc/To0VrbE3LsWmCFvN2kme9IX1I2PSTzVN0vUg0k3p5z6st8U3j0KM6sLK/Vypaj3WC4aSZcMNFn6YFgiVwAFDjr999ZUodYpt02lD+P6O719EBV/dr1/p6nl9duFMLM2ZBOVm4ha/cd902jZ3eKW0vjYPmEdoOtecrL89XWrVtJS0sjJSWFnJwcbDYbIlJlbdPll19OcnIy7733HiNHjmTq1KlFaoH79u0jIyOj1PHHjx8nPj6e8PBw5s6dy+7du4v23XHHHTz55JNcd911TJgwAYCRI0fy6quvFtWurnR9DwwbNoxp0yy1ynXr1rFmje97ppoI1mkQPWCAVxwrN9nZ2dxzzz0cO3aMsLAwunbtyuTJk7n55pv54x//SE5ODtHR0bz22ms8//zzfPrpp1x99dWlrtGqVSvGjh3L66+/zsMPP+w120IC90qyI8+KXrVKDLRFdZpzurXgjZ+3k1/oxKnQ3ZdNJw2Ad+ec+jDf5KxcyclffgGHgz033+K1xbBytBvs83pPd4pgfW427CapQ6xP6j1F5GsqiXeq6uVev2kIsevwSf722WrA+iW5UzS9+v/C7WCZCJbv8MF8NXHiRJ5++mkAxo4dy6hRo3juued44oknqjz3kUceYdy4cWzcuJGNGzcWtQVp1KgR//vf/4iPjy869rrrruOyyy6jb9++JCcn06NHDwDef/99wsPDGTduHA6Hg7POOouffvqJhx9+mPvuu49+/frhdDrp1KkTM2fO5K677uLmm2+mZ8+e9OzZk6SkJK/+PjxhHKwAkpSUxMKFC8uNx8XFlRsfN664Jn/Xrl2l9r366qs+sS/oca/MLH0L1n4MCT0DbVGdJqlDLB+MH8qMlfv43+LddIxrWPVJhqChPsw3OUuXgdNKFdaCAnKWLvONg+UH3MIWRqrdp/zT9T4aaAn8z7U9Fkj3xQ1F5CLgZazEz7dV9TkR+QDoC8xU1b+7jnsIWKeqX/jCjqpYuO0wd32wAofTSYTdhsPpkxRNS4ChSVuIbubd6xp8yscff1z0OT4+3uN3i5vHHnus1HZSUhKbN28G4N577+Xee8u33HFHteLi4jzWdXXs2JEbbrgBsFLflyxZUrRv0qRJ5Y5v0KABH374YSU/kfcxDpYhtGk3GPYsshysA2ugVf9AW1SnSeoQS2SYjf8t3s1Xq/YT0yAiVJQEDfWA6MGDkMjIorq16MGDAm1SjTlw8gAAaw+t5TcdfhNga+omqvozgIj8S1VLdhz9WkS8nnMuInbgdWA4lhz8MhH5Fjilqv1E5AcRiQGigSGq6rWmc0s+/icNtn/LqS6XENOhP5kbfyK25wUARZ+7D/oNm5f9yPpF3zItoz3xzQcw5cZB5OxYROaGn4jtdQE92jeFPUtg1y/Q8ZziyMjepVWPubejm0HOUet96w8Q2cjaZ1SADXUI42AZQpu9S+En13fQN3+GFt3NJO1j3D1Rvly1n+/WH/SJbLLBUBPcdWs5S5cRPXhQyEavVmWs4r31VgnQ3xb8jSkNppAYnxhYo+o2DUWks6ruAHC1dPBFiH4wsK3EfT4Efgs0EBEbEA44gCeA8k2BasiST/7F4PVPWh0MNiyHDa4du/5TfNCu/8As6IH1uioCyAJeK3GhXa/7Rioo9xi8e6lV22i+vw11hHrtYKmq72RHg4w626xy1wJwuFJonAXWtpmgfcqafccAH+bk12HMnON7vF0rGwhS01Mt9UCsFMHU9FTjYPmW+4F5IrIDq/dHB6BiabOa0wbYW2Lb3dT4ELAC+C/QFbCp6orKLiQit+NqhpyQkMC8efMqPLbxxi+KPrv/LEXAqdYP6/68T+JpoxnYXNvbovoS10CIPbYGARThVFQCDXIPFm1nNrXk1Use42ms9HmUegdQRz47f3qfPR1yKv0F1ndiYmI4ceJEvfkeCSZUldzc3Er/1kpSbx2sqKgojhw5QvPmzev8P1RV5ciRI0RFRQXaFO/T8Vywh1kSpLZwn0kmG4pJbNsUAJvgm5z8OoqZcwzVJTkhmQh7BAXOAsJt4SQnJFd9kqHGqOp3ItINK3gDsElV8/x4//vcn13CG3eIyESgP/CDqr7l4ZzJwGSA5ORkTUlJqfD6Sw6NgnVripyrAuzYVHG4hKTt6qSAMPb3vpO4dc8RroUUEIbzwsdo1rJJkcy32COIvvBB+O6Bou1mV/3LumiJYzyNFZ1XmIfgBAQpoS8i9gg6X3ADnc0CaaXs3Lnz/9u711g5yjqO499fodgWhCZcqnIKbZFA2pAerk3k1ki5KGhBQShvQJFQBQ0RE/GKb1QiMSGAESkgSABpGoIIKLzApAgEaEmxFyWWWziNSkVFLo0F/Pti55ThdM9lu3N5Zvf3SSY9+8wzu7/zzOx/+5zZ2WXr1q198TqSkuHXtOnTp3PoBP+A17cTrIGBAYaGhti8eXPdUSoxZcoUBgYG6o5RvJlHwQnfh4e+BZ+4ymevKjD3I61vgT/1kA9z/tGzffZqglxzbKIG9xlk2UnLWPX3VRwx4wifvarG4cAsWv8vmi+JiPhlwY+xCZiZuz2QtQEgaTGwGtgNOCAiPifpQUm3R8QOn9pZcNZlPAFM3Xg/Wz56KnvMmr/tmipg288LjlzEn3PrDj5yUesOzru39e6QWce2XmNnzH3/7XZ9xtpu6p6w5dXWv39bAwjmL/Hr9wT02+tISjp9TevbCdbkyZOZPXt23TGsCHtnXyg3Y269OfrMSfM+5MlVB1xzrBOD+wx6YlURSbcBBwBraF0DBa13sBU9wXoKODC7xmsTcA5wbpZhMnAprWuyDuS9j4/fCdgF6Oq9cwvOugy47L2G4cnTiJ8PPnLR+9fB9h/z3e5jvyfSVsHXG/Q6v440R99OsMzMzMyAI4C5UfKFgxHxjqRLgAdpTZxujoj12eqLgVsj4i1JfwSmSVoLPBAR/y4zl5kVzxMsMzMz62fraH0P1l/LfqCIeIA2n8UXEVfnfg5a38VlZg3lCZZVpPmfYhg98DuYmdl29gI2SHoS2PbhFhHx6foimVmTqWc/vjtH0mbgJVpF9B8VPewewGsVbDuRvmP1GW1du/aRbSNvN2V8O91+vL5ljW+7tqaMcX7b/SNi72IiNUNWc96kefuqiP69+HxIqd6M16cfa3pX9UbS8e3ah7+IOHW5/+PUrcrjpVOpZks1F6SbLbVc7WtORPTNAqyq8LFuqGLbifQdq89o69q1j2xrc7sR41v0GJc1vk0e4273Ty8svbqv+vH5kFK9KXOMXW+8dLkPKzteeiVbqrlSzpZqrpGL3yJYnt9UtO1E+o7VZ7R17dpHtnXzO3ar28cucozLGt+JPHaZqjqGrXtV7qt+fD6kVG/G69OPNX2HtpX0Ou3fvy5al0Lt3kUmM+tjffEWwWGSVkWEv7GxJB7f8nmMm8P7qnwe43J5fK0TKR8vqWZLNRekmy3VXCNNqjtAxW6oO0CP8/iWz2PcHN5X5fMYl8vja51I+XhJNVuquSDdbKnmep++OoNlZmZmZmZWpn47g2VmZmZmZlYaT7DMzMzMzMwK4gmWmZmZmZlZQTzBMjMzMzMzK0hfT7AkLZT0iKTrJS2sO0+vkTRJ0g8kXSvpvLrz9CJJx2bH742SHqs7j43O9aZcrjflc72xHSXpdEnLJN0l6aS68wxLuS5L2k/SPZJulnR5AnnmSLpJ0oqx2lLJlrXvKmmVpNOqztRzE6zsQHxF0roR7adIelbSxtyBGsAbwBRgqOqsTdTh+C4GBoC38fhOWCdjHBGPRMRS4D7g1jry9jPXm3K53pTP9cbG0+HzsK2IuCciLgSWAmenkouS6nJB2Q4BVkTEF4BD684TEc9HxAXjtaWSLfMNYHk3+XZYRPTUAhwHHAasy7XtBDwHzAF2AZ4B5gKTsvUzgNvrzt6EpcPxvRy4KOuzou7sTVk6GePc+uXAB+vO3m+L601S4+t6U/IY59a73vTR0uHz8BBaE/D8sk9uu58Ah6WSq6y6XFC2PYHfAw8Dn687T2677eprNzW3rGzAicA5wPnAaVU/b3amx0TESkmzRjQfBWyMiOcBJP0KWBwRG7L1/wI+UF3K5upkfIGXga1Zn3crC9lwHY7xBkn7Aa9FxOvVJjXXm3K53pTP9cbG02Gd+xGw3duxJAm4EvhtRDydSq6cQutyQWP2deCK7L5WAL+oM09ZSsy2ENiV1sRsi6QHIuJ/xaQeX8+9RXAU+9J68R02BOwr6TOSfg7cBlxXS7Le0HZ8gbuBkyVdC6ysI1gPGW2MAS6gi8JrhXO9KZfrTflcb2w8Yx0j7XwFWAScKWlpKrkqrsudjtnvgK9Kuh54se48kvbMshwq6ZujtaWSLSK+HRGXAncAy6qcXAG9dwarExFxN60XZStBRLxF68XYShQRV9SdwcbnelMu15tquN7YjoiIa4Br6s4xUsp1OSLWAWfWnWNYRLxK6xq6MdvqMFaOiLil2jQt/XIGaxMwM3d7IGuzYnh8y+cxbg7vq3J5fMvnMbbxpHqMpJoL0suWWp68lLNNSL9MsJ4CDpQ0W9IutC56u7fmTL3E41s+j3FzeF+Vy+NbPo+xjSfVYyTVXJBettTy5KWcbUJ6boIl6U7gceAgSUOSLoiId4BLgAeBPwHLI2J9nTmbyuNbPo9xc3hflcvjWz6PsY0n1WMk1VwpZkstT1OydUPZRxmamZmZmZlZl3ruDJaZmZmZmVldPMEyMzMzMzMriCdYZmZmZmZmBfEEy8zMzMzMrCCeYJmZmZmZmRXEEywzMzMzM7OCeIJllZH0rqQ1uWWWpIWS7mvTd56khyU9K+kvkr4rSdm68yVtzu5jg6QLq/9tzCx1rjlmVidJMyTdIel5SaslPS7pjNz6qyVtkjQp1+Z60wM8wbIqbYmIwdzyYrtOkqbS+sbuKyPiIGA+8DHgy7lud0XEILAQ+KGkGaUmN7Mmcs0xs1pkf6C5B1gZEXMi4nDgHGAgWz8JOAN4GTh+xOauNw3nCZal6Fzg0Yh4CCAi3qL1jd6Xj+wYEa8AzwH7SzpL0jpJz0haWWliM2sy1xwzK9rHga0Rcf1wQ0S8FBHXZjcXAuuBnwFL2t1Bvt6UG9WKtnPdAayvTJW0Jvv5hYg4Y5R+84DV+YaIeE7SbpJ2z7dLmgPMATYCNwEnR8QmSdMLTW5mTeSaY2Z1mQc8Pcb6JcCdwK9pnaWaHBFv5zuMqDfWIJ5gWZW2ZKe8i3C2pGOA/wIXRcQ/JT0K3CJpOXB3QY9jZs3lmmNmSZD0U+AYYCtwNPBJ4GsR8bqkJ4CTgeHrQ7erN3Vkth3nCZalaANwXL4h+yvOGxHxn+y687si4pJ8n4hYKmkBcCqwWtLhEfFqVaHNrLFcc8ysaOuBzw7fiIiLJe0FrKI1mZoOrM3qyzRgC+9NsLarN9YsvgbLUnQ7cIykRbDtAvRrgB+PtZGkAyLiiYj4HrAZmFl6UjPrBa45Zla0h4Epkr6Ua5uW/bsE+GJEzIqIWcBs4ERJ07Ce4AmWpeAESUPDCzAILAa+I+lZYC3wFHDdOPdzlaS1ktYBjwHPlBnazBrLNcfMShURAZwOHC/pBUlPArcCVwCnAPfn+r4J/AH4VA1RrQRq7X8zMzMzMzPrls9gmZmZmZmZFcQTLDMzMzMzs4J4gmVmZmZmZlYQT7DMzMzMzMwK4gmWmZmZmZlZQTzBMjMzMzMzK4gnWGZmZmZmZgX5Py/+S1252AAMAAAAAElFTkSuQmCC",
- "text/plain": [
- "<Figure size 864x216 with 3 Axes>"
- ]
- },
- "metadata": {
- "needs_background": "light"
- }
- }
- ],
- "metadata": {}
- }
- ],
- "metadata": {
- "kernelspec": {
- "name": "python3",
- "display_name": "Python 3.9.1 64-bit ('srenet': pyenv)"
- },
- "language_info": {
- "codemirror_mode": {
- "name": "ipython",
- "version": 3
- },
- "file_extension": ".py",
- "mimetype": "text/x-python",
- "name": "python",
- "nbconvert_exporter": "python",
- "pygments_lexer": "ipython3",
- "version": "3.9.1"
- },
- "interpreter": {
- "hash": "bc219ed3a4d84e3940fd3bdf8d6746735ad7a4a1b4666e5401f7a1755a2774f1"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
\ No newline at end of file
diff --git a/pyproject.toml b/pyproject.toml
deleted file mode 100644
index 73a7fd1..0000000
--- a/pyproject.toml
+++ /dev/null
@@ -1,3 +0,0 @@
-[build-system]
-requires = ["setuptools>=40.8.0", "wheel"]
-build-backend = "setuptools.build_meta"
\ No newline at end of file
diff --git a/requirements.txt b/requirements.txt
deleted file mode 100644
index 0731cd3..0000000
--- a/requirements.txt
+++ /dev/null
@@ -1,4 +0,0 @@
-matplotlib==3.4.2
-numpy==1.21.0
-pandas==1.2.5
-scipy==1.7.0
diff --git a/setup.cfg b/setup.cfg
deleted file mode 100644
index 0792c86..0000000
--- a/setup.cfg
+++ /dev/null
@@ -1,2 +0,0 @@
-[metadata]
-license_files = LICENSE
\ No newline at end of file
diff --git a/setup.py b/setup.py
deleted file mode 100644
index 13e34a7..0000000
--- a/setup.py
+++ /dev/null
@@ -1,32 +0,0 @@
-from setuptools import setup, find_packages
-import pathlib
-
-here = pathlib.Path(__file__).parent.resolve()
-long_description = (here / 'README.md').read_text(encoding='utf-8')
-
-PACKAGE_NAME = 'srenet'
-
-setup(
- name=PACKAGE_NAME,
- version='1.0.0',
- description='Screen & Relax : A fast solving procedure for the Elastic-Net problem by safe identification of the solution support.',
- long_description=long_description,
- long_description_content_type='text/markdown',
- author='Theo Guyard',
- author_email='theo.guyard@insa-rennes.fr',
-
- classifiers=[
- 'Development Status :: 1 - Planning',
- 'Intended Audience :: Science/Research',
- 'Topic :: Scientific/Engineering :: Mathematics',
- 'License :: OSI Approved :: MIT License',
- 'Programming Language :: Python :: 3.9',
- ],
-
- keywords='optimization, compressive sensing, screening',
-
- packages=find_packages(),
- include_package_data=True,
- python_requires='>=3.9, <4',
-
-)
\ No newline at end of file
diff --git a/srenet/__init__.py b/srenet/__init__.py
deleted file mode 100644
index c2aad6e..0000000
--- a/srenet/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-__author__ = 'Theo Guyard'
diff --git a/srenet/algebra.py b/srenet/algebra.py
deleted file mode 100644
index 3148448..0000000
--- a/srenet/algebra.py
+++ /dev/null
@@ -1,95 +0,0 @@
-import numpy as np
-
-def coherence(A):
- """
- Return the mutual coherence of a matrix A.
- """
- if A.shape[1] > 0:
- gram = A.T @ A
- n = gram.shape[0]
- for i in range(n):
- gram[i,i] = 0.
- return np.max(np.abs(gram))
- else:
- return 0
-
-def update_data(Ginv,B,Aw,b,yw,tw,i,A,y,S,Sbar,l1,l2,AtA,gram,Aty_l1):
- """
- Update the reduced problem data when a new indice is relaxed.
- """
-
- n = Ginv.shape[0] # Curr. number of relaxed atoms
- Si = np.sum(S[:i]) # Position of i in S
- Sbari = np.sum(Sbar[:i]) # Position of i in Sbar
- perm = np.insert(np.arange(n), Si, -1) # Perm. last <-> ith row/column
-
- # ------------------------------------------------------------------------ #
- # Values computed to update problem data #
- # Value Expression Dimensions #
- # ------------------------------------------------------------------------ #
- g = AtA[S,i] # |S| x 0
- q = gram[i,i] # 1 x 0
- Ginv_g = Ginv @ g # |S| x 0
- d = 1 / (q - (g @ Ginv_g)) # 1 x 0
- Ginv_g_d = d * Ginv_g # |S| x 0
- d_Ginv_g_Ginv_g = np.outer(Ginv_g_d,Ginv_g) # |S| x |S|
- Ginv_11 = Ginv + d_Ginv_g_Ginv_g # |S| x |S|
- Ginv_12 = - Ginv_g_d # |S| x |S|
- Bred = np.delete(B,Sbari,axis=1) # |S|-1 x |S|-1
- d_B_g_AtA = d * (B.T @ g + AtA[Sbar,i]) # |Sbar| x 0
- B2 = np.delete(d_B_g_AtA,Sbari) # |Sbar|-1 x 0
- B1 = Bred + np.outer(Ginv_g,B2.T) # |Sbar|-1 x |S|-1
- b1 = b # |S| x 0
- b1 += Aty_l1[S] @ d_Ginv_g_Ginv_g # |S| x 0
- b1 -= Aty_l1[i] * Ginv_g_d # |S| x 0
- b2 = - Aty_l1[S] @ Ginv_g_d + Aty_l1[i] * d # 1 x 0
- Aw1 = np.outer(A[:,S] @ Ginv_g - A[:,i], B2) # |S| x |S|
- yw1 = A[:,S] @ b1 # m x 0
- yw2 = b2 * A[:,i] # 1 x 0
- tw1 = B1.T @ (l1 + l2*b1) # |S| x 0
- tw2 = - B2 * (l1 + l2*b2) # 1 x 0
-
- # ----------------------------- FLOPS count ------------------------------ #
- # Above computatations
- # Ginv_g : |S|^2/2 (symmetric matrix)
- # d : 2 + |S|
- # Ginv_g_d : |S|
- # d_Ginv_g_Ginv_g : |S|^2/2 (symmetric outer product)
- # Ginv_11 : |S|^2/2 (sum of symmetric matrices)
- # Ginv_12 : |S|
- # d_B_g_AtA : |S||Sbar| + |Sbar|
- # B1 : |S||Sbar|
- # b1 : |S| + 4
- # b2 : |S| + 3
- # Aw1 : |S|^2 + m + |S|
- # yw1 : m|S|
- # yw2 : 1
- # tw1 : 3|S|
- # tw2 : 4
- # TOTAL : (5/2)|S|^2 + 2|S||Sbar| + m|S| + 9|S| + |Sbar| + 14
- # ------------------------------------------------------------------------ #
-
- # Update of Ginv
- Ginv_new = np.vstack((
- np.hstack((np.reshape(Ginv_11,(n,n)), np.reshape(Ginv_12,(n,1)))),
- np.hstack((np.reshape(Ginv_12,(1,n)), np.reshape(d,(1,1))))
- ))
- Ginv_new = Ginv_new[:,perm][perm,:]
-
- # Update of B
- B_new = np.vstack((B1, -B2.T))
- B_new = B_new[perm,:]
-
- # Update of b
- b_new = np.hstack((np.reshape(b1,(n,)), np.reshape(b2,(1,))))
- b_new = b_new[perm]
-
- # Update of Aw
- Aw_new = np.delete(Aw,Sbari,axis=1)
- Aw_new += Aw1
-
- # Update of yw and tw
- yw = y - yw1 - yw2
- tw = l1 + tw1 + tw2
-
- return Ginv_new, B_new, Aw_new, b_new, yw, tw
\ No newline at end of file
diff --git a/srenet/data.py b/srenet/data.py
deleted file mode 100644
index 5b5a991..0000000
--- a/srenet/data.py
+++ /dev/null
@@ -1,88 +0,0 @@
-import numpy as np
-from scipy.fftpack import dct
-
-DIC_GENERATIONS = ['gaussian', 'uniform', 'dct', 'toeplitz']
-OBS_GENERATIONS = ['gaussian', 'posgaussian']
-
-def _normalize(x):
- if len(x.shape) == 2:
- for i in range(x.shape[1]):
- x[:, i] /= np.linalg.norm(x[:, i], 2)
- elif len(x.shape) == 1:
- x /= np.linalg.norm(x,2)
- else:
- raise ValueError("Require x with 1 or 2 dimensions.")
- return x
-
-def _gaussian_dic(m, n):
- A = np.random.randn(n, m).T
- return A
-
-def _uniform_dic(m, n):
- A = np.random.rand(n, m).T
- return A
-
-def _dct_dic(m, n):
- A = dct(np.eye(max(m,n)))
- indices = np.random.permutation(n)
- A = A[indices,:]
- A = A[:n,:m].T
- return A
-
-def _toeplitz_dic(m, n):
- ranget = np.linspace(-10, 10, m)
- offset = 3.
- rangemu = np.linspace(np.min(ranget)+offset, np.max(ranget)-offset, n)
- A = np.zeros((m, n))
- for j in range(n):
- A[:,j] = np.exp(-.5 * (ranget - rangemu[j])**2)
- return A
-
-def _gaussian_obs(m,pos=False):
- y = np.random.randn(m)
- if pos:
- y = np.maximum(y,0)
- return y
-
-def sample_data(m,n,dic_gen,obs_gen):
- """
- Sample a normalized (m,n) matrix and a normalized (m,) vector using given
- generation methods.
-
- Parameters
- ----------
- m: int
- Dimension of y.
- n: int
- Number of columns in A.
- dic_gen: {'gaussian', 'uniform', 'dct', 'toeplitz'}
- Generation method for A.
- dic_gen: {'gaussian', 'posgaussian'}
- Generation method for y.
- """
-
- if dic_gen not in DIC_GENERATIONS:
- raise ValueError(f"Values allowed for parameter dic_gen are " \
- " {DIC_GENERATIONS}.")
- if obs_gen not in OBS_GENERATIONS:
- raise ValueError(f"Values allowed for parameter obs_gen are " \
- " {OBS_GENERATIONS}.")
-
- if dic_gen == 'gaussian':
- A = _gaussian_dic(m, n)
- elif dic_gen == 'uniform':
- A = _uniform_dic(m, n)
- elif dic_gen == 'dct':
- A = _dct_dic(m, n)
- elif dic_gen == 'toeplitz':
- A = _toeplitz_dic(m, n)
-
- if obs_gen == 'gaussian':
- y = _gaussian_obs(m,False)
- elif obs_gen == 'posgaussian':
- y = _gaussian_obs(m,True)
-
- A = _normalize(A)
- y = _normalize(y)
-
- return A, y
diff --git a/srenet/enet.py b/srenet/enet.py
deleted file mode 100644
index 6a232f6..0000000
--- a/srenet/enet.py
+++ /dev/null
@@ -1,19 +0,0 @@
-import numpy as np
-
-def primal_value(x, A, y, l1, l2):
- x = np.maximum(x,0)
- pv = 0.5 * np.linalg.norm(y-A@x,2)**2 + l1 * np.sum(x) + \
- (l2/2) * np.linalg.norm(x,2)**2
- return pv
-
-def dual_value(u, A, y, l1, l2):
- dv = 0.5 * (np.linalg.norm(y,2)**2 - np.linalg.norm(y-u,2)**2 - \
- (1./l2) * np.linalg.norm(np.maximum(A.T @ u-l1,0),2)**2
- )
- return dv
-
-def duality_gap(x, u, A, y, l1, l2):
- x = np.maximum(x,0)
- pv = primal_value(x, A, y, l1, l2)
- dv = dual_value(u, A, y, l1, l2)
- return np.abs(pv-dv)
\ No newline at end of file
diff --git a/srenet/projgrad.py b/srenet/projgrad.py
deleted file mode 100644
index 91b339f..0000000
--- a/srenet/projgrad.py
+++ /dev/null
@@ -1,316 +0,0 @@
-import time
-import numpy as np
-from srenet.algebra import coherence, update_data
-from srenet.enet import primal_value, duality_gap
-from srenet.utils import (
- StopCrit,
- ProjgradParams,
- ProjgradMonitor,
- validate_inputs,
- print_header,
- print_trace,
- print_footer,
-)
-
-def projgrad(x0, A, y, l1, l2, params=ProjgradParams()):
- """
- Solve the non-negative Elastic-Net problem
- min_{x>=0} (1/2) ||y-Ax||_2^2 + l1 ||x||_1 + (l2/2) ||x||_2^2
- using a projected-gradient algorithm. Can be overloaded with screening and
- relaxing tests.
-
- Parameters
- ----------
- x0 : np.ndarray
- Initial point with shape (n,).
- A : np.ndarray
- Dictionary with shape (m,n).
- y : np.ndarray
- Target observation with shape (m,).
- l1 : float, positive
- L1 regularization parameter.
- l2 : float, positive
- L2 regularization parameter.
- params : ProjgradParams
- Optional parameters for the algorithm (see `utils.ProjgradParams`).
-
- Returns
- -------
- xopt : np.ndarray
- Solution of the problem with shape (n,).
- monitor : ProjgradMonitor
- Monitoring values (see `utils.ProjgradMonitor`).
- """
-
- validate_inputs(x0, A, y, l1, l2)
- params.validate()
-
- if params.verbosity:
- print_header(A, l1, l2, params)
-
- monitor = ProjgradMonitor()
-
- ########################### Initialization #############################
-
- # /!\ No FLOPs are counted for the initialization as all values can be
- # /!\ computed off-line
-
- start_time = time.time()
- it = 0
-
- # Precomputed values
- m, n = A.shape
- ninit = n
- sqrt_l2 = np.sqrt(l2)
- Aty = A.T @ y
- Aty_l1 = Aty - l1
- AtA = A.T @ A
- gram = AtA + l2 * np.eye(n)
- L = np.max(np.abs(np.linalg.eigvals(gram)))
-
- # Iterates
- x = np.maximum(x0,0)
- u = y - A @ x
- xold = np.copy(x)
- zold = np.copy(x)
- told = 1
-
- # Best gap recorded so far for screening and relaxing tests
- best_x = np.copy(x)
- best_Atu = np.copy(A.T @ u)
- best_gap = duality_gap(x,u,A,y,l1,l2)
-
- # Screened and relaxed coefficients
- idx = np.arange(n) # Initial indices
- S = np.zeros(n,dtype=bool) # Indices relaxed
- S[params.Sinit] = True # Indices relaxed at initialization
- Sbar = np.invert(S) # Indices not fixed yet
- sc = set() # Screened coefficients
- rc = set() # Relaxed coefficients
- rc.update(idx[S])
-
- # Data of the transformed problem
- Ginv = np.linalg.inv(gram[:,S][S,:]) # (A[:,S].T @ A[:,S] + l2 * I)^{-1}
- b = Ginv @ Aty_l1[S] # Ginv @ Aty_l1[S]
- B = - Ginv @ A[:,S].T @ A[:,Sbar] # - Ginv @ AtA[S,Sbar]
- yw = y - A[:,S] @ b # y - A[:,S] @ b
- Aw = A[:,Sbar] + A[:,S] @ B # A[:,Sbar] + A[:,S] @ B
- tw = l1 + B.T @ (l1 + l2 * b) # l1 + B.T @ (l1 + l2 * b)
-
- # Monitoring
- monitor.obj.append(primal_value(x,A,y,l1,l2))
- monitor.gap.append(duality_gap(x,u,A,y,l1,l2))
- monitor.flops.append(0)
- monitor.sc_size.append(len(sc))
- monitor.rc_size.append(len(rc))
- if params.coherence:
- monitor.coherence.append(coherence(Aw))
- else:
- monitor.coherence.append(np.nan)
-
-
- ############################# Main loop ################################
-
- while True:
-
- it += 1
- flops = 0 # FLOPs of the current iteration
- S_size = np.sum(S) # ie. card(S)
- Sbar_size = np.sum(Sbar) # ie. card(Sbar)
-
- # Iterates update
- Awz = Aw @ zold[Sbar]
- u = yw - Awz
- Awtu = Aw.T @ u
- BBz = B.T @ (B @ zold[Sbar])
- Awtu_tw = Awtu - tw
- grad = Awtu_tw - l2 * (zold[Sbar] + BBz)
- x[Sbar] = np.maximum(zold[Sbar] + grad/L, 0)
- x[S] = np.maximum(b + B @ x[Sbar],0)
- t = (1 + np.sqrt(1 + 4 * told ** 2)) / 2
- z = x + ((told - 1) / t) * (x - xold)
-
- # ------------------------- FLOPs count ------------------------------ #
- # Iterates update
- # Awz : 2m|Sbar|
- # u : m
- # Awtu : 2m|Sbar|
- # BBz : 4|Sbar||S|
- # Awtu_tw : |Sbar|
- # grad : 2|Sbar| + |S|
- # x[Sbar] : 2|Sbar|
- # x[S] : 2|S|^2 + |S|
- # t : 6
- # z : 3n + 2
- flops += 4*m*Sbar_size + 4*Sbar_size*S_size + 2*S_size**2 + 3*n + m + \
- 5*Sbar_size + 2*S_size + 8
- # -------------------------------------------------------------------- #
-
- # Gap computation
- zpos = np.maximum(zold[Sbar],0)
- Atu = np.zeros(n)
- Atu[S] = A[:,S].T @ u
- Atu[Sbar] = Awtu - B.T @ Atu[S]
- y_u = y - u
- v = np.maximum(Atu-l1, 0)
- vpos = v[v>0]
- pv = 0.5*u@u + tw@zpos + (l2/2)*(zpos@BBz + zpos@zpos + b@b) + l1*np.sum(b)
- dv = 0.5 * (1 - (y_u@y_u) - (1/l2)*(vpos@vpos))
- gap = np.abs(pv-dv)
-
- # Update best primal and dual points recorded so far
- if gap < best_gap:
- best_gap = gap
- best_x = np.copy(np.maximum(zold,0))
- best_Atu = np.copy(Atu)
-
- # ------------------------- FLOPs count ------------------------------ #
- # GAP computation
- # Atu[S] : 2m|S|
- # Atu[Sbar] : n + 2|S|m
- # y_u : m
- # v : n
- # pv : 2m + 6|Sbar| + 3|S| + 5
- # dv : 2m + 2|vpos| + 4
- # gap : 1
- flops += 4*m*S_size + 2*n + 5*m + 5*Sbar_size + 3*S_size + 10
- # -------------------------------------------------------------------- #
-
-
- if params.screen:
-
- # Screening test
- radius = np.sqrt(2 * best_gap)
- screen = np.where(best_Atu <= l1 - radius)[0]
- sc.update(idx[screen])
-
- # Dimension reduction
- if len(screen) > 0:
- screen_Sbar = np.where(best_Atu[Sbar] <= l1 - radius)[0]
- A = np.delete(A, screen, axis=1)
- gram = np.delete(gram, screen, axis=0)
- gram = np.delete(gram, screen, axis=1)
- AtA = np.delete(AtA, screen, axis=0)
- AtA = np.delete(AtA, screen, axis=1)
- Atu = np.delete(Atu, screen)
- Aty = np.delete(Aty, screen)
- Aty_l1 = np.delete(Aty_l1, screen)
- x = np.delete(x, screen)
- z = np.delete(z, screen)
- best_x = np.delete(best_x, screen)
- best_Atu = np.delete(best_Atu, screen)
- idx = np.delete(idx, screen)
- S = np.delete(S, screen)
- Sbar = np.delete(Sbar, screen)
- Aw = np.delete(Aw, screen_Sbar, axis=1)
- B = np.delete(B, screen_Sbar, axis=1)
- tw = np.delete(tw, screen_Sbar)
- _, n = A.shape
-
- # ----------------------- FLOPs count ---------------------------- #
- # Screening
- # radius : 2
- # screen : 1
- flops += 3
- # ---------------------------------------------------------------- #
-
- if params.relax:
-
- # Relaxing test
- radius = np.sqrt(2 * best_gap) / sqrt_l2
- relax_Atu = np.where(best_Atu > l1 + radius * sqrt_l2)[0]
- relax_x = np.where(best_x > radius)[0]
- relax = np.unique(np.hstack((relax_Atu, relax_x)))
- rc.update(idx[relax])
-
- # ----------------------- FLOPs count ---------------------------- #
- # Relaxing
- # radius : 2
- # relax_Atu : 2
- flops += 4
- # ---------------------------------------------------------------- #
-
- # Update problem data
- if len(relax) > 0:
- for i in relax[Sbar[relax]]:
- Ginv, B, Aw, b, yw, tw = update_data(Ginv, B, Aw, b, yw, tw,
- i, A, y, S, Sbar, l1,
- l2, AtA, gram, Aty_l1
- )
-
- # -------------------- FLOPs count ----------------------- #
- # Update problem data
- # See algebra.update_data for FLOPs count detail
- flops += 1.5*S_size**2 + S_size*Sbar_size + m*S_size + \
- 9*S_size + Sbar_size + 14
- # -------------------------------------------------------- #
-
- S[i] = True
- Sbar[i] = False
- S_size += 1
- Sbar_size += 1
-
-
- if len(sc)+len(rc) == ninit:
-
- # Closed form for the optimizer
- x[S] = b
- x[Sbar] = 0
-
- # /!\ Following operations are only used for monitoring purpose
- # /!\ but are not needed in practice because convergence is
- # /!\ already certified. Thus, no FLOPs are counted.
- monitor.obj.append(primal_value(x,A,y,l1,l2))
- monitor.gap.append(duality_gap(x,y-A@x,A,y,l1,l2))
- monitor.flops.append(monitor.flops[-1] + flops)
- monitor.sc_size.append(len(sc))
- monitor.rc_size.append(len(rc))
- monitor.stopcrit = StopCrit.CONVERGED
- if params.coherence:
- monitor.coherence.append(coherence(Aw))
- else:
- monitor.coherence.append(np.nan)
- if params.verbosity and it % params.showevery == 0:
- print_trace(it, monitor)
- break
-
-
- # Displays and monitoring
- monitor.obj.append(pv)
- monitor.gap.append(gap)
- monitor.flops.append(monitor.flops[-1] + flops)
- monitor.sc_size.append(len(sc))
- monitor.rc_size.append(len(rc))
- if params.coherence:
- monitor.coherence.append(coherence(Aw))
- else:
- monitor.coherence.append(np.nan)
-
- if params.verbosity and it % params.showevery == 0:
- print_trace(it, monitor)
-
- if gap <= params.tol:
- monitor.stopcrit = StopCrit.CONVERGED
- break
- elif it+1 >= params.maxit:
- monitor.stopcrit = StopCrit.INTERRUPTED
- break
-
- xold = np.copy(x)
- zold = np.copy(z)
- told = np.copy(t)
-
-
- ############################## Postprocessing ##############################
-
- monitor.sc = sc
- monitor.rc = rc
- monitor.solvetime = time.time() - start_time
-
- xopt = np.zeros(ninit)
- xopt[idx] = x
-
- if params.verbosity:
- print_footer(monitor)
-
- return xopt, monitor
diff --git a/srenet/utils.py b/srenet/utils.py
deleted file mode 100644
index cc92b97..0000000
--- a/srenet/utils.py
+++ /dev/null
@@ -1,130 +0,0 @@
-import warnings
-import numpy as np
-from dataclasses import dataclass, field
-from enum import Enum
-
-_WIDTH = 75
-
-class StopCrit(Enum):
- NOT_OPTIMIZED = -1
- INTERRUPTED = 0
- CONVERGED = 1
-
- @property
- def has_converged(self):
- return (self.value > 0)
-
-@dataclass
-class TypeValidator:
- def validate_type(self):
- for f_name, f_def in self.__dataclass_fields__.items():
- f_type = type(getattr(self, f_name))
- if not (f_type == list or f_type == np.ndarray):
- if not issubclass(f_type,f_def.type):
- raise ValueError(f"Expected {f_def.type} for {f_name}")
-
-@dataclass
-class ProjgradParams(TypeValidator):
- """
- Parameters for the method projgrad in srenet.projgrad.
- """
- maxit : int = 10_000 # Maximum number of iterations
- tol : float = 1e-15 # Tolerance required on the gap
- verbosity : bool = True # Toogle verbosity of the method
- showevery : int = 10 # Show iter info every `showevery`
- screen : bool = False # Toogle screening tests
- relax : bool = False # Toogle relaxing tests
- coherence : bool = False # Compute coherence every iterations
- Sinit : list[int] = field(default_factory=list) # Init. knowledge of S
- def validate(self):
- self.validate_type()
-
-@dataclass
-class ProjgradMonitor():
- """
- Monitoring values of the method projgrad in srenet.projgrad.
- """
- obj : list[float] = field(default_factory=list) # Objective trace
- gap : list[float] = field(default_factory=list) # Gap trace
- flops : list[int] = field(default_factory=list) # FLOPs trace
- stopcrit : StopCrit = StopCrit.NOT_OPTIMIZED # Stopping criterion
- solvetime : float = -1 # Solvetime
- sc : set[int] = field(default_factory=set) # Screened coefs
- rc : set[int] = field(default_factory=set) # Relaxed coefs
- sc_size : list[int] = field(default_factory=list) # Nb screened coefs over iters.
- rc_size : list[int] = field(default_factory=list) # Nb relaxed coefs over iters.
- coherence : list[float] = field(default_factory=list) # Coherence trace
-
-def validate_inputs(x0, A, y, l1, l2):
- if not issubclass(type(x0),np.ndarray):
- raise ValueError("Expected np.ndarray for x0")
- if not issubclass(type(A),np.ndarray):
- raise ValueError("Expected np.ndarray for A")
- if not issubclass(type(y),np.ndarray):
- raise ValueError("Expected np.ndarray for y")
- if y.shape[0] != A.shape[0]:
- raise ValueError("Dimension missmatch between A and y")
- if x0.shape[0] != A.shape[1]:
- raise ValueError("Dimension missmatch between A and x0")
- if not issubclass(type(l1),float):
- raise ValueError("Expected float for l1")
- if not issubclass(type(l2),float):
- raise ValueError("Expected float for l2")
- if l1 <= 0:
- raise ValueError("Expected a positive value for l1")
- if l2 <= 0:
- raise ValueError("Expected a positive value for l2")
- if not np.allclose(np.linalg.norm(y,2),1):
- raise ValueError("Expected y to be normalized")
- for i in range(A.shape[1]):
- if not np.allclose(np.linalg.norm(A[:,i],2),1):
- raise ValueError("Expected columns of A to be normalized")
-
-def print_header(A, l1, l2, params):
- print('\n'+'*'*_WIDTH)
- print(("Projected-Gradient for the Elastic-Net problem\n").center(_WIDTH))
- print("Parameters")
- print(" dim : {:d}x{:d}".format(A.shape[0],A.shape[1]))
- print(" l1 : {:.2e}".format(l1))
- print(" l2 : {:.2e}".format(l2))
- print(" maxit : {:d}".format(params.maxit))
- print(" tol : {:.2e}".format(params.maxit))
- print(" screen : {}".format(params.screen))
- print(" relax : {}".format(params.relax))
-
- print()
-
- HEADER = "{}|{}{}{}{}{}".format(
- 'Iter'.ljust(7),
- 'Obj'.rjust(10),
- 'Gap'.rjust(10),
- 'Flops'.rjust(10),
- 'Screen'.rjust(7),
- 'Relax'.rjust(7),
- )
-
- print(HEADER)
- print("-"*len(HEADER))
-
-def print_trace(it, monitor):
- TRACE = "{}|{}{}{}{}{}".format(
- '{:d}'.format(it).ljust(7),
- '{:.2e}'.format(monitor.obj[it]).rjust(10),
- '{:.2e}'.format(monitor.gap[it]).rjust(10),
- '{:.2e}'.format(monitor.flops[it]).rjust(10),
- '{:d}'.format(monitor.sc_size[it]).rjust(7),
- '{:d}'.format(monitor.rc_size[it]).rjust(7),
- )
- print(TRACE)
-
-def print_footer(monitor):
- print()
- if not monitor.stopcrit.has_converged:
- warnings.warn("Algorithm did not converge")
- print("Stop crit : {}".format(monitor.stopcrit.name))
- print("Obj value : {:.5e}".format(monitor.obj[-1]))
- print("Final gap : {:.5e}".format(monitor.gap[-1]))
- print("Flops : {:.5e}".format(monitor.flops[-1]))
- print("Solve time : {:.2f} seconds".format(monitor.solvetime))
- print()
- print('*'*_WIDTH + '\n')
\ No newline at end of file
diff --git a/tests/convergence.py b/tests/convergence.py
deleted file mode 100644
index 4051d6d..0000000
--- a/tests/convergence.py
+++ /dev/null
@@ -1,64 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-from srenet.data import sample_data
-from srenet.projgrad import projgrad
-from srenet.utils import ProjgradParams
-
-m, n = 100, 300 # Size of A
-gen_dic = 'gaussian' # Sampling method for A
-gen_obs = 'gaussian' # Sampling method for y
-l1ratio = 0.1 # l1-regularizer / max(A.T @ y)
-l2ratio = 0.1 # l2-regularizer / max(A.T @ y)
-
-x0 = np.zeros(n)
-A, y = sample_data(m, n, gen_dic, gen_obs)
-lmax = np.max(A.T @ y)
-l1 = l1ratio * lmax
-l2 = l2ratio * lmax
-
-x, mnt = projgrad(x0, A, y, l1, l2, ProjgradParams(verbosity=False,screen=False,relax=False))
-x_s, mnt_s = projgrad(x0, A, y, l1, l2, ProjgradParams(verbosity=False,screen=True,relax=False))
-x_r, mnt_r = projgrad(x0, A, y, l1, l2, ProjgradParams(verbosity=False,screen=False,relax=True))
-x_sr, mnt_sr = projgrad(x0, A, y, l1, l2, ProjgradParams(verbosity=False,screen=True,relax=True))
-
-res = {
- 'aPG' : {'x': x, 'mnt': mnt},
- 'aPGs': {'x': x_s, 'mnt': mnt_s},
- 'aPGr': {'x': x_r, 'mnt': mnt_r},
- 'S\&R' : {'x': x_sr, 'mnt': mnt_sr},
-}
-
-print("----------------")
-print("Convergence test")
-for k, v in res.items():
- print(k)
- print(" Stop crit : {}".format(v['mnt'].stopcrit))
- print(" Last Gap : {:.2e}".format(v['mnt'].gap[-1]))
- print(" Last Obj : {:.2e}".format(v['mnt'].obj[-1]))
- print(" Iters : {}".format(len(v['mnt'].gap)))
- print(" FLOPs : {:.2e}".format(v['mnt'].flops[-1]))
-
-plt.rcParams["font.size"] = 12
-plt.rcParams["axes.labelsize"] = 12
-plt.rcParams["axes.titlesize"] = 18
-plt.rcParams["lines.linewidth"] = 2
-plt.rcParams["lines.markersize"] = 5
-plt.rcParams["xtick.labelsize"] = 12
-plt.rcParams["ytick.labelsize"] = 12
-plt.rcParams["text.usetex"] = True
-plt.rcParams["font.family"] = "serif"
-plt.rcParams["font.serif"] = ["Helvetica"]
-
-fig, ax = plt.subplots(1,1,figsize=(4,3))
-for k,v in res.items():
- ax.plot(v['mnt'].flops, v['mnt'].obj-v['mnt'].obj[-1], label=k, marker='.')
-ax.set_xlabel(r"FLOPs")
-ax.set_ylabel(r"$f(x)-f(x^{\star})$")
-ax.set_xscale('log')
-ax.set_yscale('log')
-ax.grid(b=True, which='major', axis='y')
-ax.legend()
-plt.tight_layout()
-plt.show()
-
-print("----------------")
\ No newline at end of file
diff --git a/tests/identification.py b/tests/identification.py
deleted file mode 100644
index 0aefba1..0000000
--- a/tests/identification.py
+++ /dev/null
@@ -1,48 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-from srenet.data import sample_data
-from srenet.projgrad import projgrad
-from srenet.utils import ProjgradParams
-
-m, n = 100, 300 # Size of A
-gen_dic = 'gaussian' # Sampling method for A
-gen_obs = 'gaussian' # Sampling method for y
-l1ratio = 0.1 # l1-regularizer / max(A.T @ y)
-l2ratio = 0.1 # l2-regularizer / max(A.T @ y)
-
-x0 = np.zeros(n)
-A, y = sample_data(m, n, gen_dic, gen_obs)
-lmax = np.max(A.T @ y)
-l1 = l1ratio * lmax
-l2 = l2ratio * lmax
-
-params = ProjgradParams(verbosity=False, screen=True, relax=True)
-x, mnt = projgrad(x0, A, y, l1, l2, params)
-I_pos = np.sum(np.abs(x) > 1e-8)
-I_null = np.sum(np.abs(x) <= 1e-8)
-
-print("-------------------")
-print("Identification test\n")
-print("aPGsr")
-print(" Stop crit : {}".format(mnt.stopcrit))
-print(" Last Gap : {:.2e}".format(mnt.gap[-1]))
-print(" Last Obj : {:.2e}".format(mnt.obj[-1]))
-print(" Iters : {}".format(len(mnt.gap)))
-print(" FLOPS : {:.2e}".format(mnt.flops[-1]))
-print("Sparisty of the solution : {}/{} ({:.2%})".format(I_pos, n, I_pos/n))
-
-fig, ax = plt.subplots(1,1,figsize=(4,3))
-
-ax.plot(mnt.gap, mnt.sc_size/I_null, label="% Screened", marker='.')
-ax.plot(mnt.gap, mnt.rc_size/I_pos, label="% Relaxed", marker='.')
-ax.yaxis.set_major_formatter(lambda x,pos: r"{0:.0f}%".format(x * 100))
-ax.invert_xaxis()
-ax.set_xlabel("GAP")
-ax.set_ylabel("Indices identified")
-ax.set_xscale('log')
-ax.grid(b=True, which='major', axis='y')
-ax.legend()
-plt.tight_layout()
-plt.show()
-
-print("-------------------")
\ No newline at end of file
--
GitLab