用java调用Union-Find-API解决蒙特卡洛模拟问题

渗透问题和蒙特卡洛模拟问题

Percolation渗透问题

问题介绍

引发：给定由随机分布的绝缘材料和金属材料构成的组合系统：金属材料占多大比例才能使组合系统成为电导体？ 给定一个表面有水的多孔渗水地形（或下面有油），水将在什么条件下能够通过底部排出（或油渗透到表面）？科学家们已经定义了一个称为渗透（percolation）的抽象过程来模拟这种情况
模型：我们使用N×N网格点来模型一个渗透系统，每个格点或是open格点或是blocked格点，一个full site是一个open格点，它可以通过一连串的邻近（左，右，上，下）open格点连通到顶行的一个open格点，如果在底行中有一个full site格点，则称系统是渗透的（对于绝缘/金属材料的例子，open格点对应于金属材料，渗透系统有一条从顶行到底行的金属路径，且full sites格点导电。对于多孔物质示例，open格点对应于空格，水可能流过，从而渗透系统使水充满open格点，自顶向下流动。）。

问题： 在一个著名的科学问题中，研究人员对以下问题感兴趣：如果将格点以空置概率p独立地设置为open格点（因此以概率1-p被设置为blocked格点），系统渗透的概率是多少？ 当p = 0时，系统不会渗出; 当p=1时，系统渗透。 下图显示了20×20随机网格（左）和100×100随机网格（右）的格点空置概率p与渗滤概率。

当N足够大时，存在阈值p*，使得当p <p*，随机N N网格几乎不会渗透，并且当p> p时，随机N N网格几乎总是渗透。 尚未得出用于确定渗滤阈值p的数学解。你的任务是编写一个计算机程序来估计p*。

public class Percolation {
   public Percolation(int N)           // create N-by-N grid, with all sites blocked
   public void open(int i, int j)         // open site (row i, column j) if it is not already
   public boolean isOpen(int i, int j)    // is site (row i, column j) open?
   public boolean isFull(int i, int j)     // is site (row i, column j) full?
   public boolean percolates()          // does the system percolate?
   public static void main(String[] args)  // test client, optional
}

约定行i列j下标在1和N之间，其中(1, 1)为左上格点位置：如果open(), isOpen(), or isFull()不在这个规定的范围，则抛出IndexOutOfBoundsException例外。如果N ≤ 0，构造函数应该抛出IllegalArgumentException例外。构造函数应该与N2成正比。所有方法应该为常量时间加上常量次调用合并-查找方法union(), find(), connected(), and count()。

Union-Find API

type UF interface {
	NewUF(int)             //初始化整数个触点
	Union(int, int)          //在两个触点之间增加一条连接
	Find(int) int            //查找该触点所在分量的标识符
	Connected(int, int) bool //判断两个触点是否相连
	Count() int              //连通分量的数量
}

我在这里就不给大家详细讲述什么叫做Union-Find了，就是我们熟知的并查集，感兴趣的朋友们可以自行去了解了解什么叫并查集，还有带权并查集，我后面也会大致提一提戴权的并查集。

问题分析与解决

拿到这个问题，我们的第一反应就是建立一个N*N的网状模型，其次，我们应该可以想到，建立一个二维的布尔数组，其中true表示open,false 表示blocked，但是我们查阅并查集的接口会发现，Union-Find是建立在一维数组的基础里进行维护的，这时我们就要想到，数组不管是几维的，在内存中它都是连续存储的，所以我们就很快可以想到，把这个问题映射到一维数组会更好。于是我们在写Percolation这个类的时候，能自然而然的想到几个参数：
1.逻辑上二维数组的长度
2.真实存储时一维数组的长度
3.一个Union-Find变量用来维护联通关系
4.一个boolean类型的一维数组

import edu.princeton.cs.algs4.StdIn;
import edu.princeton.cs.algs4.StdOut;
import edu.princeton.cs.algs4.UF;
public class Percolation{
		private int gridLength;//表示逻辑上二维数组的长度
		private int virtualLength;//表示虚拟成一维数组时的一维数组长度
		private UF uf;//Union-Find变量
		private boolean[] grid;//存放元素的一维数组
}

写完了这几个属性，接下来肯定是要实现要求我们的各种方法，而对于面向对象的语言，一个好的习惯就是先写他的构造方法
对于构造有如下几点要注意
1.参数应该是NN网格的N
2.如果N不合法，要抛出一个异常
3.程序的重点是在于实现一个连通分量，怎么实现呢，这里我是让boolean类型数组的长度为NxN+2,为什么这么设立呢?第一，是因为我们问题要求，我们NN网格里，第一个元素下标不是（0，0）,而是（1,1）。第二，是我们可以让grid[0]表示这个网状图的虚拟顶部,让最后一项grid[N*N+1]做这个网状图的虚拟底部，这个大家脑海里必须得转过来这个思维，然后我们把虚拟顶部和虚拟底部全都设置为true,对于一个系统而言，判断它是否渗透，只需要判断顶部和底部是否处于一个连通分量

public Percolation(int N) {
		if (N <= 0) {
			throw new IllegalArgumentException("长度不能小于1");
		}
		gridLength = N;// 网络长度为N
		virtualLength = N * N;// 即将二维数组的问题映射到一维数组上
		grid = new boolean[virtualLength+2];
		grid[0]=true;//上部分也是表示联通的
		grid[virtualLength+1]=true;//这个就表示联通出去的
		uf = new UF(virtualLength + 2);// 0是虚拟顶部，virtualTop+1是虚拟底部，中间的N*N代表着网络图的内容
	}

接下来，我们就可以实现别的函数了，但是因为题目要求我们要抛出每一次位置不合法导致的异常，所以我们可以把位置非法的异常封装在一个方法里面

private void vaildate(int row, int col) {
		if (row <= 0 || row > gridLength)
			throw new IndexOutOfBoundsException("行索引越界");
		if (col <= 0 || col > gridLength)
			throw new IndexOutOfBoundsException("列索引越界");
	}

紧接着，这一步也相当重要，就是我们要找出二维数组和一维数组之间映射关系，学习过数据结构的朋友们应该都知道，对于二维数组向一维数组映射而言，公式应该是
DDLocation（i，j）=DLocation(i * n+j);
但是题目要求我们下标从（0，0）开始，所以应该改写为
DDLocation（i，j）=DLocation(（i-1） * n+（j-1）);
但是这还没结束，因为在物理的存储空间中，我们的第一位用来存放我们的虚拟顶部，所以我们逻辑上的第一个元素应该在物理上的第二位，因此映射公式应该为
DDLocation（i，j）=DLocation(（i-1） * n+j);

private int getsite(int row,int col){
	return (i-1)*gridLength+col;
}

接下来我们就要写题目要求的所有方法，第一个是open方法，open方法的功能是让一个Blocked的点成为Open的点，我们要注意，这不仅仅是一个简单的赋值，而是要判断该点的上下左右是否也是open是的话，我们要将其放到同一个连通分量里，而对于第一行来说，他的上方就是我们的虚拟顶部，对于最后一行，他的下方就是我们的虚拟底部

（PS:实在懒得用电脑作图了，各位别骂我字丑~~~~）

public void open(int row,int col){
	//第一步先检查是否在表格中，位置是否合法
	vaildate(row,col);
	//如果有效就获取当前的位置
	int site=getsite(row,col);
	//按照流程图上所述，先判断
	if(grid[site]){
		return;
	}
	grid[site]=true;
	//接下来看他的前后左右，先把第一行和最后一行这两个比较特殊的拿出来
	if(row==1)
	{
		uf.union(site,0);//和虚拟顶部相连
	}
	if(row==gridLength)
	{
		uf.union(site,virtual+1);//与虚拟底部相连
	}
	int other;//接着判断其他的各点
	if (row > 1)// up
		{
			other = geisite(row - 1, col);
			if (grid[other]) {
				uf.union(site, other);
			}
		}
		if (row < gridLength) {//down
			other = getsite(row + 1, col);
			if (grid[other]) {
				uf.union(site, other);
			}
		}
		if(col>1) {//left
			other=getsite(row,col-1);
			if(grid[other]) {
				uf.union(site, other);
			}
		}
		if(col<gridLength) {//right
			other=getsite(row,col+1);
			if(grid[other]) {
				uf.union(site, other);
			}					
		}
}

接下来就是题目要求我们的Isopen()方法和isFul(),ispercolate()方法，对于isopen我们只需要返回当前结点的值就可以了，对于isfull我们只需要返回它是否与虚拟顶部有链接便是,对于ispercolate()我们只需要判断虚拟底部和虚拟顶部是否在一个来连通分量上就是了

public boolean isOpen(int row, int col) {
		vaildate(row,col);
		return grid[getsite(row,col)];
	}
public boolean isFull(int row, int col) {
		vaildate(row,col);
		return uf.connected(0,geisite(row,col));//判断是否和顶部存在于一个连通分量
	}	
public boolean ifPrecolate()
{
	return uf.connected(0,virtual+1);
}

这样，我们一个Percolation渗透类就做好啦，我们用一个简单的小例子测试下

public static void main(String[] args) {
		Percolation percolation = new Percolation(3);
        percolation.open(2, 2);
        percolation.open(3, 2);
        for (int i = 1; i <= 3; i++) {
            for(int j = 1; j <= 3; j++) {
                StdOut.println("" + percolation.isFull(i, j) + " " + percolation.isOpen(i, j));
            }
        }
        StdOut.println("percolation is " + percolation.percolates());
	}

结果如图，表示没有问题

蒙特卡洛模拟问题

题目描述和要求

蒙特卡洛模拟（Monte Carlo simulation）：要估计渗透阈值，考虑以下计算实验：
• 初始化所有格点为blocked。
• 重复以下操作直到系统渗出：

1. 在所有blocked的格点之间随机均匀选择一个格点 (row i, column j)。
2. 设置这个格点(row i, column j)为open格点。
   • open格点的比例提供了系统渗透时渗透阈值的一个估计。
   例如，如果在20×20的网格中，根据以下快照的open格点数，那么对渗滤阈值的估计是204/400 = 0.51，因为当第204个格点被open时系统渗透。
   
   通过重复该计算实验T次并对结果求平均值，我们获得了更准确的渗滤阈值估计。 令xt是第t次计算实验中open格点所占比例。 样本均值提供渗滤阈值的一个估计值； 样本标准差测量阈值的灵敏性。
   
   假设T足够大（例如至少30），以下为渗滤阈值提供95％置信区间
   

我们创建数据类型PercolationStats来执行一系列计算实验，包含以下API。

public class PercolationStats {
   public PercolationStats(int N, int T)    // perform T independent computational experiments on an N-by-N grid
   public double mean()                // sample mean of percolation threshold
   public double stddev()               // sample standard deviation of percolation threshold
   public double confidenceLo()         // returns lower bound of the 95% confidence interval
   public double confidenceHi()         // returns upper bound of the 95% confidence interval
   public static void main(String[] args)    // test client, described below
}

在N ≤ 0或T ≤ 0时，构造函数应该抛出java.lang.IllegalArgumentException例外。
此外，还包括一个main( )方法，它取两个命令行参数N和T，在N×N网格上进行T次独立的计算实验（上面讨论），并打印出均值、标准差和95％ 渗透阈值的置信区间。 使用标准库中的标准随机数生成随机数； 使用标准统计库来计算样本均值和标准差。

问题分析和实现

对于这个问题，我们首先得有数学的一点知识，这个类他想要实现的功能就是，对于一个NN的网格而言，重复对他进行T次蒙特卡洛模拟，然后计算出它的标准差，均值，上下置信区间
所以问题的实现要比第一个渗透实验简单多了，首先还是和分析第一个问题一样，我们先想明白他都需要哪些参数
1.NN的网格的长度；
2.重复的次数；
3.实验的均值，方差，上下置信区间
4，用一个数组存放每次单独实验的数据,数组长度为重复的次数

public class PercolationStats {
	private int gridLength;
	private double[] opreationtimes;
	private double ResultMean;//均值u
	private double ResultStddev;//标准差
	private double ResultConfidenceHi;//置信区间上界
	private double ResultConfidenceLo;//置信区间下界
	}

接下俩，我们首先应该想到的是完成构造方法

public PercolationStats(int N,int T) {//在N*N的网格里做T次独立实验
		if(N<=0) {
			throw new IllegalArgumentException("长度不能小于1");
		}
		if(T<=0) {
			throw new IllegalArgumentException("追踪次数不能小于1");
		}
		gridLength = N;
		if(gridLength==1)
		{
			ResultMean=1;
			ResultStddev=Double.NaN;
			ResultConfidenceHi=Double.NaN;
			ResultConfidenceLo=Double.NaN;
		}
		else{
		}
}

我们发现，我们写到else,按理说就应该给存放每次实验数据的数组赋值了，但是很我们还没有一个可以获取单个实验数据的方法，所以我们先停在这里，写单个实验的步骤

private double oneTrials(int length) {//一次测试得到的结果
		int openedCount=0;//计数，看到第几个open时系统渗透
		Percolation percolation=new Percolation(length);
		while(!percolation.percolates()) {//当系统渗透时结束循环
			int row=StdRandom.uniform(length)+1;
			int col=StdRandom.uniform(length)+1;
			if(!percolation.isOpen(row, col)) {
				percolation.open(row, col);
				openedCount++;
			}
		}
		return (double)openedCount/(length*length);
	}

这个逻辑就很简单了，如果刚才大家了解了Percolation渗透的构造，很容易明白，我就不多讲了
既然得到了一次的数据，那我们就接着完善我们构造函数中没写完的else部分

else{
	TrialsResult=new double[T];
			for(int i =0;i<T;i++)
			{
				operationtimes[i]=oneTrials(gridLength);
			}
			ResultMean=StdStats.mean(TrialsResult);
			ResultStddev = StdStats.stddev(TrialsResult);
			double diff=(1.96*ResultStddev)/Math.sqrt(T);
			ResultConfidenceLo = ResultMean - diff;
	        ResultConfidenceHi = ResultMean + diff;
}

里面的数学逻辑我就不赘述了，就是简单的一个置信区间计算，知道次数和每次的概率P是很容易计算的
接着我们按照题目要求把剩余的几个方法补全

public double mean() {//渗透阈值的样本平均值
		return 	ResultMean;
	}
	public double stddev() {//
		return ResultStddev;
	}
	public double ConfidenceLo() {
		return ResultConfidenceLo;
	}
	public double ConfidenceHi() {
		return ResultConfidenceHi;
	}

最后我们再测试一下

public static void main(String[] args) {
		int length=200;
		int trials=100;
		PercolationStats percolations = new PercolationStats(length, trials);
		StdOut.println("mean                    = " + percolations.mean());
		StdOut.println("stddev                  = " + percolations.stddev());
		StdOut.println("95% confidence interval = "
                + percolations.ConfidenceLo() + ", "
                + percolations.ConfidenceHi());
	}

第一次结果是
第二次结果是

第三次结果是

我们和实验题目给我们的几个标准输入后的数字作一比较
第一个

第二个

第三个

得出结论我们的程序设计没问题，可以获得正确估计~

结尾

我这次只是用了个Union-fInd方法，其实这并不是最好的方法，最好的方法是带权的WeightedQuickUnionUF，这个要更加强大，但是万变不离其宗，我想如果大家能理解常规的Union-Find的话，像写出weightedQuickUnionUF根本不成问题
最后祝大家学习进步~生活愉快哦，如果楼主哪里讲的不对或者大家有什么不懂得地方或着更好得思路欢迎在底下评论