function [state, cprobs] = sw_allall_bm(state, weights, biases, iters)
%SW_ALLALL_BM Swendsen-Wang algorithm for fully-connect Boltzmann machine
%
%     [state, cprobs] = sw_allall_bm(state, weights, biases, iters)
%
% Takes "iters" steps of the S-W Markov chain from state. The Boltzmann machine
% is a probability distribution over a vector of {0,1} variables. The
% distribution is a reparameterization of the classic {-1,+1} Ising model.
%
%     p(s) \propto \exp( -E(s;W,b) )
%
% Inputs:
%        state DxN Markov chain state: N independent D-dimensional vectors.
%                  Should be {0,1} variables. Use Ising functions for {+1,-1} vars.
%      weights DxD Every off-diagonal element weights(i,j) = W_{ij} contributes
%                      - W_{ij} s_i s_j
%                  to the energy. This means there are two terms for every (i,j)
%                  pair. To avoid "double counting" only put weight parameters
%                  into either the upper or the lower diagonal of the weights
%                  matrix, or put each weight parameter in both halves, but
%                  remember to divide them by two.
%       biases Dx1 or DxN
%                  The biases define the final part of the energy.
%                      -E = \sum_{i=1}^D \sum_{j=1}^D W_{ij} s_i s_j + \sum_i b_i s_i
%                  Each Markov chain can optionally have its own bias vector,
%                  which is useful in applications where the "bias" is an input
%                  from another layer of the model, different for each case.
%        iters 1x1 Number of Markov chain steps to take.
%
% Outputs:
%        state DxN Markov chain state after an update that leaves the
%                  Boltzmann machine 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

% Convert to Ising representation
ising_J = (weights+weights')/4;
ising_h = bsxfun(@plus, biases/2, sum(ising_J, 2));
ising_state = 2*state - 1;

% Run Ising code
[ising_state, cprobs] = sw_allall_ising(ising_state, ising_J, ising_h, iters);

% Convert back
state = (ising_state + 1) / 2;