什么是渗流算法？

什么是渗流(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());
 }
}