{ "cells": [ { "cell_type": "markdown", "id": "a58b49d5034e4ab1", "metadata": { "collapsed": false }, "source": [ "Author: Julia Jurkowska" ] }, { "cell_type": "markdown", "id": "1f85e355ff43d7ab", "metadata": { "collapsed": false }, "source": [ "## Simulating and Localizing Multi-Source EEG Activity with MVPURE_py\n", "\n", "This tutorial walks through a complete pipeline for simulating brain activity from multiple sources and evaluating source localization accuracy using **MVPURE_py** framework.\n", "\n", "#### Introduction\n", "This tutorial provides a comprehensive workflow for:\n", "1. **Simulating realistic multi-source brain activity** using autoregressive models with controlled spatial and temporal characteristics.\n", "2. **Forward modeling** to project cortical sources to sensor space through lead field matrices.\n", "3. **Data processing** following standard EEG pipelines to prepare signals for inverse modeling.\n", "4. **Source localization** using **MVPURE_py** and LCMV-NAI.\n", "5. **Quantitative evaluation** of localization accuracy.\n", "\n", "Controlled simulations are essential for validating inverse solutions because they provide:\n", "- **Ground truth**: we know exactly where and when activity occurs.\n", "- **Parameter control**: we can systematically vary SNR, source configuration and coupling.\n", "\n", "This tutorial uses the **MVPURE_py** package, which extends MNE-Python with tools for multivariate source analysis.\n", "\n", "#### Tutorial overview\n", "We will:\n", "- Load anatomical and forward model data for ``sample_subject`` (downloadable from [Figshare](https://figshare.com/articles/dataset/Sample_subject_data_/30102451?file=57853861))\n", "- Define regions of interest (ROIs) using Desikan-Killiany parcellation\n", "- Generate synthetic source time courses with configurable temporal dynamics\n", "- Project cortical activity to EEG sensor space via the forward solution\n", "- Add controlled sensor noise and apply standard preprocessing\n", "- Localize and reconstruct sources using **MVPURE_py** and LCMV beamforming\n", "- Quantify localization error" ] }, { "cell_type": "code", "execution_count": 1, "id": "initial_id", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:11:08.283867Z", "start_time": "2026-06-27T19:11:04.381913Z" }, "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using pyvistaqt 3d backend.\n" ] } ], "source": [ "import mne\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "import os\n", "\n", "mne.viz.set_3d_backend('pyvistaqt')\n", "\n", "from mvpure_py import viz, utils, simulation" ] }, { "cell_type": "markdown", "id": "7f78c9ba99489892", "metadata": { "collapsed": false }, "source": [ "### Simulation Configuration\n", "\n", "#### Anatomical labeling and source selection\n", "We define simulation parameters that control both the spatial distribution and temporal characteristics of simulated brain activity. The workflow distinguishes between two types of sources:\n", "\n", "1. **Task-related activity** (post-stimulus): Simulated cortical responses following stimulus onset\n", "2. **Background activity**: Ongoing oscillatory and coupled activity present throughout the recording\n", "\n", "Source locations are specified using the Desikan-Killiany ``aparc`` parcellation atlas. \n", "\n", "#### Temporal dynamics parameters\n", "The simulation uses Multivariate AutoRegressive (MVAR) models to generate temporally structured source activity. Key parameters include:\n", "\n", "- **MVAR order** (``ORDER``, ``ORDER_BG``): Determines temporal memory of the autoregressive process. Higher orders capture longer-range temporal dependencies but increase model complexity.\n", "\n", "- **Coupling strength** (``COUPLING``, ``COUPLING_BG``): Controls cross-source interactions in the MVAR model. Values near 0 produce independent sources; values near 1 create strong interdependencies. Physiologically, this reflects effective connectivity between regions.\n", "\n", "- **Innovation noise** (``NOISE``, ``NOISE_BG``): Amplitude of the white noise innovations driving the MVAR process. This controls the stochasticity vs. deterministic structure in the generated signals.\n", "\n", "- **Temporal smoothing** (``GAUSSIAN_FILTER_SIGMA``, ``GAUSSIAN_FILTER_SIGMA_BG``): Standard deviation of Gaussian kernel applied to smooth generated time courses. This removes high-frequency artifacts and produces more physiologically realistic signals.\n", "\n", "#### Signal amplitude and SNR control\n", "\n", "- **Target standard deviation** (``TARGET_STD``, ``TARGET_STD_BG``): Sets the amplitude of source activity in Ampere-meters (Am).\n", "\n", "- **ERP scaling factor** (``ERP_FACTOR``): Amplitude multiplier for event-related potential (ERP) components. This controls the strength of phase-locked, stimulus-evoked responses relative to ongoing activity.\n", "\n", "- **SNR** (``SNR_DB``): Signal-to-noise ratio in decibels on source level.\n", "\n", "- **Sensor noise factor** (``NOISE_FACTOR_SENSORS``): Additional amplitude scaling for sensor-level noise, allowing independent control of source-level and measurement noise.\n", "\n", "#### Dataset size\n", "\n", "- **Number of epochs** (``N_EPOCHS``): Determines statistical power for covariance estimation. More epochs provide more stable covariance estimates, improving beamformer performance.\n", "\n", "- **Vertices per region** (``N_VERTICES_PER_ACTIVITY_LABEL``, ``N_VERTICES_PER_BG_LABEL``): Controls spatial extent of activity within each ROI.\n", "\n", "These parameters collectively define a realistic but controlled simulation environment. For methodological studies, parameters can be systematically varied to characterize performance across different scenarios.\n", "\n", "For more details see [paper](https://arxiv.org/abs/2509.14118)." ] }, { "cell_type": "code", "execution_count": 2, "id": "31a5f858aa733e51", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:11:29.550454Z", "start_time": "2026-06-27T19:11:29.541086Z" }, "collapsed": false }, "outputs": [], "source": [ "labels_to_use = [\n", " \"cuneus-lh\",\n", " \"cuneus-rh\",\n", " \"lateraloccipital-lh\",\n", " \"lateraloccipital-rh\",\n", " \"inferiorparietal-lh\",\n", " \"inferiorparietal-rh\",\n", " \"superiorparietal-lh\",\n", " \"superiorparietal-rh\",\n", " \"superiortemporal-lh\",\n", " \"superiortemporal-rh\",\n", " \"supramarginal-lh\",\n", " \"supramarginal-rh\",\n", " \"superiorfrontal-lh\",\n", " \"superiorfrontal-rh\",\n", " \"insula-lh\",\n", " \"insula-rh\",\n", " \"caudalanteriorcingulate-lh\",\n", " \"caudalanteriorcingulate-rh\"\n", "]\n", "\n", "labels_to_use_poststimuli = [\n", " 'lateraloccipital-rh',\n", " 'superiorfrontal-rh',\n", " 'caudalanteriorcingulate-lh',\n", " 'caudalmiddlefrontal-lh',\n", " 'lateraloccipital-lh'\n", "]" ] }, { "cell_type": "code", "execution_count": 3, "id": "b0ed363b01d6ead7", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:11:31.616998Z", "start_time": "2026-06-27T19:11:31.610839Z" }, "collapsed": false }, "outputs": [], "source": [ "N_VERTICES_PER_BG_LABEL = 1\n", "N_VERTICES_PER_ACTIVITY_LABEL = 1\n", "N_EPOCHS = 100\n", "SNR_DB = 1.0\n", "ERP_FACTOR = 15e-9\n", "ORDER_BG = 7\n", "ORDER = 3\n", "TARGET_STD_BG = 15e-9\n", "TARGET_STD = 15e-9\n", "COUPLING_BG = 0.5\n", "COUPLING = 0.8\n", "NOISE_BG = 1.0\n", "NOISE = 1.0\n", "GAUSSIAN_FILTER_SIGMA_BG = 3\n", "GAUSSIAN_FILTER_SIGMA = 3\n", "N_DOMINANT_EIGVALS = 2\n", "NOISE_FACTOR_SENSORS = 0.3" ] }, { "cell_type": "markdown", "id": "8ee4a8da5957a1b9", "metadata": { "collapsed": false }, "source": [ "### Loading real data for sample subject\n", "We load:\n", "- **Epoched EEG data**: Used here only to extract metadata (sampling rate, time windows, channel names)\n", "- **Forward solution**: Contains the lead field matrix and source space geometry\n", "\n", "The forward solution has already been converted to **fixed orientation** (surface-normal constraint), which reduces degrees of freedom from 3 to 1 per source location." ] }, { "cell_type": "code", "execution_count": 4, "id": "542d1222a48348e7", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:12:01.148204Z", "start_time": "2026-06-27T19:11:59.672515Z" }, "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Reading /Volumes/UMK/oddball/subjects/sample_subject/_eeg/_pre/sample_subject_oddball-epo.fif ...\n", " Found the data of interest:\n", " t = -199.22 ... 800.78 ms\n", " 0 CTF compensation matrices available\n", "Not setting metadata\n", "621 matching events found\n", "No baseline correction applied\n", "0 projection items activated\n", "Reading forward solution from /Volumes/UMK/oddball/subjects/sample_subject/forward/sample_subject_ico4-fwd.fif...\n", " Reading a source space...\n", " [done]\n", " Reading a source space...\n", " [done]\n", " 2 source spaces read\n", " Desired named matrix (kind = 3523 (FIFF_MNE_FORWARD_SOLUTION_GRAD)) not available\n", " Read EEG forward solution (5124 sources, 128 channels, free orientations)\n", " Source spaces transformed to the forward solution coordinate frame\n", " No patch info available. The standard source space normals will be employed in the rotation to the local surface coordinates....\n", " Changing to fixed-orientation forward solution with surface-based source orientations...\n", " [done]\n" ] } ], "source": [ "# Reading real data for sample subject\n", "subject = \"sample_subject\"\n", "subjects_dir = \"subjects\"\n", "mne.set_config('SUBJECTS_DIR', subjects_dir)\n", "\n", "epoched = mne.read_epochs(os.path.join(subjects_dir, subject, \"_eeg\", \"_pre\", f\"{subject}_oddball-epo.fif\"))\n", "\n", "# Read forward\n", "forward_path = os.path.join(subjects_dir, subject, \"forward\", f\"{subject}_ico4-fwd.fif\")\n", "fwd_vector = mne.read_forward_solution(forward_path)\n", "fwd = mne.convert_forward_solution(\n", " fwd_vector,\n", " surf_ori=True,\n", " force_fixed=True,\n", " use_cps=True\n", ")\n", "\n", "# Leadfield matrix\n", "leadfield = fwd[\"sol\"][\"data\"]\n", "\n", "# Source positions extracted from forward model\n", "src = fwd[\"src\"]" ] }, { "cell_type": "markdown", "id": "19977cb321655058", "metadata": { "collapsed": false }, "source": [ "We extract key timing parameters:\n", "- sampling frequency (``sfreq``)\n", "- time window (``tmin`` to ``tmax``)\n", "- time vector (``times``) with its mask corresponding to sensory processing activity during oddball paradigm (``post_mask``)" ] }, { "cell_type": "code", "execution_count": 5, "id": "c42594da2dee0e90", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:12:10.005982Z", "start_time": "2026-06-27T19:12:09.994818Z" }, "collapsed": false }, "outputs": [], "source": [ "sfreq = epoched.info['sfreq']\n", "tmin, tmax = epoched.tmin, 0.25\n", "tmin_sensory = 0.05\n", "tmax_sensory = 0.2\n", "times = np.arange(tmin, tmax, 1/sfreq)\n", "n_times = len(times)\n", "\n", "post_mask = (times >= tmin_sensory) & (times <= tmax_sensory)" ] }, { "cell_type": "markdown", "id": "5ec3b515eb2e8042", "metadata": { "collapsed": false }, "source": [ "### Generating source locations\n", "Rather than activating entire anatomical regions uniformly, we randomly select individual vertices within each label. \n", "For this simulation:\n", "- **Background sources**: Selected from the full ``labels_to_use`` set (18 bilateral regions)\n", "- **Task-related sources**: Selected from ``labels_to_use_poststimuli`` (6 posterior regions)\n", "\n", "This configuration mimics an oddball paradigm where visual stimuli engage posterior sensory areas while a broader network maintains sustained attention." ] }, { "cell_type": "code", "execution_count": 6, "id": "2cf5b58c05cbe04e", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:13:45.479045Z", "start_time": "2026-06-27T19:13:43.101856Z" }, "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Reading labels from parcellation...\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n", "Reading labels from parcellation...\n", " read 1 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/lh.aparc.annot\n", " read 0 labels from /Volumes/UMK/oddball/subjects/sample_subject/label/rh.aparc.annot\n" ] } ], "source": [ "labels_info = simulation.get_random_vertices(\n", " n_vertices_per_label_bg=N_VERTICES_PER_BG_LABEL,\n", " n_vertices_per_poststimuli_label=N_VERTICES_PER_ACTIVITY_LABEL,\n", " noise_labels=labels_to_use,\n", " poststimuli_labels=labels_to_use_poststimuli,\n", " subject=subject,\n", " subjects_dir=subjects_dir,\n", " src=src,\n", ")\n", "\n", "noise_vertices, poststimuli_vertices = simulation.split_vertices(labels_info, labels_to_use_poststimuli)\n", "\n", "# Mapping to leadfield space\n", "leadfield_subset_indices_noise, lh_vert_to_lf_noise, rh_vert_to_lf_noise = utils.translation.transform_vertices_to_leadfield_indices(\n", " vertices=noise_vertices,\n", " src=src,\n", " hemi=\"both\",\n", " include_mapping = True\n", ")\n", "\n", "leadfield_subset_poststimuli_indices, lh_vert_to_lf_stimuli, rh_vert_to_lf_stimuli = utils.translation.transform_vertices_to_leadfield_indices(\n", " vertices=poststimuli_vertices,\n", " src=src,\n", " hemi=\"both\",\n", " include_mapping=True\n", ")\n", "\n", "leadfield_subset_indices = np.unique(np.concatenate((leadfield_subset_indices_noise, leadfield_subset_poststimuli_indices)))\n", "\n", "leadfield_to_label_mapping = simulation.assign_label_to_leadfield_index(\n", " labels_info=labels_info,\n", " lh_vert_to_lf=lh_vert_to_lf_noise | lh_vert_to_lf_stimuli,\n", " rh_vert_to_lf=rh_vert_to_lf_noise | rh_vert_to_lf_stimuli\n", ")" ] }, { "cell_type": "markdown", "id": "eb4ec2495a42bf3b", "metadata": { "collapsed": false }, "source": [ "### Simulating source activity\n", "Now we generate synthetic source time series.\n", "Our simulation generates two activity components:\n", "\n", "**1. Background activity (entire epoch)**\n", "- Higher MVAR order (``ORDER_BG`` = 7): Captures slower oscillatory dynamics\n", "- Moderate coupling (``COUPLING_BG`` = 0.5): Weak to moderate inter-regional interactions\n", "- Active throughout the epoch: Represents ongoing brain state\n", "\n", "**2. Task-related activity (post-stimulus)**\n", "- Lower MVAR order (``ORDER`` = 3): Transient, event-related dynamics \n", "- Strong coupling (``COUPLING`` = 0.8): Coordinated network response\n", "- Restricted to post-stimulus period: Evoked sensory processing\n", "- Includes ERP component: Phase-locked response to stimulus onset\n", "\n", "After MVAR generation:\n", "1. **Gaussian smoothing**: Reduces high-frequency noise, mimics dendritic integration\n", "2. **Amplitude scaling**: Normalizes to target RMS amplitude (``TARGET_STD``)\n", "3. **ERP addition**: Superimposes stereotyped waveform for evoked components\n", "\n", "This yields ``n_epochs × n_sources × n_times`` array of source time courses with realistic spatiotemporal structure." ] }, { "cell_type": "code", "execution_count": 7, "id": "7c53faca3ec83d9c", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:13:58.842533Z", "start_time": "2026-06-27T19:13:57.002647Z" }, "collapsed": false }, "outputs": [], "source": [ "X_epochs = simulation.simulate_source_epochs(\n", " n_epochs=N_EPOCHS,\n", " lf_subset_indices=leadfield_subset_indices,\n", " n_times=n_times,\n", " bg_lf_subset_indices=leadfield_subset_indices_noise,\n", " poststimuli_mask=post_mask,\n", " poststimuli_lf_subset_indices=leadfield_subset_poststimuli_indices,\n", " lf_to_label=leadfield_to_label_mapping,\n", " sfreq=sfreq,\n", " erp_factor=ERP_FACTOR,\n", " snr_db=SNR_DB,\n", " order_bg=ORDER_BG,\n", " order=ORDER,\n", " target_std_bg=TARGET_STD_BG,\n", " target_std=TARGET_STD,\n", " coupling_bg=COUPLING_BG,\n", " coupling=COUPLING,\n", " noise_bg=NOISE_BG,\n", " noise=NOISE,\n", " gaussian_filter_sigma=GAUSSIAN_FILTER_SIGMA,\n", " gaussian_filter_sigma_bg=GAUSSIAN_FILTER_SIGMA_BG,\n", " n_dominant_eigvals=N_DOMINANT_EIGVALS,\n", " seed=42\n", ")" ] }, { "cell_type": "markdown", "id": "db26fd1d6b3ff539", "metadata": { "collapsed": false }, "source": [ "### Projecting to sensor space\n", "\n", "Source activity is projected to EEG electrodes using the forward solution. The simulated sensor data is packaged into an MNE ``Epochs`` object. This allows the simulated data to be processed identically to real recordings." ] }, { "cell_type": "code", "execution_count": 8, "id": "752a6c0591d740bf", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:14:09.269669Z", "start_time": "2026-06-27T19:14:06.735035Z" }, "collapsed": false }, "outputs": [], "source": [ "sim_epochs = simulation.simulate_sensor_epochs(\n", " X_epochs=X_epochs,\n", " leadfield=leadfield,\n", " lf_subset_indices=leadfield_subset_indices,\n", " src=src,\n", " tmin=tmin,\n", " sfreq=sfreq,\n", " noise_factor=NOISE_FACTOR_SENSORS,\n", " info=epoched.info\n", ")\n", "# sim_epochs.plot()" ] }, { "cell_type": "markdown", "id": "2f0260b5c3ddb092", "metadata": { "collapsed": false }, "source": [ "### Preprocessing\n", "\n", "We apply basic preprocessing:" ] }, { "cell_type": "code", "execution_count": 9, "id": "492d53650000231f", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:14:21.724234Z", "start_time": "2026-06-27T19:14:17.379239Z" }, "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Setting up band-pass filter from 1 - 45 Hz\n", "\n", "IIR filter parameters\n", "---------------------\n", "Butterworth bandpass zero-phase (two-pass forward and reverse) non-causal filter:\n", "- Filter order 16 (effective, after forward-backward)\n", "- Cutoffs at 1.00, 45.00 Hz: -6.02, -6.02 dB\n" ] } ], "source": [ "sim_epochs = sim_epochs.filter(l_freq=1, h_freq=45, method=\"iir\")\n", "# sim_epochs.plot()" ] }, { "cell_type": "markdown", "id": "fd95da1bc3ef9d0b", "metadata": { "collapsed": false }, "source": [ "Then re-reference the data:" ] }, { "cell_type": "code", "execution_count": 10, "id": "34cbc0c551ce6449", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:14:26.154199Z", "start_time": "2026-06-27T19:14:26.040960Z" }, "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "EEG channel type selected for re-referencing\n", "Adding average EEG reference projection.\n", "1 projection items deactivated\n", "Average reference projection was added, but has not been applied yet. Use the apply_proj method to apply it.\n", "Created an SSP operator (subspace dimension = 1)\n", "1 projection items activated\n", "SSP projectors applied...\n" ] }, { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", "\n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", "\n", "\n", "\n", " \n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", "\n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "
\n", " \n", " \n", " General\n", "
MNE object typeEpochsArray
Measurement date2018-10-31 at 11:19:57 UTC
ParticipantUnknown
ExperimenterUnknown
\n", " \n", " \n", " Acquisition\n", "
Total number of events100
Events counts\n", " \n", " 1: 100\n", " \n", " \n", "
Time range-0.199 – 0.246 s
Baselineoff
Sampling frequency256.00 Hz
Time points115
MetadataNo metadata set
\n", " \n", " \n", " Channels\n", "
EEG\n", " \n", "\n", " \n", "
Head & sensor digitization131 points
\n", " \n", " \n", " Filters\n", "
Highpass1.00 Hz
Lowpass40.00 Hz
Projections\n", " \n", " Average EEG reference (on)\n", " \n", " \n", "
" ], "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sim_epochs_ref = sim_epochs.set_eeg_reference('average', projection=True)\n", "sim_epochs_ref.apply_proj()" ] }, { "cell_type": "markdown", "id": "c48be39f6e227ee6", "metadata": { "collapsed": false }, "source": [ "This mimics standard EEG preprocessing pipelines." ] }, { "cell_type": "markdown", "id": "8eb958f77a319a5a", "metadata": { "collapsed": false }, "source": [ "### Creating ground truth source estimates\n", "\n", "We construct source estimated (``stc``) corresponding to simulated activity:" ] }, { "cell_type": "code", "execution_count": 11, "id": "25a61bd76e00e52e", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:14:36.236842Z", "start_time": "2026-06-27T19:14:36.225061Z" }, "collapsed": false }, "outputs": [], "source": [ "sim_stc = simulation.add_simulated_epochs_to_stc(\n", " X_epochs=X_epochs,\n", " src=src,\n", " n_times=n_times,\n", " lf_subset_indices=leadfield_subset_indices,\n", " tmin=tmin,\n", " sfreq=sfreq,\n", " subject=subject\n", ")\n", "\n", "# sim_stc.plot(hemi=\"both\")" ] }, { "cell_type": "markdown", "id": "1c8ef3320c03d010", "metadata": { "collapsed": false }, "source": [ "This serves as **ground truth** for later evaluation." ] }, { "cell_type": "markdown", "id": "cb225e96e5ab828e", "metadata": { "collapsed": false }, "source": [ "### Localization Pipeline" ] }, { "cell_type": "markdown", "id": "8a4134cc1695f40b", "metadata": { "collapsed": false }, "source": [ "#### Covariance estimation\n", "We compute noise and data covariance matrices as these are essential inputs for beamforming." ] }, { "cell_type": "code", "execution_count": 12, "id": "9ca8f8b5cba0d052", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:15:54.028748Z", "start_time": "2026-06-27T19:15:53.812470Z" }, "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Created an SSP operator (subspace dimension = 1)\n", " Setting small EEG eigenvalues to zero (without PCA)\n", "Reducing data rank from 128 -> 127\n", "Estimating covariance using EMPIRICAL\n", "Done.\n", "Number of samples used : 5200\n", "[done]\n", " Created an SSP operator (subspace dimension = 1)\n", " Setting small EEG eigenvalues to zero (without PCA)\n", "Reducing data rank from 128 -> 127\n", "Estimating covariance using EMPIRICAL\n", "Done.\n", "Number of samples used : 3900\n", "[done]\n" ] } ], "source": [ "noise_cov = mne.compute_covariance(\n", " sim_epochs,\n", " tmin=-0.2,\n", " tmax=0,\n", " method=\"empirical\"\n", ")\n", "\n", "data_cov = mne.compute_covariance(\n", " sim_epochs,\n", " tmin=0.05,\n", " tmax=0.2,\n", " method=\"empirical\"\n", ")\n", "\n", "sim_evoked = sim_epochs.average()\n", "# Subset signal for given time range\n", "signal_sim = sim_evoked.copy().crop(\n", " tmin=0.05,\n", " tmax=0.2\n", ")" ] }, { "cell_type": "code", "execution_count": 13, "id": "b6d8cd6fed25e65b", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:15:56.255403Z", "start_time": "2026-06-27T19:15:56.237916Z" }, "collapsed": false }, "outputs": [], "source": [ "# Optional visualization\n", "\n", "# signal_sim.plot()\n", "# plt.show()" ] }, { "cell_type": "markdown", "id": "12de118138ba2752", "metadata": { "collapsed": false }, "source": [ "#### Eigenvalue spectrum visualization\n", "The eigenvalue spectrum of the whitened data covariance reveals the signal subspace structure:\n", "- Large eigenvalues: Capture coherent brain signals\n", "- Small eigenvalues: Dominated by noise" ] }, { "cell_type": "code", "execution_count": 14, "id": "a99d47b0bdde9a8", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:16:05.528476Z", "start_time": "2026-06-27T19:16:05.317382Z" }, "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8gAAAIUCAYAAADG2QXiAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAARdtJREFUeJzt3QmYnFWZP+zTWehASHdWCIEQAi6YILtAUGFUEBhkUQcVg0ZgRCGIyLiADktEDMjo5QwwCDKCY1gUBQFHmAFk+ZCdEAQCyCbwDwmBhHQngSSQ1Hedt6i2uvNWd1Uvtd73dRXd9VZV19vVD53+1TnnOU2ZTCYTAAAAoMENqvQJAAAAQDUQkAEAAEBABgAAgCwBGQAAAARkAAAAyBKQAQAAQEAGAACALAEZAAAABGQAAADIEpABAABAQAaAzi688MKw8847h6FDh4Yzzjij0qcDAJSRgAwAeTbbbLMkGH/605+u9KkAAGUmIANAnkMPPTQcfPDBYeTIkQPy9S+77LLQ1NQU/va3v4V6VO/fHwD1TUAGoOrCVe4yZMiQsPnmm4cvfelLYcGCBan3HTZs2Hq3Rf/wD/8QtttuuzKePQBQ6wRkAKrO97///fCrX/0q/OxnPwsHHHBAmDNnTth7773DqlWr1rvv6tWrw9lnnx1qxRe+8IXw5ptvhkmTJlX6VACALgRkAKpODMVHHHFE+Od//udwySWXhG9+85vh2WefDddff/16991xxx3Dz3/+8/Dyyy/3+HU/9KEPdRqhzr/867/+ayiHwYMHJ6Pe8TkBgOoiIANQ9T784Q8nH2NI7uq73/1uWLt2bVGjyHfddVfIZDKplx/84Ad9Ps841fuoo44Km266aWhubg5Tp04Nv/jFL3pco3v77beHXXfdNQnO22yzTbjooouSRmFpIbqY58g99plnnkmmp8f11K2treHII48Mb7zxRsf9fvvb3yb3u+OOO9Z7nngO8bbHHnssuf7CCy+E4447Lrz3ve8NG264YRgzZkw47LDDelxrHJ9/q622Wu94X76/5cuXhxNPPDH5uvE+m2yySdh3333D3Llzuz0XAOjJkB7vAQAVlgtho0aNWu+2yZMnhy9+8YvJKPLJJ58cJkyY0Kfnevvtt5NLDN3xY5zWHbd8iiO/3XnllVfCHnvskYS+448/PowbNy7ceOON4eijjw7t7e1JoEvz8MMPh/333z/pnj1r1qzkeeMU8/j4vj7HZz7zmeT1mT17dhIe42h8DJPnnHNOcvuBBx4YNt544/Cb3/wmmcKe79e//nUSTnPruB944IFw9913h8997nNhiy22SH4mcUusuNZ7/vz5YaONNir5te7t9/fVr341CffxPlOmTAlLlixJ3vx44oknki26AKDXMgBQJS699NJM/Kfplltuybz66quZl156KfPb3/42M27cuExzc3Nyvet9H3jggcyzzz6bGTJkSOaEE07ouH3vvffOTJ06teRzOP3005Ovm3+Jz9WTo48+OrPZZptlXnvttU7HP/e5z2VaW1szb7zxRqfzfv7555PrBx10UGajjTbKLFiwoOMxTz/9dPL9dP1nutjnyH0PRx11VKf7ffKTn8yMGTOm07HDDz88s8kmm2TefvvtjmMLFy7MDBo0KPP973+/41jua+e75557kuf57//+745jXb+/GTNmZCZNmrTeY3Pn2JvvL34+c+bM9b4mAPSVKdYAVJ199tknGT2cOHFi+Kd/+qcwfPjwZP1xHLlMs/XWWyfNry6++OKwcOHCPj13nPrbdfp1nCbcnXif3/3ud+Gggw5KPn/ttdc6Lvvtt19oa2tLnf4bR4tvueWWZGup/JHvd73rXck67L4+Rxxp7TpVPY62xtHYnM9+9rNh8eLFyTTvnDg6u27duuS2nDitOuett95Kvk48zzh9uz+mNpfy/cXnvO+++4padw4ApRCQAag6F1xwQbj55puToPaP//iPSUiKa027E5tsxSnRleho/eqrr4Zly5YlAT0G+/xLXPcbxRDaVTwWO1rHoNlV12O9eY4tt9yy0/XcFPXXX3+941ic3h3XJ8cp1Tnx89j87D3veU/HsXiep512WvKmRfxZjB07NnnueE4xvPZVKd/fj370o2RtdDyX3XbbLXlT47nnnuvzOQCANcgAVJ0YemLTqiiOrsbu05///OfDU089layZLTSKHDtfx4AV1yKXUxxtjeLzz5gxI/U+22+/fdmfo9C66ThCmxPDbnyNr7322vCf//mfyTrgP//5z+GHP/xhp8d87WtfC5deemmyDnjatGlJqI5rheOa5Ny5pSnUrTuOnvf2+4trq+NoeDzn//u//wvnnntusq76mmuuWW/kHQBKISADUNViyItNpj7ykY+E888/v9vwG0eR457JuSZU5RJHOUeMGJGEvjg9vFixYVbsXB27TXfV9Vhvn6MYcSr1L3/5y3Drrbcmja5igM6fXh3F0fwYXH/84x93HIsNzOKob3fiqHXafWJX7L58f7GpWeyqHS9xZDk25zrrrLMEZAD6xBRrAKpe7JQcR5V/+tOfJqGskLhFUhyBjFsULVq0qKwh/tOf/nSyhja3LVLX6cOFHhfD4O9///tO62ljOI7dm/vjOYoRz2H06NHJ1Op4ia917H7d9fnzR56j8847b72R4LSfSZyC/Ze//KXjWFwnHkd/e/P9xefrOqU7vtEQ13CvXr26hO8aANZnBBmAmvCtb30r2Xc37iPctflUvu9973vhV7/6VTIdO25TVC5x7fNtt90Wdt999/DlL3852X5o6dKlSWOp2Igrfp4mrp+N04Q/+MEPhmOPPTYJgHGkPG6vNG/evH55jp7Ebaw+9alPhauuuiqsXLky/Nu//dt69/nEJz6RvK5xanV83nvuuSd5zrgfcnfiFOzvfOc74ZOf/GQ44YQTkn2Y4/ZQcX1z1+ZexXx/cQ/k2KwtNm/bYYcdkin38ba4DVX+6DYA9IaADEBNiAEujkbG8BbDUyGxuVUcRY5Thstp0003Dffff3+yh3FcCxvX88bwGEN6d1O+d9lll2S0+Jvf/GY49dRTk8ZT8WvEqc5PPvlkvzxHMeKU6rhPclwzHNf4dvXv//7vySjv5Zdfnozix0Afg2nsMN2deH5xtPikk04K3/72tzv2ZX766afXC8jFfH9xv+U4rTq+qRDvE9cux595vG98gwEA+qIp7vXUp68AAPS72Djr8ccfT4IkAFAe1iADQIXFLZTyxVD8xz/+MVl7DQCUjxFkAKiw2JH5S1/6UrJVVezuHNfoxoZTDz/8cHj3u99d6dMDgIZhDTIAVNj+++8frrzyyqTzdtyXOO4zHPchFo4BoLyMIAMAAIA1yAAAAJAlIAMAAICADAAAABVq0rVu3brw8ssvhxEjRoSmpqZyPz0AAAANJpPJhOXLl4cJEyaEQYMGVU9AjuF44sSJ5X5aAAAAGtxLL70Utthii+oJyHHkOHdiLS0todLvIrS1tYXW1laj2XRQFxSiNkijLkijLkijLkijLsqjvb09GajN5dGqCci5H3oMx9UQkOMlnodiJEddUIjaII26II26II26II26KK+eXmNNugAAAEBABgAAgCwBGQAAAARkAAAAyBKQAQAAQEAGAACALAEZAAAABGQAAADIEpABAABAQAYAAIAsARkAAAAEZAAAAMgSkAEAACCEMKTSJ1CN1q7LhPufXxoWL18VNhkxLOw2eXQYPKip0qcFAADAABKQu7jpsYVh1g3zw8K2VR3HNmsdFk4/aErYf7vNKnpuAAAADBxTrLuE42PnzO0UjqNFbauS4/F2AAAA6pOAnDetOo4cZ1Juyx2Lt8f7AQAA0OABee3ateHUU08NkydPDhtuuGHYZpttwplnnhkymdoPjXHNcdeR43zxO4y3x/sBAADQ4GuQzznnnHDhhReGX/7yl2Hq1KnhwQcfDEceeWRobW0NJ5xwQqhlsSFXf94PAACAOg7Id999dzjkkEPCgQcemFzfaqutwpVXXhnuv//+UOtit+r+vB8AAAB1PMV6zz33DLfeemv461//mlx/5JFHwl133RUOOOCAUOviVk6xW3WhzZzi8Xh7vB8AAAANPoJ88sknh/b29rDtttuGwYMHJ2uSzzrrrDB9+vSCj1m9enVyyYmPj+K65UqvXc6dQ7wMamoKp33ifWHm5XOTMJx/ZrnQHG+P2yFX+rwpX11APrVBGnVBGnVBGnVBGnVRHsW+viUF5N/85jfh8ssvD1dccUWyBnnevHnhxBNPDBMmTAgzZsxIfczs2bPDrFmz1jve1tZW8SKIz79ixYrk86ampjBt4kbh/MO2DRff+Vx4bcWajvuN3XiDcMxeWye3x/OmvnWtC8hRG6RRF6RRF6RRF6RRF+WRG6jtSVOmhJQ6ceLEZBR55syZHcd+8IMfhDlz5oQnn3yy6BHk+HWWLVsWWlpaQiXFbz0G3thkLL8Y41ZOsVv1q8tXhXEjstOqB8ehYxpCoboAtUEadUEadUEadUEadVEeMYeOHDkyea27y6EljSC/8cYbYdCgzsuW41TrdevWFXxMc3Nzcukq/vCroQBy55F/LkMGN4U93zW2oudF9dUFRGqDNOqCNOqCNOqCNOpi4BX72pYUkA866KBkzfGWW26ZTLF++OGHw09+8pNw1FFH9fY8AQAAoCqUFJDPO++8cOqpp4bjjjsuLF68OFl7/JWvfCWcdtppA3eGAAAAUG0BecSIEeGnP/1pcgEAAICG3QcZAAAA6pWADAAAAAIyAAAAZAnIAAAAICADAABAloAMAAAAAjIAAABkCcgAAAAgIAMAAECWgAwAAAACMgAAAGQJyAAAACAgAwAAQJaADAAAAAIyAAAAZAnIAAAAICADAABAloAMAAAAAjIAAABkCcgAAAAgIAMAAECWgAwAAAACMgAAAGQJyAAAACAgAwAAQJaADAAAAAIyAAAAZBlBBgAAAAEZAAAAsgRkAAAAEJABAAAgS0AGAAAAARkAAACyBGQAAAAQkAEAACBLQAYAAAABGQAAALIEZAAAABCQAQAAIEtABgAAAAEZAAAAsgRkAAAAEJABAACgFwF5q622Ck1NTetdZs6cWcqXAQAAgKozpJQ7P/DAA2Ht2rUd1x977LGw7777hsMOO2wgzg0AAACqMyCPGzeu0/Wzzz47bLPNNmHvvffu7/MCAACA6g3I+dasWRPmzJkTTjrppGSadSGrV69OLjnt7e3Jx0wmk1wqKXcOlT4Pqou6oBC1QRp1QRp1QRp1QRp1UR7Fvr69Dsi///3vw7Jly8KXvvSlbu83e/bsMGvWrPWOt7W1VbwI4vOvWLEi+by7kE9jURcUojZIoy5Ioy5Ioy5Ioy7KIzdQ25OmTC9T6n777Rc22GCDcMMNN3R7v7QR5IkTJybhuqWlJVRS/NZjUG9tbVWMdFAXFKI2SKMuSKMuSKMuSKMuyiPm0JEjRyavdXc5tFcjyC+88EK45ZZbwjXXXNPjfZubm5NLV7kO2JWW340bctQFhagN0qgL0qgL0qgL0qiLgVfsa9urfZAvvfTSsMkmm4QDDzywNw8HAACAqlNyQF63bl0SkGfMmBGGDOn1EmYAAACo7YAcp1a/+OKL4aijjhqYMwIAAIAKKHkI+OMf/3jFu08DAABAf+vVGmQAAACoNwIyAAAACMgAAACQJSADAACAgAwAAABZAjIAAAAIyAAAAJAlIAMAAICADAAAAFkCMgAAAAjIAAAAkCUgAwAAgIAMAAAAWQIyAAAACMgAAACQJSADAACAgAwAAABZAjIAAAAIyAAAAJAlIAMAAICADAAAAFkCMgAAAAjIAAAAkCUgAwAAgIAMAAAAWQIyAAAACMgAAACQJSADAACAgAwAAABZAjIAAAAIyAAAAJAlIAMAAICADAAAAFkCMgAAAAjIAAAAkCUgAwAAgIAMAAAAWQIyAAAACMgAAACQJSADAACAgAwAAABZAjIAAAD0JiAvWLAgHHHEEWHMmDFhww03DO9///vDgw8+ODBnBwAAAGUypJQ7v/766+GDH/xg+MhHPhJuvPHGMG7cuPD000+HUaNGDdwZAgAAQLUF5HPOOSdMnDgxXHrppR3HJk+ePBDnBQAAANUbkK+//vqw3377hcMOOyzccccdYfPNNw/HHXdc+PKXv1zwMatXr04uOe3t7cnHTCaTXCopdw6VPg+qi7qgELVBGnVBGnVBGnVBGnVRHsW+viUF5Oeeey5ceOGF4aSTTgrf/e53wwMPPBBOOOGEsMEGG4QZM2akPmb27Nlh1qxZ6x1va2ureBHE51+xYkXyeVNTU0XPheqhLihEbZBGXZBGXZBGXZBGXZRHbqC2J02ZElJqDMK77rpruPvuuzuOxYAcg/I999xT9AhynKa9bNmy0NLSEiopfusxqLe2tipGOqgLClEbpFEXpFEXpFEXpFEX5RFz6MiRI5PXurscWtII8mabbRamTJnS6dj73ve+8Lvf/a7gY5qbm5NLV/GHXw0FkDuPajgXqoe6oBC1QRp1QRp1QRp1QRp1MfCKfW1L2uYpdrB+6qmnOh3761//GiZNmlTa2QEAAECVKSkgf+Mb3wj33ntv+OEPfxieeeaZcMUVV4SLL744zJw5c+DOEAAAAKotIH/gAx8I1157bbjyyivDdtttF84888zw05/+NEyfPn3gzhAAAADKoKQ1yNEnPvGJ5AIAAAANO4IMAAAA9UpABgAAAAEZAAAAsgRkAAAAEJABAAAgS0AGAAAAARkAAACyBGQAAAAQkAEAACBLQAYAAAABGQAAALIEZAAAABCQAQAAIEtABgAAAAEZAAAAsgRkAAAAEJABAAAgS0AGAAAAARkAAACyBGQAAAAQkAEAACBLQAYAAAABGQAAALIEZAAAABCQAQAAIEtABgAAAAEZAAAAsgRkAAAAEJABAAAgS0AGAAAAARkAAACyBGQAAAAQkAEAACBLQAYAAAABGQAAALIEZAAAABCQAQAAIEtABgAAAAEZAAAAsgRkAAAAEJABAAAgS0AGAACAUgPyGWecEZqamjpdtt1224E7OwAAACiTIaU+YOrUqeGWW275+xcYUvKXAAAAgKpTcrqNgXj8+PEDczYAAABQK2uQn3766TBhwoSw9dZbh+nTp4cXX3xxYM4MAAAAqnUEeffddw+XXXZZeO973xsWLlwYZs2aFT784Q+Hxx57LIwYMSL1MatXr04uOe3t7cnHTCaTXCopdw6VPg+qi7qgELVBGnVBGnVBGnVBGnVRHsW+viUF5AMOOKDj8+233z4JzJMmTQq/+c1vwtFHH536mNmzZydBuqu2traKF0F8/hUrViSfx4ZjEKkLClEbpFEXpFEXpFEXpFEX5ZEbqO1JnzpsjRw5MrznPe8JzzzzTMH7nHLKKeGkk07qdGITJ04Mra2toaWlJVRSLqDHc1GM5KgLClEbpFEXpFEXpFEXpFEX5VHsa9ungBzf6Xj22WfDF77whYL3aW5uTi5pJ1gNBZC/ZRXkqAsKURukURekURekURekURcDr9jXtqQmXd/85jfDHXfcEf72t7+Fu+++O3zyk58MgwcPDocffnhvzxMAAACqQkkjyP/v//2/JAwvWbIkjBs3LnzoQx8K9957b/I5AAAANExAvuqqqwbuTAAAAKCW9kEGAACAeiQgAwAAgIAMAAAAWQIyAAAACMgAAACQJSADAACAgAwAAABZAjIAAAAIyAAAAJAlIAMAAICADAAAAFkCMgAAAAjIAAAAkCUgAwAAgIAMAAAAWUPe+UiR1q7LhPufXxoWL18VNhkxLOw2eXQYPKip0qcFAABAHwnIJbjpsYVh1g3zw8K2VR3HNmsdFk4/aErYf7vNKnpuAAAA9I0p1iWE42PnzO0UjqNFbauS4/F2AAAAapeAXOS06jhynEm5LXcs3h7vBwAAQG0SkIsQ1xx3HTnOF2NxvD3eDwAAgNokIBchNuTqz/sBAABQfQTkIsRu1f15PwAAAKqPgFyEuJVT7FZdaDOneDzeHu8HAABAbRKQixD3OY5bOUVdQ3LuerzdfsgAAAC1S0AuUtzn+MIjdg7jWztPo47X43H7IAMAANS2IZU+gVoSQ/C+U8Yn3apjQ6645jhOqzZyDAAAUPsE5BLFMDxtmzGVPg0AAAD6mSnWAAAAICADAABAlinW/WTtuoy1yQAAADVMQO4HNz22MMy6YX5Y2Laq41jcFzlu/aS7NQAAQG0wxbofwvGxc+Z2CsfRorZVyfF4OwAAANVPQO7jtOo4cpxJuS13LN4e7wcAAEB1E5D7IK457jpynC/G4nh7vB8AAADVTUDug9iQqz/vBwAAQOUIyH0Qu1X35/0AAACoHAG5D+JWTrFbdaHNnOLxeHu8HwAAANVNQO6DuM9x3Mop6hqSc9fj7fZDBgAAqH4Cch/FfY4vPGLnML618zTqeD0etw8yAABAbRhS6ROoBzEE7ztlfNKtOjbkimuO47RqI8cAAAC1Q0DuJzEMT9tmTKVPAwAAgF4SkMtg7bqM0WUAAIAqJyAPsJseWxhm3TA/LGz7+17IsbN1bN5lfTIAAECdNOk6++yzQ1NTUzjxxBP774zqLBwfO2dup3AcLWpblRyPtwMAAFDjAfmBBx4IF110Udh+++3794zqaFp1HDnOpNyWOxZvj/cDAACgRgPyihUrwvTp08PPf/7zMGrUqP4/qzoQ1xx3HTnOF2NxvD3eDwAAgBpdgzxz5sxw4IEHhn322Sf84Ac/6Pa+q1evTi457e3tycdMJpNcKil3DgNxHovb3wxNqePH69+v0q8D5asLapvaII26II26II26II26KI9iX9+SA/JVV10V5s6dm0yxLsbs2bPDrFmz1jve1tZW8SKIzx9Hw6O4lro/jRq6Nmw+vLj7xdeC6jGQdUFtUxukURekURekURekURflkRuo7deA/NJLL4Wvf/3r4eabbw7Dhg0r6jGnnHJKOOmkkzqd2MSJE0Nra2toaWkJlZQL6PFc+rsY9xzREtbd8HTSkCvtbYD4bONbh4U93zfRlk9VZiDrgtqmNkijLkijLkijLkijLsqj2Ne2pID80EMPhcWLF4edd96549jatWvDnXfeGc4///xkKvXgwYM7Paa5uTm5pJ1gNRRA7jz6+1yGDG4Kpx00NelWHeWH5KZ3rsfbhwzuUyNxaqwuqH1qgzTqgjTqgjTqgjTqYuAV+9qWlM4+9rGPhUcffTTMmzev47LrrrsmDbvi513DcaOL+xxfeMTOyUhxvng9HrcPMgAAQPUoaQR5xIgRYbvttut0bPjw4WHMmDHrHScrhuB9p4xPulUvXr4qbDJiWNht8mjTqgEAAOqhizWliWF42jZjUm+L+yALzwAAAHUQkG+//fb+OZMGdNNjC8OsG+Z32i95s9Zh4fSDpph+DQAAUGY6RFUwHMcGXvnhOIpdr+PxeDsAAADlIyBXQJxWHUeO07Z/yh2Lt8f7AQAAUB4CcgXENcddR47zxVgcb4/3AwAAoDw06aqA2JCrGH9+5lXNuwAAAMpEQK6AGHiLcf5tz3Z8rnkXAADAwDLFugLiaHAMvKWMB2veBQAAMLAE5AqIU6XjaHBUbEjWvAsAAGBgCcgVEqdKX3jEzmF8a3HTrfObd9377JJwz7NLwnXzFiQfBWYAAIC+swa5wiF53ynjk27VsRnX06+sCOff9kyPj5t5xdyw7M23Oq5bnwwAANB3RpCrYLr1tG3GhEN23Dx88F1ji3pMfjiOrE8GAADoOwG5xpt3RdYnAwAA9J2AXOPNu7quT47TtQEAACidgFwjzbtGbji0qMfHtcwAAACUTpOuGmjetcmIYWFdJhOmX3Jfj4+N9wUAAKB0AnKVN+/KiWuL4/rk2JArbZVxnJIdR53jOmYAAABKZ4p1HaxPzl2Pt8f7AQAAUDoBuQ7WJ8fr8Xi8PY403/PsknDdvAXJR12tAQAAimOKdR2sT47TquPIcdwHOW71FLtZ58Rp2XFkOT4OAACAwgTkOlifHMVwfOycueutT45rluPx3AgzAAAA6UyxrgNxGnUcOU6bTJ07Fm833RoAAKAwAbkOxOnW+dOqu4qxON4e7wcAAEA6AbkOxLXI/Xk/AACARmQNch2IjbqKMXZ4c9LZumtzLwAAAATkuhCDbuxWHRtypa0yjhG4daOh4V+ufiQsatfhGgAAII0p1nUgjgLHoBt1HQ+O12NoXvbGW53CcX6H69gBGwAAoNEJyHUijgLHrZzGt3aebr1pS3MYudHQ1MfkRpvPuP7x8OdnXgvXzVuQTMHW7RoAAGhEpljXWUjed8r4pFt1bp3xukwmTL/kvoKPiVF4UfvqTvcx9RoAAGhERpDrcLr1tG3GhEN23Dz5+NqK1SV/DVOvAQCARiQg17liO1zny02wnnXDfNOtAQCAhiEgN0iH61I3c4qxeGHbqmS6NgAAQCMQkBu4w3Ux4lpmAACARiAgN3CH64Gaog0AAFCLdLFu0A7XY4c3h3+5+pHwSvuqjjXH+eJocwzUcYo2AABAIxCQG7DDdc4ZB09JulXHMJwfknNTsePU7PgYAACARmCKdQMrNPU6Xo/H4+2xi/U9zy4J181bkHzU1RoAAKhXRpAbXNep13HNcZxWHUeO4z7Icaun2M06J3bEjiPL8XEAAAD1REBmvanXUQzHcfp11/HiRW2rkuO5EWYAAIB6YYo164nTqOPIcdpk6tyxeLvp1gAAQD0xgsx64nTr/GnVXcVYHG+/99klYdCgpvWmZgMAANQiAZn1xMBbjJlXzA3L3nyr47r1yQAAQC0zxZr1xNHgYuSH4/z1yXH9MgAAQK0RkFlPnCodR4NLnSxtfTIAANAwAfnCCy8M22+/fWhpaUku06ZNCzfeeOPAnR0VEdcRx6nSUW9CclyfHNcxAwAA1G1A3mKLLcLZZ58dHnroofDggw+Gj370o+GQQw4Jjz/++MCdIRUR1xHHrZzGt3aebj1yw6FFPf7Pz7warpu3INzz7BKjyQAAQP016TrooIM6XT/rrLOSUeV77703TJ06tb/PjSoIyftOGZ+MBuc6Va/LZML0S+7r8bHn3/Zsx+eadwEAAHXdxXrt2rXh6quvDitXrkymWheyevXq5JLT3t6efMxkMsmlknLnUOnzqGZx16Y9th7dcT2OBk9obU4achX7qr3S9mY4bs5D4YLpO9dESFYXFKI2SKMuSKMuSKMuSKMuyqPY17fkgPzoo48mgXjVqlVh4403Dtdee22YMiW7XjXN7Nmzw6xZs9Y73tbWVvEiiM+/YsWK5POmJvv3FutfP75VmP3HJ5PPi/0Jxlf3olseC7ttvmHV75WsLihEbZBGXZBGXZBGXZBGXZRHbqC2J02ZElPqmjVrwosvvpgE3N/+9rfhkksuCXfccUfBkJw2gjxx4sSwbNmypNFXJcVvPX4fra2tirFEcSunM/8wP2nIVYpfHbV7GDSoKby6fFUYN2JY0jG72gKzuqAQtUEadUEadUEadUEadVEeMYeOHDkyea27y6ElB+Su9tlnn7DNNtuEiy66qOgTiz/8nk6sHBRj38Tp1rn1yU+/siKcf9szPT4mNvnK3z85f31y/tfbpILhWV1QiNogjbogjbogjbogjbooj2JzaK/XIOesW7eu0wgxjSOG12nbjEk+j92qiwnI+eE4imuZj50zNxyz1+Rw/SMLO41Ia+4FAABU7TZPp5xySrjzzjvD3/72t2Qtcrx+++23h+nTpw/cGVIT4mhvDLS92Tc5Xi668/n1pmvnwnOczg0AAFBVAXnx4sXhi1/8Ynjve98bPvaxj4UHHngg/O///m/Yd999B+4MqZnR5DjaG/XXhKHc3P9ZN8wPa95el4xS21sZAAAYKCVNsf6v//qvATsRal+cCn3hETsngTZ/NLjruuNSxBgcv9Yes28NS1eu6Thu+jUAANDf+rwGGfLFwLrvlPGdmm2ty2TC9Evu69PXzQ/H+dOvYyAXkgEAgP4gIDOgzbuiOB06jvjGUNtfE6Mz70zljqPVMZBX21ZRAABAna9BhmpZn5w//TqOVgMAAPSVgExZ1yePbx3W6XgcWf7KXpOT4Nzb8ByncgMAAPSVKdZUdH1y3B4qjjDvtOWo9Zp7jR4+NCxd2XNzr9eWr066W+d/PQAAgFIJyFR0fXJ34XmXSaPC3ufe1u3a5ZiFz/yfJzqu624NAAD0linWVF14PmTHzZOPGwwZ1OPa5a7bIee6W9/02MKkOZi9kwEAgGIZQaYm91aOI8dpeTfX3frkax4NZ1w/Pyxq//tjjC4DAADdEZCpel2nX8c1x/nTqtNC8rI34trlzuuX7Z0MAAB0R0Cm5tYuxynTfd07+aPbbhoeeuH19ZqFAQAAjUtApubEQNvXvZP3mH1rWLpyTcdx068BAABNuqg5cbQ3Btq+jPfmh+Ouzb0AAIDGJCBTc+JU6J66W5cq1+8rTr/W7RoAABqTgExNd7ce39p5uvX4luYwcqOhvQrOuenXv/zz8+GOpxbbGgoAABqMNcjUTXfrXLOtm+cvSqZLx5Dcm3j7gz8+ETYfHsKClSGMb93Q2mQAAGgQRpCpi+7Wh+y4efIxXi80ujx6+NCSv37+2uQ4mhxHlWMXbaPLAABQf4wg0zCjy7tMGhX2Pve2JPRmStwa6uRrHg1nXD8/LGpf1XGbztcAAFBfjCDTMKPLGwwZ1KvmXjEkL3vjrU7hONL5GgAA6ouATEMpNP26r52v17y9zvRrAACocaZYExp9+vVry1eHM//niV59rVzn6z1m39ppb+Xc9Ou0JmJxZBsAAKg+AjINPf06iqO9l9z1fElrk7vKD8dR/FpfnTM32XIqTs/OsW4ZAACqlynWNLwYlnuzNrk7uaCdH44j65YBAKB6CcjQzdrk8S3NyShwfwdn65YBAKD6mGINXdYm3/fckvDqkqVh3JjRYfetx4Sb5y9KRn1jSO6PCNvTumXTrwEAoDKMIEPK2uS937tJ8jFeLzS6PHr40D49V9q6ZdOvAQCgcowgQy86X8eO1LtMGhX2Pve2PjX3yhe/RtM706/jc+l2DQAA5WUEGUocXT5kx82TjxsMGTQgzb3i9OsYxAEAgPISkKEPCk2/jo29+hKc4yg1AABQXqZYwwBMv95t8uikuVecLh1HhPPXLS9d2XnrpzRjhzcnna3zv54p1wAAMLAEZOjH6dd9XbccI3DrRkPDv1z9SFjU/vdgrcM1AAAMPFOsoUrWLee2kVr2xludwnGkwzUAAAw8ARmqZN3ypi3NHWuXu8qNNscp22veXpdMv75u3oLk49p1/dFDGwAAMMUaKiBt+vW6TCZMv+S+Hjtc7zH71k57KJt+DQAA/cMIMlTJ9OvXVqwu6nH54Tgy/RoAAPqHEWSoEnEUuTcy76xfPuP6x8OIYUOToK3zNQAAlE5AhioRA22cLl2ow3V34v0Xta/uNEXb1GsAACiNKdZQJeJob6EO172RP/U6NvLS2AsAALpnBBmqsMN17FYdG3LljB4+NCxd+Vavpl6ffM2j4Yzr5xfcVzmG5fxmYaZmAwDQqARkqDJpHa53mTQq7H3ubSVPv87tqxzCW6mjy8fsNTlc/8jCTmHc1GwAABqVKdZQAx2uNxgyqF+nX2feuVx05/OdwnGkKzYAAI1KQIYam349vrV33a6LlRuhjtO8rVUGAKCRmGINNTz9euzw5vAvVz8SXmkvvfN1d+LXiiPL8XniCDYAADSCkkaQZ8+eHT7wgQ+EESNGhE022SQceuih4amnnhq4swO6nX79wXePDWcc3H9Tr7v68zOv6nwNAEDDKCkg33HHHWHmzJnh3nvvDTfffHN46623wsc//vGwcuXKgTtDoFdTr8e3NIeRGw3tU3A+/7Znw9evmhcO//m94UPn/KljXbJtowAACI0+xfqmm27qdP2yyy5LRpIfeuihsNdee/X3uQF96Hwdt2u6ef6ipOFWDMl9jbDFdr62bRQAAA25BrmtrS35OHr06IL3Wb16dXLJaW9vTz5mMpnkUkm5c6j0eVBdarUuYgbdY+vO/y/uN3V8+M/pO4Uz/zB/vUB70PabhZ//f88n10v5Ti++87nkY37kfaXtzXDcnIfClz88Odzwl/XD86mfqI9to2q1NhhY6oI06oI06oI06qI8in19mzK9/EmsW7cuHHzwwWHZsmXhrrvuKni/M844I8yaNWu94y+88EJoaWkJlRS/9RUrVoSNN944NDUZ4aJ+6yKO6j7+cnt4feXqMGp4c5g6oSUZ1b372deSwPvaijUD9ty5V/CUf9w27LnN2FDL6rE26Dt1QRp1QRp1QRp1UR5xoHbSpEnJIG93ObTXAfnYY48NN954YxKOt9hii5JGkCdOnJgE62oIyPEFam1tVYw0bF3kpkS/unxVePqVFeGC25/p9+eIr2JcI33bNz8SHnrh9eS5xtXg9OtGqw2Koy5Ioy5Ioy5Ioy7KI+bQkSNH9hiQezXF+vjjjw9/+MMfwp133tltOI6am5uTS1fxh18NBZA7j2o4F6pHI9XFkMFNYc93ZUd2Y8Ot829/tt+fI74L93Lb6jDt7D+FpSvXpK5drhWNVBsUT12QRl2QRl2QRl0MvGJf2yGlvrvxta99LVx77bXh9ttvD5MnT+7t+QFVKI7oxtAaG3INxCqY/HCc3/jrgs/vlEz91tgLAIBKKikgxy2errjiinDdddcleyEvWrQoOR6nA2y44YYDdY5AmcRQGkd0+6vzdU9yX//4Kx8O+TtF6YoNAEAllLQGudCw9KWXXhq+9KUvFT33OwbqnuZ+l4P5/qRRFyHZ73jWDet3vj54h83CxXeW3vm6VLlXvactpcpNbZBGXZBGXZBGXZBGXZRHsTm05CnWQOPuqxxHbnfactSAh+fc4y965+ulTcu+8Iida2rtMgAAdb4PMlC/Yhiets2YPofn0cOHhqUr3+q384rhOb63Gp/no9tumnTFNv0aAID+ICADAxqed5k0Kux97m392vgrfp0YwveYfWvNd8UGAKB6DKr0CQD1GZ4P2XHz5OMGQwYloTVqKlNX7LiGGgAASiUgAwMujujGNcPjW4d1Ot7fs6FzI9Rx+nXsfg0AAKUwxRooi7Tp16+vXBNmXjE3ub2/p19f9ufnw9gRzdYmAwBQNAEZqOja5QsH7TwgXbHP/J8nOn09a5MBAOiJgAxUVDm6YudvDVXouQAAQEAG6r4rdm5rqJOveTSccf38sKh9Verocly3LDwDADQuARmoufAcA20cEY7RtZSQvOyNOPL8Vuro8jF7TQ7XP7JwvanepmYDADQOXayBuumK3RuZdy4X3fl8p3Ac2TYKAKCxGEEGalLX6devLV/dqTFXf8hNzY7roD+67abhwb8tDa8uWRrGjXk77L71mGR027RsAID6ISADdTH9OgbVS+56vqS1yaVsG7XH7FvD6ytXh82Hh7BgZQjjWzdMOm13Ny1beAYAqC0CMlAXYvDszdrkYi1duSb5ujkxFMdp2V1Z0wwAULusQQbqfm3y+JbmMHKjoZ0C7kApdk1zHF2+59kl4bp5C5KP8ToAAJVlBBloiH2Vb56/aMBGl0M/bzcFAEBlGEEG6nZt8iE7bp58jNcLjS7HYPqVvSYnwbVcI8xxu6n8cJw/uvzHv7xsZBkAoEKMIAOh0UeXY4DeactRSbfq/GnRo4cPDUtXdt43eaDkYvDxVz4c8jOxpl8AAOUjIAMN2/m6p/C8y6RRYe9zb+v3ztjd6TpgrOkXAED5mGINUGBq9gZDBiUBNKrUOK2mXwAA5WMEGaAbubXLcfr1orY3O43eFtoHOR6/+J0toDJV0PTL1GwAgOIIyAA9yE2/vu+5JeHVJUvDuDGjw+5bZ5t/fXv/9xW9pjluN7Xq7XWh7Y23+i0455p+hfBWyVOzC63HBgBoVE2ZTKas8/Da29tDa2traGtrCy0tLaGS4rcezyOeT1OTPwrJUhf0V22kjdzmtptKvl6ojNxWV3Fv6Gy4zrKmuXf8ziCNuiCNuiCNuqiuHGoNMkAVbDdVzoHbXDDPD8eRNc0AQKMzxRqgzNI6Zr++ck2YeUVlR5aLXdMMAFCvjCADVMHo8j9unz6yHIPpV/aanATXcgwy59Y054fjyOgyANAIjCADVPHIcrmbfhVidBkAaAQCMkAVjiwXG55zTb9yjbcGUk8dsy/4/E5h1PBmXbEBgJolIAPUcHjO36e5kvsxR8df+XDIn21tL2YAoNYIyAANNjU7t71Tf486d12KXMxezKZlAwDVREAGaNCp2QO9pjn3NS56ZxQ7LTzH0e9C4R4AoNwEZIAGnZpdyTXNxTT96i44m7INAAwEARmgQZWypnkgOmZ31/Trq3PmdkwF7xqco7Q116ZsAwB9JSADUHWjy1F+OM4PzmlM2QYA+oOADECfRpdj9uzaoGsgdPcU9mkGAPqDgAxAn0aXX1+5Jsy8IjuyW4ac3Ot9mrsbXbamGQCIBGQA+jy6fOGgyu7FHPowuhzPsbttqIRnAGgcAjIAZd+LeSCafvVmdDmeU3fbUPW0h3MMz/c9tyS8umRpGDfm7bD71mOEZwCoYQIyABXbi7kcTb8Gag/nXHhe1PZm2Hx4CAtWhjC+dUPrnQGghgnIAFRV06/c9k5dg3P+9UqG6q7huSklPF/w+Z3CqOHNpmUDQI0RkAGomO5Gl9eblt3NPsjlnrJdSO65j7/y4U6dva1pBoDaICADUJWjy93tZ1ztU7a7bntV7JpmAKCyBGQAampNc6lTtqthdLmYNc3dbUMFAFRpQL7zzjvDueeeGx566KGwcOHCcO2114ZDDz10YM4OAEpQqw3ButuGytRsAKjigLxy5cqwww47hKOOOip86lOfGpizAoBeKmV0ubt9kMu5h3OhbaiKmZpt1BkAKhiQDzjggOQCALWku3XN397/fUXv4dw1PFdyavZX58zt6Pqdf35GnQGgd6xBBiA0+rrmUvdwzg/PcR/knJg9uzboGii5p8kPx5GGYABQxQF59erVySWnvb09+ZjJZJJLJeXOodLnQXVRFxSiNhpTDL17bD2607FYA/tNHR/2ed+m4b7nloTXliwNY8eMTsLq166cm71PqKyL73wu+Zg/XvxK25vhuDkPhQum/70h2KvLV4VxRpf7nd8XpFEXpFEX5VHs6zvgAXn27Nlh1qxZ6x1va2ureBHE51+xYkXyeVOTPwrIUhcUojZIM2XskLBi2EZh442HhKamoeH8w7ZNwulrK9Z03GfsxhuEvd8zNlwz9+XkeiX/9YuV+x83/SVcdPNj4bWVnc/xmL22DrtPHhMef7k9vL5ydRg1vDlMndAiOPeC3xekURekURflkRuo7UlTpg8pNf4Ae+pinTaCPHHixLBs2bLQ0tISKil+6zGot7a2KkY6qAsKURsUWxe5tb9dR2dvemxhOPMP1bcNVZTr8J22pvnUT3Re02zUuWd+X5BGXZBGXZRHzKEjR45MXuvucuiAjyA3Nzcnl67iD78aCiB3HtVwLlQPdUEhaoNi6mLI4Kaw57vGrne/A94/IXx86mYFt6GKKr1X8+tvvN1pYvbCttXhuMsfDsfstazbNc0agq3P7wvSqAvSqIuBV+xrW3JAjsP/zzzzTMf1559/PsybNy+MHj06bLnllqV+OQBoKL3Zhqpc2031ppN2MQ3BCoVnoRqAalNyQH7wwQfDRz7ykY7rJ510UvJxxowZ4bLLLuvfswOABlFMx+z8AJqbCp2bGh1qLDwX2n9al20AKqlPa5B7O/c7zq/vae53OZjvTxp1QSFqg0rWRdpoa5yaXY2jzr2Ve/XiaHqhNwuiWhh59vuCNOqCNOqiunKofZABoIanZpcy6lwtDcEKiecU/zQ8+ZpHwxnXzw+L2tcfXY7S3hSIt3UXqgGgGAIyANRZcO4uPOcaglVyanZ34jllu2j/vZN2bsr2V99pZNZV7ra0DtyaiAFQCgEZAOpUKQ3BqmFNc3cyRdyWH477o4kYAI1HQAaABtPd6HI9rWkeyA7cANQnARkAGlB/rGlu5PBsvTNAfRKQAYBer2luxPBczHrn+55bEl5dsjSMG/N22H3rMcIzQI0QkAGAioTnrqOz1d5lu5T1zova3gybDw9hwcoQxrduaMo2QI2wD7I9x+hCXVCI2iCNuuhZoVBYaG/nGDSj/D9Q8puHVWsjsXxNIdMRkHM7PPd2yrZQXT/8viCNuigP+yADAFU98lxKl+3x3eyDXCsduHszZbu7fZ/jawVA/zKC7N0aulAXFKI2SKMuBkZ3o6aFRp6raR10/ghy5p0R5NIeX/h8c18tvpFg5Lm2+H1BGnVRHkaQAYC6G3UudFu9NRHrad/n+Cf0ydc8Gs64fn5Y1L7+6HJk5BmgdAIyAFAXGqkDdzyn7JTst1KnbKfJNRK74PM7hVHDm0sedTYiDTQCARkAqHv9FZ5rZb1zd7cdf+XDYV3eHYsZde7uNntCA/XEGmTz/elCXVCI2iCNuqhvvV3vHPq4Brmcelrv3NNtPe0JLTz/nd8XpFEX5WENMgBAH/V2vXPcB7lWpmxn+nhbT3tCWwcN1BIBGQCgn6ds3/fckvDqkqVh3JjRYfetx5Q8ZbvW9n0uZVurGJ6768Bt1BmoJAEZAGAAwnPb2CGdpkwWGnlOm7Ld3b7P41uaw6q314W2N96qqeBcTAfuONLe3ahzb5uICd1AsaxBNt+fLtQFhagN0qgL+qMuerPvcxyJTZ6rTkaeC8m9et1N2R6IBmMDEar9viCNuqiuHCogK0a6UBcUojZIoy6oVF3c9NjCkkJhzHb53avrwUA1GOtpH+lC4bmnUO33BWnURXkIyEVQjKRRFxSiNkijLqhkXZQy8vz6yjVh5hWljTrX44h0d3oK1d2NZBczPTxtbXpkCnhj8+9IeQjIRVCMpFEXFKI2SKMuqKW6KHXUubvbqn1P6GrQNVTH7ua57b/Gt244oKPV1I5q/X1RbwTkIihG0qgLClEbpFEX1Fpd9Gejq2L2hBae/64pb3/s7LVC9+v7aLXwXDuq+fdFPRGQi6AYSaMuKERtkEZd0Oh1USiMpY1W13IH7v4OyJmOGNxfXzurp/2njUhXn0b6fVFJxeZQ2zwBADBge0IX6sBtanb5958eiBFp22tRb4wge7eGLtQFhagN0qgL0qiL3q2FLhTUCk3ZrrUGYwM5gtzf6ml7rWrn90V5mGJdBMVIGnVBIWqDNOqCNOqiZ6VO9S1Hg7GBDtW1FJDrZXutWgjcfl+Uh4BcBMVIGnVBIWqDNOqCNOqidhuMdRfU+tp8rF4Ccq1sr9VT4K4Wfl+Uh4BcBMVIGnVBIWqDNOqCNOqifgN3X6aHhwJdrKtlCngjBe4Lj9i5V1O9B2JE2u+L8hCQi6AYSaMuKERtkEZdkEZd1LfeTg8vZR9kW2UNnPh/ZOtGQ8OwIYPDovbSZhD0dzOzePy+55aEV5csDePGjA67bz2m6qaA1wsBuQj+8SKNuqAQtUEadUEadUFX3QWhgRitFqpL15uR/L40M8v9HNPeOOlL93DSCchF8I8XadQFhagN0qgL0qgL+qsu+quZmfBcvYE7f216LnL3tnt4d9PGG127gNwz/3iRRl1QiNogjbogjbqgknUx0CPStba9VrUrtnlbXzqEG5EORefQIWU9KwAAYEDFYDNtmzHrHY8hqdAI47f3f1/q8Z22HLVeqB7fzUjm+CrcXqteZIq4LT8cR4vaVoVj58zt9Yj0/r3cXquWGUH27i5dqAsKURukURekURfUU13U2vZatRa4K7X910Btr7V/FW2hlc8U6zr+JcXAUhcUojZIoy5Ioy5Ioy7Kt71WwRHuluaw6u11oe2Nt6omONfL/thNeVtoVWNINsUaAAComSng/T09PEp7TBzFjtOOTfXuX5l3Xrf4pkR83Wt1urWADAAAVL1C4bnUwB3DdhzlLGX9dH83M6tXmRCS1ye+KVHoZ1LtBGQAAKChdDciHQ10M7Ou+yDXW+BevPzv32utEZABAICG05up3r2ZAt5T4L7vuSXh1SVLw7gxo8PuW4/pVeAu1CG8UjYZMSzUKgEZAACggoG7beyQTs3behO4C3UIL+eIdNM7QT6eT60SkAEAAKpMb9ZWFwrVvZ0CfnGJ22tF8WvWaoOuSEAGAACoA/09BXynlGBdqGFZLnBX4xZPpRCQAQAA6lw5t9eqZYN686ALLrggbLXVVmHYsGFh9913D/fff3//nxkAAAAVNfid8HzIjpsnH3MhuNDxhgvIv/71r8NJJ50UTj/99DB37tywww47hP322y8sXrx4YM4QAAAAqjEg/+QnPwlf/vKXw5FHHhmmTJkSfvazn4WNNtoo/OIXvxiYMwQAAIBqW4O8Zs2a8NBDD4VTTjml49igQYPCPvvsE+65557Ux6xevTq55LS3tycfM5lMcqmk3DlU+jyoLuqCQtQGadQFadQFadQFadRFeRT7+pYUkF977bWwdu3asOmmm3Y6Hq8/+eSTqY+ZPXt2mDVr1nrH29raKl4E8flXrFiRfJ7bcwzUBYWoDdKoC9KoC9KoC9Koi/LIDdRWvIt1HG2Oa5bzT2zixInJRtgtLS2hknIBPX9TblAXFKI2SKMuSKMuSKMuSKMuyqPY17akgDx27NgwePDg8Morr3Q6Hq+PHz8+9THNzc3JJe0Eq6EAcudRDedC9VAXFKI2SKMuSKMuSKMuSKMuBl6xr21JTbo22GCDsMsuu4Rbb72149i6deuS69OmTSv9LAEAAKBKlDzFOk6XnjFjRth1113DbrvtFn7605+GlStXJl2tAQAAoGEC8mc/+9nw6quvhtNOOy0sWrQo7LjjjuGmm25ar3EXAAAA1JJeNek6/vjjkwsAAADUi5LWIAMAAEC9EpABAABAQAYAAIA+rEHuj42w29vby/3UqecSz8OeY+RTFxSiNkijLkijLkijLkijLsojlz9zebRqAvLy5cuTjxMnTiz3UwMAANDAli9fHlpbWwve3pTpKUL3s3Xr1oWXX345jBgxouLvkMR3EWJQf+mll0JLS0tFz4XqoS4oRG2QRl2QRl2QRl2QRl2UR4y9MRxPmDAhDBo0qHpGkOPJbLHFFqGaxEJUjHSlLihEbZBGXZBGXZBGXZBGXQy87kaOczTpAgAAAAEZAAAAsho6IDc3N4fTTz89+Qg56oJC1AZp1AVp1AVp1AVp1EV1KXuTLgAAAKhGDT2CDAAAADkCMgAAAAjIAAAAkCUgAwAAQKMH5AsuuCBstdVWYdiwYWH33XcP999/f6VPiTKaPXt2+MAHPhBGjBgRNtlkk3DooYeGp556qtN9Vq1aFWbOnBnGjBkTNt544/DpT386vPLKKxU7Z8rr7LPPDk1NTeHEE0/sOKYmGteCBQvCEUcckfzsN9xww/D+978/PPjggx23x56Xp512Wthss82S2/fZZ5/w9NNPV/ScGVhr164Np556apg8eXLyM99mm23CmWeemdRCjrpoDHfeeWc46KCDwoQJE5J/N37/+993ur2YOli6dGmYPn16aGlpCSNHjgxHH310WLFiRZm/E8pVF2+99Vb4zne+k/xbMnz48OQ+X/ziF8PLL7/c6Wuoi/Jr2ID861//Opx00klJS/W5c+eGHXbYIey3335h8eLFlT41yuSOO+5Igs69994bbr755uQX1cc//vGwcuXKjvt84xvfCDfccEO4+uqrk/vHX1qf+tSnKnrelMcDDzwQLrroorD99tt3Oq4mGtPrr78ePvjBD4ahQ4eGG2+8McyfPz/8+Mc/DqNGjeq4z49+9KPwH//xH+FnP/tZuO+++5I/eOK/K/FNFerTOeecEy688MJw/vnnhyeeeCK5HuvgvPPO67iPumgM8W+H+LdkHHxJU0wdxBD0+OOPJ3+T/OEPf0jC1THHHFPG74Jy1sUbb7yRZJD4Jlv8eM011yQDNQcffHCn+6mLCsg0qN122y0zc+bMjutr167NTJgwITN79uyKnheVs3jx4viWf+aOO+5Iri9btiwzdOjQzNVXX91xnyeeeCK5zz333FPBM2WgLV++PPPud787c/PNN2f23nvvzNe//vXkuJpoXN/5zncyH/rQhwrevm7dusz48eMz5557bsexWC/Nzc2ZK6+8skxnSbkdeOCBmaOOOqrTsU996lOZ6dOnJ5+ri8YU/0249tprO64XUwfz589PHvfAAw903OfGG2/MNDU1ZRYsWFDm74By1EWa+++/P7nfCy+8kFxXF5XRkCPIa9asCQ899FAyvSVn0KBByfV77rmnoudG5bS1tSUfR48enXyMNRJHlfPrZNtttw1bbrmlOqlzcWbBgQce2OlnH6mJxnX99deHXXfdNRx22GHJkoyddtop/PznP++4/fnnnw+LFi3qVButra3J8h21Ub/23HPPcOutt4a//vWvyfVHHnkk3HXXXeGAAw5IrqsLiq2D+DFOn42/Z3Li/ePfp3HEmcb5WzROxY61EKmLyhgSGtBrr72WrBvadNNNOx2P15988smKnReVs27dumSdaZxCud122yXH4j9mG2ywQccvqfw6ibdRn6666qpkqlOcYt2Vmmhczz33XDKVNi7N+e53v5vUxwknnJDUw4wZMzp+/mn/rqiN+nXyySeH9vb25I2ywYMHJ39bnHXWWcmUyEhdUGwdxI/xzbd8Q4YMSd60VyuNIU63j2uSDz/88GS9caQuKqMhAzKkjRg+9thjyTv/NK6XXnopfP3rX0/W+cTmfZD/Jlp8B/+HP/xhcj2OIMffGXE9YQzINKbf/OY34fLLLw9XXHFFmDp1apg3b17yZmtstqMugGLF2Wmf+cxnkmZu8c1YKqshp1iPHTs2eae3a+fZeH38+PEVOy8q4/jjj0+aHtx2221hiy226DgeayFOx1+2bFmn+6uT+hWnUMdGfTvvvHPyDm28xEZcsbFK/Dy+268mGlPsPDtlypROx973vveFF198Mfk89/P370pj+da3vpWMIn/uc59LOtF+4QtfSBr5xV0SInVBsXUQP3ZtFPv2228nHYzVSmOE4xdeeCF5gz43ehypi8poyIAcp8Ttsssuybqh/NGBeH3atGkVPTfKJ75LF8PxtddeG/70pz8l23TkizUSO9bm10nsLhj/IFYn9eljH/tYePTRR5NRoNwljhrG6ZK5z9VEY4rLL7puAxfXnU6aNCn5PP7+iH+s5NdGnHob14ipjfoVu9DGtYD54hvw8W+KSF1QbB3Ej/HN1/hGbU782yTWUlyrTH2H47jl1y233JJsI5hPXVRIpkFdddVVSffAyy67LOkQd8wxx2RGjhyZWbRoUaVPjTI59thjM62trZnbb789s3Dhwo7LG2+80XGfr371q5ktt9wy86c//Snz4IMPZqZNm5ZcaBz5XawjNdGYYmfRIUOGZM4666zM008/nbn88sszG220UWbOnDkd9zn77LOTf0euu+66zF/+8pfMIYcckpk8eXLmzTffrOi5M3BmzJiR2XzzzTN/+MMfMs8//3zmmmuuyYwdOzbz7W9/u+M+6qJxdj94+OGHk0v88/onP/lJ8nmuG3ExdbD//vtndtppp8x9992Xueuuu5LdFA4//PAKflcMZF2sWbMmc/DBB2e22GKLzLx58zr9Lbp69eqOr6Euyq9hA3J03nnnJX/obrDBBsm2T/fee2+lT4kyir+o0i6XXnppx33iP1zHHXdcZtSoUckfw5/85CeTX1w0bkBWE43rhhtuyGy33XbJm6vbbrtt5uKLL+50e9zK5dRTT81suummyX0+9rGPZZ566qmKnS8Dr729Pfn9EP+WGDZsWGbrrbfOfO973+v0x626aAy33XZb6t8U8U2UYutgyZIlSfDZeOONMy0tLZkjjzwyCVjUZ13EN9UK/S0aH5ejLsqvKf6nUqPXAAAAUC0acg0yAAAAdCUgAwAAgIAMAAAAWQIyAAAACMgAAACQJSADAACAgAwAAABZAjIAAAAIyAAAAJAlIAMAAICADAAAAFkCMgAAAIEQ/n8NM6aXmwkKoAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "viz.plot_RN_eigenvalues(\n", " R=data_cov.data,\n", " N=noise_cov.data\n", ")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "a0d7c0ef36f63f57", "metadata": { "collapsed": false }, "source": [ "#### Localization across ranks\n", "\n", "We systematically evaluate beamformer performance for different rank values (1 to ``n_sources``). For each rank, we:\n", "\n", "1. Construct rank-**k** localizer.\n", "2. Compute localization metrics:\n", " - number of correctly localized sources,\n", " - distance error \n" ] }, { "cell_type": "code", "execution_count": 15, "id": "42008dea5f9ccc5d", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:16:26.669413Z", "start_time": "2026-06-27T19:16:16.754866Z" }, "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[36mCalculating activity index for localizer: mai_mvp\u001b[0m\n", "\n", "[Activity Index Result]\n", " Selected indices (index_max): [1913, 2178, 1248, 2783, 3290]\n", " Index max values: [3.847179 4.40560165 4.56822079 4.71270861 4.83707129]\n", " Rank parameter (r): 1\n", "\n", "\u001b[36mCalculating activity index for localizer: mai_mvp\u001b[0m\n", "\n", "[Activity Index Result]\n", " Selected indices (index_max): [1913, 2178, 2783, 3290, 4359]\n", " Index max values: [3.847179 4.74318399 5.28592942 5.5143203 5.78097535]\n", " Rank parameter (r): 2\n", "\n", "\u001b[36mCalculating activity index for localizer: mai_mvp\u001b[0m\n", "\n", "[Activity Index Result]\n", " Selected indices (index_max): [1913, 2178, 2783, 27, 3290]\n", " Index max values: [3.847179 4.74318399 5.53240386 6.07573153 6.36922522]\n", " Rank parameter (r): 3\n", "\n", "\u001b[36mCalculating activity index for localizer: mai_mvp\u001b[0m\n", "\n", "[Activity Index Result]\n", " Selected indices (index_max): [1913, 2178, 2783, 27, 4530]\n", " Index max values: [3.847179 4.74318399 5.53240386 6.29491264 6.84870992]\n", " Rank parameter (r): 4\n", "\n", "\u001b[36mCalculating activity index for localizer: mai_mvp\u001b[0m\n", "\n", "[Activity Index Result]\n", " Selected indices (index_max): [1913, 2178, 2783, 27, 4530]\n", " Index max values: [3.847179 4.74318399 5.53240386 6.29491264 7.04857874]\n", " Rank parameter (r): 5\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGwCAYAAACHJU4LAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAMStJREFUeJzt3Qd4VFX6x/E3oSR0BBOKQOg1gAZdBWkqZUEpS1ERBQIL6x8EQgSliJRHaSpSFhGVIizCQgBZWCWyKCC49BYEQYoSmgSWEnrJ/J/3rMkmhDJ3coe5Yb6f55lnZu7M3Dk5ZJhfznnvuQEul8slAAAADhTo6wYAAADcDkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4VlbJxJKSkuTo0aOSJ08eCQgI8HVzAACAG3QJt8TERClatKgEBgbev0FFQ0rx4sV93QwAAOCB+Ph4KVas2P0bVHQkJfkHzZs3r6+bI5l9dCohIUFCQkLumm5xe/SjfehL+9CX9qAf7XPu3Dkz0JD8PX7fBpXk6R4NKQSVjH8AL1++bPqRD6Dn6Ef70Jf2oS/tQT/az52yDXoaAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAAA4FkEFAADcv0Hlxo0bsm3bNjl9+rQ9LQIAAPA0qERFRcnUqVNTQkq9evUkIiJCihcvLitXrrS6OwCAF924kSRxq36Vnat/Ndd6H9bRj76T1eoLYmJi5OWXXza3lyxZIgcPHpSffvpJZs2aJYMGDZK1a9e6va+hQ4fKsGHD0myrUKGC2R8AIGN+WLhbPukdK6eOJkrxGnkkfnOiFCyaR7qNbyy1WlXydfMyDfoxk42onDx5UgoXLmxuf/XVV9K2bVspX768dO7cWeLi4iw3oEqVKnLs2LGUy5o1ayzvAwCQ/st1RJsYOXk4Mc32k0cSzXZ9HHdHP2bCoFKoUCHZtWuXmfZZtmyZNGzY0Gy/ePGiZMmSxXIDsmbNaoJP8uXBBx+0vA8AwP/otISOAIjrFg/+vu2TqG+YvrgL+jGTTv1ERkbK888/L0WKFJGAgABp0KCB2b5+/XqpWLGi5Qb8/PPPUrRoUQkODpaaNWvKyJEjpUSJErd87pUrV8wl2blz58x1UlKSucBz2n8ul4t+zCD60T70pee0jkKnKQJ+/1NUrwMC/nud7NSRc+Z5VeuF+aydTkc/eo+Vz7XloKJ1JeHh4RIfH2+mfYKCgsx2HU3p37+/pX09/vjjMmPGDFOXotM+Wq9Sp04d2blzp+TJkyfd8zXE3FzTohISEuTy5ctWfxTc9Etz9uxZ88UQGMhR656iH+1DX3ouIeGEqaVIpl+uD5bNYW67XGmfd+LEf7cjPfrRexIT006l2RpUVJs2bcx16nDQsWNHy/tp0qRJyu1q1aqZ4BIWFibz5s2TLl26pHv+gAEDJDo6Os2Iih5tFBISInnz5vXgJ0HqLwUdIdO+5EvBc/SjfehLz4WEXDIFn8mSRwDitySKK9UfsiEhoRIaGuqDFmYO9KP36CyK14KK1qaMGDFCPv74Y/ntt99k7969Urp0aRk8eLCULFnylgHDXfnz5zeFufv27bvl4zp6kzyCk5r+J8Z/ZBmnXwr0ZcbRj/ahLz0TXjfMHJWiBZ/JtRQ6AqBfruYLVkcGiuU1z6Nvb49+9B4r/WW5Z999910zXTNmzBjJnj17ynadDvrss88kI86fPy/79+839S8AAM9kyRJoDp01Am568Pf73cY1Ms/D7dGPzmC5d2fOnCmffPKJtG/fPs1RPtWrV7e8/knfvn1l1apV8ssvv8gPP/wgf/rTn8w+27VrZ7VZAIBUdH2PgTFt5MGH0tb76QiAbmf9D/fQj75neernyJEjUrZs2VvOJ1+7ds3Svg4fPmxCyalTp8w8dO3atWXdunXmNgAgY/RL9PEWFcxRKVrwqbUUOk3BCIA19GMmCyqVK1eW77//3hS93rxi7SOPPGJpX3PnzrX69gAAC/TLVA+d1aNStOCTWgrP0I+ZKKi8/fbb5ggfHVnRUZSFCxfKnj17zJTQ0qVLvdNKAADglyxHwhYtWphz/PzrX/+SXLlymeCye/dusy15lVoAAAA7eLSOii7Ktnz5clsaAAAAYNuIysaNG81y+TfTbZs2bbK6OwAAAPuCSo8ePczy+TfTmhV9DAAAwGdBRc+cHBERkW67HvGjjwEAAPgsqOgS9rp0/s30pIJZs3pU8gIAAGBPUGnUqJE5OaCe1TTZmTNnZODAgRz1AwAAbGV5COS9996TevXqmQXfkhd427ZtmxQqVEhmzZplb+sAAIBfsxxUihUrJjt27JDZs2fL9u3bJUeOHBIZGWmWws+WLZt3WgkAAPySpaCi5/KpWLGiWYG2W7du3msVAACA1RoVHTG5fPmy91oDAACQ0XVURo8eLdevX7f6UgAAAO/WqOjKtCtWrJBvvvlGqlatas73k5qepBAAAMAnQSV//vzSunVrW94cAADA1qAyffp0qy8BAAC4NzUqAAAAjh1RKVWqlAQEBNz28QMHDmS0TQAAAJ4FlaioqHRrq2zdulWWLVsm/fr1s7o7AAAA+4JK7969b7l90qRJsmnTJqu7AwAA8H6NSpMmTWTBggV27Q4AAMC+oBITEyMFChSwa3cAAADWp370jMmpi2ldLpccP35cEhIS5KOPPrK7fQAAwI9ZDiotW7ZMcz8wMFBCQkKkfv365oSFAAAAPgsqQ4YMse3NAQAAbA0q6saNG/Lll1/K7t27zf0qVapI8+bNJUuWLJ7sDgAAwJ6gsm/fPmnatKkcOXJEKlSoYLaNHDlSihcvLv/85z+lTJkyVncJAABgz1E/vXr1MmEkPj5etmzZYi6HDh0yK9bqYwAAAD4bUVm1apWsW7cuzaHIBQsWlFGjRsmTTz5pW8MAAAAsj6gEBQVJYmJiuu3nz5+X7Nmz29UuAAAA60Hlueeek27dusn69evNGip60RGWV1991RTUAgAA+CyoTJgwwdSo1KxZU4KDg81Fp3zKli0r48ePt61hAAAAlmtU8ufPL4sXLzZH/yQfnlypUiUTVAAAAHy+jorSYKIXXVMlLi5OTp8+LQ888ICtjQMAAP7N8tRPVFSUTJ061dzWkFKvXj2JiIgw66isXLnSG20EAAB+KtCTsyRXr17d3F6yZIkcOHBAfvrpJ+nTp48MGjTIG20EAAB+ynJQOXnypBQuXNjc/uqrr+T555+X8uXLS+fOnc0UEAAAgM+CSqFChWTXrl1m2mfZsmXSsGFDs/3ixYuc6wcAAPi2mDYyMtKMohQpUkQCAgKkQYMGZruuq1KxYkV7WwcAAPya5aAydOhQCQ8PN+f6adu2rVmpVuloSv/+/b3RRgAA4Kc8Ojy5TZs26bZ17NjRjvYAAAB4XqMCAABwrxBUAACAYxFUAACAYxFUAABA5i6mPXfunNs7zJs3b0baAwAAYC2o6BmTdc0Ud+hCcAAAAPcsqHz33Xcpt3/55RezXkqnTp2kZs2aZtu///1v+fzzz2XkyJG2NAoAAMDtoKJnSE42fPhwGTt2rLRr1y5lW/PmzaVq1aryySefsJ4KAADwXTGtjp48+uij6bbrtg0bNtjVLgB+7MaNJIlb9avsXP2rudb7APyT5ZVpixcvLp9++qmMGTMmzfbPPvvMPAYAGfHDwt3ySe9YOXU0UYrXyCPxmxOlYNE80m18Y6nVqpKvmwfA6SMqH374oUycONFM9fz5z382l2rVqplt+pinRo0aZQp2o6KiPN4HgMwfUka0iZGThxPTbD95JNFs18cB+BfLQaVp06ayd+9eadasmfznP/8xF72t2/QxT2zcuFGmTJliAg8A/6TTOzqSIq5bPPj7tk+ivmEaCPAzHp2UUKd4RowYYUsDzp8/L+3btzfTSe+8884dn3vlyhVzuXl9l6SkJHOB57T/XC4X/ZhB9KPntB5Fp3sCfv/zSa91VYTk++rUkXPmeVXrhfmsnZkRv5f2oB/tY6UPPQoq33//vRkBOXDggMyfP18eeughmTVrlpQqVUpq165taV89evSQZ599Vho0aHDXoKKHPw8bNizd9oSEBLl8+bLlnwNpf2nOnj1rPoSBgSxY7Cn60XMJCSdMTUoyDSkPls1hbrtcaZ934sR/t8M9/F7ag360T2Ji2uldW4PKggUL5JVXXjGjIFu2bEkZ4dB/PB1l+eqrr9ze19y5c80+dOrHHQMGDJDo6Og0Iyo6uhMSEsKKuDZ8ALVGSPuSD6Dn6EfPhYRcMoWzyZJHUuK3JIor1R9fISGhEhoa6oMWZl78XtqDfrRPcHCw94KKjnp8/PHH0qFDBxM0kj355JN3HRFJLT4+Xnr37i3Lly93u8FBQUHmcjP9heGXJuP0A0hfZhz96JnwumHm6B4tnE2uSdGRFA0pJqjoCEuxvOZ59K11/F7ag360h5X+s9zTe/bskbp166bbni9fPjlz5ozb+9m8ebOcOHFCIiIiJGvWrOayatUqmTBhgrnNUvyAf8mSJdAcgmzcfMaO3+93G9fIPA+A/7D8iS9cuLDs27cv3fY1a9ZI6dKl3d7PM888I3FxcbJt27aUiy4ap1NKejtLlixWmwYgk9N1UgbGtJEHH/pfrYrSkRTdzjoqgP+xPPXTtWtXM2Uzbdo0MwR29OhRs1pt3759ZfDgwW7vJ0+ePBIeHp5mW65cuaRgwYLptgPwHxpGHm9RwRzdo4WzWpOi0z2MpAD+yXJQ0RMSakGRjohcvHjRTANp3YgGlZ49e3qnlQD8ioYSPQRZj+7RwlnqAQD/ZTmo6CjKoEGDpF+/fmYKSNdBqVy5suTOnTvDjVm5cmWG9wEAAO4flv9MmTlzpuzevVuyZ89uAsof/vAHE1J0HRN9DAAAwGdBpVOnTiac6Hoqqek6KpGRkbY1DAAAwKOJX10dVhd9Gzp0qP0tAgAAyEhQefnll+Xbb781y+i3adNGLl265MluAAAA7A0qWkyrnnjiCVm/fr0pqK1Vq5b88ssvVncFAABgb1DRkzElK1GihPzwww9SsmRJadiwodVdAQAA2BtUhgwZkuZQ5Jw5c8qiRYukT58+t1xaHwAA4J6to6JB5XYFtgAAAPc8qPzjH/+QJk2aSLZs2cztO9WvNGvWzM72AQAAP+ZWUGnZsqUcP37cLGWtt+8UVDjrMQAAuKdBRc/tc6vbAAAA3sSZvgAAQOYeUZkwYYLbO+zVq1dG2gMAAGAtqHz44YfuPM3UqBBUAADAPQ0qBw8etO0NAQAA3EWNCgAAuH8WfFOHDx8266kcOnRIrl69muaxsWPH2tU2AADg5ywHlRUrVkjz5s2ldOnS8tNPP0l4eLg5IaGeAygiIsI7rQQAAH7J8tTPgAEDpG/fvhIXFyfBwcGyYMECiY+Pl3r16knbtm2900oAAOCXLAeV3bt3S4cOHcztrFmzyqVLl8xJCocPHy6jR4/2RhsBAICfshxUcuXKlVKXUqRIEdm/f3/KYydPnrS3dQAAwK9ZrlF54oknZM2aNVKpUiVp2rSpvP7662YaaOHCheYxAAAAnwUVParn/Pnz5vawYcPM7b///e9Srlw5jvgBAAC+DSp6tE/qaaCPP/7Y3hYBAAB4WqOyceNGWb9+fbrtum3Tpk1WdwcAAGBfUOnRo4c5HPlmR44cMY8BAAD4LKjs2rXrlgu7PfLII+YxAAAAnwWVoKAg+e2339JtP3bsmFlXBQAAwGdBpVGjRmZ12rNnz6ZsO3PmjAwcOFAaNmxoW8MAAAAsD4G8//77UrduXQkLCzPTPWrbtm1SqFAhmTVrljfaCAAA/JTloPLQQw/Jjh07ZPbs2bJ9+3bJkSOHREZGSrt27SRbtmzeaSUAAPBLHhWV6Pop3bp1s781AAAAGQ0qen6fcePGmRMUqsqVK0vv3r2lTJkynuwOAADAnmLa2NhYE0w2bNgg1apVMxdd7K1KlSqyfPlyq7sDAACwb0Slf//+0qdPHxk1alS67W+++SZH/gAAAN+NqOh0T5cuXdJt79y5Mwu+AQAA3waVkJAQczjyzXRbaGioXe0CAACwPvXTtWtXc8TPgQMHpFatWmbb2rVrZfTo0RIdHe2NNgIAAD9lOagMHjxY8uTJIx988IFZoVYVLVpUhg4dKr169fJGGwEAgJ+yHFQCAgJMMa1eEhMTzTYNLgAAAHbL0FkECSgAAMDnQUXP6aMjKe7YsmVLRtsEAADgflBp2bKlO08DAAC490FlyJAh9r4rAACAN9ZRAQAAuFcIKgAAwLEIKgAAwLEIKgAA4P4JKt999513WgIAAJDRoPLHP/5RypQpI++8847Ex8dbfTkc6MaNJIlb9avsXP2rudb7sI5+BAAHBJUjR47Ia6+9JjExMVK6dGlp3LixzJs3T65evWr5zSdPnizVqlWTvHnzmkvNmjXl66+/trwfeO6HhbulS8kJ8laDv8nC9/9trvW+bof76EcAcEhQefDBB815frZt2ybr16+X8uXLS/fu3c2JCfWkhNu3b3d7X8WKFZNRo0bJ5s2bZdOmTfL0009LixYt5Mcff7TaLHhAv0RHtImRk4f/e86mZCePJJrtfMm6h34EAIcW00ZERJgzKOsIy/nz52XatGlSo0YNqVOnjltho1mzZtK0aVMpV66cCTzvvvuu5M6dW9atW5eRZsENOi3xSe9YEdctHvx92ydR3zB9cRf0IwA48KSE165dk8WLF5tgsnz5cnn00Uflr3/9q7Rr104SEhLkrbfekrZt28quXbvc3ueNGzdk/vz5cuHCBTMFdCtXrlwxl2Tnzp0z10lJSeYC92kdxamjiRLwe1TVaz2dU/J9derIOfO8qvXCfNZOp6MfvUc/0y6Xi8+2DehLe9CP9rHSh5aDSs+ePWXOnDnmH+uVV16RMWPGSHh4eMrjuXLlkvfff99MBbkjLi7OBJPLly+b0ZRFixZJ5cqVb/nckSNHyrBhw9Jt13Ckr4f7EhJOSPEa/zv7tX65Plg2h7ntcqV93okT/92O9OhH7/5HdvbsWfN/TWAgKylkBH1pD/rRPomJaafKbQ0qOkoyceJEadWqlQQFBd22jsXdw5grVKhg6l30H18LdDt27CirVq26ZVjRaabo6Og0IyrFixeXkJAQU4wL94WEXJL4zf/7RUkeAYjfkiiuVEE3JCRUQkNDfdDCzIF+9O6Xgp61XT/ffClkDH1pD/rRPsHBwd4LKitWrLj7TrNmlXr16rm1v+zZs0vZsmXNba1v2bhxo4wfP16mTJmS7rkajG4VjvQXhl8aa8LrhknBonlMwWdyLYWOAOiXq/mC1ZGBYnnN8+jb26MfvUu/FPh824O+tAf9aA8r/edWUPnHP/7h9g6bN28uGU2sqetQ4B1ZsgRKt/GNzVEp+mWaxu/3u41rZJ6H26MfAcC73AoqLVu2dDtpalGsu3Qqp0mTJlKiRAkzX/XFF1/IypUrJTY21u19wHO1WlWSgTFtzFErWhCaTEcA9MtVH8fd0Y8A4OOg4q0K5xMnTkiHDh3k2LFjki9fPrP4m4aUhg0beuX9kJ5+iT7eooI5KkULPrWWQqcpGAGwhn4EAO+wXKMyc+ZMeeGFF9LViujKtHPnzjXBw11Tp061+vbwAv0y1UNn9agULfhk7tUz9CMA2M/y/6SRkZHmCJ2b6dSNPgYAAOCzoKLHj2stys0OHz5spm8AAADu+dTPI488YgKKXp555hlzCHIyLaA9ePCgObMyAADAPQ8qyUf+6OJsesZkXUU29VooJUuWlNatW9vWMAAAALeDypAhQ8y1BpIXX3zxtqvSAgAA+KxGRZe211GVm61fv142bdpkV7sAAACsB5UePXpIfHx8uu1HjhwxjwEAAPgsqOhJCSMiIm5ZbKuPAQAA+CyoaG3Kb7/9lm67ri6b+kggAACAex5UGjVqZM7Rk3rRtzNnzsjAgQNZ+h4AANjK8hDI+++/L3Xr1pWwsDAz3aO0uLZQoUIya9Yse1sHAAD8muWg8tBDD8mOHTtk9uzZsn37dsmRI4dZOr9du3aSLVs277QSAAD4JY+KSnLlyiXdunWzvzUAAACpeHR6V53iqV27thQtWlR+/fVXs+3DDz+UxYsXe7I7AAAAe4LK5MmTJTo6Wpo0aSKnT5825/lRDzzwgIwbN87q7gAAAOwLKhMnTpRPP/1UBg0alOZw5EcffVTi4uKs7g4AAMC+oKJnSU4+2ufm9VUuXLhgdXcAAAD2BZVSpUrd8lw/y5Ytk0qVKlndHQAAgH1H/Wh9ip7T5/Lly+JyuWTDhg0yZ84cGTlypHz22WdWdwcAAGBfUPnzn/9s1k5566235OLFi/LSSy+Zo3/Gjx8vL774otXdAQAA2BNUrl+/Ll988YU0btxY2rdvb4LK+fPnJTQ01MpuAAAA7K9R0aN8Xn31VTPto3LmzElIAQAAzimm/cMf/iBbt271TmsAAAAyUqPSvXt3ef311+Xw4cNSo0YNs5x+atWqVbO6SwAAAHuCSnLBbK9evVK2BQQEmCOA9Dp5pVoAAIB7HlR0wTcAAADHBZVr167J008/LUuXLmVxNwAA4Kxi2mzZsqUc8QMAAOC4o350VdrRo0ebNVUAAAAcVaOyceNGWbFihXzzzTdStWrVdEf9LFy40M72AQAAP2Y5qOTPn19at27tndYAAABkJKhMnz7d6ksAAADuTVBJlpCQIHv27DG3K1SoICEhIZ7uCgAAwJ5i2gsXLkjnzp2lSJEiUrduXXPRsyd36dLFnKQQAADAZ0ElOjpaVq1aJUuWLJEzZ86Yy+LFi802XVofAADAZ1M/CxYskJiYGKlfv37KtqZNm0qOHDnk+eefl8mTJ9vWOAAA4N8sj6jo9E6hQoXSbQ8NDWXqBwAA+Dao1KxZU4YMGZJmhdpLly7JsGHDzGMAAAA+m/oZP368NG7cWIoVKybVq1c327Zv3y7BwcESGxtrW8MAAAAsB5Xw8HD5+eefZfbs2fLTTz+Zbe3atZP27dubOhUAAACfrqOSM2dO6dq1q22NAAAAsKVGZeTIkTJt2rR023WbnqwQAADAZ0FlypQpUrFixXTbq1SpIh9//LFd7QIAALAeVI4fP25Wpb2ZLqF/7Ngxu9oFAABgPagUL15c1q5dm267btOl9AEAAHxWTKtFtFFRUXLt2jV5+umnzbYVK1bIG2+8wRL6AADAt0GlX79+curUKenevbtcvXrVbNM1VN58800ZMGCAva0DAAB+zXJQCQgIMEf3DB48WHbv3m3WTilXrpwEBQV5p4UAAMBvebSOisqdO7c89thj9rYGAAAgI8W0AAAA9wpBBQAAOJZPg4qucqvTR3ny5JHQ0FBp2bKl7Nmzx5dNAgAAmS2oREREyOnTp83t4cOHy8WLF21581WrVkmPHj1k3bp1snz5cnPIc6NGjeTChQu27B8AAPhBMa0e3aPh4YEHHpBhw4bJq6++ak5MmFHLli1Lc3/GjBlmZGXz5s1St27dDO8fAAD4QVB5+OGHJTIyUmrXri0ul0vef/99c9TPrbz99tseN+bs2bPmukCBArd8/MqVK+aS7Ny5c+Y6KSnJXOA57T/9t6UfM4Z+tA99aR/60h70o32s9GGAS3v9LrRuZMiQIbJ//37ZsmWLVK5cWbJmzXrLNVb0cU8b3bx5czlz5oysWbPmls8ZOnSoGdG52d69e02dCzyn/a9BMV++fBIYSI21p+hH+9CX9qEv7UE/2icxMVHKly9v+jNv3rwZDyqp6T+OnphQp2js9H//93/y9ddfm5BSrFgxt0dU9NxDWj9ztx8Ud/8AJiQkmJNL8gH0HP1oH/rSPvSlPehH++j3t5aTuBNULC/45o0hr9dee02WLl0qq1evvm1IUbr67a1WwNVfGH5pMk5HxOjLjKMf7UNf2oe+tAf9aA8r/efRyrQ6BTRu3DhTZKt0Kqh3795SpkwZS/vRwZyePXvKokWLZOXKlVKqVClPmgMAAO5TliNhbGysCSYbNmyQatWqmcv69eulSpUq5hBjK/TQ5L/97W/yxRdfmBoTnVLSy6VLl6w2CwAA3Icsj6j0799f+vTpI6NGjUq3Xc+g3LBhQ7f3NXnyZHNdv379NNunT58unTp1sto0AADg70FFp3vmzZuXbnvnzp3NdJAVFut4AQCAn7E89aPVztu2bUu3XbfZfSQQAADwb5ZHVLp27SrdunWTAwcOSK1atcy2tWvXyujRoyU6OtobbQQAAH7KclAZPHiwKXz94IMPZMCAAWZb0aJFzWJsvXr18kYbAQCAn8rqyTHkWkyrF11ZTrEqLAAA8AaP1lFJRkABAADexNJ6AADAsQgqAADAsQgqAADg/ggq165dk2eeeUZ+/vln77UIAADAk6CSLVs22bFjh5WXAAAA3Lupn5dfflmmTp3q+TsCAAB46/Dk69evy7Rp0+Rf//qX1KhRQ3LlypXm8bFjx1rdJQAAgD1BZefOnRIREWFu7927N91icAAAAD4LKt99951tbw4AAOCVw5P37dsnsbGxcunSJXPf5XJ5uisAAAB7gsqpU6fMIcrly5eXpk2byrFjx8z2Ll26yOuvv251dwAAAPYFFT0ZoR6mfOjQIcmZM2fK9hdeeEGWLVtmdXcAAAD21ah88803ZsqnWLFiabaXK1dOfv31V6u7AwAAsG9E5cKFC2lGUpL95z//kaCgIKu7AwAAsC+o1KlTR2bOnJnmkOSkpCQZM2aMPPXUU1Z3BwAAYN/UjwYSLabdtGmTXL16Vd544w358ccfzYjK2rVrre4OAADAvhGV8PBws9Bb7dq1pUWLFmYqqFWrVrJ161YpU6aM1d0BAADYN6Ki8uXLJ4MGDfLkpQAAAN4NKqdPnzYnJty9e7e5X7lyZYmMjJQCBQp4sjsAAAB7pn5Wr14tJUuWlAkTJpjAohe9XapUKfMYAACAz0ZUevToYRZ3mzx5smTJksVsu3HjhnTv3t08FhcXZ1vjAACAfwv05Bw/ulR+ckhRejs6Oto8BgAA4LOgEhERkVKbkppuq169ul3tAgAAcG/qZ8eOHSm3e/XqJb179zajJ0888YTZtm7dOpk0aZKMGjXKey0FAAB+x62g8vDDD5sVaF0uV8o2XejtZi+99JKpXwEAALhnQeXgwYO2vBkAAIDtQSUsLMzSTgEAAHy24NvRo0dlzZo1cuLECXNCwtS0hgUAAMAnQWXGjBnyl7/8RbJnzy4FCxY0tSvJ9DZBBQAA+CyoDB48WN5++20ZMGCABAZaProZAADAbZaTxsWLF+XFF18kpAAAAK+znDa6dOki8+fP905rAAAAMjL1M3LkSHnuuedk2bJlUrVqVcmWLVuax8eOHWt1lwAAAPYFldjYWKlQoYK5f3MxLQAAgM+CygcffCDTpk2TTp062dYIAAAAW2pUgoKC5Mknn7T6MgAAAO8HFT0h4cSJE62/EwAAgLenfjZs2CDffvutLF26VKpUqZKumHbhwoVWdwkAAGBPUMmfP7+0atXK6ssAAAC8H1SmT59u/V0AAAA8wPKyAADg/hlRKVWq1B3XSzlw4EBG2wQAAOBZUImKikpz/9q1a7J161azUm2/fv2s7g4AAMC+oKKHJ9/KpEmTZNOmTVZ3BwAA4P0alSZNmsiCBQvs2h0AAIB9QSUmJkYKFChg1+4AAACsT/088sgjaYppXS6XHD9+XBISEuSjjz6ytK/Vq1fLe++9J5s3b5Zjx47JokWLpGXLllabBAAA7lOWg8rNQSIwMFBCQkKkfv36UrFiRUv7unDhglSvXl06d+7MInIAACDjQWXIkCFiZ12LXtx15coVc0l27tw5c52UlGQu8Jz2n46O0Y8ZQz/ah760D31pD/rRPlb60HJQ8aWRI0fKsGHD0m3XaafLly/7pE330y/N2bNnzYdQR8ngGfrRPvSlfehLe9CP9klMTLQ/qOg/yp0WelP6+PXr18VbBgwYINHR0WlGVIoXL26mnvLmzeu19/WXD6D++2lf8gH0HP1oH/rSPvSlPehH+wQHB9sfVLTQ9Xb+/e9/y4QJE7w+HBYUFGQuN9NfGH5pMk4/gPRlxtGP9qEv7UNf2oN+tIeV/nM7qLRo0SLdtj179kj//v1lyZIl0r59exk+fLj7rQQAALgLjyLh0aNHpWvXrlK1alUz1bNt2zb5/PPPJSwszJPdAQAAZLyYVouIRowYIRMnTpSHH35YVqxYIXXq1BFPnT9/Xvbt25dy/+DBgyb06MJxJUqU8Hi/AADAz4LKmDFjZPTo0VK4cGGZM2fOLaeCrNJzAz311FMp95MLZTt27CgzZszI8P4BAICfBBWtRcmRI4eULVvWTPPo5VYWLlzo9pvrInF6mBcAAECGgkqHDh3uengyAACAT4IKUzEAAOBe40BwAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWAQVAADgWFklE3O5XOb63Llzvm5KppeUlCSJiYkSHBwsgYHkV0/Rj/ahL+1DX9qDfrRP8vd28vf4fRtU9BdGFS9e3NdNAQAAHnyP58uX747PCXC5E2ccnG6PHj0qefLkkYCAAF83J9OnWw188fHxkjdvXl83J9OiH+1DX9qHvrQH/WgfjR4aUooWLXrX0alMPaKiP1yxYsV83Yz7in74+ABmHP1oH/rSPvSlPehHe9xtJCUZk2wAAMCxCCoAAMCxCCowgoKCZMiQIeYanqMf7UNf2oe+tAf96BuZupgWAADc3xhRAQAAjkVQAQAAjkVQAQAAjkVQAQAAjkVQ8XOrV6+WZs2amdUBdXXfL7/80tdNypRGjhwpjz32mFklOTQ0VFq2bCl79uzxdbMypcmTJ0u1atVSFtWqWbOmfP31175uVqY3atQo8xmPiorydVMynaFDh5q+S32pWLGir5vlNwgqfu7ChQtSvXp1mTRpkq+bkqmtWrVKevToIevWrZPly5fLtWvXpFGjRqZ/YY2uNq1fqps3b5ZNmzbJ008/LS1atJAff/zR103LtDZu3ChTpkwxARCeqVKlihw7dizlsmbNGl83yW9k6iX0kXFNmjQxF2TMsmXL0tyfMWOGGVnRL9u6dev6rF2ZkY7wpfbuu++aURYNgfplAWvOnz8v7du3l08//VTeeecdXzcn08qaNasULlzY183wS4yoAF5w9uxZc12gQAFfNyVTu3HjhsydO9eMTOkUEKzTkb5nn31WGjRo4OumZGo///yzmSIvXbq0CX6HDh3ydZP8BiMqgBfO6q11AE8++aSEh4f7ujmZUlxcnAkmly9flty5c8uiRYukcuXKvm5WpqMhb8uWLWbqB557/PHHzShphQoVzLTPsGHDpE6dOrJz505TlwbvIqgAXvgLVv8DYw7bc/qFsG3bNjMyFRMTIx07djR1QIQV98XHx0vv3r1NzVRwcLCvm5OppZ4e1zofDS5hYWEyb9486dKli0/b5g8IKoCNXnvtNVm6dKk5mkqLQuGZ7NmzS9myZc3tGjVqmBGB8ePHm4JQuEfro06cOCERERFpptL0d/Ovf/2rXLlyRbJkyeLTNmZW+fPnl/Lly8u+fft83RS/QFABbKCnzOrZs6eZoli5cqWUKlXK102676bT9IsV7nvmmWfMFFpqkZGR5rDaN998k5CSwQLl/fv3yyuvvOLrpvgFgoqf0w9c6r8KDh48aIbctQi0RIkSPm1bZpvu+eKLL2Tx4sVmzvr48eNme758+SRHjhy+bl6mMmDAADPUrr9/iYmJpl81/MXGxvq6aZmK/h7eXCOVK1cuKViwILVTFvXt29ccjabTPUePHjVnUNag165dO183zS8QVPycrlPx1FNPpdyPjo4211oToMVjcI8ePqvq16+fZvv06dOlU6dOPmpV5qTTFR06dDBFixr0tCZAQ0rDhg193TT4qcOHD5tQcurUKQkJCZHatWubw+X1NrwvwKVj1gAAAA7EOioAAMCxCCoAAMCxCCoAAMCxCCoAAMCxCCoAAMCxCCoAAMCxCCoAAMCxCCoAAMCxCCoA3BYQECBffvmlr5sBwI8QVIBMTpfo1wChl2zZspkTIr7xxhty+fJlXzftvu7zli1b+roZgF/gXD/AfeCPf/yjOa/QtWvXZPPmzeZcTRpcRo8eLfcj/Tk1lGV2V69elezZs/u6GYCjMaIC3AeCgoKkcOHCUrx4cfOXfoMGDWT58uUpj+vJ1PSkag899JDkzJlTqlatKnPmzEmzDz2hYq9evcxojJ49W/c3dOjQO76vnkW2SJEismPHjls+rq9/+OGHZcqUKaZt+t7PP/+8nD17NuU5GzduNCccfPDBB81JCOvVqydbtmxJsx8NXXrix+bNm5szAL/77rty48YN6dKlixlB0jNUV6hQQcaPH3/LkY8RI0ZIoUKFJH/+/DJ8+HC5fv269OvXz/ycxYoVMyEvtfj4eNNOfb4+p0WLFvLLL7+k/Eyff/65OVN28kiWnt35bq9L3R5tf9GiRU2bAdwZQQW4z+zcuVN++OGHNH+p6zRQjRo15J///Kd5vFu3bvLKK6/Ihg0b0rxWv4A1CKxfv17GjBljvtRTB55kei7Tnj17ysyZM+X77783Zzi+nX379sm8efNkyZIlsmzZMtm6dat079495fHExEQzArRmzRpzRtpy5cpJ06ZNzfbUNCD86U9/kri4OOncubMkJSWZkDF//nzZtWuXvP322zJw4EDzXql9++23cvToUVm9erWMHTvWhKvnnntOHnjgAfNzvvrqq/KXv/zFnCE3ebSmcePGkidPHvOzrV27VnLnzm1GrXQEpG/fviaM6H09w7NeatWqddfXJVuxYoXs2bPH9OvSpUst/dsCfknPngwg8+rYsaMrS5Ysrly5crmCgoL0bOiuwMBAV0xMzB1f9+yzz7pef/31lPv16tVz1a5dO81zHnvsMdebb76Zcl/3PX/+fNdLL73kqlSpkuvw4cN3fI8hQ4aYtqV+3tdff23ad+zYsVu+5saNG648efK4lixZkuZ9o6KiXHfTo0cPV+vWrdP0TVhYmNlnsgoVKrjq1KmTcv/69eum7+bMmWPuz5o1yzwnKSkp5TlXrlxx5ciRwxUbG5uy3xYtWqR5b3dfV6hQIbMdgHuoUQHuA0899ZSZGrlw4YJ8+OGHkjVrVmndunXK4zpNotMfOtpw5MgR8xf+lStXzFRMajePjOi0zokTJ9Js69Onj5lq0tEPna65mxIlSpgpp2Q1a9Y0oyE6qqDTS7/99pu89dZbZvpE30vbevHiRTl06FCa/Tz66KPp9j1p0iSZNm2aee6lS5fMz6VTTalVqVJFAgP/N3isU0Dh4eEp97NkySIFCxZM+Tm3b99uRoF0ZCQ1HZXav3//bX9Od1+n027UpQDuI6gA9wGdrilbtqy5rV/c1atXl6lTp5oaDvXee++Z+o1x48aZL0p9flRUVJopCXVzgarWX2ioSE3rSbS+JTY2Vtq3b5/htuu0j9bQaPvCwsJMCNIwc3PbtM2pzZ0710zDfPDBB+b5GhD059TpnLv9THf6Oc+fP2+myWbPnp2urSEhIbf9Odx93c0/B4A7I6gA9xkdPdBajejoaHnppZdMoanWS2hh58svv2yeo1/Ke/fulcqVK1vevxa0NmvWzOxbRyNefPHFOz5fRzu0RkSLR5WOxGgbkwtJtW0fffSRqUtJLkg9efLkXduhr9PakNT1Lnca8XBXRESE/P3vf5fQ0FDJmzfvLZ+jIyI68mP1dQCso5gWuA+1bdvWhAidGlFaoKrFm1pku3v3blM8qlMuntKi1lmzZklkZKTExMTc8bnBwcFm1ESnRrTIVI8s0mJUnfZJbpvuS9uloyE6SqPh6m70dZs2bTIjOxq6Bg8ebI4gyih9f53S0mCn7T148KCZltJ2JxfclixZ0hzppNNXGqq0kNad1wGwjqAC3Ie0RuW1114zR+5o3YrWgOhf/HpUih6GrCEhowuWtWnTxhwlpEcPLVy48LbP0ympVq1amRGTRo0amToYHUFJplNUp0+fNu3TfekXu45K3I2GLd3vCy+8II8//riZPko9uuIprdvRI4S0tkb3X6lSJTOFprUmySMlXbt2NSNCWjej0zo6uuPO6wBYF6AVtR68DgDuSg8p1iX3t23b5uumAMikGFEBAACORVABAACOxdQPAABwLEZUAACAYxFUAACAYxFUAACAYxFUAACAYxFUAACAYxFUAACAYxFUAACAYxFUAACAONX/A4+MhzjygLbaAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGwCAYAAACzXI8XAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAANK1JREFUeJzt3QlUldXex/E/TuAADgjiAGSZs1nXyovlkEOmXYfUm9qgmUsbzHlKy8pK8XbfTCqz7Jba4FCm3qxXKWdN7c0pp7I0CwcUnABN0OC8679vcDkqeoYHznng+1nrLA7Pec5h8wc8P/fez94BDofDIQAAADZUzNcNAAAA8BRBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2BZBBgAA2FYJKeSysrLk6NGjEhwcLAEBAb5uDgAAcIEuc5eWlibVqlWTYsWKFd0goyEmMjLS180AAAAeOHTokNSoUaPoBhntickuREhIiK+bY/vereTkZAkLC7tqOsbVUUfrUEvrUEtrUEfrpKammo6I7PfxIhtksoeTNMQQZLz/A01PTzd15A/Uc9TROtTSOtTSGtTReteaFkKVAQCAbfk0yMyYMUNuuummnN6SmJgYWbZsWc7jrVq1Mkks9+3xxx/3ZZMBAIAf8enQkk7emTJlitx4441mdvKcOXOkS5cusn37dmnQoIE5Z8CAAfLiiy/mPKdMmTI+bDEAAPAnPg0ynTp1cvp80qRJppdm8+bNOUFGg0tERISPWggAAPyZ30z2zczMlE8//VTOnTtnhpiyffzxx/LRRx+ZMKPBZ8KECVftlcnIyDC33LOesydg6Q2e0/ppzxl19A51tA61tA61tAZ1tI6rNfR5kNm1a5cJLjrLu1y5crJ48WKpX7++eeyBBx6Q6OhosxjOzp07ZezYsbJv3z5ZtGhRnq8XGxsrEydOvOy4Xg6nXwPe/VKlpKSYP1Jm43uOOlqHWlqHWlqDOlpHF8NzRYBDq+1DFy5ckISEBPODX7hwofzrX/+StWvX5oSZ3FatWiVt2rSR/fv3yw033OByj4xeh3769Gkuv/YS6yNYgzpah1pah1pagzpaR9+/K1asaPLB1d6/fd4jU6pUKalVq5a536RJE/nuu+8kLi5O3nnnncvObdq0qfl4tSATGBhobpfSXyh+qbynV45RS+9RR+tQS+tQS2tQR2u4Wr9i/phmc/eo5LZjxw7zsWrVqgXcKgAA4I982iMzbtw46dChg0RFRZmxsLlz58qaNWskPj5eDhw4YD7v2LGjhIaGmjkyw4cPlxYtWpi1Z1CwMjOzZPe63yQ5OUnCws5LwxbRUry43+VgAEAR49Mgk5SUJH369JHExEQpX768CSgaYtq1a2f2RlqxYoVMmzbNXMmk81y6d+8uzz77rC+bXCRtXPSDzBwaLyePpklkk2A5tDVNQqsFy8C49tKsWz1fNw8AUIT5NMi89957eT6mwUUn/cL3IWZyj4UiDpGAXB0wJ46kmePjF/YgzAAAfIaxAVx1OEl7YjTEXObPYzOHfWXOAwDAFwgyyNOe9Qly4vBVruN3iJw4lGrOAwDAFwgyyNPpxLOWngcAgNUIMshTxarlLD0PAACrEWSQpwbNo6RyjWCRgDxOCBCpHBlizgMAwBcIMsiTrhOjl1gbl4aZPz8fOO1u1pMBAPgM70C4Kr20Wi+xrlw92Ol45RohXHoNAPA5n++1BP+nYaVplzq5VvYNZ2VfAIBfIMjAJRpaGrWMlqSk0hIeHs5maAAAv8C7EQAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2CDAAAsC2fBpkZM2bITTfdJCEhIeYWExMjy5Yty3k8PT1dBg0aJKGhoVKuXDnp3r27HD9+3JdNBgAAfsSnQaZGjRoyZcoU2bp1q2zZskVat24tXbp0kT179pjHhw8fLkuXLpVPP/1U1q5dK0ePHpVu3br5sskAAMCPlPDlF+/UqZPT55MmTTK9NJs3bzYh57333pO5c+eagKNmzZol9erVM4//9a9/9VGrAQCAv/BpkMktMzPT9LycO3fODDFpL83Fixelbdu2OefUrVtXoqKiZNOmTXkGmYyMDHPLlpqaaj5mZWWZGzyn9XM4HNTRS9TROtTSOtTSGtTROq7W0OdBZteuXSa46HwYnQezePFiqV+/vuzYsUNKlSolFSpUcDq/SpUqcuzYsTxfLzY2ViZOnHjZ8eTkZPM14N0vVUpKivkjLVaMeeKeoo7WoZbWoZbWoI7WSUtLs0eQqVOnjgkt+oNfuHCh9O3b18yH8dS4ceNkxIgRTj0ykZGREhYWZiYUw7s/0ICAAFNL/kA9Rx2tQy2tQy2tQR2tExQUZI8go70utWrVMvebNGki3333ncTFxUnPnj3lwoULcubMGadeGb1qKSIiIs/XCwwMNLdL6S8Uv1Te0z9Qauk96mgdamkdamkN6mgNV+tXzB/TrM5x0VBTsmRJWblyZc5j+/btk4SEBDMUBQAA4NMeGR0G6tChg5nAq2NheoXSmjVrJD4+XsqXLy/9+/c3w0SVKlUyw0KDBw82IYYrlgAAgM+DTFJSkvTp00cSExNNcNHF8TTEtGvXzjz+2muvma4lXQhPe2nat28vb731Fj85AADg+yCj68Rca6LP9OnTzQ0AAMDv58gAAAC4iiADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsiyADAABsq4QrJ33++eduv3C7du2kdOnSnrQJAADAuiDTtWtXcUdAQID8/PPPcv3117v1PAAAgHwZWjp27JhkZWW5dCtTpoxbjQAAAMi3INO3b1+3hokeeughCQkJ8ahBAAAAlg4tzZo1S9wxY8YMt84HAADwBFctAQCAwt0jk1t6erq88cYbsnr1aklKSjJzYnLbtm2ble0DAACwLsj0799fvvrqK+nRo4fcfvvt5golAAAAWwSZL774Qv73f/9X7rjjjvxpEQAAQH7NkalevboEBwe7+zQAAADfB5lXX31Vxo4dK7/99pv1rQEAAMjPoaVbb73VTPjVVXt14buSJUs6PX7q1Cl3XxIAAKBggkzv3r3lyJEjMnnyZKlSpQqTfQEAgH2CzMaNG2XTpk3SuHHj/GkRAABAfs2RqVu3rpw/f97dpwEAAPg+yEyZMkVGjhwpa9askZMnT0pqaqrTzR2xsbFy2223maugwsPDzS7b+/btczqnVatWZvgq9+3xxx93t9kAAKAQcnto6Z577jEf27Rp43Tc4XCYkJGZmenya61du1YGDRpkwswff/wh48ePl7vvvlv27t0rZcuWzTlvwIAB8uKLL+Z8zu7aAADAoyCjWxNYZfny5U6fz5492/TMbN26VVq0aOEUXCIiIlx6zYyMDHPLlt1LpFspXLqdAtyj9dPASh29Qx2tQy2tQy2tQR2t42oN3Q4yLVu2lPySkpJiPlaqVMnp+McffywfffSRCTOdOnWSCRMm5Nkro8NVEydOvOx4cnKyuWwc3v1S6c9I/0iLFWO/UU9RR+tQS+tQS2tQR+ukpaW5dF6AQ6vtJg0EO3fuvOKmkZ07dxZP6Ovoc8+cOSMbNmzIOT5z5kyJjo6WatWqma+pi/HpHk+LFi1yuUcmMjJSTp8+LSEhIR61Df/9GWkgDAsL4w/UC9TROtTSOtTSGtTROvr+XbFiRRMMr/b+XcKT4aA+ffrIiRMnLnvM3Tkyuelcmd27dzuFGDVw4MCc+40aNZKqVaua+TkHDhyQG2644bLXCQwMNLdL6S8Uv1Te058xtfQedbQOtbQOtbQGdbSGq/Vzu8qDBw+Wv//975KYmJgz7yT75mmIeeqpp8xmlDr/pkaNGlc9t2nTpubj/v37PfpaAACg8HC7R+b48eMyYsQIs6qvt3RUS4PR4sWLzeXcNWvWvOZzduzYYT5qzwwAACja3A4yPXr0MKHjSsM6ngwnzZ07V/7973+btWSOHTtmjpcvX15Kly5tho/08Y4dO0poaKiZIzN8+HBzRdNNN93k9dcHAABFLMi8+eabZmhp/fr1Zs7KpZtGDhkyxOXXmjFjRs6id7nNmjVLHnnkESlVqpSsWLFCpk2bJufOnTOTdrt37y7PPvusu80GAACFkNtBZt68efLVV19JUFCQ6ZnJvWmk3ncnyFzrgikNLrpoHgAAgCVB5plnnjHrtDz99NPMyAYAAD7ldhK5cOGC9OzZkxADAAB8zu000rdvX1mwYEH+tAYAACA/h5Z0rZhXXnlF4uPjzZVDl072nTp1qrsvCQAAUDBBZteuXXLLLbeY+7oSb265J/4CAAAU6t2vAQAAvMGMXQAAULiDTLdu3cwulK568MEHzc7YAAAAPh9a0i0EdFtyVxe5W7p0qbz00ksSHh7ubfsAAAC8CzIaTmrXru3KqQAAAP4VZDyZ4Fu9enVP2gMAAGBtkGnZsqXrrwgAAFBAuGoJAADYFkEGAADYFkEGAAAUjSCjVy8lJCRIenp6/rUIAAAgv4JMrVq15NChQ+48DQAAwPdBplixYnLjjTfKyZMn86c1AAAA+TlHZsqUKTJ69OjLdr4GAADw+92v+/TpI7///rs0btxYSpUqJaVLl3Z6/NSpU1a2DyhUMjOzZPe63yQ5OUnCws5LwxbRUrw4c+4BoMCCzLRp0zz+YkBRtnHRDzJzaLycPJomkU2C5dDWNAmtFiwD49pLs271fN08ACgaQaZv37750xKgkIeYyT0WijhEAnJ1wJw4kmaOj1/YgzADAAURZFRmZqYsWbJEfvjhB/N5gwYNpHPnzlK8eHFPXg4o9MNJ2hOjIeYyeixAZOawr6RplzoMMwFAfgeZ/fv3S8eOHeXIkSNSp04dcyw2NlYiIyPlyy+/lBtuuMHdlwQKtT3rE+TE4bS8T3CInDiUas67qdV1Bdk0ALA9t//7N2TIEBNWdC2Zbdu2mZsuklezZk3zGABnpxPPWnoeAMCLHpm1a9fK5s2bpVKlSjnHQkNDzWXZd9xxh7svBxR6FauWs/Q8AIAXPTKBgYGSlnZ5N/nZs2fN5dgAnDVoHiWVawSbuTBXFCBSOTLEnAcAyOcg87e//U0GDhwo3377rdmyQG/aQ/P444+bCb8AnOkEXr3E2rg0zPz5+cBpdzPRFwA84Pa/nK+//rqZIxMTEyNBQUHmpkNKugdTXFycJ20ACj29tFovsa5cPdjpeOUaIVx6DQAFNUdGe19SU1Nl/vz55qql7Muv69WrZ4IMgLxpWNFLrP+7sm84K/sCQEEHGQ0se/bsMZtHEl4A92hoadQyWpKSSkt4eLjZiBUA4Dl2vwYAALbF7tcAAMC22P0aAADYFrtfAwCAohFkLl68aFb2nTBhgtmSAAAAwDZzZEqWLCmfffZZ/rUGAAAgPyf7du3aVZYsWeLu0wAAAHw/R0Yvv37xxRflm2++kSZNmkjZsmWdHndnB+zY2FhZtGiR/Pjjj2bScLNmzeQf//iH1KlTJ+ec9PR0GTlypFmELyMjQ9q3by9vvfWWVKlSxd2mAwCAQibAoavcueFqc2MCAgLkl19+cfm17rnnHunVq5fcdttt8scff8j48ePNZd179+7NCUhPPPGEfPnllzJ79mwpX768PPXUU2Y9Gw1SrtCViPV5KSkpEhIS4nLbcLmsrCxJSkpiITcvUUfrUEvrUEtrUEfruPr+7XaPzMGDB8Uqy5cvd/pcw4r+8Ldu3SotWrQwjX/vvfdk7ty50rp1a3POrFmzzJYIulHlX//6V8vaAgAA7MftIJPtwoULJtToBpIlSnj8Mk40uKhKlSqZjxpo9Eqptm3b5pxTt25diYqKkk2bNl0xyOjwk95yJ7rslKw3eE7rpx141NE71NE61NI61NIa1NE6rtbQ7QSii+ENHjxY5syZYz7/6aef5PrrrzfHqlevLk8//bTHDR42bJjZSbthw4bm2LFjx8yiexUqVHA6V+fH6GN5zbuZOHHiZceTk5PNfBt4Tn9GGjb1j5QuU89RR+tQS+tQS2tQR+ukpaXlT5AZN26cfP/997JmzRozxyWb9pq88MILHgeZQYMGmfkxGzZs8Oj5uds3YsQIpx6ZyMhICQsLY46MBX+gOg9Ka8kfqOeoo3WopXWopTWoo3WCgoLyJ8jopdcLFiwwwzr6w8rWoEEDOXDggHhCJ/B+8cUXsm7dOqlRo0bO8YiICDOEdebMGademePHj5vHriQwMNDcLqW/UPxSeU9/5tTSe9TROtTSOtTSGtTRGq7Wz+0q6xCNTsi91Llz55yCjSu0601DzOLFi2XVqlWXXRGll3frInwrV67MObZv3z5JSEiQmJgYd5sOAAAKGbeDzK233mouh86WHV7+9a9/uR0udDjpo48+MlclBQcHm3kvejt//rx5XC+76t+/vxkqWr16tZn8269fP/N1uGIJAAC4PbQ0efJk6dChg1nrRdd+iYuLM/c3btxo9mFyx4wZM8zHVq1aOR3XS6wfeeQRc/+1114z3Uvdu3d3WhAPAADA7R6ZO++8U3bs2GFCTKNGjeSrr74yQ016ObQOBbk7tHSlW3aIyZ7sM336dDl16pQZvtKVgPOaHwMAAIoWjxaA0bVj3n33XetbAwAA4AamVAMAANsiyAAAANsiyAAAANsiyAAAgKIXZPbv3y/x8fE5a77o1UYAAAB+HWROnjxp9lWqXbu2dOzYURITE81xXbhu5MiR+dFGAAAAa4LM8OHDpUSJEmabgDJlyuQc79mzpyxfvtzdlwMAACi4dWR0ATwdUsq9uaO68cYb5bfffvO8JQAAAPndI6Or6+buicmmK+9eaddpAAAAvwkyzZs3lw8++MBp08isrCx55ZVX5K677rK6fQAAANYNLWlgadOmjWzZskUuXLggY8aMkT179pgemW+++cbdlwMAACi4HpmGDRvKTz/9ZDaP7NKlixlq6tatm2zfvt3swQQAAODXm0aWL19ennnmGetbAwAAkJ89MrNmzZJPP/30suN6bM6cOe6+HAAAQMEFmdjYWKlcufJlx8PDw2Xy5MmetwQAACC/g4wuhFezZs3LjkdHR5vHAAAA/DbIaM/Lzp07Lzv+/fffS2hoqFXtAgAAsD7I9O7dW4YMGSKrV6+WzMxMc1u1apUMHTpUevXq5e7LAQAAFNxVSy+99JL8+uuvZi0Z3XNJ6YJ4ffr0YY4MAADw7yBTqlQpWbBggQk0OpxUunRpadSokZkjAwAA4PfryKjatWubGwAAgG2CjM6JmT17tqxcuVKSkpLMsFJuOl8GAADAL4OMTurVIHPvvfea7Qp000gAAABbBJn58+fLJ598Ih07dsyfFgEAAOTX5dc62bdWrVruPg0AAMD3QWbkyJESFxcnDofD+tYAAADk59DShg0bzGJ4y5YtkwYNGkjJkiWdHl+0aJG7LwkAAFAwQaZChQpy3333efbVAAAAfBlkZs2aZeXXBwCPZGZmye51v0lycpKEhZ2Xhi2ipXhxt0fLARTVBfEAwFc2LvpBZg6Nl5NH0ySySbAc2pomodWCZWBce2nWrZ6vmwfA34PMwoULzSXYCQkJcuHCBafHtm3bZlXbAOCKIWZyj4UiDpGAXB0wJ46kmePjF/YgzABFiNv9sK+//rr069dPqlSpItu3b5fbb79dQkND5ZdffpEOHTrkTysB4M/hJO2J0RBzmT+PzRz2lTkPQNHgdpB56623ZObMmfLGG2+YNWXGjBkjX3/9tQwZMkRSUlLyp5UAICJ71ifIicNpeZ/gEDlxKNWcB6BocDvI6HBSs2bNzH3d+Tot7T//qDz88MMyb94861sIAH86nXjW0vMAFMEgExERIadOnTL3o6KiZPPmzeb+wYMHWSQPQL6qWLWcpecBKIJBpnXr1vL555+b+zpXZvjw4dKuXTvp2bMn68sAyFcNmkdJ5RrBInntVRsgUjkyxJwHoGhw+6olnR+TlfWfiXSDBg0yE303btwonTt3lsceeyw/2ggAhq4To5dYm6uWLg0zf34+cNrdrCcDFCFu/7UfPnxYihcvnvN5r169zJVMTz31lBw7dsyt11q3bp106tRJqlWrJgEBAbJkyRKnxx955BFzPPftnnvucbfJAAoRvbRaL7GuXD3Y6XjlGiFceg0UQW73yNSsWVMSExMlPDzc6bjOm9HHMjMzXX6tc+fOSePGjeXRRx+Vbt26XfEcDS65VxMODAx0t8kAChkNK0271Mm1sm84K/sCRZTbQUYn9GrPyKXOnj0rQUFBbr2WrjtzrbVnNLjoBGMAyE1DS6OW0ZKUVNr8x6pYMUIMUBS5HGRGjBhhPmqImTBhgpQpUybnMe2F+fbbb+Xmm2+2vIFr1qwx/0hVrFjRTDR++eWXzbycvGRkZJhbttTUVPNR5/Vkz+2BZ7R+GmSpo3eoo3WopXWopTWoo3VcraHLQUZX8VX6A9q1a5dZDC+b3tcholGjRomVdFhJh5x0yOrAgQMyfvx404OzadMmp3k6ucXGxsrEiRMvO56cnCzp6emWtq8o/lLpoof6O8D/fj1HHa1DLa1DLa1BHa2TvU7dtQQ43Fz8RS+5jouLk5CQEE/bduWGBATI4sWLpWvXrnmeo9sg3HDDDbJixQpp06aNyz0ykZGRcvr0acvbXBT/QDUQhoWF8QfqBepoHWppHWppDepoHX3/1tEYDYZXe/92e45M7om32V9o1apVUrduXXPLT9dff71UrlxZ9u/fn2eQ0Tk1V5oQrL9Q/FJZEzippfeoo3WopXWopTWoozVcrZ/bVb7//vvlzTffNPfPnz8vt956qznWqFEj+eyzzyQ/6aXfJ0+elKpVq+br1wEAAPbgdpDRtV+aN29u7utQkI5MnTlzxqwloxNx3aFXOu3YscPcsrc50Pu6n5M+Nnr0aLMFwq+//iorV66ULl26SK1ataR9+/buNhsAABRCbgcZHauqVKmSub98+XLp3r27uYLp3nvvlZ9//tmt19qyZYvccsst5pZ9ZZTef+6558xk3p07d5oVg2vXri39+/eXJk2ayPr161lLBgAAeDZHRifO6lVDGmY0yMyfP98c18m07q4j06pVq6tuNBkfH+9u8wAAQBHidpAZNmyYPPjgg1KuXDmJjo42YSR7yEnnyQAAAPhtkHnyySfl9ttvl0OHDpldr7NnFesVRe7OkQEAACjQIKP0SiW95aZzZAAAAPwuyOgk3JdeeknKli2bs1VBXqZOnWpV2wAAALwPMro9wcWLF3Pu5+VKm0kCAAD4NMisXr36ivcBAAB8ifWTAQBA4e6R0R2oXbVo0SJv2gMAAGBtj0z58uVzbroDpW4XoKvyZtu6das5po8DAAD4VY9M7h2vx44dazaJfPvtt802AiozM9OsL3O1bbYBAAB8Pkfm/fffl1GjRuWEGKX39bJsfQwAAMBvg8wff/whP/7442XH9VhWVpZV7QIAALB+Zd9+/fqZnagPHDhgtipQ3377rUyZMsU8BgAA4LdB5n/+538kIiJCXn31VUlMTDTHqlatKqNHj5aRI0fmRxsBAACsCTK6SeSYMWPMLTU11Rxjki8AALDNppHZCDAAAMCXWNkXAADYFkEGAADYFkEGAAAU7iBTqVIlOXHihLn/6KOPSlpaWn63CwAAwJogc+HChZwrlObMmSPp6emuPA0AAMD3Vy3FxMRI165dpUmTJuJwOGTIkCFSunTpK57LNgUAAMCvgsxHH30kr732mlnNNyAgQFJSUuiVAQAA9ggyVapUMVsQqJo1a8qHH34ooaGh+d02AAAAaxfEO3jwoLtPAQAA8J/Lr9euXSudOnWSWrVqmVvnzp1l/fr11rcOAADAyiCj82Xatm0rZcqUMZN+syf+tmnTRubOnevuywEAABTc0NKkSZPklVdekeHDh+cc0zAzdepUeemll+SBBx7wvDUAAAD52SPzyy+/mGGlS+nwEvNnAACAXweZyMhIWbly5WXHV6xYYR4DAADw26GlkSNHmqGkHTt2SLNmzcyxb775RmbPni1xcXH50UYAAABrgswTTzwhERER8uqrr8onn3xijtWrV08WLFggXbp0cfflAAAACi7IqPvuu8/cAAAAbLeODAAAgD8gyAAAANsiyAAAANsiyAAAANsiyAAAgKJz1VJmZqZZM0YXxUtKSpKsrCynx1etWmVl+wAAAKzrkRk6dKi5aaBp2LChNG7c2OnmjnXr1pntDqpVqyYBAQGyZMkSp8cdDoc899xzUrVqVbMxpW5W+fPPP7vbZAAAUEi53SMzf/58sxBex44dvf7i586dM+Hn0UcflW7dul32uG5O+frrr8ucOXOkZs2aMmHCBGnfvr3s3btXgoKCvP76AACgiAWZUqVKSa1atSz54h06dDC3K9HemGnTpsmzzz6bs2LwBx98IFWqVDE9N7169bri8zIyMswtW2pqqvmoQ2CXDoPBPVo//blQR+9QR+tQS+tQS2tQR+u4WkOP9lrSPZXefPNNMxyUX3Qn7WPHjpnhpGzly5eXpk2byqZNm/IMMrGxsTJx4sTLjicnJ0t6enq+tbeo/FKlpKSYP9JixZgn7inqaB1qaR1qaQ3qaJ20tLT8CTIbNmyQ1atXy7Jly6RBgwZSsmRJp8cXLVokVtAQo7QHJjf9PPuxKxk3bpyMGDHCqUdGd+UOCwuTkJAQS9pWlP9ANbxqLfkD9Rx1tA61tA61tAZ1tI6rU0jcDjIVKlTw632WAgMDze1S+gvFL5X39A+UWnqPOlqHWlqHWlqDOlrD1fq5HWRmzZolBUF32FbHjx83Vy1l089vvvnmAmkDAADwb34bF/UqJQ0zul5N7mGib7/9VmJiYnzaNgAA4B/c7pFRCxcuNJdgJyQkyIULF5we27Ztm8uvc/bsWdm/f7/TBN8dO3ZIpUqVJCoqSoYNGyYvv/yy3HjjjTmXX+uaM127dvWk2QAAoKj3yOi6Lv369TOTbrdv3y633367hIaGyi+//JLnpdR52bJli9xyyy3mpnSSrt7XRfDUmDFjZPDgwTJw4EC57bbbTPBZvnw5a8gAAAAjwKHXiLmhbt268vzzz0vv3r0lODhYvv/+e7n++utN+Dh16pS5LNuf6HCUXratl8Nx1ZL3s/F1W4rw8HAmsXmBOlqHWlqHWlqDOhb8+7fbVdbhpGbNmpn7um1A9nXeDz/8sMybN8+bNgMAALjF7SCjE3C150XpPJbNmzfnzG9xs3MHAACgYINM69at5fPPPzf3da7M8OHDpV27dtKzZ0+/Xl8GAAAUPm5ftTRz5syc/Q8GDRpkJvpu3LhROnfuLI899lh+tBEAAMCaIHPpaoW651Fe+x4BAADkJ4+mVK9fv14eeughszDdkSNHzLEPP/zQ7MMEAADgt0Hms88+k/bt25srlnQdmYyMDHNcL4+aPHlyfrQRAADAmiCjK+2+/fbb8u677zrtfH3HHXe4taovAABAgQeZffv2SYsWLS47rovWnDlzxusGAQAA5Os6Mrn3R8qm82N0hV8AAAC/DTIDBgyQoUOHml2oAwIC5OjRo/Lxxx/LqFGj5IknnsifVgIAAFhx+fXTTz9t1pFp06aN/P7772aYKTAw0AQZ3eARAADAb4OM9sI888wzMnr0aDPEpDtS169fX8qVK5c/LQQAALAqyGQrVaqUCTAAAAB+H2QeffRRl857//33vWkPAACA9UFm9uzZEh0dLbfccgu7XAMAAHsFGb0iad68eXLw4EGz67VuUVCpUqX8bR0AAIAVl19Pnz5dEhMTZcyYMbJ06VKJjIyU+++/X+Lj4+mhAQAA/r+OjF5m3bt3b/n6669l79690qBBA3nyySfluuuuM1cvAQAA+P3u1+aJxYqZS7G1NyYzM9PaVgEAAFgdZHSna50n065dO6ldu7bs2rVL3nzzTUlISGAdGQAA4L+TfXUIaf78+WZujF6KrYGmcuXK+ds6AAAAK4LM22+/LVFRUWZjyLVr15rblSxatMjVlwQAACiYINOnTx8zJwYAAMCWC+IBAAAUiquWAAAAfI0gAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbMuvg8wLL7wgAQEBTre6dev6ulkAAMBPlBA/16BBA1mxYkXO5yVK+H2TAQBAAfH7VKDBJSIiwtfNAAAAfsjvg8zPP/8s1apVk6CgIImJiZHY2FiJiorK8/yMjAxzy5aammo+ZmVlmRs8p/VzOBzU0UvU0TrU0jrU0hrU0Tqu1tCvg0zTpk1l9uzZUqdOHUlMTJSJEydK8+bNZffu3RIcHHzF52jQ0fMulZycLOnp6QXQ6sL9S5WSkmL+SIsV8+vpVX6NOlqHWlqHWlqDOlonLS3NpfMCHFptmzhz5oxER0fL1KlTpX///i73yERGRsrp06clJCSkAFtbOP9ANRCGhYXxB+oF6mgdamkdamkN6mgdff+uWLGiCYZXe//26x6ZS1WoUEFq164t+/fvz/OcwMBAc7uU/kLxS+U9vXKMWnqPOlqHWlqHWlqDOlrD1frZqspnz56VAwcOSNWqVX3dFAAA4Af8OsiMGjVK1q5dK7/++qts3LhR7rvvPilevLj07t3b100DAAB+wK+Hlg4fPmxCy8mTJ81445133imbN2829wEAAPw6yMyfP9/XTQAAAH7Mr4eWAAAAroYgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbIsgAwAAbKuErxsAAPCtzMws2b3uN0lOTpKwsPPSsEW0FC/O/3PdRR19wxZBZvr06fLPf/5Tjh07Jo0bN5Y33nhDbr/9dl83CwBsb+OiH2Tm0Hg5eTRNIpsEy6GtaRJaLVgGxrWXZt3q+bp5tkEdfcfvo+KCBQtkxIgR8vzzz8u2bdtMkGnfvr0kJSX5umkAYPs338k9FsqJw2lOx08cSTPH9XFcG3X0Lb8PMlOnTpUBAwZIv379pH79+vL2229LmTJl5P333/d10wDA1sMg2oMgjis8+OexmcO+Muchb9TR9/x6aOnChQuydetWGTduXM6xYsWKSdu2bWXTpk1XfE5GRoa5ZUtNTTUfs7KyzA2e0/o5HA7q6CXqaB1q6Tmdy6HDIAF//ndWPwYE/OdjtpNHUs15jVpG+6yd/o465h9X/679OsicOHFCMjMzpUqVKk7H9fMff/zxis+JjY2ViRMnXnY8OTlZ0tPT862tReWXKiUlxbxxaKCEZ6ijdail53RCqs7lyKZvvpVrlTb3HQ7n85KS/nMcl6OO+SctzXmozpZBxhPae6NzanL3yERGRkpYWJiEhIT4tG2F4U0jICDA1JI3Dc9RR+tQS8/pVTU6ITVbdg/CoW1p4sj1H+GwsHAJDw/3QQvtgTrmn6CgIPsHmcqVK0vx4sXl+PHjTsf184iIiCs+JzAw0Nwupf/I8Q+d9/RNg1p6jzpah1p6Ri8N1qtqdEJq9lwO7UHQN1/zBqw9CzVCzHnUNm/UMf+4Wi+/rmqpUqWkSZMmsnLlSqf/gennMTExPm0bANiZrm+ilwYbAZc8+OfnA6fdzToo10Adfc/vK6vDRO+++67MmTNHfvjhB3niiSfk3Llz5iomAIDndH2T8Qt7SOXq/53jobQHQY+z/olrqKNv+fXQkurZs6eZqPvcc8+ZBfFuvvlmWb58+WUTgAEA7tM32aZd6uRakTacFWk9QB19J8Ch0/0LMZ3sW758eXNlA5N9vaPDeroQoU5YY6zXc9TROtTSOtTSGtSx4N+/qTIAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtv9+iwFvZCxfrCoHwfsXKtLQ0s7U6K1Z6jjpah1pah1pagzpaJ/t9+1obEBT6IKO/UCoyMtLXTQEAAB68j+tWBUV2ryVNx0ePHpXg4GAJCLh0j3W4m441EB46dIh9q7xAHa1DLa1DLa1BHa2j8URDTLVq1a7au1Xoe2T0m69Ro4avm1Go6B8nf6Deo47WoZbWoZbWoI7WuFpPTDYG8AAAgG0RZAAAgG0RZOCywMBAef75581HeI46WodaWodaWoM6FrxCP9kXAAAUXvTIAAAA2yLIAAAA2yLIAAAA2yLIAAAA2yLI4JrWrVsnnTp1Mqsr6urIS5Ys8XWTbCk2NlZuu+02s8p0eHi4dO3aVfbt2+frZtnSjBkz5KabbspZdCwmJkaWLVvm62bZ3pQpU8zf+LBhw3zdFNt54YUXTO1y3+rWrevrZhUJBBlc07lz56Rx48Yyffp0XzfF1tauXSuDBg2SzZs3y9dffy0XL16Uu+++29QX7tHVuvVNd+vWrbJlyxZp3bq1dOnSRfbs2ePrptnWd999J++8844JiPBMgwYNJDExMee2YcMGXzepSCj0WxTAex06dDA3eGf58uVOn8+ePdv0zOibcYsWLXzWLjvSHsLcJk2aZHppNCTqmwncc/bsWXnwwQfl3XfflZdfftnXzbGtEiVKSEREhK+bUeTQIwP4SEpKivlYqVIlXzfF1jIzM2X+/PmmZ0uHmOA+7Sm89957pW3btr5uiq39/PPPZgj++uuvN8EwISHB100qEuiRAXy0K7vOQ7jjjjukYcOGvm6OLe3atcsEl/T0dClXrpwsXrxY6tev7+tm2Y6GwG3btpmhJXiuadOmppe1Tp06Zlhp4sSJ0rx5c9m9e7eZF4f8Q5ABfPQ/YP0HjjF0z+kbxo4dO0zP1sKFC6Vv375mHhJhxnWHDh2SoUOHmjlbQUFBvm6OreUeftd5RhpsoqOj5ZNPPpH+/fv7tG2FHUEGKGBPPfWUfPHFF+ZqMJ20Cs+UKlVKatWqZe43adLE9CjExcWZCatwjc7PSkpKkr/85S9OQ3X6u/nmm29KRkaGFC9e3KdttKsKFSpI7dq1Zf/+/b5uSqFHkAEKiG5rNnjwYDMEsmbNGqlZs6avm1Tohuv0jReua9OmjRmiy61fv37msuGxY8cSYrycQH3gwAF5+OGHfd2UQo8gA5f+IHP/r+LgwYOmS18nqUZFRfm0bXYbTpo7d678+9//NmPmx44dM8fLly8vpUuX9nXzbGXcuHGmK19//9LS0kxdNRzGx8f7umm2or+Hl87RKlu2rISGhjJ3y02jRo0yV9PpcNLRo0fNDtgaBHv37u3rphV6BBlck67Tcdddd+V8PmLECPNR5yTo5Da4Ri8PVq1atXI6PmvWLHnkkUd81Cp70uGQPn36mEmVGgR1ToKGmHbt2vm6aSiiDh8+bELLyZMnJSwsTO68806zHIDeR/4KcGh/NwAAgA2xjgwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAALAtggwAywQEBMiSJUt83QwARQhBBijkdPsDDRh6K1mypNmscsyYMZKenu7rphXqmnft2tXXzQCKBPZaAoqAe+65x+zpdPHiRdm6davZJ0uDzT/+8Q8pjPT71NBmdxcuXJBSpUr5uhmAX6NHBigCAgMDJSIiQiIjI01PQdu2beXrr7/OeVw3utMN76pXry5lypSRRo0aybx585xeQze7HDJkiOnN0Z3P9fVeeOGFq35d3QG4atWqsnPnzis+rs+/+eab5Z133jFt0699//33S0pKSs453333ndkMsnLlymaDyJYtW8q2bducXkdDmW7K2blzZ7N786RJkyQzM1P69+9veqB0d/E6depIXFzcFXtOJk+eLFWqVJEKFSrIiy++KH/88YeMHj3afJ81atQwITC3Q4cOmXbq+XpOly5d5Ndff835nubMmWN2Oc/uCdOdua/1vNzt0fZXq1bNtBnA1RFkgCJm9+7dsnHjRqf/6eswU5MmTeTLL780jw8cOFAefvhh+b//+z+n5+obtAaFb7/9Vl555RXzpp87EGXTvWgHDx4sH3zwgaxfv97sTp2X/fv3yyeffCJLly6V5cuXy/bt2+XJJ5/MeTwtLc30IG3YsMHsJnzjjTdKx44dzfHcNEDcd999smvXLnn00UclKyvLhJBPP/1U9u7dK88995yMHz/efK3cVq1aJUePHpV169bJ1KlTTfj629/+JhUrVjTf5+OPPy6PPfaY2d04u7enffv2EhwcbL63b775RsqVK2d6vbQHZdSoUSas6Oe6O7femjVrds3nZVu5cqXs27fP1PWLL75w62cLFEm6+zWAwqtv376O4sWLO8qWLesIDAzU3e4dxYoVcyxcuPCqz7v33nsdI0eOzPm8ZcuWjjvvvNPpnNtuu80xduzYnM/1tT/99FPHAw884KhXr57j8OHDV/0azz//vGlb7vOWLVtm2peYmHjF52RmZjqCg4MdS5cudfq6w4YNc1zLoEGDHN27d3eqTXR0tHnNbHXq1HE0b9485/M//vjD1G7evHnm8w8//NCck5WVlXNORkaGo3Tp0o74+Pic1+3SpYvT13b1eVWqVDHHAbiGOTJAEXDXXXeZoZdz587Ja6+9JiVKlJDu3bvnPK7DMDq8or0VR44cMT0EGRkZZqgnt0t7VnTYKCkpyenY8OHDzVCW9p7ocNC1REVFmSGtbDExMaY3RXsldPjq+PHj8uyzz5rhGf1a2tbff/9dEhISnF7n1ltvvey1p0+fLu+//7459/z58+b70qGs3Bo0aCDFiv23c1qHmBo2bJjzefHixSU0NDTn+/z+++9NL5L2rOSmvVoHDhzI8/t09Xk6rMe8GMB1BBmgCNDhoFq1apn7+sbeuHFjee+998wcEvXPf/7TzB+ZNm2aeSPV84cNG+Y05KEunUCr8z80dOSm81l0fk18fLw8+OCDXrddh5V0Do+2Lzo62oQkDTuXtk3bnNv8+fPNMM+rr75qztcAod+nDhdd63u62vd59uxZMwz38ccfX9bWsLCwPL8PV5936fcB4OoIMkARo70POldkxIgR8sADD5iJsDpfQyeePvTQQ+YcfdP+6aefpH79+m6/vk647dSpk3lt7c3o1avXVc/X3hKdo6KTW5X25Ggbsye6atveeustMy8me8LsiRMnrtkOfZ7OTck93+ZqPSau+stf/iILFiyQ8PBwCQkJueI52qOiPUfuPg+A+5jsCxRBf//7303I0KEXpRNodXKpTgL+4YcfzORWHdLxlE66/fDDD6Vfv36ycOHCq54bFBRkel106EUnweqVUTpZVoeVstumr6Xt0t4U7eXR8HUt+rwtW7aYniENZRMmTDBXQHlLv74OmWnw0/YePHjQDHtpu7MnBF933XXmSi0dHtPQpRN9XXkeAPcRZIAiSOfIPPXUU+bKI503o3NQtMdAr6rRy6w1RHi7oFuPHj3MVU569dOiRYvyPE+HvLp162Z6XO6++24zD0d7YLLpENjp06dN+/S19I1fezWuRcOYvm7Pnj2ladOmZngqd++Mp3TekF7hpHN79PXr1atnhuh0rkt2T8uAAQNMj5LO29FhI+0dcuV5ANwXoDN+PXgeAHhNL5nWLQ127Njh66YAsCl6ZAAAgG0RZAAAgG0xtAQAAGyLHhkAAGBbBBkAAGBbBBkAAGBbBBkAAGBbBBkAAGBbBBkAAGBbBBkAAGBbBBkAACB29f88IuBjJEAiPQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "locs_results_summary = simulation.evaluate_localization_for_each_rank(\n", " subject=subject,\n", " subjects_dir=subjects_dir,\n", " n_sources=len(labels_to_use_poststimuli)*N_VERTICES_PER_ACTIVITY_LABEL,\n", " R=data_cov.data,\n", " N=noise_cov.data,\n", " forward=fwd,\n", " true_vertices=poststimuli_vertices,\n", " localizer_to_use=\"mai_mvp\",\n", " plot_sum_error_by_rank=False,\n", " plot_correct_sources_by_rank=True,\n", " show_progress=False\n", ")" ] }, { "cell_type": "markdown", "id": "7a0920d216f0b102", "metadata": { "collapsed": false }, "source": [ "### Standard LCMV Beamformer (for comparison)\n", "\n", "For comparison with the rank-reduced MVPURE approach, we can also apply MNE's standard LCMV beamformer using ``mne.beamformer.make_lcmv``. Strongest sources are identified based on peak NAI value." ] }, { "cell_type": "code", "execution_count": 16, "id": "4d6a17674ec75782", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:16:50.848905Z", "start_time": "2026-06-27T19:16:50.659548Z" }, "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Computing rank from covariance with rank=None\n", " Using tolerance 6e-13 (2.2e-16 eps * 128 dim * 21 max singular value)\n", " Estimated rank (eeg): 127\n", " EEG: rank 127 computed from 128 data channels with 1 projector\n", "Computing rank from covariance with rank=None\n", " Using tolerance 3.3e-13 (2.2e-16 eps * 128 dim * 11 max singular value)\n", " Estimated rank (eeg): 127\n", " EEG: rank 127 computed from 128 data channels with 1 projector\n", "Making LCMV beamformer with rank {'eeg': 127}\n", "Computing inverse operator with 128 channels.\n", " 128 out of 128 channels remain after picking\n", "Selected 128 channels\n", "Whitening the forward solution.\n", " Created an SSP operator (subspace dimension = 1)\n", "Computing rank from covariance with rank={'eeg': 127}\n", " Setting small EEG eigenvalues to zero (without PCA)\n", "Creating the source covariance matrix\n", "Adjusting source covariance matrix.\n", "Computing beamformer filters for 5124 sources\n", "Filter computation complete\n" ] } ], "source": [ "(\n", " lcmv_top_vertices,\n", " lcmv_n_correctly_localized,\n", " lcmv_error_info\n", ") = simulation.compare_with_strongest_sources_lcmv(\n", " signal=sim_evoked,\n", " n_sources=len(labels_to_use_poststimuli)*N_VERTICES_PER_ACTIVITY_LABEL,\n", " true_vertices=poststimuli_vertices,\n", " forward=fwd,\n", " R=data_cov,\n", " N=noise_cov,\n", " reg=0.05,\n", " pick_ori=None\n", ")" ] }, { "cell_type": "code", "execution_count": 17, "id": "9dcbe97ea560b663", "metadata": { "ExecuteTime": { "end_time": "2026-06-27T19:16:53.174915Z", "start_time": "2026-06-27T19:16:53.161942Z" }, "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of correctly localized sources using LCMV-NAI: 1\n", "Mean distance error using LCMV-NAI: 58.1 mm\n" ] } ], "source": [ "print(f\"Number of correctly localized sources using LCMV-NAI: {lcmv_n_correctly_localized}\")\n", "print(f\"Mean distance error using LCMV-NAI: {np.round(lcmv_error_info['mean_error'], 2)} mm\")" ] }, { "cell_type": "markdown", "id": "7375d0d8872131a1", "metadata": { "collapsed": false }, "source": [ "This tutorial has shown a complete pipeline for:\n", "\n", "1. **Simulating realistic EEG data** with known ground truth sources\n", "2. **Controlling experimental parameters** (SNR, coupling, spatial distribution)\n", "3. **Preprocessing** sensor data using standard pipelines\n", "4. **Applying localization algorithm** with the MVPURE framework\n", "5. **Quantitatively evaluating** localization performance.\n", "\n", "To further explore source localization: test different SNRs, coupling strengths, or source configurations." ] } ], "metadata": { "kernelspec": { "display_name": "mvpure-tools-venv", "language": "python", "name": "mvpure-tools-venv" }, "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.10.18" } }, "nbformat": 4, "nbformat_minor": 5 }