diff --git a/Two-Classes Classification (BNCI) _CSP.ipynb b/Two-Classes Classification (BNCI) _CSP.ipynb new file mode 100644 index 0000000..3220307 --- /dev/null +++ b/Two-Classes Classification (BNCI) _CSP.ipynb @@ -0,0 +1,532 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import matplotlib.pylab as plt\n", + "import numpy as np\n", + "import seaborn as sns\n", + "sns.set(context=\"paper\", style=\"white\")\n", + "import scipy.io as sio\n", + "%matplotlib inline\n", + "from mne.decoding import CSP\n", + "from mne import filter\n", + "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n", + "from sklearn.model_selection import ShuffleSplit" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## data file names and function used to extract data" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# data filenames\n", + "trainingFileList = ['BBCIData/S14T.mat',\n", + " 'BBCIData/S13T.mat',\n", + " 'BBCIData/S12T.mat',\n", + " 'BBCIData/S11T.mat',\n", + " 'BBCIData/S10T.mat',\n", + " 'BBCIData/S09T.mat',\n", + " 'BBCIData/S08T.mat',\n", + " 'BBCIData/S07T.mat',\n", + " 'BBCIData/S06T.mat',\n", + " 'BBCIData/S05T.mat',\n", + " 'BBCIData/S04T.mat',\n", + " 'BBCIData/S03T.mat',\n", + " 'BBCIData/S02T.mat',\n", + " 'BBCIData/S01T.mat']\n", + "\n", + "\n", + "validationFileList = ['BBCIData/S14E.mat',\n", + " 'BBCIData/S13E.mat',\n", + " 'BBCIData/S12E.mat',\n", + " 'BBCIData/S11E.mat',\n", + " 'BBCIData/S10E.mat',\n", + " 'BBCIData/S09E.mat',\n", + " 'BBCIData/S08E.mat',\n", + " 'BBCIData/S07E.mat',\n", + " 'BBCIData/S06E.mat',\n", + " 'BBCIData/S05E.mat',\n", + " 'BBCIData/S04E.mat',\n", + " 'BBCIData/S03E.mat',\n", + " 'BBCIData/S02E.mat',\n", + " 'BBCIData/S01E.mat']\n", + "\n", + "\n", + "def extract_data(data_folder, training_fn, validation_fn):\n", + " \"\"\"\n", + " extract data from the matfile and organize them into X and y;\n", + " adapted from awesome_bci examples.\n", + " note: trainin data and validation data were loaded together\n", + " :param data_folder: where data is stored\n", + " :param training_fn: filename\n", + " :param validation_fn:\n", + " :return:\n", + " \"\"\"\n", + "\n", + " # prepare data containers\n", + " y = []\n", + " X = []\n", + " tmp_trainingFileList = [training_fn]\n", + " tmp_validationFileList = [validation_fn]\n", + " for i in range(len(tmp_trainingFileList)):\n", + " # read file\n", + " d1T = sio.loadmat(data_folder + tmp_trainingFileList[i])\n", + " d1E = sio.loadmat(data_folder + tmp_validationFileList[i])\n", + "\n", + " samplingRate = d1T['data'][0][0][0][0][3][0][0]\n", + " trialLength = 5 * samplingRate\n", + "\n", + " # run through all training runs\n", + " for run in range(5):\n", + " y.append(d1T['data'][0][run][0][0][2][0]) # labels\n", + " timestamps = d1T['data'][0][run][0][0][1][0] # timestamps\n", + " rawData = d1T['data'][0][run][0][0][0].transpose() # chan x data\n", + "\n", + " # parse out data based on timestamps\n", + " for start in timestamps:\n", + " end = start + trialLength\n", + " X.append(rawData[:, start:end]) # 15 x 2560\n", + "\n", + " # run through all validation runs (we do not discriminate at this point)\n", + " for run in range(3):\n", + " y.append(d1E['data'][0][run][0][0][2][0]) # labels\n", + " timestamps = d1E['data'][0][run][0][0][1][0] # timestamps\n", + " rawData = d1E['data'][0][run][0][0][0].transpose() # chan x data\n", + "\n", + " # parse out data based on timestamps\n", + " for start in timestamps:\n", + " end = start + trialLength\n", + " X.append(rawData[:, start:end]) # 15 x 2560\n", + "\n", + " del rawData\n", + " del d1T\n", + " del d1E\n", + "\n", + " # arrange data into numpy arrays\n", + " # also torch expect float32 for samples\n", + " # and int64 for labels {0,1}\n", + " X = np.array(X).astype(np.float32)\n", + " y = (np.array(y).flatten() - 1).astype(np.int64)\n", + " sr = samplingRate\n", + "\n", + " return(X, y, sr)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### other parameters" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "data_folder = \"./\"\n", + "cue_time = 2.5 # cue starting time\n", + "signal_delay = 1.0 # take signal from 1 sec after cue onset\n", + "signal_dur = 1.5 # signal duration \n", + "signal_st = cue_time + signal_delay\n", + "signal_et = signal_st + signal_dur" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## load data from one subject\n", + "\n", + "### subject 11 gives best result (> 95%), subject 4 gives the worst result (60%). \n", + "### bad channels greatly influence the classification results. \n", + "e.g., the 2nd channel of subject 6 is bad. After removing it, the classification accuracy increases from 50% to 90%." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "subject_id = 6\n", + "trfn = trainingFileList[subject_id]\n", + "vlfn = validationFileList[subject_id]" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "((160, 15, 2560), (160,), 512)\n" + ] + } + ], + "source": [ + "X, y, sr = extract_data(data_folder, trfn, vlfn)\n", + "print(X.shape, y.shape, sr)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### visualize the raw response to check for potential outliers" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWEAAAD2CAYAAAAH6zaYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJztnW/MLNdd3z9nZvf5d+99fOPEDnGg\nhZL2lJg4JsQlUW3ixFaDKhWiImjVNAQ1UilCKHVrFLWVKvEOolZWS/oC9UVLGyUtoCqIVJDGpCa2\nYsdWDaIkcFoCWG1qE8z1vdf3zz7P7pzTF3POzNm5++w+3mf37j57vx9rvbNnZs78zp89z9yzM58x\nIQSEEEKshmLVAQghxK2MBmEhhFghGoSFEGKFaBAWQogVokFYCCFWSG9ZGVtrt4H7gBeBalnHEUJs\nBCXwJuA559zBvJlYa28H9o+x6WXn3IV5j7NIljYIUw/ATy4xfyHE5vEA8NQ8O1prb98/d/bPLr96\n5Tibv2Ktfcs6DMTLHIRfBPiug2/nPNscEAjADoYBgR0M1/BN2kVGbFNQEdimpKQ+fT6gwmAY4TlD\nj4pAiaGgnkvxwJD6WucSQw84INDH4OM2BTAEPAFPYEhglzLu0+4/InCOEp8VwhPYpmBIoAD6GA7w\neKCHwQChSQ/NPrsUXMWzS8F1PD1MPH57nCGeAYEehhD3S1dt52kpjwID8Vgpr5RvP9v+ChXbsXz9\nWPY+MIh1t4PhChWHeLYomvxTffUo6GMYEdiKZRwQOMRTxhgMsEuBj++XY37bFPQw9DGco+BlRs2c\n10Gs/z4Fw1g3I+pOeDXuW3XqoIj9YITnPD0uU7EdYyC2/yGBM3FfME2bemAr9oMhgSqr671Y5gEV\nO5RsxfJWQBXLvU0ReyhN3R4QmjiquK6HuaFfjoBBbIfrVJRxmwLYpmTUlJemrtO+qa3LGKuJMeX9\noIxlT7HtYLiKZ4+Ca3h2Yhu23426HkvgMhXnKDlLyeX4j1RPoIjtth2/mykmYj5n4/dgj4KLVPTi\ntgBX8JRAgWGAZwtDEfczsc+nVkvf0wGBEZ4+BYEhz25/DeK4MSf7l1+9wn/4Nz/LG+94w5Eb/cmf\nvsyP/MTHXkd9xrzRg3AFsMMWe1l33qXAxC9goGq+xNdDwY6pv0g7lE2nNXEQHgbPrskHYdMMwr3Y\nHevGNRSxYVPHqgeWegCo988HYdPsX3ewHn6sEIGdbBDeivl76i9PGoS3MJSxw1dxUPBU7MXu34+D\nRPojsEePQzzgxwbQdOw8LeWRBuGtLK+U71bcviIwYsROMwgbenGQJv4xqAe/ESZ4tk3BLmXTFvUf\nmoKtOFBux6MaPEU2CBfAXvyDtUfBkBFF8+WvB4CzoeCqKZtBuK6RwBYFJZ5dSkZxYKyo2GnK1dZB\nEfvBEM8efYbxj3XCx3z3Yl4mDlCHsX624x+KwziopHrda4alil3K5g9tXX/1fjudQbgeRD1nYhyj\nuK4/NgibOAjXe+zEus0H4Z34B3iHsqmPYRxcC2i27w7CeT8o4x/fFNsuBVXsb6nfHTUIHzJijx5n\nKBnFQTjF349/pME3fd5gOIx1PIx9+4ARfYpmEK6ymKFim4Iy7mfidyu1WvqegmcUTwR8+4f1xFOX\nb3zD7bz5m+44eoPgj163ApY5CAshxM3H+/o1bf0aoUFYCLFRBAJhytluO+G1HmgQFkJsFtWofk1b\nv0ZoEBZCbBbeg58ytazpCCGEWCLBT//xTT/MCSHEEgkzfpjTICyEEMsjBD/9hzkNwkIIsUSqasYP\nc+tlUZg5CFtr7wI+C7wVOOucG8X0R4AfdM7dv9wQhRDiNeCrGT/MrdcgfByL2gXgIeCZlBDlPPcu\nKyghhJib9MPctNcaMXMQds4NnHOvdJI/AvzCckISQogTEEJ719yk15o9V/M1+4SttX3gQefcF5YQ\njxBCnIxTdiY8zw9zHwI+ddyN+yQbVhKY1MaoQWMFq4UpW6ZoxDTJgGaojV4ApUlCmloMMqS1iuXi\nmUm3JCbNThKMmChNOcQ3kh8fBSO5GCVJcmoJSRKV1Pl5/JhpZAiNOKaKkp5eFKgUGEIjDzJRXdTG\n2ZX3JFJayoNYJ4le/FTLYVqrXLKgDWPdpm2SBS6l7ZhWYlRG8U36HOIxk3ms3i9ZwuraSHXpY7m2\nKQjUt4V64MCERsqT6iRk9TmKy8OsfmrLnG/KlMpcZnXZ1k9NEuaUTSu15fUxHpPVVyuUMfQoxvLs\nN3Vhmv4ArWkvqajKbJ/Uknl/rI9jMhMZ+BAoTdGIhoZNOUMjCPKkPpjnXDOI4qjaEtea/Oo+2n5v\nktWtFSfRtAuxnfJ+ntq0bPJq90uk8tZSnjYu37y3lrTUlu3xTTPQ1O1U101thavzOGCBZ6enzB0x\nz5M1LPDj1tpfB+621v7kgmMSQoi5CX5E8MMpr1N223Kcfvg14O3A54B/6pz7WFz3lHPu55YbohBC\nvAY27WYN59wQePiIdbo8TQixXui2ZSGEWCES+AghxAo54ZmwtfZHgA9T/wb6QeBR4J3A8865j8Zt\nHuumzYseeS+E2CyST3ja6wistW8G3uOce8g59yDwRuo7hR8Atqy191lr39FNO0m4OhMWQmwW6WaN\naeuP5v1Aaa39DeCrwO8Dn4/rHgfeTf0c127ac/OGqzNhIcRmMe1uuVnXENdnvlvOuYeAa8BtwOW4\n7hJwPr66aXOjQVgIsVGEUM18TeES8Jtx+QvU9/bsx8/7wMW4TTdtbjQICyE2Cz/DHeGnTkd8Cbgn\nLt9LfbPhQ/Hzw9Qis6cnpM2NBmEhxGZxAneEc+63gevW2ieA+4B/AQystU8ClXPuWefc8920k4Sr\nH+aEEJuFn/G05Rm3LTvnHu0k3XAJ2kkvS8vRICyE2Cz8jKsjpk9H3HRuwnREa1/KzVMAW/HwRbau\naAxSyb5konurrjjfsS35uK6fFaXIDFHjUSSjWzJ9tYax/D1EI1W9XJufkhmsiss9TPPfMFqsik75\nfHac5PfqR6OcyY6XW7kKWttX+rzVyTfPH2qjVjJwJUbkpqtkb2stXym29HmYLad6Saau0QTDW8p9\n1PHWlVm9pZ8/khmsn9VZKmde932KphVT30hGMLIYUzzDaMHLy5uOGTKDXG3ZSy494nZt+08y2KX6\nrOI+yf6XjlFFa5iJ+ydrWN7vUlmSka1nTGz/to1TXx3FeNLeyWKX2iy37tUmtzzv1m6XjHGtia3t\n+yNCY1DLzWrJHDdsPGt1n0t5pLpKfSPVwyG+KU9qq2SeK5uWrvffiiVJx00mtYLWsrgwbgGVpRBC\nrC+nTGWpQVgIsVlsmkVNCCFOFSHMcEes15ywBmEhxGYxww8xdd0K0CAshNgsNCcshBCrZMZ0xCKf\nZ7cANAgLITYLnQkLIcQK0SAshBArpKrq17T1a4QGYSHEZnEyqftN5ziPvL8L+CzwVuAs8N3AY9R3\ngz7nnHtkqREKIcRr4ZQ9bfk47ogL1O7M5Mx8AXhffNz9ndbaty0rOCGEeM2c7MkaN52ZZ8LOuQG1\nOzN9filbPaR1pgghxOoJYfqUw2mbjjgKa+09wB3Oua8uMB4hhDgZt4LK0lp7O/AJ4COzth3gGRK4\nRkUVlX89DNsUXIsn0QVwNYwY4BlQMYwSvmH8XBG4Fp8LdYjnMBMP9qOa7wDPNapGvQg0x0vLQwLX\n4/6jqArMFZhpmy0KhlHxZ4BzUcwXCI1WchjVhQWwS8FOVp5+1DZWBHYaXWdeH54iUzSWcbloYvBN\ner2uLlOtAmyVl2n/fqYnLDEMQr3Gx7pOxxrF117c/lqoGnliEes2SUOTIrHEsEXBgIpRbMthrL/D\nqLrcwjRtlbSdfQz7oWAQy+JplY7nKBud5jYFt9GLukma+j/Esx0FlAWGQ3xWl3U+tfqyredrUTQ6\nrkNN+XquZ/9o24n7HmbyySLW84CqiSXvGz0M1/HsZH1gGDWqJh6zVaXW+/SzWMjyPMDTp2jaqI9h\nh5KdqKlsdZDt8hlKerGv71KwTRHzD00b9WKZ6v6TYq/1lNvU7XGVEcOoN/XAAaHplwbYDQUDAldi\neybNZaqfs5RN3ziIfQBaxWhqy1QPBhjE3porRQdUHOKjxnSB+BmPu58hdb/ZvOYzYWttD/gk8Ghn\nakIIIVaPD4RpZ7trdiZ8nKsj+sCvAW8HPgd8kfrZSx+P88T/xDn39DKDFEKIY7NpN2s454bUTxTN\n+enlhCOEECdEKkshhFghIUyfctAgLIQQS2TTpiOEEOJUIXeEEEKskE1zRwghxKnCz5gTPm2XqAkh\nxKnilAl8NAgLITYLnQkLIcTqCJUnjI7+8S1UOhMWQojloekIIYRYIZqOGMfEVy/ansr4nkxLwJg1\nqqI1hZUYehS1jcoUTZrJ8k8msYJxw1iBaWxggRsrPdm3QjStJfNWAY1LzGPoNTHVpUk5mRh3HVNN\nP8Y6ykxkyejm4z6J5O4K0TOVjp1sbXnEITOspfKlGNL6VCsh1lU5drTcbGVifCHWpWnyS22Q7Fop\nPx/LZrL9kkXLQNNmJtZpMmIlU1dqpyrm1ZrVQqyXOt/c7AXJhtea5tq2m1SPsBNbwjf1NF4HKd4i\n9sHxMkOV2dlSekHbh6E25CUzWJHllaxhKfa2bcbbOPWV1J55P0pOu7x8reHNjFndUr8aTqjrtg+Z\nsXyaeg0BTLcP1t9RgJFpPydjXS/Gb7K663fqN6/vUeyTPourS+p/qR8sDN2sIYQQK0S3LQshxAqR\nwEcIIVZHGFXTr46Ysm4VaBAWQmwW+mFOCCFWiOaEhRBiheg6YSGEWCGeGdMRs7Ow1j4C/KBz7n5r\n7WPAO4HnnXMfjetvSJuXuZ62LIQQ60qID/qc9pqGtXYbuDcuvwM465x7ANiy1t43Ke0k8WoQFkJs\nFlUFoymv2VL3jwC/EJffBXw+Lj8OvPuItLnRICyE2CzS1RHTXkcQny7/oHPuCzHpPHA5Ll+Knyel\nzc1xHnl/F/BZ4K3Up+CjRc6HCCHEQjnZJWofAj6Vfb4E7MflfeAi9Z363bS5Oc6Z8AXgIeAZmDxH\ncpIAhBBikYQQZr6mYIEft9b+OnA38Abq8Q/gYepx8OkJaXMzcxB2zg2cc69kSQudDxFCiIUSmD4V\nMWUMds59zDn3fufc9wFfcc79NDCw1j4JVM65Z51zz3fTThLuPJeonQf+MC5fov5rIYQQa0EYecLo\n6OvQpq3Lcc7dH99vmHJd5DTsPD/MTZojOZJ+1EGWmRavVU6aRnPXMyaq81pNYE4ZFYJFzG8SZZZ/\nyicn6QZzfWPeHCHTIgKNljLXMI6ybUKjwmy1gD5um7arxvKvMbEsvpOeuFF7mcp0ozowL4+Jy+WE\n8uf1HuJ2fVM0ZTbxv7wu8v1N896m5jpSOsujWF8lSYc5rivMdaO5uDO1bTpOEZdMp+w+2z5fl7c/\ntG3ox/Zt2yzFlsqc1vtm29TWmQryCFXmeF9qj9UlKSHpbN/dtso+p341SdOZJKPpu9TvtGb+nSqB\nbVPv3Qtt/2m1pq2uNdfOmhhDQa3NnBRHNVaXScFqmjJ7Ws1rqqui2W+BnOCHuVUwzyC80PkQIYRY\nKGmUP+q1XmPw7EHYWtu31j4OvB34HNBngfMhQgixSE56s8bNZuacsHNuSH3Gm/Pl5YQjhBAnRBY1\nIYRYIWnaYdr6NUKDsBBiowhVIIyOPtsNlc6EhRBiecya99V0hBBCLBFNRwghxOo4Zc/51CAshNgw\ndCYshBCrI4wgTL6ptlm/TmgQFkJsFJqOEEKIFXLKnvOpQVgIsVmctkF46Y83KmgNSmTLwzg7npvF\nki2rwEwwbtUmp2SK6hqsurTmsTbP9Cozx1RuP4MbzWOJ3LyWG6ommcFSfGn71gw2bnhLseX7FtAY\n51JcXbtWbqCqfSTJ9Naat5KnZNSpd+I2IRqxclJ+ef2160xTlmQ2K8fqp4072bO8aU15KT3FmuIY\nZZazbllzs1iybY3XdWtyS+/J9mWaOjeNeS83syVjWl7XeTsQy3ejja61jOXGsBR/Momlvpwb01LM\ndPLL6dZ7vl1uLuv229RuVVa+vNeErP6a75hp+3UyrCVrWjG2743/fu9+R+t822Mbxvtct9xdw9po\nwjHmJpjZrzVCZ8JCiI1Cc8JCCLFCwsjgp5zthkpnwkIIsTRCMIRpg7CmI4QQYnkEP+M64TX7YU6D\nsBBiowjBELzOhIUQYiWEMP3HN/0wJ4QQSyR403lM7Y3r1wkNwkKIjcJX06+O8BqEhRBiicyYE9bN\nGkIIsURCYPp0xCbMCVtr94BfAs4Al4Afds4dLDIwIYSYh5mXqK3ZIDyvO+L7gC875x4Eno2fhRBi\n5fhgZr7WiXmnI74GfE9cPg/82WLCEUKIk+F9gZ9yfukpmDJbcdOZ90z4fwPvttZ+BXgn8KXFhSSE\nEPOTrhOe9lon5h2EPwz8qnPubuC/An931kHG1YytxnCyqrFOyzWXMK427O4ziXSsPM9ESfvHsH1v\ntYYm/tfVCvYyReW4xtI0GsUb48i3a3WcSSPYajFv1P5VsR5Mtl2+Pukwc9XmpEbNlaJdbWCuiEyq\nyUntklpk1MQ5rrosMzVhXpau9jPpGlP6MGo1kx5xFFt/kt6wewJTNn1jXIGZjp3vN0kjmddb9xhJ\nUdnqScd1nzDer0Jn3zaeuh+OQmjUozmpvvP+nZPnn3SnuZqyW9Zc5ZnK01W0Ju1mt8xtmer3to5q\nHWg51uZkfa79rpSZ7hNa/Wk5FlP9ubyhRReAr6+OOOrFhlyiZoALcfll4LbFhCOEECfDB4OfMrhP\nW7cK5h2EPwX8Z2vth4Ah8LcWF5IQQsxPCDPumNuEQdg5dxF4/4JjEUKIE1NfJzxl/U2L5HjoZg0h\nxEZRBUMVjv65qzIbcCYshBBryxpeATENDcJCiI1i1g0Zm/LDnBBCrCUzH2+kQVgIIZbHSc+ErbXf\nAzxGfUn4c865R6y1PwX8APAC8KPOueGktHninfdmDSGEWEvCMV4zeAF4n3PufuBOa+17gPfGz78D\nfMBae2c3bd54dSYshNgovC+mXh3hzfRzT+fcS9nHIXA38ET8/DjwQeDqhLRfmideDcJCiI3CM35r\n+aT1x8Faew9wB3Ax2+0StbTsPHC5kzYXmo4QQmwUATPzNQtr7e3AJ4CPUA+y+3HVPvWgPCltLjQI\nCyE2Cg/4MOU1Y39rbQ/4JPBonJp4DnhPXP0w8MwRaXOx9EE4mbGSTStZyJKVqcgMS0W2Llm2UoCm\n+Tz+V2zcjJbsZPnxW3NZP+6dLFRHTdCnRkqWqdpkFk1YE8xVyTyWytvNZ1Ilh2ybSXEkO1dukJsU\nYzp+HkNujktmrtQGeYym857aoczaJS9r2q5Hbps7uoy+U1e5gaw1g4Umhnz7ZDcrsthSXeU2thBN\nX21dTG7/ZCDL691n9VUx/s/Y3HTmO/Xa5j1eT4Zxs1mXwpixemzL2h4nb6/c0Ja/t/3qxr6WSPF3\n26UXvwe9pp3z72Vbln7crms5a+uyjT239+XpeYzT+smiCdQCn6NexzgT/iHgPuDj1tongG8Hvmit\nfQq4F/iMc+4b3bR549WcsBBio6iAaspAW83Y3zn3aeDTneSngZ/tbPez3bR50CAshNgoZs376mYN\nIYRYIou6OuJmoUFYCLFR5L8bHLV+ndAgLITYKDQdIYQQK8Qz/TFymo4QQoglUmFmXB2hM2EhhFga\nmhMWQogV4s2Mpy3r8UZCCLE8ZukqN+ZM2Fr7I8CHgRL4oHPu6wuLSggh5uSWuE7YWvtm4D3OuYcW\nHI8QQpwIb2YMwus1GzH3mfD7gdJa+xvAV4F/6JybdUu2EEIsHT/j6oh1e9DnvFKjNwJb8Uz4GvVz\nloQQYuV4M/u1Tsw7CF8CfjMufwH4jlkHKTsaQRPTcqVhri5MJJVeyDR/JuaVK/26ekafKQ5zfWJS\nIsKNur20TYq5oFYTjmsjW/1k0kPmc1DjysFWgZn2zd/TcfIYcq1nXgdJkxhodZr5u4n/5XF3FaHd\nsgZCbJdc25jqKdWRoftDRq6crLL9cpI2cjSW+41lzmMGGsViL2un1O5dPWPSOo4IjDrtlNotvUxW\nX22fS+0Xmj7araN0rNRXuxrJvLd2FZ9FTMvLaLK2SPnnOtfu+JDrMrs6yKR2Td+lvJ1TmXNVaq5o\nzVWexLoZxldFW99V9p62Se8HMcdR7Ot5fad6TvWe1+m4/nK8/y0Cf4zXOjHvIPwl4J64fC/wR4sJ\nRwghTsYCHvR5U5lrEHbO/TZwPQqP7wN+eZFBCSHEvIQZUxFhzaYj5r5EzTn36CIDEUKIRTCKr2nr\n1wndrCGE2CiCmXGzxqacCQshxDoid4QQQqyQW+KOOSGEWFduGXeEEEKsI7fKbctCCLGWVEy/AmLd\n/AoahIUQG4WmI4QQYoVoOkIIIVaIro4QQogVctqmI+YV+Lwmkr2q7Dii8s9lx0pVAH2KMTNUvne+\n7KMNbJIpbLK1atwulfZJRqtkdcptUF1Dl+l8nmSCSheN+2ihSiTjV8jS0zY+s8WldUUW6yhLy8uZ\nx5AMWmm7G8s4TteklvCdesj3y2PLSYasbp0lO1pr+xo3fnnaPjDpTGXczlWTLGhF8/nGMqRYUgy+\nE3eeZzL5ddNTrL1OXzIxBjPh2MkCl46XLHe9xtSWly007Z0MZ3m/gNZG52NeuQEtjOU1busrsuVk\nNTOd7XtZSqqDZH/L66JLMvXl7Zl/B5NJr2sPHGZ9I1kRFzlDEBqr2+TXoq1tJ0VnwkKIjaJi+hUQ\nujpCCCGWiOaEhRBihQSmXwGxXpMRGoSFEBtGmvudtn6d0CAshNgoTtvVERqEhRAbRXrm4LT164QG\nYSHExrFew+x0NAgLITYKXR0hhBArRD/MCSHECjnpD3PW2seAdwLPO+c+urDAjuBEty1bax+x1j61\nqGCEEOKk+GO8jsJa+w7grHPuAWDLWnvfksOdfxC21m4D9y4wFiGEODG5h2PSa8Z0xLuAz8flx4F3\nLzvek5wJfwT4hUUFIoQQiyCJk45+TeU8cDkuX4qfl8pcg7C1tg886Jz7woLjEUKIExGO8ZrCJWA/\nLu8DF5cUZsO8Z8IfAj513I0nFbrA0I+6v6SWTCS9Hdm6pBSclGfKIykFkw6S+NnQqvVMpt3zjYqv\njSktl9k2uVrST9g2xZDiynWCSUOZKxWTLrBVBdYqTpPpHkNcnxSQuZqyyParOrWbSwFN53Mea1fN\nWWWxpXwC47rNLkXcN/81OsWW6izln68z3PgLdauQrNP7WdytUjNkytP6fdjE3ZYlbds968nVqemC\n/Vwtmfe5Oo8WT6sgTcdJZQzZtiPG+x6d5bw1Rlnsbdw37pvrF1PfqRgXMqZ+n2tWk9LUZPn2Yv/P\n67efKTZLDFvx1c/6Y9qmjO/9mM92qJf7nXUpnzJbzrcpYz2nz/l3+6RMPwueOR3xNPBQXH4YeGZh\ngR3BvIOwBX7cWvvrwN3W2p9cYExCCDE3rcd78mvaEOycex4YWGufBCrn3LNLDne+S9Sccx9Ly9ba\np5xzP7e4kIQQYn7yf9kdtX4aN+OytJwTXyfsnLt/EYEIIcQiCEx/eoaerCGEEEtEty0LIcQKCQR8\nmHImbHQmLIQQS0M+YSGEWCES+AghxAo56dURNxsNwkKIjSLMOBPW1RFCCLFE0p2e09avExqEhRAb\nhS5RE0KIFRJCIEy7RG3NzoU1CAshNorTdnXEiZ6scRwOovlqgOcA3/xTwBN4lREFtQnrMK4/jFVU\nERhmaelz2jZ3LhXAdTxXqfBx3xLT/ApaAQcxr2tUDGN+yZQ2iuYpT/2o7B6GQ3y0RxXsUjQGqWR/\nSg1pMJyhZIeCChrjVHnEa9TslwxvdXpuAxtQAa1J7OCIbtOa2tr8k2GriBa0RHKsAuxQsE3BQajb\nY4uCHgVDfGOZC5mxrJ/Fl0xinsAg2x5ggGdA1ZRpN9R14uM6Hw1gKeZU12cpx36x7mM4wLNNQQH0\nKQgEtmOMKa5U1x7YjTm3JrQ2rl7sC9eomn134tq8blNsQ0K0hLVt3falQBmPnbdlsv31YwRDAqMQ\nmro0MYZB/A4MY12MYj1vUdR1FsvcoxgrQ1rfp2BIYI9yrJ+l/HcoCDG+krrv1/0Stim4hudVRlzD\n0wt13tfxHMSypPp5NW53GL+xQ0Lz3Uv9dYuCK8Y3361h/P4M43f1MLOWXaViOFaPZPXRmt8WwQml\n7jcdnQkLITaK03YmrEFYCLFRaE5YCCFWiK6OEEKIlTJdZbluVwprEBZCbBRVCJhw9PmublsWQogl\noh/mhBBihejJGkIIsUJCmCF11yAshBDLQ1J3IYRYIfUdeBs+J2yt/R7gMepL7p5zzj2y0KiEEGJO\nquBh6tUR63Wl8LzuiBeA98XH3d9prX3bAmMSQoi5SVL3o14bMSfsnHsp+ziEaJwRQogVc0tJ3a21\n9wB3OOe+uqB4hBDiRNwy7ghr7e3AJ4AfnrbdXfS5SMWVMGKEp28Ml8OQPgWvhkPuKHaoCFz0B+yZ\nHgbDOWPYp+QwavGgnnzejrMnHhoNYIHhOhUHwXPG1MXZoeQQH7WAhv1QUJlaHbhFwSGec5SNvq8X\n86kInKVkSKBH4CAqDStaRWLSBA6jnrBWKxZRJ1krG5Oi7zDU5TWZcnEUVYYvhQPOmB6B0Gj9Uvn6\nFGzF5auMmro8R483hJIhcN143hAKXjRDPIE9epwNBf/HHDIMnsKUUSNa/xPsMKoh9ym5wKiuO2MY\nhIoBFWdNDzAMqLidHn0MX+eAt4QdrptaNdiPP3dcoWIn6j1fZkhJj8thxIvVFQCGoeJz/+q9sH87\ng5//Rfwg8PWv7OMO9vm9rbr9HFc5wHPVD3lreRu/e/in/OX+G9gxZVNeT+BiGHEpHPJHBy9zrTrg\nxasX2N/a5erwgCp4+kXJteEB/bLHHbv7GFO307ft3snAj/AE7uydxRN4aXiZ0hT8+f55tih5xQ/Y\nL7YZYCjNFhfDEKjVlJgenoI9Cl7vC75RVNwW+tzp6/yHpqQMMDAQkyhD3S/P+vr9lWKLMsDQ0LxX\nBG7zhgNTf/m2Q/1PyZGBKya2LcLhAAAJ8UlEQVSwHwzXDFwxtTJy3xu8gbsHFb+zU/dNDxyYwHYw\nbGM4IPB6b/h/ZcE3+YIXisAeBRepqAicp0cBuHCFv2DOUIXAjjGMTPqnueEcBd9cFVwu4I/NIX8u\nbPEVDumZuk+8GoacN1tUsd+foWz+eT/A40avYDDsFn0O/Ijdog/AGdPnG9VVSlOwa/pRD2rombp/\nnzU93l7t8LVycQOjn3F9xKb8MNcDPgk82pmaEEKIlXLazoTn/WHuh4D7gI9ba5+w1r57gTEJIcTc\neDzVlFf7aIn1YN4f5j4NfHrBsQghxInxBMy0M2GzXmfCullDCLFRhBlzwus2HaFBWAixUfjA9DPh\n9RqDNQgLITYLnQkLIcQK8cHPMPj45T9m/jWgQVgIsVH4EKbPOYSgQVgIIZbFrOmIdbtxWYOwEGKj\nCMeYjpgHa+3fB/5e/PivnXOfijeu/Tvg24DPOud+ZlLatHzX6axcCCFOTKB9ztxki9rc/Dfn3LuA\nB4B/HNO+H/j9aJS831r7TUekHYkGYSHERpFuW572mgfn3B/HxVF8AbwL+Hxc/u/AXzki7Ug0HSGE\n2Ciq4AnBHLl+2jXEx+QfAL8Sl88Dl+Pypfh5UtqRLH0QvsCIPbbpRwPTGQqGps82BbebPiMCfQrK\nYpc9Sg7ifd2D+I+G/FQ9zeQUwDAazkYEtijYMSV9TGY4g0MCBYFXTPqnSJ3LkMB1PIOY0qfgKqNo\nQjONua2KFreCvEENIeZfxXUDPL24bxltbLfRYxjNbSnWIYFdehzgOW9qy9Q1KgoMfWorXC8eaxjj\nTHmmbf9vtGuVGK6aETsUDAm8SsXA1Oa4njFjeZUYdii4SsWFWLYQ660wtXXuGhV9DFsUXKbCALuU\nfMNUzfZXqZryHcayn6NX29lMj/3e+ejkgn/26O+yS8GQNzXmr2vbh4zw7FCyT7+u63Kba1S8betO\ntinwsc63KagI3G76nDM97tr9FvoU+LPJfNeWMdn0SgxFLO9hbNsSw0H0Zn3r9lmq2N88gdvKuv9t\nU9DH8PrYR1Mf245t9qdFxYjAy2bEhXI8z1QfB3j2TN1/e0Udy4AKE9cXpu6/fQy+pLGJpX4xSkY/\nY2Lf9fQp+JOi7gN/tAuGiiGBvWgJLI1p+tcrZV1fg6LuNwNG0cZX518RuNPsALBves0/yQ2GAsM1\nPL9X1nW2R8nLpuJ19PHU/7TfMgVnaA13V6nqOqTHHiXf0bu9+T4Oi7pd+rGeXt97XazT2joItRGx\njN+d3yoPOGA4afiYCx/C1IsjZg3CcfrgP3WSX3LO/e34VKG/Dnwgpl8C9uPyPvAHR6Qdic6EhRAb\nxexnZ0xfG82QD3bTrbVvBv4l8P3OufQgi6eBh4BngfdSO3V6E9KORHPCQogNY9Z88NzTEf8ceCPw\nX6I9chf4VeA7rbVPAU875148Iu1IdCYshNgo0hTKonHO/dgRqz7Y2W7YTZuGBmEhxEYRZswJ1+uO\n/uHuZqNBWAixUVTe46fcj1FfOFEevcFNRoOwEGKjCM2VUJNZn3PgGg3CQoiN4njTEeuDBmEhxEbh\nQ8BPvU745sVyHDQICyE2ilnXCa/ZGKxBWAixWVQ+UE09FV6vWeG5B2Fr7WPAO4HnnXMfXVxIQghx\nEmbfM7dOzHXHnLX2HcBZ59wDwJa19r7FhiWEEPMRZtwxt24D9Lxnwrmq7XHg3cBznW1KgOscAlCR\n7mIpGFDho2xlBAyj6AMKDmMFDZsLSQIHsdp8TKulJa1PLol0RhDfQ8ytFaWkfHtx34DhIDvWAZ4i\n6kySwMdn+SRqmU+IxyJuV+ebjjPEU1JFsUwba3o/IJCecjWgFvIkz2mS9RhqoY/J0lI5y3isWnZT\nl7eidesN4/4hNoKPdXc9HqvexmOoZTGeggGeUZLNxPRJ7ZCENSMMPtZVoOA6SfRDzLNeP4pCHh/L\nOsITKOP/TVPGUawTH2U2adljGMa2qI9J046j7D21cxL4JFFMapMUV5XlE+Ln+hhm7E6rItZLav9U\njqKTZ0ndtw9jjz5spEtwEOu4ivt5YIRpYkl9dhjjIMZuYvuMYh2U2T6j+B2q5U75d2G8P/jYdqmv\ntNKlgut4CgqGUdzTj/smGUI/6yM+xl0RYu0S+20txepT4YHr+Ob7OGralqbPQCveIvaPIn4H69gO\nU9YnvoC3KKCckkuxZrKGeQfh88AfxuVLwN0TtnkTwDPbUwVCQgiR8ybga3Puexl45c137b3uGNu+\nQqubXCnzDsJdVdvFCds8R22gfxGaP7JCCDGJknoA7v6L+tg45y5Ya99COzZN47Jz7sK8x1ok8w7C\nTwM/Bvwi8DDw77sbOOcOgKfmjkwIcasx7xlwQxxY12JwPS5zzY44554HBtbaJ4HKOffsYsMSQohb\nAzPv85aEEEKcnDX7nVAIIW4tlnbH3Gm9mcNa+63Al4HfAw6dc3/NWvtTwA8ALwA/6pwbTkpbVcyT\nsNbeBXwWeCv1Nd2jSW1y3LR1oFsm4JvptFXc7lS0V3xe2WPUV4E955x75Lixn6LyXAJ+K27yN+OP\nZx8EfoJ67vbvOOcuT0pbQRFWwlLOhDfgZo7PO+cejAPwncB7nXP3A78DfGBS2iqDPYIL1M+5egYm\nt8lx01ZVgAmMlSnStBXAKWuvF4D3xbjutNa+h2PEforK8zbgf8b2eTAOwH3qpxV/L/AfgR+blLai\n+FfCsqYjJt3McZp4r7X2SWvtI9RnhE/E9FSWSWlrhXNu4Jx7JUua1CbHTVsLJpQJxtsKTlF7Oede\ncs4N4sch9fX2T8TP02I/LeWpgO+I7fMz1loD/EXqgXlEG/uktFuGZQ3C52kvhL4UP58WXgT+EvVT\nUh+m7vDdspzG8k2K+bhp68pYW1lr7+EUlinGfQf19fanvo1SeZxzX6UeYL8XeB3wNziF5Vk2y5oT\nPs7NHGtJvL75AMBa+1nqzvHmuDqV5RL1fGSetu5MapPqmGlryYS2+k4mt83atpe19nbgE8APA9/N\n8WI/LeVJ1+1irf0M8F3Ar3Bj/zq148UiWNaZ8NPUc3dQn00+M2XbtcJaey77+FeBPwDeEz+nsjw3\nIW3dmdQmx01bSya01deY3DZr2V7W2h7wSeBR59xLHD/2U1Eea+0Za22yOKT2+V/Uj4MvaWOflHbL\nsJRB+JTfzPGAtfZ/WGu/BHzdOfdl4IvW2qeAe4HPOOe+0U1bYbwTsdb2rbWPA28HPgf06bTJpHZa\n57abUKZ/1G2rSW2zxu31Q8B9wMettU8A384xYj9F5bkHeM5a+0XgW4Bfjldx/FvgSeDDwM9PSltB\n7CtDN2sIIcQK0c0aQgixQjQICyHECtEgLIQQK0SDsBBCrBANwkIIsUI0CAshxArRICyEECtEg7AQ\nQqyQ/w/x0bP21Y90wwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.pcolormesh(X[0]);\n", + "plt.colorbar(); plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAEDCAYAAACWDNcwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJztvXuwbFl93/dZe+/u87iPuTMwMwxI\nMZIcrZiXJtJMJGwRRkAZopJFYkeobAmhhKRccqwCFCIspSq2yxXZ1qOILFIVK5WK5BBkWzIWCEuD\nBdRIYMEwKUyIwFnBQJBEGGCYmTv33vPo3nuv/LEevbtPn+4+9/TZ59x7vp+pnu7ea+31+61H/+4+\na6/13cZ7jxBCiH4oTtsBIYQ4TyjoCiFEjyjoCiFEjyjoCiFEjyjoCiFEj1Sn7UAXa+0G8CDwZaA5\nZXeEEGebErgPeMw5t3+zhVhr7wIur5D1GefckzdrJ3Gmgi4h4H74tJ0QQtxSvAz4yM2caK296/Kl\ni19/5tr1VbI/Za3908cNvGct6H4Z4PrXDb41AIzaGoDChO9b5ZC9dhyOxZM2yyHXxntTBW1VQwBK\nU7DbjGL+UEaLZ9yEcjfKAQAGk20km6UpMEwfGxQFVRGa7Ua0uVFWlKYE4Mm9awBc2timbpopGwC1\nb3LZAY+JNUlplSnxtNFXOnUNvjS+jWWYfG43fyo5nTssKpq2nbJhjMn+ecJa7cvDbQCuj/dy++X2\nB5q4pju1XVmUDIrpGapx2+bPZSctnZO4MNjMPiX7jW+5MrwYfdgFYFCUc/v76mgHgEuDLQCujXcp\nYpumejy9f53taiOWnf5wMoxmfBmWZa7bMPZt6/2BPtgbh3F0eWOb3XrELN0xCrDXjvnGrWcDcLUO\n/tZtjYl+jppQr9IU+dwijqO6rWlju2zG8gZFSRlbIbXJuK1Ja+3v3bwTgCdGz3ChDPV+anQ912ej\nDHXbi3YNhibWLpW7VW2wU6ffUii3Kiq2YnnXxjs5rTu+ACpTHPgdXhpsRpsj6Pz+ks3dej/bBWho\nsy9pXAyLKo/51I9VUWEKz9adDcS4cZNcfubadf7R//j3uffuZx+a6Stfe4If+a/edifhivi2CroN\ngG8nQbed/FYA8MZMjqWDpugci/kKcyC/iQO79Z4Yb/DxGJj8OeUvCpNtpGOtN/gU/NIxTMwLdR0G\nSltNzpnYgBSTUn7vJz/WlOYLQ+unByiYif/tpD753Kn80+d6Ou2Zz+3UKbVZp81T+3Xbuo0/7lRv\n4yd2Z+uX0hPNbP+UE59StrYF2mLKF8/8/s79V3b6IrkSy2i69Yjlmanxk+o/qVvqW+9n23RSB98e\nLKPrXnccFb6a8om2ADNTR2PysW5abj6T/sktc8BODvjWkPY3lR1bpiinbXgO2A19Gj6lcURbTJ+T\n2mTWZz89vpKfB9qlLKK7Zur3l2zmNk391Pk9pPp7DN5P92O2GTj2VOS9z76L5z3n7sMz+PbwtCNy\n1oKuEEL0T9tOXzHMS18TCrpCiHOPx+MXXM161ieXoKArhBBNHV6L0teEgq4QQrQt8yfqO+lrQkFX\nCCF8u/hmmW6kCSHEGvFLbqQp6AohxPrwvl18I01BVwgh1kjTLLmRtj5VghMJutba5wLvA14AXHTO\n1fH4W4C/5Jz77pOwK4QQN0XbLLmRtr6ge1IqY08CrwQ+lg5EMZv7T8ieEELcPOlG2qLXmjiRoOuc\n23POPTVz+I3Ar56EPSGEOBbeT3alzXut8VmSvejpWmsHwEPOuQ/1YU8IIY5Ej1e6fd1Iez3wrlUz\nD4oSj2HcNlldKP37MG6brI6VBTTwWWkpHavjHExZFh1FrwmzT0H2+Kw0tVePox8tm1FtK4lwGGNy\n2V2Sn4Oo5DRqJupPbSdtv4l+JmEVfE5Pvo/bOvswjOVV5YBxR/1slioKnOyM9/M5iaIj5jOPZmap\nTGEMg6h2tdMGFaiRb7NKW6rrRjFglmFZsTPen/K963Nqu2ExYN+P4rHYPh01qW57Xx8F5arNqNbW\n+DYrhSXRlQKT2zspV3nvs7hKalvvfbaR2qwwRUe5Lfk5zmUnFa1r7W5ur9nx1m3Hrk/PRHWx3WY/\nH6tm2qIsC8ZRPSu1aIvP6bWpc72r2C+pvLptcj2+PgoKdzv1XlY6S/UviyIrxnXHXvI5CTCNo3pZ\nl4EpD2yDvTHeZ7MaHKhHPif6lNvCFHn8Zp/KYqr9Eqnvt6ugIHd9vDsRdkoiO96v9eqzT+2Fvp4c\nYYEfs9Y+DLzQWvvjPdkVQoil+LbGt+MFrzO+DThOJ/wO8G3A+4Gfds69LaZ9xDn3SydhVwghbopb\nfXOEc24MvOqQNC0XE0KcLbQNWAgheuQYgjfW2m3g14ELwFXgdYue2aanAQshxPFWL7wGeNQ59xDw\n8fj9UBR0hRAi6ekueh3O5whXuQBXgK8vyqygK4QQx9sc8VngpdbaTwMPAH+wKLOCrhBCLAq4y9bw\nwhuA33LOvRD4F8APL8qsG2lCiHOP9w3eH34jbVEa4TnQ6bHsTwB3LMqsoCuEEK1fsiNt4fTCu4B/\nYq19PTAGfnBRZgVdIYQ4xjpd59zTwKtXNaWgK4QQ7ZIVCmd9G7AQQtxSHG964UicyaB7o97Pm0OS\nClJNONBViRpHtbGtasg4KR0V0wpOo6Y+oArWdspIalWNbye2Yv7CmKyMlMpv2jafm973mKiWjVN5\nbZt96SpiJQWupMK014wPqHc1vs1KTLv1KPud/MoqX+UAmK7bjfEesJnLBig3iuxDeg91S3UK7ztR\nnWvU1OyYiYpVqn/yOatfNU1WfUrqXIUx7EWfk/pVt8/2o0+79f4B1Snv62wj1XtQlNm/vVjHC8UG\ndaxHsrXXjLMPe/s3st2qGE+1GcBu7IOtwUY+ltolq9Q1TW6/CzFf6s/WT8ZAO+fHmMbAjfE+d21c\nDvWI6mBVMfnJpfYpTUEZjxUxX+jZpD4X7FamPKCGtu/HxGw8a+tSrus491tq48lCpdRno6Y+MPa6\nn9PvYewbNs3GVNqoGedxkOrBkI7S22iqPXfrUS4vKZbtN+Pc94m6bWjKaCP67r2flDMO5Y6KmrKE\nC+sKYdoGLIQQPdKjtKOCrhBC3OoqY0IIcUvh/ZLphdt8TlcIIXplmb7CYu2FI6GgK4QQmtMVQog+\nWTK9gKYXhBBifehKVwghekRBVwgheqRpwmtR+ppQ0BVCCL9kG/BZXzJmrX0u8D7gBcBF4DuAtxM2\nLD7mnHvLSdgVQoibosdtwCf15IgngVcCH4vfvwi8Ij5+/R5r7YtPyK4QQhyd4z054kicyJWuc24P\n2LPWpu+Pd5LHzKq0CCHEaeL94imEsz69cBjW2pcAdzvnPtOnXSGEWMjtKO1orb0LeAfwumV5L1Qb\n+CJIuV2Pgn5JVu9itcluE2XjoixdVZQMy2oq3zBK6JWmyBJxQQoxyB+OO/KNIW0jf04ycl15xrKz\nDbCIsnpJoq5um4msYJS0u2f7jmx32JHzKzdDvu0qyC9eHGx3ykl+DtltgvxgkgmsijLL4I079cmy\niLGMOzYvcN/WXQB88fpXczslycRUx9b77DOxXbaridRhaoMLg82Yv52S/0ukdq+iJOGN8R4Xh1tT\ntjAmSyBuVsPsb5YEjOzXY+7dvBOAUfu13HZ3b4dHTiX792xeyfVOfbpXj9mswufLw/A07Cd2r+b0\nNp5bmCLXO8kQlkXBwJfZXuJCMT37ltI2ysFE2rFzBZTGSraJ51s37wHgqWYXgKvNbu6rVF7r/WTs\nxTFd+4bNcmMq33Y56Z8bTfhdVKbK5V0uQl89VVxnqwjtvB/7r24bNsphbNtxPLfM56a+uFBtZenH\n2sTfSPSp68tdW5dyG3R/a3UR/4iNQ+XK8GLwo36KzSiRmX57g6LkRpTZvBj9HLV1Li+Vv1UNs8xm\n6uONcoApPGv7o/l2EzG31lbAO4G3zkw1CCHE6dN6/KKr2bN+pWutHQC/A3wb8H7g94EHgZ+N87w/\n5Zz76EnYFkKII3Orb45wzo2BV80c/tsnYUsIIY6NpB2FEKJHvF88haCgK4QQa+RWn14QQohbCmkv\nCCFEj9zq2gtCCHFL0S6Z0z3rS8aEEOKW4piCN9ba1wB/I30Ffsw595vz8iroCiHEMa90nXMPAw8D\nWGsfBT5wWF4FXSHEucc3Lb4+/GaZb1ZbvWCt/WbgK86564flUdAVQoj16en+ReCfL8qgoCuEEOu7\nkfYXCIH3UM5k0G2j9lHTWcJRx3VybeWz+lBSi+qqICXlKNNRDEv596NSkTEmn1PNqG8B7ETlo0FR\nUkad97ajDJXKvjGKSk9lOWUPgpJSHdWUuspVyZfGN1P+AlmVDCZqZYmKaUWublndejRNy9f2rk61\nT1kUFD6qi8XBUxqT2zeVM09NKim6FZisqJXeq8JM+iDWe1hO2qdbXpLLTzY3ykFun1xeWbIfFbBS\nuVVRsTO6MXXsRr2ffTadfkv4pH5VVln5izkXKlUZ2qwyJb5I7TJRqcs2Yj326t0DZRRz7E+qbLK6\n2JN1qMOorfM5+7Ftq6LKT/hO9sdtkxXAUv6dqDyX2iDkG+dzrkZbo7ZmowhqXF01tNQuvqNM52eW\nQtW+yce6bTxbz9R3XaqiojSh/1ozXUb3N9JV88sqfp3fcsqXlOsKU1KY9al8zWUNmyOstc8BRs65\nry/KdyaDrhBC9Mp6tgG/FnjPskwKukIIsQbBG+fcP1zFlIKuEOLc4+tm8eqFBWlHRUFXCCG0I00I\nIXpE0o5CCNEj61unuxQFXSGEaFkyvbA+Uwq6Qohzj1/yYMqFD608Igq6QgjRNFAXi9PXhIKuEELc\n6qsXrLXPBd4HvAC46JyrrbVvBx4APuGce9NJ2BVCiJuix6C74Hr6WDwJvBL4GIC19tsJwfdlwNBa\n++AJ2RVCiCPjvV/6WhcnEnSdc3vOuac6h74L+N34+QPAS0/CrhBC3BSeydXuvNf6Yu6JXenOcgV4\nJn6+Gr8LIcSZwNft0te66OtG2lXgcvx8GXh6lZP8nH9eupf63fRZ6bn854CZSAyajtRgkmpMx1rv\ncxndPyXSuVMyhZGuhGEqp4gSda1vJ1KITMsGdm20c+T1WnyWQMzShHPsduvTLTvVYxDlHgemZEwz\nVY8Wn8uu475yM0cqs+vTbBvPkzUsTZHrm9rC47OvE9m+Kv+Tn2QCK8q5ko2pLUwZjtV+IvOXfE11\nnSWlp3vPZVFQtIf3ZWrHrmxmaqdmzgJ5M2fMdPti1wepwySRmbyCaSlPY6b7ufVNtjeOEofdS6TW\nH7ybPort0viJLGTyq9tXuW39ZGxOyvV53Ka2K80cWdG2zdKYKb/pyDLOjtuqU0Z3/BYz46JLmX9L\nzVSbpvwLVDWPzm0wpzvLRwlzvACvIs71CiHEmcATN0gc8jrr0wvW2oG19gPAtwHvBwbAnrX2w0Dj\nnPv4SdgVQoibIW2OWPRaFycyveCcGxOuaLs8ehK2hBDi2Nzq63SFEOKWIk0jLEpfEwq6Qohzj288\nvl6gvdDoSlcIIdbHsnlbTS8IIcQa0fSCEEL0xxqeS7kyCrpCCKErXSGE6A9fg1+ww62zCfLYKOgK\nIc49ml4QQoge6fG5lAq6Qghx7oNuUBIL70kdadypdVY1Mgev+ZMiFR1hpFkFo6KjbtSlnfkboqsq\nlcpog2MHfZ5REpunSJXq1s3XVZPq1iGrbfmJ+lUzZ2Kp6CidpfJn62FMQTFjo8AcUJjyXWWtGVWn\n3K4dSlMcUA8rTEET1bOKIip8eX9Q7csUzK5FL4w5oHBlMFmNK9kYNfWUSlsinZvarutzV2Ur5Rsk\nVTVjKHyqR5HzzPbVqmQFO2PYb4PKWFJG6467pK42LKrc9i0TP8ZNOHczjuVmjnLduG0oy+BfUjI7\nTA3tYH8fHMetb3JbdZemZrt+YjepjCUK01HbS/3d8WVW+cwbP1elLZfXVcSb6QN/yO/wpvFmyaTu\nYkkza+2PAG8gRJ4fcs596bC8ZzLoCiFEnxxnTtda+zzg5c65Vx6ea4KCrhDi3ONrQ7vgatY3C690\nXw2U1toPAp8B3uycO/TxwUf6u8lae+Eo+YUQ4lbAe7P0tYB7gWG80t0BXrso80pB11r7vdbax4CP\nW2sra+3/umplhBDirJNupC16LeAq8Hvx84eAP7Mo86pXuv8t8N3AV51zNfCnVjxPCCHOPN4bfLvg\ntfhK9w+Al8TP9wNfWJR51aDbOOf2mTy0Yv4DqYQQ4hYkLYZY9DoM59wngV1r7SPAg8BvLLK16o20\nd1trfx34d6y1/xvwz1Y8Twghzjy+Nfg5y0i76Ytwzr11VVsrBV3n3P9grX0/8ELg/3bO/eGqBoQQ\n4qzTNotXL7RLgu5RWBh0rbV/cdY28K3W2m91zr17bV4IIcRpEud0F6Wvi2VXui+O7w8A28DH4+er\ngIKuEOK2wHsWTy/0JXjjnPvbANba96TdFtZaA7z3qIastdvArwMXCEH7dfHmnBBCnCq+XbILeI1B\nd9XVC99orX1O/HwP8LybsPUa4FHn3EOEK+bX3EQZQgixdlpvlr7WxaqrF94M/Ga8Wt0B3nITtj4H\nfGf8fAX4+k2UIYQQa6dtC9oF16AtBQtmH47EqqsXfh/4rmPa+izwUmvtp4GvAm87ZnlCCLEWwpzu\ngnToN+haa78T+EVgi3Cl+xPOuY8e0dYbgN9yzv2ctfatwA8D/2hexjZKOxpj5rZEMSvVaMosETcl\nwch8ybh5ZbXeZ+m5VRkUkz0iyU461pVXTGmFMRP/OraGRTV1rGnbLGXXreusvN2wrKbk75LdNu5Z\n3KyG+bxmpn22q02e3r8OQBV9rqIfVVFTmHLKZlkWbJQDAPaj5ODlwQX2mv2pcsuyoPZB62Oz3Mi+\n7tR7AGyUk3pVpornBFuNb3JbjIpBtn9xsAnAxcF2LHeQy031vzTcyjKCG+UwNlCo57T9YZZUvFBt\nARPZxW59t6tN2liPJPeY7G6UQ8axjG7/JPup/I1yQBXbcSu2RXdcTMZskaUd0/gZtxOJxVRuaSrS\nD2KjGOY0M2esJInFroxlGofD2Lbe+yybWXQkJbsynSF/lf1KbTDYLHMblEWRbQ7LasbniW9JKnLe\nb3OZfGaqY1Wc0L6sJet0wRxRqeZwVi3mF4D/1Dn3bcDrgJ+/CVsGeDJ+fgK44ybKEEKItdPnnO7K\nsds59yfx/Y9v0ta7gNfFrXI/BPzvN1mOEEKslWOqjB2JVW+kPWKt/U3gUcLc7iNHNeSce5qgOymE\nEGeKleZ018SqQfftwDcC3wo87Jz712v0QQghTpXGGxp/+B/+jen/SvfXnHN/Hvjk2iwLIcRZYc2P\nXFvEqkH3/7TWfj/wUeKzBZ1zTy4+RQghbg2W3Sxr17VejNWD7gPxlTZFeOAVa/NCCCFOkWU3yxYv\nJzsaqwbdvxk3SAB53a4QQtwW9Hmlu+qSsb818/2/XpsHQghxyvgVXutimZ7uXwb+CvASa21SFhss\nOEUIIW452rZYuHqhXbJj7igsm174bcJD195E2AYMMAIeX5sHQghxyrTxtSh9XSwM3865q865LwKf\nAv4ofn6coKMghBC3BR6z9LUuVr1mfr1zzgPE99evzQMhhDhlWqD1C15rtLXq6oUL1tqBc25srR0C\nF9fow6EYDE07Xd1x20yUmuJ77eucL6kbJcWwFp+PdRW5uqpcKX9SPUr5aSfqXb6Tr+sLBGWodE7K\nX3Rm3ttOubP1adqWsZ8+1jJRfxpSHTiW2G/GWU0qqUDVbZM/79RBAawclAcUsK6Nd7Iq1W49AmDU\njuN7ndW+Uv5x2zBq6mw3tUlWyiIpUxn265C+X45yGan99uK5lwbbWd0rlVsYw05ULUvqZcOi4un9\nG1PnDsuKnXFIvzzczj7N9t9uPZr0ZWSn3sv1TWmVKbMyWhXbYL8ZT9TIospZsr9T7+U2KGZU7QC2\nB0FRbK8esVPu53MgqGklBbc8Vn2bx00aU61vDvje/b7fhjqMmzordHX7ePbcrqpd6udUv65dY8yB\nsT5q61y/1AZP79/IimNVVIkbFgP2Ytsmkhpa19aB31m0kfxM6a2Z5EvHUjmmMJg17mbwmIUrFE5j\nydg/AD5mrf0EcD9hW7AQQtwWNECzILA2h6YcnVVFzN9lrX0Y+Bbgc9qNJoS4nVg2b9v7nK619puB\ndwD/HfCMtfZNa/NACCFOmXaF17pY9UbaLwM/A1xwztXA96/RByGEOFU8iwPuOjdHHEXE/A/XaFcI\nIc4MfS4ZW/VG2h9Za/86cMla+18Cn1+bB0IIccq0QLsgri6bXrDWPp/wkId/A4yiFO5cVg26fxV4\nYyy0Av7aiucJIcSZp8EsWb2w0pXu7zrnfnhZpmXaC3d1vv7T+AK4xOQhk0IIcUuT5nQXpa/A91hr\nPwy82zl36LLaZVe6/2yOPYP0dIUQtxGtWbw5ol3+uJ4vEx5ntg+8x1r7Qefcp+ZlXBh0nXPfA2Ct\nfRD4WcJOtD8L/J1lHgghxK3CMvnGZVe6zrl9QsDFWvs+4EUEzZoDrLp64ReAHwCuOefGwIMrnjeF\ntfZHrLUftNY+Yq193s2UIYQQ6+a463SttZc6X/8c8LnD8q56I61xzj1hrU3fjywuGYPsy51zrzzq\nuUIIcZK0Zom04/L7aC+z1v4dwtXuh51zjx6WcdWg+5i19ueBe6y1/z3wsRXP6/JqoLTWfhD4DPBm\n59w6tzQLIcRN0S5ZvbDscT3Oud8m6I8vZaUrVufcTwIfAH4F+FfOuZ9a5bwZ7gWG8Up3B3jtTZQh\nhBBrpzXLX+ti1StdnHMPAw8fw9ZV4Pfi5w8Rni68kLs37+CPr38NmEi6daUVvT/42c/Ivc1+T8eS\nzGO3vHlyeOn8JGHYzZNl/fAT+cYk3VgUWXIv2ShNkevRtZ/K7srftXPk7w7UyUxLWKbz0ufkiz+k\njHambbq2Wua0o0lmzfRxoCyKif05bZ7q6OdIaaY6NN6zUQS5wGt+59B6Fx3Jz9ZPJD1Tm3alPNs5\nkpzdcg6zARO5w9k0Y0yWz+zmKWbucDe+5WIZZCH3m1E8t+BCFaQfr4938rFipo+qooJmWibxYrWZ\nP++PJmnJryTHOeXznHp3x7KZc1c+p3dkSFN9c93ahraspvL5TnnJ1t3DywA8M75xwE43X5fZ8eO9\np4z2TZbSbCnWuDn3zDw5Ys38AfCS+Pl+4As92hZCiEPp88GUvQVd59wngV1r7SOE1Q+/0ZdtIYRY\nhF8ytbDg6exHZuXphXXgnHtrn/aEEGIV6vhalL4ueg26QghxFvFmyeaIW/VKVwghziJr0l5YCQVd\nIcS5p8/VCwq6Qohzz3G1F46Cgq4Q4tyzhm3AK6OgK4Q49zQsXqHQ+yPYhRDidkbTC0II0SOaXhBC\niB7R6gUhhOiRcz+94GlpveFPbjxxQO2q8W1WJkoKT41vs7JUek959tvxRJGqqygWj1W+zN/Tuanc\nrrJSUgDrKiClY8aYbG9KDa2YzrfHOJedfPHe52OpjBafFauST3XbTKk+pfypHNMeVCqrO+2TP8cy\nymKieDabVjcNdTFd39bPUwWb+GM6W3Zyejvpk0V9ltIa3zJq6ymfxkVDm8oponqYmUiGdNuniv21\nX4+zzxtzfE0+DDvHRk2wm5S6vPcH2qdbr3l91sShUXXGym4b1MCqWG5piuzLRhkU1SpTgilj3Sbt\nmHxJ43DU1tmHpJAW7E6PC4OZUhdLPqdjqa9a7ymTKtgh6n3he6e+nXGb1cVM7FvfTPUlwFf2nwq+\nN3VW3Zv9DcCkjUdNDfFnV/qJct3Yh/TUt6Up4s9rWv3sZvH4A7FmNn1dnMmgK4QQfdKweIWCVi8I\nIcQa0ZyuEEL0iGfxCoXbfk5XCCH6pF0yp7so7ago6Aohzj3nfvWCEEL0SY2nXhBaF6UdFQVdIYRg\nvVezi1DQFUKce7R6QQghekQ30oQQokf6vJHW2yPYE9bat1hrP9K3XSGEOIx2hdcyVo1tvV7pWms3\ngPv7tCmEEMto8TTHmF44Smzr+0r3jcCv9mxTCCEWEq5m/YLXUlaObb0FXWvtAHjIOfehvmwKIcQq\n+BVeh3HU2Nbn9MLrgXetkrHxntZ77tu+iy9e+2o41pE1PCD32JFl/JbL9wHw5d0nD6QlmcLCzJG+\n8222keT1jDFTUonAge8QZO6MSXJ5oQzjzQGJvNo3NDOSgHvNOMv0jTtll2b638OufF/dkZ6sm/C5\nKKftd8/pSjumtmvbZkqmMtivc/4ksbiItiN/mOrftG0ud78ZT9U1pCe5xHoi3xjf95sxfiPkTVKL\npSkoitAWhkkdU/+NO7KUKT3VuzAml5NsDcvJkE9+GjORrUz19nEMduvx7K3LAGyUQ+4YXgDgYrkJ\nwOP7T3Gx2go2oiTjtfEOV6qQL7V7gckyhXdVFwG40e7nMdCVrdwqhrkNUlruozKkjX3DRsy3UQSp\nyO+4/E18Yf8JAO4cXsr12qn3ALg02I7l7R4YZ9vlJkUVfLk6vgEEWcqUL7Vx2ZE9ner7NL46kpcp\nLUlAJuq2mfRZlO1s2pbRvKeVxWHd/S2HIb8eacdjrl5YObZBv9MLFvgxa+3DwAuttT/eo20hhDgU\nz+KbaEtWLxwptvV2peuce1v6bK39iHPul/qyLYQQi2iW3EhblHbU2HYq63Sdc999GnaFEGIePv63\nKH0VVolt2hwhhDj3aBuwEEL0iMcfuLE8lW5Wu9JdBQVdIcS5R3q6QgjRIxK8EUKIHjnO6oWjoqAr\nhDj3+CVXuquuXlgFBV0hxLknzOkuCrrrQ0FXCHHu0ZIxIYToEe/9Aa2UqXRNLwghxPo496sXtqsN\nfGG4Or7BZhWUk8Zt0Oa5NNjmBrsAVIOgMLRRDqmjctPX9q4CsBVVmChhtxnFfKGsum2y+lRlYhnD\nQf6XLqlKlaZgswrlJHWjzWqYFaFSPmNMVmFKSmXP2rqUbWxGu8YUOd92tZl9H7XTalyDYsDV0fVc\nNsDFajOrUyW7G+WAqpi0QfIzqVwRzLJRDrPdeephyc87hkH1aqfey2pcuS2KiTZSV8Xr0nAr1nED\ngOvjHTY2Lhw8J7ZLOrZRDql60IPvAAAPKElEQVSin02s10Yz4FmDqIq1GexulRtZSSyde+fGpG0v\nDEM7Fpjs84VB8OnaeCcratWx3l3FrKqjlJV+CVdiG+w3IwbRv9R/X90NY+ve7SFPx/55yl/LRVzz\nO7luEBTkLhahXR4fX2WWp5pQxrCostpWl9RXF2LbVqbMPjdNm+udVPESn77xJ1wZhHo8WQf/vPd5\nzO234ffQ+Iki3Fa1kY/txN/LOLbxRjHMbZveN6oBgzz2QvtcGmxnJbN9Qv8l5bW9wSj7nhTfqqLM\n4yuNvarYzb/dIv42S1Pk33dSWdushpjCwzxFspvguCLmR+FMBl0hhOiTc3+lK4QQfaI5XSGE6BGt\nXhBCiF5ZLO24zpW6CrpCiHNP4z3GH349q23AQgixRnQjTQghemRdT45YBQVdIcS5p/v057npCrpC\nCLE+JGIuhBA90uIxt9ucrrX2O4G3E5a8Peace0tftoUQYhGNb2Hh6oX1rdQtlmdZG18EXhEfUXyP\ntfbFPdoWQohDSSLmh71uyTld59zjna9joOnLthBCLOK2FjG31r4EuNs595m+bQshxDxuW+0Fa+1d\nwDuA1y3K99CFb2GLIZ9rnuFF5RUAbsQL43sZcmcbpOyux2fR39kavnQxzLncvx/SYhYeGzb5kvqb\n21DdB0Z7PBl1D/+/KsjH/alxzZUySNq9+I0hnx/VlA9+OwDlA68B4L0vfTufCip4XIsl79EyiNJ8\nP7AXjr17s+Lvvfc/A6D58G8BYO59DuWD3xvK3g8ygMVdz+WvPfC2kC92rMfzdJTfS9KOz7T7DOJs\n0OUoF7jja0ZR8u5KOrY9kbrb9+Hz2DdsmlCnLEHpW3b9OH6eHlDlpmEvSuhdLoMc4F475mKUGOza\nTOc+0YT6jIeXGURJvjviuePOXNlOrNfQlNxXBjm/r7VBqvPr9XXurcKxZ1dBHrLF873b3wLA59og\nhVhguG8YxsU3VXcA8Pn6ab6hvAzAF+qnAfhzF76JT+0/Hn0JEoOFMWxthb5/Xhls/L/NNcrYtqPY\nZgUmS2mmeidpRxP/g8U/xsKYLOl4rQ51LE2R5USTPOTlwYXcL6m80hTc6JyTfEuymul9vxlnKcnr\nUVZxnnynMYY21qeYIyPpU7mdtCTFmOyHfMG/UVNnycvki8EckM1Mvhhj8lhONL49UG+DyTetNpKM\npW/ZKAYHfFkn7ZL1C7fqjbQKeCfw1pmpBiGEOFWOe6VrrX0R8MuEadN/C/znzrm5J/V5I+0HgAeB\nn7XWPmKtfWmPtoUQ4lBaWpoFr3b56gXnnPuzzrmXxe8PHJaxzxtpvwb8Wl/2hBBiVVo8ZtGVrll8\npeucG3e+7gN/fFhebY4QQpx7/JI53VVupFlrvx/4GeCzwNcPy9fn9IIQQpxJWg9t1F+Y/1pehnPu\nvc65FwF/AnzfYfkUdIUQ5x6/wn+LsNZudL4+A/HpuXPQ9IIQ4tzT+naJ4k277Ar1Ndban4ifPwv8\ny8MyKugKIc49rfew4EYa3i8Mus659wDvWcWWgq4Q4tyz7EaanpEmhBBrxK8wvbAuFHSFEOeeZYI3\nBzdO3zwKukKIc0/YBrwww9psKegKIc49jW/x/vDr2UW71Y7KmQy6//IZh28N47bm0TYoI42aoFZ0\nx8Z2/pwEKi4Nt7g6CipXTdvmYwD1jYa6CWWURbj/2HbELTaqiVJSOnf/F8KOPoPhjg0XvQo7mPeb\nMWbHTNlK5wP8z9G3zesD3vmyfw1MVJ02ygFf3fllAJ5/+V4Arte77DfjqfoYY7g+CopRSVXpWVuX\ncr5U/61qmB+mt1WFZYJfuv4EVzaCetZezH/P1hX2m6DutVPvZxuJ3XFIu7IZzhs3NYMyDI06tv9+\nPaYqyyn729VkaeIw5m98y7VRWKKY/Bg1NUW0d30c6nXf9l18vvlKSI9KVKO2Zm8z+Py1vaAUdnl4\ngQ8886loI/TVpeFWtvF/pf5rWz4VbaR+/jftH3NxEJTOnvBXc73TuV3/66QoViVltIan924AcPfW\nZbq0vsl/iib1qe5YqIom1yf1/bAI7VOZMrdFsp/UurqUpsjnpPxT6mZxvLV4Rm1os63BMJ+fnnTQ\ndOYik9pe8qluGwZFNZXPdFS80rHun93p2Lip2a1H2VeAUVVP2mNmDnS/HtOU4VgaP4OizOMw1af2\nTXZ0GJXFRu0YU5ipcxvfUpSwvaY//NslV7q3fdAVQog+Wb79QUFXCCHWyJIrXQVdIYRYHy3rFSpf\nhIKuEOLcs2z1Qkhbz/yxgq4Q4tzTtC3tgv0PYWFDuRZbCrpCiHNPeAT74WhzhBBCrJHVphfWg4Ku\nEOLcs0yofMnTeo6Egq4Q4tyzbJ3uOtc1KOgKIc49TetpFl7qrm9Wt9ega619O+HRxJ9wzr2pT9tC\nCHE4y/ekrYvenpFmrf124GJ8LvzQWvtgX7aFEGIRHh9vph3yukV3pH0X8Lvx8weAlwKPzeQpAUwR\nKlgYKE1KiKIfxeSxGemOoik8ZVpCZybHIPyrkvIX+cTpc1M+n21NikrpiYLJXxpTokR+2s+iDK9w\nju/4GQU+ivbw+hhyvrJTnyS6kfJ3jxHLq0pDUfqpelC0k3PLiY1E1CHptIXPn9P2x259cvt02ibn\n95O+6B5L9lKaL9osDJ1tGPBRLCbVgaKlqpKQzaTcVE7yyZtJnVI/ezPtQ6r3rH+hzpO2mq1HOtb9\nPtsX3kBRTPo5+daaZlLfmM93xnJIS1I0Ezw+252kmyyKlNquKDvlRFum8FOfZ0ldX5RgzLTP3bp1\nx1Y7U15VmUkf5ALnn5tszY6fouiMw06+vOW26IyPTno6t5hcMh57AW1RdPr3kPR10WfQvQJ8Pn6+\nCrxwTp77AC4+Kz06wzBxMb17Dl6gj9k8UJXuqrvZ1pwTLaFTbvfBntM/iAtT5xZzPi9q0jHfwGb8\nfB2A7fia9evKlA8Ao86RZKvrW1Dv+oY7tzrHki872erm3BWHSZ2q7hyrZ/J02zCV223jUf60ndPD\nsW5NLuS069nqsJNe8wwAE12vHb7xyjbTtFyaOmvWv0QBjA8c3Z7rf2KUfdrONsKx+3JNRgesT4+F\nSb33eGoqV7dVU2lNp+3IxybtMi99uu1G8f/h/QIwIqiqzbZclzAWUhuk90l7TfpthxE7uezwfoGD\n7OaW7Z4L8Oy546fLpN0nvux0yhrNKSdzH/C5eQkr8Azw1POeu33nCnmfivmPRZ9B9yqT39Jl4Ok5\neR4DXgZ8mdloJ4QQ05SEgDv7F/PKOOeetNb+abr/zh/OM865J2/WVqLPoPtR4K8C/xR4FfArsxmc\nc/vAR3r0SQhxa3OzV7iZGEiPHUxXpbcbac65TwB71toPA41z7uN92RZCiLOC8evc3yaEEGIhvV3p\nCiGEOGM70k5r84S19juBtxNumz7mnHuLtfa/AV4LfBH4UefcwdvgJ+PLW4C/5Jz77lNsjx8B3kC4\nUfFDwFv79MNauw38OuFm+VXgdcDf68sHa+1zgfcBLyCsLa/n9cVJ9s+sD8B3MDNGY74THafz2iIe\nz+M0fu+tLWJ/TI1R59yXbpXNV2fmSveUN098EXhFHED3WGtfDnxP/P4p4D/uwwlr7QZwf/x8Ku1h\nrX0e8HLn3Cudcw8B956CH68BHo32Pw78jZ59eBJ4JfAxmN8XPfTPlA8cHKMvttbew8mP01k/psZp\n/N5rW8yO0Rhwb5nNV2cm6DJ/80QvOOced87txa9jwhriR07BlzcCvxo/n1Z7vBoorbUftNb+UrTb\ntx+fY7Ik9AphMXVvPjjn9pxz3QW28/riRPtn1oc5Y7QhXNU9clI+zPMj0h2n0HNbMDNGrbXlSfuw\nTs5S0L3CZOHx1fi9V6y1LwHuJqwh7tUXa+0AeMg596F46LTa415g6Jx7JWF1+h2n4MdngZdaaz9N\nCCz1KfjQZV5fnEr/pDHqnPvMafgwZ5xyCn7MjtHXnoIPN81ZCrqrbJ44May1dwHvIPwrfhq+vB54\nV+f7abXHVeD34ucPEbbJ9e3HG4Dfcs69EPgXwOAUfOgyry9675+ZMXqYXyfN7Dg9DT9mx+ifOQUf\nbpqzFHQ/Spi3gbB54mML8q4Va20FvBN4q3PuccIOl5f37IsFfsxa+zBheuPZnE57/AHwkvj5fsKf\n9n37YZgsVn8ivp/K2IjMG5u9jtc5YxTOwDi11v44/f92Z8foF07Bh5vmTK3Ttdb+IvDtwCedcz/e\no92/DPwD4NPx0E8B/yHwF4A/ItwVPrgB/uT8+UhcvXBa7fHzhD/rnwD+CvBzffphrb0C/BPCtvsx\n8IPA3+zLh/gn9O8QVgx8AvhpQjtM2T/J/pnjw+8Df53OGHXOfdRa+zZOcJzOawvn3KMx7SOd1Qt9\ntsVPAz9AZ4w650an9Xs5Kmcq6AohxO3OWZpeEEKI2x4FXSGE6BEFXSGE6BEFXSGE6BEFXSGE6BEF\nXXHTWGufb6393kPSnmOt/bsLzn3IWvuOk/MOrLWPWGuffZI2hDgqCrriODwfOBB0rbVl1Ar4qf5d\nOh5xH78QJ8aZknYUtxxvAv4Da+2LCI9iejdhEfv91tr/AniHc+77rLU/AXwfYT/8/+Sc++V5hVlr\nnw/8Y8IDTF8E/Ixz7h9ba38llvV/WGvfSlgQ/whhA8XnCbuT/i5BZevfA37SOffbsdi/Za19AfBV\nggRgY639OeBBwvj/Cefcx6POw+8Qdji9am0tJMQMutIVx+EXgX8e5fUcQRnsnc652aD1D51zryAo\nQb3ZWrto3N0L/CjwCoKO7yLuIeg0pF1zryeIn3S1VD8YbX8FeK219j8Cyigb+Z8Afz/mO8x3IdaK\ngq5YJ7vOuU/OOf6D1trfJ0jvPRdY9LjrzzjnRs65JyA/Ob67bdLM5gW+BPw/zrnd+PmuTp6Pd97/\nXcIV9KuttY8QhNKTSMphvguxVhR0xXEYMT1F1RyS76cJf7L/eYL6kzkkH0wH2MRTwDfEz//+IXkP\nC8wPdN7/LfAZ4L3x6vwhJrqrh/kuxFrRnK44Dn9IUJr6DeAnF+R7P/CvYv5ZQexV+F+Ad1lrfxTY\nPeK5D1lr3wx8DXhPfNTLy+KVridcAb/tJnwS4qaQ4I0QQvSIpheEEKJHFHSFEKJHFHSFEKJHFHSF\nEKJHFHSFEKJHFHSFEKJHFHSFEKJHFHSFEKJH/n9edGRDA31LgQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "X_power = np.std(X, axis=-1) # calculate the standard deviation for each electrode and each trial\n", + "plt.pcolormesh(np.log(X_power.T))\n", + "plt.xlabel(\"trial number\"); plt.ylabel(\"electrode\");\n", + "plt.colorbar(); plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### If the above picture is mostly black with a few bright spots, it indicates there are outliers. There may be bad electrodes or trials." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# electrode 1 has unnormal variations and indicate a bad electrode (electrode number start from 0)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWIAAAEGCAYAAABfOZ82AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJztvXmUZVd93/s5wx1q7OpWjxqQkIU2\nRgOKkDBgCELSAp5fAo4dIImN7RXyYjsvDsOTzQO/Fdus9Zx4iBXbOM5z1nsJCYHYOHbAxEAAR0aA\nJpsZwUZDa1aPNQ93OOfs98dv71O3qrvure6qe2911+/T66577jl7+O3h/HrXPnt/T+ScQ1EURRke\n8bANUBRF2e2oI1YURRky6ogVRVGGjDpiRVGUIaOOWFEUZcioI1YURRkyaT8SNcZcCnwSeAkwbq3N\n/Pl3Az9qrX11P/JVFEW5EOnXiHgauAO4P5wwxtSAm/qUn6IoygVLX0bE1toG0DDGdJ5+B/Ah4AMb\nxfPO+lbgeSDvh22Kolw0JMAR4CFrbfN8EzHG7AMmNxF03lo7fb75dKMvjng9xpgKcJu19t8YYzZ0\nxIgTvncQNimKctHwGuCL5xPRGLNvcmL89PzC4maCzxhjrumHMx6IIwbeDnxkE+GeB7i0eSmpS8md\n42R7HoBqXAHgqsoUT7XnAKjFYv7VyR6+3jwOQIFs2b6ytg+ASSo8ns0CUIkSAFaKFnPZMgCHqnsk\n/Sil6q8f83mOJzVGfL7HWpLnRDLCgXQUgO+sSJ77KxPsSeoAfGnaAnDt5GWsFC3JoyL/2cZETPt8\nJ1MJn7mCaiTlON2WznCgMk7DZQC0XO7tS4iIAFj26dajlJq3r1m0AWhTlOVo+HOH0glWfHrT2RIA\naRSzmMsgIvd53DpyBQDfbB7nBdW9YlO+XIZfziXfWZ/GSFJlIhkBIPG2zecNAnVvm8Ox0HEe4MX1\ngywEm33+y0WL2ypHAPhqIW22N67xdDZf1gHANckUX1p5EoDrRiT8wyvHGImrANxcOwzAl1ae4sr6\nJQAs+LLWopTjvn1jMZl96XhZtgOV8dKmRiF1llEAcLwxA8D1Y5fzVOvMezHx9l3t6+65bJHf2iNt\n+7ElOTdDTsWHf8atADAapYwgccf99zGaZV++Gulvl+Yxe3M593BVvo/TpoHY+c+vkD76e08f5BVN\nmXX8o+pCWZ7D8ZjYVUg/i4nKegn999p4km/nUvdN32cOpWMYJO6X81MAFM6V7XakMgHIvfb15jGf\ntuR/Y+0QAI9ns6S+fpYKyXMyGeGov4deOCLhFotmaUurkPT3JyPMO2mfYO+hdIwWbR6vPQ3eb5wn\nk/MLi/zH3/s1Dh3Yv2Gg4ydP8RP/+3v3IiPnC9YRG+AmY8zPANcZY37OWvu7ZwmXA6QupeIqxM4R\nOWm82H/XXJXEidmp/x6hWl7Hd96qk+5ejzrCE9IqoIjXpFGhQoW1eQU7gNKOxKVUndzwIQ05J+Fy\nP6ESu5TId6TUX0uIOtL2t6PLqfhbM1yruAqZlwDJXVTaHhxx7G+AhFX7Mq8ZUlB0lHO1LjIXrckj\nJiHyjhAfrk61DBPKs1p3cZlvKHcUJST+P8OEtel3xi1wRMXaxxFVV6XiQvar5RqjVtYBrGtvX64R\nqqUNtY62iKOO6/5cpSyHbwvS0pYo8nXr0rJsIbxzMUmwzztiV0RlfUYd5QzErPbRUP5DkaQXytUg\no0JIRxxdjQp1b/uovyVrzpH7fMciibvHxez1bTVGaLMI56Q8R9LUp1Fjrz9XcY2yPKF9Q9+Licp6\nCeUepVbeE6EPVl2VcW9/uJY7R+FW6wPkXivvV++IQ56JSzv6ZUeeRbx6jLRTOA7p11yVii93sLfq\nqvhqhG2Yxjy0fx+XHT6wcQBXbDWLrvRr1UQF+BTwUuAzwPutte/11764gRNWFEUZDkUhn27X+0i/\nHta1gTs3uKZL1xRF2VE4HK7LqNfRX5XKQU1NKIqi7FzyTD7drvcRdcSKoihFAUWXqeYLcWpCURTl\ngsIV3R/IXYgP6xRFUS4oXI+HdeqIFUVR+otzRfeHdeqIFUVR+kye93hY11/FBXXEiqIoRd7jYZ06\nYkVRlP6iD+sURVGGjHM9HtZ139BhjLke+ANku/WjiMrk/cB3gJa19vXd4qsjVhRF2fqI2FprXwVg\njPn3wH7gs9baH99M9vqqJEVRlKA10e3TBS/rEGgiWsmvM8bc699M1JUdOSIeiyrUqDLrmrS8FGFQ\nHpspGix6ScVFP39+dTpFyytZBam9aS/fOJqsylsGYiIK/6dGSDdzOU0ndXm6JTKJrTTjCi9pGOQm\nkyhiOl85w+Z5L+03VRO5wJn2Ipm3vZHKuf3pWCkdWPH/B7ZcxoqX+BtNROHqVL7EiaZIEe7zEoN7\nq3tKec6xWMLlHfvfJ7104NHmKS5Jx8tygshHNr1AVTgXExF79bGVTGwP++krUcIlPo/HWiJ72Cqy\nsq5C+EO1qTK9YMm+dIwnGxLnwIjYnruCRixlXPJxD0Q1WtGq/CXA/mSMZa84FtTIKlHMIwvPSn4j\n0hZLRYuZppd3HJXwSRSXfeBoLteaPl1JT7p6w2Us+/4z4eVMR6IKS5GEHfNqaTP5yqrMYyrlsO1n\nyjRyP0IKEqaSX9tfD+p4EX86L4peNloo6/2SSBTJToU+mk4y7ftA4tukwHHa97MTqbTFYpqwN5F6\nsU6kSGeKBlO+rX7/8csB+FZ+khvTg8CqjOpYVOFEIemFvrfi2qXNi76+n4sbpSpd4tXP9kVVWqy9\nX46uHONgbQqA2VLWslrWxbjvj6GPjsRVpnOxOdzTWZxT83FC38tdUdr30lRkKf86n2bE13NokyXX\nps32bTt2RYYr2l2v98IY8ybgV4FHgK8C1yJO+ePGmM9ba7+xUVwdESuKorgeo+FNPKyz1n7CWns9\n8AzwQ9baJf++zk8C13eLq45YURQlzBF3+3TBv+YtMA9rhus/CDzWLf6OnJpQFEUZKFsX/XmjMeY9\n/vgRIDfG/DUyNXGvtfaBbpHVESuKomxx1YS19uPAx9ed/vPNZq+OWFEURfWIFUVRhswWN3RsFXXE\niqIoF+M76xRFUS4knMtxbuOHdd2ubQfqiBVFUYoeUxOFTk0oiqL0F1VfUxRFGTJFj1UTm9jivBX6\n4oiNMZci2/peAowDLwPuBgrgIWttTxEMRVGUgTHkqYl+bXGeBu5A9DgBngRut9a+GjhojLmhT/kq\niqKcO1vc4rxV+jIittY2gIYxJvw+1nG5DfT3EaSiKMq5sJuWrxljbgQOWGsf7hbukeZJoiKhEiUs\nZiLdN1ssAiJXuZyJjOGyl1Q8WdvPQlskBeteSrIaS9G+nR1nti1xiw7ZyEYmsoOnvcRgI2+X8ojL\nbUk3JuK4T2fB27EQrdD2e9KX/DmHKyUvZ5tLZRqjFUn7ZFskEFsu59llkYjcOykSjMdac1TixJen\nUdqS+eUywc4Cx3xb0j5ezEga1QkWIrEh8X/cPLFwHES1kRMNkdIcm6hyKpM6CHWRRAlLbYm72JZ8\nvzsyLfa25njU12Oo1+WsyaivqyUffiFbYSX38pH+WjVKOb4s9o2nIoVYOMdyLnV6amVO7Bxd5KS3\nKaRxyi0wNSLpHG2dBmAqGWXJt8dxJN2rxg6x0BLbn2rJuenmPPVUJBX/qvE4AHlR8FQqZQpyjxER\nx5YkTjwmko7HooSTrbk19biQr3B8RcIVY3ITjvj054smjdzLR7pVSdQgjXk6lTp7YuUE9bEjUo5Y\n4u6nQhBl3eOlIie99CbAWLgaVVn00pw1b9NBlzATSR7j/tZ9rmjj/GaDf3SZ1Ofi8UM858/N+PZb\niSqlxOayr4vZ9lIphZpEkkeOYyVIY/r74UTS4CX4tvfSoqdXFsq+ecLXSz56pKyDZ5tS76lP95nm\nNE3fziHM8cZMaXtgrrVEo+77QyK2t1zGkpfBPNaQdGeSRYgL4jrbg+vhiC/EEfHZMMbsAz4IvHVQ\neSqKomwK53qsmrgIlq8ZY1Lgw8Bd66YpFEVRhs/FqDVhjKkAnwJeCnwG+AJwK/Drft74fdba+/qR\nt6IoyjlzMc4R+/c33bnu9K/0Iy9FUZSt02NqgotgakJRFGVHczGOiBVFUS4o1BEriqIMmTyXT7fr\nfUQdsaIoigrDK4qiDBlVX1MURRkyOkesKIoyZJzrPv2gUxOKoih9Rt/QoSiKMmQuRmF4RVGUC4rC\n4bqNenVErCiK0mf0Yd2ZvLh2kJqr4jr2d1e8LvD11QM8nonO7ojXcZ2KquytighvLZZzR1L5PRFV\nOVYZB+AFiZz7RvM48/HK2vCVKaqRaMFWvIZqw2XUvD7sjNeOzZ1j3GvLtpHGmc2WGfXnTq7MA/DG\nS67nRC56qkcSyb/AMTYp4a5NpwC4Od1Py5fztBPd3RdEIzzmvPZwJt97kjp5ZQ8Ap3I5dyidKLVd\nG07+dLph6irucpcB8H/VvgdANUrYk4wAorEcCGXfV5sE4Psrl/gyFqS+Lq4dvRSApmsz7/Vpg65z\n7gou8fW+Px0D4LvLz3PtpOSfep3l3BWlHvDh0X2SXpEx6cVkx3zdnWov8ouJlONdFbHpYDzK5IHr\nAUqd3J/ND/CBKQl3xNfJyXSOS6oS53+dEJs/tWB5QU3yW/Y6utUopZ5IuYOudCVO2VMR+w+k0lZR\nFFEbq5THAJM+zOXpJC1f3y3/J6vDMeb7yOWp2NGq5/z0Vz8g+f/8Pwbgnk8dpOHTuzQV2/Ic6r6r\nj/j7/blKyosj0QC+piF2vuzq5yhyiXvfU17nuHaQhu+Hl/wzeRHD/vcf5aaG2PW9+n4A5l2LK2Ox\n/xmvobw3GS3LsT8RfeyXu3HyVNKb8RrAI1FK4u0q76tLDAuF9NdD/tyeqMpcLmkX/uHWG5JDAHw4\nXuBIVfr8rL8v9iSjPNU4CcCL6hLueLrAQd8GQfv4+solPJGL1vL+mrT3lZW9tKI2j3KUbWE3yGAq\niqLsaJzrPv2gjlhRFKXP6NSEoijKkFGtCUVRlCGjWhOKoihDpugxR6zL1xRFUfrMFkV/jDHXA38A\n5MCjwD8Efgu4BfiKtfad3eLH52atoijKRUgYEXf7dMdaa19lrX2N//1yYNz/rhpjbu0WWR2xoii7\nHpcXuCzf+JN3HxH793QGmsAdwGf9788Br+wWXx2xoihKmJro9umBMeZNxphvAYeACjDvL80BU93i\nqiNWFEXZ+tQE1tpPWGuvB54BMmDSX5oEZrvFVUesKIoSNnR0+3TBGFPr+DkPOGR6AuBO4P5u8fuy\nasIYcynwSeAlyIR1Zoy5m00+QVQURRkoW9/i/EZjzHv88SPAPwbuNsbcC3zNWvtgt8j9Wr42jfxv\n8KcAxpib8U8QjTG/b4y51Vr7UJ/yVhRFOTe2KPpjrf048PF1pzc94OyLI7bWNoCGMSacegVnPkFU\nR6woyo4grI7odr2fDGpDxxTwuD+eA67rFjijIKZgrmiSe4m/lpdgzHE0nVRK4qe4x6PVYgTpzFVJ\ny6iUb3zKh6nEKYmXuhyPZWoniSJSLxF5tHkagImkTt3LO7b9/5Z7kxFGfH7faRyXcOkIVS/TmcSS\n7mzRLCUBD3qJwYSYZSerXJaRMiy4rJTTPF14CcEYnmzL3L7z/xNPJSPkHfKTsCp9GewCkYr8d+k0\nACNO5CXH4kpZ3hAnieLVuvVSjhVf/iSKORBLeke9HZUoKfMN3/W4wkSoP98WU5Uxqr5+Kr4NknhV\nejPIYV6ZTnK6aHib8rIev9kQmcORdE7Siyp8r3UKWJXtvL9WsLIkEol1n0cnoZ72VicY9bZkkdgc\nEZW21tNqWbeh3+yJ5JxLHEtRa005nm9M+7JGZR0EWdEoOvNxSyVKWHnfzwDw8f9xGIBH645RJ3Ee\ni0RGch9pmc4eJ+mciDJWvM21itjE44cpfLiH6nLtpGtR9Xk//y8fAOCpaIoX+3Iv4eU6Xc4Ka9sv\npyArj6X803FB08+Hhr5aieLyYVLoq6Fvd7IvqjDm+0Pot4u+DHvS0bIPBonVapRQ8xKo4R6ejdKy\nTau+T4+SlPdzSKPAlXKs28Iu2Vk3xzk8QVQURRkoQ5bBHNSqifs4hyeIiqIoA2Ub1hFvhb44YmNM\nxRjzOeClwGeQxc0N/wQx7/UEUVEUZaAU9FhH3N/s+/Wwro2MfDt5oB95KYqibBXX4+WhXV8sug2o\n+pqiKEqeQ9ZlgkCF4RVFUfrMLlk1oSiKsnNRR6woijJcnHPl2ueNrvcTdcSKoiiOHuuI+5u9OmJF\nUXY9Litw2cZr1Lpd2w7UESuKougcsaIoypBxdN+0oVMTiqIo/UU3dCiKogwbnZpQFEUZMgXdpyYu\nRK2JrZIj+qgFjijovfrv3LlS07bitX8zXKlTGmh5jdsiSql5rdOg+Rp36MnW/TXnHG2vndrymr0O\nqHldpEpH+kG7NY4kvZS41EutJ6KZu1i0yjyChmvaobEU9I1XXEbm0wtrFdsULOeiVRs0eOOOOhjx\nGq4V4rJ/TESVMu1gy3hSB2AqqjKLaOuGuihw1H06S5nklfj0Kx12Bs3XjOKMOk7jhKisA/meTOql\nXvS4T3/FtUsN2tzX8T4qtKJgvdgWx3Xmgg1Bu5aIJV8XFa/5fJJWWRch31FfVrF5laA9HMJ3ajOX\n30Tl9Ymgi+uapUbxpC9HO1/Vfw5xQ/kzt7oFNujpRsDj/0Ps+lpltU/FyPVjXn86jkfKfubNYImc\nWa8HfNr3gVqalGWbo+nLujpSe+KEvCh4oZ7R9v0h9Km0Q7e57MtO2rUz3EpUlDrVob/uISVxa8u9\nUrSYTEVnO9xr1TX3QbUsL4jud+LrKmhip1FS3kOddRYY8+6pQVH272MdZXBn0YA+X1zucFmXqYlc\nR8SKoij9pcccsU5NKIqi9BudmlAURRkuW3x36JZRR6woiqIjYkVRlOHiMnBR9+v9RB2xoii7Hp2a\nUBRFGTK93g/a53eHqiNWFEVRR6woijJsXNRjkrjLtW1AHbGiKLsenSNWFEUZMi6LKLqMel2+g0bE\nxpgxa+3S+WRkjBkFPgaMAXPAW621zfNJS1EUZTtxLsJ1c8R9nprYlGqGMeaHjDEPAQ8aY1JjzL8/\nj7zeCDxgrb0NeND/VhRFGTrhYV23Tz/ZrHzRLwKvBk5YazPgyvPI6zFkNAwwBZw+jzQURVG2Heci\nXNHls0Me1uXW2qYxJkxZJ11Dn51HgFcaY74NnADeu1HAzBVEriBzeSldGaQVc4pSXjI/y/tLgnxk\nO5UwMZBGa6UQY6IzJB2hQzrTz8wvFU0KJtaEabuilB5clViEFS/NWfMymEtFk3ZIrzN+IVt0Rr3s\nX4OclUJkIIOk4krRppFLelUvy1iPUuYLKVve8d9zkDSs+LK1XFbKGIauUycpbQ0kRLT9cea3DbU7\n6nMqqq4J38hbrKceVcr2CbaPxlXm2/NiezLm0y2oRc6XV9IZJabqbQ/CVlEU0fbHQYozJWK+LbNh\ndS8HeSJepu3rO0iMJlFc1su8l48Mcqli62qXDe0yipSxGiWlPGioxwJXSpWOsCodCpTymLAqrZpE\n8RmvXK9ECbYl/Wem1pC8iEs5yLlMZDAPxXWavhzVILtKXLZ329fj6QSahHqU8KeLBlf46w9XpTyL\nboWCEYBSYrUSxTR8ucN90CqyNX0JYJ6srPvQfyIiWv5H6FsL7WXGOqRHQZxCkNus+XAzBPlPR7OQ\n/ENfGo0qZfi2rzuRvxXGoyCDmTMWrXVVbVeUbb8dONf9gdxOeVj3J8aYjwEvMMb8J+C/nkdePwn8\nmbX2N4wxdwE/DvzH80hHURRlW3FFhKPLHHHRfURsjPkB4G5k3PWQtfbdxpg54Ks+yI9Ya6c3ir8p\nR2yt/dfGmM8A1wHftdZ+azPx1hEBwZBTwJ7zSENRFGXbKfLuqyaKHo4YeBK43VrbMMb8Z2PMDcA3\n/TOxnnR1xMaYH1lvD3CtMeZaa+2fbCaDDj4C/KEx5u1AG3jbOcZXFEXpD36OuNv1blhrj3X8bCMv\nGvp+Y8y9wJeA91lrN5zg6DUivsF/3wKMIqsdbkGWn52TI7bWzgJvOJc4iqIog8A5uk9NbHKO2Bhz\nI3DAWvuwMeZFwAzwb4G/DXxio3hdHbG19ld84h+31t7hj6NuCSqKolxouKLHDudNOGJjzD7gg8Bb\nAcKcsDHmvwF/gy5+c7PL164wxhz2xweByzYZT1EUZcdTuKjnpxvGmBT4MHCXtfaYMWbMGBOW6vwg\nsnx3Qza7auJdwH/zu+OWgXdvMp6iKMqOpyhiii7j0oKYLjMXAG8BbgV+3RgD8D7g94wxi8BR4Je6\nRd7sqokvAK/YTFhFUZQLDZkj7nIdujpia+1HgY+uO33zZvPflCP2a+R+GxhBRsTvsdbet9lMFEVR\ndjQ91hFDtPmJ3PNgs0n/K+DvWmtfikxE/2b/TFIURRksW50j3iqb9vHW2mf899P9M0dRFGXwBPW1\nbp9+stmHdff4JRgPIHPF9/TNIkVRlAGzqTniPrJZR3w3cAVwLfBpa+1Xe4RXFEW5YMhdRO42niDI\no50xIv6otfb1wNf6aYyiKMpQ6KG+1m8264i/box5E3AfsoeabkpCiqIoFxK9HsgVPRYRb5XNOuJb\n/Cds5HDA7X2xCNE8jSmIiEhjvzmlQ3q0GovZQVd2nISa16rNgm6x/++t7Tq1Y1eJ/Z8aQQ81aBHD\nWr3Z9drDck3iTng9VtcRfzSt+/wLKj7toPlai+Myv2aHDuzBeFRskP/jmC9aZRmDzm9CxIgvY7Bv\nfzxCzWv6NoJGsytY8JrMV1emAJgiYcnnW/Xf16V7ubfxDAD7KqKZe8Tr8z6f1Kn78k7ENUk3rnBl\nKoJ5T6Vi72uTg1i8VrBPd4yU6WgZgO+Lx8v6+S5zvs4kvVEXsT+q+DykrAtkXNGQclxRFz3dUWKu\nGj0EwMsrBwF4UZ6K8glQ9134uvphlr2u8tWRXIzrR7jFiQ1f9TZdGdWZr0j93JzsA+A4q1rLE15q\n+yXxHmZ9u4z6c4fqe8v0TyeNNfWZELHi8z9RiM7wC9IJDrUl7WudlGclWu1boX/USUo94sOFnDsR\nr/bNoD18iUtoeH9wJdLP5qMWdW9f6N9JFFHz/T70y/EopeXTORJJ3EaS08oWxYYOW4JW90gs/eGg\nS9jnb49Lq1InUxPXcMqXc9RrBY+7mMOJXA9tsce3TzVKyP19EzSNGy5bzaujHpN1Tq/tHGM+3J5E\n6rFge+dte74qaYc44l/ymzqAcl2xoijKRcGwR8SbXb72y+t+/x/bbIeiKMrQcJv49JNeesR/H/gH\nwI3GmKAcVOmzTYqiKAOlKOKuqyaKs7xabTvpNTXx58CXgXciW5wBWsCxDWMoiqJcYBTQ9Q14fX6J\nc/epCWvtnLX2SeAbwFP++Bjy/jlFUZSLAkfU89NPNjvefnt4zYf/fnv/TFIURRksBfI28Q0/fc5/\ns454zBhTATDGVIHx/pmkKIoyWBwRRZfPTlm+9jvA/caYrwA3IVueFUVRLgpyKNc5b3S9n2xWGP4j\nxphPA98HPKa76hRFuZjoNQ+8I+aIjTFXIy/F++fAvDHmnX21SlEUZYAUm/j0k83OEf8B8KvAmLU2\nA97UP5MURVEGi6O7E94pMphYa7/lX4p33hhjfgJZ+pYAP2atfXZLCSqKomwDw56a2KwjfsoY80+B\nCWPM/wY8fq4ZGWMuA15rrb3jXOMqiqL0kwIouvjafk9NbNYR/zTwDuQNHSnwT84jrzcAiTHm88DD\nwLustf1+GKkoitKTnKjHqokhPqwzxuwzxuwDJoA/Qh7W/aH/fa4cAqp+RLwMvPk80lAURdl2dvoc\n8X89iw0R56dHPAf8pT/+C0Tf+Kw453DOUYkSlr22bmC2aNIqROs06PLOkbGUiz7sUibh84r8MbFM\nRsOtDd92eZlG0GttuazUDW7moiGbRDFZ0PllVe+35Wtkwec5mtRoFm1Ju8jLuIGQx7LLWfRxKlW5\nvpK3OVmIpm/4P7fhMubbS2vKveIyZnLRf419yOeLpVInecrrBs+3l9jr9YK/0zoFwET1MHNO7Jvz\naTxYtEqt4+ONGQCerottJ7JFrqiOrQk/ly0zk4mm76mmaAtno0Wpxxu+0yjhREuuP1GZ9Gk0yjo4\n1poF4JaxvaUO8Ilc0o2iiK/UpRzfLSSNI/Eo35p7UsKNyLm9lXGeXj4JwE0TVwLwXHuu1KJe9HrD\nzzZnWBw5DEifAph2DZ5snFoTbk8ywnwhZX/G1+Oz+SLTmbTBddWDa+rpW7V5jrXnfLnlFspdUWpX\nX1s7AMCj7Vn+qnYEgG8jZZwgZczrBy8VUv4mBbmP+1ws/XKFgsVC7Gt7XeeZyJUjp6NO2mWmWGHM\n6zrPeb3q2aLJYix1H8I3Kco6eNqJLbO+bQHmXMuXJynbKtTniSinXUl8vUidfGX+KAfqonc9kYpG\n8PHKKE9mUi+Fz+uKilybyxskvr+V91nRLu/JoOHccnmp69z05VomK9MLbbI/HScjY7sooqir1GXR\n51cl9dKaeJ219nbgvYifmABeDzx4Hnl9GbjRH98EHD2PNBRFUbadYctgbnb52r8C3gIsWGvbwK3n\nmpG19mvAijHmHh//j881DUVRlH4w7HXEm31Yl1trT3UsXzsvcU5r7V3nE09RFKWfFFEPGcz+zkxs\n2hE/ZIz5TeCgMeb/Bu7vo02KoigDpeixamJHvDzUWvsLxpg3IlrED1tr/7yvVimKogyQC2VEjLX2\n08Cn+2iLoijKUBj2Gzo27YgVRVEuVnqtjBj2OmJFUZSLHtdjasLtlKkJRVGUi5XMf7pd74Yx5geQ\nF2YUwEPW2ncbY34e2UH8JPBTfunvWenvO6IVRVEuAFzU+9ODJ4HbrbWvRlaXvRZ4nf/9DeCHu0XW\nEbGiKLueoDXR7Xo3rLXHOn62geuAe/zvzwE/Bnxso/jqiBVF2fVs16oJY8yNwAFgtiPaHDDVLZ5O\nTSiKsuvZDq0Jr1T5QUQyeA6Y9JcmEce8IeqIFUXZ9RRR7083jDEp8GHgLj9N8RDwWn/5TnrsRt6R\nUxPhz4Sfae/jF528MHqhLdJA8AkIAAAV+klEQVR9reqqlnyQ0+s8Ltxarfm2W5X/C1KDOQXVOC2P\nA8teljDQKrJS2jD1Mo+LebOUjwzfBY5WkO4rViU064nI+AVZwdEOWc/wP2CbgtRLZgbZQ+dcKdMZ\n5DXDeYCWf7l33aXkUVHaEOqh7euiHdI4yx9WOcWa+gNYdNlZjwGyIoOkVpYNRI4vlG0kXi1rqzjz\n4XCQQMy8TTmO3JcnyJTGLuKFuYS7N5J6jONR8kLsDOnWo5RmLsdBSnI5b5ZtulR0ypiKfW3O7Cud\nUqWBdkedLGYijRnXvHyjt6MaJUwlIjUaZBmTKC7lSYO0Y8O1ebWXDv2ulz2tEvOqhoR7MJE+VSEq\n46Q+jUkvlRmuA9zRbJJEEu571TPtvart+0BclI4jlLvuHNm6PBquXfbr0BaSn9h62q224xHfLmGn\n72K7wZ6qtFvDy8a2K0Upy9r26wzec+gEAG97BprF2j6Vu6KUjQ33ZhRFZ/TLtis44OVJg+zoUtEi\nj7ZPBjOn+8qITbzB4i2ImNmve02e9wFfMMZ8EXgK+NfdIu9IR6woijJItrqhw1r7UeCj607fB/za\nZvJXR6woyq7ngtGaUBRFuVhRrQlFUZQho1oTiqIoQ8bhygfeG13vJ+qIFUXZ9eR0XxmxiVUTW0Id\nsaIoux6dI1YURRkyju4rI3SOWFEUpc8UPeaIu13bDtQRK4qy69FVE4qiKEMmY3UL+EbX+4k6YkVR\nFPo/6u3GwNXXjDHv9kIYiqIoO4JiE59+MtARsTGmBtw0yDwVRVF6MeyHdYMeEb8D+NCA81QURenK\ndgjDb4WBjYiNMRXgNmvtvzHGfKBb2OWiReIKfoUnSl3eQMNl5bmlbAWA2XqLZa8dG66F79limSWv\nARy2KeauoOk1VCdT0ZVt5C1WgrZtJuHHKqv/T83lklezaBN5UdYFn381Tktd1UYmaTjniL0G73wh\ntj2dJyy1JU74U2elaLPg0y61WV1RauUGOxfcqn2hHO2oUurIVoNecrtBuy7pTLcXpJ5qGXO52NDw\nmr7VOC3tX2iteDsl/aW8wUwqdRB0lnNXrGr7ejtXXLvUjg35A2T++pLPa6FoljYHjd8ZMubc2vJM\nZ8s8U3PeJgk3nbTK9IIG8VhcLfNa8OVazhqQ1qWMuZQ7c1mpl7zstXVzV5Rt0KyMl+mcas0DcInX\nGW4WWdm/ThYrZdxQT7P58ppzK3mrbO/Tvj1bRca3ff865Ft8j4uZ8Xfd5Ym8wOGwC4rDMOLfUplF\ncCQZL+MAPJrUafouWY0ape2zTtqqQPKqRAnLPsGgYd1gtS5CW2VFThJLgkGXOyE648FUg5xTXju5\nlfs+WhQs+7oP/WLatUpt6dBXf/n4YanD1hOlPvdcWzScR7y+NaxqcZ9uLYBv3nHfzjkFp30Zjzdn\nABhN67h4+/a77aYNHW8HPjLA/BRFUTZFgSvF6Te63k8GOTVhgJ81xnwauM4Y83MDzFtRFGVDZETs\nunz6y8BGxNba94ZjY8wXrbW/O6i8FUVRurErN3RYa189jHwVRVHOxrBXTeiGDkVRdj2O7g/kLsoR\nsaIoyk4i7/Gwrtu17UAdsaIoux7n/3W73k/UESuKsuvZTeuIFUVRdiQOR+G6jIgjHREriqL0lV25\nfE1RFGUnocvXFEVRhoyumlAURRkyrseIWFdNKIqi9BmZI+7miPvLjnTErSIjdo5/ll7NL7W+BqxK\nIOauOKPCVoo2K14u8ndHbgbgt3gOENnHINcXiIlK+cJSxrBYldKMvJxhRMS0lzQMcoJLWZPCrZXf\na+YtKrFI/DWydnm+lmRr4s4XDZbaIue35OUCT7bmS8nLTjtHklWpR4DFvEnsxRLnvDzjVDLKbLZU\nlgmglWflcaizBdcuZTyDZGGW57S8TGUo73QuaTXyNieyxbIO1hPqPyty5jKRg6xG0pWaxao05vPt\nWUDqOJxb9uU/7VqlTdO+DDOtBRp1Sfu0l6UciavU0oqvT/leLFalMYPE6XK7WZYjlLsap5zMRBJz\noS127q1OlCOfIPEZu7h8Yn7K29LqkFs94eVEX7PvxQBcFY/zuvgSAO7w4X8tzXmFE1nLI1659Z56\nkx++4Wmp01MiI1nfl3PqMZGr/JHrpPzf+XJMPZVI4+NSniRx3HpyytssbbFnb4NHTu4FYL4ukp8L\n6QRXRCMA3FibA+ATf/wzfOjNfwzA34kOAvBYkmELqdOXR5LG16spY5HUaaiTlxZ1Kk7S/lwq/XKK\nlLqX56zF0s6jlVUJy9B/54tWmU7oZ2HZV7NoledC/1xoL1P16c37vt8sWiKF2YFzjjyRMjZ827Zd\nThQ7L/y5dXT5mqIoypBxzpUDprNe16kJRVGU/qKrJhRFUYbMVoXhjTGXAp8EXgKMA5cDDwDfAVrW\n2td3i6+OWFGUXc82jIingTuAP+0491lr7Y9vJv9BvzxUURRlxxHmiLt9umGtbVhrZ9adfp0x5l5j\nzLt75a+OWFGUXU+xic858jxwLfA64E5jzI3dAuvUhKIoSg8ZzHNdSWytbQJNAGPMJ4HrgW9sFF5H\nxIqi7Hpy58r17mf/nJsjNsZMdPz8QeCxbuF1RKwoyq5nqw/rjDEV4FPAS4HPAF8wxrwJGRXfa619\noFt8dcSKoux6tvqGDmttG7hz3elf2Wz+6ogVRdn1ONdDGF43dCiKovSXXSMMb4z5AeBuZCXIQ9ba\nnmvrFEVRBkGBIxriFudBrpp4ErjdWvtq4KAx5oYB5q0oirIh3VdMrKoH9ouBjYittcc6fraBfKOw\niqIog2TXCcP7HSYHrLUPbxTmRdVLqFPl026BwzXRTl1MRa/0luohvhGdAmCyLlqrL4jHmBk9BMC/\nS04D8MJYtFzjdC+PZbLz8HAyDsBs0Sg1cPenYwBcXtlTbmN8si3hR+Iq11SmfBzRO35hdV+pp/qc\n17pNo5jxWPRZmz7cbRMv4lguGrhXJJLHGCn1fVLlL0NsuWbsWp6Tdd80vF7rFdEIf9F8Ruz3Wq83\nVA9wqpA6OOa1io+kE0wmoh17dSzpNabaHErlOIi1Xh2Ps6cqcZ7NxWaHI6lK3Z7yerevrV4GwDeT\nWQ4noz68XBtNaqvazamUMY0TbhiVON/n8/+r9klumbxayus1mnPnaHjt35FY7LgmGuWA17Rdrogd\nz1aXeMeEtO0zXCPlLkaYH19ZE/fNyRFOjUk5rh85AsB34pQ9qdj8ssoBAB5sHefmqvSRmZrYvDeq\nluns9Rq3sKqpe1tF0nvcLXMolaWgl/t+88nZb0t9Tt3AZ/MTAPwZXtM4i/jLROrnhVVp7yfzBUb/\n3t8E4P73PApA++iqwvP/99fSt16YJoz6AVc8J3kWETzutZlf1RA7p6IVLh2Rfvs/ff9JXURlnWb0\nO3/0P/O2tpz70Ij0qYbLuT7eA8BTvr/N5g0akdj/okS0lGdwPBXL9VNtyevqyhjf1xIDv1aT/nZJ\nfbKs78MViXsbe3kolf7ynP9j+5WZhH+kfoApX9/T/r7Ym4xw3N9Dr6leCsDXk2lemEp6+5D+M+5i\nTng7p0f3A3BVZS+tqMV3eZztYFcJwxtj9gEfBN46yHwVRVG6sWv0iI0xKfBh4K510xSKoihDpeix\nbuJielj3FuBW4NeNMfcYY145wLwVRVE2ZKvqa1tlkA/rPgp8dFD5KYqibJbCL2DbCNfnt9bphg5F\nUXY9BY6o2xxxdJGMiBVFUXYqrscc8UXzsE5RFGWnUji6j4j7vH5NHbGiKLseHREriqIMmcIVPVR/\nir4uMVNHrCjKrqdwrvv8g3PqiBVFUfpJr6mJfm9yVkesKMqux21iaqKfqCNWFGXX00v0Z+OtHtuD\nOmJFUXY9so25a4C+5r8jHfG0a1BxOSezRRZzkX483ZwH4DvVGeZykUUMUpYT9Qqz/vj55rQkMiqy\negtFk/lMZPcyLzPZKrLyeKVoAzDnGiznIv93rCEymEkcl1KJFT9V/0R7hjRKAGjkIq1Y4Hi+mAWg\nmUt6Dy4/zUQqcZe9NOYV6STfnH8SgNfvFdnOh4sFnvdSgLnfRvlUlPLMishBprHkdV11PwuutSbf\npaRNy5cj/OH0vflnGdkj8oHTbS9hWY054UMs+PpMo6QUMplpSf5P1KQOF4smy14qMuyxn20vlbKH\nofx7k1HaPt3TSLn3p2N8c/lZAF41dhUAJ91KKTMZ2nPCxZyKJM6JQtpzLl/hL6cvB+Bxnpc6q9Sw\nsyIJun9U5BE/PzbCUiZxjrakvR2ubL/7nWhKtYuMGV9nC07yalJwrCltNT4iMpxJFJUSn8eR8G0K\nHvbluHzClHkALJOXbdX29d/wbQwwE0ndzecN7nuPBeAbdTl3WQa1QtJ5zkldXMkoK2HI5b+ncsfz\nkaT5pJcLPTZ9kDn/xGg2kbI+1p6hXpX2+G5D6ufpygyn472+3I3SrqNO7oMaEn6xaFJNxAUcd5Je\nEa2WM7RZp8zmjA93dO4YL5gUidElX+/Xjk5QdRKn6vvI46nU08zKEpmvs1MtuZebabuUeV3w8uRL\nRYvnCrHzEi/b+UTU5CondbCYNXxZpimiHGpsC7krcG7jcW+3NcbbwY50xIqiKIOk6DEiVkesKIrS\nZ1zPLRvqiBVFUfpMjxGxOmJFUZT+UtB/8fduqCNWFGXX02vVhFzr3yI2dcSKoux68qKg6LJnQxZU\nJH3LXx2xoii7HtfjHRy6oUNRFKXPbG5qon+oI1YUZddTOEfRdR1xf/NXR6woyq6n1zrifq+nUEes\nKMquJy8cedchcX9nidURK4qibGJvXTeMMZcCnwReAoxbazNjzN3ALcBXrLXv7Ba/n6LzZ2CMudsY\nc68x5rcHma+iKEo3HM4/sNvg09tJTwN3APcDGGNuRhzya4CqMebWbpEHNiLuNMwY8/vGmFuttQ+t\nC5YAZFEGQBHluFgWlXgRMrIoE9UlwEVyrRm1ynBBnCPzyl55lOHivEyvTNcfh7xy8vJ6FLvyu+3T\nCS9KKaKc3M/cl+nhVvP3cV1cUIS0/V81zahVlmMhEsWqFm3yUF7f2Hm01gaABq2yTCHfdtQm86pV\nK0h6SUKZHt6mJZo0vZJXiJtHq50r5NHqqLOWDx/qx8XFGelmURvn1bMafo1lK8rK642gZBa1icMC\nIH9tIWqy0nE95Dvr6yUnK21PkmiNne2oXR6XfSHOKXwehc+qiPKy3G2fHlFR2tAu6zMq2yrY3OrI\nI5xb2xahj/r6ifKyfkKeeZRxOvbqb77PzEZQ9cch3HyUnDEicpFjxSvazfnyVIEFfxzqroiy0r4Z\n/+dzO2ozE/qXLyN0roJNyvoJddD0cVeABmvvjSWaTPt8Q/g0jcp6DH1/mSZNf70zroTJz+g/RZRL\ne7Bax3mUlXmEuA0yFn3+riNuaHu2YYFvHMu90+16N6y1DaBhjAmnXgF81h9/DnglsN7flUSu3+sy\nPMaYfwKcstb+kTHmR4HLrLW/sy7Mq4F7B2KQoigXC6+x1n7xfCIaY/YBjwJ7NxF8BrjGWjvdJb17\ngDuBX0CmJD5tjLkTeJW19gMbxRvkHPEU8Lg/ngOuO0uYh4DXAM8D+VmuK4qiBBLgCF1Gmr2w1k4b\nY64BJjcRfL6bE17HXEeak8Bst8CDdMQ9DbPWNoHz+p9NUZRdyWNbTcA718062M1yH/DTwB8hI+T/\n0C3wIB/W3YdMZoMYdv8A81YURekbxpiKMeZzwEuBzwAVZM74XiC31j7YLf7A5ogB/GqJm4GvWWt/\nbmAZK4qi7GAG6ogVRVGUMxnoOmJFURTlTHbUzrpz2Ymyzfn+AHA3ItT/kLX23caYnwfeDDwJ/JS1\ntt0tjW205d3Aj1prXz3E+vgJ4CeRp9I/Btw1SDuMMaPAx4Ax5CHvW4F/OSgbNrtLqp/ts94G4GWs\n66M+XF/76dnqwp8v+6n/PbC68O2xpo9aa58d1v2yHeyYEfG57kTZZp4Ebved6qAx5rXA6/zvbwA/\nPAgjjDE14CZ/PJT6MMZcBrzWWnuHtfY24NAQ7Hgj8IDP/0Hg/xywDT13SQ2gfdbYwJl99AZjzEH6\n30/X27Gmn/rfA62L9X3UO+Fh+o8ts2McMWffiTIQrLXH/M4YgDayxvmeIdjyDuBD/nhY9fEGIDHG\nfN4Y87s+30Hb8RgyGgZZf+4GaYO1tmGtnek4dba26Gv7rLfhLH00R0Z/9/TLhrPZ4enspzDgumBd\nHzXGJP22od/sJEc8Bcz74zn/e6AYY24EDiBrnAdqizGmAtxmrf0Lf2pY9XEIqFpr7wCWgT1DsOMR\n4JXGmG8jziYbgg2dnK0thtI+oY9aax8ehg1n6acMwY71ffTNQ7BhW9lJjvicdqJsN36r4weR/+2H\nYcvbgY90/B5WfcwBf+mP/wJ5S8yg7fhJ4M+stdcB/x1Zkzm0vsHZ22Lg7bOuj25kV79Z30+HYcf6\nPvr9Q7BhW9lJjnhoGz6MMSnwYeAua+0xZMvkawdsiwF+1hjzaWRqZD/DqY8vAzf645uQaYFB2xGx\nutPplP8e5mags/XNgfbXs/RR2AH91Bjzcwz+3l3fR48OwYZtZUetIx7Whg9jzN8Hfgf4tj/1PuBv\nAn8beAp5Gt0aoD1f9KsmhlUfv4lMCZwC/gHwG4O0wxgzBfwhUEPmQ98G/NKgbPB/fn8KWanwFeD9\nSD2syb+f7XMWG74A/FM6+qi19j5jzHvpYz89W11Yax/w177YsWpikHXxfuAtdPRRa23rQt4wtqMc\nsaIoym5kJ01NKIqi7ErUESuKogwZdcSKoihDRh2xoijKkFFHrCiKMmTUESvnjTHmKmPMD21w7bAx\n5l90iXubMeaD/bNO3h9mjNnfzzwUZTtQR6xshauAMxyxMSbx2gjvG7xJW8PrFijKQNlRMpjKBcc7\ngZcbY65H3s/1J8jC+5uMMf8I+KC19m8ZY94D/C1k//+/tdb+wdkSM8ZcBfwX5CWz1wO/aq39L8aY\n/+DT+itjzF3IIv57kE0fjyO7rP4Foj72YuAXrLV/7pP9ZWPMS4ATiFxiboz5DeBWpP+/x1r7oNe1\n+BSyU+vObashRdkEOiJWtsJvA3/qpQgtopj2YWvtekf2/1hrb0cUst5ljOnW7w4BPwXcjuggd+Mg\noksRdv+9HRGA6dSi/bzP+zjwZmPM/wIkXmLz7wC/5sNtZLui9B11xMp2smKt/dpZzr/NGPMFRKbw\nUmBvlzQetta2rLWngLo/17n9M1ofFngW+J61dsUf7+sI82DH94uQkfYbjDH3IOLzQShmI9sVpe+o\nI1a2Qou101v5BuHej/y5/3pEFSvaIBysdbqBGeByf/w3Ngi7kbO+peP7UeBh4BN+FH8bq7q1G9mu\nKH1H54iVrfAtRIHrj4Ff6BLuM8CXfPj1IuOb4f8FPmKM+Slg5Rzj3maMeRdwEvi4f83Oa/yI2CEj\n5feeh02Ksm2o6I+iKMqQ0akJRVGUIaOOWFEUZcioI1YURRky6ogVRVGGjDpiRVGUIaOOWFEUZcio\nI1YURRky/z83GioaiNR9FAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "electrodes2del = []\n", + "# electrodes2del = [9, 11, 6]\n", + "# electrodes2del = [7], subject 3\n", + "electrodes2del = [1] # subject 6\n", + "if len(electrodes2del) != 0: \n", + " X = np.delete(X, electrodes2del, axis=1)\n", + " X_power = np.std(X, axis=-1) # calculate the standard deviation for each electrode and each trial\n", + " plt.pcolormesh(X_power.T)\n", + " plt.xlabel(\"trial number\"); plt.ylabel(\"electrode\");\n", + " plt.colorbar(); plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## separate EEG data into different frequency bands\n", + "\n", + "7 to 30 Hz band pass $\\mu$ & $\\beta$\n", + "\n", + "band pass the filter and downsampling the signal to 100 Hz" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "low_freq = 7.0\n", + "high_freq = 30.0\n", + "down_ratio = 5.0" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Setting up band-pass filter from 7 - 30 Hz\n", + "l_trans_bandwidth chosen to be 2.0 Hz\n", + "h_trans_bandwidth chosen to be 7.5 Hz\n", + "Filter length of 845 samples (1.650 sec) selected\n" + ] + } + ], + "source": [ + "filtered_X = filter.filter_data(np.asarray(X, dtype=np.float64), sr, \n", + " l_freq=low_freq, h_freq=high_freq, fir_design=\"firwin\")\n", + "filtered_X = filter.resample(filtered_X, down=down_ratio)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## use CSP & LDA for decoding overall view\n", + "\n", + "use the EEG response from 3.5 to 4.5 seconds for training" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "st_idx, et_idx = int(signal_st * sr / down_ratio), int(signal_et * sr / down_ratio) # only use data after 3.5 seconds\n", + "csp = CSP(n_components=6, reg=0.001, log=None, norm_trace=False, transform_into=\"csp_space\")\n", + "csp = CSP(n_components=6, reg=0.001, log=None, norm_trace=False, transform_into=\"average_power\")\n", + "tf_X = csp.fit_transform(filtered_X[:, :, st_idx:et_idx], y)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "clf = LinearDiscriminantAnalysis()\n", + "mappings = clf.fit_transform(tf_X, y)\n", + "est_y = clf.predict(tf_X)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false, + "scrolled": true + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEDCAYAAADUT6SnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xu0XGWZ5/HvIQkEIQeCmISAYwDN\nw5AAaQFR2ki4BVRA7W5ox6xuMHS3w1wUWNrqDNPE7izvDbY6Oq1LxNWCaFw6BFQuEolBCGSkgU7S\nPOEmCiEkSEi4BXI580dVHetU6rKr9u3de/8+a7Fyqk5R+z3v3u+z3/seGhkZQUREym2PvBMgIiLp\nU7AXEakABXsRkQpQsBcRqQAFexGRChifdwLaMbO9gOOBp4CdOSdHRKQIxgEHAavc/ZXWXwYZ7KkF\n+hV5J0JEpIDmAne0vhlqsH8K4JprrmHatGl5p0VEJHgbNmxgwYIFUI+frUIN9jsBpk2bxiGHHJJ3\nWkREiqRt17cGaEVEKkDBXkSkAhTsRUQqQMFeRKQCQh2gFRGRuk9/+tOsXr2aI488kssuu2yg71DN\nXkQkRZu3buPXDz7N5q3bBvr/16xZw0svvcS1117L9u3beeCBBwb6HtXsRURSsnTFI3z7hjXs2DnC\n+HFDfPDsWZwz9/C+vuO+++7jxBNPBODEE0/kvvvu4+ijj+47LarZi4ik4Nmt20YDPcCOnSN8+4Y1\nfdfwn3/+efbdd18AJk2axNatWwdKj4K9iEgKHlu/ZTTQN+zYOcKj67f09T2TJk3ihRdeAOCFF15g\neHh4oPQo2IuIpOCw6fsxftzQmPfGjxvisOn79fU9c+bMYeXKlQDceeedzJkzZ6D0KNiLiKRg8vBE\nPnj2rNGAP37cEAvPns3k4Yl9fc+sWbPYc889+cAHPsC4ceMG6q8HDdCKiKTmnLmHM/eYg3l0/RYO\nm75f34G+YdDpls1yCfZm9pfA+dT2X17g7k/mkQ4RkbRNHp7IsQMG+SRlHuzN7GDgJHc/Netji4hU\nVR41+zOAcWZ2G7AWuNjd9TQqEZEU5TFAOxXYs16zfwl4Tw5pEBGplDyC/RZgef3nZcB/zCENIiKV\nkkewvxNozB2aAzyWQxpERArh6aef5n3vex9HHXUUO3bsGPh7Mg/27n4f8LKZ3U7tweI/zDoNIiJZ\nWbNxHWs2rhv4/99///25+uqrB15M1ZDL1Et3/2gexxURydqS1TcCMOuUSwf6//faay/22muv2OnQ\nCloRkRSs2biORcuuYO2mh1i76SEWLbsiVg0/LgV7EZEUzJoyk3NnnzX6+tzZZzFryszc0qPtEkRE\nUrJm4zouP/mS0Z/zDPaq2YuIpOS8em1+1pSZnNdUy+/H9u3bueCCC3jwwQe58MILuf/++wf6HtXs\nRUQCNmHCBK6++urY36OavYhIBSjYi4hUgIK9iEgFKNiLiFSAgr2ISAUo2IuIVICCvYhIBSjYi4hU\ngIK9iEgFKNiLiFSAgr2ISAUo2IuIVICCvYhIBSjYi4hUgIK9iEgFKNiLiFSAgr2ISAUo2IuIVICC\nvYhIBSjYi4hUgIK9iEgFKNiLiFSAgn1O1mxcx5qN6/JOhohURC7B3swuMbM78jh2KJasvpElq2/M\nOxkiUhGZB3sz2wuYk/VxQ7Fm4zoWLbuCtZseYu2mh1i07ArV8EUkdXnU7C8EvpPDcYMwa8pMzp19\n1ujrc2efxawpM3NMkYhUQabB3swmAPPcfVmWxw3Nmo3ruPzkS7j85EtUqxeRTIzP+Hh/AVyb8TGD\nc15TzV61ehHJQtbdOAZcZGY3AbPM7L9nfHwRkUrKtGbv7h9v/Gxmd7j7V7I8vohIVeU2z97d357X\nsdOiufMiEiotqkqQ5s6LSKgU7BOgufMiEjoF+wRo7ryIhC7rqZel1Zg73/hZwV5EQqJgnxDNnReR\nkKkbR0SkAhTsRUQqQMFeJDBaryFpULAXCYzWa0gaFOxFAqH1GpImBXuRQGi9hqRJUy9FAqL1Gvlr\ntKbKlvcK9iIpGDRgaL1G/hrjJbNOuTTnlCRL3TgiKdAga/GUfcxEwV4kQWUPGGVW9jETBXuJRXPC\nxyp7wCi7Mj8fumefvZmd6+5Lml6/291/km6ypCiy7N8sysBZ2oOsRcmHIirzmEmUAdqLgCVNrxcC\nCvYVt2bjOpasvpG1mx4CYNGyK1KvxRZl4CztgFGUfJCwdAz2ZnYR8F+AQ83sAWAIGAHuyyhtErBZ\nU2bC7LP41C+uBNLtrsjjxhKiLPJBrYby6hjs3f3rwNfN7MPu/uUM0xQsFYSxspoTnuWNJWRZ5INa\nDeUVpRvne2b2X4EDqNXucfe/TzVVgVJBGCvL/k0tNqpJKx/Ueiq/KMH+BuD7wNqU0xIsFYT8lXng\nrB9p5YNaT+UXJdhvdPcrU09JQtLoasmiIKiLSPJW5dZTFcpflGC/xcwWA/9KbYAWd/9RqqmKIa2u\nlrQLgrqIwlKFwt+qXauhn3wocp5VofxFCfYP1/+dXf93BAgu2Kfd1ZJW87lTutM4lkRXhcIfRT/5\nUMQ8q1IXbc8VtO7+KeB7wB3APwBfSztRgyjqysVO6dbeKvnQdgc1/eRDkfOsqHFjEFFW0P4v4Gjg\ncOBY4BrgjJTTNZCi9jk2p/sXj901UE0j1CZ0qOnqpIgDlXmPUxUxz5oVNW70K0o3zmnufpKZ/cLd\nR8xsQuqpGlBRZ2y0pnvNxnV9F5xQm9ChpqubohX+EMapipZnzYoaN/o1NDIy0vUDZrYcOJPaFgnv\nBG5y95MHPaCZnQBcCewCVrn7JW0+MwN47LbbbuOQQw4Z9FCpS6vW+oPVN44ZIGu+GNulobklcOTr\n3hREzSqEdBWtVdGvEPJYwvHEE09w6qmnAhzq7r9p/X2UXS8XAcupdeUsA+IuqHocOMXd3w5MMbOj\nYn5fbtLqVz+vXmBnTZnZNdBDtD7HPHamDKEvtOzjHp3yWDuRSjs9u3Hc/RfAW8zsQHd/Ju4B3X1D\n08vtwM6435m10EbwezWh8+pKyatpH9r5SVO7PG4+362tm7K3dqSzKN04RwAXMna7hIVxD2xmRwOf\ncfd3t/ndDALvxmnuV7/85EuCLDxJNfOLGCDanZ8i/h39aHe+X9z+MvtM2JtF9Rv9omVX1P4NcAxF\nN6Z4enXjRNobh1rXzW+TSpSZHQB8FTgvqe/MWhEGpJKaJVGWQdYi/h39aD3fL25/mcefewKAj920\nGIaGRl+H2NppPT9lP19Zi1Kz/7G7vy+pA5rZeGApsMjd7+nwmRkEXrPP2qC1nH4Ge9sdswwDgFH+\njrLUIpvP9y8eu4tf/mYlwOiNL8nWaFJ51np+3rDfwWNuTEe+7k0cf/AxzJj8+sKfnzQlUbPfaWbf\nZex2CVfESNO5wPHA580M4JPufleM76uEQWs5caaVFX3+dEOUv6MsT9xqPt+trRsg0dZoUnnWen4u\neHOtwd98vpasvpFVT94fZC2/KBWFKDX781vfc/fvpJYiVLNvlnftOk7LICSd/o488jfkfvMo0siz\n1vMDteD5m82/4yfrlvHMS88mdqyk9TqfWd0MetXsewZ7GB1MnQmsc/cHkk5km+PNIIBgH8oduwiD\nwUUWJX+TuBbyvnEnoXlKZxZ51vieEK//qOczq5t77G4cM7scOA5YCVxgZr9298uTTmgeel2MoQwQ\nxR0M3rx1G4+u38Jh0/dj8vDENJIYhEGDS5T8TeJaKEO3WCMfjpwyM5M8a3x/iJMhep3P0KYAR+mz\nP72+AAoAM/sVUIpg3+liDO0kxel3X7riEb59wxp27Bxh/LghPnj2LM6Ze/iYz4TSgokrjXGNpK+F\ndoGrCPnfmg/A6MK/tPMs5O0Mut2IQru5R1lBu4eZTQMws4Mi/j9B67RLX+O/EFZ/JuHZrdtGAz3A\njp0jfPuGNWzeum3M54q+0jTNXReTvhbarY4uQv73kw9lKT9R9Frt3rgZXH7yJW2vySxXO0ep2X8M\nuN7M9gZeAj6abpLS1+mO2+hba6w8DLHp2I/H1m8ZDfQNO3aO8Oj6LRw7PDG4FsygBq1BRa1Rp3Ut\nFC3/q7IxWpJ6tUqy7CqONECbtTgDtFH7p1vnIz/z4u8LPXDWzuat21i4+JYxAX/8uCGuumz+aN6E\nOvjVr0FmDSU5cBZnvCCrgU7JRpTzlcZgfRIDtCcBnwcmAtuAj7v77QOnKEX99E8nsa1w6CYPT+SD\nZ88akycLz5495iaYdA0sr8DUT79uGjXqQWtoWQ0OS3ainK88+vOjdONcAbzb3TeY2VTgp9QeYhKU\nTv3TUw55mUn77Dmakd0GZcvY7Dxn7uHMPebgjq2duM8dbdVvYMrj5pBkQYt748hycFjS1e/5yjrm\nRAn2v2vsVOnuT5vZE6mmaECd+qeXrLmRSa/ZE+qr8DqdiJBH/OOaPDyRY/uYcjlITXLQwFT0HTnT\nrKGFNpsjrjy7o7I4dr/nK+uYEyXY72Nmq4B7gT8CNpnZlwHc/cNpJq4fh03fj/HjhkYD/h6Tfs+E\nQx7hdy8+Cy/WgsrxBx8zGoyKXnDSEKcm2e+FnnetNcmClmYNrUwtztYbe5bBP6tKRcjnK8p2CSd1\n+p27L088RQw+QNvaZ/+u+cPc+vvvA4xOfSrD0v80xRmw7XeQtCyDw9Jdp8HIRgBOc2VpGVYtR5XU\ndgknMHY/+58mm8zdjjeDBGbj3Prbnyu49ynLvXDaHasoM0+Kks5QNN/Yz5/zZ6x68v5UA3Dz+alK\npSKJ2Tg/AZ4B1tffGqE2SBuk5v7pMvfDpyXLPGt3rKLMPCnCQHRIWrs3zk15LOJ79y/l5Vd2cNk7\nLg66ayVLUfrs93b33Xa+FEnSIH34eQTQog1Eh6L1xv6D1TemtlDt679awsZXa/NI/voH/8D8Gacx\na/bM0WNXVZStD75vZh8xs5PM7B1m9o7UUyWVM8gS+zy2Geg3nWlu5VBkvbYZGNRBE/8D6/9t+ujr\n7U8czs9ueX63LULylOUWCc2i1OzPBl4F3lh/PQL8MrUUSWxF7TKI2tzOeyZPP90CZZs+GbrH1m9h\nZN/f88q/Hw/AHsPPsuP5145uETKIpMtTXq28KMF+grtrZLNAQu4y6FZw2q1qbvfZvANov+Ma6jPO\nzmHT94MNM9lVn4K96/nXMn7cUO39ASW5VXOelZQowf5xM/srxj6W8N5UU1USWdew876Youin4HT7\nbJECaMgTBYraCuwkyhYhUSVdnvKupESq2QN/XP8PagF/YWopKpGsa9h5X0zd9FNwonw25ABaJCG3\nAgfVa4uQqNIoT3lWUiLvemlmB7r7Mymnp3GsGQTwWMJB5bmQI+RnxvYz37kqc6PzUqXFRnGEXJ5a\nJTHP/izgc8CTZnYw8El3X5p0QsuktUZw5qHzC7kNQNK0H3o4Qm4FhiTk8tSvKN04/xM4wd1fMLNJ\nwK2Agn0PP1q1kh3+FnbuGuEL62/m/OP23m275bwUYRvidp+tyrN0s1KmG2rZxh7SECXYA7zY8q90\n8ezWbdy7fP/RTdleff61fPupNcw95uAgglQR+2mjPKtA+lOmWmsRr+msRVlU9S/AXWb2z8Cd9dfS\nRbfHAeapqAt8oj5LV6onj2s6zUVRaX53z2Dv7l+jtrDqW8DZ7v6/U0lJiTS2W24Wd65vEgZZpRqC\nUG+eMtbmrdv49YNPZ3oTHvSajhNU01y5neZ3Rxmg/QTwOXffZGZDZvYJd/9sKqkpiSTn+iYtiX7a\nrPtHW59VAGHcPOUP8uxmG+SazvIBPXl/d0OUPvvTG8Hd3UfM7HRAwb6HpOb6Ji2Jftqs+0dDvnlK\n5262rMaosnr+cNGfShYl2E80s/3d/Tkz2x/YO9EUlFi/jwMMXZ4rdEO9eUr3brbQrv+4QbXITyWL\nEuwvB+4ws6eAqcAliaZABlL0B3UPomw3z7IoWjdbnKCa5gymtGdH9Qz27v5zMzsKOBB4xt2jLbnt\nwsyuBI4D7nX3j8T9vioq+oO6pTyK1s1Wpimn/Yg0z74e4DclcUAzezOwr7vPNbOvm9nx7r4qie+u\ngrw3O6tqQZGaTgvbsu5m0yKq/kVdVJWkt1JbhQvwc+BtQGbBvuirMPPuSpHq6jXjJstuNi2i6l9f\nwd7M9gH+xN3jLKzaH3i0/vMWYFaM7+pLWVZhqitFspb3jJuGvFu2RRZlnv2ewFnAucABwB0xj7kF\nGK7/PAw8F/P7IgnlYk1CmbpSit7SqopQZtyoZTu4jsHezOYD/wl4A7WNz6a6+ykJHPMu4EPAD4DT\ngKsT+M6eQrlY5Q/SbmnpRpKckGbcqGU7mG7bJSylVgt/r7t/CUhkDXT9KVfbzGwFsNPd70nie3sJ\ndQuDqkp7v5ulKx5h4eJbWPTNlSxcfAtLVzySyPdmIY9tB3ppzLhplKHx44Z4/+nGo+u3ZJ7OtB5W\nXnbdunEOAv4MuM7Mngemmtle7v5K3IPmMd2y3fSwxsV6WP33kp00W1pRu+xCrPmHPK7UPOPm4d89\nx3W3epDplPY6Bnt33wx8E/immR0EvB+41cxedffTskpgknSxpi9qAE2zWyDKjSTEoFqEcaXJwxM5\nFFh81d1Bp1N2F2WLY4Ad7n6lu78D+Js0E5S2ycMTOXT6fqOBHrRlblL66Tpp1y0QdyFOo/vjgOGJ\nXbvsQt0yuSi7exYlnTJWtwHaYeArwInAM2Z2ILXB1f+WUdpSo8Ha5A1SK01yIU5rTf34I6exau2G\ntis6Qz3/IQ2CdlOUdMpY3frs/x5Y5e7nN94ws4uAxcCH005YmnSxJm/QAJrEQpx2N5pVazdwxcUn\n8ezWbbvdSEI9/yFsOxClGy6EdEr/ugX7o9z94uY33P3rZrYs5TSlThdrd4MMXOYZQDvdaJ7duo1j\nj5i62+dDPv957u7ZzziGdiEtnm7BvtOGZ7vSSEjWtJdHe4MOXOYZQAe50YQcrPLY3XOQbjjtQlos\n3YL9EWb25Zb3hgBLMT2Z0l4eY8WdDZJXAB30RtN8/nu1ZkKcppmkUMcxJDndgv2CDu//MI2ElFWR\n9vJIosDnVduLc6Pp1ZoJcZpm0kIdx5DkdJtnv7zd+2b2NaDt72R3RdrLo+gFfpAbTa/WTBHmvich\n5HEMScYgWxy/MfFUlFxR9vKoYoHv1ZqpUvdGyOMYEl8e+9lXTpF2qaxage/Vmil6a6dfeQ66ln1c\nJG/dFlXdwO4zcoaAY1JNkeSuV4EvU6Hs1ZoJsbVTpvxvqMK4SN661ewLv1K2X2UsRO3E+TvjFsoQ\n87hXayak1k4Zg2JVxkXy1m2A9vEsE5K3MhaiduL8ne0K5VVLV7PfPnty9Btf17NghpzHvVozIcwp\nzyMoZnFzrtK4SJ6iboRWaqFujJW0uH9nu0K5cxd88Zp7e258VpU8TlPWG5Bl9UwAPWsiGwr2ZFOI\nQnggRdy/s12hbP6ebsFbOyXGl2VQTOLmHPWa77QDKpB7mSkTzcYh/RkXoXRfxP07WwcrW3Vrekc9\ndoh9+qHIcrA4btdKu2u+n3GRFfc/ycLFt3QsM7pO+qdgT/RCNMgFFtLgUxLBolEo7394E1/63r+y\nc1e0G0eUY4dyUwxZVoPFcSoG7a75b12/mquWrmHnrs7ntjEu0qvM6DoZjIJ9Xa9CNOgFFtrgUxLB\nYvLwROa9+fVsffHVvm4c3Y4d0k0xdFkMFsepGLS75neNACPRzm23MnMolPo6SbPFomDfpFMhihOI\nQlyUk1SwGOTG0enYod0UZfCKQbtrvtWgXX6Plvg6SbvFogHaCOIMLqbx+L2QTB6eyLFHTI3991R5\nRkYIg/edDHJ+W6/5cXvAHkPRz227MvP+041H12/p+cjJospitppq9hHErZ2HtCinnRAGu0JcqZqF\nsvY/txtwHbTL7+HfPTf6zOhej5wsqixatgr2ESQRiJLuZ00qQKcdbPpJZ9o3xRBuas3Kukiqofma\nH7TL71Bg8VV3j8mjbo+cLKosunsV7LtoLhgh1c6TCtBpB5tB0pnW4GNIN7WGrMcp8m5FDHJu+33k\nZCeh3ehbZdGyLX2wH/QkdyoYeQ8CJRmg0ww2Ic2uCfGmBtkO3od0PvqRRB7lfZOLKu0KZakHaAdd\n7h3C0v5Og3ZJrkRNc1A0pBWzaaYlzrWS5eB9pzy4/+FNwQ4OQ/w8CqEs9yOpCQ/tlLZmH6cmk0Xz\nuluLo1tNJMnaYJpNx5CmnKaZlrjXSp6LpPYYYnRhXFlrvJrS+welrdnHqc2lPQ2wW4ujV00k6drg\nOXMP56rL5rPor9/KVZfN5+y5h8X860glnaGmJYlrJc3aXPMxWqdDwtDoCuii13g7tYQ7nZ8DhicG\n3aJJQ6Y1ezP7G2Bh/eWX3f3atI4VpzaXZo23V4sjSk0k6dpgWoOiIQ1qp5WWIk0Zbc6DF156lS9e\nc++Y3xe1xtutJdzu/LzlyGlc+qXlwffhJy3rbpxb3P0bZjYBWAmkFuzjFsK0gkOvYB71JhXC/upR\nhJTOLG9qoc7+aOTB5q3bgulmiyNKd23z+Tlg0kQu/aflwQxUZ3mdZBrs3f039R931P9LVdyA3Ss4\nDHKiegXzItUU5Q+ar5V+d3zMQ1mus6h98o3z8+sHnw6mDz/rWUJ5DdD+Z+D6LA4U2rztKIUspO4P\nCH+Ockii7Pj456fP5E2vn5x7foZ2nQ2i3+7aUCYO5DEVNpVgb2bTgOta3t7g7u83sxOAdwHvTePY\nWYh7oqIUstabVF4BtyhzlEMRZcfHa25ygCDyM6RutkH020JJo0VThAV1kFKwd/cNwLzW983sYOAf\ngXPcfWcax85CtznLk16zZ6ST3k8hyyvg5r0Qp4gtiig7PjaklZ9FzLc4+m2hJNmiKcKCutHvT+2b\n2/s7YCrwIzMDeKe7v5xxGmLLcs5yngE3zznKRW1RtNYcx+0BIyND7BppH/yTyM/m4N664VhR8i2u\nflsoSbRo4pTNPMZMsh6g/VCWx0tLpwLdOmc5iYCcZ8DNq38z7xZFFN1qz912fGwVNz+bb4qtN5YQ\n861MirKgrqG0K2jTltWc5TwHlPKasRH6qscorY5OOz62btcbJz9bb4o7dwGEm29lk0TZzHLMRME+\nhizmLOc9RS6PGRuhzJhoZ9BWR+NaOfaIqcw/4Q2J5Ge7m2KrUPKtTJpbdWk9uzoNCvYJSDsg5z1F\nLusZG1HzM49ClESrI6n87DR2NDQ0NDp2VMS58yFr16q76rL5iT+7Og0K9glJOyAXfYpcv9J6AHxc\nIbU6Ot0U337M9CBqkmXTrVXXbm/90MaeFOwTVLWAnLZO+ZlnIcq7W61Vp5tiEa7DULo3ouq3VRfa\n2JOCvRRO3oUo7261VkWsZITUvRFVUVfrNpR2i2Mpr7S3oI4ii22Jy6poDxRp6Her7JC2+QbV7Euj\naE3iOELrSpH+5N0yiyPP1bpxKdiXQBGbxHGFVIikP6F1b/Qrj9W6SVA3TsEVtUmcBHWlFFNo3RtJ\n6/TUrLypZl9wnZrEN9/9OGec8IbSFCApl7K2zEJuZatmX3DtBisBrrnpwd2ebysSkrK1zKK2svOq\n+SvYF1xrk7hZlbp0khZqU1zSkcT57jbw3LB0xSMsXHwLi765MvPKmLpxSqDRJL757se55qYHx/yu\nKLMcQhJyU1ySl9T57jXwnPeKWtXsS2Ly8ETOOOENuc8/L7oqD3hXUZLnu9fAc5Saf5pUsy8RzT+P\nr8hzwKV/SZ/vbgPPeU85VbAvmbLOcshK3gVSspXG+e40rz7vypiCfQmFsoijiPIukJKtrM93npUx\nBXuRFmodVUvW5zuvypiCvUgbZWkdVWnPpDjKcr67UbAPhAqlJE1TSKWZgn0AVCglaXnP6ZbwaJ59\nzjSvW3oZZHVn3nO6i6ysq6dVs8+Z5nVLN4O2+jSFdDBlbmWrZp+zEJ66JGGK0+or+zbCaSh7K1s1\n+5xpXrd0ErfVpymk/Sl7K1vBPgAqlNJOEl0xVZhSmJSyd32pGycQZdvbW+JTV0y2OuU3UIoB21xq\n9mZ2PfBv7n5ZHscXKQq1+rLVmt8r7n+ShYtvKcWAbeY1ezM7Gtg76+OKFJVafdlq5PcIlGrANo9u\nnA8DX8vhuCIikZVtrUKmwd7MjgA2Ac9leVwRkX6VbVp0Kn32ZjYNuK7l7Q3AVuDvgCPSOK6ISFLK\nNi06lWDv7huAea3vm9nNwNXAAcBrzexWd1+eRhpEROIq0wB5prNx3P0MADObB5ymQC8ioSvLWoVc\npl66++3A7XkcW0SkirSoSkSkAhTsRUQqQMFeRKQCFOxFRCpAwV5EpAIU7EVEKkDBXkSkAhTsRUQq\nQMFeRKQCFOxFRCpAwV5EpAIU7EVEKiCXjdAiGAewYcOGvNMhIlIITfFyXLvfhxrsDwJYsGBB3ukQ\nESmag4BHWt8MNdivAuYCTwE7c06LiEgRjKMW6Fe1++XQyMhIu/dFRKRENEArIlIBCvYiIhWgYC8i\nUgGhDtAOzMyuBI4D7nX3j+SclunAjcCRwL7uviOU9JnZCcCVwC5glbtfYmYfA94DPA5c4O7bc0jX\nbOAb1AbmHwYWAlcQQJ4BmNklwJ+6+9sDOpczgLuBfwdedff5IZzLpvT9JXA+tQHEBcBHyTnfzOxM\n4BONl8BFwJsIIM/M7DXAEmAfYAtwHvBZYuZZqWr2ZvZmakF1LrCnmR2fc5KeBU4FVkJw6XscOMXd\n3w5MMbOTgJPrrx8A3ptTutzdT6znEcBbCCTPzGwvYE7955DOJcCt7j6vHuinEMa5xMwOBk5y91Pd\nfR4wlQDyzd1vqufXPOC3wL0EkmfAmcDd9bTdQ+2mFDvPShXsgbcCt9Z//jnwthzTgrtvc/fNTW8F\nkz533+Du2+ovtwOzgNvrr3NLW0tt6hVqN8sg8gy4EPhO/edgzmXdyWa2ot7yOI4AzmXdGcA4M7vN\nzL5ST0sw+WZmhwFPA7MJJ88eoVarB9gfGCGBPCtbsN8f2Fr/eUv9dUiCS5+ZHQ28DniOQNJmZueY\n2WpqtcAJIaTLzCYA89x9Wf2tkM7lU8BM4GTgNGrBPpS0TQX2dPdTgZeA/QgnbQB/AvyYsM7nQ8Db\nzGwNtXO5gwTSVrZgvwUYrv+BC7sHAAADfElEQVQ8TC2AhSSo9JnZAcBXqdVYg0mbuy9199nAE9Qu\n9BDS9RfAtU2vQ8qvV9z9RXffQW2M6JFQ0kYtn5bXf14GDBFO2gDOBpYS0PmkNr5xg7vPAn5CrcIT\nO21lC/Z3UWv2Q62GszLHtLQTTPrMbDzwXeCj7r6B2qq7k/JOW71fvGErtSZsCHlmwEVmdhO1Lq8D\nA0kXZjap6eUfUxvYzv1c1t0JHF3/eQ7hnE/MbBq1Ae3fE8j1XzdEbbwP4Jn6v7HzrFTB3t3vBbaZ\n2Qpgp7vfk2d6zGyCmf0cOAa4mdodOpT0nQscD3zezG4HDgd+aWZ3UCuU/zendJ1pZsvNbDm1LoDP\nEkCeufvH3f0Mdz8TWOPunwohXXVzzezXZnYn8KS7300Y5xJ3vw94uX6NHQ98kXDy7T3A9QDuvpFA\n8oxaC/K8ep4tAL5CAnmm7RJERCqgVDV7ERFpT8FeRKQCFOxFRCpAwV5EpAIU7EVEKkDBXkrBzGaY\n2bs6/G6amX2my/87z8y+ml7qwMxuN7MD0zyGSDcK9lIWM4Ddgr2ZjavvA/TJ7JMUj5m1fXC0yCBK\nt8WxVNZHgLfUt0j+EPAj4GfAHDP7K+Cr7n6WmV0KnEVtf5H/4+7faPdl9W2DrwMepbZJ1qfd/Toz\nu7r+Xf/PzD5KbYXj7cD36589GvgMtV0TjwD+1t1/Wv/aRWZ2JLARWODuO83sC9QWG40HLnX3e+p7\novyM2uKe0xLLIak01eylLP4J+HF921qntmvgd929NVj+s7ufQm3XyovNrFsZmApcAJxCbQ/2bqZQ\n29PkA8AXqO2l8x5qN6GG2+rHfhp4j5m9ExhX38r2fcDn6p/rlHaRgSnYS1m9XF+q3+rPzeyX1LaM\nnQ5M7vIda939VXd/BphYf695yflQ62eBJ4F17v5y/ecDmj5zT9O/b6LWYjijvix+CX/Y7KpT2kUG\npmAvZfEqY7sld3b43P+g1jUyn9rugUMdPgdjA3vDZuCQ+s9/1OGznW4IxzX9+zCwFlja9BCNxj7l\nndIuMjD12UtZrAZmmdkPgb/t8rmbgV/VP7+5y+c6+RZwrZldALzc5/87z8wuBjYB19cfUzm3XrMf\noVbj//gAaRLpSRuhiYhUgLpxREQqQMFeRKQCFOxFRCpAwV5EpAIU7EVEKkDBXkSkAhTsRUQqQMFe\nRKQC/j9O65JevEf9KwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.95\n" + ] + } + ], + "source": [ + "plt.plot(mappings[y==0], \"o\")\n", + "plt.plot(mappings[y==1], \"*\")\n", + "plt.legend([\"0\", \"1\"]);\n", + "plt.ylim(-5, 7); plt.xlabel(\"trial number\"); plt.ylabel(\"LDA component\")\n", + "plt.show()\n", + "print(np.sum(est_y == y) * 1.0 / len(y))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## use cross validation to measure the classification accuracy" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def classification_acc(X, y, signal_dur=1.0, signal_delay=1.0, num_splits=10):\n", + " # get the index in which data are used for training and testing: 1 seconds after cue onset\n", + " signal_st = cue_time + signal_delay\n", + " signal_et = signal_st + signal_dur\n", + " st_idx, et_idx = int(signal_st * sr / down_ratio), int(signal_et * sr / down_ratio)\n", + " # accuracy list\n", + " acc_list = []\n", + " rs = ShuffleSplit(n_splits=num_splits, test_size=0.2, random_state=7) # split the data\n", + "# pdb.set_trace()\n", + " for train_idx, test_idx in rs.split(X):\n", + " # get training and testing data\n", + " tr_X, tr_y = X[train_idx], y[train_idx]\n", + " test_X, test_y = X[test_idx], y[test_idx]\n", + " # train CSP based on the train data, transform the test data using the CSP from training data\n", + " csp = CSP(n_components=6, reg=0.001, log=None, norm_trace=False, transform_into=\"average_power\")\n", + " train_tf_X = csp.fit_transform(tr_X[:, :, st_idx:et_idx], tr_y)\n", + " test_tf_X = csp.transform(test_X[:, :, st_idx:et_idx])\n", + " # train the LDA on the train data\n", + " clf = LinearDiscriminantAnalysis()\n", + " clf.fit(train_tf_X, tr_y)\n", + " est_y = clf.predict(test_tf_X)\n", + " acc = np.sum(est_y == test_y) * 1.0 / len(test_y)\n", + " acc_list.append(acc)\n", + " return(np.mean(acc_list), acc_list)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Average classification accuracy on testing set is: 0.906250\n" + ] + } + ], + "source": [ + "avg_acc, acc_list = classification_acc(filtered_X, y, signal_dur=1.5, num_splits=10)\n", + "print(\"Average classification accuracy on testing set is: %f\" % (avg_acc))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.14" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +}