渗滤：percolation渗滤是一个由绝缘体和金属随机分布的复杂系统。那么它的金属分布在什么情况下会导致它是一个导体。科学家定义了一个抽象的被称作percolation的过程来为这些情况建模。
模型：这个模型被定义为一个n*n的方格来表示。方格site有两种状态：阻塞block和打开open。在方格中，阻塞用黑色表示，打开用白色表示。当一个方格被称为full时，表示这个方格是打开的并且通过一系列与之相邻（上下左右）的打开的方格与最上层的open的方格连接。当一个系统是percolates时，表示该系统最上层的打开的方格与最底层的打开的方格连接。也就是说存在最底层的方格为full。
image.png
image.png
图片中黑色为阻塞，白色为打开，蓝色为full。
问题：假设这些格子都是独立的，它们被打开的几率为p。那么系统为渗滤的概率是多少。显而易见，当p为0时，渗滤的概率为0；p为1时，渗滤的概率为1.下图展现了格子打开几率p在2020与100100格子下与系统为渗滤的概率的关系。
image.png
image.png
当n足够大时，存在渗滤阈值p，当p小于p时，系统几乎不渗滤；当p大于p时，系统几乎渗滤。我们的任务是估计出这个渗滤阈值p。
当网格中所有方格为阻塞状态，随机打开一个方格后，判断该系统是否渗滤，如果没有，则继续打开，直到系统渗滤。那么p就为打开的方格数除以所有的方格数。进行大量多次实验就可以估计p的值。
蒙特卡洛模拟：用来估计渗滤阈值。设定网格中所有方格为阻塞状态，随机打开一个方格后，判断该系统是否渗滤，如果没有，则继续打开，直到系统渗滤。那么p就为打开的方格数除以所有的方格数。进行大量多次实验就可以估计p的值。
这一段主要是一些数值计算，包括平均值，标准差，低95%置信区间的端点，高95%置信区间的端点。可以使用课程提供的StdStats来计算。
代码
Percalation.java
import edu.princeton.cs.algs4.StdOut;
import edu.princeton.cs.algs4.WeightedQuickUnionUF;
public class Percolation {
// n*n grid
private int n ;
//status of each site
private boolean[] eachStatus ;
//number of open site
private int openNum;
//UF alg with virtual site
private WeightedQuickUnionUF uf1;
//UF alg with only top site due to backwash
private WeightedQuickUnionUF uf2;
// create n-by-n grid, with all sites blocked
public Percolation(int n){
if(n<=0) throw new IllegalArgumentException("illegal value of n!");
this.n = n;
//plus two virtual site
eachStatus = new boolean[(n+1)*(n+2)];
//all sites are blocked
for(int i=0; i< eachStatus.length;i++) eachStatus[i]=false;
// grid with top site and bottom site
uf1 = new WeightedQuickUnionUF((n+1)*(n+2));
// grid with only bottom site
uf2 = new WeightedQuickUnionUF((n+1)*(n+2));
}
//translate 2 dimentional coordinate to 1 dimentional pattern
private int oneDimention(int row, int col){
return row*(n+1)+col;
}
// open site (row, col) if it is not open already
public void open(int row, int col){
if(row<1 || row>n || col<1 || col>n) throw new IllegalArgumentException("illegal row or col!");
if(!isOpen(row,col)) {
eachStatus[oneDimention(row,col)]=true;
openNum++;
int temp1 = oneDimention(row,col);
//if neighbor could be connected?
if(row==1){
uf1.union(0,temp1);
uf2.union(0,temp1);
}
if(row==n) {
uf1.union((n+1)*(n+1),temp1);
}
int[] dx = {1,-1,0,0};
int[] dy = {0,0,1,-1};
for(int i=0; i<4;i++){
int tempX = row+dx[i];
int tempY = col+dy[i];
int temp2 = oneDimention(tempX,tempY);
if (eachStatus[temp2]){
uf1.union(temp1,temp2);
uf2.union(temp1,temp2);
}
}
}
}
//is site (row, col) open?
public boolean isOpen(int row, int col){
if(row<1 || row>n || col<1 || col>n) throw new IllegalArgumentException("illegal row or col!");
return eachStatus[oneDimention(row,col)];
}
//is site (row, col) full?
public boolean isFull(int row, int col){
if(row<1 || row>n || col<1 || col>n) throw new IllegalArgumentException("illegal row or col!");
return uf2.connected(0,oneDimention(row,col));
}
// number of open sites
public int numberOfOpenSites() {
return openNum;
}
// does the system percolate?
public boolean percolates() {
return uf1.connected(0,(n+1)*(n+1));
}
// test client (optional)
public static void main(String[] args) {
Percolation p = new Percolation(3);
p.open(1, 2);
p.open(2, 2);
p.open(3, 2);
StdOut.println(p.isOpen(1,1));
StdOut.println(p.percolates());
}
}
PercolationStats.java
import edu.princeton.cs.algs4.StdOut;
import edu.princeton.cs.algs4.StdRandom;
import edu.princeton.cs.algs4.StdStats;
public class PercolationStats {
//trial times
private int trialNum;
//threshold P
private double[] preP;
public PercolationStats(int n,int trialNum) {
if (n<1 || trialNum<1) throw new IllegalArgumentException("Illegal n or trialNum,please check");
this.trialNum = trialNum;
preP = new double[trialNum];
for(int i=0;i<trialNum;i++) {
Percolation p = new Percolation(n);
while(!p.percolates()) {
int row = StdRandom.uniform(n)+1;
int col = StdRandom.uniform(n)+1;
p.open(row,col);
if (p.percolates()) break;
}
preP[i] = (double)p.numberOfOpenSites()/(n*n);
}
}
// mean of p
public double mean() {
return StdStats.mean(preP);
}
// standard deviation
public double stddev() {
return StdStats.stddev(preP);
}
//confidence interval:low
public double confidenceLo() {
return mean()-1.96*stddev()/Math.sqrt(trialNum);
}
//confidence interval:high
public double confidenceHi() {
return mean()+1.96*stddev()/Math.sqrt(trialNum);
}
public static void main(String[] args) {
int n =10,trialNum=1000;
PercolationStats ps = new PercolationStats(n,trialNum);
StdOut.println("grid size :" + n+"*"+n);
StdOut.println("trial times :" + trialNum);
StdOut.println("mean of p :"+ ps.mean());
StdOut.println("standard deviation :"+ps.stddev());
StdOut.println("confidence interval low :"+ps.confidenceLo());
StdOut.println("confidence interval high :"+ps.confidenceHi());
}
}