coursera算法课 Programming Assignment 1:Percolation

这次算法作业的任务是通过蒙特卡洛模拟来计算 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) {

    }


}