exp5 figure4&S3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TempAttV exp5 (visual - passive&tracking)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load & preprocess data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "clear\n",
    "clc\n",
    "% close all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "% define Paths: replace XXX by your own path\n",
    "origDir = '/XXX/exp5/';\n",
    "addpath(genpath(origDir)); cd(origDir)\n",
    "\n",
    "% define data\n",
    "SUJETS = 1:20;\n",
    "% SUJETS = []; % to remove the unfinished ones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "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(['TA18V_' subject '_Listen_1*.mat']);\n",
    "        load([ifile(end).name]) \n",
    "        IN{isuj,1}  = X;\n",
    "        Dif(isuj,1) = cfg.ndis;\n",
    "    subject = ['S' sprintf( '%02d', isuj ) ];\n",
    "        ifile = dir(['TA18V_' subject '_Tracking_1*.mat']);\n",
    "        load([ifile(end).name]) \n",
    "        IN{isuj,2}  = X;\n",
    "        Dif(isuj,2) = cfg.ndis;\n",
    "    clear X subject ifile\n",
    "    \n",
    "    if Dif(isuj,1) ~= Dif(isuj,2)\n",
    "    fprintf('Problem! : difficulty level was not matched across sessions')\n",
    "    end\n",
    "end\n",
    "Dif = mean(Dif,2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "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": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "% extract data\n",
    "OUT{1} = NaN(SUJETS(end), length(cfg.fbeat), 2); OUT{2} = OUT{1}; OUT{3} = OUT{2};\n",
    "for isuj = SUJETS\n",
    "    for itapp = 1:2\n",
    "        for ibeat = 1:length(beats)\n",
    "            y = IN{isuj,itapp}.beat == beats(ibeat);\n",
    "            out{1} = IN{isuj,itapp}.iscor(y);\n",
    "            out{2} = IN{isuj,itapp}.rt(y);\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,itapp) = nanmean(out{1});\n",
    "                OUT{2}(isuj,ibeat,itapp) = nanmean(out{2});\n",
    "                OUT{3}(isuj,ibeat,itapp) = nanmean(out{2}) ./nanmean(out{1});\n",
    "            clear out y rtmin rtmax\n",
    "        end\n",
    "    end\n",
    "end\n",
    "clear IN  isuj ibeat itapp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "% Load tapping data\n",
    "load([origDir 'exp5_tapping.mat'], 'TAPP'); % TAPP\n",
    "load([origDir 'exp5_lag.mat']); % LAG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "% remove phantom (NaN) participants\n",
    "OUT{1}(isnan(OUT{1}(:,1,1)),:,:) = [];\n",
    "OUT{2}(isnan(OUT{2}(:,1,1)),:,:) = [];\n",
    "OUT{3}(isnan(OUT{3}(:,1,1)),:,:) = [];\n",
    "TAPP(isnan(nanmean(TAPP(:,1,1,:),4)), :,:,:) = [];\n",
    "Dif(isnan(Dif)) = [];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "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": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "% Fits per participant\n",
    "xVec = (beats(end):0.01:beats(1));\n",
    "Fit = [];\n",
    "% Fit perceptual (visual) data\n",
    "for j0 = 1:length(OUT{i0})\n",
    "    for k0 = 1:2\n",
    "        x = OUT{i0}(j0,:,k0);\n",
    "        x = polyfit(log(beats), x, fit_order);\n",
    "        if fit_order == 3\n",
    "            Fit(j0,k0) = 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,k0) = exp(-(x(2)./(2.*x(1)))); % order 2\n",
    "        end\n",
    "    end\n",
    "end\n",
    "\n",
    "% Fit motor tapping data\n",
    "z = nanmean(TAPP(:,:,2,:)*100, 4);\n",
    "for j0 = 1:size(z,1)\n",
    "    x = z(j0,:);\n",
    "    x = polyfit(log(beats), x, fit_order);\n",
    "    if fit_order == 3\n",
    "        Fit(j0,3) = 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,3) = exp(-(x(2)./(2.*x(1)))); % order 2\n",
    "    end        \n",
    "end\n",
    "Fit = real(Fit);\n",
    "nsuj = length(Fit);\n",
    "clear j0 k0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ans =\n",
      "\n",
      "    0.2462    2.9685\n",
      "\n"
     ]
    }
   ],
   "source": [
    "% detect potential outliers (4std)\n",
    "for j0 = 1:size(Fit,2)\n",
    "    olr = find((Fit(:,j0)>(mean(Fit(:,j0))+4*std(Fit(:,j0)))));\n",
    "end\n",
    "% control range of Fit fq (should be between [0, 3.2])\n",
    "[min(Fit(:)), max(Fit(:))]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAB3RJTUUH4wsVDiAyG94NGAAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAyMS1Ob3YtMjAxOSAxNTozMjo0OS03EgUAACAASURBVHic7N17PFT5/zjw94z7NfdbxVCNS0ko0QqpVUm1uqcU2lt9Pt21RTe2bbfoouimTxtJkdp0xyaiVCyFQoRRDIZyZ5jb74/z3fOYn1vSMHN4Pf/Yh3nP+5zzmvaYl/f7vC8kHo+HAAAAAGEjCzsAAAAAACFISAAAAEQEJCQAAAAiARISAAAAkQAJCQAAgEiAhAQAAEAkQEICAAAgEiAhAQAAEAmQkAAAAIgESEgAAABEAiQkAAAAIgESEgAAAJEACQkAAIBIgIQEAABAJEBCAgAAIBIgIQEAABAJkJAAAACIBEhIAAAARAIkJAAAACIBEhIAAACRAAkJAACASICEBAAAQCRAQgIAACASICEBAAAQCZCQAAAAiARISAAAAEQCJCSEEPr06dOff/4p7CgAAGBYI/F4PGHHIGRtbW3ff/99cXHx8+fPhR0LAAAMX8O9hfTy5csFCxb8888/wg4EAACGu+GekCIjIy0sLHx8fIQdCAAADHfDPSEdPHjw8OHDI0aMEHYgAHyV4OBgYYcAwNcSF3YAX6C9vf3du3dFRUVFRUUdHR3jx48fP368gYGBmJhYv88pJSX12Tru7u5WVlYbN27s91XA0DMQd+PXCAkJgVsUEB0xEhKHw7l27VpwcPDHjx8RQhISEmQyub29HSEkIyOzcOHCzZs3q6ioDNDV09PTIyIiBujkgHCEezcCMIQRICE9evQoMDCwo6PD1dXV0dFRV1dXTU2Nx+PR6fSSkpL8/PyoqKh79+5t2LDB3d1dQkJC2PGCoQzuRgAGEE+0cblcR0fHO3fucDicnuqw2ezbt2+7uLhkZGT07yp//fXX1KlTe3qXSqX277RgiBmcu7F/4C4FQwAB5iFxOBysX57JZEpISPD30XM4nNLSUgqFIi4ujn0eMrk/wzRu3rwZEBDw7Nmzbt81NDR8+/Zt/4IHQwx+N/bua+7G/oG7FAwBBEhIuFmzZnl7e8+ZMwcvSU9Pd3d3j4uL09fXH7jrwq86wN25c6eurq6XCsuWLZOWlh60eHBwl4IhgADPkOLj42NiYhBCNTU1Z8+evX79OlbO4/FKS0sRQsrKysKMDwwn58+f7/q9TyaTuVwu9vPChQuFkpAAGAIIkJBsbW2joqLa2tpYLBaTyWxubsbKSSSSrq6uh4eHkpKScCMEw8e1a9c4HA5C6NSpU1FRUUeOHJk2bRqLxXr27NnZs2fl5eVhThsA/UaAhCQnJ3fx4kWEkI+Pj6urq5WVlbAjAsMX3vqJjY1du3ato6MjVujk5KSkpOTu7l5aWjqgHcgADGEESEi4P/74Q9ghAIAQQlwul8VidXr+Wl9fL6x4ABgaiJSQOBxOUlJSRUUFi8XCXjKZzOzs7GPHjkGvHRhMZDL5m2++iYiIoFAoVlZWzc3Nqamply9f1tfXh+YRAP1GpIQUFBQUGhraqVBbW1tYi7WA4Wzfvn2bN2/euXMnXmJgYHDu3DkhhgQA0REpId2+fdvGxsbf39/NzW3r1q1jx47dvXv3lClTFBQUhB0aGHZUVFTCwsLy8/MLCgqYTOb48eNNTEw+uzRiXl5eeno6lUq1sbEhkUifvUpKSgqDwViyZAle8tdffxkbGxsbG3/tBwBA9BBmtW8Oh1NdXb148WI9Pb0JEyYUFhZOmjQpMDAwOjq6sbFR2NGBYQfrQM7MzKyvr2cymenp6WfPnvXy8urlSdLOnTtXrVqVmpq6YcOGHTt29OUqUVFRHR0d+MsHDx74+PhkZWUJ4AMAIHoI00ISExNTUlJqaWlBCBkYGGRkZCCE9PX1uVxuSUnJpEmThB0gGF6+tAP57du3t27dioyMtLS0fPTo0fr167ds2TJq1KheLsFms58/f75r1y7sZXV19f79+wUSPACiiTAtJISQubn56dOnHz16ZGJikpeX9+DBg9DQUHFx8TFjxgg7NDDsYB3ICQkJampqBw8ejI6OHjt2rIODQ08dyDo6OpcuXbK0tEQItba2IoSw+Uy9yMrKUlNT09XVRQjxeLxdu3bNmjVr5MiRgv4oAIgKwrSQEEK7d+/28vJKSEjw9/c/c+bMli1bEEKzZ8+GZ0hgkGEdyN7e3ngH8pIlSwIDAxcvXrxt2zZFRcWuhygoKFhZWdXW1h46dCgpKWnlypV6enrYW6mpqffv38fXekAIeXh4GBsbp6Sk2NraYiXh4eHv378PCQmZP39+T1EFBwfDlkiA0IiUkEaNGhUfH19bWyslJRUdHf3o0SNNTU1ra2thxwWGnX53ILe0tCgrKysrK2dmZr5//15XVzcrK6u0tNTZ2Tk7O1tRUVFCQkJPT2/s2LEIoSdPnmzevBkhVFRUFBQUdOHCBTk5uV6igmwEiI5IXXYIofz8/D/++GPBggWRkZFsNruqqkrYEYFhqn8dyHp6ert37753715VVRW2KqO2tvaaNWva2tqWLl06YsSIZcuWTZs2TUJCora2tri4GPt7y9fXd/To0YWFhVevXm1ubs7IyHjy5MkgfU4ABhGREtKTJ08WL1786NGjioqKpqamqqqqnTt3HjhwQNhxgeFo9+7d0tLSCQkJTk5OFAply5Ytp0+fnjlzZk8dyFhzB3tuJCUlpaOjU1ZWhhDS1tZGCKWlpWlqamZkZOBjwVNTUy0sLGRkZBBCKioqsrKysbGxsbGxra2tOTk56enpg/Q5ARhEROqyCwgIGDt27LVr1/773/+SSKQNGzawWKzIyEhfX1+YGwsG2cePHxMSEmpqavrYgWxoaJiUlHT+/PlVq1bdvXu3oKDAw8MDe4tOp2Nju9vb20tKSrD+utTU1OnTp2MV+OfbOjo6rlu3btWqVQP30QAQFsK0kNhsdlFRkbu7O/Y3I0KIRCItW7assbHxzZs3wo0NDEObN2/esGGDuro6QkhBQWHhwoW9P85UV1f38/O7e/fu5MmTDx48uGHDBldXV+ytixcvzp07FyGkqamZmZmJEOJyuU+fPsUTEj8SidSXGbUAEBFhWkji4uJKSkpFRUX8hWlpaQghWPAfDDIul8tkMsXFv+zXZ9GiRYsWLaqsrFRTU5OQkMDLf/nlF+zl1q1bsWRDIpGuXr1qYGDQ9SSJiYlfFzsAooswCQkh5OTkdO3aNRKJVF1draSkdPLkyaioqLFjx2ITNQAYNGQyefny5WfPnl2zZo2pqSn/2r6rV6/GG/Hdwh4a8cOTE97zTCKRus1GAAxtREpIO3furKqqCg8PRwhhTSUKhXL8+PHeezCKiorS0tJ0dXXt7Oy6fdRUVlaWkpIyefJkWB8M9N2DBw8QQi9evHjx4gV/+aJFi3pPSACAnpA67eki+vLy8oqKilpbW/X09KZMmcLf9dFVaGjo0aNHp06dWlBQYG1tffLkyU4VAgMDL126NGXKlOzs7GnTpgUHB3c9iaGhYdddq8Ewx+Fw+Key4nq/IQcO3KVgCCBSCwkhlJGRERERgTePWCyWvb19T5WrqqpCQkICAwMXLFhQUFCwcOHCly9fmpub4xUaGxsjIyP9/PwWL1786tWr5cuXl5SUQFcJ6AmPx4uIiHBzcxMXFxcTE+tpbCeHw7l3756ZmRm+FgMAoC8IM8oOIXT//v3Vq1cnJCRIS0srKys/fvz4xx9/PHv2bE/18/LyyGSys7MzQsjIyMjQ0DAlJaVTHXwpZUlJSYTQp0+fBix8MBQ8ePBg7ty5V69e7fZWaWtri4iImDVr1md7kgEAXRGphXTixAkTE5M///xTWVkZIdTc3Lx3796TJ0+6ubl1u3pYeXm5uro6PhRKS0urpqaGv4KiouK2bdv27t177969nJwcJycn/vYTAJ1gg98SEhKOHDly4MABMzMzXV1dbW1tEolUUlJSUlJSVlZGJpN//PFHLy8vaWlpYccLAMEQpoXE4XDev3/v5uaGZSOEkLy8/MaNGzkczrt377o9pK6ujv9LQUZGBlt8DNfS0nL//n1lZWUdHR1lZWU6nV5dXd3tqbp9tgSGJycnp/v37+/du1dDQ+P169fnz58/e/bsu3fvjIyMtm/fHh8fv2HDBshGAPQDYVpIYmJiVCo1ISFh6dKleGFhYaGsrGxPw75VVFSam5vxl01NTYaGhvwVkpOT8/PzExMTdXR0sMdR0dHRW7du7XoqWLYS8BMXF1+5cuXKlSsRQmw2m8Viwcg6AL4eARJSTEwM1l+vp6cXHx/v7Oxsb2+voqKSn5//999/f/PNNz2tHjZy5EgGg9HY2KioqMjj8YqLi2fOnMlfoaSkRENDQ0dHByEkISFBpVLfv38/CJ8IDCXi4uJfOkMWANAtAvwiRUZG5ufn4y+Li4uLi4vxl0lJSc3NzVJSUl0PtLW1VVJS+vXXX7dt23blypWmpiZHR0eEUHR0tL6+vpWVlbW1dUhISGRk5OLFi1NSUtLT0y9cuDAInwgAAEBXBJiH9NmNNXtZWfXJkyfbt2+vr6/X0dHx9vaeN28eQsjOzm7u3Lk+Pj4IofDw8OPHjzOZTDKZvGLFin379nU9CczwAKIP7lIwBBAgIWVkZJibmwuxVwR+1YHog7u0q+Tk5MePH9vb2zs4OAg7FtAnot5lx+PxPDw8pKSkrK2tp0+fPn369FGjRgk7KACAqPP39/fz88N+9vDw2L9/P4VC+dKTQEobZARoIbW0tDx//jw1NTU1NbW8vJxCoWCZycrKanCGNsHfngDD4/HWr1/fx8o7duzoffdYwYK7lB9/NsJQKJTS0tL+ncfPz2///v0CCw70TNRbSAghOTm5mTNnYgPkaDRaSkpKamrq9evXORzO5MmTp0+fvmjRIv7llgEYICQSSU1N7bMPNTGwaaQQdcpGCCEajZacnAwNHRFHgITEj0KhUCiUNWvWtLe30+n06OjoGzduTJgwwcrKStihgWHht99+E3YI4DNoNFq35Y8fP+57QqLRaDNmzMBOFRYWhhCCRtIgIFJCamxsDAgImDVrlp2d3a5du+Li4rS0tI4cOWJpaSns0Pqjo6Ojtra2srISe6mqqqqqqtrTnCogaphMpoSEBH8ziMPhlJaWUigUmJY0BOjr6+M/02g0rMnVU05iMBgIIQ0NjUEJbSgj0m/O77//Hhsba2lpGR8ff//+fTs7OwaDsWnTptTUVDKZMGsg4d6+fZuTk8NisaSlpdlsNofDsbCwMDU1FXZcoE9cXFy8vb3nzJmDl2RmZrq7u8fFxfF/l4HBR6FQHBwckpOTO5WvXbu2j2fAmkRdC8eNG/f69WuEUHV1NZaEumVqaoptIjp69GjspYaGhqamJmSszyJSQkpJSdm8ebOrq+v27dtlZGRCQkKKi4tdXV1pNBrh9oyg0+m1tbUKCgr4yB8ajVZbW0un07GVI4Boio+Pj4mJQQjV1NScPXv2+vXrWDmPx8OemeNrLQIhunjxIt7hhpf0cZQdg8G4du1a13IajXbhwgVZWVm8pL29HZuS/+rVKxMTE2zHAIRQbm4udp7CwsKPHz+eOXMG6//A05KpqWmnVWMAhjAJicPhfPr0ady4cRwO58mTJ9bW1lJSUlgHF5vNFnZ0X6ajo6OysvLTp0/8uWfkyJE0Gu3jx49qamr4nQ1Eja2tbVRUVFtbG4vFYjKZ+GKJJBJJV1fXw8MDxteIAgqFkpSU5O/vHxYW5uDgsH///t6fHjEYjOrq6tzc3NevX+fm5vKvBcNv3bp1WEbp1Nzx9PS0t7f38PDAz4b+bUXduHHjzJkztra2WDn2VmJiYlBQkIaGxsyZM01NTaFfBEeYhCQmJjZy5Mhbt26Vl5fX19c7ODg0NTWdOXNGTU1t3Lhxwo7ui7W3t7e0tMjJyeElEhISzc3Nnz59gmwkyuTk5C5evIgQ8vHxcXV1hdE0IotCoaxdu5ZCofQyi4jBYDx8+BBLQvzllpaWhYWFnSp7eHi4ubl1ex7slsBhuUpDQ4NGo2HLno0dO/bixYsMBiM3Nzc3Nxf7gcFgXL169erVq3ibCTITAeYh4WJjY318fLhcLpVKjY6OPn78eGRk5KZNm37++ecBva7AZ3h0dHRkZWVVVFTo6OjgOamlpYVOp1MoFFNTU8hJoo/D4SQlJVVUVLBYLOwlk8nMzs4+duyYUBpJMA+p77A8lJiYyGAwNDQ08CEJ06dPt7S0xLJCWFiYp6cnfoiDg0NSUtIXXaXTGbBGG95tiCcnLDNhhVibadasWcP2aROREhJC6P3799XV1aamptLS0q9evRoxYsQgPEAeiF/1jx8/5uTkYJ2QWAmNRlNRUTEwMIBnSIRw9OjR0NDQToXa2tp37twRylBJSEifheWAK1euoH971ZSUlCwtLTU0NLpt+tBoNKzTrx8TY5OTk2fMmNGpsKeshiXIq1evYi+Hc4OJAF12UVFRbDZ79erVly9fbmxsRAhlZGTwV/D09CTcbjQKCgpqamp0Oj0nJ0dcXBwbQ6yvrw/ZiChu375tY2Pj7+/v5ua2devWsWPH7t69e8qUKTBwX9Twt4ewEg0NjZUrV3724Q0267F/F338+HHXwq4D//B43NzcZs2ahTWYEhMTExMTsV5ENze3YTX8gQAtpO+++665ufnhw4ezZ8/udspbWlqaqqrqwAUwQH97NjU1NTU18c9D0tHRgc46QuBwOOPHjw8MDJw/f/5PP/2kp6fn6+ubl5e3ePHiFy9eKCoqDn5I0ELqqlPLo7293cPDQ1NTs+9f8f1ey67r2kUYbKZa78fyh/3x40cxMbGNGzcuW7bsiwIgKAK0kGJjY7Ef4uPju02fJBJpcCMSDAUFBQUFBWgSEZGYmJiSklJLSwtCyMDAAGuy6+vrc7nckpKSSZMmCTvA4e7KlSudmkQzZ87saUhCLxwcHPq32pCenl635X1pcuENpuDg4KNHj7a1tWVmZgYEBJw+fXrID6IhQELiR9DcA4Yec3Pz06dPa2homJiYhIeHP3jwoLCwUFxcfDAXVAWddH0YMzhjBDoN+0b/ZrJOfXRfNCwiLy8PX6cKy0kODg6///67m5vbEB7yQICEdOfOnbq6ul4qLFu2TFpaetDiAQAhtHv3bi8vr4SEBH9//zNnzmzZsgUhNHv2bHiGJBSdUpGSktLcuXP70SQSFAqFcvHixfDwcKzjDnv5RY0tf3//TiVtbW1nzpzJzc3tX2uPEAjwDGnBggVdO8fJZDKXy8V+Tk9PHzFixMAFAL3zoFs8Hq+2tlZdXb2pqenRo0eamprW1tbCCmZI3qXYIxw9PT3+xkcnXVtFIjUQoN8bWHTbG6SqqmpjY4MQMjU13bJly9BrKhGghXTt2jVswf9Tp05FRUUdOXJk2rRpLBbr2bNnZ8+elZeXH9BsBEBP8vPz//e//717987Z2VldXb2qqkrYEQ0dNBrN09MT7/J6/Phx1x32OqUiU1NTbOzc4Eb6Gfb29n5+fvb29l96IIVC6TqGa/78+YcPH/bx8cnNzd2+fbtwW4EDgQAtJNy0adOWL1++efNmvCQ9Pf2zy1kWFRWlpaXp6ura2dl1u0XNmzdvMjMz9fT0vvnmm27XaR6Sf3uCr/TkyZMffvhBSkpKTExsxYoVsrKyJ0+eXL169d69ewV4lZSUFAaDsWTJEoQQjUZLTU2dOHGimZlZ15pD7C6dMWNGpwcw/JN4OqUirAtriDUXPD09u67xmpSU5ODgwP/xh1hTiTCLZHO5XBaL1Sl91tfX935UaGioi4tLYmLizp07t27d2rXCpUuXli9fnpycvGvXLmdn5z7uvQZAQEDA2LFjnz17NmnSJBKJtGHDhvXr19++fVuwt1BUVFRHRwdCyNvb28XF5dGjR99///28efMIt37jFwkLC+s6ZSc5OTk5OZnBYFy5cmXdunX41/GFCxeG0jcyrmuLEB/yh3VLXrhwQUNDIzc318fHB5vtOxTwiGPz5s0WFhY3b96sqKh4+/bt//73PwcHh9mzZ/dUv7Ky0tTU9NatWzweLz8/n0qlZmVl8VdoaGiYNGnSn3/+yePxGhsbzc3NscqdUKlUQX8UQGwsFsvIyCg6OprH43l5eQUGBvJ4vIqKCiqVmp2dLcCrmJubl5WV0Wg0IyOjzMxMHo9XX18/ceLEu3fvdqo8lO7STkvD4fbs2ePyLx8fn+rqamFHOrBKS0svXrzo4eHh5+eXlJTUtUJ1dXVkZCT2DxIZGTnoAQoeAZ4h4fbt27d58+adO3fiJQYGBufOneupfl5eHplMdnZ2RggZGRkZGhqmpKSYm5vjFbD/xx4eHvX19ZKSksnJyfxrywPQE3FxcSUlpaKiIv7CtLQ0hJAAn2hmZWWpqanp6uq+ffvW29vbwsICIaSoqKiiovL69et58+YJ6kKipqfJOmlpabKyshoaGlu2bBG1Z0UDgUKheHh49DKgA5+x5OPjg63vcOHChUEMUPCIlJBUVFTCwsLy8/MLCgqYTOb48eNNTEyw/Ui6VV5erq6ujj8W0tLSqqmp4a9QWVmpoqKybt067Ktk+fLlsEsx6CMnJ6dr166RSCRsK7aTJ09GRUWNHTtWV1f3S0+Vmpp6//59fNQoQsjDw8PY2DglJQXbtsDQ0NDQ0BB7Ky4ujk6n29nZdT1PcHDwxo0b+/uBREi3O+ypqqpSKBTRTEVd5yENJg0NjT/++AObC7xu3Tr+QYZhYWHh4eE0Go1CoXx2Dw5RQKSEhBASExObMGHChAkT+lK5rq6Of36SjIwMNrUe9/Hjx4qKiilTpsTFxT179szPz8/U1BR7gNzJkPlVB4Kyc+fOqqqq8PBwhBDWVKJQKMePH//SudtZWVmlpaXOzs7Z2dmKiooSEhJ6enpjx45FCD158oR/CA9CKCUlZdeuXWvWrMHG/nYyZG5RbNYO/w57FArl4MGDQ2xEmQBhTSUNDY2rV6+GhYVVV1e7ubnxL19Eo9GSk5OxDkBhBvpZwu4zHECXLl1ycHDAX3p6eh46dIi/QkhIiJGR0cePH7GXrq6u27Zt63qeodQ7DwTrzZs3sbGxV65cefr0aUdHRz/OQKfTeTxefHx8VVVVbGwsm83GymtqaiZMmNDa2orXjImJMTY2xp5XdTXE7tLIyEhHR0czMzMqlXr48GFhh0MYDx8+xB4p7dmzp+u3PYVCEXaAn0GwFtIXGTlyJIPBaGxsVFRU5PF4xcXFnabL6ejokEgkfKVwHo8H0+zBFzExMTExMfmaM2hrayOE0tLSnJycMjIy5s+fj5WnpqZaWFjgN+fFixcDAgJ27dq1du3ar4xZxOXm5gYFBTEYDFlZ2alTp/7+++/CjohIsE0rfHx87ty50/VdrJ0kyh13hBn23Q+2trZKSkq//vornU4/evRoU1OTo6MjQig6Ojo9PR0h5OTkJCcnd+jQoebm5tu3b+fl5XXbNQ9AtzIyMjZt2jR37ty5c+euX7++2x0H+oJOp2Nju9vb20tKSrDC1NTU6dOnYz9nZ2cfPnx4wYIFxsbG6enp6enp3S57T3QMBsPX19fX1xfbN+/333+HbNQP2CMlTU3Nbt/t94Yag4NIE2MRQhkZGREREXiX/YoVK3qfAv3kyZPt27fX19fr6Oh4e3tjA5Ps7Ozmzp3r4+ODEHrx4sWOHTuqq6tlZWW///77//znP11PMsSmHAKBuH///tatW0kkkrGxsYyMzKtXrzgcztatW/uxf/HBgwft7OymT59+5MiR0aNHL1++nMvl2tjYXLp0CRvL4Ovre+PGDf5DFi9e3OnLmuh36ZUrV7CpRf1emRvwo9Fo3S4XIOJf+ERKSAL8CvgiRP9VBwNh9uzZsrKyf/75p7KyMkKoubl579698fHxz58//9L9kFgsloSEBEKIw+GQSCQymczj8UpLSw0MDPp+EuLepQwGIygoCNuPbubMmdgyteDrdd2TCVvoQTjR9JEwH2B9IScnp+++++7Tp0/Yy6ampi1bthgbGzc0NAzodYfY42Lw9dhstpGR0bVr1/gLi4uLqVQqNn118BH0LsXndXp5eeXk5Ag7nKEGW2IcG0bf7dRaUUOYZ0gcDuf9+/dubm7YH6QIIXl5+Y0bN3I4nHfv3gk3NjDciImJUanUhIQE/sLCwkJZWdl+zEMannJzc+fPn491061cufLChQsiOMGoj7pdd04UeHh4JCUllZaWEqBthBAi0Dwk/Ctg6dKleCF8BYBBFhMT8+nTJ4SQnp5efHy8s7Ozvb29iopKfn7+33///c0338BAzc/C1qNLTExEQ25tUPCVCJCQ4CsAiI7IyMj8/Hz8ZXFxcXFxMf4yKSmpubm5l9VDQGJiYlBQEBpagxd6WnwPfCkCDGr47rvv+L8CukpLS1NVVR24AIj7uBgI3GcX8+52i5NBIPp3Kf/gBVNTUxjSDboiQAup04DXroT1FQCGoa43W1hY2Lt373777TehxEMU/A0j0VyPDogCAiQkyDdAlOXl5WVlZQk7CtEFo7pB3xEgIeF4PN7Tp0+LioqYTCZ/uYeHB77CCgBAdEDDCHwRIiWkU6dOBQcHdy1fuXIlJCQARAo0jEA/ECkhXb161crK6tChQ1paWvyL/JPJhJlNBYaeb7/91sjISNhRiJbc3FxfX180bBpGwt0PaSghTELicDgNDQ3ffffdyJEjhR0LACgsLKyhoQF/eeLECYSQhISEsrLytGnT9PT0hBeaMEHDCHwNwiQkMTExU1PTlJSUxYsXCzsWAFB8fDw2loFMJpNIJP7h4JKSksHBwYSYGC9Yw61hhIN5SIJCmISEEHJycjp06NC6devGjx/P/9AIBjWAwTd//vySkhJfX187OzsxMbG0tLTDhw/b2Nj89NNP/v7+u3fvfvr0qbBjHFT8iy/AHCPQPwSYGItzcnIqKyvrWg4TY8Eg43K5U6dO9fT03LBhA14YERFx5MiRrKysgoKCRYsWJSYmjho1atBCEuJdymAwfHx8GAwGQmjlypVDY/EFIBRE1ek9mwAAIABJREFUaiHFxcV1Ww6DGsDgI5FI7e3t/CWtra1MJpPJZEpLSyOE2Gy2kEIbVPg+RrAqHfh6op6QeDxeRESEm5ubuLh4L4mHw+Hcu3fPzMxs2D5MBoOJTCbb2tpevnxZXl7e3t5eTU0tMzPz4sWLCxYsYLPZJ06cGDlypIhvzfn1+McvdG0Y0Wg0bFvbYfgsDfSbqCckhNCDBw8iIiK8vLxmz56toqLS6d22trbr16//+eefCKHw8HBhBAiGoz179tTU1Bw5cuTIkSNYib6+vo+PT0pKSnx8PLYf8RDW+4zXsLAwT09P7GcKheLh4bF//34hRAmIhhjPkBISEo4cOVJeXm5mZqarq6utrU0ikUpKSkpKSsrKyshk8o8//ujl5YV1lQgcPEMC3eJwODk5OW/evGEymWPGjMFGNzAYDA6Ho62tPcjBDNpd2nvDCCGUnJw8Y8YM/hIKhYLtFDcI4QkFzEMSFAK0kBBCTk5Ojo6OMTExz58/f/369d27d7lcroGBgYmJyZIlS+bMmaOpqdnTsUVFRWlpabq6utj3RU/V/vrrL2NjY2Nj44H5BGCo4fF4z549wxeyKigoKCgoQEN9zGdfBnb7+/t3KqHRaP7+/kM4IQFBIUYLqRM2m81isfryax8aGnr06NGpU6cWFBRYW1ufPHmy22oPHjzYsmXLvn37Vq1a1fVdaCGBrkJCQrpdyOrFixdKSkqDH88g3KX84xd6Gditr6+PPT3iR6FQSktLBzQ8MAQQo4XUibi4uLj45yOvqqoKCQkJDAxcsGBBQUHBwoULX758aW5u3qladXU1dHCDLzWsFrLCB3YrKSl5eHjMnDmzl8oUCqVrQoLmEegLQiakPsrLyyOTyc7OzgghIyMjQ0PDlJSUTgmJx+Pt2rVr1qxZz58/H+TwmpubS0pKcnJysJcGBgZjxozppe8RiI5htZDVlw7sXrt2bXJyctfCAQoPDCVDOSGVl5erq6vjbSktLa2amppOdcLDw9+/fx8SEjJ//vxBDu/vv/9OSkpiMpmysrIsFuvp06dOTk6urq6DHAboh2GykNVnxy90y8PDo6yszM/PDy8Z2iMagAARKSHdunXL2tq6722Iuro6/nF3MjIyLS0t/BWKioqCgoIuXLggJyfX+6mCg4M3btz4pQH3Iicnp6ysTFlZeeLEifwlOTk5eAkQZUN+IauvWZhu//79eDvJwcFhyE/JAoJCpIQUEhKya9euadOmubq6fvvtt1JSUr3XV1FRaW5uxl82NTUZGhryV/D19R09enRhYWFhYWFzc3NGRoaenp6trW3XUwk2GzU3N+fk5NDpdP5tC6hUanZ2dnFxsYGBgby8vAAvBwYC1ov15MmTJ0+e8JcvW7aM6AmJwWDgC9P1e8VubPqRgCMDQx2REtKVK1fu3LkTGxu7fft2BQUFZ2dnV1fXroMUcCNHjmQwGI2NjYqKijwer7i4uNPDWBUVlfr6+tjYWIRQa2trTk6Orq5utwlJgGg02o4dOz58+NDe3p6dnc3/1ocPH+Tl5aHXjhCG6kJWubm5QUFBDAZjuK3Y/TVgHpKgECkhqaure3l5eXl5FRQU3Lp16++//46OjqZQKIsXL160aJGamlqn+ra2tkpKSr/++uu2bduuXLnS1NTk6OiIEIqOjtbX17eysjp37hxe2dHRcd26dd0O+xa4efPmpaam0ul0CQmJhISERYsWSUhItLe3S0tLGxgYNDc3QwtJ9BE98XSrjwO7ARggREpIOG1tbT09PT09vQ8fPtDp9FOnTp08eXL37t0rV67kryYpKXn48OHt27ffuXNHR0fnwIED2Pz5U6dOzZ0718rKir8yiUTiH7w7cLCuDC0traioqI8fP4qJiWloaHA4HDabbW9vP3HiRMhGoiwqKorNZq9evfry5cuNjY1dK3h6ehKxy65/4xcABvZDEhQiJaSOjo6kpKTbt28/fvyYxWKZmpr6+fnNmzcPIXTgwAF/f/9vv/22UzvJ1tb2xYsXnc6TkpLS9eRYj/ngoNPpCKFRo0ZxuVyEUF1dXU1NjZycnKSkJIxoEHFRUVHNzc2rV6+OiIjoOtsGIbR8+XLCJSS8YaShofHHH3/Ait1AWIiUkBYtWlRUVKSqquru7r5o0aJx48bhby1cuPD27du1tbVdO+5ETUdHR2VlZUtLi42Njbm5eVRUlImJCbbfqKmpaUdHh6SkpLBjBD3CnjgihOLj47td5USw7eyUlBQGg7FkyRK8RLBrXEHDCIgUIiUkJyenLVu2zJgxA1+SjsPhlJaWUiiUSZMm/f3336NHjxZuhH3U3t7e0tKCJ9S5c+cihF6+fNnW1gbZiEBevXp15cqV/Px8Doczfvz4+fPn29vbC/YSUVFR/KNsHjx44OPjs2/fPoEkpGG74zgQWUR6MHv79m02m82/QGpmZua8efOwwWm6urqD8xDo65HJZDk5Of5JUS0tLfLy8vLy8h0dHUIMDPTdnTt3Vq1ade/ePTk5uXHjxmVnZ//888+XLl0S4CXYbPbz58/xhCTANa6wgd1YNpo5c+aFCxcgGwFRQIAWUnx8fExMDEKopqbm7Nmz169fx8p5PB62XKOysrIw4/tCkpKS48aNa2tro9PpeCOppqZGU1NTW1sbWkhEERERMXHixNDQUEVFRYQQl8vdtWtXUFDQypUrJSQkBHKJrKwsNTU1XV1d9OVrXGELbGNTUzttRwQDu4HIIkBCsrW1jYqKamtrY7FYTCYTn+tKIpF0dXU9PDyEsrjy11BQUFBTU6PT6Tk5OTweLysrS0JCQl9fX0dHR9ihgT7hcDivX78ODAzEshFCiEwm//DDD7du3crLyzMzM/uis6Wmpt6/fx8b4YLx8PAwNjZOSUnBm0dftMZVpx2J/Pz8wsLCsL/eYGD3QIB5SIJCgIQkJyeHjar08fFxdXXtNFybiCQlJSkUiqqqamVlJZvNNjU1VVVVhWxEIGJiYmPHjr179y42yBOTmpoqJiamr6//RafKysoqLS11dnbOzs5WVFSUkJDQ09MbO3YsQujJkyebN29GfV7jCl/gqtsdiQICAt6+fctgMBCMXwCiipD7IQ2yAd1phkSC/wWEdOLEidOnT1MoFDs7O2Vl5WfPnmVlZenp6WEtGBKJ5OXl1ZcO2MrKSm1t7YSEBDMzs+fPn7u4uGBPSWtra2fMmJGeni4jI7N06VImk4mlkOPHj0+bNm3JkiWdlhThv0u7fZiqqqpqY2MDA7uBKCNAC+nmzZt1dXW9VFi5ciWxZn50dHTU1tZWVlZiL2k0mqqqqoKCgnCjAl/k3r17CCEajcY/G6m4uDgoKAj72c3NrS8JCZusnZaW5uTklJGRgffIpaamWlhYYDe2QNa4kpGRgYYREHEESEhhYWHY5tA9WbhwIbES0tu3b3NyclgslrS0tLi4+OPHjy0sLODZMrH0tJYdru9rC9HpdGx0ZXt7e0lJCdZfl5qaOn36dKzCl65x5eDg0HVHogMHDkA2AiKOAAkpJiYGmzfaE/49JkQfnU6vra1VUFDA1uQnkUjKysq1tbV0Oh0eIxEImUx++fKlQOYhXbx4EZuLpqmpmZmZOXbsWC6X+/Tp059++qlr5b6scXXx4sVOj7IcHBzgkTsQfUR6gHHjxg3+7SRwy5cvH9CcJMBnSB0dHbm5uTQaTUdHB3tAPXny5GfPntFoNCqVamhoCMO+ieLOnTs7d+5ECJmammpqaubn55eXl/v4+KxZs+ZLT8VisbCR4hwOh0QikclkbEqDgYFB30/S6S6l0Wjh4eH4jkSCmsAEwIAiUkKaNWvWhw8fOhWOGDHi4cOH+OjbgSDYhPTPP/+8e/du0qRJWMnkyZP/+eefly9fUigUgc/zBwNn2bJlZDK50zykhw8fvnjxQlDzkL7IgA69AWBwEGmlhjt37mT9Ky0t7cqVK0ZGRg4ODgOajQQOlmkYArB5SO7u7p3mIbW0tOTl5Qk3NjD4PD09w8LChB3FUECkhCQjIyP3L1VVVUtLy19//fXWrVu1tbXCDq2vsGUaVFRUsAW/MbBMA+Hg85D4C/s3DwkAgCPAoIbPKi8vF/FFvvknKra3t3/8+LGmpobJZLa0tHC5XD8/Py0tLRMTk64jv6HrX2TNnDnz9OnTs2fP5p+HRKFQIiMj0ZfMQwJDAOyHJChEeoYUExPT1NSE/czlcpuamv7666+mpqZnz54N6LDvr+ydx5YUc3Bw4C989+7dmzdvsOW9JSQkSCSSgoLCuHHj+AfaYUdBThJNTk5OZWVlvVTIyMgYzM5keIYEhgAitZDOnz/f6StAQ0Pj0KFDoj8JqVNeodPpz58/NzAwYDKZTCYTIVRfXy8jIzNt2rRZs2bB4G9CEOA8JAAAhkgJ6e7du/wTkiQlJfm3oiCKjo6O9+/fFxcXl5SUYAuLSUlJsVisxsZGLpc7YcIENTU16OoRfZBvABA4IiUkSUlJBoNRUVFhbm6OEHrw4IG9vb2srGzvRxUVFaWlpenq6trZ2XWbwGg0Wmpq6sSJE790keZ+6+joePv27cePH1VUVLDpsQihly9flpaW5ubmWlhYDE4Y4GtwOJykpKSKigoWi4W9ZDKZ2dnZx44dI9zy8wCICCIlpMzMzB9//HHixInYI8TDhw/7+/v7+/vPnj27p0NCQ0OPHj06derUgoICa2vrkydPdqrg7e0dFxc3ZcqUkydPamho3Lp1S1x8wP9NWltbGQxGZWUl9i0mKSkpKyurpKRUU1NTVVX16dMnFRWVgY4BfKWgoKDQ0NBOhdra2kRstQMgIojU7RAeHq6jo7Nt2zb8pbm5+YkTJ3qqX1VVFRISEhgYeOnSpUuXLsXHx798+ZK/QllZ2b179y5dunTx4sWHDx+Wl5fHx8cP7GdACBvF0NLSUlNTw2AwGAxGWVlZYWFhe3s71gkJ2YgQbt++bWNjk5CQoKamdvDgwejo6LFjxzo4OMAiucMQzEMSFMIkJA6Hk5qa+sMPP+CLkOrp6e3YsaO4uLiwsLDbQ/Ly8shksrOzM0LIyMjI0NAwJSWFvwKTyfT29sa6yBQVFVVUVF6/fj3AnwPR6fSKigo2my0lJdXe3t7Q0NDQ0IAt/q2goKClpdXt8khApHA4nOrq6sWLF+vp6U2YMKGwsHDSpEmBgYHR0dGNjY3Cjg4AoiJMl52YmJiSklJmZuaCBQvwQmxWfE+Pl8vLy9XV1fEuOC0trZqaGv4KhoaGhoaG2M9xcXF0Ot3Ozm4ggudfeplOp+fn59fU1LS1tbW1tYmLi/N4PDKZjD2NSEhIKC4uxo/qNFgciAjsbsSW2zAwMMjIyEAI6evrc7nckpISfF0oMEzAPCRBIUxCQgjZ2dndvHlTV1fX2tpaUlKyoKDg2LFjFAplzJgx3davq6vjX3RVRkaGf8EefikpKbt27VqzZo2NjU23FfC9OL8eg8Goq6uTkpJSU1NrbGwUFxdns9lMJpNEIklJSY0YMUIgVwEDzdzc/PTp0xoaGiYmJuHh4Q8ePCgsLBQXF+/pbgQAfBaREtLu3bvr6uoCAgLwkjFjxhw/fryn1fhVVFT4u7+amprw9hC/69ev79u3z8vLy9vbu6dLf2U2wuchNTU1Xbt27dGjRxwOh8fjSUlJ1dfXs1ismpoadXX1KVOmfPfddxMmTIBh36Jv9+7dXl5eCQkJ/v7+Z86c2bJlC0Jo9uzZ8AwJgH4jUkKSlJQ8efIk1uXV2Nioq6s7adKkXgY1jRw5ksFgNDY2Kioq8ni84uLimTNndqpz8eLFgICAXbt2rV27doDDRwghKSmpUaNGYQPqOjo6DAwMRo4cWVdXRyKRjIyMVFRUmpqaIBsRgpqa2v3797HGbnR09KNHj2bOnFlVVcVmswdhoCYAQxIBfnOioqLYbPbq1asvX77M/8S4srLyxYsXCCFPT89uF2uwtbVVUlL69ddft23bduXKlaamJkdHR4RQdHS0vr6+lZVVdnb24cOHFy5caGxsnJ6ejhDS0NDAJwYNEDk5uZEjR5aVldFotDdv3rS2tiKEZGRkFBQUDAwMsDW/ISeJPhcXF29v7zlz5iCEFBQUFi5cmJ6e7u7uHhcXB+urAtA/xEhIzc3Nq1evjoiIoNFoXSssX76824QkKSl5+PDh7du337lzR0dH58CBA9ra2gihU6dOzZ0718rKKjo6msfjxcbGxsbGYocsXrz4999/H7jPIikpaWxsnJubKyYm1tzc3NTUhPU3crncnJwcHR2dRYsWQTYSZfHx8TExMQihmpqas2fPXr9+HSvHttRDCCkrKwszPgCIjAAJSUFBYdOmTQghPz8/GxubrqvB9rKjs62tLdaK4ocP/v79998HNP10q7m5uaOjo7W1VVpaWkJCQlFRkcVikUgkDodTXl5eVFQEa9mJMltb26ioqLa2NhaLxWQy8YeUJBJJV1fXw8MDlmkYhjw9Pe3t7WGT+K9HgIT08uVLKysrR0fHvXv34p0kxEKj0fCR34WFhbm5udXV1e3t7QoKCjweT1ZWlsPhkMnkDx8+3Lp1C8+4NBptoPsPwZeSk5PDxvj6+Pi4urpaWVkJOyIAhg4CbD+xYsWKV69ejRkz5v3791paWl3/Ar1w4YIob2FOo9E8PT3xl83NzVVVVdgAdEVFxaamJnl5eTabjRBqa2uTlZUdN24cXvnixYuQk0BfwPYTYAggQAvp2LFjISEhdXV1ZWVlCgoKXffi66XLThRQKJSkpCTs546Ojri4uFu3bmGLGI0cOfLZs2cmJiZtbW1kMplEIk2fPn3Lli3Q7QMAGIYIkJB27tzp6enp6Ojo5ua2devWKVOmCDuiryIvL6+pqTlmzJiysrIPHz5wuVxsdANCyNLScsaMGZCNAADDEwHWsnv58uWbN28QQgwG4+PHj8IO52uJi4tLSUlpamoymcyqqqr29vaioqKKioq2tjYNDQ17e3thBwgAAMJBgBbShAkTTp06FRcXV11dffTo0QsXLnSqMNDPkARLUlJyzJgxWlpaLBaLRqMlJSWNGzdOXFycQqFMmTKFTqfDKDsCaW5u3rdvH4PBCAgIgP9xAHwlAiQkoj9D4icpKSknJ5eQkFBbWysmJiYvLy8uLq6hocHhcFpbW2/cuCEvLw87xhJFW1vb/PnzAwMD8/Lyli9f7uLisnPnTmEHBQCBESAh6ejoYLOFvLy8PDw8BmhB7kGjoKAwYsSIDx8+MJlMZWXl6dOns1isESNGjBgx4v37921tbZCNRFlubm5JScm0adPU1dX9/Pyw4SqTJ09es2ZNdHT0nTt35s+fL+wYwWCDeUiCQoCEhPvzzz+FHYIAUCgUf3//Z8+e5eXlqaury8jIjBgxQl1dvaOjg06njxo1CpYOEmWmpqampqbV1dU3btyIjY01NTX99ttvNTU1EULLly/fs2cPJCQA+o0ACamntexwPa1lJ7JUVFSsra1lZWU/ffqEzzqqqKjQ1NTU1taGbCT6NDU1Fy9e/M8//6xevZpOp587d+758+dOTk5Lly7FFvMVdoBgUMF+SIJCgImx3333XXNz88OHD2fPnt3tWnZpaWmqqqoDF8BATDns6Oh4+/ZtTk4Oj8cTFxdnMpkSEhITJ07E98MFoq+ysjI3N9fJyQkvycrK+vvvv9+9e+fq6optVTxoYGIsGAII0ELCVz6Nj4/vNn0SaFADTlJSkkKhqKqqVlZWYiWqqqowTItYtLW1g4KCWlpaXF1dsZK7d+/u27dPuFEBQFwESEj82traPnz4oKKioqamRsQ8xE9BQUFBQQGSEKEdPnx4z549+/btU1NTa2pqOnTo0ABdqKioKC0tTVdX187Orpc9wAAgNAJ02SGEOBzOuXPnYmJiKisrsYBlZGRcXFzWr18/cuTIgb46dIaA3rW2ttbU1IwaNWqAUkVoaOjRo0enTp1aUFBgbW198uTJrnXgLgVDAAFaSCwWy93d/eXLl1paWk5OToaGhnV1dQUFBTExMY8ePYqMjIT90IBwycrK6unpDdDJq6qqQkJCAgMDFyxYUFBQsHDhwpcvX5qbmw/Q5QAQIgIkpJs3b758+dLNzW3//v385U+ePFm/fn1QUNCJEyeEFRsYVng83vr16/tYeceOHWPGjPn6i+bl5ZHJZGyIhJGRkaGhYUpKSu8JKTg4eOPGjV9/aZEi3A/V+9UHdB7SV35wYt0MBFjLLjExkUKh7N27t1O5ra3tokWL8N32ABhoJBJJTU1NuW8E1X1XXl6urq4uLv5/fztqaWnV1NT0fkhISMiXXiU4OFiANT9bp5cKPb3Vjw8lQH25et//Dft36f6dv1PkAxSkoBDgGdLSpUvV1NTOnDnT9a2IiIjffvstKytLTk5u4AKA3nkgRCdOnHj48OGdO3ewl5s3byaTycePH+9Uzd3dPT09fdCjAwRjZWUVEREh7Ch6RIAuO2yJnW7fwnZqYDKZX5OQYPwSEGUqKir4RukIoaamJkNDw67VRPlbBoA+IkCX3YAKDQ11cXFJTEzcuXPn1q1bhR0OAJ2NHDmSwWBga5TweLzi4uJRo0YJOygABgQBWkgIocePH7u7u3ctr62t/ZrTwvglIPpsbW2VlJR+/fXXbdu2XblypampydHRUdhBATAgCJCQKBQKm83uKfcYGBjgz3u/VD/GLwEwyCQlJQ8fPrx9+/Y7d+7o6OgcOHBAW1tb2EEBMCAIkJAGblhIP8YvCRyNRqNQKIN8UUAstra2L168EHYUAAw4AiSkgVNXVyctLY2/lJGRaWlp6VrNysqq28fIhPDf//6XQLMQAADD2bBOSDB+CRBdX8aI/vXXX8bGxsbGxoMcm2B9+vQpNjbWy8tLWNei0Wj5+fn8JRoaGpaWloK9NI1GS01NnThxopmZWae3uFzuo0eP3r9/P3nyZFNT024X88zLy0tPT6dSqTY2Nl0rlJWVpaSkTJ48WXRvBt4wlpiYaGJi0tDQwOPxuFyunZ3d5cuXhR0UAH117tw5KpXq7u4+ZcqUjRs3dlvn/v37VCqV6Dd2a2urm5vb1KlThXityMhIm39ZW1tTqdQNGzYI9tLbt28fP368h4fH5MmTnZ2dWSwW/haXy123bp2hoeHKlStNTU19fX27Hv7LL79MmjTJy8vLzMxs+/btnd4NCAiYMGGCp6enhYXFf//7X8FGLijDOiG1t7dPmzZt+/btFRUVgYGB5ubmdDpd2EEB0CeVlZWmpqa3bt3i8Xj5+flUKjUrK6tTnaqqqilTphA9IWVlZc2aNYtKpQ5CQurjta5duzZ+/Pi8vDwBXppGoxkZGWVmZvJ4vPr6+okTJ969exd/9+HDh1Qq9eHDh1hNQ0PDxMRE/sMLCgoMDQ3/+ecfHo+XmJhIpVI/fPiAv9vQ0GBmZnb9+nUej/fy5UsqlVpcXCzA4AVlWM9DwsYvpaamzpgx4969ezB+CRBIt2NE+SvweLxdu3bNmjVrEFbEH1CRkZEWFhY+Pj4icq3KyspDhw5t3LhRsB1fTCbT29vbwsICIaSoqKiiovL69Wv83YKCAiUlJQcHB4SQnp6eoaEh/7sIIR0dnUuXLmFdiK2trQghDofDX6GjowP7AduT+tOnTwIMXlCG9TMkBOOXAGF9doxoeHj4+/fvQ0JC5s+fL4wABebgwYNSUlI3b94UkWsdPXpUXV39+++/F+ylDQ0N8WfYcXFxdDrdzs4Of3f8+PFNTU11dXVqamrt7e0fPnzotH22goKClZVVbW3toUOHkpKSVq5cyb8CvaKi4rZt2/bu3Xvv3r2cnBwnJyfRnN8yrFtIABBX72NEi4qKgoKCAgICBnSZx8EhJSUlOteqr6+Pj49ftWrVwC0zlpKSsmvXrjVr1tjY2OCF06dPV1JS+umnny5duvTzzz+3tbUxmcyux7a0tGBr+2ZmZr5//56//P79+8rKyjo6OsrKynQ6vbq6eoDi/xqQkAAgpK5jRLW0tPCXvr6+o0ePLiwsvHr1anNzc0ZGxpMnT4QR5lBz8+ZNcXFxfNN6gbt+/frPP//s7u6+e/du/nIxMbGzZ8+qqKjExsaampp+++23urq6XQ/X09PbvXv3vXv3qqqqrl+/jpcnJyfn5+fHxMT89ttv9+/fr6ysjI6OHqCP8DWGe5cdAASFr3GnqKjI4/GKi4tnzpyJv6uiolJfXx8bG4sQam1tzcnJ0dXVtbW1FV68Q0R8fPycOXPk5eUH4uQXL14MCAjYtWvX2rVrO73FYDDKyspCQ0NJJBKPx5s5c6a1tTV/hSdPnsTExBw7dkxMTExKSkpHR6esrAx/t6SkRENDQ0dHByEkISFBpVL520+iA1pIABASvsYdnU4/evQovsZddHR0enr6uXPnov+loaGxbt26bdu2CTtkosL+SbGfi4uLx48fPxBXyc7OPnz48IIFC4yNjdPT09PT07GnRNjVxcTE9u7de+3aNSaTefr06fb29jlz5vAfbmhomJSUdP78+aampqtXrxYUFGAjILDDra2tq6qqIiMjmUxmQkJCenr6smXLBuJTfCVISAAQUk9jRE+dOpWYmMhfk0QidTuJknAG81PwXwv/J62trW1sbKRSqQNxxejoaB6PFxsb6/6v0NBQ/Oqqqqq//PJLSEjIpEmT4uLijh07pqKiwn+4urq6n5/f3bt3J0+efPDgwQ0bNmD9itjhU6ZM8fX1DQwMnDRp0pYtW1asWMH/gEp0EGCDPgBAT7hcbk1NjaamprADAYOBy+XW1dWpqqr2UqeyslJNTU1CQqLrWxwOp7q6WlVVdTDHiXwRSEgAAABEAnTZAQAAEAmQkAAAAIgESEgAAABEAiQkAAAAIgESEgAAAJEACQkAIEwcDodGozU2Ngo7EEGqqqr6+PFj/45lsVhFRUWCjYcoICEBADozMTEx/JeFhcXSpUtjYmK+5oR5eXlXr17tWv4+0pZAAAAgAElEQVTkyRNLS8vZs2d7e3t/zflFCpPJXLp06bt37xBCM2bM+OmnnzpV6LYQx+Pxfv755+Tk5AENUjTBWnYAgM54PJ6Zmdm8efM4HE59ff3jx4/37NkjIyPj4uLSvxNu2bLFwsJi5cqVncojIiLExcXPnDmjr6//1VGLisjISGyLP/TvltydKvB4PC6X29PhkpKSmzZtCgoKsre3HxpLbPQdtJAAAN2gUqlr16718vLatm3blStXpKSkOq1IhBDqaVp936fb19fXT5gwwdHRsZeE9PVX6cshgloioLm5OTQ01N3d/WtO4uLiUlVV9eDBA4GERCCQkAAAnyErKysnJ4f/Ud/Y2Ojr6/vNN9+Ym5uvXr36zZs3WHlpaen69eu/+eabCRMmzJ0799KlS1j5xo0bKyoqHj582GlBz+XLl+fn52dnZ7u4uFRWVrq4uJw5c8bGxmbevHnYU6Vur4IQqq6u3rJli5WVlaOj45UrV+bPn5+SksJgMFxcXJKSkrA6VVVVLi4u2C66PZ2qvr7excUlJiZm5cqVEydO/Pbbb/FNGVpbW/ft2+fo6GhjY+Pr61tcXLxnzx7+RbhLSkrw8/O7dOkSi8X65ptv+vIPe/XqVZcu2Gy2mJiYo6PjiRMn+nKSoQS67AAA3aDT6Y8ePUIINTY2Pnz4sLW11c3NDSHE4/F++OGHd+/erVixQlNTMyIiws3N7cGDB2pqaqtWrZKXl9+0aROZTL527drBgwctLS3Hjx+/bNmyzMxMCoXSaVcFLy+vI0eOyMnJff/995qamtimgi4uLm1tbRQKZfny5V2voqOjw2Qyf/755+rq6p9//rmhoSEgIKCtra2hoQEbC4APjujo6CgqKmpoaOgpYB0dHeyQPXv2LFiwwNXVNSoqat++fRYWFgYGBv/5z39evXq1ZMkSKpV69uzZoqKixYsXx8TEZGZmYtuE//XXX+Xl5diO4/yePXs2bdo0/qXk6HR6VFQUfx18K0UzMzN8z993796FhYVNnToV2/rP3t7+xo0bVVVV/NtcDX08AAD4/xkZGVH/fytWrKitreXxeCkpKVQq9dy5c1jNmpoaY2PjAwcO0Gi0X375JSMjAyv/559/qFTq7du3sZfffvvtzp07u15o6dKlP/30E4/Hq6yspFKpW7duxcp7ugqPx0tNTaVSqSkpKdhbwcHB2IXKy8upVGpsbCxWXlZWhpX3cioGg0GlUnfs2IG9lZ2dTaVSb926lZOTw3+q1NTUH374oa2tzczMbM+ePTwej81m29ra4gfys7KyOnz4MP7S3t6e2p3vv/+e/6iKigpbW1tnZ2csg/J4vPz8fP6POUxACwkA0I158+ZhKaSpqSkvL++PP/5Yvnz5w4cPCwsLEUL379/Ht6AlkUjv3r3D9yq9e/duaWlpXl4eQojNZn/RRSdMmID90NNVEEKvX79GCE2ePBkrnz59enBwcC/n7OVUnS6KtUWYTCYWPL5Bg62tLba34ezZs+Pi4vbs2ZOens5gMBYvXtzpWrW1tfX19dg+IDg7O7tz587xl2A7V+EaGxt/+OEHLpd7/vx5RUVFrBDbTK+wsHD69Om9fLohBhISAKAbsrKy2K4WWlpa48aNq6ioOHHixJs3b5hMJkLI0tIS3wTB2tpaXV2dTqcvXbq0paXF0tLSyspqzpw5fn5+vC8cKTBixAjsh56ugv5NcnV1dTIyMgghZWVl/jNwOBzsh7a2ts+eCoOdByFEJv/fM3Vsb3is64zf4sWLY2Njk5OTExISRo0aZWVl1akCdvVOn5pEIuFn7qq9vX39+vUVFRWRkZFYEur24wwTkJAAAJ9XW1uLEJKRkcGGw40ePdrDwwMhxGQyDx06pKWlde/evdra2ujo6EmTJiGE4uLiEEL8g5t7GejcVU9XQQgZGBgghB4/fowNIn/69Cl2iKSkJEKooaEBe1lcXPzZU/WEQqEghF6/fm1vb48QSk9P37Bhw4kTJ6ZNm6anp3f79u1nz555eXl1HZOtqampqqpaUlLSx4/J5XJ37Njx8uXLkJCQThvRYvEP0O60IgsSEgCgG4WFheHh4QghFotVUFBw9+5dMzMzfX19HR0dHR2d8+fPk8nkSZMmhYWF3b9/f/78+Vjj5vHjx+rq6vn5+f7+/gih1tZW7GzS0tK5ubnx8fGzZ8/uy9UdHR27vQpCaNasWbq6umfPnpWUlJSQkMC2VUUIqaioKCkpXbp0ady4cS0tLQcPHvzsqXoyffr00aNH+/n5rV+/ftSoUUePHtXS0rKxsSGRSIsWLTp+/Dj2Q7fHjh8/vu8J6ejRo/i/CT6q3sLCQllZGTsJ3p04TEBCAgB0RiKRsrOzs7OzEUISEhIaGhrLli3btGkTiUSSlpYODQ319fXFvvHHjh27e/duS0tLLpc7f/78sLCw06dPKygo/Oc///nf//6XlZW1Zs0ahNDChQuPHTu2adOmp0+fqqmp8V8Ia2fw/xch1NNVEEKSkpLnzp3bunWrr68vmUy2tram0+kIITExsV9//XXPnj3r1q2TlpbesGEDljl6OVWni+IhSUpKnjlz5pdfftm7d6+YmJiNjc3GjRuxbjdXV9cTJ05MnTq1a/caZsKECVeuXGGz2djwuW73j8cLX7x4gRCKj4+Pj4/H3w0PD7e2tn79+vWoUaPwPszhQsiDKgAAxNTQ0FBZWdmpsL29/cOHD1wut2v99vb2xsZGgVwFU11d3dDQgI+mwwo5HM779+87Ojq+6FQ9+fTp06dPn/hLSkpKqFTqvXv3ejrkw4cPEyZMuHHjxhddqJOGhoZJkyadP3/+a05CRNBCAgD0h6KiIj4kDCcpKTlq1Khu60tKSmKPeb7+KhgNDQ2EUH19PX8hmUwePXr0l56qJ/wjJqqrq1+8eHHt2jV1dXUnJ6eeDhk1apSHh8f//vc/V1fXfi/8c/nyZTU1tU7TtoYDWKkBAEBgJBJJTEys64g4gauqqtqxY8ebN28CAgLw2azdwmbsdl1pqY+YTGZERMTOnTv5Z9cOEySegFZwAgCAoa26ulpFRaUveaKmpkZCQkJJSakfV2Gz2eXl5dhIv+EGEhIAAACRAM+QABgo7e3t7969KyoqKioq6ujoGD9+/Pjx4w0MDAahfwkAIoKEBIDgcTica9euBQcHY9uGSkhIkMnk9vZ2hJCMjMzChQs3b96soqIi7DABEC3QZQeAgD169CgwMLCjo2POnDmOjo66urpqamo8Ho9Op5eUlOTn50dFRTU1NW3YsMHd3X0YPrgGoCeQkAAQJB6PN2vWrK1btzo7O/e0ghmHw7l//35oaOj+/fvxRUIBAJCQABAwDofTl6dE2EzAXpbdBGC4gYQEwIBjsVjPnz9/8eKFmZnZ9OnTpaWlhR0RAKIIEhIAAwtb5K28vFxfX7+srExeXv7Ro0fw6AiArqC7AAABKy8v53/59OnTT58+vXjxIjY2NjMzU0lJ6cGDB8KKDQBRBgkJAEHi8XizZ8/et29fVVUVViItLd3Q0HDz5s3s7Oy7d+9WVlbCgG8AugVddgAI2KtXr4KCgjIzM1esWPHTTz+pqqoeOXLk5s2bHz9+VFNTmzdvnq+vr7BjBEAUQUICYEC8ePHi+PHjBQUF7u7u33//vYKCQmVlpY6OTr9XgAZgyIOEBMAASk5OPnHixPv37z09PT08POTl5YUdEQCiCxISAALG4/Gys7MTExNlZWWdnJwMDAzi4+NPnDjx6dOnH374YdWqVTIyMsKOEQBRBAkJAAE7ceJEaGiokZFRW1tbaWlpcHDwrFmzOBzO7du3Q0JCmExmQEDAN998I+wwARA5kJAAELApU6YcOXLE3t4eIfTnn38+fvw4PDwce4vNZsfExMjJyS1YsKCnw2k0Wmpq6sSJE83MzPDCoqKitLQ0XV1dOzs7WCwcDFWw2jcAAqaoqHjjxg0Wi8Vms//++28DAwP8LXFx8ZUrV/ZyrLe3d1xc3JQpU06ePKmhoXHr1i1xcfHQ0NCjR49OnTq1oKDA2tr65MmTA/8hABACaCEBIGAZGRknTpzIysqSlZW1tbX95ZdfdHR0+nJgWVnZnDlzIiMjLSwsGhoa7Ozsfv/9d0tLSycnp99++23BggUFBQULFy6MiooyNzcf6E8BwOCDFhIAAjZlypTLly+3tLRIS0t/Ufcak8n09va2sLBACCkqKqqoqLx+/VpGRoZMJjs7OyOEjIyMDA0NU1JSICGBIQlWagBAkHg83qVLl9hstpycXC/ZCBvjUFZWxl9oaGi4bt067Oe4uDg6nW5nZ1deXq6uri4u/n9/O2ppadXU1Axc/AAI0XBJSMHBwcIOAQwXDx48mDt37tWrVz99+tT13ba2toiIiFmzZh0/frynSbIpKSm7du1as2aNjY1NXV0d/+rgMjIyLS0tXQ9xd3eHmxwQ3XDpsgsJCdm4caOwowBDH4lEunr1akJCwpEjRw4cOGBmZqarq6utrU0ikUpKSkpKSsrKyshk8o8//ujl5dXtPhTXr1/ft2+fl5eXt7c3QkhFRaW5uRl/t6mpydDQsOtR6enpERERA/e5ABgEwyUhATCYnJycHB0dY2Jinj9//vr167t373K5XAMDAxMTkyVLlsyZM0dTU7PbA/9fe+ceF3P+Pf7XNEnlNpIkqdBOdFFC8tmhlJUvpU2LjVqVDWs367ZaYYnNUotSSKtV0sU95FIf3ZRraCPVNsWUWw0rlW5qZn5/vL77/s6vmTJNM++5dJ5/eMyceb9frzM57zmvy3mdc+zYsZCQkJ9//nnp0qVYMmLECDabXV9fP3DgQB6PV1FR4ejoSOJXAQDy6C1RdiYmJn///bestQB6Ke3t7W1tbZ9M0FBYWLho0SJXV1d3d3cs0dHR0dPTmzFjxtSpU9etW5eYmJiYmHj58uXhw4d3uBcsHFACYIYEAFJHVVWViErogpMnT/J4vJSUlJSUFCxxd3fftWvXnj171q9ff+nSJT09vZ07dwp6IwBQDmCGBADKAFg4oATADAlQVLKzs3Nycuzs7Ozt7WWtCwAAEkDew77fvXv3559/flLIZDLj4uKysrI4HA6J2gGyJCcnZ/v27UFBQdnZ2bLWBQAACSDXDqm5udnf3z86OrprYXR0tLOzc0ZGRkBAwNq1a0lXE5ANeG6E50my1gUAAAkgv0t2BQUFGzdurKqqGjx4cBfC6urqyMjI0NBQItNXQUEBJFbpDdjb2+fk5Mjz9IjH4928eZPJZLa0tPDLvb29oSQSAAgivw4Jp5hcsmRJVFRUF8Li4mLI9AXIJwcPHhSaPcHDwwMcEkAOirXVKr9LdsHBwXv27Bk0aFDXQtEzfUFiFYBkkpKSbGxsMjMzi4uLS/ig0WiyVg3oLeCtVkVZ1iZphtTa2lpeXs5kMplM5sePH83MzMzMzEaPHt1F9sm+ffuKIhQx0xdCCFIHSRbFGnmRD4fDqaur+/LLL0eMGCFrXQBAMZC6Q+JwOKdOnYqIiPjnn38QQn369FFRUWltbUUIaWhouLq6/vjjj1paWmK3L2KmL0Di4JHX9u3bwSEJhUqlWlhY3Lhxg0i7AABA10jXIWVmZoaGhn78+NHNzc3BwcHAwEBbW5vH47169erp06clJSXJycmXL19etWqVl5dXnz59xOgCMn3JCiLIDSZJnTFr1qzdu3cvW7bMzMyMf9MIghoAQChSdEg8Hi84OHjt2rVz5sxRUfm/zSoKhaKvr6+vrz99+vRvv/32ypUr0dHR48ePnzRpkhi9MBgMGo22Y8cOnOmroaHBwcFBcl8C6BQc5IZnSOCQhJKUlIQQysvLy8vL45cvXLgQHBIACCJFh0ShUNLT07uumEmlUl1cXJydnbvIYCS0ZgwhVFNTg0xfgHxy7do1oXL+8RkAAATSXbIT9EZtbW137ty5e/eupaXltGnTcDwChULprFKZm5ubm5tb10IGg3H37l2JKg4AEkBFRSU/Pz8qKqqsrKy5uZlOp8+fP/+rr76StV4AIKeQeg6Jy+V++eWXL168GDVqVEJCQv/+/TMzM8XbOgIA+Sc7O3vFihUqKipWVlaampr5+fkPHjx4//79t99+K2vVAEAeke7SwYsXL/jf3rx58927d3fv3k1JSXnw4AGNRrt69apUFQCkB4vFYrFYstZCrtm3b9+YMWNu376dlJQUExNz7949Z2fnQ4cOQcZFABCKFB0Sj8dzcnL65ZdfqqursURdXb2uru78+fOFhYWpqamvX7/uScA3IEOCgoJGjRoVGxuLEIqNjQXPJEh7e3t5ebmnpydxDFZdXd3Pz6+xsbGkpES2ugG9BBaLhXNrKcoTKkWHRKFQEhISqqqqvvjii+Dg4Ldv306aNMnHxyciImLhwoV79uyZP38+g8GQngKAlPDx8dm+fTvxlsVizZgxQ3bqyCmqqqr9+vXrUKMIv9XU1JSRUkAvIjs7e9SoUdghxcbGEq/lGenuIVlZWcXGxt69e3f//v0zZ8708vJavnz5+vXrX79+raen11kgAyDn4IkRP3ggBsHfHZgxY8b58+dpNJqdnd2AAQNycnJOnDhhaGhoZGQka9UA5afDMJHFYvn4+Dx79kxW+ogCGeGnU6ZMSU5ODgsLy8vLc3BwOHTo0KBBg8AbKSidDbIUJVkWmWzcuHH8+PFRUVEeHh7Ozs6hoaEUCmXfvn0Q9g1IG8FRI+JbwZNbpDtD4vF4hYWFGRkZmpqas2bNOnfuXFpaWnh4eHx8vJ+f35IlS+B4oMLR2eje0NCQXEUUAG1t7ePHj9+/fx9XoDA0NGQwGPypFwFASlRWVspaBXGQrkM6cOBAdHT02LFjm5ubDxw4EBERMXv27C+++OLixYuRkZHHjh0LCQn5/PPPpaoDIFmMjIxwxqAOcm9vbxloI/eoqKjY2NjY2NjIWhGgd2FnZydULufLxdJ1SCdOnDh06BD+0/z555/x8fEzZ86kUqlubm4uLi6nT5/GGVd7IQqdKvvYsWMzZswg4naMjIyOHTsmU43ki+Tk5Pb2dk9PzxMnTtTX1wte4OPjA2sDgFTBCb06DBy3b9/eqx3SwIEDz54929bW1t7e/t///nf06NH/17GqqoeHh1R7l2eILHAIIYXzSUZGRllZWdinxsbGent7y+QryG1Ia3Jy8ocPHzw9PePj44XqtmjRInBIgLTx9PR89epVWVkZQsjIyMjb23vbtm2yVuoTSNch7d69Ozw8fPXq1ZqamgwG4/vvv5dqdwoEkSpbQTOTYvuW4Tr1jBkziNEf3r+Vn1laSkoKfnHhwgVVVVX+BFoUCqW1tbXrBI8A0HMeP36ckpJCp9MtLCzOnj2rEN4ISTvKbvLkySdOnMjPz797925YWJienp5Uu1MgFNQPSZbs7OygoCAxwn5iY2M73BUbGys0rEi2ODs7X79+ncLHvXv3xo8f//z5c1mrBigzbDY7MDAQIbRr1y4LCwtZq9MNpFt+Ij4+fvHixf369eviMg6Hc/nyZUtLSwjT6m2IXeIvKChIUBgXFycngRVpaWmnT59GCL158yYqKurMmTNYzuPx8CmQwYMHd93Cu3fvUlJSfH198dtnz57hhReMpaWlrq6uVFQHFB82m71p0yaE0Jo1aywsLM6dOydrjbqBdJfsrl69Gh8f7+vr6+TkJJglqLm5+cyZM3/++SdCKC4uTqqaAHKIZEv8yc9OEoPBSE5Obm5ubmtra2lpISoaUygUAwMDb29vIpmQUJqbm/39/SsqKgiHFBcXd+7cuf79++O3O3bsAIcEdEZYWBibzXZ0dFTEUqXSrYeUlJSUnp7++++/79y509LS0sDAYPjw4RQK5enTp0+fPq2srFRRUVm+fLmvry8czlA4srOze3jITuwSf0ZGRoLuR36WQPv164c3tDZt2uTm5tatmO+CgoKNGzdWVVXxz6JKS0t9fX3XrFkjeV0B+aa74biBgYGPHz+2sLBQUGuRevmJWbNmOTg4nD59+s6dO0VFRampqVwud/To0aampl999dXs2bOHDRsmbR0AaZCTkyOrU9/Hjh0bNWpUB+HSpUtlokwX/Prrr1lZWXFxcW1tbQghDofT0tJSWFi4b9++ziZJCQkJ1tbWS5YsiYqKwhIul4sdUklJiaqqqrGxMWQ56T10a1mb8Ea7du2SvmpSgYx6SDjCGwd5t7e3t7W19fKYVxaLFRQUhDfh5WehqbvY2dnhFKudHcGTHkZGRs+ePSPOQuGDUPIzQyIICwuLjo7uIBw+fHgXUXbBwcF9+/Y9f/48Iamqqmpubg4MDMQLgLa2tlFRUUKfoIiICH9/f+KtQp91AzCiL2srgTdC5OSy40dVVbWXeyOcgpcICYuNjfXx8ZGpRmJib2+/bdu2bdu2yeT3Dvsk7BFldRDqk1y8eHHq1Knp6ena2trBwcEnT540Nja2t7cfMGBAZ7f07du3g6SpqYnBYERGRhYUFBw+fPjevXuCTg7D743Qv4NryDGo0BDnW7v+fwwLC1MCb4TId0jd5d27dzjqgYDJZMbFxWVlZRFVzp49e5bGB1F+ST4RdD+CQcyAEsDhcGpqatzd3Q0NDc3NzcvKyqysrEJDQ0+ePCk0fUNnmJqaxsTE2NraqqioODg4TJo0qaCgQJQbicE1WJdyExYWlpGRoaOjo+jeCMm5Q8LhRvzjwejoaGdn54yMjICAgLVr12JhXFzcTz/9FPQvRUVFMtL303RWZVW8YazY53gAEqBSqTQarbGxESE0evTohw8fIoRGjRrF5XKfPn0qejs3b97ctGlTe3s7fltXV/fZZ5+JcqOIg2spAcZJDuvXr8feKCYmRugFeGmd/HV18SBjD0k8BMONqqurIyMjQ0ND582bV1pa6urqWlBQMGHChF4bgyT2OR6AHCZMmHDo0CEdHR1TU9O4uLirV6+WlZWpqqqOGTNG9EbGjBlz6dKlkSNHenp6Xrx48e+//xZ6DEveUOjkWIpCYGBgWVlZ1yt1inUGn6QZUktLC7HChuFwOOXl5cS4TxAcboRPeGGKi4tVVFTmzJmDEBo7dqyJicmNGzdwDJKpqWlJSQmTyeTxeCKqJJMRHE6VLSiXyfhFtmPY7OzsGTNm4E0gxY3s6JrNmzerq6unp6fPmjXLyMhozZo1hw4dcnR07GIPiYAIpdPV1d22bduVK1cmT54cGhq6YcOGCRMmSFlxCUAsGPZkfgbTrM7AuRiUY9+IH5IckrOz83//+19+yYMHD+bOndtFDpXg4OA9e/YMGjSIkLx48WLo0KGqqv87q9PV1X3z5g0Rg/T11187Ozt7e3s3NzcLbTAiIoL/LR7BkW/ughmlZLUnL8NN79jY2A7J6BRi1N9d9PX109LS1q9f37dv35MnT4aEhMTFxR04cOCTN7q5ud2+fZt4u2DBgtTU1L///ruwsNDPz0+aKksMiQzMIS5DKGw2W2miGDog3SW7nuRQEQw3qq2t5T8/q6Gh0djYiGOQ/Pz8bGxssrOzv//+++jo6B9//FGwwQ4xSLJKb2pvb//s2TMi7Nve3l5+soKShtDIDvIDlKUUGN1F+YmXL18+fPgQyk8AYsNms5ctW4YQUj5vhKTtkHqYQ6UDWlpaRAsIoYaGBhMTExyDhCXdikHCaQJkshqAz80YGRn1zu2fzoork62H1DbhlKD8BJvN7uwjHR0dMjXpCcp3EisjIyMsLAwh5OHhsXjxYlmrI3mk65B6kkNFkBEjRrDZ7Pr6+oEDB/J4vIqKCkdHx5s3b6ampu7cuRMv5dXV1U2ZMkUy2gPSobMSYUqzk0SUnwgJCbG0tBTc1yQh1UJ2draPjw/+k2ZnZy9dulTwz469zuPHj2tqathsNpvNxi86XKajo9NBSPikYcOG4dcWFhY6OjrEWzlByQIrEhMTk5KSEEJr1qxRxDx1okBSlJ0YOVQEYTAYNBptx44d69atS0xMbGhocHBwoFAoihiD1HM6iyCXfzpzSErwk9GBH3/80dTU9NChQyT3iwNGOrzFi+RsNvv69etFRUWPHz/uuhHCtXz8+JF4jT0T4Z+IFxkZGfx3WVhYNDc3y7x8oqJXHeMfVURERFhaWuLQGMWqKNEtSHJIYuRQISCGk2pqanv27Fm/fv2lS5f09PR27tw5fPhwhNC2bdvi4uLCw8PV1dUVJQaph/Dnd8ArYApRfQuDVyw7bCPhNUwZaSQVuFxuS0sLEYNDJoJbdCwWa9q0aYKDPx0dHRqNNnHixI8fP06cOBHxTXq6APshYjqFfdvz58/LysqwZO/evS9evMAXx8bGlpeXu7u741kUmT+mMlyW7zkdRhX//PPP7du3CcevrJD0tOAcKkFBQYsXL167dq2xsfHmzZsnT578yfhXNzc3Nzc34i2Dwbh7926HaxYsWLBgwQLJKy2v4HET8ZbFYm3fvt3Q0FD0akD8Iy+Z+DMcWIhDHLH+clLKSIKoqKgsWrQoKirqm2++sbCw4HcGnp6eUt1DEjopYbFYVlZW2CXgf8X2DdhjEX6Lf/mIzWbHxsampqbyX5+Xl8fhcPCPKV7ZMzc374kCCoTY+1iCKz3Nzc1BQUEKNPQUAzIcEs6hsmHDBiKHyldffRUaGuru7r5u3bqBAweSoIO80ZOELkJrRwUFBYn+m84/8hLDn0kEPE8KCgrCR5GUkqtXryKE7t6922EUNX/+fPKDGiwsLGJiYqS9x6Ojo3P48GFBeVlZ2fLlyx8/fow3qx4/fpyUlESs71lYWCjrpojYgTNCfx+ys7PBIfWUDjlU8vPzEV8OFSsrKxJ0EES2+bZ7spLQQ4WF5nLtlj8DROTq1atcLldQ3qdPH6n26+3tLRjKuHHjRhlGHBAxytgbYc+E1/oyMjIyMjISExPxZY6OjnjmRAxWxJvEs1isuLg4Rcypj48ZaWhoCB6pVLJlbUFIWrKTSA4VScFisYjKBZPOEacAACAASURBVAih2NhYIyMjkscdEq/d0MM9W8V6YhUFKpUqyi6pxNm2bRteCyUkZB6+Flo+kfgl1dHRIYqZduacdHR0Lly48M8//+Bb8CSexWKJfmKPxWLxV8ySyTMuHkQo3ZAhQ4h9OAI5LPolWUhySJs3b/b19U1PTw8KCjp8+DDOO+fk5CRKDhWJQ2yfEODkgyQfjxW7O/xz00EouqUq/SBLfuBwOFlZWS9fvuxJcKkYGBkZZWVlZWdnx8XFZWdnb9++nczfYtHtszPnlJ2dTXgjgtjY2KVLl4r41MjJyetuQbgihJCHh0dERAT/uBkpWlY68SDJIeEcKm/fvsU5VDIzM4cNG2Zra0tO7x0QulaWk5OjKP/Z9vb2eFWakHRrhXrp0qWC2zZKvJEjQ3oSXNpDjIyMvL29KysryY8xw/YZGxvbrfKJ/M4pNjaWP3MSQVhY2F9//SVKtJ7gt2axWGI84+QcrU1MTHzw4EFZWRlCSEdHhwjsluGoQlaQV36ipKTkt99+mzdvXkJCQnt7u6yqFnW2NqVYa1bbtm179uwZ3vXprqXiDBH8ElxqT6IKAgiJVaBPOdi2bVtWVpbY5RM7m8Q3NDTgVAWBgYHLli3DdYAEL5PgsyzVZHpsNjsxMXHZsmVJSUllZWXYFcXExBC+Fo8qFGWgLBFIckh5eXnu7u6ZmZkvX75saGiorq4OCAjYuXMnOb3zI1f5tnuCkZGR2ItvHfxZVlaWBBUDMJIq0Keg9GRlWKht29vbJyUlEbMHNpuNnZOLi4ugcxLau6GhodgqSRbsijw9PZOSkthsto6OjoeHR0xMjLKGGooOSUt2ISEhxsbGp06d+uGHHygUyqpVq9ra2hISEgIDA8nf9RVc48YjEZLVkC098WeSoofFTOW8Fqp8BpcqBHgPjD8qAYckEMt6OMvR48ePccoJ7Jywf8KreQwGo8M8yd7eXubPOJvNvnLlSm5uLj4+PHLkSITQ4sWLwQ8RkOGQ2tvbmUxmUFAQcfaCQqEsXLjw8OHDT548GT9+PAk68GNvb5+VlUUUnsBvSdYBQP/Gvtvb24s3Pe3J7SwWi4j7l97RYLkKLlUsjIyMiKT4gk+ojo4OsY0k1DkhhCwtLV+8eIGDI8aNGyfGM85vIUKTAYreQnZ2dm1tbUVFBaG/o6PjzJkz5Sr1nzxAhkNSVVWl0WhMJpNfeOvWLYQQf7kjMsHxKvigQ69aopUriOLK4v0X9OR2Hx8fYnYlRlSxiMg2uFTOZ5CfhJjEd/3/25lzQgiNHDmyrKysrKyMw+G4uLgQl4kSFsFfrwsfFBElNIOAzWYvWrSIaCE7O/vu3btWVlbjxo2DKVEXkLRkN2vWrFOnTlEolJqaGhqNduDAgeTkZGNjYwMDA3IUAOSQHoaxin17bGys4C91t6KKReTJkyfp6elv3ryRSXCpbNO4ycQddnBOCKHt27fjeAH077YT+jcVLJEnAv3rooi3ghbCYrF8fHxwglpBCEeIEMLxcjj1HP81zc3Nr1+/xgNxoDNIckgBAQHV1dU45w2eKhkZGe3fv5+EPPwAICtCQ0PXr1/v4ODg7u7OYDBcXV3J7F3ih6+7hcyzmmIHM2zYMITQzJkzv/vuO4QQ/zlc7LH4XRRxY3p6umCDeCaND5AR6WWfP38uWEoUISS0FjbO0C/zvVt5hiSHpKmpeeTIkeLiYiaT2dTUZGhoOHnyZGknUFF6iDUrmdyuuHT2iyDxX4ojR46kpKRcunQpLS1NR0fH1dV1/vz5o0ePlmwvnSHbc5SydYeCYP9EHHVC/05riH/5X2hqagptJDc3l/8jHR0d7I34k/LhXpQ7SaP0IDU3vqmpqampKZk9KjeyWvJSdPDmhGAUlsQd0pgxY9avX79u3bp79+5duHAhOTn5jz/+mDBhwp9//tnZT57SIA/W1fWyIV7fE5Sz2ezi4mL+BMQYe3v7wMBA4l78orNqHULdsBihrYq+EdhdyHNI+fn58fHxxHrd119/LSdDJ6C3gaOK+fOySDXRGYVCmTBhQn19fUNDQ3p6+qNHjz5+/Kj0DkkeEC8OEzuqDslQEELbtm0T3cVif9zBl4gRNSPzlU+SIckhXblyZe3atRQKZdy4cRoaGjk5OZmZmWvXrl25ciU5CgAAP9gndRZVLCl4PN79+/cvXrx47dq1+vp6Y2PjjRs3urq6SjWRHUDQkzhMnHwPJ0Cyt7cXo4DksWPHcLpxvG/UrSA9Anlb+ZQ2JDmk8PBwU1PTP//8c/DgwQihDx8+bN269cCBA4sXL5ZVPaTeNhcGOiBiVHFPWLRoUWFh4YABA+bOnTt//nxLS0spdaSs9HCbsyfLhsSkGZ8MEWM5l78FsTMAycPKJ5mQVKCvqqpqx44d2BshhPr37+/v73/lypXy8nJra+su7n337l1KSoqvry8hYTKZt27dMjAwmD59OpHlQaiwa3p4KhMAPsm4ceO8vLy++OILdXX1bt0ootkrPb3t5xggI5cdlUql0+kdIinLyso0NTW7PofU3Nzs7+/Pny85Ojra2dk5IyMjICBg7dq1XQg/CR58dWtdGAC6RVBQkIuLS3e9kYhmD5AMnFEhAenOkE6fPv3u3TuEkKGhYVpa2pw5c+zs7LS0tEpKSv773/9+/vnnXRxZLygo2LhxY1VVFTGvqq6ujoyMDA0NnTdvXmlpqaura0FBwfDhwwWFEyZM+KRuMPgC5BARzV4UCwd6Tq89GiErpOuQEhISSkpKiLcVFRVENieEUFZW1ocPH4QeK8P3WltbL1myJCoqCkuKi4tVVFTmzJmDEBo7dqyJicmNGzcsLCwEhfC4AgqKiGYPFk4O/MNWHo8nU116BdJ1SAcPHhwxYgSHw+nsgi5Ww4ODg/v27Xv+/HlC8uLFi6FDh6qq/q/Ourq6b968ESqUkPoAQDYimr2MtAMA6SLFPSQejzdr1ixra2t/f/9Tp069fv2aKkAXtwvOnGpra/mX4zU0NBobG4UKhTYYERHRWV8UCkXoAjHIlVsupbP0Xl5emZmZCKHbt2/jHDOiI6LZC723CwsHAIVAig6JQqHcu3cvNDRUR0fn6NGjjo6OTk5Ov/76a05OTnNzsxgNamlpffjwgXjb0NCgq6srVCj0dn9//85a5vF4QufjIFduuZQcUkFBwZMnTxBCW7duffjwYQ9bk4iFA4BCIN0lu379+hHJo1gs1o0bN3Jzc8+cOcPhcCZNmjRt2rT58+eLfkhwxIgRbDa7vr5+4MCBPB6voqLC0dFx+PDhgkJpficA+ATm5uYHDx68du1aTU3N3r17Y2JiOlwQExMj+vE7oWYvaZUBaQFhEd2CvNRB+BziN99809raev/+/dzc3LNnz5qbm9vY2IjYAoPBoNFoO3bsWLduXWJiYkNDg4ODw5AhQwSFUv0iANA1+/bti4yMrK2traysHDBggLa2docLuhVALNTsJaovIEUgmrdbkOSQYmNj6+rq+CUaGhpz586tqKgYNmxY17XuiadXTU1tz54969evv3Tpkp6e3s6dO4cPH44QEioEAFmhp6e3a9cuhJCvr6+3t/f06dPFaOSTZg+QDIVCgUA7aUPSn9jDwwMvpquoqFAoFP64OzU1tYiICGkPIkxMTP7++2+pdgEoHLhGAD4fLaUu8vPzo6KiysrKmpub6XT6/Pnzv/rqK2l0BBYubcAhkQAZmRoQQi4uLjQaLSQk5NatW3fu3AkPD9fT03N3d09PT588efLmzZvJUQMAyCQ7O9vT0/PWrVv6+vqWlpZFRUWbN28+evSorPUCxAG8EQmQ4fO5XO6UKVN8fHxWrVpFCOPj43///feHDx+WlpbOnz8/IyNDX19fejrA+BEQJDs7OycnR7xs0KIwb9689vb2xMREHLnT0tKyefPmrKys/Px8ieejAwsHlACSZkgUCqW1tZVf0tTU1NLS0tLSgs9YtLe3k6MJABDY29tLL5lhe3t7eXm5p6cnEUeqrq7u5+fX2NjIn74EAAACMoIaVFRUGAzGiRMn+vfvb2dnp62t/eDBg2PHjuHxY3h4+IgRI6DOPKBkqKqq9uvXr8OsBb+F6nwAIBSSouy2bNny5s2b33///ffff8eSUaNGbdq06caNG2lpaZs2bSJHDQAgkxkzZpw/f55Go9nZ2Q0YMCAnJ+fEiROGhoYw/AIAoZAXN8LhcB49evTkyZOWlpYxY8bgsi5sNpvD4ZAQxgor7AD5vH37ds2aNfn5+YRET08vIiLC3Nxc4n2BhQNKAEkOicfj3bx5k8lktrS08Mu9vb01NDRIUAAeV0AmcLnc+/fvY8s3NDRkMBjdLY8kImDh0gbCvkmApCW7gwcPCs386OHhQY5DAgCZoKKiYmNjI3o6EgDozZDkkJKSkmxsbHbv3q2rq8ufN0VFhaQwPwAAgJ4A0yMSIMMhcTicurq6L7/8csSIESR0BwAAACgiZExQqFSqhYXFjRs3SOgLAOSHlpaWDtUpORxOeXk5nLoDAKGQtGQ3a9as3bt3L1u2zMzMjH/TiLSgBgAgH2dn5w0bNsyePZuQPHjwwMvL69q1a6NGjZKhYgAgn5C3h4QQysvLy8vL45cvXLgQHBKgZKSlpZ0+fRoh9ObNm6ioqDNnzmA5j8d79uwZQmjw4MGy1A8A5BWSHNK1a9eEyiGoAVA+GAxGcnJyc3NzW1tbS0sLUe+VQqEYGBh4e3uLXpQSAHoVJDkkcDxA76Ffv37Hjh1DCG3atMnNzQ1ivpUDOIdEAtJ1SMnJye3t7Z6enidOnKivrxe8wMfHB5bsAGXlt99+k7UKAKBISN0hffjwwdPTMz4+nsViCV6waNEicEiAsiLzBCWABIHpEQmQNwkV2hH/IVmpAolVAPKJjIwUmqDk7t27Et9GAgsHlACS9pAQQn/99VdiYmJJSQmHwzEzM3NxcbGzs+tuI1wuNzMzs6qqatKkSRYWFtifPXv2rKysjLjG0tJSV1dXkqoDgFhAghIA6BYkOaRLly4FBAQghCwsLIYNG1ZYWHj58uVNmzZ98803ojfC4/GWL1+el5dnbW0dFhbm4uISHByMEIqLizt37lz//v3xZTt27ACHBMgcSFACAN2FJIcUHx8/fvz46OjogQMHIoS4XO7PP/8cFhbm4eHRp08fERvJzMzMzc09dOiQo6NjZWWlk5OTo6Ojg4NDaWmpr6/vmjVrpPkNAKB7EAlK3N3dZa0LACgGJOWyKyoqCg0Nxd4IIaSiouLn53fhwoXi4mJLS0sR2yktLaXRaLjgtKGhoYmJSVFRkb29PXZIJSUlqqqqxsbGpO1LAUDXQIISAOgWZDgkKpVqbGycmpo6d+5cQpibm0ulUruVQMXMzKyhoaG2tlZbW7u1tfX58+csFquqqqq5uTkwMBAfQrS1tY2KioKnHZAHIEGJMgHnkEiApCU7R0fHQ4cOOTk5TZ8+ffDgwbdv33748KGRkVFCQgJCiEKh+Pr6qqmpdd3ItGnTaDTaihUrXF1ds7KympubW1pampqaGAyGn5+fjY1Ndnb2999/Hx0d/eOPPwreHhER4e/vL5WvBwDCgAQlANAtSPL5s2bNqqys7OKC/Px8YkGvCx49ehQREfHPP/8wGAwWi6Wnp/fzzz/zX+Dl5UWlUmNjYzvcCEGxgEwoLi4+evRoeXn5nDlzhg4dSqVSv/zyy+42IkocKVg4oATIOJcdgShjRjabXVlZGR0djefOjo6Otra2N2/eTE1N3blzp6qqKkKorq5uypQpklEaAHpGXl6en59f3759qVRqQ0MDh8M5cODA48ePt27d2q12II4U6CWQl8uuoKCgh+eQqFTq1q1bm5qaXF1dY2JiWltbZ8+e/fHjx0uXLo0cOdLT0/PixYt///13UFCQlL4FAHSLkJAQY2PjU6dO/fDDDxQKZdWqVW1tbQkJCYGBgVQqVfR2II4U6CWQtJZ96dKlJUuWXL58uV+/fp999llhYeHKlSuPHz/erUaGDBmycePGyMhIKyura9eu7du3T0tLS1dXd9u2bVeuXJk8eXJoaOiGDRsmTJggpW8BAKLT3t7OZDK9vLyI+AUKhbJw4cL6+vonT56I3g6Xyy0tLTU1NS0pKWEymbCvDigxinQOCSG0ePHixYsXdxAuWLBgwYIFklQXAHqMqqoqjUZjMpn8wlu3biGEBg0aJHo7EEcK9B7ImCHhc0heXl4dziE1NjYWFxeToAAAyIRZs2adOnVq165dNTU1r1+/PnDgwL59+4yNjQ0MDERvBMeRRkZGFhQUHD58+N69e9HR0UKvFJo3DwAUCEU6hwQAikVAQEB1dXVcXBxCCE+VjIyM9u/f362z26ampjExMfi1g4PDpEmTCgoKhF4JpxqkCpxDIgFFOocEAIqFpqbmkSNHiouLmUxmU1OToaHh5MmTu7VGjRCCOFKg90CSQ7p8+TJCiMVi8VdFqqioCAsLw68XL14MDglQDjqrS1lbW/vXX3+hbtalHDNmDMSRygMwPSIBkiahXC636wukfXYdjg0CpPHll19++PDh+vXrTk5OQutS3rp1a8iQIaI3ePr06bi4OCaTqa6u/sMPP/j5+QleAxYOKAG9ZVUUHldAJpBWlxIsHFACSFqy43A4WVlZL1++bGtrw29bWloKCwv37dsn8dKZACA/vHnz5uXLl/hs3NWrV+3s7DQ1NWWtFADIKSQ5pLCwMMFY1eHDh3frvDoAKBYPHjxYvnz5+PHjjx07hhDas2dPUFBQUFCQk5OTrFUDAHmEpEwNFy9enDp1anp6ura2dnBw8MmTJ42Nje3t7QcMGECOAgBAPnFxcXp6euvWrSPeTpgwITw8XLZaAYDcQtLB2JqaGnd3d0NDQ3Nz87KyMisrq9DQ0JMnT/LHIAGAMsHhcHJzc/38/CwsLLDE0NDwp59+qqio4E/dDSgKUPmTBMhwSFQqlUajNTY2IoRGjx798OFDhNCoUaO4XO7Tp09JUAAAyAeb/YMHD/iFODUJ1EMCAKGQtIc0YcKEQ4cO6ejomJqaxsXFXb16taysTFVVdcyYMeQoAADkM3369PPnzxsYGNja2qqpqZWWlu7bt8/IyAjMXhHpJQHJsoUkh7R582ZfX9/09PSgoKDDhw/jRPpOTk6whwQoMZs3b66trQ0JCSEkY8aM6W7qIADoPUj3HNLHjx9VVFRwyhMej/f27duhQ4c2NDRkZmYOGzbM1tZWel13AE5pALLi1atXJSUl9fX1BgYGVlZWUoosBQsHlADpOiQXFxcLC4tdu3ZJrwsRgccVII3OUgcRdCt1kIiAhQNKAElLdgDQe0hOTv7w4YOnp2d8fLzQ1EGLFi2CgkYAIAg4JACQMAMGDFi9ejVCaPv27VOnThVchIA9JAAQitQdUmlp6eHDhzv71NvbG4aKgJJRUFBgY2Pj4OCwdevWDRs2zJ49W9YaARIA6iGRgNT3kLo+A9jdtMdcLjczM7OqqmrSpEkWFhbESJPJZN66dcvAwGD69OlCN41hhR0gja+//vqvv/4aM2ZMVVWVrq6uYLbGmJgYonqypAALlzbgkEhA6jOkOXPmBAQEdPbp4MGDRW+Kx+MtX748Ly/P2to6LCzMxcUlODgYIRQdHb13794pU6aUlpba2toeOHBAAnoDgLjs27cvMjKytra2srJywIAB2traHS6AJTtFBLwRCUjdIWloaOjq6kqkqczMzNzc3EOHDjk6OlZWVjo5OTk6OpqamkZGRoaGhs6bN6+0tNTV1bWgoAAnVwYAmRAQEODj4+Pg4LB48eK1a9dOnjxZ1hoBgGKgSClMSktLaTSavb09QsjQ0NDExKSoqKi4uFhFRWXOnDkIobFjx5qYmNy4cUPGigK9m4KCgidPniCE2Gz2P//8I2t1AEBhkO4MydHRcfjw4ZJqzczMrKGhoba2Vltbu7W19fnz5ywWi0ajDR06FJ+9RQjp6uq+efNGUj0CgBiYm5sfPHjw2rVrNTU1e/fujYmJ6XCBNPaQAEAJkK5DwimCJMW0adNoNNqKFStcXV2zsrKam5tbWlpqa2vV1dWJazQ0NHAWV0EiIiL8/f0lqA8ACAX2kABAPBTpHBKVSo2KioqIiEhJSWEwGAMGDNDT09PS0vrw4QNxTUNDg4mJidDbwRsB5KCnp4ezk/j6+np7e0+fPl3WGgGAYqBIDonNZldWVkZHR+P4S0dHR1tbW11dXTabXV9fP3DgQB6PV1FR4ejoKGtNAQAhhP78809ZqwBIDAj7JgFFckhUKnXr1q1NTU2urq4xMTGtra2zZ8/u378/jUbbsWPHunXrEhMTGxoaHBwcZK0pAPwvxcXFR48eLS8vnzNnztChQ6lU6pdffilrpQBATlGkKLshQ4Zs3LgxMjLSysrq2rVr+/bt09LSUlNT27NnT25u7owZMy5fvrxz504JhlEAQE/Iy8tzd3fPzMx8+fJlQ0NDdXV1QEDAzp07Za0XIA4wPSKB3jIJhXPsAPnMmzePx+OdOnXqhx9+GDdu3Pr168PDwxMSEu7cuSPxIhRg4YASoEgzJABQINrb25lMppeXF5GtkUKhLFy4sL6+Hp9SAgCgA4q0hwQACoSqqiqNRmMymfzCW7duIYQGDRrU3dY+ma0RAJQAcEgAIC1mzZp16tQpCoVSU1NDo9EOHDiQnJxsbGxsYGDQrXYgWyPQSwCHBADSIiAgoLq6Oi4uDiGEp0pGRkb79+/v1sHY6upqyNYI9BJ63R5SRESErFWQCjL8XgrdtVSV19TUPHLkyPnz50NCQrZv337s2LHU1FQ6nd6tRsTI1qiURi7bLxURESGr/Bo9/OKKZQy9ziFFRkaKcZfo/6miXNn1NeJ9Kt73kggidi2NB4PoWuzG+ZWXhob5+fl79+79/fff9+7dGxkZeeHChe628OLFi+5maxTDGMi08K4vkEMLF7F3Kf3099DIO2gu5/6pt4R9e3l53bt3T9ZaAHKNjY1NfHy8BBvMzs5esWKFioqKlZWVpqZmfn5+a2vrTz/99O2334reSHh4+PXr1y9duoTf/vjjjyoqKvv37+9wGVg4IAoSN3LJ0lv2kOT5/wBQVvbt2zdmzJjExERcNLalpWXz5s2HDh3y8fERPVJOxGyNYOGAEtDrluwAgBza29vLy8s9PT2JEubq6up+fn6NjY0lJSWitzNixAicrREhhLM16uvrS0VjAJA14JAAQCqoqqr269evQ/YE/FZTU1P0dhgMBs7W+OrVq71790K2RkCJ6S1LdgBAPjNmzDh//jyNRrOzsxswYEBOTs6JEycMDQ2NjIxEbwRna1y/fv2lS5f09PQgWyOgxPSWoAYAIJ+3b9+uWbMmPz+fkOjp6UVERJibm8tQKwCQW8AhAYAU4XK59+/fZzKZLS0thoaGDAaDv8AxAAD8ULdv3y5rHaQIk8lMTU2tra01MDBQURG+YXbu3DmE0NChQ8lVTZK8e/cuKSmJnNP7QvtisVi3b98u56O+vl5PT0+yXbNYLBz9rKur2+EjLpebkZGRnZ1NpVJ1dHSEnmEsLi6+fPlyU1OTvr6+4AWVlZWXLl1SU1OTlCV8/PixsbFRXV19xIgR48ePV1NTGzNmjBhZ7Lqml1g4AiOXSyOXPDzl5ciRI3Q63cvLa/Lkyf7+/kKvuXLlCp1OP3HiBMm6SZCmpqbFixdPmTJFhn0lJCRM/RdbW1s6nb5q1SrJdr1+/XozMzNvb+9JkybNmTOnra2N+IjL5S5btszExMTDw8PCwiIwMFDw9o0bN1pZWfn6+lpaWq5fv77DpyEhIebm5j4+PtbW1j/88EPPtb106ZKdnV10dDQh8fX1HTduXEhICJfL7Xn7mF5i4Twwcrk0cmmgtA7p9evXFhYWFy5c4PF4JSUldDr94cOHHa6prq6ePHmyQj+uDx8+nDlzJp1OJ+FZFbGvU6dOmZmZFRcXS7BrFos1duzYBw8e8Hi89+/fjx8/PjU1lfj0+vXrdDr9+vXr+EoTE5OMjAz+20tLS01MTO7fv8/j8TIyMuh0+vPnz4lP6+rqLC0tz5w5w+PxCgoK6HR6RUVFT7TNy8szMTGZMGHC1atXCWFqaur06dPpdPqBAwd60jhBL7FwHhg5j8eTPyOXEkob9v3JDGA8Hu/nn3+eOXPmiBEjZKSjBEhISLC2tt60aZOc9PX69evdu3f7+/uPGzdOgl23tLRs2LDB2toaITRw4EAtLa2ioiLi09LSUhqNZm9vjxAyNDQ0MTHh/xQhpKend/z48YkTJyKEmpqaEEIcDof/go8fP+IXampqCKF37971RNuoqChtbe0rV67Mnj2bEM6dOzctLW3ChAkxMTHNzc09aR/TSywcgZEjhOTPyKWE0oZ9fzIDWFxcXFVVVWRkpIuLiywUlAzBwcF9+/Y9f/68nPS1d+/eoUOHdis1jiiYmJgQ6QmuXbv26tWr6dOnE5+amZk1NDTU1tZqa2u3trY+f/6cxWLx3z5gwAAbG5u3b9/u3r07KyvLw8PD0NCQ+HTgwIHr1q3bunXr5cuXHz16NGvWrB5uVJSWls6ZM0dwD0BdXd3Nze2XX35hsVg9/y3rJRaOwMgRQvJn5FJCaWdItbW1/OFMGhoajY2NxFsmkxkWFhYSEtKvXz9ZaCcx+vbtKz99vX//Pi0tbcmSJdKrIHfjxo2ff/75m2++mTp1KiGcNm0ajUZbsWLF8ePHV65c2dzc3NLSInhvY2Pj4MGDBw8e/ODBg6qqKn75lStXBg8erKenN3jw4FevXtXU1PRESRUVFSI7QweMjY170jI/vcTCERg5Qkj+jFxKKK1DEswAxj9iDQwMHDlyZFlZWVJS0ocPH/Lz8/Py8mShplJx/vx5VVVVNzc3KbV/5syZlStXenl5bd68mV9OpVKjoqK0tLRSUlIsLCy++OILoRXwDA0NN2/efPny5erq6jNnzhDy7Ozs1N20owAADI5JREFUkpKS06dP//rrr1euXHn9+vXJkyd7oufkyZNzc3OFfnTnzp2BAwcKTUbXXcDCZQIYuVRRWofUdQYwLS0tTU3NlJSUlJSUpqamR48eQabknpOWljZ79uz+/ftLo/Fjx45t3bo1ICBgw4YNHT5is9mVlZXR0dHnzp1bu3ZtUVFRh2c1Ly/vxx9/xEvqffv21dPTq6ysJD59+vSpjo4Ojt/t06cPnU7nH1qKweeff15cXHzgwIEOi/jp6elHjx6dNGlSZ/HZ3QIsXCaAkUsVpd1DIjKArVu3LjExkcgAdvLkyVGjRh05coS40sHBYdmyZUuWLJGdsgoM/nva2NgghCoqKpydnaXRS2Fh4Z49e1xdXceNG4d/WHV0dIyMjHDvY8aM2bp1a1NTk6ura0xMTGtrK380AULIxMQkKyvrjz/+WLJkSWpqamlpqbe3N6G8ra1tZGRkQkKCu7v7jRs37t27FxMT0xNtFy1aVFRUdPDgwYyMjMmTJ+vr67969aqsrOzOnTt6enqrV6/uSeMEYOGkAUZOHrIN8pMqubm5NjY2dDrd3t6eCKCcNm3arl27+C9zcHBISEiQhYIS49y5c7a2tjLpi/h7vnnzhk6n3717Vxqdbtq0if7/s2nTJv7eExISGAyGiYmJs7PznTt3BFs4e/bs3Llz6XS6mZlZWFhYB+VjY2MtLS1NTEzGjRsXFBTUc4W5XG5ERISzs/PYsWOxwnZ2dj///HNdXV3PGyfoPRbOAyOXPyOXBkqeOojL5b5582bYsGGyVgSQOlwut7a2dsiQIV1c8/r1a21t7T59+gh+xOFwampqhgwZItkt9MbGxvLy8pEjR2ppaUmwWQKw8F6FfBq5BFFyhwQAAAAoCkob1AAAAAAoFuCQAAAAALkAHBIAAAAgF4BDAgAAAOQCcEgAAACAXAAOSd7hcDgsFgsfyFcaqqur//nnH/HubWtrYzKZktUHkCFg4R3ozRYODkkcTE1NTf7F2tp6wYIFp0+f7kmDxcXFSUlJgvK8vLyJEyc6OTkJphJRXFpaWhYsWFBeXo4QmjFjxooVKzpcIFRIwOPxVq5cmZ2dLVUlezlg4T0BLFxslDZ1kFTh8XiWlpZz587lcDjv37/PycnZsmWLhoaG2DlF1qxZY21t7eHh0UEeHx+vqqp6+PDhUaNG9VhreSEhIQEXQEP/JgrpcAGPx+NyuZ3drqamtnr16rCwMDs7O6FVnIGeAxbeE8DCxQZmSGJCp9OXLl3q6+uLM4n17ds3IyOjwzWdHToW/TDy+/fvzc3NHRwcunhce96LKLdI6gD1hw8foqOjvby8etKIs7NzdXX11atXJaISIBSwcPEAC+8J4JAkgKamZr9+/YghT319fWBg4Oeffz5hwgRPT88nT55g+bNnz7777rvPP//c3Nz8f/7nf44fP47l/v7+L1++vH79+sKFC/mbXbRoUUlJSWFhobOz8+vXr52dnQ8fPjx16tS5c+fiNXehvSCEampq1qxZY2Nj4+DgkJiY6OLicuPGDTab7ezsnJWVha+prq52dnbGNUY7a+r9+/fOzs6nT5/28PAYP378F198QaSsb2pq+uWXXxwcHKZOnRoYGFhRUbFly5alS5cSOjx9+pRon5/jx4+3tbV9/vnnovxhk5KSnAVob2+nUqkODg7h4eGiNAL0HLBwsHBygCU7MXn16lVmZiZCqL6+/vr1601NTYsXL0YI8Xg8Pz+/8vLyr7/+etiwYfHx8YsXL7569aq2tvaSJUv69++/evVqFRWVU6dOBQcHT5w40czMbOHChQ8ePDAyMuI3d4SQr6/v77//3q9fv2+//XbYsGG45Jqzs3Nzc7ORkdGiRYsEe9HT02tpaVm5cmVNTc3KlSvr6upCQkKam5vr6urwTimxdfzx40cmk4lzfQpVWE9PD9+yZcuWefPmubm5JScn//LLL9bW1qNHj/7+++//+uuvr776ik6nR0VFMZlMd3f306dPP3jwABdRPnfu3IsXL3A9Zn5u3779n//8hz/R1qtXr5KTk/mvIQrNWVpaEhVRy8vLY2Njp0yZgguj2dnZnT17trq6WrAwKyARwMLBwmWA9PO3KiFECmeCr7/++u3btzwe78aNG3Q6/ciRI/jKN2/ejBs3bufOnSwWa+PGjfn5+Vh+//59Op1+8eJF/PaLL74ICAgQ7GjBggUrVqzg8XivX7+m0+lr167F8s564fF4ubm5dDr9xo0b+KOIiAjc0YsXL+h0ekpKCpZXVlZieRdNsdlsOp3+008/4Y8KCwvpdPqFCxcePXrE31Rubq6fn19zc7OlpeWWLVt4PF57ezuDwSBu5MfGxmbPnj3EWzs7O7owvv32W/67Xr58yWAw5syZQ2TLLikp4f+agGQBCwcLlwkwQxKTuXPn4gesoaGhuLj4t99+W7Ro0fXr18vKyhBCV65cIQp0UiiU8vJyopJjamrqs2fPiouLEULt7e3d6tTc3By/6KwXhFBRURFCaNKkSVg+bdq0iIiILtrsoqkOneKRWktLC1aeqK/MYDAYDAZCyMnJ6dq1a1u2bLl37x6bzXZ3d+/Q19u3b9+/fz98+HB+4fTp0/mL9yCEcF0fgvr6ej8/Py6X+8cffwwcOBALcamxsrKyadOmdfHtALEBCwcLJx9wSGKiqamJc/7r6up+9tlnL1++DA8Pf/LkCS50P3HiRCJFvK2t7dChQ1+9erVgwYLGxsaJEyfa2NjMnj17+/btvG7uow4aNAi/6KwX9O9PQG1trYaGBkJo8ODB/C0QNUybm5s/2RQGt4MQIuqc4srZeGGBH3d395SUlOzs7PT0dH19fVzQTLD3Dt+aQqF0UUG1tbX1u+++e/nyZUJCAn5EhX4dQOKAhYOFkw84JMnw9u1bhJCGhgYOFho5ciQu19jS0rJ7925dXd3Lly+/ffv25MmTVlZWCKFr164hhPhDP7sIAxWks14QQqNHj0YI5eTk4BDbmzdv4lvU1NQQQnV1dfhtRUXFJ5vqDCMjI4RQUVGRnZ0dQujevXurVq0KDw//z3/+Y2hoePHixdu3b/v6+gpGrA4bNmzIkCFPnz4V8WtyudyffvqpoKAgMjLSzMyM/yOsfwchID3AwsHCSQAckpiUlZXFxcUhhNra2kpLS1NTUy0tLUeNGqWnp6enp/fHH3+oqKhYWVnFxsZeuXLFxcUFD/1ycnKGDh1aUlISFBSEEGpqasKtqaurP378OC0tzcnJSZTeHRwchPaCEJo5c6aBgUFUVJSamlqfPn2io6PxLVpaWjQa7fjx45999lljY2NwcPAnm+qMadOmjRw5cvv27d99952+vv7evXt1dXWnTp1KoVDmz5+/f/9+/ELovWZmZqI/rnv37iX+JkTMsbW19eDBg3EjxGILIHHAwsHCyQcckjhQKJTCwsLCwkKEUJ8+fXR0dBYuXLh69WoKhaKurh4dHR0YGIifB2Nj482bN0+cOJHL5bq4uMTGxh46dGjAgAHff//90aNHHz58+M033yCEXF1d9+3bt3r16ps3b2pra/N3hEdh/P8ihDrrBSGkpqZ25MiRtWvXBgYGqqio2Nravnr1CiFEpVJ37NixZcuWZcuWqaurr1q1Cj9XXTTVoVNCJTU1tcOHD2/cuHHr1q1UKnXq1Kn+/v54UcLNzS08PHzKlCmCiw8Yc3PzxMTE9vZ2HFxEfMEOXWDh3bt3EUJpaWlpaWnEp3Fxcba2tkVFRfr6+sQKDyBZwMLBwmWDLCMqlJq6urrXr193ELa2tj5//pzL5Qpe39raWl9fL5FeMDU1NXV1dUSsERZyOJyqqqqPHz92q6nOePfu3bt37/glT58+pdPply9f7uyW58+fm5ubnz17tlsddaCurs7KyuqPP/7oSSNADwELFwpYeE+AGZK0GDhwIBEwQ6Cmpqavry/0ejU1NbwI3vNeMDo6Ogih9+/f8wtVVFRGjhzZ3aY6g38/uaam5u7du6dOnRo6dOisWbM6u0VfX9/b2/vo0aNubm5ip0U5ceKEtrZ2h0MtAMmAhQsFLLwnQKYGJYdCoVCpVMF4IYlTXV39008/PXnyJCQkhDjrJxR8nlEwD42ItLS0xMfHBwQE8J89BHotYOHKBIUnoQxOAFBTU6OlpSXKU/TmzZs+ffrQaDQxemlvb3/x4gWOgwIAMgELlzbgkAAAAAC5AJbsAAAAALkAHBIAAAAgF4BDAgAAAOQCcEgAAACAXAAOCQAAAJALwCEBAAAAcgE4JAAAAEAuAIcEAAAAyAXgkAAAAAC54P8B5LdatsZpeosAAAAASUVORK5CYII=\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",
    "    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, 1] )\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",
    "% 2. beat-tapping Lag\n",
    "subplot(2,2, 2); hold on\n",
    "x = nanmean(LAG,4);    \n",
    "    x = bsxfun(@plus,bsxfun(@minus, x, mean(x,2)), mean(x(:)));\n",
    "    errorbar(beats, mean(x), std(x)/sqrt(size(x,1)), ...\n",
    "    '.','color', 'k', 'MarkerSize',15, 'Linewidth', 1.5,'Capsize', 1)\n",
    "x = mean(nanmean(LAG(:,:,1,:)), 4);\n",
    "    x = polyfit(log(beats), x, fit_order);\n",
    "    freqCVmin = 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', 'k', 'Linewidth', 1.5);\n",
    "    h.Color(4) = 0.7;\n",
    "\n",
    "% + plot max fit\n",
    "[~,b] = max(x);\n",
    "    maxLAG = max(x);\n",
    "h = line([xVec(b), xVec(b)], [.0, max(x)]);\n",
    "    set(h, 'color', 'k', 'LineStyle', ':')\n",
    "hold off\n",
    "    set(gca, 'FontSize', 10, 'FontName', 'Arial')\n",
    "    ylabel( '|\\phi tapping-beat| (rad)' )        \n",
    "    xlabel( 'Beat frequency (Hz)' )\n",
    "    xlim( [.38, 4] )\n",
    "    xticks( round(beats([end:-2:1 1]),1) )\n",
    "    set(gca, 'XScale', 'log')\n",
    "    xlabel( 'Beat frequency (Hz)' )\n",
    "\tylim( [0, 2.7] )\n",
    "    yticks( [0 pi/4 pi/2 3*pi/4] )\n",
    "    yticklabels({'0','\\pi/4','\\pi/2', '3\\pi/4'})\n",
    "    set(gca,'Layer','top','Box','off','TickLength',[.01 .01])\n",
    "\n",
    "        \n",
    "% 3. mean: actual/expected (tapspeed)\n",
    "subplot(2,2, 3); hold on\n",
    "x = nanmean(TAPP,4);\n",
    "x = 1e3./x(:,:,1); % in Hz\n",
    "    x = bsxfun(@rdivide, x, beats)*100; % actual/expected\n",
    "    x = bsxfun(@plus,bsxfun(@minus, x, mean(x,2)), mean(x(:)));\n",
    "    errorbar(beats, mean(x), std(x)/sqrt(size(x,1)), ...\n",
    "    '.','color', 'k', 'MarkerSize',15, 'Linewidth', 1.5,'Capsize', 1)\n",
    "        hold on; plot(beats, 100*ones(1,length(beats)),'k:'); hold off\n",
    "    set(gca, 'FontSize', 10, 'FontName', 'Arial')\n",
    "    ylabel( 'Tapping/beat (%)' )        \n",
    "    xlabel( 'Beat frequency (Hz)' )\n",
    "    set(gca,'Layer','top','Box','off','TickLength',[.01 .01])\n",
    "    ylim( [88, 112] )\n",
    "        xticks( round(beats([end:-2:1 1]),1) )\n",
    "        set(gca, 'XScale', 'log')\n",
    "        xlim( [.38, 4] )\n",
    "   \n",
    "% 4. coefficient of variation (std/mean*100)\n",
    "subplot(2,2, 4); hold on\n",
    "x = nanmean(TAPP,4);    \n",
    "x = x(:,:,2)*100; % in Hz\n",
    "    x = bsxfun(@plus,bsxfun(@minus, x, mean(x,2)), mean(x(:)));\n",
    "    errorbar(beats, mean(x), std(x)/sqrt(size(x,1)), ...\n",
    "    '.','color', 'k', 'MarkerSize',15, 'Linewidth', 1.5,'Capsize', 1)\n",
    "x = mean(nanmean(TAPP(:,:,2,:)*100, 4)); % listen\n",
    "    x = polyfit(log(beats), x, fit_order);\n",
    "    freqCVmin = 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', 'k', 'Linewidth', 1.5);\n",
    "    h.Color(4) = 0.7;    \n",
    "\n",
    "% + plot max fit\n",
    "[~,b] = min(x);\n",
    "    minCV = min(x);\n",
    "h = line([xVec(b), xVec(b)], [.0, min(x)]);\n",
    "    set(h, 'color', 'k', 'LineStyle', ':')\n",
    "hold off\n",
    "    set(gca, 'FontSize', 10, 'FontName', 'Arial')\n",
    "    xlim( [.38, 4] )\n",
    "        xticks( round(beats([end:-2:1 1]),1) )\n",
    "        set(gca, 'XScale', 'log')\n",
    "        xlabel( 'Beat frequency (Hz)' )\n",
    "    ylim( [0, 20] )\n",
    "        ylabel( 'Coefficient of variation (%)' ) % Tapping precision fn beat(ms/Hz)\n",
    "    set(gca,'Layer','top','Box','off','TickLength',[.01 .01])            \n",
    "clear x h   \n",
    "\n",
    "% Print figure\n",
    "cd(origDir)\n",
    "% print('figure4b&S3a','-dpng','-r600')\n",
    "f.PaperSize = [46.81 33.1];\n",
    "f.PaperPosition(1:2) = [0 0];\n",
    "print -dpdf figure4b&S3a.pdf\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "min CV from indiv Fit:  mean = 1.3951 Hz, med = 1.3317 Hz, std = 0.7672 Hz \n",
      "min CV from group averaged data:  11.45 %, 1.73 Hz \n"
     ]
    }
   ],
   "source": [
    "fprintf('min CV from indiv Fit:  mean = %1.4f Hz, med = %1.4f Hz, std = %1.4f Hz \\n', mean(Fit(:,3)), quantile(Fit(:,3),[.5]), std(Fit(:,3)))\n",
    "fprintf('min CV from group averaged data:  %1.2f %%, %1.2f Hz \\n', mean(minCV), freqCVmin)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "% load tapping CV from auditory exp2\n",
    "A = load('../exp2/exp2_values.mat');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit beat frequency : p= 0.0371,  t= -2.1804,  df= 30.5832,  bf= 1.9270\n",
      "   df without Welch-Satterthwaite correction : 38\n"
     ]
    }
   ],
   "source": [
    "%% Stat compare audition vs vision on average (across beats)\n",
    "x = nanmean(A.TAPP(:,:,2,:)*100, 4); x(:,8) = [];\n",
    "y = nanmean(TAPP(:,:,2,:)*100, 4); y(:,8) = [];\n",
    "X = mean(x,2); Y = mean(y,2);\n",
    "[h,p,ci,stats] = ttest2(X,Y, 'Vartype', 'unequal'); % unpaired welsh ttest\n",
    "\n",
    "% unpaired Bayes Factor\n",
    "    T = stats.tstat;\n",
    "    df = numel(X)+numel(Y)-2;\n",
    "    N = numel(X)*numel(Y)/(numel(X)+numel(Y));\n",
    "    r = sqrt(2)/2; % Ã  voir !!! .7 or 1...\n",
    "numerator = (1+T.^2/df).^(-(df+1)/2);\n",
    "fun  = @(g) ( ((1+N.*g.*r.^2).^-0.5) .* (1+T.^2./((1+N.*g.*r.^2).*df)).^(-(df+1)/2) .* (2*pi).^(-1/2) .* g.^(-3/2).*exp(-1./(2*g))  );\n",
    "bf10 = numerator/integral(fun,0,inf); % Integrate over g\n",
    "bf10 = 1./bf10; % Return BF10\n",
    "% bf10 = 2*(1-pValue)*bf10; % if one-tail\n",
    "\n",
    "fprintf('Fit beat frequency : p= %6.4f,  t= %6.4f,  df= %6.4f,  bf= %6.4f\\n', p, stats.tstat, stats.df, bf10) \n",
    "fprintf('   df without Welch-Satterthwaite correction : %1.0f\\n',df);\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "min CV from auditory indiv Fit:  mean = 1.4179 Hz, med = 1.4072 Hz, std = 0.4439 Hz \n",
      "min CV from visual indiv Fit:    mean = 1.3951 Hz, med = 1.3317 Hz, std = 0.7672 Hz \n",
      "Fit beat frequency : p= 0.9092,  t= -0.1150,  df= 30.4387,  bf= 0.3104\n",
      "   df without Welch-Satterthwaite correction : 38\n"
     ]
    }
   ],
   "source": [
    "%% Stat compare beat frequency corresponding to minimal CV in audition vs vision\n",
    "fprintf('min CV from auditory indiv Fit:  mean = %1.4f Hz, med = %1.4f Hz, std = %1.4f Hz \\n', mean(A.Fit(:,3)), quantile(A.Fit(:,3),[.5]), std(A.Fit(:,3)))\n",
    "fprintf('min CV from visual indiv Fit:    mean = %1.4f Hz, med = %1.4f Hz, std = %1.4f Hz \\n', mean(Fit(:,3)), quantile(Fit(:,3),[.5]), std(Fit(:,3)))\n",
    "X = Fit(:,3); Y = A.Fit(:,3);\n",
    "[h,p,ci,stats] = ttest2(X,Y, 'VarType', 'unequal'); % unpaired welsh t test\n",
    "\n",
    "% unpaired Bayes Factor\n",
    "    T = stats.tstat;\n",
    "    df = numel(X)+numel(Y)-2;\n",
    "    N = numel(X)*numel(Y)/(numel(X)+numel(Y));\n",
    "    r = sqrt(2)/2; % Ã  voir !!! .7 or 1...\n",
    "numerator = (1+T.^2/df).^(-(df+1)/2);\n",
    "fun  = @(g) ( ((1+N.*g.*r.^2).^-0.5) .* (1+T.^2./((1+N.*g.*r.^2).*df)).^(-(df+1)/2) .* (2*pi).^(-1/2) .* g.^(-3/2).*exp(-1./(2*g))  );\n",
    "bf10 = numerator/integral(fun,0,inf);\n",
    "bf10 = 1./bf10; % Return BF10\n",
    "% bf10 = 2*(1-pValue)*bf10; % if one-tail\n",
    "\n",
    "fprintf('Fit beat frequency : p= %6.4f,  t= %6.4f,  df= %6.4f,  bf= %6.4f\\n', p, stats.tstat, stats.df, bf10)\n",
    "fprintf('   df without Welch-Satterthwaite correction : %1.0f\\n',length(Fit)+length(A.Fit)-2);\n",
    "clear x y  h p ci stats bf\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parametric repeated-measure Anova on Proportion Correct \n",
      "Beat                 F(  7.00,133.00) =   0.882,\tp = 0.4763\n",
      "\n",
      "3nd order fit: r2=0.78, p=0.007 \n"
     ]
    }
   ],
   "source": [
    "% Anova CV tapping in visual exp5 between conditions\n",
    "x = nanmean(TAPP(:,:,2,:)*100, 4);\n",
    "fprintf('Parametric repeated-measure Anova on %s \\n', varname{i0})\n",
    "[~,F,~,~] = repanova(x, 8, {'Beat'});\n",
    "\n",
    "% compute fit quality: À VÉRIFIER\n",
    "x = mean(x);\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",
    "clear x y r p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAB3RJTUUH4wsVDigYCLxOxgAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAyMS1Ob3YtMjAxOSAxNTo0MDoyNB313aUAACAASURBVHic7J15fFNV+v+fLDdpmqa9Tdu0TZekLC1L07CWOmKLVHHDMurMaBUQFcSvA18cxQXRH+JXeTkyKJsOMjgzUqEyOIoFRGYs0oAVUJY2sS20pQmlaZu26W1u0uy5vz8uXkOallKaNC3n/fKluSfnnnMS8H7ynPMsLIqiAIFAIBCIoYY91AtAIBAIBAIACRICgUAgQgQkSAgEAoEICZAgIRAIBCIkQIKEQCAQiJAACRICgUAgQgIkSAgEAoEICZAgIRAIBCIkQIKEQCAQiJAACRICgUAgQgIkSAgEAoEICZAgIRAIBCIkQIKEQCAQiJAACRICgUAgQgIkSAgEAoEICZAgIRAIBCIkQIKEQCAQiJAACRICgUAgQgIkSAgEAoEICZAgIRAIBCIkQIKEQCAQiJAACRICgUAgQgIkSAgEAoEICZAgIRAIBCIkQIKEQCAQiJAACRICgUAgQgIkSAgEAoEICZAgIRAIBCIkQIKEQCAQiJAACRICgUAgQoKRLEhbtmwZ6iUgEAgEor+wKIoa6jUEioyMjPPnzw/1KhAIBALRL7hDvYBeqa2tLS8vT01Nzc3N5XA4Pu96PJ4jR45cunRp2rRpCoWCxWINySIRCAQCMViE6Jbd9u3b586dW1pa+vLLL//pT3/yeZeiqKeffnrZsmXffvvt/PnzX3vttSFZJAKBQCAGkVC0kFpaWrZu3bp+/fqCgoKampp58+adPXt28uTJTIcjR44cO3bsww8/zM/P1+l0d911V35+/uzZs4dwzQgEAoG4QULRQqqqqmKz2ffeey8AjBs3LiMjQ6VSeXeoqanBcXzWrFkAIJPJMjIyNBrNkCwVgUAgEINFKArS5cuX4+LiuNwr1ltCQkJbW5t3h4kTJ5Ik2dnZCQB2u72xsVGr1fodCjnaIRAIxHAhFLfsOjs7w8LCmEuBQGCxWLw73HbbbTiOL126dN68ed99953VarXZbH6HWr58eWDXikAgEIhBIhgWks1mc7vd3i1ut7uurs7lcvntLxaLzWYzc0mSZEJCgncHDoezbds2sVi8b98+hUJx5513pqamBmLlCAQCgQgawRCkuXPn/ve///VuOX369H333dfY2Oi3f1JSksFgMJlMAEBRVH19fXJysncHg8Gg0+m2b9/+xRdf/OlPf9JoNEiQEAgEYrgTwC27w4cP7927FwDa2tq2bdv2+eef0+0URTU0NABAdHS03xtnzpyJ4/ibb775/PPP7969myRJ2oNuz549aWlp2dnZHA7n9ddf7+7unjdv3scff2y32+++++7AfZDecBKEVacDALKiAgDo1wxcHAcALCqKi+MCmQzDcYFMNijzkmo1qdGIMjNFCsWgDIhAIBChQAAzNVgslmXLllmt1srKytTUVBzHr0zJYvH5/NmzZy9cuLC3e48fP/7CCy8QBCGVSleuXHnfffcBQG5u7j333LNq1SoA2L1791//+te2traxY8e+9tprM2bM6DlIgDI1OAmCrKiw6nTOri4XQTDtlNPJwjBHe7ubJAVpaUwL04GL41hUlEAmEymV2C/fxgDQFxc3FxeLMjMTCwuRJiEQiBFDMFIHrVq16oEHHsjOzg70RD4MuiAZy8pMlZXeIsTF8cisLAAQyOXcqCgMxztKS0mNRr5ihZMgAKC7ro6FYd11dW6LxduEom8U5+Vd7xocBoN240ZSowGAmPx8+YoVg/DBEAgEIgQIUi47g8HQ1NREB7ceOnQoLy8vPDw80JMOoiB5SxGtJVwcj1Qqr2sQeouPtq7oFoFMdl0GE6lWX1i92rsFaRICgRgxBEOQTp8+/fTTT2dlZf3jH/8AgFmzZtlstrVr1951110BnXdQBMmq0xnLymgJGbBZ44OPMtHD9keW1IsXOwwGn8b0t99GG3cIBGIEEIw4pE8++UQqlT7//PPM5TvvvLNp06ZAC9IN4iQIQ0kJoxni3NzrNYl6A8NxDMcjlUonQRhVKrKioqO01FRZ2bcskWp1TzUCAFKjQYKECBEcDkd7e3tzczN9GRMTExMTIxKJhnDYnveKRCK73d7HaAH6FH2vatCnGI4EXJDcbvexY8fWrl2r+OWhKZPJXnzxxXvuuefChQvp6emBXsDAMFVUGFUqF0EMllXkFwzH4wsKxLm5ZEWFUaUyqlS0LPmdjh8f73cQvkQSiLUhhjtD8sg7f/58ZWWl1Wp1Op0dHR1ut3vcuHEzZszo/9RarXbt2rU+je3t7RaLxWq1MsGLUqlULBb3Z0D6Xu8WFotFUVQfo/W8RSwWX3P9a9askcvl0L9vnv6inE5nWFiYy+Vyu91TpkxR3PS/LAMuSBwOB8fx06dPFxQUMI1VVVUAwGaHYuIiAGjauZM2jMS5uQGSIm8wHBfn5YmUSkaWrDqdOC+vp5s4TyLpaSTxkCAh/BH8R55er29vbxeJRAKBoLa2lsVidXd3V1ZWtra2zpw5MyEhoT/qSKtRntf/dx0dHW1tbQaDwWw2u91uPp9PkqTdbhcKhdOmTet7SfS9breb8fK9fPlye3u7x+OJjY2tq6sjSXL06NHt7e18Pn/ixIk8Hs/pdMbHx3vfQhAE3T8xMbG3icrKytauXUufSvh883a7PS0tLSoqivnsTqeT/qJoAQMArVbb3t6u1+ulUukATLoRQzC27HJzc7/88svU1NScnBwej1dTU/Pee+/J5fLRo0cHYXZvrhnBQ++hWXU6Lo7HFxQMVuRQf2Bkid4nbNq500cOeRKJtLBQu2mT913I8xvhF0Yb/D7yAjGjw+Fobm42Go08Ho8gCKFQmJiY6HK59Ho9i8X66aefHA5HREREf9QxLy9v0aJFzLBqtfrs2bMpKSkej4devMvlqqioGDNmzO233y6TyXqzSOh7tVqtVCoVCoUA4HQ6q6ury8vLo6Oj4+LipFLpqVOnMjIyGhoaYmJioqKiRCJRREQERVHh4eHd3d30gOHh4SaTacKECWKxuKOjw2cWZtllZWXg75tXqVSlpaUSiSQ+Pp7Wp4iICC6X6/0HkZSUpNVqOzo6YmNje/6SEIlEtNFpNptxHB/B5lQwBGn16tWdnZ3vvvsu0zJ69Oj3338/+FX1SI2mubjYnJkJAD2f47QMAIBAJkvqPUaqbxi374HdjuG4pKCgN1MpJj9fpFB4u31LCwsHNhFiBMNoQ2+PPB6PF4h57XZ7V1eXUCg0mUxxcXEAwOVyu7u7MQwzGo0xMTEKhQLDMLhOdTSbzXq9XiQS0WPSw4aFhTU1NZlMJofD0YctaLfbLRYLrUY0DoejpaWFw+GMHz++trbWarUy4pGSkuJyuQwGg16vj4qKcrlcfD6fHpDFYjkcDj6f37fF2fObb2tr4/F4FEWlpKTIZDIMw+rq6pqamlgs1tixY5kbMQwzm81Go7G9vd1Hz86dO6fVamNiYhISEpYuXbp///5A/7YYQoIhSDweb/PmzXq9vrq62mQypaamTpo0qWcR2CAgysw0Z2aSGk1ED0cAq07XWlICwdqm6wPGVKL9HVpLSrxPlXgSSYRCQWo0WpA7JNnyIVwoIgR44oknejZ+9NFHPR/EzCOPViO/N/YTemOqJ2w2Ozw8vKury2q17ty5k240mUwul6u7uzs6OvrIkSN0o9PpdDgcIpGo577T0aNH867+H5DD4XC53K6uLiZDmM1mi4yM7Orqamtri4mJ6c0WjI2NZbPZQqHQ+6twOp3h4eFOp5PFYhEEwWazeTxeeHi42Wzu7OxUKBTnzp0zmUzR0dHjx49nBmxubuZyuWPHjmUUqKcqHD16dNGiRZ2dnQ6Hg9mgIwjCarU6HI7y8nL6mMpqtXZ0dHA4HObbAIAXX3wxIiKCz+dfunTJW8+cTiebzW5paUlISGCO1oLw22KoCFK276qqqh07dtTV1d17770ej6exsfG3v/1tcKb2RqRQkBoN2aN4EmMbSQoKbtCVLiY/PyY//0ZGoKH9HbCoqI7SUtrZIWnBAtoBT5SZqcmcW6aJuQVkt9z4TIjhTF4vP556PogtFktERERERITD4eDxeL3dOGB4PN7YsWNJktRqtQKBIC0tLSwszGg0yuVyk8lkNBrHjh0bGxuL/ZK7pLGxUSwW93Rr8iklw+PxxowZU1lZqdVqrVarQCAAgM7OToFAEBERER4e7vMEh6uf12PHjrVarXq9njZHMAwLDw+XSCR2u93pdGIYxufzHQ5HWFhYeHh4RESEx+Oh/R184mFsNpvNZouNjfU7C60Kcrl8zpw5Op3OaDTSawOAy5cv0580MjKS+fgNDQ20HnsbUvHx8QkJCXq93ueXhMfjIQjC4XDMnTt3+/bt0OO3xUgiGIJ0/PjxJUuW8Pl8DodDkqTb7d68ebNarX799deDMPs1YdQoaeHCQB8aXW8aOtpUaioqchFEU1ERbSqJFIoOTYxW04zUCMGctfjg8yCGXx55iYmJ9FOstxtvBJFIlJCQEBMTo9VqhUJhdHR0QkJCWlpaY2Oj3W6fNGlSSkoK3dNisej1erlcrlAofJ6q9EmMN2KxOCsrq6Wl5fTp00lJSXa7ncvlCoXCiRMn+n2Cez+vRSJRbGysXq+vrKzkcrk2m43FYk2aNKmuru706dMGg6G9vd1kMolEovT09Li4OKfTGRERYTabXS5XXV0dh8Ox2+0sFsvj8Xi7OYA/VZDL5YsWLero6KisrKQFmD6yqqqqoo/M6S1Hi8Wi0+lYLJbVauVyufSqSJJMS0uTSqUtLS0+vyQcDgc9b1ZWFi1gPr8tBvcPcWgJhiC9++67Y8aM+de//rVs2TIWi/Xss886nc5du3a9+uqrwdy4I9Vq7aZNtJdaR2lpbH4+TyIJtBoRhLOry6XVWgGgq8t5ubzFXGUIS6mNmIBxRSIAwHEuAERFYTjOxXFMJhP4jIDheNKCBcypUlWV+d1DV4pxGAyOQV8wYmTQ80GMYRj9yAvcpDweTy6XR0ZGnjhxora2FgCEQiFJkkKhUCwWexct81HHaw47YcIEp9P5/fff22w2DofD4XBEIlFUVJTfJ7jP81oul8fExHi7PISFhaWkpHz//feNjY1sNttkMsXHx0dERMTFxdXV1cXFxeE4LpPJGKcG2k2cx+N1dXUxe3G9qYLPN9/V1cVisSwWC3MA1tbWlpycHB8fLxQKvVcllUppK7OnSZecnOx2u2Uy2UcffXS9397wIuCZGlwul0KhWLt27R/+8Iennnpq/PjxK1eu1Ov1t99++969e7OysgI3tXemhp5Jd3gSSerTT7cdPgyDqkaMAul0Vp3OiuNcgvi17JPN2OU2NLeSvGhpdHSS/2TnOE5nxeMqlZHe+mSqqDj713/p9fYGSDsBOXRjZqZo3boQjeVCDC0kSZIk2fORF/ypIyMjSZKsrq6mKIpRx6ysLL9+YmvXrv3nP//JHAgx2O12AHA4fv0RxuPx+Hy+3W7v7u7mcDjMR3M4HG63Ozw8nM/n97ZCejStVut0OoVCIYZhGIZRFEVRFJ/PZ7FYPgPSahoWFtbbLFqtdtGiRWvWrPH5+E6nkyCI1tZW+iSs789OD0v7aDDfFYvFor3s+vPtDXcCbiFxuVwcx+mfSwzl5eUAwPzWCALNxcU+LQ6D4dK2bYK0tEFRI4JwVlSQfkUIx7kymSAqCgOAcAeY9x3YT2ZMnhifW5BEEE66j05nBQCCcOl0VoJwEYRLp4OKCpK2n7KyIuVygUyp3KbvyofSNGiQgKEU8i0g1GhItZpUKEZgRALiBqFdBobEEavn1HSZTR919Hvv448/fr3nW1qtlpYWLpfrcrkwDJPL5T0lzS8Oh4MkSR9nbr8Dwi8C1tsss2bN8vvx/f4y6G09fk26nnFII8+/jiYYuezWrFmzb9++hx9+uLy8PCMjQyaTffbZZ9HR0QcOHAio57e3heQ3CxxXJEpft+4G1aiszKjTWbu6nIwI4Tj3ioRcvf9Gl42gX9NBRT3dH2iJooWN1ifmLZJ01dRYhGBJg4sK0FhASGtSYWFiYeHI/NuJQPSHQbcF/Q5ot9uHyuK8eQjGGdLLL7/c0tLyySefAABtKsnl8iDHIfH95TiIys4esBrRJpFKZWRaaB3KdJ0lNRp5nm8ckrcaAYDDYNAXF4sUCp88CziO0f9WKiPpWWhxqqggeTw2AFhA2ACjAEABmnwovQhpEol8YB8BgRgZDLot6HdA2kVipJomIUIwBKm2tvYvf/lLY2NjbW1td3e3TCabPn065lW5LgjE5Of39PYWZWYOYCgfKaJ1KC/vSiKsjlL/d/ndM2wvLe07shXHMVqccnPFXV2ud9+92Nhos4BQAwoAUIBmFDRM8KgBZg3ggyAQCERIEYwtu1mzZk2YMOHDDz8M9EQ+eG/Z0RZJR+mvciHKzExft+66BvSRIplMoFSKaFPmmpz2SuXHcL3VjAwGx6uvnmec6yRgmBdxVCzG5HPzk+6+vf/jIBAIRAgScAvJ4/HYbDYuN0gRuL3Bk0jkK1aIMjPp1D6JhYXXm3SnrMxYWWmiD4pkMkFBgYTeXrvBVV1Xf4mEt25dhlpNlpZ2aDRk7IRRVbyEjKaDpt2HT5zonPbEPT29xv1yveFQCAQCEQQCrhNsNvvhhx/etm3bwoULFQqFd2TZ/Pnz6bjroBGTn283GHru3fUNQThLSgy0fwGOcxcsSBqAFB2FvFngG/FXAZOud0NaIuHl58cYDHaNhlQqRffcM6rmtLhld1H4xVPfvu/Cptyalye+pizROf0AJWZFIBChRDAMl0OHDgHAyZMnT5486d3+4IMPBlmQBkBZmZHeo8Nxbm6uuJ8bdD2pgEkAQGvSVzAvCggdyEcbBPdc/1BqNalWm+nXOI7l5Musyj9Wb/07XDzTforcqctTKkW5ueLeVNNhMJjVavrFwD4LAoFABIJgCNKuXbvonBw+7UH2awAAUq2mn8X9wdswys0VM24LA4MAvAxmVcCkPCgDgAqYRAB+i2QgPqMaDanRkJmZoszMK+FHgthoxXNPRf79E/ElLbRCRUWeTmf1drVg6CgtZQpYdJSWkmp1xrp1qKISAoEIBYJRIu/3v//9s88+i/UgCFP7QGdWFWVmXtO/Tqezbtmio6NcFy5MukE1AoAVK+QAgAOBAzEPvlLCOQDIz4+9xm3+yMwUFRYmFhYmesfDYjgue/Lx5Cljp4mbpjp/IAiXSmXcskVbUWFi+tDJk7yHon09BvaJEAgEYnAJXaeG2tra8vLy1NTU3NxcvynvdDqdSqWaNm0akyX+mogyM6Gw8JqH+cw23QAMo97qISkUosLCRH3xUTloASBewn97RbpkQBaSQiHym5qBrqUEJSXpOu0YWfgJyNHprLSRR+/gdXilu2cg+20yIhAIREAJUaeG7du3b9iwYcaMGTU1NTk5OZs3b/bpsH79+p07d06fPn3jxo2/+c1vtmzZ0p/FiBSKax7j79zZRBtGBQXx/XRa6w8SCa+wUPoD/OZoMQDAb/KnBiLfD61JrV98YdNV3aXEamQzVSojndMoKyvSbwyww2BwGAxo1w6BQAw5wYhDomuE9GwvLy+PiYnp2d7S0jJnzpy33nqroKCgpqZm3rx5n3322eTJk5kOJpMpNzf39ddff+ihh86dO/fwww8fOnRo1KhRPuN4xyHRqNVXTl/8igFzaDRgV7prQi8AAHpbw6DgJAi6DrpIqeTn3sPETtka6u9q/6cQLD79U3Z8PjBbDYFAIAaRIHnZeTyenu29HSNVVVWx2ex7770XAMaNG5eRkaFSqbwFCbyS/tK5pIxGY09B6olGQxYXNxcWQk8xIAhnUVETQbhkMsHChUn9+FgDobfdtsGFtpMMJSVkRQUA5BUUyOWCkpLWC7z4b+DudDivgF8d3z+BxzcjNUIgECFAMASJw+GcOXNm27ZtFy5csFqt6enpDz744O9+97ve+l++fDkuLo45dkpISGhra/PuEBkZ+fzzz7/++usHDx6srKycM2eOj1z1RmamKDPTrFabMzOvypCt01l37myCwfCmCxFoTWoqKiIrKrCoKFle3oIFST/+GLFti80FXAtE5MAJAvBzoBydP3WoF4tAIBAAwfGyO3r06Pz588vLy5OTk5VKpUajWb169Y4dO3rr39nZGRYWxlwKBAKL5apdJovF8vXXX0dHR0ul0ujoaL1e39ra6ncon7MlhUKkUETQbtNMI6NGBQWS4KjRpk3a0tKOQM9CV/bj4rhRpTKWleE4duedsc+uGFsD436CaSVQ8D3c2pF5P+3+h0AgEENOMCyk9957b/To0bt376Y9Gmw22+rVqz/88MMnnnjCr/ucWCw2m83MJUmSGRkZ3h2OHj1aXV1dWloqlUqdTmdeXt6ePXv+9Kc/9Rxq+fLl3pcGg4MJKaVh1GjhwqRBdGEIEWhNaioqMqpUACDOy8vPj1EoRN8Xn/i+VBg/QZZ310D8zhEIBCIQBNxCokvTz58/n/GvCwsLW7JkicViqa6u9ntLUlKSwWAwmUwAQFFUfX19cnKyd4eLFy9KJBI6DzyGYenp6ZcuXbrmSoqL9YsXq2nbqLS0o7hYP1RqtGKFPD/fjzdHIPCxkwBAIuElSDAbhDUTGB2rxNQJRCAQiCEk4ILE5XKFQqGPtxt9GR4e7veWmTNn4jj+5ptv6vX6DRs2kCQ5e/ZsANizZ8+pU6cAICcnp6WlZdeuXTab7T//+c+pU6f+8Ic/9L2M0tKO4uJm5tJgcBQXN3/wgQ4GW428UyGECLQmAYCpspLWJC3ILkHqWGU8Xdy2qKjJuxIgAoFADAmcN954I9Bz1NXVHThwwGaz8fl8s9n85Zdfbt++PSYmZtmyZX5r9HE4nIyMjL/97W8fffRRS0vLa6+9NmXKFAB49tlnORzObbfdlpSUFBkZuWHDhq1btx4+fPiRRx6ZP39+z3G2bt3KbNmtW1dvsbh9Olgs7pUrRw2ubWRtaHAYDHhOziCOeeNwwsJESmWnSuXs6uInJHCj8Ph43rRpUbfdJg4L45w/b9HprDabRy4faZuWCARiGBGMOKT29vbnnnvuxx9/ZFqkUumWLVsyB1Qfr/94xyEtXqxmyggxSCS8HTtuonTXpooKQ0kJF8fjCwq8S+WWlRlLSzswjNV/J8O+I7oQCARiAATDqSE2Nnbnzp0//fRTbW2tzWaTyWQzZ8709qMLAhIJv6cg3WwP00ilEgAMJSWtJSVJCxZgv5zq5eWJlUpRUVGTSmWsrDT1Jyi4j4guBAKBGBgBP0NyOBwmk4nNZmdnZz/22GMzZsxQKpVBViMAKCxM7Nk4e3aQPAtCh0ilUpyb6yKIpqIiJ0Ew7TiOLViQpFSK6CMl5OaAQCCCT2AF6cCBA3PmzNmzZw/T8v777+fl5a1fvz4IW4Xe0LlNmUuJhLdihXyoft0HJw6pN8R5eYwmebfjOJabK87NFdOaVFZm7G0Exnu+p9GJQCAQAyaAgvT999+vXLnSZDKlpKQwjQ8++GBcXNyOHTu2bt0auKn9Ulgo3bFDQftb5+fHBM3xOgQRKZUipdJFEK0lJd7tOI7l5V3RJJXK6FeTSks7GO95tZosLtYHadEIBGKkE0BB2rZtW2xs7Ndff3333Xczjffdd9/hw4cnT5788ccfW63BdjWWSHihkEU0mHFIfsFwXJybK5DJyIoKH00CgLw88fLlMhzn9tQktZrctEnLXNLe80No7SEQiJFEAAWppqYmPz8/ISHBpz0sLOyBBx6wWq1arTZws/vFu/h3gAjBOCS/0MnuuDhOVlTQwUne0EdKtCZt2aJl2o8c8aM9yEhCIBCDQgAFic1me1c/8mbMmDGBm7cPehb/vplhAmaZJA7e0JokkwkIwsVkc2ht9XNohE6SEAjEoBBAQZo+ffqxY8f8vnXixInIyEifDHVBwG/x78ElJj+/Z7nYkAXDcdny5QBgqqy09qhZheNYQYGEcXMgCKdCETEUy0QgEDcFARSkW2+9taqqavPmzW73VSkS/vOf/+zYsWPatGlsdjByjXujUIgKC6UodMYbeu+OdnDwdgSn8XZzKCpqGj1a2HMElC8cEbo4HKDXw+nTV/7RaoEkr30XYogIYKYGj8fz+uuvf/755+PGjZs+fXpycrJer79w4cKJEyekUukHH3wwfvz4AE1N07NiLKI3jGVlRpWKi+Pyq/OjM5SVGVUqI45z29udTU02ZpuusDCxsFAaxJUiAozDAe3t0PxL4seYGIiJAZG/33AOBzQ3Q3U1tLcDAERFwejREB8PdjtcugQEAe3tEBUFOA4pKVcGMZvh4kWorLwywqhRV24JHGq1s7ra7XaDQAAuF8fhwBQKUNxE+VmGF4FNHURR1AcffHD48OG6ujq6aGxiYuItt9yyatWqyMjIwM1LE7KCtGmTNjNTFGp+57QmCWSypIUL/XaoqDCVlBhwnJuQEKbX20pLO5AajUD6eIL7KFBn53e7dtVWV3djmJvF4rndeqtVZ7MlhIVls9lCHo/N4XR6PPVOZ73Ldc5mUwP8FuA2kYiHYd1stjgi4qG7746eMwceeODK1D5aGBkJFPWrQdOHNPaGXu+uq3PY7VRaGt3Aamjg8fmcMWNAiv7ehiKBTR3EYrGWLVu2bNkyi8VSV1eXkpIiFo+EeqwjEpFSadXprDpda0lJfEFBzw5KZaRMJtiyRUcQ5p6ZahEjAZ3OU1/vaW+nn9dUbKy7o4NdW8tpaoK4OKirc1VW2k0mT3g45XCwdTrs4sUHH3yQd8stAEB1d2NGI6bV2o1GB493Tq3OVCq5JlNYSgpn8mTu2LEcLtdZU9Ntt7uzssBs/vuGDbqOjohTp7CICMjJAZEIzp+/SgsvXmRxOC6p9Mrl+fPXZ9w4HNDc7DKZKK84SCopyXXxIqejA2JjgTf0ESAIH4KRyw4AhEKhUqkMzlyhT2geuvhUPRfn5fXsg+PY8uWyOjZebwAAIABJREFUoqKmpiZb8FeICCwOB5w4YauocIvFpZ99Njs3l6PVcrhcW3MzFhVF4Th18aKrtdU9YYJ78mTQ6ai2NoFIBCTpamjwCIUcs9l96ZK7vZ0KD6eioy18vismxsXhUAZDmFrNtlo5LpfdZvNMmgSXL3vq60cThEAqtVVWsjs7OXw+jBnj7ux0RUdfsWba2916PcvjYclkJSdPTp06NcluZ3d2cvT6/hs3HqeTslhA6HXwyeNRZrO7q4uD1CgkCbZbwYhnuMQh+cW7cpKposJvH9odPCVFIJXyPZ7grg8RUJqbXV1drrAw57RpX9bUtMTFuXg8Z2cnZbG4kpIosdiDYY7oaBePBwRBdXW5PB5TWJhVKPRERbnS0lw47m5v97S3ezCMslpNAoEtLc2mUFjDwx1Op9tgcDY2Ug4HxWKxL16ElhYXn+8QiZxstq2lxXPhAlRWukwmSiIBAHA6ob3dKRC4MYzV2XnmxInTp09TSUkukwk6OsDR3zADFovFEgrBYvm1yWJhRURwBIL+D4IIJkGykBDDBQzHkxYubNq506hSYTjuXaWCAcexP/4xtaio6eLF7rIyYz8rViBCGnqDq7ubiosDqxUAXvzTn+QUNZGiumy2yh07eGz2VB6vWyDo4vG67PZUFmtUeLiAxeKQpMfppLhcd1QU2+Phu91gtzswzIlhAMCiKPB43CyWi8vlURQLgFNZ6ezqcsbHt8bGXjhxItlqdXZ1tZSWtlmthEDwI4cDABjAeDZbyudz3O4Wp/M0RU3Jyblu44bHY40Zw7XbqcuXPb9EmLAMBkwshsREtF8XmiBBGmRi8vNj8vOHehU3hEAmE+fmGlUqnyoV3tB2El2xoqvLWVAQSEcpRFDwAIBQyMZxT2srAKxfvz7GYuEcOcKbOJF1662s1la3Wm3j8dweD1siYXk8znPnfj5xguV2s1ksj8sFNhuXxWIJhW4ej8VihXk8AMA2GjlhYRw+n4NhlFLJt1rtJ07YIyNdcjkAzJ0yJTU5mZWWxq2rE4SFuZRK9+jRIBSC00lduODS69keD5aUtKa4GOBq46afciIScaKjPS0tcO4ccLlgtXI4HHZCAvJoCFmQICH8QB8gGVWqpqKi3hzBGU2qqCABAGnScIfFYrHi4jgU5bFaH8/JwS5d4nR0sNlsLoa5JBIQCtkkidXUAIvF4vMpp5PD5/MsFszhYF26xGlr41it3ISEMJK0mkwckSiuqyusooJtMvEiIjg4zpZKOZmZHK3WZbdzGxo8FDW+vV0kErEjI93JydTFiyAWc6OirlgzGMZKTubqdBynk5JI5j74IAzMuOHxQC7HYmIwby92pEYhTDAEiaKo77//nq7O592+aNEigQDVzA5Rrul0B79oUkmJAWnSsIfZ4DIYsMmTc1JSAIBVW4tFRrolEuDxKAB2XBx25gzbbPY4ncBisY3GLqvVqVBE0Z3NZn5KCttoxKqrWQ0N4RaLwGRiCwRYejpr3DgsJQVGjYK0tDC7nXX0KMbhNHk8iePHuydNYjU1cePiOKNGQXS0p67uijVjNnNcLhaAS6+fzueD1crp7ByIcSMSgUiERGi4EIwS5lu3bt2yZUvP9pMnT/aW7G5QQHFIN4iTIJqKilwEIc7N9et0R0MQzpISg05nVSpFSJOGMQ7HFcdrNvvKBhdFYWFhTpeLbmF1drK7urCwMIiIAJcLSPLzL7746scfbSxWGEW53O6fbTYDQBLAXQCJYWEUl9sM0E5RZrf7Z5tNDQAACQB3c7njIyOt8fELn3wSFwi4bW3CMWPYs2eDSAQk+WscUlgYCIXQ8Us+X9q4QWc/I5pgCNKtt946atSod955JyEhgcViMe2BTh2EBOnGcRKEbssWAJAUFET27riPNGmEQJJXSQIdiGq3g8+WF60KPTv35y2HA86edarVboIALhfsdk5UFJadDVOmBOtDIkKXgAuS2+1WKpVr16596KGHAjpRT0JWkIYXVp2uaedOLo7HFxT4dbqjQZqE6C99yBXi5ibgcUgcDkehUKhUqkBPFCIM6zgkvwhksj6yrzLQqcFlMkFFBVlS0hrMFSKGGfShztSpV/6Ry5EaIWiC4dQwZ86cd95556mnnpo4caK3F0PfTg21tbXl5eWpqam5ubkcDsf7La1WW11d7d0ikUimTp066CtH0EQqlS6C8Ha6I9VqUqMRZWaKvFK50JqEfBwQCMTACMYZ0pw5c3Q9au0AQHl5eUyM/3OU7du3b9iwYcaMGTU1NTk5OZs3b/Z+d/fu3Vu3bqVfUxRlNBrvuOOODz74wGcQtGU3iDgJwlBSYtXpREplfEGBvri4ubg4sbBQWljo0xPt3SEQiIERDEHy9JJhpjenhpaWljlz5rz11lsFBQU1NTXz5s377LPPJk+e7Lfz3r17165du3fv3p7FLJAgDS7eTnc2vb43QQKkSQgEYkAEI5cdm80+ffr0kiVL8vLysrOz58+f/8UXX/ThYldVVcVms++9914AGDduXEZGRm9HUM3Nze+8887y5csDXVoJAVeXPCd7yXRHg86TEAjEAAiGIB09enT+/Pnl5eXJyclKpVKj0axevXrHjh299b98+XJcXByXe+V8KyEhoa2tzW/PDRs2xMXFLV68OCDrDhibNmlLSzuu3S/0oDPdAUAf3g00SJMQCMT1EgxBeu+990aPHv3DDz8UFxd//PHHp06dmjt37ocffuhT2pyhs7MzLCyMuRQIBBbvfL2/QBDE4cOHH3vsMR+XB2/8BuQibgQ6011/etKahOPcigqyrMwY6IUhEIjhTsAFyeVy1dXVzZ8/n0nKEBYWtmTJEovF4uMpxyAWi81mM3NJkmRCQkLPbl9++SWXy32AKTfpj+W95GEbWlaskA+LqFi/OAyG9sOH7Xo9AHSUlurpxJe9QOcWwnGuSmX00SS1miwu1qvVZG/3IhCIm42ACxKXyxUKhT7OBfRleHi431uSkpIMBoPJZAIAiqLq6+uTk5N7djt8+PDdd98dERERgFUPnJEXh+SDevFiUqOhXzsMhubiYlKt7qN/b5qk0ZDFxc0aDRIkBAJxhWBs2d1+++1ffvnl+++/f+bMmdra2h07drz//vsymUwul/vtP3PmTBzH33zzTb1ev2HDBpIkZ8+eDQB79uw5deoU062+vn7ixIlBWD+Cwa/WXlOA/WpSZqYoM1OkVpuRkYRAIGiCIUgvvfRSVlbWtm3bCgsL586du379ehaL9d577/XmaMfj8f785z8fO3bs9ttvP3jw4P/93/8lJiYCwAcffFBaWkr3aW9vN5lM6enpQVj/dRGTny9fsWKoVxFseisvy9BTkxQKkUIRodGQyEhCIBA0wcjUEBsbu3Pnzp9++omuQCGTyWbOnOntttCTmTNnnjx50qfR2/k7NjYWxRgFH0erH385h8HQR3lZBu+afgAwfnyEWm3uoz8CgbjZCGBg7GeffeZyuebPn//pp5/SB0I+PPHEEwGth4QCYwcdUq2+sHq1T2P46NGYWMzF8d7Ky3pDEM6//KWhocFKki66RSLhrVuXIZGgbGYIxM1OAAXpt7/9rdls/vbbb++66y6tVtuzQx+pgwaFkBWkYVR+oid00iDmUpSZmb5unbGszKhScXG8t/Ky3rz66gWfbTqJhLdjh6K3/ggE4iYhgFt2+/bto1989dVXXC7XO1qIxWLZ7fY+4ocQIYu0sDA2P1+7cSOdXDV93TroX3lZGrXaz6GRweBQq0mFQhTAdSMQiJAnGE4Nc+fO/fbbb1lenDp1Kisrq7GxMQizhyDDOg4JAHgSSYRCAQARv6T6xnBcUlDAxXGyosJYVtbHvQaD47raEQjEzUMALaTDhw/v3bsXANra2rZt2/b555/T7RRFNTQ0AEB0dHTgZh8qOkpLSY3mJnS0ozPdNRUVGVUqAOit5HlvZ0XIPEIgEAG0kGbOnOl2u81ms9PptNls5l/o7u5OTU1dvXo1fq0DcMTwgsm+aqqstPorOAIACoUoM9NXe0QiLo/H8tsfgUDcPASj/MTzzz+/dOnSjIyMQE/kQ8g6NYwA+qiHZKqoMJSU9OF0ZzA4Nm7UMidJEglPJhPgOHfBgiQcxwK+dAQCEaoE4wzpzJkzm0Z0Np2bEFFmZmJhoSgzs+dbkUqlODfXRRBNRUV+k4JLJLx169J37FDQB2n5+TG5uWKCcBUVNQV83QgEIoQJuCB5PB6bzcbUkkCMDEQKhbSw0Lt+uTfivDxGk3obQSLhMedJSqWI1qSdO5EmIRA3LwHXCTab/fDDD2/btm3hwoUKhcL73Gj+/PkBDYwNWYZ1HFI/uaYjuFpNMpkacBxTKkU6nVWns5aUtKIiswjEzUkwDJdDhw4BwMmTJ32yAT344IM3pyDdDNCO4E1FRWRFBRYV1dPpjs5iR6dYhV+KJ5WUGCoqyKgoLC9PPBSrRiAQQ0kwnBrcbrfH4+nZjmGBPcFGTg1DjpMgmoqKXAQhzs310SQ6QjYzU+Tt8E0QzqKiJoJw5eaKkSYhEDcbwRCkoWJIBOmmjUPqDSdB6LZs4eJ4fEFB39lXaZAmIRA3LcHwsgOAH3/88amnnrrtttumTZv26KOPMkGyiBEPvXfnIojWkpLegpO8oZOCA0BlpQkVPkcgbiqCYSEdPXp06dKlbDZ70qRJ4eHhP/74o91uf/HFFxcvXhzQedGWXehwzeAkHwjCuWWLDse5BQXxMhk6aEQgbgqCIUgFBQUul2v37t20i53NZlu9evV33333448/BjS/KhKkkILJCN5PTaqoMJWUGJAmXRcOB7S3Q3PzlcuYGIiJAZG/rEz97zkkOBzQ3AzV1dDeDgAQFQWjR0NKSgitEBEIAi5ILpcrKyvrtddee/TRR5nGmpqaefPm/fvf/870F1k5WCBBCjWuq0oFeGkSSuLQT9RqqK52ut1ugQBcLnA4OAoF5jdarLeefoWKz/fTePr0TxcvQn29kG4cPRri4y1CYX+KLvI7OqCjgweAdXUBABYe7sQwZ0wMiEQOHs8OAFqtqLISb2/nsNlOqTTeauUmJvLuuYfb2wqRUI0MAu72zeVyhUKhjzDQl+Hh4YGePTS5GeKQ/CJSKgHAqFJpt2zpjyYplZEAUFJiKCpqQpp0TfR66Ox0R0e70tKu/MpsaKA6O9l6PUcq7W/PlhY4c8bZ0eHm88HlAoGAM20aFh5+lXqdP89xOCpffHF7ZORtPB6fx3NQFGY2Q1hYRXz8D9dcJ0nKHY6Uzk7C44kTicZ6POywsBY+n7DbgcdrEIm0dnuMxTKOw0kSiUxdXVWjRj04fvystjZWfT07Jobd0eG7mN5E1y9Iz0KZYMQh3X777V9++SWO43l5eSKRqKys7NNPP5XJZHK5PAizI0IHDMf7XzmJRqmMJAiXSmUsKmpavlwe+DUOV+g9LpPJlZLy655HUhJ18aKro4MTGws83rV71tRwGhrcly/b4uJctCA1NWEnTrBTUsBHvfbvpzIzf3v//bdOmeKmG8+c4fD5D952G5aV1dciq6vh1CmP0+msqTl79uz5xx+/VyBw83jclBQOSbLYbJ5IxKmqgpoaO5drT093Hzv2RXt72y23uC0WV1MTt7KSjWH9Et3eOH++v3oWCOlCctg3wRCkl156Sa/Xb9u2bdu2bXSLVCp977332Owg+fiFGitWyId6CUMG7XRnKCkhKyoAoD+aRDt/q1TGLVu0SJP6wOn0WCyUUAgAsHTp0jVr1kilUrOZ6upy83ic3nrS8HjQ1UXV1LgpypGS4ho9mvql3anRcCwWzp13/tpZLKYaGthOZ3Ry8q/xhRMmuE+dstfXY6NGQUSEn+U5HFBWBl995Wpvd/B41MWLYrN5cmsrW6Fwd3S4u7o4CQnU0aNuo9HT0uJua3NERDhcLhZBiAQCE4aB0wkGgzssjBMff23R7Q1v0/D+++/fv39/H3rWf+nqP4EYcyQRDEGKjY3duXPnTz/9VFtba7PZZDLZzJkzw8LCgjB18EFxSNfkmkkceoISC/UHFoslFLIsFhAKQa/XHzhwIDZW1twcdv68q6qqw9tCamwUm82YTmfDsCt1EZ1OXmOj6NIlMY/HffLJZGZMqZQqL3ezWBAR8asG8HhAklaxOEEi+bUxLAzMZqqjwxkR4X9n9fvvobzc3t5ut9sdPB5FUeB2YxcvusPDuXy+myQ9IhHn4kW32eyKjqZsNhef7zQaWRZLhMHQWVp6ork5giS5Op0zJoadmqr3HlmrlURGWkePNvScdNGiRcxrb9NQr78yQm961v/9z/4TiDFHGEHKeVpRUbF3797q6mq32z1x4kQMw/L68RhCjFT6Wc2PwTuxEAAgTeoJjwdjxrDsdu7ly1RGhgcA9Hr9zz9bMYwdHt6m1XZc3T2mqyvx/HkQi1vp687OaAwTGo0pTiebolIArjwxHQ6WxwNuN7u72xMVdaXRYoHmZr3Nxm5vF8fGUr+MwIqIYIWFYWazHwtJp4P6erfdbk9Lc9jt9uRkd3d3NwBps0kuX+YlJrI5HPbZsyybzZGaSmVmUvX17KYmLovlcLn4Z88ajEaP3R4VEWEhSTtBUHp9q7eUWq3xFGW8fFkP4Ft32FuQwMs07Or69XvraUT2f/+z/wRizJFHMARp//79L7/8MgAoFIr4+PiKioqDBw+uWrVq4cKFQZg9yMTk58fk5w/1KoYBtsZGNpfrIsnr0qSioiaU7K43RCKIjua0tHjOnQObLX3GjP+Nj48fP97PjpDDcWXviM12c7lgtQKHwxkzBuvooM6etet09vHjr5wMGQxsmQwTCrmXLzsZQSJJ1j33JFZUdJ4/z46NvdKzoYGdkBCWleVHjehncVeXk6I8kZFut5tlsbDi48mWFnd3N7ulxZ2SwgsPh6YmcDrZ6enOyEgqOprT1IQZjay2tpSoqN+PH5+SmsqfM4crlcLPP7tJ0kmLLgA0NLCEQr5czu6PkcEYkQwWC0REsAQCjsNxlST43dX0u//ZfwIx5ggjGIJUVFSUlZW1ffv2yMhIAPB4PK+88srGjRsLCwv7SGdXW1tbXl6empqam5vrN1zp559/Pn36tEwmu/XWW1F5i2EHqdG0HToUlpICAP3XpAULkoqKmlQqI/xytoRg4PFALoeYGKy5Gdu8+XWpVBoTA34f09496Ra6p9HIamriXbjgpigWlws2G6utjT1qFCaTsdvbPefOAaNe0dGm6OgOk4n3zTdODAOLhYVh3LQ0bm8eDS6Xx2oFgYCy2yEmhmUwsEiSYzZzMYwymTCdjhUdze7sZGEYKyoKMAySkjxCIauhga/RUFyua9YsYWYm0JmnGNFlFpOQ0C818jYiRSLYv38/ABgMLLEYS0z0NVC89z9pepOu/hOIMUcYAX+Ou91ujUazfv16Wo0AgM1mL1my5KuvvqqqqlIqlX7v2r59+4YNG2bMmFFTU5OTk7N582afDjt37nz33Xezs7Orq6tFItGhQ4cCGmOLGEQcBoN240ZSowHaTuLzhRkZ16VJW7bokCb5RSQCkQikUpg69RpPaKanN2IxjB/Pdjp5RiPF5YLdDlIpW6lkjx0LJHmVehGE4eTJMj5f0Np6Zfs0Pt6s17f+859+5qJPrZqaMLM5ksXixsd3AVBGI+lwRHg8rXy+1WSy/PwzYbWGd3XxDhxgJSRc2WBsbxebzWcnTYq/775M5nntV0r7//3Qema1JhuNvepZz/1P6F26+kkgxhx5BFyQOBzOmDFjDhw4cN999zGNx44d43A4aWlpfm9paWnZunXr+vXrCwoK6BDas2fPTp48melgMpnef//9F1544YknniBJMi8v7+DBgwX98NcKEW7aOCQaRo1oPHa7s7NTOG5c/zVp+XLZli26ykqTXC5ASRwGER4PJkyAlBTM2y9ZKgUez1e95syZ9cMPZY2NB5mWxkZobOxj7JiurhSCSGaxRC0trvBwkiRxHGe53ZqwsItWa1tbWweLldDdPbGiQnzxIgngoii+yUTJ5VFvvfWw9/Par5T2/wP2U8+89z+v1xTrjUCMOcIIxk5Xfn7+hx9+eNddd+Xm5kZHR//www9nzpyRy+W7du0CABaL9eSTT/K8/sZVVVWx2ex7770XAMaNG5eRkaFSqbwF6bvvvqMoatGiRQRB8Hi8o0eP3rQxtsMOUq32ViMaW2OjtLCw6+zZ6/VxoJ3ukCYNIv183Mvl8n/84x/9H5Y+tTpzxnnpklunc7tckJjISkjAFArIzLxyyuVwwNmzUFnpZOyz6Gj2jBncwX1e9/MD9l+6+k8gxhxhBEOQDh48CABarVar1TKN9fX1GzdupF8/+uij3oJ0+fLluLg45lgoISGhra3Ne8Dm5maxWPzUU0+Vl5cDwMMPP7xmzRq/U2/ZsmV5/7LUBJObOQ7JYfDjmwsAHrv9uvzumCQOJSWtKIlD6MM8i+vrMYsF2tshKgri4iAh4dcncm/22VBxI6ZYMMccSQRDkL755pu+O/hEyHZ2dnpHKQkEAou3WwxAR0dHU1PT9OnTv/nmmx9++OGNN95QKBS/+93veo4cfDVCcUh9w5NI/LaLFIrr9QX3TuKANCn06c+zGD2vb3KCkSuBfS18+ovFYrP51xSNJEkmJCR4d8BxnM1mv/zyy3K5vLCwcOLEiT/8cO0MWohQQKRQiHpk1BVlZtJCRWsSF8eNKpWxrOyao+XliXNzxQThKipqCshyEQhEEAnFAn1JSUkGg8FkMgEARVH19fXJycneHaRSKYvFEgiunBxQFCUKmWxQMfn5yDzqG/lzz3lrEk8iSSwsZC77r0lqNVlcrBeLMVqTtmzRBm7NCAQiCARDkI4ePTp//vzy8vLk5GSlUqnRaFavXr1jx47e+s+cORPH8TfffFOv12/YsIEkydmzZwPAnj17Tp06BQBz5swRCoXvvPOO2WwuKSmpqqrKzc0NwgdBDAo8iUT+3HNM+HBMfr7o6tBNb01qLSnpbRyNhiwubtZoSKVSpFSKCMJVUtIa2KUjEIhAEqIF+o4fP/7CCy8QBCGVSleuXEm7jOfm5t5zzz2rVq0CgJMnT7744outra3h4eGLFy/+4x//2HMQVA8plGHc7USZmSJ/2SWdBNFUVOQiiN7qJ/1Q/EN5cflvCn9zS+EtBOEsKTHodFalUoQSCyEQwxRUoG8IuMnjkPqPkyCMKhVZUeG3zqy+uLi5uDixsFBaWAgABOH8+98vWyzu3FwxCphFIIYjAd+yQwX6EAMGw3Fxbq44N9dFEE1FRU6C6KMzjmNPPpmM41yVylhWZgzaIhEIxGCBCvQNATdzHNL1guE47f9tVKmaioois7IYd3CDwUH/m3ESvt5kd2o1qdGQmZkihSJUnGIQiJuZYDg1PP3001OmTNm2bVthYeHcuXPXr1/PYrFGaoG+jtJS7aZNQ72KkYY4L4+2k4wqVUdpqcHgWLxYXVraDgClpe2LF6tpcYJfNAkAKitN17STGLeIQK8fgUD0h2BYSIsXL54+fXpRUdHNUKAPESDEeXkipbKpqKizvPz038st5K9hswaDY+NG7bp16fQlSnaHQAxTAi5IHo/HZrPZ7fbs7Ozs7OxATzfkoHpIgYN2B698/2MP2fwH+FcqXAKASVABAGWaWWo1yey89TPZnRx0eVAuh98AoNwACMTQE/BNMzab/fDDDx8+fHjhwoXr16//mxdWqzXQsyNGGBiOO3N/Hw2dtBoBAA7ELCibB1/59FQqIwsKJHRwEkE4/Y4mA+0sKJOBNqBrDk22bNky1EsYHIL/QYbkq7uRSYfRn3UwTnEOHToEACdPntyxY8dfvOju7g7C7IgRRhy0xUGbT6MctPIeuqJURjKJhXpqksHg0KhJ+gVzBHXzsHXr1uu9pZ/PtWt267tDb+/21j6AD3KD9GfGQdcAZtIBjOy94BAXp2AExrrdbo/H07O9j3KxgwKKQxqR9OY2kv72234DbMvKjCqVUSYTFBRImASspaUdmzZp8+DoLCg7CnnVkjkrVsiHna/dqVOnFixYMNSrQAwnsrOzi4qKhnoVvRIMpwYOh4PKuSIGi97yhbd88QU3MlJAl7n2gnb+ppOCL18uBwCDwbFpk9a7D92yY4cfPQtlsrOzQ/MnFwIxMAK7ZWcymfbv3//BBx8cPnzY7XYHdK5hxIoVcmQeDRh+vJ/MQGw+n8VmN+3caSwr6xk/yyQFpxOwqtV+/Lxvzo07BCKkCKCFpNfrH3/88UuXrhw+K5XK3bt3M2X3RiqoHlKg4Ukk8hUrfHbt5CtWUG63UaUyqlSmykrv+FkapVIEACqV8b33GhIT+TgQeVA2Cc4B46cHs4L3GRAIhD8CKA+7du26dOnS888/P3Xq1L179+7bt+/zzz9/5JFHAjcj4iaBThCuLy7uKC0VZWamr1tHt4uUSrKiwq8s4TimVIp0OqtOZ7140ToPvmKcIGg/PQCQSKYOxadBIBBXCKAgVVRU5OTkLF26FAAmTpx48ODBqqqqwE0XIqA4pODAk0jow6QIL0cGOs+QSKk0lJQ4u7oYWRIplRiOM8FJtacbWyDBxytvjqQ6yB8BgUD4EMAzJLPZHP/Ldr9AIMjMzGxpaQncdAgEDYbjSQsXxhcUCGQyOttQU1FRy7//baqooDUpNZHdAGlquCrTvMNgINXqoVozAoGAgFpIbrfb+8RIKBQivwZE0BDIZEkLFzIFLMwEYa6qMqpUAplsBugug0ADCgBQgIa5xa+7xMigtra2vLw8NTU1Nze3N5fXL774Yvz48ePHjw/y2m4Eo9G4b9++J598MvizaLXa6uqrrGqJRDJ16uDs+mq12mPHjmVlZSmVSp+3PB7PkSNHLl26NG3aNIVCwWKxvN+tqqo6depUenr6Lbfc4vMWAOh0OpVKNW3atJD9Ux7hLgahCYpDChoYjscXFIhzc8mKCqtOZ9XVtiYxAAAfBklEQVTpSIIIs9vvg58vQ3IHxDRCSgo0AoAoM7M3h/Lhzvbt2zds2DBjxoyampqcnJzNmzf37HPo0KFVq1b9v//3/0L2UdUTq9W6fPny+vr6gApSb7OUl5czAacURRmNxjvuuGNQBGnlypXffPPN9OnTN2/eLJFIvvrqK+aXPUVRTz/99PHjx6dMmbJx48b777//7bffZm58+eWX//Of/9Bv3XHHHX/5y1+8h12/fv3OnTunT5++cePG3/zmN6EZIRtYQSorK2MC92pqagDAO47vgw8+iIyMDOgCECMYUWYmFBaK+lHjkalh4SQIWpmcRmPy5ct8sPPAAQA8iSSxsDDgKx4KWlpatm7dun79+oKCArow5tmzZydPnuzdp7W1dc2aNUO1woFx9uzZl1566dKlS9HR0UMyy6OPPsoUHd27d+/atWuXLVt24zPqdLqDBw/u2rVrypQpXV1dubm5hw8fpktmA8CRI0eOHTv24Ycf5ufn63S6u+66Kz8/f/bs2QBw/vz5r776ateuXVOnTj1y5Mj//M//PPfcc8nJyfSNJpNp165db7zxxkMPPXTu3LmHH3744sWLo0aNuvEFDy4BFCS5XO5yudrb2+nL2NhYAGAuASAISSJCE1QPaVAQKRR+UzP0AaNMkoICh8HQduiQj5/eyKOqqorNZt97770AMG7cuIyMDJVK5S1IFEW98sord9xxx4kTJ4ZumdcN/ch+7LHHtm3bNrSzNDc3v/POO8uXLx8U49Jms61cuXLKlCkAEBkZKRaLNRoNI0g1NTU4js+aNQsAZDJZRkaGRqOhBUkqle7cuZM20eisbD5HJA7HlTA7Ho8HAEaj8eYSpNA0CQMNikMaFmA4juF41+nTcLWf3sjj8uXLcXFxzJ5PQkJCW9tVmQA/+eSTS5cubd269f777x+KBQ6Qt99+m8/nf/nll0M+y4YNG+Li4hYvXjwoM2ZkZGRkZNCvv/nmG71en5uby7w7ceJEkiQ7OztjY2PtdntjY6NWq6XfEolE2dnZ7e3t77zzznfffVdYWCjzyloSGRn5/PPPv/766wcPHqysrJwzZ46PlRwijMASeQgEgqGzs9O79phAILBYLMxlbW3txo0b3333XaFQOBSrGzh8Pj8UZiEI4vDhw4899tigZ0dTqVSvvPLKwoULb7nlFqbxtttuw3F86dKlO3fufOaZZ6xWq81m877LYrFER0dHR0efPn2aSUpAt3/99dfR0dFSqTQ6Olqv17e2tg7uggcFJEiDTEx+PjKPEKGDWCw2m83MJUmSCQkJzOWrr76akpJy4cKF4uJis9n8448/Hj9+fCiWOVz58ssvuVzuAw88MLjDfv75588888yCBQtWr17t3c7hcLZt2yYWi/ft26dQKO68887U1FTvDjKZbPXq1QcPHmxpafn888+Z9qNHj1ZXV+/du/ett976+uuvm5ub9+zZM7hrHhSQlx0CMZJJSkoyGAwmkykyMpKiqPr6+nyvwG2xWEwQxL59+wCgu7u7srIyNTV15syZQ7feYcbhw4fvvvvuiIiIQRzzH//4x7vvvvvKK688/vjjPm8ZDAadTrd9+3YWi0VRVH5+fk5ODv3W8ePH9+7d+95773E4HD6fL5VKdTodc+PFixclEolUKgUADMPS09O97afQIXQFqe/IiYaGhgsXLjCXSqXS+3cfAoGgmTlzJo7jb7755vPPP797926SJOkz8D179qSlpX300UdMz9mzZz/11FOPPfbY0C12eEB/dXT96/r6+rlz5w7i4BUVFX/+85/nzZs3fvz4U6dOAYBEIpHL5fSko0ePfv3117u7u+fNm/fxxx/b7fa7776bvjEjI+O7777729/+9thjjx04cKCmpmbRokXManNycrZu3bpr166HHnpIpVKdOnXq448/HsRlDxpUSPLRRx+lp6cvWLBg+vTpy5cv79lhzZo1CoXill/473//27NPenp64Fc6EDZubPj22/ahXkVQ6GygvnuD+ucsaqOc+u6NoV6NL6bKyqbdu02VlUO9kMBy7Nix7Ozs9PT0WbNmHThwgG687bbb1q1b591t9uzZu3btGooFDpwvvvgiJycnyLMwX11bW1t6evrJkycHca5Vq1alX82qVau8J921a9fMmTMzMjLmzp174sQJ73v//e9/33fffenp6RMnTty4caPPav/5z38qlcqMjIzx48evXbt2ENc8iASjQN/10tLSMmfOnLfeeouJnPjss898fEIeeeSRnJyc5557ro9xUIG+IYbQwqa0q1pwOaxoGKLV3NR4PJ62trb4kZuK4qbC4/F0dnbGxPh/gDQ3N8fGxvotf+p2u1tbW2NiYoLjEjIAQtGpwW/khHcHj8dTU1MzYcKE6urq2tra/mhqSPmgD2I9pCH/XH0toGytbwuhhaM9GgM0e+BvH0aw2WykRiMGNpvdmxoBQGJiYm/FuDkcjlQqDVk1gtAUpGtGTly6dMlqtb766quPPPLI3LlzFy1aZLVa+x7Tu6p8P7mupxXTubcC2/0ZsI8Ovb01gM81uPS1AO1R+r9b1F7/8+iOBmL2gUmLz+JvHn1CIEKTUNyy27Rp07fffrt//376csWKFWw2+/3332c6VFVVbdiwYcmSJdnZ2UePHv3jH//4zDPPrOjhbL1gwQL6VBCB6A/Z2dlFRUVDvQoE4uYlFL3sekZOMKHLNBMmTGBcRGbPnj1t2rSzZ8/2HAc9XIaYr56Ac//0bcx7A2YNs7RpCAQiOITilh0TOQEAFEXV19czKQJpvv/++1WrVrlcLvqyq6tr7NixQ7BQRN/krQFcflULLkdqhEAgeiMUBYmJnNDr9Rs2bPCOnKC34EaPHr1///7t27ebTKZPP/30/PnztAcEIrTA5fD4d5D3BshnAS6HvDeQix0CgeiDUDxDAoDjx4+/8MILBEFIpdKVK1fSyW5zc3PvueeeVatWAcDevXs/+eST2trasLCwZcuWLVmyZKiXjEAgEIgbIkQFCYFAjHjcbndjY6NYLB5JddFaWlowDOvDLbsPnE6nVqu9mQ8gQnHLDoFADBUTJkzI+IUpU6b8/ve/37t3740MWFVVVVxc3LP9+PHjU6dOveuuu1auXHkj44cUNpvt97//fV1dHQDcfvvtS5cu9engt5GBoqhnnnnm6NGjAV1kKBOKXnY3SN9J8Gi++OKL8ePHD6Nqzb1hNBr37dsX0BLO15xOq9VWV1d7t0gkkkGp5ewzy7Fjx7KyspRKpc9bHo/nyJEjly5dmjZtmkKhYLFYfkeoqqo6depUenr6Lbfc0rOPTqdTqVTTpk0bAX8rbgSKopRK5X333ed2uwmCKCsre+211wQCwYAztj333HNTpkwp7FGQt6ioiMvl/vWvf01LS/N743Bk165d6enpM2bMgF+ysvl0oCjK4/H0djuPx/vf//3fjRs35uXl9fbXeIQzJAmLAsc1k+BRFPX111+np6d/+umnQV7boNPd3f3oo4/OmDFjaKfbtWsXk1QwJycnPT392WefHdypX3jhhYkTJy5atGjatGn33nuv0+lk3vJ4PE899VRGRkZhYaFCoXj11Vf9jvDSSy9NmjTpySefVCqVL7zwgs+77777bmZm5hNPPDFlypRly5YN7uKHF+PGjVu9ejVzaTabFQrFc88959PN4/H4vb1n+5133vnyyy/37PmHP/zh8ccf73sx/Z/lmvRxywBG8wtJktnZ2d999x19mZeXt2TJEp8+eXl5ixcv7mMQl8s1Y8aMgwcPDsqShh0jSpCam5sVCsVXX31FUVR1dXV6evqZM2d8+rS0tEyfPn0ECNKZM2fuuOMO+udY6Ez3r3/9a+LEiVVVVYM4tVarHTdu3OnTpymKIv5/e/ceFFX5P3D8AZI0y8S8xWDa7SBCguDQipi6aZq7ZkaiRVlRllZaZpoQ+tWapnK0ciojsVIK856VmTSQipVpmVGYDgpegUVMWRNFBc/vj+fXmR1uaUJ+xPfrD8c9u3ues8DwZs+ec57S0q5du1pXCDVNMyMjwzCMjIwM/cjAwMDMzMwqa9ixY0dgYODPP/9smmZmZqZhGPv377fudbvdoaGhy5YtM01z69athmHk5eXV4/ZfXKoE6cyZMzabbdy4cfqm2+1OSEiIiooKDQ2Ni4vLycnRy/Pz80ePHh0VFdWlS5eBAwcuWLBAL3/66ae7dOkSERExbNgwz1FiY2NvueWWsLAwh8NRWFjocDjmzJljs9kGDRq0e/fu2kYxTdPlcj3zzDPdu3fv27dvWlqa0+lcv359cXGxw+H49ttv9WOKioocDsf69evr2OAjR444HI4lS5aMGDEiJCSkX79+ixYt0neVlZVNmTKlb9++NpstISFh165dL7744siRI61tyMvLs9bv6d133+3WrdupU6f0zbqDtHDhQkc1+i+thISEO+644+y/ZY1Jo9plV+NF8Dyvymqa5uTJk/v16/fjjz9euM2sH2lpaeHh4XFxccnJyUKGKyoqeu2118aOHVu/e73Ky8uff/758PBwpVSLFi1atWqVk5OjD7xUSu3YsaNly5Z9+vRRSnXs2DEwMDAnJ0efJ2Dx9/dPTU3VexGPHz+ulKqsrPR8wKlTp/R/fH19lVKHDx++4YYb6vElXFwKCwu//fZbpdTRo0czMjL0O2OllGmao0aN2rVr14gRI9q1a/fxxx/ff//9X3/9devWrePi4q688spx48Z5e3svWbLklVdeiYiICA4Ojo2N3bJlS6dOnapM7RMfHz9z5szmzZs/9thj7dq10xPXOp3OEydOdOrUafjw4dVH8ff3Ly8vHz16dHFx8ejRo91u94wZM06cOOF2u0+fPr1z50595qJS6tSpUzt37nS73bVtsL+/v35KUlLSXXfdNXTo0EWLFk2dOjU8PPyGG2546qmnfv3113vvvdcwjOTk5J07d8bExCxdunTLli36R2jFihUHDhzQP5CeNm7cGBUV5XkducLCwkWLFnk+xpquNzQ01Lo62q5du+bPn3/rrbfqjxh69+69fPlyl8t1KU6pc4GDWK8WLFjQr18/6+aoUaM8/9YzTfOjjz6y2+3Hjh3r27fvxf4Oqby83DTNFStW/DfvkM5muAkTJgwYMKCioqLhNkPvbv3hhx+sJWvXrg0KCiopKdEb2a1bt/Hjx9f43JKSkgkTJoSHh//vf/+rcldKSkpQUNAjjzwSERHx9NNPN+hLEK5z585Vpj8YMWLEoUOHTNPMysoyDOP999/XjywpKQkKCnr55Zf37NkzadKkn376SS//+eefDcP44osv9M3adtkNGzbsiSeeME2zqKjIMAzru1bbKKZpbtiwwTCMrKwsfdfbb7+tBzpw4IBhGCtXrtTL9+7dq5fXsaqDBw8ahjFx4kR9V3Z2tmEYn3/++W+//ea5qg0bNowaNerEiROhoaFJSUmmaVZUVERHR1tP9BQZGfn6669bN3v37m3UpMouu4KCgujo6EGDBumCmn/v3bFe5iWlUb1DOnLkSNOmTa2bzZo1s/4eUUrpv8I++OCD5s2bX4itq2f/8SV7/3G40tLS9PT0SZMm1XYgyfnLysqaPHnyyJEje/ToYS3s1atXy5Ytn3jiiSFDhqxdu/bEiRPl5eU1Pr2srMzPz8/Pz2/Lli379u2z5n4uKytbvXq1n5+fv7///v37CwsLi4uL9dyalyaHw6ET8tdff/3xxx+vvvrq8OHDMzIy9JSYq1evtqY59/Ly2rVrlzVt9qpVq3bv3v3HH38opawLqZylkJAQ/Z/aRlFK5eTkKKW6d++ul/fq1avu6+HWsaoqg+r3IuXl5XrjrR+w6OhoPX/ugAED1qxZk5SUtHnz5oMHD8bExFQZ69ChQ6Wlpddee63nwttuu81zCkSlVJX37kePHh01atSZM2dSUlKsY9/1z15ubm6vXr3qeHWNUqMKUt0XwUtMTOzQoUNubm5ubu6xY8d++umnjh07Mltzffnss88uu+yyoUOHNtD6ly1bNnXq1Pj4+CpHCfv4+CQnJ7/99tsrV66Mjo6+6qqramuJ/r35/PPPR0dHL1u27LnnntPL161bt3379szMTL0np3fv3osXLx4/fnwDvRD5rrjiCj1XRfv27W+++eaCgoLZs2dv27ZNlz4iIsI6ycZms7Vp06awsHDYsGFlZWURERGRkZEDBw6cNm2aeY4nOF599dX6P7WNov6O3JEjR5o1a6aU8vPz81yDtRvWuvZ/HavS9HqUUt7e/38CjP4FUv2PqpiYmJUrV65bt+6bb74JCAjQ08VWH73Kq/by8rLWXN3JkyfHjBlTUFCQlpZW/Ye2yl7lS0SjCpJ1EbwWLVqYppmXl3f77bdb97Zq1aq0tHTlypVKqePHj//222/XXXcdQaov6enpAwcOvPLKKxti5R999NGMGTMmT55c5aMIpdTBgwf37t07d+5cLy8v0zRvv/12m81W5THffffd0qVL33jjDR8fn8svv9zf33/v3r3Wvfn5+W3bttW/EZo0aWIYxr59+xriVVykDh06pJRq1qyZPj67Q4cOem7s8vLy1157rX379l999dWhQ4cWL14cFhamlFqzZo1SyvPg5joOdK6utlGUUvqDvfXr1+uDyL///nv9FP3Jn9vt1jfz8vL+cVW16dSpk1IqJyend+/eSqnNmzc/+eSTs2fPjoqK6tix4xdffLFx48b4+Pjqx2S3a9fummuuyc/PP8uXeebMmYkTJ27duvWdd94JDg72vEtvf5WFl4hGFSTrInjPPffcwoULPS+Cd/3113u+d7bb7Y8++mhcXNyF29jGQH9h9V+LeXl5//pUlbplZ2e//vrrQ4YMCQoK0hczbNu2badOnfToN95445QpU44fPz5kyJAPPvjg5MmTAwcOrLKGwMDAtWvXpqSkxMXFrVq1aseOHfo3lF6DzWZ755130tLSYmJisrKyNm/ebF1L/tKUm5u7YMECpdTp06d37NixatWq0NDQ66+/3t/f39/fPyUlxdvbOywsbP78+atXrx48eLB+c7N+/fo2bdps3759+vTp6u+DR5RSTZs2/f3339PT0wcMGHA2o9vt9hpHUUr169fvuuuuS05O9vX1bdKkydy5c/VTWrVq1bJly9TU1JtvvrmsrOyVV175x1XVplevXh06dJg2bdqYMWMCAgJmzZrVvn17feLaPffc8+abb+r/1Pjc4ODgsw/SrFmzrK9JZmamXhgeHu7n56dXYu1OvLRc0E+w6t+GDRsiIyMNw+jTp491cLA1q7zFbrenpaVdiA2sZytWrLDZbBdqOOsLW1JSYhjGpk2bGmLQhISEKh8LJyQkeI6elpYWHR0dGBjodDp//PHHGleyfPlyh8NhGEZwcPBbb71VZfvnz58fGhoaGBgYFBQ0ffr0hngVF4ugoCDr6xwcHNy3b98pU6boY0ZM08zNzdWHnxmGMWjQoNTUVNM0KysrJ0yYEBYWZhhGRETEhx9+GBUV9cwzz+inzJs3r0uXLoZhWCvRYmNjR48ebZqmy+UyDGP58uXWXTWOouXl5d11112GYXTu3Pnhhx+2jp5Ys2ZN9+7dDcPo2rVrcnJyYGDgl19+Wceq9E/s0qVLPW8uWbJEP+Xuu+82DCMoKCg+Pn7r1q36MS6Xq3PnznWcO/XWW29FRkZaJ8n16dPn8ccfr/KYPn366GPBY2Jiqh/vsHHjRtM0p02bZrfbz/L71cg0tiCZpllZWelyuS70VuA/VVlZqY8Eq1thYaF1mkgVFRUVBQUF+mBC1M3tdhcVFVVZePLkyf3799d4kunJkyePHj1aL6NoxcXFbrfbOppOL6ysrNy3b1+N3986VlWbw4cPHz582HNJfn6+YRh1nLK6f//+kJAQz7L+C263OywsLCUl5XxWcvHi4qoALkr79u3r37//zJkz694Ld/6Ki4s3bdq0ZMmSPXv2rFu3zjp/qLpZs2ZlZmZ+9dVX//rCP3PmzPnss89Wr17teT7TpYOLqwK4KHl5efn4+DTcaQYWl8s1ceLEbdu2zZgxo44aKaX0GbvWZ0Lnqry8/OOPP37hhRcuzRoppp8AgH9UXFzcqlWrs+lESUlJkyZNWrZs+S9GqaioOHDggD7S79JEkAAAIrDLDgAgAkECAIhAkOSqrKzcs2ePdQ3jxsHlcv3555//7rn6Cs31uz0A5CBI54YJns8HEzwDqEOjunTQf8BkgufzwATPAOrAO6RzZhjGQw89FB8fr6+Yd/nll1c/7aC2YxfP/pjG0tLSkJAQu91eR5DOf5SzeUp9HYd57NixuXPnPvjgg+ezEqfT6XK5vv7663rZJACiEKTzcsUVVzRv3tz6o/7o0aOJiYk9e/bs1q3bAw88sG3bNr189+7dY8aM6dmzZ0hIyJ133pmamqqXjx07tqCgICMjIzY21nO1w4cP3759e3Z2ttPpLCoqcjqd7733Xo8ePRwOh/5UqcZRlFLFxcXPPvtsZGSk3W5fuHDh4MGDs7KyDh486HQ6165dqx/jcrmcTmdWVlYdG1xaWup0OpcuXXrfffd17dq1f//+ixcv1ncdP3586tSpdru9R48eiYmJeXl5SUlJnhfhzs/Pt9bvKTU19fTp0z179jybL+ynn37qrKaiosLHx8dut8+ePftsVgLg4sIuu3PGBM9M8AygQfzXF8+7yDHBMxM8A2ggvEM6Z0zwzATPABoCQTpnTPDMBM8AGgJBOl9M8MwEzwDqBUE6Z0zwzATPABoCQTo3Xl5e2dnZ2dnZSqkmTZq0bds2NjZ23LhxXl5eTZs2nTt3bmJiov6Nf9NNN7344osRERFnzpwZPHjw/Pnz58yZc9VVVz311FPz5s375ZdfRo4cqZQaMmTIG2+8MW7cuO+//75169aeA+n3GZ7/KqVqG0Up5evr+/77748fPz4xMdHb29tmsxUWFiqlfHx8XnrppaSkpEcffbRp06ZPPvmkLkcdq6oyqLVJvr6+77333qRJk6ZMmeLj49OjR4+xY8fq3W5Dhw6dPXv2rbfeWn33mhYSErJw4cKKigp9+Jz1AqsMoRdu2rRJKZWenp6enm7du2DBApvNlpOTExAQYO3DBNB4XOCDKhojJniuERM8A6gb75DqX4sWLaxDwiy+vr4BAQE1Pt7X11d/zHP+o2ht27ZVSpWWlnou9Pb27tChw7muqjaeR0xYEzy3adPmjjvuqO0pAQEBDz/88Lx584YOHfqvL/zzySeftG7dusppWwAaB67U0GgxwTOAiwszxqIeMMEzgPNHkAAAIrDLDgAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAIBAkAIAJBAgCIQJAAACIQJACACAQJACACQQIAiECQAAAiECQAgAgECQAgAkECAIhAkAAAIhAkAIAI/wdP22MAy0USNAAAAABJRU5ErkJggg==\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. plot performances\n",
    "subplot(2,2, 1); hold on\n",
    "x = OUT{i0}(:,:,1); % listen\n",
    "    x = bsxfun(@plus,bsxfun(@minus, x, mean(x,2)), mean(x(:)));\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",
    "x = OUT{i0}(:,:,2); % tracking\n",
    "    x = bsxfun(@plus,bsxfun(@minus, x, mean(x,2)), mean(x(:)));\n",
    "    errorbar(beats, mean(x), std(x)/sqrt(size(x,1)), ...\n",
    "    'r.', 'color', [.8 .3 .3], 'MarkerSize',15, 'Linewidth', 1.5,'Capsize', 1)\n",
    "\n",
    "% + plot staircase\n",
    "plot(2, .5, '.', 'color', [1 .5 .0], 'Linewidth', 2, 'MarkerSize', 15) %!\n",
    "        \n",
    "% Listen: plot Fit & max fit\n",
    "x = mean(OUT{i0}(:,:,1),1); % listen\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",
    "[~,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",
    "  \n",
    "% Tracking: plot Fit & max fit    \n",
    "x = mean(OUT{i0}(:,:,2),1); % tracking\n",
    "    x = polyfit(log(beats), x, fit_order);\n",
    "    freqtrack = 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', [.8 .3 .3], 'Linewidth', 1.5);\n",
    "    h.Color(4) = 0.7;     \n",
    "[~,b] = max(x);\n",
    "maxtrack = max(x);\n",
    "    h = line([xVec(b), xVec(b)], [.0, max(x)]);\n",
    "    set(h, 'color', [.8 .3 .3], 'LineStyle', ':')  \n",
    "    \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( [.5, .92] ) %!\n",
    "yticks( .5:.1:1 )       \n",
    "    xticks( round(beats([end:-2:1 1]),1) )\n",
    "    set(gca, 'XScale', 'log')\n",
    "    xlim( [.38, 4] )\n",
    "    \n",
    "\n",
    "% 2. boxplot individual preferred Fq\n",
    "subplot(2,2, 2); hold on\n",
    "boxplot(Fit,'orientation', 'horizontal', 'boxstyle','outline', 'colors','k','OutlierSize',4,'Symbol','','Widths',0.2)\n",
    "scatter(Fit(:,1),ones(nsuj,1).*(1+(rand(nsuj,1)-0.5)/10),  ...\n",
    "    'MarkerFaceColor',[.3 .3 .8],'MarkerEdgeColor','b',...\n",
    "    'MarkerFaceAlpha',.2,'MarkerEdgeAlpha',.2)\n",
    "scatter(Fit(:,2),ones(nsuj,1).*(2+(rand(nsuj,1)-0.5)/10),  ...\n",
    "    'MarkerFaceColor',[.8 .3 .3],'MarkerEdgeColor','r',...\n",
    "    'MarkerFaceAlpha',.2,'MarkerEdgeAlpha',.2)\n",
    "scatter(Fit(:,3),ones(nsuj,1).*(3+(rand(nsuj,1)-0.5)/10),  ...\n",
    "    'MarkerFaceColor',[.1 .1 .1],'MarkerEdgeColor','k',...\n",
    "    'MarkerFaceAlpha',.2,'MarkerEdgeAlpha',.2)  \n",
    "set(gca, 'FontSize', 10, 'FontName', 'Arial')\n",
    "set(gca,'Layer','top','Box','off','TickLength',[.01 .01])\n",
    "xlabel( 'Beat frequency (Hz)' )\n",
    "    xticks( round(beats([end:-2:1 1]),1) )\n",
    "    set(gca, 'XScale', 'log')\n",
    "    xlim( [.3, 4] )\n",
    "yticks( [] )\n",
    "set(gca, 'YColor', [1 1 1])\n",
    "\n",
    "% Print figure\n",
    "cd(origDir)\n",
    "% print('figure4a&c','-dpng','-r600')\n",
    "f.PaperSize = [46.81 33.1];\n",
    "f.PaperPosition(1:2) = [0 0];\n",
    "print -dpdf figure4a&c.pdf\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "optimal beat frequency passive:  0.85 %, 0.73 Hz \n",
      "optimal beat frequency tracking: 0.89 %, 0.63 Hz\n"
     ]
    }
   ],
   "source": [
    "fprintf('optimal beat frequency passive:  %1.2f %%, %1.2f Hz \\n', mean(maxpass), freqpass)\n",
    "fprintf('optimal beat frequency tracking: %1.2f %%, %1.2f Hz\\n', mean(maxtrack), freqtrack)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statistics visual vs. tracking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Parametric repeated-measure Anova on Proportion Correct \n",
      "Condition            F(  1.00, 19.00) =   0.565,\tp = 0.4616\n",
      "Beat                 F(  7.00,133.00) =  62.874,\tp = 0.0000\n",
      "Condition x Beat     F(  7.00,133.00) =   2.794,\tp = 0.0286\n",
      "\n",
      "     beat 3.82hz: p= 0.7335,  t= -0.3455,  df= 19.0000,  bf= 0.2381\n",
      "     beat 2.90hz: p= 0.8107,  t= -0.2428,  df= 19.0000,  bf= 0.2306\n",
      "     beat 2.20hz: p= 0.0150,  t= -2.6729,  df= 19.0000,  bf= 5.4418\n",
      "     beat 1.66hz: p= 0.0385,  t= -2.2241,  df= 19.0000,  bf= 2.2614\n",
      "     beat 1.26hz: p= 0.8728,  t= 0.1623,  df= 19.0000,  bf= 0.2267\n",
      "     beat 0.96hz: p= 0.1923,  t= 1.3519,  df= 19.0000,  bf= 0.5602\n",
      "     beat 0.72hz: p= 0.7173,  t= 0.3674,  df= 19.0000,  bf= 0.2400\n",
      "     beat 0.42hz: p= 0.0715,  t= 1.9087,  df= 19.0000,  bf= 1.2921\n"
     ]
    }
   ],
   "source": [
    "% Anova on performance\n",
    "fprintf('\\nParametric repeated-measure Anova on %s \\n', varname{i0})\n",
    "[~,F,~,~] = repanova( reshape(OUT{i0}, length(SUJETS),[]), [2,8], {'Condition','Beat'});\n",
    "\n",
    "% Post-hoc analysis: tracking vs. vision (beat-by-beat)\n",
    "fprintf('\\n')\n",
    "for j0 = 1:length(beats)\n",
    "    [~, p, ~, S] = ttest(OUT{1}(:,j0,2), OUT{1}(:,j0,1));\n",
    "    bf = bftest(OUT{1}(:,j0,2)-OUT{1}(:,j0,1),'bic'); \n",
    "    fprintf('     beat %1.2fhz: p= %6.4f,  t= %6.4f,  df= %6.4f,  bf= %6.4f\\n', beats(j0), p, S.tstat, S.df, bf)\n",
    "        pval(j0) = p;\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "3nd order fit passive session: r2=0.95, p=0.000 \n",
      "3nd order fitactive session: r2=0.86, p=0.002 \n",
      "\n",
      "Tapping  :  mean = 1.40 , med = 1.33 , std = 0.767  \n",
      "Tracking :  mean = 0.65 , med = 0.62 , std = 0.171  \n",
      "Passive  :  mean = 0.83 , med = 0.73 , std = 0.462  \n",
      "\n",
      "Passive vs  Tracking:  t = 1.47 ,df = 19.00 , p = 0.16, bf = 0.6586 \n",
      "Passive vs  Tapping:   t = -2.93 ,df = 19.00 , p = 0.01, bf = 9.2782 \n",
      "Tracking vs  Tapping:  t = -4.32 ,df = 19.00 , p = 0.00, bf = 209.7942 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "% optimal beat frequency\n",
    "x = mean(OUT{i0}(:,:,1),1); % passive\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 passive session: r2=%1.2f, p=%1.3f \\n', fit_order, r*r, p)\n",
    "\n",
    "x = mean(OUT{i0}(:,:,2),1); % tracking\n",
    "    y = polyfit(log(beats), x, fit_order);\n",
    "    y = polyval(y, log(beats));\n",
    "[r,p] = corr(x', y', 'type', 'Spearman');\n",
    "fprintf('%1.0fnd order fitactive session: r2=%1.2f, p=%1.3f \\n', fit_order, r*r, p)\n",
    "\n",
    "% mean values of optimal frequency\n",
    " fprintf('\\nTapping  :  mean = %1.2f , med = %1.2f , std = %1.3f  \\n', mean(Fit(:,3)), quantile(Fit(:,3),[.5]), std(Fit(:,3)))\n",
    " fprintf('Tracking :  mean = %1.2f , med = %1.2f , std = %1.3f  \\n', mean(Fit(:,2)), quantile(Fit(:,2),[.5]), std(Fit(:,2)))\n",
    " fprintf('Passive  :  mean = %1.2f , med = %1.2f , std = %1.3f  \\n', mean(Fit(:,1)), quantile(Fit(:,1),[.5]), std(Fit(:,1)))\n",
    " fprintf('\\n')    \n",
    "  \n",
    "% compare optimal beat frequency btwn tapping/passive/tracking  \n",
    "[h,p,ci,stats] = ttest(Fit(:,1),Fit(:,2));\n",
    "    bf = bftest(Fit(:,1)-Fit(:,2),'bic'); \n",
    "     fprintf('Passive vs  Tracking:  t = %1.2f ,df = %1.2f , p = %1.2f, bf = %6.4f \\n',stats.tstat,stats.df, p, bf)\n",
    "[h,p,ci,stats] = ttest(Fit(:,1),Fit(:,3));\n",
    "    bf = bftest(Fit(:,1)-Fit(:,3),'bic'); \n",
    "     fprintf('Passive vs  Tapping:   t = %1.2f ,df = %1.2f , p = %1.2f, bf = %6.4f \\n',stats.tstat,stats.df, p, bf)\n",
    "[h,p,ci,stats] = ttest(Fit(:,2),Fit(:,3));\n",
    "    bf = bftest(Fit(:,2)-Fit(:,3),'bic'); \n",
    "     fprintf('Tracking vs  Tapping:  t = %1.2f ,df = %1.2f , p = %1.2f, bf = %6.4f \\n',stats.tstat,stats.df, p, bf)\n",
    " fprintf('\\n')    \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "passive2 vs  passive1:  t = -0.07  df= 33.0325, , p = 0.95, bf = 0.2911 \n",
      "   df without Welch-Satterthwaite correction : 46\n"
     ]
    }
   ],
   "source": [
    "% load data from exp4 to compare experiments\n",
    "V = load('../exp4/exp4_values.mat', 'Fit'); % Fit\n",
    "\n",
    "% compare visual experiments: exp4 vs. exp5\n",
    "X = Fit(:,1); Y = V.Fit;\n",
    "[h,p,ci,stats] = ttest2(X,Y, 'Vartype', 'unequal'); % unpaired welsh ttest\n",
    "\n",
    "% unpaired Bayes Factor\n",
    "    T = stats.tstat;\n",
    "    df = numel(X)+numel(Y)-2;\n",
    "    N = numel(X)*numel(Y)/(numel(X)+numel(Y));\n",
    "    r = sqrt(2)/2; % Ã  voir !!! .7 or 1...\n",
    "numerator = (1+T.^2/df).^(-(df+1)/2);\n",
    "fun  = @(g) ( ((1+N.*g.*r.^2).^-0.5) .* (1+T.^2./((1+N.*g.*r.^2).*df)).^(-(df+1)/2) .* (2*pi).^(-1/2) .* g.^(-3/2).*exp(-1./(2*g))  );\n",
    "bf10 = numerator/integral(fun,0,inf); % Integrate over g\n",
    "bf10 = 1./bf10; % Return BF10\n",
    "% bf10 = 2*(1-pValue)*bf10; % if one-tail\n",
    "\n",
    "fprintf('passive2 vs  passive1:  t = %1.2f  df= %6.4f, , p = %1.2f, bf = %6.4f \\n',stats.tstat,stats.df, p, bf10)\n",
    "fprintf('   df without Welch-Satterthwaite correction : %1.0f\\n',length(Fit)+length(V.Fit)-2);\n",
    "clear h p ci stats\n"
   ]
  },
  {
   "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
}