1function [CC,KAPPA,tsd]=findclassifier(D,TRIG,cl,MODE1,t0,MODE)
2% FINDCLASSIFIER
3%   identifies and validates a classifier of a BCI systems [1-3].
4%   Several evaluation criteria are obtained [4]. Several Cross-validation
5%   procedures are supported.
6%
7% By default, a Trial-based Leave-One-Out-Method is used for Crossvalidation
8%	MODE.Segments = class_times;
9%	MODE.WIN   = t_ref;
10%       [CC,Q,TSD] = findclassifier(D,TRIG,Class,MODE,t_ref,TYPE);
11% Also this will work but its use is discouraged (it might become obsolete).
12%       [CC,Q,TSD] = findclassifier(D,TRIG,Class,class_times,t_ref,TYPE);
13% An K-fold cross-validation can be applied in this way:
14%       ng = floor([0:length(Class)-1]'/length(Class)*K);
15%       [CC,Q,TSD] = findclassifier(D,TRIG,[Class,ng],...);
16%
17% D 	data, each row is one time point
18% TRIG	trigger time points
19% Class class information
20% class_times	classification times, combinations of times must be in one row
21% t_ref	reference time for Class 0 (optional)
22% TYPE  determines the type of classifier (see HELP TEST_SC for complete list)
23%       bthe default method is 'LD3'
24%
25% CC 	contains the classifier and the validation results
26% Q  	is a list of classification quality for each time segment (as defined by 'class_times')
27% TSD 	returns the discrimination
28%
29% Example:
30%  [CC,Q,TSD]=findclassifier(d,find(trig>0.5)-257,~mod(1:80,2),reshape(1:14*128,16,14*8)');
31%
32% see also: TRAIN_SC, TEST_SC, BCI4EVAL
33%
34% Reference(s):
35% [1] Schlögl A, Neuper C, Pfurtscheller G
36% 	Estimating the mutual information of an EEG-based Brain-Computer-Interface
37%  	Biomedizinische Technik 47(1-2): 3-8, 2002.
38% [2] Schlögl A, Keinrath C, Scherer R, Pfurtscheller G,
39%	Information transfer of an EEG-based Bran-computer interface.
40%	Proceedings of the 1st International IEEE EMBS Conference on Neural Engineering, Capri, Italy, Mar 20-22, 2003
41% [3] Schlögl A, Lee FY, Bischof H, Pfurtscheller G
42%	Characterization of Four-Class Motor Imagery EEG Data for the BCI-Competition 2005.
43%	Journal of neural engineering 2 (2005) 4, S. L14-L22
44% [4] Schlögl A, Kronegg J, Huggins JE, Mason SG;
45%	Evaluation criteria in BCI research.
46%	(Eds.) G. Dornhege, J.R. Millan, T. Hinterberger, D.J. McFarland, K.-R.Müller;
47%	Towards Brain-Computer Interfacing, MIT Press, p327-342, 2007
48
49%   $Id$
50%   Copyright (C) 1999-2006 by Alois Schloegl <alois.schloegl@gmail.com>
51%   This is part of the BIOSIG-toolbox http://biosig.sf.net/
52
53
54% BioSig is free software; you can redistribute it and/or
55% modify it under the terms of the GNU General Public License
56% as published by the Free Software Foundation; either version 3
57% of the  License, or (at your option) any later version.
58%
59% This program is distributed in the hope that it will be useful,
60% but WITHOUT ANY WARRANTY; without even the implied warranty of
61% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
62% GNU General Public License for more details.
63%
64% You should have received a copy of the GNU General Public License
65% along with this program; if not, write to the Free Software
66% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
67
68if ~isempty(t0)
69	warning('arg5 (t_ref) should be empty. Use MODE.WIN=t_ref instead.')
70end;
71
72CC = []; Q = [];tsd=[];md=[];
73if isstruct(MODE1)
74	CC.T = MODE1;
75	T  = MODE1.Segments;
76	t0 = MODE1.WIN;
77else
78	T = MODE1;
79end;
80
81if nargin<6,
82        MODE.TYPE='LD3';
83elseif ischar(MODE);
84        tmp = MODE;
85        clear MODE
86        MODE.TYPE = tmp;
87end;
88if nargin>4,
89        if isempty(t0),
90                t0=logical(ones(size(T,1),1));
91        end;
92end;
93tmp=cl;tmp(isnan(tmp))=0;
94if any(rem(tmp,1) & ~isnan(cl)),
95        fprintf(2,'Error %s: class information is not integer\n',mfilename);
96        return;
97end;
98if length(TRIG)~=length(cl);
99        fprintf(2,'Warning FINDCLASSIFIER: number of Triggers do not match class information');
100end;
101
102if size(cl,1)~=length(TRIG);
103        fprintf(2,'Warning FINDCLASSIFIER: Classlabels must be a column vector');
104        if length(TRIG)==size(cl,2),
105                cl = cl';
106        end;
107end;
108
109TRIG = TRIG(:);
110tmp  = ~any(isnan(cl),2);
111TRIG = TRIG(tmp);
112cl   = cl(tmp);
113
114if size(cl,2)>1,
115        cl2 = cl(:,2);          % 2nd column contains the group definition, ( Leave-One (Group) - Out )
116        cl  = cl(:,1);
117else
118        cl2 = [1:length(cl)]';  % each trial is a group (important for cross-validation); Trial-based LOOM
119end;
120[CL,i,cl] = unique(cl);
121CL2 = unique(cl2);
122%[CL,iCL] = sort(CL);
123
124% add sufficient NaNs at the beginning and the end
125tmp = min(TRIG)+min(min(T))-1;
126if tmp<0,
127        TRIG = TRIG - tmp;
128        D = [repmat(nan,[-tmp,size(D,2)]);D];
129end;
130tmp = max(TRIG)+max(max(T))-size(D,1);
131if tmp>0,
132        D = [D;repmat(nan,[tmp,size(D,2)])];
133end;
134
135% estimate classification result for all time segments in T - without crossvalidation
136CMX = zeros([size(T,1),length(CL)*[1,1]]);
137KAPPA = repmat(NaN,size(t0));
138for k = 1:size(T,1),
139        if t0(k),
140                c = []; d = [];
141                for k1 = 1:length(CL),
142                        %t = perm(TRIG(cl==CL(k1)),T(k,:));
143                        t = perm(TRIG(cl==k1),T(k,:));
144                        d = [d; D(t(:),:)];
145                        %c = [c; repmat(CL(k1),prod(size(t)),1)];
146                        c = [c; repmat(k1,prod(size(t)),1)];
147                end;
148                cc{k} = train_sc(d,c,MODE);
149                r     = test_sc(cc{k},d,MODE,c);
150                KAPPA(k)  = r.kappa;
151        end;
152%        fprintf(1,'search for segment: %i-%i kappa=%4.2f\n',T(k,1),T(k,end),KAPPA(k));
153end;
154[maxQ,TI] = max(KAPPA.*(t0~=0)); %d{K},
155%[maxQ,TI] = max(q.*(t0~=0)); %d{K},
156CC = cc{TI};
157CC.KAPPA = KAPPA;
158CC.TI = TI;
159CC.TC = T(TI,:);
160if isstruct(MODE1)
161	CC.T = MODE1;
162end;
163
164if isnan(maxQ)
165	fprintf(2,'ERROR FINDCLASSIFIER: no valid classifier available.\n');
166	return;
167end;
168
169%% cross-validation with jackknife (group-based leave-one-out-method)
170nc  = max(max(T))-min(min(T))+1;
171
172M = length(CL);
173tsd = repmat(nan,[nc*length(TRIG),M]);
174tt  = tsd(:,1);
175IX  = find(~isnan(cl(:)'));
176
177for l = 1:length(CL2);        % XV based on "Leave-One(group)-Out-Method"
178        ix = find(cl2==CL2(l));         % identify members of l-th group
179        t  = perm(TRIG(ix), T(CC.TI,:));        % get samples of test set
180
181%        fprintf(1,'\nX-V (%i/%i):',l, length(CL2));
182        % decremental learning
183        if 0, ~isempty(strfind(CC.datatype,'statistical')),
184                c  = repmat(cl(cl2==CL2(l))', size(T,2),1);     % classlabels of test set
185                cc = untrain_sc(CC,c(:),D(t(:),:));             % untraining test set
186
187        elseif 1, %~isempty(strfind(CC.datatype,'statistical')),
188                t  = perm(TRIG(cl2~=CL2(l)), T(CC.TI,:));       % samples of training set
189                c  = repmat(cl(cl2~=CL2(l))', size(T,2),1);     % classlabels of training set
190                cc = train_sc(D(t(:),:),c(:),MODE);             % train classifier
191
192        end;
193
194        t  = perm(TRIG(cl2==CL2(l)), min(min(T)):max(max(T)));  % samples of evaluation set (l-th group)
195        ix = perm((find(cl2==CL2(l))-1)*nc, 1:nc);              % save to ...
196        if any(~isnan(tt(ix(:))))
197                fprintf(2,'WARNING FINDCLASSIFIER#: overlapping segments %i\n',sum(~isnan(tt(ix(:)))));
198        end;
199        tt(ix(:)) = t;
200
201        if ~isempty(strfind(CC.datatype,'svm:lib:1vs1')) | ~isempty(strfind(CC.datatype,'svm:lib:rbf')),
202                c = repmat(cl(cl2==CL2(l))', size(t,1), 1);     % classlabels of test set
203                r = test_sc(cc,D(t(:),:),MODE,c(:));            % evaluation of l-th group
204        else
205                r = test_sc(cc,D(t(:),:),MODE);                 % evaluation of l-th group
206        end;
207        tsd(ix(:),1:M) = r.output(:,1:M);                                % save results of l-th group
208end;
209
210%CC.TSD  = bci4eval(tsd, (0:length(cl)-1)'*nc, cl, 1, nc);
211CC.TSD  = bci4eval(tsd, (0:length(cl)-1)'*nc, CL(cl), 1, nc);
212CC.TSD.T = CC.TSD.T - 1 + min(T(:));
213