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