算法实验报告1
[在此处输入文章标题]
算法分析与设计实验报告
实验一渗透问题（Percolation）
- 1. 实验目的
使用合并-查找（union-find）数据结构，编写程序通过蒙特卡罗模拟（Monte Carlo simulation） 来估计渗透阈值。
- 2. 实验环境
Java8 Eclipse algs4.jar包
- 3. 实验内容
我们使用 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*。
图1 实验内容示例
图2 可能性P的变化规律
- 4. 实验程序设计
Percolation数据类型。模型化一个Percolation系统，创建含有以下API的数据类型Percolation。
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
}
为了应用union-find的两种算法，我们这里初始化两个节点，一个标号为0，表示最上面的一个节点，另外一个标号为N*N+1表示最底下的一个节点，一开始的时候，我们将第一行的所有节点与0号相连，而最后一行的所有元素和N*N+1号节点相连，为了判断是否渗漏，只要判断0号节点与N*N+1号节点是否相连就行，这样简化了判断第一行是否和最后一行相连这个问题。
而一旦出现“open”的元素，我们判断他四个周围的节点，一旦有一个周围节点是open 的，我们就将此节点与他的周围节点连通。
而为了使用蒙特卡洛模拟，我们使用以下的算法分析方法：
1）初始化所有格点为 blocked。
2）重复以下操作直到系统渗出：
3）在所有 blocked 的格点之间随机均匀选择一个格点 (row i, column j)。
4）设置这个格点(row i, column j)为 open 格点open 格点的比例提供了系统渗透时渗透阈值的一个估计。
5）通过重复该计算实验 T 次并对结果求平均值，我们获得了更准确的渗滤阈值估计。 令 xt是第 t 次计算实验中 open 格点所占比例。样本均值μ提供渗滤阈值的一个估计值；样本标准差σ测量阈值的灵敏性，计算以下值
还有95%的置信区间
下面进入我们的核心程序：
1） Quickfind算法的实现
//并查集代码1
public class QuickFindUF {
private int[] id;
private int count;
public QuickFindUF(int N)
{
count = N;
id = new int[N];
for(int i=0;i<N;i++)
{
id[i] = i;
}
}
public int count()
{
return count;
}
public int find(int p)
{
return id[p];
}
public boolean connected(int p,int q)
{
return find(p) == find(q);
}
public void union(int p,int q)
{
int pid = find(p);
int qid = find(q);
if(pid == qid)
return;
for(int i =0;i<id.length;i++)
{
if(id[i] == pid)
id[i] = qid;
}
count--;
}
}
2） WeightedQuickUnion算法的实现
//并查集代码2
public class WeightedQuickUnionUF {
private int[] id;
private int[] sz;
private int count;
public WeightedQuickUnionUF(int N)
{
count = N;
id = new int[N];
for(int i=0;i<N;i++)
{
id[i] = i;
}
sz = new int[N];
for(int i=0;i<N;i++)
{
sz[i] = 1;
}
}
public int count()
{
return count;
}
public int find(int p)
{
while(p!=id[p])
p = id[p];
return p;
}
public boolean connected(int p,int q)//判断连通性
{
return find(p) == find(q);
}
public void union(int p,int q)//合并！
{
int pid = find(p);
int qid = find(q);
if(pid == qid)
return;
if(sz[pid]<sz[qid])
{
id[pid] = qid;
sz[qid]+= sz[pid];
}
else
{
id[qid] = pid;
sz[pid] += sz[qid];
}
count--;
}
}
3） Percolation类的实现，这里我给出了一份代码，表示的quick-find算法实现的代码，而weightedquickunion算法的代码以注释的形式，把与quick-find的代码不同的地方写出来，因为两份具体实现的代码非常相似，仅仅是调用了不同的并查集算法，也就是一个使用的是uf1实例，一个使用的是uf2实例
public class Percolation {
private int N;
private QuickFindUF uf1;
//private WeightedQuickUnionUF uf2;
private int grid[][];
public Percolation(int N)//0为最上方的节点，N*N+1自然是最下方的节点
{
this.N = N;
grid = new int[1+N][1+N];
uf1 = new QuickFindUF(N*N+2);
//uf2 = new WeightedQuickUnionUF(N*N+2);
for(int i=1;i<=N;i++)
{
for(int j=1;j<=N;j++)
{
grid[i][j] = 0;//一开始的时候都是block的状态
}
}
for(int i=1;i<=N;i++)//初始化，使得最上方与最下方的节点分别和第一行，最后一行连起来
{
uf1.union(0, i);//第一行
//uf2.union(0, i);//第一行
uf1.union(N*N+1, N*N-i);//最后一行
//uf2.union(N*N+1, N*N-i);//最后一行
}
}
public void open(int row,int col) {
grid[row][col] = 1;//置grid为1
int recentkey = (row-1)*N+col;
int key;//下面判断是否超出边界，以及算key,key代表Unionfind算法中的id值
if(row-1>=1&&row-1<=N&&col>=1&&col<=N)
{
key = (row-1-1)*N+col;
if(grid[row-1][col]==1)
{
uf1.union(recentkey, key);//union操作，将两个grid连通起来
//uf2.union(recentkey, key);
}
}
if(row+1>=1&&row+1<=N&&col>=1&&col<=N)
{
key = (row+1-1)*N+col;
if(grid[row+1][col]==1)
{
uf1.union(recentkey, key);
//uf2.union(recentkey, key);
}
}
if(row>=1&&row<=N&&col+1>=1&&col+1<=N)
{
key = (row-1)*N+col+1;
if(grid[row][col+1]==1)
{
uf1.union(recentkey, key);
//uf2.union(recentkey, key);
}
}
if(row>=1&&row<=N&&col-1>=1&&col-1<=N)
{
key = (row-1)*N+col-1;
if(grid[row][col-1]==1)
{
uf1.union(recentkey, key);
//uf2.union(recentkey, key);
}
}
}
public boolean isOpen(int row,int col)
{
return grid[row][col] == 1;
}
public boolean isFull(int row,int col)//是否与最上方的0号节点连通
{
return uf1.connected(0, (row-1)*N+col)&&isOpen(row,col);
}
public boolean percolates()//仅判断第一和最后即可
{
return uf1.connected(0, N*N+1);
}
}
- 5. 实验具体代码以及运行结果截图
全部的代码如下
public class Percolation {
private int N;
private QuickFindUF uf1;
//private WeightedQuickUnionUF uf2;
private int grid[][];
public Percolation(int N)//0为最上方的节点，N*N+1自然是最下方的节点
{
this.N = N;
grid = new int[1+N][1+N];
uf1 = new QuickFindUF(N*N+2);
//uf2 = new WeightedQuickUnionUF(N*N+2);
for(int i=1;i<=N;i++)
{
for(int j=1;j<=N;j++)
{
grid[i][j] = 0;//一开始的时候都是block的状态
}
}
for(int i=1;i<=N;i++)//初始化，使得最上方与最下方的节点分别和第一行，最后一行连起来
{
uf1.union(0, i);//第一行
//uf2.union(0, i);//第一行
uf1.union(N*N+1, N*N-i);//最后一行
//uf2.union(N*N+1, N*N-i);//最后一行
}
}
public void open(int row,int col) {
grid[row][col] = 1;//置grid为1
int recentkey = (row-1)*N+col;
int key;//下面判断是否超出边界，以及算key,key代表Unionfind算法中的id值
if(row-1>=1&&row-1<=N&&col>=1&&col<=N)
{
key = (row-1-1)*N+col;
if(grid[row-1][col]==1)
{
uf1.union(recentkey, key);//union操作，将两个grid连通起来
//uf2.union(recentkey, key);
}
}
if(row+1>=1&&row+1<=N&&col>=1&&col<=N)
{
key = (row+1-1)*N+col;
if(grid[row+1][col]==1)
{
uf1.union(recentkey, key);
//uf2.union(recentkey, key);
}
}
if(row>=1&&row<=N&&col+1>=1&&col+1<=N)
{
key = (row-1)*N+col+1;
if(grid[row][col+1]==1)
{
uf1.union(recentkey, key);
//uf2.union(recentkey, key);
}
}
if(row>=1&&row<=N&&col-1>=1&&col-1<=N)
{
key = (row-1)*N+col-1;
if(grid[row][col-1]==1)
{
uf1.union(recentkey, key);
//uf2.union(recentkey, key);
}
}
}
public boolean isOpen(int row,int col)
{
return grid[row][col] == 1;
}
public boolean isFull(int row,int col)//是否与最上方的0号节点连通
{
return uf1.connected(0, (row-1)*N+col)&&isOpen(row,col);
}
public boolean percolates()//仅判断第一和最后即可
{
return uf1.connected(0, N*N+1);
}
}
//并查集代码1
public class QuickFindUF {
private int[] id;
private int count;
public QuickFindUF(int N)
{
count = N;
id = new int[N];
for(int i=0;i<N;i++)
{
id[i] = i;
}
}
public int count()
{
return count;
}
public int find(int p)
{
return id[p];
}
public boolean connected(int p,int q)
{
return find(p) == find(q);
}
public void union(int p,int q)
{
int pid = find(p);
int qid = find(q);
if(pid == qid)
return;
for(int i =0;i<id.length;i++)
{
if(id[i] == pid)
id[i] = qid;
}
count--;
}
}
//并查集代码2
public class WeightedQuickUnionUF {
private int[] id;
private int[] sz;
private int count;
public WeightedQuickUnionUF(int N)
{
count = N;
id = new int[N];
for(int i=0;i<N;i++)
{
id[i] = i;
}
sz = new int[N];
for(int i=0;i<N;i++)
{
sz[i] = 1;
}
}
public int count()
{
return count;
}
public int find(int p)
{
while(p!=id[p])
p = id[p];
return p;
}
public boolean connected(int p,int q)
{
return find(p) == find(q);
}
public void union(int p,int q)
{
int pid = find(p);
int qid = find(q);
if(pid == qid)
return;
if(sz[pid]<sz[qid])
{
id[pid] = qid;
sz[qid]+= sz[pid];
}
else
{
id[qid] = pid;
sz[pid] += sz[qid];
}
count--;
}
}
import edu.princeton.cs.algs4.StdRandom;
//统计类，使用加权并查集,这里仅仅需要改变其中的并查集类就行了,跟quickfind的统计类大致形似
public class PercolationStatsWeighted {
private int T_number;
private int matrixLength;
private double[] threshold;
private double mean; //平均值
private double stddev; //标准偏差
private double confidenceLow; // 最低置信度
private double confidenceHigh; // 高置信度
public PercolationStatsWeighted(int N, int T)
{// perform T independent computational experiments on an N-by-N grid
matrixLength = N;
T_number = T;
mean = 0;
threshold = new double[T];
for(int i=0;i<T;i++)
{
PercolationWeighted percolation = new PercolationWeighted(N);
int count = 0;
do {
int row = StdRandom.uniform(N)+1;
int col = StdRandom.uniform(N)+1;
if(percolation.isOpen(row, col))
{
continue;
}
else
{
++count;
percolation.open(row, col);
}
}while(!percolation.percolates());
threshold[i] = (double)1.0*count/(N*N);
}
}
public double mean() // sample mean of percolation threshold
{
for(int i=0;i<threshold.length;i++)
{
mean+=threshold[i];
}
return mean=mean/threshold.length;
}
public double stddev()
{// sample standard deviation of percolation threshold
double result = 0;
for(int i=0;i<threshold.length;i++)
{
result+=(threshold[i]-mean)*(threshold[i]-mean);
}
result/=(T_number-1);
return stddev=Math.sqrt(result);//赋值并且返回结果
}
public double confidenceLo()
{// returns lower bound of the 95% confidence interval
return confidenceLow=(mean-stddev*1.96/Math.sqrt(T_number));//赋值并返回结果
}
public double confidenceHi()
{// returns upper bound of the 95% confidence interval
return confidenceHigh=(mean+stddev*1.96/Math.sqrt(T_number));//赋值并返回结果
}
}
package percolation;
import edu.princeton.cs.algs4.StdRandom;
public class PercolationStats {
private int T_number;
private int matrixLength;
private double[] threshold;
private double mean; //平均值
private double stddev; //标准偏差
private double confidenceLow; // 最低置信度
private double confidenceHigh; // 高置信度
public PercolationStats(int N, int T)
{// perform T independent computational experiments on an N-by-N grid
matrixLength = N;
T_number = T;
mean = 0;
threshold = new double[T];
for(int i=0;i<T;i++)
{
Percolation percolation = new Percolation(N);
int count = 0;
do {
int row = StdRandom.uniform(N)+1;
int col = StdRandom.uniform(N)+1;
if(percolation.isOpen(row, col))
{
continue;
}
else
{
++count;
percolation.open(row, col);
}
}while(!percolation.percolates());
threshold[i] = (double)1.0*count/(N*N);
}
}
public double mean() // sample mean of percolation threshold
{
for(int i=0;i<threshold.length;i++)
{
mean+=threshold[i];
}
return mean=mean/threshold.length;
}
public double stddev()
{// sample standard deviation of percolation threshold
double result = 0;
for(int i=0;i<threshold.length;i++)
{
result+=(threshold[i]-mean)*(threshold[i]-mean);
}
result/=(T_number-1);
return stddev=Math.sqrt(result);//赋值并且返回结果
}
public double confidenceLo()
{// returns lower bound of the 95% confidence interval
return confidenceLow=(mean-stddev*1.96/Math.sqrt(T_number));//赋值并返回结果
}
public double confidenceHi()
{// returns upper bound of the 95% confidence interval
return confidenceHigh=(mean+stddev*1.96/Math.sqrt(T_number));//赋值并返回结果
}
public static void main(String[] args) {
// TODO 自动生成的方法存根
for(int x=10;x<=2000;x*=2)
{
long starttime = System.currentTimeMillis();
PercolationStats quickfind = new PercolationStats(x,30);
long endtime = System.currentTimeMillis();
long realtime = endtime - starttime;
System.out.println("使用quickfind算法实现，所使用的n为 "+x+",t是30");
System.out.println("mean() "+quickfind.mean());
System.out.println("stddev() "+quickfind.stddev());
System.out.println("confidenceLo() "+quickfind.confidenceLo());
System.out.println("confidenceHi() "+quickfind.confidenceHi());
System.out.println("运行时间"+realtime+"毫秒");
starttime = System.currentTimeMillis();
PercolationStatsWeighted Weightedqucikunion = new PercolationStatsWeighted(x,30);
endtime = System.currentTimeMillis();
realtime = endtime - starttime;
System.out.println("使用WeightedquickUnion算法实现，所使用的n为 "+x+",t是30");
System.out.println("mean() "+Weightedqucikunion.mean());
System.out.println("stddev() "+Weightedqucikunion.stddev());
System.out.println("confidenceLo() "+Weightedqucikunion.confidenceLo());
System.out.println("confidenceHi() "+Weightedqucikunion.confidenceHi());
System.out.println("运行时间"+realtime+"毫秒");
System.out.println("---------------------分割线-----------------------------");
}
}
}
实验运行结果截图如下
- 6. 实验结果分析
N=10的时候，我们进行如下分析
|
算法时间（ms） |
n |
2n |
4n |
8n |
16n |
32n |
64n |
|
Quickfind |
15 |
7 |
46 |
259 |
4379 |
62345 |
1203317 |
|
Weightedquickunion |
5 |
3 |
4 |
5 |
46 |
241 |
1350 |
我们都知道Quickfind应该是N2的算法，下面，我们对两边取lg对数，然后进行线性回归
此时lgT = lg（3.9192）+2.912lgN + lg30
其中30 = t
所以得到 T = 3.9192*t*N^(2.9192)，此时，我们发现，在这种条件下
总时间T与N^(2.9192)成了正比。
而对Weigtedquickunion，我们直接进行线性回归
在matlab中，输入
得到 T = -136.7414+2.0561*N，所以WeightedQuickUnion此时近似为O（2N）的算法，而如果算上t，那么
T = （-136.7414+2.0561*N）/30*t
