#include
#include
/* Swendsen-Wang helper code for binary models with pairwise correlations (Known
* as Ising models, Boltzmann machines, binary potts models, etc.)
*
* See the calling Matlab code for the bigger picture of the S-W algorithm. This
* C code works out the clusters given the bonds provided by the Matlab code. It
* then works out the probability of flipping each cluster given biases passed
* by Matlab, makes a decision for each node and passes out both the decision
* and the conditional probability. The Matlab code can then turn this flip
* information into actual settings of the Ising model.
*
* Here the graph of variables is assumed to be all-all connected. And the
* implementation is fairly naive. If many of the connections are weak, it would
* be possible to reduce the computational complexity of each update
* dramatically with a little pre-processing. Note that many of the connections
* better be weak, or percolation will be above a phase transition, the system
* will all lock together and the sampler will do nothing interesting!
*
* bonds and flips are stored as doubles throughout because Matlab makes
* dealing with non-doubles difficult and I couldn't be bothered with
* conversions
*/
int
findroot(int *ptr, int i)
/* Adapted from:
* A fast Monte Carlo algorithm for site or bond percolation
* M. E. J. Newman and R. M. Ziff
* arXiv:cond-mat/0101295 v2 8th April 2001
*
* See that for an explanation of ptr[]
*/
{
if (ptr[i]<0)
return i;
return ptr[i] = findroot(ptr, ptr[i]);
}
void
percolate(double* bonds, int *ptr, double* biases, int nn)
/* This routine populates ptr[], from which the clustering implied by the bonds
* can be colored more easily. The biases of the nodes are destructively accumulated:
* the root of each cluster has the total bias in biases, any other entries are
* arbitrary values.
*/
{
int i, j;
int base_idx;
int l, m;
for (i=0; i