%% === Core Analysis Steps ===


%% Training and cross-validation
% load the training data
traindata = io_loadset('data:/imag.set');

% define approach
myapproach = 'CSP';

% train model and run cross-validation
[trainloss,lastmodel,laststats] = bci_train('Data',traindata,'Approach',myapproach,'TargetMarkers',{'S  1','S  2'}); 
disp(['training mis-classification rate: ' num2str(trainloss*100,3) '%']);

% visualize model
bci_visualize(lastmodel);



%% Testing
% load test set
testdata = io_loadset('data:/imag2.set');

% calculate predictions
[prediction,loss,teststats,targets] = bci_predict(lastmodel,testdata);
disp(['test mis-classification rate: ' num2str(loss*100,3) '%']);
disp(['  predicted classes: ',num2str(round(prediction{2}*prediction{3})')]);  % class probabilities * class values
disp(['  true classes     : ',num2str(round(targets)')]);

% annotate data with BCI-derived outputs
testdata_bci = bci_annotate(lastmodel,testdata);

% take a quick look at the data
testdata_bci.data = diag(1./std(testdata_bci.data,[],2))*testdata_bci.data;
vis_artifacts(testdata_bci);



%% Basic Online Processing

% stream dataset in background ("playback")
run_readdataset('laststream',testdata);

% start a prediction background task, updating a figure
run_writevisualization(lastmodel,'laststream','bar(y);ylim([0 1])');

% run until user clicks into the figure
pause(5);
waitforbuttonpress;
onl_clear;
close(gcf);


%% Inspect model and statistics

% most important stats properties
laststats
laststats.mcr_mu, laststats.mcr_std  % mean and standard deviation of mis-classification rate
laststats.TP, laststats.FP           % true-positive and false-positive rate
[laststats.per_fold.mcr]             % error rate in each fold
vertcat(laststats.per_fold.targ)     % target vales in each fold (here chronological)
laststats.per_fold(1).pred           % predictions on first fold
laststats.per_fold(1).pred{2}        % probability of class 1 and class 2 for each trial in 1st fold
tmp = vertcat(laststats.per_fold.pred); vertcat(tmp{:,2}) % probabilities for both classes across trials

% some model properties
lastmodel
lastmodel.args                       % arguments passed in to the BCI paradigm's calibrate() function
% feature extraction properties
lastmodel.featuremodel               % model returned by the paradigm's feature_adapt function, if any (feature-extraction model)
lastmodel.featuremodel.filters       % in case of CSP: [Channels x Filters] array of spatial filters 
lastmodel.featuremodel.chanlocs      % channel locations at feature-extraction stage, e.g., for topo plotting
figure;topoplot_grid(lastmodel.featuremodel.filters,lastmodel.featuremodel.chanlocs) % spatial filters of the model
figure;topoplot_grid(lastmodel.featuremodel.patterns,lastmodel.featuremodel.chanlocs) % spatial patterns (forward projections) of the model
lastmodel.conditioningmodel          % rarely-used step between feature extraction and machine learning (e.g., dimensionality reduction)
% machine learning properties
lastmodel.predictivemodel.args       % arguments passed into the machine learning function
lastmodel.predictivemodel.model      % model produced by the machine learning function (here: LDA)
lastmodel.predictivemodel.model.w    % vector of linear classifier weights, one element per feature (some also contain an extra element for the bias)
figure;bar(lastmodel.predictivemodel.model.w) % linear weights for each spatially filtered signal
lastmodel.predictivemodel.model.b    % bias value for the classifier (not supported by all classifiers)
% signal processing properties
lastmodel.tracking                   % information for online application of the model
lastmodel.tracking.filter_graph      % cell array of filter chains/trees, one per input stream of the paradigm's prediction function (almost always just 1 stream)
lastmodel.tracking.filter_graph{1}   % recursive data structure describing filter operations yielding the paradigm's first input stream during online processing
exp_fullform(lastmodel.tracking.filter_graph{1}) % string version of 1st filter chain/tree 
stage=utl_find_filter(lastmodel.tracking.filter_graph{1},'flt_fir') % properties of the FIR filter stage used in this model
stage.parts                          % cell array of arguments used in FIR filter application
tmp = hlp_varargin2struct(stage.parts); % struct version of FIR input arguments
tmp.fspec                               % frequency bands used in FIR filter
tmp.signal                              % input signal to the FIR stage (here: resampling)
figure;plot(tmp.state.b)                % filter kernel used by FIR filter (here: minimum-phase bandpass)
hlp_varargin2struct(getfield(utl_find_filter(lastmodel,'set_makepos'),'parts')) % arguments to epoch extraction
% input data to reproduce the model
lastmodel.source_data                % collection of source data sets used to compute the model
lastmodel.source_data{1}.streams{1}  % expression to reproduce the first stream of the first data set
exp_eval(lastmodel.source_data{1}.streams{1}) % input data used to compute the model, reloaded from disk
getfield(utl_find_filter(lastmodel.source_data{1}.streams{1},'io_loadset'),'parts') % arguments to load command
lastmodel.control_options            % options used for bci_train
lastmodel.control_options.approach   % approach and parameters used to compute the model


%% === Defining Customized Approaches ===
% set up some fixed inputs to bci_train for use in below lines
params = {'Data',traindata,'TargetMarkers',{'S  1','S  2'}};

%% using a different epoch size
app = {'CSP','SignalProcessing',{'EpochExtraction',[0.5 2.5]}};
fprintf('Error rate with shorter epoch: %.3f%%\n',bci_train('Approach',app,params{:})*100);

%% using a different frequency band
app = {'CSP','SignalProcessing',{'FIRFilter',[6 7 15 17]}};
fprintf('Error rate with alpha band: %.3f%%\n',bci_train('Approach',app,params{:})*100);

%% using a different number of CSP patterns for feature extration (in the prediction function)
app = {'CSP','Prediction',{'FeatureExtraction',{'PatternPairs',2}}};
fprintf('Error rate with 2 pattern pairs: %.3f%%\n',bci_train('Approach',app,params{:})*100);

%% using a different classifier (in the prediction function)
app = {'CSP','Prediction',{'MachineLearning',{'Learner','logreg'}}};
fprintf('Error rate with logistic regression (Variational Bayesian): %.3f%%\n',bci_train('Approach',app,params{:})*100);

%% using a sparse classifier variant
app = {'CSP','Prediction',{'MachineLearning',{'Learner',{'logreg','Variant','vb-ard'}}}};
fprintf('Error rate with sparse logistic regression: %.3f%%\n',bci_train('Approach',app,params{:})*100);

%% searching over the optimal number of patterns (note performs a nested cross-validation)
app = {'CSP','Prediction',{'FeatureExtraction',{'PatternPairs',search(1,2,3)}}};
fprintf('Error rate with optimized pattern number: %.3f%%\n',bci_train('Approach',app,params{:})*100);

%% using a different BCI paradigm altogether, also visualize model
app = {'SpecCSP'};
[loss,model] = bci_train('Approach',app,params{:});
fprintf('Error rate with Spec-CSP: %.3f%%\n',loss*100);
bci_visualize(model);

%% using Spec-CSP with a sparse non-linear classifier (note: displays "all folds" 6x, takes a while)
app = {'SpecCSP','Prediction',{'MachineLearning',{'Learner',{'hkl'}}}};
fprintf('Error rate with Spec-CSP and HKL: %.3f%%\n',bci_train('Approach',app,params{:})*100);

%% using a heavily customized approach that uses the (log-) power spectral density 
%  in standard bands of independent components and a sparse classifier
app = {'Spectralmeans', ...
    'SignalProcessing',{ ...
        'ICA',{...
            'Variant','robust_sphere', ...
            'DataCleaning','highpass', ...
            'TransformData',true}}, ...
    'Prediction',{ ...
        'FeatureExtraction',{ ...
            'FreqWindows',[0.5 3; 4 7; 8 13; 14 30; 31 42], ...
            'LogTransform',true}, ...
        'MachineLearning',{'Learner',{'logreg','Variant','lars'}}}};
fprintf('Error rate with ICA-bandpower and LARS: %.3f%%\n',bci_train('Approach',app,params{:})*100);




%% === Customized Analysis Steps ===


%% Data curation example I: using EEGLAB tools
% load subset of channels, evaluate it for editing
smalldata = exp_eval(io_loadset('data:/imag.set','channels',{'C3','C4','Cz','FCz','Oz'}));
% relabel events to 'left'/'right'
[smalldata.event.type] = celldeal(strrep(strrep({smalldata.event.type},'S  1','left'),'S  2','right'));
% use middle third of the data for training and last third of the data for testing
first_part = pop_select(smalldata,'time',[smalldata.xmax*1/3 smalldata.xmax*2/3]);
second_part = pop_select(smalldata,'time',[smalldata.xmax*2/3+1 smalldata.xmax]);

% train model and run cross-validation on first part, then test on second part
[trainloss,lastmodel] = bci_train('Data',first_part,'Approach','SpecCSP','TargetMarkers',{'left','right'}); 
[prediction,testloss] = bci_predict(lastmodel,second_part);
disp(['training mis-classification rate: ' num2str(trainloss*100,3) '%']);
disp(['test mis-classification rate: ' num2str(testloss*100,3) '%']);


%% Data curation example II: predicting time-on-task

% mass-insert 300 markers at regular intervals, excluding first and final minute of data
tmpdata = exp_eval(set_insert_markers('Signal',traindata, ...
    'SegmentSpec',{60,smalldata.xmax-60},'Event','myevent','Count',300));
% manually set each marker's .target value to the current time, remove all other markers (including 'boundary')
[tmpdata.event.target] = celldeal(num2cell([tmpdata.event.latency]/tmpdata.srate));
tmpdata.event(~strcmp({tmpdata.event.type},'myevent')) = [];


% define approach that can perform regression
myapproach = {'Spectralmeans', ...
    'SignalProcessing',{'EpochExtraction',[-2 2]},...
    'Prediction',{'MachineLearning',{'Learner',{'logreg','Variant','lars','Regression',true}}}};

% evaluate, quantify performance using Mean Absolute Error (mae) metric
[trainloss,lastmodel,laststats] = bci_train('Data',tmpdata,'Approach',myapproach,'Metric','mae'); 
fprintf('training error: %.1f seconds (+/-%.1fs); total duration: %.1fs\n',trainloss,laststats.mae_std,tmpdata.xmax);

% plot the predictions versus the true time on task
figure; hold; xlabel('marker #'); ylabel('time on task (s)');
plot(vertcat(laststats.per_fold.pred));
plot(vertcat(laststats.per_fold.targ),'r--');


%% Different Cross-Validation Schemes (using previous regression approach)

%% use 10-fold chronological (=blockwise) cross-validation (default is 5-fold chronological)
[trainloss,lastmodel,laststats] = bci_train('Data',tmpdata,'Approach',myapproach, ...
    'Metric','mae','EvaluationScheme',{'chron',10}); 
fprintf('error with 10x chronological CV: %.1f seconds (+/-%.1fs); total duration: %.1fs\n',trainloss,laststats.mae_std,tmpdata.xmax);
figure; hold; xlabel('marker #'); ylabel('time on task (s)'); plot(vertcat(laststats.per_fold.pred)); plot(vertcat(laststats.per_fold.targ),'r--');

%% use 10-fold randomized cross-validation
[trainloss,lastmodel,laststats] = bci_train('Data',tmpdata,'Approach',myapproach, ...
    'Metric','mae','EvaluationScheme',10); 
fprintf('error with 10x randomized CV: %.1f seconds (+/-%.1fs); total duration: %.1fs\n',trainloss,laststats.mae_std,tmpdata.xmax);
% need to sort the trials for plotting
[targets,order] = sort(vertcat(laststats.per_fold.targ));
predictions = vertcat(laststats.per_fold.pred);
predictions = predictions(order);
figure; hold; xlabel('marker #'); ylabel('time on task (s)'); plot(predictions); plot(targets,'r--');

%% use 10-fold chronological CV with 50 trials margin between train and test (omitted from training set)
[trainloss,lastmodel,laststats] = bci_train('Data',tmpdata,'Approach',myapproach, ...
    'Metric','mae','EvaluationScheme',{'chron',10,50}); 
fprintf('error with 10x chronological CV: %.1f seconds (+/-%.1fs); total duration: %.1fs\n',trainloss,laststats.mae_std,tmpdata.xmax);
figure; hold; xlabel('marker #'); ylabel('time on task (s)'); plot(vertcat(laststats.per_fold.pred)); plot(vertcat(laststats.per_fold.targ),'r--');


%% Custom Prediction Parameters

%% using mean absolute error for prediction, predicting on training session
[prediction,testloss,teststats,targets] = bci_predict('Model',lastmodel,'Data',tmpdata,'Metric','mae');
fprintf('error of whole-data trained classifier on training set (illegal): %.1f seconds\n',testloss);
figure; hold; xlabel('marker #'); ylabel('time on task (s)'); plot(prediction); plot(targets,'r--');

%% on second session
% curate test set
tmpdata2 = exp_eval(testdata);
tmpdata2 = exp_eval(set_insert_markers('Signal',tmpdata2,'SegmentSpec',{60,tmpdata2.xmax-60},'Event','myevent2','Count',300));
[tmpdata2.event.target] = celldeal(num2cell([tmpdata2.event.latency]/tmpdata2.srate));
tmpdata2.event(~strcmp({tmpdata2.event.type},'myevent2')) = [];

% perform prediction
[prediction,testloss,teststats,targets] = bci_predict('Model',lastmodel,'Data',tmpdata2,'Metric','mae');
fprintf('error on second session: %.1f seconds\n',testloss);
figure; hold; xlabel('marker #'); ylabel('time on task (s)'); plot(prediction); plot(targets,'r--');


%% Corollary: Attempt to distinguish data from first vs. second session

% merge both sessions, remove every second marker (to reduce training time)
both_sessions = exp_eval(set_merge(tmpdata,tmpdata2));
both_sessions.event = both_sessions.event(1:2:end);

% use Spec-CSP, 10x chronological with 25 trials margin, markers myevent vs myevent2
[trainloss,lastmodel,laststats] = bci_train('Data',both_sessions,'Approach','SpecCSP', ...
    'EvaluationScheme',{'chron',10,25},'TargetMarkers',{'myevent','myevent2'}); 
disp(['error rate at session classification: ' num2str(trainloss*100,3) '%']);
% display session classifier
bci_visualize(lastmodel);




%% === Custom Online Scripting ===


%% Using custom main loop for online processing

% train a model on training set (using Spec-CSP again, left vs. right)
[loss,model] = bci_train('Data',io_loadset('data:/imag.set'),'Approach','SpecCSP','TargetMarkers',{'S  1','S  2'}); 

%% load test data for real-time playback of recording
testdata = exp_eval(io_loadset('data:/imag2.set'));
pos = 0; t_start = tic;  % start time for main loop

% define new online stream and set up new online predictor
onl_newstream('mystream', 'srate',testdata.srate, 'chanlocs',testdata.chanlocs);
onl_newpredictor('mypredictor',lastmodel);
while true
    % get current position into data set for playback
    last_pos = pos;
    pos = 1+round(toc(t_start)*testdata.srate);
    if pos>testdata.pnts
        break; end
    
    % extract new chunk from testdata
    chunk = testdata.data(:,last_pos+1:pos);

    % append chunk to online stream and query classifier
    onl_append('mystream',chunk);
    y = onl_predict('mypredictor');
    
    % visualize classifier output
    if iscell(y)
        bar(y{2}); 
        ylim([0 1]); 
        drawnow; 
    end
end


%% Simulate online processing in non-real-time

%% call onl_predict at 3.5 seconds after each marker
[predictions,latencies] = onl_simulate(testdata,model,'offset',3.5,'markers',{'S  1','S  2'});
% get known target labels and manually calculate the error rate
targets = strcmp({testdata.event.type},'S  1') + strcmp({testdata.event.type},'S  2')*2;
fprintf('pseudo-online error rate: %.1f%%\n',100*mean(argmax(predictions')~=targets(targets~=0)));

%% call onl_predict at a rate of 5 Hz, covering the first 5 minutes of data
[predictions,latencies] = onl_simulate(testdata,model,'sampling_rate',5,'interval',[0 300]);
figure('Position',[0 0 1200 400]); hold; plot(latencies,predictions(:,2));



%% === Batch Analysis ===

%% define data and approaches
% define data to train on (and perform cross-validation on)
trainsets = {'data:/imag.set','data:/imag2.set'};
% define data to test on (here we test on the other session for lack of data)
testsets = {'data:/imag2.set','data:/imag.set'};
% define approaches
approaches.CSP_7to30 = {'CSP' 'SignalProcessing',{'FIRFilter',[6 8 28 32]}};
approaches.CSP_7to30_iir = {'CSP' 'SignalProcessing',{'FIRFilter','off','IIRFilter',[6 8 28 32]}};
approaches.CSP_alpha = {'CSP' 'SignalProcessing',{'FIRFilter',[6 8 14 15]}};
approaches.BP_channels = {'Bandpower','SignalProcessing',{'SurfaceLaplacian','off'}};
approaches.BP_surflap = {'Bandpower','SignalProcessing',{'SurfaceLaplacian','on'}};

%% Run basic batch analysis 
results = bci_batchtrain('Datasets',trainsets,'Approaches',approaches,'TargetMarkers',{'S  1','S  2'})
% investigate some results
[results.BP_surflap.loss]           % average losses per data set for approach "BP_surflap"
[results.CSP_7to30_iir.model]       % models per data set for CSP approach
[results.BP_channels.stats]         % all stats for BP_channels approach
BP_losses = [results.BP_surflap(2).stats.per_fold.mcr] % mis-classification rates on all 5 folds of second data set, for BP-surflap approach
CSP_losses = [results.CSP_7to30(2).stats.per_fold.mcr]
ttest(CSP_losses,BP_losses,0.01)    % paired-samples t-test at 1% significance level between CSP and BP
signrank(CSP_losses,BP_losses)      % p-value for paired-samples Wilcoxon signed rank test between CSP and BP


%% Run with sets to predict on, too
results = bci_batchtrain('Datasets',trainsets,'PredictSets',testsets,'Approaches',approaches,'TargetMarkers',{'S  1','S  2'})

%% Run with a directory of data sets 
% (and run a second time to observe speedup due to cached intermediate results)
datapath = env_translatepath('bcilab:/userdata/*.set')
results = bci_batchtrain('Datasets',datapath,'Approaches',approaches,'TargetMarkers',{'S  1','S  2'})

%% Run with custom arguments to io_loadset, bci_train, and bci_predict, respectively
results = bci_batchtrain('Datasets',trainsets,'PredictSets',testsets,'Approaches',approaches,'TargetMarkers',{'S  1','S  2'}, ...
    'LoadArguments',{'channels',{'C3','C4','Cz','FCz','Oz'}}, ...         % using small subset of channels
    'TrainArguments',{'EvaluationScheme',{'chron',3}, 'Metric','auc'},... % using negative area under curve and 3x chronological cross-validation
    'PredictArguments',{'Metric','auc'})                                  % also using negative area under curve

%% Run again, reusing existing results for given approaches and data sets
% (note that the results from the previous run, using a non-standard metric are looked up)
results = bci_batchtrain('Datasets',datapath,'Approaches',approaches,'TargetMarkers',{'S  1','S  2'}, ...
    'ReuseExisting',true)

%% Run again, this time declaring a custom "study tag" that allows to keep analyses separated on disk
% (note that this time no existing results are present at first run and the correct metric is used)
results = bci_batchtrain('Datasets',datapath,'Approaches',approaches,'TargetMarkers',{'S  1','S  2'}, ...
    'ReuseExisting',true,'StudyTag','ImagBenchmark_1.1')

%% Specify a custom storage path pattern (default is in the home directory)
results = bci_batchtrain('Datasets',datapath,'Approaches',approaches,'TargetMarkers',{'S  1','S  2'}, ...
    'ReuseExisting',true,'StudyTag','ImagBenchmark_1.1', ...
    'StoragePattern','bcilab:/results/batchtrain/%caller-%study/%approach-%set.mat')


%% Run with the default error handler (none)
approaches.buggy = {'sfdsdfsdf'};
results = bci_batchtrain('Datasets',datapath,'Approaches',approaches,'TargetMarkers',{'S  1','S  2'}, ...
    'ReuseExisting',true,'StudyTag','ImagBenchmark_1.1', ...
    'StoragePattern','bcilab:/results/batchtrain/%caller-%study/%approach-%set.mat')

%% Run using a custom error handler (note: this one will be set by default in future versions)
approaches.buggy = {'sfdsdfsdf'};
results = bci_batchtrain('Datasets',datapath,'Approaches',approaches,'TargetMarkers',{'S  1','S  2'}, ...
    'ReuseExisting',true,'StudyTag','ImagBenchmark_1.1','ErrorHandler','hlp_handleerror', ...
    'StoragePattern','bcilab:/results/batchtrain/%caller-%study/%approach-%set.mat')


%% Run on a pre-configured cluster 
% (WARNING: If you are on UNIX/linux, this will by default attempt to spawn 
% multiple MATLAB workers on your own machine, see also Settings/Cluster Settings/Edit...)

% acquire the cluster (same as pressing the cluster switch in the GUI)
env_acquire_cluster;

% run computation remotely if possible
results = bci_batchtrain('Datasets',datapath,'Approaches',approaches,'TargetMarkers',{'S  1','S  2'}, ...
    'ReuseExisting',true,'StudyTag','ImagBenchmark_1.2','ErrorHandler','hlp_handleerror', ...
    'StoragePattern','bcilab:/results/batchtrain/%caller-%study/%approach-%set.mat')

% release the cluster (workers eventually shut themselves down; same as releasing the cluster switch)
env_release_cluster;
