From cee6af6aef6b08eaf81578a2ab34616b9af91467 Mon Sep 17 00:00:00 2001 From: nicolas barascud <10333715+nbara@users.noreply.github.com> Date: Wed, 7 Apr 2021 16:14:57 +0200 Subject: [PATCH 1/5] [FIX] Detrending error with high-order polynomials (#36) --- meegkit/detrend.py | 33 ++++++++++++++++----------------- tests/test_detrend.py | 32 +++++++++++++++++++++++++------- 2 files changed, 41 insertions(+), 24 deletions(-) diff --git a/meegkit/detrend.py b/meegkit/detrend.py index b829a1a5..de2f2778 100644 --- a/meegkit/detrend.py +++ b/meegkit/detrend.py @@ -80,8 +80,8 @@ def detrend(x, order, w=None, basis='polynomials', threshold=3, n_iter=4, elif basis == 'sinusoids': r = np.zeros((n_times, order * 2)) for i, o in enumerate(range(1, order + 1)): - r[:, 2 * i] = np.sin[2 * np.pi * o * lin / 2] - r[:, 2 * i + 1] = np.cos[2 * np.pi * o * lin / 2] + r[:, 2 * i] = np.sin(2 * np.pi * o * lin / 2) + r[:, 2 * i + 1] = np.cos(2 * np.pi * o * lin / 2) else: raise ValueError('!') @@ -129,7 +129,8 @@ def regress(x, r, w=None, threshold=1e-7, return_mean=False): Weight to apply to `x`. `w` is either a matrix of same size as `x`, or a column vector to be applied to each column of `x`. threshold : float - PCA threshold (default=1e-7). + PCA threshold (default=1e-7). Dimensions of x with eigenvalue lower + than this value will be discarded. return_mean : bool If True, also return the signal mean prior to regression. @@ -145,6 +146,7 @@ def regress(x, r, w=None, threshold=1e-7, return_mean=False): w = _check_weights(w, x) n_times = x.shape[0] n_chans = x.shape[1] + n_regs = r.shape[1] r = unfold(r) x = unfold(x) if r.shape[0] != x.shape[0]: @@ -178,9 +180,9 @@ def regress(x, r, w=None, threshold=1e-7, return_mean=False): else: yy = demean(x, w) * w rr = demean(r, w) * w - V, _ = pca(rr.T.dot(rr), thresh=threshold) - rrr = rr.dot(V) - b = mrdivide(yy.T.dot(rrr), rrr.T.dot(rrr)) + V, _ = pca(rr.T @ rr, thresh=threshold) + rr = rr @ V + b = mrdivide(yy.T @ rr, rr.T @ rr) z = demean(r, w).dot(V).dot(b.T) z = z + mn @@ -189,25 +191,22 @@ def regress(x, r, w=None, threshold=1e-7, return_mean=False): if w.shape[1] != x.shape[1]: raise ValueError('!') z = np.zeros(x.shape) - b = np.zeros((n_chans, n_chans)) + b = np.zeros((n_chans, n_regs)) for i in range(n_chans): - if sum(w[:, i]) == 0: - print('weights all zero for channel {}'.format(i)) - c = np.zeros(x.shape[1], 1) + if not np.any(w[:, i]): + print(f'weights are all zero for channel {i}') else: wc = w[:, i][:, None] # channel-specific weight - yy = demean(x[:, i], wc) * wc + xx = demean(x[:, i], wc) * wc + # remove channel-specific-weighted mean from regressor r = demean(r, wc) rr = r * wc - V, _ = pca(rr.T.dot(rr), thresh=threshold) + V, _ = pca(rr.T @ rr, thresh=threshold) rr = rr.dot(V) - c = mrdivide(yy.T.dot(rr), rr.T.dot(rr)) + b[i, :V.shape[1]] = mrdivide(xx.T @ rr, rr.T @ rr) - z[:, i] = r.dot(V.dot(c.T)).flatten() - z[:, i] += mn[:, i] - b[i] = c - b = b[:, :V.shape[1]] + z[:, i] = np.squeeze(r @ (V @ b[i, :V.shape[1]].T)) + mn[:, i] return b, z diff --git a/tests/test_detrend.py b/tests/test_detrend.py index b4bd7761..477e6cb6 100644 --- a/tests/test_detrend.py +++ b/tests/test_detrend.py @@ -39,14 +39,15 @@ def test_regress(): w[:, 1] == .8 [b, z] = regress(y, x, w) assert z.shape == (1000, 2) - assert b.shape == (2, 1) + assert b.shape == (2, 3) def test_detrend(show=False): """Test detrending.""" # basic - x = np.arange(100)[:, None] - x = x + np.random.randn(*x.shape) + x = np.arange(100)[:, None] # trend + source = np.random.randn(*x.shape) + x = x + source y, _, _ = detrend(x, 1) assert y.shape == x.shape @@ -57,20 +58,37 @@ def test_detrend(show=False): assert y.shape == x.shape - # weights + # test weights trend = np.linspace(0, 100, 1000)[:, None] data = 3 * np.random.randn(*trend.shape) + data[:100, :] = 100 x = trend + data - x[:100, :] = 100 w = np.ones(x.shape) w[:100, :] = 0 + + # Without weights – detrending fails y, _, _ = detrend(x, 3, None, show=show) + + # With weights – detrending works yy, _, _ = detrend(x, 3, w, show=show) assert y.shape == x.shape assert yy.shape == x.shape + assert np.all(np.abs(yy[100:] - data[100:]) < 1.) + + # detrend higher-dimensional data + x = np.cumsum(np.random.randn(1000, 16) + 0.1, axis=0) + y, _, _ = detrend(x, 1, show=False) + + # detrend higher-dimensional data with order 3 polynomial + x = np.cumsum(np.random.randn(1000, 16) + 0.1, axis=0) + y, _, _ = detrend(x, 3, basis='polynomials', show=True) + + # detrend with sinusoids + x = np.random.randn(1000, 2) + x += 2 * np.sin(2 * np.pi * np.arange(1000) / 200)[:, None] + y, _, _ = detrend(x, 5, basis='sinusoids', show=True) - # assert_almost_equal(yy[100:], data[100:], decimal=1) def test_ringing(): """Test reduce_ringing function.""" @@ -87,4 +105,4 @@ def test_ringing(): if __name__ == '__main__': import pytest pytest.main([__file__]) - # test_detrend(True) + # test_detrend(False) From 52419ad1bc814a29b50ffeda6d12cc52787a9172 Mon Sep 17 00:00:00 2001 From: gferraro2019 <47594930+gferraro2019@users.noreply.github.com> Date: Wed, 7 Apr 2021 20:25:11 +0200 Subject: [PATCH 2/5] Add Task Related Component Analysis (TRCA) (#32) * add the Python translation of Masaki Nakanishi's TRCA and Filterbank Matlab tutorial * cleaning and editing scripts * cleaning and editing scripts * deleting lines scripts * cleaning and editing scripts * cleaning and editing scripts * cleaning and editing scripts * cleaning and editing utils * cleaning and editing utils * cleaning and editing utils * last attempt cleaning and editing utils * start refactor - actually integrate TRCA code into python package - regroup trca utils - update README - fix PEP errors, improve docstrings all around - draft example * continue refactor + lots of cleanup * add basic Class implementation + unit test * fix doc * comments * Update example_trca.ipynb * use tuples for filterbank spec Co-authored-by: Giuseppe Co-authored-by: nbara <10333715+nbara@users.noreply.github.com> --- README.md | 19 +++ doc/index.rst | 1 + doc/modules/meegkit.utils.rst | 12 +- examples/example_asr.py | 1 - examples/example_trca.ipynb | 258 ++++++++++++++++++++++++++++++++++ examples/example_trca.py | 143 +++++++++++++++++++ meegkit/__init__.py | 6 +- meegkit/trca.py | 218 ++++++++++++++++++++++++++++ meegkit/utils/stats.py | 6 +- meegkit/utils/trca.py | 137 ++++++++++++++++++ setup.cfg | 2 +- tests/data/trcadata.mat | Bin 0 -> 20987516 bytes tests/test_trca.py | 107 ++++++++++++++ 13 files changed, 901 insertions(+), 9 deletions(-) create mode 100644 examples/example_trca.ipynb create mode 100644 examples/example_trca.py create mode 100644 meegkit/trca.py create mode 100644 meegkit/utils/trca.py create mode 100644 tests/data/trcadata.mat create mode 100644 tests/test_trca.py diff --git a/README.md b/README.md index 82a095ad..27d6f260 100644 --- a/README.md +++ b/README.md @@ -91,3 +91,22 @@ If you use this, you should cite the following article: [1] Cohen, M. X., & Gulbinaite, R. (2017). Rhythmic entrainment source separation: Optimizing analyses of neural responses to rhythmic sensory stimulation. Neuroimage, 147, 43-56. ``` + +### 4. Task-Related Component Analysis (TRCA) + +This code is based on the [Matlab implementation from Masaki Nakanishi](https://github.com/mnakanishi/TRCA-SSVEP), and was adapted to python by [Giuseppe Ferraro](mailto:giuseppe.ferraro@isae-supaero.fr) + +If you use this, you should cite the following articles: + +```sql +[1] M. Nakanishi, Y. Wang, X. Chen, Y.-T. Wang, X. Gao, and T.-P. Jung, + "Enhancing detection of SSVEPs for a high-speed brain speller using + task-related component analysis", IEEE Trans. Biomed. Eng, 65(1): 104-112, + 2018. +[2] X. Chen, Y. Wang, S. Gao, T. -P. Jung and X. Gao, "Filter bank + canonical correlation analysis for implementing a high-speed SSVEP-based + brain-computer interface", J. Neural Eng., 12: 046008, 2015. +[3] X. Chen, Y. Wang, M. Nakanishi, X. Gao, T. -P. Jung, S. Gao, + "High-speed spelling with a noninvasive brain-computer interface", + Proc. Int. Natl. Acad. Sci. U. S. A, 112(44): E6058-6067, 2015. +``` diff --git a/doc/index.rst b/doc/index.rst index f733c353..e4a89444 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -26,6 +26,7 @@ Contents ~meegkit.ress ~meegkit.sns ~meegkit.star + ~meegkit.trca ~meegkit.tspca ~meegkit.utils diff --git a/doc/modules/meegkit.utils.rst b/doc/modules/meegkit.utils.rst index a1919394..9797de75 100644 --- a/doc/modules/meegkit.utils.rst +++ b/doc/modules/meegkit.utils.rst @@ -71,9 +71,19 @@ Signal ---- - Statistics ---------- .. automodule:: meegkit.utils.stats .. autosummary:: + + +| + +---- + +TRCA utilities +-------------- +.. automodule:: meegkit.utils.trca + + .. autosummary:: diff --git a/examples/example_asr.py b/examples/example_asr.py index 77c70752..d8634d6c 100644 --- a/examples/example_asr.py +++ b/examples/example_asr.py @@ -11,7 +11,6 @@ import matplotlib.pyplot as plt from meegkit.asr import ASR -from meegkit.utils.asr import yulewalk_filter from meegkit.utils.matrix import sliding_window # THIS_FOLDER = os.path.dirname(os.path.abspath(__file__)) diff --git a/examples/example_trca.ipynb b/examples/example_trca.ipynb new file mode 100644 index 00000000..78382686 --- /dev/null +++ b/examples/example_trca.ipynb @@ -0,0 +1,258 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "# Task-related component analysis (TRCA)-based SSVEP detection\n", + "\n", + "Sample code for the task-related component analysis (TRCA)-based steady\n", + "-state visual evoked potential (SSVEP) detection method [1]_. The filter\n", + "bank analysis [2, 3]_ can also be combined to the TRCA-based algorithm.\n", + "\n", + "Uses meegkit.trca.TRCA()\n", + "\n", + "References:\n", + "\n", + ".. [1] M. Nakanishi, Y. Wang, X. Chen, Y.-T. Wang, X. Gao, and T.-P. Jung,\n", + " \"Enhancing detection of SSVEPs for a high-speed brain speller using\n", + " task-related component analysis\", IEEE Trans. Biomed. Eng, 65(1): 104-112,\n", + " 2018.\n", + ".. [2] X. Chen, Y. Wang, S. Gao, T. -P. Jung and X. Gao, \"Filter bank\n", + " canonical correlation analysis for implementing a high-speed SSVEP-based\n", + " brain-computer interface\", J. Neural Eng., 12: 046008, 2015.\n", + ".. [3] X. Chen, Y. Wang, M. Nakanishi, X. Gao, T. -P. Jung, S. Gao,\n", + " \"High-speed spelling with a noninvasive brain-computer interface\",\n", + " Proc. Int. Natl. Acad. Sci. U. S. A, 112(44): E6058-6067, 2015.\n", + "\n", + "This code is based on the Matlab implementation from\n", + "https://github.com/mnakanishi/TRCA-SSVEP\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Author: Giuseppe Ferraro \n", + "import os\n", + "import time\n", + "\n", + "import numpy as np\n", + "import scipy.io\n", + "from meegkit.trca import TRCA\n", + "from meegkit.utils.trca import itr, normfit, round_half_up\n", + "\n", + "t = time.time()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Parameters\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Results of the ensemble TRCA-based method:\n\n" + ] + } + ], + "source": [ + "len_gaze_s = 0.5 # data length for target identification [s]\n", + "len_delay_s = 0.13 # visual latency being considered in the analysis [s]\n", + "n_bands = 5 # number of sub-bands in filter bank analysis\n", + "is_ensemble = True # True = ensemble TRCA method; False = TRCA method\n", + "alpha_ci = 0.05 # 100*(1-alpha_ci): confidence interval for accuracy\n", + "sfreq = 250 # sampling rate [Hz]\n", + "len_shift_s = 0.5 # duration for gaze shifting [s]\n", + "list_freqs = np.concatenate(\n", + " [[x + 8 for x in range(8)],\n", + " [x + 8.2 for x in range(8)],\n", + " [x + 8.4 for x in range(8)],\n", + " [x + 8.6 for x in range(8)],\n", + " [x + 8.8 for x in range(8)]]) # list of stimulus frequencies\n", + "n_targets = len(list_freqs) # The number of stimuli\n", + "\n", + "# Preparing useful variables (DONT'T need to modify)\n", + "len_gaze_smpl = round_half_up(len_gaze_s * sfreq) # data length [samples]\n", + "len_delay_smpl = round_half_up(len_delay_s * sfreq) # visual latency [samples]\n", + "len_sel_s = len_gaze_s + len_shift_s # selection time [s]\n", + "ci = 100 * (1 - alpha_ci) # confidence interval\n", + "\n", + "# Performing the TRCA-based SSVEP detection algorithm\n", + "print('Results of the ensemble TRCA-based method:\\n')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Load data\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "path = os.path.join('..', 'tests', 'data', 'trcadata.mat')\n", + "mat = scipy.io.loadmat(path)\n", + "eeg = mat[\"eeg\"]\n", + "\n", + "n_trials = eeg.shape[0]\n", + "n_chans = eeg.shape[1]\n", + "n_samples = eeg.shape[2]\n", + "n_blocks = eeg.shape[3]\n", + "\n", + "# Convert dummy Matlab format to (sample, channels, trials) and construct\n", + "# vector of labels\n", + "eeg = np.reshape(eeg.transpose([2, 1, 3, 0]),\n", + " (n_samples, n_chans, n_trials * n_blocks))\n", + "labels = np.array([x for x in range(n_targets)] * n_blocks)\n", + "\n", + "crop_data = np.arange(len_delay_smpl, len_delay_smpl + len_gaze_smpl)\n", + "eeg = eeg[crop_data]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## TRCA classification\n", + "Estimate classification performance with a Leave-One-Block-Out\n", + "cross-validation approach.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Block 0: accuracy = 97.5, \tITR = 301.3\n", + "Block 1: accuracy = 100.0, \tITR = 319.3\n", + "Block 2: accuracy = 95.0, \tITR = 286.3\n", + "Block 3: accuracy = 95.0, \tITR = 286.3\n", + "Block 4: accuracy = 95.0, \tITR = 286.3\n", + "Block 5: accuracy = 100.0, \tITR = 319.3\n", + "\n", + "Mean accuracy = 97.1%\t(95% CI: 97.0-97.1%)\n", + "Mean ITR = 299.8\t(95% CI: 299.4-300.2%)\n", + "\n", + "Elapsed time: 13.8 seconds\n" + ] + } + ], + "source": [ + "# We use the filterbank specification described in [2]_.\n", + "filterbank = [[(6, 90), (4, 100)], # passband freqs, stopband freqs (Wp, Ws)\n", + " [(14, 90), (10, 100)],\n", + " [(22, 90), (16, 100)],\n", + " [(30, 90), (24, 100)],\n", + " [(38, 90), (32, 100)],\n", + " [(46, 90), (40, 100)],\n", + " [(54, 90), (48, 100)]]\n", + "trca = TRCA(sfreq, filterbank, is_ensemble)\n", + "\n", + "accs = np.zeros(n_blocks)\n", + "itrs = np.zeros(n_blocks)\n", + "for i in range(n_blocks):\n", + "\n", + " # Training stage\n", + " traindata = eeg.copy()\n", + "\n", + " # Select all folds except one for training\n", + " traindata = np.concatenate(\n", + " (traindata[..., :i * n_trials],\n", + " traindata[..., (i + 1) * n_trials:]), 2)\n", + " y_train = np.concatenate(\n", + " (labels[:i * n_trials], labels[(i + 1) * n_trials:]), 0)\n", + "\n", + " # Construction of the spatial filter and the reference signals\n", + " trca.fit(traindata, y_train)\n", + "\n", + " # Test stage\n", + " testdata = eeg[..., i * n_trials:(i + 1) * n_trials]\n", + " y_test = labels[i * n_trials:(i + 1) * n_trials]\n", + " estimated = trca.predict(testdata)\n", + "\n", + " # Evaluation of the performance for this fold (accuracy and ITR)\n", + " is_correct = estimated == y_test\n", + " accs[i] = np.mean(is_correct) * 100\n", + " itrs[i] = itr(n_targets, np.mean(is_correct), len_sel_s)\n", + " print(f\"Block {i}: accuracy = {accs[i]:.1f}, \\tITR = {itrs[i]:.1f}\")\n", + "\n", + "# Mean accuracy and ITR computation\n", + "mu, _, muci, _ = normfit(accs, alpha_ci)\n", + "print()\n", + "print(f\"Mean accuracy = {mu:.1f}%\\t({ci:.0f}% CI: {muci[0]:.1f}-{muci[1]:.1f}%)\") # noqa\n", + "\n", + "mu, _, muci, _ = normfit(itrs, alpha_ci)\n", + "print(f\"Mean ITR = {mu:.1f}\\t({ci:.0f}% CI: {muci[0]:.1f}-{muci[1]:.1f}%)\")\n", + "if is_ensemble:\n", + " ensemble = 'ensemble TRCA-based method'\n", + "else:\n", + " ensemble = 'TRCA-based method'\n", + "\n", + "print(f\"\\nElapsed time: {time.time()-t:.1f} seconds\")" + ] + } + ], + "metadata": { + "kernelspec": { + "name": "python388jvsc74a57bd0d64e410d98a0dc7c6b3fb09ececfc32281268599ac952adfc85e199a2f396698", + "display_name": "Python 3.8.8 64-bit ('miniconda3': virtualenv)" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.8-final" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file diff --git a/examples/example_trca.py b/examples/example_trca.py new file mode 100644 index 00000000..f3f93df7 --- /dev/null +++ b/examples/example_trca.py @@ -0,0 +1,143 @@ +""" +Task-related component analysis (TRCA)-based SSVEP detection +============================================================ + +Sample code for the task-related component analysis (TRCA)-based steady +-state visual evoked potential (SSVEP) detection method [1]_. The filter +bank analysis [2, 3]_ can also be combined to the TRCA-based algorithm. + +Uses meegkit.trca.TRCA() + +References: + +.. [1] M. Nakanishi, Y. Wang, X. Chen, Y.-T. Wang, X. Gao, and T.-P. Jung, + "Enhancing detection of SSVEPs for a high-speed brain speller using + task-related component analysis", IEEE Trans. Biomed. Eng, 65(1): 104-112, + 2018. +.. [2] X. Chen, Y. Wang, S. Gao, T. -P. Jung and X. Gao, "Filter bank + canonical correlation analysis for implementing a high-speed SSVEP-based + brain-computer interface", J. Neural Eng., 12: 046008, 2015. +.. [3] X. Chen, Y. Wang, M. Nakanishi, X. Gao, T. -P. Jung, S. Gao, + "High-speed spelling with a noninvasive brain-computer interface", + Proc. Int. Natl. Acad. Sci. U. S. A, 112(44): E6058-6067, 2015. + +This code is based on the Matlab implementation from +https://github.com/mnakanishi/TRCA-SSVEP + +""" +# Author: Giuseppe Ferraro +import os +import time + +import numpy as np +import scipy.io +from meegkit.trca import TRCA +from meegkit.utils.trca import itr, normfit, round_half_up + +t = time.time() + +############################################################################### +# Parameters +# ----------------------------------------------------------------------------- +len_gaze_s = 0.5 # data length for target identification [s] +len_delay_s = 0.13 # visual latency being considered in the analysis [s] +n_bands = 5 # number of sub-bands in filter bank analysis +is_ensemble = True # True = ensemble TRCA method; False = TRCA method +alpha_ci = 0.05 # 100*(1-alpha_ci): confidence interval for accuracy +sfreq = 250 # sampling rate [Hz] +len_shift_s = 0.5 # duration for gaze shifting [s] +list_freqs = np.concatenate( + [[x + 8 for x in range(8)], + [x + 8.2 for x in range(8)], + [x + 8.4 for x in range(8)], + [x + 8.6 for x in range(8)], + [x + 8.8 for x in range(8)]]) # list of stimulus frequencies +n_targets = len(list_freqs) # The number of stimuli + +# Preparing useful variables (DONT'T need to modify) +len_gaze_smpl = round_half_up(len_gaze_s * sfreq) # data length [samples] +len_delay_smpl = round_half_up(len_delay_s * sfreq) # visual latency [samples] +len_sel_s = len_gaze_s + len_shift_s # selection time [s] +ci = 100 * (1 - alpha_ci) # confidence interval + +# Performing the TRCA-based SSVEP detection algorithm +print('Results of the ensemble TRCA-based method:\n') + +############################################################################### +# Load data +# ----------------------------------------------------------------------------- +path = os.path.join('..', 'tests', 'data', 'trcadata.mat') +mat = scipy.io.loadmat(path) +eeg = mat["eeg"] + +n_trials = eeg.shape[0] +n_chans = eeg.shape[1] +n_samples = eeg.shape[2] +n_blocks = eeg.shape[3] + +# Convert dummy Matlab format to (sample, channels, trials) and construct +# vector of labels +eeg = np.reshape(eeg.transpose([2, 1, 3, 0]), + (n_samples, n_chans, n_trials * n_blocks)) +labels = np.array([x for x in range(n_targets)] * n_blocks) + +crop_data = np.arange(len_delay_smpl, len_delay_smpl + len_gaze_smpl) +eeg = eeg[crop_data] + +############################################################################### +# TRCA classification +# ----------------------------------------------------------------------------- +# Estimate classification performance with a Leave-One-Block-Out +# cross-validation approach. + +# We use the filterbank specification described in [2]_. +filterbank = [[(6, 90), (4, 100)], # passband, stopband freqs [(Wp), (Ws)] + [(14, 90), (10, 100)], + [(22, 90), (16, 100)], + [(30, 90), (24, 100)], + [(38, 90), (32, 100)], + [(46, 90), (40, 100)], + [(54, 90), (48, 100)]] +trca = TRCA(sfreq, filterbank, is_ensemble) + +accs = np.zeros(n_blocks) +itrs = np.zeros(n_blocks) +for i in range(n_blocks): + + # Training stage + traindata = eeg.copy() + + # Select all folds except one for training + traindata = np.concatenate( + (traindata[..., :i * n_trials], + traindata[..., (i + 1) * n_trials:]), 2) + y_train = np.concatenate( + (labels[:i * n_trials], labels[(i + 1) * n_trials:]), 0) + + # Construction of the spatial filter and the reference signals + trca.fit(traindata, y_train) + + # Test stage + testdata = eeg[..., i * n_trials:(i + 1) * n_trials] + y_test = labels[i * n_trials:(i + 1) * n_trials] + estimated = trca.predict(testdata) + + # Evaluation of the performance for this fold (accuracy and ITR) + is_correct = estimated == y_test + accs[i] = np.mean(is_correct) * 100 + itrs[i] = itr(n_targets, np.mean(is_correct), len_sel_s) + print(f"Block {i}: accuracy = {accs[i]:.1f}, \tITR = {itrs[i]:.1f}") + +# Mean accuracy and ITR computation +mu, _, muci, _ = normfit(accs, alpha_ci) +print() +print(f"Mean accuracy = {mu:.1f}%\t({ci:.0f}% CI: {muci[0]:.1f}-{muci[1]:.1f}%)") # noqa + +mu, _, muci, _ = normfit(itrs, alpha_ci) +print(f"Mean ITR = {mu:.1f}\t({ci:.0f}% CI: {muci[0]:.1f}-{muci[1]:.1f}%)") +if is_ensemble: + ensemble = 'ensemble TRCA-based method' +else: + ensemble = 'TRCA-based method' + +print(f"\nElapsed time: {time.time()-t:.1f} seconds") diff --git a/meegkit/__init__.py b/meegkit/__init__.py index 795795bd..162a3c54 100644 --- a/meegkit/__init__.py +++ b/meegkit/__init__.py @@ -1,7 +1,7 @@ """M/EEG denoising utilities in python.""" __version__ = '0.1.1' -from . import asr, cca, detrend, dss, sns, star, ress, tspca, utils +from . import asr, cca, detrend, dss, sns, star, ress, trca, tspca, utils -__all__ = ['asr', 'cca', 'detrend', 'dss', 'ress', 'sns', 'star', 'tspca', - 'utils'] +__all__ = ['asr', 'cca', 'detrend', 'dss', 'ress', 'sns', 'star', 'trca', + 'tspca', 'utils'] diff --git a/meegkit/trca.py b/meegkit/trca.py new file mode 100644 index 00000000..2f3c2e75 --- /dev/null +++ b/meegkit/trca.py @@ -0,0 +1,218 @@ +"""Task-Related Component Analysis (TRCA).""" +# Author: Giuseppe Ferraro +import numpy as np +import scipy.linalg as linalg + +from .utils.trca import bandpass +from .utils import theshapeof + + +def trca(X): + """Task-related component analysis (TRCA). + + This function implements the method described in [1]_. + + Parameters + ---------- + X : array, shape=(n_samples, n_chans[, n_trials]) + Training data. + + Returns + ------- + W : array, shape=(n_chans,) + Weight coefficients for electrodes which can be used as a spatial + filter. + + References + ---------- + .. [1] M. Nakanishi, Y. Wang, X. Chen, Y. -T. Wang, X. Gao, and T.-P. Jung, + "Enhancing detection of SSVEPs for a high-speed brain speller using + task-related component analysis", IEEE Trans. Biomed. Eng, + 65(1):104-112, 2018. + + """ + n_samples, n_chans, n_trials = theshapeof(X) + + S = np.zeros((n_chans, n_chans)) + for trial_i in range(n_trials - 1): + x1 = np.squeeze(X[..., trial_i]) + + # Mean centering for the selected trial + x1 -= np.mean(x1, 0) + + # Select a second trial that is different + for trial_j in range(trial_i + 1, n_trials): + x2 = np.squeeze(X[..., trial_j]) + + # Mean centering for the selected trial + x2 -= np.mean(x2, 0) + + # Compute empirical covariance between the two selected trials and + # sum it + S = S + x1.T @ x2 + x2.T @ x1 + + # Reshape to have all the data as a sequence + UX = np.zeros((n_chans, n_samples * n_trials)) + for trial in range(n_trials): + UX[:, trial * n_samples:(trial + 1) * n_samples] = X[..., trial].T + + # Mean centering + UX -= np.mean(UX, 1)[:, None] + + # Compute empirical variance of all data (to be bounded) + Q = np.dot(UX, UX.T) + + # Compute eigenvalues and vectors + lambdas, W = linalg.eig(S, Q, left=True, right=False) + + # Select the eigenvector corresponding to the biggest eigenvalue + W_best = W[:, np.argmax(lambdas)] + + return W_best + + +class TRCA: + """Task-Related Component Analysis (TRCA). + + Parameters + ---------- + sfreq : float + Sampling rate. + filterbank : list[[2-tuple, 2-tuple]] + Filterbank frequencies. Each list element is itself a list of passband + `Wp` and stopband `Ws` edges frequencies `[Wp, Ws]`. For example, this + creates 3 bands, starting at 6, 14, and 22 hz respectively:: + + [[(6, 90), (4, 100)], + [(14, 90), (10, 100)], + [(22, 90), (16, 100)]] + + See :func:`scipy.signal.cheb1ord()` for more information on how to + specify the `Wp` and `Ws`. + ensemble: bool + If True, perform the ensemble TRCA analysis (default=False). + + Attributes + ---------- + traindata : array, shape=(n_bands, n_chans, n_trials) + Reference (training) data decomposed into sub-band components by the + filter bank analysis. + y_train : array, shape=(n_trials) + Labels associated with the train data. + coef_ : array, shape=(n_chans, n_chans) + Weight coefficients for electrodes which can be used as a spatial + filter. + classes : list + Classes. + n_bands : int + Number of sub-bands + + """ + + def __init__(self, sfreq, filterbank, ensemble=False): + self.sfreq = sfreq + self.ensemble = ensemble + self.filterbank = filterbank + self.n_bands = len(self.filterbank) + self.coef_ = None + + def fit(self, X, y): + """Training stage of the TRCA-based SSVEP detection. + + Parameters + ---------- + X : array, shape=(n_samples, n_chans[, n_trials]) + Training EEG data. + y : array, shape=(trials,) + True label corresponding to each trial of the data array. + + """ + n_samples, n_chans, _ = theshapeof(X) + classes = np.unique(y) + + trains = np.zeros((len(classes), self.n_bands, n_samples, n_chans)) + + W = np.zeros((self.n_bands, len(classes), n_chans)) + + for class_i in classes: + # Select data with a specific label + eeg_tmp = X[..., y == class_i] + for fb_i in range(self.n_bands): + # Filter the signal with fb_i + eeg_tmp = bandpass(eeg_tmp, self.sfreq, + Wp=self.filterbank[fb_i][0], + Ws=self.filterbank[fb_i][1]) + if (eeg_tmp.ndim == 3): + # Compute mean of the signal across trials + trains[class_i, fb_i] = np.mean(eeg_tmp, -1) + else: + trains[class_i, fb_i] = eeg_tmp + # Find the spatial filter for the corresponding filtered signal + # and label + w_best = trca(eeg_tmp) + W[fb_i, class_i, :] = w_best # Store the spatial filter + + self.trains = trains + self.coef_ = W + self.classes = classes + + return self + + def predict(self, X): + """Test phase of the TRCA-based SSVEP detection. + + Parameters + ---------- + X: array, shape=(n_samples, n_chans[, n_trials]) + Test data. + model: dict + Fitted model to be used in testing phase. + + Returns + ------- + pred: np.array, shape (trials) + The target estimated by the method. + + """ + if self.coef_ is None: + raise RuntimeError('TRCA is not fitted') + + # Alpha coefficients for the fusion of filterbank analysis + fb_coefs = [(x + 1)**(-1.25) + 0.25 for x in range(self.n_bands)] + _, _, n_trials = theshapeof(X) + + r = np.zeros((self.n_bands, len(self.classes))) + pred = np.zeros((n_trials), 'int') # To store predictions + + for trial in range(n_trials): + test_tmp = X[..., trial] # Pick a trial to be analysed + for fb_i in range(self.n_bands): + + # Filterbank on testdata + testdata = bandpass(test_tmp, self.sfreq, + Wp=self.filterbank[fb_i][0], + Ws=self.filterbank[fb_i][1]) + + for class_i in self.classes: + # Retrieve reference signal for class i + # (shape: n_chans, n_samples) + traindata = np.squeeze(self.trains[class_i, fb_i]) + if self.ensemble: + # Shape of (# of channel, # of class) + w = np.squeeze(self.coef_[fb_i]).T + else: + # Shape of (# of channel) + w = np.squeeze(self.coef_[fb_i, class_i]) + + # Compute 2D correlation of spatially filtered test data + # with ref + r_tmp = np.corrcoef((testdata @ w).flatten(), + (traindata @ w).flatten()) + r[fb_i, class_i] = r_tmp[0, 1] + + rho = np.dot(fb_coefs, r) # Fusion for the filterbank analysis + + tau = np.argmax(rho) # Retrieving the index of the max + pred[trial] = int(tau) + + return pred diff --git a/meegkit/utils/stats.py b/meegkit/utils/stats.py index 1fe1a6e7..cc694857 100644 --- a/meegkit/utils/stats.py +++ b/meegkit/utils/stats.py @@ -33,7 +33,7 @@ def robust_mean(X, axis=0, percentile=[5, 95]): return m -def rolling_corr(X, y, window=None, fs=1, step=1, axis=0): +def rolling_corr(X, y, window=None, sfreq=1, step=1, axis=0): """Calculate rolling correlation between some data and a reference signal. Parameters @@ -44,7 +44,7 @@ def rolling_corr(X, y, window=None, fs=1, step=1, axis=0): Reference signal. window : int Number of timepoints for to include for each correlation calculation. - fs: int + sfreq: int Sampling frequency (default=1). step : int If > 1, only compute correlations every `step` samples. @@ -83,7 +83,7 @@ def rolling_corr(X, y, window=None, fs=1, step=1, axis=0): corr = corr.squeeze(-1) # Times relative to end of window - t_corr = (timebins + window) / float(fs) + t_corr = (timebins + window) / float(sfreq) assert len(t_corr) == corr.shape[0] diff --git a/meegkit/utils/trca.py b/meegkit/utils/trca.py new file mode 100644 index 00000000..c596dd2a --- /dev/null +++ b/meegkit/utils/trca.py @@ -0,0 +1,137 @@ +"""TRCA utils.""" +import numpy as np +import scipy + +from scipy.signal import filtfilt, cheb1ord, cheby1 +from scipy.stats import chi2, t + + +def round_half_up(num, decimals=0): + """Round half up round the last decimal of the number. + + The rules are: + from 0 to 4 rounds down + from 5 to 9 rounds up + + Parameters + ---------- + num : float + Number to round + decimals : number of decimals + + Returns + ------- + num rounded + """ + multiplier = 10 ** decimals + return int(np.floor(num * multiplier + 0.5) / multiplier) + + +def normfit(data, ci=0.95): + """Compute the mean, std and confidence interval for them. + + Parameters + ---------- + data: array, shape=() + Input data. + ci : float + Confidence interval (default=0.95). + + Returns + ------- + m : mean + sigma : std deviation + [m - h, m + h] : confidence interval of the mean + [sigmaCI_lower, sigmaCI_upper] : confidence interval of the std + """ + arr = 1.0 * np.array(data) + num = len(arr) + avg, std_err = np.mean(arr), scipy.stats.sem(arr) + h_int = std_err * t.ppf((1 + ci) / 2., num - 1) + var = np.var(data, ddof=1) + var_ci_upper = var * (num - 1) / (chi2.ppf((1 - ci) / 2, num - 1)) + var_ci_lower = var * (num - 1) / (chi2.ppf(1 - (1 - ci) / 2, num - 1)) + sigma = np.sqrt(var) + sigma_ci_lower = np.sqrt(var_ci_lower) + sigma_ci_upper = np.sqrt(var_ci_upper) + + return avg, sigma, [avg - h_int, avg + + h_int], [sigma_ci_lower, sigma_ci_upper] + + +def itr(n, p, t): + """Compute information transfer rate (ITR). + + Inputs + ------ + n : int + Number of targets. + p : float + Target identification accuracy (0 <= p <= 1). + t : float + Average time for a selection (s). + + Returns + ------- + itr : float + Information transfer rate [bits/min] + + References + ---------- + .. [1] M. Cheng, X. Gao, S. Gao, and D. Xu, + "Design and Implementation of a Brain-Computer Interface With High + Transfer Rates", IEEE Trans. Biomed. Eng. 49, 1181-1186, 2002. + + """ + itr = 0 + + if (p < 0 or 1 < p): + raise ValueError('Accuracy need to be between 0 and 1.') + elif (p < 1 / n): + raise ValueError('ITR might be incorrect because accuracy < chance') + itr = 0 + elif (p == 1): + itr = np.log2(n) * 60 / t + else: + itr = (np.log2(n) + p * np.log2(p) + (1 - p) * + np.log2((1 - p) / (n - 1))) * 60 / t + + return itr + + +def bandpass(eeg, sfreq, Wp, Ws): + """Filter bank design for decomposing EEG data into sub-band components. + + Parameters + ---------- + eeg : np.array, shape=(n_samples, n_chans[, n_trials]) + Training data. + sfreq : int + Sampling frequency of the data. + Wp : 2-tuple + Passband for Chebyshev filter. + Ws : 2-tuple + Stopband for Chebyshev filter. + + Returns + ------- + y: np.array, shape=(n_trials, n_chans, n_samples) + Sub-band components decomposed by a filter bank. + + See Also + -------- + scipy.signal.cheb1ord : + Chebyshev type I filter order selection. + + """ + # Chebyshev type I filter order selection. + N, Wn = cheb1ord(Wp, Ws, 3, 40, fs=sfreq) + + # Chebyshev type I filter design + B, A = cheby1(N, 0.5, Wn, btype="bandpass", fs=sfreq) + + # the arguments 'axis=0, padtype='odd', padlen=3*(max(len(B),len(A))-1)' + # correspond to Matlab filtfilt : https://dsp.stackexchange.com/a/47945 + y = filtfilt(B, A, eeg, axis=0, padtype='odd', + padlen=3 * (max(len(B), len(A)) - 1)) + return y diff --git a/setup.cfg b/setup.cfg index 4cd146d5..300b0b8e 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,5 +1,5 @@ [flake8] -exclude = __init__.py,*externals*,constants.py,fixes.py,*examples* +exclude = __init__.py,*externals*,constants.py,fixes.py ignore = E241,E305,W504 [pydocstyle] diff --git a/tests/data/trcadata.mat b/tests/data/trcadata.mat new file mode 100644 index 0000000000000000000000000000000000000000..a5f2d4e084587f6eeade02ecbeb32ab509e37783 GIT binary patch literal 20987516 zcma%?Lzg8Cj6l0=+qP}nwr!iQY}@Rz)n(hZZQJwBV)m0&PJTd6l6xa3q$(>UOvpyh zNGKKe<=k_*dUgbx(6X zd)m*>&fNJ}ODxvoLhQaAzII2u;11Z2TKoUS3YGr;7MGd)3C8qrBuYPx2-r72wvwakHW1t2MNDDk|iKBZ~pQZ{uG>SD2q8ZFz z8h*|cEKaER-EoDy)z4j+Zet78*tX2oVKoouq)Ij*7|JA?Dm50d$%Xz>sQ5k(%^^fP z|LGY0;5f0Q4Pm>mHLhAsa^>&Y9Gb5Q@PWR6x`T3wh(g91W^qd?F_BikTLOOX>cEC_ zbp5Dj5lurAb8BE$jy{_|o>ft}{eQ(*h!3DzX_TKUS8wZyIJp08GRS4q#3}|cp z#+DZ_7|m{UfcmsRD3OY-dc(TOG&cZo!bHelHe~+@ROrcJ)m_vanqAg0YcdAhMUao| z!p87Y9zJ{`$$=!3LhR3Rd?t}4>>b9RkL0m(>NLJo2b4y{#iTfz z-N+gf-|RYQdU+3t&t{3gt@p1l0rdY6&7+zOMH8xanB5XAo@5Y3SHA{+Que zV{^Cn2Yk)cjaA1x6)Y)xU6iEDtS#Nty!vWYGxF0ueqr_!ZV@qmdypG}C=h+>B|@Hp zaF%U$#x&V;Iw=lK0=JVR_-O%$+KY6N>U`Di6otm-|;dhUBCAVl97tYAUYAPEJvE zXb6C6qQzM;I8^qc|13@r6oG0U#a6sUtISFNwfi$_q4ZXsxj}Nw5L;`xg^y`~6cm(P z+RMc#e)Kaw9Cgk|kx_NVinX&ZF?fIqAV5KuZht1pxQoI%2o7~6+5H*5iELph%ZfG9 zZMGoE)sfMxm^13+Edg*!y*pX-kx^6HaB}f4patd##PL7Iw-*kMZygGJi`lmN0D|qO z$@3Z+9F^{ddGSlc*d}3517-jw!F=-|Z^gX}@DJy{h>p$gMPx}qx^dAx8f5EpQ7U*% zorRFz{CcKBSAPTbxwTXu{}6;MLz6?3)%Uy}Uc+Q24#Yf=i6!7i5iA-(0jXeE+pCsy0<{0o;{Ou;)lnDxqc^cJ zySy|KZ;u6PebtQLFCWKe%6yF)cy~e>OD0F$rFfZQc&yh)dW&VA_sbfc+G~V)u)Q6R z25`#qS!IvnV`c8efAaL%xw*}i;3iE;PrPF}&tvRj?W==RfE)yik2G(m@Rby-9zGjo zp;m6`bW)aAz(vj#Fy)j;+RxBe*BGWFZyV0f=U;`>WUF#oVnHx(zBi{S24rHvO+;Wb z8SM69+Bbm_aFhm`Me7~Uzmh(EeyrPcDf1R|J}i9|_EJW$@xxupZ}*WV7yIAiCW`y- z2E}_TVg=OZw$4pRf`a0z@GQ}?ApzE}(|O=4gvC^1!f2Buu16Uphw!xfX4oAs@8p>Z z{RN!2Ru{kePP^UCl!(!+IuOj7W2K2|_fk`2K*jHt-02~kq=pqa5W0z#aMy7&bnJ{y zCE`Y0UP1Gc=XMGY@Is-a`yK`u-48yY$)68ARUjaXelUX?Rvy)CKaPdEtb@>?TD2(n zDO`xd(hbbm4=fqe8$1+gbr2i;V%7QOkOLfwh5k0a$nk!WWSTEXinq=n@>4ZK7Jr0{ zo$$(IQIZ{qkELVDzSQgX^R4%JoKF)^X0lc5ruP@FreQ>0cbl$f4Tj*OekCwnmP%OZ z7Q-Vm3~KMm|wjzC;NUqF0-ceznHh@(dqY%n5# zuPb*f2EfIM4ZlkpZo9hLN2Nn;<=)oYPCW`Yb=}2kxC2uHN^7&y8C^aPZnDF?j~rN*PtOZpEx}1^b$ut{eF8L-$eu4 zeD~Z~adpk1yWe~&aRYA*r*>-%$%{Fcm7{TQB)_fiEEY=eVI=^MJGH%!!}9AiZItc_ z0eyC2=p0a!?JaJ-vbAC%n{tGi*$g@L?-dN*Op{`fM(Y~H)YTL9mT`2%P*?LNnP^`y_T_+3C&0( ztPSDjyyST!;$C!UoP*<7Xv5@J!8m}ta<5zDS0Ne;sp=igczXP7FXJc#X%JT5Cz$(W zo-|suTb5rJK}FO~yz-hjex-gnFpx7*C+Ms$?(<$zhM4$GpHgj`EYLrON+0nuQOGmR zIjcICduz*Y%>pohTHCtrJ#o?v%FRH2aC|Nhh;y9FuBAsp$o1Pb-E-<&{!{uTy+EYS zsf%;&L7{_soqHtO+ZK6OZMfA!n1oR->`%AG%fMa8gnPIxCQOWQi$V`{;nYZ8OCvckMyDeCldJU1|wTmzn8mcg^`jPpn}b zQ(x!#K7;?gPLX4H*23fxE(Uoa&)S0;nKudy_Z+fPj*=nnEtZ%_gQjorR3EW+61%XX zp_3}4)7GJ|w8Pmg*EOTwFt+LosX#>i4*ByHcnn6b`Ncg7N=K9IJzF|#G?6HDM{fnY z`DoJ~M9v@Bly~OJ>*j;fnn52m=8Px^t!}}&o5nOqtnw&bF$ub&(;bg_PzYKgbXkh_ z&Xdq-cW%C{#Qzo{PkwP-Oq8h9kwg+B^VZN>jG?QghN#wdpD{`XeV~!KrR64$s>h>D zmAN<^J=W0#BOBiq_`M&QBL~c&Yz~89N3%76Pqux|?-PKaCh^oV9)DL9QqB-Jh)UI# zjPeoAmjT0>MkR%z3L5xmsa`H#mfkS-B2|sKlRq^>}A6-By~Kw)@5+ zzL(@;uGBFN0o7?>9u~%GqB!p5MA8ElS$oM}~+-fwU*{>{q$0 zf(-d{;@y)eCjRp1)Amg}9>LSH%+DQ+FmV(P*+OccpwokFSY7}OjnCk)Qq~ZBa`-!; zg{E;Yf^1w#z&|JAWkatZ3xvhit7I2L!YC3o8Z6MBwWwJ%oh>!&CiYg)T9F~pgrb%1 zZ?^d+9W-~-Z98bk>gvCH3@TSc<$L{89~Br(ouy$DgJFa!)IeGVu@I-~%J4>xacYA6 zK|DS>nJ5&nQHp6LfugU99)-rn;Y>_P+c~{46z=>*fNLbZFLJwNKW$-<0%Vk@wlnw| z*(VNErHfBoUtaOvfbdX!uRJ*Q50o6W6Ah`DW;Xb%D(&}C3U{DLTQ@xl3)1oIqAIq3)k7{RK;M5)*u!G-H#sPT> z!5}cP6gSSj#UCjMJ)_wl#KEAo0?{>CL`#OV0`H?!VHbvsDjp5Q1tf0NKG}P&&88=d|VX51RFb(y+7S6iX%)nO;jg z>yjdZjMaSi^xn<#EBwH*CoGB z69f@M3(tF+-G8$9w$xF_lG7cSjJ-kSWMU3p z)ihMO@}j;Vqu&?>APVKfA5UwhJeknfb({X`r-ucUgv{+6(O+>+KK_2`^&ZFy*m;#o zhai$gl(l{=1Hbg~2i2bq#jNf6U-JBu41$hsz~bOWb~fP3cV>=;@CI;VuTs&Qm=n;; zAUjb%n0+VOjtxgGkmh<=QQQaeOuzTI+~KOb2dB|ck0U_f9PAjU#v^RzRFk|81|j*;VKmp zTG`L0X_=R;_hFaZ%=W>$%4pfWB?vpmaMIy^Q? z^pD|WHPLsIgv-VG#H{TkjIy28tbmDzB1(FRo2n&7ewf+ZO!Nx{NvgCqh@E>w5?|~n z*pe2m4APUwJN%0HG57Ewv(j0WkO#Z&WI!lP%%o{?@2a=>DsH@F#BhBXdVNrUB<_#O zJM-AVm7H@alPbh#iE^w|tVqQjnIMsQYJ)dd}u>WjX7<&sH2;#KGN!9 zPN&oAbw%CeaTzDnmqIRodQVYgQa9ux7|`)e$Kdb4*@)yv7AN77=~9Lvj_BD@CYlz3 zH1=^LN$IgYAd{|HCL1U(=U-&B4bROFq%gd$V0%_~3 zsZv<*K2qG=n9@S(K&z(Qs^sn!2Ib;O;)wMHoqN_Wf0#MYRWxC$U;LpaI4J*GV92G5 ze52%(aM3nX!ZVWyxt5m_UF^;1>!U=N6~MXD5tC`l|JFyy>*ST}@B&nau_9Q1NFUBK z9*Fr2LF#ARG2$Cw0Z3oU))&{>@MxLvwl>9tiEfrt`%#P?;Vlo`!y| zys`T6KrxqZu%?y|hWvYJB$@SraJi$iI=4k4$PBdt@BYSMI;V?pOFU&!o&c@gT7TA{ z^^EEs`5}rs+pOf-&p+5t$p;(S11*Xr`C7Y|uuO`&4Y|>~%7X9aS=SOPSVTKv>Yjb# z;COq^n<*bm@iZC!q+z*9!c%NF0&e1#<>oKfRq6-W`URWXc~038qU%@C@b-};R+lC( zxn(Ix;i?~~DwUuU8i)u`#bfr5ItygqTXa`kE8qM~Sp;UvuVA^1B*eowKEu_H;LQwn z3RW|r>2C^P$2}`-?0<)xZVzU};aO`d5bK}v{+PUiyg*k^nOJ33Hl9v{B`Wf6b3<#6 zbiEsVXW1cx9)JN}VbE|gEbN(@cInj_qcUOn@3er&H_e}hmnn;1l|uENA;Rm|m5aI& z%}^XD71)&-62Sp23}{UbRFs#MQVX|r)s}^tX1hBku!>l2wpk^I2=>UT4-#ZCU_F}c z^Pj0KeX8()on!;#1jK4K0JjlhzC)(;nB7M>!wJEH)0zF?*{ou&9n8CtMC(_O0kZ(| zTUFh`U-_^lP(Pq+9Y0JzhbNs)9$QInS?LMiO@|(Qy*<_zioUJeq`3!c^$5yz+W1H_ z;N5iV(sU>Ph`*%X7jKTlGwOVLk4ntE<@6|Ku*iWEMvKtbduqOh0)Cds?0A~`N&;`! z#w79{43Q48fzUzSK7$P7F?-z0br$YJ)%&#k{GKN!LY#R^WF*a7@i=Ujm#=cQd#3q7 zfpcO`1#6SqN6YpHTH%l0lkVWL2MHPkX&Iv5HuI zL0?hY&WF#efxFMF(&{F`&&9XQzjCyD^Vzxp5F39ZsWQ7;?NN5)=x*?H6>Q_V;*ak& z<}z;IG>%%?Uipxg_|@c2dMMzRlAig^+^|Od`?cQRcrd-Dv>L-=_d^o1@nERpseqrj zdE*3F--eP6@Ny#7$tJnN^aJSawb(hY~+-KH)692$Aj$)Q-?xS-*af=PyBxYoxINoPEI>3!PhiU z?-tE_JFB(3Hc+;A&f{Ua0&B71xdl1Pj+w*CJ6p20UGw4VUbzX+-3&l=iySD+P}qh} z{R7L#piqB*$n7wh68cll>&J{7u(M2OY2~9##D{zHD2H`@!VUH`wK>ygU}=?_Y;PFlO!YXjWa-h<*m!N}sCveBbwN#+UUjB{U0Urg^;| z|MFUfqiP^XGzN)&k6)>6Dt+b8I-NgvnyOmXU5(*P!bBcoa81#WA%Z9j9jS((`m$L5 zrZit~9thxm`I4q1aQQM#d+m#c9x=e}^Grk}xaiYQ=URanBw(-O<@W{*cOPa%up@`9 z)p*LTrG~g~{3dGn(fvV+2bqyiwZE)Q4U4$pHE(Od%5r|G=hW2G9{c)wW$~Yro;Pe( zgm2as#x;1Cb-0~>J|#NkDqy)Z+%EAes_#f}Y57`(ezEk%mt9%8cPUix1wc^Se@3`w zDzcETk!~Q#{Ke|0_%TEu_5l471e`>2>?w2%t&J?Svd;UtcDQqg^ z)za{~Jx8L6kHGAlTRT)R9?bOkc0>g53oJGeVM7UtcFBp!w}USD89HDM5P@SQQc>uW zMig5WyA9D3MUY2|s~vLt-ua&Ee$IkpX|}pt#95gNF|4r=w3k8_O)=!VqUux#-sWK% z*L4N@-F7SO>BDfc_O8o776=ay9Tgf^09^&#)Q&Ik90YKzM$;J$wheCC2lmNZIbzCu z&7Jgy(y{YjGV<^SjiS@65KLKs>5Kwq*q(9A9hoS0wtRr2U!?afI-9m{^|XuW-9_(A zExp~?av*ba7=5J*TJy%fig7U&#eLxEW6)s+9NvZ079_)(mKpPa zn=B`Kr3Y(#l(HHEIm9-u56>BpT;p}Tpw874E#$8g1-mo#cnxllK*mNMz^JNAE#;?E zKMq92RPSbMz`4%_pG-|sk1$&~axA$TY8px0W3=JVd?yR?V5drymack$BlS{l|*0NOT!op?%c>d68t- zP4b+IMcK)E?sEh<93SnWXqsd!II3rCZIb{t?4Ke%cq|6L{rZ?Q6Vrv~q_wG)U7Hy2 zOXrK}`fHc#w~@WHkee&}EThx|om(xbu$Zhyq;<6b&{LfI9W5$2f!CDUI z)m(azvDizZu@X_Y7WXPIDlOEHM(=3lcCoE?xVd__2aYxtag4#J2rhoiW%ifv2-r0r zL|^o`k7BD*q!Sgr;nVhFmZAb>Efzuocz1n1AmytrSx5vVM}?o?5r(wiKbNSpf(WRQBR9LUUNP!Hlpf6pAc({n1#?T@y-=o8#e^ zb?l4-%qsNuwzaJM9iWK;Q_~{?JUSkSry<9etrEx^(>ynR*OvdQmK56(#D+x#1=uiXYAS>c|CM9Ij26i|{^5dtQ&Ka?F z`L-Lt)?O>P_{ed@Yiu_y{b;;{{1W`f#D9Nb(>-{X?J&cy`QlvSykcB%eauLa{S_o( z<`cq{pep_h+v0CxR~?Y6})`MBHqLlvx``**{F9I0~T9%5iz=AVf}@Kz;P#3kcPvf||xg7GYANs(xvpokKqFR{NmB zA{#OaC9ZIltfbYSR<;n~LJrF#7+PT>vnmsy@@mTHPi#MVrWAUjrjykV4M&I#L&MZd z*a8mP#s6t(F%P=7$k(FPLjfvW>kN)YIXIJI9kqYloH%>)hCpDFM+WoOS+NUZNxDHgUM_1+FginJl6D4Z$$tR z*j12X#Y_%L!(YC?IdO`KOyzXJgWFI3nu?0HUCvM)jN_2nHYR*xddjp8p~UV043Kw8 zLjL$d-4Qa>7p2Q<*;ND`u2IATiR2aA*3a4dJ}UA-iI>!R^O}(LH=gyyn*7%oJ$D>J zUj>*lyg&Emexk6dojoM=8^bV%@R;@{dj!KAY2?10SL}5^b3n~g2B7Q(*4BQUy8B)> zctUd%jsoWvfuG}3edfB1Jja8L9u&KwIrI%a&EUkwi8R)d|j(9y_<=_ zH1#b`blsHOyv97d43rrDzoEU&gu51c!C8QLIPy zpcoPvf!Q>D2*GmO<1S}C+-CQ$jd@|5^6Ps$GCa7-^R8sDvoo_I!BSKU3rjIWv!TF}&ri5CcXwU&H%ib|`Yse+lg=jD4KedRc0%`!m zOa1c(+cR&?xpQxg)S03*2C(zP?PZkdSV5c$9iw4s3H}hJ{A~Vv=dm{$H z9@IY=SwlKzK)J?Y&G5$ze;`Auj?}xS%=vT&ikk69))TFzqA*fQ;RIhtiIZY`^`%P- z1WE?JH|I?OeV*mjgNimeK{`%^_6lnVq{H06z_my?XDtr{H zUo5o7-XZpRIK=Saog>AT0eE8NnRGwhc1$Y~ga=a+tV(-pqdzMw!Yb~eJsypyZQ`m> z_J;IFGrw>WoiXUz=Fa=(R|=Lbhaf&3L34xzoDRTinGPBoY3HN5TP1k8-@l3?B=eRU zzX^|Fw;urMQHILEjPvo@`nFMg?ykF5Rd54ZhbZ}3nN|~jH`x8uKR`3-O!(^GSaz?; zVE2Ir6jh{8*uRn*&Y;K_ssB-c|+liRs1 zx~MtKuGbKC+c5#FD0(4w-C}6|8iJb54b>?Fo1Rq;ASLmoVHy8zJZxBTKZbXH%lbDl zP2}u9o9LaKz4rkjdl1EgC;mz*C8+cg$r=vFR5K0|D;vvduog5sH$NgR=weOaex6<6 zE)2W)u)-n2k>ls%9isu6aG{g;Y<)Ot;J5kn@9Lb80o3idUB)OcDz6(|UQ{2~NeD-H zr3I7gq4Dn>C<`Wdd-?Xo0XTf6-*#e~QcQPECz<^%xBt1rxpqziZ(A2)t2-(Rc(?`g zA|PeWaQ2@B$$=5y8%8(nL}#eFZhonbp57#Is`sY6bS2&1h59+b++5YL*w0>!+7I;c zTX3A!;WxzYdQZy!FWn*gE>?;T{9@ZHZjI!; zsi9oMjMuWkD$^jK4jZnFHGo7=kJu7`k*C$;uy2Qo)#A7WxUDS&sv50W?|=2Z>y$L4 zukD7qmA-t6#VHB!k{Y>|KowH;ziXc+(V^z8PKV0YyXced-LSlVL<*Zqum07{hB{m4 zIkaXtP(zOD!Cdi@mNDDTO8^#RA?>J@Q@PWTqBtK^W+pgbMBKA#- zXcm~U|IVWi0yhwlY@kXPbO877sqQ2T35toT&5Dc&%7%IO4k6GWo2q`vvx}RBMVfpf2n$&koXlK|Ux!4}FN7xeDe_IELa;bzS`bh!?uy-2L}qH@s& zZFyv#Z-grTZ=2ZO9;Jdw$=>x#zAhgl`A1vgkbsBxn=#xoMFcZJLqS2-iz4{^Ds(x! zNa{&EA>g?Xjm5R$jR{RZCrCRaw)j#*yFod7YDc3KRDJ*S54Or`i?;D<)o=;ce7Nva zXIK{6=@OfECUWWk5WPUFiBTey?EGCHS4B;#t9dTtI|V!8-1};m{-5#Y1`l$wAxoLD zSVoz_hrB7m;`L%G0OhSN4E9Y2>G-SDsH}>M2b=n0H%m~}H&y1BiD1Q42;(zj4@=Kx z^1D8|69mppowi9wr9tsW+y3tpBRjh59T?ypB7B;6NEZ|t_9)>Llk7>gl2P-G8Lwtr zL|kOzD|0%Gw!q4fT;$^ty;qa5-|X2b!mO7qWqh(|4t{QeJ7 ziNQ7ZQG!Fj3s1bL5~nwXMVvZOFzhWJO=lJ%kWP;9)$q1>Ln90+T~!Iu6{GS+2Mc{L zssfz^D|xNaqd&IS6h!7c51h`RAB}W9`_UDacRmhhvatNGEgJZ9TXvkytk(QhQQ1e- z5AaXpLsH%h6mm{F794t2e!l@Fyq?@g1>`^WS=-(=24fU%SGr$hds9d4&Y6p4>a^#h z_iCds&H+39-KcXT5G#f%51NVKwbZ`!kYHqCM&1&;&N6EnX>LRFp~rX#ZG{#ZX?kdq z?|y?H3D1qY+MF;vE~SexaYC<{9?m#BP_WXp-bK;y9oAG$a}>pd*Eb>*E}^G8L@#kB zeXtvhA^gOx$w1b-KNs78D9Sg7rRfM#g0J%b+#$H-0G8s}*PMv_k+FL`_*ny%`fz)T ziEQN*qB(T_iuBrrhUzq(T{b3(_f0z3!m14KsS1iVe>PlPZMi(OA4`Y zPH98e%=Q*QY$X4ZDSjY{2`Wcq_L;mWIf1M>YG({7jMEXxGD3*kuKg0HIrZDj9S$(G z0{J?ct|k9+AVnORB=%#h+)2pLZoN2EI1axDd#W7&$3Z%v`H5KC66-ym9|xqSj&91> zXpEG~&raZM{0L**<(J~IZw|mKH_wc+8|e z@fubPZ}5lTHk(g3AFO-58v>qV#|(RXRC!c1EQ6jutM}Q3;p+jTBOb&Si3CF;n8&Oz z@Q5m4&RbL5`pu!sFZXK1V_9CF*;*(>tWsa1HwB0x-6#I&?zLzFdhh#-Bw*+X?cY}j zn^x`tM=luQ%GCAv<2au_3KKz^p_*XaqhbF1Z^D?K4I#77DdUj>$%1m=EQRxu(DTgs z<2)seeU2+8c)11ev{8gNcP{zltl&O+K0Y}Zd<3~k7iXzpRVHq*rA*idJr72`XZ0bI=|AlfcQIK(yvT zZQ-`}f+`}X4wI{O87yi1gEum(Q;vMoJ>+jZECzvQ>n%V;(EKJ0Jkffp54j=czyx$} zwTk?*sZ)S`N_Td@ry_ zNbc!K3ju!;Yt;VbvBTRqt6*To9@Q)UzaWKdYj-XJAuz)5^!@S^Med%%KHs$AdOOl+YG4 zs+Ffs<&i2LJDJN2N9J1IoLFO=XcV3^0Bp~o3kzFIm28vdwP!=UpT|B4eKNoCZPtp0 zhZ^6Q)SI{aRC3{g){9B@c6%YiIuZc+h5W5mdD=8CA*x=!Ov#ZFdB)!}Gz8YW)%6MA zRWYy{P2vuHLOrMuNTSky>9XdNPQh(nBV7JbplvVs344?1n-TxX2ZZ(e@v<|m1Gl-Y0yGfAN!ND!yqfVTFMHdo1tll^Rue^LkR~(c7T!pm{jQaA*CFI0jH-> z-cFG+42O5Mbaul{^^{JX$L7_HaOtt-^cBX&P~@W9{Gf*t+8oa_!`J%g&}4Yr-J$ChLkWFGMLm7P5D=YZy{-i zB_Qw^lY8Vevy21c-ZG4LMC51*Al!oh!) z7_4b%J2_p2D^SA4i@G))!hJ`C7x6?S^PdE<&%$%O_;Hx zq?hrX46j@Q41GlM0g=XE6uT;>&(}uV7!PJ$=x(<=9=z*qU8aU2630&&j;=Dk$oD$42sBxggQm z@wE>)?ycf&-U3y6PzDzX-&EJOzTG0CM=;NH)Z8IRcg|niZ z*d==gdnm-25I0wHQu zQSk(uS5$%P79z(Z_dtB%`X!x7l;l{(ARM2$qo>J#KLtu1&B_Fv{QPh<#S!oKCT!#) zV7s$ia%n?e!W`-RmrMIe^5&a2J6Qi>&Hw7R+F%3RsOYyw6mh&t&mymRSfUv`O13MM zmV^!7#KJn-$0cX$cxk%Sak0)K|LJ1mk35Nc_-Sx{V=PAH-hE?92ovi%ylHsU2gCjz zCKFho_2*lX`*=d2AXZ;SDA07(-}sISC}%BL-Us~Vbsdp)(H99+n2I19)dHVx*BlkH z^f!ZaX%hK-!n_-!LyIP!qRLcx5e4)|?0mT5;qVstG>Qws>jU5>ZnaZ_R?WKn^7NqI zR@kyEHH$$tJ*V`43x}h<>DnmxP($`6MWTZ}QZ)-B#T4C4$RIYI!jQWwt3d_tO-At< zOg+E}`+c%DwTCWP=hLHtZWn`c$`jY09@h%M7km~B zUH3>K8vkt*eS5{mjl|=d)F>YXmS*Ru;kak6xorb5q{E;n*IT0Y*_lFT|Jj|(aSH;a z(q@KEY7G<^vvjOp-2uybEC5F(ZA# zh9gnmh436hHz&-8eyxFY>3$Xbs4N*vc_H7%L7xAr_J{6#B3!feRzjJlAqY>!%NbIW zfavjFU9TB$IJxIie4QEVmC5ns!$$dT8tRm<$6Wi=YwwwFg#z};k4LT*!Lu3+>-*j9 zRLjSF9-&mg&;fPR!c)rNT)ZS?B@K!DB)h>e=q7HilV=z9mqth-a`@9%cZ$N@qdp;{ zHBgebUJp)uJZh%@EGsSS|Hc1RjW%Ng9HIlAx+g1;j0Q-4!I|w7rNrSoIMO_S< zo$`0GV?PpGRp_alx6K`15hi_cN_Q?$+leYOSd?x602b9f4~7CB-aWkJ6F6DB90@lvPT(D5wR?)=kvqHE^NwkcW|6yF z;Ced}bgKNs%5zKq>5XExwQ|HCy!qB-Q6u_uNNW%*K8@_ZHFx6)U_na}tp>+DMXxXM z+0|0N__*_ro!Xv^t3Dg~rm13^9JsTowtUzpK|YwqVNb**378TVdV*j;H+Ax7@-}>B zT|r(NEtqgFNak{_jUz(x%&qyS*!^xVHAXD%7>G8K?pD)OnRbUz&09HW=BH`z8Mg?S z22Hptm-e*k12XjnyX{Qkx)MJcZ~nfoS$kR?o)jK|)qSnldudf6dHa|7yOi$5ZyeGA zQ5F#+##^qQ36{zt68c@GEC+_<*>Td(`OgHJ19kh^_JRfpf2}BLw*(b(YG=}lod)(j ziKFXNY=Wuf*ly_B@9U20acR*>1{k%@mUH@hUr6@l-0>uUb;<%%qWR7-@QY$w^0i+- z=XD|^!H*QpF4OTpZ9G*11UEUczcu|X>Mgg!zuP7wtadGcLk_J?w|V{NzY$ym-71oW z=fsI?odx-yDl~M?Q*+x9E{=G>~yGV`l)np0qP2tWu{%$4yLPUJWOmKrWg zwpJ!dtg8GpW3e7p_v?DQkg2b)RgWEoW{oD+Tor#uW~L~p(dT!rM@KmO-HD1j%m9@4 z8iMS4BgqMdpv+tcgB#f`nNcGLO~^j?Bfkgy$*5Cy*wSLXWN%$vH&7-4N=|HwXQ!4- z@bW))0kw)o8>PBdZ}W(+$dLO>MBHC`nb3l zrDO{*yf;gYp_2XFG8tLGgTA_fSC94^4~x_4;Tkk6fJ@x8$6DU~iM5b$kelY_iyCLc zi=Hq^6l*@j6?r&;j%Yn~e%mfL`XXrxYzNJRTE2zi=$)#DZsj@z{^tr=+?z-w0elHK zAcLC!n(q}=LS!X{{wWfaJnq!Z)AKGr8hLF$RJ#&{wZ=s}Ru(@Q{vDvtNf0;vome6e zawX+*D)2H22BufZjS{yd@suu@y-^Wnmee|gLrG>q^R3lc9V<@=h-ELt`PBpM;z`v# zx&&v6Z_gE;Ei+Ymb?Uc;DFW@lteHTw*^M8*u&rnNNd;~%7En+|gpe)}y@}tD4|&)` zp;B7SQ2YsATv{68XQgtU8URA2(un8s35;B&8I-yX!B|P@If$y zAv%r9WY0)sAUM8x$H?aUs@ux&VaD+;y?(QC4Tt~I_|^}U`a!&_)2+imQzk<@E$Wm}DE+eVLS6gpcs}hzc z8l~4d^7UjhA5w!Bzpjg7Ih`a)jt|n0vm0;>R=IO`4*=Ds@AXJ~MCe z)5Mia0;(pfFYBv7S5V8xvL2QRV?VapU0DhiL7S^p_3pk|@CImSKpGr&PL4v012Y9t z{x#e5y}~DBO)(i-r-^s{h@m*HNI1cbsViTS_hc-#!~XIM@fKxQ2J_eRA^437HkR$* zwvf+r&N4}W&OndfNX&Bbn-J12+jt(mWWKrtM4qo36b4!F>R;i+F%^DQYPZijdLQbi z4hN^(p1P2{odX*R{;xQN-kXR6$t2;RPqDJ|Q?JCU%*i66%q+_qtux^4DBL%_FjxU@ zG;_X^UbhLT0+`uB<2-KjeS)3i%scG`AY|t4Hw{ifmPU*D7{S1~Uh(UAmov&rQhj;7 zm?-Z>JWdf9OwIeZ$EV0uQObKs>}BmVp|(N!R&@}hr|YO(&gW6`aKBsrNyoT7WM5}s zpmpQ1kNwa4?O;+(9o_C-#9~!A04TolgnnM+$e@qfhX2<%hOa{JpI zm843Esgbt=Avjy+hDi{v8P^M*>}GBa1*8xGf;*|=ptPz+s`zUZN0rM@t8vb1=7554bT zI_*ymQNKp}%gG-x$e7m}on#mSu7XxGrpIfjb*$XFtPZAhi@iUNQG%2rFKGkLlu#Dg zHP+Wm8gm~O_L!-#4>7aRB38Dbbf!2fHwn&7$`x}B&qgd#DVK1%USoXsGCUVs-3qfd zT~_sLrJpY05M=>>xjsh7H*MA;)~i)bNKD36s|n7JEJ`u`#;Md;Q|3qz2z*o?{`m`I zEEZh5t9g#D%vj#WOYu(z6#lx^t2mVuA&|xOpB=q&0T$lsVFon}19iOU>Xq$)bddXL zFa-p6f9typB@G(<(<7Y!as?aWO3F8~Wb6n#_=Ca&?JIsT0BZ*!W6 z(ns8%iiyIA-ManuKz7)kR=0gjS{`L`HoQ42v4ltVmCa`=h(bF4N#pErDznn{MgFNzPSc=1BYK1aU_DFzS%T}D@X8K5m#lixW{hj&w#chQP?QK8Q5 zrDiG<6T0Sd>cB8RWHzYyl|+i6?EI9-Jnv1wD_QVrXsib-zF6|{MTY?L4$JFqXO~1h zE*2lwzYC#LH?ACgG_i)4riP{YX}b^gg%UR zrlFrXztm$MP=QuzM2#=eL6@>}_)2A0v9f{xB)Ul2=vS4c!yghKiatSAaxT~-TW>rAY<_|2T*c3B8@EQ{1nB*Tr#TA6klJMeVt z_0%ww5t(?^a=MX@OdtDK9uOy>+h0tcI$RJ0pN=2GRYT$^%zZ3-gMs z{E34T?Fo<~+;-ANPNAIKW9k=AW*g zwmvNmzprcj_g0mJK79On%zGaJo|bJ>)Yl!ww%tk1cv9Dd_YUrSO{0sVFJ0*OpWC94 zSuRt2kB0)ajmL-l{_sJ&%O}=1>|^2~{Xu-vH~Cm+B^emoqmKCDtk8NW*!=fa z&W(I*r_SZP!z~LE!;|ef;BT`o0sLHXBugA#72J^k6S zpi&gXI))A`H>P7vaidB~;gVpNmU%lQP8_~c*$>q07lZ4>Hve&F4)jeDbL8jg5Th6N)QH}+7Aszj4>=H{kzO_DGJ*4pRZ$6!zdZR2EQ8Q zN+7y~(0b%=VU$~mrFX8d!cU$-j;I}M$YuY&nrc1@#NFO>A>V!!TluFn!0N*f-hJZn z+p`GZVLWG<^-L8ORbMLa9H4{h#V?(M6LffTMEg>E01YjllF`#X;dl zjBV$bG)nD@@wxS04mP#ja{VUOgAXp>Fy-VWqEa{O*Nd^DX#acvQn~9xSio?byBvuJ z{6u$YT+Cg;q~+fd2iTOLIe{0yYP^E)S9zL!Sy2>FRtC+`MFx!!+~OijG~uO)Q)BLS z4MfNfTMlHhA_L99<{4ZVWtC2(lz&;po-{psd}vA+h0qRdNgAVogF^qdabGe>w7K6PibJ>!Fk z_zZ_@hH9`{VIH{Ytt<#`{w2z(A`2mmcu~a%;^@@UX7-1NMWEJmvCh3)415S3EEl?^ z(EX%i0hKO_AbIjoQv|6EAN?e3v1~*^@GfXLutFH_MK`Uzdp3cw?#(ja3P z@b>CHK2V@K#;^sgV@rxp((SKyV49)JPwOA3Bh@J_X?JZZ#N^5LxbjiK@Xb)3=G9f4 z-mGlxpff0$ItcHqC*V|=uL3cm(|OQzn6z?|J&w5-!2c2&Qt z_|+CUDAtSJoc2ipo~hZs{%0x&ThtD3&XpEHac0G>hN}dKyxDm3*@-@FAiz zeOJc?H~#6AnNR=AXg*k=z7nL*38D`M{YEar%lIpQo%s3ECG3T=;I}cWS=>P_@s0jd zRyZ{Eo-WU$4lpQ5=lDlJZ=@`=zNAXR+1g_Y?Ba4Lc57bNnv4QGl3B@eq^YABZmta7 zARZ7z0lMai94P4g*(dvyyKv6&`P{f^UG$Dk@`RL<99*Q*LWHHL@K?ItZ%tEetoPbbKsXPiFRPqmvyulN{xcSie$vQPTISE`p*~zyCm~mka}l@q zHvB*xVFAWD#hZRg63`D~oRQFNWP-VnA@NSb^ddZ4^;3pJ$F!fL<2; zfc-cL(T}v15zSk52rL4rRFa+U*tkM-GNY{xRH~z ziPjFaHC&fV=F7R?qVO{KOyK<%I--{E#z#FA%skqsuqBdjd8 zlOfEyC~b?=F`ZD6LnXmT!~mwflGfe@q5i; z$bN2=FMpj36xACj>VqPBzzVh8hh;o&MyOKeR{?j(VB;-wuwj&Y#<0>3tZ$%)bpD#ORGD*&(gTzYuS z2JZ@XHtpHkhG9qEUK<)yL`+=+X_bw0+0SPhgxC=f>&`97q=PJo?MJ=E^KoIs4*my{ zeCP1SQaDlN@CWZ&5U)P2N`NxCE-^=UAw=50GwpM-DkKf-cuXoOfMvyD)%%C#;r{L} zYax3vWLo6uwO57%#?*UQJbti)_V2vRF{Tu$1EMK7c4^HX>5 z0Kqu@V909)WY4jsS|?N;p85jQhN+JnyZdAW-w;4i?1f*gwg{rTsN`H@%wo2){;$1_ zHY4X;wV?J`CZ6$Oj+9lw5299!MoL_r*kk<(k*+`?WKjDApJgaRm4raCz%dd;=acMa zV|bxDJTFRsR}D#73vk6+F%T))xqp_(fhv^F7~)?!QC^oFdF&ez`Rges4PD?sx8#N7 zQkHq3z-!ng?xFzP>rQt(Y%7MG=9N;kHp_yO@RC^kKRWzmG!w#bC9pfua*C5Fj`XES zyGtwj@i~*_*>d?$c-mvu8Ig24inQ5#zfnUGSTktJ-(&w_X;g2e$D9)AmBv2%9fbn0 znSIBdk)8!?%a@&h1vd#pkU4*{d+<0ek6#R`wkH9RktNZ`v527$;pv$!OZZh@H;tX0 z0>~>u>3PE=K`?#b`E64c8{)Yg`@;Ug3T~>Ln-S&l7mt7Uo9%{P2ezFUxcD)X4GG`8 zLaWtegPpH*oZQ}XLHMPHZ%uq8@O`hLyF!}4-yUbE{V1Km+j8hn^H@Yd+xe*A*{MZb z{{HOUn`Hzvc4!>}cCO&u3qOW&TVY|*zwpl!#{M?Vrj~kGJiPK@8Y5e4_h3H4e z#__ppHObdaMd6k8(x~tl6HnO9nqZkI4x!Q}KFkI|@KW`x&*7v&ozYSo>gPh$KBvBh z>WV>+(2hfQACthLCuP-c{XO;|fIV@`hL1k>&fI+7l?JXouH`HFz1X(Vgoih(Wzj5k z_+iB-X(ZJ5>P@T15cY!QRmnka5imFLKKhAQ5z-XQIeicj^0|DAY@bR3U(Y?GQ`WL5 z``ue}u%klbxA=<|!U|{!CVI7v*^u_W1t|^@c~tZ)OugXr7#7o-EH}PW5H&EDP3dW5 zH1(B+{W|&-^H%OP;l56@O)zcMAjks%QxTTF*i35{F#k;|(M!ph(e2@iqO^*wb-^??1i?pr5PX&ndQxqR6BV ziIu(Iv2ymPpPMmGG`F%%>)D+S+$nY6L$@0;aI7IYb6iap(3u4D5cNT9mvQdpv^I5Q zbaXdbP-X>-?HXmJdlqp5A?ZQq26rCMy2(E`5(K&W{w32g5|sOtOGxb{!j3=eM}BSO zzHeKYjZQfQxX$&c=S~pew_)M%m4_?XbKlRwlf25{JmszMg`I-5e=y>{x`;y()7(IM zqt_Z2jb${-P?3^5WslP?2H1Aj>m1(T+PyTrn1M-AT^ET zIDv_u+n2KKl92>N>{yc>xkUoQC0V!f`wYM%k0ed>5>d+cTTOKyQsC)RFA*ay0iSn{ zxyH<`V*?8-)j?wtkSl!B#=lYm`gmXZAFCHYC-T!cU8VkFAsIa1gc;%}Om5`XS~v%^ zP3A^xzTtr1CE#12Mh9;CkC@APyy!*vX~&brBv3qTBSa0MBV^ci;A%1#1Sk*4PkS(s z)ZgR9FQ@wPMtjEH3*PLAl&kIO7rlV-l0})zEe!PPh1ou9&SCsxlB|0`{QzFTZ)#ti zNd)x~!@qnZllbz{v*|T>}o?;N3f!PgV+bo786(j?n^Nllz>w*1ao`Sjef#E z+x!f%Q;4;y%Nt{OuRg3tYuJw3_&DBLR7G?w? z+V+;Xy&xBSpdX1HEbqVsVw+F=M-f8hx1L&aT^Pl_gl9cIc!3H<*LG6OhJ?|zSGEB+ z4|76RDQEnmISq`8i$BJG)&s-vN9E1@e5ifmq`de96&i@H7dvkYz}Iu>l!=U<3*Py_sgp`iX(??_G&eq7JL> zRzU_?gPgyvF|jnq&nJa(5%i3HkLBlP6$Df41FqTC7%jmdv->kUw1iiG_|(#cRUZG; z6yvJ|qe9i?`wvM0O?ZRn*Z<+P^(_rZ5l_%?|ncON-0G4q*BhL@=vgn8{M=j&O(<6~|o>ez}wR79xQPoGw7 zbxqsdJ5L>$eVS*Bek@~V4izDN7#U;)g|x8G^hMHS!BdPoI3BIY{JtND@;(ZX81 zJ$V4r+wZvj#Ha!akiWTX-5`f{jI_IpHH~6tcmEa$C>z47<9%;>dd^{f4B3)%G$l|8 z)qEx*NPt^hwG`IF-0)_XomF#_7Gke7pTAVihYmYCB&yPz@so}E@ZM7ZwNGAY>Gu|Z z6Gs1@MFhOV;#b3>XfL?Y(!(7W){jVl{iuIuFCw9@DcF2r-HIyAb}Ujxjek9x47c%%|GjDb$2XImWPrqj!B<%1Kb%t{nMv^d&a- zMd@9V_h0-;d3*73ej*A;p#F*JAH-dx^INBUmax&HdB-s#6Ytf(M)sTGhYb4F3~V18 zI&>tJz48$|XkA=-s1Y%P2T1Xtp3&LJ(Zbh1UL=ZOqr*d#+u)T;P^pVz5t*iC6UcZY%DY!tYln^0_+mLi~%% z`J?PO_C9BG#d^jHE_XLA+~1u5a+%J<*#b=bteo|535^18Y%dxYxJ=_sCfm+>T#yF6 zOu4WccRqCBhF6bU9S!bK@{b#d2t(&%p?FuFPORlu25w=&4ds6VRTA1{!2qx4Z{66h z31aEIYTuNRM=jf0E1iz6njRpl+>wXA-}KzXU_Mkss|Lk=QV^cWQ^7dEflzMMna&MA z8Lc?Vb}CgHr9SyWAJHVD7kV#VeCAU@{LXLctyWgCD-q}BHJk{DjcKH4;3bH-{Qmo( zN}t6!?r`}Risj%3-mcCzcZwhu<|UK;-s6}bYmUU`S_15eI4iR3`v*^dp#ND}UjSx9 z$_^F;i^7=6rrCKZH8l0yNu=Cj6934ua|@?aGL|Mt^O<&U!PUQeJ-pbg2m^1mY)oHL z&|bgsL?Vp})RCR!;OPlmGn#EdFPj3!e|9@Op^_o@%a+a;WEyI+&w3=pqYDa^X|+N5 zOicIwjZe3CZ}3`c_iNJ$3Y>lX=y$LvFO+_gn6nU-2b0lm^R=%V{d||(S9%Q#Jh!@( zHea-HZ;_q7Db{pU`T3;XUKIu))ze}TUPPq6=~w2s$%aqt4Z9jP$$^qHYF}z7&f<=z ztp?9u)<$}jLZ>oDDCofPl6{w|R1vTK(Mk6i6c{@B<|GA*qywX`VMKUm7233RR}OslpskjsO0ci4F+M|`4f?~7j>E;D zSizb1Kfg+lQDF3yk|--#cvHE&g>StN^V?*cnRbo_(c8ZN5q4g}PjE)}s5dm?(`_pu zv*H?vP?M`2KK>Y2@3|Rt@-q!>tvWtD%)oop!CA5C}69tYRtAAiK)QVeQ^-a)I)`V%HA#q1ZVWhoq zFp9ldh9Re%bI%iQjd^H1Td>xXX(%+la9DP-n+(;|LCfHO1tH4}fABxLEK zD}}_GZyuW4Q-QH#*A7J}D58g)SMT;nuj8u@uTNBM`2W9(!vPwPHs-k7m!~3)1;xK+ z#y$0C;%+NHv-wU`;j!saKZiH`;)w2%{WqCJFt?u0@0yf_5ADkK^}k7=<>no>&`pPc z@{eq@d5UNi1-=AN5-1cM@;+xkg8u;k0RR6C*M}pOeG~`qYj3V`FPF=`M)uCGJU=QS zWmF=YR3uT!NF_>Eg;F6aLX=49MfDs>QTB|8jLeV~GG2efIp1?W-=w=Go?nyZ@CV#N zJjRL)aM{_v#lnpT4R+CDId$0~yZ{+$UJwKKt{XAG+Ft;S;HsJ#9ju} z4hc}^JcZG!Z@}aM3O~~Q&`e4t38P#Icl!OldffE-UglJfHM~qB(E)rUAg(`8<@nPD zT$A--D1D+G|H38keCPoM885x6nd6$p-=+oom@%oMlSl9N9QLNuHK@@r2?eu4{Iee))Nbh?-5xUxy*KA&E zfm^}kf7g{c5Z`@&UL|u`2oKMZoARQ2H)wf+0hR>=FoZn42QR6@Obji_s44G1b~VfvFGno3CG#g=)_T zlL|}-;oPm=P%i?Z?VPLk9*p9=V0HO5E(jOG3wjtb1mMrHr#dx%#L%O*a@zXi?Ql+` zG0D$<8JiHgys!~K0og^9-ioccC1|1m+f05qx~slor=JwWPp`%+iYh=zM5{*S!&Q7x zJE;F1O$ynbGGv)6ltrV9$3x;$rGW98g}-_xJ8~+*2d4j!z_sg(9eMvUhM5J=J8x%4 z`xHyrIj$?B$$|sz@3h&`i~Q5f1>4lo(eZVO&l~(eShgy;Mz6<{4(l}CcuYVlnf(s} zOvmtw0@Z6P<4d?>?ww+b(-O$7Dm?3-2osvCUXqnI=Z0q@k+<$0-UhiYbCiuHM)*&f z?-|=KG8A7tqg8W49yt!r*vj0V#tX9*OFoGTpxd99GzKUHq&}@?%=C5*?-UP7TbL9< z4~uA5TMf8jkgii*FTISvSMWS#C_aT9r_Pt2FcE+lT6v#BH7BapuAQ;>&<6R09xunM z8+f~Z_p|ft(r~Z$>1w~IH1r=i;r{o#Ae2?Y;ZA})e3B?1_aZ4kt>e|Y1sXpZnj|H8 z*oq-7={Hw}O-0dt#`LC%!>phu>?V*n*MfO-{c5kW7J)LRk7^FXKk!tg4-R&2{OF@Z zmcKVy2Av~Vog^He!p?BGJ^$3Dga(*s*dH?v_$t29L!t4*s`GkOwEHslbFKJm@E%^& zZk+#G$chUY_-@;^E1MssWt5Br&$MBgd`v;Az0elMqmy>$ z65)`b@4~Jpsz{WYcDLsHUyPGIn&59gi)|XOvwxRt!JSK%cLo&@;9X~lqgK2Cy7&5e zo5{u|-YJ@tR_CD%%y&D6IKOS(8|zm4=XV*QR>nDzq}+#{NLN+-8m0v=!Y(%n5*P8t zN%NMH85wxaLbW+?KnCb4tgXpjEFd4ZSMhzcJS5MohFHW>;MwsrA%{mO$VNn}?6%ny zo+55nZECrK&&88{$sgDu?i{DXPw5W)*F;#qqPZaSp6PQ-e%OZ}JiuYwGFOEgrk`OZ zC-9+$wiTzyY)0_Bqg*N&ML?0GsooY`%ed*Nc2Dt59&i&5o3yw+h7I&<&<<#d!F};T z!dGJkNFp1STqHBV(&RPL>ku{={EV^EtA(J)oH_L7_ zfTQ<<{Uw6%?TF0J3Wp)wLUgs@l{6Q`%zxunaURFj&jo2!El_}TpSNAxmJ3brD(Jk7 z76Yo(_+a*70=(BfUaA+gfcbZRntnOiiT(MfbWC z6Dx|@EMPAdMV@D>$RP)%r1#Zh6Szp4!kfo-%a}Az>rtbIX`G?{W;)@$C`_dtY_RM6 zg=>>sz1sw3(fZqrdY9|M@b;KXFsA}Dm|k&GW)W8a<&^ahZ!K2%rEm97=70cbHzy_A zB+Dbi&|DVBe_z*s&5dYQ;rmZw(x7oD=KAOOulW8}+uCSldGPfK zkqObDLS}Q#HA5l~TpUT)FwtQHPWBZ?>4y@?FXri7$_)ibY8?Dn@KOY*qUn#}j0|)a z&gP7N5rQ>?E0LvBi`bQpu7LNv{di1SkAf)+4{SOM%YRQ5hA`%ohT9u*NH6(!mTn0n z*tR?~@@nh=A`F&mBaFm-f9+vs-U z|ei_p%_>h#vJ%T43IMa95SQI$N8EvPA7V*Gb z@TVnypWjo@V#c0IBLFa9K#o&2+D!2whfLnc<5jd$rmRLY67l5 zKTt`)aizCHGiDAW>>7@?x$qH-=P`Pz%)<*|RYj))6iD#no}TZyx64@JNqe`Y{oiqv z(0*03QWy;1_4%18{=}1~l=2i>NYE~s>%Yf`6{tOvL++vV_~4JTyey!IOnI71hY>Sq zEJ(dGFB3yi&B4w=yiCyc{jbgyYfd<~J2YP_V-~k6VRhix+FOaSWYgnihV$$jG z!xGRvEsAaP17ujNht5QVb22s89L?9Rk}n^Vk%=zIVV|Vp{46}&z?Rp_|4Zj zs@yz^MIF1F$nD4rv-};M3Ik)9p6jvdBYTI&r&`}yKOe(c1H+)G3Q}Y9X~N%zxrC1SRrtI;(Jd)XbC@kUPX%yBOq$Z z&xo!?ZuqIy?-dZw1IY{prCztCV6S+1qjWnFB$>D8jgqCnze?S5l|c&qH=T{#sANzV z;V*9t6NN`IZA;OOTYu=^qRQzM5+pI#EtEuYAYwqo!8N)RXmy=%X(Wij&!cyqd(INk zj=L)d7FU?y2i8%^to;+GUD=Lgm4ITh1xUPZw$8 zLAT`ttk#15VcBlS+)@L@;6mB{hnX#8blAN{(Pxwgl|8Tx>vm*C(~)h0j6PrShWj^f zVy~vL9Sv{AG~MX1(^k|o`KlmnhXLm^K8)!2A*bF&Mlza=!^T(*1wf~Z&q!=h85L55 zbo9(AkbPz^wFC))VM5p;`K2Xn*8leIn7%*wb~O*}JYFJ-ukfhPY2-t+x|hDT^NV;7 zrQug$HQ?8kdnHN9@=#zAJ`qsP3kwz!F_OfO? zUBXO*Mvrbvt>9|&jh*LH6v6M>$rlYbcu?q<&F`kUa=`E@|NA>}aTqF8C>x(&aMA~Z*Ym9|JrLXW|G)tZkqd>oDXxlu`hii(fTP&bD0vhv-`_3Xf8 zh~FkwesVx^)LM7A8VN|5>vhW2(un(p^f4`r3vKW|cXe?P1nu!BiOL1s(A;sO(8sD5 zdlW78@wFE(EIBdh&@%txerncc@m;D&=Va$mYZ$_uuOvKnk5oZf(Lo>YOY_5Kr))me zWF8O7uXcB(&f*hSU&*kaFB`{xJP)@{el3ck z3QWYN%cOyI_*})hdPzUc4v%e zv%>{zO|cgjTCoCM#_XGuG;~DH!c?o9hA3I$F;1KUFiY*W)mZq0ekNmTp%U6Px}TdxpBw7^Do%+etD%K2DN{-1WZ2JR&M=oWg9$vS zHny>x#v{-4`1adu;CAsgMvv0F1_; zgcmRqkyM%5`t{T~+;7#P{Ub{!UYV(G!J*6r0h#w5>{z<6U6*}U3^pa8-jckeXu8!? z`;(6(^%-=00LW$wqiX;5r_kc6EKioDUQa>qDXc z3Lefw{#!_)0pIWY8XCpiuPe$rl|8i zX5@nKE{2N_1%==~m*I)b*Gy=lvmcok@Pl%Imh7VdEhw@)nm#2AU-!o=E;SNCD3yF@adZ&RycqCZG?5FwJdqUt!TJ{uI&3$NUnc{- z&&EdXN;f8;GN@n5E&=9p8@9coqxh5ACkF}Llekq%g2#uLdi?O)MN-i*K_J-J+I8+7 z!Z|f|I2_}pqi{Ysn9vhL7eAlt!4?QG_em&tjX@DobFa)&#`r*=Yd<+uMG#3>pzJ&4 zR4C$SPofN~p-2A}*-KxQ28^iF_mM>oW%VU>_Rom}t1IsZz786C99;fA-c%Tz^8&a= z_saqM9+hf?u{NyD<XXF2y;r z0&AlE@t}ih@L88Huwd)%#+0$x$L`al(zn#(Z-QOT#&weaJP-C5W1#ZAp7Tl7;GOt zf4b$8I;tqVncFd~83mG$%%LLgPCD$l6a=^!Dv<_Qyk=oAn;%p5E z1T#s!^Ot8wXIGf=hO#9f&SAOyF7q@_RV>hoDdR*jIcdxX^@zx(S<{L|=Qnm>SE$C^ zhb=65$^pWR~AbQhhJyn&r4W0?&{(p)YA?i#? zUN1=v^n&)>)a&O5p>6tIQeMlL3d4@G-j-5I=77%?8XtrUq+trLTwoa*!xZz+cE@`T2z@>{ejYaZy0$3ka22Js!V{yD$Bz!3>tegAzk(;R6ECuOZ(` z+?i{naS_kq5gfUs^$e>l9d%8A!wC-Bg?T$Jih)Nl*Px>v=kEI^B8R2s;&Lnw-ZluQcr z4Gf9kYZ*QCv_Kfy4#{y)vt>b6f1gwG9SP{YRAT<&(j3Nb9kI?cMMpyG?v@oE%Q#=` z*;IZHIb=4(RGscfhi4nES=S#4fcFWR@{>^;_{YAZ7u2^LJLcoV-66t({L?BvC~-b1efJS1+G^^_(#ysyYrh7yPImM)yCu0(DJ|fpAwTlMzb$umg ze##@G*A@SMxc|dsT#QcTJo+F5b6y?xl7)qE;^Cb;)^Gz#2(OQ}3KYqh8q4e?LH>tZ4np5% zvGVX8VpH>KAn`(WD7AwZ|T1dcU7qdFiJx8lYq2)3Svkk*~y{voih4ex>>XLz9>lS+!q}&)VGzZ zB@8<*0i1?yF%*3+Xvuhh0#!G#nLlqVm)+G-L&D&f`1W~xA^O1XCOZujmD56fbwdAB?mLS$3CQk?t@Zr#9$fc1>88OXFNm&oPVCl|LVs?Y2p!#WUFYM0PsBbhSU2dn zf3HItrN24G|9Ivv79A{El@r2%n6G6dtmxBW*xmMOi{Ci@YkQ>Z-h36bG8i%#e@g%| zn$K;#80p1&T9udvv29@FIhVXN+f~}Y$^Xnr? z{B&6JB3P*4b67-)Xkqj4Jbus2YmXnx8a5ahSlwa90`a&)%*WH+c)^HRD_Nc!I?CtP z4LnHD)qa@m2;UqYn8V| z#eQP`@?He|j@|k1{TEhrr~Xh<^lvUO7^7H85MU;<@9;n%7(A>?&h z@P67YMUcMh?T4gau24H7WwPge?TcztPYs zwA+$ZLK5;S_vUxoG-9Alzu{`71k(HeLun>NDE^dlVysUHtyEqv#{yFnfFYNGHnW?S5ch0^gj!DayRHMM3CU^-KJ|~ zBRT}I>Yh^`RE5*(dnHpPNXT)#>^!qK;Iw%Il{gk=6!a%O>!{LSOfa0QRNa6EyM_Nf z72di-#|CB|lYJ$D7{73qzkveXBeuuG(PCH3> zR`%qgetA7!<2+jr^h@92oyDU&%sN$Cn#08YwHg_oC^5liq!y3TeXVDfq7&7w*Ii=^xkAwGWO?&59~S|$dz z&(0X1-#&shJ&_`pG3ek%R~E=|+`t)koA?>m72xEBsdKe$|L~yrA2*IXU&8`3G_`~} z#33c>k1J2CG<^9m;cMQ{3C5_Y@y5{we0XM$%wq~a%4kw=Fw&F-gE$RAt2=KoeRWa) zKeDnQs(W}Klp+8gYzA{x-d%VHhnh>u6Ef00V~uKcI8eBlEce4vG9urc|3c&ufy#%E zWIIlaffa?-_xe#fJk0Ci{o0@i+U4eybnyvnHMh|EaG(gpe0LbpwH5<@zN|#0y8@`x zK2{=2l^dDteeCpQn>0At&*TSbOG1#Pe+yp1gKj9rviBDY!~D;s>D*H^IQFR|W3qem z?)&xc!4GIaVhOldGS`b8mJ`u*I7Wrk+b$VJ!$cTEM$gnSDb%7n_=~_y2FtOmu)sz zvd~cnnN>RosNY>AB-Wl8&2z$hz%eQc*16UgUo8fg4m_TSnk6E$c;>oYf>en86r~jF zKm(0vB1t^69rsw=p4T`kg6sv3znf|%qdnS3$@S@U)RB6DIrF*{i)Jn1rQFNMA_M;zI3B*Dlt zo610!H#Wp;LlH>aF9w}f6XCBsL?9H6^FQcSL0LA1q1v<&Jm0@z;@O!gY|hx?)BLmo zR2pazI~+)RY7W1LvfWu6zy`aO$w{*#C`WKHclO6L%SDboSXDmkaCc8<=3w?_Z+lEh${Z*0=U02=yz$Pm$x}SK|ahj#iNMzKMd8-F#`~D;i*Hu~ax`K!Zc|9WN;q zc@VoI5L#$Fgx%MdGrhW93hiW3>6`aahIegkYoptx;r(`&YwIWE;5Ku6;b1&J5ZA0Y zI-AO|Gx-CBFUGadbNiHpt|D<%Gv0exIYSne@r~-&f0@Iws|GKqwMoFQNne%TzYMrE zqjkK_pAJv2yqnolD1#aw$*R4^#6hRyi%D_cGF~1&{ip4L5IPmJhI3X(!H`C3s>N#w zFb^BoGnJ76mzK1X+&d)DwXLDD^PSwVHN*Evd=wLEdwXRnHBlVp1;lvjL}^2|{M(DQ z*}^E`&2b8M-x`i*?dy1MPllbdsQG@UGD^sQFjd_m22n@P%stW*KrJ~vd(ZdMVdR5T zbDuj0{4>x$F>WFZulA{F76wYgxA!sI^i_CK*aMaTA3asDzIk+o*hT~VfL4qmE)M3> zBcAkpQM53^?vi^?8M&XRp}R|nL0f(9Khni1TxbKCd1(uQPmW}YLXH%$t;Fkf@r+}~ zWsLLP>(tRlyK7l;K2)U6(|Tbqr#iZB%Hfxs$pzDdD@Uqd5rLIiZU>`T5#B!J>i2jf z3g=Q9Yf4t7ki(Jw*h8`c@XNxZsDFtUB6malffx$n{88BP$bk%m=mQe@2~u$Mc0`5u zm^5%7vG3&ZprWgt=^`XOao{q4B1Vsz#J3r?yf_gpg;Zyh?QDDGU@a^&n>C#ug*lW3 zvos6AL3d_Hf*v1ODEk{a$Y>&)%y~WvCIsy9TE#Q>c~MkVQr^;_0{GB}q#RU;K(XxQ zxUxkUK1oFUcv~kAT@NDk5L-LW>b#p=|C0ga-+DWY6t;ootde@mpG7PuR9gs+mnxB{hwz#_g|{n5!u9O7cuF679=>Tlfq z7C{15w^t6D7qWrD#Hi&`FgG}@9Di1Lhzrb_5B5@;g%Q&u&xP0L1))C3#;PTQ74B># z)aA`@U@wzae7uSXV4!@3R?xM97cgG98tSmX2SeinuhJ-BoL|AbU^#{TCt9{|rB)1R zhSlhzMgUgxUxWY2`)l}pHs{vTgCy7*F1TYyQyu1~JH%x{;R5SB0p7ShH%Kjlqmc<6 z(4WSa&Vzl}L7tExzw8NYsXIga^d4S->{eqBGB^D4KXkTbKLO6z#@&>BLxLC*BSWcs z(`OY92XAcd&lr8b<6J*KtgPnBE?YA|VvZp^;YfzWI_>?L&s z76s9pni14;I=Y>x;w-Szji;SD)-}wk1!B7;bZ!0kUvok}6#Krvpor zUe9};L4#dpEu90+WVm=mY_2U|82P=9PB`7i0P(N;16?K;aK+9Pf8UqNP$Bu>HjI4% zd);xqiF9}xzr48Bv(cIcKJi5wtS?B&`I~&cUHw12@XFb7zbY~`?cuxXpr`;LlHG(w zUuJlR16Neg1=MqMe}7F7n2)=R_!S(fp${8!z)I#r=>ZP;Rt+v1JUPWm~Rl^ z^W}RA(H^4c&nOg^}T67SY!&f|~w2M~K! z>2Si>_-obB5ME+J3f>c_1Y50re%Ui^;DM1(yie73<5)Ha-8D%Dt=(fvj$I(YQq6}t zvdb!_IrEwQ8YTu?ey5w=HS_zfAw^X=8&O zX6I$wBPt+TSiYV!EC{=Kxvt215m8XB_c;zS z74(8s?UEYEsKMZ!nRc@jEM0o_*XR)!85Ay?Nl8$eIIXZ{o*Q+yxVi8DIf8xDQDb|JcjFBzWlIOW#o%VMeY6;L7LPAP zsx~8JXykioy{g6ytzO0@9?p?=Srz4A9__ z%f{6Rff+nzS6G-|Jq371q~6>)DF^5JNxPo(zQIq-ZGv6O(j>(Ff<*Om^!fJe+Cg{qV2M@R2QAQB8OSLr#K*RUbJ)DQ)?jYBdQKTl4fA zeA&@StF_|$*(@m6)mW98#)d{y%fG8M$b*E}i$T3GF)(0iHh=wC4!Iulc+A1c3y$wA zgu|+dNbOtFJF^ptV8XW;O(U)24t2p!%<;m|#dRV6%V8NH`B??T$Sz`$Ll>T0H<5#K zlll154GOST{~#Z#7DaGkHr-B42o`S@Gv0_wK%j}@{m*-7;9PQ4s)C&t;-V@|XDd}; z)X1>)u^ky)B{I+U-zmp4pv#2cS`}@%SAAOZy9gY(DOJ6&E`xrmDNGwXbYPWxgCf)1 zx4`o_CGcEd$A>*i(y|WYw#*xRcHv8%PM~N zplS(Iblt|HYC4B&<}2EmWy?X=O)K+@=7ZRQ+d5@wR#Nb8K{>mPCevT(Cpckkx`E_8oEr`5l53Fo>bP25^Qft(MH{J$c5akiFKZLJ3~kY_W; zq;_TV&lFHVzrknA z&0W!qP_>TuJB5eynwMOmYJ=t~7g_&43Ct6P)mLA#p@p&f*33P0;LoRYed{2iKS=Z| z`XvY-z8ky<+9LpgA$?!V6I5XU=J#j*!FC!F@m)-)*+YelPe)uAwYd@A+Cr$k_TqWe!M_#e(wr#>&+Lqky^Oy^fp*kF-e zihjnK8N~-?hTl1<3jZBB7IH_K-JLd+D@w>VAWH0r>X?iT?XCDp1EyUkn3uZ?SBM+)> z$W~*W)O1>n2PcwsOV_An=Yl|Ar>@E`>v(eRqv0gOea!uEcS0x+8NTbr zM$(I?@E;ropKY~k@WF&KyGcGy(33HCPouZvofj};+lpm8oxrY<`C1Y_lu-Lba#&E} z;^PzPmBO$;*iXVtrW@B8`Ve!-oE5l3amTVNB9I(3<>z8Ij+L$L=QYgYhR0!&-S$uU z0f}g{G;i|o9`<^_%L;%W^|@fuZ~kr_iXi|Sm#VIpBLS}I4Wx<~3&2?l)$y~I z1Ux=>Kk?6*AGp5z1)mff1-N*4p~*I!0$F$0t2iae&@xa${;IQrZA`vBLLCqU2_~)B z@XuQ4Ki4R0o(u})5{&LY4yL2tF_T-hK3tHYaj1nV+l}RklzF9G;{?99!r~=SYH-O- z&bE2eANAOa#+Xt?xW>qdj#cO>6qNf(@igrkdWuS&v`|C^5FKZCDtNy4kOOV`PyPi z5JTOTVJpB22F7dq=MPDO3+-x1lrtZcOjWEWzhgz|4u(tF8lqq`D|B^NCkYxjlFgXM zbn>%-Sin^H10JiO=Ez?30z*9D)7N#o<=evvT?OIBgV{SqzgW?Wn9BQKJNY46F_k*hqm9o0 zRIs`{L;=1-<6Omu=;%faC2q`x0G>0RR;8PIyRL{#XA7}F>c!^6*fABj^qx0bV1OO% z%H(C1Nmc@XvEt5~ssvQ@SBYHjK8f$tev@@|ObIa`{H-{kvB`hHpx*^q5GIt*4QKO$ zRrrnq(Va4&Bp)%z9@31}6uxs7Jt&C0o{#j4wJ{*T;8KuR`A@8qeq)PdiZbxu)(RCj zWrI%jH-HvCxkUfsbUcwWVsN?@|PXOyAV)ags0;7az0OAOtV&Oi2945QWU%s+o7;{g~0O z%!ZrU3Q!m3AJ93@1yY2xc-}{{5V`1g+EQf$vowBKs_nhWyS1cB&hlwB3JDUl?sq?tL+YrxbFFJG{R>)>G(t@ThvuUxuoBb7QR{Jee z22|Q+WIsG2f{LMbrji~TaGNd5R~IXxupzG#5-~DRarWhj&fQxTB%?~MG4K%sg1Oqme<8x~`vX@Q!B!ErykOCJIrtfG zf9K`JJkANcS;pmlfs$}q{mPCLmnFg3da&bdr2uGucXo>i##wuiZm}#zNPwcxeu_kXUlc=UXivunKDH z=B#KaIyOHqP8>eAJ=JK@UBlhytmt1;7$9K$>qOcvEjWB)=MOarBBV+xR{7W}fI{lN zt3`%Fur?z{wX@y8vb%T+G!rz@hiaqq38!RX%ZjZs(L)&izVs=^9a+IsefQlXj(k81 z+2^zPC0md5sdN<+n_& zP@WJR=a!c*w3xz^wq4KE1UeM$KROvwpa=^SVthEycWfd0PjV=)B&^kl9Mjo8hc*7Z zsEHYiqk|7T%1{21LnpmvF4oI)z;Ce?-s(gmJhS1Zn#RxLrH8h6Zwzw4UDFxK-vu++ z(I(~W`)O}zcvn}9O{t*idR5mWq8(p2JFmmussN*3IOHR?QXoYr z?3tQ25%TV2=6bGjA!9@Nv&kmBkjQ3qlg~y6&23~;PbiClb|~hRR7yvKoVeW+4HkHJ z+~_v7oehPH?|A>yfDN3(^k*yr6yd+iPsQ&Z@U8kl#h6{sFlDu;66$JHwa<)8}$1t_Df7wZa<9 z3ee%oqno2Vgqyv3*kol!LJdV8Gqozh$exnC+s&R1`d?4p*-|%&>9{T5Tjyton-2M< z!8KI$@Y>5*}e%8L7rQ33@Y3vl;nGPoDgq=^vDG<;*?I~y`3fkkpo-FTRhQy$Xe*J-e_#5@3TJJY`qF!a` zVACH8h-SI_n5lC>M5h+QO5elKiu_?RRFv@UF0_8!7zg2>W6Ivn6PHWx@gwTk=9!Pe;q0(=jWZka zdT_j;Tr)(}UI)D>t8a>Y%Zg@6m4#7? zTc28;BA6yTQ`jn`4851Go_FWn{JR%3zZv=nLLSA5NS0s3-tIe_Hk_=DZ+ zih1w!R$AXF2_&J%Pdj!oV8&DAMpm5w)W`*!UDzc8gf#Pd@5AfZ{{R30|Njh^hd?y8Mj6?eQTARx zf5P*e*E#2XHrU@T5+v1QCb_#l=6Eo`R0ajx39AAlHshYP2N}L*Fq69Mmob3`nvs*f zJXERrjn( zN$Dxf9x8#UaYJtFZ#9tjA+ooelYt?^iH(;Nn=y05%vx(xZ0@>A{TSpcMkdo&6t#ktg{#zc?VPdP8 zyz3qX;OoO`@a8jV3QTzMNVP-%tQvg%EybecP{Gekk?UwNH?l8y-vp^N7(4pf)3#L( z6c_cK^b$A`G5QM$KU|8>%V(tB>ms5Op%smPj%;Xrdm3}`9>6y2AU>Cg0S!;x|F-H#RXc=nKf zaKvjS{M9=4;r$zJ;P1ZD9qT~@+Ge`-^k+G^mv&xSX)_5-hx&_mB@be;AD#zs^QobP zH@jcOepH2yZ2J406$a$`xTjQo{f$f9^bQa5WP|i3$%60a25{=XH!0(FQpkFLWYy=B z8X#N4l8rmc30_BwTf$lxu=F+hOuZ5rv=5QSF02=W)YC3^bplkNTI>6Xz;7gUa7wjd zT2CC#a*WRe9hHHhJ-8X8P7z#mMfRQ9&I>p0IA;-BdBN=aV(o)yB4Cx_ryhjJ10hlU zMvFNIxSqqW_$7>DI@W!>8dHF+QVqYf#WG-71c6)N+XljyXW1k3D{Si?JoC9piiH^Jc*M1i^&&IjES*=NG95`qcTMnBpJd$Fy)1r`!SpPellTD zr!lv6SqJnjsWrCQQn0|(_Jo%*VEEnb#P47c@J}>!7p;{+k0KT_>aS^|A8uq@3rQQ9Q#h1Xn|YIrb#a;W8a!XJTbkpmK~Zrp?8P_KxX_ zX)P>bSAUxowb@I<>Xgi{#hrBc-Ya)W!)yqTGkdtxwnI&9!E3;k~E{j^ts6P#MemJ$E)L#gk9$aD~a!9Q&94`G2$ zu=0SV?=)6a{`UUY^fA2eMDtuyU^6CJs^$Lam?}D`b<0lJjEVAJ!KH=2M9^|pN{$he z2SLh(4=4YVfbH##3ql8Gak4*YIr-iUZpolca?47h|Ehm?+pR9(cG_*CX@w-n*n1hi zNQyxWjW7j=pm(vh=G%3#VvA^57@EByPQAR66`GLG4)z-l~5{dN!rqoqJz zcV}7nbkr_=w+nkMcW6W#Cm~=*m=j1HGO(2}q+7UE2Z)>fOLEnvkk(nT@%1~%$XC4W zgJ||RF06Q{`Mx?Ec84tyuRmVE)WP$C`&J=z>0AWY&Z$*w&&8tG!gw8)D0{dl#oQPP zS_!PG@(7@96V*xFU*wUignuXXlsXjdSa%24SjJ}4S5{eZYzXQ;sN)}?3}r~;fSYo-O83E3my2*3w9tF4`SU3%{xKYQyZo6OCCNRUBd>Y|@v!kYp#sXUZvHd?SxNF%DQ@9dp83FJn`HAUobCKz|UdzVBff^^;x;B-Ku8AVZ(@DAKd*{%{SO5-^02^3D{@P^qB>|u{!wlGOMuzCh5Dqo z@}Ml#-L`3;1adfGvh3ucinc%X5Lw-z13D_w^LsP;@Dl&ATAd42xU>IB9Ud+Z%}Ta% z>)e;|c;T&j7lv7&xZmE{#Y75)(QGbeP07Kq^;ylc!8G`6s_1%Xvl7(ftasvf72$@r zyG3Cc7n1i{apu@S1^4O-@1bBe_-T4vX-T9aUbp356@F<{Fv{vFa2&z6UJRNUnIeN< zRC8t$E{^8)ZmJXv&taz@m}REfGvHu0OKGw|9_lu(=TaVGq7nSOz?q{&bm;Fj;n;Nw zNRpC!7FQ!deYe|KbnYCc#%2s1&=H5-CELdLOUc3SG<4?*ML^taI}C z=qN5XF07{M#zxTPZOQ1-K!t~HY;wp}f(~yA_Q6aH4llm{HASF-%}&nLo-je^b5ECvK;sw( zqs=)wiw=Ba7iKFKD5#_7#Pfz&6<9i%C^Vl(fdcx=D)BxB7}gDXiFB$!?i0@W4fogj z?F6x7D2Wbs;ZGKa4r>F+qlj$zT@vZV6gGcnR^TJr#Yqk8MN!$h{h9VI!ce)|q!&K> z2N%_BAK7$61XB5}Hf9ZUW1>Hvy-zOa!Gz&R;kvvb{CJ1!L^PKm>|m@2{9+oY@z4H>D`pr$vP@f3Ia}80NTWR2L#@1+RlmKa$$aKK~0*nyXcP#qU7qmP_)@xNv0EQMbo2PQ@T}!R zQElYL%-7r~;<}!K|lxCF4yC_y z6awFy7C+_;2*67jP)?H=!VA-y3I$@tkxHI$jTARO+Mwn-$wipN51uA?zuPGV#h;oA zWUF|AvMJ;UX}dhS-BE3Hzn26%IfK3l+TF#zR+b+mS<@hxsukIoYd(A~e5G zo`z1vEAe+HXu_rw0hec7h``;gZ%zs#!xEaXH(TSqZ`|1yBHCS8xV)3_grf*ps_gXc005 z!RC1SK3O@LwR{BqRbKjPGEGab>ds)nClE?S|qNMza$8y!YjpYWfRyy zam>f1|A>eyI~oHFlxE(tSy9>=S#=-@D1YBW(S16N)x?9cwB4~#aVi9@$M z_$nRKGaV!&`J8h;k-zxShzYlA=wmWkdJ+_U^ocZ#S8vjImCFsiv&%!t72+^WvX