calculating percolation probability

well this is a nice application of unionfind algorithm. but i think I’m doing wrong as it’s still very slow for just 100 by 100 matrix! the tests are taken with 10 by 10 and amazingly it correct results!

#include <functional>
#include <numeric>
#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <cstring>
#include <cassert>
#include <vector>
#include <list>
#include <map>
#include <set>
#include <deque>
#include <queue>
#include <bitset>
#include <sstream>
#include <algorithm>
#include "unionfind.h"
double monte_carlo_simulation()
{
int nmax=10;
unionfind uf(nmax*nmax + 2);
for (int i = 1; i <= nmax; i++) {
uf.unionof(0,i);
uf.unionof(nmax*nmax -nmax +i,nmax*nmax +1);
}
bool percolate_board[nmax*nmax +2];
percolate_board[0] = true;
percolate_board[nmax*nmax +1] = true;
for (int i = 1; i < nmax*nmax +1; i++)
percolate_board[i] = false;
int calls = 0;
while(not uf.isconnected(0,nmax*nmax +1))
{
int tmp = rand()%(nmax*nmax -calls) + 1;
calls++;
int count = 0;
int skip_count = 0;
while (count < tmp)
{
if(percolate_board[skip_count + count + 1])
skip_count++;
else
count++;
}
tmp += skip_count;
percolate_board[tmp] = true;
if (tmp < nmax) {
if (percolate_board[tmp-1]) uf.unionof(tmp,tmp-1);
if (percolate_board[tmp+1]) uf.unionof(tmp,tmp+1);
if (percolate_board[tmp+nmax]) uf.unionof(tmp,tmp+nmax);
continue;
}
if (tmp == nmax) {
if (percolate_board[tmp-1]) uf.unionof(tmp,tmp-1);
if (percolate_board[tmp+nmax]) uf.unionof(tmp,tmp+nmax);
continue;
}
if (tmp > nmax*nmax -nmax +1) {
if (percolate_board[tmp-1]) uf.unionof(tmp,tmp-1);
if (percolate_board[tmp+1]) uf.unionof(tmp,tmp+1);
if (percolate_board[tmp-nmax]) uf.unionof(tmp,tmp-nmax);
continue;
}
if (tmp == nmax*nmax -nmax +1) {
if (percolate_board[tmp+1]) uf.unionof(tmp,tmp+1);
if (percolate_board[tmp-nmax]) uf.unionof(tmp,tmp-nmax);
continue;
}
if ((tmp-1)%nmax == 0) {
if (percolate_board[tmp+1]) uf.unionof(tmp,tmp+1);
if (percolate_board[tmp+nmax]) uf.unionof(tmp,tmp+nmax);
if (percolate_board[tmp-nmax]) uf.unionof(tmp,tmp-nmax);
continue;
}
if (tmp%nmax == 0) {
if (percolate_board[tmp-1]) uf.unionof(tmp,tmp+1);
if (percolate_board[tmp+nmax]) uf.unionof(tmp,tmp+nmax);
if (percolate_board[tmp-nmax]) uf.unionof(tmp,tmp-nmax);
continue;
}
if (percolate_board[tmp-1]) uf.unionof(tmp,tmp-1);
if (percolate_board[tmp+1]) uf.unionof(tmp,tmp+1);
if (percolate_board[tmp+nmax]) uf.unionof(tmp,tmp+nmax);
if (percolate_board[tmp-nmax]) uf.unionof(tmp,tmp-nmax);
}
int percolation_count = 0;
for (int i = 1; i < nmax*nmax +1; i++) {
if(percolate_board[i])
percolation_count++;
}
return percolation_count;
}
int main(int argc, const char *argv[])
{
double total;
srand(time(0));
for (int i = 0; i < 10000; ++i)//10000 experiments
{
total += monte_carlo_simulation()/100;
}
std::cout << double(total)/10000 << std::endl ;
return 0;
}
#include <map>
#include <set>
#include <list>
#include <cmath>
#include <ctime>
#include <deque>
#include <queue>
#include <stack>
#include <bitset>
#include <cstdio>
#include <limits>
#include <vector>
#include <cstdlib>
#include <numeric>
#include <sstream>
#include <iostream>
#include <algorithm>
#include <iterator>
class unionfind
{
private:
//long long int maxsize=100000;
long long int *id;
long long int *wt;
long long int n;
// recursive and has backtracking
long long int findroot(long long int& n);
public:
// new long long int [n]
unionfind(long long int );
// delete
~unionfind();
// through findroot
bool isconnected(long long int , long long int );
// using findroot and making union by taking weights into consideration
void unionof(long long int ,long long int );
// just for testing
void printufstate();
};
unionfind::unionfind(long long int n=100000)
{
this->n = n;
id = new long long int [n];
wt = new long long int [n];
for (long long int i = 0; i < n; ++i)
id[i]=i;
for (long long int i = 0; i < n; ++i)
wt[i]=1;
}
unionfind::~unionfind()
{
delete[] id;
delete[] wt;
}
long long int unionfind::findroot(long long int& n)
{
if (n==id[n])
return n;
n=findroot(id[n]);
return n;
}
bool unionfind::isconnected(long long int p, long long int q)
{
long long int tmp1=id[findroot(p)];
long long int tmp2=id[findroot(q)];
if (tmp1 == tmp2)
return true;
else
return false;
}
void unionfind::unionof(long long int p,long long int q)
{
long long int rootp=id[findroot(p)];
long long int rootq=id[findroot(q)];
if (rootp!=rootq)
{
if (wt[rootq]>wt[rootp]) // checking their weights
{
id[rootp]=rootq;
wt[rootq]+=wt[rootp];
wt[rootp]=0;//it's useless though
}
else
{
id[rootq]=rootp;
wt[rootp]+=wt[rootq];
wt[rootq]=0;
}
}
}
void unionfind::printufstate()
{
for (long long int i = 0; i < this->n; ++i)
std::cout << id[i] << " ";
std::cout << '\n';
}
view raw unionfind.h hosted with ❤ by GitHub

Leave a comment