% S = sw(W,T,q[,M,init]) Swendsen-Wang Monte Carlo sampling
%
% Produces a sequence of M states S in a Potts model with the Swendsen-Wang
% algorithm. Reference:
% 
%   Robert H. Swendsen and Jian-Sheng Wang: "Nonuniversal Critical Dynamics
%   in Monte Carlo Simulations", PRL 58(2):86-88 (Jan. 2, 1987).
%
% In:
%   W: NxN sparse matrix containing the edge values (affinities in
%      clustering, interactions in the Potts model); Wij >= 0. Note:
%      W is assumed triangular (i.e., if W is symmetric use triu(W,1)).
%   T: temperature; T > 0.
%   q: number of values for the Potts spins (labels).
%   M: number of sample states to produce. Default: 1.
%   init: initial state, one of: 'rnd' (random), '1' (constant) or an Nx1
%      vector containing the initial state. Default: '1'.
% Out:
%   S: NxM list of M sample states, each an Nx1 vector with values in 1..q.
%
% Any non-mandatory argument can be given the value [] to force it to take
% its default value.

% Copyright (c) 2004 by Miguel A. Carreira-Perpinan

function S = sw(W,T,q,M,init)

% ---------- Argument defaults ----------
if ~exist('M','var') | isempty(M) M = 1; end;
if ~exist('init','var') | isempty(init) init = '1'; end;
% ---------- End of "argument defaults" ----------

N = size(W,1);
[I,J,p] = find(W);	% List of edges
p = exp(-p/T);		% Probability to remove each edge

S = ones(N,M+1);

% Initial state
if isstr(init)
  if strcmp(init,'rnd')
    S(:,1) = floor(q*rand(N,1))+1;
  end
else
  S(:,1) = init;
end
% Ensure we use labels 1,2,...
q2 = sort(unique(S(:,1))); q3 = zeros(q,1); q3(q2) = 1:length(q2);
S(:,1) = q3(S(:,1));

for m=2:M+1
  
  % Build SW-graph: for each edge i~j connecting sites with the same label,
  % remove such edge with probability exp(-Wij/T).
  tmp = find(S(I,m-1) == S(J,m-1));	% Edges connecting same-label sites
  tmp = tmp(rand(size(tmp)) >= p(tmp));	% Surviving edges
  II = I(tmp); JJ = J(tmp);
  
  % Build SW-clusters: connected components of the SW-graph.
  [K, comp] = conncomp(N,II,JJ);
  
  % Generate new configuration S: assign a random label in 1..q to every
  % SW-cluster.
% $$$   for k=1:K
% $$$     S(comp(k).c,m) = floor(q*rand)+1;
% $$$   end
  % Ensure we use labels 1,2,...
  q1 = floor(q*rand(K,1))+1; q2 = sort(unique(q1));
  q3 = zeros(q,1); q3(q2) = 1:length(q2);
  for k=1:K
    S(comp(k).c,m) = q3(q1(k));
  end
  
end

