FNIRS_SVM_Group_analysis.ipynb

{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "EhPm_3tHt4mQ",
        "outputId": "b926a420-d6ad-44b7-951f-52d0d162065f"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Mounted at /content/drive\n"
          ]
        }
      ],
      "source": [
        "## for google colab\n",
        "\n",
        "# If encounter google drive issue, remount\n",
        "# !rm -rf /content/drive\n",
        "\n",
        "from google.colab import drive\n",
        "drive.mount('/content/drive')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "11eC2IiTuLL6"
      },
      "outputs": [],
      "source": [
        "%%capture\n",
        "!pip install mne ipywidgets ipyevents nibabel trame trame-vtk trame-vuetify pyvista nilearn PyWavelets seaborn statsmodels mne_nirs snirf optuna sktime scikit-learn numpy pandas optuna"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q8463p6_t4mT"
      },
      "source": [
        "## 1. Load Library"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "3bjepRqht4mU",
        "outputId": "500660e5-1672-4230-b823-f7698e094d2e"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "INFO:pysnirf2:Library loaded by process 327\n"
          ]
        }
      ],
      "source": [
        "# Import common libraries\n",
        "from collections import defaultdict\n",
        "from copy import deepcopy\n",
        "from itertools import compress\n",
        "from pprint import pprint\n",
        "\n",
        "# Import Plotting Library\n",
        "import matplotlib.pyplot as plt\n",
        "import pandas as pd\n",
        "import seaborn as sns\n",
        "\n",
        "# Import StatsModels\n",
        "import statsmodels.formula.api as smf\n",
        "from mne import Epochs, events_from_annotations, set_log_level\n",
        "from mne.preprocessing.nirs import (\n",
        "    beer_lambert_law,\n",
        "    optical_density,\n",
        "    scalp_coupling_index,\n",
        "    temporal_derivative_distribution_repair,\n",
        ")\n",
        "\n",
        "# Import MNE processing\n",
        "from mne.viz import plot_compare_evokeds\n",
        "\n",
        "\n",
        "# Import MNE-NIRS processing\n",
        "from mne_nirs.channels import get_long_channels, picks_pair_to_idx\n",
        "from mne_nirs.datasets import fnirs_motor_group\n",
        "from mne_nirs.signal_enhancement import enhance_negative_correlation\n",
        "\n",
        "# Set general parameters\n",
        "set_log_level(\"WARNING\")  # Don't show info, as it is repetitive for many subjects\n",
        "from itertools import compress\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import os\n",
        "import mne\n",
        "import snirf\n",
        "\n",
        "from scipy import signal"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jdmfGB71t4mU"
      },
      "source": [
        "## 2. Define function"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ThZGa0b4t4mV"
      },
      "source": [
        "#### 2.0 Remove close channels"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Zkbin2zlt4mV"
      },
      "outputs": [],
      "source": [
        "def remove_close_channels(raw_intensity):\n",
        "# Remove channels that are too close to each other (less than 1cm)\n",
        "    picks = mne.pick_types(raw_intensity.info, meg=False, fnirs=True) # [0 ... 111]\n",
        "    dists = mne.preprocessing.nirs.source_detector_distances(\n",
        "        raw_intensity.info, picks=picks\n",
        "    )\n",
        "    raw_intensity.pick(picks[dists > 0.01]) # Number of channels with distance > 1 cm: 96\n",
        "\n",
        "    return raw_intensity"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YzFo5Ctxt4mW"
      },
      "source": [
        "#### 2.1 Annotation function"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "XHeqs7m6t4mW"
      },
      "outputs": [],
      "source": [
        "def process_annotations(raw_intensity):\n",
        "\n",
        "    # Define wanted stimuli\n",
        "    wanted_stimulus = [\"2\", \"3\"]\n",
        "    unwanted = np.nonzero(~np.isin(raw_intensity.annotations.description, wanted_stimulus))\n",
        "    # The ~ operator is used to negate the result of np.isin\n",
        "    raw_intensity.annotations.delete(unwanted)\n",
        "    raw_intensity.annotations.rename(\n",
        "        {\"2\": \"Reappraise\", \"3\": \"view negative\"}\n",
        "    )\n",
        "\n",
        "    return raw_intensity\n",
        "\n",
        "# Usage:\n",
        "# raw_intensity = process_annotations(raw_intensity)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0Eva31tat4mW"
      },
      "source": [
        "#### 2.2 Wavelet motion Correction function"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "EGTuhG-Ct4mW"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import pywt\n",
        "import mne\n",
        "\n",
        "def hmr_motion_correct_wavelet_standard(raw_od, iqr=1.5, turnon=1, wavename='db2'):\n",
        "    if iqr < 0 or turnon == 0:\n",
        "        return raw_od\n",
        "\n",
        "    data = raw_od.get_data()\n",
        "    dod_wavelet = np.zeros_like(data)\n",
        "    signal_length = data.shape[1]\n",
        "    n_levels = max(int(np.ceil(np.log2(signal_length))) - 4, 4) # decomposition level, L = 4 is the lowest wavelet scale used in the analysis\n",
        "\n",
        "    # Identify bad channels\n",
        "    bads = raw_od.info['bads']\n",
        "    bad_indices = [raw_od.ch_names.index(bad) for bad in bads if bad in raw_od.ch_names]\n",
        "\n",
        "    # Determine good channels\n",
        "    all_indices = set(range(data.shape[0]))\n",
        "    bad_indices_set = set(bad_indices)\n",
        "    good_indices = sorted(all_indices - bad_indices_set)\n",
        "\n",
        "    if not good_indices:\n",
        "        print(\"All channels are marked as bad. No motion correction applied.\")\n",
        "        return raw_od.copy()\n",
        "\n",
        "\n",
        "    for idx, channel_no in enumerate(good_indices):\n",
        "        # Remove DC component\n",
        "        dc_val = np.mean(data[channel_no, :])\n",
        "        data_no_dc = data[channel_no, :] - dc_val\n",
        "\n",
        "        # Normalize\n",
        "        med_dev = np.median(np.abs(data_no_dc))\n",
        "        if med_dev == 0:\n",
        "            norm_coef = 1\n",
        "            data_normalized = data_no_dc\n",
        "        else:\n",
        "            norm_coef = 1 / (1.4826 * med_dev)\n",
        "            data_normalized = data_no_dc * norm_coef\n",
        "\n",
        "        # Wavelet Transform\n",
        "        coeffs = pywt.wavedec(data_normalized, wavename, mode='periodization', level=n_levels)\n",
        "\n",
        "        # Thresholding (Artifact Removal)\n",
        "        for i in range(1, len(coeffs)):\n",
        "            c = coeffs[i]\n",
        "            q1, q3 = np.percentile(c, [25, 75])\n",
        "            iqr_val = q3 - q1\n",
        "            threshold = iqr * iqr_val\n",
        "            coeffs[i] = np.where((c > (q3 + threshold)) | (c < (q1 - threshold)), 0, c) # np.where(condition, value if true, value if false)\n",
        "\n",
        "        # Inverse Wavelet Transform\n",
        "        reconstructed = pywt.waverec(coeffs, wavename, mode='periodization')\n",
        "\n",
        "        # Denormalize and add DC component\n",
        "        reconstructed = reconstructed / norm_coef + dc_val\n",
        "\n",
        "        # Truncate to original signal length\n",
        "        dod_wavelet[channel_no, :] = reconstructed[:signal_length]\n",
        "\n",
        "    # For bad channels, retain the original data\n",
        "    if bad_indices:\n",
        "        dod_wavelet[bad_indices, :] = data[bad_indices, :]\n",
        "        print(f\"Retained original data for {len(bad_indices)} bad channels.\")\n",
        "\n",
        "    raw_od_corrected = mne.io.RawArray(dod_wavelet, raw_od.info)\n",
        "    return raw_od_corrected"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-qj5fBMxt4mX"
      },
      "source": [
        "#### 2.3 hmr_Bandpass_Filter function\n",
        "(Auxiliary --> accelerometer/gyroscope data).\n",
        "Both optical density and auxiliary data uses the same bandpass filter function"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Xjys4sr3t4mX"
      },
      "outputs": [],
      "source": [
        "def hmr_Bandpass_Filter(data, hpf=0.01, lpf=0.4):\n",
        "    fs = data.info[\"sfreq\"] # sampling frequency = 5.09 Hz\n",
        "    if hpf > (fs/2) or lpf > (fs/2): # if bandpass filter > nyquist frequency\n",
        "        print(\"Warning: Bandpass filter is larger than the Nyquist frequency\")\n",
        "        return data\n",
        "\n",
        "    # Get the data array\n",
        "    x = data.get_data()\n",
        "    filtered_data = x.copy()\n",
        "\n",
        "    # low-pass filter\n",
        "    lpf_norm = lpf / (fs/2) #  The normalized lpf_norm = 0.157 indicates that the low-pass cutoff is at 15.7% of the Nyquist frequency.\n",
        "    if lpf_norm > 0:\n",
        "        # Create filter\n",
        "        filter_order = 3\n",
        "        # Design an 3rd-order digital Butterworth filter and return the filter coefficients.\n",
        "        b, a = signal.butter(filter_order, lpf_norm, 'low')  # butter(N, Wn, btype='low')\n",
        "            # N: The order of the filter\n",
        "            # For digital filters, if fs is not specified, Wn units are normalized from 0 to 1, where 1 is the Nyquist frequency (Wn is thus in half cycles / sample and defined as 2*critical frequencies / fs). If fs is specified, Wn is in the same units as fs.\n",
        "            # btype {‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’}, optional, The type of filter. Default is ‘lowpass’.\n",
        "\n",
        "        for ch in range(x.shape[0]):\n",
        "            filtered_data[ch, :] = signal.filtfilt(b, a, filtered_data[ch, :]) # b : The numerator coefficient vector; a : The denominator coefficient vector\n",
        "\n",
        "    # High-pass filter\n",
        "    hpf_norm = hpf / (fs/2)\n",
        "    if hpf_norm > 0:\n",
        "        filter_order_h = 5\n",
        "        b_h, a_h = signal.butter(filter_order_h,hpf_norm, 'high')\n",
        "        for ch in range(x.shape[0]):\n",
        "            filtered_data[ch, :] = signal.filtfilt(b_h, a_h, filtered_data[ch, :])\n",
        "\n",
        "    data_filtered = data.copy()\n",
        "    data_filtered._data = filtered_data\n",
        "\n",
        "    return data_filtered\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YP7c-yE4t4mX"
      },
      "source": [
        "#### 2.4 Define MotionCorrectCbsi function\n",
        "\n",
        "This function assumes that the true psychological signals of HbO and HbR are negatively correlated\n",
        "\n",
        "Page 5 of https://doi.org/10.3390/s23083979"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "g-7VC372t4mY"
      },
      "outputs": [],
      "source": [
        "def hmrR_MotionCorrectCbsi(raw_haemo):\n",
        "\n",
        "    hbo_channels = [ch for ch in raw_haemo.ch_names if 'hbo' in ch.lower()]\n",
        "    hbr_channels = [ch for ch in raw_haemo.ch_names if 'hbr' in ch.lower()]\n",
        "\n",
        "    hbo_indices = [idx for idx, ch in enumerate(raw_haemo.ch_names) if 'hbo' in ch.lower()] # [0-47]\n",
        "    hbr_indices = [idx for idx, ch in enumerate(raw_haemo.ch_names) if 'hbr' in ch.lower()] #[48-95]\n",
        "\n",
        "    haemo_data = raw_haemo.get_data() # (96, 2396)\n",
        "    n_timepoints = haemo_data.shape[1] # 2396\n",
        "    n_channels = len(hbo_channels) # 48\n",
        "\n",
        "    data_dc = np.zeros((n_timepoints, 3 ,n_channels)) # (2396, 3, 48)\n",
        "    data_dc[:,0,:] = haemo_data[hbo_indices,:].T\n",
        "    data_dc[:,1,:] = haemo_data[hbr_indices,:].T\n",
        "\n",
        "    # avoid changing original data\n",
        "    dc = data_dc.copy()\n",
        "\n",
        "    # Copy data for correction\n",
        "    dcCbsi = dc.copy()\n",
        "\n",
        "    for idx_ch in range(n_channels):\n",
        "\n",
        "        # Mean-center the HbO and HbR signals\n",
        "        dc_oxy = dc[:, 0, idx_ch] - np.mean(dc[:, 0, idx_ch])\n",
        "        dc_deoxy = dc[:, 1, idx_ch] - np.mean(dc[:, 1, idx_ch])\n",
        "\n",
        "        # Calculate standard deviations\n",
        "        sd_oxy = np.std(dc_oxy, ddof=0) #Delta Degrees of Freedom = 0: divides by n\n",
        "        sd_deoxy = np.std(dc_deoxy, ddof=0)\n",
        "\n",
        "        # Compute alpha, avoiding division by zero\n",
        "        if sd_deoxy != 0:\n",
        "            alfa = sd_oxy / sd_deoxy\n",
        "        else:\n",
        "            alfa = 0\n",
        "\n",
        "        # Apply CBSI correction\n",
        "        dcCbsi[:, 0, idx_ch] = 0.5 * (dc_oxy - alfa * dc_deoxy)\n",
        "        if alfa != 0:\n",
        "            dcCbsi[:, 1, idx_ch] = - (1 / alfa) * dcCbsi[:, 0, idx_ch]\n",
        "        else:\n",
        "            dcCbsi[:, 1, idx_ch] = 0  # Handle division by zero\n",
        "\n",
        "        dcCbsi[:, 2, idx_ch] = dcCbsi[:, 0, idx_ch] + dcCbsi[:, 1, idx_ch]\n",
        "\n",
        "    raw_copy = raw_haemo.copy()\n",
        "\n",
        "    # Step 4: Extract and replace HbO and HbR data\n",
        "    corrected_hbo = dcCbsi[:, 0, :].T  # Shape: (n_channels, n_timepoints)\n",
        "    corrected_hbr = dcCbsi[:, 1, :].T  # Shape: (n_channels, n_timepoints)\n",
        "\n",
        "    # Replace data in-place\n",
        "    raw_copy._data[hbo_indices, :] = corrected_hbo\n",
        "    raw_copy._data[hbr_indices, :] = corrected_hbr\n",
        "\n",
        "    # Step 5: Add HbT channels\n",
        "    corrected_hbt = dcCbsi[:, 2, :].T  # Shape: (n_channels, n_timepoints)\n",
        "\n",
        "    # Define HbT channel names\n",
        "    hbt_channels = [ch.replace('hbo', 'hbt') for ch in hbo_channels]\n",
        "\n",
        "    # Define channel types for HbT\n",
        "    hbt_channel_types = ['misc'] * n_channels  # 'misc' since 'hbt' is not standard in MNE\n",
        "\n",
        "    # Create Info object for HbT channels\n",
        "    hbt_info = mne.create_info(ch_names=hbt_channels, sfreq=raw_copy.info['sfreq'], ch_types=hbt_channel_types)\n",
        "\n",
        "    # Create RawArray for HbT channels\n",
        "    raw_hbt = mne.io.RawArray(corrected_hbt, hbt_info)\n",
        "\n",
        "    # Add HbT channels to the Raw copy\n",
        "    raw_copy.add_channels([raw_hbt], force_update_info=True)\n",
        "\n",
        "    return raw_copy"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yM9IFOo3t4mY"
      },
      "source": [
        "#### 2.5 Define get ROI channel function\n",
        "\n",
        "take away 'hbr' data to reduce redundancy"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "YlKvv7j-t4mY"
      },
      "outputs": [],
      "source": [
        "def get_ROI_channels(raw_copy):\n",
        "    # Get channel names from your fNIRS dataset\n",
        "    channel_names = raw_copy.ch_names\n",
        "\n",
        "    # Define the selection criteria\n",
        "    selection_criteria = ['S2', 'S4', 'S5', 'S6']\n",
        "\n",
        "    # Create a list to store selected channel names\n",
        "    selected_channels = []\n",
        "\n",
        "    for channel in channel_names:\n",
        "        if any(criterion in channel for criterion in selection_criteria) and ('hbo' in channel or 'hbr' in channel or 'hbt' in channel): #HbO channels: 760nm HbR channels: 850nm\n",
        "            selected_channels.append(channel)\n",
        "\n",
        "    # Now, get the indices of these selected channels\n",
        "    picks = [channel_names.index(name) for name in selected_channels]\n",
        "\n",
        "    return picks"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7LpOWhJht4mY"
      },
      "source": [
        "#### 2.6 Define get evoked metadata for Reappraise / view negative"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "jGOJhHu8t4mY"
      },
      "outputs": [],
      "source": [
        "def get_evoked_metadata(raw_copy, events, event_dict, picks):\n",
        "    # Rejection criteria for bad epochs\n",
        "    reject_criteria = dict(hbo=70e-6)\n",
        "\n",
        "    epochs = mne.Epochs(raw_copy,\n",
        "                       events,\n",
        "                       event_id=event_dict,\n",
        "                       tmin=-2,\n",
        "                       tmax=20,\n",
        "                       picks=picks,\n",
        "                       baseline=(-2, 0),\n",
        "                       reject=reject_criteria,\n",
        "                       preload=True,\n",
        "                       verbose=False)\n",
        "\n",
        "    # Get rejected epochs information\n",
        "    drop_log = epochs.drop_log\n",
        "    dropped_indices = [i for i, drops in enumerate(drop_log) if drops]\n",
        "\n",
        "    print(f\"\\nNumber of rejected epochs: {len(dropped_indices)}\")\n",
        "    print(\"\\n\")\n",
        "    print(\"\\nRejected epochs event types:\")\n",
        "    for i in dropped_indices:\n",
        "        event_type = events[i][2]\n",
        "        print(f\"Epoch {i}: Event type {event_type}\")\n",
        "\n",
        "\n",
        "    # Separate epochs by condition and channel type\n",
        "    epochs_dict = {\n",
        "        \"Reappraise\": {\n",
        "            \"HbO\": np.mean(epochs[\"Reappraise\"].average(picks=\"hbo\").data, axis = 1),  # average over time\n",
        "        },\n",
        "        \"view_negative\": {\n",
        "            \"HbO\": np.mean(epochs[\"view negative\"].average(picks=\"hbo\").data, axis=1),  # average over time\n",
        "        }\n",
        "    }\n",
        "\n",
        "    # Get channel names for each type\n",
        "    hbo_ch_names = [ch for ch in epochs.ch_names if 'hbo' in ch.lower()]\n",
        "\n",
        "\n",
        "    # Remove frequency encoding from channel names\n",
        "    hbo_ch_names = [ch[:-4] for ch in hbo_ch_names]\n",
        "\n",
        "\n",
        "    # Create metadata DataFrames\n",
        "    metadata_df_reappraise = pd.DataFrame({\n",
        "        'epoch': range(len(epochs[\"Reappraise\"])),\n",
        "        'condition': 'Reappraise',\n",
        "        'times': [epochs.times] * len(epochs[\"Reappraise\"]),\n",
        "        'hbo_ch_names': [hbo_ch_names] * len(epochs[\"Reappraise\"]),\n",
        "    })\n",
        "\n",
        "    metadata_df_view_negative = pd.DataFrame({\n",
        "        'epoch': range(len(epochs[\"view negative\"])),\n",
        "        'condition': 'view_negative',\n",
        "        'times': [epochs.times] * len(epochs[\"view negative\"]),\n",
        "        'hbo_ch_names': [hbo_ch_names] * len(epochs[\"view negative\"]),\n",
        "    })\n",
        "\n",
        "    return epochs_dict, metadata_df_reappraise, metadata_df_view_negative, epochs.times"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sQhf4AI9t4mY"
      },
      "source": [
        "#### 2.7 Define interpolation function"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "_23Htii7t4mY"
      },
      "outputs": [],
      "source": [
        "# identifies bad fNIRS channels via scalp coupling index, interpolates them\n",
        "def interpolate_bad_channels(raw_od, sci_value):\n",
        "    # Make a copy for preprocessing\n",
        "    raw_od_copy = raw_od.copy()\n",
        "\n",
        "    # Set montage\n",
        "    montage = raw_od_copy.get_montage()\n",
        "    raw_od_copy.set_montage(montage)\n",
        "\n",
        "    # Calculate scalp coupling index\n",
        "    sci = mne.preprocessing.nirs.scalp_coupling_index(raw_od_copy)\n",
        "\n",
        "    # Mark bad channels based on SCI\n",
        "    raw_od_copy.info[\"bads\"] = list(compress(raw_od_copy.ch_names, sci < sci_value))\n",
        "    # print(f\"Bad channels before interpolation: {raw_od_copy.info['bads']}\") #####################################################\n",
        "    # Interpolate bad channels\n",
        "    raw_od_copy.interpolate_bads(\n",
        "        reset_bads=False,      # Keep bad channels in the bads list after interpolation\n",
        "        mode='accurate',       # Use high-quality interpolation\n",
        "        method=\"nearest\"       # Specify interpolation method for fNIRS channels\n",
        "    )\n",
        "    # print(f\"Bad channels after interpolation: {raw_od.info['bads']}\") #########################################################\n",
        "    # Create a RawArray with processed data and original info\n",
        "    processed_data = raw_od_copy.get_data()\n",
        "    original_info = raw_od.info.copy()\n",
        "\n",
        "    # Create new Raw object\n",
        "    raw_od = mne.io.RawArray(processed_data, original_info)\n",
        "\n",
        "    return raw_od"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "N9zYNQmGt4mZ"
      },
      "source": [
        "## 3. Subject Level Analysis\n",
        "\n",
        "Each person has 12 metadata, every metadata is a 2D matrix of shape {39 channels [HbO, HbR, HbT], 113 timepoints [-2 to 20s]}"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "FpLBaWAft4mZ"
      },
      "outputs": [],
      "source": [
        "import os\n",
        "import mne\n",
        "import numpy as np\n",
        "from glob import glob\n",
        "import mne_nirs\n",
        "import logging\n",
        "\n",
        "import warnings\n",
        "\n",
        "# Suppress specific warnings\n",
        "warnings.filterwarnings(\"ignore\", message=\"Fiducial point nasion not found, assuming identity unknown to head transformation\", category=RuntimeWarning)\n",
        "warnings.filterwarnings(\"ignore\", message=\"No bad channels to interpolate. Doing nothing...\", category=RuntimeWarning)\n",
        "\n",
        "base_dir = r\"/content/drive/MyDrive/fNIRS data/neurofeedback\"\n",
        "participant_id = 1\n",
        "\n",
        "def subject_wise_analysis_1_to_6(base_dir, participant_id, ROI=True):\n",
        "    participant_dir = os.path.join(base_dir, f\"p{participant_id:02d}\")\n",
        "\n",
        "    # Find all directories that have at least one underscore in their name\n",
        "    date_dirs = glob(os.path.join(participant_dir, \"????-??-??_???\"))\n",
        "\n",
        "    if not date_dirs:\n",
        "        print(f\"No data directories found for participant {participant_id}\")\n",
        "        return None, None, None\n",
        "\n",
        "    # Initialize dictionary to store epochs\n",
        "    all_epochs_dict = {\n",
        "\n",
        "            \"Reappraise\": {\n",
        "                \"HbO\": []\n",
        "            },\n",
        "            \"view_negative\": {\n",
        "                \"HbO\": []\n",
        "            }\n",
        "\n",
        "    }\n",
        "\n",
        "    run_numbers = [run[-3:] for run in date_dirs]\n",
        "    # print(\"All run numbers:\", run_numbers)\n",
        "\n",
        "    # Get only run 001\n",
        "    sorted_runs = sorted(run_numbers)\n",
        "    first_run = [sorted_runs[0]]  # Get only first run (001)\n",
        "\n",
        "    selected_runs = [run for run in date_dirs if run[-3:] in first_run]\n",
        "\n",
        "    for run in selected_runs:\n",
        "        run_num = run[-3:]\n",
        "        print(\"\\nRun\", run_num)\n",
        "        snirf_files = glob(os.path.join(run, \"*.snirf\"))\n",
        "\n",
        "        if not snirf_files:\n",
        "            print(f\"No SNIRF file found in {run}\")\n",
        "            continue\n",
        "\n",
        "        fnirs_snirf = snirf_files[0]\n",
        "\n",
        "        # Load raw intensity info\n",
        "        raw_intensity = mne.io.read_raw_snirf(fnirs_snirf, preload=True)\n",
        "        raw_intensity.load_data()\n",
        "\n",
        "        # Process the data\n",
        "        raw_intensity = process_annotations(raw_intensity)\n",
        "        raw_intensity = remove_close_channels(raw_intensity)\n",
        "        raw_od = mne.preprocessing.nirs.optical_density(raw_intensity)\n",
        "        raw_od = interpolate_bad_channels(raw_od, sci_value=0.9)\n",
        "        raw_od_wavelet = hmr_motion_correct_wavelet_standard(raw_od, iqr=1.5, turnon=1, wavename='db2')\n",
        "        raw_od_wavelet_filtered = hmr_Bandpass_Filter(raw_od_wavelet, hpf=0.01, lpf=0.4)\n",
        "        raw_haemo = mne.preprocessing.nirs.beer_lambert_law(raw_od_wavelet_filtered, ppf=6.0)\n",
        "        raw_copy = hmrR_MotionCorrectCbsi(raw_haemo)\n",
        "\n",
        "        # Get channels\n",
        "        if ROI:\n",
        "            picks = get_ROI_channels(raw_copy)\n",
        "            hbo_ch_names = [ch for ch in np.array(raw_copy.ch_names)[picks] if 'hbo' in ch.lower()]\n",
        "        else:\n",
        "            picks = list(range(len(raw_copy.ch_names)))\n",
        "            hbo_ch_names = [ch for ch in raw_copy.ch_names if 'hbo' in ch.lower()]\n",
        "\n",
        "        # Get epochs\n",
        "        events, event_dict = mne.events_from_annotations(raw_intensity)\n",
        "        epochs_dict, metadata_df_reappraise, metadata_df_view_negative, epochs_times = get_evoked_metadata(\n",
        "            raw_copy, events, event_dict, picks)\n",
        "\n",
        "        # Store epochs in the dictionary\n",
        "        all_epochs_dict[\"Reappraise\"][\"HbO\"] = epochs_dict[\"Reappraise\"][\"HbO\"]\n",
        "        all_epochs_dict[\"view_negative\"][\"HbO\"] = epochs_dict[\"view_negative\"][\"HbO\"]\n",
        "\n",
        "    # Convert lists to numpy arrays and print shapes\n",
        "    for condition in all_epochs_dict:\n",
        "          all_epochs_dict[condition][\"HbO\"] = np.array(all_epochs_dict[condition][\"HbO\"])\n",
        "          array_shape = all_epochs_dict[condition][\"HbO\"].shape\n",
        "          print(f\"Participant {participant_id} - {condition}: {array_shape}\")\n",
        "\n",
        "    return all_epochs_dict, epochs_times, hbo_ch_names"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ffgV8f-Kt4mZ",
        "outputId": "d67ff7e5-7129-4bd5-90ed-2bfc29613ab7"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 1 - Reappraise: (13,)\n",
            "Participant 1 - view_negative: (13,)\n",
            "{'Reappraise': {'HbO': array([ 6.42470800e-08,  7.40716326e-08,  1.48861541e-08,  3.96452970e-08,\n",
            "       -2.53012675e-09,  3.93469066e-08,  2.12135430e-08,  4.93735923e-08,\n",
            "        3.75904438e-08,  1.22373408e-08,  4.97176675e-08,  4.24024541e-08,\n",
            "        1.52321513e-08])}, 'view_negative': {'HbO': array([-8.11709026e-08, -9.95629926e-08, -6.05969031e-08,  1.15628102e-08,\n",
            "        2.29510115e-08,  3.35584298e-08,  6.56851344e-09, -1.47226080e-08,\n",
            "       -1.04553863e-08, -2.29109197e-08,  2.75062621e-08, -7.46029532e-09,\n",
            "       -1.69376182e-08])}}\n"
          ]
        }
      ],
      "source": [
        "all_epochs_dict, epochs_time, hbo_ch_names =subject_wise_analysis_1_to_6(base_dir, participant_id, ROI=True)\n",
        "print(all_epochs_dict)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FJ9xtAWEt4ma"
      },
      "source": [
        "## 4. get all participants' data\n",
        "\n",
        "1. If a channel is bad in any run, it might be unreliable in other runs too\n",
        "2. For group analysis, we need consistent channel layouts across all runs\n",
        "3. So, we remove channels marked as bad in any run from all runs"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "_iwm2Jeht4ma",
        "collapsed": true,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "e72ecf0a-19d5-4c99-85fb-1597edfd37d7"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 1 - Reappraise: (13,)\n",
            "Participant 1 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 3 - Reappraise: (13,)\n",
            "Participant 3 - view_negative: (13,)\n",
            "\n",
            "Run 007\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 6 - Reappraise: (13,)\n",
            "Participant 6 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 7 - Reappraise: (13,)\n",
            "Participant 7 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 10 - Reappraise: (13,)\n",
            "Participant 10 - view_negative: (13,)\n",
            "\n",
            "Run 007\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 11 - Reappraise: (13,)\n",
            "Participant 11 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 13 - Reappraise: (13,)\n",
            "Participant 13 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 15 - Reappraise: (13,)\n",
            "Participant 15 - view_negative: (13,)\n",
            "\n",
            "Run 002\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 16 - Reappraise: (13,)\n",
            "Participant 16 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 17 - Reappraise: (13,)\n",
            "Participant 17 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 20 - Reappraise: (13,)\n",
            "Participant 20 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 21 - Reappraise: (13,)\n",
            "Participant 21 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 24 - Reappraise: (13,)\n",
            "Participant 24 - view_negative: (13,)\n",
            "\n",
            "Run 007\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 28 - Reappraise: (13,)\n",
            "Participant 28 - view_negative: (13,)\n",
            "\n",
            "Run 007\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 30 - Reappraise: (13,)\n",
            "Participant 30 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 31 - Reappraise: (13,)\n",
            "Participant 31 - view_negative: (13,)\n",
            "\n",
            "Run 006\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 33 - Reappraise: (13,)\n",
            "Participant 33 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 37 - Reappraise: (13,)\n",
            "Participant 37 - view_negative: (13,)\n",
            "\n",
            "Run 007\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 38 - Reappraise: (13,)\n",
            "Participant 38 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 40 - Reappraise: (13,)\n",
            "Participant 40 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 41 - Reappraise: (13,)\n",
            "Participant 41 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 43 - Reappraise: (13,)\n",
            "Participant 43 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 45 - Reappraise: (13,)\n",
            "Participant 45 - view_negative: (13,)\n",
            "\n",
            "Run 002\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 46 - Reappraise: (13,)\n",
            "Participant 46 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 48 - Reappraise: (13,)\n",
            "Participant 48 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 50 - Reappraise: (13,)\n",
            "Participant 50 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 55 - Reappraise: (13,)\n",
            "Participant 55 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 56 - Reappraise: (13,)\n",
            "Participant 56 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 2 - Reappraise: (13,)\n",
            "Participant 2 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 5 - Reappraise: (13,)\n",
            "Participant 5 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 9 - Reappraise: (13,)\n",
            "Participant 9 - view_negative: (13,)\n",
            "\n",
            "Run 003\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 12 - Reappraise: (13,)\n",
            "Participant 12 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 14 - Reappraise: (13,)\n",
            "Participant 14 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 18 - Reappraise: (13,)\n",
            "Participant 18 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 19 - Reappraise: (13,)\n",
            "Participant 19 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 22 - Reappraise: (13,)\n",
            "Participant 22 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 23 - Reappraise: (13,)\n",
            "Participant 23 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 25 - Reappraise: (13,)\n",
            "Participant 25 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 26 - Reappraise: (13,)\n",
            "Participant 26 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 27 - Reappraise: (13,)\n",
            "Participant 27 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 29 - Reappraise: (13,)\n",
            "Participant 29 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 34 - Reappraise: (13,)\n",
            "Participant 34 - view_negative: (13,)\n",
            "\n",
            "Run 007\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 35 - Reappraise: (13,)\n",
            "Participant 35 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 36 - Reappraise: (13,)\n",
            "Participant 36 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 39 - Reappraise: (13,)\n",
            "Participant 39 - view_negative: (13,)\n",
            "\n",
            "Run 007\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 42 - Reappraise: (13,)\n",
            "Participant 42 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 44 - Reappraise: (13,)\n",
            "Participant 44 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 47 - Reappraise: (13,)\n",
            "Participant 47 - view_negative: (13,)\n",
            "\n",
            "Run 007\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 49 - Reappraise: (13,)\n",
            "Participant 49 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 51 - Reappraise: (13,)\n",
            "Participant 51 - view_negative: (13,)\n",
            "\n",
            "Run 001\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 52 - Reappraise: (13,)\n",
            "Participant 52 - view_negative: (13,)\n",
            "\n",
            "Run 007\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 53 - Reappraise: (13,)\n",
            "Participant 53 - view_negative: (13,)\n",
            "\n",
            "Run 013\n",
            "\n",
            "Number of rejected epochs: 0\n",
            "\n",
            "\n",
            "\n",
            "Rejected epochs event types:\n",
            "Participant 54 - Reappraise: (13,)\n",
            "Participant 54 - view_negative: (13,)\n",
            "X shape: (106, 13)\n",
            "y shape: (106,)\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "from scipy import stats\n",
        "\n",
        "def process_fnirs_data():\n",
        "\n",
        "    # We reject participant's data if the negativity reaponse pattern happens to be \"invariant responding\", which is a clear indicator of low-quality data\n",
        "\n",
        "    base_dir = r\"/content/drive/MyDrive/fNIRS data/neurofeedback\" # google colab: r\"/content/drive/MyDrive/fNIRS data/neurofeedback\"\n",
        "    participant_ids= [1,3,6,7,10,11,13,15,16,17,20,21,24,28,30,31,33,37,38,40,41,43,45,46,48,50,55,56] # google colab: [1,3,6,7,8,10,11,13,15,16,17,20,21,24,28,30,31,33,37,38,40,41,43,45]\n",
        "      # reject p8 as 51/63 responses are 1s (79.7%), 13/64 responses are 2s\n",
        "    base_dir_sham = r\"/content/drive/MyDrive/fNIRS data/sham\"\n",
        "    participants_id_sham = [2,5,9,12,14,18,19,22,23,25,26,27,29,34,35,36,39,42,44,47,49,51,52,53,54]\n",
        "      # reject p4 as 56/63 responses are 1s (88.9%), 7/64 responses are 2s\n",
        "      # reject p32 as 60/64 responses are 1s (93.75%), 4/64 responses are 2s\n",
        "\n",
        "    # Initialize empty lists for features (X) and labels (y)\n",
        "    X = []  # Will store all HbO data\n",
        "    y = []  # Will store corresponding labels (0 for view_negative, 1 for Reappraise)\n",
        "\n",
        "    # Process neurofeedback group\n",
        "    for participant_id in participant_ids:\n",
        "        participant_epochs,_,_ = subject_wise_analysis_1_to_6(base_dir, participant_id, ROI=True)\n",
        "\n",
        "        # Add Reappraise data\n",
        "        if participant_epochs['Reappraise']['HbO'].size > 0:\n",
        "            X.append(participant_epochs['Reappraise']['HbO'])\n",
        "            y.append(1)  # Label 1 for Reappraise\n",
        "\n",
        "        # Add view_negative data\n",
        "        if participant_epochs['view_negative']['HbO'].size > 0:\n",
        "            X.append(participant_epochs['view_negative']['HbO'])\n",
        "            y.append(0)  # Label 0 for view_negative\n",
        "\n",
        "    # Process sham group\n",
        "    for participant_id in participants_id_sham:\n",
        "        participant_epochs,_,_ = subject_wise_analysis_1_to_6(base_dir_sham, participant_id, ROI=True)\n",
        "\n",
        "        # Add Reappraise data\n",
        "        if participant_epochs['Reappraise']['HbO'].size > 0:\n",
        "            X.append(participant_epochs['Reappraise']['HbO'])\n",
        "            y.append(1)\n",
        "\n",
        "        # Add view_negative data\n",
        "        if participant_epochs['view_negative']['HbO'].size > 0:\n",
        "            X.append(participant_epochs['view_negative']['HbO'])\n",
        "            y.append(0)\n",
        "\n",
        "    # Convert lists to numpy arrays\n",
        "    X = np.array(X)\n",
        "    y = np.array(y)\n",
        "\n",
        "    # Verify the shapes\n",
        "    print(\"X shape:\", X.shape)  # Should be (120, 13)\n",
        "    print(\"y shape:\", y.shape)  # Should be (120,)\n",
        "    return X, y\n",
        "\n",
        "X,y = process_fnirs_data()"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Create training dataset"
      ],
      "metadata": {
        "id": "GAY-OpB9J2TD"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from sklearn.model_selection import train_test_split\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "\n",
        "scaler = StandardScaler()\n",
        "X_scaled= scaler.fit_transform(X)\n",
        "# Split data into training and test sets\n",
        "X_train, X_test, y_train, y_test = train_test_split(\n",
        "    X_scaled,\n",
        "    y,\n",
        "    test_size=0.165,\n",
        "    random_state=42,\n",
        "    shuffle=True,\n",
        "    stratify=y\n",
        ")\n",
        "print(\"X_train shape:\", X_train.shape)\n",
        "print(\"X_test shape:\", X_test.shape)\n",
        "print(\"y_train shape:\", y_train.shape)\n",
        "print(\"y_test shape:\", y_test.shape)"
      ],
      "metadata": {
        "id": "ERialGPjKTfo",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "c4a63a40-3e32-43a2-d0cf-e564094d5b98"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "X_train shape: (88, 13)\n",
            "X_test shape: (18, 13)\n",
            "y_train shape: (88,)\n",
            "y_test shape: (18,)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## RFE (Recursive Feature Elimination)"
      ],
      "metadata": {
        "id": "IWSKtZhvmo32"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from sklearn.feature_selection import RFE\n",
        "from sklearn.svm import SVC\n",
        "\n",
        "def apply_rfe(X_train, y_train, X_test, n_features_to_select):\n",
        "    # Define the base classifier for RFE\n",
        "    base_clf = SVC(kernel='linear', C=1.0)\n",
        "    rfe = RFE(estimator=base_clf, n_features_to_select=n_features_to_select)\n",
        "\n",
        "    # Fit RFE on the training data and transform both train and test\n",
        "    X_train_rfe = rfe.fit_transform(X_train, y_train)\n",
        "    X_test_rfe = rfe.transform(X_test)  # Use transform (not fit_transform) for test data\n",
        "\n",
        "    print(\"Selected Features:\", rfe.support_)\n",
        "    print(\"Training data shape after RFE:\", X_train_rfe.shape)\n",
        "    print(\"Test data shape after RFE:\", X_test_rfe.shape)\n",
        "\n",
        "    return X_train_rfe, X_test_rfe, rfe\n",
        "\n",
        "# Apply RFE to both train and test data\n",
        "X_train_rfe, X_test_rfe, rfe = apply_rfe(X_train, y_train, X_test, n_features_to_select=4)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "DbZTjIEWmEId",
        "outputId": "e0223b63-04d6-473c-ab12-a6bae906d1a7"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Selected Features: [False False  True False  True False False  True False False False False\n",
            "  True]\n",
            "Training data shape after RFE: (88, 4)\n",
            "Test data shape after RFE: (18, 4)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Define hyperparameter function"
      ],
      "metadata": {
        "id": "XjVP0-E4e8dM"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from sklearn.decomposition import PCA\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "from sklearn.metrics import log_loss, accuracy_score\n",
        "\n",
        "def optimize_hyperparameters(X_train, y_train, n_trials=100):\n",
        "    def objective(trial):\n",
        "        # hyperparameter 'C': penalty for misclassification / margin violation\n",
        "        # smaller C allows for more margin violation, vice versa\n",
        "        svc_c = trial.suggest_float('C', 0.1954,0.195401)\n",
        "\n",
        "        # Initialize the classifier with a fixed linear kernel\n",
        "        clf = SVC(C=svc_c, kernel='linear', probability=True, random_state=42)\n",
        "\n",
        "        # Define cross-validation strategy\n",
        "        cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)\n",
        "\n",
        "        fold_losses = []\n",
        "\n",
        "        for train_idx, val_idx in cv.split(X_train, y_train):\n",
        "            # Split data for training and validation\n",
        "            X_train_fold, X_val_fold = X_train[train_idx], X_train[val_idx]\n",
        "            y_train_fold, y_val_fold = y_train[train_idx], y_train[val_idx]\n",
        "\n",
        "            # Train the model on training fold\n",
        "            clf.fit(X_train_fold, y_train_fold)\n",
        "\n",
        "            # Get probability predictions and class predictions\n",
        "            y_pred_proba = clf.predict_proba(X_val_fold)[:, 1]\n",
        "\n",
        "            # Calculate metrics for this fold\n",
        "            fold_loss = log_loss(y_val_fold, y_pred_proba)\n",
        "\n",
        "            # Append results\n",
        "            fold_losses.append(fold_loss)\n",
        "\n",
        "        # Return mean log_loss\n",
        "        return np.mean(fold_losses)\n",
        "\n",
        "    # Create an Optuna study to minimize the log_loss\n",
        "    study = optuna.create_study(direction='minimize')\n",
        "    study.optimize(objective, n_trials=n_trials)\n",
        "\n",
        "    # Output the results\n",
        "    print(\"Number of finished trials: \", len(study.trials))\n",
        "    print(\"Best trial:\")\n",
        "    trial = study.best_trial\n",
        "\n",
        "    print(f\"  Value (log_loss): {trial.value}\")\n",
        "    print(\"  Params: \")\n",
        "    for key, value in trial.params.items():\n",
        "        print(f\"    {key}: {value}\")\n",
        "\n",
        "    return study.best_trial.params"
      ],
      "metadata": {
        "id": "aQpWDYRDIyUs"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Model Evaluation"
      ],
      "metadata": {
        "id": "LtRswfK4mw_a"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from sklearn.metrics import accuracy_score, hinge_loss, confusion_matrix\n",
        "import optuna\n",
        "from sklearn.svm import SVC\n",
        "from sklearn.model_selection import cross_val_score, StratifiedKFold\n",
        "optuna.logging.set_verbosity(optuna.logging.WARNING)\n",
        "def main():\n",
        "    # Optimize hyperparameters using Optuna\n",
        "    best_params = optimize_hyperparameters(X_train_rfe, y_train, n_trials=200)\n",
        "    if best_params is None:\n",
        "        print(\"No optimal parameters found.\")\n",
        "        return\n",
        "\n",
        "        # Train the final model with the best hyperparameters\n",
        "    clf = SVC(\n",
        "        C=best_params['C'],\n",
        "        kernel='linear',\n",
        "        probability=True,\n",
        "        random_state=42\n",
        "    )\n",
        "    clf.fit(X_train_rfe, y_train)\n",
        "\n",
        "    # Evaluate on the test set\n",
        "    y_pred = clf.predict(X_test_rfe)\n",
        "\n",
        "    # Calculate accuracy\n",
        "    test_accuracy = accuracy_score(y_test, y_pred)\n",
        "    print(f\"Test Accuracy: {test_accuracy:.4f}\")\n",
        "\n",
        "    # Calculate loss (using hinge loss for SVM)\n",
        "    test_loss = hinge_loss(y_test, clf.decision_function(X_test_rfe))\n",
        "    # print(f\"Test Loss: {test_loss:.4f}\")\n",
        "\n",
        "    # Generate and print confusion matrix\n",
        "    cm = confusion_matrix(y_test, y_pred)\n",
        "    print(\"\\nConfusion Matrix:\")\n",
        "    print(cm)\n",
        "\n",
        "if __name__ == \"__main__\":\n",
        "    main()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "VQclKr-wb7s2",
        "outputId": "72c1be2f-e1e5-49fc-ef00-f4917556c1c0"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Number of finished trials:  200\n",
            "Best trial:\n",
            "  Value (log_loss): 0.6981508971494995\n",
            "  Params: \n",
            "    C: 0.19540099998697116\n",
            "Test Accuracy: 0.7778\n",
            "\n",
            "Confusion Matrix:\n",
            "[[5 4]\n",
            " [0 9]]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Test of Significance"
      ],
      "metadata": {
        "id": "Myvw9Cy5SZAa"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Chi-Square Test"
      ],
      "metadata": {
        "id": "sHLWvVDIjhL0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from scipy.stats import chi2_contingency\n",
        "\n",
        "# Your confusion matrix\n",
        "confusion_matrix = [[5, 4],\n",
        "                   [0, 9]]\n",
        "\n",
        "# Perform chi-square test\n",
        "chi2, p_value, dof, expected = chi2_contingency(confusion_matrix)\n",
        "print('chi-square:', chi2)\n",
        "print('p-value:', p_value)\n",
        "print('degree of freedom:', dof)\n",
        "print('expected value table:', expected)\n",
        "if p_value >= 0.05:\n",
        "    print('H0 is accepted')\n",
        "else:\n",
        "    print('H0 is rejected')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "1Xc3lsRZjg9K",
        "outputId": "7749e3c2-3188-49eb-b217-e1f6afda920a"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "chi-square: 4.430769230769231\n",
            "p-value: 0.035296596841407936\n",
            "degree of freedom: 1\n",
            "expected value table: [[2.5 6.5]\n",
            " [2.5 6.5]]\n",
            "H0 is rejected\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "import numpy as np\n",
        "\n",
        "# Confusion matrix\n",
        "confusion_matrix = np.array([[5, 4],\n",
        "                             [0, 9]])\n",
        "\n",
        "# Labels for the axes\n",
        "labels = ['Reappraise', 'View Negative']\n",
        "\n",
        "# Plot the confusion matrix\n",
        "plt.figure(figsize=(5, 4))  # Smaller overall box\n",
        "sns.heatmap(confusion_matrix, annot=True, fmt='d', cmap='Blues',\n",
        "            xticklabels=labels, yticklabels=labels, cbar=False, annot_kws={\"size\": 16})  # Larger annotations\n",
        "\n",
        "# Add labels, title, and adjust font sizes\n",
        "plt.xlabel('Predicted', fontsize=16)  # Larger font size for axes labels\n",
        "plt.ylabel('Actual', fontsize=16)\n",
        "plt.title('Confusion Matrix', fontsize=18)  # Larger title font\n",
        "plt.xticks(fontsize=14)  # Larger tick labels\n",
        "plt.yticks(fontsize=14)\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 407
        },
        "id": "nF-cLnec8Vxs",
        "outputId": "e290438c-2443-441a-f95b-7568165f8d39"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 500x400 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAGGCAYAAAC0W8IbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVw0lEQVR4nO3dd1gUV9sG8HvoSAelCFIEu8beK1hibCQSWxRRrHlfEzWK3SgxiZpimjEmRrFFjRpbrIkKilijYmxYQcSCShWkc74//JjXlQVhKTvg/bsur7Bnzpx5Zsnuw5k554wkhBAgIiIiRdLRdgBERERUMCZqIiIiBWOiJiIiUjAmaiIiIgVjoiYiIlIwJmoiIiIFY6ImIiJSMCZqIiIiBWOiJiIiUjAmaqJycvHiRQwcOBAODg7Q09ODJElo0qSJ1uIJCQmBJEmQJElrMZB6UVFR8u8mKipK2+GQljFRU4WSk5ODzZs3Y/jw4ahduzYsLS1hYGAAW1tbdOjQATNnzsSlS5e0HWY+kZGRaN++PbZs2YKHDx/CwsICdnZ2qFq1qrZDq5DykpgkSahXr94r6585c0ZlnxEjRpRqPOHh4Zg/fz6+/fbbUm2XCAD0tB0AUVGdPHkSfn5+uH79ulymr68PMzMzxMXFISwsDGFhYVi0aBH69++PjRs3wsDAQIsR/8/PP/+Mp0+fwsPDAyEhIXB0dNR2SKhSpQrq1Kmj7TBKLCIiAidOnEDbtm0LrLNq1aoyjSE8PByBgYFwcXHBpEmTStyevr6+/LvR19cvcXtUsbFHTRXCn3/+iS5duuD69euwsbHBwoULcf36dWRmZiIuLg6ZmZk4c+YMZsyYAXNzc2zbtg3Pnj3TdtiyixcvAgC8vb0VkaQBoFWrVoiIiEBERIS2Q9GYq6srACAoKKjAOunp6di0aRMkSYKLi0s5RVYyjo6O8u9GKf+/kPYwUZPi3bhxA8OGDUNGRgbq16+P8PBwzJgxA7Vq1ZLr6OrqokWLFli4cCEiIyPh7e2txYjzy/ujwdTUVMuRVC7Dhw+HJEn4/fffC/zDbNu2bUhMTETnzp3lxE5UkTBRk+LNmTMHycnJMDIywvbt2+Hk5FRofWtra+zYsQMWFhb5tj18+BABAQFo0KABTExMYGJiggYNGmDatGmIjY1V297LA3tiY2MxceJEuLm5wcjICHZ2dhg8eLDanqmrqyskSUJISAgAIDAwUOVeaV75/PnzIUkSunTpUuB5vWrw16lTpzB06FA5LhMTE7i4uKBz585YsGABYmJiitWeNt6v4nJzc0Pnzp2RnJyMP/74Q22dvMveI0eOLLStZ8+eYePGjRg+fDiaNGmCatWqwdDQENWrV8fbb7+Nffv2qd1PkiS57Tt37qj8fiVJwvz58+W6I0aMkO+RCyHw66+/okOHDrCxsYEkSVi9ejWAggeTxcXFwcnJCZIk4e2331YbT3Z2Ntq3bw9JkvDGG28gPT290POmCkAQKdjDhw+Fjo6OACBGjRpVorZCQkKEpaWlACAACBMTE2FiYiK/trKyEqGhofn2i4yMlOvs3r1b2NraCgCiSpUqwtDQUN5mbm4uwsPDVfZt0aKFsLOzE/r6+vIx7ezs5H9hYWFCCCHmzZsnAIjOnTsXGH9wcLB8rJetXr1aSJIkbzc0NBTm5ubyawAiKCioyO1p6/0qqhfPac2aNQKA8PT0zFcvKipKSJIkzMzMRGpqqujcubMAIPz8/PLVDQoKktuVJElYWFiIKlWqqLyHU6ZMybefnZ2d/F7r6Oio/H7t7OzEl19+Kdf18/MTAMTw4cOFj4+PvI+VlZXQ0dGRf0cvvoeRkZEqxwsJCZE/E0uXLs0Xz+zZswUAYWxsLC5fvly8N5YUiYmaFG3jxo0qX/qaio6OlpNO/fr1xbFjx+RtR48eFXXq1BEAhLW1tYiJiVHZ98UvTSsrK9G+fXtx5swZIYQQWVlZ4u+//xYODg4CgOjYsaPa4+cliHnz5qndXpJEnZqaKszMzAQAMWzYMHHz5k15W0pKivjnn39EQECA2LNnT5HaU8L79SovJuq885ckSdy+fVul3vz58wUAMXr0aCGEKDRR79ixQ0ydOlUcO3ZMpKamyuX3798XgYGB8h9bO3fuzLdvXpJ3cXEpNO68RG1qair09PTEV199JZKSkoQQQjx9+lTcv39fCFF4ohZCiLlz5woAwsjISPz7779yeXBwsJzEly9fXmgsVHEwUZOizZkzR/7CunfvnsbtjB8/Xk4cDx48yLf97t27cq/ov//9r8q2F78069atK549e5Zv/127dsl17t69m297WSbqU6dOyT3erKysAvcvantCaP/9epWXrxKMHj1aABAff/yxXCc3N1e4uroKAPKVi8IS9at8+eWXAoDo2rVrvm3FTdQAxPfff19gvVcl6uzsbNG+fXv5D6lnz56JJ0+eCEdHRwFA9O/fv7inRwrGe9SkaHFxcfLP1tbWGrUhhMDmzZsBAOPHj4e9vX2+Ok5OThg/fjwAYNOmTQW2NWXKFBgbG+crf+utt+SpYHkjvMuLpaUlAMgj4EuqIr5f/v7+AIA1a9ZACAEACA4ORlRUFOrUqYN27dqV+Bi9e/cGAJw4cQI5OTklasvKygrjxo3TeH9dXV1s2LABVlZWuHLlCiZOnAh/f3/cu3cPNWrUwK+//lqi+EhZmKip0ouMjER8fDwAoFu3bgXW6969O4DnfxxERkaqrdO6dWu15Xp6eqhWrRoAyMcqL+7u7qhbty6ysrLQunVrLF68GOHh4Ronk4r4frVt2xZ169bFnTt3cOjQIQBFH0T2otjYWMybNw9t27aFjY2NvIKcJEmoX78+gOeDzhISEkoUb8uWLUs8x9/Z2RkrVqwAAKxYsQK7du2Crq4u1q9fDysrqxK1TcrCRE2KZmNjI/+s6Rf6o0eP5J8Lm5P64mjyF/d5kZmZWYH76+k9Xz8oKyuruCGWiK6uLjZt2gQ3NzfcuXMHM2bMQNOmTWFubo7u3bvjp59+Ktac8or6fuUl5KCgICQnJ2Pbtm3Q1dXF8OHDi7T/iRMnULduXXzyySc4efIk4uPjYWxsDFtb23yryKWmppYoVltb2xLtn8fHxwc+Pj7y66lTp6JTp06l0jYpBxM1KVqDBg3kn8+fP6/FSJStcePGiIiIwB9//IGxY8eiYcOGSEtLw8GDB/Gf//wHdevWLfdL8uXN19cXurq62L59O5YvX460tDT07NkTDg4Or9w3OzsbQ4YMQWJiIpo0aYK9e/ciOTkZT58+RWxsLB4+fIiTJ0/K9fMur2tKV1e3RPvniYqKwsGDB+XXYWFhJb4sT8rDRE2K5unpCR2d5/+bbt++XaM2Xuy9vDyX+EUvbiutHk9R5fUuC5vzmpSUVGgbBgYG6N+/P37++WdcvHgRjx8/xvLly2FtbY27d+/Cz8+vSLFUhPdLHQcHB/Ts2RNpaWmYO3cugKJf9j5x4gTu3LkDXV1d7N69G2+99Va+qwEPHz4s9ZhLIu+Pi6SkJNSuXRuGhoY4duwYFixYoO3QqJQxUZOi2dnZyZf2NmzYoLLO96vk9Xrc3NzkgWh59y/VyeuZ2NjYwM3NTdOQNZJ3T/Hu3bsF1jl16lSx2rSxscG4ceOwePFiAM+vSBRlsFlFeL8KkjeoLDMzE1WrVkW/fv2KtF/e+16tWrUCL/e/2HN9Wd4fkyXtaRfHvHnzcPLkSVSpUgU7duyQf8+ffvopjh07Vm5xUNljoibF+/TTT2Fqaoq0tDT0798f9+7dK7R+QkICfHx85B6oJEkYNGgQgOcPx1DXM7p//z5+/vlnAMCQIUNK+QxerXHjxnIc6hLyo0eP5IFDL8vIyCi07RdHXecllMJUhPerIH379kVAQACmTJmCb7/9tsgPtMhbxS42NlbtimsxMTH4/vvvC9zf3NwcAJCYmFj8oDUQHByMRYsWAQC++eYb1KtXDxMnTkTv3r2Rk5ODoUOHlnjAGykHEzUpXu3atbFu3ToYGBjg8uXLaNKkCRYvXoybN2/KdXJycnD+/Hl8/PHHqFmzJrZt26bSxqxZs2BpaYn4+Hh069YNx48fl7eFhYWhW7duSExMhLW1NWbMmFFu55anXbt28gMj/Pz88M8//0AIgdzcXISEhKBLly7Izc1Vu++mTZvQvn17/Pzzz7h9+7ZcnpOTgwMHDsjn07Zt2yKPBlb6+1UQfX19fPHFF/jqq68wdOjQIu/XoUMHmJiYQAiBgQMHyldu8t7DLl26FLrUasOGDQEAycnJ8tS2shIXFwdfX1/k5uaif//+GDt2rLwtKCgIDg4OiI6OxpgxY8o0DipH2pvCTVQ8x44dEx4eHipLOhoYGAhra2t5NSb8//KPQ4YMEZmZmSr7h4SECAsLiwKXxLS0tBRHjx7Nd9xXLT6Rx8XFRe1SnUK8esETIYTYv3+/vPoV/n/JTSMjIwFA1KpVS2WVthe9uPQl/n/5UBsbG5X3pHr16uLq1asq+xVlCVFtvV+vktd+cfctbMGTn376SeV9NDU1ld//qlWrqizSou68unbtKm83MzMTLi4uwsXFRXzzzTdynbwFT1614Eph72G/fv0EAFGjRg0RHx+fb9+///5bXk72l19+KcK7QkrHHjVVGO3bt0dERAQ2btyIoUOHwsPDA0ZGRnj69Cmsra3RoUMHzJ49G1evXsWGDRvyXfbs3Lkzrl69iilTpqBevXrIzc2FEAL16tXD1KlTcfXqVXTs2FFLZwe8+eabCA0NRZ8+fWBlZYWcnBzUqFEDM2bMwNmzZ9UuPAIA/fr1w9q1azFy5Eg0btwYFhYWSEpKgpmZGVq1aoUFCxbg8uXLqFu3brHiUfr7VdrGjx+PPXv2oEuXLjA1NUV2djYcHR3xwQcf4MKFC2jUqFGh+2/duhWTJ09G7dq1kZWVhTt37uDOnTulejn8xx9/xK5du6Cjo1PgfOlu3bohICAAADBp0iRcvXq11I5P2iEJUY6jH4iIiKhY2KMmIiJSMCZqIiIiBWOiJiIiUjAmaiIiIgVjoiYiIlIwJmoiIiIFY6ImIiJSMD1tB6BtQ9aGazsEogqrfnVTbYdAVGHN7eZRpHrsURMRESkYEzUREZGCMVETEREpGBM1ERGRgjFRExERKRgTNRERkYIxURMRESkYEzUREZGCMVETEREpGBM1ERGRgjFRExERKRgTNRERkYIxURMRESkYEzUREZGCMVETEREpGBM1ERGRgjFRExERKRgTNRERkYIxURMRESkYEzUREZGCMVETEREpWIVJ1Nu3b8fAgQPxxhtvwMPDQy6PiIjAF198gXv37mkxOiIiorKhp+0AXiU3NxdDhgzB1q1bAQDGxsZIS0uTt1tZWWH27NnIycnBzJkztRUmERFRmVB8j/qbb77Bli1bMG7cOCQkJGDq1Kkq2+3s7NCxY0fs2bNHSxESERGVHcUn6tWrV6Nly5ZYtmwZzM3NIUlSvjoeHh6IjIzUQnRERERlS/GJ+ubNm+jYsWOhdWxsbBAXF1dOEREREZUfxSdqY2NjJCUlFVrnzp07sLS0LJ+AiIiIypHiE3XTpk1x4MABpKenq90eHx+P/fv3o02bNuUcGRERUdlTfKL+8MMPERMTAx8fH8TExKhsu3XrFt555x0kJSXhww8/1FKEREREZUfx07O8vb0xffp0LF68GC4uLjAxMQEA2NraIi4uDkIIzJ07F15eXlqOlIiIqPQpvkcNAAsXLsSBAwfQp08fVKlSBbq6usjNzUXPnj2xb98+BAYGajtEIiKiMqH4HnWe7t27o3v37toOg4iIqFxViB41ERHR60rxifrixYtYtWoVkpOT5bK0tDS8//77cHR0hLu7O5YvX67FCImIiMqO4hP1p59+irlz58LMzEwumzVrFn7++Wc8ffoUMTEx+O9//4u///5bi1ESERGVDcUn6tOnT8PT01NeOjQ7OxtBQUFo1aoVHj16hMjISFSrVg3fffedliMlIiIqfYpP1I8fP0aNGjXk12fOnEFycjLGjx8PIyMjVK9eHd7e3rhw4YIWoyQiIiobik/Uenp6yMjIkF+HhIRAkiR4enrKZTY2Nnjy5Ik2wiMiIipTik/Urq6uCA4Oll9v2bIFbm5ucHFxkcvu3bsHGxsbbYRHRERUphSfqH19fXHhwgW0bt0anTp1woULF/Dee++p1Pn3339Rq1YtLUVIRERUdhSfqCdMmIABAwbgn3/+wbFjx/DWW29h1qxZ8vbLly/jwoULXEKUiIgqJcWvTGZoaIjff/8dycnJkCRJZZoWANjZ2eH8+fNwdXXVToBERERlSPGJOo+5ubna8qpVq6Jq1arlHA0REVH5UPylbyIioteZ4nrUNWvWhCRJOHjwINzc3FCzZs0i7SdJEm7dulXG0REREZUvxSXq3NxceRUyda8LIoQoy7CIiIi0QnGJOioqqtDXRERErxPeoyYiIlIwJmoiIiIFU9yl74LExMQgODgY9+/fV1n7O48kSZg7d64WIiMiIio7FSJRBwQE4LvvvkNOTo5cJoSQB5nl/cxETURElY3iL32vWLECX3/9NTw9PbF161YIIeDn54eNGzdi/Pjx0NPTw4ABA3D48GFth0pERFTqFN+j/uWXX+Dq6op9+/ZBR+f53xWurq4YNGgQBg0ahIEDB6J79+4YMGCAliMlIiIqfYrvUUdERKBnz55ykgaA7Oxs+efOnTujd+/e+Oqrr7QRHhERUZlSfKIGAEtLS/lnExMTxMXFqWyvU6cOLl++XM5RERERlT3FJ2pHR0fExMTIr93d3XHq1CmVOpcuXYKJiUl5h0ZERFTmFJ+o27dvj5MnT8qvvb29cf78eYwbNw579uzBzJkzsW/fPnTq1EmLURIREZUNxQ8m8/X1xf3793Hnzh24uLggICAAu3fvxooVK/Drr79CCAFXV1d8+eWX2g6ViIio1EmiAj7NIisrCzt37sStW7fg4uKCvn37anzpe8ja8NINjug1Ur+6qbZDIKqw5nbzKFI9xfeojx49CnNzczRp0kQu09fXx7vvvqu9oIiIiMqJ4u9Re3p64pdfftF2GERERFqh+B61ra0tjIyMtB0GlZPx7ZzR2cO60DrD119AVm6Fu2NDVO7ObV+FKwf/AAA07uOLRm8N1nJEpAnFJ+ru3bsjJCREZW1vqvwiHqUgNjlT7bbcijesgqjcPb59BVcPbQckCeBnpkJTfKJetGgR2rZti7Fjx2Lx4sWwti68t0WVQ/CNeBy9Fa/tMIgqpOzMdBxf9w2MLaxg41Ibdy+c0HZIVAKKT9TDhg2DpaUlVq1ahfXr18PNzQ12dnb5eteSJOHQoUNaipKISDnO71yDp4/uw/P9+bhzLlTb4VAJKT5Rh4SEyD9nZGQgIiICERER+erxsjgREfDw+r+4duRP1GzdFY4NWzJRVwKKT9S5ubnaDoG0oIG9KZwtjWCkr4OUjBzcfPIM4feSkc1BZEQFykpPw8n138LIzBLN3x2j7XColCg+UdPrqZN7/rEICc+y8PPxaFy4/1QLEREp37ntK5ESF4vOY+fAsIqZtsOhUlLhEvXt27eRlJQECwsL1KxZU9vhUCm7k5CG1adjcPlBCp6kZsJATwcuVkbwaWyPOrammOrphs8P3sbV2BRth0qkKPevnsONY/vg0rwTajRuq+1wqBQpfsETAEhKSsLEiRNhZWWFWrVqoUWLFqhVqxasrKwwadIkJCUlaTtEKiX7rj7GgYgniElKR3p2LpLTs3HxQQrm77+JM9FJ0NPVwfCWjtoOk0hRMtNScXL9dzA0tUDLgeO1HQ6VMsX3qB89eoSOHTvixo0bsLS0ROfOnWFnZ4fY2FiEh4fj+++/x759+xAaGgpbW9tC28rIyEBGRoZKWU5WJnT1DcryFKiUbL3wAC2dLeBqbQzrKvqIf5al7ZCIFOGfrb/gWeITdPSfASNTC22HQ6VM8Yl65syZuHHjBmbMmIHZs2erPHwjNTUVn376KRYvXoxZs2bh119/LbSthQsXIjAwUKWswdvj0Ogd/gVaEdxL+t8fWTZM1ESyuxdOQNLRxbXQPbgWukdlW/LDuwCAmyf+woNr4TA2t0JH/+naCJM0pPinZ9na2uKNN97AwYMHC6zj5eWFS5cu4dGjR4W2pa5HPXpLBHvUFYSlsR5+GtAQABCwKwIxielajoj49Cxl+H3qQGSlpRaprom1Ld5ZEFTGEVFRVJqnZ6WmpqJNmzaF1mnbti1Onz79yrYMDQ1haGioUsYkXXG0dbUCADzLzMGDJCZpojyDvtpc4Lbja5fg9qlDXOu7AlP8YLKGDRsiKiqq0DpRUVFo2LBh+QREZcbFyhjNncyh89LaNRKALh7WGNzUAQCwP+IxchR9HYiIqPQovkc9a9YsDBo0CCNGjEC3bt3ybf/rr7+wdetWbN26VQvRUWmqZmqAKZ5uSMnIRmR8GpLSsmBioAsnS2NUM31+5SMsMgF/XHio5UiJiMqP4hN1UlISevTogTfffBPdu3dHhw4d5FHfoaGhOHjwIPr06YOEhASsXbtWZd/hw4drKWrSxJ2ENOy98gg1baqgurkh6lQzASQgKS0bJ6MSceRWHMLvcbETInq9KH4wmY6ODiRJwqvCfHGt77xHYubk5Lyy/SFrw0saItFri4PJiDRXaQaTBQVxdCIREb2+FJ+o/fz8tB0CERGR1ih+1DcREdHrTPE96jznzp3DmjVrcP78efmhHM2aNcPw4cPRrFkzbYdHRERUJipEog4ICMA333yT79nUx44dw9KlS/HRRx/hiy++0FJ0REREZUfxl76XLl2Kr7/+GrVq1cK6desQFRWFtLQ0REVFYe3atfDw8MDXX3+NZcuWaTtUIiKiUqf46Vn169dHamoqLl26BDOz/A9CT0pKQqNGjWBqaoorV64Uu31OzyLSHKdnEWmuqNOzFN+jjoyMhI+Pj9okDQAWFhbw8fFBZGRkOUdGRERU9hSfqF/1jOk8dnZ2ZRwJERFR+VN8oh4yZAj++OMPpKSkqN2enJyMP/74A0OGDCnnyIiIiMqe4hN1YGAgmjRpglatWmHTpk2IiYlBVlYWYmJisHHjRrRp0wbNmjVDYGCgtkMlIiIqdYqfnlWlShUAz9fvHjp0aL7tQghcu3YNxsbGKuWSJCE7O7tcYiQiIiorik/UHTt2VHngBhER0etE8Yk6JCRE2yEQERFpjeLvURMREb3OmKiJiIgUTPGXvgEgJycHmzdvxsGDB3H//n1kZGTkqyNJEg4dOqSF6IiIiMqO4hN1amoqevTogZMnT0IIAUmS8OKqp3mvOeCMiIgqI8Vf+v70009x4sQJBAYG4smTJxBCYP78+Xjw4AF+//131KxZEwMGDFDbyyYiIqroFJ+ot23bhjZt2mDOnDmwtraWy+3s7DBgwAAEBwfj4MGD+PLLL7UYJRERUdlQfKKOjo5GmzZt5Nc6OjoqvWcnJyf07t0ba9as0UZ4REREZUrxidrExAQ6Ov8L08LCAg8ePFCpY29vj+jo6PIOjYiIqMwpPlG7uLioJOGGDRvi8OHDcq9aCIFDhw7BwcFBWyESERGVGcUn6q5duyI4OFhet9vPzw/R0dFo27YtAgIC0KFDB4SHh8PHx0fLkRIREZU+xU/PGjNmDGxsbPD48WM4ODjA398f58+fx7JlyxAeHg4A8PHxwfz587UaJxERUVmQxIuTkiuQx48f4/bt23BxcYG9vb3G7QxZG156QRG9ZupXN9V2CEQV1txuHkWqp/gedUGqVauGatWqaTsMIiKiMlVhEvXDhw+xbds2REREIDU1FStXrgTwvGcdGRmJRo0a5XsmNRERUUVXIRL1smXLMGXKFHmktyRJcqJ+9OgR2rZti+XLl2PMmDHaDJOIiKjUKX7U959//okJEyagUaNG2LVrF95//32V7Q0aNMAbb7yBHTt2aCdAIiKiMqT4HvWXX34JZ2dnBAcHw8TEBGfPns1Xp1GjRggNDdVCdERERGVL8T3q8PBw9O7dGyYmJgXWcXR0RGxsbDlGRUREVD4Un6hzc3Ohr69faJ1Hjx7B0NCwnCIiIiIqP4pP1HXq1Cn0snZ2djaOHj2KRo0alWNURERE5aNI96i9vLw0PoAkSTh06JDG+w8dOhRTp05FYGAg5s2bp7ItJycHU6dOxe3btzF9+nSNj0FERKRURVqZ7MWnVxX7AJKEnJwcjffPyspCjx49cPToUbi7u8PIyAiXL1+Gj48P/vnnH0RFRaFHjx7Yt28fJEkqdvtcmYxIc1yZjEhzpboyWXBwcImCKQl9fX0cOHAAgYGBWL58ORISEgAAW7duhbm5OaZPn47AwECNkjQREZHSVai1voUQuHbtGuLj42Fubo569epBV1cXkZGRCAwMxOrVq4vdJnvURJpjj5pIc5VyrW9JklC3bl35dXR0NBYsWIC1a9ciOztbo0RNRESkZIod9X3s2DF4enrC3Nwc1tbW8Pb2xrVr1wAAz549w0cffYTatWtj5cqVqFatGr7//nstR0xERFT6StSjfvDgAXbu3Ilr164hOTkZ6q6iv7gud1GdPXsW3bp1Q2Zmplz2559/4p9//kFoaCj69euHK1euoHr16pg+fTrGjh3LedRERFQpaZyof/jhBwQEBCArK0suy0vUeQO7hBAaJeovvvgCmZmZWLhwIUaNGgUAWLFiBWbPno2OHTsiNjYWc+bMwaxZs2BkZKTpKRARESmeRoPJDh06hO7du8Pc3Bz//e9/ceTIEZw4cQLLly/H9evXsW3bNkRFRWHixIlo3Lgx/Pz8itW+k5MT6tati4MHD6qUd+3aFSEhIfjyyy/x0UcfFTdstTiYjEhzHExGpLmiDibT6B71d999B0mScODAAXz22WeoVasWAGDMmDH48ssvceXKFfj5+WHVqlXo2LFjsdt/9OgRmjdvnq88r6y4iZ+IiKii0ihRnz59Gs2aNUPr1q3Vbjc0NMRPP/0EIyMjfPLJJ8VuPzs7W+1DOPLKbGxsit0mERFRRaRRok5ISIC7u7v8Ou+hGWlpaXKZoaEhOnbsWKLlQ4mIiF53Gg0ms7a2RmpqqvzaysoKwPN5zXXq1JHLc3JyEBcXp1Fg69evx8mTJ1XKbt68CQDo1atXvvqSJGHPnj0aHYuIiEipNErUzs7OuHv3rvy6YcOGEEJg9+7dcqJOSUlBaGgonJycNArs5s2bcmJ+2f79+/OVcQlRIiKqjDRK1J07d8Y333yD2NhY2NnZoXfv3jAxMcGsWbPw8OFDODs7Y82aNYiPj8fgwYOL3X5kZKQmYREREVU6GiXqAQMG4Pz58wgPD8ebb74Ja2trLFmyBOPHj8eSJUsAPJ9D7erqisDAwGK37+LioklYRERElU6pPpTj3Llz2LJlC+Lj41GvXj2MHDkSFhYWpdV8meA8aiLNcR41kea08lCOZs2aoVmzZqXZJBER0WtNsQ/lICIiIiZqIiIiRdPo0reurm6R60qShOzsbE0OQ0RE9NrTKFEXZ/xZKY5VIyIieu1odOk7NzdX7b+cnBzcvn0b33//PaysrDBv3jzk5uaWdsxERESvjVId9S1JElxdXTFhwgQ0bNgQ3bp1Q8OGDeHj41OahyEiInptlNlgsi5duqBp06byAihERERUfGU66rtmzZq4ePFiWR6CiIioUivTRH3jxg0OJiMiIiqBMknU2dnZ+OyzzxAeHo6mTZuWxSGIiIheCxoNJvPy8ipw29OnT3H79m0kJiZCR0cHs2bN0jg4IiKi151GiTokJOSVdWrVqoVFixahZ8+emhyCiIiIoGGiDg4OLnCbgYEBHB0d4ezsrHFQRERE9JxGibpz586lHQcRERGpodHzqNeuXQsPDw+0a9eu0HonT57E9evXMXz4cI0DLGvpXIacSGNWLSdoOwSiCivt/NIi1dNo1PeIESPw66+/vrLeypUrMXLkSE0OQURERCjjedScQ01ERFQyZZqoHz16hCpVqpRae/Hx8bh7926ptUdERKR0RR5MdvToUZXXDx8+zFeWJzs7G5cvX8Zff/2FRo0alSjApKQkfPzxx9i0aROePHmi8nzrU6dOITAwEAsWLEDz5s1LdBwiIiIlKvJgMh0dHUiSBOD5Je28nwsjhMCvv/4Kf39/jYKLj49Hu3btcP36dTRr1gzp6em4evUqcnJyAABpaWmwt7fHqFGjNH74BweTEWmOg8mINFfUwWRF7lF36tRJTs5HjhyBra0t6tatq7augYEBnJyc4OPjg169ehX1EPnMnz8f169fx6ZNmzBw4EAEBgbik08+kbcbGxujc+fOOHz4sMbHICIiUrIiJ+oXVyPT0dHBW2+9hVWrVpVFTLJdu3ahT58+GDhwYIF1XF1dcfz48TKNg4iISFs0XpnM3t6+tGPJ58GDBxg8eHChdQwNDZGamlrmsRAREWmDolcms7GxeeUo74iICDg4OJRLPEREROVNo+lZ+/fvh5eXV6H3hg8dOgQvLy/8/fffGgfXqVMn7Ny5EzExMWq3X7lyBfv370e3bt00PgYREZGSaZSog4KCcPr0abRs2bLAOq1atcKpU6ewevVqTWPD7NmzkZOTg/bt2+O3337DkydPAABXr17FypUr4eXlBUNDQwQEBGh8DCIiIiXTaK1vd3d3ODg44NixY4XW69ChA2JjY3Hjxg2NA9y1axd8fX2RkpIC4H9Tw4QQMDMzw8aNG0s0spzTs4g0x+lZRJor9elZL3rw4AFat279yno1atRAeHi4JoeQ9evXD5GRkVizZg1OnTqF+Ph4mJubo3Xr1hg5ciSqVq1aovaJiIiUTKNEbWBggKdPn76yXkpKCnR0Sr5KqbW1NSZPnlzidoiIiCoajbJorVq1EBYWhmfPnhVY59mzZwgLC0PNmjU1Dm7lypVITk7WeH8iIqKKTqNE3bdvXyQmJmLChAlqn5AlhMAHH3yApKQkeHt7axzcmDFjYG9vj8GDB2P37t3y0qFERESvC40GkyUmJqJhw4Z48OABmjZtCn9/f3k50YiICKxatQrnz5+Hvb09Ll68CGtra42CW7JkCdatW4cLFy5AkiTY2Nhg8ODBGDZsGFq1aqVRmy/jYDIizXEwGZHmijqYTKNEDQAXLlxA3759ERMTk+8BHUIIODk5YdeuXWjSpIkmzau4fPky1q5di40bN8rH8/DwgK+vL4YOHQo3NzeN22aiJtIcEzWR5so8UQPP70OvWLECBw4cwJ07dwAAzs7O6NmzJ0aPHg0TExNNm1ZLCIHg4GCsW7cO27dvR3JyMiRJQrt27RAaGqpRm0zURJpjoibSXLkk6sLExcVh7dq1WLVqFS5evFjq7aenp+Pbb7/FvHnzkJ2drfH9ayZqIs0xURNprkznURdECIH9+/dj5cqV2L17N7KyskqzeQBAUlISNm/ejPXr1yMsLAy5ubmwsLAo9eMQEREpQakk6sjISKxatQqrV6/G/fv35ZHgzZo1w/Dhw0vcflZWFvbs2YN169Zh7969yMjIgJ6eHnr27AlfX98SjSwnIiJSMo0TdUZGBrZu3YqVK1fi6NGjEELIy3tOmzYNw4cPR/369UsU3LFjx7B+/Xps3boVCQkJEEKgefPm8PX1xZAhQ1CtWrUStU9ERKR0xb5HffbsWaxcuRKbNm1CUlIShBDQ09NDr1698O+//+LOnTulNt85b1UzZ2dnDB06FL6+vvI0sNLCe9REmuM9aiLNleo96oSEBKxfvx4rV66UB4YJIVC3bl34+/tj+PDhsLW1RceOHeXR36XB398fvr6+5fb8ayIiIqUpUqJ2cHBAVlYWhBAwNTXFoEGD4O/vj7Zt25ZpcL/++muZtk9ERKR0RUrUmZmZkCQJTk5OWLduHXu4RERE5aRIibpRo0a4ePEiYmJi4OXlhUaNGsHf3x9Dhw6FjY1NqQXj5eUFSZKwZs0aODk5wcvLq0j7SZKEQ4cOlVocRERESlHkwWT//PMPfv31V2zatEleEUxfXx/e3t7w9/dHjx490KlTJxw/flzjwWQ6OjqQJAlXr15F7dq1i/yITEmSuOAJkRZwMBmR5spsZbK0tDRs3rwZK1euxLFjx543IklwdHREWloa4uPjK9RTrpioiTTHRE2kuaIm6mI/5tLY2Bh+fn44evQorl27hmnTpsHOzg4xMTGIi4sDALRv3x6//PILkpKSits8ERERvUCj51HnqVWrFhYtWoS7d+9ix44d6NOnD3R0dHDixAm8//77cHBwwODBgzVu39/fH7t27Sq0zu7du+Hv76/xMYiIiJSsRIk6j66uLvr164ddu3bh7t27+Oyzz+Du7o709HRs2bJF43ZXr16N8PDwQutcuHABa9as0fgYRERESlYqifpF9vb2mDlzJq5fv47g4GAMGzastA+hIj09HXp6pfpsESIiIsUo0wzXuXPnEs+5liRJbbkQAnfv3sW+fftQvXr1Eh2DiIhIqUq9R11SOjo60NXVha6uLgBg/vz58usX/+np6cHNzQ3nzp0r0X1wIiIiJVPcNeNOnTrJveijR4/C2dkZrq6u+erp6urC2toaXl5eGDNmTDlHSUREVD4Ul6hDQkLkn3V0dDBy5Eh8/PHH2guIiIhIixSXqF+Um5ur7RCIiIi0SnH3qImIiOh/FN2jBoCcnBxs3rwZBw8exP3795GRkZGvDh/KQURElZWiE3Vqaip69OiBkydPQggBSZLw4tLkea8LmsJFRERU0Sn60venn36KEydOIDAwEE+ePIEQAvPnz8eDBw/w+++/o2bNmhgwYIDaXjYREVFloOhEvW3bNrRp0wZz5syBtbW1XG5nZ4cBAwYgODgYBw8exJdffqnFKImIiMqOohN1dHQ02rRpI7/W0dFR6T07OTmhd+/eXOubiIgqLUUnahMTE+jo/C9ECwsLPHjwQKWOvb09oqOjyzs0IiKicqHoRO3i4qKShBs2bIjDhw/LvWohBA4dOgQHBwdthUhERFSmFJ2ou3btiuDgYGRnZwMA/Pz8EB0djbZt2yIgIAAdOnRAeHg4fHx8tBwpERFR2VD09KwxY8bAxsYGjx8/hoODA/z9/XH+/HksW7ZMfk61j48P5s+fr9U4iYiIyookXpyYXEE8fvwYt2/fhouLC+zt7UvUVnp2KQVF9BqyajlB2yEQVVhp55cWqZ6ie9QFqVatGqpVq6btMIiIiMqcou9RExERve4U3aOuWbPmK+vo6OjA3NwcderUwTvvvIOBAweWQ2RERETlQ9H3qF1dXZGdnY379+8DAPT09FC1alU8efJEHglevXp1JCcnIyUlBZIkoWvXrti9ezcMDAyKdAzeoybSHO9RE2muqPeoFX3pOzw8HA4ODvDy8sLx48eRkZEhP0Hr+PHj6Nq1K6pXr47o6Ghcv34dvXr1wqFDh/D1119rO3QiIqJSoege9bhx43DixAmEh4errFCWJycnB02bNkW7du2wfPlypKeno379+jAzM8OFCxeKdAz2qIk0xx41keYqRY96586d6NWrl9okDQC6urro1asXdu7cCQAwMjKCl5cXbt68WZ5hEhERlRlFJ+rk5GQkJycXWicpKQlJSUny66pVq5Z1WEREROVG0Ym6fv362LhxI27fvq12++3bt7Fp0ybUr19fLouOjuYcayIiqjQUPT1r1qxZePfdd9GkSROMHj0a7du3h62tLR49eoSwsDCsXLkSKSkpmDVrFgAgMzMTf/31F3r06KHlyImIiEqHogeTAcCqVaswadIkefpVHiEETE1NsWTJEowePRoAkJiYiCNHjqBBgwbw8PAoUvscTEakOQ4mI9JcUQeTKT5RA8/vQ+/cuRMXLlxAcnIyzM3N0bhxY3h7e8PCwqJEbTNRK9NfB/bh940bcO1aBLKysuBcwxm9+vTFsOEjoK+vr+3w6P8xUSuLk50lpozojh7t68PRzhJPUzNw/updLNsYgv3HLms7PHpJpUrUZYmJWnm+WPgZflu/Fnp6emjZqg2qVKmC06dP4mlyMpo2a47lK1bByMhI22ESmKiVpHl9Z+z88b+wsTTBg8dJ+OdSFKwtTNCykSsM9PXw2c978enyvdoOk15Q6R7KkZKSguvXryM1NRUdO3bUdjhURg4fOojf1q9FlSpVsGrNetSr3wAAkJAQjzH+fjh/7ix+/OE7TAmYruVIiZTD0EAPG78aDRtLE2w5cBZj561HekYWgOcJfMfS/2D2uF44fv42Dp+K0HK0VFyKHvUNAFFRUfD29oaVlRVatmwJT09PeVtYWBjq16+PkJAQ7QVIperXX5YDAPxHj5WTNABYWVlj1px5AIBNG9bj6dOnWomPSIm8vRqjhoM1EpKf4YNPN8lJGgDOXonGwhX7AACzxvbUVohUAopO1NHR0WjTpg327t0Lb29vtG3bFi9eqW/dujWePHmCjRs3ajFKKi2xsbG4fOkiAOCt3n3ybW/WvAXs7R2QmZmJY0ePlHd4RIrVvIELAOD81WgkpaTl23741DUAQNsmNWFnY1ausVHJKTpRz5s3DwkJCThy5Ai2bt2K7t27q2zX09NDx44dERYWpqUIqTRFXL0CALCwsISTUw21deo3bKhSl4gAE2NDAEB8Yqra7XEJKQCeP22wST31ny1SLkUn6gMHDuCdd95Bu3btCqzj4uKCe/fulWNUVFbu3YsBANg7OBRYx97eXqUuEQGP45/fCnJzUr8y44vlrtW5emNFo+hEHR8fD1dX10LrCCGQkZFRPgFRmXqW+rw3YGxsXGCdKlVMAAApKep7DkSvo5Az1wEATevVQOM6Tvm2j363g/yzmSlnTFQ0ik7UdnZ2uHHjRqF1Ll68CGdn53KKiIhIeY6cuY7Qszego6ODrd+OQ69ODWFuagRXRxssnPwOhvZphcys53NRRW6ulqOl4lL09Kzu3btj3bp1+Pfff/HGG2/k2x4aGorDhw9j0qRJRWovIyMjX+9b6BrC0NCwNMKlEqpi8ry3nJaWfzBMnmfPnvekTU1NyiUmoopiaMBKbPp6DNo1dccf341X2fbD+sNo19QdzRu4ID7pmZYiJE0pOlHPmTMHW7duRadOnRAQECA/vnLfvn04fvw4lixZgqpVqyIgIKBI7S1cuBCBgYEqZbPnzsOcj+eXduikgerVHQEAsQ8fFFjn4cOHz+s6OpZLTEQVxeOEFHT1/wZereuiS6vasLYwwaP4ZOwOuYhzV6Jx+6/PAACXb97XcqRUXIpfmezUqVMYPHgw7ty5A0mSIISQ/+vs7IytW7eiRYsWRWqLPWpli334ED26dgYA7DlwUO3I7ze7dsHDhw+w6Iuv1U7hovLFlckqBjenqrjy53w8SUhBzR6zkZWdo+2QCJVoZbLWrVvjxo0b+PPPP3Hq1CnEx8fD3NwcrVu3hre3NwwMDIrclqFh/qTMJUSVw87eHg0aNsLlSxexb89ujBn3vsr2c2f/wcOHD2BgYIAOnTprKUqiimfS8K4AgFXbwpikKyDFJ2rg+Xzpd955B++88462Q6EyNnrseEz+8L9Y9esv6NCxk7w6WWJiAj7/9Plti8HvDYOZGRdtIHpR3Zr2uBebiKep6XKZrq4OPvLrhtE+7XEz+hEW/3pAixGSpipEoqbXh1fXbnhvmC82rF+HYUMGoXWbNjA2roJTp07gaXIymjRthv9+MFHbYRIpzqj+7THKpz3OX72L+48SYWigh5aN3GBf1Rw3ox+h9/ileJaeqe0wSQOKu0e9du1ajfYbPny4Rvvx0rcyHdi/9/ljLiOuIjs7G041nNG7T1/4Dh8B/WLc7qCyxXvUytG1TV2MG9gRTerWQDVrM2RkZuPGnVhs+/s8fvr9qMr636QMFfYxlzo6OpAk6ZX18gaVAYAkScjO1izjMlETaY6JmkhzFXYw2eeff16kRJ2dnY3Vq1fj1q1b5RAVERGRdiguUc+YMeOVdbZu3Yo5c+bg1q1bsLCwwPTpfDYxERFVTopL1IUJDQ3FtGnTcPr0aRgYGGDy5MmYPXs2rK2ttR0aERFRmagQifrSpUuYMWMG9u3bB0mSMGzYMCxYsIBrfBMRUaWn6EQdExODuXPnYv369cjJyUHPnj2xaNEitet+ExERVUaKTNSJiYn47LPP8OOPPyI9PR0tW7bE4sWL0aVLF22HRkREVK4Ul6gXL16MxYsXIykpCe7u7vj888/x7rvvajssIiIirVDsPGp3d3eMGDECenpF+1ti2rRpGh2P86iJNMd51ESaq9ALnuTJe0rWq0iShJwczRaaZ6Im0hwTNZHmKuyCJ0FBQdoOgYiISDEUl6j9/Py0HQIREZFi6Ly6ChEREWkLEzUREZGCMVETEREpGBM1ERGRgjFRExERKRgTNRERkYIxURMRESmYohP1ggULEBYWhuxsLh9GRESvJ8UtIfqivHW/jY2N0bZtW3h6esLT0xOtWrWCrq5uqRyDS4gSaY5LiBJprsIuIfqi4OBgBAcH4/DhwwgNDcWhQ4cgSRJMTEzQrl07OXG3aNFCZY1wIiKiykLRPeoXpaen49ixYwgODsaRI0dw5swZ+ZK4mZkZEhMTNWuXPWoijbFHTaS5Cvv0rKK4evUqdu3ahSVLluDx48d8ehaRljBRE2muUlz6znPz5k0cPnwYwcHBCAkJwaNHjyCEgIeHB95++214enpqO0QiIqIyoehEPXz4cAQHB+P+/fsQQsDFxQVvvfUWPD094eXlBUdHR22HSEREVKYUnajXr18PSZLQtWtXfPLJJ2jTpo22QyIiIipXih4q3a9fP1haWuLgwYNo3749GjZsiIkTJ2LHjh0aDx4jIiKqSBQ/mEwIgfDwcBw+fBghISEIDQ1FcnIydHR00LhxY/kyeK9evTRqn4PJiDTHwWREmqu0o75zc3Nx7tw57N+/Hz/88AOePHkCSZI0Xr2MiZpIc0zURJqrVKO+gecJ+syZM/IiKGFhYXj27BkAQE+vwpwGERFRsSg6w507d05OzKGhoUhJSYEQAnp6emjevLm8MlmHDh20HSoREVGZUHSibtGiBYDna343bdpUTswdO3aEqamplqMjIiIqe4pO1JMmTYKnpyc6deoECwsLbYdDRERU7hSdqJcsWaLtEIiIiLRK0Yk6T2ZmJg4ePIiIiAikpqZi7ty5AJ4/qCM5ORlVq1bl07OIiKhSUnx227VrF5ydndG3b19MnToV8+fPl7f9+++/cHBwwKZNm7QXIBERURlSdKIOCwvDu+++C0NDQ3z33Xd47733VLa3atUKHh4e+OOPP7QUIRERUdlS9KXvBQsWwNLSEmfPnkXVqlURFxeXr06LFi1w6tQpLURHRERU9hTdoz516hS8vb1RtWrVAuvUqFEDDx8+LMeoiIiIyo+iE3VGRgbMzc0LrZOYmMiBZEREVGkpOsPVrFkTZ86cKbTOiRMnULdu3XKKiIiIqHwpOlH7+PggLCwMQUFBard/9dVXuHTpEgYNGlTOkREREZUPRT89KyUlBW3atMHVq1fh5eWFjIwMhIWFYcqUKThx4gSOHz+OJk2a4Pjx4zA0NNToGHx6FpHm+PQsIs1VmsdcJiQkYMKECdi8eTNycnLkckmSMHDgQCxbtgxWVlYat89ETaQ5JmoizVWaRJ0nLi4OZ86cQXx8PMzNzdGyZUvY2dmVuF0maiLNMVETaa7SPY/axsYGPXv21HYYRERE5UrRg8mIiIhed4rqUfv7+0OSJHz++eews7ODv79/kfaTJAkrV64s4+iIiIjKn6LuUevo6ECSJFy9ehW1a9cu8kImkiSpDDQrDt6jJtIc71ETaa5C3qOOjIwEADg6Oqq8JiIiel0pKlG7uLggIyMDenp68msiIqLXmeIGkzk4OGDChAk4d+6ctkMhIiLSOsUl6vT0dCxbtgwtW7ZEs2bN8OOPPyIxMVHbYREREWmF4hJ1bGwsfvrpJ7Ro0QLh4eH48MMPUb16dbz33ns4dOiQtsMjIiIqV4oa9f2yK1euYNWqVfjtt98QGxsLSZLg7OwMf39/jBgxAjVq1CjxMTjqm0hzHPVNpLlKtYRoTk4O/vzzT6xatQr79+9HdnY2dHR00L17d4waNQre3t7Q19fXqG0maiLNMVETaa5SJeoXxcbGYu3atQgKCkJERAQkSYKNjQ0ePXqkUXtM1ESaY6Im0lxRE7Xi7lG/ip2dHQICAvD777+jffv2EEIgLi5O22ERERGVCUXNo36Vp0+fYsOGDVi5ciXOnj0LIQRMTEwwcOBAbYdGRERUJipEog4ODsaqVauwfft2pKWlQQiBNm3aYNSoURg0aBBMTU21HSIREVGZUGyijomJQVBQEFavXo2oqCgIIVCtWjWMHz8eo0aNQr169bQdIhERUZlTXKL+/fffsWrVKhw+fBg5OTnQ0dHBm2++KY/uzltelIiI6HWguKw3ZMgQAICbmxtGjhyJESNGwMnJSctRERERaYciE/WoUaPg5eWl7VCIiIi0TnGJ+rffftN2CERERIpR4RY8oddHRkYGFi5ciJkzZ8LQ0FDb4RBVKPz8VB5M1KRYycnJsLCwQFJSEszNzbUdDlGFws9P5VHhViYjIiJ6nTBRExERKRgTNRERkYIxUZNiGRoaYt68eRwIQ6QBfn4qDw4mIyIiUjD2qImIiBSMiZqIiEjBmKiJXuLq6gpXV1dth0HlICoqCpIkYcSIEdoO5bUnSRK6dOmi7TAUiYm6ksn74nnxn76+PhwdHTFw4ED8888/2g6RqNS99957kCQJGzduLLRecnIyqlSpAktLS6SlpZVTdJpbvXq1/DleuHCh2jqLFi2CJElYvXp1+QZXTF26dIEkSdoOo0JS3FrfVDrc3d0xbNgwAEBqairOnj2LLVu2YMeOHTh48CA6deqk5QiV69ChQ9oOgYpp1KhR2LhxI1atWiU/gU+djRs3Ii0tDX5+fjA2NoajoyOuXr0KCwuLcoxWM4sXL8a4ceNgbW2t7VDKxNWrV1GlShVth6FITNSVlIeHB+bPn69StmjRIsycORNz587FkSNHtBNYBeDu7q7tEKiYvLy84ObmhsOHDyM6OhrOzs5q661atQrA88QOAPr6+qhbt265xakpd3d33Lp1C5999hm+/vprbYdTJirC70FbeOn7NZL35XT27FmV8szMTCxZsgTNmjWDiYkJzMzM0LFjR+zatStfG9evX8e0adPQrFkz2NjYwMjICLVr18aMGTOQkpKSr37e5a709HTMmDEDzs7OMDIyQr169fDDDz/g5dmBeZf6Vq9ejZ07d6JVq1aoUqUKqlWrBn9/f8TGxuY7Rt69rXv37mH48OGwt7eHjo4OQkJCAADBwcHw9/dHnTp1YGpqClNTU7Ro0QK//PKL2vdJ3T3q9PR0fP3112jcuDEsLCxgYmICV1dXDBw4EBcuXMjXxs6dO9G1a1dYWVnByMgIDRs2xFdffYWcnBy1x6SSkSQJI0eORG5uLoKCgtTWuXz5Mk6fPo033ngDLVq0AFD4PeqnT59i3rx5aNCgAYyNjWFpaYk333wTx44dU6k3efJkSJKU77bS22+/DUmS5CtbeUJCQiBJEgIDA4t8fiNGjICHhwd+/PFHREdHF3m/yMhIjB49Gs7OzjA0NISDgwNGjBiBO3fuqK2/bds2tGjRAsbGxrCzs8OYMWOQkJCg9jNRnO8CSZLkzsGLt+VefN9fvkc9atQoSJKEo0ePqo11yZIlkCQJK1asUCn/999/MXjwYDg4OMDAwAAuLi744IMPEBcXV8R3TYEEVSqRkZECgHjzzTfzbXv06JEAICwsLOSy9PR00aVLFwFANGnSRHzwwQdi/PjxokaNGgKA+OGHH1TaWLhwobC2thY+Pj5i8uTJYuLEiaJ169YCgGjTpo3IzMxUqd+5c2cBQPTt21c4OTmJiRMniokTJwonJycBQHz00Ucq9YOCggQA0adPH6Gvry+GDBkiZs6cKTw9PQUA4eHhIeLj41X2ASAaNmwoatSoIRo3biwmTpwoxo0bJ86ePSuEEOLNN98U7u7uYujQoWL69Oli3LhxwsXFRe3xhRDCxcVFuLi4qJQNHDhQABBvvPGGmDhxopg2bZoYMmSIsLe3FytWrFCpO2PGDAFAODo6Cn9/fzF58mTRokULAUC8++676n9xVGLR0dFCR0dHuLq6itzc3Hzbp0yZIgCI7777Ti7L+7z4+fmp1I2LixMNGjQQAET79u3FpEmThL+/v7CxsRF6enpi+/btct2dO3cKAOKLL76Qy3JycoSVlZUAIKpXr67S9scffywAiCNHjrzynPI+DwsXLhS///67ACCGDx+uUmfhwoUCgAgKClIpP3nypLCwsBB6enri7bffFgEBAWLAgAFCT09P2Nrailu3bqnUX7lypQAgzM3NxdixY0VAQIBo0KCBaN68uahevXq+z0RxvgvmzZsnf+bmzZsn/3vxfQQgOnfuLL8ODg4WAMSYMWPUvjdNmjQRhoaGIiEhQS7buXOnMDQ0FMbGxmLw4MEiICBA9O7dWwAQtWrVyvfdUVEwUVcyhSXqzz//XAAQvXv3lstmzZolAIi5c+eqfLklJyeLFi1aCAMDA3Hv3j25PCYmRmRkZORrOzAwUAAQ69evVynPS9R16tQRiYmJcnliYqKoU6eOkCRJnDlzRi7P+2ICIPbv36/SVl4CnDBhgkp5Xv2RI0eK7OzsfLHdvn07X1lWVpbo3r270NXVFXfu3FHZ9nKiTkxMFJIkiebNm+drPzs7W+WL4q+//pLf/5SUFLk8NzdXjB8/XgAQW7duzRcPlY6ePXsKAOLgwYMq5VlZWcLOzk4YGhqKuLg4ubygRP3ee+8JAPn+CIuNjRU1atQQ1apVE2lpaUIIIRISEoSOjo5466235Hpnz54VAETXrl0FAHHt2jV5W8eOHYWxsbHaz9HLXkzUubm5omXLlkJHR0dcuHBBrqMuUWdmZgpXV1dhZmYmzp07p9JmaGio0NXVFX369JHLEhIShKmpqTAxMRHXr19Xed+8vLwEgHyJWtPvgoK8nKhzc3OFs7OzsLKyEunp6Sp1L168mO8P3ydPnghzc3Ph6OgooqKiVOpv3LhR7XdHRcFEXcnkffG4u7vLf7VOnTpV7pHa2dmJK1euCCH+91e/u7u72h7Irl271Paq1YmLixMAxIgRI1TK8z6cL39ohRBi3bp1+T48eV9M3bp1y1f/6dOnwtLSUpibm4ucnBy5HIAwMDAQjx8/fmWcL/rjjz8EALF69WqV8pcTdVJSktyzUvc+vahfv34CQL7kL8T/Er6Pj0+x4qSi27JliwAg3nvvPZXy7du3CwBi4MCBKuXqEvXjx4+Frq6u8PLyUnuM77//XgAQf/75p1zWvHlzYWpqKrKysoQQQnz11VcCgAgLCxMAxE8//SSEEOLZs2fCwMCgwLZf9mKiFkKIw4cPCwAqfxSoS9Tbtm0TAMQnn3yitt3+/fsLHR0dkZSUJIQQYvXq1QKA+PDDD/PVPX78uNpEXZBXfRcU5OVELYQQM2fOFADEH3/8oVI+bdo0AUDs2LFDLluyZIkAINauXau2/WbNmomqVasW6RyUhoPJKqlbt27luwdmb2+P0NBQeHh4AACuXbuGhIQEVK9eXe39ssePHwMAIiIi5DIhBIKCgrB69WpcunQJSUlJyM3Nlbffv39fbTwdO3YssOz8+fNFqm9qaoomTZogJCQEt2/fls8DANzc3FC1alW1x3769Cm++uor7NixA7du3UJqaqrK9oJizmNubo5evXph7969aNasGQYMGIAuXbqgZcuW0NfXV6l78uRJmJiYyIOWXmZsbKzyflLp8vb2RrVq1bB9+3YkJSXJo7lfHkRWmDNnziAnJwcZGRn5BmQCwI0bNwA8/1z06dMHAODp6YmzZ8/izJkzaNu2LYKDg1GvXj20a9cOLi4uCA4Oxvjx4xEWFobMzEx4enpqdH6enp7o2bMn9u3bhyNHjqBz585q6508eRLA88+4unN4+PAhcnNzcf36dbRo0UIeZ9GhQ4d8dVu3bg09vfypQtPvguLw9fXFwoULsW7dOvTv3x8AkJubiw0bNsDGxga9evWS6+ad86lTp3Dr1q18baWnp+PJkyd48uRJgd8VSsVEXUm9+eab2L9/P4DnCXfNmjWYPn06+vXrh9OnT8PU1BTx8fEAng+yuXz5coFtvZjYPvzwQyxduhQ1atRAv3794ODgIC/6HxgYiIyMDLVt2NnZFViWlJRUpPqF7VNQ/czMTHTp0gXnzp1D06ZN4evrCxsbG+jp6SEqKgpr1qwpMOYXbdmyBZ9//jk2bNiA2bNnA3iewEeOHInPP/9cnlYSHx+P7OzsQgcKvfyHApUefX19+Pr6YsmSJdiwYQPef/99PHz4EPv27YOzszO6dev2yjbyPhdhYWEICwsrsN6Lv0dPT0989dVXCA4ORqtWrRAaGgpfX1952969ewE8H9iYV6apRYsW4a+//sK0adNw6tSpQs/ht99+K7StvHNITk4GANja2uaro6OjozaxafpdUBz16tVD8+bNsXfvXiQkJMDKygohISGIiYnBf/7zH5U/lPPO+ccffyy0zdTU1AqXqHnpu5Ip7B71nDlzBAAxceJEIYQQ//77rwBQ5EuxsbGxQpIk8cYbb4jU1FSVbQ8ePFB76Srvcpe6S8FRUVHyJeU8eZf6AgMD1caQN/Dtxo0bcpm64+bZvHmzACBGjRqVb1vefat58+aplKsbTPai27dvi5UrV4qWLVsKAGLs2LHyNmtra2FjY1PgvlT2Ll++LACI5s2bCyGE+OKLL9T+noVQf+k775bPlClTinzM5ORkoaenJ7p16yZOnz6tMhZhzZo1AoC4dOmSaNeunTAxMck36LIgL1/6zuPr6ysAiM2bN6u99P3RRx/luzxfmMmTJ8vtvSwnJ0fo6empfCZK8l1QkII+x99++60AIJYvXy6EEGLEiBECgDhx4oRKvf79+wsA4uLFi68424qH07NeI7NmzUL16tWxbNkyREVFoV69ejA3N8c///yDrKysV+5/+/ZtCCHQrVu3fAsThIaGFrqvuu15ZU2bNi1S/ZSUFISHh8Pc3Bw1a9Z8ZbwA5Etg3t7eRTpGUbi5ucHf3x9HjhyBqampyjS21q1bIy4uTr48SuWvfv36aNOmDc6ePYt///0XQUFB8vStomjZsiUkScKJEyeKfEwzMzM0b94cYWFh2L9/PyRJknvNXl5eAIA///wTZ86cQbt27fLdMimuBQsWwNDQELNnz0Z2dna+7a1btwaAIp9D48aNAUDtFYTTp0/nO4Ym3wW6uroAUOwpikOGDIGenh7Wr1+PtLQ0bNu2DR4eHmjTpo1KveKec0XCRP0aMTY2xvTp05GVlYUFCxZAT08P77//Pu7cuYOpU6eqTdaXLl3Co0ePAAAuLi4AgOPHj6vci4qJicHMmTMLPfaCBQtULlcnJSXh008/hSRJ8PPzy1f/4MGDOHDggErZZ599hsTERAwfPhw6OkX7Xzcv5pfnvh45ciTf/MuCPH78GJcuXcpXnpCQgIyMDBgZGcllH374IQDA399f7bzNhw8f4urVq0U6Lmku7170f/7zH1y9ehXdunWT/194FXt7ewwcOBDHjx/Hl19+mW+uP/D8PuizZ89Uyjw9PZGWloYffvgBjRs3llcQc3JygoeHB5YsWYKsrKwSXfbO4+Ligv/85z+4ceOG2qVDvb294ezsjCVLlqidh5yVlaXymfD29oapqSlWrlypcn83Ozsbc+fOVXt8oHjfBXnvx927d4t2kv/P1tYWPXr0QFhYGL799lskJyfnm5sOACNHjoSZmRlmz56t9lbes2fP5PvYFY6We/RUygq79C2EEGlpaaJ69epCT09P3Lx5U6Snp4vu3bvLI8VHjhwppk+fLoYNGyYaN26c7xKTj4+PACCaNm0qpk6dKnx9fYWVlZV49913C73cVVrzqN3d3dXOoy7o0vfTp0+Fq6urACB69eolpk2bJry9vYWurq4c86sufZ8/f14AEI0bNxa+vr5ixowZYsyYMcLW1lYAEEuXLlXZf+7cuQKAsLS0FIMHDxbTp08Xo0ePFl26dBG6urr5LmNS6UtOThYmJiby1L1NmzaprVfYPOomTZoIAKJRo0Zi7Nix8tz5WrVqCQDiwYMHKvscOHBAPt7kyZNVto0ZM0be9vIl28IUdOlbiOfTkSwsLOR2X55Hffr0aWFjYyMACC8vLzFx4kQxadIk8c4774iqVauKOnXqqNT/5ZdfBPB8nYVx48aJadOmiYYNG8rzqN3c3FTqF/e7YNmyZQKAaNasmZg9e7ZYsGCB2LVrl7y9sM9x3m0qfX19AUDcvHlTbb3du3cLY2NjoaurK3r37i2mTJkiJkyYIPr06SPMzMwK/F5UOibqSuZViVoIIX744QcBQPj6+gohns8F/vnnn0X79u2Fubm5MDQ0FM7OzqJnz57ip59+UpkP/PTpUzFlyhTh6uoqDA0NRa1atcSCBQtEZmZmoYk6LS1NTJs2TdSoUUMYGBiIOnXqiO+//z7fdKe8L6agoCCxY8cO0bJlS2FsbCxsbGzEiBEj8n05ClH4B1yI5/eUfXx8RLVq1USVKlVEy5YtxaZNm+QFFV6VqBMSEsT8+fNFp06dhIODgzAwMBDVq1cXPXv2FPv27VN7zL///lv07dtXVKtWTejr6wt7e3vRtm1bsWDBAhEdHV1grFR6Ro4cKQAIa2vrfPNw8xSUqIV4PpXqiy++EM2bNxcmJibC2NhYuLm5ibffflusXbtWnoqVJzU1VU4kL98b3rBhgwCgMoWrKApL1EL8b20EdYlaiOdznSdOnChq1aolDA0Nhbm5uahXr54YPXq0OHToUL76W7ZsEU2bNhWGhobC1tZWjB49WsTFxQlTU1PRuHFjlbrF/S7IysoS06ZNE87OzkJPTy/f+17Y5/jZs2fC3NxcABBt27Yt7C0TERERYtSoUcLFxUUYGBgIKysr0ahRI/Hhhx+K06dPF7qvUklCqLmuQ1RKunTpgiNHjqi9fKjO6tWrMXLkSAQFBfHRg0QKcPPmTdSqVQsDBw7E77//ru1wXku8R01ERPKYixelpaVh8uTJAJ6vXU7awXnURESEI0eOYNSoUejRowecnZ3x5MkTHD58GFFRUfDy8sKgQYO0HeJri4maiIjQoEEDdO/eHWFhYdixYweA54/LXbBgAaZOnVrkmRZU+niPmoiISMH4JxIREZGCMVETEREpGBM1ERGRgjFRExERKRgTNRERkYIxURNpgaurKyRJUvlnaGgIZ2dnDBo0SOMne5W2+fPnQ5IkzJ8/X6V89erVkCSpQq8eV9C5ESkNEzWRFrVv3x5+fn7w8/PDW2+9hdzcXGzevBmdO3fGkiVLtB1eucj7oyUqKkrboRApEhc8IdKi0aNHq/RK09PTMW7cOKxduxbTpk1Dnz59ULt2be0FWIB33nkHbdq0gYWFhbZDIar02KMmUhAjIyP8+OOPMDExQU5ODrZt26btkNSysLBA3bp14eDgoO1QiCo9JmoihTE1NUWdOnUAQL4cnHcfGwCCgoLQtm1bWFhY5LtkfP/+fXz00UeoV68eqlSpAjMzM7Rs2RJLly5Fdna22uOlpaVh/vz5qFWrFgwNDeHg4AA/Pz9ER0cXGOOr7lHfu3cPAQEBaNSoEczMzGBiYoLatWtjxIgROH78uEobd+7cAQC4ubmp3LMPCQlRabO8zo1IaXjpm0iBkpOTAQCGhoYq5R988AGWLVuGdu3aoXfv3rh9+7acwI8ePYq3334bCQkJcHV1Rffu3ZGRkYHTp0/jgw8+wJ9//ondu3dDX19fbu/Zs2fo2rUrTp48CRMTE/To0QPGxsY4cOAA9uzZg969exc79kOHDuHdd99FYmIibG1t0bVrVxgYGCAqKgobNmwAALRr1w4eHh7w8/PD1q1bkZqaCh8fH5iamsrt2Nvbyz8r5dyItEKrT8Mmek25uLgIACIoKCjftgsXLggdHR0BQKxatUoIIQQAAUCYm5uLEydO5NvnwYMHwsbGRkiSJJYtWyZycnLkbU+ePBFeXl4CgAgMDFTZb+rUqQKAqFu3rrh3755cnpqaKry9veXjzps3T2W/oKAgAUD4+fmplEdHRwsLCwsBQMyYMUNkZGSobI+NjRWhoaFq34vIyEi171V5nxuR0jBRE2mBukSdmJgo9uzZI9zd3QUAUb16dZGSkiKE+F+i/uSTT9S2N336dAFATJgwQe32mJgYoa+vL6pVqyZyc3OFEEI8e/ZMmJmZCQBi3759+fZ58OCBMDIyKlainjRpkgAg+vbtW8R34tWJurzPjUhpeI+aSItGjhwp35O1tLRE7969cevWLbi7u2Pv3r0wMTFRqf/uu++qbWfPnj0AUOAzgx0dHVGrVi08fvwYN27cAACcO3cOT58+RdWqVdGzZ898+9jb26NHjx7FOp/9+/cDAMaOHVus/QqjlHMj0hbeoybSovbt28PDwwMAYGBgAFtbW7Rp0wY9e/aEnl7+j6erq6vadm7fvg0A6Nix4yuP+fjxY9SuXRsxMTGFtgk8H+BVHHkDw+rWrVus/QqjlHMj0hYmaiItenke9asYGxurLc/NzQXwvMf9ci/8ZTY2NkU+nhJU5nMjKgomaqJKoEaNGrhx4wamT5+OFi1aFGkfR0dHACh0RbDirhbm7OyMa9euISIiQr5SUFJKOTcibeE9aqJK4K233gIAbN68ucj7NG/eHKampnjy5An++uuvfNtjY2PVlhcm737wihUriryPgYEBABQ4F1op50akLUzURJVAQEAALC0tsWTJEnz99dfIzMzMVycyMhLr16+XXxsbG8uDviZPnowHDx7I29LS0vD+++8jLS2tWHF89NFHMDMzw65duzBnzhxkZWWpbH/06BGOHTumUubk5AQAuHz5sqLPjUhrtD3snOh1VNg8anXw/9OzCnPkyBFRtWpVAUDY2toKLy8vMXToUNGnTx95ylfr1q1V9klJSRGtWrUSAISpqano27evGDBggLC3txc2NjZi+PDhxZqeJYQQBw4ckKdG2dnZibffflsMGDBAtGrVSujr6+fbZ+nSpfLx+/fvL0aNGiVGjRolIiIitHJuRErDe9RElUSnTp1w+fJlLF26FHv27MGZM2eQkZEBW1tbODs7Y9iwYfDx8VHZx8TEBMHBwVi0aBE2bNiAAwcOwMrKCt26dcOnn36K1atXFzuOHj164NKlS1iyZAn279+P/fv3Q09PD9WrV4evry/GjBmjUv/999/H06dPsX79euzduxfp6ekAgGHDhslLqSrl3Ii0QRJCCG0HQUREROrxHjUREZGCMVETEREpGBM1ERGRgjFRExERKRgTNRERkYIxURMRESkYEzUREZGCMVETEREpGBM1ERGRgjFRExERKRgTNRERkYIxURMRESkYEzUREZGC/R+Tf4nhRag8pAAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        }
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [
        "q8463p6_t4mT"
      ],
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.11.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}