adapt_apcluster.m

function [labels,NC,labelid,i,Sp,Slam,NCfix,netsim,dpsim,expref,idx] = ...
    adapt_apcluster(data,dtype,pvalues,folds,adapt,varargin)

adapt = adapt+1;            % 0: turning to original AP
if adapt < 2
if strcmp(dtype, 'euclidean') || (isnumeric(dtype) && dtype == 1)
   [Dist, dmax] = similarity_euclid(data,2);
elseif strcmp(dtype, 'correlation') || (isnumeric(dtype) && dtype == 2)
   Dist = 1-(1+similarity_pearson(data'))/2;
   dmax = 1;
end

nrow = size(Dist,1);
nap = nrow*nrow-nrow;
s = zeros(nap,3);
j=1;

for i=1:nrow
   for k = [1:i-1,i+1:nrow]
     s(j,1) = i;
     s(j,2) = k; 
     s(j,3) = -Dist(i,k);
     j = j+1;
   end;
 end;
Dist = [];
else
  s = data;
  data =[];
  %pmin = -min(s(:,3));
end

pfixed = 0;
dn = find(s(:,3)>-realmax);
pmedian = median(s(dn,3));                     % Set preference to median similarity
pstep = folds*pmedian;                        % decreasing step of p
if length(pvalues) < 1
    pvalues = pmedian*0.5;
else
    pfixed = 1;
end

% Handle arguments to function
if nargin<2
    error('Too few input arguments');
else
    maxits=500; convits=50; lam=0.5; plt=0; details=0; nonoise=0;
    i=1;
    while i<=length(varargin)
        if strcmp(varargin{i},'plot')
            plt=1; i=i+1;
        elseif strcmp(varargin{i},'details')
            details=1; i=i+1;
		elseif strcmp(varargin{i},'sparse')
			[idx,netsim,dpsim,expref]=apcluster_sparse(s,pvalues,varargin{:});
			return;
        elseif strcmp(varargin{i},'nonoise')
            nonoise=1; i=i+1;
        elseif strcmp(varargin{i},'maxits')
            maxits=varargin{i+1};
            i=i+2;
            if maxits<=0 
                error('maxits must be a positive integer'); 
            end;
        elseif strcmp(varargin{i},'convits')
            convits=varargin{i+1};
            i=i+2;
            if convits<=0
                error('convits must be a positive integer'); 
            end;
        elseif strcmp(varargin{i},'dampfact')
            lam=varargin{i+1};
            i=i+2;
            if (lam<0.5)||(lam>=1)
                error('dampfact must be >= 0.5 and < 1');
            end;
        else i=i+1;
        end;
    end;
end;

if lam>0.9
    fprintf('\n*** Warning: Large damping factor in use. Turn on plotting\n');
    fprintf('    to monitor the net similarity. The algorithm will\n');
    fprintf('    change decisions slowly, so consider using a larger value\n');
    fprintf('    of convits.\n\n');
end;

% Check that standard arguments are consistent in size
if length(size(s))~=2 
    error('s should be a 2D matrix');
elseif length(size(pvalues))>2 
    error('pvalues should be a vector or a scalar');
elseif size(s,2)==3
    tmp=max(max(s(:,1)),max(s(:,2)));
    if length(pvalues)==1 
        N=tmp; 
    else
        N=length(pvalues);
    end;
    if tmp>N
        error('data point index exceeds number of data points');
    elseif min(min(s(:,1)),min(s(:,2)))<=0
        error('data point indices must be >= 1');
    end;
elseif size(s,1)==size(s,2)
    N=size(s,1);
    if (length(pvalues)~=N) && (length(pvalues)~=1)
        error('pvalues should be scalar or a vector of size N');
    end;
else error('s must have 3 columns or be square'); 
end;

% Construct similarity matrix
if N>3000
    fprintf('\n*** Warning: Large memory request. Consider activating\n');
    fprintf('    the sparse version of APCLUSTER.\n\n');
end;
if size(s,2)==3
    S=-Inf*ones(N,N); 
    for j=1:size(s,1)
        S(s(j,1),s(j,2))=s(j,3); 
    end;
else
    S=s;
end;
s = [];

% In case user did not remove degeneracies from the input similarities,
% avoid degenerate solutions by adding a small amount of noise to the
% input similarities
if ~nonoise
    rns=randn('state'); 
    randn('state',0);
    S=S+(eps*S+realmin*100).*rand(N,N);
    randn('state',rns);
end;

% Place preferences on the diagonal of S
if length(pvalues)==1 
    for i=1:N 
        S(i,i)=pvalues; 
    end;
else
    for i=1:N 
        S(i,i)=pvalues(i); 
    end;
end;

% Allocate space for messages, etc
dS=diag(S); 
A=zeros(N,N); 
R=zeros(N,N); 
t=1;
if plt 
    netsim=zeros(1,maxits+1); 
end;
if details
    idx=zeros(N,maxits+1);
    netsim=zeros(1,maxits+1); 
    dpsim=zeros(1,maxits+1); 
    expref=zeros(1,maxits+1); 
end;



% Initialization
dn=0; i=0;
stoptimes = max([maxits/10 2000]); % Stop if K's unchanging times >= stoptimes
if pfixed
    stoptimes = convits;
end
Hstop = zeros(N,stoptimes);  % recording unchanging times for stop condition
Hconvits = zeros(N,convits); % recording convergence at each K
Tdelay = 10;                            % delay time for detecting variation of K
Hdelay = Tdelay;                    % counting delay
Hconverg = 0;                      % whether K examplars convergence
nhalf = round(0.3*convits);
Hdelay2 = Tdelay;
Hconvhalf = zeros(N,nhalf);  % recording K at half convergence
Hsavehalf = 0;
Hn1 = 0; Hn2 = 0;
Kmean=0; Kdown=0; 
Wstart = max([100 round(convits/2)]); % start oscillation monitoring
wsize = 40;                  % monitoring window size for detecting oscillations
Kunchange=0;                 % recording variation of K
Kocil = ones(1,wsize);              % recording K's non-oscillations
Noscil = wsize+10;               % statistic of non-oscillations, initial value
Svib = 0;                               % counting potential oscillations
Hvib = 10; Tvib = 2;                    % counting oscillations, initial values
Hguid = 1;                            % label to start solution recording
Sprefer = pvalues;                     % asigned to diag of S for realizing p decreasing
astep = pstep;                      % adjustable step instead of fixed pstep
Kset = []; Kold = 0;              % recording K
Kfix = 0; nKfix = 0;             % used for speeding up
Kmax = 0; nfix = 0;

% Execute parallel affinity propagation updates
while ~dn
    i=i+1;

    % Compute responsibilities
    AS=A+S; 
    [Y,I]=max(AS,[],2); 
    for k=1:N 
        AS(k,I(k))= -realmax;
    end;
    [Y2,I2]=max(AS,[],2);
    AS = [];
    Rold=R;
    R=S-repmat(Y,[1,N]);
    for k=1:N 
        R(k,I(k))=S(k,I(k))-Y2(k); 
    end;
%     R=(1-lam)*R+lam*Rold;
    Rold = lam*Rold;
    R = (1-lam)*R;
    R = R+Rold;
    Rold = [];

    % Compute availabilities
    Rp=max(R,0); 
    for k=1:N 
        Rp(k,k)=R(k,k); 
    end; 
    Aold=A;
%     A=repmat(sum(Rp,1),[N,1])-Rp;
    A = sum(Rp,1);
    A = repmat(A,[N,1]);
    A = A-Rp;
    Rp = [];
    dA=diag(A);
    A=min(A,0);
    for k=1:N
        A(k,k)=dA(k);
    end;
%     A=(1-lam)*A+lam*Aold;
    Aold = lam*Aold;
    A = (1-lam)*A;
    A = A+Aold;
    Aold = [];

    % Check for convergence
    E = ((diag(A)+diag(R))>0);
    Hconvits(:,mod(i-1,convits)+1) = E;
    K=sum(E);
    Kset(i) = K;
    newp = Sprefer(1);
    newlam = lam;
    Hstop(:,mod(i-1,stoptimes)+1) = E;
    Hconvhalf(:,mod(i-1,nhalf)+1) = E;
    if mod(i,100) == 1 || i == maxits
        fprintf('** running at iteration %d, K = %d\n', i,K);
    end
    
    Hsave = 0; Hsave1 = 0; Hsave2 = 0; Hsave3 = 0;
    if i>=Wstart || i>=maxits
        se = sum(Hconvits,2);
        se1 = sum(se==convits);
        se2 = sum(se==0);
        unconverged = (se1+se2) ~= N;
        Hconverg = ~unconverged; 
        se = sum(Hstop,2);
        se1 = sum(se==stoptimes);
        se2 = sum(se==0);
        if (se1+se2) == N || i == maxits
            dn=1;
            if (se1+se2) == N
               Hsave1 = 1;
            end
        end
        se = sum(Hconvhalf,2);
        se1 = sum(se==nhalf);
        se2 = sum(se==0);
        Hsavehalf = (se1+se2) == N;
        Hsavehalf = Hsavehalf && Hguid == 2;
    end
    
 if adapt   
    if i > 5
      Kmean(i) = mean(Kset(i-5:i));        % covering 6 points
      Kdown(i) = Kmean(i)-Kmean(i-1) < 0;  % covering 7 points
      if Hguid == 2
          Kdown(i) = Kdown(i) && K <= Kold;
      end
      Kunchange(i) = sum(abs(Kset(i)-Kset(i-5:i-1)));
      Kocil(:,mod(i-1,wsize)+1) = Kdown(i) || Kunchange(i) == 0;
      Noscil = sum(Kocil);                     % non- oscillations
    end
    
    % reducing parameter pvalues to yield smaller NC when unchanging
    if Hconverg
      Hdelay = Hdelay+1;
      if Hdelay >= Tdelay
        Hsave1 = 1;        % starting pvalues reduction & result saving
        Hdelay = 0;
        Hn1 = Hn1+1;
        if K == Kfix
            nKfix = nKfix+1;
        else
            nKfix = 0;
        end
        Kfix = K;        
        stepfold = sqrt(K+50)/10;
        pstep = folds*pmedian/stepfold;
        if nKfix > 1
            astep = nKfix*pstep;                 % speeding up
        else
            astep = pstep;
        end
      end
    elseif Hsavehalf
        Hdelay2 = Hdelay2+1;
        if Hdelay2 >= Tdelay
            Hsave2 = 1;        % starting pvalues reduction & result saving
            Hdelay2 = 0;
            Hn2 = Hn2+1;
        end
    end
    if ~Hconverg
        Hn1 = 0;
        Hdelay = 0;
    end
    if ~Hsavehalf
        Hn2 = 0;
        Hdelay2 = 0;
    end
 
    if ((K == 1 || K == 2) && Hsave1)   % avoiding influence of oscillation
       dn = 1;
       unconverged = 0;
    end

    if Hguid == 1 && Hsave1
        Hguid = 2;                               % starting guidance
        labels = zeros(N,K);
        labelid = zeros(N,K);
        NC = zeros(1,K);
        NCfix = zeros(1,K);
        Sp = zeros(1,K);
        Slam = zeros(1,K);
        Kmax = K;
        stepfold = sqrt(Kmax+50)/10;
        pstep = folds*pmedian/stepfold;
    end
      
    if Hsave1                                     % K >= Kold unchanging, K rise
      Svib = 0;
      if ~pfixed
          Sprefer = Sprefer+astep;           % reducing pvalues
          if length(pvalues)==1
              for k=1:N
                  S(k,k) = Sprefer;
              end
          else
              for k=1:N
                  S(k,k) = Sprefer(k);
              end
          end
      end
      
    else
       Svib = Svib+1;                             % counting potential oscillations
       HSvib = (Svib > wsize && Noscil < 0.66*wsize) || Svib > 150;
       HSvib = HSvib && i > Wstart;
       if HSvib                          % = & ! should be 2/3, otherwise oscillations
          Hvib = Hvib+1;
          if Hvib > 10
            lam = max([0.7 lam]);
          elseif Hvib >= 1
            if Tvib >= 3
                if lam >= 0.9
                    lam = min([0.98 0.025+lam]);
                    if lam >= 0.95 && mod(i,9) == 2
                        rns=randn('state');
                        randn('state',0);
                        S=S+(eps*S+realmin*1000).*rand(N,N);
                        randn('state',rns);
                        fprintf(' # A small amount of noise is added\n');
                    end
                end
            else
            lam = min([0.9 0.05+lam]);
            end
            if lam >= 0.85
                Tvib = Tvib+1;
                if ~pfixed
                    if Hguid == 2 && Kold
                       if Kmax <1
                          stepfold = 2;
                       else
                          stepfold = max([3/(sqrt(Kmax)/10+0.4) 1]);
                       end
                       Kvar = 2*sqrt(std(Kset(i-49:i)));
                       astep = min([0.8*Kvar+0.2*Tvib stepfold])*pstep;
                       %0.5*((Kvar+min([Tvib 4]))*pstep); max([sqrt(abs(K-Kold))])
                    else
                       astep = min([Tvib 2])*pstep;
                    end
                    Sprefer = Sprefer+astep;  % escaping oscillations
                    if length(pvalues)==1
                        for k=1:N
                            S(k,k) = Sprefer;
                        end
                    else
                        for k=1:N
                            S(k,k) = Sprefer(k);
                        end
                    end
                    fprintf(' # Escaping oscillation turns on\n');
                end
            end
          end
          Hvib = 0;
          Svib = 0;
          fprintf(' # Damping factor is increased to %g\n', lam);
       else
           Tvib = max([Tvib-0.002 0.98]);
           if lam > 0.9 && Tvib < 1
              lam = max([lam-0.0001 0.5]);
           end
       end
       
    end
 end
 
   if Hguid >= 2 && (K < Kold && K > 1 || Kmean(i) == Kmean(i-1))
      Hsave3 = 1;              % catching decreasing K
   end
   if Hsave1 || Hsave2 || Hsave3
       Hsave = 1;
   end
   Kold = K;

    % Handle plotting and storage of details, if requested
    if plt || details || Hsave
        if K==0
            tmpnetsim=nan; tmpdpsim=nan; tmpexpref=nan; tmpidx=nan;
        else
            I=find(E); 
            [tmp c]=max(S(:,I),[],2); 
            c(I)=1:K;
            if Hsave || dn == 1
               if Hsave1
                  nfix = Hn1*Tdelay+convits;
               elseif Hsave2
                  nfix = Hn2*Tdelay+nhalf;
               elseif Kmean(i) == Kmean(i-1)
                  nfix = 6;
                  if Kmean(i) == Kmean(i-5)
                    nfix = 10;
                  end
               else
                  nfix = 1;
               end
               
               if (K <= Kmax && nfix > NCfix(K)) || (K > Kmax && nfix >= 10)
                   NCfix(K) = nfix;
                   labels(:,K) = c;
                   labelid(:,K) = I(c);
                   NC(K) = K;
                   Sp(:,K) = newp; 
                   Slam(:,K) = newlam;
               end
               if K > Kmax
                   Kmax = K;
                   if length(NCfix) < K
                       NCfix(K) = 0;
                   end
               end
            else
              tmpidx=I(c);
              tmpnetsim=sum(S((tmpidx-1)*N+[1:N]'));
              tmpexpref=sum(dS(I)); 
              tmpdpsim=tmpnetsim-tmpexpref;
            end
        end
    end;
    if details
        netsim(i)=tmpnetsim; dpsim(i)=tmpdpsim; expref(i)=tmpexpref;
        idx(:,i)=tmpidx;
    end;
    if plt
        netsim(i)=tmpnetsim;
        figure(1); 
        tmp=1:i; 
        tmpi=find(~isnan(netsim(1:i)));
        subplot(2,1,1)
        plot(tmp(tmpi),netsim(tmpi),'r-');
        xlabel('# Iterations');
        ylabel('Fitness (net similarity)');
        subplot(2,1,2)
        tmp=max([1 i-3*wsize]):i; 
        plot(tmp,Kset(tmp),'b-');
        xlabel('# Iterations');
        ylabel('Number of clusters');
        drawnow; 
    end;
end;


fprintf(' # Programs run over at K= %d\n', K);
I=find(diag(A+R)>0); 
K=length(I);                                      % Identify exemplars
if K>0
    [tmp c]=max(S(:,I),[],2); 
    c(I)=1:K;                                       % Identify clusters
    % Refine the final set of exemplars and clusters and return results
    for k=1:K 
        ii=find(c==k);
        [y j]=max(sum(S(ii,ii),1)); 
        I(k)=ii(j(1)); 
    end
    [tmp c]=max(S(:,I),[],2); 
    c(I)=1:K; 
    tmpidx=I(c);
    tmpnetsim=sum(S((tmpidx-1)*N+[1:N]')); 
    tmpexpref=sum(dS(I));
    labels(:,K) = c;
    labelid(:,K) = tmpidx;
    NC(K) = K;
    NCfix(K) = nfix;
    Sp(:,K) = newp;
    Slam(:,K) = newlam;
else
    tmpidx=nan*ones(N,1); 
    tmpnetsim=nan; 
    tmpexpref=nan;
end;
if details
    netsim(i+1)=tmpnetsim; netsim=netsim(1:i+1);
    dpsim(i+1)=tmpnetsim-tmpexpref; dpsim=dpsim(1:i+1);
    expref(i+1)=tmpexpref; expref=expref(1:i+1);
    idx(:,i+1)=tmpidx; idx=idx(:,1:i+1);
else
    netsim=tmpnetsim; 
    dpsim=tmpnetsim-tmpexpref;
    expref=tmpexpref; 
    idx=tmpidx;
end;

if length(NC) > 1
   NC(1) = 0;
end
if length(NC) < 1
  Sp = [];
  Slam = [];
  NC = 0;
  NCfix = 0;
  labels = ones(nrow,1);
  labelid = ones(nrow,1);
else
  S = find(NC);
  Sp = Sp(:,S);
  Slam = Slam(:,S);
  NC = NC(:,S);
  labels = labels(:,S);
  labelid = labelid(:,S);
  NCfix = NCfix(:,S);
end

if plt || details
    fprintf('\nNumber of identified clusters: %d\n',K);
    fprintf('Fitness (net similarity): %f\n',tmpnetsim);
    fprintf('  Similarities of data points to exemplars: %f\n',dpsim(end));
    fprintf('  Preferences of selected exemplars: %f\n',tmpexpref);
    fprintf('Number of iterations: %d\n\n',i);
end;
if unconverged
    fprintf('\n*** Warning: Algorithm did not converge at K = %g ! \n', NC(1));
    fprintf('    The similarities may contain degeneracies - add noise to\n');
    fprintf('    the similarities to remove degeneracies. To monitor the net\n');
    fprintf('    similarity, activate plotting. Also, consider increasing\n');
    fprintf('    maxits and if necessary dampfact.\n\n');
end;