From 4ec69159ebaea27969ad8cbd6af99cb7da2becc1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Eneko=20Uru=C3=B1uela?= Date: Mon, 17 Jan 2022 18:52:14 +0100 Subject: [PATCH 1/4] Set up docs --- docs/Makefile | 23 +++++ docs/_templates/function.rst | 12 +++ docs/_templates/module.rst | 40 ++++++++ docs/conf.py | 176 ++++++++++++++++++++++++++++++++++ docs/index.rst | 20 ++++ docs/make.bat | 36 +++++++ docs/sphinxext/github_link.py | 85 ++++++++++++++++ 7 files changed, 392 insertions(+) create mode 100644 docs/Makefile create mode 100644 docs/_templates/function.rst create mode 100644 docs/_templates/module.rst create mode 100644 docs/conf.py create mode 100644 docs/index.rst create mode 100644 docs/make.bat create mode 100644 docs/sphinxext/github_link.py diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000..50c0f47 --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,23 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line. +SPHINXOPTS = +SPHINXBUILD = python -msphinx +SPHINXPROJ = maPCA +SOURCEDIR = . +BUILDDIR = _build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +clean: + rm -r _build generated + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) \ No newline at end of file diff --git a/docs/_templates/function.rst b/docs/_templates/function.rst new file mode 100644 index 0000000..ed57231 --- /dev/null +++ b/docs/_templates/function.rst @@ -0,0 +1,12 @@ +:mod:`{{module}}`.{{objname}} +{{ underline }}==================== + +.. currentmodule:: {{ module }} + +.. autofunction:: {{ objname }} + +.. .. include:: {{module}}.{{objname}}.examples + +.. raw:: html + +
\ No newline at end of file diff --git a/docs/_templates/module.rst b/docs/_templates/module.rst new file mode 100644 index 0000000..941dfd6 --- /dev/null +++ b/docs/_templates/module.rst @@ -0,0 +1,40 @@ +{{ fullname }} +{{ underline }} + +.. automodule:: {{ fullname }} + + {% block functions %} + {% if functions %} + .. rubric:: Functions + + .. autosummary:: + :toctree: + {% for item in functions %} + {{ item }} + {%- endfor %} + {% endif %} + {% endblock %} + + {% block classes %} + {% if classes %} + .. rubric:: Classes + + .. autosummary:: + :toctree: + {% for item in classes %} + {{ item }} + {%- endfor %} + {% endif %} + {% endblock %} + + {% block exceptions %} + {% if exceptions %} + .. rubric:: Exceptions + + .. autosummary:: + :toctree: + {% for item in exceptions %} + {{ item }} + {%- endfor %} + {% endif %} + {% endblock %} \ No newline at end of file diff --git a/docs/conf.py b/docs/conf.py new file mode 100644 index 0000000..0ebf632 --- /dev/null +++ b/docs/conf.py @@ -0,0 +1,176 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# +# mapca documentation build configuration file, created by +# sphinx-quickstart +# +# This file is execfile()d with the current directory set to its +# containing dir. +# +# Note that not all possible configuration values are present in this +# autogenerated file. +# +# All configuration values have a default; values that are commented out +# serve to show the default. + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +import os +import sys + +sys.path.insert(0, os.path.abspath("sphinxext")) +sys.path.insert(0, os.path.abspath(os.path.pardir)) + +from github_link import make_linkcode_resolve + +# -- General configuration ------------------------------------------------ + +# If your documentation needs a minimal Sphinx version, state it here. + +# needs_sphinx = '1.0' + +# generate autosummary even if no references +autosummary_generate = True +add_module_names = False + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + "matplotlib.sphinxext.plot_directive", + "sphinx.ext.autodoc", + "sphinx.ext.autosummary", + "sphinx.ext.doctest", + "sphinx.ext.ifconfig", + "sphinx.ext.intersphinx", + "sphinx.ext.linkcode", + "sphinx.ext.napoleon", + "sphinx.ext.todo", + "sphinxarg.ext", +] + +from distutils.version import LooseVersion + +import sphinx + +if LooseVersion(sphinx.__version__) < LooseVersion("1.4"): + extensions.append("sphinx.ext.pngmath") +else: + extensions.append("sphinx.ext.imgmath") + +# Add any paths that contain templates here, relative to this directory. +templates_path = ["_templates"] + +# source_suffix = ['.rst', '.md'] +source_suffix = ".rst" + +# The master toctree document. +master_doc = "index" + +# General information about the project. +from datetime import datetime # access current time and date + +project = "maPCA" +copyright = "2021-" + datetime.today().strftime("%Y") + ", maPCA developers" +author = "maPCA developers" + +# The version info for the project you're documenting, acts as replacement for +# |version| and |release|, also used in various other places throughout the +# built documents. +# +# The short X.Y version. +import mapca + +version = mapca.__version__ +# The full version, including alpha/beta/rc tags. +release = mapca.__version__ + +# The language for content autogenerated by Sphinx. Refer to documentation +# for a list of supported languages. +# +# This is also used if you do content translation via gettext catalogs. +# Usually you set "language" from the command line for these cases. +language = None + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This patterns also effect to html_static_path and html_extra_path +exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"] + +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = "sphinx" + +# If true, `todo` and `todoList` produce output, else they produce nothing. +todo_include_todos = False + +# ----------------------------------------------------------------------------- +# Napoleon settings +# ----------------------------------------------------------------------------- +napoleon_google_docstring = True +napoleon_numpy_docstring = True +napoleon_include_init_with_doc = False +napoleon_include_private_with_doc = False +napoleon_include_special_with_doc = True +napoleon_use_admonition_for_examples = False +napoleon_use_admonition_for_notes = True +napoleon_use_admonition_for_references = False +napoleon_use_ivar = True +napoleon_use_param = False +napoleon_use_keyword = True +napoleon_use_rtype = False + +# -- Options for HTML output ---------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +# installing theme package +import sphinx_rtd_theme + +html_theme = "sphinx_rtd_theme" + +# Theme options are theme-specific and customize the look and feel of a theme +# further. For a list of options available for each theme, see the +# documentation. + +html_theme_options = { + "includehidden": False, +} + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ["_static"] + + +# https://github.com/rtfd/sphinx_rtd_theme/issues/117 +def setup(app): + app.add_css_file("theme_overrides.css") + app.add_js_file("https://cdn.rawgit.com/chrisfilo/zenodo.js/v0.1/zenodo.js") + + +html_favicon = "_static/mapca_favicon.png" + + +# -- Options for HTMLHelp output ------------------------------------------ + +# Output file base name for HTML help builder. +htmlhelp_basename = "mapcadoc" + +# The following is used by sphinx.ext.linkcode to provide links to github +linkcode_resolve = make_linkcode_resolve( + "mapca", "https://github.com/me-ica/" "mapca/blob/{revision}/" "{package}/{path}#L{lineno}" +) + +# Example configuration for intersphinx: refer to the Python standard library. +intersphinx_mapping = { + "python": ("https://docs.python.org/3/", None), + "numpy": ("https://numpy.org/doc/stable/", None), + "scipy": ("https://docs.scipy.org/doc/scipy/reference/", None), + "matplotlib": ("https://matplotlib.org/", None), + "nibabel": ("https://nipy.org/nibabel/", None), + "pandas": ("https://pandas.pydata.org/pandas-docs/stable/", None), + "nilearn": ("http://nilearn.github.io/", None), +} diff --git a/docs/index.rst b/docs/index.rst new file mode 100644 index 0000000..cd1ce8e --- /dev/null +++ b/docs/index.rst @@ -0,0 +1,20 @@ +.. maPCA documentation master file, created by + sphinx-quickstart on Mon Jan 17 18:00:49 2022. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +Welcome to maPCA's documentation! +================================= + +.. toctree:: + :maxdepth: 2 + :caption: Contents: + + + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/docs/make.bat b/docs/make.bat new file mode 100644 index 0000000..44b76d2 --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,36 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=python -msphinx +) +set SOURCEDIR=. +set BUILDDIR=_build +set SPHINXPROJ=meica + +if "%1" == "" goto help + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The Sphinx module was not found. Make sure you have Sphinx installed, + echo.then set the SPHINXBUILD environment variable to point to the full + echo.path of the 'sphinx-build' executable. Alternatively you may add the + echo.Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.http://sphinx-doc.org/ + exit /b 1 +) + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% + +:end +popd \ No newline at end of file diff --git a/docs/sphinxext/github_link.py b/docs/sphinxext/github_link.py new file mode 100644 index 0000000..81df6d3 --- /dev/null +++ b/docs/sphinxext/github_link.py @@ -0,0 +1,85 @@ +""" +This script comes from scikit-learn: +https://github.com/scikit-learn/scikit-learn/blob/main/doc/sphinxext/github_link.py +""" +import inspect +import os +import subprocess +import sys +from functools import partial +from operator import attrgetter + +REVISION_CMD = "git rev-parse --short HEAD" + + +def _get_git_revision(): + try: + revision = subprocess.check_output(REVISION_CMD.split()).strip() + except (subprocess.CalledProcessError, OSError): + print("Failed to execute git to get revision") + return None + return revision.decode("utf-8") + + +def _linkcode_resolve(domain, info, package, url_fmt, revision): + """Determine a link to online source for a class/method/function + + This is called by sphinx.ext.linkcode + + An example with a long-untouched module that everyone has + >>> _linkcode_resolve('py', {'module': 'tty', + ... 'fullname': 'setraw'}, + ... package='tty', + ... url_fmt='http://hg.python.org/cpython/file/' + ... '{revision}/Lib/{package}/{path}#L{lineno}', + ... revision='xxxx') + 'http://hg.python.org/cpython/file/xxxx/Lib/tty/tty.py#L18' + """ + + if revision is None: + return + if domain not in ("py", "pyx"): + return + if not info.get("module") or not info.get("fullname"): + return + + class_name = info["fullname"].split(".")[0] + if type(class_name) != str: + # Python 2 only + class_name = class_name.encode("utf-8") + module = __import__(info["module"], fromlist=[class_name]) + obj = attrgetter(info["fullname"])(module) + + try: + fn = inspect.getsourcefile(obj) + except Exception: + fn = None + if not fn: + try: + fn = inspect.getsourcefile(sys.modules[obj.__module__]) + except Exception: + fn = None + if not fn: + return + + fn = os.path.relpath(fn, start=os.path.dirname(__import__(package).__file__)) + try: + lineno = inspect.getsourcelines(obj)[1] + except Exception: + lineno = "" + return url_fmt.format(revision=revision, package=package, path=fn, lineno=lineno) + + +def make_linkcode_resolve(package, url_fmt): + """Returns a linkcode_resolve function for the given URL format + + revision is a git commit reference (hash or name) + + package is the name of the root module of the package + + url_fmt is along the lines of ('https://github.com/USER/PROJECT/' + 'blob/{revision}/{package}/' + '{path}#L{lineno}') + """ + revision = _get_git_revision() + return partial(_linkcode_resolve, revision=revision, package=package, url_fmt=url_fmt) From 1e8a8db9858c4923aa77c4b12058ef1c2a6b9ed8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Eneko=20Uru=C3=B1uela?= Date: Tue, 25 Jan 2022 17:11:01 +0100 Subject: [PATCH 2/4] First, silly draft of notebook --- .gitignore | 3 + docs/Makefile | 23 -- docs/_templates/function.rst | 12 - docs/_templates/module.rst | 40 --- docs/conf.py | 176 ---------- docs/index.rst | 20 -- docs/make.bat | 36 -- docs/sphinxext/github_link.py | 85 ----- notebooks/pipeline.ipynb | 632 ++++++++++++++++++++++++++++++++++ 9 files changed, 635 insertions(+), 392 deletions(-) delete mode 100644 docs/Makefile delete mode 100644 docs/_templates/function.rst delete mode 100644 docs/_templates/module.rst delete mode 100644 docs/conf.py delete mode 100644 docs/index.rst delete mode 100644 docs/make.bat delete mode 100644 docs/sphinxext/github_link.py create mode 100644 notebooks/pipeline.ipynb diff --git a/.gitignore b/.gitignore index b6e4761..11f1e4f 100644 --- a/.gitignore +++ b/.gitignore @@ -127,3 +127,6 @@ dmypy.json # Pyre type checker .pyre/ + +# Nilearn cache +nilearn_cache \ No newline at end of file diff --git a/docs/Makefile b/docs/Makefile deleted file mode 100644 index 50c0f47..0000000 --- a/docs/Makefile +++ /dev/null @@ -1,23 +0,0 @@ -# Minimal makefile for Sphinx documentation -# - -# You can set these variables from the command line. -SPHINXOPTS = -SPHINXBUILD = python -msphinx -SPHINXPROJ = maPCA -SOURCEDIR = . -BUILDDIR = _build - -# Put it first so that "make" without argument is like "make help". -help: - @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) - -clean: - rm -r _build generated - -.PHONY: help Makefile - -# Catch-all target: route all unknown targets to Sphinx using the new -# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). -%: Makefile - @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) \ No newline at end of file diff --git a/docs/_templates/function.rst b/docs/_templates/function.rst deleted file mode 100644 index ed57231..0000000 --- a/docs/_templates/function.rst +++ /dev/null @@ -1,12 +0,0 @@ -:mod:`{{module}}`.{{objname}} -{{ underline }}==================== - -.. currentmodule:: {{ module }} - -.. autofunction:: {{ objname }} - -.. .. include:: {{module}}.{{objname}}.examples - -.. raw:: html - -
\ No newline at end of file diff --git a/docs/_templates/module.rst b/docs/_templates/module.rst deleted file mode 100644 index 941dfd6..0000000 --- a/docs/_templates/module.rst +++ /dev/null @@ -1,40 +0,0 @@ -{{ fullname }} -{{ underline }} - -.. automodule:: {{ fullname }} - - {% block functions %} - {% if functions %} - .. rubric:: Functions - - .. autosummary:: - :toctree: - {% for item in functions %} - {{ item }} - {%- endfor %} - {% endif %} - {% endblock %} - - {% block classes %} - {% if classes %} - .. rubric:: Classes - - .. autosummary:: - :toctree: - {% for item in classes %} - {{ item }} - {%- endfor %} - {% endif %} - {% endblock %} - - {% block exceptions %} - {% if exceptions %} - .. rubric:: Exceptions - - .. autosummary:: - :toctree: - {% for item in exceptions %} - {{ item }} - {%- endfor %} - {% endif %} - {% endblock %} \ No newline at end of file diff --git a/docs/conf.py b/docs/conf.py deleted file mode 100644 index 0ebf632..0000000 --- a/docs/conf.py +++ /dev/null @@ -1,176 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -# -# mapca documentation build configuration file, created by -# sphinx-quickstart -# -# This file is execfile()d with the current directory set to its -# containing dir. -# -# Note that not all possible configuration values are present in this -# autogenerated file. -# -# All configuration values have a default; values that are commented out -# serve to show the default. - -# If extensions (or modules to document with autodoc) are in another directory, -# add these directories to sys.path here. If the directory is relative to the -# documentation root, use os.path.abspath to make it absolute, like shown here. -# -import os -import sys - -sys.path.insert(0, os.path.abspath("sphinxext")) -sys.path.insert(0, os.path.abspath(os.path.pardir)) - -from github_link import make_linkcode_resolve - -# -- General configuration ------------------------------------------------ - -# If your documentation needs a minimal Sphinx version, state it here. - -# needs_sphinx = '1.0' - -# generate autosummary even if no references -autosummary_generate = True -add_module_names = False - -# Add any Sphinx extension module names here, as strings. They can be -# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom -# ones. -extensions = [ - "matplotlib.sphinxext.plot_directive", - "sphinx.ext.autodoc", - "sphinx.ext.autosummary", - "sphinx.ext.doctest", - "sphinx.ext.ifconfig", - "sphinx.ext.intersphinx", - "sphinx.ext.linkcode", - "sphinx.ext.napoleon", - "sphinx.ext.todo", - "sphinxarg.ext", -] - -from distutils.version import LooseVersion - -import sphinx - -if LooseVersion(sphinx.__version__) < LooseVersion("1.4"): - extensions.append("sphinx.ext.pngmath") -else: - extensions.append("sphinx.ext.imgmath") - -# Add any paths that contain templates here, relative to this directory. -templates_path = ["_templates"] - -# source_suffix = ['.rst', '.md'] -source_suffix = ".rst" - -# The master toctree document. -master_doc = "index" - -# General information about the project. -from datetime import datetime # access current time and date - -project = "maPCA" -copyright = "2021-" + datetime.today().strftime("%Y") + ", maPCA developers" -author = "maPCA developers" - -# The version info for the project you're documenting, acts as replacement for -# |version| and |release|, also used in various other places throughout the -# built documents. -# -# The short X.Y version. -import mapca - -version = mapca.__version__ -# The full version, including alpha/beta/rc tags. -release = mapca.__version__ - -# The language for content autogenerated by Sphinx. Refer to documentation -# for a list of supported languages. -# -# This is also used if you do content translation via gettext catalogs. -# Usually you set "language" from the command line for these cases. -language = None - -# List of patterns, relative to source directory, that match files and -# directories to ignore when looking for source files. -# This patterns also effect to html_static_path and html_extra_path -exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"] - -# The name of the Pygments (syntax highlighting) style to use. -pygments_style = "sphinx" - -# If true, `todo` and `todoList` produce output, else they produce nothing. -todo_include_todos = False - -# ----------------------------------------------------------------------------- -# Napoleon settings -# ----------------------------------------------------------------------------- -napoleon_google_docstring = True -napoleon_numpy_docstring = True -napoleon_include_init_with_doc = False -napoleon_include_private_with_doc = False -napoleon_include_special_with_doc = True -napoleon_use_admonition_for_examples = False -napoleon_use_admonition_for_notes = True -napoleon_use_admonition_for_references = False -napoleon_use_ivar = True -napoleon_use_param = False -napoleon_use_keyword = True -napoleon_use_rtype = False - -# -- Options for HTML output ---------------------------------------------- - -# The theme to use for HTML and HTML Help pages. See the documentation for -# a list of builtin themes. -# -# installing theme package -import sphinx_rtd_theme - -html_theme = "sphinx_rtd_theme" - -# Theme options are theme-specific and customize the look and feel of a theme -# further. For a list of options available for each theme, see the -# documentation. - -html_theme_options = { - "includehidden": False, -} - -# Add any paths that contain custom static files (such as style sheets) here, -# relative to this directory. They are copied after the builtin static files, -# so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ["_static"] - - -# https://github.com/rtfd/sphinx_rtd_theme/issues/117 -def setup(app): - app.add_css_file("theme_overrides.css") - app.add_js_file("https://cdn.rawgit.com/chrisfilo/zenodo.js/v0.1/zenodo.js") - - -html_favicon = "_static/mapca_favicon.png" - - -# -- Options for HTMLHelp output ------------------------------------------ - -# Output file base name for HTML help builder. -htmlhelp_basename = "mapcadoc" - -# The following is used by sphinx.ext.linkcode to provide links to github -linkcode_resolve = make_linkcode_resolve( - "mapca", "https://github.com/me-ica/" "mapca/blob/{revision}/" "{package}/{path}#L{lineno}" -) - -# Example configuration for intersphinx: refer to the Python standard library. -intersphinx_mapping = { - "python": ("https://docs.python.org/3/", None), - "numpy": ("https://numpy.org/doc/stable/", None), - "scipy": ("https://docs.scipy.org/doc/scipy/reference/", None), - "matplotlib": ("https://matplotlib.org/", None), - "nibabel": ("https://nipy.org/nibabel/", None), - "pandas": ("https://pandas.pydata.org/pandas-docs/stable/", None), - "nilearn": ("http://nilearn.github.io/", None), -} diff --git a/docs/index.rst b/docs/index.rst deleted file mode 100644 index cd1ce8e..0000000 --- a/docs/index.rst +++ /dev/null @@ -1,20 +0,0 @@ -.. maPCA documentation master file, created by - sphinx-quickstart on Mon Jan 17 18:00:49 2022. - You can adapt this file completely to your liking, but it should at least - contain the root `toctree` directive. - -Welcome to maPCA's documentation! -================================= - -.. toctree:: - :maxdepth: 2 - :caption: Contents: - - - -Indices and tables -================== - -* :ref:`genindex` -* :ref:`modindex` -* :ref:`search` diff --git a/docs/make.bat b/docs/make.bat deleted file mode 100644 index 44b76d2..0000000 --- a/docs/make.bat +++ /dev/null @@ -1,36 +0,0 @@ -@ECHO OFF - -pushd %~dp0 - -REM Command file for Sphinx documentation - -if "%SPHINXBUILD%" == "" ( - set SPHINXBUILD=python -msphinx -) -set SOURCEDIR=. -set BUILDDIR=_build -set SPHINXPROJ=meica - -if "%1" == "" goto help - -%SPHINXBUILD% >NUL 2>NUL -if errorlevel 9009 ( - echo. - echo.The Sphinx module was not found. Make sure you have Sphinx installed, - echo.then set the SPHINXBUILD environment variable to point to the full - echo.path of the 'sphinx-build' executable. Alternatively you may add the - echo.Sphinx directory to PATH. - echo. - echo.If you don't have Sphinx installed, grab it from - echo.http://sphinx-doc.org/ - exit /b 1 -) - -%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% -goto end - -:help -%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% - -:end -popd \ No newline at end of file diff --git a/docs/sphinxext/github_link.py b/docs/sphinxext/github_link.py deleted file mode 100644 index 81df6d3..0000000 --- a/docs/sphinxext/github_link.py +++ /dev/null @@ -1,85 +0,0 @@ -""" -This script comes from scikit-learn: -https://github.com/scikit-learn/scikit-learn/blob/main/doc/sphinxext/github_link.py -""" -import inspect -import os -import subprocess -import sys -from functools import partial -from operator import attrgetter - -REVISION_CMD = "git rev-parse --short HEAD" - - -def _get_git_revision(): - try: - revision = subprocess.check_output(REVISION_CMD.split()).strip() - except (subprocess.CalledProcessError, OSError): - print("Failed to execute git to get revision") - return None - return revision.decode("utf-8") - - -def _linkcode_resolve(domain, info, package, url_fmt, revision): - """Determine a link to online source for a class/method/function - - This is called by sphinx.ext.linkcode - - An example with a long-untouched module that everyone has - >>> _linkcode_resolve('py', {'module': 'tty', - ... 'fullname': 'setraw'}, - ... package='tty', - ... url_fmt='http://hg.python.org/cpython/file/' - ... '{revision}/Lib/{package}/{path}#L{lineno}', - ... revision='xxxx') - 'http://hg.python.org/cpython/file/xxxx/Lib/tty/tty.py#L18' - """ - - if revision is None: - return - if domain not in ("py", "pyx"): - return - if not info.get("module") or not info.get("fullname"): - return - - class_name = info["fullname"].split(".")[0] - if type(class_name) != str: - # Python 2 only - class_name = class_name.encode("utf-8") - module = __import__(info["module"], fromlist=[class_name]) - obj = attrgetter(info["fullname"])(module) - - try: - fn = inspect.getsourcefile(obj) - except Exception: - fn = None - if not fn: - try: - fn = inspect.getsourcefile(sys.modules[obj.__module__]) - except Exception: - fn = None - if not fn: - return - - fn = os.path.relpath(fn, start=os.path.dirname(__import__(package).__file__)) - try: - lineno = inspect.getsourcelines(obj)[1] - except Exception: - lineno = "" - return url_fmt.format(revision=revision, package=package, path=fn, lineno=lineno) - - -def make_linkcode_resolve(package, url_fmt): - """Returns a linkcode_resolve function for the given URL format - - revision is a git commit reference (hash or name) - - package is the name of the root module of the package - - url_fmt is along the lines of ('https://github.com/USER/PROJECT/' - 'blob/{revision}/{package}/' - '{path}#L{lineno}') - """ - revision = _get_git_revision() - return partial(_linkcode_resolve, revision=revision, package=package, url_fmt=url_fmt) diff --git a/notebooks/pipeline.ipynb b/notebooks/pipeline.ipynb new file mode 100644 index 0000000..88a2fbe --- /dev/null +++ b/notebooks/pipeline.ipynb @@ -0,0 +1,632 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# maPCA walkthrough\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Import libraries" + ] + }, + { + "cell_type": "code", + "execution_count": 99, + "metadata": { + "execution": { + "iopub.execute_input": "2022-01-25T15:53:51.147898Z", + "iopub.status.busy": "2022-01-25T15:53:51.147398Z", + "iopub.status.idle": "2022-01-25T15:53:51.168751Z", + "shell.execute_reply": "2022-01-25T15:53:51.163109Z", + "shell.execute_reply.started": "2022-01-25T15:53:51.147833Z" + } + }, + "outputs": [], + "source": [ + "from mapca import mapca, utils\n", + "import numpy as np\n", + "from nilearn import datasets, input_data, masking, image, plotting\n", + "import matplotlib.pyplot as plt\n", + "from nilearn.plotting import plot_roi, plot_epi, show\n", + "import nibabel as nib\n", + "from nilearn._utils import check_niimg_3d, check_niimg_4d\n", + "from scipy.stats import kurtosis\n", + "from sklearn.decomposition import PCA\n", + "from sklearn.preprocessing import StandardScaler" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Load 4D data, reshape into space x time and normalize" + ] + }, + { + "cell_type": "code", + "execution_count": 104, + "metadata": { + "execution": { + "iopub.execute_input": "2022-01-25T15:55:33.607242Z", + "iopub.status.busy": "2022-01-25T15:55:33.606137Z", + "iopub.status.idle": "2022-01-25T15:55:35.554404Z", + "shell.execute_reply": "2022-01-25T15:55:35.553469Z", + "shell.execute_reply.started": "2022-01-25T15:55:33.607193Z" + } + }, + "outputs": [], + "source": [ + "# Download data\n", + "dataset = datasets.fetch_development_fmri(n_subjects=1)\n", + "\n", + "# Generate masker\n", + "brain_masker = input_data.NiftiMasker(\n", + " detrend=True, mask_strategy='epi',\n", + " low_pass=0.1, high_pass=0.01, t_r=2,\n", + " memory='nilearn_cache', memory_level=1, verbose=0)\n", + "\n", + "# Get 4D data\n", + "img = check_niimg_4d(dataset.func[0])\n", + "mask = check_niimg_3d(brain_masker.fit(dataset.func[0]).mask_img_)\n", + "data = img.get_fdata()\n", + "mask = mask.get_fdata()\n", + "\n", + "# Get shape and reshape\n", + "[n_x, n_y, n_z, n_timepoints] = data.shape\n", + "data_nib_V = np.reshape(data, (n_x * n_y * n_z, n_timepoints), order=\"F\")\n", + "mask_vec = np.reshape(mask, n_x * n_y * n_z, order=\"F\")\n", + "X = data_nib_V[mask_vec == 1, :]\n", + "\n", + "n_samples = np.sum(mask_vec)\n", + "\n", + "# Normalize time-series\n", + "scaler = StandardScaler(with_mean=True, with_std=True)\n", + "X = scaler.fit_transform(X.T).T" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "execution": { + "iopub.execute_input": "2022-01-25T15:05:42.185314Z", + "iopub.status.busy": "2022-01-25T15:05:42.184711Z", + "iopub.status.idle": "2022-01-25T15:05:42.236272Z", + "shell.execute_reply": "2022-01-25T15:05:42.213135Z", + "shell.execute_reply.started": "2022-01-25T15:05:42.185260Z" + } + }, + "source": [ + "## Step 1: Perform SVD on original data" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": { + "execution": { + "iopub.execute_input": "2022-01-25T15:08:14.604696Z", + "iopub.status.busy": "2022-01-25T15:08:14.604202Z", + "iopub.status.idle": "2022-01-25T15:08:16.070047Z", + "shell.execute_reply": "2022-01-25T15:08:16.048190Z", + "shell.execute_reply.started": "2022-01-25T15:08:14.604639Z" + } + }, + "outputs": [], + "source": [ + "# V = eigenvectors\n", + "V, eigenvalues = utils._icatb_svd(X, n_timepoints)\n", + "\n", + "# Reordering of values from highest to lowest\n", + "eigenvalues = eigenvalues[::-1]\n", + "dataN = np.dot(X, V[:, ::-1])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 2: Apply kurtosis threshold of Gaussianity to SVD data to get first estimate of subsampling depth" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": { + "execution": { + "iopub.execute_input": "2022-01-25T15:12:24.621362Z", + "iopub.status.busy": "2022-01-25T15:12:24.620713Z", + "iopub.status.idle": "2022-01-25T15:12:24.647967Z", + "shell.execute_reply": "2022-01-25T15:12:24.646372Z", + "shell.execute_reply.started": "2022-01-25T15:12:24.621289Z" + } + }, + "outputs": [], + "source": [ + "def apply_kurtosis(dataN, eigenvalues):\n", + " # Using 12 gaussian components from middle, top and bottom gaussian\n", + " # components to determine the subsampling depth.\n", + " # Final subsampling depth is determined using median\n", + " kurt = kurtosis(dataN, axis=0, fisher=True)\n", + " kurt[kurt < 0] = 0\n", + " kurt = np.expand_dims(kurt, 1)\n", + "\n", + " kurt[eigenvalues > np.mean(eigenvalues)] = 1000\n", + " idx_gauss = np.where(\n", + " ((kurt[:, 0] < 0.3) & (kurt[:, 0] > 0) & (eigenvalues > np.finfo(float).eps)) == 1\n", + " )[\n", + " 0\n", + " ]\n", + " idx = np.array(idx_gauss[:]).T\n", + " dfs = np.sum(eigenvalues > np.finfo(float).eps) # degrees of freedom\n", + " minTp = 12\n", + "\n", + " if len(idx) >= minTp:\n", + " middle = int(np.round(len(idx) / 2))\n", + " idx = np.hstack([idx[0:4], idx[middle - 1 : middle + 3], idx[-4:]])\n", + " else:\n", + " minTp = np.min([minTp, dfs])\n", + " idx = np.arange(dfs - minTp, dfs)\n", + "\n", + " idx = np.unique(idx)\n", + " \n", + " return idx" + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "metadata": { + "execution": { + "iopub.execute_input": "2022-01-25T15:12:37.071679Z", + "iopub.status.busy": "2022-01-25T15:12:37.071243Z", + "iopub.status.idle": "2022-01-25T15:12:37.227203Z", + "shell.execute_reply": "2022-01-25T15:12:37.224985Z", + "shell.execute_reply.started": "2022-01-25T15:12:37.071627Z" + } + }, + "outputs": [], + "source": [ + "idx = apply_kurtosis(dataN, eigenvalues)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 3: Calculate the subsampling depth by increasing its value until the entropy rate of the subsampled sequence reaches the upper bound" + ] + }, + { + "cell_type": "code", + "execution_count": 53, + "metadata": { + "execution": { + "iopub.execute_input": "2022-01-25T15:24:38.960093Z", + "iopub.status.busy": "2022-01-25T15:24:38.959800Z", + "iopub.status.idle": "2022-01-25T15:24:38.972372Z", + "shell.execute_reply": "2022-01-25T15:24:38.970012Z", + "shell.execute_reply.started": "2022-01-25T15:24:38.960065Z" + } + }, + "outputs": [], + "source": [ + "def estimate_subsampling_depth(dataN, mask_ND):\n", + " # Estimate the subsampling depth for effectively i.i.d. samples\n", + " sub_depth = len(idx)\n", + " sub_iid_sp = np.zeros((sub_depth,))\n", + " for i in range(sub_depth):\n", + " x_single = np.zeros(n_x * n_y * n_z)\n", + " x_single[mask_vec == 1] = dataN[:, idx[i]]\n", + " x_single = np.reshape(x_single, (n_x, n_y, n_z), order=\"F\")\n", + " sub_iid_sp[i] = utils._est_indp_sp(x_single)[0] + 1\n", + " if i > 6:\n", + " tmp_sub_sp = sub_iid_sp[0:i]\n", + " tmp_sub_median = np.round(np.median(tmp_sub_sp))\n", + " if np.sum(tmp_sub_sp == tmp_sub_median) > 6:\n", + " sub_iid_sp = tmp_sub_sp\n", + " break\n", + " dim_n = x_single.ndim\n", + "\n", + " sub_iid_sp_median = int(np.round(np.median(sub_iid_sp)))\n", + " if np.floor(np.power(n_samples / n_timepoints, 1 / dim_n)) < sub_iid_sp_median:\n", + " sub_iid_sp_median = int(np.floor(np.power(n_samples / n_timepoints, 1 / dim_n)))\n", + " \n", + " N = np.round(n_samples / np.power(sub_iid_sp_median, dim_n))\n", + " \n", + " return sub_iid_sp_median, N" + ] + }, + { + "cell_type": "code", + "execution_count": 54, + "metadata": { + "execution": { + "iopub.execute_input": "2022-01-25T15:24:39.519332Z", + "iopub.status.busy": "2022-01-25T15:24:39.518694Z", + "iopub.status.idle": "2022-01-25T15:24:42.161757Z", + "shell.execute_reply": "2022-01-25T15:24:42.160013Z", + "shell.execute_reply.started": "2022-01-25T15:24:39.519278Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Effective number of i.i.d. samples 3032\n" + ] + } + ], + "source": [ + "mask_ND = np.reshape(mask_vec, (n_x, n_y, n_z), order=\"F\")\n", + "sub_iid_sp_median, N = estimate_subsampling_depth(dataN, mask_ND )\n", + "\n", + "print(\"Effective number of i.i.d. samples %d\" % N)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 4: Generate subsampled data" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "metadata": { + "execution": { + "iopub.execute_input": "2022-01-25T15:25:30.255454Z", + "iopub.status.busy": "2022-01-25T15:25:30.254886Z", + "iopub.status.idle": "2022-01-25T15:25:30.273395Z", + "shell.execute_reply": "2022-01-25T15:25:30.269014Z", + "shell.execute_reply.started": "2022-01-25T15:25:30.255380Z" + } + }, + "outputs": [], + "source": [ + "def generate_subsampled_data(mask_ND, sub_iid_sp_median, n_timepoints, mask_shape):\n", + " n_x, n_y, n_z = mask_shape\n", + " mask_s = utils._subsampling(mask_ND, sub_iid_sp_median)\n", + " mask_s_1d = np.reshape(mask_s, np.prod(mask_s.shape), order=\"F\")\n", + " dat = np.zeros((int(np.sum(mask_s_1d)), n_timepoints))\n", + " for i_vol in range(n_timepoints):\n", + " x_single = np.zeros(n_x * n_y * n_z)\n", + " x_single[mask_vec == 1] = X[:, i_vol]\n", + " x_single = np.reshape(x_single, (n_x, n_y, n_z), order=\"F\")\n", + " dat0 = utils._subsampling(x_single, sub_iid_sp_median)\n", + " dat0 = np.reshape(dat0, np.prod(dat0.shape), order=\"F\")\n", + " dat[:, i_vol] = dat0[mask_s_1d == 1]\n", + " \n", + " return dat" + ] + }, + { + "cell_type": "code", + "execution_count": 58, + "metadata": { + "execution": { + "iopub.execute_input": "2022-01-25T15:25:30.971632Z", + "iopub.status.busy": "2022-01-25T15:25:30.971184Z", + "iopub.status.idle": "2022-01-25T15:25:31.294542Z", + "shell.execute_reply": "2022-01-25T15:25:31.293798Z", + "shell.execute_reply.started": "2022-01-25T15:25:30.971585Z" + } + }, + "outputs": [], + "source": [ + "dat = generate_subsampled_data(mask_ND, sub_iid_sp_median, n_timepoints, (n_x, n_y, n_z))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 5: Calculate SVD of covariance matrix of subsampled data" + ] + }, + { + "cell_type": "code", + "execution_count": 60, + "metadata": { + "execution": { + "iopub.execute_input": "2022-01-25T15:27:01.965930Z", + "iopub.status.busy": "2022-01-25T15:27:01.965556Z", + "iopub.status.idle": "2022-01-25T15:27:02.641985Z", + "shell.execute_reply": "2022-01-25T15:27:02.640663Z", + "shell.execute_reply.started": "2022-01-25T15:27:01.965887Z" + } + }, + "outputs": [], + "source": [ + "# Perform Variance Normalization\n", + "temp_scaler = StandardScaler(with_mean=True, with_std=True)\n", + "dat = temp_scaler.fit_transform(dat.T).T\n", + "\n", + "V, eigenvalues = utils._icatb_svd(dat, n_timepoints)\n", + "eigenvalues = eigenvalues[::-1]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 6: Make eigenspectrum adjustment" + ] + }, + { + "cell_type": "code", + "execution_count": 61, + "metadata": { + "execution": { + "iopub.execute_input": "2022-01-25T15:27:42.446795Z", + "iopub.status.busy": "2022-01-25T15:27:42.446281Z", + "iopub.status.idle": "2022-01-25T15:27:42.516496Z", + "shell.execute_reply": "2022-01-25T15:27:42.514602Z", + "shell.execute_reply.started": "2022-01-25T15:27:42.446741Z" + } + }, + "outputs": [], + "source": [ + "# Make eigen spectrum adjustment\n", + "eigenvalues = utils._eigensp_adj(eigenvalues, N, eigenvalues.shape[0])\n", + "# (completed)\n", + "if np.sum(np.imag(eigenvalues)):\n", + " raise ValueError(\"Invalid eigenvalue found for the subsampled data.\")\n", + "\n", + "# Correction on the ill-conditioned results (when tdim is large,\n", + "# some least significant eigenvalues become small negative numbers)\n", + "if eigenvalues[np.real(eigenvalues) <= np.finfo(float).eps].shape[0] > 0:\n", + " eigenvalues[np.real(eigenvalues) <= np.finfo(float).eps] = np.min(\n", + " eigenvalues[np.real(eigenvalues) >= np.finfo(float).eps]\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 7: Estimate dimensionality with AIC, KIC and MDL" + ] + }, + { + "cell_type": "code", + "execution_count": 81, + "metadata": { + "execution": { + "iopub.execute_input": "2022-01-25T15:41:21.924817Z", + "iopub.status.busy": "2022-01-25T15:41:21.924354Z", + "iopub.status.idle": "2022-01-25T15:41:21.942970Z", + "shell.execute_reply": "2022-01-25T15:41:21.940437Z", + "shell.execute_reply.started": "2022-01-25T15:41:21.924765Z" + } + }, + "outputs": [], + "source": [ + "def estimate_dimensionality(n_timepoints, eigenvalues, N):\n", + " p = n_timepoints\n", + " aic = np.zeros(p - 1)\n", + " kic = np.zeros(p - 1)\n", + " mdl = np.zeros(p - 1)\n", + "\n", + " for k_idx, k in enumerate(np.arange(1, p)):\n", + " LH = np.log(np.prod(np.power(eigenvalues[k:], 1 / (p - k))) / np.mean(eigenvalues[k:]))\n", + " mlh = 0.5 * N * (p - k) * LH\n", + " df = 1 + 0.5 * k * (2 * p - k + 1)\n", + " aic[k_idx] = (-2 * mlh) + (2 * df)\n", + " kic[k_idx] = (-2 * mlh) + (3 * df)\n", + " mdl[k_idx] = -mlh + (0.5 * df * np.log(N))\n", + "\n", + " itc = np.row_stack([aic, kic, mdl])\n", + "\n", + " # AIC\n", + " criteria_idx = 0\n", + " dlap = np.diff(itc[criteria_idx, :])\n", + " a = np.where(dlap > 0)[0] + 1 # Plus 1 to\n", + " if a.size == 0:\n", + " n_aic = itc[criteria_idx, :].shape[0]\n", + " else:\n", + " n_aic = a[0]\n", + "\n", + " # KIC\n", + " criteria_idx = 1\n", + " dlap = np.diff(itc[criteria_idx, :])\n", + " a = np.where(dlap > 0)[0] + 1 # Plus 1 to\n", + " if a.size == 0:\n", + " n_kic = itc[criteria_idx, :].shape[0]\n", + " else:\n", + " n_kic = a[0]\n", + " \n", + " # MDL\n", + " criteria_idx = 2\n", + " dlap = np.diff(itc[criteria_idx, :])\n", + " a = np.where(dlap > 0)[0] + 1 # Plus 1 to\n", + " if a.size == 0:\n", + " n_mdl = itc[criteria_idx, :].shape[0]\n", + " else:\n", + " n_mdl = a[0]\n", + " \n", + "\n", + " return np.array([aic, kic, mdl]), np.array([n_aic, n_kic, n_mdl])" + ] + }, + { + "cell_type": "code", + "execution_count": 91, + "metadata": { + "execution": { + "iopub.execute_input": "2022-01-25T15:47:07.481070Z", + "iopub.status.busy": "2022-01-25T15:47:07.472318Z", + "iopub.status.idle": "2022-01-25T15:47:08.096438Z", + "shell.execute_reply": "2022-01-25T15:47:08.094187Z", + "shell.execute_reply.started": "2022-01-25T15:47:07.480073Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0, 0.5, 'A.U.')" + ] + }, + "execution_count": 91, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABJ4AAAJNCAYAAABwab9RAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAADMHUlEQVR4nOz9d5gc1Z33f3/OdPfknGc0ozzKEQmByHERGYMx2MYGjAkOi73evQ342d+a9Qab3b3W4b531ws2AmyDyDnIJBGMjJCIyjmMwuQcO5znj2qNAgKrYapruvv9uq66qrqqTp1vg81Inzl1jrHWCgAAAAAAABhuaV4XAAAAAAAAgORE8AQAAAAAAABXEDwBAAAAAADAFQRPAAAAAAAAcAXBEwAAAAAAAFxB8AQAAAAAAABX+L0uIN5KS0vt2LFjvS4DAAAAAAAgaaxatarZWlt2+PmUC57Gjh2rlStXel0GAAAAAABA0jDG7DjSeV61AwAAAAAAgCsIngAAAAAAAOAKgicAAAAAAAC4IuXmeAIAAAAAAHBDMBhUfX29+vv7vS7FNZmZmaqpqVEgEDiq+wmeAAAAAAAAhkF9fb3y8vI0duxYGWO8LmfYWWvV0tKi+vp6jRs37qja8KodAAAAAADAMOjv71dJSUlShk6SZIxRSUlJTCO6CJ4AAAAAAACGSbKGTvvF+v0IngAAAAAAAJLIE088IWOM1q9fL0navn27ZsyYMXR9xYoVOuWUUzR58mTNnTtX3/zmN9Xb2+tKLQRPAAAAAAAASeSBBx7QSSedpAceeOBj1xoaGnT55Zfrjjvu0IYNG/Tee+9p0aJF6urqcqUWgicAAAAAAIAk0d3drTfffFO//e1vtWTJko9d/6//+i9dffXVWrhw4dC5L37xi6qoqHClHoInAAAAAACAJPHkk09q0aJFmjRpkkpKSrRq1apDrq9evVrz5s2LWz3+uPUEAAAAAACQIv7x6TVau6dzWJ85rTpfP75w+qfe88ADD+h73/ueJOnKK6/UAw88oO9+97vDWkcsCJ4AAAAAAACSQGtrq1555RV99NFHMsYoHA7LGKPvfOc7Q/dMnz5dq1at0sUXXxyXmgieAAAAAAAAhtlfGpnkhkceeURf+9rX9L//+79D50499VTt2rVr6PN3v/tdLViwQOeff76OO+44SdJjjz2mE0880ZV5npjjCQAAAAAAIAk88MAD+sIXvnDIucsuu0w//elPhz5XVFRoyZIl+ru/+ztNnjxZU6dO1dKlS5WXl+dKTcZa68qDR6r58+fblStXel0GAAAAAABIMuvWrdPUqVO9LsN1R/qexphV1tr5h9/LiCcAAAAAAAC4guAJAAAAAAAAriB4AgAAAAAAgCsIngAAAAAAAOAKgicAAAAAAAC4guAJAAAAAAAAriB4AgAAAAAASBK5ublDx88995wmTZqkHTt26Pbbb9d//Md/DF37j//4D02ZMkVz5szRscceq/vuu8+VevyuPBUYwa594VpJ0uJFi+Pb8eLzowU8e9RNrvjf5ZKkB29c6EZFAAAAAIAk9fLLL+vmm2/W0qVLNWbMmEOu/frXv9aLL76oFStWKD8/X52dnXr88cddqYPgCSnnP0/7T286/lLs6fH/XDXPhUIAAAAAAMns9ddf1/XXX6/nnntOEyZM+Nj1f/3Xf9WyZcuUn58vScrPz9fVV1/tSi0ET0g5RZlF3nScUxJzk+KcdBcKAQAAAAAkq4GBAV1yySVatmyZpkyZ8rHrnZ2d6urq0vjx4+NSD8ETUs4Tm5+QJF0y8ZL4dvzeH5z93K8edZOHV+6SJF0+v9aNigAAAAAAbnn+VmnfR8P7zMqZ0rk/+9RbAoGATjjhBP32t7/VL3/5y+Ht/zNgcnGknCc3P6knNz8Z/47fv9/ZYvDIqno9sqrepYIAAAAAAMkmLS1NDz30kFasWKF//dd//dj1/Px85ebmauvWrXGphxFPAAAAAAAAw+0vjExyU3Z2tp599lmdfPLJqqio0HXXXXfI9dtuu03f+c539OCDDyo/P1/d3d167LHH9PWvf33YayF4AgAAAAAASDLFxcV64YUXdMopp6isrOyQa9/61rfU3d2tY489VoFAQIFAQH/7t3/rSh0ETwAAAAAAAEmiu7t76Li2tlbbtm2TJF100UVD540x+uEPf6gf/vCHrtfDHE8AAAAAAABwBSOekHL++6z/9qbjrz4cc5N7rl3gQiEAAAAAAMQHwRNSTpY/y5uO07NjbpKV7nOhEAAAAAAA4oNX7ZBylqxfoiXrl8S/4xV3OVsMfrd8u363fLs79QAAAAAA4DKCJ6ScpduXaun2pfHveM0TzhaDZz7cq2c+3OtKOQAAAAAAuI3gCQAAAAAAAK4geAIAAAAAAEgSxhhdddVVQ59DoZDKysp0wQUXSJLuuecelZWVae7cuaqrq9M555yjt956a+j+a665Ro888siw1UPwBAAAAAAAkCRycnK0evVq9fX1SZJefPFFjRo16pB7rrjiCr333nvatGmTbr31Vl166aVat26dK/UQPAEAAAAAACSR8847T88++6wk6YEHHtCXv/zlT7z39NNP1w033KA777zTlVoInpByFi9arMWLFse/42ufdbYYPHjjQj1440KXCgIAAAAAJKMrr7xSS5YsUX9/vz788EMdd9xxn3r/Mccco/Xr17tSi9+VpwIAAAAAAKSwO1bcofWtwxvmTCmeolsW3PIX75s1a5a2b9+uBx54QOedd95fvN9aOxzlHRHBE1LOPavvkSRdM+Oa+Hb8p185+xNvPuomd76+RZJ0wykT3KgIAAAAAJCkLrroIv3d3/2dli1bppaWlk+997333tPUqVNdqYPgCSnntfrXJHkQPG1c6uxjCJ5eXtcoieAJAAAAABLN0YxMctM3vvENFRYWaubMmVq2bNkn3vfaa6/pzjvv1KuvvupKHQRPAAAAAAAASaampkY333zkgQ8PPvig3nzzTfX29mrcuHF69NFHDxnxdOONN+r73/++JKm2tlbLly//zHUQPAEAAAAAACSJ7u7uj5077bTTdNppp0mSrrnmGl1zzTWf2P6ee+4Z1npY1Q4AAAAAAACuYMQTUk6GP8ObjgOZMTfJDPhcKAQAAAAAgPggeELK+fVZv/am46sejbnJvd9Y4EIhAAAAAADEB6/aAQAAAAAADBNrrdcluCrW70fwhJTz6w9+rV9/4MGop9f+zdli8KuXN+lXL29yqSAAAAAAwHDKzMxUS0tL0oZP1lq1tLQoM/Pop5LhVTuknLf3vi1Jumn2TfHteOtrzv7UHx51kz9tbpYk3XxmnRsVAQAAAACGUU1Njerr69XU1OR1Ka7JzMxUTU3NUd9P8AQAAAAAADAMAoGAxo0b53UZI4qrr9oZY/7GGLPGGLPaGPOAMSbTGDPOGPO2MWazMeZBY0x69N6M6OfN0etjD3rObdHzG4wx5xx0flH03GZjzK1ufhcAAAAAAADExrXgyRgzStLNkuZba2dI8km6UtIdkn5urZ0oqU3SddEm10lqi57/efQ+GWOmRdtNl7RI0n8bY3zGGJ+k/5J0rqRpkr4cvRcAAAAAAAAjgNuTi/slZRlj/JKyJe2VdIakR6LX75V0SfT44uhnRa+faYwx0fNLrLUD1tptkjZLWhDdNltrt1prByUtid4LfKrCjEIVZhTGv+PsImeLQVF2uoqy010qCAAAAAAAd7k2x5O1drcx5j8k7ZTUJ+mPklZJarfWhqK31UsaFT0eJWlXtG3IGNMhqSR6/s8HPfrgNrsOO3+cC18FSebnp//cm46v+H3MTX79tXkuFAIAAAAAQHy4+apdkZwRSOMkVUvKkfOqXNwZY24wxqw0xqxM5pnlAQAAAAAARhI3X7U7S9I2a22TtTYo6TFJJ0oqjL56J0k1knZHj3dLqpWk6PUCSS0Hnz+szSed/xhr7Z3W2vnW2vllZWXD8d2QwH6x6hf6xapfxL/jl253thjc8cJ63fHCelfKAQAAAADAba69aifnFbvjjTHZcl61O1PSSkmvSvqinDmZrpb0ZPT+p6Kfl0evv2KttcaYpyTdb4z5Tzkjp+okrZBkJNUZY8bJCZyulPQVF78PksQHTR940/Gud2Ju8u6ONhcKAQAAAAAgPtyc4+ltY8wjkt6VFJL0nqQ7JT0raYkx5p+j534bbfJbSb8zxmyW1ConSJK1do0x5iFJa6PP+Y61NixJxpjvSloqZ8W8u621a9z6PgAAAAAAAIiNmyOeZK39saQfH3Z6q5wV6Q6/t1/S5Z/wnH+R9C9HOP+cpOc+f6UAAAAAAAAYbm7O8QQAAAAAAIAU5uqIJ2Akqsip8Kbj/OqYm1QVZLpQCAAAAAAA8WGstV7XEFfz58+3K1eu9LoMAAAAAACApGGMWWWtnX/4eV61AwAAAAAAgCsInpBy7lhxh+5YcUf8O37+VmeLwT8+vUb/+DSLNQIAAAAAEhNzPCHlrG9d703H+z6KucnaPZ0uFAIAAAAAQHww4gkAAAAAAACuIHgCAAAAAACAKwieAAAAAAAA4ArmeEpEA11S5x6pdJJkjNfVJJwx+WO86bhkQsxNxpfluFAIAAAAAADxYay1XtcQV/Pnz7crV670uozPJfLGz5X28u3SbbuljFyvywEAAAAAACnOGLPKWjv/8PO8apeAltVHJEn9HY0eVwIAAAAAAPDJCJ4SUCCvTJLU2brX40oS0+1v3a7b37o9/h0/dbOzxeC2xz7UbY996FJBAAAAAAC4izmeElBGQYUkqbe1weNKEtOOzh3edNyyJeYmW5t6XCgEAAAAAID4YMRTAsoqLJck9XcQPAEAAAAAgJGL4CkB5ZVUSZKC3c0eVwIAAAAAAPDJCJ4SUFFhsQasXxGCJwAAAAAAMIIxx1MCys8KqEH5Mr0ET5/FlOIp3nRcOTPmJtOq810oBAAAAACA+CB4SkDGGHWYAgX6W70uJSHdsuAWbzo+92cxN/nxhdNdKAQAAAAAgPjgVbsE1e0vUMZgm9dlAAAAAAAAfCKCpwTVFyhSdojg6bO49Y1bdesbt8a/40evd7YYfH/Je/r+kvdcKggAAAAAAHfxql2CGswoVl5fh9dlJKSGngZvOu7cE3OTvR39LhQCAAAAAEB8MOIpQUWySpSjPik04HUpAAAAAAAAR0TwlKiySyRJoa5GjwsBAAAAAAA4MoKnBOXLK5MkdbV69NoYAAAAAADAX8AcTwkqPb9cktTTtk9FHteSaGaXzfam49pjY25yzBj+7QIAAAAAEhfBU4LKKqiQJPW1MeIpVt+f931vOj7r9pib3LJoyvDXAQAAAABAnPCqXYLKLXGCp2BXk8eVAAAAAAAAHBnBU4IqLC5TyKYp3E3wFKu/efVv9Dev/k38O37wKmeLwU2/W6WbfrfKpYIAAAAAAHAXr9olqKKcTLUpT+pp8bqUhNM+0O5Nx71tMTdp6x10oRAAAAAAAOKDEU8JKuBLU7vJl7+f4AkAAAAAAIxMBE8JrMtXqPTBVq/LAAAAAAAAOCKCpwTW5y9SVrDd6zIAAAAAAACOiDmeEthAepHyutu9LiPhHFd1nDcdjz815iYnTix1oRAAAAAAAOKD4CmBhbNKlNfVLYWDki/gdTkJ46bZN3nT8ak/jLnJzWfWuVAIAAAAAADxwat2Ccxmlzj7XiYYBwAAAAAAIw/BUwJLyy2TJPW2NXhcSWK56aWbdNNLHox6+v1lzhaDq+9eoavvXuFSQQAAAAAAuIvgKYEF8pzgqbt1n8eVJJaB0IAGQgPx7zjY72wx6A+G1R8Mu1QQAAAAAADuInhKYJlFFZKk3nZGPAEAAAAAgJGH4CmB5RY6wdNgZ6PHlQAAAAAAAHwcwVMCyy+uUMQahbqavC4FAAAAAADgY/xeF4DPriQ/S+3KkXpY1S4Wp9ac6k3Hk86JucmZU8tdKAQAAAAAgPggeEpgWQGf9ilfaX0ET7G4ZsY13nR84s0xN7nhlAkuFAIAAAAAQHzwql0CM8aoM61A6QOtXpcCAAAAAADwMQRPCa43UKjMYJvXZSSUa1+4Vte+cG38O158vrPF4Ir/Xa4r/ne5SwUBAAAAAOAugqcE1x8oUk6o3esyAAAAAAAAPobgKcGFMouVZ7ukSMTrUgAAAAAAAA5B8JTgIlml8iki9bd7XQoAAAAAAMAhCJ4SXFpuqSRpsLPB40oAAAAAAAAO5fe6AHw+/rwySVJXS4NKKqd6XE1iOGfsOd50PP2SmJtcMKtq+OsAAAAAACBOCJ4SXEZBhSSpp32fSjyuJVFcOeVKbzpecH3MTb62cOzw1wEAAAAAQJzwql2Cyy50gqeBjkaPK0kcfaE+9YX64t/xYK+zxaBvMKy+wbBLBQEAAAAA4C5GPCW4vBIneAp1NXlcSeL49kvfliQtXrQ4vh3/4XJnf+2zR93kmsUrJEkP3rjQjYoAAAAAAHAVI54SXEl+njptlmw3wRMAAAAAABhZCJ4SXEFWQK3Kl+lr8boUAAAAAACAQ7gWPBljJhtj3j9o6zTGfN8YU2yMedEYsym6L4reb4wxvzLGbDbGfGiMOeagZ10dvX+TMebqg87PM8Z8FG3zK2OMcev7jFRpaUadpkCB/lavSwEAAAAAADiEa8GTtXaDtXaOtXaOpHmSeiU9LulWSS9ba+skvRz9LEnnSqqLbjdI+h9JMsYUS/qxpOMkLZD04/1hVfSe6w9qt8it7zOS9fgLlTHY5nUZAAAAAAAAh4jX5OJnStpird1hjLlY0mnR8/dKWibpFkkXS7rPWmsl/dkYU2iMqYre+6K1tlWSjDEvSlpkjFkmKd9a++fo+fskXSLp+Th9pxGjL1ConIEtXpeRMC6eeLE3Hc/5SsxNvjivxoVCAAAAAACIj3gFT1dKeiB6XGGt3Rs93iepIno8StKug9rUR8992vn6I5xPOcGMYuX1dUjWSqn3tmHMLpl4iTcdz/1qzE0un1/rQiEAAAAAAMSH65OLG2PSJV0k6eHDr0VHN9k41HCDMWalMWZlU1Pyrf4WzipRQCFpoNPrUhJCW3+b2vo9eDWxp8XZYtDaM6jWnkGXCgIAAAAAwF3xWNXuXEnvWmsbop8boq/QKbpvjJ7fLeng4R010XOfdr7mCOc/xlp7p7V2vrV2fllZ2ef8OiOPySmVJEW6mz2uJDH8YNkP9INlP4h/xw993dli8K3fr9K3fr/KpYIAAAAAAHBXPIKnL+vAa3aS9JSk/SvTXS3pyYPOfz26ut3xkjqir+QtlfRXxpii6KTifyVpafRapzHm+Ohqdl8/6FkpxZfnhGndrfs8rgQAAAAAAOAAV+d4MsbkSDpb0o0Hnf6ZpIeMMddJ2iHpS9Hzz0k6T9JmOSvgXStJ1tpWY8w/SXonet9P9k80Lunbku6RlCVnUvGUm1hckjLyyyVJ3W0Nyve4FgAAAAAAgP1cDZ6stT2SSg471yJnlbvD77WSvvMJz7lb0t1HOL9S0oxhKTaBZRU687P3dzT8hTsBAAAAAADiJx6v2sFluSWVkqRgZ/JNnA4AAAAAABKXqyOeEB/FBYXqtRmKdBM8HY0rJl/hTcfHfiPmJlcdP8aFQgAAAAAAiA+CpyRQlBNQk/JkelnV7mgsGrfIm45nXBZzkwtnV7tQCAAAAAAA8cGrdkkgw+9Tuwrk72/9yzdD+3r2aV+PBysAdtQ7Wwz2tPdpT3ufSwUBAAAAAOAuRjwliW5fgfIHCZ6Oxm1v3CZJWrxocXw7fiy6uOO1zx51k7958H1J0oM3LnShIAAAAAAA3MWIpyTRGyhSdrDd6zIAAAAAAACGEDwlicGMYuWGO7wuAwAAAAAAYAjBU5IIZxYrUwPSYI/XpQAAAAAAAEgieEoaNqfU2fc0eVwJAAAAAACAg8nFk4Qv1wmeBjoalVk01ttiRrirp1/tTccnfDfmJtefPN6FQgAAAAAAiA+CpyQRyC+XJHW1NihzrLe1jHSn1Z7mTceTz425yVnTKlwoBAAAAACA+OBVuySRVeAEFH3tDR5XMvJt69imbR3b4t9x8yZni8GWpm5taep2qSAAAAAAANzFiKckkVPsBE+DnY0eVzLy/WT5TyRJixctjm/HT3/f2V/77FE3+dFjH0mSHrxxoQsFAQAAAADgLkY8JYmiwmINWL/CXUwuDgAAAAAARgaCpyRRlJuhVuVLvS1elwIAAAAAACCJ4Clp5GX41aY8pfW3el0KAAAAAACAJIKnpGGMUVdaodIJngAAAAAAwAjB5OJJpMdfqKxgbKumpaIbZt3gTcen/F3MTf76jDoXCgEAAAAAID4InpLIQEaRcnravS5jxFtY7dEKcRNOj7nJSXWlLhQCAAAAAEB88KpdEglllijH9kqhAa9LGdHWt67X+tb18e9474fOFoM1ezq0Zk+HSwUBAAAAAOAuRjwlEZtV4hz0tkj51d4WM4LdseIOSdLiRYvj2/ELtzn7a5896iY/eXqtJOnBGz0apQUAAAAAwOfAiKckYnLLJUnBzgaPKwEAAAAAACB4Sir+ggpJUk/rHo8rAQAAAAAAIHhKKpmFVZKk3haCJwAAAAAA4D2CpySSVzpKkjTQvtfjSgAAAAAAAJhcPKmUFBao02YrzBxPn+p7x3zPm47P/IeYm/xw0WQXCgEAAAAAID4InpJIWV6GGmyBTA/B06eZUz7Hm45HHxdzk3ljil0oBAAAAACA+OBVuySSk+5TiylSoLfJ61JGtPcb39f7je/Hv+OdbztbDFbtaNWqHa0uFQQAAAAAgLsInpKIMUZd/mJlDjR7XcqI9st3f6lfvvvL+Hf88k+cLQb/9sIG/dsLG1wqCAAAAAAAdxE8JZm+9BLlhRghAwAAAAAAvEfwlGQGs8qUZXulwR6vSwEAAAAAACmO4CnJRHIqnIPuRm8LAQAAAAAAKY/gKcn48p3gKdS5z+NKAAAAAABAqvN7XQCGV3phlSSpu2WPCsd6W8tIdcuCW7zpeNFPY27yDxdOc6EQAAAAAADig+ApyeQUV0uSelv3qNDbUkasKcVTvOm4albMTaZXF7hQCAAAAAAA8cGrdkkmv6RSYWs02L7X61JGrOV7lmv5nuXx73jLq84Wgzc3NevNTc0uFQQAAAAAgLsY8ZRkyvKz1aIChTsbvC5lxLrzwzslSQurF8a349f/w9lPOP2om/zfVzZJkk6qK3WjIgAAAAAAXMWIpyRTlpehJlugtB6CJwAAAAAA4C2CpySTGfCpzRQpvb/J61IAAAAAAECKI3hKQt2BYmUPtnhdBgAAAAAASHEET0moL6NMeaE2KRLxuhQAAAAAAJDCmFw8CYWyy+TvCUn97VJ2sdfljDj/sPAfvOn4wl/E3ORfL505/HUAAAAAABAnBE/JKLdCapLU3UDwdATjCsZ503FpXcxNJpTlulAIAAAAAADxwat2SciXXylJGmzf63ElI9OyXcu0bNey+He84Xlni8FLaxv00lpWKAQAAAAAJCZGPCWhzMIqSVJXy26VeFzLSHTvmnslSafVnhbfjt/6f85+8rlH3eSuN7ZKks6aVuFGRQAAAAAAuIoRT0kop7RaktTfyognAAAAAADgHYKnJFRcWKI+m65gB8ETAAAAAADwDsFTEirLz1STLZDtZm4gAAAAAADgHYKnJFSSm64mFcrX0+R1KQAAAAAAIIUxuXgSCvjS1J5WpJoBgqcj+enJP/Wm40v/N+YmP79izvDXAQAAAABAnBA8Jane9BLlDK73uowRqTKn0puOC2piblJdmOVCIQAAAAAAxAev2iWp/swy5UY6pdCg16WMOC9se0EvbHsh/h2vftTZYvD0B3v09Ad7XCoIAAAAAAB3MeIpSYWzy6VOST1NUsEor8sZUR7c8KAkadG4RfHt+J27nf2My466ye//vEOSdOHsajcqAgAAAADAVYx4SlImt9w5YGU7AAAAAADgEYKnJBUoqJIk9bft9bgSAAAAAACQqlwNnowxhcaYR4wx640x64wxC40xxcaYF40xm6L7oui9xhjzK2PMZmPMh8aYYw56ztXR+zcZY64+6Pw8Y8xH0Ta/MsYYN79PIskscoKn7pbdHlcCAAAAAABSldsjnn4p6QVr7RRJsyWtk3SrpJettXWSXo5+lqRzJdVFtxsk/Y8kGWOKJf1Y0nGSFkj68f6wKnrP9Qe1i/OkPSNXXkl0xFM7E1MDAAAAAABvuDa5uDGmQNIpkq6RJGvtoKRBY8zFkk6L3navpGWSbpF0saT7rLVW0p+jo6Wqove+aK1tjT73RUmLjDHLJOVba/8cPX+fpEskPe/Wd0okpYX5arO5Cncwx9Ph/vO0//Sm4y/dF3OT/7lqnguFAAAAAAAQH26uajdOUpOkxcaY2ZJWSfqepApr7f6Jh/ZJqogej5K066D29dFzn3a+/gjnIaksL0NNtkAZTC7+MUWZRX/5JjfklMTcpDgn3YVCAAAAAACIDzdftfNLOkbS/1hr50rq0YHX6iRJ0dFN1sUaJEnGmBuMMSuNMSubmprc7m5EKMpOV7MK5e9Lje8biyc2P6EnNj8R/47f+4OzxeDhlbv08Mpdf/lGAAAAAABGIDeDp3pJ9dbat6OfH5ETRDVEX6FTdN8Yvb5bUu1B7Wui5z7tfM0Rzn+MtfZOa+18a+38srKyz/WlEoUvzajTV6zMgWavSxlxntz8pJ7c/GT8O37/fmeLwSOr6vXIqvq/fCMAAAAAACOQa8GTtXafpF3GmMnRU2dKWivpKUn7V6a7WtL+BOApSV+Prm53vKSO6Ct5SyX9lTGmKDqp+F9JWhq91mmMOT66mt3XD3oWJPWmlyg32CpZ1weVAQAAAAAAfIybczxJ0l9L+oMxJl3SVknXygm7HjLGXCdph6QvRe99TtJ5kjZL6o3eK2ttqzHmnyS9E73vJ/snGpf0bUn3SMqSM6k4E4sfZCCrXBkD/dJgt5SR53U5AAAAAAAgxbgaPFlr35c0/wiXzjzCvVbSdz7hOXdLuvsI51dKmvH5qkxekexyqV1SdyPBEwAAAAAAiDs353iCx9LynQUDbdc+jysBAAAAAACpyO1X7eChjMIqSVJv617ljPW2lpHkv8/6b286/urDMTe559oFLhQCAAAAAEB8EDwlsayi/cFTvXI8rmUkyfJnedNxenbMTbLSfS4UAgAAAABAfPCqXRIrKClX0Po00M6rdgdbsn6JlqxfEv+OV9zlbDH43fLt+t3y7e7UAwAAAACAywieklh5fpaaVaBIZ4PXpYwoS7cv1dLtS+Pf8ZonnC0Gz3y4V898uNeVcgAAAAAAcBvBUxIrzc1Qky2Q6SF4AgAAAAAA8UfwlMQKsgJqUaECfc1elwIAAAAAAFIQwVMSM8aoy1+srEGCJwAAAAAAEH8ET0muL6NMuaE2KRL2uhQAAAAAAJBi/F4XAHcFs8rk64tIva1SbpnX5YwIixct9qbja5+NucmDNy50oRAAAAAAAOKDEU/JLrfc2XczwTgAAAAAAIgvgqck58uvlCSFuwie9rtn9T26Z/U98e/4T79ythjc+foW3fn6FpcKAgAAAADAXQRPSS6jsEqS1NtS73ElI8dr9a/ptfrX4t/xxqXOFoOX1zXq5XWNLhUEAAAAAIC7CJ6SXE6JEzz1te31uBIAAAAAAJBqmFw8yRUXFavbZmqwneAJAAAAAJC8wpGwBsIDGgwPHrI/eBsMD6o/3H/ItcHwoEKRkEKRkIKR4CH7kA0NXTv8nqHjg+6J2IjCNqyIjQxtn/b5komX6NYFt3r9j85VBE9JrjQ3Q022QH4mFwcAAAAAeCwYCao32Ku+UN+Bfaj3Mx0fHiqFIqHPXV+aSVMgLSB/mt/ZjH/oeP/5Q66n+ZVhMoaOfcanNJPm7NPSlKY053OaT0Zm6Pr+bU75nM//D3WEI3hKcmV5GVqrQo3qafK6lBEjw5/hTceBzJibZAZ8LhQCAAAAALEbDA+qa7BL3cFudQ92qyvY5eyPcK47GD1/8HGwWwPhgaPuz2/8ygpkKdufrSx/lrIDzr4ks0Q1uTXK8mcpw5ehdF+6MnwZyvBnOPvouUxf5tC1/Z8Pvv/ge/YHR2mGGYmGG8FTkstJ96nVFGp8/z6vSxkxfn3Wr73p+KpHY25y7zcWuFAIAAAAgFRlrVVfqE8dAx3qGOxQx0CH2gfa1THQoc7BTuf8J5wbjAz+xedn+7OVm56rvECectNzVZBZoJq8mqFzOYEc5QRyhkKkbH/2Icf7A6Zsf7YCvkAc/onAbQRPSc4Yo55AibIH13hdCgAAAABgmA2EB9TW36bW/tah/f6trb9NbQNt6hzoPCRoCkaCn/i8TF+m8jPyVZhRqIKMAo3NH6uCjALlp+crL90Jk3IDuc7x/n30XG4gV7403trAoQieUkB/Rpmye7qlYP9net0r2fz6A2fE002zb4pvx6/9m7M/9YdH3eRXL2+SJN18Zp0bFQEAAAAYYUKR0IHwqK9VrQPOvm3g46FSa3+reoI9R3yOP82v4oxiFWYWqjCjUOMLxys/PV8FGQVDoVJBeoHyMw6cy0/PV6afvzNieBE8pYBQdpnUI6mnUSoc7XU5nnt779uSPAietr7m7GMInv60uVkSwRMAAACQ6HqDvWrua1ZTX5Oa+prU0teipt6Djvua1NzXrLb+NlnZj7X3G7+KMotUnFmsoswijSodpZLMkkPO7f9clFmkvECejDEefFPgUARPKcDkVUhNkroaCJ4AAAAAYBj1Bnu1r3efGnoa1NjbeEiQ1NTbpJZ+J2DqDfV+rK3f+FWSVaKyrDJV51ZrVtkslWWVqTSr9GOhUn56PkESEhLBUwrw5VdJkkIde+Sv9bgYAAAAAEgA1lp1DnZqX88+NfY2qqG3wdl6DuwbexvVFez6WNucQI7KsspUklWiqcVTdfKok1WaVaqybCdUKs0qVVlWmQoyClhFDUmP4CkFZJQ4aVNv0w7le1wLAAAAAIwEnYOd2tu9V3u692hfbzRc2h8qRYOl/nD/IW2MjEqzSlWRXaGxBWN1XNVxqsipUEW2s5Vnl6s0q1TZgWyPvhUw8hA8pYD8ogr124AGWnZ6XcqIUJhR6E3H2UUxNynKTnehEAAAACC5WWvV0t/iBEs9e7S3e692d+/W3p4Dn7uD3Ye08Ru/yrLLVJFdoSnFU3RqzalOoBQNlipzKlWSVaJAWsCjbwUkJoKnFFCWn6ndtlS5HfVelzIi/Pz0n3vT8RW/j7nJr782z4VCAAAAgMQWjoTV2Nt4IEzq3nPIfm/PXg2EBw5pkxvIVVVulapzqjW/Yr6qc6qHPu8PlXjtDRh+BE8poKogS5ttsSZ37va6FAAAAAA4Kj3BHtV31Ttbd712de1SfbfzeXf3boUioUPuL84sVnVOteqK6nRqzalDoVJ1rhMw5acz8QjgBYKnFFCWl6E3VarZfWu9LmVE+MWqX0iSvj/v+/Ht+KXbnf1Ztx91kzteWC9JumXRlOGvBwAAAPBQxEbU2Nuo+q5DQ6X9+9b+1kPuz0vPU21erSYXTdaZo89UTV6NRuWMUlVulapyqpTpz/TomwD4NARPKcCXZtSZXq6cwTekcFDypfY7yR80feBNx7veibnJuzvaXCgEAAAAiI9wJKy9PXu1s3OndnTt0I7OHdrZuVP13fXa3bVbg5HBoXvTTJqqcqpUk1ej02tPV01ejWrzalWTV6Oa3BoVZBR4+E0AfFYETyliILtKaV0RqWuvVDja63IAAAAAJAlrrZr6mrSjc8fHtl1duxSMBIfuzfJnaUz+GE0snKjTak5zQqW8GtXm1qoyt5KJu4EkRPCUIiL5o6QuSR27CZ4AAAAAxMRaq/aB9o8FSzu7dmpH5w71hfqG7k1PS9fo/NEamz9Wp9aeqjF5YzQm39lKs0pljPHwmwCIN4KnFOEvqpV2S7ajXvxnHgAAAMCRhCNh7enZo20d27S1fau2dR7Ydwx0DN3nMz7V5NVodN5oza+YPxQsjckfo8qcSlaHAzCE4ClFZJWOkST1Ne9Utse1eK0ip8KbjvOrY25SVcAEiQAAABh+A+EBbe/Yrm2d27StfZu2dmzV1o6t2tG5QwPhgaH7ijOLNb5gvM4Zc47GFowdCpeqc6t5LQ7AUSF4ShHlpSXqtNkKNu9I+eDpZyf/zJuOL7sr5ia/uHKuC4UAAAAgVXQMdGhbxzZnBFM0XNrWsU27u3crYiOSJCOjUbmjNK5gnBZWLdT4wvEaXzBe4wrGMaE3gM+N4ClFVBZkabctUWl7vdelAAAAABhmvcFebe3Yqk1tm7S5fbOztW1WY1/j0D3paekaUzBG00qm6YLxF2hcwTiNLxivMfljlOlnpD0AdxA8pYjqgkyttiUq79rtdSmeu2PFHZKkWxbcEt+On7/V2Z979COu/vHpNZKkH1843Y2KAAAAkGCC4aC2dW7T5jYnXNrUvkmb2zZrd/duWVlJUoYvQ+MLxuv46uM1oXCCJhRM0LiCcRqVO0q+NJ/H3wBAqiF4ShGluRnapxId37fK61I8t751vTcd7/so5iZr93S6UAgAAABGunAkrPruem1ui4ZL0RFMOzp3KGRDkpwJvsfmj9X00um6eOLFqius08SiiarJrSFgAjBiEDyliLQ0o66MCmWH2qVgnxTI8rokAAAAAJI6Bzu1oXWDNrZt1IbWDdrQtkFb2rcMTfJtZFSTV6OJhRN1xugzNLFwoiYWTdTY/LFK96V7XD0AfDqCpxQymFMldUjq2C2VTvS6HAAAACClRGxE9V312tC2YShg2tC6QXt79g7dU5xZrElFk3TF5CtUV1SnusI6jSsYp+xAqi8RBCBRETylEJtf4wRPnfUETwAAAICLeoO92ti28ZBRTJvaNqk31CtJSjNpGps/VnPK5uhLk7+kyUWTNbl4ssqyymSM8bh6ABg+BE8pJFBUK+2SbEe9UvlH2Zj8Md50XDIh5ibjy3JcKAQAAADDqaWvRWtb1mpd6zqtb12vjW0btbNz59Bk37mBXE0qmqSLJ16sKcVTNLlosiYUTmAlOQApgeApheSUjZYk9TbvVCrHGbefcLs3HV/0q5ib/PTSWS4UAgAAgM/CWqvG3kata12ndS3rtLZlrda2rlVjb+PQPTW5NZpSPEXnjz9/aBRTdU41o5gApCyCpxRSUVygJpsvX/OOlA6eAAAAgL/EWqu9PXu1rmWd1rSsGQqbWvpbJDkTfo8rGKdjK4/V1OKpmlYyTVOKpygvPc/jygFgZCF4SiFVBZnaY0s1qr3e61I8dftbtzv7eI98eupmZx/DyKfbHvtQEiOfAAAA3GSt1a6uXVrbunZoJNO61nXqGOiQJPmMTxMKJ+ikUSdpaokTMk0umsyE3wBwFAieUkhVQZZW2RKN6d7jdSme2tG5w5uOW7bE3GRrU48LhQAAAKQua60aehu0pnmNPmr+SKtbVmtt81p1BbskSf40v+oK63TW6LM0rWSaphZPVV1RHfMxAcBnRPCUQkpy0tWgEmX1rfG6FAAAACAuOgY6tKZ5jVa3rNZHzR9pTfMaNfU1SZL8xq9JxZN07rhzNa1kmqaVTNPEwokK+AIeVw0AyYPgKYWkpRl1Z1YoI9gr9XdImQVelwQAAAAMm/5Qv9a3rndGMjWv1urm1drZtXPo+tj8sTq+6nhNL52umaUzNbl4sjJ8GR5WDADJj+ApxQzmVEvtkjrqCZ4AAACQsEKRkLa0b9Hq5uhIppY12tS2SWEbliRVZFdoRukMfaHuC5pROkPTS6Yz8TcAeIDgKcWYgppo8LRbqpjudTmemFI8xZuOK2fG3GRadb4LhQAAACSe1v5Wfdj0oT5o+kAfNH2g1c2r1RfqkyTlp+drRukMXTfzOs0omaEZpTNUll3mccUAAIngKeUEimqlHVKko15pXhfjkVsW3OJNx+f+LOYmP74wNcNBAACQ2kKRkDa2bdQHTR8MhU27unZJcuZlmlI8RV+Y+AXNKpulmaUzVZtXK2OMx1UDAI6E4CnF5JfXKGTTNNC8QzleFwMAAABIau5rPmQ009qWtUOjmcqyyjS7bLYun3S5ZpfN1rSSaawwBwAJhOApxVQW5qpBRcpq2ZmywdOtb9wqSfrZybGPQPpcHr3e2V9211E3+f6S9yRJv7hyrhsVAQAAxF0wEtTG1o1DIdMHTR9od/duSZI/za+pxVN1Wd1lml02W7PKZqkqp4rRTACQwAieUkxVQab22BKN76j3uhTPNPQ0eNNx556Ym+zt6HehEAAAgPjpHuzWh00f6t3Gd/Ve43v6sOlD9YedP+OUZ5VrdvlsfXnKlzW7bLamlkxllTkASDIETymmqiBTf7IlmtydusETAAAA3NPQ06D3Gt/Tu43v6v3G97WhbYMiNqI0k6YpxVP0xUlf1Ozy2ZpTNkeVOZVelwsAcBnBU4opzklXgylVdv8qyVqJYcsAAAD4jCI2oq3tW4dGM73X+N7Qa3NZ/izNLputG2fdqLnlczWrbJZyAqk62QMApC5XgydjzHZJXZLCkkLW2vnGmGJJD0oaK2m7pC9Za9uM8+L2LyWdJ6lX0jXW2nejz7la0t9HH/vP1tp7o+fnSbpHUpak5yR9z1pr3fxOic4Yo97MCvkHB6WeZimXZWYBAABwdAbDg1rTskbvNhwImjoHOyVJpVmlmls+V1dNvUpzK+ZqctFk+dP4PTcApLp4/CQ43VrbfNDnWyW9bK39mTHm1ujnWySdK6kuuh0n6X8kHRcNqn4sab4kK2mVMeYpa21b9J7rJb0tJ3haJOn5OHynhBbMqZYGJXXsSsngaXbZbG86rj025ibHjClyoRAAAICj0xvs1fuN72tlw0qtalil1c2rNRgZlCSNKxins8ecrbnlc3VM+TGqyathEnAAwMd48SuIiyWdFj2+V9IyOcHTxZLui45Y+rMxptAYUxW990VrbaskGWNelLTIGLNMUr619s/R8/dJukQET3+RKaiV2iR17pZGHeN1OXH3/Xnf96bjs26Puckti6YMfx0AAACfoHuwW+82vusETftWaW3LWoVsSD7j0/SS6frK1K9obvlczSmfo+LMYq/LBQAkALeDJyvpj8YYK+l/rbV3Sqqw1u6NXt8nqSJ6PErSroPa1kfPfdr5+iOcx1+QXlIrbZci7fVK87oYAAAAeKZjoEPvNjhB08qGlVrful4RG5E/za+ZpTN17YxrNb9ivuaUz1F2INvrcgEACcjt4Okka+1uY0y5pBeNMesPvmittdFQylXGmBsk3SBJo0ePdru7Ea+wtEr9NqBw806l4vSOf/Pq30iSfn76z+Pb8YNXOfsrfn/UTW763SpJ0q+/Ns+NigAAQIpp72/XqoZVQ0HThtYNsrJKT0vXzLKZun7m9ZpfOV+zy2Yry5/ldbkAgCTgavBkrd0d3TcaYx6XtEBSgzGmylq7N/oqXWP09t2Sag9qXhM9t1sHXs3bf35Z9HzNEe4/Uh13SrpTkubPn5/yk49XFWRpry1WYeuulAye2gfavem4ty3mJm29gy4UAgAAUkVrf6ve2feOVu5zgqbN7ZslSRm+DM0pm6NvzfmW5lfM16yyWcrwZXhcLQAgGbkWPBljciSlWWu7osd/Jeknkp6SdLWkn0X3T0abPCXpu8aYJXImF++IhlNLJf2rMWb/LMt/Jek2a22rMabTGHO8nMnFvy7p/7r1fZJJVWGm9thSFXbs+ss3AwAAIGF0DXZpVcMqvb33bb29721tatskScryZ2lO2RydO+5cza+YrxmlM5TuS/e4WgBAKnBzxFOFpMejK1v4Jd1vrX3BGPOOpIeMMddJ2iHpS9H7n5N0nqTNknolXStJ0YDpnyS9E73vJ/snGpf0bUn3SMqSM6k4E4sfheqCLK1TiWb3bPC6FAAAAHwO/aF+vdf4nlbsW6EVe1doTcsahW3YGdFUPkc3z71Zx1Yeq+ml0xVIC3hdLgAgBbkWPFlrt0r62Lr11toWSWce4byV9J1PeNbdku4+wvmVkmZ87mJTTGF2QI2mRNkDTVI4JPm8WNwQAAAAsQpGglrdvFpv731bK/at0PuN7ysYCcpnfJpZOlPXzbxOx1Uep9nls3l1DgAwIpA4pCBjjPqyqpQ2EJG690kFNX+5URI5ruo4bzoef2rMTU6cWOpCIQAAIFFEbEQbWjcMvTq3qmGV+kJ9MjKaUjxFX5nyFS2oWqB5FfOUE0jF2TsBACMdwVOKCuVWSwOSOupTLni6afZN3nR86g9jbnLzmXUuFAIAAEYqa612de3SW3ve0tt739Y7De+oY6BDkjSuYJwumnCRjq86XvMr5qsws9DbYgEAOAoET6mqoEZqkRM8AQAAwDMdAx16e+/bWr53uZbvWa7d3c5CzVU5VTq99nQtqFyg46qOU3l2uceVAgAQO4KnFJVVMlraKkU6divN62Li7KaXnBFPvz7r1/Ht+PeXOfurHj3qJlffvUKSdO83FrhREQAA8EAwHNQHTR8MBU1rWtYoYiPKDeRqQeUCXTP9Gp1QfYJq82oVXagHAICERfCUokpKS9Vps+Rr3qFUmw1gIDTgTcfB/pib9AfDLhQCAADiyVqrbZ3btHyPEzS9s+8d9YZ6hyYEv3HWjTqh+gTNKJ0hfxp/PAcAJBd+sqWoqoJM7bUlqmjd5XUpAAAASaetv01v731bb+15S8v3Lte+nn2SpNF5o3XhhAu1sHqhFlQuUF56nseVAgDgLoKnFFVVkKU9tkSVnczxBAAA8HkFI0F90PiB/rTnT3prz1ta17JOVlZ56Xk6vup43TDrBi2sWqiavNRa1AUAAIKnFFVdkKX3bInSe97zuhQAAICEtK9nn/60+096c/eb+vPeP6s72C2/8WtW2Sx9Z853tLB6oaaXTJcvzed1qQAAeOYzBU/GmBustXcOdzGIn/wsv5rSSpUVbHPmHgpkel1S3Jxac6o3HU86J+YmZ05l9RoAAEaKYDio9xrf05u739Qbu9/Q5vbNkqSK7AqdM/YcnTzqZB1XdZxy03M9rhQAgJHjs454YnmNBGeMUV92ldQvqXO3VDLB65Li5poZ13jT8Yk3x9zkhlNS598LAAAj0Z7uPXpz95t6c/ebenvv2+oN9cqf5te88nm6eN7FOmnUSZpQOIHV5wAA+ASfKXiy1v7vcBeC+IvkVjvBU0d9SgVPAAAAn2QwPKhVDauGwqatHVslSdU51bpg/AU6adRJWlC1QDmBVFsXGACAz+YTgydjzA8OO2UlNUt601q7zdWqEBemsNb5N9qRWhOMX/vCtZKkxYsWx7fjxedHC3j2qJtc8b/LJUkP3rjQjYoAAICk3d279Ub9G3pz95tasW+F+kJ9CqQFNL9ivi6tu1QnjzpZ4wrGMaoJAIDP4NNGPB1pbdexkv5/xpjbrbVL3CkJ8ZJVOkaRTUZq3a40r4sBAACIk1AkpA+aPtDr9a/r9frXh+Zqqsmt0cUTnNfnjq08VtmBbI8rBQAg8X1i8GSt/ccjnTfGFEt6SRLBU4KrKM7XHpWouGmr+GMVAABIZh0DHXpz95t6vf51vbn7TXUOdspv/JpXMU9fmP8FnVxzssbmj2VUEwAAwyzmOZ6sta2Gn8hJobIgUzsj5Sps2ep1KQAAAMPKWqst7Vv0Wv1rer3+db3f9L4iNqLizGKdVnuaTq05VQurFyov/UiD/AEAwHCJOXgyxpwuqc2FWhBn1QVZes+W65jO1V6XAgAA8LkNhAf0zr539Nqu1/TG7je0u3u3JGlK8RR9c+Y3dWrNqZpROkNphkkGAACIl0+bXPwjOROKH6xY0h5JX3ezKMTHqKIsPWkrlDmwTBroljJyvS4pLs4Ze443HU+/JOYmF8yqGv46AABIIg09DXpj9xt6rf41vb33bfWF+pTpy9TxVcfrupnX6eRRJ6syp9LrMgEASFmfNuLpgsM+W0kt1toeF+tBHOVm+NWaXu38m23fIVVM97qkuLhyypXedLzg+pibfG3h2OGvAwCABGat1ca2jXpl1yt6deerWte6TpJUlVOliyZcpFNrTtWxlccq05/pcaUAAED69MnFd8SzEHgjmD9G6pDUui1lgqe+UJ8kKcufFd+OB3udffrRT+XeNxiWJGWl+9yoCACAhBCMBLWqYZVe3fmqlu1apj09e2RkNKtslr53zPd0as2pmlg4kYnBAQAYgWKe40mSjDHPWGsPHxGFBOQrHecET23bvS4lbr790rclSYsXLY5vx3+43Nlf++xRN7lm8QpJ0oM3LnSjIgAARqzuwW69uedNvbrzVb2x+w11DXYpw5ehhVULdePsG3VKzSkqzSr1ukwAAPAXfJbJxU+S1OhCLfBAaVmlOjbnKK91m5hmEwAAeGlfzz4t27VMr+56VSv2rVAoElJRRpHOqD1Dp48+XQurFio7cPQjhwEAgPeOKngyxsyV9BVJl0vaJulRN4tC/IwuztZOW6a6pi1iJgQAABBP++drenXXq3p116ta27JWkjQmf4yumnqVTq89XbPLZsuXxivnAAAkqk9b1W6SpC9Ht2ZJD0oy1trT41Qb4qC2OFs7bIUmtm7zuhQAAJACQpGQ3m14dyhs2t29W0ZGM8tm6nvHfE9n1J6hcQXjmK8JAIAk8WkjntZLekPSBdbazZJkjPmbuFSFuBldnK0PbbnSu1dJkbDEbxQBAMAwGwgP6O29b+vFHS9q2a5lah9oV3pauhZWL9T1M6/XqbWnMl8TAABJ6tOCp0slXSnpVWPMC5KWSOJXT0mmqiBTu1Qhnw1JnbulwtFel+S6iyde7E3Hc74Sc5MvzqtxoRAAANzXE+zRG7vf0Ms7XtYbu99QT7BHeYE8nVp7qs4cfaZOqD6B+ZoAAEgBnxg8WWufkPSEMSZH0sWSvi+p3BjzP5Iet9b+MS4VwlV+X5p6c0ZLA3JWtkuB4OmSiZd40/Hcr8bc5PL5tS4UAgCAO9r727Wsfple3vGy3trzlgYjgyrOLNa5487VWaPP0oLKBQr4Al6XCQAA4ugvTi5ure2RdL+k+40xRXImGL9FEsFTkrBFY6R9klq3SeNO8boc17X1t0mSijKL4ttxT4uzzyk56iatPYOSpOKcdDcqAgDgc2vsbdQrO1/RSztf0sp9KxW2YVXlVOlLk7+ks8acpTllc5gcHACAFHZUq9rtZ61tk3RndEOSyCkbo9A+n/xt270uJS5+sOwHkqTFixbHt+OHvu7sr332qJt86/erJEkP3rjQjYoAAPhMdnXt0ss7XtZLO1/SB00fSJLG5o/VN2Z8Q2eOOVPTiqcxOTgAAJAUY/CE5FRTkqddkVLVtmzlfxAAAOCINrdt1os7X9TLO17WhrYNkqSpxVP13Tnf1dljztb4wvEeVwgAAEYicgZodHG2dtoKVTYTPAEAgAM2t23WH3f8UX/c/kdt6dgiI6O55XP1f+b/H50x+gzV5LEIBgAA+HTkDFBtUbY+tOU6of0dr0sBAAAe2x82Ld2+VFs7tsrIaF7FPP1oyo901uizVJZd5nWJAAAggRA8QaOLs/W0rVAg2CH1tUlZcZ50GwAAeGpz22Yt3bFUf9z+x0PCpiunXEnYBAAAPheCJ6gwO6Amf5XzoW170gdPV0y+wpuOj/1GzE2uOn6MC4UAAFKdtVab2w+8Rndw2PTlKV/WWWPOUmlWqddlAgCAJEDwBBljFMwfI3XJCZ6q53pdkqsWjVvkTcczLou5yYWzq10oBACQivaHTUu3L9Ufd/xR2zq2ychofuV8wiYAAOAagidIkvylY53gqXWb16W4bl/PPklSZU5lfDvuqHf2BUc/Eeue9j5JUnVhlhsVAQBSwKa2TYeETWkmTfMq5ukrU75C2AQAAFxH8ARJUkVpqZq3FqikbbuM18W47LY3bpMkLV60OL4dP3ajs7/22aNu8jcPvi9JevDGhS4UBABIVrs6d+n57c/r+W3Pa3P7ZsImAADgGYInSHImGN9py5TftEXpXhcDAABi1tDToKXbl+r5bc9rdctqSdLc8rn60XE/0tljziZsAgAAniB4giSppjhbO2yFprdt97oUAABwlNr62/Tijhf1/LbntaphlaysphZP1Q/m/UCLxi5SVW6V1yUCAIAUR/AESc6Ipw9tudK7l0uhQcnPuCcAAEai7sFuvbLrFT237Tn9ec+fFbZhjSsYp2/N+ZYWjV2kcQXjvC4RAABgCMETJEmjCrO001bIKCJ17JJKJnhdEgAAiOoP9ev1+tf1/Lbn9Xr96xqMDKo6p1pXT79a5407T5OKJsmYZJ+lEQAAJCKCJ0iSMgM+9WTXSEFJbduSOni6evrV3nR8wndjbnL9yeNdKAQAkAiC4aCW712u57c9r1d2vqLeUK9Ks0p1+eTLtWjsIs0um03YBAAARjyCJwyxheOkJkmt27wuxVWn1Z7mTceTz425yVnTKlwoBAAwUkVsRO83vq9ntz6rpTuWqmOgQ/np+Tp33Lk6d9y5ml8xX740n9dlAgAASbJWshEpEnb2NnzQceSga+HD7jvoWmaBlJ/cczISPGFIXukoDTSlKyPJJxjf1uEEa3GfA6N5k7MvrTvqJluauiVJE8py3agIADBCbG3fqme2PqPntj2n3d27lenL1OmjT9f5487XCdUnKOALeF0iAADuCoekUJ8z53Co39nC+48HDmzh/ccHnz/s3vCgFA5KkVB0v/845Bwfcm3/uU+4ZsNOwHSkcEn283/veddIF/7y8z9nBCN4wpDaklztiJRpQss2JfPvUn+y/CeSpMWLFse346e/7+yvffaom/zosY8kSQ/euNCFggAAXmrqbdJz257Ts1uf1brWdUozaVpYtVDfmfMdnTn6TGUHsr0uEQCAAyIRKdgjDXRLgz3SYNeB42BvdOs7sB/siX7uO+h674Fzh1+PBD9/jWkByZ8p+QLOlhaQfH5nn+Y/cLz/Wnr2J1/z+Z3zxiel+SSTdmBL8znnh47NYZ9juFac/NOrEDxhyOiSLO205RrTsjWpgycAALzSE+zRSzte0rNbn9Xb+95WxEY0vWS6bjn2Fi0at0ilWaVelwgASCaRiDTYLQ10Sv2dUn/HgeOBTufa4P4wqeug425pIPp5sNs5F+yJrW9/phTIjm5Z0S1bysiTcisOPbf/Hn9mdEuP7jMkX4az92cc4dxB9/oypLQ0d/454nMheMKQ2qJsfWgr5Ot43RlKyISlAAB8bsFIUG/tfkvPbH1Gy3YtU3+4X6NyR+n6mdfr/PHnx//VbwBA4oiEnbCor03qa3f2Ax3Ouf3h0dD+sHP7Q6ajeR0sPVdKz3H2GblSep6UVxU9znXCosOvp+dEj3OkQI4zemgoQMoiBMIQgicMGV2crWdsufyhXqmnWcot87okAAASkrVWHzR9oGe2PqOl25eqfaBdhRmFunjixbpg/AWsSAcAqcRaZ+RQX5uz9bcfGiR94ucOJ2T6NMYnZeZLGfnRfYFUONqZsHro3MH7ggPXMvKc4CiQQ0gEVxE8YUhZXob2pFU6H9q2ETwBABCjHZ079PSWp/Xs1mdV312vDF+GTq89XReMv4BJwgEgWQz2Sr0tB22tB477Wo98Pjz4yc9LC0hZhVJWkZRZKOVWSmVTDnzOKjr0embBgSApPYc3VTDiETxhiDFG4fwxUo+ktu1S7QKvS3LFDbNu8KbjU/4u5iZ/fcbRr4AHAPBGx0CHlm5fqqe2PKUPmj5QmknTgsoFumn2TTpz9JnKTWdlUgAY0YJ9Unej1NPkbN2NUk+j1N0k9TYfFCJFg6RQ3yc8yDjhUHaJsxWOkarnOsdZRVJ28ZGDJMIjJDmCJxzCVzLOCZ5at3ldimsWVnu0QtyE02NuclIdk8wCwEgUioT01p639OTmJ7Vs1zINRgY1sXCifjDvBzp//Pkqzy73ukQASF3WOnMbdTc5AdJQmNR0hICpyZk8+0gy8p3QKKdUyq+WKmc64dH+YCnroOPsEidMSmOZJuBwBE84RHVpoRp2FKu8bZuSNXNf37pekjSleEp8O977obOvmnXUTdbscd7pnl5d4EZFAIAYbWjdoCe3PKlntz6r1v5WFWYU6ouTvqiLJl6kacXTmLcJANw02CN17ZO6G6SuvVJXw0HBUtOBUUo9TVJ44AgPME5wlFPuTCsy6pgDxznlUk7ZoceBzLh/RSAZETzhELXF2dpuy1XcvFXJOgvFHSvukCQtXrQ4vh2/cJuzv/bZo27yk6fXSpIevNGjUVoAADX3NevZrc/q6S1Pa0PbBvnT/Dq15lRdNOEinTzqZOZtAoDPa6DLCZG69kZDpX2HHUfDpoHOj7dNCzghUU6plFsulU09cDwUKkXDpOwSycdfgYF44/91OMTo4mztjJTrmLYNXpcCAIBnBsIDenXXq3p6y9P60+4/KWzDmlEyQz867kc6d+y5Ksws9LpEABj5wiGpe5/UuUfq3B3d74mOVtp3YAv2fLytP1PKq3Qm2q6YJk08U8qtkPKqpLwK53xepTNPEqNNgRGN4AmHqC3O0ge2QoHe151J9gJZXpcEAEBcWGv1QdMHemrLU3ph+wvqGuxSeXa5rp5+tS6acJEmFE7wukQAGDmC/VLXngNh0tB2UMDU0yjZyKHt/FlOYJRXJVXNliadc2iglFflfM4sIFACkgTBEw5RW5StnTY6IWrbDqk8zvMgAQAQZ/t69umpLU/pqS1PaUfnDmX6MnXmmDN10YSLdFzlcfIxUSyAVBMOOaOSOuqljl3RrV7q2H0gbOpt+Xi7jAJnEu78KmeUUv6o6OdRTqCUX80IJSAFuR48GWN8klZK2m2tvcAYM07SEkklklZJ+pq1dtAYkyHpPknzJLVIusJauz36jNskXScpLOlma+3S6PlFkn4pySfpN9ban7n9fZJdToZfHVmjnH/SbdsIngAASWkwPKhXd72qxzc/ruV7litiI5pXMU/XzbhOZ485W7npuV6XCADu6e88KEzav6+X2qPHXXs+PlIpu+RAiFRz7IHj/GopLxo2ZeR5830AjGjxGPH0PUnrJOVHP98h6efW2iXGmF/LCZT+J7pvs9ZONMZcGb3vCmPMNElXSpouqVrSS8aYSdFn/ZeksyXVS3rHGPOUtXZtHL5TUrNF46RmSa3bvC7FFd875nvedHzmP8Tc5IeLJrtQCACkrvWt6/XE5if0zNZn1DHQoYrsCn1z5jd1yYRLVJtf63V5APD5WeuMRmrbIbVvl9p3HhoqddRLAx2HtkkLSAWjpIJaadzJzr6gJrrVOtfSczz5OgASn6vBkzGmRtL5kv5F0g+Ms8bwGZK+Er3lXkm3ywmeLo4eS9Ijkv5f9P6LJS2x1g5I2maM2SxpQfS+zdbardG+lkTvJXj6nAqLK9XdnK3ctuQMnuaUz/Gm49HHxdxk3phiFwoBgNTSMdChZ7c+qyc2P6F1resUSAvozNFn6gsTv6DjqniVDkACGuiW2ndEw6WD99udoGmw+9D7s4qcEKlojDT2xIOCpeg+t0JKS/PkqwBIfm6PePqFpB9K2j/mskRSu7U2FP1cL2lU9HiUpF2SZK0NGWM6ovePkvTng555cJtdh52P/W/2+JjRJTnaHqnQtJYtSsYfP+83vi/JgwBq59vOPoYAatWOVkkEUAAQq3AkrLf3vq3HNz+ul3e+rGAkqKnFU3Xbgtt0/vjzVZBR4HWJAPDJQoPOK3Bt248cMB0+v1J6rlQ4RioaK4071QmYCsdE96N5BQ6Ap1wLnowxF0hqtNauMsac5lY/R1nLDZJukKTRo0d7WUpCGF2crW22QpObkzN4+uW7v5QkLV60OL4dv/wTZ3/ts0fd5N9e2CBJevDGhW5UBABJZ1fXLj25+Uk9ueVJ7evZp/z0fF0+6XJdMvESTS2Z6nV5AOCwVuppklq2OPOqHh4sde6RZA/cnxaQCmudMGnqhQeFSmOdsCm7mAm7AYxYbo54OlHSRcaY8yRlypnj6ZeSCo0x/uiopxpJu6P375ZUK6neGOOXVCBnkvH95/c7uM0nnT+EtfZOSXdK0vz58+2R7sEBNcVZWmkr5e98RwoHJV/A65IAAPhEfaE+vbTjJT2x+Qmt2LdCRkYnjDpBfzv/b3V67enK8GV4XSKAVBSJSN37pNatR9i2HfY6nHEm6S4cI4075aBgKbrPq5J4LRhAgnIteLLW3ibpNkmKjnj6O2vtV40xD0v6opyV7a6W9GS0yVPRz8uj11+x1lpjzFOS7jfG/KecycXrJK2QZCTVRVfJ2y1nAvL9c0fhcxhdnK1HI5UyNuy8I14yweuSAAA4hLVWa1rW6LFNj+n5bc+rO9itmtwa/fXcv9ZFEy5SZU6l1yUCSAWRiNS5+5PDpVDfgXvTAs7opOLx0piTnH3xeKl4nDPPkp+QHEByiseqdoe7RdISY8w/S3pP0m+j538r6XfRycNb5QRJstauMcY8JGfS8JCk71hrw5JkjPmupKWSfJLuttauies3SVJVBVnalVblfGjZQvAEABgxuga79OzWZ/XIxke0oW2DMn2Z+quxf6VLJl6ieRXzlGaS8SVxAJ4Kh6TO+gOBUstB4VLbdik8cOBeX4YTJBWPlyacceC4eLwzkTejlgCkoLgET9baZZKWRY+36sCqdAff0y/p8k9o/y9yVsY7/Pxzkp4bxlITwvaO7Xqv8T19oe4Lrjzfl2YUyh8r9Upq3eJKHwAAHC1rrd5vel+PbHxEf9z+R/WH+zW1eKr+/ri/13njz1NeOpPmAvicrHUm7G7eJLVskpo3Ss2bneO2HVIkeOBef5YTJJXWSZPOORAslUyQ8qpZHQ4ADuPFiCd8Tq/Vv6b/WPkfOrX2VBVnurPaWUFplXp2ZSundasrz/fSLQtu8abjRT+Nuck/XDjNhUIAIDG097frqS1P6bFNj2lLxxZl+7N14YQLddmkyzS9ZLrX5QFIRKFBZzLv5mi41LL5wHF/+4H7fOlS8QSpfKozmffQa3ETpLxKJvIGgBgQPCWgSUWTJEmb2jbpuKrjXOljTGmutu+o1LSWLUq2H6tTiqd403HVrJibTK9muW8AqSViI3pn3zt6dOOjemnnSwpGgppVOkv/eMI/atHYRcoOZHtdIoCRzlqppzkaLG2KjmLa7Hxu2yE5s3Y4ciudkUvTvyCVTnKOSyZKhaN5LQ4AhgnBUwKqK6qT5HLwVJKtrZFyTWnZomT7kbt8z3JJ0sLqhfHteMurzn7C6Ufd5M1NzZKkk+pK3agIAEaM5r5mPbn5ST266VHt6tqlvPQ8XT7pcl026bKhX7gAwCFCg860EEOvx206cNzfceA+f6YzUqlyljTjMqmkTiqd6Owz872rHwBSBMFTAirNKlVxZrE2tm10rY+xJTn6yFYqreMdKRyUfAHX+oq3Oz+8U5IHwdPr/+HsYwie/u8rmyQRPAFITuFIWMv3LtejGx/Vsl3LFLIhzauYp2/N/pbOHnO2Mv2ZXpcIYCQIDTiBUtN6qWmD1LTO2bdsOXT0Ul61EyjN+KIzcqm0zgmXCmqZdwkAPETwlKDqCuu0qW2Ta88fXZKtpyOVMjbsDEkunehaXwCA1NLQ06DHNj+mxzc9rr09e1WUUaSrpl2lS+su1biCcV6XB8Arwb5owHRQuNS03lk9zkace4zPmWupbLI09SJnv//1uAwWGgCAkYjgKUHVFdXp0U2PKmIjriwdXVOUpZ2qdD60biV4AgB8LhEb0fI9y/XQhof0Wv1rCtuwjq86Xj+Y/wOdUXuG0n3pXpcIIF4Ge535lg4PmNq2HxowlUyQyqdJ0y91AqbyqU7A5M/wtHwAQGwInhLUpKJJ6gv1qb6rXqPzRw/78zP8PvXnj5X65bw7DwDAZ9Da36onNj+hhzc8rPruehVnFuvq6Vfri3VfVG1+rdflAXDTYE80VDo8YNohyTr3pPmdMKlyljTzSwcCpuIJkp9AGgCSAcFTgto/wfjGto2uBE+SVFhSpZ492cpp3erK8wEAyclaq1UNq/TQxof00g5nZbp5FfN08zE368zRZzK6CUg24ZAzQr5xjdSwVmpcKzWscUYwDQVMAeeVuOq50uyvOAFT2RRnVFMSzSUKAPg4gqcENaFwgoyMNrVt0lljznKlj9GlOdqxu1LTWpJrxNM/LPwHbzq+8BcxN/nXS2cOfx0A4JLOwU49veVpPbzhYW3p2KK8QJ6+NPlLunzS5ZpQOMHr8gB8XtZK3Y1Sw+pouLTWCZuaNkihfucek+aMVqqaJc3+sjN6qWyKMy+Tj796AEAq4r/+CSrLn6XavFptandvgvGxJdnaHK7Q5JYt8rnWS/x5NnFtaV3MTSaU5bpQCAAMrzXNa/TQxof0/Lbn1Rfq04ySGfrJCT/RonGLlOXP8ro8AJ/FYI/UuP6gkGmNs+9tOXBPboUzB9Ox33T2FdOckCnA/+8BAAcQPCWwuiJ3V7YbU5KjDbZCaR1vS6HBpHnPftmuZZKk02pPi2/HG5539pPPPeomL61tkCSdNa3CjYoA4DPrDfbq+W3P66GND2lty1pl+bN03rjzdPnkyzW9ZLrX5QE4WtZKnbulfR9J+1ZL+z50jg9+TS6Q7YxcmnyeVDE9GjJNl3JKvawcAJAgCJ4SWF1RnV7d9ar6Qn2u/EZ5bEmOlkYqZWxEat+ZNCvb3bvmXkkeBE9v/T9nH0PwdNcbzvxaBE8ARopNbZv00IaH9MzWZ9Qd7NbEwon60XE/0gXjL1BeOkuZAyNaOOisJrfvo+gWDZn62g7cUzxeqpzpvCZXEQ2YCsdKacO/ijIAIDUQPCWwSUWTFLERbW3fqumlw//b5dHF2dpuK50PrVuSJngCAMQmGA7qpZ0vacn6JXq38V0F0gI6Z+w5+tLkL2lO2RwZY7wuEcDh+juiI5gOCpma1kvhQee6P9MZuTT1IidoqpzlBE0ZBMgAgOFF8JTA6goPrGznRvCUle5TT+5oKShnpRIAQEpp7G3Uwxsf1iMbH1FzX7Nqcmv0g3k/0CUTL1FRZpHX5QGQnFflOnYdFDBFQ6b2nQfuyS51Jvue8K1owDRDKpnIZN8AgLjgp00Cq82rVaYv09UJxgtKqtTTkKOcJFvZDgBwZNZarWxYqSXrl+iVna8obMM6adRJunLKlTpp1ElKM7xuA3gmHHRGLe2NviLXEJ2Tqb8jeoNxAqVR86V510ZHMs10JgFnZCIAwCMETwnMl+bT+MLxrk4wPrY0RzsbKjWVEU8AkNR6g716ZuszemD9A9rcvln56fn66tSv6orJV6g2v9br8oDUExp0VpHb+4G0931pz/vOynLhAed6INuZf2n6pYe+Kpee42XVAAB8DMFTgqsrrNObu9907fljSnK0KVSuyS1blCy/4/7pyT/1puNL/zfmJj+/Ys7w1wEAB9nWsU0PbnhQT25+Ut3Bbk0tnqqfnPATLRq3yJWFKwAcQWjACZX2vu8ETXved0Kn/fMxZRQ4r8otuF6qnuuETCUTpDSfl1UDAHBUCJ4S3KSiSXpyy5Nq6WtRSVbJsD9/TEm2ttoKmY63nd+8+dOHvY94q8yp9KbjgpqYm1QX8pc+AMMvFAnptfrXtGT9Ev1575/lT/PrnLHn6MrJV2p22WwmCwfcFOw7EDLted/ZN66TIiHnemahVDVbOv5bUtUc57hoHKvKAQASFsFTgqsrciYY39S+yZXgaWxJjl6JVMrYiDNJZRKsbPfCthckSYvGLYpvx6sfdfYzLjvqJk9/sEeSdOHsajcqApBiWvtb9dimx/TQhoe0t2evKrIr9Ndz/1qX1l2q0qxSr8sDkk9owJmHafe7h4ZMNuxczyqWqudIJ5ztBEzVc6TCMczHBABIKgRPCW4oeGrbpOOrjh/2548uydZ2Gx0h1LolKYKnBzc8KMmD4Omdu519DMHT7/+8QxLBE4DP56Omj/TA+gf0wvYXFIwEdVzlcfrhsT/UabWnyZ/GHwWAYREJS80boyHTu86+YfWB1+WyS5zX5CYtcgKmqtlSQS0hEwAg6fGnzQRXmlWq4sxi1yYYz88MqD2rVopIYoJxAEgYwXBQS3cs1f3r7tdHzR8p25+ty+ou05VTrtSEwglelwckNmultu0HAqY97zlzMw12O9fT85xw6bibpFHHSKPmETIBAFIWwVMSqCusc3Vlu4KSSvU05yinZYtrfQAAhkdzX7Me3viwHtrwkJr7mjU2f6xuW3CbLp54sXICrHYFfCZd+w4dybTnPamv1bnmy3BWlZvzFan6GCdoKqljTiYAAKIInpJAXVGdHtn4iMKRsHwurG4ytjRXO5urNLWV4AkARqq1LWv1h3V/0PPbnlcwEtRJo07SV6d+VSdUn6A0w1+AgaM22OPMx7R7pVT/jhM0de52rpk0qWyqNOV8J2CqPkYqn5YUi68AAOAWgqckMKlokvrD/arvrteY/DHD/vzRxdnaFCrTlJatYoA4AIwcwUhQL+98Wfevu1/vNb6nLH+Wvjjpi/rylC9rXME4r8sDRr5IRGrZJNXvD5lWSg1rD0z+XTRWGn2886pc9TFS1SwpnZGDAADEguApCRw8wbgbwdPY0mxts5VSx9tSaDDhf6v3n6f9pzcdf+m+mJv8z1XzXCgEQKJr62/To5se1ZL1S9TQ26Ca3Br98Ngf6pKJlygvPc/r8oCRq6flwEim+pXOaKaBDudaRr4ziunkH0ij5ks186UcVnsEAODzInhKAhMKJ8jIaFPbJp015qxhf/6Ykhy9EamQsRGpfWfCr2xXlFnkTcc5JTE3Kc5J7JAPwPDa0LpB96+/X89ufVYD4QEdV3Wc/v74v9fJo0525VVrIKGFBqV9Hx0YyVS/Umrb5lwzaVLFdGnGpU7AVHMs8zIBAOASgqckkOXPUm1erTa1uzPB+NiSHG23lc6H1i0JHzw9sfkJSdIlEy+Jb8fv/cHZz/3qUTd5eOUuSdLl82vdqAhAAghHwlq2a5n+sP4PemffO8r0ZeqiCRfpK1O+oolFif3fY2DYWCu174i+MrfSCZr2fiCFB53redVSzTxp/rXOaKbqObwyBwBAnBA8JYm6ojptbNvoyrOLsgNqSq9xPiTBynZPbn5SkgfB0/v3O/sYgqdHVtVLIngCUlHnYKce2/iYHlj/gPb07FFVTpV+MO8HurTuUhVkFHhdHuCtwV5p9ypp19sH5mfqbXau+bOk6rnScTc5o5lGzZcKRnlbLwAAKYzgKUlMKpqkV3a+or5Qn7L8WcP6bGOMCkoq1NOWo5zWrcP6bADAoXZ27tTv1/1eT2x+Qn2hPs2rmKf/c+z/0Wm1p8mfxo9tpKiO3U7ItGuFtOvPzit0kZBzrXSSNOkcZwLwmmOdVeZ8/H8FAICRgp/KSaKuqE5WVlvbt2p66fRhf/6Y0lztaq/SlNbEH/EEACONtVarGlbpvrX3admuZfKl+XTeuPN01dSrNLVkqtflAfEVDkkNq6NBUzRs6nBePVcg2wmYTvy+VHucVHuslOXR3I0AAOCoEDwlibpCZ2W7jW0bXQmexpZka9O6Mk1u3Soz7E8HgNQUDAf1wvYX9Lu1v9O61nUqzCjU9bOu15WTr1RZdpnX5QHx0dfuvCq3P2iqXyUFe5xr+aOcgGnhd6XaBVLlTMkX8LRcAAAQG4KnJFGbV6tMX6ZrE4yPKc7RVlsptb/trBLjZ7U1APisOgY69PDGh/XAugfU2NeocQXj9A8L/0EXjL9g2F+XBkYUa6XWrQdCpp1vS03rJVnJ+KTKGdLcq5yQafTxUkGN1xUDAIDPieApSfjSfBpfON61CcbHlGRreaRSxkacVWNK61zpJx7++6z/9qbjrz4cc5N7rl3gQiEAvLK9Y7t+v+73emrLU+oL9Wlh1ULdfsLtOnHUiUozLOOOJBTsl/a+fyBk2vX2gUnAMwqcgGnGZc5+1DwpI9fTcgEAwPAjeEoik4om6fX611159tjSHO2wFc6H1q0JHTx5NpogPTvmJlnpPhcKARBP1lq9s+8d3bf2Pr1W/5oCaQFdMP4CXTXtKk0qmuR1ecDw6u9w5mTa8Zaz7XlXCg8614rHS3V/dWA0U+lkKY3AFQCAZEfwlETqCuv0xOYn1NzXrNKs0mF9dnlehvb4q50PLYk9wfiS9UskSVdOuTK+Ha+4y9kvuP6om/xu+XZJ0tcWjh3+egC4ajA8qOe3Pa/frf2dNrRtUHFmsb41+1v60uQvDft/owHPdDdJO986EDQ1rJZsRErzS1VzpONulEYvlGoWSLnMWwYAQCoieEoiU4qnSJLWt67XSaNOGtZnG2NUUFSp3u4cZSf4ynZLty+V5EHwtOYJZx9D8PTMh3slETwBiaStv00PbXhISzYsUXNfsyYWTtQ/nvCPOn/8+crwZXhdHvDZWSu175R2Lpd2/EnasVxqic4t6c+SauZLp/xQGrNQqjlWSs/xtl4AADAiEDwlkf1Lbq9pXjPswZMkjSnN0faeGk1r2jDszwaARLezc6fuW3ufntz8pPrD/Tpx1In6l6n/ooXVC2UM64EiAVkrNW04aETTcqmz3rmWWeCMZJp7lTTmRKlqNguPAACAIyJ4SiJ56Xkamz9Wq1tWu/L8saU5WrO5WlMbPxJ/hQIAxwdNH+jeNffqpR0vyZ/m14UTLtTXp31dEwoneF0aEJtwSNr3oRMy7VzubL0tzrXcCmnMCdLo7zn78mnMzwQAAI4KwVOSmV46Xe/sfceVZ48uztb68CiZ3lelnmYphzlKAKSmiI1o2a5lunfNvXq38V3lpefpmzO/qS9P+bLKspnHBgki2C/tXhUNmt5yJgUf7HauFY2TJi1yRjWNOcGZGJyRewAA4DMgeEoyM0pm6Nmtz6qxt1Hl2eXD+uyxJTlaamucD43rpHEnD+vzAWCkGwgP6OktT+veNfdqe+d2jcodpVsX3KovTPyCsgOxr1wJxFVoQKpfKW1/Q9r+phM0hQeca+XTpdlXRkc1nSDlV3lbKwAASBoET0lmeul0Sc48T+Wjhzd4GlOSrY2RaPDUtD5hg6fFixZ70/G1z8bc5MEbF7pQCIBYtfe368END+r+9fertb9V00qm6d9P+XedNeYs+dP4UYoRKjTgjGja/qa07XWp/h0p1C/JSFWznMUuxpwojT5eyi72uloAAJCk+NNykplSPEVpJk2rW1br9NGnD+uzqwuz1OorVr8vV5mN64b12QAwEu3q2qXfrf2dntj8hPpCfTp51Mm6dsa1ml8xnwnDMfKEBqU97zojmra94YxoCvVJMlLlDGn+ddLYk5xV57KKvK4WAACkCIKnJJPlz9KEwgla07Jm2J/tSzMaXZyj+sExmti0ftifHy/3rL5HknTNjGvi2/GffuXsT7z5qJvc+foWSdINpzBJMRBPq5tXa/HqxXpp50tKM2m6YPwFunra1ZpYNNHr0oADwkFpz3vOaKbtb0q73paCvc61ihnSvGuiQdMJjGgCAACeIXhKQjNKZujVXa/KWjvsv5GvK8/T+p2jNLFxpbPMcgL+xv+1+tckeRA8bVzq7GMInl5e1yiJ4AmIh4iN6I36N7R4zWKtalilvECerp1+rb4y9SvDPmce8JlEwtLe952gadsb0s4/S8Ee51r5dGnu15ygaexJBE0AAGDEIHhKQtNLpuvxzY9rT88ejcodNazPrqvI1XsbKnWBv1XqaZJy+csYgMQWDAf1zNZndM+ae7S1Y6uqcqr0w2N/qEvrLlVOIMfr8pDKrJWaN0pbX5O2vea8Qtff4VwrmyrN+Yoz3+KYE1lpFgAAjFgET0loRukMSc6rIsMdPE0sz9VDkegzG9cRPAFIWD3BHj2y8RHdt/Y+NfY2akrxFN1x8h06e+zZCqQFvC4Pqaqj/kDQtPU1qXufc75wtDT1Imn8adK4U/j5CwAAEgbBUxKqK6qTP82vNS1rdM7Yc4b12RPLcw9d2W78qcP6fABwW2t/q/6w7g96YP0D6hrs0oLKBfqnE/5JC6sXMmE44q+3NfrqXDRoanXm9lN2qRMwjT9VGneqVDzO2zoBAAA+I4KnJJTuS9fkosla0zz8E4xPKMtVsylUvz8/YVe2y/BneNNxIDPmJpkBnwuFAKmpvqte9665V49vflyD4UGdOfpMfWPGNzSzbKbXpSGVDPZIO5ZL25Y5QdO+jyRZKT3XeWXu2OucoKl8mpSW5nW1AAAAnxvBU5KaXjJdz217ThEbUZoZvj+4ZgZ8zsp2kTGamKDB06/P+rU3HV/1aMxN7v3GAhcKAVLLhtYNunv13Vq6famMMbpowkW6Zvo1GlfACBLEQWhQ2r3qwIim+nekSFDypUs1C6TTf+QETaOOkXy84gkAAJIPwVOSmlE6Qw9tfEg7OncM+1+u6spztWH3KE1sWp6wK9sBSG7WWq1sWKm7V9+tN3e/qWx/tr427Wu6aupVqsip8Lo8JLNIRGpYfSBo2vFWdOU5I1XNlhZ+2wmaRi+U0rO9rhYAAMB1BE9JalrJNEnSmpY1wx48TSzP06rNFTrf1yF17ZPyq4b1+W779QfOiKebZt8U345f+zdnf+oPj7rJr17eJEm6+cw6NyoCkk7ERvTqrld19+q79WHThyrOLNbNc2/WlyZ/SQUZBV6Xh2TVuVfa+qq05RVpy6tSb7NzvnRSdOW5U6SxJ0nZxd7WCQAA4AGCpyQ1oXCCMn2ZWtO8RheMv2BYn11XnqtHwjWST1LTuoQLnt7e+7YkD4Knra85+xiCpz9tdv7yQvAEfLpgOKhntj6jxWsWa1vHNtXk1ujvj/t7XTzxYmX6Y59fDfhUg73OSKYtrziBU+Na53xOuTTxTGn86c6k4PnV3tYJAAAwAhA8JSl/ml9TiqdoTcvwTzBeV3HQynaN66UJZwx7HwBwNHqDvXp448O6b+19auxt1JTiKfq3U/5NZ485W/40fsRhmOx/fW7LK862c7kUHpR8GdKYhdLsK52fheXTmRAcAADgMPypPInNKJ2hRzY+olAkNKx/AZtQlqsWFagvUKispsScYBxAYusY6ND96+7X79f9Xp2DnVpQuUA/OeEnOqH6BBnmncNw6NrnvDa3f1RTT5Nzvny6tOAGacLp0ugTmKcJAADgLyB4SmLTSqapP9yvrR1bNalo0rA9NyfDr1GFWdptxmhi4/phey4A/CVNvU26b+19emjDQ+oN9er02tP1zZnf1KyyWV6XhkQX7Dvw+tyWV6XG6Ijh7FJnNNOEM6TxpyXc6+UAAABecy14MsZkSnpdUka0n0estT82xoyTtERSiaRVkr5mrR00xmRIuk/SPEktkq6w1m6PPus2SddJCku62Vq7NHp+kaRfyplt6DfW2p+59X0S0YzSGZKkNc1rhjV4kqSJ5bla31CjiU1vJNzKdoUZhd50nF0Uc5Oi7HQXCgESz+7u3Vq8erEe3/S4QjakRWMX6bqZ1w37f9uQQqyVGtYceH1ux1tSeEDypTsrzp31j07YVDGD1+cAAAA+BzdHPA1IOsNa222MCUh60xjzvKQfSPq5tXaJMebXcgKl/4nu26y1E40xV0q6Q9IVxphpkq6UNF1StaSXjDH7/6bxX5LOllQv6R1jzFPW2rUufqeEMiZ/jHIDuVrTskZfqPvCsD67rjxXq7ZV6AJfp9S5RyoYNazPd9PPT/+5Nx1f8fuYm/z6a/NcKARIHFvbt+q3q3+rZ7c+K2OMLp5wsa6bcZ1q82u9Lg2JqKdZ2vzygdfnuhuc82VTpWO/6QRNY3h9DgAAYDi5FjxZa62k7ujHQHSzks6Q9JXo+Xsl3S4neLo4eixJj0j6f8aZqONiSUustQOSthljNktaEL1vs7V2qyQZY5ZE7yV4ikozaZpWMk2rm1cP+7PrKnL1WGjUgZXtEih4AjDyrW1Zq9989Bu9tOMlZfgy9OUpX9bV069WZU6l16UhkUTC0u53pc0vSptelPa8J8lK2SXOynMTznDmamL1OQAAANe4OseTMcYn53W6iXJGJ22R1G6tDUVvqZe0P7EYJWmXJFlrQ8aYDjmv442S9OeDHntwm12HnT/Oha+R0KaXTNfv1v1Og+FBpfuG77WtieV52mij/xoa10sTzxq2Z7vtF6t+IUn6/rzvx7fjl2539mfdftRN7njBmUPrlkVThr8eYARa1bBKd310l/60+0/KC+TpmzO/qaumXaXizGKvS0Oi6G6StrzsBE1bXpH6WiWTJo2aL53+I+fnVdUcXp8DAACIE1eDJ2ttWNIcY0yhpMclefK3Z2PMDZJukKTRo0d7UYJnppdOVygS0qa2TZpeOn3YnjuxPFdtyldvoFjZCbay3QdNH3jT8a53Ym7y7o42FwoBRhZrrf6050+668O79G7juyrOLNb3jvmerph8hfLS87wuDyNdJCztXuUETZv3j2qSlFMmTTrHCZomnCFlE14CAAB4IS6r2llr240xr0paKKnQGOOPjnqqkbQ7ettuSbWS6o0xfkkFciYZ339+v4PbfNL5w/u/U9KdkjR//nw7LF8qQUwvccKmNS1rhjV4KsgKqDwvQ3vSWdkOwGcTsRG9vPNl3fXhXVrXuk4V2RW6dcGturTuUmX5s7wuDyNZd6MzV9Pm/aOa2pxRTTXHSqf/vVR3llQ5m1FNAAAAI4Cbq9qVSQpGQ6csOZOA3yHpVUlflLOy3dWSnow2eSr6eXn0+ivWWmuMeUrS/caY/5QzuXidpBWSjKS66Cp5u+VMQL5/7ihEjcodpcKMQq1uXq0vTf7SsD67riJX65trNLFpWcKtbAfAO6FISM9te06/+eg32taxTWPyx+gnJ/xEF4y/QAFfwOvyMBJFwlL9ygNzNe193zmfUy5NOtcJmsafzqgmAACAEcjNEU9Vku6NzvOUJukha+0zxpi1kpYYY/5Z0nuSfhu9/7eSfhedPLxVTpAka+0aY8xDciYND0n6TvQVPhljvitpqZwpru+21q5x8fskJGOMppdM15qW4f9HU1eep1U7K3RBWpfUUS8VssoUgE8WDAf19NanddeHd6m+u16Tiibp30/5d5095mz50nxel4eRprtR2vzSgbma+tsPG9V0tlQ5i1FNAAAAI5ybq9p9KGnuEc5v1YFV6Q4+3y/p8k941r9I+pcjnH9O0nOfu9gkN710un770W/VF+ob1tdXJpbn6slgtZQhqWl9wgRPFTkV3nT8GVZNqirIdKEQIL4Gw4N6YvMT+s1Hv9Henr2aVjJNvzr2Vzqt9jQZRkpiv/0r0G1aKm36o7Q3Oh9fTrk0+TxGNQEAACSouMzxBG9NL5musA1rXcs6HVNxzLA9t648VxttjfOhca3z2+cE8LOTf+ZNx5fdFXOTX1z5sewWSBj9oX49uulR3b36bjX2NmpW2Sz9f8f/fzpp1EkETnD0dzqjmTZGw6be5uiopgXSGX8vTWRUEwAAQKIjeEoBc8rnSJLebXx3eIOnijx1KFe96aXKZoJxAFG9wV49vPFhLV69WC39LTqm/Bj984n/rOOrjidwgtSyxQmaNr4g7fiTFAlJmYXOLy8mLWIFOgAAgCRD8JQCijOLNbFwolbsXaFvzvzm8D03J13FOenOynZN64btuW67Y8UdkqRbFtwS346fv9XZn3v0I67+8Wlnbq4fXzh8KxICbukJ9mjJ+iW6d829ahto03FVx+nfZ/27jq081uvS4KVwUNr5Zydo2rhUatnknC+bIi38jhM21SyQfPyRBAAAIBnxp7wUsaBygR7f/LiC4eCwrho1sTxXG9prNLHpRSkSSYjXIda3ejQ6a99HMTdZu6fThUKA4dU52Kn7192v36/7vToGOnTiqBN106ybhkZbIgX1tjqTgm98Qdr8sjTQIfnSpbEnSQuul+r+Sioe53WVAAAAiAOCpxSxoHKB7l9/vz5q/mjY53latbdC56tX6tgpFY0dtmcDGNk6Bjr0u7W/0/3r7ldXsEun1ZymG2bdoJllM70uDfFmrdS47sCopvoVko04E4NPu9AZ1TT+NCkjz+tKAQAAYhaJWIWtVcRaRSJy9tYqYiVZ57NVdG8le9jnA+c/fm9+pl/l+cm9qBTBU4qYXzlfRkYr9q0Y9uDp6YEqZ2W7xvUET0AKaO1v1X1r7tMD6x9Qb6hXZ40+SzfMukFTS6Z6XRriKdgvbX/zQNjUsdM5XzVbOuX/SJPOkarmJsRIWAAA8PlZaxWKWA2EIhoIhjUYjmggGNFAKKLBUEQDoXB0H9FgOKJQ2CoUiSgYtgqFIwqGo8dD56xzLhK9NxxRMLL/Xufawc/Y/zkYiSgcsUObtXJCo8PCo3DEDgVIh9974Jpzn5u+vGC0fnppcv/iluApRRRkFGhS0SS9s+8d3TT7pmF7bl1FnjbtX9muaZ00edGwPRvAyNLc16x7Vt+jhzY+pP5Qv84Ze46un3W9JhVN8ro0xEvXvujE4Eulra9KwV7JnyVNOF065W+dV+jyq72uEgCAlBeJBkB9wbCzDYbVH3S2/Z/79n8eDKsvGBm63h8MO+HRQYHR/s9HCpGGQqZQRHaYMxpjpIAvTYE0I78vTQGfUcCXJr/PKJDm7P1paQr4999jlBvwy5dm5DNGaUN7Kc0Y+dKM0oyJHku+NCNjnHsOXHPOD7U1Oug5h94jOc815sDeGCNzyHnJ6NB7Dr53THH28P5DG4EInlLIsZXH6qEND2kgPKAMX8awPLOuPFedylFvRlnCrGw3Jn+MNx2XTIi5yfiyHBcKAWLT0NOgxWsW65GNjygYCeq8cefp+pnXa3zheK9Lg9v2v0K34Tln273KOZ9fI835ivMK3diTpECWt3UCAJCAQuGIeoNh9Q6E1TMYUu9AWL2DIfUOHvjcs//zgLPff70/GhQdHCrtD5mc8CjymWrKDKQpw+9Thj9NGYE0pfuin6PHBVkB51wgzbnHf+D+9Ojn9I+d8w1d238u4Ev7WIAUiAZLfl+a/GnO5/3hDhIbwVMKWVC5QL9f93t92PThsK0yVZaXobxMv+rTx2tSw+pheabbbj/hdm86vuhXMTf56aWzXCgEODp7u/fqt6t/q8c2PSZrrS6YcIGun3m9RueP9ro0uCkcknYulzY8L214Vmrb7pyvPkY64++lSedKFdOdX0ECAJBCQuGIegbC6hoIqnsgpO7+kLqi+4M/9w2G1DMYVu9AdD8YUs/AgX1f0AmSBkJHHw7504yy033KTvcrO8OnrIBPmQFnX5Sdrqx0n7ICaUPnMgO+6DlnywikOcfpB65/7B5/mtIIeuACgqcUMq9yntJMmt7Z986wBU/GGNWV52pN7zhNanxYCvbxm28gwdV31es3H/1GT255UpJ0ycRLdN2M61STV+NxZXDNQJez+tyG56VNS6W+NmcVunGnSid+zwmb8qu8rhIAgM8kHLHq6g+qsy+kroGguvoPhEUHgqPgx4Okw8KlvmD4qPpzAqJoSJTuU06GX3mZflXmZyo7w6ecaHiUE72ene5XTkZ0n+5TdsZh+3S/0v3MmYjERfCUQvLT8zWleIpW7Fuhb+vbw/bcuvI8vbW2Rl+wYalhjVQzf9ie7Ybb37rd2cd75NNTNzv7GEY+3fbYh5IY+YT42N29W3d9eJee3Pyk0kyaLp90ub4x4xuqzKn0ujS4oXOvtPF5af1z0rbXpPCglFUk1Z0jTTlPmnAGq9ABAEaEcMSquz+kzv6gOvqC6oyGSJ1Dx0F19h/8OXTI+e6B0F/sw5dmlJvhV240JMrN8Ks4J12ji7OHPudmBJSb6Vdehl+5+88d9jkn3c+oIeAwBE8pZkHlAv1h3R/UH+pXpn94lmysq8jV4pW1UqakPe+N+OBpR+cObzpu2RJzk61NPS4UAhxqT/ce3fXRXXpi0xMyxujyyZfruhnXqSKnwuvSMJyslRrXOnM1rX9O2vOuc75onLTgBmnyuVLt8ZKPPxoAAIZfJGLV1R9Se9+g2nuDau8Lqr13UB19Qedzb1DtfYOHhEZd0TCp6y8ER8ZIeRl+5WcFlJ8ZUH6WX6OLsw/57OwDyjticBRQZiBNhtfIAVfwp8sUc2zlsbpnzT16v+l9HV91/LA8c2J5rnarVMGMIgX2vj8szwTgvr3de3XXR3fp8c2Py8joi5O+qOtmXscIp2QSCTvzNa1/1tnao8H7qPnSmf8gTT5PKpvCfE0AgKMWCkfU2R9Se++g2vuC6ogGRvvDo45ooNQWDZc69t/XF/zUFc/2B0cFWZ8eHOVnHhow5WcFlMsoI2BEI3hKMceUHyOf8WnF3hXDFjzVVeRJMmrKm6rqPR8MyzMBuGdfzz795qPf6NFNj0qSLqu7TN+c+U0Cp2QR7Je2LpPWP+3M2dTbIvkypPGnSSf/wFmJLo9/1wAAZxRSZ39QrT2DausdVEu3s2/tCaq1Z0CtPUHnfM+g2qL3dPV/+uij/Ey/CrPTVZQdUEF2usYUZ6swO6DCLOdzYVbA+ZwdUEFWenQfUMDHHEZAsiJ4SjG56bmaVjJN7+x7Z9ieWV2QqZx0n7b4Jqq66X7nLz2B4XmND8Dw2R84PbbpMVlZXTrxUn1z5jdVlcuk0Qmvv0Pa9KK07mlp80vSYLeUkS9NOkeacoE08SwpI9frKgEALusPhodCotaDtoPDo4PPtfUGFY4ceRhShj9NJTnpKs5NV1F2usaWZKsoO30oRCrMTlfBQceFWc6IJB8jjwAchuApBR1beazuW3OfeoO9yg5kf+7nGWM0sTxXK4OjdXIkFJ1gfN4wVOqOKcVTvOm4cmbMTaZV57tQCFJNQ0+Dfrv6t3pk4yOy1uqSukt0/czrVZ1b7XVp+Dy6G6Ov0D0jbX1NigSl3App5uXS1AuksadI/nSvqwQAfA7hiFV776CauwfV0j2g5p5BNXcNqKVnQM1dg86+e1DN3QNq6R78xFXX0oxUlJ2uopx0Feeka0JZropy0lWSs/9cQMU5GSrOTldRTkAlORnKSvfF+dsCSFYETyloQeUC3b36br3f+L5OGHXCsDxzalW+Xlxdpb+RpL3vjejg6ZYFt3jT8bk/i7nJjy+c7kIhSBWNvY367UdO4BSxEV088WJdP+t6jcod5XVp+KxatzlB07pnpF1vS7LO5ODH3yRNuVCqOVZK41UFABjJ+oPhoaBoaH9IkLT/2qBaewZ0pAFJvjSjkpx0leRmqDQ3XeNKc1Sck66S3PRoeHRQqJSdroKsAHMgAfAMwVMKmls+V37j14p9K4YteJo+qkBL3ilUpLBIaXveH5ZnAvhsmnqbdPfqu/XwxocVioScwGnm9arJq/G6NMTKWqlhtRM0rX/GOZacEZSn3eaMbCqfxuTgAOCxUDiilp5BNXUNqKlrQI1d/QcdO/vmbmd0UvcnrNCWk+4bCpJqi7M1d3SRSnOdAKk0L0MlOc610twMgiQACYXgKQVlB7I1o3TGsM7zNL06X5JRa8E0lY7wle1ufeNWSdLPTo59BNLn8uj1zv6yu466yfeXvCdJ+sWVc92oCEmmua9Zv/3ot0OB00UTLtL1s65XbV6t16UhFpGwtGuFEzStf0Zq2y7JSKMXSuf8qzTlfKlorMdFAkDys9aqZzCsxs7+QwKkpu4BNXY6eydc6ldLz+ARV2wryAqoLC9DZbkZmllTOBQc7d+X5GY4wVIur7YBSF4ETynq2Mpjdffqu9UT7FFOIOdzP29qZb7SjLTVP1Gle/8woicYb+hp8Kbjzj0xN9nb0e9CIUg27f3tunvN3Xpg3QMKRoK6YPwFunHWjarNJ3BKGKEBadvrzuTgG56TepokX7qzEt1JP5AmnyfllnldJQAkBWutOvtC2tfZr4bOfu3r7FdjZ/+BYOmgkOlIcyYFfEZluRkqy8vQqMIszaktVFlehsrznHP7j0tzM5QZIEwCAIKnFLWgaoHu+ugurWpYpVNqTvncz8tK92lCWa5WBcdoQSQkNa6RRo3ceZ6AZNA12KX71t6n3639nXqDvTp//Pm6afZNGpM/xuvScDQGe6UtL0trn5Q2LpUGOqX0XKnur5xX6CaeLWWywAAAxKI/GFZDZ78aOgeGAqV9Hf1q6BpQQ0e/GrqcsKk/GPlY24NHJ80dXTgULpXnZ6gsNzO65zU3AIgVwVOKml02W4G0gN7Z986wBE+S87rdH7dU6VuStOd9gifAJb3BXt2//n4tXr1YnYOdOnvM2fr27G9rYtFEr0vDXzLQLW36oxM2bfqjFOyVsoqlaRdLUy+Sxp8q+TO8rhIARpxwxKq5e8AZoXRwkDQ0YskJmjr6gh9rmxlIU2V+psrzMzW7plAV+RmqyM9URX6mKgsyVZHnhEqMTgIAdxA8pagsf5Zmls7Uin0rhu2Z06sL9MT7+YoUFipthM/zBCSigfCAHtrwkH7z0W/U2t+qU2pO0XfmfEfTSqZ5XRo+TX+nM6Jp7RPS5pekUL+UUy7N/rITOI05UfLx4xhA6gqFI2rqHtCedidU2tvRp70H7fd1OK/BhQ9b3s2X5rzyVlGQqTEl2TpufPFQoFSRnzEUNuVn+mVYhAEAPMOfdFPYgqoFuvPDO9U52Kn89M//Osf0Uc4E4x2F///27ju+6vLu//jryt7Jyd4JBAibEAKIiIKi4kBcdbZ1tI7aWu3W1o7bX+/7tq21Wqt1Vax3Hdi6cIsDtdYFgbA3JCGDANl7Xb8/vocQFDSBnJyM9/Px4PEN3/Mdn4PHk5x3rutzTcA1gFe2mxI3xTs3Tpve61NyM1weKEQGm7aONp7f9jwPrnmQisYKZibN5Hs53yMnPsfbpcmRNFXB5tdgw1JnOl1HK4QnQe4VTtiUfhz46DfrIjL0HSlUKq9pprSmyRm9VNvM5zIlQgJ8SYoMIikymNmjYkmKDDo4SskdLMWEBeKrKW8iIgOegqdhbEbiDB4oeICV5SuZlz7vmK83ISkSgJ3+o3CVPOk0yx2AU0Zunnazd248/ze9PuVnC8b2fR0yaLR3tvPyjpd5oOABSupLyInL4X9P+F9mJM3wdmlyOA37YfMrzjS6Hcuhsx0i02DGtc40utTp4OPj7SpFRPpMZ6dlb30LJdVNlFUfPlQ63EilYH9fkqKCSHaHSsmRQSRGBpMUFdQVNmmUkojI0KHgaRibEjeFEL8QPij5oE+Cp8gQf1JdweS3Z5Lb2QZ71kNKbh9UKjK8dNpO3tz1Jvetvo9dtbsYHzOe2467jdnJs/VD+EBTXwGbXnbCpp0fgO0AVybM+q4zsik5F/TfTEQGqZb2DsqqmympbqKkqsnZur8urXHCptaOQ5t0HwiVkiKDukYqJUUGO9uoIJIigokIVqgkIjKcKHgaxgJ8A5idMpv3it+j87hOfMyx/yZ+YnIk75Ql822AstUDMnj6wbs/AOBP8/7Uvzde8nVne/E/enzK9f+3EoAHvqFG7cOBtZZ3i9/lL6v/wtaqrYyKGsXdc+/m5PST9QP6QFJbBhtfcsKmov+A7YSYUXDCzU7YlDhZYZOIDAo1TW1dgVJpt1DpQMC0t67lkOONgYTwIFJcwUxOjeKMicGkuIJJiQoiOSpYoZKIiByWgqdhbl7aPJYVLmPj/o1MiJ1wzNebkBzBH9eHYaOiMGUFfVBh36tuqfbOjRuren1KVWOrBwqRgcZay39K/8O9q+5l/f71ZERk8Ls5v+P0zNPxVR+ggaFmt9OvacOLUPwJYCFuHJz4Uydsih+nsElEBhRrnWlwu6ua2F3lHqXUfcRSdRN1Le2HnBPg50NKVDApUcGcnB1PctSBYCmYVFcwCRFBBPhpyrCIiPSOgqdhbk7KHHyMD+8Wv9s3wZO7wXitawKRA7jBuMhAsaJ8Bfeuupf8inySQ5O5/fjbWZi1ED8fvT17Xc1uWP8CrH8eSlY4+xImwbxfwPhzIC7bq+WJyPBmraW6sY3dVU0UVzWyu6qR4soDXzexu6qR5rZDp8FFBvuTEhVMekwIs7JinJDJFewETFHBxIYFaLSSiIj0OX2yGeaigqKYGj+Vd4vf5XtTv3fM15uY7DQYLwoczaTiJ6G9FfwCjvm6IkPNmr1r+Muqv/BR2UfEB8dz28zbOH/0+fj7+nu7tOGtttQZ1bT+effIJiBpCpzya2dkU0yWd+sTkWGlvqWd4spGiisbuwKm4konVNpd1UT950YsRQT5kRYdwqi4MOaOiSMtOoS06GBSokJIcQUTFqgf/UVEpP/pu48wL20ed664k5L6ElLCUo7pWvERQcSGBbK6PZNJHa1QsQGSc/qmUJEhYFPlJu5bdR/Ldy8nOiian+T9hIuyLyLIL8jbpQ1fdXsOhk1FHwHWGdl0yq9g/LkKm0TEY5rbOpxAqdI9Yqnraydkqm5sO+T4kABf0lwhpLqCOW5kDKmuYFJdTriU6gohMli/vBARkYFHwZMwN20ud664k+XFy7l83OXHfL0JyRG8U5XMN8BpMD7AgqeZSTO9c+ORJ/X6lNmjYj1QiHjDjpod3LfqPt4sfJPwgHBuyr2Jy8ZeRoh/iLdLG57q98LGF52pdLv+DViIHw/zfg4TzoPY0d6uUESGiOrGVgr3N1JY2UjR/gaKKhsp3N9IUWUj5bXNWHvw2AA/n64waXJqJGnRTsiU5gohLToEV4i/psKJiMigo+BJyIjIYETkiD4LniamRPDgtjBsRCSmdDUMsAXZrp9yvXdufNJPe33K90/Rh9/BrrS+lAcKHuDF7S8S5BvEtZOv5YoJVxAREOHt0oafhv2w6SVY9xzs+sBZjS52DJz0Mydsih/r7QpFZBDq6LSU1TRRVNlIUVfA1EhhZQNF+xupbT50OlxceCAZ0SHMGhlDekwIGTEhpEeHkOoKIS4sEB8fBUsiIjK0KHgSwJlu9/j6x6lrrSM8IPyYrjUhOZL2TmiInkBY2eq+KVBkkNnftJ+H1z7MM5ufwWC4fNzlfHvSt4kOivZ2acNLYyVsesWZRrdjOdgOiM6COT9yh03jtRqdiHyl5raOrpFKhfsbKK48GDDtrmqiteNgE28/H0OqK5j0mFCmprm6gqV09zYkQD9+i4jI8KLvfAI4wdOj6x7lw5IPWTBiwTFda0KyM5KjOCibcYVPDLgG49e/5Yx4emD+A/17439c4Gy//myPT7ni0U8B+PvVMzxRkXhAbWstj617jH9s/AetHa2cO+pcrp9yPYmhid4ubfhoqobNrzph0/Z3obMNXJkw+yYnbEqcpLBJRL6gsbWdwv2N7NrXwM79DRTua3S2+xvYU9tyyLHhgX6kx4QwNimc0yYkkh59cORSclQwvhq1JCIi0kXBkwAwKXYS0UHRvFP8zjEHT2muEMID/SjoyGRcRyvs3eisCjVAtLS3fPVBntDW3OtTmts6PFCIeEJTexNPbnySR9c9Sm1rLWdknsENOTeQGZnp7dKGh+Za2PyaO2x6GzpaITIdjvsOTDwfknIUNokITa0dFFY2OOHSPmf00s59Dew6TLgUGxbIiNgQ5oyOI8M9YikjJpSM6BCi1GtJRESkxxQ8CQC+Pr6cmHoibxe+TVtnG/4+R78qio+PYXxyBO/VpXAJQOnqARU8ifSlto42nt36LA+ueZB9Tfs4MfVEbpx6I2Oj1S/I41obYctrTs+mrcugowUiUmDGtTDhfEjJVdgkMgw1t3VQuL+xK1DqCpf2Oc28u4sNCyAzJpQ5o+PIjAkhMzaUzJhQMmJCCA/SCnEiIiJ9QcGTdJmbNpcXtr1A/p78Y175bUJyJE9+WoUNicCUrQau6JMaRQaKjs4OXt35Kvetvo+S+hKmJUzjrrl3MTV+qrdLG9raW50RTWv/5YxwamuA8CTIu9qZRpc6HXx8vF2liHjYgZ5LTqDUwC73FLld+xsoqzk0XIoJDSAzNpTZo2K7wqURsQqXRERE+ouCJ+kyK2kWgb6BLC9e3gfBUwTNbZammImElK7qmwJFBgBrLe8UvcO9q+5le812xkWP45fzf8nxycdr2oWndHbArn/Dun/BhqXQXA3BLpj8NZh4IWQcDz6+3q5SRPpYW0cnxe5wacdep+/SLnfQVFbbjLUHj40JDSAjJoRZWTGMiAklIzbUvQ0hQuGSiIiIVyl4ki4h/iEcl3Qc7xa/y0+n//SYPkRPTIkEoDh0ItnbH4WWeggM66tSj8lJqSd558ZjTu/1KaeMi/dAIXK0Pir9iD/n/5l1+9eRGZHJH0/6I/Mz5uNjNMKmz1kLJSudkU3rn4f6cvAPhbFnwaQLYeS8AbVogYgcHWstlQ2t7NjXwI699ezY28D2vQ3s2FdP0f5G2jsPpkvR7nDpuJExZLpHLDkjl0KJDFa4JCIiMlApeJJDzE2by3u732Nb9TZGu0Yf9XWy4kIJ9PMh34wju7Mddn8KWSf3YaVH78qJV3rnxrO/3+tTrj0xywOFSG8V7C3g3vx7+aT8E5JCk7j9+NtZmLUQPx+9hfa5PRuckU3rnoWqXeAbAKNPc8Km0adDQIi3KxSRo3Cg79KOvfXs2NfAdnfItGNvPbXN7V3HBfj6kBkbwpj4cBZMSGRkXBgj40LJig0jMkThkoiIyGCkT01yiAOjgZYXLz+m4MnP14exieG8WRvIpcYXCv8zYIInkZ7aUrWFe1fdy/Li5UQHRXPLjFv42pivEeCrkTZ9qnKnEzStexYqNoDxgZFz4cSfwrizISjS2xWKSA9Yaymvbe4KlLbvdZp679hXz+6qpkOmxiVGBDEyLpSFU5IPCZdSXMH4+mjasoiIyFCi4EkOERcSx6TYSbxb/C7XTL7mmK41ISWSlwtKsclTMIX/6aMKj91Vr18FwOIFi/v3xovPchfwSo9PufjBjwBYct0sT1QkR1BcW8x9Bffx6o5XCfMP4/tTv8/l4y4nxF+jbfpMXbkzhW7tv6BkhbMv7Tg4804Yfy6ExXm1PBE5stb2Tnbtb2BbRX3Xn+1769m5r4HG1o6u40ICfBkRG0pOmovzp6Y64VJcGCNiQwkN1I+gIiIiw4W+68sXzE2by72r7mVv417iQo7+w9+E5Aie/KSIuoQZRKxZDG3N4B/Uh5WK9K09DXt4aM1DPLf1Ofx8/Lh64tVcNfEqIgM14qZPNFbCxqVO2LTr34CFxEkw/79g4vkQle7tCkWkm/qWdrYfCJf2ugOminoKKxvp6NZ7KSUqmKz4MKZnRpMVF9o1gikxIkiLLoiIiIiCJ/miA8HT8t3L+dqYrx31dSYkOx/WtwZPZlpHi9MoOHN2X5Up0mdqWmr429q/8eSmJ+mwHVw45kKunXztMQWv4tZSD5tfc/o2bXsbOtsgOgtO+qmzIl3cGG9XKDKsWWvZV9/aFS5td49e2lZRT1lNc9dxfj6GzNhQxiSEc+akJEbFhzEqXqOXRERE5KvpJwX5gtFRo0kNS+XtorePKXgamxiOr4/ho/YxTMM4fZ4UPMkA0tzezBMbn+Bv6/5GfWs9Z488mxtybiA1PNXbpQ1u7S2w7S1nZNOW16GtEcKTYeZ1TpPwpBzQKAiRftXZaSmpbjpketyBUUw1TW1dx4UE+JIVF8ZxI2MYFR9GVpwTMGXEhODvqxU8RUREpPcUPMkXGGM4NeNU/m/D/1HTUnPU04yC/H0ZFRdGfgWQMAEK/w38pE9rFTka7Z3tvLjtRe4vuJ+KxgpOTD2R70/9PtnR2d4ubfDq7ICd7zsjmza+BM01EBwNUy5xRjalzwIffWgV8bQDAdOWPXVs3lPH1j31bNlTx/a99TS3dXYdFxMaQFZ8GGdNTmKUO1zKig8jKSIIHzX3FhERkT6k4EkO69SMU1m8fjHv7X6Pc7LOOerrTEyJZPnmCuy04zGr/gEdbeDr3eWQT8883Ts3nnBur085e3JS39cxjFlreafoHe5ZdQ87a3YyOW4yv5vzO/IS87xd2uBkLZTkw9pnYN1z0FABAWEw9mxnZNPIuV7//11kqLLWUlbTzBZ3uOSETHVsrag/pMF3UmQQoxPCu0YwjYoPY1RcGK5Qrc4pIiIi/UPBkxzWxNiJJIYmsmzXsmMKnqZnung2fzcV0XkktD0EZQWQ6t0P+ZeMvcQ7N57R+1UCvzErs+/rGKY+K/+Mu/PvZs3eNYyIHMHd8+7m5LST1fj2aOzfDmuegbX/hMrt4BsAo0+DSV+DMaeDf7C3KxQZMqy17K1vYUu5M3Jpa0Udm8udsKmupb3ruNiwQLITw7goL43sxHDGJIQxKj6cyGCFvyIiIuJdCp7ksIwxzE+fzzObn6G+tZ6wgLCjuk5epguAjzvGsAig8EOvB09N7U0ABPv184fj1kZnGxDS41Oa3L+1Dg7w9URFw8Lmys3ck38PH5R8QHxIPP91/H9xTtY5+Pno7a9X6vbA+uecwKk0HzCQeQKccDOMOweCo7xcoMjgV9nQypY9dd3+OGFTdePBHkyuEH9GJ4Rz7tQUxiSEMSYhnDEJ4RrBJCIiIgOWPnnJEZ2acSr/2PgP3t/9PmeOPPOorpEVF4YrxJ8Py3xYFDsGdn0Is2/q40p754a3bgBg8YLF/XvjJ9yN2q96pcenXLn4UwCWXDfLExUNaSX1Jdy36j5e3vEyYQFh/GDaD7hs7GUE+QV5u7TBo7kWNr3shE073wPbCYmT4bTfwoTzITLF2xWKDEo1TW1s7RYsHfizr76165jwID/GJIRzxsTErnBpTEI4sWEBGqkpIiIig4qCJzminPgc4oLjeKvoraMOnowxTMuIZsWuKsg+HtY97zQh9tEIHvGMyuZKHl7zMEs2L8HH+HDlxCv51sRvHXWT/GGnvdW9It0zsPk1aG+GqAw44Ycw+SKIUwN2kZ5qbutgW0U9m8rr2Fxeyyb3FLny2uauY0ICfBmdEM687HgnXHJPk0uMCFLAJCIiIkOCgic5Ih/jw8npJ7N0+1Ia2xoJ8e/5FLHupme6eGvjHupOnEH4ysdgz3pImty3xcqw19jWyOMbHuex9Y/R1N7EeaPO4/op15MYmujt0ga+zk4o+sgJm9a/AM3VEBIDU7/h9G1KmwH6ACxyRNZadlc1dQVMG8udPkw79zXQ0WkBCPDzYVRcGMdnxTA6IbxrmlxKVLBWkRMREZEhTcGTfKnTMk5jyeYlfFj6IadmnHpU18jLjAZgJeOZC06fJwVP0kfaOtt4dsuzPFDwAPub93NK+il8f+r3GRk10tulDXx71jvT6NY9CzXF4B8CY8+CSRdB1jytSCdyGDVNbWwur2OTewTTprJatuypp75bo++06GCyEyI4Y2Ii2YnhjE0MJzMmFD9fHy9WLiIiIuIdCp7kS+Um5OIKdLGscNlRB08TUyII9PPh3xWBzI3KcIKn477Tx5XKcNNpO3lz15v8edWfKa4rZlrCNO6edzc58TneLm1gqy6Gdf+CNf+EivVgfGHUKXDKryD7TAg8uoUERIaa1vZOduyrZ3N5HRvLDk6VK6s5OE0uIsiPsUkRnJ+b4g6YIshODCcsUD9eiYiIiBzgsZ+MjDFpwONAAmCBh6y19xhjooElQCawC7jIWltlnEYG9wBnAo3AldbafPe1rgBuc1/6t9bav7v3TwMeA4KBV4GbrLXWU89pOPLz8ePk9JN5bedrtHS0EOgb2OtrBPr5MiU1is8KqyBjNmx9A6z12tSdRaMWeeW+5FzW61MunJbqgUIGv49KP+JPK//ExsqNjHaN5r5T7mNOyhz1QzmSpipY/7wTNhX9x9mXOgPOvBMmnAehsd6tT8SLrLXsrWthfVktm7oFTNv31tPW4fxI4e9ryIoLY+aIaLITIxibGM7YpHD1YRIRERHpAU/+Sq4d+JG1Nt8YEw6sNMYsA64E3rbW3mGMuQW4BfgZcAYw2v1nJvBXYKY7qPo1kIcTYK00xiy11la5j7kG+AQneFoAvObB5zQsnZpxKs9ufZaPSj9ibtrco7pGXqaLh97fQeuM4wgoeBL2bob4sX1baA+dO+pcr9yXqZf3+pSv5aV5oJDBa/3+9dy98m4+LvuY5NBk/ueE/+HMEWfiq2b1X9TeCtuWQcFTsOUN6GiF2DEw7zaYdCFEj/B2hSL9rr2jkx37GthQWsvGslo2lNWyobSW/Q0HV5NLjgwiOzGcudnxjEsKJzsxnJGxYQT4aZqciIiIyNHwWPBkrS0Dytxf1xljNgIpwCJwWv0AfweW4wRPi4DH3SOWPjbGRBljktzHLrPWVgK4w6sFxpjlQIS19mP3/seBc1Hw1OdmJM0gIiCCZYXLjjp4mp4Zzf3Lt7POfxK54Ey381LwVNVcBYAryNW/N27Y72xDY3p8SqX7w1B0aIAnKho0iuuKuTf/Xl7b9RpRgVH8dPpPuTj7YgJ8h/e/yxdYCyX5Tti07lloqoTQOJj+bZh8MSRNUZNwGTbqmtvYVF7HhlInXNpQVsvmPXW0tncCEODrw5jEME4eG8/45AjGJ0UwNjGCyBD1NhMRERHpS/3ShMAYkwlMxRmZlOAOpQDKcabigRNKFXc7bbd735ft332Y/dLH/H38mZs2l3eL36Wtow3/o2g4nJvuwhj4cF8YueFJUPgfmP4tD1T71X64/IcALF6wuH9v/Mw3ne1Vr/T4lO/8YyUAS66b5YmKBrzq5moeXPMgT29+Gn8ff66ZdA1XTbyK8IBwb5c2sFQXwZolUPA07N8GfkFOv6Ypl0LWyeCrfjMydFlrKa1p7gqYDoxkKqps7DrGFeLP+OQIrpiVwfjkCMYlRZAVF4a/mn2LiIiIeJzHP40YY8KAZ4GbrbW13XshWGutMcbjPZmMMdcC1wKkp6d7+nZD0mkZp7F0+1I+Kf+EE1JO6PX5kSH+ZCeE81lRNWQc74x48mKfJxnYWjpaeHLjkzy85mEa2hs4b9R53JBzA/Eh8d4ubeBoroUNLzphU+G/nX0ZJ8Dsm2D8IgiK9G59Ih7Q2t7J1oo6d8BUx4ayGjaU1lLb7KwoZwxkxoQyKSWSi6enMS4pnPFJkSREBKoXk4iIiIiXeDR4Msb444ROT1hrn3Pv3mOMSbLWlrmn0lW495cA3RvapLr3lXBwat6B/cvd+1MPc/wXWGsfAh4CyMvLU/PxozAreRah/qG8VfjWUQVP4PR5emFVKZ1nzcZn3bNQtROiteS9HNRpO3llxyvcu+peyhrKmJMyhx9M+wGjXaO9XdrA0NEOO951ptJtegXamyFmlNO3afJF4MrwdoUifaaqofVgHyZ3L6buDb+D/X3JTgzn7CnJjE9yRjGNTQwnVCvKiYiIiAwonlzVzgB/AzZaa+/q9tBS4ArgDvf2xW77v2eMeRqnuXiNO5x6A/gfY8yBhjynAbdaayuNMbXGmONwpvB9E7jXU89nuAvwDeDE1BN5p+gdbjvuNvx8ev/SycuI5h8fF7EjdDKjwJlup+BJ3D4u+5i7VtzFxsqNjIsex/+b/f+YmTTT22V5n7VQvgYKlsDaf0JDBQS7YOrXnal0KdM0clAGNWstZTXNrCupYV1pLetLathYVktpTXPXMfHhgYxPjuDksfGMS4pgfHIEmTGh+ProtS8iIiIy0Hny14KzgW8Aa40xq937fo4TOD1jjPkWUAhc5H7sVeBMYBvQCFwF4A6Y/h/wmfu42w80GgduAB4DgnGaiquxuAedlnEar+18jZV7Vh5VIJCX6WSHH1bHMiokxgmepn69r8uUQWZr1VbuWnkX/y75N0mhSfzvnP/lzBFn4mOGee+V2lInaCp4Gio2gI8/ZC+AyZfA6NPAT43VZfCx1lJU2ci6klrWldawrqSG9aW1XQsp+BjIigtjxojoroBpXFIEsWGBXq5cRERERI6WJ1e1+zdwpF9FnnKY4y3w3SNc61Hg0cPsXwFMPIYypRdmp8wm2C+YN3a9cVTBU0pUMEmRQXxWWMUVGcfDzg+80ufp4uyL+/V+XaZf3etTvn7c0J06tadhD/cX3M8L214g1C+UH077IZeNu4xA32H8AbO1ATa+7Eyl27EcsJA6A876I0w4H0KivV2hSI91dFp27mtgvTtgOhA21bn7Mfn7GsYkhHPquAQmpkQwISWScYkRBAf4erlyEREREelLaoQgPRbsF8z89Pm8tvM1fpz3Y0L8Q3p1vjGGvMxoPttZiT39VMzGl6BkJaTmeajiw1swYkG/3q/LxAt6fcrCKckeKMS7GtoaeHTdozy+/nHabTuXj7ucayddS1RQlLdL847ODtj5vrMq3Yal0NYAUelw4k9gyiUQk+XtCkW+UltHJ9sq6rtGMK0rqWFDWS2NrR0ABPj5MC4pgnOmJDMpJZKJKZGMTggj0E8hk4iIiMhQp+BJeuWCMRfw0o6XWFa4jEWjFvX6/OmZLl4qKKU05XRS/IJg9ZP9HjyVN5QDkBia2K/3pWa3s41M/fLjuimtbgIgOSrYExX1q7bONp7d8ix/Lfgrlc2VnJF5Bjfm3khaeNpXnzwUVWx0ptGteQbqSiEwAiZd4PRtSjsOfIb5VEMZsFraO9hSXt81VW5dSQ0by+tobe8EICTAlwnJEVyUl8bElEgmpkSQFReGv69e0yIiIiLDkYIn6ZXc+FwyIzJ5duuzRxU85WU4U4U+K+sgZezZsO5ZWPC/4Nd/06tu/eBWABYvWNxv9wTgueuc7VWv9PiUHyxZDcCS62Z5oKD+Ya3lneJ3uHvl3eyq3cW0hGn85eS/MClukrdL63/1e2Hdv5ypdGUFYHxh1Hw4/b8h+wzwH/wBowwtTa0dbCirPWS63JY9dbR3OivLRQT5MTElkiuPz2RCcgQTUyLV9FtEREREDqHgSXrFGMMFoy/gjyv/yI7qHYyM6t2qdNmJ4YQH+vHZrkrOzbnM+RC++TWYcK5nChavKthbwF0r7iK/Ip8RkSP487w/MzdtLmY4rcLW1gSbXnGm0m17G2wHJOXAgjtg4oUQFuftCkWAAyFTDWt217B2dw1rS2rYvrced8ZEdGgAE1MimZsd54xkSo4kLTp4eP3/LCIiIiK9puBJem1h1kLuWXUPz259lp9M/0mvzvX1MUzNcLFiVxUsmgvhyc7oDwVPQ0pxbTF359/Nm4VvEhMUwy+P+yXnjz4fP59h8pbT2QmFH8Kap2H9i9BaBxEpcPyNTt+m+HHerlCGuea2DjaV17F2d7UTNJXUsGVPXVfIFB8eyOTUSM6clNQ1XS4xIkghk4iIiIj02jD5FCh9KSY4hnlp83hp+0vclHsTAb69W9Z9eoaLPy7bQk1zJ5GTL4L/3Av1FRAW76GKpb9UNVfx4JoHWbJ5Cf4+/lw/5XqunHAlof6h3i6tf+zd7PRtWvtPqCmGgDAYvwgmXwyZJ4CPGilL/2vr6GRzeR1rS9yjmUqq2VxeR1uHkzLFhAYwKTWS08YnMDk1ikmpkSREBHm5ahEREREZKhQ8yVG5YPQFLCtcxjvF77Ags3erxOVlOn2eVhZVcnLOZfDh3U6D5eO/54FKpT80tzfzxMYneGTtIzS2N3LeqPO4IecG4kOGQZhYv9fpVVbwFJStBuMDWSfD/N9A9pkQ0LvVH0WORXtHJ9v21ndNl1tTUsPGstquxt+Rwf5MTo3kmjkjmZwayaTUKJIjNZJJRERERDxHwZMclVnJs0gOTea5Lc/1OnjKSYvCz8ewYlcVJ48dC8m5zof2fgqerphwRb/c5wuO4vldM6d3PbT6W6ft5OUdL3PvqnspbyjnpNSTuDn3Zka5Rnm7NM9qa4LNr0LBEtj2ltO3KXEynP4/Tt+m8ARvVyjDQGenZce+BtaWVHcFTetLa2lq6wAgLNCPiSkRXHl8JpNSIpmSGqWeTCIiIiLS7xQ8yVHxMT6cO/pc7l99P7vrdpMantrjc4MDfJmYEun0eQLIuQxe/TGUrYGkyR6q+KC5aXM9fo/Dyj6j16fMHz9wA4yPSj/irpV3salyE+NjxvPfs/+bGUkzvF2W53R2QtF/nJB0w1JoqVXfJuk31lpKqpsoKK6hYHc1BcXVrC+tpb6lHYBgf18mpkRw6Yx090imSEbEhOKj1eVERERExMsUPMlRO2/UeTxQ8ADPb3ueG6fe2KtzZ46IZvGHu6hpaiNy4gXw+q3OB/p+CJ521uwEYETkCI/f6xD7tjrb2NE9PmX73noAsuLCPFHRUdlcuZk/rfwTH5Z+SHJoMnfMuYMzRpyBj/HxdmmesXeL0yR8zTMH+zaNOwemXAyZc9S3STyiurGVgt01FBQ7IVPB7mr21bcCEODrw7jkCM7PTWFSSiSTU6PIigvFz3eI/j8oIiIiIoOagic5aomhicxOns0LW1/gO1O+06sVy86YlMSD7+/gjXXlXDQ9DbIXOB/sT70dfP09WDXc/tHtACxesNij9/mCl252tle90uNTfv7cWgCWXDfLAwX1Tll9Gfetvo+l25cSFhDGj6b9iEvHXUqgb6C3S+t7DfsO9m0qXXWwb9Mpv4axZ0LAMGmWLv2iua2D9aU1h4xm2rW/EQBjYFRcGCeNiScnLZIpaVGMTYwgwE8hk4iIiIgMDgqe5JhcMOYCbn73Zj4s+ZCT0k7q8XlTUiPJiAlhaUGpEzxNuQw2vuT0yzmKKWniOTUtNTyy9hGe3PgkFss3xn+DaydfS2RgpLdL61ttTbD5NVjj7tvU2Q6Jk+C0/4ZJF0J4orcrlCGgo9OyraKeguJqVrtDps3ldbR3OivMJUUGMSU1iounpzMlLZJJKZGEB3k2jBcRERER8SQFT3JMTkw9kZigGP619V+9Cp6MMSyaksxf3t1GRV0z8aNPhZBYWP2kgqcB4sBKdX9b9zfqW+tZmLWQ7+Z8l+SwZG+X1nc6O6HoI3ffphedvk3hyTDruzD5EkgY7+0KZRCz1lJa09w1XW51cTVrS2pobHWaf4cH+ZGTFsV1J41kSmoUU9KiSIgI8nLVIiIiIiJ9S8GTHBN/H3/OHXUuj61/jIrGCuJD4nt87jk5yfz5nW28sqaMq2aPgMkXwWePQGMlhER7sGr5Mu2d7SzdvpT7Vt9HRWMFc1LmcFPuTWRHZ3u7tL6zbysUHOjbVAT+oTB+kfo2yTGpbW5zAqYiJ2T6fF+m8ckRXJSXxpQ0Z4W5TDX/FhEREZFhQMGTHLPzR5/P39b9jRe2vcC1k6/t8Xmj4sMZnxTBi6tLneBpyqXw8f1Ob50Z13iwYjkcay3Li5dzT/49bK/ZzqTYSdwx5w6mJ073dml9o6tv09NQmu/0bRo5D075JYw9S32bpFfaOzrZsqee1cXVrCqqYlVxNdv31mOt05cpKy6MudnxTEmLYkpqpPoyiYiIiMiwpeBJjll6RDozk2ayZPMSrpxwJQG+AT0+95ycZO54bROF+xvISJoMCROd6XYeDJ56E471qRN/3OtTbjy55yvgHYtVFav408o/sapiFZkRmdw19y7mp8/HmEE+GqOtGba8BgVLYNsy9W2So1ZR28yq4mpWFVWzuriKNbsPTpmLDg1galoUi6YkMzXdxeS0SCLUl0lEREREBFDwJH3k6olXc92y61i6fSkXjrmwx+ctnOIETy8VlPK9k0dDzmXwxs+hfK0TEHjArGQvrRCXNa/Xp5wwOtYDhRy0vXo79+Tfw7vF7xIbHMuvZv2Kc0edi7/PIP7QbC0Ufez0bVr/ArTUQHgSHHcDTLkEEiZ4u0IZ4A6sMreqqJpV7qlzJdVNAPj7GsYnOVPmpqZHkZMWRXp0yOAPaUVEREREPMRYa71dQ7/Ky8uzK1as8HYZQ461lktfuZSalhpeOu8l/Hx6nml+7YH/UN3Yxps/OBHTVAV3jYfJX4Nz7vVIrZsqNwEwNnqsR65/RGVrnG3S5B6fsr60BoAJyX27glx5Qzl/LfgrL2x7gRC/EK6eeDWXj7ucEP+QPr1Pv6rc6axIV/AUVO1y+jaNW+iETSNOVN8mOSxrLUWVjU7I5J4yt7GslrYO53tjSlRwV8A0Nd3FhOQIgvz1WhIRERER+TxjzEprbd7n92vEk/QJYwzXTL6Gm9+9mTd2vcFZI8/q8bnn5KTwyxfWsam8jnFJ0U6T8TXPwPz/8kiT8d99+jsAFi9Y3OfX/lKv3+psr3qlx6fc/tIGAJZc1zejtGpaanh03aM8sfEJOmwHl429jGsnX4sryNUn1+93zTXOqKaCp5zV6TBOyDT3Vhh7NgSGebtCGWAONABf5W4Avrq4msoGpwF4SIAvk1Mj+fackUxNiyInPYr4cK0yJyIiIiJyLBQ8SZ+ZlzaPUVGjeGTtI5wx4gx8TM8a6Z41KYn/WrqeF1eXMi4pAmZcC/l/h1X/gNnf93DVw0NLRwtPbXyKh9c+TF1rHWeNPIvv5nyX1PBUb5fWex3tsP0dJ2za/Cq0N0PsGDjl105oGTkIn5N4xJc1AAcYHR/G/HHx5KS5mJoexZiEcHy1ypyIiIiISJ9S8CR9xsf48K1J3+LWD25lefFyTk4/uUfnRYcGcMLoWF4qKOWnp2fjkzgRMmbDZ4/ArO9qitQx6Ojs4KUdL3Hf6vsobyhndspsbs69uf+nGfaF8rXOinRrnoGGCgiOhtxvOlPpknOdpcRkWKtqaGVVcRUrC6vIL6ymYHf1IQ3Ac9QAXERERESk3yl4kj61IHMB9626j4fXPMy8tHk9bri7KCeZHywpIL+oirzMaGdVu39eCVuXQfYCzxY9BFlreX/3+9ydfzfbqrcxIWYCv539W2YmzfR2ab1TtwfW/tMJnPasBR9/GHM6TLkURp8Gfj1fQVGGls5Oy9aKevKL3EFTURU79jYA4OvjNAD/2rRUpqY7o5nUAFxERERExDsUPEmf8vPx4+pJV3P7R7fzUdlHHJ98fI/OO3V8IoF+a3lxdakTPI09G8KT4dMHFTz10uqK1fxp5Z/Ir8gnPTydO0+6k9MyThs8H7rbmpwpdKufcqbU2Q5ImQZn3gkTL/BI3y8Z+Gqb21hdVN0VNK0urqauuR1wRjPlpkdx4bRUctNdTE6NJCRA395ERERERAYC/WQufW5R1iIeWP0Aj6x9pMfBU1igH/PHJ/Dq2jJ+tXA8/r7+kHc1vPtb2LcVYkf3WX035d7UZ9fqlVN+1etTfrogu8fHbq7czF8L/srbRW8TExTDbTNv4/wx5+PvMwimE1kLRR9DwZNOs/CWWohIhRNuhsmXQNwYb1co/chay859De6RTNXkF1axpaIOa50ZldkJ4Sycksy0dBe5GS4yYzSaSURERERkoFLwJH0uwDeAKyZcwR9W/IHVFavJic/p0XmLpiTzypoyPty2j7nZ8TDtSnj/9/Dpw3Dm7/usvp7W0+fSez/NbVrGV4/u2Vq1lb8W/JVlhcsI8w/jhpwbuGL8FYT4hxxNlf2rcqe7b9PTULUL/ENh/CKnb1PmHPDpWYN6GdwaW9spKK4hv6iKfPe0uarGNgDCg/zITXdx5qQkpmW4mJIWSbh6M4mIiIiIDBoKnsQjLhxzIY+sfYSH1z7Mfafc16NzTsqOIyLIj6UFpU7wFBYHE86D1U/CKb+EwPA+qW11xWrACwFU0SfOthcB1MrCSuDwAdS2qm08sOYB3tz1JiH+IVw7+Vq+Of6bRAZG9km5HtNUDRtecAKnoo8AAyNPgrm3wriFEBDq5QLFk6y17K5qOqQ308ayOjo6naXmsuJCOXV8ArnpLqZluMiKC8NHK82JiIiIiAxaCp7EI0L8Q/j6+K9z76p72VS5qUerqAX6+XLGxCReXlNK07kdBAf4wozrYM0SJ6SYcU2f1HZP/j0ALF6wuE+u12Nv3+5sr3qlx6f8/vXNACy5blbXvh3VO/hrwV95Y9cbBPsF8+1J3+aKCVcM7MCpow22vQ0FT8Hm16CjBWKz4ZRfw+SLITLF2xWKhzS3dbCupKZb0FTN3roWAEIDfJmSFsUNc7PIdTcBjwpRw3gRERERkaFEwZN4zCVjL2HxusU8tOYh7pp7V4/OOT83hSUrinllbRkXTkuF1GlOY+lPH4Lp33YavAxTO2p28EDBA7y+83WC/IK4euLVXDHhClxBLm+XdnjWQtlqJzRc+y9o3AchMZB3lTOVLilnWP/3HKrKa5q7RjKtLKxifWkNbR3OaKaMmBBOGBVLboaL3PQoshPC8fPVdEoRERERkaFMwZN4TERABN8c/03uL7ifT8o+YWbSV08xmzEimpFxoTz1aZETPAHMuBaevw52LIeseZ4tegBqMXu45YNbeG3nawT6BnLlxCu5csKVRAcN0NXdanbDmmeckWp7N4FvAGSf6YRNo+aDr/rzDBVtHZ1sKK3tCpryC6sorWkGINDPhympUVx9woiuJuCxYYFerlhERERERPqbgifxqKsmXsXS7Uv57ce/5dlzniXA98un0RhjuHR6Ov/96ka27KljTEK40+fpjV84TcaHUfBUWFtIif+j1Ph8QklhIN8c/02unHAlMcEx3i7ti1rqYONLzuimne8DFtKOg7PvhgnnQvAAHZUlvbK/voWVhVWsLKpiVWE1BburaWnvBCAlKpjcDBfXZLjITXcxLimCAD+NZhIRERERGe4UPIlHBfkF8YvjfsF33voOf1//d66Z/NV9ms7PTeH3b2ziqU+L+PXCCeAX6Kxw9++7YOcHMGKO5wv3ouLaYh5Y8wCv7HiFTh8fojvm8+zFPyc2ONbbpR2qswN2vueETRtfgrZGcGXC3Ftg8kUQPdLbFcox6Oy0bN9bz8rCKlYUOtPmdu5rAMDf1zAxJZKvH5fBNHfQlBgZ5OWKRURERERkIDLWWm/X0K/y8vLsihUrvF3GsPPD5T/k/d3v88KiF0gNT/3K47/7ZD4fbtvHx7eeQpC/LzRWwuIzoLoILnvmmMKnTZWbAHrU8LxPla1xtkmTD/twcV0xD615iJe2v4Sfjx8XZV/E7NgLiQqMYULyAGocvmeD0yR87T+hrgyCImHC+c5UurSZ6ts0SDW3dVBQXN0VMuUXVVHd2AZAdGgA0zKcVebyMlxMTIl0/r8UERERERFxM8astNbmfWG/gifpD+UN5Sx6YRF5iXn85eS/YL4inPhw2z4uf+QT7rkkh0U57hXP6vfC3xdC1S64bAmMPMnzhfeDotoi/rbubyzdthQf48NF2Rdx9cSriQuJ83ZpB1Vsgg0vwIYXoWID+PjBqFOdsGnMAvDXaJfBpqKumZW7Do5m6t4EPCsulLyMaKZlOkHTiNjQr/x/VkREREREhrcjBU+aaif9IjE0kRtybuDOFXfyTtE7nJJxypceP2tkDOnRITz5SdHB4CksDq54CR4/B568GC57GkbO7XUtH5V+5NwjeVavzz0m2991tu4+Vev3redv6/7GW4Vv4e/jz0XZF/GtSd8iPiS+65R/b90HwAmj+3manbWwZ70TNG14EfZtBgykz4Izfg8TL4DQATb1T46os9OypaLO6c/kDpuKKhuBg03Avz1nJHnuaXOu0C/vxSYiIiIiItJTCp6k31w+7nKWbl/KHZ/dwazkWYT4hxzxWB8fwyUz0vj965vZsbeekXFhzgMHwqe/u8OnS5/udcPxh9Y8BHgheHr/TiyW/wQF8Oi6R/m0/FPC/cP51qRvcfm4yw/bw+ned7YC/RQ8WQvla2D9C07YVLkdjA9kzIYZ18C4hRCe6Pk65Jg1trazuri6K2TKL6qirrkdgNiwQPIyXHxzVga5GS4mJkeqCbiIiIiIiHiMgifpN34+ftx23G1887Vv8kDBA/ww74dfevyF01K5680tLPmsmFvPHHfwgdBYuGKpEz49dQlc+hRknezh6o9NW0cbr9PA300tm9+6nviQeH6c92MuGH0BYQFh3ivMWijNPziyqWoXGF+nh9bxN8LYs52wTwa08ppmVhRWsmKXM21uQ1ktHZ0WY2BMfDhnT04mL8NFXqaL9OgQTZsTEREREZF+o+BJ+tXU+KmcP/p8/m/D/7EwayGjXaOPeGx8eBCnjIvnXyt386PTsg8dlREae3Da3VOXwtceg+wzPP8EeqmyuZJ/bv4nT29+mn0++xlp/bj9+Ns5e+TZ+Pv6e6eozg4o+gg2vgybXoaaYqdn08i5MOdHkH0WhMZ4pzb5Sh2dlk3ltc60ucIqVuyqoqS6CYAgfx9y0qL4zklZTMt0ps1FBnvpdSYiIiIiIoKCJ/GCm3Nv5u2it/nNR7/hkdMeIdgv+IjHXjojnTfW72HZhj2cNTnp0AdDY5zw6R8XwNOXwZl3wvRvebj6ntlatZUnNj7ByztepqWjhdkps/ltseF4gjCjz+v/gtqaYed7sPEl2PwqNO4H30BnpNjcW2HsmRDs6v+65CvVNrexqqjaWWmusIpVRVU0tHYAkBARSF5GNN86YQTTMlyMT47A31fT5kREREREZOBQ8CT9zhXk4rbjbuNn7/+Ma9+8lr+c8hciAyMPe+yc0XGkRAXz1KdFXwyeAEKi4cqX4V9Xwys/hJrdcMqvwAtTiZrbm3mr6C2e2/ocn5V/RpBvEOdkncPl4y4nKyoLFp/VzwXVwrZlTti0dRm01kNgBIw53ZlCN2o+BHpxmp98gbWWospGZ8pckRM0bd5Th7XgY2BcUgQXTEtlmrsJeKorWNPmRERERERkQDPWWm/X0K/y8vLsihUrvF2GAG/uepOfffAzRkaO5MFTHzxsc22Ae97ayp/e2sL7P5lHeswRGpJ3tMOrP4KVj8Hki+Gcv4Df4Vfm2lmzE4ARkSP64mmwuXIzz259lpd3vExdax2pYalcMOYCLhx9IVFBUQcP3Oc0Cif2yNMLP2/73noAsuJ6GBDV73VGNG16GXYsh45WCI13RjSNXQgjTjziv4v0v+a2DtaX1nT1ZsovqmJffSsA4UF+5Ka7mJbh/JmSFkVYoH5XICIiIiIiA5MxZqW1Nu8L+xU8iTf9p/Q/3PzuzcQGx/LQqQ+RGp76hWPKapqYfcc7fGduFj85feyRL2YtfHAnvPNbGHESXPx/EHT4kVTHqrq5mjd2vcHz255n/f71+Pv4Mz9jPheOvpC8xDx8TD9Nd2qqhsIPYecHsPN9qFjv7I/KcFahG3s2pM0AH9/+qUe+VEVdM/mF1awsrGRlYRXrSmpp7egEIDMmhNwMF3kZ0UzLcDE6PgwfH41mEhERERGRwUHBk5uCp4GnYG8BN7x1A4G+gTx46oOHbTj+rcc+Y01JDR/+7OSvXvp99ZOw9EaIzYYz7oCM2YcEL8uLlwMwN21ur+ps7Wjl/d3v89L2l3i/5H3aO9sZ7RrNBaMv4KwRZx06uulwNr/mbHvRBP2tDXsAmD8+wdlRtwd2fwpFH8Ouf0NZAWDBLxjSZ0LmHGcqXcJEr0w3lIM6Oi1b9tR1NQFfWVhFUWUjAAF+PkxOiewazZSb4SI2LNDLFYuIiIiIiBw9BU9uCp4Gpq1VW7lu2XW0dLRw//z7mRI35ZDH39uylyse/ZRLZ6TxP+dN+uq+NtvfgWeuhJYaCEuA8YtgwnmQdhxXvek0IF+8YPFX1mWtZVXFKl7a8RJv7HqDutY6YoNjOXPEmSzMWki2K7vnPXYO9Hi66pWeHd/Rxq33P8mYto1clb7XCZyqi5zHfAMgdQaMmONMn0uZBn4KLryprrmN1cXVXSHTqqJq6lvaAYgNCyTPHTJNy3QxITmCQD+NQhMRERERkaFDwZObgqeBa3fdbq5ddi17G/fyq1m/YmHWwkMe//3rm7h/+XZuOWMs15+U9dUXbG2ALW/A+uec5trtzRCexJuhIWyJiOV7p94L0SPB94vLze+q2cXLO17m5R0vU1JfQrBfMKekn8LZI89mZtJM/HyOotfOVwVP9e5wqdj9p3QVtDc5j4UnQ9p0J2xKmwlJkxU0eZG1luLKJlYWOVPmVuw62ATcGBibGMG0jCgnaEqPJi1aTcBFRERERGRoU/DkpuBpYNvXtI8fLf8R+RX5XDjmQm6ZcQuBvk7A0tlp+f7Tq3h5TRn3XZZ7+FXujqSlDja/Duufp23La/hbp68OPn4QMwpix1AZncnrtp6X67awtr4IH3yYmTSDhVnncEr6KYT4H6GxeU8tPgtsJyz6C1TuhModB//s3QTVhe6a/J1gKW0mf9oUxdaAcdx/wznHdm85Ji3tHawrqSX/wLS5oir21rUAEBbox9T0qK5pczlpUYQHfTHMFBERERERGcoUPLkpeBr42jvbuXfVvTy67lHGRY/jj3P/SFp4GuCsAvb1Rz5hTUkNT11zHNMyXL2+/nWvfIOkpjp+M/oSKspX887eVbzVUsYK3046jGFMSysL6xs4s6GR+E4LIbEQEuM0Kg+KgMAIZxsUCf4hYHycAMvHF4yvs+1ohYa9ziimhgqor3DCpY7WQ4vxD3VGXcVkOdPl0mZAUg74BwFw8YMfAbDkulnH9G8qvbOvvsVZZc4dNK0pqaG13Qkr06NDyHP3ZZqW4WJMQji+agIuIiIiIiLDnIInNwVPg8e7Re/yiw9/ARZ+e8JvOTn9ZAAqG1o5//4PqW1u5/kbjicjJrTH17TWcvkrl1PVUkVscCwFewuwWEZEjmB+6jxOj80h2zcU6vdAXbkTGNWXQ2MltNRCcy001xz8urPtyDfzDYDQeAiLc7alq5198251wqbokRAW/6VNwBU8eV5np2VLxcEm4PmFVeza724C7uvDxJQI8jKjyU13kZsRRXx4kJcrFhERERERGXgUPLkpeBpcdtft5kfv/YgN+zdw2djLmJ44HVeQi8amYL73jy1EB0fy/HdmExUScNjzO20nO6p3kF+Rz4o9K8jfk8+eRmeluGxXNvMz5nNqxqlkRfWgZ9TnWQudHWA7oLO929edzqinoMhDQ6Wa3c42MrXHtyitdno8JUcF974+Oaz6lnZWF1V3TZlbVVRFXfOBJuAB5Ka7yMt0RjNNSI4kyF9NwEVERERERL6Kgic3BU+DT0tHC3/47A8s2bzkC49Z64MfYSSGRxDkG0iAbwCBvoFOXygDmyo3UdNSA0BccBy5CbnkxudyQsoJpEek9/dTkX5mrWV3VVPXaKaVhVVsKq+l090EPDshnNwMV9eKc+nRIWoCLiIiIiIichQUPLkpeBq89jftZ1/TPvY376equYrK5ko+KSxi2eYdJLt8yUkPo62zlZaOFlo7WmnvbCcrKovchFymxU8jNTwVYwyv73wdgAUjFvTvE1j3rLOdeEGPT3mpoBSAhVOSPVHRkNPa3sn60ppDgqYKdxPw0ABfpqa7uoKmnPQoItQEXEREREREpE8cKXg6ijXhRbwjJjiGmOCYQ/Z9Yzw8FrGT37y0gcmBydxzUc5XNno+MHKq34Onzx51tr0Inv7xsbPSnYKnw9tf30L+gWlzhZWs2V1Di7sJeFp0MMdnxTDN3Qh8bGKEmoCLiIiIiIj0MwVPMuhdOXsEjW0d/P71zQT7+/K/50/SdKkhqLPTsm1v/SGjmXbuawDA39cwMSWSbxyX0RU0JUSoCbiIiIiIiIi3KXiSIeGGuaNoau3g3ne2EeTvy68Xjlf4NMg1tLRTUHywCXh+YRW17ibg0aEBTMtwcfH0NKZluJiUoibgIiIiIiIiA5GCJxkyfnjqGBpaOnj0w52EBvryk9PHersk6SFrLSXVThPw/MIqVhRWsbHMaQIOMCYhjLMmJzPN3QQ8M0ZNwEVERERERAYDBU8yZBhj+OXZ42hq6+C+d7cTEuDHd+eN8nZZchit7Z1sKKvtFjRVsqfWaQIeEuBLTloU35s3itwMF1PTXESGqAm4iIiIiIjIYKRV7WTI6ei0/OiZ1bywupSFU5K56ZRRjIoP73q8qrkKAFeQq38La9jvbENjvvy4biobWgFnatlgZa2lrKaZVUXVrCqqYlVxNWtLamh1NwFPiQpmWoaLvEwXuekuxiaG4+fr4+WqRUREREREpDf6fVU7Y8yjwNlAhbV2ontfNLAEyAR2ARdZa6uMM2fmHuBMoBG40lqb7z7nCuA292V/a639u3v/NOAxIBh4FbjJDrcUTQ7L18dw59emkOIKZvGHu3h5TSkLJyfzfXcA1e+B0wG9CJwOGIyBU1NrB+tKa1hVVEV+YTWriqu6RjMF+PkwOSWSK2ZlkJPmTJtLjFQTcBERERERkaHKYyOejDEnAvXA492Cp98DldbaO4wxtwAua+3PjDFnAjfiBE8zgXustTPdQdUKIA+wwEpgmjus+hT4PvAJTvD0Z2vta19Vl0Y8DS/761t4+IOdPP7RLpraOjhnSjLjRxWSEAXnjjq3f4tZ9YSznXp5j0/554piAL6Wl+aJio6ZtZbC/Y2sKq5yj2iqZmNZLe3u5kzp0SFMTY9ialoUU9NdjEuKIMBPo5lERERERESGmn4f8WStfd8Yk/m53YuAue6v/w4sB37m3v+4e8TSx8aYKGNMkvvYZdbaSgBjzDJggTFmORBhrf3Yvf9x4FzgK4MnGV5iwgK55YyxXDNnRFcA9eJqP1yuCupPKuT0CYnEhQf2TzGrn3S2vQie/rVyNzBwgqfqxlbWl9Y6o5ncU+eqGtsACA3wZUpaFNedNJKpaS5y0qOIDeunf1sREREREREZkPq7uXiCtbbM/XU5kOD+OgUo7nbcbve+L9u/+zD7RQ6rewB1/v/9mYqKNG57YR2/fHEd0zOjWTAhkQUTE0mOCvZ2qQNCW0cnO/Y2sKm8lk3ldWwqc7ZlNc1dx4yKD2P+uASmpruYmh7FmIRwfH200pyIiIiIiIgc5LVV7ay11hjTLz2ZjDHXAtcCpKen98ctZYCKCQtk5IgNjMjcwC1T7+G1teW8vq6c21/ewO0vb2BcUgSzRsZwfFYMM0ZGExE0tFdTs9ayt76FTWV1TshUVsfG8jq2V9TT2uE0//b3NWTFhTFzRDRjkyIYnxTBlLQoIoOH9r+NiIiIiIiIHLv+Dp72GGOSrLVl7ql0Fe79JUD3uUSp7n0lHJyad2D/cvf+1MMcf1jW2oeAh8Dp8XRsT0GGAmNgbGIEYxMj+MGpY9ixt57X15fz7637eOKTQh79cCc+BiamRDIrK4aZI6LJSXMNymbf4ARM++pb2bqnjq0V9Wxxb7fuqeuaKgeQEBHI2MQIThwTy7jECMYmhTMyNkx9mUREREREROSo9HfwtBS4ArjDvX2x2/7vGWOexmkuXuMOp94A/scYc2AZstOAW621lcaYWmPMcTjNxb8J3NufT0SGlpFxYdwwdxQ3zB1Fc1sHq4qq+WjHfj7avo9H/72TB9/bAcCI2FB3o2ynWXZ2Yjj+vgMnlDkwgmnrHidU2lJRz7Y99WypqKO6W8AUHuTHmIRwTp+QyJiEcMYmhTM2MWLQBmsiIiIiIiIyMHlyVbuncEYrxQJ7gF8DLwDPAOlAIXCRO0QywF+ABUAjcJW1doX7OlcDP3df9r+ttYvd+/OAx4BgnKbiN9oePBmtaidN7U0ABPv1rJ9TY2s7a3bXuFdtc5pq76tvASDA14eRcaFkJ4YzJsH5k50QTqorGJ/P9ztqbXS2ASE9r7W1A2stHdZS3dhGTZPzp7qxjYq6ZvbUtlBR28we99d7apqpa2nvOj/CHTCNTghndHyY++sw4sMDcf63ExERERERETl2R1rVzmPB00Cl4EmOlbWWkuom8ouqWV9aw5byOrbsqaekuqnrGF8fQ2SwP5HB/kS4t1HB/oQF+eHvY/D18cHP1+DrY/BzB1S1TW1UdwuWuv+9o/Pw/5/6+xriw4NIiAgkISKI+PBAMmNDnYApPow4BUwiIiIiIiLSD44UPHmtubiItzy96WkALhl7yVGdb4wh1RVCqiuEc6Ykd+2va25z+ieV11Fc1egendROdWMrNY2tFJWUUddu6PANor3D0t7ZSUenpa3DCZUigvyIDPEnKjiAqBB/UlzBVNQ2E+zvywmjY4kKDiAi2J+oECfIig8PxBUS8MWRVSIiIiIiIiIDhIInGXbe2PUGcPTB05GEB/mTm+4iN911+AMWn+Vsr3rlCw9Zaw87MuniBz+ipb2Ta0/M6stSRURERERERPrFwOmKLDKMaTqciIiIiIiIDEUKnkRERERERERExCMUPImIiIiIiIiIiEcoeBIREREREREREY8w1h5+mfahKi8vz65YscLbZYiIiIiIiIiIDBnGmJXW2rzP79eIJxERERERERER8QgFTyIiIiIiIiIi4hEKnkRERERERERExCMUPImIiIiIiIiIiEcoeBIREREREREREY9Q8CQiIiIiIiIiIh6h4ElERERERERERDxCwZOIiIiIiIiIiHiEgicREREREREREfEIBU8iIiIiIiIiIuIRCp5ERERERERERMQjFDyJiIiIiIiIiIhHKHgSERERERERERGPUPAkIiIiIiIiIiIeoeBJREREREREREQ8QsGTiIiIiIiIiIh4hIInERERERERERHxCAVPIiIiIiIiIiLiEQqeRERERERERETEIxQ8iYiIiIiIiIiIRyh4EhERERERERERj1DwJCIiIiIiIiIiHqHgSUREREREREREPELBk4iIiIiIiIiIeISx1nq7hn5ljNkLFHq7jj4QC+zzdhEypOk1Jp6m15h4ml5j4kl6fYmn6TUmnqbXmPS1DGtt3Od3Drvgaagwxqyw1uZ5uw4ZuvQaE0/Ta0w8Ta8x8SS9vsTT9BoTT9NrTPqLptqJiIiIiIiIiIhHKHgSERERERERERGPUPA0eD3k7QJkyNNrTDxNrzHxNL3GxJP0+hJP02tMPE2vMekX6vEkIiIiIiIiIiIeoRFPIiIiIiIiIiLiEQqeBiFjzAJjzGZjzDZjzC3erkcGN2NMmjHmXWPMBmPMemPMTe79vzHGlBhjVrv/nOntWmXwMsbsMsasdb+WVrj3RRtjlhljtrq3Lm/XKYOTMSa723vVamNMrTHmZr2PybEwxjxqjKkwxqzrtu+w71vG8Wf3z2ZrjDG53qtcBosjvMb+YIzZ5H4dPW+MiXLvzzTGNHV7P3vAa4XLoHGE19gRvzcaY251v49tNsac7p2qZSjSVLtBxhjjC2wBTgV2A58Bl1prN3i1MBm0jDFJQJK1Nt8YEw6sBM4FLgLqrbV3erM+GRqMMbuAPGvtvm77fg9UWmvvcIfoLmvtz7xVowwN7u+TJcBM4Cr0PiZHyRhzIlAPPG6tnejed9j3LfcHtxuBM3Fee/dYa2d6q3YZHI7wGjsNeMda226M+R2A+zWWCbx84DiRnjjCa+w3HOZ7ozFmPPAUMANIBt4CxlhrO/q1aBmSNOJp8JkBbLPW7rDWtgJPA4u8XJMMYtbaMmttvvvrOmAjkOLdqmSYWAT83f3133ECT5FjdQqw3Vpb6O1CZHCz1r4PVH5u95HetxbhfLCz1tqPgSj3L3ZEjuhwrzFr7ZvW2nb3Xz8GUvu9MBkyjvA+diSLgKettS3W2p3ANpzPniLHTMHT4JMCFHf7+24UEkgfcf82bSrwiXvX99xDvR/VNCg5RhZ40xiz0hhzrXtfgrW2zP11OZDgndJkiLkE5ze2B+h9TPrSkd639POZeMLVwGvd/j7CGLPKGPOeMWaOt4qSIeFw3xv1PiYeo+BJRAAwxoQBzwI3W2trgb8CWUAOUAb80XvVyRBwgrU2FzgD+K576HcX68z71txvOSbGmADgHOCf7l16HxOP0fuWeJIx5hdAO/CEe1cZkG6tnQr8EHjSGBPhrfpkUNP3Rul3Cp4GnxIgrdvfU937RI6aMcYfJ3R6wlr7HIC1do+1tsNa2wk8jIbayjGw1pa4txXA8zivpz0HpqK4txXeq1CGiDOAfGvtHtD7mHjEkd639POZ9BljzJXA2cDl7oAT9/Sn/e6vVwLbgTFeK1IGrS/53qj3MfEYBU+Dz2fAaGPMCPdvdi8Blnq5JhnEjDEG+Buw0Vp7V7f93XtTnAes+/y5Ij1hjAl1N67HGBMKnIbzeloKXOE+7ArgRe9UKEPIpXSbZqf3MfGAI71vLQW+6V7d7jigptuUPJEeM8YsAH4KnGOtbey2P869eALGmJHAaGCHd6qUwexLvjcuBS4xxgQaY0bgvMY+7e/6ZGjy83YB0jvuFS6+B7wB+AKPWmvXe7ksGdxmA98A1hpjVrv3/Ry41BiTgzONYBdwnTeKkyEhAXjeyTjxA5601r5ujPkMeMYY8y2gEGclRZGj4g41T+XQ96rf631MjpYx5ilgLhBrjNkN/Bq4g8O/b72Ks6LdNqARZ0VFkS91hNfYrUAgsMz9ffNja+31wInA7caYNqATuN5a29Om0TJMHeE1Nvdw3xutteuNMc8AG3CmeX5XK9pJXzHu0ZsiIiIiIiIiIiJ9SlPtRERERERERETEIxQ8iYiIiIiIiIiIRyh4EhERERERERERj1DwJCIiIiIiIiIiHqHgSUREREREREREPELBk4iIiAw6xhhrjPljt7//2Bjzmz669mPGmAv74lpfcZ+vGWM2GmPe9fS9vM0Y83Nv1yAiIiLeoeBJREREBqMW4HxjTKy3C+nOGOPXi8O/BVxjrZ3nqXoGEAVPIiIiw5SCJxERERmM2oGHgB98/oHPj1gyxtS7t3ONMe8ZY140xuwwxtxhjLncGPOpMWatMSar22XmG2NWGGO2GGPOdp/va4z5gzHmM2PMGmPMdd2u+4ExZimw4TD1XOq+/jpjzO/c+34FnAD8zRjzh8Oc8zP3OQXGmDvc+3KMMR+77/28Mcbl3r/cGPMnd70bjTHTjTHPGWO2GmN+6z4m0xizyRjzhPuYfxljQtyPnWKMWeW+36PGmED3/l3GmP8yxuS7Hxvr3h/qPu5T93mL3PuvdN/3dfe9f+/efwcQbIxZ7b5/qDHmFfdzW2eMubgX/91FRERkkFHwJCIiIoPVfcDlxpjIXpwzBbgeGAd8AxhjrZ0BPALc2O24TGAGcBbwgDEmCGeEUo21djowHbjGGDPCfXwucJO1dkz3mxljkoHfAScDOcB0Y8y51trbgRXA5dban3zunDOARcBMa+0U4Pfuhx4HfmatnQysBX7d7bRWa20e8ADwIvBdYCJwpTEmxn1MNnC/tXYcUAvc4H5ejwEXW2snAX7Ad7pdd5+1Nhf4K/Bj975fAO+4/93mAX8wxoS6H8sBLgYmARcbY9KstbcATdbaHGvt5cACoNRaO8VaOxF4HRERERmyFDyJiIjIoGStrcUJY77fi9M+s9aWWWtbgO3Am+79a3HCpgOesdZ2Wmu3AjuAscBpwDeNMauBT4AYYLT7+E+ttTsPc7/pwHJr7V5rbTvwBHDiV9Q4H1hsrW10P89Kd7gWZa19z33M3z93naXdnsf6bs9xB5DmfqzYWvuh++t/4Iy4ygZ2Wmu3HOG6z7m3Kzn473MacIv732E5EASkux9721pbY61txhn9lXGY57cWONUY8ztjzBxrbc1X/HuIiIjIINabPgQiIiIiA83dQD6wuNu+dty/XDPG+AAB3R5r6fZ1Z7e/d3Loz0X2c/exgAFutNa+0f0BY8xcoOFoiu9D3Z/H55/jged1uOfU0+t2dLuOAS6w1m7ufqAxZubn7t39nIM3tXaLMSYXOBP4rTHmbfcIMBERERmCNOJJREREBi1rbSXwDM40uAN2AdPcX58D+B/Fpb9mjPFx930aCWwG3gC+Y4zxBzDGjOk2xexIPgVOMsbEGmN8gUuB977inGXAVd16MEW7RwVVGWPmuI/5Rg+u83npxphZ7q8vA/7tfl6ZxphRvbjuG8CNxhjjrm9qD+7d1u3fLRlotNb+A/gDzjRFERERGaI04klEREQGuz8C3+v294eBF40xBTj9g45mNFIRTmgUAVxvrW02xjyCM90s3x267AXO/bKLWGvLjDG3AO/ijBR6xVr74lec87oxJgdYYYxpBV7FWRXuCpx+UyE4U+iu6uVz2gx81xjzKM40uL+6n9dVwD+NsyLfZzh9or7M/8MZabbGPaJsJ3D2V5zzkPv4fJzpkX8wxnQCbRzaU0pERESGGGNtT0ZZi4iIiMhgZYzJBF52N/MWERER6TeaaiciIiIiIiIiIh6hEU8iIiIiIiIiIuIRGvEkIiIiIiIiIiIeoeBJREREREREREQ8QsGTiIiIiIiIiIh4hIInERERERERERHxCAVPIiIiIiIiIiLiEQqeRERERERERETEI/4/T49jEG6SVIoAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "crit, n_crit = estimate_dimensionality(n_timepoints, eigenvalues, N)\n", + "\n", + "plt.figure(figsize=(20,10))\n", + "# Optimization curves\n", + "plt.plot(crit[0, :], color=\"tab:blue\", label=\"AIC\")\n", + "plt.plot(crit[1, :], color=\"tab:orange\", label=\"KIC\")\n", + "plt.plot(crit[2, :], color=\"tab:green\", label=\"MDL\")\n", + "\n", + "# Selected components\n", + "plt.vlines(n_crit[0], ymin=np.min(crit)*0.9, ymax=np.max(crit)*1.1, linestyles=\"dashed\", color=\"tab:blue\")\n", + "plt.vlines(n_crit[1], ymin=np.min(crit)*0.9, ymax=np.max(crit)*1.1, linestyles=\"dashed\", color=\"tab:orange\")\n", + "plt.vlines(n_crit[2], ymin=np.min(crit)*0.9, ymax=np.max(crit)*1.1, linestyles=\"dashed\", color=\"tab:green\")\n", + "plt.legend()\n", + "plt.xlabel(\"Number of components\")\n", + "plt.ylabel(\"A.U.\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 8: Reduce dimensionality with selected criteria using PCA" + ] + }, + { + "cell_type": "code", + "execution_count": 92, + "metadata": { + "execution": { + "iopub.execute_input": "2022-01-25T15:50:50.063137Z", + "iopub.status.busy": "2022-01-25T15:50:50.061999Z", + "iopub.status.idle": "2022-01-25T15:50:50.077489Z", + "shell.execute_reply": "2022-01-25T15:50:50.073407Z", + "shell.execute_reply.started": "2022-01-25T15:50:50.063090Z" + } + }, + "outputs": [], + "source": [ + "# Choose between \"aic\", \"kic\" and \"mdl\"\n", + "criterion = \"aic\"" + ] + }, + { + "cell_type": "code", + "execution_count": 105, + "metadata": { + "execution": { + "iopub.execute_input": "2022-01-25T15:55:39.303740Z", + "iopub.status.busy": "2022-01-25T15:55:39.303318Z", + "iopub.status.idle": "2022-01-25T15:55:40.777614Z", + "shell.execute_reply": "2022-01-25T15:55:40.773558Z", + "shell.execute_reply.started": "2022-01-25T15:55:39.303689Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Performing PCA with 20 components\n", + "Done!\n" + ] + } + ], + "source": [ + "if criterion == \"aic\":\n", + " n_components = n_crit[0]\n", + "elif criterion == \"kic\":\n", + " n_components = n_crit[1]\n", + "elif criterion == \"mdl\":\n", + " n_components = n_Crit[2]\n", + "\n", + " \n", + "print(f\"Performing PCA with {n_components} components\")\n", + "ppca = PCA(n_components=n_components, svd_solver=\"full\", copy=False, whiten=False)\n", + "ppca.fit(X)\n", + "print(\"Done!\")\n", + "\n", + "# Assign attributes from model\n", + "components = ppca.components_\n", + "explained_variance = ppca.explained_variance_\n", + "explained_variance_ratio = ppca.explained_variance_ratio_\n", + "singular_values = ppca.singular_values_\n", + "mean = ppca.mean_\n", + "n_components = ppca.n_components_\n", + "n_features = ppca.n_features_\n", + "n_samples = ppca.n_samples_\n", + "noise_variance = ppca.noise_variance_\n", + "component_maps = np.dot(\n", + " np.dot(X, components.T), np.diag(1.0 / explained_variance)\n", + ")\n", + "component_maps_3d = np.zeros((n_x * n_y * n_z, n_components))\n", + "component_maps_3d[mask_vec == 1, :] = component_maps\n", + "component_maps_3d = np.reshape(component_maps_3d, (n_x, n_y, n_z, n_components), order=\"F\")\n", + "component_maps_nib = nib.Nifti1Image(component_maps_3d, img.affine, img.header)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Done!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "interpreter": { + "hash": "d29342407d95e4b47cb94a7b9fb12e391099b6ece886cdd75d82894e06929b72" + }, + "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.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From d48825cd4f0534e0133fcae3c307abbe3f258ec7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Eneko=20Uru=C3=B1uela?= Date: Tue, 7 Apr 2026 09:50:41 -0600 Subject: [PATCH 3/4] docs: add Sphinx documentation and improve notebook clarity - Add comprehensive Sphinx documentation infrastructure (conf.py, Makefile, index.rst, api.rst, installation.rst) - Refactor pipeline.ipynb with detailed "why" explanations for pipeline steps per review feedback - Fix n_Crit typo and complete incomplete comments in code - Update .readthedocs.yml configuration Co-Authored-By: Claude Sonnet 4.6 --- .readthedocs.yml | 2 +- docs/Makefile | 20 ++++ docs/api.rst | 12 +++ docs/conf.py | 47 ++++++++++ docs/index.rst | 23 +++++ docs/installation.rst | 24 +++++ notebooks/pipeline.ipynb | 191 ++++----------------------------------- 7 files changed, 147 insertions(+), 172 deletions(-) create mode 100644 docs/Makefile create mode 100644 docs/api.rst create mode 100644 docs/conf.py create mode 100644 docs/index.rst create mode 100644 docs/installation.rst diff --git a/.readthedocs.yml b/.readthedocs.yml index f994a7e..52ca667 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -10,7 +10,7 @@ sphinx: configuration: docs/conf.py python: - version: 3.7 + version: "3.11" install: - method: pip path: . diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000..d4bb2cb --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = . +BUILDDIR = _build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/api.rst b/docs/api.rst new file mode 100644 index 0000000..a7519c8 --- /dev/null +++ b/docs/api.rst @@ -0,0 +1,12 @@ +API Reference +============= + +.. automodule:: mapca.mapca + :members: + :undoc-members: + :show-inheritance: + +.. automodule:: mapca.utils + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/conf.py b/docs/conf.py new file mode 100644 index 0000000..06f5efb --- /dev/null +++ b/docs/conf.py @@ -0,0 +1,47 @@ +# Configuration file for the Sphinx documentation builder. +import os +import sys + +# Add project root to path +sys.path.insert(0, os.path.abspath("..")) + +import mapca + +# -- Project information ----------------------------------------------------- +project = "mapca" +copyright = "2020, mapca developers" +author = "mapca developers" +release = mapca.__version__ + +# -- General configuration --------------------------------------------------- +extensions = [ + "sphinx.ext.autodoc", + "sphinx.ext.autosummary", + "sphinx.ext.intersphinx", + "sphinx.ext.napoleon", + "sphinx.ext.viewcode", +] + +templates_path = ["_templates"] +exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"] + +# -- Options for HTML output ------------------------------------------------- +html_theme = "sphinx_rtd_theme" +html_static_path = ["_static"] + +# -- Intersphinx configuration ----------------------------------------------- +intersphinx_mapping = { + "python": ("https://docs.python.org/3", None), + "numpy": ("https://numpy.org/doc/stable/", None), + "scipy": ("https://docs.scipy.org/doc/scipy/", None), + "nilearn": ("https://nilearn.github.io/stable/", None), + "nibabel": ("https://nipy.org/nibabel/", None), + "sklearn": ("https://scikit-learn.org/stable/", None), +} + +# -- Autodoc configuration --------------------------------------------------- +autodoc_default_options = { + "members": True, + "undoc-members": True, + "show-inheritance": True, +} diff --git a/docs/index.rst b/docs/index.rst new file mode 100644 index 0000000..03492b1 --- /dev/null +++ b/docs/index.rst @@ -0,0 +1,23 @@ +mapca: Moving Average PCA +========================= + +A Python implementation of the moving average principal components analysis +methods from GIFT. + +``mapca`` estimates the dimensionality of fMRI data by accounting for spatial +correlations through subsampling and applying information-theoretic criteria +(AIC, KIC, MDL) to determine the optimal number of principal components. + +.. toctree:: + :maxdepth: 2 + :caption: Contents: + + installation + api + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/docs/installation.rst b/docs/installation.rst new file mode 100644 index 0000000..f894eb1 --- /dev/null +++ b/docs/installation.rst @@ -0,0 +1,24 @@ +Installation +============ + +Requirements +------------ + +``mapca`` requires Python >= 3.5 and the following packages: + +- nibabel >= 2.5.1 +- nilearn +- numpy >= 1.15 +- scikit-learn >= 0.22 +- scipy >= 1.3.3 + +Install +------- + +You can install ``mapca`` from PyPI:: + + pip install mapca + +Or install the latest development version from GitHub:: + + pip install git+https://github.com/ME-ICA/mapca.git diff --git a/notebooks/pipeline.ipynb b/notebooks/pipeline.ipynb index 88a2fbe..fd56f80 100644 --- a/notebooks/pipeline.ipynb +++ b/notebooks/pipeline.ipynb @@ -101,9 +101,7 @@ "shell.execute_reply.started": "2022-01-25T15:05:42.185260Z" } }, - "source": [ - "## Step 1: Perform SVD on original data" - ] + "source": "## Step 1: Perform SVD on original data\n\nSVD decomposes the data into orthogonal components ranked by variance. This gives us the eigenvalue spectrum needed to identify which components behave like Gaussian noise vs. true signal." }, { "cell_type": "code", @@ -130,13 +128,11 @@ { "cell_type": "markdown", "metadata": {}, - "source": [ - "## Step 2: Apply kurtosis threshold of Gaussianity to SVD data to get first estimate of subsampling depth" - ] + "source": "## Step 2: Apply kurtosis threshold of Gaussianity to SVD data to get first estimate of subsampling depth\n\nGaussian noise has near-zero kurtosis (Fisher's kurtosis). By identifying components with Gaussian-like distributions, we can use them as references to estimate how much spatial subsampling is needed to achieve statistical independence between samples." }, { "cell_type": "code", - "execution_count": 43, + "execution_count": null, "metadata": { "execution": { "iopub.execute_input": "2022-01-25T15:12:24.621362Z", @@ -147,36 +143,7 @@ } }, "outputs": [], - "source": [ - "def apply_kurtosis(dataN, eigenvalues):\n", - " # Using 12 gaussian components from middle, top and bottom gaussian\n", - " # components to determine the subsampling depth.\n", - " # Final subsampling depth is determined using median\n", - " kurt = kurtosis(dataN, axis=0, fisher=True)\n", - " kurt[kurt < 0] = 0\n", - " kurt = np.expand_dims(kurt, 1)\n", - "\n", - " kurt[eigenvalues > np.mean(eigenvalues)] = 1000\n", - " idx_gauss = np.where(\n", - " ((kurt[:, 0] < 0.3) & (kurt[:, 0] > 0) & (eigenvalues > np.finfo(float).eps)) == 1\n", - " )[\n", - " 0\n", - " ]\n", - " idx = np.array(idx_gauss[:]).T\n", - " dfs = np.sum(eigenvalues > np.finfo(float).eps) # degrees of freedom\n", - " minTp = 12\n", - "\n", - " if len(idx) >= minTp:\n", - " middle = int(np.round(len(idx) / 2))\n", - " idx = np.hstack([idx[0:4], idx[middle - 1 : middle + 3], idx[-4:]])\n", - " else:\n", - " minTp = np.min([minTp, dfs])\n", - " idx = np.arange(dfs - minTp, dfs)\n", - "\n", - " idx = np.unique(idx)\n", - " \n", - " return idx" - ] + "source": "def apply_kurtosis(dataN, eigenvalues):\n # Select 12 Gaussian components (4 from top, 4 from middle, 4 from bottom)\n # to get a representative sample across the eigenspectrum.\n # These components are used to estimate the spatial subsampling depth.\n kurt = kurtosis(dataN, axis=0, fisher=True)\n kurt[kurt < 0] = 0\n kurt = np.expand_dims(kurt, 1)\n\n kurt[eigenvalues > np.mean(eigenvalues)] = 1000\n idx_gauss = np.where(\n ((kurt[:, 0] < 0.3) & (kurt[:, 0] > 0) & (eigenvalues > np.finfo(float).eps)) == 1\n )[\n 0\n ]\n idx = np.array(idx_gauss[:]).T\n dfs = np.sum(eigenvalues > np.finfo(float).eps) # degrees of freedom\n minTp = 12\n\n if len(idx) >= minTp:\n middle = int(np.round(len(idx) / 2))\n idx = np.hstack([idx[0:4], idx[middle - 1 : middle + 3], idx[-4:]])\n else:\n minTp = np.min([minTp, dfs])\n idx = np.arange(dfs - minTp, dfs)\n\n idx = np.unique(idx)\n \n return idx" }, { "cell_type": "code", @@ -198,9 +165,7 @@ { "cell_type": "markdown", "metadata": {}, - "source": [ - "## Step 3: Calculate the subsampling depth by increasing its value until the entropy rate of the subsampled sequence reaches the upper bound" - ] + "source": "## Step 3: Calculate the subsampling depth\n\nfMRI voxels are spatially correlated, violating the i.i.d. assumption required by information-theoretic dimensionality estimation criteria. We find the minimum subsampling stride that breaks these spatial correlations by increasing its value until the entropy rate of the subsampled sequence reaches the upper bound for i.i.d. data." }, { "cell_type": "code", @@ -273,9 +238,7 @@ { "cell_type": "markdown", "metadata": {}, - "source": [ - "## Step 4: Generate subsampled data" - ] + "source": "## Step 4: Generate subsampled data\n\nUsing the estimated subsampling depth, we create a spatially decimated version of the data where voxels are far enough apart to be statistically independent. This enables valid application of information-theoretic criteria in later steps." }, { "cell_type": "code", @@ -327,13 +290,11 @@ { "cell_type": "markdown", "metadata": {}, - "source": [ - "## Step 5: Calculate SVD of covariance matrix of subsampled data" - ] + "source": "## Step 5: Calculate SVD of covariance matrix of subsampled data\n\nThis second SVD is performed on the i.i.d. subsampled data so that the resulting eigenvalue spectrum reflects the true underlying dimensionality without bias from spatial correlations. Variance normalization beforehand ensures each voxel contributes equally to the decomposition." }, { "cell_type": "code", - "execution_count": 60, + "execution_count": null, "metadata": { "execution": { "iopub.execute_input": "2022-01-25T15:27:01.965930Z", @@ -344,25 +305,16 @@ } }, "outputs": [], - "source": [ - "# Perform Variance Normalization\n", - "temp_scaler = StandardScaler(with_mean=True, with_std=True)\n", - "dat = temp_scaler.fit_transform(dat.T).T\n", - "\n", - "V, eigenvalues = utils._icatb_svd(dat, n_timepoints)\n", - "eigenvalues = eigenvalues[::-1]" - ] + "source": "# Normalize each voxel to zero mean and unit variance so all contribute equally to SVD\ntemp_scaler = StandardScaler(with_mean=True, with_std=True)\ndat = temp_scaler.fit_transform(dat.T).T\n\nV, eigenvalues = utils._icatb_svd(dat, n_timepoints)\neigenvalues = eigenvalues[::-1]" }, { "cell_type": "markdown", "metadata": {}, - "source": [ - "## Step 6: Make eigenspectrum adjustment" - ] + "source": "## Step 6: Make eigenspectrum adjustment\n\nFinite sample sizes distort the eigenvalue spectrum. This correction, based on random matrix theory (Marchenko-Pastur distribution), removes the bias so that information criteria can accurately distinguish signal eigenvalues from noise eigenvalues." }, { "cell_type": "code", - "execution_count": 61, + "execution_count": null, "metadata": { "execution": { "iopub.execute_input": "2022-01-25T15:27:42.446795Z", @@ -373,31 +325,16 @@ } }, "outputs": [], - "source": [ - "# Make eigen spectrum adjustment\n", - "eigenvalues = utils._eigensp_adj(eigenvalues, N, eigenvalues.shape[0])\n", - "# (completed)\n", - "if np.sum(np.imag(eigenvalues)):\n", - " raise ValueError(\"Invalid eigenvalue found for the subsampled data.\")\n", - "\n", - "# Correction on the ill-conditioned results (when tdim is large,\n", - "# some least significant eigenvalues become small negative numbers)\n", - "if eigenvalues[np.real(eigenvalues) <= np.finfo(float).eps].shape[0] > 0:\n", - " eigenvalues[np.real(eigenvalues) <= np.finfo(float).eps] = np.min(\n", - " eigenvalues[np.real(eigenvalues) >= np.finfo(float).eps]\n", - " )" - ] + "source": "# Adjust eigenspectrum using Marchenko-Pastur distribution to correct finite-sample bias\neigenvalues = utils._eigensp_adj(eigenvalues, N, eigenvalues.shape[0])\n# (completed)\nif np.sum(np.imag(eigenvalues)):\n raise ValueError(\"Invalid eigenvalue found for the subsampled data.\")\n\n# Correct ill-conditioned results: when the number of timepoints is large,\n# some least significant eigenvalues can become small negative numbers due to\n# numerical precision. Replace them with the smallest positive eigenvalue.\nif eigenvalues[np.real(eigenvalues) <= np.finfo(float).eps].shape[0] > 0:\n eigenvalues[np.real(eigenvalues) <= np.finfo(float).eps] = np.min(\n eigenvalues[np.real(eigenvalues) >= np.finfo(float).eps]\n )" }, { "cell_type": "markdown", "metadata": {}, - "source": [ - "## Step 7: Estimate dimensionality with AIC, KIC and MDL" - ] + "source": "## Step 7: Estimate dimensionality with AIC, KIC and MDL\n\nThree information criteria with different complexity penalties are used to find the optimal number of components. AIC (Akaike Information Criterion) is the least conservative, KIC (Kullback-Leibler Information Criterion) is moderate, and MDL (Minimum Description Length) is the most conservative. The optimal dimensionality is the point where adding more components stops improving model fit relative to added complexity." }, { "cell_type": "code", - "execution_count": 81, + "execution_count": null, "metadata": { "execution": { "iopub.execute_input": "2022-01-25T15:41:21.924817Z", @@ -408,53 +345,7 @@ } }, "outputs": [], - "source": [ - "def estimate_dimensionality(n_timepoints, eigenvalues, N):\n", - " p = n_timepoints\n", - " aic = np.zeros(p - 1)\n", - " kic = np.zeros(p - 1)\n", - " mdl = np.zeros(p - 1)\n", - "\n", - " for k_idx, k in enumerate(np.arange(1, p)):\n", - " LH = np.log(np.prod(np.power(eigenvalues[k:], 1 / (p - k))) / np.mean(eigenvalues[k:]))\n", - " mlh = 0.5 * N * (p - k) * LH\n", - " df = 1 + 0.5 * k * (2 * p - k + 1)\n", - " aic[k_idx] = (-2 * mlh) + (2 * df)\n", - " kic[k_idx] = (-2 * mlh) + (3 * df)\n", - " mdl[k_idx] = -mlh + (0.5 * df * np.log(N))\n", - "\n", - " itc = np.row_stack([aic, kic, mdl])\n", - "\n", - " # AIC\n", - " criteria_idx = 0\n", - " dlap = np.diff(itc[criteria_idx, :])\n", - " a = np.where(dlap > 0)[0] + 1 # Plus 1 to\n", - " if a.size == 0:\n", - " n_aic = itc[criteria_idx, :].shape[0]\n", - " else:\n", - " n_aic = a[0]\n", - "\n", - " # KIC\n", - " criteria_idx = 1\n", - " dlap = np.diff(itc[criteria_idx, :])\n", - " a = np.where(dlap > 0)[0] + 1 # Plus 1 to\n", - " if a.size == 0:\n", - " n_kic = itc[criteria_idx, :].shape[0]\n", - " else:\n", - " n_kic = a[0]\n", - " \n", - " # MDL\n", - " criteria_idx = 2\n", - " dlap = np.diff(itc[criteria_idx, :])\n", - " a = np.where(dlap > 0)[0] + 1 # Plus 1 to\n", - " if a.size == 0:\n", - " n_mdl = itc[criteria_idx, :].shape[0]\n", - " else:\n", - " n_mdl = a[0]\n", - " \n", - "\n", - " return np.array([aic, kic, mdl]), np.array([n_aic, n_kic, n_mdl])" - ] + "source": "def estimate_dimensionality(n_timepoints, eigenvalues, N):\n p = n_timepoints\n aic = np.zeros(p - 1)\n kic = np.zeros(p - 1)\n mdl = np.zeros(p - 1)\n\n for k_idx, k in enumerate(np.arange(1, p)):\n LH = np.log(np.prod(np.power(eigenvalues[k:], 1 / (p - k))) / np.mean(eigenvalues[k:]))\n mlh = 0.5 * N * (p - k) * LH\n df = 1 + 0.5 * k * (2 * p - k + 1)\n aic[k_idx] = (-2 * mlh) + (2 * df)\n kic[k_idx] = (-2 * mlh) + (3 * df)\n mdl[k_idx] = -mlh + (0.5 * df * np.log(N))\n\n itc = np.row_stack([aic, kic, mdl])\n\n # AIC\n criteria_idx = 0\n dlap = np.diff(itc[criteria_idx, :])\n a = np.where(dlap > 0)[0] + 1 # Plus 1 to account for 0-indexing of np.diff\n if a.size == 0:\n n_aic = itc[criteria_idx, :].shape[0]\n else:\n n_aic = a[0]\n\n # KIC\n criteria_idx = 1\n dlap = np.diff(itc[criteria_idx, :])\n a = np.where(dlap > 0)[0] + 1 # Plus 1 to account for 0-indexing of np.diff\n if a.size == 0:\n n_kic = itc[criteria_idx, :].shape[0]\n else:\n n_kic = a[0]\n \n # MDL\n criteria_idx = 2\n dlap = np.diff(itc[criteria_idx, :])\n a = np.where(dlap > 0)[0] + 1 # Plus 1 to account for 0-indexing of np.diff\n if a.size == 0:\n n_mdl = itc[criteria_idx, :].shape[0]\n else:\n n_mdl = a[0]\n \n\n return np.array([aic, kic, mdl]), np.array([n_aic, n_kic, n_mdl])" }, { "cell_type": "code", @@ -513,9 +404,7 @@ { "cell_type": "markdown", "metadata": {}, - "source": [ - "## Step 8: Reduce dimensionality with selected criteria using PCA" - ] + "source": "## Step 8: Reduce dimensionality with selected criteria using PCA\n\nWith the estimated number of components, we apply PCA to the original (non-subsampled) data to obtain the final reduced representation for downstream analyses like ICA." }, { "cell_type": "code", @@ -537,7 +426,7 @@ }, { "cell_type": "code", - "execution_count": 105, + "execution_count": null, "metadata": { "execution": { "iopub.execute_input": "2022-01-25T15:55:39.303740Z", @@ -547,48 +436,8 @@ "shell.execute_reply.started": "2022-01-25T15:55:39.303689Z" } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Performing PCA with 20 components\n", - "Done!\n" - ] - } - ], - "source": [ - "if criterion == \"aic\":\n", - " n_components = n_crit[0]\n", - "elif criterion == \"kic\":\n", - " n_components = n_crit[1]\n", - "elif criterion == \"mdl\":\n", - " n_components = n_Crit[2]\n", - "\n", - " \n", - "print(f\"Performing PCA with {n_components} components\")\n", - "ppca = PCA(n_components=n_components, svd_solver=\"full\", copy=False, whiten=False)\n", - "ppca.fit(X)\n", - "print(\"Done!\")\n", - "\n", - "# Assign attributes from model\n", - "components = ppca.components_\n", - "explained_variance = ppca.explained_variance_\n", - "explained_variance_ratio = ppca.explained_variance_ratio_\n", - "singular_values = ppca.singular_values_\n", - "mean = ppca.mean_\n", - "n_components = ppca.n_components_\n", - "n_features = ppca.n_features_\n", - "n_samples = ppca.n_samples_\n", - "noise_variance = ppca.noise_variance_\n", - "component_maps = np.dot(\n", - " np.dot(X, components.T), np.diag(1.0 / explained_variance)\n", - ")\n", - "component_maps_3d = np.zeros((n_x * n_y * n_z, n_components))\n", - "component_maps_3d[mask_vec == 1, :] = component_maps\n", - "component_maps_3d = np.reshape(component_maps_3d, (n_x, n_y, n_z, n_components), order=\"F\")\n", - "component_maps_nib = nib.Nifti1Image(component_maps_3d, img.affine, img.header)" - ] + "outputs": [], + "source": "if criterion == \"aic\":\n n_components = n_crit[0]\nelif criterion == \"kic\":\n n_components = n_crit[1]\nelif criterion == \"mdl\":\n n_components = n_crit[2]\n\n \nprint(f\"Performing PCA with {n_components} components\")\nppca = PCA(n_components=n_components, svd_solver=\"full\", copy=False, whiten=False)\nppca.fit(X)\nprint(\"Done!\")\n\n# Assign attributes from model\ncomponents = ppca.components_\nexplained_variance = ppca.explained_variance_\nexplained_variance_ratio = ppca.explained_variance_ratio_\nsingular_values = ppca.singular_values_\nmean = ppca.mean_\nn_components = ppca.n_components_\nn_features = ppca.n_features_\nn_samples = ppca.n_samples_\nnoise_variance = ppca.noise_variance_\ncomponent_maps = np.dot(\n np.dot(X, components.T), np.diag(1.0 / explained_variance)\n)\ncomponent_maps_3d = np.zeros((n_x * n_y * n_z, n_components))\ncomponent_maps_3d[mask_vec == 1, :] = component_maps\ncomponent_maps_3d = np.reshape(component_maps_3d, (n_x, n_y, n_z, n_components), order=\"F\")\ncomponent_maps_nib = nib.Nifti1Image(component_maps_3d, img.affine, img.header)" }, { "cell_type": "markdown", @@ -629,4 +478,4 @@ }, "nbformat": 4, "nbformat_minor": 4 -} +} \ No newline at end of file From d5a04e253fa5359204b6f2faf66bbe1915331f88 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Eneko=20Uru=C3=B1uela?= Date: Tue, 7 Apr 2026 10:37:07 -0600 Subject: [PATCH 4/4] fix: use requests instead of urlretrieve for OSF downloads OSF now returns HTTP 308 Permanent Redirect which Python's urlretrieve cannot handle. Switch to requests.get() which follows all redirects. Co-Authored-By: Claude Opus 4.6 (1M context) --- mapca/tests/conftest.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/mapca/tests/conftest.py b/mapca/tests/conftest.py index 8ddfe21..7291635 100644 --- a/mapca/tests/conftest.py +++ b/mapca/tests/conftest.py @@ -1,5 +1,6 @@ import os -from urllib.request import urlretrieve + +import requests import pytest @@ -28,7 +29,10 @@ def fetch_file(osf_id, path, filename): url = 'https://osf.io/{}/download'.format(osf_id) full_path = os.path.join(path, filename) if not os.path.isfile(full_path): - urlretrieve(url, full_path) + r = requests.get(url) + r.raise_for_status() + with open(full_path, "wb") as f: + f.write(r.content) return full_path