function [state, cprobs] = sw_allall_ising(state, ising_J, ising_h, iters)
%SW_ALLALL_ISING Swendsen-Wang algorithm for fully-connect Ising model
%
% [state, cprobs] = sw_allall_ising(state, ising_J, ising_h, iters)
%
% Inputs:
% state DxN N D-dimensional {+1,-1} vectors
% ising_J DxD This matrix should be symmetric, contributing energy:
% - 0.5 \sum_i \sum_j J_{ij} s_i s_j
% That is both J_{ij} and J_{ji} contain the coupling between
% s_i and s_j. The 0.5 takes care of double-counting.
% ising_h Dx1 or DxN
% field terms \sum_i h_i s_i
% There can be one common h, or one for each vector.
% iters 1x1 Number of Markov chain steps to take.
%
% Outputs:
% state DxN Markov chain state after an update that leaves the
% Ising model distribution(s) invariant.
% cprobs DxN Marginal probabilities that each site would have been set to
% the "+1" state conditioned on the bond variables in the last update.
% Iain Murray, August 2007, November 2009
potts_J = 2*ising_J;
if isvector(ising_h)
ising_h = repmat(ising_h, 1, size(state, 2));
end
bond_probs = 1 - exp(-abs(potts_J));
flips = zeros(size(state));
flipprobs = zeros(size(state));
for it = 1:iters
for nn = 1:size(state,2)
% Find bond variables
bonds_allowed = ((state(:,nn)*state(:,nn)').*potts_J > 0);
bonds = (rand(size(bond_probs)) < bond_probs) .* bonds_allowed;
% Work out which spins to flip (and with what probability):
flip_bias = 2*ising_h(:,nn).*state(:,nn); % bonus associated with sticking rather than flipping
[flips(:,nn), flipprobs(:,nn)] = sw_allall_flip_conditionals(bonds, flip_bias);
end
flips = 2*flips - 1;
if it == iters
% This book-keeping has to be done before updating the state
cprobs = flipprobs;
downs = (state < 0);
cprobs(downs) = 1 - cprobs(downs);
end
% New state
state = state.*flips;
end