这次算法作业的任务是通过蒙特卡洛模拟来计算 percolation threshold即渗滤阈值。
渗滤:percolation渗滤是一个由绝缘体和金属随机分布的复杂系统。那么它的金属分布在什么情况下会导致它是一个导体。科学家定义了一个抽象的被称作percolation的过程来为这些情况建模。
模型:这个模型被定义为一个n*n的方格来表示。方格site有两种状态:阻塞block和打开open。在方格中,阻塞用黑色表示,打开用白色表示。当一个方格被称为full时,表示这个方格是打开的并且通过一系列与之相邻(上下左右)的打开的方格与最上层的open的方格连接。当一个系统是percolates时,表示该系统最上层的打开的方格与最底层的打开的方格连接。也就是说存在最底层的方格为full。
图片中黑色为阻塞,白色为打开,蓝色为full。
问题:假设这些格子都是独立的,它们被打开的几率为p。那么系统为渗滤的概率是多少。显而易见,当p为0时,渗滤的概率为0;p为1时,渗滤的概率为1.下图展现了格子打开几率p在20*20与100*100格子下与系统为渗滤的概率的关系。
当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) {
}
}