AP聚类Matlab程序

function [idx a r output] = affprop(s, options)
m = size(s,1);

defaultopt = struct( 'MaxIter', 200,...
                     'StallIter', 20, ...
                     'FunValCheck','off', ...
                     'Dampening', .5, ...
                     'OutputFcn',[]);

if nargin < 2
    options = [];
end;

maxiter   = optimget(options,'MaxIter',   defaultopt,'fast');
stalliter = optimget(options,'StallIter', defaultopt,'fast');
lambda    = optimget(options,'Dampening', defaultopt,'fast');
outputfcn = optimget(options,'OutputFcn', defaultopt,'fast');

if isempty(outputfcn)
    haveoutputfcn = false;
else
    haveoutputfcn = true;
end


% useful indices
ri = (1:m)';
di = 1:(m+1):m*m;         %index to diagonal elements

% initialize responsibility matrix
r = zeros(m);

% availability matrix,a, evidence from other datapoints how an exemplar
% point i will make
a      = zeros(m);
a_prev = zeros(m);

%initialize exemplar vector
idx_prev = ri;

iter = 1;
converged = false;
ccount = 0;

while ~converged
    %% Part I. Update or initialize r, the responsibility matrix
    % update (or initialize) r, the responsibility matrix, such that element
    % r(i,j) indicates how good an exemple point x(j,:) would make for
    % point x(i,:)

    %  given as = a + s,
    %  set d(i,j) = max( as(i,k) ), for k not equal to j
    as = a + s;

    %rmax is the row maximum for all elements
    [rmax k] = max(as,[],2);
    kp = sub2ind( [m m], ri, k );

    %rmax2 is the row maximum for when j=k would be the row maximum
    as(kp) = nan;
    rmax2  = max(as,[],2);

    d = repmat( rmax, 1, m );
    d(kp) = rmax2;

    %     r = s - d;                                                 %eq(1)
    r = lambda*r + (1-lambda)*(s-d);       % dampended


    %% Part II. Update the availability matrix [ eq2 and eq3 ]
    %  The availability matrix is evidence that point k is an exemplar based on
    %  positive responsibiliities sent to it from the other points
    % a(i,j) = min[ 0, r(j,j) + sum( max(0, r(k,j) ) ],
    % for all k not in (i,j)

    % calculate c = max(r,0)
    c = r;
    c(c<0) = 0;     % only consider positive elements

    % calculate the sum( max(0, r(k,j) ), where k not in (i,j) this is done
    % by setting c(j,j) to zero. Then sum the columns and form a new
    % matrix,cs, by replication. Finally, subtracting c(i,j) gives
    % the desired matrix, c (reuse storage)
    c(di) = 0;                   % set diagonal elements to zero
    cs = repmat(sum(c),m,1);     % sum columns
    c = cs - c;

    %     a = min(0, repmat( diag(r)', m, 1) + c );                   %eq(2)
    %     a(di) = cs(di);                                             %eq(3)
    %
    a = lambda*a_prev + (1-lambda)*(min(0, repmat( diag(r)', m, 1) + c ));
    a(di) = lambda*a_prev(di) + (1-lambda)*cs(di);

    a_prev = a;

    if haveoutputfcn
        outputfcn(a,r);
    end

    iter = iter+1;
    if iter > maxiter
        converged = true;
    end

    [rmax idx] = max( a+r,[],2);
    if isequal(idx,idx_prev)
        ccount = ccount + 1;
        if ccount == stalliter
            converged = true;
        end
    else
        idx_prev = idx;
        ccount = 0;
    end




end

output.iter      = iter;
output.stalliter = ccount;

M = [5 0; 0 5; 10 10];
idx = repmat( 1:3, 50, 1 );
x = M(idx(:),:) + randn(150,2);

% generate similarity matrix
m     = size(x,1);
s     = zeros(m);

o     = nchoosek(1:150,2);      % set up all possible pairwise comparisons 
xx    = x(o(:,1),:)';           % point 1
xy    = x(o(:,2),:)';           % point 2
d2    = (xx - xy).^2;           % distance squared
d     = -sqrt(sum(d2));          % distance

k     = sub2ind([m m], o(:,1), o(:,2) );    % prepare to make square 
s(k)  = d;             
s     = s + s';
di = 1:(m+1):m*m;         %index to diagonal elements

s(di) = min(d);

%% clustering
options.StallIter = 10;
options.OutputFcn = @(a,r) affprop_plot(a,r,x,'k.');

figure
ex       = affprop(s, options );

function affprop_plot( a, r, x, cc )
%AFFPROP_PLOT throw away function to illustrate affinity propagation
% 
%Example
%   for example see
%See also affprop_demo, affprop

if nargin < 4
    cc = '.';
end

m = size(x,1);
subplot(2,2,1);
cla
plot(x(:,1),x(:,2), cc, 'markersize', 3 );

[rmax k] = max( a+r,[],2);


% color examplars
i = (1:m)';
hold on;
plot(x(i==k,1),x(i==k,2), 'r.', 'markersize', 7 );

j = i~=k;

% quiver( x(j,1),              x(j,2), ...
%         x(k(j),1) - x(j,1),  x(k(j),2) - x(j,2), ...
%         0 );

quiver( x(j,1),              x(j,2), ...
    x(k(j),1) - x(j,1),  x(k(j),2) - x(j,2), ...
    0, '.' );

subplot(2,2,2)
colormap(copper)
imagesc(a);
title('availability')

subplot(2,2,3);
colormap(copper)
imagesc(r);
title('responsibility')

subplot(2,2,4);
colormap(copper)
imagesc(a+r);
title('availability+responsibility');
pause(0.1);
drawnow