Programming Assignment 1: Percolation

通過蒙特卡洛模擬方法來估計滲流閾值。

Percolation. 給一個有隨機分布的絕緣和金屬材料的組成的復合系統。例如我們想知道哪些部分必須是金屬材料才能讓這個復合系統是一個電導體。或者在一個多孔的地形，在表面有水或者油，在什么情況下水或者油能夠從最表面滲透到最底層。科學家把這種過程的模型叫做Percolation。

The model. 在Assignment中，用一個NxN的格子表示percolation系統，每一個格子是打開或者關閉，打開是白色關閉是黑色。如果一個格子是full，首先他必須是打開額，然后表示從最頂上通過相連(4方向)的打開的格子可以滲透到這個位置。當一個系統是percolates，表示能從最頂層滲透到最底層，也就是說，最底層存在打開的格子是full。

 

The Problem. 研究人員對一下的問題感興趣，如果每一個格子是獨立的，並且被打開的概率為p，那么系統percolates的概率是多少？p=0，percolates概率為0，p=100，percolates的概率為100。下圖是20x20和100x100格子的概率p的分布：

當N足夠大時, 有一個閾值P*, 使得當p < p*時候，任意的N*N網格，幾乎不能被滲透, 並且當p > p*, 基本能夠被滲透。p*沒有准確的數值解。任務是寫一個計算估計p*的算法。

Percolation data type.

public class Percolation {
   public Percolation(int N)              // create N-by-N grid, with all sites blocked
   public void open(int i, int j)         // open site (row i, column j) if it is not already
   public boolean isOpen(int i, int j)    // is site (row i, column j) open?
   public boolean isFull(int i, int j)    // is site (row i, column j) full?
   public boolean percolates()            // does the system percolate?
}

其中網格的左上角為(1, 1)，右下角為(N, N)。

使用boolean[] matrix來標記網格的打開和關閉。

isOpen表示當前網格是否是打開，直接返回matrix當前位置的值即可。

思考open，isFull，percolates的使用並查集如何實現和復雜度。

open要做什么？思路比較簡單，將一個關閉的網格位置打開，並且和周圍4方向上已經打開的格子連通在一起。

isFull判斷當前點是否能夠滲透到，最簡單的想法，枚舉第一層每個打開的格子，並查集查詢判斷與當前格子是否在一個集合中，如果在那么就是成功的。復雜度枚舉最差N, 再加上每次需要進行一次1個Find()的判斷。

percolates判斷系統是否能夠被滲透，最簡單的想法，枚舉第一層和最后一層打開的格子，判斷是否在一個集合中，如果在那么滲透成功，復雜度N^2，再加上並查集Find()判斷。

最簡單的想法太暴力，復雜度太高，那么我們開始優化

1. 加入虛擬節點，首尾加入兩個虛擬節點，首虛擬節點與第一層所有打開的元素相關聯，尾虛擬節點和最后一層所有打開的節點相關聯，那么對於isFull操作，如果首虛擬節點與當前點在一個集合，那么成功，對於percolates操作如果首尾虛擬節點在一個集合中，那么系統是滲透成功的。這樣所有的這些操作中都不會有枚舉的for循環出現了，但是也會出現一個問題，backwash，因為有一些格子雖然從上面不能滲透到，但和尾虛擬節點連接后，他們也和首虛擬節點在一個集合了。

2. 如何避免backwash？想了個比較方法，又加入了一個並查集，現在兩個並查集，一個只負責維護首虛擬節點，當進行isFull判斷是否只考慮這個並查集，另一個並查集進行首尾虛擬節點的維護，用在percolates的判斷中，這樣backwash的問題也解決了，但是犧牲了空間。

下面是實現代碼：

public class Percolation {
    
    private boolean[] matrix;
    private int row, col;
    private WeightedQuickUnionUF wquUF;
    private WeightedQuickUnionUF wquUFTop;
    private boolean alreadyPercolates;
    
    public Percolation(int N) {
        if (N < 1) throw new IllegalArgumentException("Illeagal Argument");
        wquUF = new WeightedQuickUnionUF(N*N+2);
        wquUFTop = new WeightedQuickUnionUF(N*N+1);
        alreadyPercolates = false;
        row = N;
        col = N;
        matrix = new boolean[N*N+1];
    }
    
    private void validate(int i, int j) {
        if (i < 1 || i > row) 
            throw new IndexOutOfBoundsException("row index i out of bounds");
        if (j < 1 || j > col) 
            throw new IndexOutOfBoundsException("col index j out of bounds");        
    }
    
    public void open(int i, int j) {
        validate(i, j);
        int curIdx = (i-1)*col + j; 
        matrix[curIdx] = true;
        if (i == 1) {
            wquUF.union(curIdx, 0);
            wquUFTop.union(curIdx, 0);
        }
        if (i == row) {
            wquUF.union(curIdx, row*col+1);
        }
        
        int[] dx = {1, -1, 0, 0};
        int[] dy = {0, 0, 1, -1};
        for (int dir = 0; dir < 4; dir++) {
            int posX = i + dx[dir];
            int posY = j + dy[dir];
            if (posX <= row && posX >= 1 
                    && posY <= row && posY >= 1 
                    && isOpen(posX, posY)) {
                wquUF.union(curIdx, (posX-1)*col+posY);
                wquUFTop.union(curIdx, (posX-1)*col+posY);
            }
        }
    }
    
    public boolean isOpen(int i, int j) {
        validate(i, j);
        return matrix[(i-1)*col + j];
    }
    
    public boolean isFull(int i, int j) {
        validate(i, j);
        int curIdx = (i-1)*col+j;
        if (wquUFTop.find(curIdx) == wquUFTop.find(0)) return true;
        return false;
    }
    
    public boolean percolates() {
        if (alreadyPercolates) return true;
        if (wquUF.find(0) == wquUF.find(row*col+1)) {
            alreadyPercolates = true;
            return true;
        } 
        return false;
    }
    
    public static void main(String[] args) {
        Percolation perc = new Percolation(2);
        perc.open(1, 1);
        perc.open(1, 2);
        perc.open(2, 1);
        System.out.println(perc.percolates());
    }
    
}

Monte Carlo simulation. 估計percolation的閾值，初始化時候格子都是關閉的，隨機尋找一個關閉的位置打開，直到系統可以滲透為止，打開的格子比上總格子數就是閾值。

這一部分就是一些數值計算比較簡單，不贅述了。

public class PercolationStats {
    private double[] x;
    private int expTimes;
    
    public PercolationStats(int N, int T) {
        if (N <= 0 || T <= 0) 
            throw new IllegalArgumentException("Illeagal Argument");
        x = new double[T+1];
        expTimes = T;
        for (int i = 1; i <= T; ++i) {
            Percolation perc = new Percolation(N);
            boolean[] isEmptySiteLine = new boolean[N+1];
            int numOfLine = 0;
            while (true) {    
                int posX, posY;
                do {
                    posX = StdRandom.uniform(N)+1;
                    posY = StdRandom.uniform(N)+1;
                } while (perc.isOpen(posX, posY));
                perc.open(posX, posY);
                x[i] += 1;
                if (!isEmptySiteLine[posX]) {
                    isEmptySiteLine[posX] = true;
                    numOfLine++;
                }
                if (numOfLine == N) {
                    if (perc.percolates())
                        break;
                }
            }
            x[i] = x[i] / (double) (N*N);
        }
    }
    
    public double mean() {
        double mu = 0.0;
        for (int i = 1; i <= expTimes; ++i) {
            mu += x[i];
        }
        return mu / (double) expTimes;
    }
    
    public double stddev() {
        if (expTimes == 1)
            return Double.NaN;
        double sigma = 0.0;
        double mu = mean();
        for (int i = 1; i <= expTimes; ++i) {
            sigma += (x[i]-mu)*(x[i]-mu);
        }
        return Math.sqrt(sigma / (double) (expTimes-1));
    }
    
    public double confidenceLo() {
        double mu = mean();
        double sigma = stddev();
        return mu - 1.96*sigma / Math.sqrt(expTimes);
    }
    
    public double confidenceHi() {
        double mu = mean();
        double sigma = stddev();
        return mu + 1.96*sigma / Math.sqrt(expTimes);
    }
    
    public static void main(String[] args) {
        int N = Integer.parseInt(args[0]);
        int T = Integer.parseInt(args[1]);
        PercolationStats percStats = new PercolationStats(N, T);
        StdOut.printf("mean = %f\n", percStats.mean());
        StdOut.printf("stddev = %f\n", percStats.stddev());
        StdOut.printf("95%% confidence interval = %f, %f\n", 
                      percStats.confidenceLo(), percStats.confidenceHi());
        
    }
    
    
}