...
 
Commits (13)
# ---------------------------------
# Finds module CSV
# Adds library to target
# Adds include path
# ---------------------------------
GET_PROPERTY(OV_PRINTED GLOBAL PROPERTY OV_TRIED_OpenViBEModuleCSV)
SET(INCLUDED_OV_SDK_COMPONENTS CSV)
INCLUDE(AddOpenViBESDKComponents)
SET_PROPERTY(GLOBAL PROPERTY OV_TRIED_OpenViBEModuleCSV "Yes")
# Add all the subdirs as projects of the named branch
OV_ADD_PROJECTS("CONTRIB_APPLICATIONS_DEMOS")
PROJECT(p300-speller-unsupervised)
INSTALL(DIRECTORY share/ DESTINATION ${DIST_DATADIR}/openvibe/scenarios/bci-examples/p300-speller-unsupervised/ PATTERN ".svn" EXCLUDE REGEX "signals-test/.*csv" EXCLUDE)
Unsupervised P300 speller
=========================
Status: Very experimental at the moment. Not well tested.
Relying almost fully on the work of David Hübner,
Thibault Verhoeven, Pieter-Jan Kindermans, Michael
Tangermann and others(*). Many thanks to the original authors.
For original code, see https://github.com/DavidHuebner/Unsupervised-BCI-Matlab
*) These original authors are not responsible for any issues introduced by the OpenViBE integration
What is it
----------
These scenarios provide a P300 speller without explicit training phase. Simply focus on
the letters of your choice (but not the hash marks) and the used machine learning technique
should catch on after a while. Typically you should notice the predictions start
to improve after 5 letters or so, but it may also take longer.
Most people will be interested to run the scenario 'unsupervised-p300-online-1.xml'. The
other scenarios are intended solely for testing purposes and may require supplementary files.
Requirements
------------
Matlab at the moment.
Note that P300 can be very time sensitive. Its recommended to minimize background processes on the
computer. On Windows, always run OpenViBE release build outside Visual Studio, and set the Power
Options to 'High Performance'.
Design notes
------------
* The flash groups and sequences (type 1 or type 2) are read from CSV files. To allow the user
to choose the flash timings, these files are used like queues: flash 1 gets the group
and sequence of line 1 in those files, flash 2 gets associated with line 2, and so on. This
is currently done inside the "P300 Speller Stimulator II" box.
* This speller is 'multiletter' meaning that it can correct previous letters. This is why
you see the previously spelled text change during the spelling. Each prediction from the
classifier is expected to be a stimulation set with markers coding
'segmentStart','row','col','row','col',...'segmentStop'
which is then mapped to multiple letters. The last letter is highlighted in red in the grid.
See P300_Process.m for the actual event id's used.
* The timeline and target letter sequence is generated by "P300 Speller Stimulator II"
box, except in those test scripts that rely purely on canned CSV data.
* The speller is currently limited to a buffer of 15 last letters to keep the computation
delays acceptable.
* Many parameters are presently hardcoded. The code could use lots of polish, and the
matlab parts should be ported to C++.
* Software tagging is currently used. The P300 Speller Visualisation II forwards the
flash event markers to Acquisition Server (using TCP Tagging) after each flash render.
Glossary
--------
What David & co. call 'stimuli' in their code are called 'Flash Groups' in the OV scenarios;
The set of flash time event markers are called 'Timeline'. 'Sequence' is called 'Group types'.
How to test & debug
-------------------
* See the scenarios named like unsupervised-p300-test*xml. Required data files can be imported
from the format used by David Huber & co. ( https://zenodo.org/record/192684 ) by the script
'matlab-test/dump_data_to_csv.m'. Using the S1.mat file from the aforementioned url should get
the 'epoched' test scenarios to print 'franzy_jagt_im'... etc.
- J.T. Lindgren 01.Dec.2017
%
% CSV export datasets downloaded from
%
% https://zenodo.org/record/192684
%
% for initial testing of the P300 in OV
%
clear;
% Change to the current folder and add required paths to the Matlab repository
addpath('../matlab/helpers')
haveNonEpoched = false;
haveEpoched = false;
if(exist('S1_cnt_mrk.mat','file'))
haveNonEpoched = true;
end
if(exist('S1.mat','file'))
haveEpoched = true;
end
assert(haveEpoched || haveNonEpoched);
if(haveNonEpoched)
rawdata = load('S1_cnt_mrk.mat');
sampleTimes = (0:(size(rawdata.cnt.x,1)-1))'/rawdata.cnt.fs;
featTime = rawdata.mrk.time'/1000; % orig time is in milliseconds
stims = [0,32769,0];t=0;
for i=1:size(rawdata.mrk.time,2);
stims=[stims;featTime(i),32779,0; featTime(i)+0.1,32780,0];
end
stims = [stims;featTime(end)+2,32770,0];
exportcsv20('../signals-test/01-raw-data.csv',rawdata.cnt.x, sampleTimes, false, rawdata.cnt.fs);
exportcsv('../signals-test/02-flashes.csv',stims); % 20's csv reader has issues atm for stims
exportcsv20('../signals-test/03-groups.csv',rawdata.mrk.stimuli', featTime, true);
exportcsv20('../signals-test/04-sequence.csv',rawdata.mrk.sequence', featTime, true);
if(haveEpoched)
load('S1.mat');
fv = proc_jumpingMeans(epo,[50 120; 121 200; 201 280;281 380;381 530; 531 700]);
data = reshape(fv.x,31*6,12852)';
data = data(1:size(featTime,1),:);
exportcsv20('../signals-test/01-epoched-data.csv',data,featTime, false);
end
end
if(haveEpoched && ~haveNonEpoched)
load(fullfile(your_data_path,'S1.mat'));
load(fullfile(your_data_path,'sequence.mat'));
fv = proc_jumpingMeans(epo,[50 120; 121 200; 201 280;281 380;381 530; 531 700]);
data = reshape(fv.x,31*6,12852)';
stimuli = epo.stimuli;
% We lack the original flash onsets, make up
flashDuration = 0.1;
featTime = (0:(size(data,1)-1))'*flashDuration;
stims = [0,32769,0];t=0;
for i=1:size(rawdata.mrk.time,2);
stims=[stims;featTime(i),32779,0; featTime(i)+0.1,32780,0];
end
stims = [stims;featTime(end)+2,32770,0];
exportcsv20('../signals-test/01-epoched-data.csv',data,featTime,false);
exportcsv('../signals-test/02-flashes.csv',stims);
exportcsv20('../signals-test/03-stimuli.csv',stimuli',featTime,true);
exportcsv20('../signals-test/04-sequence.csv',sequence,featTime,true);
end
% verify
if(0)
% it appears in the epoched files, the first 4284 rows of sequence+stimuli correspond to the 4284 rows of
% the rawdata file
rawdata = load('S1_cnt_mrk.mat');
S1 = load('S1.mat'); % epo
seq = load('sequence.mat');
assert(max(max(abs(rawdata.mrk.stimuli - S1.epo.stimuli(:,1:4284))))==0);
assert(max(max(abs(rawdata.mrk.sequence' - seq.sequence(1:4284))))==0);
end
function exportcsv( fn, mat, freq )
% Can be read as feature vectors or signal in openvibe. For signal, provide sampling frequency 'freq'.
h=fopen(fn,'w');
fprintf(h,'Time(s)');
for i=1:size(mat,2)-1
fprintf(h,';f%03d',i);
end
if(nargin==3)
fprintf(h,';Freq\n');
for i=1:size(mat,2)
fprintf(h,'%.4f;',mat(1,i));
end
fprintf(h,'%d\n',freq);
startIdx = 2;
else
fprintf(h,'\n');
startIdx = 1;
end
fclose(h);
dlmwrite(fn,mat(startIdx:end,:),'-append','delimiter',';','precision','%.4f');
end
function exportcsv20( fn, mat, sampleTime, isInt, freq, stims)
% Can be read as feature vectors in openvibe 2.0
epochSize = 32;
if(nargin<4)
isInt = false;
end
if(nargin<5)
freq = 0;
end
if(nargin<6)
stims = [];
end
h=fopen(fn,'w');
if(freq>0)
fprintf(h,'Time:%dHz,Epoch',freq);
numEpochs = floor(size(mat,1)/epochSize);
endTime = repmat((0:(numEpochs-1)),[epochSize 1]);
endTime = endTime(:);
diffSize = size(mat,1)-size(endTime,1);
if(diffSize>0)
endTime = [endTime;repmat(numEpochs,[diffSize 1])];
end
timePattern = '%.4f,%d';
else
fprintf(h,'Time:%d,End Time',size(mat,2));
% Assumes equispaced vectors in time
duration=sampleTime(2)-sampleTime(1);
endTime = sampleTime + duration;
timePattern = '%.4f,%.4f';
end
for i=1:size(mat,2)
fprintf(h,',f%03d',i);
end
fprintf(h,',Event Id,Event Date,Event Duration\n');
if(isInt)
pattern = ',%d';
else
pattern = ',%.4f';
end
% can't use dlmwrite as we need 'empty' events... ouch slowness
stimCnt=1;
for i=1:size(mat,1)
fprintf(h,timePattern,sampleTime(i),endTime(i));
for j=1:size(mat,2)
fprintf(h,pattern,mat(i,j));
end
if(isempty(stims) || stimCnt>size(stims,1))
fprintf(h,',,,\n');
else
if(stims(stimCnt,1)<=sampleTime(i))
fprintf(h,',%d,%.4f,0\n',stims(stimCnt,2),stims(stimCnt,1));
stimCnt = stimCnt + 1;
else
fprintf(h,',,,\n');
end
end
end
fclose(h);
end
function likelihoods = classifier_compute_log_likelihoods(classifier,projection)
% Compute the LOG!! likelihood based on the projection
% Originally written by Pieter-Jan Kindermans
%
% Commented by
% david.huebner@blbt.uni-freiburg.de, 11-2015
%
% In this function the probability p(x|w,a) = N(x^T*w|y(a),beta) is
% computed (see Page 6 in the paper "True Zero-Training Brain-Computer
% Interfaces"). Here, beta is chosen to be 1.
%
% classifier.label are the rescaled target and non-target class labels.
% One can show that when they are chosen as N/N+ and N/N- respectively,
% that this leads to an equivalence between least square regression and
% linear discriminant analysis in terms of the projection w. For instance,
% in the paradigm of the MIX scenario, there are N+ = 8 targets and N-=26
% non-targets contained in every 34 symbols.
%
% For more information, see:
% Bishop's Machine learning book, Springer, 2006, Chapter 4.1.5
%
% The function gauss_log_lik then computes the log likelihood of the projections
% given the two different class means and assuming a unit variance.
%
% The return value is a cell with 1xlength(projection) entries where each
% entry contains one vector for the target probabilities and one for
% non-target
likelihoods={};
for c_i =1:length(projection)
likelihoods{end+1}=[gauss_log_lik(projection{c_i},classifier.label(1)),...
gauss_log_lik(projection{c_i},classifier.label(2))];
end
end
\ No newline at end of file
function projection = classifier_compute_projection(classifier,data)
%% Computes w^T*x for each trial
projection = {};
for d_i = 1:length(data)
projection{end+1}=data{d_i}*classifier.w;
end
end
\ No newline at end of file
function [probs,data_log_likelihood] = classifier_compute_trial_probs(classifier,likelihoods,stimuli)
%% Compute the probabilities for the symbols given the log likelihoods of the projections, the stimuli and the classifier
% likelihoods -- cell array, with per cell (Sx2) likelihoods,
% first column, contains likelihood for positive stimulus
% second column, contains likelihood for negative
% stimulus
%
% stimuli -- cell_array with per cell (nr_commands X S) matrix,
% BINARY ARRAY to indicate whether command a was intensified
% during stimulus s
%
% By Pieter-Jan Kindermans, 2014
%
% Modified and commented by david.huebner@blbt.uni-freiburg.de, 2017
data_log_likelihood = 0;
probs = zeros(length(likelihoods),classifier.nr_commands);
% Loop over each trial (spelled character)
for c_i =1:length(likelihoods)
cur_stimuli = stimuli{c_i};
cur_likelihoods = likelihoods{c_i};
% Sum the log-likelihood for each symbol and add the prior
for a_i = 1:classifier.nr_commands
index = cur_stimuli(a_i,:);
probs(c_i,a_i) = probs(c_i,a_i)+sum(cur_likelihoods(index,1)); % Target
probs(c_i,a_i) = probs(c_i,a_i)+sum(cur_likelihoods(~index,2)); % Non-Targets
probs(c_i,a_i) = probs(c_i,a_i)-log(classifier.nr_commands); % Prior: 1/#Classes
end
% Normalize across selection. Effectively, this is computing
% exp(probs(d_i,:))/sum(exp(probs(d_i,:))
normalizing_constant = logsumexp(probs(c_i,:));
probs(c_i,:) = exp(probs(c_i,:)-normalizing_constant);
% ==> probs are now rescaled to be in [0,1] and to sum up to 1
% Update log-likelihood with logsumexp(probs(c_i,:)).
data_log_likelihood = data_log_likelihood+normalizing_constant;
end
end
function classifier=classifier_mix_add_letter(classifier,data,stimuli)
%% Add letter to the classifier
% parameters:
% classifier -- the classifier struct to which the data has to be added
% data -- cell_array with per cell (n_ept x feat_dim) matrix
% containing the features
% stimuli -- cell_array with per cell (nr_commands X n_ept) matrix,
% BINARY ARRAY to indicate whether command a was intensified
% during stimulus s
%
% By Pieter-Jan Kindermans, 2014
%
% Modified and commented by david.huebner@blbt.uni-freiburg.de, 2017
% Loop over each character
for c_i = 1:length(data)
cur_data = data{c_i};
cur_stimuli = stimuli{c_i};
% Store data and stimuli
classifier.data{end+1} = cur_data;
classifier.stimuli{end+1}=cur_stimuli;
% Cache XTY and XTX
classifier.XTX_list{end+1} = cur_data'*cur_data;
classifier.XTX = classifier.XTX + cur_data'*cur_data;
classifier.X2TX2 = classifier.X2TX2 + (cur_data.^2)'*(cur_data.^2);
temp_XTYp = zeros(classifier.feat_dim,classifier.nr_commands);
temp_XTYn = zeros(classifier.feat_dim,classifier.nr_commands);
% Loop over all possible selections
for a_i =1:classifier.nr_commands
index = cur_stimuli(a_i,:);
nr_pos = sum(index == 1);
nr_neg = sum(index == 0);
temp_XTYp(:,a_i) = temp_XTYp(:,a_i) + cur_data(index==1,:)'*(classifier.label(1)*ones(nr_pos,1));
temp_XTYn(:,a_i) = temp_XTYn(:,a_i) + cur_data(index==0,:)'*(classifier.label(2)*ones(nr_neg,1));
end
classifier.XTYp{end+1} = temp_XTYp;
classifier.XTYn{end+1} = temp_XTYn;
end
end
\ No newline at end of file
function [classifier] = classifier_mix_expectation(classifier)
%% Execute the expectation step for the EM algorithm
% During this step, we compute the probability for each option symbol per
% trial p(c_t|X_t,w)
%
% parameters:
% classifier -- the classifier struct should contain the data and the stimuli
%
% returns
% classifier --
% .probs: matrix [nr_trials x nr_commands] where entry(i,j)
% describes the probability that symbol j was target in trial i.
% .data_log_likelihood: A measure of how good the classifier
% is separating the data
%
% Internal variables are:
% projection: w^T*X for each epoch x
% likelihood: The target and non-target likelihood of each epoch x
% (It is computed based on the projection)
%
% --
% written by Pieter-Jan Kindermans
%
% Updated by david.huebner@blbt.uni-freiburg.de
projection = classifier_compute_projection(classifier,classifier.data);
likelihoods = classifier_compute_log_likelihoods(classifier,projection);
[probs,data_log_likelihood] = classifier_compute_trial_probs(classifier,likelihoods,classifier.stimuli);
classifier.probs = probs;
classifier.data_log_likelihood = data_log_likelihood;
end
\ No newline at end of file
function classifier = classifier_mix_maximization(classifier)
%% Execute the maximization step for the EM algorithm
%
% Parameters:
% classifier -- the classifier struct should contain the classifier.probs matrix.
%
% returns:
% classifier -- classifier.em_pos and classifier.em_neg are updated.
%
% written by Pieter-Jan Kindermans
%
% Rewritten by David Hübner
% Compute total number of epochs (N) and number of trials (nr_trials)
N=0;
nr_trials = length(classifier.data);
for i=1:nr_trials
N = N+size(classifier.data{i},1);
end
xTyp = zeros('like',classifier.w);
xTyn = zeros('like',classifier.w);
% Add contributions for each selection(trial) that has to be made
for c_i = 1:nr_trials
cur_probs = classifier.probs(c_i,:);
cur_XTYp = classifier.XTYp{c_i};
cur_XTYn = classifier.XTYn{c_i};
% Weigh the contribution per option with the probability of that option.
for a_i = 1:classifier.nr_commands
xTyp = xTyp + cur_probs(a_i)*cur_XTYp(:,a_i);
xTyn = xTyn + cur_probs(a_i)*cur_XTYn(:,a_i);
end
end
classifier.em_pos = xTyp/N;
classifier.em_neg = -xTyn/N;
end
# Unsupervised-BCI-Matlab
Unsupervised Learning for Event-Related Potential based Brain-Computer Interfaces.
## Overview
The provided code is a standalone Matlab implementation of the MIX method proposed by Thibault Verhoeven et al. in
> Verhoeven, T., Hübner, D., Tangermann, M., Müller, K.-R., Dambre, J., & Kindermans, P.-J. (2017). Improving zero-training brain-computer interfaces by mixing model estimators. Journal of Neural Engineering, 14(3), 036021.
This method allows a user to control an event-related potential (ERP) speller without a calibration phase. It continously learns from unlabeled data, i.e. when the user's intentions are unknown. The decoding performance is almost as good as a comparable supervised method with the same amount of data and full label access. Hence, this plug & play system is a powerful alternative to traditional supervised methods. For a visual ERP speller, about 3 minutes of unlabelled data are sufficient to perfectly decode the attended symbols.
## Get Started
Two simple steps are necessary to run this code. First, you have to download the example data set (__online_study_1-7.zip__ & __sequence.mat__) from Zenodo https://zenodo.org/record/192684 and extract it. Second, you have to clone this Github repository. Then you can open `example_script.m`, replace `your_data_path` with the data directory where you extracted the data and run it. The functions were tested with Matlab 2014b. The Matlab output should read something like:
-------Start simulating online experiment!--------
The classifier starts from a random initalization and is retrained after each trial.
The binary target-vs-non target area under the curve (AUC) on all data
until the current trial is then reported, if label information are available.
Additionally, the mixing coefficient for the target and non-target classes are reported.
-----------------------------------------------------
Trial 1. AUC: 64.663%. Gamma_pos: 0.370, Gamma_neg: 0.713. Runtime: 0.081s
Trial 2. AUC: 70.162%. Gamma_pos: 0.271, Gamma_neg: 0.476. Runtime: 0.085s
Trial 3. AUC: 78.980%. Gamma_pos: 0.245, Gamma_neg: 0.379. Runtime: 0.230s
Trial 4. AUC: 76.728%. Gamma_pos: 0.230, Gamma_neg: 0.329. Runtime: 0.137s
Trial 5. AUC: 82.918%. Gamma_pos: 0.226, Gamma_neg: 0.311. Runtime: 0.198s
Trial 6. AUC: 87.837%. Gamma_pos: 0.217, Gamma_neg: 0.272. Runtime: 0.267s
Trial 7. AUC: 86.641%. Gamma_pos: 0.217, Gamma_neg: 0.249. Runtime: 0.295s
Trial 8. AUC: 90.805%. Gamma_pos: 0.200, Gamma_neg: 0.227. Runtime: 0.293s
Trial 9. AUC: 92.076%. Gamma_pos: 0.196, Gamma_neg: 0.210. Runtime: 0.326s
Trial 10. AUC: 93.401%. Gamma_pos: 0.211, Gamma_neg: 0.218. Runtime: 0.362s
Trial 11. AUC: 97.160%. Gamma_pos: 0.215, Gamma_neg: 0.215. Runtime: 0.381s
Trial 12. AUC: 97.817%. Gamma_pos: 0.224, Gamma_neg: 0.219. Runtime: 0.666s
Trial 13. AUC: 97.800%. Gamma_pos: 0.236, Gamma_neg: 0.228. Runtime: 0.447s
Trial 14. AUC: 97.988%. Gamma_pos: 0.237, Gamma_neg: 0.224. Runtime: 0.495s
Trial 15. AUC: 97.893%. Gamma_pos: 0.244, Gamma_neg: 0.227. Runtime: 0.543s
Trial 16. AUC: 97.764%. Gamma_pos: 0.260, Gamma_neg: 0.244. Runtime: 0.572s
Trial 17. AUC: 98.000%. Gamma_pos: 0.258, Gamma_neg: 0.241. Runtime: 0.636s
Trial 18. AUC: 98.327%. Gamma_pos: 0.265, Gamma_neg: 0.247. Runtime: 0.696s
Trial 19. AUC: 98.387%. Gamma_pos: 0.263, Gamma_neg: 0.244. Runtime: 0.634s
Trial 20. AUC: 98.411%. Gamma_pos: 0.268, Gamma_neg: 0.247. Runtime: 0.949s
Trial 21. AUC: 98.346%. Gamma_pos: 0.264, Gamma_neg: 0.240. Runtime: 0.732s
Trial 22. AUC: 98.241%. Gamma_pos: 0.267, Gamma_neg: 0.241. Runtime: 0.717s
Trial 23. AUC: 98.181%. Gamma_pos: 0.270, Gamma_neg: 0.247. Runtime: 0.796s
-----------------------------------------------------
Simulation completed.
Find the final classifier in C.classifier.
Statistics are saved in C.statistics.
-----------------------------------------------------
Final sentence: franzy jagt im komplett
Of course, you can also use the algorithm with your own data. Just make sure that the data is in the right format and satisfies the prerequisites of the algorithms (see below).
## Background
The proposed MIX method is an analytical combination of two other unsupervised learning method. Depending on the variance of the individual estimators, it assigns a higher or lower weight to the invididual classifier.
The first one is an __Expectation-Maximization (EM)__ algorithm for a Gaussian mixture model. It was initially proposed for brain-computer interfaces by Pieter-Jan Kindermans et al. in 2012:
> Kindermans, P. J., Verstraeten, D., & Schrauwen, B. (2012). A Bayesian model for exploiting application constraints to enable unsupervised training of a P300-based BCI. PloS one, 7(4), e33758.
The EM heavily utilizes that ERP paradigms introduce a certain structure in the data: Knowing the attended symbol uniquely determines for each highlighting event (e.g. row or column intensification in the classical matrix speller by Farwell & Donchin), whether it was a target or non-target. Hence, we can define the attended symbol as the latent variable. The algorithm then proceeds by alternatively estimating the value of the latent variable [E-step] and maximizes the parameter given these estimates [M-step]. It is known to work well in practice, even though no theoretical guarantees exist regarding convergence. The initialisation is random.
The second unsupervised learning method is __Learning from Label Proportions (LLP)__. It was first applied to BCI by Hübner et al. in 2017:
> Hübner, D., Verhoeven, T., Schmid, K., Müller, K.-R., Tangermann, M., & Kindermans, P.-J. (2017). Learning from label proportions in brain-computer interfaces: Online unsupervised learning with guarantees. PloS one, 12(4), e0175856.
It relies on a simple general principle by Quadrianto et al. (2009), that allows inferring information about class means when the unlabelled data can be divided into at least two sub-groups that have different, but known class proportions. Then, the class means can be found by solving a simple linear system of equations.
> Quadrianto, N., Smola, A. J., Caetano, T. S., & Le, Q. V. (2009). Estimating labels from label proportions. Journal of Machine Learning Research, 2349-2374.
To enable this approach in BCI, a modification of the paradigm is necessary. Additional '#'-symbols are introduced to the spelling grid which take the role of visual blanks and as such, are never attended by the user. Now, two sequences S1 and S2 are introduced where S1 highlights only highlights ordinary symbols and S2 also highlights the visual blank '#' symbols. This leads to a higher probability of highlighting the attended character in S1 than in S2. These different probabilities can be quantified by (known) target and non-target ratios. S1 and S2 are hence, a linear combination of target and non-target events. The target and non-target class means can then be reconstructed by averaging the responses for S1 and S2 and solving the linear system, see the schematic overview:
![Schematic overview of LLP](readme_LLP.png)
Other ERP paradigms naturally have these different target to non-target ratios, e.g., in paradigms with a two-step selection procedure, where the number of items differs in the first and second step. LLP is determinstic. It was shown to outperform EM when only unlabeled data from a few trials is avaialable, but falls behind when more data is available.
___
@David Hübner, 10.8.2017.
Contact: david.huebner@blbt.uni-freiburg.de
%% Example script
% This function simulates the training of an unsupervised classifier on
% visual ERP data
%
% To begin the example, load data from the Zenodo repository
%
% https://zenodo.org/record/192684
%
% and extract it to a local folder <your_data_path>;. Also, make sure that
% you cloned this whole GitHub repository.
%
% https://github.com/DavidHuebner/Unsupervised-BCI-Matlab
%
% You can then run this example.
%
% Please report any problems to david.huebner@blbt.uni-freiburg.de
%
% @David Hübner, Jul, 2017
% !! Set your data path !!
your_data_path = 'data';
% Change to the current folder and add required paths to the Matlab repository
addpath('EM')
addpath('helpers')
% Load stimuli and epoch data for subject 1 (out of 13)
load(fullfile(your_data_path,'S1.mat'));
% Load the sequence data for LLP indicating whether epoch k was part of
% sequence 1 or sequence 2. It is the same for all subjects.
load(fullfile(your_data_path,'sequence.mat'));
% Extract a feature vector where the features are the mean amplitudes in
% the 6 given intervals
fv = proc_jumpingMeans(epo,[50 120; 121 200; 201 280;281 380;381 530; 531 700]);
% Bring feature matrix in the shape [N * feat_dim]
% where feat_dim = n_channels * n_time_intervals
data = reshape(fv.x,31*6,12852)';
% Next, look at "stimuli" which encodes which symbols are highlighted for
% each visual highlighting event (stimulus)
% Shape: [feat_dim * N] with entry(i,j)=1 if symbol i was highlighted
% during epoch j and 0 else.
stimuli = epo.stimuli;
% Denote the number of epochs per trial, i.e. the number of highlighting
% events which occur when writing a single character
nr_ept = 68;
% Decide about the classifier:
% gamma = -1; MIX, gamma = 0; EM, gamma = +1; LLP
gamma = -1;
% Mixing matrix for the LLP method describing the target and non-target
% ratio for the different subgroups in the data.
% !! Note that only two subgroups are supported at the moment !!
A = [3/8 5/8; 2/18 16/18];
% [OPTIONAL] Matrix y contains the label information for performance assessment
% Shape: [N] where entry at position k is 1 if epoch k was a target and 0
% if it was a non-target
y = epo.y(1,:);
% [OPTIONAL]: Set seed. The classifier is randomly initialized and a
% fixed seed ensures reproducibility of the results
rng(1234);
% [OPTIONAL]: Which trials should be analyzed If no values "[]" are given,
% then the whole data set is analyzed.
trials = 1:23;
%%%%%%%%%%%%%%%%%%%%%%%%%%% Run %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C = simulate_MIX([],data,stimuli,sequence,y,nr_ept,A,gamma,trials);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Obtain and display spelled symbols
[probs,select]=max(C.classifier.probs,[],2);
disp(['Final sentence: ' convert_position_to_letter(select)])
function classifier = clear_fields(classifier,feat_dim)
classifier.data = {};
classifier.stimuli = {};
classifier.XTX_list = {};
classifier.XTYp = {};
classifier.XTYn = {};
classifier.XTX = zeros(feat_dim,feat_dim);
classifier.X2TX2 = zeros(feat_dim,feat_dim);
end
\ No newline at end of file
function [EM_var_pos, EM_var_neg] = compute_EM_variance(C,data_norm,mean_pos,mean_neg)
% Computes the variance of the EM-estimator.
% Returns one variance value for each of the two classes
%
% !!! Important: This function assumes that the data is whitened !!!
%
% Originally by Thibault Verhoeven, 2017
[N,feat_dim] = size(data_norm);
prec = eye(feat_dim);
dist_p = data_norm*C.classifier.Whiten - repmat(mean_pos'*C.classifier.Whiten,N,1);
dist_n = data_norm*C.classifier.Whiten - repmat(mean_neg'*C.classifier.Whiten,N,1);
p1 = -0.5*sum(times(dist_p,dist_p),2)*(8.0/34);
p2 = -0.5*sum(times(dist_n,dist_n),2)*(26.0/34);
H11 = zeros(feat_dim,feat_dim);
H22 = zeros(feat_dim,feat_dim);
for i=1:N
distmat = dist_p(i,:)'*dist_p(i,:);
temp = (p1(i)+p2(i))*p1(i)*(distmat-prec);
temp = temp - p1(i)*p1(i)*distmat;
H11 = H11 + temp/((p1(i)+p2(i))*(p1(i)+p2(i)));
distmat = dist_n(i,:)'*dist_n(i,:);
temp = (p1(i)+p2(i))*p2(i)*(distmat-prec);
temp = temp - p2(i)*p2(i)*distmat;
H22 = H22 + temp/((p1(i)+p2(i))*(p1(i)+p2(i)));
end
EM_var_pos = trace(inv(-1.0*H11));
EM_var_neg = trace(inv(-1.0*H22));
end
\ No newline at end of file
function letter = convert_position_to_letter( pos )
% Converts a position in a matrix to the corresponding letter using the
% following matrix configuration:
%
% A B # C # D E
% # F G H I J #
% K L M # N O P
% Q # R S T U #
% V W X Y # Z ␣
% # . # , ! ? <-
A = ['a','b','#','c','#','d','e','#','f','g','h','i','j','#','k','l','m',...
'#','n','o','p','q','#','r','s','t','u','#','v','w','x','y','#','z',...
' ','#','.','#',',','!','?','<'];
letter = A(pos);
end
function log_lik = gauss_log_lik(X,mu)
% Computes the log-likelihood for data points X of a normal distribution
% with mean mu and unit variance(!)
log_lik = -log(sqrt(2*pi))-0.5*((X-mu).^2);
end
\ No newline at end of file
function C = init_classifier(w,feat_dim,nr_commands, nr_ept, nr_tpt)
% The C struct has two fields: "classifier" which is used for storing the
% current data, performing the EM steps, caluclating the projection and
% predicting the attended character and the "statistics" field which is
% storing useful information from all trials
C.classifier = struct();
C.classifier.w = w;
C.classifier.nr_commands = nr_commands;
C.classifier.feat_dim = feat_dim;
% classifier.label are the rescaled target and non-target class labels.
% One can show that when they are chosen as N/N+ and N/N- respectively,
% that this leads to an equivalence between least square regression and
% linear discriminant analysis in terms of the projection w. For instance,
% in the paradigm of the MIX scenario, there are N+ = 16 targets and N-=52
% non-targets contained in every 68 symbols.
%
% For more information, see:
% Bishop's Machine learning book, Springer, 2006, Chapter 4.1.5
C.classifier.label = [nr_ept/nr_tpt -nr_ept/(nr_ept-nr_tpt)];
C.statistics.w_init = w;
C.statistics.gamma_pos = {};
C.statistics.gamma_neg = {};
C.statistics.auc = [];
C.statistics.data_log_likelihood = [];
C.statistics.w = {};
C.statistics.projection = {};
C.statistics.probs = {};
end
\ No newline at end of file
function out = logsumexp(x)
% By Pieter-Jan Kindermans
%
% Comment by david.huebner@blbt.uni-freiburg.de:
% There is an equivalence trick used here in order to avoid overflow of the values,
% using the identity: log(sum(exp(x))) = x_max+log(sum(exp(x-x_max))
% see https://hips.seas.harvard.edu/blog/2013/01/09/computing-log-sum-exp/
% for a more detailed explanation
x_max = max(x);
x_temp = x-x_max;
x_exp = exp(x_temp);
x_expsum = sum(x_exp);
out = log(x_expsum)+x_max;
end
\ No newline at end of file
function loss= loss_rocArea(label, out, varargin)
%LOSS_ROCAREA - Loss function: Area over the ROC curve
%
% Synopsis:
% LOSS= loss_rocArea(LABEL, OUT)
%
% IN LABEL - matrix of true class labels, size [nClasses nSamples]
% OUT - matrix (or vector for 2-class problems) of classifier outputs
%
% OUT LOSS - loss value (area over roc curve)
%
% Note: This loss function is for 2-class problems only.
%
% Benjamin Blankertz
%
% Taken from the BBCI Toolbox distributed with a MIT License
% https://github.com/bbci/bbci_public
%
% Copyright (c) 2015 The BBCI Group
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
% SOFTWARE.
if size(label,1)~=2,
error('roc works only for 2-class problems');
end
N= sum(label, 2);
lind= [1:size(label,1)]*label;
%resort the samples such that class 2 comes first.
%this makes ties count against the classifier, otherwise
%loss_rocArea(y, ones(size(y))) could result in a loss<1.
[so,si]= sort(-lind);
lind= lind(:,si);
out= out(:,si);
[so,si]= sort(out);
lo= lind(si);
idx2= find(lo==2);
ncl1= cumsum(lo==1);
roc= ncl1(idx2)/N(1);
% area over the roc curve
loss= 1 - sum(roc)/N(2);
end
\ No newline at end of file
function dat= proc_jumpingMeans(dat, nSamples, nMeans)
%PROC_JUMPINGMEANS - compute the mean in specified time ivals of epochs
%
%Synopsis:
% dat= proc_jumpingMeans(dat, nSamples, <nMeans>)
% dat= proc_jumpingMeans(dat, intervals)
%
%Arguments:
% dat - data structure of continuous or epoched data
% nSamples - number of samples from which the mean is calculated
% if nSamples is a matrix means over the rows are given
% back. nMeans is ignored then.
% nMeans - number of intervals from which the mean is calculated
% intervals - each row contains an interval (format: [start end] in ms)
% in which the mean is to be calculated
% (requires field 't' in dat)
%
%Returns:
% dat - updated data structure
%
% Taken from the BBCI Toolbox distributed with a MIT License
% https://github.com/bbci/bbci_public
%
% Copyright (c) 2015 The BBCI Group
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
% SOFTWARE.
if nargin==0
dat=[];return
end
[T, nChans, nMotos]= size(dat.x);
if length(nSamples)==1,
if ~exist('nMeans','var'), nMeans= floor(T/nSamples); end
dat.x =permute(mean(reshape(dat.x((T-nMeans*nSamples+1):T,:,:), ...
[nSamples,nMeans,nChans,nMotos]),1),[2 3 4 1]);
if isfield(dat, 'fs'),
dat.fs= dat.fs/nSamples;
end
if isfield(dat, 't'),
dat.t = mean(dat.t(reshape((T-nMeans*nSamples+1):T,nSamples,nMeans)));
end
elseif size(nSamples,1)==1 && size(nSamples,2)~=2,
intervals= nSamples([1:end-1; 2:end]');
dat= proc_jumpingMeans(dat, intervals);
else
nMeans= size(nSamples ,1);
da = zeros(nMeans, nChans, nMotos);
for i = 1:size(nSamples,1),
if any(isnan(nSamples(i,:))),
da(i,:,:) = NaN;
else
I = find(dat.t>=nSamples(i,1) & dat.t<=nSamples(i,2));
da(i,:,:) = mean(dat.x(I,:,:),1);
end
end
dat.x = da;
dat= rmfield(dat, 'fs');
dat.t = mean(nSamples, 2)';
end
function C = simulate_MIX(C,data,stimuli,sequence,y,nr_ept,A,gamma,trials,verbose)
%% This function simulates/trains an unsupervised mixing classifier
%
% Terminology:
% The series of flashes which is used to spell one character is called a
% trial. An epoch is called the segmented interval around a highlighting event.
% Furthermore, the following names are used:
%
% nr_commands: -- Number of different possible commands (e.g. number of
% symbols the user could spell in an ERP paradigm
% N -- Total number of epochs
% nr_ept -- Number of epochs per trial
% nr_tpt -- Number of targets per trial
% nr_trials -- Number of trials
% feat_dim -- Feature dimensionality
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input / Data Format:
%
% C -- Pre-estimated classifier, can be empty ([])
% data -- Array of size (N X feat_dim) feature matrix
% stimuli -- Array of size (nr_commands X N) matrix, BINARY ARRAY
% to indicate whether command was intensified during stimulus s_i
% sequence -- ARRAY of length N indicating whether a
% stimulus k belongs to sequence 1 or 2
% y -- [OPTIONAL]. BOOLEAN ARRAY of length N indicating whether
% stimulus k was a target (TRUE) or non-target (FALSE)
% Only used for statistics, not for training the classifier
% nr_ept -- Number of epochs (highlighting events) to spell one character
% A -- Mixing matrix for LLP
% gamma -- Value deciding about the decoder:
% -1 : MIX method
% 0 : EM-Method
% +1 : LLP-Method
% trials -- Which trials should be analysed, e.g. trials = 10:20
% analyses the spelling attempts for characters 10 to 20
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Background:
% The Expectation-Maximization (EM) Algorithm was implemented by Pieter-Jan Kindermans:
%
% Kindermans, P. J., Verstraeten, D., & Schrauwen, B. (2012).
% A bayesian model for exploiting application constraints to enable
% unsupervised training of a P300-based BCI. PloS one, 7(4), e33758.
%
% The learning from label proportion (LLP) idea is based on work by David Hübner et al.:
%
% Hübner, D., Verhoeven, T., Schmid, K., Müller, K. R.,
% Tangermann, M., & Kindermans, P. J. (2017).
% Learning from label proportions in brain-computer interfaces:
% Online unsupervised learning with guarantees. PloS one, 12(4), e0175856.
%
% The MIX-Method was implemented by Thibault Verhoeven:
%
% Verhoeven, T., Hübner, D., Tangermann, M., Müller, K. R.,
% Dambre, J., & Kindermans, P. J. (2017).
% Improving zero-training brain-computer interfaces by mixing model estimators.
% Journal of neural engineering, 14(3), 036021.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This function was created by
% David Hübner <david.huebner@blbt.uni-freiburg.de>
%
% Jul, 2017
if(nargin<10)
verbose = true;
end
% Adjustable parameters
nr_em_steps = 5; % Number of EM-steps [Default = 5]
% Internal parameters
nr_tpt = 16; % Number of targets per trial
A_inv = inv(A);
A_inv2 = A_inv.^2;
nr_commands = size(stimuli,1);
[nr_total,feat_dim] = size(data);
% Select specified trials in case this option was chosen
if ~isempty(trials)
idx = [];
for t=trials
idx = [idx (t-1)*nr_ept+1:t*nr_ept];
end
data = data(idx,:);
if ~isempty(y)
y = y(:,idx);
end
stimuli = stimuli(:,idx);
nr_trials = length(trials);
else
assert( rem(nr_total,nr_ept)==0, 'Number of data points divided by epochs per trial is not an integer');
nr_trials = nr_total/nr_ept;
trials = 1:nr_trials;
end
if(verbose)
disp('-------Start simulated online experiment!--------')
disp('The classifier starts from a random initalization and is retrained after each trial.');
disp('The binary target-vs-non target area under the curve (AUC) on all data ');
disp('until the current trial is then reported, if label information are available.');
disp('Additionally, the mixing coefficient for the target and non-target classes are reported.')
fprintf('-----------------------------------------------------\n\n');
end
% Init classifier
start_idx = 1;
if(isempty(C))
C = init_classifier(randn(feat_dim,1),feat_dim,nr_commands,nr_ept,nr_tpt);
else
% If we get classifier passed in, assume the caller is already accumulating data
% in trials and that we just want to integrate everything so far
start_idx = nr_trials;
end
% Loop over new trials
for c_i = start_idx:nr_trials
if(verbose)
fprintf('Trial %3d. ', trials(c_i));
end
tic();
N = c_i*nr_ept;
% Clear stored data fields (this is necessary becasuse the data is
% normalized on all data until the current trial)
C.classifier = clear_fields(C.classifier,feat_dim);
% Perform global normalization on all data until current trial
current_data = data(1:N,:);
current_stim = stimuli(:,1:N);
current_sequence = sequence(1:N,1);
global_mu = mean(current_data,1);
global_std = std(current_data,1,1);
data_norm = current_data-repmat(global_mu,N,1);
data_norm = data_norm./repmat(global_std,N,1);
% Split data into trials
data_list = {};
stimuli_list = {};
for j = 1:c_i
data_list{end+1} = data_norm((j-1)*nr_ept+1:j*nr_ept,:);
stimuli_list{end+1} = current_stim(:,(j-1)*nr_ept+1:j*nr_ept);
end
% Add all letters
C.classifier = classifier_mix_add_letter(C.classifier,data_list,stimuli_list);
% -----Determine covariance matrix with Ledoit-Wolf shrinkage ------ %
% as implemented for BCI by B. Blankertz et al., NeuroImage, 2010.
% Yields an analytical shrinkage parameter lamb and the regularized
% covariance matrix C.classifier.Sigma
nu = mean(diag(C.classifier.XTX));
T = nu*eye(feat_dim,feat_dim);
numerator = sum(sum(C.classifier.X2TX2-1.0*((C.classifier.XTX.^2)/N)));
denominator = sum(sum((C.classifier.XTX-T).^2));
lamb = (N/(N-1))*numerator/denominator;
lamb = max([0,min([1,lamb])]);
C.classifier.Sigma = ((1-lamb)*C.classifier.XTX+lamb*T)/(N-1);
% Determine whitening matrix
[V,D] = eig(C.classifier.Sigma);
C.classifier.Whiten = V*(diag(diag(D).^(-0.5)));
% ----------- Estimate Means with LLP ----------- %
N1 = sum(current_sequence==1);
N2 = sum(current_sequence==2);
% Averaging sequences
average_O1 = mean(data_norm(current_sequence==1,:),1);
average_O2 = mean(data_norm(current_sequence==2,:),1);
% Solution of the linear system
llp_mean_pos =(A_inv(1,1)*average_O1+A_inv(1,2)*average_O2)';
llp_mean_neg =(A_inv(2,1)*average_O1+A_inv(2,2)*average_O2)';
% Calculate variance for mixing it with EM
var_O1 = var(data_norm(current_sequence==1,:)*C.classifier.Whiten,1,1);
var_O2 = var(data_norm(current_sequence==2,:)*C.classifier.Whiten,1,1);
llp_var_pos = sum((A_inv2(1,1)/N1)*var_O1 + (A_inv2(1,2)/N2)*var_O2);
llp_var_neg = sum((A_inv2(2,1)/N1)*var_O1 + (A_inv2(2,2)/N2)*var_O2);
% ----------- Estimate Means with EM ----------- %
% EM Iterations
for em_it = 1:nr_em_steps
C.classifier = classifier_mix_expectation(C.classifier);
C.classifier = classifier_mix_maximization(C.classifier);
end
C.classifier = classifier_mix_expectation(C.classifier);
%Extract EM means
em_mean_pos = C.classifier.em_pos;
em_mean_neg = C.classifier.em_neg;
% ----------- Mix Mean Estimations ----------- %
% gamma = -1 ==> MIX means
% gamma = 0 ==> EM-algorithm
% gamma = 1 ==> LLP-algorithm
if gamma == -1
[EM_var_pos, EM_var_neg] = compute_EM_variance(C,data_norm,em_mean_pos,em_mean_neg);
pos_est_diff = C.classifier.Whiten'*(em_mean_pos - llp_mean_pos);
neg_est_diff = C.classifier.Whiten'*(em_mean_neg - llp_mean_neg);
gamma_pos = max([0,min([1,0.5*((EM_var_pos-llp_var_pos)/dot(pos_est_diff',pos_est_diff)+1)])]);
gamma_neg = max([0,min([1,0.5*((EM_var_neg-llp_var_neg)/dot(neg_est_diff',neg_est_diff)+1)])]);
else
gamma_pos = gamma;
gamma_neg = gamma;
end
%Mix EM means with LLP based on mixing coefficients gamma
mean_pos = (1.0-gamma_pos)*em_mean_pos+gamma_pos*llp_mean_pos;
mean_neg = (1.0-gamma_neg)*em_mean_neg+gamma_neg*llp_mean_neg;
%Update projections with new means
C.classifier.w = C.classifier.Sigma\(mean_pos-mean_neg); %Sigma^(-1)*(mu_pos-mu_neg)
%******Evaluate******
projection = classifier_compute_projection(C.classifier,C.classifier.data); % x^T*w
% Compute AUC if labels are given
if ~isempty(y)
label = [y; 1-y];
total_projection = vertcat(projection{:})';
auc=loss_rocArea([y(1:length(total_projection)); ...
1-y(1:length(total_projection))],total_projection);
C.statistics.auc = [C.statistics.auc auc];
if(verbose)
fprintf('AUC: %.3f%%. ',100*auc);
end
end
if(verbose)
fprintf('Gamma_pos: %.3f, Gamma_neg: %.3f. Runtime: %.3fs\n', gamma_pos, gamma_neg, toc());
end
% Store informative results in C.statistics
C.statistics.projection{end+1} = projection;
C.statistics.data_log_likelihood = [C.statistics.data_log_likelihood C.classifier.data_log_likelihood];
C.statistics.w{end+1} = C.classifier.w;
C.statistics.probs{end+1} = C.classifier.probs;
C.statistics.gamma_pos = [C.statistics.gamma_pos gamma_pos];
C.statistics.gamma_neg = [C.statistics.gamma_neg gamma_neg];
end
if(verbose)
fprintf('\n-----------------------------------------------------\n');
disp('Simulation completed.')
disp('Find the final classifier in C.classifier.');
disp('Statistics are saved in C.statistics.');
fprintf('-----------------------------------------------------\n\n');
end
function res = mysplitapply( fun, data, groups )
% mysplitapply applies function fun to rows of data specified by groups
% [length(groups) size(data,1)]
assert(length(groups)==size(data,1));
u = unique(groups);
tmp = fun(data(groups==u(1),:));
res = zeros(length(u), length(tmp));
res(1,:) = tmp;
for i=2:length(u)
res(i,:) = fun(data(groups==u(i),:));
end
end
% Initialize.m
% -------------------------------
% Author : Jussi T. Lindgren / Inria
% Date :
%
%
function box_out = p300_Initialize(box_in)
% profile on;
addpath('matlab');
addpath('matlab/helpers');
addpath('matlab/EM');
% We carry settings in this structure instead of using global variables
box_in.user_data.flashes_per_trial = box_in.settings(1).value;
box_in.user_data.debug = box_in.settings(2).value;
box_in.user_data.max_trials = 15; % max number of trials retained for building the classifier
% and how many past letters can be affected
box_in.user_data.keyboard_cols = 7; % @FIXME bad hardcoding
box_in.user_data.previous_time = 0;
box_in.user_data.gfeatures = [];
box_in.user_data.oldSelections = []; % set of letter predictions from before max_trials was reached
box_in.user_data.total_flashes = 0;
box_in.user_data.total_trials = 0;
box_in.user_data.C = []; % classifier
% fprintf(1,'%d\n', box_in.settings(1).value);
box_out = box_in;
end
\ No newline at end of file
% p300_Process.m
% -------------------------------
% Author : Jussi T. Lindgren / Inria
% Date :
%
%
function box_out = p300_Process(box_in)
% to use some global variables, we need to declare them
global OVTK_StimulationId_SegmentStart;
global OVTK_StimulationId_Label_01;
global OVTK_StimulationId_Label_08;
global OVTK_StimulationId_SegmentStop;
haveChunk = (OV_getNbPendingInputChunk(box_in,1)>0) && ...
(OV_getNbPendingInputChunk(box_in,2)>0) && ...
(OV_getNbPendingInputChunk(box_in,3)>0);
% only proceed after we have gathered one of each inputs
if(~haveChunk)
box_out = box_in;
return;
end
% As the current EM code gets slow when enough data is accumulated,
% here we always work just 'maxFeatures' last flashes worth of data.
maxFeatures = box_in.user_data.max_trials*box_in.user_data.flashes_per_trial;
[box_in, start_time, end_time, features] = OV_popInputBuffer(box_in,1);
[box_in, ~, ~, flashed] = OV_popInputBuffer(box_in,2);
[box_in, ~, ~, sequence] = OV_popInputBuffer(box_in,3);
if(size(features,2)==1 && size(features,1)>1)
% fix orientation
features = features';
flashed = flashed';
sequence = sequence';
end
if(size(features,1)>1 && size(features,2)>1)
% change to [time x chns]
features = features';
% if we get in a chunk of data, assume its an epoch, turn it into a
% feature vector by just taking in 6 equally sized groups. In
% OV, data matrix is [chns x samples].
if(0)
groupIdxs = round(linspace(1,6,size(features,1)))';
if(exist('splitapply','file'))
features = splitapply(@mean, features, groupIdxs);
else
features = mysplitapply(@mean, features, groupIdxs);
end
% change to row vector
features = features(:)';
else
baselineDuration = 200/1000; % 0.2sec baseline
a = [50 120; 121 200; 201 280;281 380;381 530; 531 700]/1000; % ms
a = floor((baselineDuration+a)*100); % -> sampleIdx @fixme assume hz=100;
baselineStop = round(baselineDuration*100);
% baseline subtraction
baseline = mean(features(1:baselineStop,:),1);
features = features - repmat(baseline,[size(features,1) 1]);
% save('E:\jl\debug.mat');
% fprintf(1,'%d %d\n', size(baseline,1), size(baseline,2));
features = [ mean(features(a(1,1):a(1,2),:),1), ...
mean(features(a(2,1):a(2,2),:),1), ...
mean(features(a(3,1):a(3,2),:),1), ...
mean(features(a(4,1):a(4,2),:),1), ...
mean(features(a(5,1):a(5,2),:),1), ...
mean(features(a(6,1):a(6,2),:),1) ...
];
end
end
if(isempty(box_in.user_data.gfeatures))
box_in.user_data.gfeatures = zeros(maxFeatures,size(features,2));
box_in.user_data.gflashed = zeros(maxFeatures,size(flashed,2));
box_in.user_data.gsequence = zeros(maxFeatures,size(sequence,2));
end
if(box_in.user_data.total_flashes >= maxFeatures)
% treat as a fifo
box_in.user_data.gfeatures = circshift(box_in.user_data.gfeatures,[-1 0]);
box_in.user_data.gflashed = circshift(box_in.user_data.gflashed,[-1 0]);
box_in.user_data.gsequence = circshift(box_in.user_data.gsequence,[-1 0]);
end
box_in.user_data.total_flashes = box_in.user_data.total_flashes + 1;
% box_in.user_data
idx = min(size(box_in.user_data.gfeatures,1),box_in.user_data.total_flashes);
box_in.user_data.gfeatures(idx,:) = features;
box_in.user_data.gflashed(idx,:) = flashed;
box_in.user_data.gsequence(idx,:) = sequence;
time_now = end_time;
if(rem(box_in.user_data.total_flashes,box_in.user_data.flashes_per_trial)==0)
% we've collected features for the whole trial, update model
box_in.user_data.total_trials = box_in.user_data.total_trials + 1;
% @fixme, hardcoded
gamma = -1;
A = [3/8 5/8; 2/18 16/18];
numTrials = min(box_in.user_data.total_trials,box_in.user_data.max_trials);
box_in.user_data.C = simulate_MIX(box_in.user_data.C, ...
box_in.user_data.gfeatures, logical(box_in.user_data.gflashed'), ...
box_in.user_data.gsequence, ...
[],box_in.user_data.flashes_per_trial,A,gamma,1:numTrials,false);
% predict which letter was selected
[~,select]=max(box_in.user_data.C.classifier.probs,[],2);
% Include predictions from before maxFeatures was reached
box_in.user_data.selection = select;
select = [box_in.user_data.oldSelections;select];
if(box_in.user_data.debug)
fprintf(1,'Trial %03d end (at %f s, %d flashes total, %d trls used, fd %d): Spelled so far: %s\n', ...
box_in.user_data.total_trials, end_time, ...
box_in.user_data.total_flashes, numTrials, size(features,2), convert_position_to_letter(select));
end
% Convert select to stimulations
select = select - 1;
rowidx = floor(select / box_in.user_data.keyboard_cols);
colidx = rem(select, box_in.user_data.keyboard_cols);
stim_set = [double(OVTK_StimulationId_SegmentStart); time_now; 0];
for q=1:length(rowidx)
stim_set = [stim_set, [rowidx(q) + double(OVTK_StimulationId_Label_01); time_now; 0]];
stim_set = [stim_set, [colidx(q) + double(OVTK_StimulationId_Label_08); time_now; 0]];
end
stim_set = [stim_set, [double(OVTK_StimulationId_SegmentStop); time_now; 0]];
box_in = OV_addOutputBuffer(box_in,1,box_in.user_data.previous_time,time_now,stim_set);
if(box_in.user_data.total_flashes >= maxFeatures)
box_in.user_data.oldSelections = [box_in.user_data.oldSelections;box_in.user_data.selection(1)];
end
else
box_in = OV_addOutputBuffer(box_in,1,box_in.user_data.previous_time,time_now,[]);
end
box_in.user_data.previous_time = time_now;
box_out = box_in;
end
% p300_Uninitialize.m
% -------------------------------
% Author : Jussi T. Lindgren / Inria
% Date :
%
function box_out = p300_Uninitialize(box_in)
rmpath('matlab/helpers');
rmpath('matlab/EM');
rmpath('matlab');
box_out = box_in;
end
\ No newline at end of file
To run test scenarios, export test signals here from matlab using
matlab-test/dump_data_to_csv.m