exp1 figure1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TempAtt: exp1 (auditory)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load & preprocess data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "clear\n",
    "clc\n",
    "% close all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "% define Paths: replace XXX by your own path\n",
    "origDir = '/XXX/exp1/'; %!\n",
    "addpath(genpath(origDir)); cd(origDir)\n",
    "\n",
    "% participants\n",
    "SUJETS = 1:31;\n",
    "SUJETS(SUJETS==9 ) = []; % to remove the unfinished participant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "% load data: 2 sessions per participant\n",
    "IN = []; Dif = NaN(SUJETS(end),1);\n",
    "for isuj = SUJETS\n",
    "    subject = ['S' sprintf( '%02d', isuj )];\n",
    "        ifile = dir(['TA18_' subject '_1*.mat']);\n",
    "        load([ifile.name])\n",
    "        IN{isuj,1}  = X;\n",
    "        Dif(isuj,1) = cfg.ndis;\n",
    "    subject = ['S' sprintf( '%02d', isuj ) 'b'];\n",
    "        ifile = dir(['TA18_' subject '_1*.mat']);\n",
    "        load([ifile.name]) \n",
    "        IN{isuj,2}  = X;\n",
    "        Dif(isuj,2) = cfg.ndis;\n",
    "    clear X subject ifile\n",
    "end\n",
    "Dif = mean(Dif,2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "% varnames & beat\n",
    "varname = {'Proportion Correct', 'Mean RT (ms)', 'Hit IE (sec)'};\n",
    "beats = cfg.fbeat(end:-1:1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "% extract data: group 2 sessions, remove trials with suspect reaction times, average trials per condition\n",
    "OUT{1} = NaN(SUJETS(end), length(cfg.fbeat)); OUT{2} = OUT{1}; OUT{3} = OUT{2};\n",
    "for isuj = SUJETS\n",
    "    for ibeat = 1:length(beats)\n",
    "        out = [];\n",
    "        for isess = 1:2 %!\n",
    "            y = IN{isuj,isess}.beat == beats(ibeat);\n",
    "            out{1}(isess,:) = IN{isuj,isess}.iscor(y);\n",
    "            out{2}(isess,:) = IN{isuj,isess}.rt(y);\n",
    "        end\n",
    "        % group 2 sessions\n",
    "            out{1} = reshape(out{1},[],1);\n",
    "            out{2} = reshape(out{2},[],1);        \n",
    "        % nan rt >|< 2sd\n",
    "            rtmin = mean(out{2}) -2*std(out{2});\n",
    "            rtmax = mean(out{2}) +2*std(out{2});\n",
    "            out{2}(out{2}>rtmax)  = NaN; out{2}(out{2}<rtmin) = NaN;\n",
    "            out{1}(isnan(out{2})) = NaN;\n",
    "        % avg trials\n",
    "            OUT{1}(isuj,ibeat) = nanmean(out{1});\n",
    "            OUT{2}(isuj,ibeat) = nanmean(out{2});\n",
    "            OUT{3}(isuj,ibeat) = nanmean(out{2}) ./nanmean(out{1});\n",
    "        clear out y rtmin rtmax\n",
    "    end\n",
    "end\n",
    "clear IN  isuj ibeat isess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "% remove the phantom (NaN) participant\n",
    "OUT{1}(isnan(OUT{1}(:,1)),:) = [];\n",
    "OUT{2}(isnan(OUT{2}(:,1)),:) = [];\n",
    "OUT{3}(isnan(OUT{3}(:,1)),:) = [];\n",
    "Dif(isnan(Dif)) = [];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "% define variable of interest\n",
    "i0 = 1; % performance (or RT, or IE, see varname)\n",
    "fit_order = 3; % fit order (n = 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "% Fits per subject - define after rstudio analysis\n",
    "xVec = (beats(end):0.01:beats(1));\n",
    "Fit = [];\n",
    "for j0 = 1:length(OUT{i0})\n",
    "    x = OUT{i0}(j0,:);\n",
    "    x = polyfit(log(beats), x, fit_order);\n",
    "    if fit_order == 3\n",
    "        Fit(j0) = exp((-x(2)-sqrt(x(2).^2-3.*x(1).*x(3)))./(3.*x(1))); % order 3\n",
    "    elseif fit_order == 2\n",
    "        Fit(j0) = exp(-(x(2)./(2.*x(1)))); % order 2\n",
    "    end\n",
    "end\n",
    "Fit = real(Fit)';\n",
    "nsuj = length(Fit);\n",
    "\n",
    "% save data for inter-experiment comparisons\n",
    "save([origDir 'exp1_values.mat'], 'Fit', 'Dif');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "% option: exclude outlier participants in the Fit\n",
    "% olr = find((Fit>(mean(Fit)+4*std(Fit)))); % detect outliers (4std)\n",
    "% [min(Fit), max(Fit)] % edit range of Fit fq\n",
    "% % exclude outlier on all data types\n",
    "%     nsuj = nsuj-numel(olr);\n",
    "%     Dif(olr) = [];\n",
    "%     Fit(olr) = [];\n",
    "%     OUT{i0}(olr,:) = [];\n",
    "% [min(Fit), max(Fit)] % new rang of Fit fq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Parametric repeated-measure Anova on Proportion Correct \n",
      "Beat                 F(  7.00,203.00) =  15.277,\tp = 0.0000\n",
      "\n",
      "3nd order fit: r2=0.86, p=0.002 \n"
     ]
    }
   ],
   "source": [
    "% statistics\n",
    "fprintf('\\nParametric repeated-measure Anova on %s \\n', varname{i0})\n",
    "[~,F,~,~] = repanova(OUT{i0}, 8, {'Beat'});\n",
    "\n",
    "% compute fit quality:\n",
    "x = mean(OUT{i0},1);\n",
    "    y = polyfit(log(beats), x, fit_order);\n",
    "    y = polyval(y, log(beats));\n",
    "[r,p] = corr(x', y', 'type', 'Spearman');\n",
    "fprintf('\\n%1.0fnd order fit: r2=%1.2f, p=%1.3f \\n', fit_order, r*r, p)\n",
    "length(x);\n",
    "clear x y F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAB3RJTUUH4wsVCzsCkvo0xQAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAyMS1Ob3YtMjAxOSAxMjo1OTowMkNvQqYAACAASURBVHic7J15XFNn9v+fe5PcAEkghE2ikKAIVYOoVepWsNJqXapFWy1WrdLazji1Vov71yr6q3WpaHFjrNaFUtw6Ki6IDirouJZaBAHZJCAJS4CQBLLe3N8f10kzgDEKIQvP+9XXvJKHGzgMJp/7POecz0EIggAQCAQCgVgb1NoBQCAQCAQCABQkCAQCgdgIUJAgEAgEYhNAQYJAIBCITQAFCQKBQCA2ARQkCAQCgdgEUJAgEAgEYhNAQYJAIBCITQAFCQKBQCA2ARQkCAQCgdgEUJAgEAgEYhNAQYJAIBCITQAFCQKBQCA2ARQkCAQCgdgEUJAgEAgEYhNAQYJAIBCITQAFCQKBQCA2ARQkCAQCgdgEUJAgEAgEYhNAQYJAIBCITQAFCQKBQCA2ARQkCAQCgdgEUJAgEAgEYhNAQYJAIBCITQAFCQKBQCA2ARQkyDOKi4uPHDly7do1HMfbvUAoFCYlJRUUFHRxYBAIpJtgWUHKz88/fPjwrVu3CIKw6A+CdJD9+/dPnjw5IyNjxYoVS5YsaXvBtm3bJk+efO3atdmzZy9atKjrI4RAII4PYTGWL18+aNCgmJiY0NDQb775xnI/CNJBxGJxSEjI2bNnCYIoKCgICgr6448/jC9oamoKDQ09deoUQRAPHjwICgoqLS21TqwQCMRxsdQO6fHjx2fPnj1w4MDBgwfj4+PPnTv39OlTC/0sSAfJz89HUXTixIkAgNdeey04ODgrK6vVNRqNhnyAYRgAoKGhoYuDhEAgDg/VQt+Xy+UePXr09ddfBwC0tLQAAJ6XmYBYnadPn3p5eVGpz/4x9OjRo66uzvgCV1fXpUuXrl279sKFCw8fPhw3btzgwYOtESkEYn3y8/Pv3bsXFBQ0YsQIBEHaXiAUCrOysoYOHdqvX7+uD8+usdQOicVihYWFSSSS2NjYdevWRUdH83g8C/0sSAdpbGx0cnIyPHV2dm5ubja+oLm5+eLFi+7u7lwu193dXSQS1dTUtP0+c+bM2bVrl8XDhUCsx4oVKz7++OMbN24sXLhw2bJlbS+A2daOYNYOSa1Wl5SUFBcXFxcXazSaAQMGDBgwoHfv3hQKxfQLm5ub3d3d3d3ds7OzKyoq/P39W10wZ86csLAw+GezLhwOR6FQGJ7K5fLg4GDjC65fv15QUJCRkcHlcrVabURExPHjx9vWPty7dy8pKakrIoZArAGZiUhOTn799devXr3697///euvv+7Vq5fhAplMlpycvH79+unTp//5558zZ84sKyvr3bt3q++za9cu+KHXLi8QJBzHT5w4sWvXrvr6egAAjUZDUVStVgMAnJ2dp06dunjxYg6H87yX83i8NWvWxMbGjh49+tSpU0uXLm11AfwIswV69uxZW1srk8lcXV0JgigtLY2MjDS+oKyszNvbm8vlAgBoNFpQUFBFRYWVgoVArIY5mYi22da2grR7924oSO1i6sju6tWrkydPPnDgQFRU1K+//nrz5s3c3Nw///wzIyPjp59++vvf/56VlTVu3Liff/5Zq9W2eu3NmzcXL15M/rXodDqXyxUKhRb8PSAdYPTo0Ww2e8OGDSKRaPv27XK5fOzYsQCA48eP37t3DwAwfPjw6urq5ORklUp1+fLle/fuzZgxw9pRQyBdzQszEYZsa0xMzNy5c2G29WV5riARBPHdd9/94x//uHLlyrJly15//XUvLy8EQVAU7dWrV3h4+BdffPHvf/973bp1p0+fzsnJafXy4ODga9eu/fTTT3K5PCUlpbCwcMyYMZb9VSCvCoZhW7ZsuXHjxltvvXXhwoWNGzf6+voCAPbs2ZORkQEAGDZs2OrVq7dt2zZo0KCvv/76o48+GjFihLWj7hC5ufKUFFFurtzagUDsj1aZiFZfMifbCgCA2db2MVESrtPpyAdKpdLw2PCl4uJirVZLEIRer8dxvO3Lf/vtt0mTJgUFBQ0YMGDnzp3t/oigoKBXqlaH2CJ29Nf89deq9977fdWqxw8fyqwdC8QuUalUQ4cO3b59u/Hi+fPnX3vttaqqKoIgNBrNiBEj4uPj277Wjt4pXYypIztDzcLkyZOvXLli/KXs7OxJkyZVVlYCAMhtU9uXT5s27fz5848fP87Ly1u8eHHnaSgE0lEEApZAwMrLk+flwU0SxFxemImA2dYOYqqoIT09/eTJkwCAurq6xMTEU6dOkesEQTx58gQA4O7u3gUhQiCdTkgIVCPIS9cPGzIRH3/88fnz5wsLC+fNmwcAOH78eEBAQFhY2PDhw3fv3p2cnDx9+vSsrKx79+4dPHiwS38lO8eUII0ePfrYsWNKpVKr1apUKkNlMIIg/v7+8+bNY7PZXRIkBNItyM2V5+XJBQJWSAjL2rE4OK9WP+zl5bV+/fqff/55x44dNBpt4cKFUVFRAIA9e/ZMmDAhLCzMkG3duHEjiqIOkG3tYhDCDNvTVatWRUVFhYWFdfqPDw4Ofvz4cad/W4hVsKO/Zm2tZufO8rw8eWSkx+LFfGuH84yUFFFKijg62jc6mmvtWByZq1evbtu2TaPRvPvuu2PHjvX39/f09CQIQiQSlZWVFRQUHDt2TC6XL1y4cM6cOTQardMDsKN3ShdjVmPs999/b+k4IJ2ISqWi0WjGxw44jj958oTP5xv8gboz5Oc++Zgst4MC0H0gCOK7775bsmTJxIkTjZPfCIL06tWLLCH+7LPPLl68uH///oEDBw4dOtSK0XY3zPp4wnH82rVrVVVVZL8RjuMqlSonJyc+Pt4uTu00Go1YLC4oKJBIJAAANze3Pn36+Pn5sViOeTAyefLk2NjYd99917CSnZ09Z86cS5cuBQQEWDEwWyAjo96gRgCA2lpNSooYnpJ1HxAEuXz5Mnm7Vl9fz2KxyA5WksbGxps3b44fP/69996bPHmyOQdIkE7ELEHauXPn/v37Wy36+vq+0DrIRnj8+PHFixefPHlCoVBQFNVoNDwe77333gsJCbF2aJ0JLEIxh5QUUXuLYihI3QfDB9fMmTNdXV0TExO9vb3JFaFQGBsbO3LkSA8PDwRB2vVOhVgOs8xVU1NTR4wYcfnyZU9Pz+++++748eOBgYFjxoyxix2GSCQqLi7W6XQCgWDs2LFjxowRCAQEQRQXF4tE7Xw22S+jR4/GcVyhUBiKUEhaWlr8/f3XrFljF9tZq1Bbq7Z2CH9RW6uxdgjdiEePHn344YeFhYXWDgQCgDk7JBzHa2pqYmNjeTyeQCAoKir64IMPtm3bNn369KVLl7q6unZBlK+MRqOpqKgoLS2lUCiG/YG3t7dSqayoqODz+Z6ensYbdruGwWAcOnQIWLIIxQHw9qa3/cQ3bI+kUi0AQChUSqU6AEBTU2tPLBIezxkAwGbT3NyobHYnJL0NRRYAgIyM+txc+aZNwd7eDvIv05ZZtGjRuXPnoqOjd+zYAd1krM6LBYlCobDZbHIeQe/eve/fvw8ACAgI0Ov1ZWVlgwYNsniMHUOj0UilUiaTaZiwQKVStVptfX29XC53GDUy5vvvv6+trX3w4AHpo5WWlhYREeHi4mLtuGyC6GjfNWtatx9RqUhqak1OznPbkiQSDZ2OsljP3i/GV7LZVDc3Go/nzGZTQ0Nf8f7MoEYkpD5t2hT0at8NYj5cLvfEiRN///vf//73v69atSo0NNTaEXVrzMohDR48eO/evd7e3v379z9y5EhaWlpRURGVSu3Tp4+l4+s4OI4jCCIWi/V6PZ1OZzKZVCoVRVFPT086na7RaBxPk7Kzsz///POBAweSG6YtW7bExcXFxcWNHz/e2qFZn5AQVnS0b0ZGPblPotPRgADnhgatWq0nLyAFhvxf8ikA4Pjx6gEDmAIBEwAgkWibm3VaLdHcjAuFyspK1aNHirw8OYtFTU2tZbOpPJ5zaKgruYsyh4yM+rYtunl58txcOcxsdQFubm6HDx9esWLFd999Z07bUHl5+Y0bNwYOHNiueun1+qtXr1ZUVAwdOjQkJARmoV4KswRpzZo1MTExly9fjouL27dv39dffw0AGD9+vO3nkDAMo9FoGo2mqqpKoVCw2WwcxxkMhqura1BQkL+/v+OpEQDgyJEjXC7XMOzjyJEjmzdv/vHHH6EgSaVaoVApkWh69KBTqYhIpPbwoPn5OQ0c6AoA4POdn6ciJrY+ZBG5XI7370/T6wmpVCeVynNy5K+gTK2AyaQuA8Ow+Ph4Lpd74MAB01fGxsZeunRp2LBhCQkJ3t7eZ8+eNW6lIAji888/v3nz5pAhQ3bu3Pnee+999913Fo7doTBLkHr16pWeni6RSOh0+vHjx69everj4zN8+HBLB9dxRCKRUqnkcrk6na66urq2trapqQlFUYFA0Lt3b9JyysHAcfzGjRtxcXGGGkIej7ds2bIJEyYUFRUFBXXTU6CcHJlQqDQctfn4YDodIRKpIyI4HWxCEghYAoEiL0/eu7fzhAleQqGS/M+gTDyeM4/nHBHx3LFhz8sVwe2RpVm8ePGAAQPIxwiCLFu2jMfj/fvf/37efapQKLxw4UJycvKQIUOamprCw8PT09MnTZpkuODq1as3btzYu3dvZGSkUCgcP358ZGQkOcwFYg7mtkkWFBQcOHCgpKRk4sSJXl5e1dXVFg2rUyDbj2QyWUhISN++fSsqKqRSKY7jzc3NI0aM4PP51g7QIpA5v+zs7ClTphgW8/PzAQDteuA6NlKpNidH/vChjCxSAAAY5CElRZSfrzD9cnMwtsVjs2lsNo3cTpE/OiurgdSnhw9lPJ5zeDinbRFESMgzp1fjRYGABYsaLM17773XamXGjBkmZn2pVKrY2NghQ4YAAFxdXTkcTl5enrEgFRYWstlssjiCx+MFBwfn5eVBQTIfswTp5s2bCxYsoNPpFApFLpfjOJ6QkJCbm7t27VpLx9dB1Gp1c3Nz3759AQB+fn7k4oMHD1gslu2fN74y4eHhp0+f9vf3Hz58OIZhhYWF8fHxfD7fLnJ+nYVBD8inbDZ14EBXwzYlN1eem9sJamQCNpsWEcEJDWWROzOhUCmVyoVCZbvneF9/zTeua/D2xmBFg+VYsmSJ6Xmhhw8fbrd+ODg4ODg4mHx86dIlkUgUHh5ufMGAAQPkcnljY6Onp6dara6srCwvL++8wB0fswRp69atgYGBJ06c+PLLLxEEWbhwoVarTU5OXr16tY33xqIoymAwmpubGQwGudLc3MxkMplMpkOWM5CsWbOmsbFx69athpU+ffrs2LGj++RXMzMbjKUoPJzTKglE7mnIIRQd/FkpKSIT2mbYMxkEkjzHCw1lGe+WSAX68cfyjIx6m7LXc1QM74WGhgaRSNS/f3/j84MXniVkZWWtXLly7ty5rYog3nzzTTab/cUXX0ydOvXatWtKpVKlUrX7HXbt2gWnmLflxYJEzuKLi4tzdn52T4cgyIwZM/bt2/fo0aOBAwdaOMJXB8Owvn37KpVKkUhEbpIAAHV1dT4+Pr6+vo6qRgAADMMSEhJEIlFBQYFMJvP39x80aJCN3zp0FkKhMjW1hjygY7OpU6b4tFtTIBCwoqNBBx2Dams1n32Wa3iakVFv4hsaNkykLJF7JuNNG/hvMgme1FmaHTt2GB7/+uuvcXFxJ06cMN9E9dSpU99++21MTExsbGyrL1EolMTExF27dp05c2b06NEsFut5iWqoRu3yYkGiUqlsNru4uNh48datWwAANzc3S8XVSbBYLE9PT5FI9PDhQyqVSrqOBgQEOGQ5gzH5+fmGnJ9er6+srHz//fetHZRlkUq1qam1QqESmJQikpCQl5Oi1auLoqN9W71k585y46e1tZoffyw/cMCUH1UrWcrKanj4UNZKliC2zKFDh7Zu3bpy5cpPPvmk7Vdra2uFQuH+/fsRBCEIIjIy0i6Kv2wHs47sxo0bd+LECQRBampq2Gx2QkLCsWPHAgMD/f39LR1fB8EwjM/ne3h4iMXP/DQ9PDwcXo3sN+dnAtOzggxndK1yRZaOp9Viba3GnOYhgyyRWyVSlubM6WmxYCGdQ05OzpYtW6ZOndqvX7979+4BALy9vfl8vmFAH4VCWbt2bUtLy9SpUw8ePKhWq409jiEvxCxBWrFiRXV19ZEjRwAA5FaJz+fbS06CrF9weBEyxn5zfibIy5OnpIgFAgX433poqVSblFRFntGFh3MsJEXmlxiY3zzEZtPIbVxGRr1UqktKqqqqsiFLPUhbjh8/ThDEmTNnzpw5Q65Mnz5906ZNhgF9Hh4ey5cv371797p16/r27RsfH992yh/EBGYJkouLyz//+c/8/Pzi4uKWlhYejzds2DBLzK2CdBz7zfmZxtDuk5fHNAiS8cbI9Bldp+PjQ293/WUzQGTFHblP0usJLpfu59d1vwXkpdi0adOmTZvarmdlZRkez5o1a9asWV0YlENhbh/S/fv3k5KSDNsjrVYbERFhycAgr4hd5/xMYNzuA/43Y2S5jZEJvL2xyEiPjIx648VXq5IwnOCRW70HD5pwnIBZJcuxdOlSQ9l3Q0MDAGDmzJnG5z2HDh2ycdtoR8UsQbp48eKSJUsQBOnXr5+zs3NmZubVq1eXLFnyt7/9zdLxQV4B+835mYlQqDx6tApYY2NkTHQ019sbM4z7EwhYHWkeYrNpc+b0NBQ7NDVpp0zx6aRIIf8DiqKGyXvu7u7kHAA4i88WQMz5M4wfP97FxeXnn38m/3IKhWLt2rXp6el37tzp4H0EnC1vCVpaWpYsWXL9+nXDCp/P37Vrl6V9gyz91yRd40aOZGu1BLDSxqgttbWajAxJSoo4Otq3U0ahG7JibDZ1zpyenTLeAmJTwM+952HWPKSKiooNGzYY5gkxmcxFixZdvHixpKSEdNGA2BTFxcU//PBDZWWlQ+b8ysuVAwYwrbgxakWntw0Zb5WSkqpgUXjnQhBEUlLSrFmzjE1R24Lj+IULF0JDQ3k8XpfFBnmxuRmFQgkKCrp8+bLxYlFRkYuLi8McATkYixcvXrFiRf/+/adOnRodHT1y5EgHUCPS6QAAwGBQ58zp2cVqtHp1UW7uc6cldTpkVik8nAMAyMpqyMxs6LIf3R1IS0ubMGFCSkoKmUBqhVKpTEpKevvtt+2lkNiRMHWPcPLkSfIPxuPx0tPTJ06cGBERweFwCgoKrly5MmrUKAe2g7Nf9Hq9SqUyffdnd5BJI3Kc69ChrjZ1imU5WzxDpQPp0Dp3LmxU6gQQBElJSbl8+fIPP/ywcePG0NBQf39/X19fBEHKysrKysqEQiGKop9//nlMTIxhqqeNo9EAiQT8t9kSeHgADw9gjx/PpnJI77//fkFBgYkX37p1y8PDoyM/Hp6lWoIdO3YkJia+8cYbISEhbDbbsD579mxDLbglsNBfMydHdvJktUiklkg0AACBgPX113zb8dchM1sAgM7KIbXi1VJKpvuIIQAAnU538uTJO3fulJSUlJeX6/X63r179+/fXyAQvPvuuz4+Fqwo6fR3Sm4uKCjQ4jju7Ax0OqDRUEJCaCGmPENsFFOChOO46Rd3vMsSCpIlGDduXLtmxh2/gTCNJf6amZkNV65IHj78n+Myb29s06ZgG9Ekg2WD5T79Dcas5msSKZMW0kjHQ6fTabVai96uGdO57xSRCJSU4Gq1JiDg2Yf5kycInY4FBlLszg/A1MHOH3/8MXjwYAc7/OkOJCcne3l5abXaVut2l0lKTa3JyZGLRK39C2prNSkpIhuxxH5ZW7xXgM2mhYayAABkmYM5BR1kH3FurkIggHPQXwyVSrXTDzqNBvz00zmhUOHsXPPfkQZAowEymQ+DIZ8+PWTMmBcPZbcdnlvUQBDEvHnzwsLCFi5cmJKS8vTp064MC9IRPvzww4ULF9LaYO24Xo6jR6vIQeByua7tV2tqutd4b0OZg1SqS02teWGZQ0gIKySEadxKDHFUfH39mpspBMF4+tTn2jXiwgUnoZBfUuIulXrk5dnZ+dNzbwoQBLl3796dO3du3Lhx4MCB9evX8/n8N99888033wwLC+uyvS3kZXGMooajR6uEQiWbTV20iC8Wa+CnKklEBIfNpqam1j58KCOfWjsiiPUJDQ3VaPqrVNqGBn1DQxGLhbu7v+blReHzsZEj7ew21NTHFoPBiIyMjIyMBACUl5dnZWXduHHj1KlTOI4PHTr0zTffnDZtmnHOHGILoCg6c+bMxMTEuXPndnFRQ2dBqhGP50zWlZF3+q2uCQlhWiM060N63+3aJYSaBAEAYBhgsRAEQRob8V69tKWlFWy20ts7wNnZ2d0dVyhQkcieMknm3kfz+Xw+nz937ly1Wi0SiY4fP/7bb78JBIKwsDCLxgd5BdLS0gAAd+/evXv3rvH6tGnTbF+QWqkRACAy0jMjo97YRVsgYHVxrr7deUjWgs2mLVrEg5rUcc6ePTt8+HCLFtRZGo0GyGRAr9drtUhpKaHVYnV1NH9/QBCEry+or9e5u1M8PYG9jCM1S5BkMtnWrVvffvvt8PDwlStXXrp0qUePHj/88MPrr79u6fggr0BaWpper2+7bvtppLZqBP5bU5eRISFlqYOWcY6BaU0i3YxI41fzx2F0Q3bv3r1y5cqRI0dGRUW98847dHr7Du42j55GI/h85NEjql5PBQB4eVEaG4FarScItKkJxzC7GTrzYqcGAMCmTZtOnTrV2NiYnp5+8eLF0aNHu7q6fvXVV+1+6kGsDoVC+fPPP//2t7+NHTt2xIgRn3zyydmzZ+1UjUi8vbHoaG5kpAew0mHdpk1BNrI9MkBqklSqe/hQ1qrG4bPPclNSxKQUZWTUG89Zhxjz66+/Llu2rK6u7ptvvhk1atS333774MEDawf10iAIotUiZWV6Gk2HojoaTfv0qV6nA3o9Sqcjzs4Ujf3ck5glSFlZWYsXL46Kivr3v//t7Oy8e/fu77//XiKRlJeXWzg8yKtw/fr12bNn37p1q1evXqGhoXl5eWvWrDlw4IC14zKFCTWCPI92NenHH8tbXUZWyXd1cPaAl5dXTExMamrq2bNnP/zww1u3bn300Ufjx4/fv3+/RCKxdnRmQeaQUBR1cQG9eyM0mjOGEfX1RG0tIhQCJpPm62s353XAHEHCcbyhoaFv3744jt+8eXP48OF0Op00DdLp2qnHhVid+Pj4Pn363L59OyUl5eDBg/fu3Zs8efLevXtf2OlsLVJTa8iaOqhGL0tbTWrXc89C5kYOg6+vL4/HI31URSLRnj17xowZk5KSYu24XgyZQyIIgsGgFBRgVVW9Hz7sn5+P/f47+uefSHU1akcVDcBMc9WePXuePXs2KSlJKpWOGTNGLpfv27fP09Ozb9++pl9bXl6elJSUk5PTSdFCXoxOpyspKZk9e7ahvs7JyWnBggXNzc2mjaAsRG6uPCVFZMKZNDOzgew3WrSI34VxOQ7GmkROLISYiUajSU9P/8c//jFq1Kh169Y1NTWtX7/+P//5z40bNyZMmBAXF2cn+yQ9k0kIBCiFQsNxjbOzks8Hvr46AHC1Wi+yq72xWUUNixYtWrVq1eXLl4OCgqZMmbJjx44zZ8589dVXpq1wY2NjL126NGzYsISEBG9v77Nnz9p7c4xdQKVSGQxGK2MS8qmLi0vXx5OXJ09JEQsECgBA2zRMTo7M4IjT9bE5DGw2bcIEr7S0utTUGhaL2raQodtWyZtm2rRpxcXFHh4ec+bMmTZtmvEd9tSpU1NTUyUSiaenpxUjNAcEQZhM8OgR0tKiwzCRp6dKIPBRqxGdjqirQx8/dnK0Krv3339/yJAhNTU1ISEhTk5OkyZNmjVrVkBAgImXCIXCCxcuJCcnDxkypKmpKTw8PD09fdKkSZ0UNsQUb7311unTp9lsdkREBIvFyszM/OWXX3g8Hp/P7/pgSA+bvDx5Xh6zlSAJhcrU1FoAwJQpPjZl4G2PDB3qRqMhqam1DAaFTkfV6r8KjsiSECvGZrO8+eaba9asCQwM9PLyAgCkpaVFRESQ922DBg26cuWKn5+ftWN8ARgGAgORP/+k3b2r0unULi6EWk1/9Ajt2ZMqEOhzcnSVlXozawVsAVOCdOzYMZ1ON3v27F9++UUmkwEA7t+/b3zB/Pnzn9fXolKpYmNjyfF9rq6uHA4nLy8PClLXsHz5cpFIlJiYmJiYSK5wudz4+HgUtcK/y5AQVrsGNoYx5HPnmjXcyHJTHszBpvqQnkdoqKtUqsvKaggOZvTs6ZSd3QSr5E3z9ttvf/755wMHDjx06BAAYMuWLXFxcXFxcePHj2cymUymfWwrFQpAoRA0Gk6l6hBEBwDAMJwgKDodQhCETmdPo9lfIEgKhWL27NlJSUntFtTNnDnzeYIUHBwcHBxMPr506ZJIJAoPD+9wtBCz8PT0PHr06O+//15cXKxSqXg83ujRo21qsotUqn0pNQIAkKomELAEAptWBetCNiRlZTXodPrISI+UFDE8rDPBkSNHuFzu0qVLDU83b978448/jh8/3rqBvRCjD2QsNxeoVERwMBAKiaYmCoVS07NnQ0uL08OHTgShRRB5ebmGPLKzyhnJS2FKkM6cOUM+SE9Pb3dKhTnjFLOyslauXDl37twRI9o3nd21a9eiRYvMCBViFhqNRqVSubq6hoWFhYWF5eXl+fj42JQaAQDIk7rwcI75g18FAlZ0tAVHPJjGjjYZoaGspiZtTo780SNYWWcKHMdv3LgRFxcX8t/BQTweb9myZRMmTCgqKgoKst2/eFxc3OHDh//7DNPpQqjUXk5OYq12EI3W08mp4V//uiCX++l0NBy/S6VmAfBsct+8efPWrVtnpajNwtwqg1cb5Xvq1Klvv/02JiYmNjb2eddANepEzp8//8MPP3z88ccLFiwgV3bs2HH79u358+fHxsbayDxmQ8vRS3nedMGUB8eAzaaRpuBV9oouHgAAIABJREFUVSprx2LTUCgUNpudnZ09ZcoUw2J+fj4A4IWH2w0NDWfOnImJiWm1Xl5e3qqW1dvbu9Mdba5fv37o0KExY8YAADQa8OABUVen6dEDr6oCJSVauVyv0yG1tRRPT+epU6PeeOOvV8XFxdmxIJ07d66xsdHEBTNmzDBx633o0KGtW7euXLnyk08+efUAIWbzn//8JzY21sXFxTgTO23atJKSkgMHDjg5OVlF+1NSRMa5H7LlCDbAWhQ2mzZlind+vkIkUkulsFnwuYSHh58+fdrf33/48OEYhhUWFsbHx/P5/D59+ph4lVKpXLRoUWlpaVtBunXr1u7du8nHBEE0NDS8/fbbFrVYI4sa1GpqTQ0REKD38aFJJKCqCunThz5gADpokOV+skUwJUg//fRT27GGKIoaHIOmTp36PEHKycnZsmXL1KlT+/Xrd+/ePQCAt7e37Z9g2jWJiYmenp6nTp3q0aOHYXHSpEmRkZHz5s07ePDgZ5991pXmqrW1mtWrHxtKkDMy6qVSnUSigQ2wXUBlpQpFEQCAUKgk7wCsHZEtsmbNmsbGxq1btxpW+vTps2PHDhNnCQ8ePFi+fHlFRYW7u3vbr86aNWvWrFnk45MnT8bFxX355ZedHnYrWCzg7k6prtYXFQEqFSiVwN2d0q8fao8jzE0J0okTJ8je/j179hw7duyHH34YOXKkVqu9fft2YmIik8l0c3N73muPHz9OEMSZM2cMiajp06dv2rSpc6OHGFNYWDhx4kRjNSJxcnKKior69ttvy8vL+/Xr12Xx7NxZbtwQU1uruXq1PjiYMWWKHZsr2wtkDYifn5NeT6Sm1pg5+Ly7gWFYQkKCSCQqKCiQyWT+/v6DBg2iUExZkZKtLB9//LGhhLVdxGLx5s2bFy1a1AXvOAwDfD7w8KCJxc/+xB4ewL4MGgyYEiTD7ufMmTOffPLJ2LFjycVx48ax2ew5c+Y8efLked1ImzZtgvLTxaAo+rzxVIGBgV0cTG1tO1P11Gr9W295wLv1LsBQA9LQoCUHn0MjjLYQBPGf//yHLEYFAFRXV5PHOfPmzXveWcJ3331Hp9NPnz5t+jtv377dy8vrs88+e94FnVvMxWIBFsteRciYFxc16PV6rVbbqspOKpVaLKTORKPRSCSSiooKqVQqkUjc3NzYbLafn5+Hhwdpx+dIDBs27MaNG0uWLGn7pTt37ri6uhoK8buAmhp1u+u1te2v2zh20YdkjKEGRCrVkkV3qak1cG/aij179uzatavtenR09PMEyZwRFVKpND09ffny5SY2W7CYq11eLEgoio4aNSopKYnP54eFhSkUihs3bvzyyy8BAQGmzRpsgcePHz98+PDJkye1tbUajYZCofj5+fn6+g4ZMiTEHk9YTTJq1Ki4uLiEhIR//OMfxu+Ey5cvHzhwYPjw4V3ZGOvj0/771tvbTkfO2Ctk0Z1QqMzJkbu50eA0P2NSUlLCwsI2b97co0cP47xRB98pp0+fplKpUVFRHQ7wJdBogEQCxM8KvIGHB/DwAHZ3121W2fe33367ePHiFStWGFZ69+79z3/+02JRdQ4ikUgikWi1Wi6Xy2azfX19RSKRi4sLQRASiUQkEnEdYItrxMyZM/Py8vbs2ZORkTFs2LBevXqJRKKioqI7d+5wudyvvvqqK4Px9saio31TUsTGiwIBi5xpZHfYUR9SW9hs2pw5Pclpfny+MzwyJcFxvKmp6f333+/Zs5NLbNLT0999990uNnp4/BgUFGhxHHd2BjodePyYEhJCs7u7brMEicPhHD58uKCgoLCwUKVSDRgwoH///jY+XVGj0Zw/fz4zM1OhUBAEQaPRnJycdDpdS0uLs7MzjUbjcDjTp08na/kdAxRF/9//+3++vr7p6enJyclkMaSvr29UVNSqVatcXV27OJ7ISE8AgEGToIeNFSELwVNTa2GBgwEKhRISEpKVlTV9+vQOfqvjx48HBASEhYWRT0tLSydPntzhAF8CkQg0NuLu7rqAgGe5lSdPiMZGVCSi2Nddt7mNsRQKRSAQCAQCi0bTuej1epVKhaKoUqkkM0ZUKrWlpYXJZMrlchMlgvYLgiBffvnll19+2dzcXFJS4ufnx+FY7Yjmv2NePTMyJNDDxuoYnO5ggYOBcePGbd68+dNPPx0wYIBx0shEUYMB4yO+PXv2TJgwgRQkiUQik8m60uhBowFiMZDJdH5+f2X6e/Ykysp09fUUO7L6BuYLkj0yaNAgNpvd3Nwsl8u9vLycnJyUSmVdXZ27u7uTk1Pfvn0dL41kgMFghIaGWjsKAADw9rafd4OjExrKItuSMjMbYDIJAECO4Lt58+bNmzeN12fMmGFakKKiooxTRFlZWYbHnp6ebds3LY1Wq29uJhiMv1YwDCgURFMTjmGmqthtDYcVJAzD+vbtq1QqhUKhq6trXV2dn59fY2Mjh8MhCKJXr16+vr6YHd05QCAdhjy427VLmJXVAJNJAIBLly61u24VX/yOgCAIg4E0NwODJjU3AyYTcXamaDT2tEOys//fXwoWi+Xp6YkgSG1tbVlZ2ZUrV4qKiiorK2k0mpubm4NVNEAg5kBqEgAgNbVGKtVaOxwrgz4Ha8f1cmAY8PNDcJx64wZaWAgKC4FYDCoqEA6H5utrT2oEzN8h3b9/Pykpqbi4GADA5/M/+uijiIgISwbWCWAYxufzPTw8SktLm5ub9+/fz2az33777R49ekA1grwUdteHZAJDMik1tRZ6OOXn5x84cKCkpGTixIleXl4UCuX999+3dlAvTXU1KCkhnj7FRSIcANDcTOHxaHw+anefc2YJ0sWLF5csWYIgSL9+/ZydnTMzM69evbpkyZK//e1vlo6vg7BYLBaLRcrP3bt3AQDDhw+3dlDdC+sO1oO0C0wmkdy8eXPBggV0Op1CocjlchzHExIScnNz165da+3QXgKRCCgUeN++uhEj0NpaCgBAJEK8vCh214QEzBSkH3/8sX///j///DPpJ6hQKNauXZuQkDBr1qyuLyaGvJBWhigGzKkd6nQcY7CegxWsG5JJ3bwzaevWrYGBgSdOnPjyyy8RBFm4cKFWq01OTl69erVpRzvbwVBix+PpGQzg7U0AAIKCQFkZkMvp9pVAAuYIEo7jFRUVGzZsMLjbMpnMRYsWXbx4saSkhBxSDrEpXsEQxXJYd7Ae5HkYdyZ1zypwnU5XXFwcFxdneFMgCDJjxox9+/Y9evRo4MCB1g3PfBymxA6YI0gUCiUoKOjy5csffvihYbGoqMjFxcXf39+SsUFekVczRCkuLr5165a/v394eHir28OOjB2Dg/VsltBQ15wcuVCo7J42d1Qqlc1mk3lxA7du3QIA2FeTIoIgdDoiFAKl8tkKjQaoVPsrsQOmBenkyZMNDQ0AAB6Pl56ePnHixIiICA6HU1BQcOXKlVGjRjmeP6kD8GqGKPv379++ffsbb7xRWFg4fPjwhIQE4692/dgxSNcwZYp3UlJVTo48NNS1Gx7cjRs37sSJEwiC1NTUsNnshISEY8eOBQYG2tGtNjmgLy+PevOm1tlZi2EAx0FTE23AACe7K7EDpgUpOTnZ+L64tLS0tLTU8PTatWsKhcLGDYSAkeG3SCSqq6s7d+6cAxt+g1cyRKmurt69e/e2bdumTJlSWFg4derUBw8eDB482HCBVcaOQboA0nq121oKrVixorq6+siRIwAAQwmx6QF9NohCAaqrEZGIcHMj6HRCrUakUsLT09phvRKmBOm3334z/WK7yPsZDL+LiopUKtWlS5cc2PCb5GUNUfLz81EUnThxIgDgtddeCw4OzsrKMhYkA105dgzSNYSGupJe4FlZDd3t4M7FxeWf//xnfn5+cXFxS0sLj8cbNmwYjWZPqqzRgIICAIB64EACRWlNTSiDAZycCLlck5PjZF++QcC0INmF3pjGYPjt6+vr7OzM5XL79evnwIbfJC9riPL06VMvLy8q9dk/hh49etTV1bX7nV84dsxRcaQ+pLYY5lPweM6hod2obva3335TKP7qSSguLib3STNnzjSMJ7Vx1GpQXa1vbtb5+YG6Oh2NpqfTCQxDHj4k3NywyEg7a/I1q+zbpsqIzUej0YjF4pqaGgCATCYjQ/Xy8hKLxb6+vjU1NZ6enp6eno5nIPSyhiiNjY3Gbz9nZ+fm5ua2l5kzdqxz52BCugY2mzZlis/Ro1VZWQ08nnP3Objbt29fZWVlq0U3N7eoqCh7ESQEAXo9kMtRhULLYum4XB0AQKVCystRuRwXi1Eez9ohvgxmCZJNlRG/FGq1urm5mUajtbS0kP/CSMNvpVKJIEhDQ4PjqREAAEXR+/fvJyYmFhUVKZXKoKCgadOmffDBB8+7nsPhGN8nyuXydmfLmjN2zFHVyMH6kNrC4zmHhrJycuTdyr7h3Llz5JQWAIBKpSovL9+wYUNwcLAdtVdiGOjRA0EQamWlZsgQnFyUSNCePRFnZ319PbCv0gazBMlCcxW7ABRFGQyGQqFwcXFRqVROTk4qlcrFxQUAQKfTmUymRqNxPE26fv36F198gaLooEGDXFxc7t+/n52dLZVKn3fU1rNnz9raWplM5urqShBEaWlpZGRk28usMnYM0mWQB3dCoTInR9ZNDu6M76cZDIaHh8eGDRtmzJixfPlyTzupCsAwMHAgcusWNTeXgqI0FAVaLdBqkcBAGpNJtLTYWSvSixXFuIyYQqHYkQUhafjt4+NDpVJdXV2VSiUAgDT8BgA4sOF3fHx8nz59bt++nZKScvDgwXv37k2ePHnv3r04jrd7/ejRo9ls9oYNG0Qi0fbt2+Vy+dixYwEAx48fv3fvnuGy0tLSAQMGdNHvAOlyyIM7AEBWVkM39119+vSptUN4CVxcgE6HKpVoURFaXo6IRIhKhVKpSK9ez1qR7AizGmM7a65i10MafldUVNTV1SmVSolE4uHh4e3tHRAQ4KiG3zqdrqSk5P/+7//YbDa54uTktGDBgvPnzxcUFLQ7YhHDsC1btnzzzTfnzp3jcrkbN2709fUF1h47Bul6DAd33aTi7uTJk3K5nHys1+vlcvm//vUvZ2fndo+sbZacHEAQGk9PwsmJcHfHdTqUxQI4Tuj1Dur23ZG5itbF2PD7yZMndXV177zzjpeXlwMbflOpVAaD0WpEGPmUPKtsl9GjR5Pms8ZYfewYpOvpVhV3P/30k1AoNF7x9vbevHmzjX+sGSMUgtJSnMVSv/kmUVeHoCjS1ITI5ahcjgQEoHZy7vgX5uaQwCvNVbQFDIbfV69e5XK57733nrUjsjhvvfXW6dOn2Wx2REQEi8XKzMz85ZdfeDwen8+3dmgQW8fQKpuV1eDwgnT+/Hnjc2wMw+yr14V0VlUotCyWnscjvL2R+nqMwwFiMVWrpbFYegyz9cRKK8wSJLueq0g6NYjFYpFIBAAoLy93VI8GA8uXLxeJRImJiYmJieQKl8uNj4+3i7+XbeLYfUit6D4ed+fOnTOuL20FgiDR0dE23ier0+kVCuDkhGZnIyoVCgBgMgmlkmCzcQ6H5jhedgRBJCUlzZo1i0qlmvggw3H8woULoaGhPFstdyedGrRabWNjo06ny8zMdGCPBhJPT8+jR4/+/vvvZOsYj8cbPXq0vfRVQGwBcjiFw3vcJScnP3r0iHyMoqihBNzA+++/b5uCVF5eXl5ertEAiYTW0kJ78IDS3Iy7umpQlFCpaAoFdeBAVKGoFon+5yVWC9dsEIIg2v0CQRCzZs2SSCQxMTHjx48nK9OMUSqVp06d+vnnnwEAR44ceTU7wuDgYItmJkQi0ePHj+VyOZ/P379/PwBg3LhxLBYrODjY8XJIx44d0+l0s2fP/uWXX2QyWdsL5s+fb9EjVkv/NSFdSU6OLDW1ls2mOvBwipSUlPj4+KVLl44aNcrT0/POnTt79uwBABw5coTsb2EYD3XoPDr4TnnrrbcM6qJWB2o041xcBri7K+n0RgBwpbKnQuGiVucBcJROrzd+IZ/Pv3btWkcitzTP3SEhCJKSknL58uUffvhh48aNoaGh/v7+vr6+CIKUlZWVlZUJhUIURT///POYmBgbufs+fPhwqxWhUNjQ0MBkMnNzc7Ozs319ff/444/GxsYePXqQhWStmDdvXhfEaSGOHTumUChmz56dlJTU7t3QzJkzbT/nB7EReDxnHs/ZgafKEgSxZcuWhQsXRkdHkytjx451dXX9+OOPa2pq+vTpY93wTHPo0KExY8ZoNODuXXDhgrahodnDQ00QhFYLCIJCo9H8/EaNHPmVQPDXkd3169fj4uKsGvWLeUEOady4cWPHjj158uSdO3fy8vLOnz+v1+t79+7dv3//Dz744N133/XxsZUj5vnz55eXl7fK29fX10skkh49egAASAV6+PChRCJhMpltIy8vL8/MzDx06FBXhdzJnDlzhnxw9uxZKpVqnJ5FEEStVttXwhZiXYynyoaGshzPT4ggCCqV2uqYrra2FgDQyiPtFWhoaDhz5kxMTEwHv48Z6CkUvZ8f6uZGIyvYXVxI/2/c7rpigTlFDVQqNTo6mryJ0Ol0Wq3WZu+yP/nkE+Mtjkaj+eOPP6qqqrhcrmHr3dzcLBKJ+Hx+SEhIq67Yw4cPZ2ZmdmXAFmLy5MmxsbHvvvuuYeXevXtz5sy5dOlSQECAFQOD2BdsNs2B/YRQFB01atTBgwd9fHyGDRvW1NR08+bNlJQUHo/XQTN7pVK5aNGi0tLSLhAkFEU8PBCpFGlsxGk0HMNAczMiFtPc3SlUqmMN6GvnairV4Alt+5BODUqlUiQS9e3bl1ysq6vz8fFxSI+G9PT0kydPAgDq6uoSExNPnTpFrhME8eTJEwCAYQg9BGImBj8hoVDpeNUNcXFxixYtWrlypWElMDBw586dHalHffDgwfLlyysqKrrg7YZh4LXXkCdPqKWlWhoN4fMJAIBKhWq1FCcnCp1uZ2oEXlaQ7A7SqUEkEj18+JBKpapUKhqNFhAQ4HgVDQCA0aNHHzt2TKlUarValUplqGdFEMTf33/evHkG7wYIxEyMJ/g5XnUDm80+fPjw48ePCwoK1Go1n89/4403OniynZycPGTIkI8//tjQdGFRMAwgCCqRUJqaqOXlAEGARoMyGKhGo6dSgYPvkOwOg1ODWCwmVzw8PBxSjQAADAaDTIAtXbr0iy++sC/7ExunW/UhtcKx25IoFEr//v379+/fWd/wu+++o9Ppp0+f7qxvaBo6HTg56X18AI9HISvUmUwCw/CmJl15OSUkxOFySACAs2fPDh8+3HbqF14Kg1ODtQPpOv74448ff/xx79691g4E4gjk5sqbm3G5XOeQB3f5+fkHDhwoKSmZOHGil5cXhUJ5//33O/IN6XS6OZd11uQwtRooFIiHhy4kRMfhPOvhkcmQBw/0LS0UuRzYlweAWUelu3fvHjNmzKeffnr+/Hm1Wm3pmCAdQa/Xq1QqO0r12QWbNgV1z+0RACAvT56WVieT6aRSXWZmg7XD6Uxu3rw5ffr0q1evVlVVyeXy6urqFStWbNy4sQt+dGdNDqPTAZcLdDpqeTm1qOjZfyUlVI2G5uKity81AmbukH799ddz586dOXPmm2++YbFYEydOjIqKGjx4sKWDg7wCKIrOnDkzMTFx7ty5ISEhxnmj2bNn22yFJMRmEQhYAoGiqUnreJukrVu3BgYGnjhx4ssvv0QQZOHChVqtNjk5efXq1XbUI+HiArRa6v37uEajxzBCo0FQFOXzKVwu6pg5JC8vr5iYmJiYmMLCwrNnz165cuX48eN8Pn/69OnTpk2zl0lW3Ye0tDQAwN27d1sZeE+bNg0KEuRlCQlh5eXJU1LE/fszAQAOU92g0+mKi4vj4uIMbwoEQWbMmLFv375Hjx4NHDjQuuGZCYYBd3dErwd6Pe7jo8MwAseBWEyj0RAUtTM1Ai9b1ODr68vj8Xg8XmVlpUgk2rNnT0JCwpo1awytzhBbIC0tra0rFwDANl25IPYCj+fMYFAcxruBSqWy2ezi4mLjxVu3bgEA3NzcOv79jYdrWw6NBpSVAb1e16OHjsvVUChAowE0GqLRaKVSB+1D0mg0165dS01NzczM1Gq1ISEh69evnzRpEgBg48aNcXFx77zzzvP2SV3Ysdw+Go2msrLyzp07OTk5crncxcWlb9++YWFhffv2dVTPbwqFYkcHDhB7gc2mjhzpfvRolcN4N4wbN+7EiRMIgtTU1LDZ7ISEhGPHjgUGBr6aM6cxUVFRUVFRnRKkadRqIJEQCKJ74w1Eo3kmPr17gz//1FVW6s2sErAdzBKkadOmFRcXe3h4zJkzZ9q0aYYmUwDA1KlTU1NTJRJJu4LUlR3Lz+Px48eHDx++e/euXC5HUZQgiNu3bxcWFn766acO7Pl9//79xMTEoqIipVIZFBQ0bdq0Dz74wNpBQeweg8GdY4yUXbFiRXV19ZEjRwAA5FaJz+fv2LGjazY3nQKCAL1eT6EQCPKXTbZOB3Q6Qqdr3zjbljF3YuzXX3/91ltvGe67cRx/8uQJn88fNGjQlStX/Pz82r6qKzuWn4dIJPr9999ramq8vLwEAoG7u3tDQ0NNTY1EIvn9998dtSfp+vXrX3zxBYqigwYNcnFxuX//fnZ2tlQq/eyzz6wdmr3SnfuQWkEa3JWUtDhAdUNxcfEPP/xQWVlZXFzc0tLC4/GGDRtmXyfbGAb8/NDCQurt25oePbQYBnAcVFXRXFyoXK79HZOYtaFLTU3V6XTGp0DZ2dmTJk2qrKxkMpn+/v7t3lCQHcurVq3qtGBfEo1GU1ZWlp2d3dTUxGAwSDs7Nzc3V1dXhUKRn5//9OlTjUZjrfAsR3x8fJ8+fW7fvp2SknLw4MF79+5Nnjx57969xsMxIZBXg/RuaG7GHaAEfPHixStWrOjfv//UqVOjo6NHjhxpX2oEAMAw4OODYBhVJqOgKIUgKDodpbmZQqPRfHzsLIEETO+QOuiNZmbHcmc1iLWLRqORyWQqlYrD4ZDmdRQKhSAImUwmFouVSqXjOdrpdLqSkpL/+7//MxR8Ozk5LViw4Pz58wUFBQKBwLrh2SmbNgVZOwQbIjSU9fChzN5LwB2jY0+jAS0twN0dUKm0xkZUpwNNTYiTE6pWE3V1oKEBtJlkZ9OY+mN00BvNzI5ly6kRAIBCoTAYDAqFUltb+/TpU51Oh+M4juNUKpVc12g0DqZJ5K/WavwX+dTFxcVKQUHsmNpaTW7u/8z5dgyDOwfq2NN7eOhHjED/+INSVaXz8NBzOLqaGv29e4SvL+Y4gmTwRlu1alVUVFRYWFhXRdU5YBgmEAjGjBnz+PHjsrIyAICTk5NWq1UqlRwOh8/n9+7d28HUiOStt946ffo0m82OiIhgsViZmZm//PILj8drNSwKAnkhKSmilJRnPpAZGfUAgOhoLjCqbsjJkYWGulozxA7gGB17Oh2CokhODtLcrGcw9IGBehoNeHoCtZqoq8NFIoodJcrN2q5+//33lo7DQrBYLIIgDLsiHMcRBMEwDMMwnU5n7egsxfLly0UiUWJiosFvmMvlxsfHd8RUH9INycioN6gRAKC2VpOSIhYIWCEhLDabFhHBOXq0Kiurwe4Eac6cOfPnzx87duz+/fv5fL5Wq211ge1nksrLyw1ToUUi7NEjl9JSXKfDcZy4c4fi4kJ57bVmHq9OKKTrdIhCocEw0O4UaVvDlCCdPn26sbHRxAXR0dE2fh+h0WikUimHw3FyctLr9TiO02g0DoejVqvlcvnTp089PT0db5Pk6el59OjR33//vbi4WKVS8Xi80aNH28iYeYgdkZIiam9RTFYbGjZJducC/uDBg7CwsLFjx3722WetRlnaBWPGjImLiyPnkavVHgD0pVJHI8hrdDoHQRoJwk2lomVlCV1czmq1zihah6LP9n/G80ttE1OCdPjw4cLCQhMXTJ061RxBsm5RP1nF0K9fP+MDqzt37ohEIocsaiDJyck5efJkQUEBjuMDBgyg0WgRERHWDgpiZ9TWvqAGldwkSaU6qVRrR32yAoFgz549ly5dqqmp2b59+8GDB1tdcPDgQVdX2932rVu3bt26dQAAjQbk5oKCArVYrCsqUuO4urGxFkXrXV0H0WgBwcFve3pSBAIsJMRuyu1MCdLJkydNFwqbc9PdZR3LJmAymTKZrLy8nJQfvV6PIIiHh4dDFjUAAM6dO7dixQoAQEhIiI+PT05OzoULF1atWjV37lxrh2avdM8+JIGAlZcnb7Xo4/PX+8VO+2Tj4+N3797d2NgoFArJGZ6tLrCjxtiWFn1lJaHV4lyunkolFAqcyVS5uREajfbpU4TDoXt62o0aAdOCZPik/u233wwldsbMnDnTpg6Cjhw5kpmZ2WqxpqamqKhIo9EwGAzyN1IoFCiKenh4VFVVGbsHlZeXO0baPykpaeDAgfv37yfv8vR6/cqVK3fu3BkdHW37h+MQ2yE62nfNmtaCNHash/FTsk/WvkrAuVzupk2bAAAxMTHz5s0LDw+3dkSvDooiNBqoqyPc3Qk2GwkOHqjRIE1NoKEBpVIBlUrweNYO8WUwq6hh3759lZWVrRbd3NyioqJsSpD4fH7bg6n6+nomk1lSUqJSqXQ6HUEQFArFy8srLCwsNDTUKnFaFBzH8/Lytm3bZjhzQFF0wYIFZ8+ezc/Pd8hfuQvonn1IISGs6GhfQ12Dtze2eDG/1TaRzaaFhrJycuQ5OTJ7ESQDP//8s7VD6BAYBvr0QR4+pKpUWioVcXJCeDxKWRnVzY3m7Iy6u9P4fDvzVzVLkM6dO2dwj1apVOXl5Rs2bAgODra1Y9aIiIhWWTuNRpObm+vu7o7jeFlZmVwuBwB4enq6u7u/+eabAoGg1Xld2w2W3UGhUAIDA8+fP0+635IAwQ/aAAAgAElEQVTcuHGDQqEEBARYMTCIPRIdzY2M9ExJEWVk1EdGerR7aBkezsnJkefkyENDXe1Ok+wdDgeEhFDFYqykhKiro7i66p2cUAxDGQwqn0/r08ee1AiYKUjGlQsMBsPDw2PDhg0zZsxYvny57Q9DUqvVVCp16NChb7zxhmHxwYMHcrnc8bJHJJGRkXv37h0/fnx4eLi7u/vt27f/+OMPPp+fnJwMAEAQJCYmxlF/d0in4+2NeXub+tdC9slmZTVkZjbMnduzywKDAAAwDPTvD+Ry+tmz1OJiLZWqp9EIV1ekb18an4/aUQcSSYdsM8iy6c4KxULgOK5UKjMyMlQqFQCAyWSS5QxMJtMhKxoAABcuXAD/26kAACgtLd25cyf5eNasWQ75i0OsRWgoKyurwb4ySQAAlUpFo9GMXToNttF2ZCnE4QAOBwwcqOdwELWaiuMAw9ABA6jBwdaO7OUx6//0kydPkoddAAC9Xi+Xy//1r385OzsH2/xvjGEYgiB//vlneXk5m80my+rodPo777zj6+vrqB/Kly5dMn0B7JCFdC5sNm3KFG+7MxOaPHlyqz6k7OzsOXPmXLp0yY7Ot0UioFbjffvqxo17Nm/iyROETkckEnvyaCAxS5B++uknoVBovOLt7b1582Yb74oFAAiFwvz8/JqaGrVarVAoSOsgvV5fX19v7dAsCNQbSNfD4zmz2VSpVGf7m6QO2kbbFBoNEIuBTKbz8/tr+lHPnkRZma6+nmJfNd/ATEE6f/68cUMShmF2MZBUo9HcvHnz/PnzdXV1DAajpaVFLpczmUxvb28cx6urqx3SpoEEDujrXLpnH9JLYUeOqx20jbY1tFp9czPBYPy1gmFAoSCamnAMs4MPamPMEiQMw2pra6uqqgYPHgwASEtLi4iIsH3raLFYXFFRoVAofHx8/Pz8FAqFTCZTq9V1dXUIggwbNsxR1QgO6INYBcMmycYdV+3dNroVCIIwGEhzMzBoUnMzYDIRZ2c7q/kGZg7oy87OnjBhQkJCAvl0y5YtY8eOTU9Pt2RgHUWj0VRUVDx9+tTJyYmc1FdfX4/juFwur6ysbGlpQVHUIafzATigzwJs2hTUbbdHubnyVuMnnge5SQIAZGXZx+y+77//3t7VCMNAYCDi6kp9+vSvD/PaWoTDofn62pkaATN3SEeOHOFyuUuXLjU83bx5848//jh+/HhLxtZRCILAMIxOpysUCiqV2rNnTwCAUqkkD+6YTKZD7pDggD5I55KXJ29rIPQ87GssBY7j165dq6qqIg2/cRxXqVQ5OTnx8fF2dGrHYgF3d0p1tf7PPwGVCpRKQKFQevSwv5pvYI4g4Th+48aNuLi4kJAQcoXH4y1btmzChAlFRUVBQbbbwY6iqLu7u0gkUqvVKpWqpKREp9M1Nze7ubk5OTkRBOGQZd9wQB+kcxEIWNHRzx688GLSuIF0t2srSLm58rw8OTnAwhKhviw7d+7cv39/q0VfX1/TOfLi4uJbt275+/uHh4e3e+WjR4+ys7N5PN6oUaO6oHwcwwCXCwiC9vDhM2Mwf39gj2oEzBEkCoXCZrOzs7OnTJliWMzPzwe2Xc2FYVi/fv3q6uqePHnCZrN79uzZ3Nwsl8s9PT0HDhzIYrEcuDEWDuiDdCIhIS+nHyY2SXl58pQUcXQ0sBFBSk1NHTFiRFxc3KxZs5YsWRIYGLhmzZphw4YZW1y2Yv/+/du3b3/jjTcKCwuHDx9uSGQYOHr06NatW8PCwgoKClgsVlpaWheUgIlEoLBQS6HgaWn/0umQyZM/cHWleXi8+IW2hlmKEh4efvr06YMHDz569Ki4uPjcuXPbt2/n8/l9+vSxdHwdgcVi9ezZk8FgaLVakUgkkUiUSiWDwWCz2T4+PmRjrLVjtAiff/75kCFDEhMTo6OjJ0+evG3bNgRB4IA+SNdAbpJAe5kkgYAlELBycxW5ueaeAVoOHMdramqmT5/O4/EEAkFRUdGgQYO2bdt2/PhxmUzW7kuqq6t37969bdu2o0ePHj16ND09/cGDB8YXyGSyHTt2fPPNNz///PPly5clEgnZpW5RRCLQ2Ii7u+uGDdP36FHbq1eNl5eusREXtTPNytYxazu5Zs2axsbGrVu3Glb69OmzY8cOGzdpJzdJn376aXJyckFBgZubG4PBcHd3Lysr43A4DtwY+9lnnw0bNiwpKQkO6INYhedtkkJCWOQmKS+PafVNEnn209zcDADo3bv3/fv3AQABAQF6vb6srGzQoEFtX5Kfn4+i6MSJEwEAr732WnBwcFZWFll7THLt2jWCIObNmyeVSjEMu379uqXPybtjHxKGYQkJCSKRqKCgQCaT+fv7Dxo0yC5akVgsFofD6d+/P0EQnp6eVCpVrVY3NjY68PxyvV6vUqnUanVYWJi9VxDZDrAP6aUwnUmyHQYPHrx3715vb+/+/fsfOXIkLS2tqKiISqU+7+zn6dOnXl5ehrRQjx496urqjC8Qi8UcDufTTz+9desWAGDmzJnkJD2L0l36kI4dO6bT6WbPnv3LL78Y72HFYvHdu3cBAPPnz7dxswaNRiMWi5VKZa9evcjxE25ubsHBwTiO19fXO2RjLIqiM2fOTExMnDt3bkhIiHGx0OzZs2387wVxGHg8Z8NYCpvVpDVr1sTExFy+fDkuLm7fvn1ff/01AGD8+PHPyyE1NjYanzQ4OzuTGywD9fX1VVVVw4YNu3Tp0u3bt9evXx8SEtJuT/quXbsWLVrUKb8F2YfU2IhotaC62hsAcP8+2tKCEoT99SG9QJAUCsXs2bOTkpKMbToNzJw50/Y/4EpLSx8+fOjn50en03U6XVNTE5fLbWlpaWhocDw1IklLSwMA3L17l7xvMDBt2jTb/3vZJt1zHlJHYLNpPJ5zTo5cKFTarCD16tUrPT1dIpHQ6fTjx49fvXrVx8dn+PDhz7uew+EYjyqVy+Wt/DzZbDaKoitWrOBwOHw+/+TJk7dv325XkDpLjcg+JLWaev++Tq/XyeUMOp0oKtKgqFPv3namRsC0ILFYrK+++goAsH79+hEjRhAE0eoCG88hAQDEYnFTU5ObmxudTi8qKiosLFQqlefPnx8wYICrq6tDln0DANLS0gzzq4yB42IhXQmZScrJkfN4zjarSQUFBQcOHCgpKZk4caKXl1d1dbWJi3v27FlbWyuTyVxdXQmCKC0tjYyMNL6Ay+UiCGK47SMIwkTBXmfBYgEcpzQ06PR6nYsLrlYDb28Cw3AKBReJ7Mxf1ZQgPXjwICwsbOzYsWvXrm3liWubCIXC69evG68UFRUJhUKRSHThwgWpVKrX61EUValUVVVVSqWylYViKwNZ+4VCodhFhg/i2Nh+JunmzZsLFiyg0+kUCkUul+M4npCQkJubu3bt2navHz16NJvN3rBhw9KlS3/99Ve5XD527FgAwPHjxwMCAsLCwsaNG7dp06bNmzcvW7bs6tWr+fn5nbUTMg2VCvh8QKfTystb6HQgENBcXYnKSvurazAlSAKBYM+ePZcuXaqpqdm+ffvBgwdbXXDw4EHbGRq7bt26+fPntxIkhUIhlUq1Wm1TUxOCIOS/PBqN1tLSkp6enpmZ2aqxtwsykBZFJpNlZmZWVFQEBga+/fbbUJYg1sW43M7X1ykjQ5KRUQ8AqK21iY6LrVu3BgYGnjhx4ssvv0QQZOHChVqtNjk5efXq1e2+dzAM27JlyzfffHPu3Dkul7tx40ZfX18AwJ49eyZMmBAWFsZgMHbv3r1s2bJjx465uLh89dVXpGJZGhpNT6MRgweD+/clAABfXwCAXdY1mBKk+Pj43bt3NzY2CoVCFovVdhafTR3Z8fn8a9euGa9oNJq7d+9evXr14sWLTk5OOp3O2dm5R48eWq2WIAiFQjFu3LjY2FgPe+wfaw+RSPTJJ59UVFSQT0NDQ3/99Vc7mjMGcTwMm6QrV+qN25JIWVq8mG+1yADQ6XTFxcVxcXGGEzYEQWbMmLFv375Hjx4NHDiw3VeNHj26VWoWAJCVlWV4/MYbbxg/7RoM/qoG7NRf1dSn1YoVK+bPnz927FiyjXnYsGFdFlangGFY375979+/T1qpkgd0GIYxGAy1Wi0WiyUSicOoEQAgOTm5oqJi6dKlr7/++smTJ8+cOXPq1KmPPvrI2nFBujXkJik9XdJqPSOjXiBgRUZa7Q1IpVLZbHZxcbHxIlmu7ebmZqWgXgVDXcPTp3+l+e3UX9VU6/6DBw8ePXoEAKitrbXTiXYYhpEjzKVSqVgsbmhokMlker0eQRAmk4kgiJ3+Xu2Sk5MzfPjwL774YujQoevXr6fRaKTDE6TjrF5dZAvmAvYIuUmSy9vp/CP3SVZk3LhxJ06c2LRpU01NjVgsTkhIiI+PDwwM/P/t3X1cVNW+x/EFyPhAGCAoTipoOaMBYqAICWpkWol2rNSsY3lMr1pHr5lmGpXVsaOeLD2aJfYqxSOZZJmGStdnTZOTkonAFSEUeZTkSWRAhn3/2DVnLg8DKsMs7PP+o9ew9561foOv+O7Ze+21evToYdvCbpQ6v2p5ucOlS13y8jr/+9/2RUVtnJ1b3/yqt889pHplZ2eXlJR07dr13LlzBoOhTZs2VVVV1dXVv/7669133z1mzJjb6RvS1atXTbfE2rdv7+vra3nIENAyGlpAtksXG5+9L1iwIC8vb+PGjUII9auSt7e3/HPQ1KXRCG9v0amT4759ZUIIna5dp06tcn7V2+ceUl05OTlpaWkajWbAgAHV1dUZGRllZWWVlZUXL17s1KmTXq+PiIiwdY3NyWg0mt8xcnJyYgGk5sJzSLfCxcXR19e57hoWTZk+3Ko6dOiwbt265OTktLS0a9eueXl5DRw4sJU+HeHsLJydhVabK4QIDLR1NTfLUiBptdp3331XCDFlypTJkycPGTKkpapqBuoCfenp6Y6Ojj4+Pl26dDlw4EBKSoqrq6tGoxk1atT06dNtXSPwRzFnjvfUqWfMt3TurLHhDSRz9957b+/evU+cOGFvb19dXd1KA+n2YFf3cdeWpNfra63c01yqqqrmzZu3du3aNm3amM9yrSiK0Wisrq5WFGXx4sWtfZy3udGjR1+5cqVXr17qj6mpqUKIPn36mA748MMPrXqJ1Xr/mrgNJCSUbNx4KSvLIITw9XW27ZfO3bt3b9261Wg0Tpo0yc/P709/+lNRUZEQomvXrlFRUdZe5q3Z/0+pqhKFhSI3V6jLOy1c+F+dOgnrP5Lb/G5mLjsTyeeyGz9+vFarVRSle/fuap0GgyErK8vBwSEsLCwgIOA2m6bB29u7urq6sPC34UzqJVbTj0II25584A9Op+vg6dnWaFRycir9/O6wYSWHDx+eO3eui4uLi4vLnDlz7r77bgcHh8WLF5eXl69evXrZsmV175dL7n//V6SkXDcajRUVd1RX2x07dt3Pz/H3FVVbk9t2Ljt17YmioqLExMTLly+rw2bUkyCdTtejR4/bLI2EEKtXr7Z1CUCD1OF22dkGWxcidu3a5eHhERcX5+zsPG/evJ07d7733nujR48WQhQXF0dHR9u6wBtjWg+pZ08lMbFACOHhUV1UZN/q5g0SlgNp+/bt6ov4+Ph6T64lH9SgLtCn3rHMyMioqakxGo09e/bs1auXttX9QwGt35AhbsnJV3NyKm1bRkZGho+PjzrL3ODBg3fu3Nm7d2911913311ZWanOVmfTGpvqj7gekhCioqIiKyvLzc3N3d1d8hwyUR+M7dix4/nz59UrV3feeWefPn28vLxsXRpaGdZDahbqFODHjhXbtoyKigrTsixqLJmmQFXHqV6/ft1Wtd2ETz/d+D//k+TgcEoIkZub27Vr1507dxoM/WpqcsPD23/22We2LvAGNBJIRqNx3bp1sbGxubm56pek9u3bR0REzJw586677mqRCm+Js7Ozs7PzPffcY+tCAEjEdFbdWk6vLXj00UcHDBju6WlwcvrtS1J5uV1eXjt3d7uRI1vZpSBLgXT9+vVJkyYlJiZ6enqOGDFCr9cXFRWlpqbGxsbu379/8+bNPXv2bLFCARviOaTbzPHjx2fOnCmEUJd8jYyMVFfey8/Pt3FlN0ijEcHBHmfPGsvKrmu1vy0688svdnp9W29v+9Z1vU5YDqSvv/46MTHx6aefrjU2+ujRozNnzly5cuWqVausXB6A28eZM2Vnzlxt/Dgr8/DwKCoqOnPmjOlH8xntPDw8zB8UkZ86b1BeXs1PP4k2bURFhXBwcPD0bH3zBgnLgbRv3z5vb++6S4OEhoY+/vjjO3bssGZhaGlpaWnHjh3r0aPHkCFD6p17/+zZsydPnvTy8ho8eDCTiOMmJCWVJSWV+fo623aOhk8//dSGvTc707xBubm/PdLbSucNEpYDSX3Kst6ThXvuuefatWvl5eVOTk5Wqw0tJyoqasWKFYMGDUpNTQ0ODv7nP/9Z64Do6Ojly5cHBQWlpKQ4Ozvv3r2bxZZs7+Bb4sJBUZwpvIcJ/+eE9zBbF9QIX1/niROFr68zw0Oa1+/zBtm6jltmKZAMBkOtNVVN1DEqBoPBciA1etINGeTl5a1Zs+Yf//jHmDFjUlNTH3vsscTExPvuu890QGlp6QcffPDyyy//5S9/KSsrGzp0aFxc3JgxY2xYM8TGB0Tmwd9e/7RBZB4UQ98U/SfbsKJG+fkRRbDEipdKo6KiIiIi9u3bt2DBgpdeesl6HeEWJScn29vbP/roo0KIPn366PX6WiuMHThwQFGUyZMnFxcXOzg4HDx4UD34FsnzJK88lagaryfz4H/SSFWcKQ69ZbN6bNEUbj+NBNKhQ4cm1Wft2rWW32g66Y6Ojo6Ojo6Pj09MTGy+stGcLl265OHhYbot5OnpqQ49MsnNzXVzc3v++eeDg4MDAgJWrFjRLHd916xZ08Qjm/JXrNFjLBzQaCUtvB5S47+ZzEOml6vP/D5FaXGmKM60aj23Hie1Phr5BHOWLtnVmhutll69elm4s13vSbf5VSDIo6ioSB3zqmrfvn25+WLIQvz666/Z2dkDBw7cs2fP8ePHFy9e7Ofn9+STT9ZqJygoSK/X31DXTT++KenV6DEWDmi0km3bGu2/OTXhN/OfkehrTJkUM9La9TT9NKLRpoQQQUFBt9gabieWAulWTl4aPemGPNzc3K5e/c9g3LKyslp/DV1cXOzt7RcsWODm5ubt7R0bG3v8+PG6gbRp06aWKBdCiMyDYuMDtTe6eIv//sUW1QDNw1qDdxs96VbdxDl1K/LXv/511qxZtq6icXfddVdBQYE6f5eiKOnp6Q8++KD5AVqt1s7OzjSRrqIozq1xavvbifcwMXSxOLT4/218rDVNEvNH1nr/7ln7b5q1AqnRk24V59QyCA0NdXFxefvtt+fOnRsTE1NWVhYeHi6E+OKLL3r27BkUFDRixIh333136dKl8+fP379/f3JycqsI2tvcsDeFEL8N+3bxFkPflH/YN1T83WuItQKp0ZNuyEOj0Sxbtuzll1/euXOnVqt95513unbtKoT48MMPH3nkkaCgICcnpzVr1syfP3/Lli0dOnSYPXu2mliwsWFvCnH7rDAJWGvF2KqqqgceeCAkJEQ96Y6JiYmLi1P/zAEAUJcVlzA/evToyy+/XFxcrNVq582bN2rUKCt1BAC4DVgxkAATy3N22HyWvCtXrmzfvn3KlCkt37WFAjIzM1NSUsy3dO7cOTAw0NrFZGZmHjlypF+/fv7+/rV21dTU7N+//+LFiwMGDPDz87O8dkNycnJCQoJOpwsJCal75IULFw4fPjxgwIC+ffs28wdAq0UgweosT5Rn81nyKioqpk6dmp6e/sMPP7Rkv40WEBMTY3roR1GUK1euDB8+/MMPP7RqMfPmzduzZ8/AgQOTkpI6d+78zTffmE4RFEWZNm3a0aNHAwICkpKSRo8evWTJkobaWbBgwXfffRcQEHDy5Mnhw4e/99575nvVR+YHDhx4+vTp+++/n8dj8RsFsKbc3Fw/P79vvvlGUZSUlBSdTnfq1CnT3pKSkv79+3/66aeKopSWlt53333qkS3m1KlTw4cP1+l0gwYNasl+b7SArVu3+vj4JCcnW7WYzMzMPn36nDx5UlGU4uLifv36ffvtt6a9e/fu1el0e/fuVY/U6/X79u2rt53U1FS9Xv/jjz8qirJv3z6dTpeVlWXaW1JS4u/v/+WXXyqKkpiYqNPp0tPTrfq50Fq0pmU/0BpZnijPSrPkNd3mzZsDAgIWLlzYkp3eaAG5ublLly6dNWuWta9uGQyGefPmBQQECCE6duzo5uaWlJRk2puamuri4jJs2DAhhJeXl16vN99rTqvVRkdHq1cXr127JoQwGo3mB1RVVakvNBqNEOLKlSvW+DhodVjVBtZlec4O0yx5x44dE0JMmDCh1mqQ1rZkyZK2bdt+/fXXLdnpjRawYsUKDw+PqVOnWrsYvV5vel5wz549OTk5Q4YMMe318fEpKysrKipyd3evrKzMysrKzMystx1nZ+egoKDCwsKlS5ceOHBg4sSJXl5epr0dO3acO3fu66+/HhcX9/PPP48YMYJJxaDiGxKsy/KcHeoseR4eHnv27HnzzTe3bNny1VdftWR5bdu2bcnubqKA4uLi+Pj4Z555piVvrR0+fPjVV1999tlnQ0JCTBvDwsJcXFymT58eHR09Y8aMiooKg8FgoZHy8nJXV1dXV9eTJ09evHjRfPuuXbtcXV21Wq2rq2tOTk6rWzgcVkIgwbrqztnh6elp+tE0S563t/fEiRN9fHyOHz9uizLl9fXXX7dp02bs2LEt1uOXX345Y8aMSZMmvfbaa+bbHRwcPv74Yzc3t+3bt/v5+T300EM9evSw0I6Xl9drr70WFxeXl5f35ZdfmrYfPHgwJSUlNjb2b3/7265du3Jzc7/44gtrfRi0Klyyg3VZnrODWfIaFR8f//DDD99xxx0t091nn322fPnyV1999bnnnqu1q6Cg4MKFC1FRUXZ2doqiPPjgg8HBwfU2cvTo0djY2Pfff9/BwaFt27ZarfbChQumvRkZGZ07d9ZqtUIIR0dHnU5n/v0Jf2R8Q4J1mSbKy8nJWbFihflEeQkJCSNGjHByclq6dOnVq1d37NiRnJxsftPiD0v95aiv09PTfXx8Wqbf06dPL1u2bMyYMX379k1ISEhISFDvEqn1ODg4vP7661u3bjUYDGvXrq2srHz44YfrbUev1x84cGD9+vVlZWWff/55amqqOhRCbSc4ODgvL2/z5s0Gg+G7775LSEgYP358y3xAyM6mY/zwh3DkyJGgoCCdTjds2DDTMOKwsLB3331XUZQffvghLCxMp9P1799/zZo1Nqnwq6++Cg4OtknX9RZg+uVcvnxZp9OdOHGiZcpYuHCh7v9buHCheT2bN28ODQ3V6/URERE//PCDhaa2bds2atQonU7n4+OzcuVKdaOpnQ0bNvj7++v1+r59+7711lvW/2RoHXgwFi2hpqbm8uXLXbp0aeiAgoICNzc3m0zTgBtSU1NTVFTUqVOnxg8VIjc3193d3dHRse4uo9GYn5/fqVMnm48rgTwIJACAFLiHBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQT8xmg0ZmZmlpaW2rqQ5pSXl/frr7/e3HuvX7+elpbWvPUAFhBIkNq9996r/11AQMC4ceNiY2NvpcHk5OTPP/+87vajR48GBgaOHDly3rx5t9K+VAwGw7hx486fPy+EeOCBB6ZPn17rgHo3miiKMmPGjIMHD1q1SMCE5xAhNUVR/P39R40aZTQai4uLDx06FBkZ2b59+4iIiJtrcM6cOQEBARMnTqy1fdOmTW3atPnoo4969ux5y1XLYvPmzerSf+L3OVlqHaAoSk1NTUNv12g0s2fPXrly5dChQy2vVg40C74hQXY6ne65556bMmXK3LlzY2Ji2rZtu2/fvlrHNPR8d9Of+y4uLvb19Q0PD7cQSLfeS1Pe0lzPql+9ejUqKmrSpEm30khEREReXt7u3bubpSTAMgIJrUmHDh2cnJxMJ/WlpaWLFi0aPHjwfffd9+c///ns2bPq9l9++WXmzJmDBw/29fV95JFHoqOj1e2zZs3Kzs7eu3dvrdk8J0yYkJKScvr06YiIiNzc3IiIiI8++igkJGTUqFHqXaV6exFC5Ofnz5kzJygoKDw8PCYmZvTo0YcPHy4oKIiIiDhw4IB6TF5eXkREhLpObkNNFRcXR0RExMbGTpw4sV+/fg899JBpRYZr16698cYb4eHhISEhixYtSk9Pj4yMNJ+KOyMjw9S+uejo6OvXrw8ePLgpv9jPP/88oo7q6moHB4fw8PBVq1Y1pRHgFnHJDrLLycnZv3+/EKK0tHTv3r3Xrl17+umnhRCKokybNu38+fNPPfVUly5dNm3a9PTTT+/evdvd3f2ZZ5654447Zs+ebW9vv3Xr1iVLlgQGBvr4+IwfP/7kyZPe3t611laYMmXKe++95+TkNHXq1C5duqSlpa1cuTIiIqKiosLb23vChAl1e9FqtQaDYcaMGfn5+TNmzCgpKVm+fHlFRUVJSYk6FsA0OKKqqiotLa2kpKShgrVarfqWyMjIMWPGjB07dsuWLW+88UZAQECvXr1efPHFn3766cknn9TpdB9//HFaWtoTTzwRGxt78uRJdY3wr7766tKlS+q64+aOHz9+//33m88jl5OTs2XLFvNjTIsl+vv7myYSPH/+/IYNGwYNGqQuCTh06NBt27bl5eWZL2QFWEWLT+cK3IA+ffrUmn/6qaeeKiwsVBTl8OHDOp1u3bp16pGXL1/u27fvO++8k5mZ+corr/z73/9Wt//44486nW7Hjh3qjw899NCCBQvqdjRu3Ljp06cripKbm6vT6V566SV1e0O9KIpy5MgRna7731gAAAm7SURBVE53+PBhddfq1avVji5duqTT6bZv365uv3DhgrrdQlMFBQU6nW7+/PnqrtOnT+t0um+++ebnn382b+rIkSPTpk2rqKjw9/ePjIxUFKW6ujo0NNT0RnNBQUHLli0z/Th06FBdfaZOnWr+ruzs7NDQ0EcffVRNUEVRUlJSzD8mYD18Q4LsRo0apUZIWVlZcnLy3//+9wkTJuzdu/fcuXNCiF27dh09elQ90s7O7vz586aFSr/99ttffvklOTlZCFFdXX1Dnfr6+qovGupFCJGUlCSEGDBggLo9LCxs9erVFtq00FStTtXvIgaDQS3etI54aGhoaGioEGLkyJF79uyJjIxMSEgoKCh44oknavVVWFhYXFzctWtX841DhgxZt26d+RZ1bSqT0tLSadOm1dTUrF+/vmPHjupGdSW9c+fOhYWFWfh0wK0jkCC7Dh06qOtWeHp69u7dOzs7e9WqVWfPnjUYDEKIwMBA01IIwcHBHh4eOTk548aNKy8vDwwMDAoKevjhhxcvXqzc4EiBO++8U33RUC/i95ArKipSV7x1dXU1b8FoNKovKioqGm1KZVo5197+t5u76urv6qUzc0888cT27dsPHjz43XffdevWLSgoqNYBau+1PrWdnZ2p5boqKytnzpyZnZ29efNmNYTq/TiA9RBIaGUKCwuFEO3bt1eHw3Xv3n3y5MlCCIPBsHTpUk9Pz7i4uMLCwi+++KJ///5CiD179gghzAc3WxjoXFdDvQghevXqJYQ4dOiQOoj8+++/V9+i0WiEECUlJeqP6enpjTbVEG9vbyFEUlLS0KFDhRAJCQkvvPDCqlWr7r//fi8vrx07dhw/fnzKlCl1x2R36dKlU6dOGRkZTfyYNTU18+fPT0xMXLNmTa0FatX6W2zVWvyREUiQ3blz5zZu3CiEuH79empq6rfffuvv79+zZ0+tVqvVatevX29vb9+/f/8NGzbs2rVr9OjR6pebQ4cOeXh4pKSkvPXWW0KIa9euqa21a9fuzJkz8fHxI0eObErv4eHh9fYihBg+fHiPHj0+/vhjjUbj6OgYFRWlvsXNzc3FxSU6Orp3797l5eVLlixptKmGhIWFde/effHixTNnzuzWrduKFSs8PT1DQkLs7Owef/zxDz74QH1R73t9fHyaHkgrVqww/U5Mo+oDAgJcXV3VRkyXEwHrIZAgNTs7u9OnT58+fVoI4ejo2Llz5/Hjx8+ePdvOzq5du3ZRUVGLFi1S/+Lfc889r732WmBgYE1NzejRozds2LB27VpnZ+cXX3zxk08+OXXq1LPPPiuEeOyxx95///3Zs2d///337u7u5h2p3zPM/yuEaKgXIYRGo1m3bt1LL720aNEie3v74ODgnJwcIYSDg8Pbb78dGRn5/PPPt2vX7oUXXlCTw0JTtTo1laTRaD766KNXXnnl9ddfd3BwCAkJmTVrlnrZbezYsatWrRo0aFDdy2sqX1/fmJiY6upqdfic6QPW6kLdeOLECSFEfHx8fHy8ae/GjRuDg4OTkpK6detmuoYJWJGNB1UAt6ykpCQ3N7fWxsrKyqysrJqamrrHV1ZWlpaWNksvqvz8/JKSEtNoOnWj0Wi8ePFiVVXVDTXVkCtXrly5csV8S0ZGhk6ni4uLa+gtWVlZvr6+27Ztu6GOaikpKenfv//69etvpRGgifiGhFavY8eOpiFhJhqNplu3bvUer9Fo1Ns8t96LqnPnzkKI4uJi84329vbdu3e/0aYaYj5iIj8//8SJE1u3bvXw8BgxYkRDb+nWrdvkyZM/+eSTsWPH3vTEP//617/c3d1rPbYFWAkzNQDNw87OzsHBoe6IuGaXl5c3f/78s2fPLl++3PQ0a73UJ3brzrTURAaDYdOmTQsWLDB/uhawHjulmSbOAtBi8vPz3dzcmpITly9fdnR0dHFxuYleqqurL126pI70A1oAgQQAkAKX7AAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFL4P8inxcm/91WeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure; set(gcf,'color','w'); colormap([.5 .5 .5; 1 1 1])\n",
    "\n",
    "% 1. Difficulty\n",
    "h = subplot(2,4, 1); hold on\n",
    "    pos = get(h, 'Position');\n",
    "    set(h, 'Position', pos +[-.08 0 0 0] )\n",
    "    boxplot(zeros(10,1), 'colors','k','Widths',.01)\n",
    "boxplot(Dif,'boxstyle','outline', 'colors','k','OutlierSize',4,'Symbol','','Widths',0.6)\n",
    "scatter(ones(size(Dif)).*(1+(rand(size(Dif))-0.5)/10),Dif, ...\n",
    "    'MarkerFaceColor',[.1 .1 .1],'MarkerEdgeColor','k',...\n",
    "    'MarkerFaceAlpha',.2,'MarkerEdgeAlpha',.2)\n",
    "    set(gca, 'FontSize', 10, 'FontName', 'Arial')\n",
    "    ylim( [0, 3] )\n",
    "    ylabel( 'Difficulty (nb dis/beat)' )\n",
    "    xlim( [.0, 2] )\n",
    "    set(gca,'XTick', [] )\n",
    "    set(gca,'Layer','top','Box','off','TickLength',[.01 .01])\n",
    "    \n",
    "    \n",
    "% 3. boxplot individual preferred Fq\n",
    "h = subplot(2,4, 4); hold on\n",
    "    pos = get(h, 'Position');\n",
    "    set(h, 'Position', pos +[.04 0 0 0] )\n",
    "    boxplot(zeros(10,1), 'colors','k','Widths',.01)\n",
    "boxplot(Fit, 'boxstyle','outline', 'colors','k','OutlierSize',4,'Symbol','','Widths',0.6)\n",
    "scatter(ones(nsuj,1).*(1+(rand(nsuj,1)-0.5)/10), Fit(:,1), ...\n",
    "    'MarkerFaceColor',[.3 .3 .8], 'MarkerEdgeColor','b',...\n",
    "    'MarkerFaceAlpha',.2,'MarkerEdgeAlpha',.2)\n",
    "    set(gca, 'FontSize', 10, 'FontName', 'Arial')\n",
    "    set(gca,'Layer','top','Box','off','TickLength',[.01 .01])\n",
    "ylabel( 'Beat frequency (Hz)' )\n",
    "    set(gca, 'YLim',[.4, 4], 'YScale', 'log' ) %!\n",
    "   yticks( [ .6 1 1.7 2.9 3.8] )\n",
    "xlim( [.0, 2] )\n",
    "    set(gca,'XTick', [] )\n",
    "\n",
    "\n",
    "% 2. plot performances\n",
    "subplot(2,2, 1.5); hold on\n",
    "x = bsxfun(@plus,bsxfun(@minus, OUT{i0}, mean(OUT{i0},2)), mean(OUT{i0}(:)));\n",
    "errorbar(beats, mean(x), std(x)/sqrt(size(x,1)), ...\n",
    "    'b.', 'color', [.3 .3 .8], 'MarkerSize',15, 'Linewidth', 1.5,'Capsize', 1)\n",
    "\n",
    "% + plot staircase\n",
    "plot(2, .6, '.', 'color', [1 .5 .0], 'Linewidth', 2, 'MarkerSize', 15) %!\n",
    "        \n",
    "% + plot Fit\n",
    "x = mean(OUT{i0},1);\n",
    "    x = polyfit(log(beats), x, fit_order);\n",
    "    freqpass = exp((-x(2)-sqrt(x(2).^2-3.*x(1).*x(3)))./(3.*x(1))); % order 3\n",
    "    x = polyval(x, log(xVec));\n",
    "    h = plot(xVec, x, 'color', [.3 .3 .8], 'Linewidth', 1.5);\n",
    "    h.Color(4) = 0.7;\n",
    "\n",
    "% + plot max fit\n",
    "[~,b] = max(x);\n",
    "maxpass = max(x);\n",
    "h = line([xVec(b), xVec(b)], [.0, max(x)]);\n",
    "    set(h, 'color', [.3 .3 .8], 'LineStyle', ':')\n",
    "hold off\n",
    "set(gca, 'FontSize', 10, 'FontName', 'Arial')\n",
    "ylabel( varname{i0} )        \n",
    "xlabel( 'Beat frequency (Hz)' )\n",
    "set(gca,'Layer','top','Box','off','TickLength',[.01 .01])\n",
    "ylim( [.6, .8] ) %!\n",
    "yticks( .5:.1:1 )       \n",
    "    xticks( round(beats([end:-2:1 1]),1) )\n",
    "    set(gca, 'XScale', 'log')\n",
    "    xlim( [.5, 4] )\n",
    "\n",
    "% Print figure\n",
    "cd(origDir)\n",
    "% print('figure1b','-dpng','-r600')\n",
    "f.PaperSize = [46.81 33.1];\n",
    "f.PaperPosition(1:2) = [0 0];\n",
    "print -dpdf figure1b.pdf\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "optimal beat frequency:  mean = 1.34 Hz, med = 1.12 Hz, std = 0.802 Hz \n",
      "difficulty level (staircase):  mean = 1.01 , med = 0.85 , std = 0.700  \n"
     ]
    }
   ],
   "source": [
    " fprintf('optimal beat frequency:  mean = %1.2f Hz, med = %1.2f Hz, std = %1.3f Hz \\n', mean(Fit), quantile(Fit,[.5]), std(Fit))\n",
    " fprintf('difficulty level (staircase):  mean = %1.2f , med = %1.2f , std = %1.3f  \\n', mean(Dif), quantile(Dif,[.5]), std(Dif))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Matlab",
   "language": "matlab",
   "name": "matlab"
  },
  "language_info": {
   "codemirror_mode": "octave",
   "file_extension": ".m",
   "help_links": [
    {
     "text": "MetaKernel Magics",
     "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
    }
   ],
   "mimetype": "text/x-octave",
   "name": "matlab",
   "version": "0.16.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}