什么是滲流(Percolation)？
滲流是物理中用來理解物理現象的一種算法，是很多物理系統的模型。
我們假設一個NxN的方形網格，其中的小網格我們叫做位(site)，每個site開放(open)的概率是 1−p1−p ，site閉合(blocked)的概率是 p−1p−1 ，如果頂端(top)和底端(bottom)的site連通且能連通出一條連通路徑的話，我們稱這個系統是滲流(percolates)的，如下圖所示。
實際問題
專業術語往往都是晦澀難懂的，我們來看一個比較切合實際的滲流問題：
將一個不透水的均質方塊分割為NxN，最上方為水源，隨機打開方塊中任意格子，重復此項操作多次，產生一條路徑使水能穿過這個方塊到達最下方。
上面兩張圖分別展示了滲流(percolates)和不滲流(not percolates)兩種情況，和我們上面所說的特質相同，使用NxN的格點，每個格點site擁有open & blocked兩種狀態。如果最上面與最下面由open site路徑連通，那么它就是滲流的。
問題分析
這個問題和之前在學習的并查集的動態連通性問題(Dynamic Connectivity)相聯系，因為當open了一個site之后，這個site需要和周邊的(上下左右四個方位)的site連接，類似于我們在并查集中做查詢合并操作。這里我們可以用并查集相關的算法來解決。
那么隨之而來的問題就是，我們如何確定最頂層和最底層是相連通的呢？這里最頂層和最底層相連通每次都是最頂層的site和最底層的site，可是這里的site是random的，是具有不確定性的。
在Coursera上的《Algorithms 4th Editon》由普林斯頓大學錄制的課程中，提到了這個問題的解題思路，我們可以在最頂層與最底層再加一行，但是那一行只擁有一個site，我們如果要確定這個系統是否是滲流的通過這兩個site即可確定。
且每個site的open & blocked的狀態是互相獨立的，那么我們需要一個單獨的容器來維護每個site的open & blocked狀態。
分析完問題，直接看代碼來理解會理解的更快。
Code
下面的代碼都可以在我的Github倉庫找到。
package net.sunzhenyu.miscellaneous.dsa.other.percolation;
import edu.princeton.cs.algs4.WeightedQuickUnionUF;
/**
* Percolation Model
*
* @author <a href="mailto:sunzhenyucn@gmail.com">SunZhenyu</a>
* @version 0.0.1
* @date 2018/4/19
* @since 1.8
*/
public class Percolation {
/**
* Initial number
*/
private final int n;
/**
* Open & Blocked status array
*/
private final boolean[] open;
/**
* The n * n grid
*/
private final WeightedQuickUnionUF uf, uf2;
/**
* Open site count cache
*/
private int count;
/**
* Constructor
*
* @param n create n * n grid
* @throws IllegalArgumentException if the initial num not between 1 and n
*/
public Percolation(int n) {
if (n <= 0) {
throw new IllegalArgumentException("The initial num must between 1 and " + n);
}
this.n = n;
open = new boolean[n * n + 2];
// Because we add 2 sites, so we need plus 2
uf = new WeightedQuickUnionUF(n * n + 2);
// Why use this?
uf2 = new WeightedQuickUnionUF(n * n + 1);
for (int i = 0; i < this.n; i++) {
// The initial status is blocked
open[i] = false;
}
}
/**
* Verify coordinate for avoid {@link IllegalArgumentException}
*
* @param row coordinate
* @param col coordinate
* @throws IllegalArgumentException if coordinate not between one and {@code n}
*/
private void outBoundsCheck(int row, int col) {
if (row <= 0 || row > n || col <= 0 || col > n) {
throw new IllegalArgumentException(
"The coordinate must be [1, n], please check your coordinate");
}
}
/**
* Transfer coordinate to array index
*
* @param row coordinate
* @param col coordinate
* @return array index of coordinate
*/
private int index(int row, int col) {
outBoundsCheck(row, col);
return (row - 1) * n + col;
}
/**
* Check specified coordinate's site is open
*
* @param row coordinate
* @param col coordinate
* @return if open return true
*/
public boolean isOpen(int row, int col) {
outBoundsCheck(row, col);
return open[index(row, col)];
}
/**
* Check specified coordinate's site is full
*
* @param row coordinate
* @param col coordinate
* @return if full return true
*/
public boolean isFull(int row, int col) {
outBoundsCheck(row, col);
if (!isOpen(row, col)) {
return false;
}
return uf2.connected(index(row, col), 0);
}
/**
* Return this percolate model's open site count
*
* @return this percolate model's open site count
*/
public int numberOfOpenSites() {
return count;
}
/**
* Return this percolate model's percolates status
*
* @return this percolate model's percolates status
*/
public boolean percolates() {
return uf.connected(0, n * n + 1);
}
/**
* Open specified coordinate's site and neighborhood's site if they are not open already
*
* @param row coordinate
* @param col coordinate
*/
public void open(int row, int col) {
outBoundsCheck(row, col);
if (isOpen(row, col)) {
return;
}
//Open this site
open[index(row, col)] = true;
// Top
if (row == 1) {
open[0] = true;
uf.union(index(row, col), 0);
uf2.union(index(row, col), 0);
}
// Bottom
if (row == n) {
open[n * n + 1] = true;
uf.union(index(row, col), n * n + 1);
}
// Up
if (row > 1 && isOpen(row - 1, col)) {
uf.union(index(row, col), index(row - 1, col));
uf2.union(index(row, col), index(row - 1, col));
}
// Down
if (row < n && isOpen(row + 1, col)) {
uf.union(index(row, col), index(row + 1, col));
uf2.union(index(row, col), index(row + 1, col));
}
// Left
if (col > 1 && isOpen(row, col - 1)) {
uf.union(index(row, col), index(row, col - 1));
uf2.union(index(row, col), index(row, col - 1));
}
// Right
if (col < n && isOpen(row, col + 1)) {
uf.union(index(row, col), index(row, col + 1));
uf2.union(index(row, col), index(row, col + 1));
}
// Increment open count
count++;
}
public static void main(String[] args) {
int n = 3;
Percolation p = new Percolation(n);
// while (!p.percolates()) {
// int i = StdRandom.uniform(1, n + 1);
// int j = StdRandom.uniform(1, n + 1);
// if (!p.isOpen(i, j)) {
// p.open(i, j);
// System.out.println("open (" + i + "," + j + "), isFull::" + p.isFull(i, j));
// }
// }
p.open(2, 1);
System.out.println(p.isOpen(1, 1));
System.out.println(p.isOpen(1, 1));
System.out.println(p.isFull(1, 1));
}
}
下面我們來說一下為什么需要兩個WeightedQuickUnionUF。
Backwash
我們使用虛擬頂點與虛擬底點的方式來解決對系統是否滲流的判斷問題，從而引出了一個新的問題，那就是Backwash。
為什么會出現這樣的情況呢？因為有的site只與最底層連通，而沒有連通至最高層的site，但在并查集的數據結構中確實顯示能夠與頂層的site連通，從而被標記成了full的狀態。
那么要解決這個問題，就需要另一個WeightedQuickUnionUF，這個不包含虛擬底點，這樣在isFull()方法中檢查的時候就不會出現這種問題了。
蒙特卡洛模擬
蒙特卡羅方法（英語：Monte Carlo method），也稱統計模擬方法，是1940年代中期由于科學技術的發展和電子計算機的發明，而提出的一種以概率統計理論為指導的數值計算方法。是指使用隨機數（或更常見的偽隨機數）來解決很多計算問題的方法。摘自Wikipedia 蒙特卡洛方法 詞條
通俗理解就是通過大量地隨機樣本，進而得到目標問題的所要計算的值。
在本題中，用來估計滲濾閾值。設定網格中所有方格為阻塞狀態，隨機打開一個方格后，判斷該系統是否滲濾，如果沒有，則繼續打開，直到系統滲濾。那么$p$就為打開的方格數除以所有的方格數。進行大量多次實驗就可以估計$p$的值。
package net.sunzhenyu.miscellaneous.dsa.other.percolation;
import edu.princeton.cs.algs4.StdOut;
import edu.princeton.cs.algs4.StdRandom;
import edu.princeton.cs.algs4.StdStats;
/**
* Percolation Monte Monte Carlo simulation
*
* @author <a href="mailto:sunzhenyucn@gmail.com">SunZhenyu</a>
* @version 0.0.1
* @date 2018/4/19
* @since 1.8
*/
public class PercolationStats {
private static final double NUM = 1.96;
private int t;
private double[] fraction;
public PercolationStats(int n, int t) { // perform t independent experiments on an n-by-n grid
if (n <= 0 || t <= 0) {
throw new IllegalArgumentException("n and t must be bigger than 0");
}
this.t = t;
fraction = new double[t];
for (int count = 0; count < t; count++) {
Percolation pr = new Percolation(n);
int openedSites = 0;
while (!pr.percolates()) {
int i = StdRandom.uniform(1, n + 1);
int j = StdRandom.uniform(1, n + 1);
if (!pr.isOpen(i, j)) {
pr.open(i, j);
openedSites++;
}
}
fraction[count] = (double) openedSites / (n * n);
}
}
public double mean() { // sample mean of percolation threshold
return StdStats.mean(fraction);
}
public double stddev() { // sample standard deviation of percolation threshold
return StdStats.stddev(fraction);
}
public double confidenceLo() { // low endpoint of 95% confidence interval
return mean() - NUM * stddev() / Math.sqrt(t);
}
public double confidenceHi() { // high endpoint of 95% confidence interval
return mean() + NUM * stddev() / Math.sqrt(t);
}
public static void main(String[] args) // test client (described below)
{
int n = Integer.parseInt(args[0]);
int t = Integer.parseInt(args[1]);
PercolationStats ps = new PercolationStats(n, t);
StdOut.println("mean = " + ps.mean());
StdOut.println("stddev = " + ps.stddev());
StdOut.println("95% confidence interval = " + ps.confidenceLo() + ", " + ps.confidenceHi());
}
}
