Last active
December 22, 2015 18:28
-
-
Save foxish/6512711 to your computer and use it in GitHub Desktop.
Percolator Problem: Monte Carlo simulation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
from pprint import pprint | |
import random | |
import time | |
NUM_ELEM = 600 | |
NUM_SD = NUM_ELEM*NUM_ELEM #single dim | |
dyns = None | |
matrix = None | |
experiment_attempts = 2 | |
class Cell: | |
def __init__(self, i): | |
self.parent = i | |
self.index = i | |
self.weight = 1 | |
def connect(self, node): | |
self.parent = node.index | |
self.weight += node.weight | |
def main(): | |
global dyns, matrix | |
ratio_sum = 0 | |
start_time = time.time() | |
for k in range(0, experiment_attempts): | |
num_times = 0 | |
# 2-d actual matrix | |
matrix = [[0 for i in range(NUM_ELEM)] for i in range(NUM_ELEM)] | |
# structure 1-d | |
dyns = [Cell(i) for i in range(NUM_SD)] | |
for j in range(1, NUM_ELEM): | |
connect(j, 0) | |
connect(NUM_SD-j-1, NUM_SD-1) | |
while(not is_percolated()): | |
randint1 = random.randint(0, NUM_ELEM - 1) | |
randint2 = random.randint(0, NUM_ELEM - 1) | |
if site_open(randint1,randint2): | |
num_times += 1 | |
ratio = float(num_times)/NUM_SD | |
print(ratio) | |
ratio_sum += ratio | |
print "Final: ", ratio_sum/experiment_attempts | |
print time.time() - start_time, "seconds" | |
def site_open(i, j): | |
if matrix[i][j] == 1: | |
return False | |
else: | |
matrix[i][j] = 1 | |
if i > 0 and matrix[i - 1][j] == 1: | |
connect(compute_offset(i, j), compute_offset(i - 1, j)) | |
if j > 0 and matrix[i][j - 1] == 1: | |
connect(compute_offset(i, j), compute_offset(i, j - 1)) | |
if i < (NUM_ELEM - 1) and matrix[i+1][j] == 1: | |
connect(compute_offset(i, j), compute_offset(i + 1, j)) | |
if j < (NUM_ELEM - 1) and matrix[i][j + 1] == 1: | |
connect(compute_offset(i, j), compute_offset(i, j + 1)) | |
return True | |
def compute_offset(i, j): | |
return (i*NUM_ELEM + j) | |
def connect(i, j): | |
var1 = root(i) | |
var2 = root(j) | |
if(var1 == var2): | |
return False | |
if(var1.weight <= var2.weight): | |
var1.connect(var2) | |
else: | |
var2.connect(var1) | |
def root(i): | |
node = dyns[i] | |
while(node.parent!=node.index): | |
node = dyns[node.parent] | |
return node | |
def is_connected(i, j): | |
return root(i) == root(j) | |
def get_dyns(): | |
ret = {} | |
i = 0 | |
for dyn in dyns: | |
ret[i] = dyn.parent.index | |
i += 1 | |
return ret | |
def is_percolated(): | |
return is_connected(0, NUM_SD - 1) | |
if __name__ == '__main__': | |
main() | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment