Daniel@0: function [a d] = mircluster(a,varargin)
Daniel@0: %   c = mircluster(a,f) clusters the segments in the audio sequence(s) 
Daniel@0: %       contained in the audio object a, along the analytic feature(s) 
Daniel@0: %       f, using the k-means strategy. Multiple analytic features have to
Daniel@0: %       be grouped into one array of cells.
Daniel@0: %       Example:
Daniel@0: %           sg = mirsegment(a);
Daniel@0: %           mircluster(sg, mirmfcc(sg))
Daniel@0: %           mircluster(sg, {mirmfcc(sg), mircentroid(sg)})
Daniel@0: %   c = mircluster(d) clusters the frame-decomposed data d into groups
Daniel@0: %       using K-means clustering.
Daniel@0: %       Example:
Daniel@0: %           cc = mirmfcc(a,'Frame');
Daniel@0: %           mircluster(cc)
Daniel@0: %   Optional argument:
Daniel@0: %       mircluster(...,n) indicates the maximal number of clusters.
Daniel@0: %           Default value: n = 2.
Daniel@0: %       mircluster(...,'Runs',r) indicates the maximal number of runs.
Daniel@0: %           Default value: r = 5.
Daniel@0: %
Daniel@0: %   Requires SOM Toolbox (included in the MIRtoolbox distribution).
Daniel@0: 
Daniel@0: 
Daniel@0:         nruns.key = 'Runs';
Daniel@0:         nruns.type = 'Integer';
Daniel@0:         nruns.default = 5;
Daniel@0:     option.nruns = nruns;
Daniel@0:      
Daniel@0:         nclust.position = 2;
Daniel@0:         nclust.type = 'Integer';
Daniel@0:         nclust.default = 2;
Daniel@0:     option.nclust = nclust;
Daniel@0: 
Daniel@0: specif.option = option;
Daniel@0: 
Daniel@0: specif.nochunk = 1;
Daniel@0:     
Daniel@0: d = a;
Daniel@0: if isa(a,'mirdesign')
Daniel@0:     if not(get(a,'Eval'))
Daniel@0:         % During bottom-up construction of the general design
Daniel@0: 
Daniel@0:         [unused option] = miroptions(@mircluster,a,specif,varargin);
Daniel@0:         type = get(a,'Type');
Daniel@0:         a = mirdesign(@mircluster,a,option,{},struct,type);  
Daniel@0:         a = set(a,'NoChunk',1);
Daniel@0:     else
Daniel@0:         % During top-down evaluation initiation
Daniel@0:         
Daniel@0:         e = evaleach(a);
Daniel@0:         if iscell(e)
Daniel@0:             e = e{1};
Daniel@0:         end
Daniel@0:         a = mircluster(e,varargin{:});
Daniel@0:     end
Daniel@0: else
Daniel@0:     if not(isa(a,'mirdata'))
Daniel@0:         mirerror('mircluster','The input should be either frame- or segment-decomposed.');
Daniel@0:     end
Daniel@0:    
Daniel@0:     if isempty(varargin) || (not(isa(varargin{1},'mirdata') || ...
Daniel@0:                              (iscell(varargin{1}) && ...
Daniel@0:                               isa(varargin{1}{1},'mirdata'))))
Daniel@0:         % mircluster version for frame-decomposed data:
Daniel@0:         % frames are clustered into groups using K-means clustering.
Daniel@0:         [unused option] = miroptions(@mircluster,a,specif,varargin);
Daniel@0:         da = get(a,'Data');
Daniel@0:         lva = length(da); % Number of audio files in the audio object.
Daniel@0:         c = cell(1,lva);
Daniel@0:         display('Clustering frames...');
Daniel@0:         if mirwaitbar
Daniel@0:             handle = waitbar(0,'Clustering frames...');
Daniel@0:         else
Daniel@0:             handle = 0;
Daniel@0:         end
Daniel@0:         for j = 1:lva           % For each audio file,...
Daniel@0:             va = [];            % Data transmitted to the kmeans_cluster function.
Daniel@0:             v = da{j};
Daniel@0:             if iscell(v)
Daniel@0:                 v = uncell(v,-Inf); %v{1};
Daniel@0:             end
Daniel@0:             if size(v,4)>1
Daniel@0:                 v(end+1:2*end,:,:,1) = v(:,:,:,2);
Daniel@0:                 v(:,:,:,2) = [];
Daniel@0:             end
Daniel@0:             % Standardization
Daniel@0:             %stv = std(v,0,2);
Daniel@0:             %stv(find(stv == 0)) = 1;
Daniel@0:             va(end+1:end+size(v,1),:,:) = v;%...
Daniel@0:                 %(v - repmat(mean(v,2),[1 size(v,2) ])) ...
Daniel@0:                 %./ repmat(stv,[1 size(v,2) ]);
Daniel@0:             if isa(a,'mirscalar')
Daniel@0:                 m = get(a,'Mode');
Daniel@0:                 if not(isempty(m))
Daniel@0:                     m = m{j};
Daniel@0:                     val = [];
Daniel@0:                     for l = 1:nseg
Daniel@0:                         vl = m{l};
Daniel@0:                         if iscell(vl)
Daniel@0:                             vl = vl{1};
Daniel@0:                         end
Daniel@0:                         val(:,l) = vl;
Daniel@0:                     end
Daniel@0:                     stv = std(val,0,2);
Daniel@0:                     stv(find(stv == 0)) = 1;
Daniel@0:                     va(end+1:end+size(val,1),:) = ...
Daniel@0:                         (val - repmat(mean(val,2),[1 size(val,2) ])) ...
Daniel@0:                             ./ repmat(stv,[1 size(val,2) ]);
Daniel@0:                 end
Daniel@0:             end
Daniel@0:             if size(va,3)>1
Daniel@0:                 mel = 1;
Daniel@0:                 va = reshape(va,size(va,2),size(va,3))';
Daniel@0:             else
Daniel@0:                 mel = 0;
Daniel@0:             end
Daniel@0:             [cc, p, err, ind] = kmeans_clusters(va',option.nclust,option.nruns);
Daniel@0:             [minind select] = min(ind);
Daniel@0:             c{j}.centr = cc{select}';
Daniel@0:             c{j}.index = p{select};
Daniel@0:             c{j}.weight = zeros(1,size(cc{select},1));
Daniel@0:             c{j}.covar = zeros(size(cc{select}'));
Daniel@0:             ii = 1;
Daniel@0:             for i = 1:size(c{j}.centr,2)
Daniel@0:                 clus = va(:,c{j}.index == ii);
Daniel@0:                 if isempty(clus)
Daniel@0:                     higher = find(c{j}.index > ii);
Daniel@0:                     c{j}.index(higher) = c{j}.index(higher)-1;
Daniel@0:                     c{j}.centr(:,ii) = [];
Daniel@0:                     c{j}.weight(ii) = [];
Daniel@0:                     c{j}.covar(:,ii) = [];
Daniel@0:                 else
Daniel@0:                     c{j}.weight(ii) = size(clus,2)/size(va,2);
Daniel@0:                     if c{j}.weight(ii) == 0
Daniel@0:                         pause
Daniel@0:                     end
Daniel@0:                     c{j}.covar(:,ii) = mean((clus'-ones(1,size(clus,1))*c{j}.centr(:,ii)).^2);
Daniel@0:                     ii = ii+1;
Daniel@0:                 end
Daniel@0:             end
Daniel@0:             if handle
Daniel@0:                 waitbar(j/lva,handle);
Daniel@0:             end
Daniel@0:         end
Daniel@0:         if handle
Daniel@0:            delete(handle)
Daniel@0:         end
Daniel@0:         a = set(a,'Clusters',c);
Daniel@0:     else
Daniel@0:         % mircluster version for segmented audio:
Daniel@0:         % segments are clustered into groups using K-means clustering.
Daniel@0:         da = varargin{1};
Daniel@0:         varargin(1) = [];
Daniel@0:         [unused option] = miroptions(@mircluster,a,specif,varargin);
Daniel@0:         display('Clustering segments...');
Daniel@0:         if isa(da,'mirdata') || (iscell(da) && isa(da{1},'mirdata'))
Daniel@0:             if not(iscell(da))
Daniel@0:                 da = {da};
Daniel@0:             end
Daniel@0:             vala = get(a,'Data');    % Data contained in the audio object a.
Daniel@0:             lva = length(vala); % Number of audio files in the audio object.
Daniel@0:             clus = cell(1,lva);
Daniel@0:             for j = 1:lva           % For each audio file,...
Daniel@0:                 va = [];            % Data transmitted to the kmeans_cluster function.
Daniel@0:                 nseg = length(vala{j}); % Number of segments in the audio file.
Daniel@0:                 for i = 1:length(da)    % For each analytic feature,...
Daniel@0:                     v = get(da{i},'Data');
Daniel@0:                     v = v{j};
Daniel@0:                     if iscell(v)
Daniel@0:                         v = uncell(v,-Inf); %v{1};
Daniel@0:                     end
Daniel@0:                     val = [];           
Daniel@0:                     if size(v,4)>1
Daniel@0:                         v(end+1:2*end,:,:,1) = v(:,:,:,2);
Daniel@0:                         v(:,:,:,2) = [];
Daniel@0:                     end
Daniel@0: 
Daniel@0:                     % Standardization
Daniel@0:                     stv = std(v,0,2);
Daniel@0:                     stv(find(stv == 0)) = 1;
Daniel@0:                     va(end+1:end+size(v,1),:) = ...
Daniel@0:                         (v - repmat(mean(v,2),[1 size(v,2) ])) ...
Daniel@0:                         ./ repmat(stv,[1 size(v,2) ]);
Daniel@0:                     if isa(da{i},'mirscalar')
Daniel@0:                         m = get(da{i},'Mode');
Daniel@0:                         if not(isempty(m))
Daniel@0:                             m = m{j};
Daniel@0:                             val = [];
Daniel@0:                             for l = 1:nseg
Daniel@0:                                 vl = m{l};
Daniel@0:                                 if iscell(vl)
Daniel@0:                                     vl = vl{1};
Daniel@0:                                 end
Daniel@0:                                 val(:,l) = vl;
Daniel@0:                             end
Daniel@0:                             stv = std(val,0,2);
Daniel@0:                             stv(find(stv == 0)) = 1;
Daniel@0:                             va(end+1:end+size(val,1),:) = ...
Daniel@0:                                 (val - repmat(mean(val,2),[1 size(val,2) ])) ...
Daniel@0:                                     ./ repmat(stv,[1 size(val,2) ]);
Daniel@0:                         end
Daniel@0:                     end
Daniel@0: 
Daniel@0:                 end
Daniel@0:                 [cc, p, err, ind] = kmeans_clusters(va',min(option.nclust,nseg),option.nruns);
Daniel@0:                 clus{j} = p{end};
Daniel@0:             end
Daniel@0:             a = set(a,'Clusters',clus);
Daniel@0:             t = get(a,'Time'); 
Daniel@0:             fp = get(a,'FramePos'); 
Daniel@0:             for j = 1:lva           % For each audio file,...
Daniel@0:                 aj = vala{j};
Daniel@0:                 tj = t{j};
Daniel@0:                 fpj = fp{j};
Daniel@0:                 clj = clus{j};
Daniel@0:                 k = 2;
Daniel@0:                 while k <= length(aj)
Daniel@0:                     if clj(k) == clj(k-1)
Daniel@0:                         aj{k-1} = [aj{k-1};aj{k}];
Daniel@0:                         aj(k) = [];
Daniel@0:                         tj{k-1} = [tj{k-1};tj{k}];
Daniel@0:                         tj(k) = [];
Daniel@0:                         fpj{k-1} = [fpj{k-1}(1);fpj{k}(2)];
Daniel@0:                         fpj(k) = [];
Daniel@0:                         clj(k) = [];
Daniel@0:                         k = k-1;
Daniel@0:                     end
Daniel@0:                     k = k+1;
Daniel@0:                 end
Daniel@0:                 vala{j} = aj;
Daniel@0:                 t{j} = tj;
Daniel@0:                 fp{j} = fpj;
Daniel@0:                 cl{j} = clj;
Daniel@0:             end
Daniel@0:             a = set(a,'Data',vala,'Time',t,'FramePos',fp,'Clusters',cl);
Daniel@0:         end
Daniel@0:     end
Daniel@0: end