1% Benchmark for the BioSig toolbox
2%  The benchmark test performs a typical data analysis task.
3%  Except for data loading, it tests the performance of the computational speed
4%  and can be used to compare the performance of different platforms
5%
6%  Requirements:
7%  Octave 2.1 or higher, or Matlab
8%  BioSig4OctMat from http:/biosig.sf.net/
9%
10%  References:
11%  [1] Alois Schloegl, BioSig - An application of Octave, 2006
12%      available online: http://arxiv.org/pdf/cs/0603001v1
13
14
15%	$Id: bench_biosig.m,v 1.7 2009/04/21 10:39:20 schloegl Exp $
16%	Copyright (C) 2005,2006 by Alois Schloegl <alois.schloegl@gmail.com>
17%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/
18
19
20% This library is free software; you can redistribute it and/or
21% modify it under the terms of the GNU Library General Public
22% License as published by the Free Software Foundation; either
23% Version 3 of the License, or (at your option) any later version.
24%
25% This library is distributed in the hope that it will be useful,
26% but WITHOUT ANY WARRANTY; without even the implied warranty of
27% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
28% Library General Public License for more details.
29%
30% You should have received a copy of the GNU Library General Public
31% License along with this library; if not, write to the
32% Free Software Foundation, Inc., 59 Temple Place - Suite 330,
33% Boston, MA  02111-1307, USA.
34
35clear
36
37if ~exist('l1.gdf','file') % download test file if not available
38        system('wget http://pub.ist.ac.at/~schloegl/bci/bci7/l1.gdf');
39end;
40
41K=0; tic; t0=cputime();
42
43%[signal,HDR]=sload({[p,'x_train'],[p,'x_test']});
44[signal,HDR]=sload('l1.gdf');
45KK = 5;
46ng = ceil([0:length(HDR.Classlabel)-1]'/length(HDR.Classlabel)*KK);
47Classifiers = {'LD2','LD3','LD4','MDA','MD2','MD3','GRB','GRB2','QDA','MQU','INQ','Cauchy','SVM','SVM11','RBF','REG'};
48Classifiers = {'LDA','LDA/GSVD','LD2','LD3','LD4','NBC','aNBC','REG','MDA','MD2','MD3','GRB','GRB2','QDA','MQU','IMQ','Cauchy','LDA/sparse','SVM:LIB','SVM:OSU'};
49
50
51K=K+1;jo{K}='sload l1.gdf'; t(K)=cputime()-t0,t1(K)=toc
52z = zscore(signal);
53
54try
55        bp = bandpower(z,HDR.SampleRate);
56        bp(bp<-10)=NaN;
57        K=K+1;jo{K}='bandpower'; t(K)=cputime()-t0,t1(K)=toc
58end;
59
60[SIGMA1,PHI1,OMEGA1,m01,m11,m21] = wackermann(z(:,1:2),HDR.SampleRate);
61[SIGMA2,PHI2,OMEGA2,m02,m12,m22] = wackermann(z(:,2:3),HDR.SampleRate);
62K=K+1;jo{K}='wackermann'; t(K)=cputime()-t0,t1(K)=toc
63
64if 1,
65[a,f,s] = barlow(z,HDR.SampleRate);
66BARLOW = [a,f,s];
67K=K+1;jo{K}='barlow'; t(K)=cputime()-t0,t1(K)=toc
68
69[A,M,C] = hjorth(z,HDR.SampleRate);
70HJORTH = [A,M,C];
71K=K+1;jo{K}='hjorth'; t(K)=cputime()-t0,t1(K)=toc
72
73uc  = 30:5:80;
74[F] = tdp(z,3,2^(-uc(7)/8));
75K=K+1;jo{K}='TDP'; t(K)=cputime()-t0,t1(K)=toc
76
77a = [];
78for ch = 1:size(signal,2),
79        fprintf(1,'.%i',ch);
80        INI.MOP = [0,3,0];
81        INI.UC  = 2^(-uc(7)/8);
82
83        X = tvaar(signal(:,ch),INI);
84        X = tvaar(signal(:,ch),X);
85
86        a = [a,X.AAR];
87        K = K+1; jo{K} = ['aar #',int2str(ch)]; t(K)=cputime()-t0,t1(K)=toc
88end
89
90%% K-fold Crossvalidation
91
92try,
93        CC1 = findclassifier(bp,HDR.EVENT.POS(HDR.EVENT.TYP==hex2dec('300'))-1,[HDR.Classlabel,ng],reshape(1:1152,16,72)',[1:72]'>24,'LD3');
94        CC1.TSD.T = CC1.TSD.T/HDR.SampleRate;
95        K=K+1; jo{K}='findclassifier bp LD3'; t(K)=cputime()-t0,t1(K)=toc
96catch,
97end;
98
99CC2 = findclassifier(BARLOW,HDR.EVENT.POS(HDR.EVENT.TYP==hex2dec('300'))-1,[HDR.Classlabel,ng],reshape(1:1152,16,72)',[1:72]'>24,'LD3');
100CC2.TSD.T = CC2.TSD.T/HDR.SampleRate;
101K=K+1; jo{K}='findclassifier barlow LD3'; t(K)=cputime()-t0,t1(K)=toc
102CC3 = findclassifier(HJORTH,HDR.EVENT.POS(HDR.EVENT.TYP==hex2dec('300'))-1,[HDR.Classlabel,ng],reshape(1:1152,16,72)',[1:72]'>24,'LD3');
103CC3.TSD.T = CC3.TSD.T/HDR.SampleRate;
104K=K+1; jo{K}='findclassifier hjorth LD3'; t(K)=cputime()-t0,t1(K)=toc
105CC4 = findclassifier(a,HDR.EVENT.POS(HDR.EVENT.TYP==hex2dec('300'))-1,[HDR.Classlabel,ng],reshape(1:1152,16,72)',[1:72]'>24,'LD3');
106CC4.TSD.T = CC4.TSD.T/HDR.SampleRate;
107K=K+1; jo{K}='findclassifier aar LD3'; t(K)=cputime()-t0,t1(K)=toc
108CC5 = findclassifier([SIGMA1,PHI1,OMEGA1,SIGMA2,PHI2,OMEGA2],HDR.EVENT.POS(HDR.EVENT.TYP==hex2dec('300'))-1,[HDR.Classlabel,ng],reshape(1:1152,16,72)',[1:72]'>24,'LD3');
109CC5.TSD.T = CC5.TSD.T/HDR.SampleRate;
110end;
111
112K=K+1; jo{K}='findclassifier Wackermann LD3'; t(K)=cputime()-t0,t1(K)=toc
113CC6 = findclassifier([SIGMA1,PHI1,OMEGA1,SIGMA2,PHI2,OMEGA2],HDR.EVENT.POS(HDR.EVENT.TYP==hex2dec('300'))-1,[HDR.Classlabel,ng],reshape(1:1152,16,72)',[1:72]'>24,'LD3');
114CC6.TSD.T = CC6.TSD.T/HDR.SampleRate;
115K=K+1; jo{K}='findclassifier TDP LD3'; t(K)=cputime()-t0,t1(K)=toc
116
117
118fprintf(1,'Version: %s\ncputime\t toc[s]\ttask\n================================\n',version);
119t0 = [0,0];
120for k=1:K,
121	fprintf(1,'%7.3f\t%7.3f\t%s\n',t(k)-t0(1),t1(k)-t0(2),jo{k});
122	t0 = [t(k),t1(k)];
123end;
124fprintf(1,'-----------------------------------\n%7.3f\t%7.3f\t%s\n',t(k),t1(k),'total');
125return;
126
127clear ACC KAP MI AUC
128TI = 36; 35;27;
129for k=[1:18], %:17,18:length(Classifiers);
130        try
131        CC = findclassifier(a,HDR.EVENT.POS(HDR.EVENT.TYP==hex2dec('300'))-1,[HDR.Classlabel,ng],(1:16:1152)',[1:72]'==TI,Classifiers{k});
132        CC.TSD.T = CC.TSD.T/HDR.SampleRate;
133        K=K+1; jo{K}=['findclassifier AAR ',Classifier{k}]; t(K)=cputime()-t0,t1(K)=toc
134        CC90{k}=CC.TSD;
135        ACC(:,k)=CC90{k}.ACC00;
136        KAP(:,k)=CC90{k}.KAP00;
137        MI(:,k)=CC90{k}.I(:,1);
138        AUC(:,k)=CC90{k}.AUC(:,1);
139        catch
140        end;
141end;
142figure(1)
143subplot(221),
144plot(CC90{k}.T,ACC);
145subplot(222),
146plot(CC90{k}.T,KAP);
147subplot(223),
148plot(CC90{k}.T,MI);
149subplot(224),
150plot(CC90{k}.T,AUC);
151legend(Classifiers);
152
153
154figure(2)
155subplot(221),
156plot(ACC(500:700,:)');
157subplot(222),
158plot(KAP(500:700,:)');
159subplot(223),
160plot(MI(500:700,:)');
161subplot(224),
162plot(AUC(500:700,:)');
163legend(Classifiers);
164
165
166CC91 = findclassifier(a,HDR.EVENT.POS(HDR.EVENT.TYP==hex2dec('300'))-1,[HDR.Classlabel,ng],(1:16:1152)',[1:72]'==TI,'sparse_lda1');
167CC91.TSD.T = CC91.TSD.T/HDR.SampleRate;
168CC92 = findclassifier(a,HDR.EVENT.POS(HDR.EVENT.TYP==hex2dec('300'))-1,[HDR.Classlabel,ng],(1:16:1152)',[1:72]'==TI,'LD3');
169CC92.TSD.T = CC92.TSD.T/HDR.SampleRate;
170CC93 = findclassifier(a,HDR.EVENT.POS(HDR.EVENT.TYP==hex2dec('300'))-1,[HDR.Classlabel,ng],(1:16:1152)',[1:72]'==TI,'LD3/GSVD');
171CC93.TSD.T = CC92.TSD.T/HDR.SampleRate;
172CC94 = findclassifier(a,HDR.EVENT.POS(HDR.EVENT.TYP==hex2dec('300'))-1,[HDR.Classlabel,ng],(1:16:1152)',[1:72]'==TI,'REG');
173CC94.TSD.T = CC94.TSD.T/HDR.SampleRate;
174CC95 = findclassifier(a,HDR.EVENT.POS(HDR.EVENT.TYP==hex2dec('300'))-1,[HDR.Classlabel,ng],(1:16:1152)',[1:72]'==TI,'SVM:LIB');
175CC95.TSD.T = CC95.TSD.T/HDR.SampleRate;
176
177
178if exist('OCTAVE_VERSION','builtin');
179        om = 'Octave';
180elseif 1,
181        om = 'Matlab';
182else
183        om = 'FreeMat';
184end;
185outfile = sprintf('bench_biosig1.75+_%s_%s_%s.mrk',computer,om,version);
186
187try
188        unix(['cat /proc/cpuinfo >"',outfile,'"'])
189catch
190end;
191fid = fopen(outfile,'a');
192fprintf(fid,'\n\nDate:\t%s\n',date);
193fprintf(fid,'Revision:\t$Id: bench_biosig.m,v 1.7 2009/04/21 10:39:20 schloegl Exp $\n');
194fprintf(fid,'Computer:\t%s\nSoftware:\t%s\nVersion:\t%s\n',computer,om,version);
195
196tmp = [diff([0,t(:)']);t(:)']';
197for k=1:K,
198        fprintf(fid,'%25s:\t%6.1f\t%6.1f\n',jo{k},tmp(k,:));
199end;
200fclose(fid);
201
202