percolation

percolation.py

class Site:
    def __init__(self, parent=-1, top=False, bottom=False):
        self.open = False
        self.parent = parent
        self.top = top
        self.bottom = bottom

    def makeOpen(self):
        self.open = True

    def isOpen(self):
        return self.open

    def percolates(self):
        return self.top and self.bottom

    def cardinality(self):
        if not isinstance(self.parent, list):
            return -self.parent

    def __str__(self):
        return '%s %s %s %s' % (self.open, self.parent, self.top, self.bottom)


class MonteCarloPercolation:
    def __init__(self, n):
        self.n = n
        self.sites = []
        for i in range(n + 1):
            row = []
            for j in range(n + 1):
                site = Site()
                if i == 1:
                    site.top = True
                elif i == n:
                    site.bottom = True
                row.append(site)
            self.sites.append(row)
        pass

    def open(self, i, j):
        self.sites[i][j].makeOpen()
        pass

    def isOpen(self, i, j):
        return self.sites[i][j].isOpen()
        pass

    def isFull(self, i, j):
        pass

    def percolates(self, i, j):
        parent = self.find(i, j)
        return self.sites[parent[0]][parent[1]].percolates()
        pass

    def rand(self):
        from random import randint
        while True:
            a, b = randint(1, self.n), randint(1, self.n)
            if not self.isOpen(a, b):
                return [a, b]
        pass

    def simulates(self):
        cnt = 0
        while True:
            cnt += 1
            x, y = self.rand()
            self.open(x, y)
            sides = [[x - 1, y], [x + 1, y], [x, y - 1], [x, y + 1]]
            for side in sides:
                if not self.isValid(*side):
                    continue
                if self.isOpen(*side) and self.find(x, y) != self.find(*side):
                    self.union([x, y], side)
                    if self.percolates(x, y):
                        return cnt / (self.n * self.n)
        pass

    def isValid(self, x, y):
        return 1 <= x <= self.n and 1 <= y <= self.n
        pass

    def find(self, x, y):
        if isinstance(self.sites[x][y].parent, list):
            self.sites[x][y].parent = self.find(*self.sites[x][y].parent)
            return self.sites[x][y].parent
        else:
            return [x, y]
        pass

    def union(self, site_a_index, site_b_index):
        pa_index = self.find(*site_a_index)
        pb_index = self.find(*site_b_index)
        parent_a = self.sites[pa_index[0]][pa_index[1]]
        parent_b = self.sites[pb_index[0]][pb_index[1]]
        if parent_a.cardinality() > parent_b.cardinality():
            parent_b.parent = pa_index
            parent_a.top |= parent_b.top
            parent_a.bottom |= parent_b.bottom
        else:
            parent_a.parent = pb_index
            parent_b.top |= parent_a.top
            parent_b.bottom |= parent_a.bottom
        pass


class PercolationStats:
    def __init__(self, n, trials):
        self.n = n
        self.trials = trials
        self.thresholds = []
        for i in range(trials):
            self.thresholds.append(MonteCarloPercolation(n).simulates())

    def mean(self):
        from statistics import mean
        return mean(self.thresholds)

    def stddev(self):
        from statistics import stdev
        return stdev(self.thresholds)

    def confidenceLo(self):
        return self.mean() - self.marginOfError()
        pass

    def confidenceHi(self):
        return self.mean() + self.marginOfError()
        pass

    def marginOfError(self):
        from math import sqrt
        return 1.96 * self.stddev() / sqrt(self.trials)

def main(n,trials):
    test = PercolationStats(n, trials)
    print("mean                    = %.6f" % test.mean())
    print("stddev                  = %.17f" % test.stddev())
    print("95%% confidence interval = %.17f, %.17f" % (test.confidenceLo(), test.confidenceHi()))


def time_wrapper(n, trials):
    import time
    start = time.clock()
    print("--------------%d : %d--------------"%(n,trials))

    main(n,trials)

    end = time.clock()
    print("time used               = %.2f s"%(end-start))
    print()

if __name__ == '__main__':
    tests = [
        [20,100],
        [20,500],
        [100,100],
        [100,300],
        [100,3000],
        [1000,3000],
        [2,100000]
    ]
    for test in tests:
        time_wrapper(*test)