FSK_src.zip_feature selection

% ========================================================================= % ================ FEATURE SELECTION FUNCTIONS ================== % ========================================================================= function res = featSelect (X,labels,method,dimSel,kertype,maxmem_kb) res = []; if nargin < 4 kertype = 'rbf'; end if nargin < 6 maxmem_kb = 100*1024; % 100 MBytes end if ~strcmp(kertype,'rbf') && ~strcmp(kertype,'linear') fprintf ('Error : kertype must be either ''rbf'' or ''linear''.\n'); return end switch upper(method) case 'FAS' % Forward Alignement Selection func = @featSelectForward; args = {}; case 'SFS' % Scaled Frobenius Selection func = @featSelectScaled; args = {0}; case 'SAS' % Scaled Alignement Selection func = @featSelectScaled; args = {1}; case 'SCSS' % Scaled Class Separability Selection func = @featSelectScaled; args = {3}; case 'KFDS' % Kernel Fisher Discriminant Selection func = @featSelectKFD; args = {}; otherwise fprintf ('Error : Unknown method %s.\n',method); res = []; return; end res.method = method; res.dimData = size(X,2); if nargin < 3 dimSel = res.dimData; end res.dimSel = dimSel; X = featNormalize(X); fprintf('* Feature selection method %s with kernel %s (max memory occupation %dkB).\n',method,kertype,round(maxmem_kb)); t0cpu=cputime; tic; [res.iSelected,res.r] = feval(func,X,labels,dimSel,kertype,maxmem_kb,args{:}); res.tcpu = cputime-t0cpu; res.telaps = toc; fprintf('Elapsed time: %f, CPU time: %f\n',res.telaps,res.tcpu); % ------------------------------------------------------------------------- function X = featNormalize(X) nbFr = size(X,1); Xmean = mean(X,1); Xvar = std(X,0,1); Xvar(Xvar==0) = 1; X = X - repmat(Xmean,nbFr,1); X = X ./ repmat(Xvar ,nbFr,1); % ------------------------------------------------------------------------- function [iSelected r] = featSelectScaled (X,labels,dimSel,kertype,maxmem_kb,critID) if strcmp(kertype,'linear') && critID == 0 % SFS with linear kernel fprintf ('WARNING : linear kernel is not advised with the SFS method. Frobenius criterion will diverge to infinity.\n'); end dim = size(X,2); ker = kerInit(kertype); ker = kerGradCrit(ker,X,labels,critID,maxmem_kb); [r iSelected] = sort(ker.w,'descend'); if dimSel < dim r = r(1:dimSel); iSelected = iSelected(1:dimSel); end ker.w = ker.w(iSelected); ker.dim = dimSel; % ------------------------------------------------------------------------- function [iSelected r] = featSelectForward(X,labels,dimSel,kertype,maxmem_kb) dim = size(X,2); % input data number of features info.dim = 0; % number of features selected so far info.testdim = []; info.bufdim = []; critID = 1; % criterion = KTA % selection sets iRemain = true(1,dim); % remaining not-selected features iSelected = zeros(1,dimSel); % selected features list r = zeros(1,dimSel); % selected features ranks % kernel initialisation ker = kerInit(kertype); ker.dim = dim; gamma = kerGamma(ker); ker.dim = 0; % blocs initialisation blocs = dataBlocs(X,labels,maxmem_kb); [blocs.preK] = deal(0); nBlocs = length(blocs); % ##### MAIN LOOP ##### for si=1:dimSel fprintf('Feature %3d/%3d: ',si,dimSel); dimevallist = find(iRemain); nbdimeval = length(dimevallist); ker.dim = ker.dim + 1; crit = zeros(1,nbdimeval); for idim=1:nbdimeval dimeval = dimevallist(idim); res = kerFrobResInit(blocs(1).nbFr); for bi=1:nBlocs tmpblocs = blocsComp(blocs,ker,2,dimeval,bi); K = blocs(bi).preK + tmpblocs.dimK{dimeval}; if strcmp(kertype,'rbf') K = exp(-gamma*K); end res = kerFrobAdd(res,blocs(bi),K); end crit(idim) = kerFrobCrit(res,critID); end % KTA maximizing choice [r(si) imax] = max(crit); idim = dimevallist(imax); % blocs update info.testdim = idim; for bi=1:nBlocs tmpblocs = blocsComp(blocs,ker,2,dimeval,bi); blocs(bi).preK = blocs(bi).preK + tmpblocs.dimK{dimeval}; end % update ... iSelected(si) = idim; iRemain(idim) = false; fprintf(' KTA=%1.4f (with dim %d)\n',r(si),idim); end % ------------------------------------------------------------------------- function [iSelected r] = featSelectKFD(X,labels,dimSel,kertype,maxmem_kb) % the kertype and maxmem_kb arguments are only here for % compatibility with the other featSelect functions. if ~strcmp(kertype,'linear') fprintf ('Warning : kernel forced to linear with KFDS method.\n'); kertype = 'linear'; end maxiter = 30; e = 1e-2; step = 0.8; dim = size(X,2); ker = kerInit(kertype); ker.w = ones(1,dim); for it=1:maxiter res = kerFisher(ker,X,labels); wh = abs(res.w); prevw = ker.w; ker.w = ker.w .* wh.^step; cvg = norm(ker.w-prevw)/norm(ker.w); fprintf ('Iteration %d (cvg = %.2f)\n',it,cvg); if abs(cvg) < e, break; end end [r iSelected] = sort(ker.w,'descend'); if dimSel < dim r = r(1:dimSel); iSelected = iSelected(1:dimSel); end % ========================================================================= % ==================== DATA BLOCS FUNCTIONS ===================== % ========================================================================= % ------------------------------------------------------------------------- function blocs = dataBlocs(X,labels,maxmem_kb) if nargin < 3, maxmem_kb = 100*1024; end % already done in dataRetrieve, but done here again for reliability [foo foo2 labels] = unique(labels); [labels isort] = sort(labels); X = X(isort,:); tmp = find(labels ~= labels(1),1)-1; nbFr = [tmp length(labels)-tmp]; dim = size(X,2); offset = [0 nbFr(1)]; maxnb = floor(sqrt(maxmem_kb*1024/8)); % each matlab double value takes 8 Bytes nBlocs = ceil(nbFr/maxnb); offsetnB = [0 nBlocs(1)]; bloc = struct('X',[],'Xi',[],'i',[],'sign',[],'coeff',[],'part',[],'K',[],'preK',[],'dimK',[]); bloc = bloc(1:0); blocs = repmat({bloc},1,5); Xi = 0; for ci=1:2 nbB = nBlocs(ci); tmp = offset(ci)+round(linspace(0,nbFr(ci),nbB+1)); for b=1:nbB Xi = Xi+1; blocs{ci}(b).X = X(tmp(b)+1:tmp(b+1),:); blocs{ci}(b).Xi = [Xi Xi]; end i = [ci ci]; [blocs{ci}.i] = deal(i); [blocs{ci}.sign] = deal(+1); [blocs{ci}.coeff] = deal(1); [blocs{ci}.part] = deal(ci); o = ones(nbB); tmp = o-tril(o); [X1 X2] = ind2sub([nbB nbB],find(tmp)); nb = size(X1,1); if nb==0, continue; end blocsXi = mat2cell(offsetnB(ci)+[X1 X2],ones(1,nb),2); [blocs{ci+2}(1:nb).Xi] = deal(blocsXi{:}); [blocs{ci+2}.i] = deal(i); [blocs{ci+2}.sign] = deal(+1); [blocs{ci+2}.coeff] = deal(2); [blocs{ci+2}.part] = deal(ci); end o = ones(nBlocs); [X1 X2] = ind2sub(nBlocs,find(o)); nb = size(X1,1); try blocsXi = mat2cell([X1 nBlocs(1)+X2],ones(1,nb),2); catch keyboard end [blocs{5}(1:nb).Xi] = deal(blocsXi{:}); [blocs{5}.i] = deal([1 2]); [blocs{5}.sign] = deal(-1); [blocs{5}.coeff] = deal(2); [blocs{5}.part] = deal(3); blocs = [blocs{:}]; blocs(1).nbFr = nbFr; [blocs.dimK] = deal(cell(1,dim)); % ------------------------------------------------------------------------- function blocs = blocsComp (blocs,ker,mode,dimlist,bloclist) % switch type % case 0, returns blocs with field K containing the Gram matrix (default) % case 1, returns blocs with field preK containing the Gram pre-matrix % case 2, returns blocs with field dimK containing the feature-wise Gram pre-matrix if nargin < 3, mode = 0; end nBlocs = length(blocs); if nargin < 5, bloclist = 1:nBlocs; end if mode==2 dim = blocsDim(blocs); if nargin < 4 dimlist = 1:dim; else dimlist = dimlist(:)'; end end for bi=bloclist X1 = blocs(blocs(bi).Xi(1)).X; X2 = blocs(blocs(bi).Xi(2)).X;