% [idx,netsim,dpsim,expref]=apclusterLeveraged(s,p,frac,sweeps)
%
% frac is the fraction of data points to be considered.
% sweeps is the number of times to run affinity propagation.
%

function [idx,netsim,dpsim,expref]=apclusterLeveraged(S,p,frac,sweeps);

% C = zeros(size(S),'uint8'); % for counting S-computations
if isa(S,'function_handle'), N=S(0); else N=size(S,1); end; M=ceil(frac*N);
s=zeros(M*N+N,3);
s(1:M*N,1)=repmat([1:N]',[M,1]);
s(M*N+(1:N),1)=(1:N)'; s(M*N+(1:N),2)=(1:N)'; s(M*N+(1:N),3)=realmin('single');
ii=randperm(N); ii=ii(1:M);
netsimbest=-1e300; netsimbest=-Inf;
for sw=1:sweeps
    for m=1:M
        s((m-1)*N+1:m*N,2)=ii(m);
        s((m-1)*N+1:m*N,3)=S(1:N,ii(m))+realmin('single');
% 		C(1:N,ii(m)) = C(1:N,ii(m))+1;  % counting S-computations
    end;
     [idx,netsim,dpsim,expref]=apclustermex(sparse(s(:,1),s(:,2),s(:,3),N,N),p);
	fprintf('.');
    if netsim > netsimbest | isinf(netsimbest), % should allow updates if netsim=-Inf;
        netsimbest=netsim;
        idxbest=idx; dpsimbest=dpsim; exprefbest=expref;
        tmp1=unique(idx)';
        tmp2=setdiff([1:N],tmp1);
        ii=[tmp1,tmp2(randperm(length(tmp2)))];
        ii=ii(1:M);
    else
        tmp1=unique(idxbest)';
        tmp2=setdiff([1:N],tmp1);
        ii=[tmp1,tmp2(randperm(length(tmp2)))];
        ii=ii(1:M);
    end;
end;
idx=idxbest; netsim=netsimbest; dpsim=dpsimbest; expref=exprefbest;