diff --git a/.gitignore b/.gitignore index 39e49b8a2..5882fd01e 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,6 @@ +**/.DS_Store docs/generated/ +docs/_generated_images/ .pytest_cache/ # Byte-compiled / optimized / DLL files diff --git a/docs/_static/physics_kundu_2017_TE_dependence.jpg b/docs/_static/physics_kundu_2017_TE_dependence.jpg deleted file mode 100644 index 265bbbcc0..000000000 Binary files a/docs/_static/physics_kundu_2017_TE_dependence.jpg and /dev/null differ diff --git a/docs/_static/physics_kundu_2017_multiple_echoes.jpg b/docs/_static/physics_kundu_2017_multiple_echoes.jpg deleted file mode 100644 index 6a1fe3a95..000000000 Binary files a/docs/_static/physics_kundu_2017_multiple_echoes.jpg and /dev/null differ diff --git a/docs/index.rst b/docs/index.rst index 695bc289b..ab078a61f 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -160,6 +160,7 @@ tedana is licensed under GNU Lesser General Public License version 2.1. installation multi-echo + multi-echo_walkthrough usage approach outputs diff --git a/docs/multi-echo_walkthrough.rst b/docs/multi-echo_walkthrough.rst new file mode 100644 index 000000000..ba266b15d --- /dev/null +++ b/docs/multi-echo_walkthrough.rst @@ -0,0 +1,314 @@ +######################## +What is multi-echo fMRI? +######################## + +.. admonition:: TL;DR + + Most echo-planar image (EPI) sequences collect a single brain image following + a radio frequency (RF) pulse, at a rate known as the repetition time (TR). + This typical approach is known as single-echo fMRI. + + In contrast, multi-echo (ME) fMRI refers to collecting data at multiple echo times, + resulting in multiple volumes with varying levels of contrast acquired per RF pulse. + Multi-echo fMRI can be used to identify and remove certain types of noise + that can't be removed from single-echo data. + +To understand what multi-echo fMRI is and why it's useful, +we will walk through a number of things. + + +******************* +The physics of fMRI +******************* + +In a typical fMRI sequence, +the protons in a brain slice are aligned perpendicularly to the main magnetic +field of the scanner (:math:`B_0`) with an "excitation pulse". +Those protons start releasing energy from the excitation pulse as they fall back in line with :math:`B_0`, +and that energy is actively measured at an "echo time" (also known as `TE`_). + +The time it takes for 37% of the protons fall back in line with :math:`B_0` +is known as the transverse relaxation time, or "T2". +The exact value of T2 varies across voxels in the brain, based on tissue type, +which is why protocols which are T2-weighted or which quantitatively measure T2 are useful for structural analyses. + +fMRI works based on the fact that differences in blood oxygenation +(indicative of gross differences in neural metabolism ostensibly driven by neural activity) +also impact observed T2, also known as T2*. +Thus, changes in the measured signal from an fMRI sequence reflect (at least in part) changes in blood oxygenation. + + +***************************** +BOLD signal and BOLD contrast +***************************** + +As mentioned above, energy is released after the excitation pulse, and this is our observed fMRI signal. +The echo time is the point at which that signal is recorded. + +Let's take a look at how fMRI signal varies as a function of echo time. + +.. image:: https://mfr.osf.io/render?url=https://osf.io/m7aw3/?direct%26mode=render%26action=download%26mode=render + :alt: physics_signal_decay.png + +As you can see, signal decays as echo time increases. +However, one very important feature of this signal decay is that the signal decays differently +depending on the level of blood oxygenation in the voxel. + +When a voxel contains more deoxygenated blood, +its signal decays more slowly than when the blood within it is more oxygenated. + +.. image:: https://mfr.osf.io/render?url=https://osf.io/ve7cf/?direct%26mode=render%26action=download%26mode=render + :alt: physics_signal_decay_activity.png + +This is the "BOLD contrast" functional neuroimagers care about. +Namely, you can compare the signal from a voxel during a cognitive task +(when you expect that voxel's brain region to be more "active") +against that voxel's signal from a different cognitive task in order to determine if there is a +meaningful different in activation between the two tasks. + +There are two relevant features to the BOLD contrast. +First, as repeatedly noted above, overall signal decays as echo time increases. +Eventually, with a long enough echo time, there is not enough signal to detect over noise +(i.e., signal-to-noise ratio is too low) +and meaningful signal cannot be observed. +Second, BOLD contrast increases as echo time increases. +That is, the relative difference in signal between active and inactive states increases with echo time, +even as the overall signal decreases. + +While the signal decay curve can be described using many models, +one of the most useful approximations (and the one we used to simulate the signals in this walkthrough) +is a monoexponential decay curve. +In this model, the signal is driven by two factors: +the "intercept" (signal at echo time = 0s), also known as S0; and the "slope" (decay rate), also known as T2*. +For more information about different models of fMRI signal decay, please see XXXXX. + +.. note:: T2* + We said above that observed T2 is T2*, and now we're saying that the slope of the decay is T2*. + That's because they're the same thing! + +In the above figure, the difference between the "inactive" and "active" signals was driven by a change in T2* from 20ms to 40ms. +Thus, the point at which the difference in signal is maximized between the two states in this voxel is ~30ms. + +And, in fact, this is the standard approach for the most common version of fMRI: single-echo fMRI. + + +************************* +What is single-echo fMRI? +************************* + +Above, we saw that the difference in BOLD signal between active and inactive voxels is maximized at a specific echo time. + +With single-echo fMRI, we have one echo time for each excitation pulse. +As such, we record one data point for each voxel, at each time point. +That data point is assumed to reflect blood oxygenation-level dependent signal, +and the single echo time is chosen to maximize BOLD contrast across the brain. + +However, a major drawback to this approach is that it ignores several features of fMRI data: + +1. T2* varies across voxels in the brain, depending on features like tissue type distribution, brain-air boundaries, etc. + As such, the "optimal" echo time _also_ varies across the brain. +2. Fluctuations in voxel activation (i.e., T2*) mean that T2* (and thus optimal echo time) + varies _within_ each voxel as well as between voxels. + Basically, "activation" is not a binary state: a voxel can be active at many different levels, and each level reflects a change in T2*. +3. With only a single data point at each time point, there is no way to characterize the signal decay curve. + Unfortunately, there are multiple factors which can influence how signal decays, and many of them are problematic. + With single-echo fMRI, we must rely on more standard signal- and image-processing techniques to remove noise. + +Taking these drawbacks into account, one might wonder why anyone would use single-echo fMRI +when there _must_ be something called multi-echo fMRI (what with that being the topic of this walkthrough). +There are a bunch of reasons, but perhaps the biggest is that fMRI is always balancing the amount of data against its utility. +You can always get more useful information from data that is closer to its original form (e.g., in k-space, before combining across coils), +but this involves increased work for the researcher, in that they must choose the appropriate tools and process all of those data, +as well as _massively_ increased storage requirements. + +Single-echo fMRI performs quite well considering its limitations. +While T2* varies across the brain, for _most_ regions the typical echo time for a given magnetic strength +(e.g., 30ms for 3T) will result in sufficient signal-to-noise ratio (SNR). +BOLD contrast can reliably be detected over noise, and in most analyses the SNR is sufficient. +With a standard univariate analysis, more data will generally swamp issues like thermal noise, +especially when those issues are not correlated with measures of interest (like cognitive tasks). + + +========================================== +Sources of fMRI signal fluctuations +========================================== + +There are, in fact, many factors that impact observed fMRI signal, but we will focus on a small number. + +1. Neurally-driven BOLD signal. + This is the signal we generally care about in fMRI. +2. Non-BOLD noise. + This noise is often driven by things like instrument noise, subject motion, and thermal noise. +3. Non-neural BOLD signal. + There are physiological sources of changes in blood oxygenation that are unrelated to neural activity, + including heart rate and breathing changes. + +Let's take a look at what single-echo data looks like over time. + +.. image:: https://mfr.osf.io/render?url=https://osf.io/g9dqc/?direct%26mode=render%26action=download%26mode=render + :alt: fluctuations_single-echo.gif + +As you can see, the single data point fluctuates over time. +Let's assume that those fluctuations reflect meaningful BOLD signal. +Nothing to be concerned about, right? + +Okay, let's check out the underlying signal decay curve we're sampling from. + +.. image:: https://mfr.osf.io/render?url=https://osf.io/5yjwx/?direct%26mode=render%26action=download%26mode=render + :alt: fluctuations_single-echo_with_curve.gif + +Everything still looks fine, right? +We know there's an underlying signal decay curve, and we're sampling that curve at a single point, at our TE. + +What if we describe the curve in terms of S0 and T2*? + +.. image:: https://mfr.osf.io/render?url=https://osf.io/6a7nv/?direct%26mode=render%26action=download%26mode=render + :alt: fluctuations_single-echo_with_curve_and_t2s_s0.gif + +Now we see that the changes in the signal are driven by changes in _both_ S0 and T2*. +Why should we care about that? +Well, we know that T2* reflects BOLD signal, but we don't really care about S0. +In fact, S0 changes are driven by non-BOLD noise. +Things like motion, thermal noise, instrument noise, etc. + +So if our observed signal is affected by both S0 and T2*, +and the S0 changes are introducing noise into our data, +is there anything we can do? + +Well, first, let's see if there's a way to tell S0-based fluctuations and T2*-based fluctuations apart. +We'll plot two signal decay curves. +One will _only_ include S0 changes and the other will only include T2* changes. + +To make sure we can _really_ see the curves, we'll also make the S0 and T2* changes roughly equivalent. +They have different scales, so we'll use the same time series of fluctuations, +scaled to have matching percent signal changes between the two values. + +.. image:: https://mfr.osf.io/render?url=https://osf.io/g29ez/?direct%26mode=render%26action=download%26mode=render + :alt: fluctuations_t2s_s0.gif + +Hey, look at that! +The curves change differently! +If you look at the whole curve, you can differentiate S0 changes from T2* changes. + +Now that we know that, what about single-echo fMRI? + +.. image:: https://mfr.osf.io/render?url=https://osf.io/mx4ku/?direct%26mode=render%26action=download%26mode=render + :alt: fluctuations_t2s_s0_single-echo.gif + +Hm... with only one data point per time point, we really can't tell whether the changes are due to S0 or T2*. + +What if... what if we had _multiple_ data points for each volume? + +*************** +Multi-echo fMRI +*************** + +Multi-echo fMRI involves defining and acquiring multiple echo times in your sequence. +Instead of sampling the decay curve at one point, you sample at multiple points. + +Typical multi-echo protocols use somewhere between three and five echoes, +though more are possible if you make certain compromises with your parameters. + +Here we have some simulated data with six echoes. + +.. image:: https://mfr.osf.io/render?url=https://osf.io/mf3ae/?direct%26mode=render%26action=download%26mode=render + :alt: fluctuations_t2s_s0_multi-echo.gif + +Now we can tell the two curves apart again! + +Okay, so what does this all mean? +Simply put, you need multiple echoes in order to differentiate S0 and T2* fluctuations in your fMRI data. + + +.. _multi-echo physics2: + +****************************** +The physics of multi-echo fMRI +****************************** + +Multi-echo fMRI data is obtained by acquiring multiple echo times (commonly called +`TEs`_) for each MRI volume during data collection. +While fMRI signal contains important neural information (termed the blood +oxygen-level dependent, or `BOLD signal`_, +it also contains "noise" (termed non-BOLD signal) caused by things like +participant motion and changes in breathing. +Because the BOLD signal is known to decay at a set rate, collecting multiple +echos allows us to assess non-BOLD. + +The image below shows the basic relationship between echo times and the image acquired at +3T (top, A) and 7T (bottom, B). Note that the earliest echo time is the brightest, as the +signal has only had a limited amount of time to decay. +In addition, the latter echo times show areas in which is the signal has decayed completely ('drop out') +due to inhomogeneity in the magnetic field. +By using the information across multiple echos these images can be combined in +an optimal manner to take advantage of the signal +in the earlier echos (see :ref:`optimal combination`). + +.. image:: https://mfr.osf.io/render?url=https://osf.io/m7aw3/?direct%26mode=render%26action=download%26mode=render + :alt: physics_signal_decay.png + +.. image:: https://mfr.osf.io/render?url=https://osf.io/m7aw3/?direct%26mode=render%26action=download%26mode=render + :alt: physics_multiple_echos.png + +In order to classify the relationship between the signal and the echo time we can consider a +single voxel at two timepoints (x and y) and the measured signal measured at three different echo times - :math:`S(TE_n)`. + +For the left column, we are observing a change that we term :math:`{\Delta}{S_0}` - that is a change +in the intercept or raw signal intensity. +A common example of this is participant movement, in which the voxel (which is at a static +location within the scanner) now contains different tissue or even an area outside of the brain. + +As we have collected three separate echos, we can compare the change in signal at each echo time, :math:`{\Delta}{S(TE_n)}`. +For :math:`{\Delta}{S_0}` we see that this produces a decaying curve. +If we compare this to the original signal, as in :math:`\frac{{\Delta}{S(TE_n)}}{S(TE_n)}` +we see that there is no echo time dependence, as the final plot is a flat line. + +In the right column, we consider changes that are related to brain activity. +For example, imagine that the two brain states here (x and y) are a baseline and task activated state respectively. +This effect is a change in in :math:`{\Delta}{R_2^*}` which is equivalent +to the inverse of :math:`{T_2^*}`. +We typically observe this change in signal amplitude occurring over volumes with +the hemodynamic response, while here we are examining the change in signal over echo times. +Again we can plot the difference in the signal between these two states as a function of echo time, +finding that the signal rises and falls. +If we compare this curve to the original signal we find +that the magnitude of the changes is dependent on the echo time. + +For a more comprehensive review of these topics and others, see `Kundu et al. (2017)`_. + +.. _TE: http://mriquestions.com/tr-and-te.html +.. _TEs: http://mriquestions.com/tr-and-te.html +.. _BOLD signal: http://www.fil.ion.ucl.ac.uk/spm/course/slides10-zurich/Kerstin_BOLD.pdf +.. _Kundu et al. (2017): https://www.sciencedirect.com/science/article/pii/S1053811917302410?via%3Dihub + + +******************* +Why use multi-echo? +******************* + +There are many potential reasons an investigator would be interested in using multi-echo EPI (ME-EPI). +Among these are the different levels of analysis ME-EPI enables. +Specifically, by collecting multi-echo data, researchers are able to: + +**Compare results across different echos**: currently, field standards are largely set using single-echo EPI. +Because multi-echo is composed of multiple single-echo time series, each of these can be analyzed separately +and compared to one another. + +**Combine the results by weighted averaging**: Rather than analyzing single-echo time series separately, +we can combine them into an "optimally combined time series". +For more information on this combination, see :ref:`optimal combination`. +Optimally combined data exhibits higher SNR and improves statistical power of analyses in regions +traditionally affected by drop-out. + +**Denoise the data based on information contained in the echos**: Collecting multi-echo data allows +access to unique denoising methods. +ICA-based denoising methods like ICA-AROMA (`Pruim et al. (2015)`_) +have been shown to significantly improve the quality of cleaned signal. +These methods, however, have comparably limited information, as they are designed to work with single-echo EPI. + +``tedana`` is an ICA-based denoising pipeline built especially for +multi-echo data. Collecting multi-echo EPI allows us to leverage all of the information available for single-echo datasets, +as well as additional information only available when looking at signal decay across multiple TEs. +We can use this information to denoise the optimally combined time series. + +.. _Pruim et al. (2015): https://www.sciencedirect.com/science/article/pii/S1053811915001822 diff --git a/examples/plot_approach_figures.ipynb b/examples/plot_approach_figures.ipynb index 29f293668..80ba7db47 100644 --- a/examples/plot_approach_figures.ipynb +++ b/examples/plot_approach_figures.ipynb @@ -982,9 +982,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python [conda env:python3]", + "display_name": "Python 3", "language": "python", - "name": "conda-env-python3-py" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -996,7 +996,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.4" + "version": "3.8.5" } }, "nbformat": 4, diff --git a/examples/plot_metric_simulations.ipynb b/examples/plot_metric_simulations.ipynb index 8a9fc50ee..b707cefca 100644 --- a/examples/plot_metric_simulations.ipynb +++ b/examples/plot_metric_simulations.ipynb @@ -672,9 +672,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python [conda env:python3]", + "display_name": "Python 3", "language": "python", - "name": "conda-env-python3-py" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -686,7 +686,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.4" + "version": "3.8.5" } }, "nbformat": 4, diff --git a/examples/plot_signal_decay_maps.ipynb b/examples/plot_signal_decay_maps.ipynb new file mode 100644 index 000000000..cc652811e --- /dev/null +++ b/examples/plot_signal_decay_maps.ipynb @@ -0,0 +1,127 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/taylor/Documents/tsalo/nilearn/nilearn/datasets/__init__.py:87: FutureWarning: Fetchers from the nilearn.datasets module will be updated in version 0.9 to return python strings instead of bytes and Pandas dataframes instead of Numpy arrays.\n", + " warn(\"Fetchers from the nilearn.datasets module will be \"\n" + ] + } + ], + "source": [ + "%matplotlib inline\n", + "import os\n", + "import os.path as op\n", + "from glob import glob\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import nibabel as nib\n", + "import numpy as np\n", + "from nilearn import image, plotting" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "data_dir = \"/Users/taylor/Documents/datasets/logans-dset/sub-pilot/func/\"\n", + "files = sorted(glob(op.join(data_dir, \"sub-pilot_task-checkerboard_run-2_*_bold.nii.gz\")))\n", + "echo_times = np.array([9.58, 21.95, 34.32, 46.69, 59.06, 71.43, 83.8, 96.17])\n", + "\n", + "img = image.index_img(files[0], 0)\n", + "data = img.get_fdata()\n", + "vmax = np.max(data)\n", + "idx = np.vstack(np.where(data > 1000))\n", + "\n", + "min_x, min_y, min_z = np.min(idx, axis=1)\n", + "max_x, max_y, max_z = np.max(idx, axis=1)\n", + "\n", + "imgs = []\n", + "for f in files:\n", + " img = image.index_img(f, 0)\n", + " data = img.get_fdata()\n", + " data = data[min_x:max_x, min_y:max_y, min_z:max_z]\n", + " img = nib.Nifti1Image(data, img.affine, img.header)\n", + " imgs.append(img)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABbkAAAD6CAYAAACMNkrcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAADXcUlEQVR4nOy9d9hdVbW2PyihhwSSQEhCekgICUVApDdFsIDY4KhgwXZUFP15bJ+Kfh712PCAnoNHKYoKCqhwFEUpQXonhEB6b6QQairg+v3hN1bu/b5jZK+9s5O8Ozz3deVist6155przjHnmnvt+TxzKzMrTAghhBBCCCGEEEIIIYRoQ7be3AUQQgghhBBCCCGEEEIIIZpFL7mFEEIIIYQQQgghhBBCtC16yS2EEEIIIYQQQgghhBCibdFLbiGEEEIIIYQQQgghhBBti15yCyGEEEIIIYQQQgghhGhb9JJbCCGEEEIIIYQQQgghRNuil9xCCCGEEEIIIYQQQggh2pZN9pK7KIqG/r33ve81M7PLL7+87rmXX355y8u7yy672L//+7/bpEmTbNWqVbZ8+XK78cYb7YQTTmgon2OPPXa9Zf/2t78dfm7w4MF28cUX26RJk2zFihX25JNP2t13320f+tCHrFu3bq24RdEk7RTLw4cPt8997nN2yy232Ny5c23NmjX25JNP2nXXXWfHHXdc+Jn+/fvbl770Jbv66qtt2rRp9vLLL1tRFDZs2LCmyjBs2DC77LLLbN68ebZmzRpbuHChXXHFFTZ06NDw/HHjxq23jrbffvumyiGq004xPmDAAPuv//ovu/fee23RokW2evVqW7Bggd1+++32vve9z7bddtu6eWy33Xb22GOPWVEUNm/evIauf9ppp9lvfvMbmzRpki1fvtxWrlxpU6dOtSuvvNIOPvjgTufvtNNO9q53vct+/etf26RJk+yFF16w5557zh544AH7zGc+o/F9E9JOcR7xs5/9rLze+sbnPffc0y644AKbPHmyrVy50pYvX24PPfRQOv9YH43m9da3vtXGjRtnzzzzjK1cudImTpxoX/jCFxTnm4h2ivFBgwat93pXXXVV+LlG5xjrY+utt7bzzjvPHn30UVu5cqU99dRTdsMNN9jhhx++3s8NGzbMfvrTn9rMmTNt1apVtnTpUrvnnnvsM5/5TMNlEI3RTjFe5Zo333xzzWdaPScnVeY+n/3sZ+2GG26wWbNm2fPPP2/PPvusTZgwwX7wgx9Y//79N7gMonnaKfbN/hlvH/vYx+y+++6zpUuX2vPPP29PPPGEXXjhhTZw4MBO5x999NF2xRVX2GOPPWbLli2zVatW2cyZM+36669v+H2MmVmfPn3sO9/5jj322GP23HPP2bJly+zBBx+0z372s7bLLru04hZFk7RbLDf7nrAV8/FRo0bZ1772Nbvuuutszpw55X1us8024fn13kX6vwEDBjRUB12d+t/+W8TXvva1TsfOO+8869mzp/3nf/6nPfPMMzV/Gz9+fM3/X3fddZ2OZeduKD179rQ777zT9ttvP5s4caL95Cc/sV122cVOO+00u+WWW+ycc86xyy67rKE8b7vtNrvttts6Hb/zzjs7HTvkkENs3LhxtuOOO9qNN95o119/ve2666725je/2X7605/a2972Njv55JObvT2xgbRTLH/jG9+wM8880x5//HH785//bMuXL7eRI0faqaeeaqeddpp98pOftB/96Ec1nznkkEPsm9/8pv3jH/+wWbNm2bPPPmu77bZbU9c/+OCD7dZbb7Vdd93Vbr75Zrvqqqts0KBBduaZZ9qpp55qxx13XHrPUT2bmb300ktNlUVUp51ifNiwYfbud7/b7rvvPrvuuuts+fLl1qtXLzvllFPs8ssvt7POOstOOukke/nll9M8vvWtb9mgQYOauv5pp51mhx56qD3wwAO2cOFCW7t2rQ0fPtxOP/10O+OMM+zDH/6wXXrppeX5Rx99tP3617+2p556ysaNG2fXXXed7bbbbnbqqafaD37wA3vrW99qJ554oq1Zs6ap8ojqtFOcd+RNb3qTffCDH7Tnn3/eunfvnp53xBFH2J/+9Cfbaaed7M9//rP94Q9/sB133NGGDx9uZ555pn3xi1+sfM1G8/rmN79pX/rSl+z555+33/3ud7Z8+XI7+uij7dvf/radeOKJdsopp2g838i0Y4yPHz/errvuuk7HJ06c2OnYhswxIn7zm9/YO97xDps8ebL9+Mc/tt13393OOOMMu/322+1tb3ub/e///m+nz5x++ul25ZVX2osvvmh/+tOfbNasWdajRw8bOXKkvfWtb7ULLrigkdsXDdJOMX7dddfZ7Nmzw7+dddZZNmzYMPvLX/5Sc7yVc/KOVJn7fOQjH7EXXnjB/v73v9vixYutW7dudtBBB9lnPvMZO+eccxruY6J1tFPsb7PNNnbLLbfYUUcdZZMmTbKrrrrK1qxZY4ceeqh98pOftLPPPtuOOOIImzRpUvmZE044wU444QS777777NZbb7UVK1bYwIED7dRTT7VTTz3VvvGNb9hXv/rVStcfNGiQ3XfffbbnnnvauHHj7C9/+YvtsMMOdtJJJ9n3vvc9e8973mOvec1rbPXq1S29b1GNdorlZt8Ttmo+/vrXv97OP/98e+mll2zatGm2atUq23HHHdPzZ8+enb5TGTt2rL3tbW+zxx57zObPn1/p+u1Esbn+zZo1qyiKohg0aFB6zuWXX14URVG8973v3WTl+s///M+iKIri2muvLbbZZpvyeJ8+fYo5c+YUK1asKPr3718pr2OPPbYoiqI4//zzK1//T3/6U1EURXH22WfXHN9pp52KiRMnFkVRFEcfffRmazf96/yvq8bye9/73uLAAw/sdPyYY44p1qxZU6xevbro27dvzd/69+9fHHXUUUX37t0LMyvGjRtXFEVRDBs2rOHrjx8/viiKojjvvPNqjh955JHFiy++WDzyyCOdPuPX29xtqn+1/7pqjHfr1q3YaqutOh3fdttti1tvvbUoiqJ4xzvekX7+2GOPLV5++eXiIx/5SFEURTFv3ryGrr/99tuHx8eMGVOsWrWqePrpp4tu3bqVxw844IDiXe96V80xMyt22WWX4sEHHyyKoig+85nPbPb2fqX+66pxzn+9e/cuFi1aVFx11VXrHZ/33HPPYunSpcWsWbOKESNGdPr7tttuW/majeZ10EEHFUVRFMuXLy+GDBlS87f//u//LoqiKD796U9v9vZ+Jf7rqjE+aNCgoiiK4vLLL6/8mWbmGNm/M888syiKorjzzjtrxvVDDjmkWL16dbF48eJil112qfnMfvvtV6xatap46KGHij333LNTno30Mf1r3b+uGuPZvx49ehQrVqwoVq9eXfTq1avmb62ck/Nf1blPNsf54Ac/WBRFUdxwww2bvf70b92/rhr7b3/724uiKIqbbrqp05z9a1/7WlEURXHppZfWHM9ir1+/fsWTTz5ZvPTSS52+w2b/fvzjHxdFURRf/epXa45vvfXWxc0331wURVGcddZZm7399G/dv64ay828J2zlfHyfffYpXv3qVxc77LBDTT2xLFX/XXnllUVRFMW555672du71f/kyR1w+umnm5nZV7/61ZrVf0uXLrULLrjAdtppJ/vABz6w0a7vEsuOK0ZWrlxpt9xyi5n9U3JDiqKwcePG2R577GGXXnqpPfnkk/bCCy/YXXfdZUcddZSZ/VMm/93vftdmz55tq1evtokTJ9rb3/72Ttfv1q2bnXvuufbQQw/Z8uXLbcWKFTZr1iy77rrr7MQTT9wYtyw2Er/4xS/CXzBvv/12u+2222z77be3I444ouZvCxYssDvvvNOef/75Dbr2kCFD7IADDrDFixfbhRdeWPO3u+66y/70pz/ZgQceaEcfffQGXcdlOOeff74dfPDB9pe//MWeeeYZW758uV177bWl/GbIkCF21VVX2ZIlS2zlypV266232v77798pvz322MO+973v2eTJk+2FF16wp59+2iZPnmyXX365DRkyZIPKKlrPiy++aP/8TaSWl156qVwNOGLEiPCz3bt3t5///Od2yy232P/8z/80df1sxfXEiRNt0qRJ1rNnz5rx+tFHHy1X/ZEXXnjBfvCDH5iZdbISUowL8tOf/tTMzD7+8Y+v97wvfelL1rt3b/voRz9q06ZN6/T3RlZRN5rXW97yFjMzu+SSS2zWrFmd8orK/973vreUob72ta+122+/3Z5//nlbsmSJXXbZZdajRw8zMzvwwAPtj3/8oy1fvtyef/55u/7668PViEOGDLH/+Z//sWnTppXWExMmTLCLL77Ydt9998r3LromrZ5j/Ou//quZmX35y1+uGdcffPBB++1vf2t77LFHpznzt771Ldtuu+3s3e9+ty1evLhTnh37hVuyXH755TZ06FC75pprbNmyZfbcc8/ZX//6V9tvv/3MzKx37972P//zP7Zw4UJbtWqV3X///aHF3C677GJf/vKX7bHHHrNnn33WnnvuOZs+fbr95je/sVe96lWV7ltsfs466yzbaaed7Pe//7099dRTNX9r1ZycNDL3yeY4V199tZl1nl9pviIi/N3GDTfc0GnOfv3115tZ53cbWewtXLjQ7r77bttmm20q21Jl71b+8Y9/2A033BBef9asWTZr1izbeeed7YILLrC5c+faypUr7ZFHHrHTTjvNzP65Qv1LX/qSTZ061VatWmXTp09P52Znn3223XXXXbZkyRJbtWqVzZ0712688UZ75zvfWekeRNegmfeErZyPT5061e6///4NVh306tXLTj/9dFu5cqVdccUVNX/bEsZxveQO6Nu3r5mZzZw5s9Pf/FijL3uHDx9uH//4x+2LX/yivf/977fhw4en5z7++ONmZvbGN76x5viOO+5oJ5xwgq1YscLuueeeTp/r2bOn3XXXXXbQQQfZVVddZb/73e/skEMOsb/+9a+2//772y233GKnnXaa/elPf7Jf/OIXNnDgQPvtb39rhx12WE0+P//5z+2iiy6ybt262RVXXGEXXXSR3X777TZ27FjZpGxB+Eu2jSUX9340e/bs8CVkvb70zne+0z7/+c/bpz/9aTv55JNtu+22W+/1Dj30ULvjjjvM7J9etffff7+97W1vs5tvvtlGjhxp999/vw0YMMCuuOIKu+GGG+zYY4+1m266yXbeeecyjx133NHuuusu++xnP2tz5syxiy++2C699FJ77LHH7LTTTrPRo0c3VRdi07P11lvbG97wBjMzmzBhQnjORRddZLvttpudc845Lb/+iBEjbOTIkbZ06VJbtGhRpc/U65OKcfHe977XTj/9dPvIRz5iy5cvX++5//Iv/2LLly+3v/71r7bvvvvaJz7xCfvc5z5nb3vb22piogqN5rW+eZRPlIcNG2aDBw/u9PdTTz3VbrjhBlu6dKn95Cc/sWnTptn73/9++8Mf/mCHHXaY3XnnnbbtttvapZdeanfddZedeuqp9qc//cm22mqrmus/8MAD9v73v98ef/xxu+iii+yXv/ylzZo1y8466yzba6+9Grp/seno16+fffjDH7YvfvGL9uEPf9jGjh0bnrehcwziP/ivWLGiHGOJW0jQb7N79+72xje+0R599FGbPHmyHXroofbpT3/aPvvZz9ob3/jG9frODx48uJTO//znP7e//e1v9trXvtZuu+02Gz58uN1777126KGH2m9/+1u7+uqr7YADDrC//OUvtvfee9fkc+ONN9o3vvENe+655+ySSy6xiy++2O677z475phj6vqIi67Dhz70ITNb9wPmxqYVc583v/nNZpbPrzRfEcTfbZxyyik1z2qzf9qvmVknP/qMPn362GGHHWarV6+2KVOmNHT9ju9WttpqKzvllFPs5ZdftltvvbXT57p162Y33XSTveENb7Drr7/efvnLX9qwYcPsd7/7nZ1wwgn229/+1v71X//VbrvtNrvkkktsl112sR//+MedXlx/85vftF/84hfWt29fu/rqq+2CCy6wm2++2fr372/veMc7Kt2D6Bo0856wlfPxVvHe977XdthhB7vmmmvs2WefDc9p53F8k3lybyhvectbwi9DZv/00OMg96lPfcp69uxZOe/x48eXvyKamS1btsz69etnQ4YMqfGGMlv3S+DIkSOrF97M3vOe99h73vOemmPXXnutfehDH+rkM/TlL3/ZjjjiCPv5z39u73znO+2JJ56wXXfd1d70pjfZtttua29/+9vDlyYHHnig/eQnP7GPfexj5YT/pptusl/+8pc2btw4u+uuu+y4444rfxn95S9/aXfccYd9/vOft7e+9a1mZrbrrrvamWeeaQ8++KAddthh9o9//KPmGlr9tOFsyljOGDhwoJ144om2YsUKu/322yvn3wjLli0zM0v9/ur1pd/+9rc1/7948WL7+Mc/br/73e/C89/4xjfau9/9brvyyivLY5dccomdc845dvfdd9sPfvAD+9a3vlX+7ctf/rJ94xvfsHPOOccuuugiM/vnQ2n48OH2wx/+sNOGUd26ddOmlxXZHDHeq1cv+8QnPmFbbbWV9enTx173utfZiBEj7Ne//rX96U9/Csv4vve9z84555yGN5uMOPHEE+2oo46y7bbbzoYMGVJ+AfzgBz8YvoCJ8F/+b7zxxvDvivGuxaaO84EDB9qFF15ov/zlL0NvYDJ48GDr06eP3X///fbDH/7QzjvvvJq/L1u2zM4+++xO/q+tysvH/2iFRo8ePcq5xMiRIzt51J566ql24oknls+mrbbayv7617/a6173Ovvzn/9sH/7wh8M+8OY3v7msl7e//e3Wq1cv+9SnPlXGvrPTTjt1mtuImM0xlp900kl20kkn1RwbN26cvfe9760Zqzd0jkGGDRtm2267rU2aNCncv8FXXu2zzz7lsYMPPti22WYbmz17tv32t7/t9FJjzpw59va3v90efPDBTvkdd9xx9n/+z/8Jx+v77rvPrr766nAu/+lPf7oct8eMGWNHHnmk/eEPfyjn8M5WW21VKh/E+tncc/LXvOY1tv/++9uUKVPCvZtaTbNzn3POOccGDBhgu+yyi40dO9Ze+9rX2uzZs+0LX/hCeL7mK12fTRn7N9xwg/3ud78r/X9vvvlmW7t2rR188MF21FFH2UUXXWT/9V//FeZ18MEHl+9ABgwYYG9+85utR48edu6553ZSPmR897vftTe96U327//+73b88cfbww8/bNttt52ddNJJ1rdvX/vgBz8YKp/79+9vDz/8sB133HG2du1aM1v3/uSaa66xGTNm2JgxY8qXhL6p4Be+8IVS7WD2T2/7+fPn25gxY2zVqlU11+jVq1elexA5Xfk9YSvn463Ef1xdn5qn3cfxzeaV0ojXzvo47bTTwnyr0tH/76c//WlRFEVx9dVXF1tvvXV5vHfv3sXs2bOLoiiK1atXV7rH0aNHF5/73OeK/fbbr9h5552LXr16Fa9//euLhx56qCiKorjjjjtCP9l+/foVf//732vKuWbNmuI73/lO0bNnz07nF0VRvPDCC538Arfeeuti7dq1RVEUnXwxzayYOXNmMXPmzPL/u3fvXhTFPz0JN2dstNu/rhrL0b/tttuuuOOOO4qiKIrPfvazdc/fEP+/KVOmFEVRFJ/85Cdrjh9++OHFiy++WBRFUdx44401fzvvvPOKN77xjUW/fv2K7bffvthnn32Kb37zm8Xq1auLl156qXj9619fc7773t9+++2drn/00UcXRVEUM2fOrOnLZlYMHDiwKIqiuOyyy8pjb3rTm4qiKIpvfvObmz2mutq/rh7jI0eOrDnv5ZdfLr773e+GPmd77LFHsWTJkk5ekkXRuCe3//v2t79dc/2FCxcWJ510UuXPf/zjHy+KoigefvjhTmVWjCvOt9pqq2LcuHHF/Pnza+YA2fj86le/uiiKonjxxReLFStWFB/72MeK3r17F3379i0++9nPFmvXri1WrlxZjBo1qm6dNJPXEUccURTFPz25O9ale2MWRVGceeaZ5fH3vve9RVEUxRVXXNGpDGeddVZRFEXx97//vdPfjjnmmKIoar02P/GJTxRFURQf+tCHNntMdbV/XTXG+/TpU3z9618vDjrooKJHjx5Fjx49iqOPPrq45ZZbiqIoiqlTpxY77bRTzWeamWNE/w4//PCiKP45J4/+Pnz48KIoimLy5MnlsXe+851lv1i2bFlx5plnFj179iwGDhxYfOc73ymKoiiWLFlS47HsvuPReL333nsXRbH+ufytt95aHhszZkxRFEXx61//erPHVFf711VjPPp32WWXFUVRFP/f//f/Vbq3DZmTb8jc55577qm5t/vuuy8sg+Yriv31xf75559fjs3OTTfdVBx22GFped033nn22WeL97znPQ3XTY8ePYrf/e53NXm9/PLLxU9+8pNiwIABaV0OHTq0099mzJhRFEVRHH/88Z3+duuttxZr166tifFly5YVM2fOLLbbbrvNHiPt8q+rxnKj7wlbOR9fXz014snt8+YJEyaEf99CxvH2CN5NaSjft2/fYs6cOWXj//CHPyx++tOfFosXLy4efvjhoiiKYuXKlRt0je7du5cD5KmnnlrztwMPPLCYPXt2cf/99xdHHnlksfPOOxf9+/cvPv/5zxcvvvhiMWXKlGLXXXet+UxR/PPlSHStefPmFcuXLw//dscddxRr166tOXb99dcXRVEU48ePL77yla8Uxx13XLHjjjtutjhph39dNZY7/tt6662L3/72t0VRFMVVV11V6TMbMqE+/vjji9WrVxdFURR/+9vfiu9+97vFVVddVaxZs6bsS3/+858r5cWXgDzuA/EFF1zQ6TPDhg0riqIofv/733f62zbbbFOWy4917969mDdvXvHyyy8Xf/nLX4pzzz23eNWrXtVpEH8l/munGN97772LT37yk8UzzzxT3H333cVuu+1Wc871119fLF++vNhrr71qjhdF8y+5/d9OO+1UHHjggcWvfvWr4uWXXy6+9KUv1f3M6aefXrz44ovFwoULwx8jFeOb7l9XjfPPfOYzRVEUxSmnnFJzPBufX/Oa1/j8vPj85z/fKb/vf//7RVEUxU9+8pO61242r5/97GdFUfzzC+lll11WfP/73y/uu+++YtWqVcUTTzxRFEVRvPOd7yzP95fcHV9Ymllx4okn1u0DP/3pT8tjAwcOLJ577rli7dq1xbXXXlt86EMfKkaPHr3Z46sr/OuqMZ7922abbcoXbB1jo1VzjGZecvtGlUVRFGeccUanz1x77bVFURTFF77whfKYv+Re33i9vrn81KlTy//feuuty3u88847i3/7t38rDj/88E6bGr8S/7VLjO+6667FCy+8EG44mf3bkDl5K+Y+u+++e/Ha1762uO+++4pnnnmm04/5mq9s3n9dNfa333774re//W3x3HPPFR/+8IeLPffcs+jevXtx8sknF1OnTi3WrFnT6X1IlMeoUaPKjf8uvvjiytcfNGhQ8eijjxZTp04tTj755KJ79+7FnnvuWXz4wx8unn/++WLRokXF4MGDO9Xl+t6fFEXR6X2MmRW//OUvi6Ioin79+pXHLrzwwqIoimL69OnFt771reL1r399+Fn96/qx3Oh7wlbOx9dXT4285P7Vr35VFEVRfOITnwj/voWM4wre6N8ee+xR/OhHPypmzZpVrFmzpliwYEFx0UUXlQ07e/bsDb7GN77xjaIoiuL73/9+TXBMnTq1WLFiRbhT+wUXXFAURVGcf/75NceLoijGjRuX1vOsWbPCv/lkicd22GGH4vzzzy9XyHhnveKKK4o99thjs8VLV/7XlWPZ/2299dblLrq/+c1vKg+GG7qT+4EHHlhce+21xZIlS4o1a9YUkyZNKs4777ziHe94R1EURfHzn/+8Uj7bb799qUrgKicfiDv2CbN1XyizFQVRv+nfv39xySWXFEuWLCnjf8mSJcXXvva1hnY/3tL+tUOMd/x3xhlnFEVRFD/60Y/KY74iNNpFvSg2/CU3//35z38uXn755eKQQw5JzznttNPKZ8w+++wTnqMY33T/umKcjxgxoli1alVx6aWXdvpbNj6PGjWqbNsoro488siiKPIXaq3K60Mf+lBx//33Fy+88ELx/PPPlyu2/vjHPxZFURTHHXdcea6/5I7qtZk+MGrUqOI3v/lN8cwzz5TlnzNnzha5i3wj/7pijNf7d8455xRFURTXXnttp7+1Yo4xevTooijyVU0HH3xwURRFce+995bHTj755KIo/rkScIcdduj0mXe/+91FUdR+EWxmvGa7dZzL9+zZs7jggguKuXPnljH+7LPPFhdddFGx8847b/Z221z/2iXGP/axjxVFURRXXnll5c80Oydv9dynR48exaJFi4rFixfXxL/mK5v3X1eN/fPPP78oiiJ8/u6///5FURTpu4ro38UXX1wURVG87W1vq3S+95uxY8d2+tsnP/nJMC4bfX/SsX7ZBltvvXXxqU99qhg/fnwZx2vXri2uu+66pr9fb+n/umosmzX2nrCV8/H11VPV9zq77bZbsWrVqmLFihVFjx49wnO2hHG8bTy5G6EVnmlLliyxc889184999ya48cff7yZmT3wwAMbXM6lS5eamdWYs48aNcpGjBhhDz30ULhT+7hx4+zTn/60HXzwwRt8/YzVq1fb17/+dfv6179uAwYMsGOOOcbe97732VlnnWWDBw+2Y445ZqNdW9TSKv+/bbfd1n7961/bO9/5Tvv1r39tZ5999ibzJB0/fry9/e1v73T861//uplV70tr1qyx559/3nbffXfbeeed7YUXXmhpOZ0FCxbYBz/4QTMzGz16tJ1wwgn28Y9/3M4//3zbeuut7atf/epGue4rlY3lO2+2brOw4447rjz2qle9yszMrrjiik67SZuZDRgwoPRB7dmzZ7oZRxVuvPFGO+WUU+zYY48NfVnf/va325VXXmlPPvmknXDCCTZ9+vSmr9UIivFNz4bE+ejRo22HHXawD3zgA512bHc8dt7ylrfY9ddfbzNmzLAXX3zRunXr1mnfDzOzp59+2sz+uSFMPTYkr5/97Gf2s5/9rNPxsWPH2ssvv2wPP/xw3es3y+TJk+3MM8+0bbbZxg444AB77Wtfa+eee65ddNFFtmLFCrvssss22rVfiWzMsTyaLzOfDZ1jzJgxw1566SUbOnSobbPNNp18uUeMGGFmZlOnTi2Pucfn6tWrbfXq1Z3ybKSPNcszzzxjn/nMZ+wzn/mMDRs2zI499lj7yEc+Yueee6717NnTzj777I127VcirY7xKp6oraLVc59nn33W7rnnHjv99NNtv/32s4ceeqj1hTbNV7oKGxr7vrnkuHHjOp07YcIEW758uQ0ePNh23333uptqm/1zfv/Rj37UjjvuuHSvJmeXXXax4447zp566il77LHHOv3dy7Qx36384x//sAsvvNAuvPBC69Onjx111FF25pln2jvf+U7bb7/9bL/99is9v8XGZVO/J2zlfLwV+IaTP//5zzfoO24jbI5xfIt8yX3eeeel5vMRP//5zytPtH3CSAP2ZnnNa15jZrW7s7r5eu/evcPP9OnTx8xskw2E8+fPtyuvvNKuuuoqmzJlih199NGVH0Biw2lFLHfr1s2uvvpqe8tb3mK/+MUv7P3vf3/ljfA2Fttuu639y7/8i61du9auvfbaSp/ZZ599bPfdd7fnnnuu3HBqY/PEE0/YE088Ydddd53NmzfP3vKWt2hC3WI25njdv39/MzN76aWXymP33HOP7bLLLuH5H/zgB23FihV21VVXmZmVm/Q2S3R9513vepf94he/sAULFtjxxx9vs2bN2qBrNYtifNOwIXE+e/Zsu+SSS8Lz3vjGN9pee+1lV199tT333HPlJo4vvvii3XHHHXbCCSfYmDFj7NZbb6353JgxY8zMKsVdK/MyMzv22GNt0KBB9r//+7/23HPPVfrMhuAv0x9++GG7++677Y477rC3vOUtesndYjbmWB7Nl9dHo3OMNWvW2N13323HHHOMHX300Z02ADzllFPMzGpif9asWTZjxgwbNmyYDR06tFPZGu0XG8qMGTNsxowZduWVV9qSJUvstNNO2yTXfSXRyhh/9atfbQceeKBNmTLF/v73v7eohDkbY+6zvjnOxkDzlc3Hhsa+v9/w9xhku+22s+7du5tZ9fcbjcTedtttZ2Zmu+66q3Xr1s1efPHFmr9v6ncrS5cutT/84Q/2hz/8wXr16mUnnniijRkzZqP+6C/WsanfE7Z6Dr2h+I+rP/3pTzfJ9TqyqcbxLfIl95AhQzbo81tttZXttNNOtmLFiprj73nPe+zss8+2u+66y6677rqav/Xq1ct69+5ty5Ytq9np9+CDDw5/3X73u99tZ5xxhq1Zs6Zm992JEyfa008/bYMGDbJzzjnHLr300vJvPXr0sM9+9rNmZnbLLbds0D1m9O7d2/r27WsTJ06sOb7zzjvbLrvsYi+++KJ+adyEbGgsb7fddvb73//e3vjGN9oll1xiH/7whzfKC+6+fftajx49bNGiRTUvLXbaaSdbvXp1zarxbbbZxi666CIbMWKE/cd//EeNYmHw4MH27LPPlr9qOr1797bLL7/czP65S3LHVVatYvTo0bZs2TJbsmRJzfE999zTzMxWrly5Ua77SmZDY/yggw6yRx99tJMyYeedd7YLL7zQzP65q7tz9dVX14y55IMf/KA9/fTT5QSA7L333rbTTjvZ3Llzy53Rt9tuOxs1apRNmDCh0/mHHHKIffSjH7WXXnrJbrzxxpq/nX322XbZZZfZnDlz7Pjjj7e5c+c2dtMbgGJ887Ahcf7oo4+GMWn2zxVIe+21l33pS1+yGTNm1PztRz/6kZ1wwgn2f//v/7V77723bNsePXrYV77yFTOz8qWGE8V5s3l1797dnn/++ZpjAwcOtEsuucTWrFljX/7ylxutisq86lWvsunTp3d6ia4433i0YiwfP358pznKCSecYJ/+9KfNzOxXv/pVzd8anWOY5TF+8cUX2zHHHGP//u//bieeeGL5ou+QQw6xM844w5YsWdJpxeCPf/xj++EPf2jf+c537MwzzyznJv379y/L/Jvf/GZDqiVl8ODBttVWW3X6YrzbbrvZ9ttv32keJTacDY1x8uEPf9jMNs6LhmhO3szcZ++997Y1a9Z0mi+Y/bP8r371q23u3Lnh6thWoPlK12FDY/+OO+6wsWPH2pe+9CW76667at4lfO1rX7Nu3brZ/fffX6PSPfTQQ0MlztChQ+1LX/qSmdXO783i2F++fLk98cQTNnr0aPvKV75S80Jt++23L+ciG+vdynbbbWeHHHKI3X333TXHt912W9t9993NTLG8Kdkc7wlbOR/fEI466igbPXq0PfbYY3bPPfe0JM96bK5xvG1ecr/lLW9Jf3WZPXu2/eIXv2jZtXbaaSdbvHix3XTTTTZjxgz7xz/+YUceeaQdccQR9sQTT9g73vGOTpPwT3ziE/a1r33Nvva1r5USSTOza6+91l566SV78MEHbf78+bbDDjvYoYceaocddpi9+OKL9pGPfMTmzJlTnr927Vo777zz7PLLL7dLLrnEzjzzTHvkkUdst912s1NPPdX22GMPu+eee2pefreS/v372/jx423ChAk2YcIEmzdvnu266672pje9yfbaay+78MILN5pNxCuFTRnLP/nJT+yNb3yjLV261BYsWBD+Unbbbbd1WkXiL5TN/mmhY2b2ne98p3xhcckll9hdd91VnvPtb3/b3ve+99n73ve+mvIff/zxdskll9jNN99s8+fPt1122cVOPvlkGz58uF1zzTXlwO4ce+yx9pOf/MTuvPNOmzlzpi1fvtwGDhxob3jDG6xnz572wAMP2Oc+97kNr5iE173udfa9733P7rnnHps6daotWbLEBgwYYKeddpq9/PLL9r3vfW+jXXtLYlPG+Fe/+lU78sgj7e6777a5c+faypUrbe+997ZTTjnFdtttN7vrrrvs29/+9gZf54orrrDjjjvOjjvuuLK/7Ljjjvboo4/ao48+ahMnTrT58+fbTjvtZPvuu6+dcMIJZmb2b//2b6Ws3eyf1imXXXaZbbPNNjZu3Dh7//vf3+lazzzzTPmCvtUoxlvHpozzZrjuuuvssssusw984AP22GOP2V/+8hfbZptt7E1vepMNGDDArr322k4vDaM4bzavSy+91AYNGmQPP/ywLV++3IYMGWKnnnqqdevWzc4666yN9mLEzOyss86yj3zkI3bnnXfajBkz7Omnn7Zhw4bZm9/8Zlu9erX953/+50a79pbEpozxCy64wEaMGGF33323zZ8/38zM9t9/fzvxxBPNzOzLX/5ypy9ljc4xzPIY/81vfmNvfetb7R3veIc98sgj9sc//tF69eplZ5xxhm2zzTb2oQ99qNOPNj/60Y/s5JNPtre//e02fvx4u+WWW6x79+72lre8xXbffXf7wQ9+YLfffnvL6ogccMAB9vvf/94eeOABmzRpki1cuND69Oljp512mm233Xb2ne98Z6Ncd0tjc4zj3bt3tzPOOMNWr15dKf9Wzckb5VWvepVdc801ds8999j06dNt8eLF1qtXL3vNa15j+++/vz3//PN21llnbTT7Q81XNi6bMva/+c1v2pvf/GZ77Wtfa5MnT7Ybb7zRVq1aZUceeaQddthhtnLlSvvUpz5V85m//e1vtmTJEnvkkUds3rx5tu2229qwYcPs5JNPtm7dutlFF11kN998c81nstj/5Cc/aTfccIN95Stfsde97nV2991324477minnHKKDR482KZNm7bRxswdd9zR7rrrLps2bZo99NBDNmfOHNthhx3sda97nY0ePdquv/56mzx58ka59iuFrv6esJXz8V69etn3v//98v/d/eHSSy8tr/sf//EfNd89nY3542rG5hzHN6lRO/81Yii/PrJNWpr9t+222xaXXHJJMXny5OKFF14oXnjhheKRRx4pvvjFLxY77rhj+BnfUKGjQfvnPve54m9/+1sxd+7cYuXKlcWqVauK6dOnF5dddlmx//77p2U4+uiji9/97nfFwoULi7Vr1xbPP/988eCDDxaf//zni+23376SyTvruerGCT169Ci+8pWvFLfccksxf/78YvXq1cXChQuLcePGFWeeeeZmi5Wu/q+rxrK37/qINhWoR8cNHrKNH0aMGFFce+21xdy5c4vVq1cXy5cvL2699dbiXe96V1jeMWPGFJdffnkxYcKEYtmyZcXatWuLp556qrj99tuLT3ziE0W3bt06faaVmyOMGjWq+MEPflA88MADxZIlS4rVq1cXs2bNKq655pri8MMP3+xxphjv/O8Nb3hD8ctf/rKYMmVK8cwzzxRr164tFi9eXNx0003Fhz70oYZ2my6KfPMl70vHHntseWzbbbct/s//+T/F3/72t2LevHnFqlWripUrVxbTpk0rfvGLXxSvfvWrO+Xjm+utj47jtWJccZ79q7IJ2TnnnFNu/rhixYrigQceKD72sY8VW221VaU4bzavs88+u7jzzjuLZcuWFWvWrCnmzp1b/PznPy9GjRoV5t3KjSdf/epXF//93/9djB8/vnjqqafKfnnZZZcV++2332aPM8V4538f+MAHij/+8Y/FrFmziueff75YvXp1MWfOnOI3v/lNcdRRR4WfaXSOUS/Gt9lmm+K8884rJkyYUKxcubJYvnx5ccMNN6x3bOzWrVvx2c9+tvzMc889V9xxxx3hnLmVG0/279+/+OY3v1nceeedxaJFi4rVq1cX8+bNK/785z8XJ5988maPM8V4/u+jH/1oURTVN5ysR9U5+fryj+Y+e++9d/G9732vuPfee4snn3yyWLt2bfHcc88V48ePL773ve8VAwYM6PQZzVcU+9m/3r17F9/73veKJ554oli1alWxZs2aYvbs2cVll11WjBw5stP55557bvHHP/6xmD17drFixYrymXD11VcXJ5100nrvLYr9sWPHFldccUUxZ86cYs2aNcXKlSuLiRMnFt/85jfDDfhatfHktttuW/zbv/1b8ec//7mYM2dOsWrVqmLJkiXFPffcU3zkIx8Jv9vqX9eN5WbeE/q/VszHfRxdH9H8pmfPnsXKlSvXu+Gk/9sSxvGt/l9CCCGEEEIIIYQQQgghhGg7tt7cBRBCCCGEEEIIIYQQQgghmkUvuYUQQgghhBBCCCGEEEK0LXrJLYQQQgghhBBCCCGEEKJt0UtuIYQQQgghhBBCCCGEEG2LXnILIYQQQgghhBBCCCGEaFv0klsIIYQQQgghhBBCCCFE26KX3EIIIYQQQgghhBBCCCHaFr3kFkIIIYQQQgghhBBCCNG26CW3EEIIIYQQQgghhBBCiLZFL7mFEEIIIYQQQgghhBBCtC16yS2EEEIIIYQQQgghhBCibdFLbiGEEEIIIYQQQgghhBBti15yCyGEEEIIIYQQQgghhGhb9JJbCCGEEEIIIYQQQgghRNuil9xCCCGEEEIIIYQQQggh2ha95BZCCCGEEEIIIYQQQgjRtugltxBCCCGEEEIIIYQQQoi2RS+5hRBCCCGEEEIIIYQQQrQteskthBBCCCGEEEIIIYQQom3ZdmNkutVWW23wuUVRhOd4OjpmZrbNNtuEafLyyy+bmdm22667/dWrV4fX7tu3b5n283ku2WGHHcp09+7dw3P8s8uXLy+PPfvss52uYWb2j3/8o9Pnt946/l2C52Z15/fNv7eajZl3V6VeDDfSHzbkehFZvPA4015mtmPWpuxfUZwxJqNYzsjuL8uD5zcSf6/EWN2YNDqWb8xrVi0HP18lvrL8ujLtUs4tmSpxurH6ycaM3yyPVpY/qxfFdddiQ9p8Y82VWk0Uc1lMtvpeFO+bD2/LKuNdvfFqc8b4xhxLG/kepFjuumxIO3Xl8bte3220bwvxShjL2z3mtZJbCCGEEEIIIYQQQgghRNuyQSu5N9WvGNF1fLVox79HK0o7Hn/ppZfMrHbl9ejRo8t0z549y/Suu+5apl988UUzM1u8eHF57LnnnivTXOHNz40dO7ZMd+vWzczMnn766fLYnDlzyvTs2bPL9LJly8q0r/D2sndMZyt2Sb1V8PzFhqsao7yr/LrT7r8ArY9Nsbozu0aVX9k8XqqoHKIV11VWclN1wHO83/HaVeKz2dXe0S+azf7KuSX9grmpqPcc6CqrO+qNY1X6W/bZZtiYq/7ExmFjqhU2Zww0O+5tijKrb3QNGlmx3OgK/3Zp4w1VDjWK5iCbj1a0db322xRKGF5nQ75T1PtcdL16x0TXo0o7bepxsBVkMd5sHxXtSb2V+12dRsbWZudm7Vo3EVrJLYQQQgghhBBCCCGEEKJt0UtuIYQQQgghhBBCCCGEEG3LRtl4shVQWh4tuXfbj47QooR58Px+/fqZmdmIESPKYwMGDCjTPXr0KNNr1qwp05MmTTIzs6VLl5bHVqxYUabXrl3b6VwzsxdeeKFMH3zwwTVlMDNbuXJlmaYFCS1P/Jq83o477hh+jrYp0eabVeRD0WaEGbSZaBfZ6YZSb2PUKjQi9cv6A8+N7Gt22mmn8pjb7XTMg3l7HlWsQ5hfVFbmy37JdPS5jCpymmalPI1sKCW63oZKWVvV26i30X68KeyINjWK8+ZoVobb6Of8/OyZ3Ih1RKNxXW9czGj2c42guO0aNGsB1lXGv41FszZpzeYhujZRG29Oy89mP9/Is6dKHqI9aEcbtlaguH3l0qwt6qYqU4Q2nqxFK7mFEEIIIYQQQgghhBBCtC16yS2EEEIIIYQQQgghhBCibdnKzBpei76xluVn+dJqJJLGd+/evUxT0rvddtuVado27LHHHmZmNmzYsPAazz33XJlesGBBmX7yySfNzGzJkiXlMVqUZHLi7bffvkzvvPPOZlZrj9KnT58yvXjx4jL9/PPPl2m3HVm2bFl57Nlnny3T2267znmGthWRfQv/nkn0I0laFdlCu0sbqrIppFuZFIYWNDzusWVmtssuu3T6O61uGAO8jtvz0K6kiv1DRFbmzEbI+yvjmtY7jPHMTqURmXwjsfpKieuqNGuL0IrrMZY8zWtn9kkcC1vRns1aMnj88/MsW7PS9WZ3jxfVadZ6oVHLkHrj7Ctth/R6vJLudXPRbIxvzOfC5rTKqnftRvroprDHEjHNxs6War3TCouqV+pzqB1otm22JNupVtupKca3TLrifKWROcOmsDbp6mgltxBCCCGEEEIIIYQQQoi2RS+5hRBCCCGEEEIIIYQQQrQtXcquhFCSTnbYYQczM9txxx3Dv++5555letCgQWWaNiAvvPCCmZn17t07vB7tEmbOnFmm3crBrSA6fo7WCn6uWSyZp7UE07wv2jM4lDHTVoWyg5UrV5bpF198sdNn+fdM2hBJ/jMLi0albO1KK+K+Eak6259x5n3ALLbCMVtnR8I4pEUJrUsYIx6fmY1Ndi9ZbDi0JeHnGO+77757TRnMamN11apVZZp9LbIuyexMsvLXi88tIX43lI0Z//Wg9RT7grcL25vxQyK7klbLcKvk532B57IPZlSJ6agcsi5pHVXGjFb0Ez7no+dvRldp32j3+WZp1OpFNEfUZq2QxLaCRp7VG9Mqodn7rvc5xe+mpdlxpBXjf6v7TqvjXVYOWzb1Yr8d7UdIo/2yWfsHxXV70sjzt937AmmFNVE7opXcQgghhBBCCCGEEEIIIdoWveQWQgghhBBCCCGEEEII0bZ09sPYiNRb+k97Bp5Lefquu+5qZuusDcxqbRh69uxZprt3716mn3rqqTLtdg/PP/98eWzOnDllmrJg2inQLsGh5QltTmgTwfzcXoL3R9sRytYp0Y9sTmhbwTro379/maZNi1tUsI5Y5kwO3Wr5WrvSCklLlgfl6TvttJOZ1cYer5FZNzCePaYy2xHafTAPL0cVaXxkadMx7TA+s7yXLl1aUwazWmsg3jevEdnvMA9Srw2rtHG7y3eapRUy28iWJ7OmYhvS7mbt2rVl2uMqK1sj5ayye3xk/cFyZrERWaxEddExnVmvNGsH0UobiVcSURxVia2ovrN2J9HY2mybbU6Lj1Zcu5F6rpK3xvL104r6btbmpFlpeBVLtU3d1o1YYcmGZ9PS6vlCI3RFi5J6eW/I3K8R+xax8WnEmrHKnLirsSEWhNHfm72O4nrT04pxsRXP7a7CK9WihGgltxBCCCGEEEIIIYQQQoi2RS+5hRBCCCGEEEIIIYQQQrQtle1KNpaMi+dSmu2WDWZme+65Z5l2iw7KeQcMGFCmaXkyceLEMs3zXe6e2YSwHPycH2eZM6k6y8H83DKEn6MtBc/lOV6Op59+Ovw75fy0VWHdeZmWL19eHqP0n3YX0bVfiWwqWaDb8JitiwdaijAu2GY8zvK5Lc4LL7wQXi+zf/B2Z/vXs2vI7iuThNF2hHjc0tqEsez3ZFZrB8Tju+22m5nV3jcteZq1lunqsqRNQbMS9UwOXE8myNjObDs8Tjl+k3pWKFXiOSOymqgiffZ7yex+GrEBIq22gNiSJGtdiUbtIJxm7bKqSPObbetmpc+bk2YtMbZEWvE83FDrpOwazcZOq9txY/WNVlxDNMfmHIs2pu1hszQyl1esth/NzvdaMfa1wjJwY1Lvu4xiv2vTaMy0Ys7aVeayEa2I0Ua+63f1PqGV3EIIIYQQQgghhBBCCCHalo2+8WQjb/y5wrNPnz5lukePHmXaV+ZxJee8efPKNI9zpTZXe3p+PMYVf1xJx40b/XN77LFHeE9cEcv75spUPz9bmcjNMrka3De45KpvXxVuVrvql6tfueFmv379zMxs3333LY9NmzatTHNlIlfCet00uqI7+kWsq//qsz6a3XAl+2WQq6l53OMo27w0q0PGs7clY4irXXmcab8Oy0YYy9Fmesyjympw3pf3L6o4GMtRPzKr7Qc+hnBlPD/Hc6MytdMvlJuTZldnRDETqWU6fo6xxDz8nCw26q2mZp9hP8j6XkS2kWCmUvK+xY2Ds1XrvK/oXljmbBV8vTKTrrxCYVOwse4/U0ZV2dipnuogW53t12n16u1GNhXMnhFZfUT9f2OOyRrjO9PspnzNKhSqfrbqua1WIdbro1Wu14oNoJrdIE10phWqgyy/RvpJK5QLzW5aurFiMjtHm2RvPlodW1Wu09Vox427RX2qfJdpRHlc5TpbyvekRmO7HcdkreQWQgghhBBCCCGEEEII0bboJbcQQgghhBBCCCGEEEKItmWDNp7MlrpHG3KZxXJV5kFbEh5/5plnyvQuu+xiZrWyVtqS0GZh4MCBZZoWCG5ZsHTp0vDatEph3i69HTlyZHls9913L9OPPPJI+Dmm/V5obdK3b98yTfn5nDlzyrTfNy1daD9B6Tslk7Qd8fvlRnzMj+WkLYpbRlSRDbeTjKEj0f21WtafybfZJr169TIzs5133rk8xk1B2Wa0M4hig32RscXPRRYLHm8d/75s2bJO55rlG2A6mTUQ49MtIiL7FLPaPpxt0OrHGcv8HPsMLSkiS4AtRZK0OcjGisi2I7OvyTaNJH5+FjP1nlHZho8ZjUjgMoshTzP+WC+8Fz7Poj6WyZ2bley38/jdCM1KxhvdxLEReV+9Z2q2yXU9Insfs9wKpxGZe70yZ/ZWjEmmo3tsVB4vui714qXZ7xQbYpXS7IZrHreNblS1pdn4bUm0wpakXltmNmz1nklVnjetsCAh9a4dnVvlelU+pz6xaWnlJqOt/u60KWwkN5VVpWynNh6NPn8jGhnj2ynOWxF37RinWskthBBCCCGEEEIIIYQQom3RS24hhBBCCCGEEEIIIYQQbUtlu5JGlqlnMu1IQk1pNi0ZaN+w6667lukRI0aYWa30dfbs2WW6Z8+eYZko+/a8n3rqqfIY00OGDCnTH/jAB8r01KlTzcxs7733Lo/ttttuZXrhwoVlevHixWWaViK0RXGGDh1appcsWVKmaa3gaR5jHVGKwPqihcXq1avNzGzRokVh2QYMGFCm3dKFadZhlV3kN5X8Z2PQShlKJhfPrBncXoN/Z99gHrTioMWM2+FQLk6LnP79+5dp9hm3yIni1KzWOijDyx3dU8fjHpMs/8qVK8tjtHTgvfBeWR/ePzJLFFJPgtmoVU27xThpteyqiq2DU8XeirBt/bOMI36Otjt81ixfvjwtT8drZLYOfm3208xehHlEf8/uNRq/s2uTRuTH7Ry3jdCK3dSb7fuRpZJZ3n5RnPFzhGMh84tsxjJbn8ziJ7IMaURiT1j+TLIflaMRux1+7pUS143S6udXvWdnK+S/2RzKn/dZPGX9rh717EV4vIq8v4rdQ73PNVJOsX6arbcNqe/oGZKNq0z73KHKmLmxYqBKvq2YQ8rKYeNTbzxu9ZhSJS7qWf9tCprtz43mXe9eFe8bl0bGmEbeZzVqZdhIOVvxuVdqXGkltxBCCCGEEEIIIYQQQoi2RS+5hRBCCCGEEEIIIYQQQrQtle1K6lFluX8kx6ItwqBBg8q0W4OY1Uq2XX4+ffr08HOUi9OKg+Xw/CiBpFXHscceW6bHjh3b6TjlALSRoJz83nvvLdO0D3GbiJ133rk8RnsUWjXwnO23397MamXvrLtIym5Wa0Xhtg2sC5c0m5n16dOnTEc2Mrx2lfbekuURVXbx9XbI5OK04uBxtxWhNY23v1ltjLOteb5LGNkGe+21V5l+9atfXaZ33333Mr3HHnuYmVnv3r2j26sp54IFC8o0r+Nxxn6bSeoj2SWtTdiHGau0Xonq4+mnny6PsZ+z/Cwf6855JUjINpbEtBG5VhWpOceeyOaH5zJ+aMsTWYb06tWrTNPaZNasWWWasRFJ0jL7hszyxM9nn+Y9MbYpsY+sULK+1KhdQD3auS/Us69oRZ3UOyd7BmTtxNiJrByYB8dIPlOiOOOcguMpPxfFe9YvM5sTLz/vqYo9T3Q8qov15RGVp91itllaYafQaN5RfDbSN7I8spgj3g8YIxwnOZYyrutJlqvcS7O2OK2Qq79S4rkRmrXhaTRWG7lGFEdVnsfRcz2bW2QWiM0SlXlD8m1kfGjFXGVLppFxokodN2JdkuXRrP1SFFPZ+E+ieXyjbOq5Qb161Hje9W1AG7EkrPe5LM55bj0LwSpsauuhrtJWRCu5hRBCCCGEEEIIIYQQQrQteskthBBCCCGEEEIIIYQQom1pyq4kWrKeLdXPZIRu4UG5OGXmtNmgJcikSZPMzGzRokXlMVoruN2CWa2Ml6xYsaLTNV7/+teX6VGjRpXpZcuWlWm3NGG+tDw48sgjw889++yzZXq33Xar+a9ZLimgvNflxJQp066E5aDMgeVwWOeUc/Laffv2LdNuhfLkk0+WxzLJcia9WN+xrkgksYv+XuUc1gnbj/2BdegyclrX0O6D/YTloGWN58frjR49ukyPGTOmTDM+PYYZIyw/rXduu+22Mk0JmcuCs13fu3fvXqbZjz3GGcu8NuOdUntexz/Lcmb2PLxOxzJ0vHZGu8RzszQi5c3qIpPW+vjGMYhxznMzi56Ifv36lel99923TNOCxJ8ZtLqiXRNhP2R8+P3y/mgNwXIy7urJd1kHfAYQP59jAc9lf2TefrzRHe+j8a5dYr+eXLAKjVg5RPYabI9sTsRxLJIqcrzK5k1Lliwp034+50f77LNPmZ49e3aZnjZtWqcym62LYR7jM4Xw+cN77JhXR7K8o/7Feua4Ec1H2iU+W0krLHmqnBu1ZRW7j8y+0NuY7c+4JhxLfbzlnJxza3534FyY5Yjk71mM0/IkivFsvM7uO8ori/F61hciphW2F9Fcvtnrsa0Zt4w5xplfM5uHc37CuXw9i5GsnIzhyG6ritVUFKtZn8rmKq9kGv3eubEsDBotRz1rnWzcjWIj+n5mVts3ovGR8ZtZXkXnNGrJ08icMOpTTGscr6WReUyj9l6ttOWpYkMYWVVmFqr8HGO7kdisV3cbYmMnuxIhhBBCCCGEEEIIIYQQYhOil9xCCCGEEEIIIYQQQggh2pbKdiWRTCtbxs50JrHdeeedO+VHmfb+++9fpp966qky7XKsoUOHlseWLl1apt2KxKxWEkCJo0uEeU8sM69HOZnL3SkxzqQxxxxzTJmOdnWnrIxplpn34jIAWjbQmuXxxx8v07QV4X17/Q8cOLBTedZ3bZfGUar/3HPPdSpbRzb1DsatpBEpEoliirJb9gFCKa1LUygzZ4xTusK8Kc31+s7sEzL7B//cM8880+mYWa3cnTJJ2o54n2GM7LXXXp3+bmY2ZcqUMh1J1nbZZZcyzX7CvKN7pGUF64vlzOwYnCrx26x8anNST3bVrOSdMZVZzjBvH2PYfhmZjVN0PT5bmDf7k8cgy0ZZ2YgRI8r0008/XabZx/w5wc+xP7I++EyJZGOMfebHNJ+PXu7dd9+9PEZ5G68XWTlkVkKN2pi8konsFszi+mR8cj7AWM7axGM461OMHcaIf442PLRo41yCsZrt8O6wf/FznHfUu6cq1kYO7a2YX2YB4+fUs05b3/FXCvX6exWLmejcrK2jdiKMZcY7557RnJXjPGOc+XEuzLTHUWY7xRjnsyCSxzeK1xnnRFmMEz+nisx6S6OV9iNm+ffAKuNVM2XKvitn3yWjGOc8g3mwn9A+ysvP+2Bcs5/wu5/nvSFWDl7+aGw3i61ZWNZ2/h7ZarI4q/c9pMq7mnrfgTJb0uiabGueG70LIdm8JitHvTxYjmyO1oqYasQCg7zS47kVNNInzOrbRTbbJlXmPPXiJLL+MYtjPuvfvN6GzE2cZu0eu8q4rZXcQgghhBBCCCGEEEIIIdqWyiu5o18pqvyCws9Fm4pxpSY3fOSv1lwFetNNN5lZ7cpX5jtjxowyffjhh5fpaPXR8uXLy2NcLcQy8Tq+EpArWLNf0bPVUP7rDH+l5CqObEWVr9rmakSuKnzhhRfKNH+J57V91R/rJVp9ZWY2efLkTtfhJm38e/bLajtvPNnsL7OMsyjG99577zLNFXZcmXz//febWW29cqMZrrJm3nvuuWeZ9l/wGBeE/YurhbwcjMlsYxjGIq/jKwTnzZtXHuMKdq4wId4nuDqV12ZfzDZs8/rnSlz2k2wjEt+ELcs3+1W13Wnkl+Vs5bGnWUdUmTBO2Ia+qoOxxtXIjBOu6uMKVC8Tz+VqEZaZx/2ajFumeT1uwsv68LjiqkD+cs4N0Fg3vgorW+HC5wvrJlppyDGEfSXbsMzrf2OuaukKtGK1X7Sqgm3GFaMchzmv8NjI6pXxwj7F+Yi3ZbYRDdudcxCfF/E+Fi9eXKapjOH1OH/wfpWNhdlmYn5NjqfZcyRb9eKxn22wzXslXqdVxumustqkVWxI3Eef5fjDOQPHqHrqkAzm7W3GcTBb9Re1O1e2chznsz/DYy7bHKzeZoPZudmKqmjTKdYt09mzrN4q1y0trs2aj+16G71xjGKa462f3+jKOObncBzM5sL1nj2M62wT9npUUew1Mj8k0UpePiuy7x/16rdRxWG7xX4r5y1VxqVoFX09Jdf68o422s6e+9HzIpsTZ/P4eu8YqtyLX7Pe5pDZ53hOppaONp7PaLeY3ZTUUyw0osqp8r0ne18QnZ/FV/RcrjJ2NfJeLbtevWs0e062cjxTnm1OVY5WcgshhBBCCCGEEEIIIYRoW/SSWwghhBBCCCGEEEIIIUTbUtmuhERL2TNzdcpLKVdyC5L+/fuXxyiBpN0ApbSeNy0buBkM7SAGDBhQpmm1MXfuXDOrlTgedthhZXrBggUW4UvxWU7aU2SSScorXcLCpfrMb7/99ivT3CTP64lSHUqPWR/ZxlCHHnqomdXaB1Cis3DhwjI9evToMu2SukY22up4vN0kOI1Y8mRpt2zINltiHVMi7hYyHqdmtdJJyre5weLw4cM7lZN2NPx7ZgvhEmLGGe2AGNe0WImsFGhZQTk/yzxnzpxOZc6uwfEj2rjWbJ21BK1b+LnM6sXLwXwjSbbZlrXZU704r7Jhm7d3JtUbOXJkmabFg0v4OH5HdgVmtRY2fDb4mMw44jWyzXT9OJ9PLDP7RCah9fP5dz4DKMllfXnfY1wyvjLLKvZJvw7rhffHPGh15HYlG7LRaLvFeb2NWqpIU6PNURlbnF9wvuKWTRzTMqkfx/XI4iezXWKsMga8zJTKEj4DGFvE+1cm22T/YdqpIuHMYtyvw3Ga+fFc9gOvjyrzknpt35VjvRX2DSSyneIYnNmT+TOaz8tszsO5BO2cPG4Zv9nGY1GaYzD7A2Of8cIyRRviZbJnEj07q8j7eY73QY7RvG/2Uc5N/H6rxOeWYt/QyJw8Ox5ZMmSbVjOOItumrC757I1sOTLrJ5aD8expxgJhORuxDOK59az5qmziS6L5Du8vk7hHc6Js7t3IHKZd4r6R+XgjtjKsb8Zn9LnMDqTKho7eVlXqO5pXR2Njx/zqfRfLytxoDEfXyMrh6WgOZJbbVdXbQLjRMnVlGilzI5YaVTZ5jJ7LbJNG7Uqi/Jp9D1bFzieajzQ6Dkc2Ro1uyO7HszGkXpk3B1rJLYQQQgghhBBCCCGEEKJt0UtuIYQQQgghhBBCCCGEEG1LZbuSSBpTZQdyyvAox3UrkT59+pTHHn/88TJN+xDKVsaMGWNmZq961avKY5MmTSrT++yzT5mm9JFl9XNop0CZCeXztDxxiSblPPwcr8fl/rQr8TT/Tqkt5Vq0XHB5Oi0saL3gkmYzs2HDhpXpV7/61WXapaeUt1GeR9kEJcJuq3HHHXeUxyg7zeRypN1kYxGZHCWT5rptAuuKknPWN/PYd999O117yZIlZZqWPEOGDCnTlN94/xg4cGB5jNLjp556KvycxzDPZfl5nFYikSQ+sukxq71v5ufx/sADD5TH2L84DvB67PNuV5JZm2S7frvdwJQpUzrdh1ljOye3E42UmW3IdvM0pbmsL8aBt4+Z2bPPPmtmtTZJzIPnDh06tEzTDsKfL9mu7DNmzAiPu7SeEns+i2g7snz58jLNMdf7DfsuZfyEFipeZtpNsb5YB3ymMI79uUTJf/ZcYv37M5j5kirx3G5jebPljCyYOFYyFgif2z6XYFywfdmmbEvvG2brnsW8D7YfY5h4HuyrfK4zxvk84PPcx1FeO7OAYMx53pyXVJFwctx2W6usXlj+qByM+0Zpl9huBt4bx3SvN45FjAX2B8a4x1/2jGT8ZZZ+UXky25HIyiGzlGI8ET77/dlRRSrPuZCXNbNty2S8PMdjm98Bqti+RGXb0mlkrpJZOUT2AIwdziN69epVpn2Oko1bmUVJZB+SxQXHq0jqzWMsc2ZREknHMwl5Fp9+X+x/kb1Px2sz7fM1ztsy65V6FgONSvTbgXr2AVXmZPXOYbxwDON47POZKvYwmdVe1GaMl8y2w8ufWUI0aqsQXY/5RVZYkVVax2tkx73P81lX5Ttjs7RbjJu1vsz1xkWOw5ybRu8qqtgXMj7qWX9UsTeNaNRirx7RfCQrW3aN6N1uNq8i0fOgUbupVqCV3EIIIYQQQgghhBBCCCHaFr3kFkIIIYQQQgghhBBCCNG2VLYriZa4Z5IvSjYox6WkyyXb2S7XlG5RRr7//vubWa3kgFJ2HmeZKTXYa6+9OpWfdhD+d7Na6aAv0Wc5M3lNJlVxuS3lnsyD0kdKGF1mTtk7bStYty7zNVtnw2C2TrrAsi1YsKBMRzJl5pdZlFAWwvpvdysHs2r3wPZj7Ls8j5KwRYsWlWlKxBnjI0eONLPaNn3yySfLNPsG7UgY47ymM2fOnDIdyXWZN9ufaV6befBe/L7Zp5YuXdqpPGa19eX3wnulLJ/yZt4f68D7Lss2ffr0MA/2wUiWH8npOtKOu1w3SyZp9TE+k+rR7oNj06GHHmpmZvPmzSuPcbxlPLMvRJYLvHYW58zPnz+0+yGR7ZBZbRx7zLBPML6y8vvnOG6yDzLmOeZyvPdnHseTqVOnlmneN6/jzxGWs4rtGNkS4rzR/uxpxh7bl2Md5yM+lnN8ZCzwGmwzzqd8vsFYoG1HNmZ5fhy/GQs8zvuK7E8oL2d8sp+wD3r/4ueyvk1oEeB2RbT1YT2TqK+xPjMpKdkS4jqjyvzL6yWzDKHVGi1kfEznmMLYyubIjAc/h+XktXmc8++ozaI+0LEcUd/O7EqYZt34ddinMssGlp/fiXzex/GB3wE4z2E5vH6zfrQlx3IV64LsOeZtxhhhvUZzBLN17URrtaxP1ZOqZ9+bM/y+snhqRNaejYOZZYunObZXiTn2Cf8ey7G7ipWhH2fZsv6V0Q7WavXarJ5tXJXPZTESvZPhfD2zt2nEYqVKjEdk1jWNfOfK4iWyceAzi2XOysG68zGdz8isT0XtXcWm4pVElfcu0Tyuir2Oj2Wcx2Y2Y1Xif1MTPQ+iv5vFcVfPyqtjHlE9cm5Whc05DmsltxBCCCGEEEIIIYQQQoi2RS+5hRBCCCGEEEIIIYQQQrQtle1K6kmfuCzeJdFmtRJb7q7skgHaW1Dqceedd5bpAw44oNN1KPmjxJUSEcrNKPt7+umnzSyW3JvVyhgoVfey8hqUblHyS8khiaS0rIN99tmnTFOe/thjj5lZrWyMViQu/Tczmzx5cpl+5JFHyrRL8ykrfeKJJ8o0JQiUXPv90lJj2rRp0e2lMZFJhLsSkcQxk6VQjkciWwWPN7NamR7biTHgkhDa1TAWHn744TLtcWFmdsghh5Rpj33KdChx9520zcz22GOPMu3xRSlxJvlivLP8nveyZcvKY5QYM8bZdydNmmRm+e7tr3rVq8o05b20afDjjFXWM8vPPLz/u0TerFZemUmYmt0BvqtRRRKZ2Sl43WU2BqxHjpcueT/wwAPLY1n/ePTRR8s0nwce0xy7eG0e57X5XHI41vNZxGcH297zyyTx2TPR+wXLRvk/+2xmCfTggw+aWe0YwXEms9RyaTbvP7KnMGvPOI6oJy3M+jPbwWOAMcRzGfu03/GYojUS251jMmOfcyVvn8x+LWsnPz/rG1Usa7xM2TX4HGTa+w/zpS0A45NjCZ9LU6ZMMbPavsN5VWaJ0Yg0sl5stFsfqCJ1zo77czeyyjGrjU/aJPmch2MKy8HxhfOEqNyZRVgmG47Kxr7ImMvk6t4nsrGbsG4iyW4mzee1WQduxcPrsb5YTsb7K4VINl1lrpLJrb2es79zXI0sodgeTLPN+MyuZ3mXlZl5+DmM8UxyznQU4yT7TsbjXo5Gv8ux/D6mM36zOUd0nUbH8UYtTboS9WKEsB1Yb17PmeVG9u4hyiuaO3bMu15sVLEdcbK2y65R7/ld5XtbNP6TLN5ZB1437HMcS+rZAFaxv9iSbEyatfuod9+sf7ZPZl0XlSezR6oXr1l+EZl9TSOxXeV6jdhMZ3nwHK9fvufiuNBIfyQb065HK7mFEEIIIYQQQgghhBBCtC16yS2EEEIIIYQQQgghhBCibamsf6u3IyyXr1PeTYkVpQR+fiZrf+CBB8LjLl3df//9y2OUVNLOg/LJqEz1ZDsd83YZZJUdXHlfrAOXKfBcylqy3eU9D1qUUKpOmS+tV2jlMGPGDDOrleLTUmK//fYr06w7l1+z/SjtzyTXze6qvLloVipC6Qnjz9uYdck+QysZxrj3jVGjRpXH2DaMfcYt697bhFIS2h1wl3jet1tHMF+2e+/evcMyR3I45styzJkzJ0z7OcOGDSuPjR49ukwz5mjZMnfu3DLtsZpJjih9p1WK11e0u3vHPKL+3E5EksBMDp5JvKM8OFbOnDmzTDNvWuO4Hcm+++5bHhs7dmyZZvvQroRjml+TfYw2DJlk0NNsy8zKIbMjiXarznao5o71fpwWJQMHDizTHEPYP5YuXVqmve/RNiCTVbLf+DOFz4stnUZ2Fc/wOGL789lJ+TXHWe8zAwYMKI+5bZhZrcWMWxeY1T77PZ3Z2NSTvNPig+1exZIqukZ2Pfa16FyO6+yvfBZxLuFplj9rK1qleJn4DKsik6xyvCtRT0pb5dyoH2R2Txwr2Sb+rKUlD+12eG5mJeJ5Z/2yXrxn30V4DR4n0fcZzh+yeb3HO2OPYy1jnGM369TrhvnyOZu11fqOrY92iGvSrKy9nvQ6s29gfEYWJHxmc77DGGDMMe9642O949nYV8XWwdPZ2J3J572++Hf2bd43ny3RfXM8zspZz5aIbCn2DfXuo8rYTaIY53jGcZwx6edw/hnZoJjVtllkXZLNJxp5TmXxmZ0f2TFkVhD1Pse4Zj/n2B3ZYTDuSfZuIXr2tPOcpCr17qUR+5ZsLM/mD35ttnFmv5PZ3dR7t1VvXG9FWzZqcRPVF2M7Gy+IH8/ivB5V7rvVcd5+b2iEEEIIIYQQQgghhBBCiP9Hy3Yy4S8C/KWXx7nyyX9N5Ao3roDiRnv8hcFXKfMXSK6SyFaR77XXXmXaN8zjL8ssB1dScGWRr4DmNfjLPz/HlYBcQeK/kvKXkOnTp5fp+fPnl+mRI0eW6TFjxnQqM1ca8HNcyc1Vv15urpziJlksE3+xjDYf4cZ+2Sr4dqbeKs1o4z2zeOV0tjKKK1V5HY9P/qLFlbGMP8YIyzR79mwzq10dyJXXjE+uJvTVGPyVk30jW73CX/99dWm2UQmvzRWOp5xyipnV9iPmwb7IcnDjQm8LbsgZrWjpWCaPYa5W4P2x3ZpdXdRVqLIaxMl+VWXdRZvsMh7YnlSR+BjCvz/00ENlmmMXN5u87777yrTHBGMjWp1iFm+iw8/NmjWrTDNO+NyKVoZnqwh5bW7S6mocrgpjn+bqbD5LuSmql2nBggXlMW60nMW53zf7T1ZfW8pmTvVWy1RZpRxt/sjxIXt2+spOtj/nOVz5yfhkDHh+nGtUaQ8fy/h3rrZjfix/tHFalfpiTFGx4fB5x3pkHTCG/bnD+sxWAEYbArIvVln5Tlq54mZjUW8jpCor2iLFU6beyTYe9jbjWMXnNtuX8JwoX8Lncr0VRywnvwPU2+Ap27SPn+PqbJ9rsAysAz4Dsw1rfX6WbSqZxZ/nV2VTt425qdOmJJqTkypKzEitkK0EZP14u7J9OS7xOwDHKI6x9TbOrNdOmWo4U6NF8ZwpBgifa9Ech3HNMZ8qSV7bnzlVNnON+nY7x2yjRM+dKt83GlFisb4ZDz4Hzd6hsK2jsbtj3tG1GxmLspXjrfj+lall6pWNcwo+C/z7YxbLjcylG23vLYVmVUn1Nl7viLdhNn+vtwEwz6kylkcqmWwczjYLrqdea0SpWEXNw3MyZVz0uXr9tFFFZSvQSm4hhBBCCCGEEEIIIYQQbYtecgshhBBCCCGEEEIIIYRoWyrbldQzMKccl7IWpmlB4rL1O++8szzGjYle85rXlGnKZ3yzRdowcHMubn5DqQ3TLqWlxIzX4JJ8LqP349kmI5S185xINjZ+/PjyGDfZHDp0aJlmHbjFA+W8lMbwczxOmxa/73322ac8RquX//3f/y3Tf//738v00Ucf3en+aJWSSXhoS9EORHINyi+iTSXNatuaFgRet7T7oNyRG0vyuG/uSCsZWo3werR/YFs/+eSTZlZrc5K1H+/F+xKl5Yxxym5ZpmiTHvYv1mP//v3LNOvAbRx4PcrrWf+UEVEO7TGe2T9MnDixTE+ZMqXTtaMNtcxyKXOzGzBsThrZ0If3yvrILJ2cSOprVmu54cd5rtvsdLweY5t5uNUTLZOYXyZ5ZL9xKC3mmEY7H8Z/tOEm4455UMrr8c9NXFlm2lkQjvH+HMzaMrMp8r7S6Iap7SyPrCctzCSCkYSRbZo9A3iOj0d8FtJajLYefLazTD53iTYuY9nMYmkk5z6McZ5L6T2JJO2ZZRWPe3/NNvvjs411zg1Y/Tj7Xya1jDYKjPpqx3tqx/G7KvU2uzKrb1fC9stkq/5c5jFa4jGuGSNRrGbPm0ziHdnpsByRRZVZvNldJtdlmnMQ/87DvDgmsM8zb84ROTeMzq1nsVJlM8JGrQ+6Eq3YiC+qi2xDsaxveJr1zbkiP8fvQNFztor1QnTfmZVPlTrwclS57+gZx2cIy8F5FI9zbuPzw2yDSfaZehu6ZdSzxGiX+Usj8V4vj2yOl40v3saMX7Y759hss8jeIbt2Ixv5VbFka2Q8y/qJlyPbpJ51wOOR5V9m3cD64nUaic92ieFG2JjPo3qWVdmzNWqfjudHz5EqG4ZGfbPZdm1248lss/hss+2orFn5643ljY5lrbCq0kpuIYQQQgghhBBCCCGEEG2LXnILIYQQQgghhBBCCCGEaFsq25UQXzbOpemUpGc7h1LCdPvtt5tZ7ZL2N7/5zWWaxykR9Pwo+aNshFIrl3R3xCVWXP5OSSytF2jb4FYos2bNKo9RbkxZ+NixY8s05cku3Xr88cfLYywHJemUvrv0kTIB3h9lt5n0wtuLsk3Kgw8//PAyfdVVV5VptzfhubR6YDko7aFcsx1kY41IKSjpymRVLtllP2H70raDsm9v9yFDhpTHMiucTFbi8nPKYdnubBtK1b0vZTvOU4bMtqbFgpcpk7TR2oj4dXh/tCLhWMJzGH8uIWW/5djEe33kkUfKtLcnLVYoz2N98L5JO8S4WWOyKx5ne0aScI47lKizDWk14nXOdmVb0WaGebu1DM+hDQOfHZnFQwTbO7NyYDki6VYmx4xkXNnO9RwXBg8eXKYju56RI0eWx2jtwz67YMGCMu3PvkzmznJkcd6uZP0ysyWK7DB4LuuK8c4Y9nkKj7EdKXlnjNCeJ7IDyiw3ovhkmbN+ksnEo+OZHQTTPp/KLHvYh9nvOOdx+zrOtxiTLAfnb16PG2JL0tXH8IwqVg715jnZGJbZ0vnzl32AbZPZGrKtfY6c2eZkdjNOFWuWTErr+THfKmNFNCcnvB6fgewHPmZnFlWZ3YrfF6/BGK8ns+Y57RbrGyJzjqw2Gdesb9antxnH7szKh2m2u8/FszlJPeuSzEolk47Xk3pnfSazMnQyWwfOwxnjbtvIMrNuo/6c0ah8vR1iOyrjhsj0oxivYsnj36n43YpzSsY15z5Rftk1SDYHjT5X5VkW9ZPsXF47Kkc0jzKrjWvWk/ftrJ7rWcqRKu8houPtEOtVadZiq4pNjscx4zl6V2WWx3lkeZJdr5F72ZjWYlFZs/rKLLf8+2P2biezm3MarZdWxLRWcgshhBBCCCGEEEIIIYRoW/SSWwghhBBCCCGEEEIIIUTb0pRdicPl/pQhcik7bUUoCZg4caKZmb3zne8Mz128eHGZ5nJ5l/VxeTsl3Zm0hEvnXR7LclISTGsISrf69etX8/mOn6MUf968eWWaUmW3aaHsjfJ0ynh5316/vA+WjRIypilx7JiXWa1szOXBZrVt8dBDD5mZ2SmnnFIeo5yH90K7B7ZhI5K0rkAkYcokFZS0sE1cGj5q1KjyGOuYuBUOr80yMA/CtqZ1jjN79uwyzf5HqRVjyu0feB9sU16P57DfeczxGOUvlCQzLrx+WR5KmjN5PaW+bvFCiyD2c/ZR5s17cdiukbzZrLb/tKNcLJI2RjYNZrX1yDrwmGAbu8WRWa1kl+OKX4fjNC2t2FaMGZ5/2GGHmZnZLbfcEp6bSWs9/vl3tmUVKwfPL5M4Z5Iuh3XIMZl9k/YMtLbweuIxxittKTgG+P02Imt+JVFP6so5D8cPthmtkhzOAfhMZn8gjI3hw4ebWa1FGMna0uOS12M8ZTJDPncyOaPDeGcMex/k3ID9IZNDMm7d/oIyeJbNZfAdjzuZ/QT7cyMy4q5GPRuGjHrS0CwWsvl+NLfjWMQYyOTcPp7RWqmK1Dy6b47pWX8mHsNVLE8iyxCOA7w/Pvf4TKINl8c2+wbPjayKzOK2yu6btFoCvSmJ7FVaIXPmGMZxiTEcjYMcRzi+sx3YPh4nHOPqWTaY1bcmyOqAZfbrVHnW87j3efZh3l82h+Z1fB7HWGaMs/6j+XSVctbr5115PG91n4xilfcfWdCYrWvr7FyO/4wBHvd5UPa9qBXWDdkzO7JKyc4l3kczO8ZsbOZx/75Dy1p+387sT+rZRW4pljytop69WkY0f8jG/cxWKbKkYj+oYl1S7+/ZHDSamzR6vchWNDuX8c/3P/5+ku9I+T60nj1oled1q8dyreQWQgghhBBCCCGEEEII0bboJbcQQgghhBBCCCGEEEKItqWyXUkkbaVEi3A5fJ8+fcLjLtOl1QhlXJT/8touAeESeu4ATNsHykl4bS8/pVa0eqAMdvz48WXa5Vi0F5k5c2aZ3nvvvcv0woULyzRtO9xageWkhcWAAQPCz7mcgnYmlC9zZ3uWnzIMvzbvm7IEHt93333L9GOPPWZmtXXIz/E488jqvx1waQclKplNBeXUkZzJbW7MamV/PJfSJo9nxjj/zhjJpH7e72hzwn7C9ovuizFEeC6l75ms2aEMh/fFuvPPUaLCMrMOeG3KSn33X5bBj5nVWgmwH7uEmOWh5QnJdhDeUmRjkYWMWa1EPWojjvWM+cz6w602KH1iHmwLPmsi2xH2iblz54bljOwLMksR9o/MqibKg+eynJk9g0OZPmW/7Fesf+9DtCshrEc+a7yt9tprr/IYnzMcQ6LdxDseb1eqyPmjcZHjN+cuPM5nsY/x06dPL4/xGUm7EtobsE28TLTpYYww5iKpX3avWT8nft+Z9UdmNeJjMudxfN5lMccx2WXuzJdjOceEyA6Gz5nIjsqs/m7vXXlMb7YfZtZx3sczWxmOP5GFDG3DWPfsG4TSbpe88lnOMSezjPLjLHMjNidMZ5JZfo5yXD9Oiy3GO8/lOM26c5n7jBkzymOce7HuWD7vM+x/2XNoSxivSZWxu57tVDb2sf0Yww5jge2YWUBwvIq+A2UWCjw3shrJZO317Eiq2DdEsZ9Zq3GezfGd8+yhQ4eamdnkyZPLYxwrli9fXqYjy4jse1dGI5YY7UAVO5roe0gW4xzTIwtC5ss5B+HzNHq2Vvm+Xy8+q1gb1Ds/e9ZF1+azh3XE+Rq/j/Kdkb/DcWtXs9rvEJmNgx+v0rcbtTHpyjRi4VFvTMue1RzLI0sd5sVzmQdjO7L+rUJ0L1n/qGKZ5vFRpR9E9cw4zyy5+Dk+5/x7dvYuJit/ZD3V6DjdbMxrJbcQQgghhBBCCCGEEEKItmWDNp7kL1xcccBfTbiqiasVopWy/IWYqxmiTQy46ilbYcxf8Vg+/9WTv0D4hnVmZk8++WSZ5gaS/iser8FfOrlag7+sPv7442XaV9VxUzXfWMos/zXRy5qtsM1WumSbCEWf42oo1sf9999vZrW/pHAlIFezs114LzSn76pEv5Yxbpjm/fCeuamQ1y3rmOcyzqJfCfnLGleEZJt6cEWE/8rMVUZc7crVs4xbL0f262m20jPaoIF9n32bscW8fQUuy8xfUtlfs5UJXn7GHs9lmaigmDRpkpnVjmksx5w5c8o0V/g0++tuVybbkI5tGP3anf0qzHMZx16/fC7w2lw1wXGWbegxwRhg+3AVKIlWULCfMo9skyQvKz/HOuB4yhUzHCMc9sfsucXjnh+vwY2dWEdU//h1uHKBdcf2yZQ47bDKNSJbjZyt5OP9e7szLjh3YbszXvx5z+c2Pxet2O6Yh/eDKhtw1dvUqMrGNhHZBrS8l0ilxP7HZyafOezbLKvnzVWBXAHIZymfE17XkeLDLN480GzLUCtU2VguUyL5cbYpxxfWIevbn+ds02xOTtg+HhtZrJJ698i+GK287phHdL2sn7DuvJ6efvrp8HO8HvtMtFqYYzDLzGcu8/ZrZpvEZv253cbsRshWKdfblJZ1nK2G9+9oVORUUcKwL/nqt0bH8UiRkymOsu+/EVl/iO6LfZvjOPPInoHR5pX8exb7jaxg39JWb1eh3njGY4zrTP3h6ksqIRkL/Fw2xvo8qcqmeI2MRVVW9vrxKhtIR6tWs2cFx/fs+4l/NlPecBzgd4hoVe4rJX6r0sgGi9F3JLPa7/WucuXcm2RqrCqbTNb7e73NdBvdGLge0diZPRuz90oc+71+WUeZepn179/TGx2/WzFf0UpuIYQQQgghhBBCCCGEEG2LXnILIYQQQgghhBBCCCGEaFsq25VwCblLVaJN48xqbRF4fPbs2WXaNzbkBmWUW1PCxCXrbsOQbeRVb6Mxs3VS7khWbGZ2/PHHl2lKUlyKcuutt1oEZW/cZINyWy8/NwTM7BsoPXPJAKUDrFtKKWizwPK7NJL1EkmNzGol7i7viOrQLN5s0ay2PtrBroS4jIN1klkisH25Uanff2SnYVbbfjzHY5j1R1iXtPCgrMo3NKOkm5t6UVZVz5KHFgaZTI1E5Wc8sa8xdlyKzvjNLJEymY3XKSU0vAbzY9143pnEL9voMpL+tJMk2O832+iXccmNVh599NFOeWXjMOs8knzRrinbxILPDo49PjZl0uHMusTjIysz44fxxbaPjjHNa0eyXvZ/9olIEm9WO/Z7mTLLIOYRbR7EPkFLAtYR4yDaoLNdNhOut9lLtmkkn52RlQPrO7JD4zW5KTX/zv7AjW7rbdqSWUBEG1BlEvwqkslIKp+VI4pbXptxxjGefZB17rHIOs/6Jc/x/OptGNvx2q2WSW4OMtlnZpfBDYQ8hqtsWEa8HbjZbbZBHMedaLNqtgctUTKbtEjunZ1LouOZfU82T/NzeIx1l9UXv+f4nC2zE8g2l4vGtGxeX0Ui3FWJ+mS9DbnM4g0MO57fMV+zfLNanwOz72S2P9lzOJKLZxL4yEYpm5+QKhY/9f4ezW953ySb+zDG3SIt66P8XLSJYTYPyb5/tBuNWFVksR9ZBmT5cjzjOfyeV+8a9Z6t2buJzB4sKk/2PK5nuZaN49n44OXP+lT2vTOKccZ1Zv9Qb/NE0uhGnO1Mo5Yt2VjnZPNen4MwLnkux6DMciza+JlkVkLR985G268Rq5TIMoffsaN5c8c054i+WTg/l82tu4otj1ZyCyGEEEIIIYQQQgghhGhb9JJbCCGEEEIIIYQQQgghRNtS2a6Ey9dd0pXJOLJd75988skyvd9++5lZrVRp+fLlZZpWB7RFcUkNr5fJaDJbAZcp0KaBkvq99tqrTLvtg9k6eQoltVOmTCnThxxySJlesGBBmd5///3L9JgxY8zMbOzYseUx3gvri3m75Jz1QhkBJXk8J5LXsS7YhpQxUGrpMvmFCxeWx2jRkckjWP8uTejK0hqWzeuKEhXKuWgxw7pnm3jsUMabWfJQJuntQLkHr8E82Nasb98FntdmmvHJePZ7ZPvSJoHnZtYlLkun9J/ydMY783bLhsw6g/XBMYYWC5FdCeHxGTNmlGnvE5E9EfPtmAfLSil2VyaS/mXyKsYl659jhbcz25vxyj7EOPDxnmMv65BtzP7G54jHMcvMNuEzJZLsZ3ZE2U7rkcSKsUFrEB6P6oDxxfqiHRGfP5G8k+NN9hx55JFHOpWJ7ZfJ3FkfbJeuPIY7UTtlcjyOY2wz3rP3A7YTP8cxgWNnJAHmuYyB7Fnj9Z1ZXWV5e7vWk052PB5ZeGSyxmyO5X2aYwLjOrPT4bPN7zuyjTGrtefiOV5PvD+OMXxe1OvbXTnWo/JmElCOiawXxpzXRWZ1wbbmnNzrPus7kRWbWe0Y5GXieJxJzaMYZ74sP+OF6cjmo56U3iy2waANT2anwGcZ+7bfbyaRZqwyP6+DTCq8JRJZJ0V/N8vrIrIy4uc4HnPu4/NUzrc5d+X4whjmuBSVqZ71AstXz86k43Hi57CPsj+zL7JMHuPsw5l1CesusiBhXbD+o/m7WWzTlc3rM9rhe2craMQ2ke3HdvJ24DObMc7YZ4xHc8PMqiPro41YF2T35WNhZrfI+VMU45n1DmHdscwe45y/ZHDM9nqK3j10PJ7RLjFeb36VWfFkY109y0m2FdM+DmVWltn3hKisUb9bH/XG8ujcjudHZcvOZSz5GJ+9r8vqnGOyf9/MvndmdtGRvWZGPRvFRuNcK7mFEEIIIYQQQgghhBBCtC16yS2EEEIIIYQQQgghhBCibWnKrsRtDyifo1SJ0ry//vWvZZqSAZdsU8LKPLgcnhJ3X/ZeRRqZ7art0gTfEdfM7KmnnirTQ4cOLdP9+vXr9DnaPlDCM2LEiLBMlC8PHDjQzMwWL15cHqNtAqWPlFB4OVgvlAOQTObp7cXPZTvQ8py9997bzMzuvffe8th73vOeMr3HHnuU6cxuYP78+WFZuxKR1QplS6xXtuljjz1WphmXHieMEV6DVgoDBgwo094nMkklj7NMjB2PZ0oEGU+LFi0Ky+xpxhn/znJQohlJnFlHtEdhTFKm69YMHBMySQ7lcJHkmvFLWVJmxeHHWS/77rtvmc4sDVj/XV0u5kSyJN4Hx2m2Mfsw69fPp4Sd+THOKYX0OmXskKVLl5ZpjjEcV9y6JJP6Mo55L5FtRybjz8ZLrzvm67tPm9XaXvFZVE/SGUnMWGazdXXGuGS7ZffldUMZKp+vfH5WkZx2VaLy8n4412B9Z9YYXoeMPcYCY5zjSmQrEO303rFMzNvHRfYT3guJZIvZGMXrZW3t4292bvTsYDqTgWb9NXveRmXms4PzN792JnOvMm9qFwmw4+WsJxftSBQvjF/WN2M8ilVeg/lyDOYYxTx8jpTNeTh+RpJl3ivjjGTjlvcr2ubwXtjvou8rPJflyOS/jHG/L47XkdUDz+U5VeYfmTw5siprB6o8izLLGr9nxkgkZTernb96W3O8zr6Dsv049/dy1LMG6VjmSOpdxbokgmXO+lQ0DrI+MysHwvvyumE/4ZyJdd7I84tkddMO43ez86l6ljyEMc65H+3wPMbZ1lmMZ7Y3DvPIrP8i+6hsnpTZQkR1l9l9ZO86PC75PMrssRi3PO7jNK+RWaTVG7+yPtzOMZ6RxX6V+/M2yqw1WOeRHV/23Ga7Rs/Ljuc7mYVy1B835DtUPauULHajOSLh56LvuWbr6pTnZs+URqgSw83GuVZyCyGEEEIIIYQQQgghhGhb9JJbCCGEEEIIIYQQQgghRNtS2a6Ey+Fdck5pM6XZlLVMmjSpTB9zzDFlesiQIWZWK91ifszj6aefLtO+ND6zaaCckDJJ5udSA7cOMau1D6G1wp577lmmDzroIDOrleVTnk7LiUMPPbRMU6LvO6vzGplly6hRozqVmdIuSgqYJqwb3/mX8oJMlkAJ0oEHHmhmtffKctCmZcqUKWWaUj2/ZmZN0BWIJBg8RjkrpRNz584t07S6OeCAA8wst8hgu5NIrpRZZLAd2H/69u1rZrmkm/2VEvzp06eb2bo4NauVoDBeeD23tDFbJwNi/2M/d4sJM7PBgweXaY+5rD9E8jaz2jqNysw45FjiY5CZ2f33329mteME5U60KOIO2pFMqKvLfyPZW2bHwvjPpLw+vrGemR/rvF7dsC15PUosWSbvb4w1xi77I4+79J5yTUq7Mvk4Y97P4bjI/sExftiwYZ3OycbCKmOyH89sGJgeNGhQmfbncSaTzizD2s2uhHjZM9sl1hVjlW3t52fWLhxnI3l8ZO9jVhtnHHsopfS+xvLTOoIxwJjydmU5SWbnxthwyS3PZZl5X5wX+bMys5UiLH80hvDvhM8t9mMvK6/N8rMeM+uAdpAAR2XMrFjYvpns3NOcM1JyndnU1LPBy+KaeFvynmgnyDkB+4+3K++P8c7PZbHoccYYyaxXWB9+HZaZ8cRnGfNjP/E6Y4zzc6znmTNndioz882s3bJ0O8Q48bJn88DMxiCyYmOM8LnKeKkX42yzyErGrLZ9/JpsX1qykeg6mRVOZucRzeNYnix2+Ll6c2v2RcJ5hN83+xy/S/F6rI/IMqqKvUE7E91fFRsitp/HH/Nie2RjqcdRNk5mNpLsB54f25TzbhLFcDaGZbYj2djmMFajOjJbN6Zn/TazheBcxeuX3xmzZ2c058iei1s6VSyYGrEr4XjFuTDj3Os6++6U2ZhEY3w2T8isPzx+svE7u+/oO0WVd2lZX4mORf3YrLYv+3F+N+ffCcsX2aJmZGNLsxaCWskthBBCCCGEEEIIIYQQom3RS24hhBBCCCGEEEIIIYQQbUtluxIuX3d5ev/+/ctje+yxR5mePXt2p3PNzA4++OAyHe3anElVouX82ZJ2SmNYZsoYXEbCsjG9cOHCMs3y+TlnnHFGeYx2LFy2T1kzZT4PPvigmdXaoFAeSgku5TAug6EEhhJ+Sm0oD3XbCrN1liyZ5CuTJo0ePdrM1lk6mJnNmjWrTLuNi1lt/VPu7vdCOXhXg/fv5d13333LY6xvSo7Y7rTAcEletvM6JSisKyeTSWVyY5bJ44SyQMbOxIkTyzRteyJJCCWEvF4Wq97vKJH7+9//XqZplcIxxOuXeWWSomwXa6/HTIqfWfy4rcW8efPKY+wD7KOsg2j34naSnrmUdOTIkeUxWjcR3jfHe7eqYX1FcsCORLu50womI3o2MNb222+/Mv3www+Xabahx0k2bnL8JuyHHj/Dhw8vj3FcZN60DHGrCT6fWHeZTI3l97hjn2cdMC7Zx9xCJbMi4fVoScN28b6V2Uh0BSI7Eo5zmUw8mueYrbMYyKSFmTWBtxntTDKbqsxyIapn2nNwvsL28zJxLGRb816ZZj/w82nXxLpjfrSs8jxYL5TyspwcV1jnnkdmF8Fxg237xBNPdPpcNudhjPP53k5juNm6uQKtMGhPltnlMO3jFccRzkEyybvHONuX8c5xMJOPexxlc6WM6BlCsjyiazMWaKnGOmI/iGwkOB4zzTzYt6PnF/NjPXMu71Yu2Rwxmx+xjfz8rjyOR3BMYZ/NrEQ4vnhbc16cScsZU95mjF9aRmX2PNG4wzJntlmN2Doyj+hemWacVZHMe93w/jIrB/Y1fjfw/Bj3zIN9g+3i9cRrZ1L1ejHeLkTtVGW+wHrxe45skTrmzfrhcYfzPrZD9v3K2zUrc5ZHVJ7seD1Lrsx6J4udyEaNeZBsruL3nY0fmUWsz8MzCybWbfZ9oVkbh01NZEfCes7mXPXsEdkOzI/1xeeBj0O8Xma/wXkM2yWK3exdZmSpVcXyMYtBJ5tXsZ4j65XMXiR7jnAs9/qoYo2W2cbVg/2jFXMTreQWQgghhBBCCCGEEEII0bZUXsnNX0v8bT5X7XI1Gzcf5MrW7FcuJ/sFlp+Lfo3gyiL+6sBfG/jrv+fNX2n4iyU/x3N81TZ/IeKKKubBTcemTp1apv3XF/4Kw3tlHszb4Wo+rghgHWUbKPk5/LWFv7BwdWa0KpgrcBctWtSpbGb5xnF+vF1Wcnu9cFM5rqRhm2a/HHus8lc2xjhXfHBlpdc3f93lJkyMcebHePH2YzyxrbNfi33zVMZWv379yjT7LVddcNWYH+fqCtYdy8RV214+9g3eN/tdtgrFxySey7rl6gaWw8cvjiW8J26QyT5VZcVFVyPaZIP3zXblijb+qhqtFOOvxoxLEm2wxfbLVj9kGzt5W2QrwHk99ptoozP+ne2arTDw2OUzjjGTlSO6b8LxJFvx6vWRrbDiqkTWqW+AOWHChPIY+0e0UnF912kHvJ45pnGVAWM1Wx3l40a2yWO2qaUfZzty3OcYyfEjOs65COMzW5nhx7PV/lm/izbl45jA8hMe9+tkKzGyjaRYJr/fbDUW5xpsT49xqiqYL8/NNsNsB6KVmtkK6myVEWOOq+EjmF8032RsVdmIj3l4H8w22o3GfF47G6uyTXWjlUhZjLAvsi/5+B9tjtSRbOWT1z/nQWzDaEUs72XOnDnhNTjeRKs8O6a7KtHqxEwNQLINtaigjfJgW0fjTjYWZSsS2Q6+8juLcbZHtMkfr8fPZX0+WmGXPd8i1bDZujE2i/F639PN1s3JM9Uc3y1EG7nzO2O2WpXjW1df0VqFaNO29aVZ995m2Ubk0SZ2ZvEm9NkKasZn9L05W8lKotjJ5gVZm0Z5Z2Ntpnr274rR85R/N8vVkv6ug9djG3KuEs3/qsxB223evT6iZ24V5V2khM/aLVN++zWzzSaz+Wik2MrmFFmZ6pWTZBsKR/006yuRoimLqUjVZ1Ybrz4WZ9fLNrr0a2fjdKQG6Xi8WbSSWwghhBBCCCGEEEIIIUTbopfcQgghhBBCCCGEEEIIIdqWpjae9GX0lHZxaT0l7q973evKNM/3pf2UkHCZOuUklBm7DJvL3nntzDSe8kNfAs9l+JnUj+e4dcJdd91VHsusHLgUn5YXLj9kmSnXmj9/fpmmtYLnR2l5Jreg9IzXdmlFJsXJJGR+TVo2TJ48Obw2pYTRffP+ujIuu8jkHNycj1KkaJOsTDLLGI+k3ozfbGMBxlyUH+Oa59J6hrYy/jnaFriFiVlt/EV9w2ydPJGWCQceeGCZ/tOf/lSmae/g5cs2AMsk0NHmKNm5HB84Jnh8ss/NmDGjTLO/snysu1ZIazYFkUyI7Urbm8y2JtqkkXWRxWi0wWK2UVNmGxXZonDMY1uNGjWqTLN83vbsB4xXlpOxxH4fWQTQumT8+PFlmuOej62s83pyfLNa6Vkk06ZULJNYetvyvlk2Xpt1wDZqN9mkxxTHCcrWI6sls9jKhrGQSRyjsZxy1Gyjs6zu3UqIYzLblJua0r7A44X3xA1VaTfHuQTl4V4+xh7HSG56SQszvw7vI+ujrA/Wr8dwZtPCPsPjXl98LjMPjjGsg6jfdWVLh+h5wzjjfDub6/IZ6PWZzQkzCyc/J7LyWB+RXRvHYF4vsxfyZwuvlz1DGONZTDl8BvK5EfXXKhtbZ1L5jnmZ1cY175v4mJVt0sQyZbHfDhuWRW0ZbfZtlkvLI2sdjg3ZJvLRJqjZpu+Z7D6a72SbmLOdIgvETM7P/LK5lKdZL/yOx3EwsquIxuWO8F4i69Eq9my8ts89sw1FSSMbdbYDVawb6lmXZPZR2bw6ss7ILGYz+zqP8cy6hu0XWSplscV7yTYDdLLN77KNr72sHAcyGyQej2KcedSzmzBb118za0bSbvPu9RE9O0m2IW9k85FtwJiNyV7XrE/WP+M82xA0itPMyoqfizarzvpxPZuxelYqHcvk9ZRZDZFs/uDxz3rOrOJ43MuafZ+tZ++yIWgltxBCCCGEEEIIIYQQQoi2RS+5hRBCCCGEEEIIIYQQQrQtle1KKBfxpfpDhw4tj1EmRYk75bHRDuJcvs50toutS1Ej+xSzWglJJGskkX2KWa2MnEvnXeJMqRvJ7qVv375l2qXALA/riMv9Kf91+STrljI1ygcoNaB0IbJyYJtQxhDtwDx69Ojy2B133FGmWaYxY8aU6VmzZpXpriz7dSKpKS1aMgnHyJEjy/Qee+xRpr3eGGes1yyOaPPSMS+zWmkhY5z17TIVWo1kO8BTru/tRLluFhfs85Swe34sG+uFdTphwoROebC/u01Ax/zYLyNrEubBes52g3fpJstJuxLWHe0B2EfbEa8DWpEwzllfjAPaXXjdZGNyZknlFjeZdQHHFVpNzJ07t0z7c2Ls2LHlMUoUKTVnv/Fxj88ZjqeMGeZBaZzHOetr//33L9MzZ84s0/7cYjnYl9h/Mose9rdIZsaxnHmwzNGYk8nGaHNBC4t2GMuJ1xXl2SSToLJPeFtldlNsDx536zb2HcL2ZbzTFscZNmxYmAeJ7CdYZtpBcF7CGGf5vXz8O5937Iv1+nZmvcA64LV9PGH78FmaSZ99vOG5fF4Tjk1ZG7UDPgZxDGadME045vlnszrmOMf5q49LmZUDz2WMRDZpHJ8YL8uXLy/THKP8nCjuzWpjJ7Nc8vjjMfYNzseiWM2eIYT9gM8yr3PeU2aDwbb1++KYxnyzZ2o72zp4HWV2T7w3nsM68rbK7GP4OR73vDM7kGw+ynHHr8lr8BnL752RxD3r29m4yjy8TLwnzvHZv6Jxk3ll1mqZJVdkD5B9F2Ff8rx5f5zDZVYHXdl+pypVLAwYwxw3vW4zCzW2H+PB25pxmFlUET4jPO/ou5VZra1lZLfKe83G8cxKxPNgearMF/xzvD/GHD+X2cV6OrtGZjXo6cy+Yksiuq8stkk0t+P5mQVG9m7L4dic2ZIwv8iWI5ub8vsSz/H7zWzUSD3LnGxuwD7NMvszICpPx+PZO1BPZ999MmuZqO3r/X195zeCVnILIYQQQgghhBBCCCGEaFv0klsIIYQQQgghhBBCCCFE21LZroS2AS5topRlzpw5ZbpPnz5lmvYVlDM52W7JlCVEsptMjk3JXiav9KX2mayF5aSViNcBl+pTbsbPUW5MuY7LCjKJAonsM6pIUAnrycuRyRwyuanLHNiu/fv3L9Pz588v0wcddFCZZp1mOwZ3JSJLkGwnatZFds8ux6WUjDEXyU/N1rUDY5Z9jfYilH+xH7iNBuOTZaY8kbL8aOdexlBmbUBbFIdtTvkLbY5YDpcX8Vz2tYxIpkqpEsvJOuK9RHYllHNSJsl7zXZU7sqwnF6/meyX9gb77rtvmeaY5nGeyQszSZRLyDjuUCZIyZfbPpjF9lSMNVp/0HImkuZTxkb5PNuVNhIc9zwPxgbTjHOWz+uO/SpLc4yPLALYv9muHC/4Obfc4ueyvsI8ojbMnt1dDS8nx2a2O2OONgU8zhiNYP+JJNyZTdKCBQvKdGYZ4uXmuMh5CY/zHr39WJ5sN3XeH/P28nOc5r2wvhgvPhYwX5aDscO+xuP+WY41kTTarHZc8X5Hm4lMgsprZzZN7QRjNrOs4JhOvC44j83iJbIb4BiX2U5xfIyuzTaI5r9mtffoMcdrc76VyaH5LPC+ls15OFayf3l+/BzHa44xzCOyNuIYzHvJpMADBw40M7NJkyZZRJVYbpf5SkeyGMks6vj8i+aYHLdI9L0za2uOc+w/Ufxx3s/xh32GY5uneR/ZmJj11+h6/Bz7Q1QO3h/Lz3rkO4KoH/CeeC/ZOO7fV2j1lsV1FXuPdiW7t8xiwduBbZbZdrB+IjuSbEzk/CSaX3CMZrqenUHWptkzO7Lwq2fXalbb571M2bhS5d2FlyP7HsUys/84mRXnlhTLLL/38+hYRzIrC3+mRhZ3HT/H4z4eZe9i2N6ZxZK3Lccutn2Wn5cpa+PsmRLVTWZDltlJRc8Rljl7pxpdJ7PXyex6vD4yu68qRHVXBa3kFkIIIYQQQgghhBBCCNG26CW3EEIIIYQQQgghhBBCiLalsl1JJAujJGzKlCll+vDDDy/TlP9PnTq1TLtchFID5pdJYyK5GaV+mUSTchG/Nm0fKF/JZBOeN5fLZ7tOU6ITySZYnkzeRomR581rRLtnm9VKQnmO55dJjTL5kx/nPdHKYfLkyWF+mSyqqxLJSBmTjz76aJmmFcHgwYPLNGXdbm/CdmQdRxYlZrGVQmZTw/KNHDmy0zmMLbcqMKuNd9o7eL/ifVC+zvhk3oxLjw1aTGT3SusPLxPvO9sVPJO9eF3z3Ex+xBj3a1NWzD7KMmU7LmfjRleDce51wGNz584t04yvQYMGlWmOIW6Zs9dee5XHKNeKbDZ4DsfybDzlGE/bFH+msDxsB7YnxywvB6+dyTsz+bjnx37APFhmxpJLjSnD5TXYN6PdsZl3Js9jOVj/fj7rIooHs9wqqN3sSpwq98YYZ725/RjHUMZnZukUSQQzWSxjbsSIEWXax1HGKq2nsn7i5WAs8L4Zk5wzsB87jHGOhZSoc6yIZO60H2I5WSae7zHKcmZtGEk+eY3Muo5jRTvMUTL8/hiHvE/eGy1meNxt5/i8z8ZxjkXePplFEtuMbcmx2ctNm79MPhs9txk3LAfzYAxHz3naBfH++B2mnsSWfZR1l8nRvf8w38zWgWXysSKzMyGNyoK7Opn1UGQ/Z1Y7Drg9FOe8mdVUZN+YWSFk4zvHRM+Dc2teL7N18Dji+JnZLNGuk88Tz5vl5DV47Xo2mdmcPLM58nOycSCzxPDy81gmr2/nsTuCscz25T1znGOM+7wys6DJLH48NpgXLToZ79F7BbN1bc1xMLOdymwMI3ivtOjktb18mSUKx2P2Sz8ne6+TWfxE70gyWzoSWQKynrMYb3frkkbKn9nWsN3cBo1jV2b/FFnxMB4iC8+OeTAGvXzZ94jMtsbJrLWy+47Ixojs/ZDfYzZuZv0mmofVsznqWCbvp9F40/HcVsd5e7yVEUIIIYQQQgghhBBCCCEC9JJbCCGEEEIIIYQQQgghRNtS2a4k2qGaS8+5gzrtSjL5RiT/qrJc3vOgHI2yHEoXKD+MlsPz2Lx588o0JcuUZvk1o91xzWqlDZQ1Zzu0OpRocud71p3L0HjtzJqF0hdKLFxakdmjZPfidU4ZU7TzeMfPZVLRrgrr09O8H0qwxo4dW6YpfaXs2yWRzCOTR0c2DtnOvYwttgnL5+fMmjWrPHbvvfeW6SOOOCIsk5eZknTKJDMpDPGyUsrO9s8k0B6fLA/7Ns9lfuwz0d9Zz5QzRTI6Xo+y/ccff7xMUzacSZu6MpHcN7OTOeSQQ8p0//79y/TChQvLtMvbWbeZbVE0lkfPBbPa+qcskWOaX5sSYKZp4cP78nbmtSOrh47Ho/Iz/vjMyeSifFZG5zLOeS8sh/ch5pXJI3ltr0fGAMcvyp0zaSk/2w5ENmNsR1omMcYpi3Uy6WomY/W6ZxtwXO/Xr1+Z5pjFsdPH9SVLlpTHOO5wzsPPeRxlck7GHNNRHFHezxjn9VinkYSRczZSr89wPhM9o81qxwfvG6wjxjjHgWxuFsVMVyOau2VzQt4/452WfVG/ziznGEf1ZOIctwiP+zjH8YzjYBbjfBY4kaTZrDaGWX6PAZafVlIcE4jHEe8jq4PIYoXH2fc5v+B983M+hrB9OK9iHlsKHu/Zd0qOAWwzzos9Nth3WFdZ3lGMZPM9zs95jscU44L9j+WPxrPsuzLTWSzWs66i/RXx7wHMi+VgPyI838/h9dgm7M+RzWdmoxmN12a19et5dOVxnETvJpjmGEDbqcgmL6u3bGyI2jKzXeBYylj1OROvx3kU24E2UF7+7PtxZN9jVvtc8/6VWUhEcyOzdX07s0mInjEdz/fyZXMLlpkx7nAOlN13u9kDViEbQxl3tE9j3fhcMbMMrTd3yZ7PJHv+ensxD45pbDfGqJcpm3dmFm3R9y/2Y94ry8y8fY7Ba0TjdMf7isZZXi97NrD8fg4/F83BOl6bNDuGt8dbGSGEEEIIIYQQQgghhBAioPJK7mgzJf5Cx9VsXD3HX6r5K7P/qhBtHGBW+8sX3/L72//I9N8s38iDx/0XEP56zY0zs9WZfu1sIwz+ipf9iuq/JnL1NlebsExcxefXzjZx4K862bW9DnhPWT1Hq2x4rN4v/yyz2bp7zFYedwWiX/YYC1xZOnDgwE7ndsTjL/vVK/tFzX9l5i9h2WZ0bL/oOH/tv/XWW8s0V3ONGTOmTHtMZasEsg0TotUTjHGOFQceeGCZ5ipxv3b2S2RUR2bxyoQsxrONSHwFPlcoZJt+ZBsmRKsRuuKqV9aB1znrkHEerZI3q42DbAWxw/riuT5WZBs0MuazzaH813zGKNULbMMDDjigTHss8RrZxmmMu2jVOVVALAd/wY82YmP8ZWqP7Nnm12G+XAUbrVo3W7cCgWqlSZMmleloczOz+Jf2rrwRTrTijs8m1km2wofUW+WarTDzduD1ss1Vs829PC75d6785DyBKxg9/rgSOlvdws9FG1qNHz++PMa44PgQrS6PNig0q7YZXqT+ysYB9iWfV3Bz6Pvvv79MZ20crT7vyjEe1WG2eozpSEmSkan+GCP+7GRcZPNw1j3j0mOb/YRl5spqjl0OV3tlG9nxnGhzSiolsnuJNo3NNghn385WmPlnGdfZysFo423OQ6k2yxQPXXE+sj6i/pf1Q9Y9V0ZGz9BoJadZvnGZtxPnvCwb45btxxj386PVyma1KlCucvXYqLK6N1O3eMzx+zjjKYvxaAVepAzreG0+CyLFHlU90XdswrqYOXNm3TKTdpirRFSZX/O7UxQDjL1M5RFtQpk9pzl+Rkoens9Y4Pyez55IQcyyZepv9kHm7ddctGhReYzxWU+NyPvIxm4SzdEitY1ZbczxuMPxg/fdLurgZqmyUSfHRdav12n2vi7bENHH/kyVyZjKNoF2sud99o7Q45FtnG1yyutxvPT7zVb889rR95LsOZK9T2Q5vPyNuG4wzb6bOW1k43OzqpwtuwcJIYQQQgghhBBCCCGE2KLRS24hhBBCCCGEEEIIIYQQbUtlu5IHH3ywTJ944olmZjZ//vzyGC1KuPSfkhnKt33JPZfnZ9LWetLjTIK72267lWlKS1yeTBkE7RloszB48OAy7efzPrINZXg9nuOyS0qMs80deNzlPNmGU5nkIaqbaOMDs9xmwduFcozRo0eX6XHjxpVp2lJQuuDtlUlIugLRJm+UlYwaNapMZ7IY1q1LkCgR5N8p3aKMw+uQ9Z1ZlDDG2X+8fJT6sZy812gDA8pjMquUbJMNh3KULOYYqx7jHDP4OfbXbANJLwdtR5gfz+Vxj8+99967PEZLnunTp4f3FUktq2xosTmJ4pXxQAl09jnWjbcLJUeZBJX17/WYjVeMO1pKsN38OPtBtrETnxl+38w32+iVeUQbVjEeMrliJG3ONkfN6pFl8s9GMjaz2riLZGF83rGctJBiPUcbAnVl2S+fo94mrB/eP+Oa90/bHq+3bPOfTKroMRdJx81qx3KO1ZyDeKwyxinJZfsuWLCgTI8YMaJT2Ziut8kay89nBP+e9VdPs29nG1Rn44M/Hynl5fX4PIg2A2d9ZfYEpIqFSlcikrRnNl2MfY47jHGPs2zzWbZTZE3AsYhtzTzY7zgee5n4XOH3C9riMPbdZofXYLxkcGzzsmaxxWuzfv3ZwbyyTd0yabGPN4zVbKNmpv06HFdIvY34mObfu9qYXk/OnlmxMN65EaTHOK0ess0To/yyDXo5rrIfUPruacYL+yK/A3A+5nHGvNhOvF62oZnfL8do1hH7YlRm9gfGJ8fV7Du732Nm35bhfYP1nNn+tPPG2FH/5P1kdpcksvnLNhXO+rvHJduJf+dcOZsv+GdpL8Uxkd/L+K7A+yi/72X9i3O0aOO/zMKTn2Of8WuyPNlcIHs/4/fIvsP+ldkichxysvcw2aacEV1tHDdrbEPYevZJZutinmMX5x2ZRYy3IduP12AbRpamZuvimMeYB5/L7Ic+/2F/zdJ8NrCPRXYl2bjI/Pxesnk4ycYL/yzHgswak3Xuc0rm2+h8u1nrKa3kFkIIIYQQQgghhBBCCNG26CW3EEIIIYQQQgghhBBCiLalsl0JccnAnDlzymOHHnpoeG60A61ZvEyey9C5PD/avTPb0TOTwUYScMoVKFHI5EGRTDuTa2Uy3Xnz5plZrSyTMuVsV/poV1PCz2Uyrqj8zK/e7t2sW8pOBwwYUKYpdaUdh0sv2kVK5vdK2fjYsWPLdGQNYlbb1h5TmSUP4zaKz0x+SVkjZTGRXQ7blzYUbLMoHhhDvEa20zHb2tO0+qG8hdejNM5lwZl0iLIf9l3m5/XLdsh2E6d8zfsg+/OgQYPK9BNPPFGmKTeN6IpSsQyva0r1hg0bVqYzeSvjw9OUJ/Fc1mkkPcv+TliOyG6Fx2idRdgXWD4nky1nMn2XyR144IHh3zmuU3rmUjfGIsvD+GKa/d5jjFYVlIJyHGYfcnk/xzXeH89lnUZ1wHvtyjZUTrYTOuuVcR3tcs9YzeYB0e7mkSVUxzJxnKV9gefBGOnXr1+Z5jjLMvv4xr7Da7D8HNdZJo9xtz4xq71XxhzL53Jg9ueobB3LEdVNJumnTJWyUq9rzk8j+xqey+uRTLbZ1fCycy7CGM/scqL5SCaLjqTtvHbWB1ivlLNGkmRem2MR5wmMB3/O83O0F4meU2a1cwkvR+/evctjzI/z8GisjGTrZrX9IftuE/XtbCxl+f38yGLPrFYGX689u3Jck6h/st+znbJYjOYq7OOZ9aXXYWbRV8Vaza/DcvIamW1TZBfJWI2+c5jVzg08b46fTLMemYeP49ncqYoNhtcTP5dZNUZzH57LOmBfa0QG35XjPbIrYd2znbJ29/vL3rdk87YoxjNLB9p9sB38sxzns+c0+4mPV+wD/DuvwXuNrBeyazN2mIc/GzN7EZLFvo830dyp43GOzT7visaJjp9rl3cnVYjiPIt5Pu/Zhh4rkVWdWW6r5McZ+4xzXpvHOQfxz2b2KCwT0/69jHmxL3FezBhkLHmcZxZNjP+oHNl70exeiJcjmm+b5Xawnua5m8omTSu5hRBCCCGEEEIIIYQQQrQteskthBBCCCGEEEIIIYQQom1pyq7El/5zCT0lRfVkV2brJAPZTuLZUnbPg8dYjswmJCpHJgOilQPL7OdQ+pNZTmRyTZcZU47MZf0sR5Q37y+TTmc7M/vxzKKExyM7GdYzr822571QAu1yiq4sFYvg/TAWKCXJrDP8s5lspp7tBfPNdrNmjBM/J9t9njJzyls8BhYsWFAeo0SI12MMR32bksrMkoLSTS9f1jdIJpHxsmZxncW+9zvK0Vj+aBdss7w/tgseryx7FuesO9pheHtntgMcNyKpJOswi7VMYlnPBomSR8a55832ZhyxzJkMzevDLUDMzJYsWVKmGdsss9c1ZfXz588v06wPjq3MwyWWlCSzrVhmHveYZ160Grv11lvD67XbuB2R2QAwLmjrwPjzeMh29qbMkDHsZFY4jNtsrIvklYwttu/w4cPLtPcNjl3sw1n/iuT2PJbFeGTtRsse1hHPzcYel7Fnccg8WL8+7+D90YKJ9cHrRTL3RneA31x4vbAdGYeUx9IeIIpxwrbJnm/ePpHFnVltO3FOyGt73px3MnZ4L0OHDi3T3n60X+I4nj0XaBfo5zDGWTbOldhHvV+ynHPnzg3LwRhmHjwnujbbjfXv9mK0Gcss7bI5osP24XO7q+Exzhjh8ziz8Yus06LvlGa18cLrRPUSfUcyqx3nIquUzJKB90I7Ko8Xzsk5N504cWKZ5n3TksFhv+S9su6I1xftKdnXGJP1pPuRtURHON9xK55onGiUdpm/RN/R+YzlOMjjxOMls+GJrDHN4vl4ZqmX9Rlvq8wqi88Y2gq6fdTChQvLY7TXnDlzZlimKMaz7+OM8Sj+OI5ndpeEefs1s++52Xd2fx5m75E2laXDpiayK+H4x3tl/ETftxnb2bgSvYPh9bI4z2xMvI9F3wc7ln/fffct0z6GZ1ZDmSUrvyv6dz5+b+H3cX4u+t7Bz2X2udl3JR9zWEeZXS/Tfn4Wz1nMk2bjXyu5hRBCCCGEEEIIIYQQQrQteskthBBCCCGEEEIIIYQQom1pyq7El55T/jFt2rQyTflKJgnwpeeZxJHL5aMdhaOdtjvmx+XyPO7loISB+WXyCJdS8T4oT2F9ULJGWZGfwzriDumZ5NyvSTlXJqfj53idSI7KY9EuyTyHMssJEyaUacroKLOLJAjtJrlhW1OKevjhh5dpSkUoBYnuNbPIofzD641tSvlXZgEU7RbMPBjLjJ0onrN+S2sGSo8pq3V5L8vGc3kvmeylY3k6ksW7S/RZfu7OnOXnEkxaFU2ePDm8Xiavb0e8b7P9li5dWqYpp6WUim0YSZsYMxz/Itkk6zPbdZpjD6/t4yzH7Exiz7byGGVsRDuum9XGaLSrPMdK1gWl5KyP6BnGGGZ+mZTcnwPsM08++aRF0J7By8F6njdvXqd8zWrHi3a2LvFnIGOBMT5mzJgyzfuPdn7P6iSzI/FYjWw4zGqfF4yzyPons0Zi3jzutk+ZvJfXYOwwxl0amdl2cJ4TxXi0q33H40zzHr1vMsYjWyyz2vb0sZ9zEc5RIqsls9q5UrvFuMci24Bj2z777FOmOZZGc7QsltlO0bOabcPnRr3xn+XP+ldms+BtzPGa544YMaJMM8Y5/tF+wWFcc47MPupjM8doli2T9JNoPsLxP7KaMltnV8G/83sEbawyCxLPr11iPapDxghtZTI7O69vxllmURJZjbCuOHZHtglZHtlYymtzvjxkyBAzq21rpseOHVumM9ssH/8yiTvvK5LrZ9YL2TOLRDZcPJZZP/pziOMHxzT2k3rjeD3rma5G9iziuMpzovchjPEq9+x9hjGSfe/ktaPnAsvDvshxl/Hn1iT83kbbSsY4+wbnUrNmzTKz2nuN7DDNYjtCWkJwLM2+MzKm/JwqlrvEYz/rf43YmLTLOJ6RWebxvqLxOXsXyPiP6oZ5ZZam2bjhbZHZCzPuGOejRo0ys9r44r3uueeeZZrPcM7Jn3jiiZoymNWOoSwnx0jvy9l7pyzWIisX9gnWV/b+sp5VcjY/yixNor9naCW3EEIIIYQQQgghhBBCiLZFL7mFEEIIIYQQQgghhBBCtC1N2ZX48nUuMZ80aVKZpqydMsLIaoOyEUoDeJzL5SPZMCU1mcSdy/J9OT8/l8kVol1Ls513KaPJdjr3a7PuuMQ/s1jx+sisSHhuZsPixylLYJ1ncmgvH4/dfPPNYTlIJBdpF0mN1xtjaMaMGWX6pJNOKtORxJp5ZJYimUzb6yiTmGUyNOYX2VAQxiSl3N5faXHA2KIklNfu06dPp+Ps71mMRxL3zA6IZNIt/yz7X2a9Eu0yzjL/7W9/C89lmuXzfpDJUbsikRx1+vTpZZoSb0pyKbfyOGddZDs1sy08jtl+HJOZzvqQtzPLxtjm82fOnDmdzmGcM1/utE75JqVl3sd43ywHxw4e9/7BfpD1aY7fkR1RFpeZ3Yqnme9DDz3UKd+O5WhnKwePHcYeLVoWLVpUpikXZOxHFlKZZUg03vDvlOFmtlFsa3+msGzZ9VzKaGZ24IEHmllsU2a2TiJsVhuL7D9RH+UzheVg2j9HGxE+A/icZB1EVlbR/MmsNsYje4LMaiyzJ2hnCXBk/cS6Z4yz3RcuXFimPeZYx4xDxgjHM2+zbF6ZSV+jeTTbmucyxqdMmVKmPRZ5bVo2jBw5MiwT48z7OeOCMc5xnGXyMZS2VCw/Y5J1wHtxMosAphmL3pd435y7RPZ3Zu09jnvZGZ/8zkXbmUGDBpVpzlW8XjJJN+F45e2Q2e1Ec0mW2Wxd38wso/hM5lzF59aMG85baLfF796MVZ/j83qMVcYL8fviWJLZrEUWRrwOx2v2KZYj+t7J/sy2zKwCo7hul1j3+GJd8v45T6R1JN9T+JhQxTopsmZljLMc7A88zv7j9ZzZErD8bi9itm685bjLOcJ+++1XpjMrkdmzZ9eUoWOZM3sbz4N9p56VqFntPfp9Rda6ZrUxHr1rYhuzb2TvAKK5Sjt974z6Y/YukDYyjIkor8zGlHHn42JmJ8M2zGyavHxV3mt4XPJeeA32TdrycF4RWVgyLmmPUs8mjTYomW1XVndeH/wenJUjsqeObJA7Xq8ejY7lWskthBBCCCGEEEIIIYQQom1paiW3v63nG3z+ssxVvv/yL/9SpvlrsP+iEW3kYlb7C130i0C2yqzK6gk/P1vFw/z4C2m0kVO0Sq5jOfgLoa/Y4K8j/DU12zgkWuXKX/xY5mzjhWiFE8uR/Rrkv1bx/rgCKCt/tLoh2qSuK+J1z7Zj2e+9994y/YY3vKFMRxu3ZRsZZGqFaIMatm+WjlaDZitB+WseV8NEq4zYb/krOvtM9CsnV6zwF/psYyXva1zpyHEg2kShI36dbHMt3l+0qoX3lykUWF/Rqp12WTVitq7MjEXWERU6xx9/fJmONpzLVCHZStJoQ71sEz22RbThX7TytWMePO5ptiX7On9F571S9eDl79u3b3mMSodoAyezdbGZbYqTrYgk3ka8p2wjEsajj9UzZ84M/x5tlGUWr95ql1WBkfKM7TF16tQyzdVRfK5Fc4bsOUoiJQrrleMU6zhaucG+wXEqW7XvczKO9VRpMN55TrSKzDc/M6u9V+bBePY4yjbuyTaGi1YGsv8xPz5Too2+qbzK4jpbjdlum/JF82LWFVcTUXXF2HGyDR+jVcw8no1h2SZlTPt4xnM5NmcqoSVLlphZbfxyDsY8GCO8js9vuFFp1qcY417nLE/Wh7N5h7dRtvk8YX6R4q7K5sCRAi5TxXU16pWXyhuujGc/8Dkp64qr0rLx2Oeb2dw724iL45wfZ3n4/ZLjKvPwvsvPcdV6tskor+3zGfb9bHP2eiu8M9VZNlfxe2E5o2dFx+t5G3L+yHpuJMbbjey7R6aSZ9378WzFNsdBEqlDmEc2n2WceTlYnuw9BWOcSquonIwd5sFzvG44L8jmFhxj/X55f/zuminmIwUqnzfRNTrmF9V59m6rHu0U69F9ZRvvcpyKvmNmG0hmm6L7OVndMp6jd3C8Jtst2xSd4+m0adPMrDaOGHfZ+4dofpNtth2pYXicZWPdZqvgo3exfG7x/rL3kN6fIkeNjtfLaDa+tZJbCCGEEEIIIYQQQgghRNuil9xCCCGEEEIIIYQQQggh2pam7EpcBhBJFs1qNxSg7COyzqD0KZNGEl/mn22IQ7gUP9oQINscgVKJaIk8l+dTFpdtwsGl/X7NTPrJ+qgnFePfM+uVaENKbs6ZbX6Wye8cSiXYFrwe825EdtMV8LqiBIXtTpko7y3aRCXa9MCsVoYXbWiRbY4TSYw75u3XzGSIlMWwX/p1KGtknGWSZEpFo00ZBgwYEJaDUjYn20Q1i2vet98v74n1zPxYvqhvZBuKckxrp80+Iupt0pNtbBXJH9kO2aaR0YYWjCnGOcegzGLA25BjLMdkSta4maSP8b179y6PZVY8jBPalTi89ogRI8p0vU3IMhletMEkP2e2rm6yjXp4Pcar501blUxyt6Xg8cWYZGyxXufPn1+mGUd+PuuKYwznDNG8iHGdyX6j8YjloFSZ57JfDh06tEz7Ndn+7A+RJYpZvFGbyyzNzAYPHlymWQcss98XradY/irSQ2839iOO5ZG0nXmzb2Tzj3rzknaRvkcbfmXl5aarHP+cbOzINnjycSzbkD3bJD6aN9EaLduojmOsX5Ptn1lQsW44f/O+TUuXzLoksnDiWMr+xTGG982x3usp20SV8Rdtypl9/9iY8t/NTXafrFeOYdFGv5klT2a5EW0gnFkuZTaD3sbRpt0d89t3333LtFvyZNZPtC7hvUTzGZ7LuM2+z/l9ce7EvphtKs778nkJ4zezZ4g2+4w2PutIFsvtGuMZjDOOc9EcL3oem+Ux7sez+Xhm7cW29O9PmSUKy8Rx3GOcf+eznnYmHPMY414Oxjg3LWQcRRsNMl8S2SN2TEebm2dWL3wWRBayWT2T6LtwO8V6dF+ZnSDjIBqrM8skjkFZ3k62uWjWV7wNM1tU5sdNUzn3cnh/tNjj3IxzfD+f84usnJGVWjTnM6u16GR+0XdTxjCvkfWViGxunb3nqfe5DK3kFkIIIYQQQgghhBBCCNG26CW3EEIIIYQQQgghhBBCiLalKbuSaKdmSgAo2eBy/mjXZi5757L9bOdqX6rOv1Pmle2KHuWRSY+zHeD9frNl+LxXyioo3YnkcrwG74syAa8b1jNlGpnVSLTbaSZBYpkohRo4cKCZ1Ur1eS5h/VPGH8lpWf6uhtdbJmdmHTNuGeN+/2ybTM4atXsWI4yzSD6bQZkK24Yx5XKzLMYziwXKgr0vsV9ShpxJa/xeeB+878z2JZLXRbswm9W2Fc/xuqFkJ7PeoQSOtLNsLLIRMauVdTNm+vbtW6a9jVhfbHuOw9F12CeYR7Yre7SrOaFNAdubckSPNR5jGzPfzP7E75ufoxyNz4DIRiCz5WFcZhJfz4P3R8kqx2/idcPrZbYA2VjdTvFNsucNxxi2Jec3Pu5F9kodPxeNTawzxjLzi+zJeJzlHDRoUJnmczmyr2H8Ml445+HnWFaPd5Zz4cKFZZqxE8kIMzk+Yfmi+GPdspx8prDOvUyZ5RDzy6SR7Rbj0TOVY1g2N+DnIrl2NgeJrpeNZ5mVD+vYY5zHMispxoBfh+MkLcc4X6FtBfuaP4c4h6H9GsfVSPZfxZqL/SSag7Bu2Q6M8eh5wmuwbMyP1JMCd2WiuUoW15n82duSn8ss0qIxnXWc2WQ+88wzZZpt5rHBGKGNH+dXjGE/zjijfQifJz179izTkQ0ZLc0YW+w/7Nt+v/w750C8NvtJVHeZLQnri3l4ObKxuEr8Rue0w9ie3RvrMLPw8zpknGUxTryd2KbZXIV5R+9O+HfGavYs9+8ZjN/o3UXHc9iWbunA+OT3l8zm0+uOzwf2DfajzB7K88hsHlh3zCOq82bjs53nL2b5u6jMrsjjnHWXjU1Me368RhbnmSWQn88453wla28f9zj+ZXOs7LufX5vXyKwMifdTxjnfd2S2u5F9Ka9Rz8bOLLefboRm369oJbcQQgghhBBCCCGEEEKItkUvuYUQQgghhBBCCCGEEEK0LU2tIXfJEyWzlPZOnjy5TE+fPj08xyVKlGBRWpLZdviyd0r+MqsOLp2nrMB3lc6khVxaH8kEWGZKzCiB43FeZ8899zSzWqkBr005D/OL5ImEMo5sV1mvD0o6KBnicdadyxv+8Ic/hOfSuoD3ynaJ7E26sl1JZGfAe5g6dWqZpqSFUipvS8YTZbeZnCaSW7OuKDmMdok2WyetoSyX/YtQsuIyycwyhHHGHX95vueRWb2w/OxLfj7/nkncMxl8lEcVub632+23314ey2wxMuuGdpRGepkzqTPrlvcdSR45pnF8YFtQ5uRtxfE9a29em23h5c/krywnrUSGDBnSKd/MCsp3fjerjTtKuqLr8b4plYxilON+thN4ZOXAPCin4/OHY5LXE3erz+SrpJ2lkB5fmWQuk4pyvPc8MqsoxjvbwfPj2Mx8sx3ZmYdfhzHEsZdMmzatTPfv39/Mats3s1bgc4n4GMkxlJ/jsy+SVzI+OfYSzhNYB15W9g3GMuugX79+nfJlv80skbK4brex3OMze+ZyjGX8RVJZ1jfbmm3JZ6fD/sV25HHmF8UU457xwnIyxl0Kn8mbeZx9l/Xhbc15LOsg69vR31lHvO/sO4UfZx60quDxvffeu0x7fWX2d43EaleOa1LPEoBjG+uYsej1xrrKvvdE4z+vl1lm1otF9p3seTJhwoQy7TJ4XoP3ylhl+ThX8brL7Nmy7zAez5xTZfP+7PuK1xOvx/ka+4w/s8zW9f/MFiazpGk3Gx7iZc/ugfGUvd+ILNWqWMFGdlXZnD+zD/EYz2xomcekSZPKtPeJ7LthZrHKe/FY5DuN7PnFcdyP8++M8ciiquO9RDGePcsi+8xs/t8I7TKOZ2TzMsK48zhv9Pujz6Mzu5/M0pTneNuzLbM25PtQjyX20chys2MekS1tZJ/SsRzRdwY+7ziu17tXloPPVD7PsneLUZxXiddWjOtayS2EEEIIIYQQQgghhBCibdFLbiGEEEIIIYQQQgghhBBtS1N2Jb4kncvbKcelRcnf//73Mv2Od7yjTPuydS57j3ZANatdpu7Hs11Usx1TI5lZJoGnlIpL7l3SW2U3Ucp/e/XqVaZdTpHtqEp5QSSv4XJ/1hFlGpE0nvfFuqPUgNIF7oi8cOFCMzO76aabymMDBw4s07StWbx4cXhtp10kNV7PrB+2B+/t3nvvLdMnnXRSp3Mi+ZhZbftREuXHeT3GCD/H+uQ5ns7kZrwv2pj4NdkfmKZEZs6cOWWaO/b6vUQ2Ima1/SuSAvO+MwkT05F8ldfOJJXsl/Pnzzczs3vuuac8Rvk9rZbYtzPJZLvgZc7GI+4y/sgjj5Tpo446qlMebMvM6oVjpLch24fPgyzOI7k944uxwbGQtjweH5Ro8nM8zucZ4z+yDOG9Mi5ZNy7vpEQ4k/JGMjWmWUeEefBePHYffvjh8hjzcCuvjvfF/hRJaruy9ZTHDuuYcUhbC8YfZbjerpQy8v6j8ZuwT2Vy9kx26W3Na7AcvBda8riFB8f6rA5mzpxZpnm+X5ufY4wzJnmOj518tmRzFB5nmTymMqsH5sHP+TyGdcG/cyxhe2c2VO1A1P8Yc5zbsV+zv/v52XM7m69EMZLZQDF2WD6PffYpjo+Ec8wBAwaYWe38g/my3RctWhTei5/DWOD4mVl5eT9hjHOux2dFPQuLepY9ZrXjg49ZrIvIxqvjtaM5ebtYUUVl5zHWPWMxsi/I5uTZ3MFjILPOYDloN8O29DjjNTJ7Ss5VBg8ebGa1zyPGMuOM/ZzX8VjlPIQxzthhmb1MnPNm3x+z79sef7y/rD9w3HDLNVqvZXYsmeVYV47niKi82dy33twws0SNxm6zdfESWYua1dY9251zEY8jjrvROG+27r2Cmdlhhx1mZrXfuZ588skynVl/8L68TNncKBsTvQ7Y5zKLRRLFWRbjmS2dP+NYR8wjm1e3syWP2br6qnIfrBvGh48bmQVGFNtm62IwslTq+LlsnurXZnxlljoLFiwo04cccoiZ1T4jsndmjGOW1eOc/Tv7rhzZt/A7PcfybP5AfOxgPybZtb1usnrOxpws70bQSm4hhBBCCCGEEEIIIYQQbUtTK7n9l9VsRQw3cPnrX/9aprn6z1eQcBUEfzXJfpGJfjHLVjfzV4Pol0X+WpH9IsMVzf6rG39Z5io5rgjjr6U8x399YfkzE/rs1z0n+/WD5zI/r99s46FsJeMFF1xgZrWbMfDvXGHw2GOPlelslVQ7EK2EzjZ55Orf4cOHl2lfgcFfrLPVPdEKNh7LNn3KVkT4+WwnxgLj01dDMT/+UsdVUlmMZ5snOPxVkr9i1ovrbDM29u1oo7dsFVi0oZuZ2TXXXNPpcyxztkFrtDlRu60eMauti2x16ezZs8s0NwjaZ599OuXB1Q8cV6LNoRgDrPPsecC297ZgTEWbW5rV9k0fw6lIYTmfeuqpTuU0i9UzUb8zy2MmUvNkqytZp7xHJ1tl6Bscm9Wu7ho3bpyZ1T7jWOfRJkIdy+H11C5xHpUz2vzNrLYdqFQZPXq0mcUrKszizbF57WxDs+z5yzJ7+3As54oW5u190WxdjA8bNqw8xv6XrWBnv/O8GePZJmXRJmqMM/49WkHc8dpep9HGcWa1Mc7VMK44qbeq0az+qr92jvFsYyWOiYwBjx3WcbaiOdrENVspx76RqSD9fObLdmKMDB06tEz76jufa3UsB1ctMf4YZ1E5s82bos2lsrGUn2P5WaeR8iyLcc4jXV2UbZJF6m1k1m4xno2T0epMs9pnr6v3+HzPFAqknnop2+w0+o7GOWg2x+f8ylWGI0aMCMvJlYL8Ps2866085LOF42a0+W8W47zXSAma1QsVlVyh7n0320yOtEsMVyVbFZ9txMY28e9rVTYkJdH3l0wxm/W7SJ3JcZzvENju/hw64ogjymPsU+yvXOEdlT97DrHfsRyRkod1x++/mULJ84g2ozSrnSfxmeRzKZ7b6MaT7Rj70XcIprMNfiOFDtuN1BvLs2tkY3x0PNrQ0qz2mUMVuI/P2VjO8Y/vGaOV1ZnCkd8NopXaHCsyFWX2/cj7dbbRMuOcz6Jo0/lG47zZ751ayS2EEEIIIYQQQgghhBCibdFLbiGEEEIIIYQQQgghhBBtS1N2Jb5cnDIBSmK5zP7RRx8t0/fff3+Zftvb3mZmtTJTyvgoJ4lkf5TA8HqUi2T2Ei55yEzlKWuMNmGiXJdyAEKJAmVoXiZKJSijyTbMi+QRzJe2KpQBUIbmsgLaU1D6wzZkW3m6d+/eYdm4KQnrJpJctYu0xsvOeGK9MUYopZoyZUqZdvlvJt3NZHjRZkTZxhVsh0h+w2uzf1GSyOMeU9mmJrxvSlOiDTCzjUqYH6/tn4s20Ox4biYb8zStHdg+tKzg2OTWBFF5zMzmzp1bptlumUS4XYjsn9hW0WYbZrXjtsc588g22mPsRmMh2zKTjEfS+2yDVT4bGK/+OW5Ox89xbM2sITydbeyXbcwcbV6VbfDG+orsgbJNpyhznzBhQpmmHM5hfVJiRjKpajsQydyzTSMJNyLy2IhspTqmI4k6+0YmC8xsgvw445P9knMQPqM8zripZNY3MosfLx/7AOurXoxnm8vwc5ktm/cr1kv2/KHths9BMol2FRu1dovxenLW7DnFcdwtEhgLVdrMr51ZlmWS36jfZZtjZzJ3f7Zz3kU4H+X16lmlZBtlM+1xlNlnZHM93penee1sg9ZoPOI12mUj4GaJYrzK/Iv15rGT2VpmbenpbNO+zL6B+DmMcZLZ6fj4nm0WncU47yX63pxt0BfJ4LPY4rksP497n8m+m/NeonljvTbpSHS8XTZXdbKN3zLbPs79/Dt/ZkvC9ou+M2bzk8xSI9oANLNhy2LHLTyyGM82m4z6XZWxNNqkls+v7Psv5+7sS6yziMwuzK/TaEy2Qwyvj3p9NOuv/F7m4172zM3iPBrLs/G7Xpxn388Ij/uziHHOuMvenxGP6WxzV8Y5+43Hbvbsy/ob5zx+zcxOiUT9NNsguErbV/17R9r7DY0QQgghhBBCCCGEEEKIVzR6yS2EEEIIIYQQQgghhBCibWnKrsSlq1xaTysLSme4xH/atGnrLvz/lvBzKT+Xy0e7M5utk27xXBJJbTviy/257J3SXVp8cDm8S91o5cGdck8//fRO1+h4HS8fy59JnXncZQKUEWTSuWzXWD+H18gkPFdddVWna7M8e++9d5mmHUQm82k3eY2XPZPjsR1ojUE5k5/vuwCb1cpbaDUQxS3rMrPRiCQ0/CzlPZmcMLKLuOuuu8LPvelNbwqPR22dyWKyPurnc1xh389sLSK5KcuTtc/NN99cpr1OszGD9Z/dV3Ssq8d9tNM025X3wtjlGO+xRhk525D1z34TWQLRFoGxm8kcPQ+eS/sGSgrZFt6ed9xxR3mM49hrX/vasMyRrIply2Ijsg2iBJ/xRVgfzM/HZJaH9cy8H3/88TLtdZ09a6tYObQbXkdVZG4cH1gv/jzfbbfdwnNJPeuSrL5Zpmg85d8Zq3yGs0x+Psdy9ttDDz00vDbz8L7G+2DZsl3WPb54f4zxqC92zNtjPJvnsJ/TdsjLkdlnZGP8lkAmc8/mm2w/Hxs45mR5M0Z8XOIzJOo7ZrXtV88WgTZ4mY2V53ffffeFfz/ggAM6XaPjOVGcUQpcL8YzS7XMkjCqg8xagM9iftfwMmXzoFcKjBfGHNspstLjmJlZvjBGonllZpHGvhZ9T2XZOG/J5gB+/JFHHgnzGDx4cN0y+T3yXlkH2Zjo8cU+wL7B+uK4weNeZ2yfbL7GeWM0/8gk7tnzPBsPuyr1rNWq4O1e5btT9L4hG4uy+U5kBZHZFbJMkV0ErfWIW2mZ5XMmj0uWn9eLnlksR2btyfrn94lo3ODfGdeZZaOXo4o9TSOWDl2dRqynsnuNvstUeYZHFp3Ru5OO+fF8T2e2gdk7OL92Zq9GC77sXaCXKbN8ir5j83PZ91XWczaPcbI4z6xmo3GtCtH5jca7VnILIYQQQgghhBBCCCGEaFv0klsIIYQQQgghhBBCCCFE29KUXYkvX6dktk+fPmX6ySefLNOUUi1cuLBMT5061czM9txzz/JYZvHBdCSDzaSMmSzH8+C5vDZl9zyHliwOl85TGkO5OD/n8sNMTpdJDl1iwM9RMsA6ogwtksux7mjNcumll5ZpyincmoTyD8q2H3vssTL99NNPh2VqNymly1HYNtylndIVxgtl0y4v7devX3mMcUYoMWHaySQ79aTXlIzw2gMHDizTUX+lfU8mU8skjPUk7tlxJ7Np4X0z9hnj0c7J7JfXXHNNmZ49e3aZ9j6TSVcz6Vm7xXVHIilSZmVB6TTHfpee0i5j9913L9PZLtYuTWUZsnhmDEaSLn6O6b59+5bpGTNmlGmP/8xqienMcsFjl3/P5JGsg0huxnjmtSkBjnbvZnko76d8n2OSly+TfGWyvXaT/UZk8j6SWXF47DPGGYesHz4DPTYyaWRmacO693Jnf+cznM9tfy5RGk4yaTtj0cd1xjjrjp+L5I58lmUxzvkK68nzY4wvWLCgTHP8Zj/28YRly+y5quz2Hv29q8IyZrY5vGfWbWTjR6s1nrt48eJOeWfzvcgyr2OZIkuebJ46ffr0Mu0xHkljeU8dr8f5iqczuXr2XcPHUo7zmQVEZgHjsF74PYnfI/j8jcavdojPVpDNhSN7EbNYts6/00qG4/uSJUvKtMcAx44slrPvnX48KyfnDpyTR89s9inGeBaLnkdmHcFnRGQ7xM+xvtg3GOORtJ/3yvEj+x5eD9ZHFvvt1iei5042fmYWpV6H/BzbjHNzzhmjuXc2d8/myn4+r53ZSfCZHX2PyqwPsnHa6yCyeTWrnYuwvvwcfsfLxvHMAsbTvA/OA5kH78Vp1Jak3efjXv4q30My+xmP88yKhzHP8c3TWf+JLFHM6sd5NHc1q32X4udnc2+mGa+8l+h7IPtYZOdmtq7uaAnFOXn0jtQsfkeQ2Ztm40I9Wx7S6jFbK7mFEEIIIYQQQgghhBBCtC16yS2EEEIIIYQQQgghhBCibWnKrsSXw8+dO7c8tv/++5dpypa4rJ0SrF/+8pdmZvbRj360PJZZK1Bq4OdkUoNMokpZgS/b59J5yiGZB5fiP/PMM2ZmNmzYsPLY8OHDw2tk0lxfip/JzbIdZilbjM5lOetJdNg+/Ny0adPKNO1IHMqzKe2hvCCS8Jitq+t2sXfwclLawTphvbEu2E633XabmZm97nWvK48xLihjoVTE24ySqUyqFEmfzNbFQCbLpPWK26qYmQ0aNMjMzEaMGFEey2KSUjHGuI8PmfVHFi9ejzxG2Uwm+YrsVPg52lRQ4h7lzbKxvbMdvdt9l2uP8+yeGNuUivGcBx980MzMRo0aVR7ba6+9yjRlv5G0NpO/ZrYIkdw+szlh21MG7nFOO5PMGoTPBl7by5rJZqNnjtm651lkYdIRjjORjJj1MmvWrDLNZ20kq2b7se8Sfi6K7XaJ98gKIXvmsp0Yi26RQOspxjjPjSy7MtulTCbJ8vkzI7JtM6t9bnMs32OPPcxsnd1YR9i+7IO8TmZ14jDGOe/wuGVekTWVWe34HdVBZuXAe43mFZm8OovbdrUoMYvLmc0TsuOLFi0ys9rxjuMSY5xzl6hdM0uYTKbr1+TfGZN8brN/+TOatoiMp0yyT/z8bG6ajeP+Od5/ZlXBeSTrJrIh5BzSv3Osr3yNEM1X2mUOE43jLHs232TdeuxwPM6ehZz7eBtnkvTsu2ZkucNzOY5zbOO1/TmTPW8IvzOwPjxWs2cP4XeUejYthM8KnhPZg/Jc3ks0F2k0Jtt5HK9HZjvFmPKxhmMV/84xhfNLz4/nsj2qzFW8rTOLEtqj8N2EW6j4nMWsNi7qWZSYrYuvzKqI8Ljnnb1TItl83J9lrIus/PXG8S3RhqcjUflZL9kzKbJgy6xJOUdhu0RxnvWlzPKNMegwfvi8Z5z7+EzLoKhsHe+FseZxlT1/suPRu9OsDrI5uZ/Pesks4aJnTRbbmY1JK8ZyreQWQgghhBBCCCGEEEII0bboJbcQQgghhBBCCCGEEEKItqUpuxJfQs4l9JQG9O/fv0xPmjSpTFNS6MdnzpwZ/p3L5SMLDC5Zz5bf8xzKZ3x5fXYuoYTdl/tzx/ko345lplzByXai571Ekq4qMlDKHSP5Gu81272dEme/DtuV8lG2PYlkOe2yK7CXPZNfsF4pTeHu7I8//riZmQ0dOrQ8tt9++5XpehKmLMajcprVSrA8XrL4pFR4yJAhZdpjjjY8jC1KsNiWUd6U9GQyHMZOJKdhHWQyZOLXZDkpa2dfZFt5G/LaPLeKlKodiSTAJJMzcRyaMmWKmdXKXGltlFmJeN78exbnjEHGj5eJsU+pL+N83333LdMeS/vss094LsdCxjGP+/mZRJifo8Te+0p0H/y7WS5D83My6yn2leg5l8nwqkgl25Vs13rWPWOAbTlv3rxOxzgPYH0yFqNdxavs6s6xyY+zHTOrK8a4f44SYMYF44VWDjzHx1HGOJ8pvNfIVi6zaWGZs+egX5v3nUk7o3OyZ8SWFtdm8T1lUvPMnsGfk6w3xjg/F0lps2tkMcBnRGSdkUlmOV/xsvbp06c8xr7BZzhjnLHqefP5xRjnvbLPeHzycxw/qtgj+XH2gaweIxl1ZjVF6s25t5T+kNVFZM+T2Zkw5iIrlKwus7kuY87zyGThhDZq3j6UuHM+we8f2Tju98hxnHHNeOfcOZqrZLGazTm8D7Jv8HrZd/1o3pjNvbN2abfYrlfe7P45vvj8gnWcxRz7jH8u+66T2RXS1tGP8xrZOM7P+TUza0CO14xhju+eN+ck2bshlqnefJywnnm+90HWZ2a1Gb3DqWKhtiXZZDZS5iwePZ3VBes8ev5m86PMaqdXr15l2tue43BmqcnvaH6cMcxzOb9gfLH8Xm7+PfrOYVbbP6LPZe8QWXesA++H2TvLzD63ntXYxrTl0UpuIYQQQgghhBBCCCGEEG1LUyu5/S0+f/HwVatmZq997WvLdLYyz3nggQfKNFeP8hc9/vLgv0xEKyrM4k1dqsBfQrKNLP1XhezXa/6aX2/TLn6O9chfXupdm/cXrYoxq/31yTcjO+CAA8pjv//978NrR6vb+OsTf1nNVh5GKwTbZeNJJ9tMIPr1uiMet75pmVntytGszby+GePZatfsuLcD65vtx5WC0ao6xjI3KmFfY37RZh+MT652ilb8kWyVDcuU/Yrp/YPX8w21zGpjnHl7Wblyje2drcSKVr624y/rJPs1lkR1nm1ow9VFjEdP89xoZYlZvjLQz89W7HH8i1ZLsx9wg0yOz9kGPn6c7c24Y/xE/TTbvIP9KtuQks9Hh/VIog1woo13zNpvfG6WbHVIplbwNmMdM832iMaHrI6zsYRx5MczxRCfI8Rjim3NVa6McZ6TxbvDlSns89Gm3/U2jO1ItFEy+zZjmflFG75mSr1Gxud2GcsbWa1bb/Mf1jef1VxJGuXNGK+yspLjtB9nPDHOsk0jPb6yDQMZf4xxxob3bR7j9bJ5gM/bM3VBtqKQ833vBywzx5XoerzOhqzua5fYXh9V7iFaOck6ZrxzXs/2izYIzVbcE47T0TyJz/psdZ9fk7HMeMlW40XfBzi2cwxmftFcN5sTZmqaaLVgtsIw+w4TrXzfUjfic5pVOrPePAb4rGS88ztO1E7ZivvoGma1cxU/J3svEs1reG3GAsuc1UsU44xDXi+LuUhRSrKNyaNrZ6uHWV/RvbT7yuxW0KhKI4rX7L1atNI5i/Ps+0C04Wk2/nHuEt1XphrO1KbEY5D3lL1fiZ4HPJY9L7L6iN5xZP2m3ibX2blEK7mFEEIIIYQQQgghhBBCvKLRS24hhBBCCCGEEEIIIYQQbUtTdiXR5lXLli0r00899VSZ5gaGs2bN6nR83Lhx5bGRI0eW6UMOOaRMczm8L9WnvCrbYJJL4ClVcVlYlQ2gIqlUJj3O5DzMO/pcPVm7WSyHoTyC16CslPU0ZswYMzObOHFieYz1v+eee5Zp1tfYsWPNzGzBggXlMbc+6UhkAUHaRYoTlTOKIbPamOM9+zm08undu3eZPuqoo8o0pWUe45mtRySjN6tta782y5NZH0REGzOZ5RvsRHmzvjIrC0puvL9Gm6iY1UqIKQeKYJ37xogdyxlJRbMN0Ug9S4d2kp5FMqJs85RsE1Cvuzlz5nQ6ZmZ26KGHhteOJEqM86xMjCW3iYjGWLNY9su8ufEu86DMmETxyr7LsrE/Rpuv8XNZnBMe93thnPvmiLyGWa2Uzdstso0R+bPfY4NzG9pQ+TPSLJ4HZBuLZpvLMC6jeKgnf+Vx5stYiGTGZvFYzr6zfPnyMs0Y5704vHYmr2S/47jheTOuWS/ZZrTRxkJVqCKR3dKIbPX47F+4cGGZZtuw/Tw2GCOZFJj9gLHjsVFlDh1J7Hltzg3Yplke0d85jtPihzEcfR9g3+F9s75YPh9PaJUVbVTf8V6a+XvVc9qJKnOt6PnGcYnjGduJ8e6xGG3Ka1bb1pkU3WMn20g327Tar8nrRZuZdcw7snvjNTiWZuO4fy77/sH+kG3W6jHOfsT7rmKlVI9X+saT0XyA4xLbNxuLPL8sxrO4jTYcZfszD5af/cvz4D3x7zye2Sr4cV6D71bYF6PvNdnzP5sbRfXBuqhi8RPFbaPzED+n3WK9Kll9RJYhjDXGeWQ1klkDZ5vfRvOAbMNelpN9zGOQMcW/Z3MJnh/FOecrvK8ozjO7xMyaiNfxOOd4kr0naXa+0uq5t1ZyCyGEEEIIIYQQQgghhGhb9JJbCCGEEEIIIYQQQgghRNvSlF2JE1kzmJndf//9Zfrkk08u05Szu3xgxIgR5bG//vWvZfqwww4r05GMl8vsXbLesUxc7h9JHympqbLLqBPZp3Q8TtkKy+pp3hPh7t7RzuqUTDDf7F4HDhxYphcvXmxmZt/5znfKY7TPyKRJLvugXQnvj/WV7T4eSYLaDdYJ5RqUm0RymV69epXHaC9w4IEHlmlK0l3el0lwsx14KceKLHn4d8YtY9xjjpKXyKLBrDYGIrk4JUCMScZ4FC9ZDGU2Qez/48ePN7PasYQxxzZkX+rZs6eZ1d53Jj/aknZ49/vKYiq7J7aLtwXlYfPnzy/TtKHaY489yrTHOcc5PkeicprVyuk9NpkH4zyzuPHz61nnmMVyTKYZJ4zRzHbE6yuSzHe8HvsKy/TII4+YmdmECRM65dsR1pf3dZaz0bhtxzivSmZ747D9KRGk9JsWbT7GZHMGjs88HkkOGRf8O+OaY6HHQyZFrzK2+v1yLK/SXz0/jgnMg32G9mq8r8cee8zMai0zSDbvaJZ2spnqSL3d60l2n36c9cr2YIz36NGjTPuzPZt7R/Nfs9qxNzqWzUE4JnqssmxZ/8pi3MvE/sUYp9SZeH3xetmciOWn5ZFb8dAyI5MTN2vl0G6x3ArqWTgxLtjWHCsjezCS2f/x2szP25IWCtF3vI74tTObK5Yjk7tHFmkcjzMbQq+nKp8j/l3TbN0zMLNAjCyTGqWe9H1L6QONSPw5dnBs5hjFuUNk75XZFmTfK/15n1ltZv3SY459I3sPk+UX9dGsz0SWIZm9GWOc13766afLtPcP/n1j2p9tKfFMWmG3xTaMLGDN4jiJbFDMasdQxnlk7ZPZ00Tfldk/su8AmSWh58G/Z/ZwxMucWbNk3z94vqez+XYr5tCtjm2t5BZCCCGEEEIIIYQQQgjRtugltxBCCCGEEEIIIYQQQoi2pSm7kkh6QfkHpYMuxzMzGzVqVJm+5557zMxswIAB5THKP2666aYyfdppp5VpX0ZPmWy2m262LN/LWkUmFS2/z3Y9zXbe5dJ/z4/yxcyKgrjshufyniiLoyUA87v44ovNrLZ9dt111zLN+j/ppJPK9BNPPGFmZjNnziyPURJB6UIkqzBrX3kNy8374T1Tkk3Jx7Jly8ysVkJOCc3DDz9cpl/3uteVaW+TTDofyXzN4rpnOzGdSdm8/JG9T0eYH8vkfYJ153YgZrV9g3IfT7O+GNdMU7LMa992221mZrZkyZLyGPsMr8fYd/kQ+2o0ZqyPejtld2UyaV0kXe14vo9N/DvjnBY9bFsfwyOLBbNc0hXBmOKzIbLiIey7GYwD5uHPFx6jjJ9509oiKnMmt2d9MHbdDoZ1lEnyWHeRtcCWTtQHM1ldJhH0OuTf+bnZs2eXacrc/XO05GAeHLvY1tF4wzlDlblG9AzInmeZRYLHHPszx03eKyXqfk3eH+djmdyZMe4yd9ZXNu9o9TjbrmN5No5XOR7FONuJz9QoFvmMZzyxfZk3z/Exm/lmcvXI9iCz4cme4Uz7+M3Y473w+cQY97yjuU/HMnPuwvP9+0o232rFHLpKHGzoNTYljcZ1ZMWWWevwOxDb0tMc+0hmP8kxyq+ZjeOZRZpfO4txkn1H8frg84FzaM5V2M89P+bLcrCemTfTHu/Zs7UVMZfl1w7xvLGIbMfMat+RsH48phjjmRVONP/n+fxcZp/A+YCfn819I3uRjnh8ZTFO25GlS5d2uhfmy/Lz/liOyNoyqs/1lT+ycd0Qy7FXOtm8Mvq+w2d8lWcu499jl+3A753MI7IxyeyhqjzP/HzeH++Fefs7KH6uyvdLXptx7nVQr2xm8TOqXj/YGGgltxBCCCGEEEIIIYQQQoi2RS+5hRBCCCGEEEIIIYQQQrQtTdmV+JJ1LkennJVSj7vvvrtMv+lNbyrTw4cPN7PaXcW53P/Xv/51me7du3eZPuKII8ysVgLgslaz3AaEEp2ozFUkvb68Ptu1NZOKRXKeTPJFWB8uw8hkDsxv2LBhZfo//uM/yrTbBowcObI8xp3ejznmmDI9d+7cMu12JZnFRSZvje6rnaU1mRSJVgOUpbvsK9tJ/MEHHyzTtDkYPXq0mdXKZNlPCOOdceltwr7IcxmTlKn4+exfmaUDz2G8e5qxmu3iG+0AzzqiHJ5pWvJcc801Zdrrac899wzLyWtHdfBKkzpGUr3svjNJdb3dyWlzxDFtyJAhZlZrWcWxPNvNPRp7GEeMuyzOfezkudlO8tm47sd5T6w7xhqfRZR/OYxnjgXss3yWumyef2eZM+sSLx/bL4v5RnaEb4e+0ugO9zw/im3GBZ+jM2bMKNM+d2GMZ/2IcxDGqscw/8725fgW2akw9iL5fEcieS6fAVmMUDLpz0T2DcqgMwsg1mNkD8d7rdeejdoakGhcbDcavedIisoYoeUGbTu8XTlfIRxjs7HZr53ZfWQy2Ehiz76TfS/htT2PbBxnmvJ3r4NevXqVx1gHnBfWK1Nmexg9b0g7jLutpl6fNcv7rddhFuOsY8Z7PVszxnX27PXrZHY6PDd6rnNeQMueyNrErL4dBP/O63Ec9xjv06dPp/KY1X4vYTlIZIOUfZesN95WifctzcqhyjOoXuwztlj3WQxEn8vscqL+E32v63i96J0LY5z9L7NSiOYw2Xc89lF+l/TvHIz7yDatY/mjuVu98Toja+MqY107U6WvNtKfs3dRkXVPZkVSxY44mptm3xOj+SvjPJsnZN9BI8tjxiXjnHMaf2bwWcbvouxvpJ7dV/QuwKy+TWR2vNWxrZXcQgghhBBCCCGEEEIIIdoWveQWQgghhBBCCCGEEEII0bY0ZVfiS/ijZexmtcvNuTT+rrvuKtNuO5LJVpnfFVdcUaZ9R9ETTjihPEYpFfPgUvxo59xsl2vKCSk/dMsTXiOSXHbMj3m4lCa7b8oLKLvxsvIYJZqUTP7qV78q0+PGjSvTLmPg/R144IFlet68eWV64sSJZdrrLpMJ1bMoIe0iJasnV852iaX0xOVPrB/a5rAdbrvttjLtO7y/5jWvKY9lUmDKFqOyZjtGc2deWqx4fpS/UELDNOVmlMT7NXm9TArDeKeEzOF9M33TTTeV6fvuuy/Mu2N5zGrbjfUfyY/YnzO6cgy3mkzC7WMZ44HxwzqfMmVKmX7mmWc6nUtbh0w+yD4UydyZ5thKSxBv+yzO2T+yneL9eLYTNcd42ud4H+PnmC+vPXXq1DI9efLkTudzHM5stKJ2a8SyYUukkV3MSbRjeMfPLViwoEy7nJttSjuxzFaJMefly6TDjFvm57HBPsBrcEzObCJ8zsA4i2x/zGr7bmSvxs/Nnz+/TLN8tEmL5P0sW2YzUC9uG7Gk2lL6QLNS4OzvtJVxOTfHPreiMqu1+GCbce7iefPvjHF/VpjV2hd6TGXjOK9RT4bM/sWxlDDGPRZZL0uWLCnTCxcuDK8dyZMbjbN643iV/LaU2I6oF+9ZvTH+onbiOL733nuXaX4XYywy/iLbDsYF4zYaN7MY5/hJOG4y3h3GOOuDzxDvg+zbTz75ZJnmOMBr1LNpycbrRmJyS45f0oqxgTD+GM9+nLHFsZbvchiLWRw5jAuey/cvHl/8O9PZ94loTM/erXC+w/vy753+vdus9nnDckTWf6TRcbeR+fiWyIY8p6K6q2JD5eewvfldM4tzpqPxtIp9s5c5e2+YvXOIxvLMdohxzveC/q6FFlN8PlUZkyO7S9IV59NayS2EEEIIIYQQQgghhBCibWlqJbf/QpL9gkL4CwPf4k+YMMHMzA466KDyGFf38BcN5n3llVeaWe2Kn9NPP71Mc1V3tqLP0/wlkdfj5/hLx9KlSzv9natLs40so1VZJPtlNdq0iXXIOvjv//7vMs0NDYmXg6b3kyZNKtPTpk0r06wvv2ajGxlFtMuvlFE7VdnYJvpFjStIs80J2O4PPfSQmdX+ov7617++TGcbOjLO/JfxbJU9f61k3PqvfMyXv3ozrlkHkSoi+/Xd+5GZWf/+/cu0xzjrZdasWWXa68WsdnO3aDOQ7BfKbJM2p94v9a80qqx4ivoCY5Hnsn49DvxZYFa7iijbRCnbhNJhPLOvLFu2rEx7fDAv/pqfreinYsH7QrSJmVn+i7rHOVeLPPbYY2Waihr2lWiVeLYRThbHnm5kU5yOebQrVcbveserKHuiDW24uts33ebfzWrjIRrXs01uOE/ghoCueOAchWoZxjv7DPPzGM/GTY6zjE/Pm/dNFQdXAGarwrwOsk1uGlm9TbaEWN4QGtmws8qmZx4jXCHEldycx/Acjr0+H8k2m2RMMg8f0zlHieLXrDY+eU60ERVjjp9jvHuZORfhdxiuBM5WuUYbIUYKsyq80uOaNLuBYbZCzduJKjJeg896rhDkmO1zh2xTPs5VmJ/HOOcZHMepeOD3iGgzzGzMZIwzD+8/nIdzHhVtAm7W2Dy6Fd8rxTqaHd89vrKV0oxJzlU4xkZqXObH74l8LvhYyXM5186U9vWUldm7FV7b+3z0/aBjfvW+Ayl+uw714pwwvtj22TPc59RZX8s2H/a44xyG4zTfy7CcLF8U50xnm1x7PPIZxnE/Wp3Oz5nF85WMaFzfHP1DK7mFEEIIIYQQQgghhBBCtC16yS2EEEIIIYQQQgghhBCibdnKzJpeM15lyTqhFMylHgMHDiyP7bvvvmV6+vTpZTrajIAy7sGDB5fpI488skzvv//+Zbpfv35l2uUBlN9w8z3KWig/dIkjN5ehPIfWC9lS/EgmTtkw64ibmcycOdPMzP785z+Xxx555JEyTVkF24VS5VGjRpmZ2Zw5c8pjtCvJcLlFlU1uGpEjtIOcp8r9ZLIRb2vGCK1iKCuhxDGyQWA/GTNmTJnmRkiMF5d6MZYp6fXNosxqrRL69u3bqTyMT262k20U6/VEGTrrIItxj+F77rknLBsl7tkGVd5HKTnKNkqLYjiL5WyjUdIO8dwoWfzX22CG8ioS2VAxTmhXMmLEiDLN8TmyccqkvJRmcTNd7wuME/ZN9rfMjsTJ4pz1FcnXHn300fIY4zyTRBOPx2gDUP69I69k6XCj85V6eTRif8L2oHUN44zSW8p6I2kk/04Y416+bLPqvfbaq0xnGwp7uRnjlFeSaLNTbpbqcxiz2j4VWQ6RLJZbPQ5vCTHeKM32CY41Ub0xzjjvIIwjn68w9vg8YTkZRz5noLSd18sseaLnE59N7F+MT5bP73v27NnlMX43iCz/OlJv87JWxPgrMa4boRHrKh5j/DLOsg2CXYrOY5lEnN87fS5OKTvHbs45Mom7lzvrA5ntod9LZINlltuR1ptbtxrF+PqpMs5H53C84xjLduc5HqPZRtaE73A8pnguvwtUsZaN/s57YoxH1n60YOH35syipBEa2TxRsdwYzc7JeS7HvMxqJNqQndfIxlPGlc+zM5tMXqPZOM/mTX4+x28+L6rY/9Wz8W2WjRnzWskthBBCCCGEEEIIIYQQom3RS24hhBBCCCGEEEIIIYQQbUtTdiXNShwj+XZmyeDWGma1UinflZR5UfbC45QFv+Y1rynTp556ak0ZzGrlNw8++GCZptRg0KBBZlYra+fO6kOHDi3TBx10UJmmrCCSu1CatmDBgjJNWYFbOFDWnllAUP4wfPjwTnlTQpzJO6Id4DOalVe2gyynUSkM8TZhm7PNKP+iXNylUpSdZG3Qp0+fMs22PuKIIzqVmfF+xx13lOlFixaVaY9bSm+mTZtWpocNG1amDz744DLN/uOxwzy4czWtUGib4tY53Mk9swnhfbEe/XxaEbH+G5GNEdmV1D8n25WdRHWUXYNj2m677VamadGzzz77mFltG1PW+/jjj5dpxuDIkSPNrDb+Fi5cWKZHjx7d6RpmsTyX90pLHcYgrzN37txOx2ibUsWewa/drC1JxpYYw60gkjhmlg2RfVVmfcQxmeMYx3WX7TLGKZWnRQLjyPsBY5IyXNoB+dymY/micTHbnZ3x7n2N5aFUnvOLaPf2jFbb6bzS470VFj718qIEl/FOCyq35GH7M8Zpx8eY8j7BmOS8ye3XOqbrWTsxlhmrjHcfv/kdJotxUi/GW8ErPa7r0WrrqsySgfMWj33GBa3XGHOMcZ9zMM7YN9iPMuvLCM4/2Gd4He9XfPYwxnk8e8Y1i6wcNg71Yj+b43Dszr7T+lyY82TOazh+RjGeWbKxnzA/zr2jOTGvEdnwmK2bE/FzLH+VWNb8Y9MTzcmjv2dkcc6xnHHOGIysNvl3jqF8t+PnZ9ae7GOM/8gKkOVnnPNeWAfe9xjPWbpZO7RGnqubKva1klsIIYQQQgghhBBCCCFE26KX3EIIIYQQQgghhBBCCCHalg2yK2lkJ9OO+GcpDeAS//79+5dpysX9HEpwKQeI5IRmtfIUz487VNNCweXkZrVyBJf/7rvvvuUxSre43J/nsJ48P+b7xBNPlGnao7DMkWSNuw9Twk85z7x588q025SwnJlsM5LoVLFyqMeWIs9pRJaRyfjYTpQZunyKsZxJutl/2A/23ntvM6uVvVMOyRinhN2lZ/369SuPUQpDCQ0tHSL5OY/RkofX4zluRcT74/WYJpE0KNt5mNSzzsh2i3+lUq8es7Eksy5x+DnKwyjpYppx7pZUtJDi2OoxZVYb/54fpZb8XGaPQiLrKVqi8HqRBJj1STkaj/NzEY1Ym2QotutTb1fxes+DrD2yHdk5r3AZJGOV+dHKgc8MJ9stns8cWjlk5XY49+L1WOZoF3mms+eZrBw2H83OaZysP2TyWY63HqO0YWDf4JyB47/nl1m7sc/Q1oH4Z1l+zsMZ1yyzj828v8wiolkpMGlkviJiNlTinp2btTvHWz+H89hs/s7nPmPbyb5HZM8IvzaPMcY5l4m+d0YWXB1p1gpQcbtxaIUlD9uGcZhZmzqZRRXjmmNpFJ+8Np8LzC8ae5lHZlvJ/uPlb7S+oudds9atYtNRb+6SzVei73ucozDmGduZpVN07azf1IvzyMKtI368lRZ1Zl0/nrWSWwghhBBCCCGEEEIIIUTbopfcQgghhBBCCCGEEEIIIdqWpuxKyg83uOy93vmUw2Sy2oEDB5pZrTw92kHUrHYJPyW9LinnuZQGZDuceplog0KynYHrSR5ItqOq20707t27PEb5JSUDU6ZMKdOzZs1abzkzGSWPR+1WJY8tWUpZRYrk51SxeaFMxaVZ2W7QmfwlkmNRNlNFQuxpysoySX0m0fTzWRfMj9eLJHDZbsNZP48k+pnVSLMxuaXFbzNkddfsLtaezvLKYjuKXcZ5JO/tmJ/3N0rRM7uVTE7m57Bv8tpZH4uukd13ZAPEc1oRl4rt6lQZP+pJuzMJcL1+kM0d6o3rkWS+Y371LFQ4fmfPgHpjQj05+/rO2VAU4zGNjN3NWvJkn4vmD1XslxopZzaO89oe75F8fn3li2xOqsR4dLyR+WSVfEVMs8/NZiXuUR5VLJnqla/Kd4ponM7mKiSba1UtW6tRjLeWRt7bVJmbR3m3+vtXlTJ5rGbxm819WjGuypKnvcliO2u36LsraaS9q1ikNhLnJHoeVLlGPdopnrWS+/9v7w5yE4aBKIBmC1fg/qdD4gxdJZpKdnESQ/vpeysrDcFqHVS+MmMAAAAAAGIJuQEAAAAAiDWtXcnex96ftcCoasuCtdzqer1ux2632za+XC7buFc6vl6j7mbd29m6Pu6/ljPWYyPllXUe6/zqbsG1RUkd13PW+dd53u/3bfx4PLZx3VG4Va7wqp3eZ57/Cfas8Za6Zus90GqZsCzt1iS9Eq1atljL09frjfy96jxa827dc3Vuy/L93ljH9ef1XqytJZ61/Zm93v7j+h01c7fmXtnvSIlv63VHy816a7RV6jXS5qRqlTbO+EzuUfI+z952JT8dO/Pe1bPWYiPXe1Z+3NNrcdZ6/TvakvRY48fMXrd77pnq6GfYnu8iI2t1dnuoTy0R/lRn7oeZLcZa1z1zzp42Qe/6zviq3xfjZv5vP8uMtol7jLzH0bZT1vbf8xfXfMuM1jhH119SWx5PcgMAAAAAEOvUk9wAAAAAAPCbPMkNAAAAAEAsITcAAAAAALGE3AAAAAAAxBJyAwAAAAAQS8gNAAAAAEAsITcAAAAAALGE3AAAAAAAxBJyAwAAAAAQS8gNAAAAAEAsITcAAAAAALGE3AAAAAAAxBJyAwAAAAAQS8gNAAAAAEAsITcAAAAAALGE3AAAAAAAxBJyAwAAAAAQS8gNAAAAAEAsITcAAAAAALGE3AAAAAAAxBJyAwAAAAAQS8gNAAAAAEAsITcAAAAAALGE3AAAAAAAxBJyAwAAAAAQ6wuSBuDBnvx1BgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.style.use(\"dark_background\")\n", + "\n", + "fig, axes = plt.subplots(figsize=(26, 4), ncols=len(files))\n", + "for i_echo, img in enumerate(imgs):\n", + " te = echo_times[i_echo]\n", + " if i_echo == 0:\n", + " data = img.get_fdata()\n", + " vmax = np.max(data)\n", + "\n", + " plotting.plot_epi(\n", + " img,\n", + " cut_coords=[0],\n", + " display_mode=\"z\",\n", + " annotate=False,\n", + " draw_cross=False,\n", + " axes=axes[i_echo],\n", + " figure=fig,\n", + " vmax=vmax,\n", + " cmap=\"gray\",\n", + " )\n", + " axes[i_echo].set_title(f\"TE={te}ms\", fontsize=20, pad=0)\n", + "\n", + "fig.savefig(\"../docs/_static/physics_multiple_echos.png\", bbox_inches=\"tight\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/plot_t2star_fluctuations.ipynb b/examples/plot_t2star_fluctuations.ipynb new file mode 100644 index 000000000..9dfa1fc16 --- /dev/null +++ b/examples/plot_t2star_fluctuations.ipynb @@ -0,0 +1,1002 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Make figures and examples for dependence metric calculation\n", + "\n", + "This notebook uses simulated T2*/S0 manipulations to show how TE-dependence is leveraged to denoise multi-echo data.\n", + "\n", + "The equation for how signal is dependent on changes in S0 and T2*:\n", + "$$S(t, TE_k) = \\bar{S}(TE_k) * (1 + \\frac{{\\Delta}{S_0}(t)}{\\bar{S}_0} - {\\Delta}{R_2^*}(t)*TE_k)$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import os\n", + "import os.path as op\n", + "\n", + "import imageio\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import seaborn as sns\n", + "from IPython.display import Image\n", + "from nilearn.glm import first_level\n", + "from scipy import signal\n", + "\n", + "sns.set_style(\"whitegrid\")\n", + "plt.rcParams.update({\n", + " \"text.usetex\": True,\n", + " \"font.family\": \"sans-serif\",\n", + " \"font.sans-serif\": [\"Helvetica\"]\n", + "})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def pred_signal(echo_times, s0, t2s):\n", + " \"\"\"\n", + " Predict multi-echo signal from S0, T2*, and echo\n", + " times (in ms) according to monoexponential decay model.\n", + " \n", + " This is meant to be a sort of inverse to the code used\n", + " in tedana.decay.fit_decay\n", + " \"\"\"\n", + " if not isinstance(t2s, np.ndarray):\n", + " t2s = np.array([t2s])\n", + "\n", + " if not isinstance(s0, np.ndarray):\n", + " s0 = np.array([s0])\n", + "\n", + " neg_tes = (-1 * echo_times)[None, :]\n", + " r2s = (1 / t2s)[:, None]\n", + " intercept = np.log(s0)[:, None]\n", + " log_data = np.dot(r2s, neg_tes) + intercept\n", + " # Removed -1 from outside exp because it messes up dt_sig2\n", + " data = np.exp(log_data).T\n", + " return data\n", + "\n", + "\n", + "def compute_metrics(data, B, tes):\n", + " tes = tes[:, None]\n", + " data = data[None, ...]\n", + " B = B[:, None]\n", + " n_echos = len(tes)\n", + " alpha = (np.abs(B)**2).sum(axis=0)\n", + " mu = np.mean(data, axis=-1)\n", + " X1 = mu.T # Model 1\n", + " X2 = np.tile(tes, (1, 1)) * mu.T # Model 2\n", + "\n", + " # S0 Model\n", + " # (S,) model coefficient map\n", + " coeffs_S0 = (B * X1).sum(axis=0) / (X1**2).sum(axis=0)\n", + " pred_S0 = X1 * np.tile(coeffs_S0, (n_echos, 1))\n", + " SSE_S0 = (B - pred_S0)**2\n", + " SSE_S0 = SSE_S0.sum(axis=0) # (S,) prediction error map\n", + " F_S0 = (alpha - SSE_S0) * (n_echos - 1) / (SSE_S0)\n", + "\n", + " # R2 Model\n", + " coeffs_R2 = (B * X2).sum(axis=0) / (X2**2).sum(axis=0)\n", + " pred_R2 = X2 * np.tile(coeffs_R2, (n_echos, 1))\n", + " SSE_R2 = (B - pred_R2)**2\n", + " SSE_R2 = SSE_R2.sum(axis=0)\n", + " F_R2 = (alpha - SSE_R2) * (n_echos - 1) / (SSE_R2)\n", + " \n", + " return F_S0, F_R2, pred_S0, pred_R2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Constants\n", + "OUT_DIR = '../docs/_generated_images/'\n", + "\n", + "# Simulate data\n", + "MULTIECHO_TES = np.array([15, 30, 45, 60, 75, 90])\n", + "SINGLEECHO_TE = np.array([30])\n", + "\n", + "# For a nice, smooth curve\n", + "FULLCURVE_TES = np.arange(0, 101, 1)\n", + "\n", + "# logan's TEs\n", + "# FULLCURVE_TES = np.array([9.58, 21.95, 34.32, 46.69, 59.06, 71.43, 83.8, 96.17])\n", + "\n", + "# dan's TEs\n", + "# FULLCURVE_TES = np.array([15.4, 29.7, 44.0, 58.3, 72.6])\n", + "\n", + "n_echoes = len(FULLCURVE_TES)\n", + "pal = sns.color_palette('cubehelix', 8)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Simulate $T_{2}^{*}$ and $S_{0}$ fluctuations" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Simulate data\n", + "# We'll convolve with HRF just for smoothness\n", + "hrf = first_level.spm_hrf(1, oversampling=1)\n", + "\n", + "N_VOLS = 21\n", + "\n", + "SCALING_FRACTION = 0.1 # used to scale standard deviation\n", + "MEAN_T2S = 35\n", + "MEAN_S0 = 16000\n", + "\n", + "# simulate the T2*/S0 time series\n", + "ts = np.random.normal(loc=0, scale=1, size=(N_VOLS+20,))\n", + "ts = signal.convolve(ts, hrf)[20:N_VOLS+20]\n", + "ts *= SCALING_FRACTION / np.std(ts)\n", + "ts -= np.mean(ts)\n", + "\n", + "t2s_ts = (ts * MEAN_T2S) + MEAN_T2S\n", + "s0_ts = (ts * MEAN_S0) + MEAN_S0" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Plot BOLD signal decay" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fullcurve_signal = pred_signal(FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), [30])\n", + "\n", + "fig, ax = plt.subplots(figsize=(14, 7.5))\n", + "ax.plot(\n", + " FULLCURVE_TES, \n", + " fullcurve_signal[:, 0], \n", + " alpha=0.5,\n", + " color=\"black\",\n", + ")\n", + "\n", + "ax.set_ylabel(\"Recorded BOLD signal\", fontsize=24)\n", + "ax.set_xlabel(\"Echo Time (ms)\", fontsize=24)\n", + "ax.set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)\n", + "ax.set_xlim(0, np.max(FULLCURVE_TES))\n", + "ax.tick_params(axis=\"both\", which=\"major\", labelsize=14)\n", + "fig.tight_layout()\n", + "\n", + "# save frame\n", + "fig.savefig(op.join(OUT_DIR, \"physics_signal_decay.png\"))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Plot BOLD signal decay and BOLD contrast" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fullcurve_signal_active = pred_signal(FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), [40])\n", + "fullcurve_signal_inactive = pred_signal(FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), [20])\n", + "\n", + "fig, ax = plt.subplots(figsize=(14, 7.5))\n", + "ax.plot(\n", + " FULLCURVE_TES,\n", + " fullcurve_signal_active[:, 0],\n", + " alpha=0.5,\n", + " color=\"red\",\n", + " label=\"High Activity\",\n", + ")\n", + "ax.plot(\n", + " FULLCURVE_TES,\n", + " fullcurve_signal_inactive[:, 0],\n", + " alpha=0.5,\n", + " color=\"blue\",\n", + " label=\"Low Activity\",\n", + ")\n", + "\n", + "ax.set_ylabel(\"Recorded BOLD signal\", fontsize=24)\n", + "ax.set_xlabel(\"Echo Time (ms)\", fontsize=24)\n", + "ax.set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)\n", + "ax.set_xlim(0, np.max(FULLCURVE_TES))\n", + "ax.legend(fontsize=20)\n", + "ax.tick_params(axis=\"both\", which=\"major\", labelsize=14)\n", + "fig.tight_layout()\n", + "\n", + "# save frame\n", + "fig.savefig(op.join(OUT_DIR, \"physics_signal_decay_activity.png\"))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Plot single-echo data resulting from $S_{0}$ and $T_{2}^{*}$ fluctuations\n", + "\n", + "This shows how fMRI data fluctuates over time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fullcurve_signal = pred_signal(FULLCURVE_TES, s0_ts, t2s_ts)\n", + "singleecho_signal = fullcurve_signal[SINGLEECHO_TE, :]\n", + "\n", + "out_file = op.join(OUT_DIR, \"fluctuations_single-echo.gif\")\n", + "if op.isfile(out_file):\n", + " os.remove(out_file)\n", + "\n", + "filenames = []\n", + "\n", + "for i_vol in range(N_VOLS):\n", + " filename = f\"se_{i_vol}.png\"\n", + " fig, axes = plt.subplots(\n", + " nrows=2, \n", + " figsize=(14, 10),\n", + " gridspec_kw={\"height_ratios\": [1, 3]}\n", + " )\n", + " \n", + " axes[0].plot(\n", + " singleecho_signal[0, :], \n", + " color=\"black\",\n", + " zorder=0,\n", + " )\n", + " axes[0].scatter(\n", + " i_vol, \n", + " singleecho_signal[:, i_vol], \n", + " color=\"orange\",\n", + " s=150,\n", + " label=\"Single-Echo Signal\",\n", + " zorder=1,\n", + " )\n", + " axes[0].set_ylabel(\"Signal\", fontsize=24)\n", + " axes[0].set_xlabel(\"Volume\", fontsize=24)\n", + " axes[0].set_xlim(0, N_VOLS - 1)\n", + " axes[0].tick_params(axis=\"both\", which=\"major\", labelsize=14)\n", + " \n", + " axes[1].scatter(\n", + " SINGLEECHO_TE, \n", + " singleecho_signal[:, i_vol], \n", + " color=\"orange\",\n", + " s=150,\n", + " alpha=1.,\n", + " label=\"Single-Echo Signal\",\n", + " )\n", + "\n", + " axes[1].set_ylabel(\"Signal\", fontsize=24)\n", + " axes[1].set_xlabel(\"Echo Time (ms)\", fontsize=24)\n", + " axes[1].set_xticks(MULTIECHO_TES)\n", + " axes[1].set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)\n", + " axes[1].set_xlim(0, np.max(FULLCURVE_TES))\n", + " axes[1].tick_params(axis=\"both\", which=\"major\", labelsize=14)\n", + " fig.tight_layout()\n", + " \n", + " # save frame\n", + " fig.savefig(filename)\n", + " plt.close(fig)\n", + " filenames.append(filename)\n", + "\n", + "# build gif\n", + "with imageio.get_writer(out_file, mode=\"I\") as writer:\n", + " for filename in filenames:\n", + " image = imageio.imread(filename)\n", + " writer.append_data(image)\n", + "\n", + "# Remove files\n", + "for filename in filenames:\n", + " os.remove(filename)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "with open(out_file, \"rb\") as file:\n", + " display(Image(file.read(), width=600))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Plot single-echo data and the curve resulting from $S_{0}$ and $T_{2}^{*}$ fluctuations\n", + "\n", + "This shows how single-echo data is a sample from a signal decay curve." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fullcurve_signal = pred_signal(FULLCURVE_TES, s0_ts, t2s_ts)\n", + "singleecho_signal = fullcurve_signal[SINGLEECHO_TE, :]\n", + "\n", + "out_file = op.join(OUT_DIR, \"fluctuations_single-echo_with_curve.gif\")\n", + "if op.isfile(out_file):\n", + " os.remove(out_file)\n", + "\n", + "filenames = []\n", + "\n", + "for i_vol in range(N_VOLS):\n", + " filename = f\"se_{i_vol}.png\"\n", + " fig, axes = plt.subplots(\n", + " nrows=2, \n", + " figsize=(14, 10),\n", + " gridspec_kw={\"height_ratios\": [1, 3]}\n", + " )\n", + " \n", + " axes[0].plot(\n", + " singleecho_signal[0, :], \n", + " color=\"black\", \n", + " zorder=0\n", + " )\n", + " axes[0].scatter(\n", + " i_vol, \n", + " singleecho_signal[:, i_vol], \n", + " color=\"orange\",\n", + " s=150,\n", + " label=\"Single-Echo Signal\",\n", + " zorder=1,\n", + " )\n", + " axes[0].set_ylabel(\"Signal\", fontsize=24)\n", + " axes[0].set_xlabel(\"Volume\", fontsize=24)\n", + " axes[0].set_xlim(0, N_VOLS - 1)\n", + " axes[0].tick_params(axis=\"both\", which=\"major\", labelsize=14)\n", + " \n", + " axes[1].plot(\n", + " FULLCURVE_TES, \n", + " fullcurve_signal[:, i_vol], \n", + " alpha=0.5,\n", + " color=\"black\",\n", + " zorder=0,\n", + " )\n", + " axes[1].scatter(\n", + " SINGLEECHO_TE, \n", + " singleecho_signal[:, i_vol], \n", + " color=\"orange\",\n", + " s=150,\n", + " alpha=1.,\n", + " label=\"Single-Echo Signal\",\n", + " zorder=1,\n", + " )\n", + "\n", + " axes[1].set_ylabel(\"Signal\", fontsize=24)\n", + " axes[1].set_xlabel(\"Echo Time (ms)\", fontsize=24)\n", + " axes[1].set_xticks(MULTIECHO_TES)\n", + " axes[1].set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)\n", + " axes[1].set_xlim(0, np.max(FULLCURVE_TES))\n", + " axes[1].tick_params(axis=\"both\", which=\"major\", labelsize=14)\n", + " fig.tight_layout()\n", + " \n", + " # save frame\n", + " fig.savefig(filename)\n", + " plt.close(fig)\n", + " filenames.append(filename)\n", + "\n", + "# build gif\n", + "with imageio.get_writer(out_file, mode=\"I\") as writer:\n", + " for filename in filenames:\n", + " image = imageio.imread(filename)\n", + " writer.append_data(image)\n", + "\n", + "# Remove files\n", + "for filename in filenames:\n", + " os.remove(filename)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "with open(out_file, \"rb\") as file:\n", + " display(Image(file.read(), width=600))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Plot single-echo data, the curve, and the $S_{0}$ and $T_{2}^{*}$ values resulting from $S_{0}$ and $T_{2}^{*}$ fluctuations\n", + "\n", + "This shows how changes in fMRI data can be driven by both S0 and T2* fluctuations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fullcurve_signal = pred_signal(FULLCURVE_TES, s0_ts, t2s_ts)\n", + "singleecho_signal = fullcurve_signal[SINGLEECHO_TE, :]\n", + "\n", + "out_file = op.join(OUT_DIR, \"fluctuations_single-echo_with_curve_and_t2s_s0.gif\")\n", + "if op.isfile(out_file):\n", + " os.remove(out_file)\n", + "\n", + "filenames = []\n", + "\n", + "for i_vol in range(N_VOLS):\n", + " filename = f\"se_{i_vol}.png\"\n", + " fig, axes = plt.subplots(\n", + " nrows=2, \n", + " figsize=(14, 10),\n", + " gridspec_kw={\"height_ratios\": [1, 3]}\n", + " )\n", + " \n", + " axes[0].plot(\n", + " singleecho_signal[0, :], \n", + " color=\"black\", \n", + " zorder=0,\n", + " )\n", + " axes[0].scatter(\n", + " i_vol, \n", + " singleecho_signal[:, i_vol], \n", + " color=\"orange\",\n", + " s=150,\n", + " label=\"Single-Echo Signal\",\n", + " zorder=1,\n", + " )\n", + " axes[0].set_ylabel(\"Signal\", fontsize=24)\n", + " axes[0].set_xlabel(\"Volume\", fontsize=24)\n", + " axes[0].set_xlim(0, N_VOLS - 1)\n", + " axes[0].tick_params(axis=\"both\", which=\"major\", labelsize=14)\n", + " \n", + " axes[1].scatter(\n", + " SINGLEECHO_TE, \n", + " singleecho_signal[:, i_vol], \n", + " color=\"orange\",\n", + " s=150,\n", + " alpha=1.,\n", + " label=\"Single-Echo Signal\",\n", + " zorder=3,\n", + " )\n", + " \n", + " axes[1].plot(\n", + " FULLCURVE_TES, \n", + " fullcurve_signal[:, i_vol], \n", + " alpha=0.5,\n", + " color=\"black\",\n", + " zorder=0,\n", + " )\n", + "\n", + " t2s_value = pred_signal(\n", + " np.array([t2s_ts[i_vol]]), \n", + " np.array([s0_ts[i_vol]]), \n", + " np.array([t2s_ts[i_vol]])\n", + " )[0]\n", + " axes[1].scatter(\n", + " t2s_ts[i_vol], \n", + " t2s_value, \n", + " marker=\"*\", \n", + " color=\"blue\", \n", + " s=500,\n", + " alpha=0.5,\n", + " label=\"$T_{2}^{*}$\",\n", + " zorder=1,\n", + " )\n", + " axes[1].scatter(\n", + " 0,\n", + " s0_ts[i_vol],\n", + " marker=\"*\", \n", + " color=\"red\", \n", + " s=500,\n", + " alpha=0.5,\n", + " label=\"$S_{0}$\",\n", + " clip_on=False,\n", + " zorder=2,\n", + " )\n", + " \n", + " axes[1].legend(loc=\"upper right\", fontsize=20)\n", + "\n", + " axes[1].set_ylabel(\"Signal\", fontsize=24)\n", + " axes[1].set_xlabel(\"Echo Time (ms)\", fontsize=24)\n", + " axes[1].set_xticks(MULTIECHO_TES)\n", + " axes[1].set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)\n", + " axes[1].set_xlim(0, np.max(FULLCURVE_TES))\n", + " axes[1].tick_params(axis=\"both\", which=\"major\", labelsize=14)\n", + " fig.tight_layout()\n", + "\n", + " # save frame\n", + " fig.savefig(filename)\n", + " plt.close(fig)\n", + " filenames.append(filename)\n", + "\n", + "# build gif\n", + "with imageio.get_writer(out_file, mode=\"I\") as writer:\n", + " for filename in filenames:\n", + " image = imageio.imread(filename)\n", + " writer.append_data(image)\n", + "\n", + "# Remove files\n", + "for filename in filenames:\n", + " os.remove(filename)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "with open(out_file, \"rb\") as file:\n", + " display(Image(file.read(), width=600))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Plot $S_{0}$ and $T_{2}^{*}$ fluctuations\n", + "\n", + "This shows how fluctuations in S0 and T2* produce different patterns in the full signal decay curves." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "s0based_fullcurve_signal = pred_signal(FULLCURVE_TES, s0_ts, np.full(N_VOLS, MEAN_T2S))\n", + "t2sbased_fullcurve_signal = pred_signal(FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), t2s_ts)\n", + "\n", + "out_file = op.join(OUT_DIR, \"fluctuations_t2s_s0.gif\")\n", + "if op.isfile(out_file):\n", + " os.remove(out_file)\n", + "\n", + "filenames = []\n", + "\n", + "for i_vol in range(N_VOLS):\n", + " filename = f\"total_{i_vol}.png\"\n", + " fig, axes = plt.subplots(\n", + " nrows=2, \n", + " figsize=(14, 10),\n", + " gridspec_kw={\"height_ratios\": [1, 3]}\n", + " )\n", + " \n", + " axes[0].plot(ts, color=\"black\")\n", + " axes[0].scatter(\n", + " i_vol, \n", + " ts[i_vol], \n", + " color=\"purple\",\n", + " s=150,\n", + " )\n", + " axes[0].set_ylabel(\"$S_0$/$T_{2}^{*}$\", fontsize=24)\n", + " axes[0].set_xlabel(\"Volume\", fontsize=24)\n", + " axes[0].set_xlim(0, N_VOLS - 1)\n", + " axes[0].tick_params(axis=\"both\", which=\"major\", labelsize=14)\n", + " \n", + " axes[1].plot(\n", + " FULLCURVE_TES, \n", + " s0based_fullcurve_signal[:, i_vol], \n", + " color=\"red\",\n", + " alpha=0.5,\n", + " label=\"$S_0$-Driven Decay Curve\",\n", + " )\n", + " axes[1].plot(\n", + " FULLCURVE_TES, \n", + " t2sbased_fullcurve_signal[:, i_vol], \n", + " color=\"blue\",\n", + " alpha=0.5,\n", + " label=\"$T_{2}^{*}$-Driven Decay Curve\",\n", + " )\n", + " axes[1].legend(loc=\"upper right\", fontsize=20)\n", + "\n", + " axes[1].set_ylabel(\"Signal\", fontsize=24)\n", + " axes[1].set_xlabel(\"Echo Time (ms)\", fontsize=24)\n", + " axes[1].set_xticks(MULTIECHO_TES)\n", + " axes[1].set_ylim(0, np.ceil(np.max(s0based_fullcurve_signal) / 1000) * 1000)\n", + " axes[1].set_xlim(0, np.max(FULLCURVE_TES))\n", + " axes[1].tick_params(axis=\"both\", which=\"major\", labelsize=14)\n", + " fig.tight_layout()\n", + " \n", + " # save frame\n", + " fig.savefig(filename)\n", + " plt.close(fig)\n", + " filenames.append(filename)\n", + "\n", + "# build gif\n", + "with imageio.get_writer(out_file, mode=\"I\") as writer:\n", + " for filename in filenames:\n", + " image = imageio.imread(filename)\n", + " writer.append_data(image)\n", + "\n", + "# Remove files\n", + "for filename in filenames:\n", + " os.remove(filename)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "with open(out_file, \"rb\") as file:\n", + " display(Image(file.read(), width=600))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Plot $S_{0}$ and $T_{2}^{*}$ fluctuations and resulting single-echo data\n", + "\n", + "This shows how single-echo data, on its own, cannot distinguish between S0 and T2* fluctuations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "s0based_fullcurve_signal = pred_signal(FULLCURVE_TES, s0_ts, np.full(N_VOLS, MEAN_T2S))\n", + "t2sbased_fullcurve_signal = pred_signal(FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), t2s_ts)\n", + "s0based_singleecho_signal = s0based_fullcurve_signal[SINGLEECHO_TE, :]\n", + "t2sbased_singleecho_signal = t2sbased_fullcurve_signal[SINGLEECHO_TE, :]\n", + "\n", + "out_file = op.join(OUT_DIR, \"fluctuations_t2s_s0_single-echo.gif\")\n", + "if op.isfile(out_file):\n", + " os.remove(out_file)\n", + "\n", + "filenames = []\n", + "\n", + "for i_vol in range(N_VOLS):\n", + " filename = f\"total_{i_vol}.png\"\n", + " fig, axes = plt.subplots(\n", + " nrows=2, \n", + " figsize=(14, 10),\n", + " gridspec_kw={\"height_ratios\": [1, 3]}\n", + " )\n", + " \n", + " axes[0].plot(\n", + " ts, \n", + " color=\"black\",\n", + " zorder=0,\n", + " )\n", + " axes[0].scatter(\n", + " i_vol, \n", + " ts[i_vol], \n", + " color=\"purple\",\n", + " s=150,\n", + " zorder=1,\n", + " )\n", + " axes[0].set_ylabel(\"$S_0$/$T_{2}^{*}$\", fontsize=24)\n", + " axes[0].set_xlabel(\"Volume\", fontsize=24)\n", + " axes[0].set_xlim(0, N_VOLS - 1)\n", + " axes[0].tick_params(axis=\"both\", which=\"major\", labelsize=14)\n", + " \n", + " axes[1].plot(\n", + " FULLCURVE_TES, \n", + " s0based_fullcurve_signal[:, i_vol], \n", + " color=\"red\",\n", + " alpha=0.5,\n", + " label=\"$S_0$-Driven Decay Curve\",\n", + " zorder=0,\n", + " )\n", + " axes[1].plot(\n", + " FULLCURVE_TES, \n", + " t2sbased_fullcurve_signal[:, i_vol], \n", + " color=\"blue\",\n", + " alpha=0.5,\n", + " label=\"$T_{2}^{*}$-Driven Decay Curve\",\n", + " zorder=1,\n", + " )\n", + " axes[1].scatter(\n", + " SINGLEECHO_TE, \n", + " s0based_singleecho_signal[:, i_vol], \n", + " color=\"red\",\n", + " s=150,\n", + " alpha=1.,\n", + " label=\"$S_0$-Driven Single-Echo Signal\",\n", + " zorder=2,\n", + " )\n", + " axes[1].scatter(\n", + " SINGLEECHO_TE, \n", + " t2sbased_singleecho_signal[:, i_vol], \n", + " color=\"blue\",\n", + " s=150,\n", + " alpha=1.,\n", + " label=\"$T_{2}^{*}$-Driven Single-Echo Signal\",\n", + " zorder=3,\n", + " )\n", + " axes[1].legend(loc=\"upper right\", fontsize=20)\n", + "\n", + " axes[1].set_ylabel(\"Signal\", fontsize=24)\n", + " axes[1].set_xlabel(\"Echo Time (ms)\", fontsize=24)\n", + " axes[1].set_xticks(MULTIECHO_TES)\n", + " axes[1].set_ylim(0, np.ceil(np.max(s0based_fullcurve_signal) / 1000) * 1000)\n", + " axes[1].set_xlim(0, np.max(FULLCURVE_TES))\n", + " axes[1].tick_params(axis=\"both\", which=\"major\", labelsize=14)\n", + " fig.tight_layout()\n", + " \n", + " # save frame\n", + " fig.savefig(filename)\n", + " plt.close(fig)\n", + " filenames.append(filename)\n", + "\n", + "# build gif\n", + "with imageio.get_writer(out_file, mode=\"I\") as writer:\n", + " for filename in filenames:\n", + " image = imageio.imread(filename)\n", + " writer.append_data(image)\n", + "\n", + "# Remove files\n", + "for filename in filenames:\n", + " os.remove(filename)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "with open(out_file, \"rb\") as file:\n", + " display(Image(file.read(), width=600))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Plot $S_{0}$ and $T_{2}^{*}$ fluctuations and resulting multi-echo data\n", + "\n", + "This shows how S0 and T2* fluctuations produce different patterns in multi-echo data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "s0based_fullcurve_signal = pred_signal(FULLCURVE_TES, s0_ts, np.full(N_VOLS, MEAN_T2S))\n", + "t2sbased_fullcurve_signal = pred_signal(FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), t2s_ts)\n", + "s0based_multiecho_signal = s0based_fullcurve_signal[MULTIECHO_TES, :]\n", + "t2sbased_multiecho_signal = t2sbased_fullcurve_signal[MULTIECHO_TES, :]\n", + "\n", + "out_file = op.join(OUT_DIR, \"fluctuations_t2s_s0_multi-echo.gif\")\n", + "if op.isfile(out_file):\n", + " os.remove(out_file)\n", + "\n", + "filenames = []\n", + "\n", + "for i_vol in range(N_VOLS):\n", + " filename = f\"total_{i_vol}.png\"\n", + " fig, axes = plt.subplots(\n", + " nrows=2, \n", + " figsize=(14, 10),\n", + " gridspec_kw={\"height_ratios\": [1, 3]}\n", + " )\n", + " \n", + " axes[0].plot(\n", + " ts, \n", + " color=\"black\",\n", + " zorder=0,\n", + " )\n", + " axes[0].scatter(\n", + " i_vol, \n", + " ts[i_vol], \n", + " color=\"purple\",\n", + " s=150,\n", + " zorder=1,\n", + " )\n", + " axes[0].set_ylabel(\"$S_0$/$T_{2}^{*}$\", fontsize=24)\n", + " axes[0].set_xlabel(\"Volume\", fontsize=24)\n", + " axes[0].set_xlim(0, N_VOLS - 1)\n", + " axes[0].tick_params(axis=\"both\", which=\"major\", labelsize=14)\n", + " \n", + " axes[1].plot(\n", + " FULLCURVE_TES, \n", + " s0based_fullcurve_signal[:, i_vol], \n", + " color=\"red\",\n", + " alpha=0.5,\n", + " label=\"$S_0$-Driven Decay Curve\",\n", + " zorder=0,\n", + " )\n", + " axes[1].plot(\n", + " FULLCURVE_TES, \n", + " t2sbased_fullcurve_signal[:, i_vol], \n", + " color=\"blue\",\n", + " alpha=0.5,\n", + " label=\"$T_{2}^{*}$-Driven Decay Curve\",\n", + " zorder=1,\n", + " )\n", + " axes[1].scatter(\n", + " MULTIECHO_TES, \n", + " s0based_multiecho_signal[:, i_vol], \n", + " color=\"red\",\n", + " s=150,\n", + " alpha=1.,\n", + " label=\"$S_0$-Driven Multi-Echo Signal\",\n", + " zorder=2,\n", + " )\n", + " axes[1].scatter(\n", + " MULTIECHO_TES, \n", + " t2sbased_multiecho_signal[:, i_vol], \n", + " color=\"blue\",\n", + " s=150,\n", + " alpha=1.,\n", + " label=\"$T_{2}^{*}$-Driven Multi-Echo Signal\",\n", + " zorder=3,\n", + " )\n", + " axes[1].legend(loc=\"upper right\", fontsize=20)\n", + "\n", + " axes[1].set_ylabel(\"Signal\", fontsize=24)\n", + " axes[1].set_xlabel(\"Echo Time (ms)\", fontsize=24)\n", + " axes[1].set_xticks(MULTIECHO_TES)\n", + " axes[1].set_ylim(0, np.ceil(np.max(s0based_fullcurve_signal) / 1000) * 1000)\n", + " axes[1].set_xlim(0, np.max(FULLCURVE_TES))\n", + " axes[1].tick_params(axis=\"both\", which=\"major\", labelsize=14)\n", + " fig.tight_layout()\n", + " \n", + " # save frame\n", + " fig.savefig(filename)\n", + " plt.close(fig)\n", + " filenames.append(filename)\n", + "\n", + "# build gif\n", + "with imageio.get_writer(out_file, mode=\"I\") as writer:\n", + " for filename in filenames:\n", + " image = imageio.imread(filename)\n", + " writer.append_data(image)\n", + "\n", + "# Remove files\n", + "for filename in filenames:\n", + " os.remove(filename)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "with open(out_file, \"rb\") as file:\n", + " display(Image(file.read(), width=600))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot $T_{2}^{*}$ against BOLD signal from single-echo data (TE=30ms)\n", + "\n", + "When the BOLD signal is driven entirely by T2* fluctuations, the BOLD signal is very similar to a scaled version of those T2* fluctuations, but there are small differences." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(16, 4))\n", + "\n", + "scalar = np.linalg.lstsq(\n", + " t2s_ts[:, None], \n", + " t2sbased_fullcurve_signal[SINGLEECHO_TE[0], :], \n", + " rcond=None\n", + ")[0]\n", + "ax.plot(\n", + " t2sbased_fullcurve_signal[SINGLEECHO_TE[0], :], \n", + " color=\"orange\", \n", + " label=f\"BOLD signal (TE={SINGLEECHO_TE[0]}ms)\",\n", + ")\n", + "ax.plot(\n", + " t2s_ts * scalar, \n", + " label=\"$T_{2}^{*}$ (scaled to signal)\", \n", + " linestyle=\"--\",\n", + ")\n", + "ax.set_xlim(0, N_VOLS - 1)\n", + "leg = ax.legend()\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot $S_{0}$ against BOLD signal from single-echo data (TE=30ms)\n", + "\n", + "When the BOLD signal is driven entirely by S0 fluctuations, the BOLD signal ends up being a scaled version of the S0 fluctuations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(16, 4))\n", + "\n", + "scalar = np.linalg.lstsq(\n", + " s0_ts[:, None], \n", + " s0based_fullcurve_signal[SINGLEECHO_TE[0], :], \n", + " rcond=None\n", + ")[0]\n", + "ax.plot(\n", + " s0based_fullcurve_signal[SINGLEECHO_TE[0], :],\n", + " color=\"orange\", \n", + " label=f\"BOLD signal (TE={SINGLEECHO_TE[0]}ms)\",\n", + ")\n", + "ax.plot(\n", + " s0_ts * scalar, \n", + " label=\"S0 (scaled to signal)\", \n", + " linestyle=\"--\",\n", + ")\n", + "ax.set_xlim(0, N_VOLS - 1)\n", + "leg = ax.legend()\n", + "fig.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}