Percolate

Hur perforerad behöver en N x N kub vara för att det ska finnas en väg från toppen till botten?

Den röda pricken i grafen anger vilken punkt som fick simulationen att "perkolera".


In [3]:
# percolate:
import numpy as np

from bokeh.plotting import figure
from bokeh.io import push_notebook, output_notebook, reset_output, show

import sys
import random

output_notebook()

p = figure(plot_width=400, plot_height=400, tools='pan,box_zoom,wheel_zoom,save')

N = 40
opened = np.zeros(N*N)
sz = np.ones(N*N+2)
uf = np.arange(N*N+2)

x = [x for r in range(N) for x in range(N)]
y = [N-y for y in range(N) for x in range(N)]
colors = ['olive' for x in range(N*N)]

last_x = []
last_y = []

p.square(x, y, size=(400 // N)*0.75, color = colors, alpha = 0.5)
p.circle(last_x, last_y, size=(400 // N) * 0.5, color = 'red', alpha = 0.5)

t = show(p, notebook_handle=True)


def findRoot(p):
    parent = uf[p]
    if parent == p:
        return p
    
    return findRoot(parent)

def union(p, q):
    # todo: verify that p and q are within boundaries
    p_parent = findRoot(p)
    q_parent = findRoot(q)
    if p_parent != q_parent:
        if sz[p_parent] > sz[q_parent]:
            uf[q_parent] = uf[p_parent]
            sz[p_parent] += sz[q_parent]
        else:
            uf[p_parent] = uf[q_parent]
            sz[q_parent] += sz[p_parent]

    #cu_uf[cu_findRoot(p)] = cu_uf[q] <- we do not use this one either
    #cu_uf[p] = cu_uf[q] <- this one is wrong!

def isConnected(p, q):
    return findRoot(p) == findRoot(q)


def openCell(row, col):
    if not isOpen(row, col):
        idx = row*N + col
        opened[idx] = 1
        
        if row == 0:
            union(idx, N*N)
        elif isOpen(row-1, col):
            union(idx, idx-N)

        if col > 0 and isOpen(row, col-1):
            union(idx, idx-1)

        if col < N-1 and isOpen(row, col+1):
            union(idx, idx+1)

        if row == N-1:
            union(idx, N*N+1)
        elif isOpen(row+1, col):
            union(idx, idx+N)

        color = 'white'
        if isConnected(idx, N*N):
            color = 'blue'
        colors[idx] = color
        
        for i in range(N*N):
            if isConnected(i, N*N):
                colors[i] = 'blue'
        push_notebook(handle = t)

def isOpen(row, col):
    return opened[row*N+col] == 1

def isFull(row, col):
    return opened[row*N+col] == 0

def numberOfOpenSites():
    return np.sum(opened)

def percolates():
    return isConnected(N*N, N*N+1)

while not percolates():
    idx = random.randint(0, N*N-1)
    openCell(idx // N, idx % N)

last_x.append(idx % N)
last_y.append(N - (idx // N))


push_notebook(handle = t)

reset_output()

print(percolates(), numberOfOpenSites(), numberOfOpenSites() / (N*N))


Loading BokehJS ...
True 995.0 0.621875

In [ ]: