這次演算法作業的任務是通過蒙特卡洛模擬來計算 percolation threshold即滲濾閾值。
滲濾：percolation滲濾是一個由絕緣體和金屬隨機分佈的複雜系統。那麼它的金屬分佈在什麼情況下會導致它是一個導體。科學家定義了一個抽象的被稱作percolation的過程來為這些情況建模。
模型：這個模型被定義為一個n*n的方格來表示。方格site有兩種狀態：阻塞block和開啟open。在方格中，阻塞用黑色表示，開啟用白色表示。當一個方格被稱為full時，表示這個方格是開啟的並且通過一系列與之相鄰（上下左右）的開啟的方格與最上層的open的方格連線。當一個系統是percolates時，表示該系統最上層的開啟的方格與最底層的開啟的方格連線。也就是說存在最底層的方格為full。

當n足夠大時，存在滲濾閾值p*，當p小於p*時，系統幾乎不滲濾；當p大於p*時，系統幾乎滲濾。我們的任務是估計出這個滲濾閾值p*。
當網格中所有方格為阻塞狀態，隨機開啟一個方格後，判斷該系統是否滲濾，如果沒有，則繼續開啟，直到系統滲濾。那麼p*就為開啟的方格數除以所有的方格數。進行大量多次實驗就可以估計p*的值。
下面是Percolation的api：

public class Percolation {
   public Percolation(int n)                // create n-by-n grid, with all sites blocked
   public    void open(int row, int col)    // open site (row, col) if it is not open already
   public boolean isOpen(int row, int col)  // is site (row, col) open?
   public boolean isFull
 
(int row, int col)  // is site (row, col) full?
   public     int numberOfOpenSites()       // number of open sites
   public boolean percolates()              // does the system percolate?

   public static void main(String[] args)   // test client (optional)
}

實現程式碼：

import edu.princeton.cs.algs4.WeightedQuickUnionUF;

public class Percolation {
    //虛擬頂點在grid[]中的index
    private final static int virtualTop = 0;
    //網格的邊長
    private int n;          
    //為加上虛擬頂點和虛擬底點和網格所有方格的數量。
    private int gridLength;
    //為虛擬底點。通常為n*n+1
    private int virtualBottom;
    //用於標記方格的開啟和關閉，grid[0]為虛擬頂點，grid[virtualBottom]為虛擬底點
    private boolean[] grid;
    //是否為滲濾，滲濾為true
    private boolean isPercolated;
    //並查集，包含虛擬頂點，虛擬底點和網格所有方格。
    private WeightedQuickUnionUF unionInstance;
    //為了防止backwash問題，創立的並查集，不包含虛擬底點
    private WeightedQuickUnionUF backWash;
    // create n-by-n grid, with all sites blocked 
    public Percolation(int n) {
        if (n < 1) {
            throw new IllegalArgumentException("Illegal Argument"); 
        }       
        this.n = n;
        gridLength = n * n + 2;
        isPercolated = false;
        virtualBottom = n * n + 1;
        unionInstance = new WeightedQuickUnionUF(gridLength);
        backWash = new WeightedQuickUnionUF(gridLength - 1);
        grid = new boolean[gridLength];
    }
    /*  首先判斷row，col是否越界,沒越界則將該方格開啟。當方格為第一行，將其與虛擬頂點連線
     * 當為最底部，與虛擬底點連線。然後看方格的上下左右，若為開啟的，則與方格連線。
     * */
    // open site (row, col) if it is not open already   
    public void open(int row, int col) {
        isValidBounds(row, col);
        if (isOpen(row, col))
            return;
        int gridId = xyTo1D(row, col);
        grid[xyTo1D(row, col)] = true;
        if (1 == row) {
            unionInstance.union(virtualTop, gridId);
            backWash.union(virtualTop, gridId);
        }
        if (n == row) {
            unionInstance.union(virtualBottom, gridId);
        }
        int[] dx = {-1, 0, 0, 1};
        int[] dy = {0, -1, 1, 0};
        for (int i = 0; i < 4; i++) {
            int posX = row + dx[i];
            int posY = col + dy[i];
            if (isPosValid(posX, posY) 
                    && isOpen(posX, posY)) {
                int posId = xyTo1D(posX, posY);
                unionInstance.union(gridId, posId);
                backWash.union(gridId, posId);
            }   
        }


    }

    //將方格的行列位置化為方格在grid中的位置
    private int xyTo1D(int row, int col) {
        int i = (row - 1)  * n + col;
        return i;
    }

    //判斷方格行列是否越界
    private boolean isPosValid(int row, int col) {
        if (row >= 1 && row <= n && col >= 1 && col <= n) {
            return true;
        }
        return false;
    }

    //判斷是否row，col越界
    private void isValidBounds(int row, int col) {
        if (row < 1 || row > n)
            throw new IndexOutOfBoundsException("Row index out of bounds");
        if (col < 1 || col > n)
            throw new IndexOutOfBoundsException("column index out of bounds");
    }   

    // is site (row, col) open?   
    public boolean isOpen(int row, int col) {
        isValidBounds(row, col);
        return (grid[xyTo1D(row, col)] == false ? false : true);
    }
    // is site (row, col) full?   
    public boolean isFull(int row, int col) {
        isValidBounds(row, col);
        return backWash.connected(virtualTop, xyTo1D(row, col));
    }
    // number of open sites   
    public int numberOfOpenSites() {
        int openNum = 0;
        for (int i = 1; i < virtualBottom; i++) {
            if (grid[i])
                openNum++;
        }
        return openNum;
    }
    // does the system percolate?   
    public boolean percolates() {
        if (isPercolated) 
            return true;
        if (unionInstance.connected(virtualTop, virtualBottom)) {
            isPercolated = true;
            return true;
        }
        return false;   
    }
    // test client (optional)   
    public static void main(String[] args) {
        Percolation perc = new Percolation(3);
        perc.open(1, 1);
        perc.open(1, 2);
        perc.open(2, 2);
        perc.open(2, 3);
        perc.open(3, 1);
        perc.open(3, 3);
        System.out.println(perc.percolates());
    }

}

兩個問題：
1.判斷系統是否滲濾就是看頂層格子是否與底層的格子連線。如果單純的檢視是否連線，那麼只能遍歷頂層的格子和底層的格子，看它們的連線情況。所以可以設立一個虛擬頂點和一個虛擬底點，虛擬頂點與最上層的格子連線。虛擬底點與最下層的格子連線。這樣只要看虛擬頂點和虛擬底點是否連線，就可以判斷系統是否滲濾了。
2.當使用虛擬頂點和虛擬底點時，會出現backwash的情況。有一些格子只與最底層連線，而沒有連線到頂層格子，但是在並查集中確實顯示能與頂層格子連線，從而被標記為full狀態
如下圖所示。

那麼要解決這個情況，需要再建立一個並查集，只包含虛擬頂點和網格所有方格，不包含虛擬底點。這樣在isFull（）可以使用它來判斷就不會出現backwash問題了。

我用boolean[] grid來標記方格的開啟和關閉。
isOpen()就是檢視方格在grid[]是否為true。
isFull()就是看方格是否與虛擬頂點連線。
percolates（）通過虛擬頂點和虛擬底點是否連線來判斷系統是否滲濾。
open（）為開啟方格。思想很簡單，當開啟一個方格，就將該方格上下左右的開啟的方格與它連線。

蒙特卡洛模擬：用來估計滲濾閾值。設定網格中所有方格為阻塞狀態，隨機開啟一個方格後，判斷該系統是否滲濾，如果沒有，則繼續開啟，直到系統滲濾。那麼p*就為開啟的方格數除以所有的方格數。進行大量多次實驗就可以估計p*的值。
這一段主要是一些數值計算，包括平均值，標準差，低95%置信區間的端點，高95%置信區間的端點。可以使用課程提供的StdStats來計算。
實現程式碼：

import edu.princeton.cs.algs4.StdOut;
import edu.princeton.cs.algs4.StdRandom;

public class PercolationStats {
    private double[] x;
    private int expTimes;

    public PercolationStats(int n, int trials) {
        if (n <= 0 || trials <= 0) 
            throw new IllegalArgumentException("Illeagal Argument");
        x = new double[trials+1];
        expTimes = trials;
        for (int i = 1; i <= trials; ++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) {

    }


}