www.gusucode.com > mbctools 工具箱 matlab 源码程序 > mbctools/getClusters.m

    function [dataInds, designInds] = getClusters(design, data, tol)
% [DATAINDS, DESIGNINDS] = GETCLUSTERS(DESIGN, DATA, TOL)

% CLUSTERS is a cellarray of clusters {[2x1 cell], [2x1 cell], .... }
% each [2x1 cell] is {[designInds]; [dataInds]}
% ADJACENCYMATRIX has numRows = numDesignPoints, numCols = numDataPoints
% ones label where adjacency occurs
% 

% [CLUSTERS, ADJACENCYMATRIX, MDESI, MDATI, UNMDESI, UNMDATI] = GETCLUSTERS(DESIGN, DATA, TOL)
% MDESI matched design indices
% MDATI matched data indices
% UNMDESI unmatched design indices
% UNMDATI unmatched data indices

%  Copyright 2000-2004 The MathWorks, Inc. and Ford Global Technologies, Inc.


designByData = [];
clusters = {};

%% ==== FIND THE EDGES FOR THE ADJACENCY MATRIX =====

%% data has records as rows, columns are factors
[numDesign, numDesignFactors] = size(design);
[numData, numDataFactors] = size(data);

%% check that tolerance has an entry for each dim and sizes are consistent
if ~(isnumeric(tol) & numDesignFactors==numDataFactors & length(tol)==numDesignFactors)
   return
end
% initialize output to be a logical array
designByData = false(numDesign, numData);
tolMtx = repmat(tol(:)', numData, 1);

%% check for "closeness" to each design point in turn
for thisDesPoint = 1:numDesign
   reppdThisDesPoint = repmat(design(thisDesPoint, :),[numData,1]);
   distToData = abs( data - reppdThisDesPoint );
   designByData(thisDesPoint,:) = all(distToData < tolMtx, 2)';
end

%% ==== FIND THE CLUSTERS =====

[numDesign, numData] = size(designByData);
%clusters = {};
dataInds = {};
designInds = {};
% The only rows worth looking at
LRowInd = any(designByData,2);

while any(LRowInd)
   %% set LDesignInds = [0 0 0 0 0]
   LDesignInds = false(numDesign,1);
   LDataInds = false(numData,1);
   %% label the starting row
   LDesignInds(min(find(LRowInd))) = true;
   
   while 1
      %% look for data adjacent to points in LDesignInds
      moreDataInds = any(designByData(LDesignInds, :),1);
      %% look for design points adjacent to points in moreDataInds
      moreDesInds = any(designByData(:, moreDataInds),2);
      
      if isequal(moreDesInds,LDesignInds) & isequal(moreDataInds,LDataInds)
         break
      else
         LDesignInds = moreDesInds;
         LDataInds = moreDataInds;
      end
      
   end
   %% remove from LRowInd the design points in this cluster
   LRowInd = LRowInd & ~LDesignInds;
   %% save the designPoints and dataPoints in this cluster (make sure
   %% designInds are horizontal)
   dataInds{end+1} = find(LDataInds);
   designInds{end+1} = find(LDesignInds)';
%   clusters{end+1} = {find(LDesignInds); find(LDataInds)};
%   out{end+1} = {LDesignInds; LDataInds};
end


% ----- find stats -----
numClusters = length(dataInds);
[numDesign, numData] = size(designByData);
%% ----- find lone gunmen -----
unmatchedDesignInds = 1:numDesign;
unmatchedDataInds = 1:numData;
unmatchedDesignInds = setdiff(unmatchedDesignInds, [designInds{:}]);
unmatchedDataInds = setdiff(unmatchedDataInds, [dataInds{:}]);
%% ----- find tribal groups -----
matchedDesignInds = setdiff([1:numDesign], unmatchedDesignInds);
matchedDataInds = setdiff([1:numData], unmatchedDataInds);


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%  ALTERNATIVE IMPLEMENTATION not checking twice any row/column
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% while any(LRowInd)
%    %% set LDesignInds = [0 0 0 0 0]
%    LTotalDesignInds = false(numDesign,1);
%    LTotalDataInds = false(numData,1);
%    %% label the starting row
%    LTotalDesignInds(min(find(LRowInd))) = true;
%    
%    %% initialize LDesignInds with a design point to start building the
%    %% cluster
%    LDesignInds = LTotalDesignInds;
%    while 1
%       %% look for data adjacent to points in LDesignInds
%       LDataInds = any(designByData(LDesignInds, :),1)';
%       %% look for design points adjacent to points in moreDataInds
%       LDesignInds = any(designByData(:, LDataInds),2);
%       
%       if ~any(LDataInds) & ~any(LDesignInds)
%          break
%       else
%          %% design ind to check next iteration are those not already
%          %% checked
%          LnextDesignInds = LTotalDesignInds & ~LDesignInds;
%          %% augment to  running total of design points
%          LTotalDesignInds = LTotalDesignInds | LDesignInds;
%          %% augment the running total of data inds
%          LDataInds = LTotalDataInds | LDataInds;
%          %% give correct name to design inds to be checked in the next
%          %% iteration
%          LDesignInds = LnextDesignInds;
%       end
%       
%    end
%    %% remove from LRowInd the design points in this cluster
%    LRowInd = LRowInd & ~LTotalDesignInds;
%    %% save the designPoints and dataPoints in this cluster
%    out{end+1} = {find(LTotalDesignInds); find(LTotalDataInds)};
% %   out{end+1} = {LDesignInds; LDataInds};
% end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  ALTERNATIVE IMPLEMENTATION USING INDICES RATHER THAN LOGICAL ARRAYS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% allDesignInds = find(LRowInd);
% out = {};
% 
% while ~isempty(allDesignInds)
%    designInds = allDesignInds(1);
%    dataInds = [];
%    
%    while 1
%       %% find dataPoints that are connected to these designPoints
%   %    [null, moreDataInds] = find(designByData(designInds, :));
%       moreDataInds = find(any(designByData(designInds, :),1));
%       moreDataInds = unique(moreDataInds);
%       %% find designPoints that are connected to these dataPoints
% %      [moreDesInds, null] = find(designByData(:, moreDataInds));
%       moreDesInds = find(any(designByData(:, moreDataInds),2));
%       moreDesInds = unique(moreDesInds);
%       
%       if isequal(moreDesInds,designInds) & isequal(moreDataInds,dataInds)
%          break
%       else
%          designInds = moreDesInds;
%          dataInds = moreDataInds;
%       end
%       
%    end
%    
%    %% kick out the designPoints in this cluster from those we are searching
%    %% over
%    allDesignInds = setdiff(allDesignInds, designInds);
%    %% save the designPoints and dataPoints in this cluster
%    out{end+1} = {designInds; dataInds};
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%