使用蒙特卡洛模拟和并查集求解渗透阈值

代码由C++语言编写

我使用了三个类，分别如下

第一个类，实现了并查集的基本操作。

#pragma once

#include<iostream>
#include<vector>
using namespace std;

class UnionFind
{
public:
	vector<int> parent; // 用于存储父节点
	vector<int> rank; // 结点所在集合的树高

	// 构造函数
	UnionFind(int size);

	// 查询根节点
	int find(int x);

	// 合并集合
	void unionSets(int x, int y);

	// 检查两个结点是否在同一个集合中
	bool connected(int x, int y);

};

#include "UnionFind.h"


// 构造函数，初始化parent和rank数组
UnionFind::UnionFind(int size) : parent(size), rank(size, 1) 
{
	for (int i = 0; i < size; ++i) 
		parent[i] = i;
}

// find函数，返回父节点的同时实现路径压缩
// 压缩：将路径上所有的节点都连接到最终的根节点上
// 递归实现
int UnionFind::find(int x)
{
	if (parent[x] != x)
		parent[x] = find(parent[x]);
	return parent[x];
}

// 合并集合
// 将集合树低的集合并到树高的集合上，避免树太高
void UnionFind::unionSets(int x, int y) 
{
    int rootX = find(x);
    int rootY = find(y);
    if (rootX != rootY) 
    {
        if (rank[rootX] < rank[rootY]) 
            parent[rootX] = rootY;
        else if (rank[rootX] > rank[rootY])
            parent[rootY] = rootX;
        else 
        {
            parent[rootX] = rootY;
            rank[rootX]++;
        }
    }
}

// 单纯检查两个结点是否在同一个集合里
bool UnionFind::connected(int x, int y)
{
	return find(x) == find(y);
}

第二个类，实现给定具体矩阵时，判断是否渗透。

#pragma once
#include<iostream>
#include<vector>
using namespace std;
# include "UnionFind.h"

class Solve
{
public:
	int m_rows; // 成员变量，给定矩阵的行和列
	int m_cols;
	vector<vector<int>> m_grid; // 拷贝一份给定矩阵

	// 利用并查集的方法将矩阵中的0划分为几个集合
	Solve(vector<vector<int>> grid);

	// 判断头结点和尾结点是否在一个集合里，也就是返回是否渗透
	bool isConnected();

};

#include "Solve.h"


Solve::Solve(vector<vector<int>> grid)
{
	m_grid = grid;
	m_rows = m_grid.size();
	m_cols = m_grid[0].size();
}


bool Solve::isConnected()
{
	// 额外两个节点，分别用于连接第一行和最后一行
	//  m_rows * m_cols连接第一行，m_rows * m_cols + 1连接最后一行
	// 给矩阵编号为
	// m_rows * m_cols
	// 0	1	 2
	// 3	4	 5
	// 6	7	 8
	// m_rows * m_cols + 1
	UnionFind uf(m_rows * m_cols + 2);
	int head = m_rows * m_cols;
	int tail = m_rows * m_cols + 1;
	// unionSets函数最终是以第二个参数的根节点作为合并后的根节点的
	// 初始化第一行
	for (int i = 0; i < m_cols; i++)
	{
		if (m_grid[0][i] == 0)
			uf.unionSets(i, head);
	}
	// 初始化最后一行
	for (int i = 0; i < m_cols; i++)
	{
		if (m_grid[m_rows - 1][i] == 0)
			uf.unionSets((m_rows - 1) * m_cols + i, tail);
	}
	for (int i = 1; i < m_rows; i++)
	{
		for (int j = 0; j < m_cols; j++)
		{
			if (m_grid[i][j] == 0)
			{
				// 将相邻的0合并为一个集合
				if (m_grid[i - 1][j] == 0)
					uf.unionSets(i * m_cols + j, (i - 1) * m_cols + j);
				if (j - 1 >= 0 && m_grid[i][j - 1] == 0)
					uf.unionSets(i * m_cols + j, i * m_cols + j - 1);
			}
		}
	}
	return uf.connected(head, tail);
}

第三个类，用于生成随机矩阵

#pragma once
#include <iostream>
#include <vector>  
#include <random>  
using namespace std;

class CreatMatrix
{
public:
	int m_rows;
	int m_cols;
	double m_prob_zero;

	CreatMatrix(int rows, int cols, double prob_zero);

	vector<vector<int>> creatMatrix();

};

#include "CreatMatrix.h"

CreatMatrix::CreatMatrix(int rows, int cols, double prob_zero)
{
    // 假设的数组大小  
    m_rows = rows;
    m_cols = cols;

    // 0出现的概率  
    m_prob_zero = prob_zero;
}

vector<vector<int>> CreatMatrix::creatMatrix()
{
    // 创建一个二维vector  
    vector<vector<int>> matrix(m_rows, vector<int>(m_cols));

    // 初始化随机数生成器  
    random_device rd;  // 用于获取随机数种子  
    mt19937 gen(rd()); // 以随机数种子初始化Mersenne Twister生成器  
    uniform_real_distribution<> dis(0.0, 1.0); // 分布范围在[0.0, 1.0)  

    // 填充二维数组  
    for (int i = 0; i < m_rows; ++i) {
        for (int j = 0; j < m_cols; ++j) {
            // 生成一个[0.0, 1.0)之间的随机数  
            double random_value = dis(gen);

            // 根据概率决定是0还是1  
            if (random_value < m_prob_zero) {
                matrix[i][j] = 0;
            }
            else {
                matrix[i][j] = 1;
            }
        }
    }

    return matrix;
}

主函数

#include "UnionFind.h"
#include "Solve.h"
#include "example.h"
#include "CreatMatrix.h"
#include <chrono>  
#include <fstream>
#include <string>

using namespace std;

int main()
{
	double path_length = 0.01; // 步长，用于for循环中0结点出现概率的递增
	int num_experiment = 1000; //对每种概率下的矩阵判断渗透的次数
	int matrix_row = 200; // 矩阵行
	int matrix_col = 200; // 矩阵列

	// 记录开始时间  
	auto start = std::chrono::high_resolution_clock::now();

	// 创建csv文件
	// CSV（逗号分隔值）文件是一种非常通用的文件格式，可以被多种程序读取，包括Python。
	
	// 定义CSV文件的路径  
	string csvFilePath = "your_path";

	// 创建一个ofstream对象来写入文件  
	ofstream csvFile(csvFilePath);

	if (!csvFile.is_open()) {
		cerr << "Unable to open file";
		return 1;
	}

	// 创建两个数组
	// 第一个是每次零出现的概率
	// 第二个是在这个概率下，实验一千次有多少次成功连通
	int documents = 1.0 / path_length; // 记录次数

	vector<double> probabilities;
	probabilities.reserve(documents);
	vector<int> counts;
	counts.reserve(documents);

	// 写入CSV头部
	csvFile << "Probability,Connected Count\n";
	for (double prob = 0; prob <= 1; prob += path_length)
	{
		// 用于记录成功渗透的次数
		int count = 0;
		for (int times = 0; times < num_experiment; times++)
		{
			CreatMatrix cm(matrix_row, matrix_col, prob);
			vector<vector<int>> grid = cm.creatMatrix();
			Solve solve(grid);
			if (solve.isConnected())
				count++;
			cout << "这是" << prob << "概率下第" << times << "次实验" << endl;
		}
		probabilities.push_back(prob);
		counts.push_back(count);
	}

	for (size_t i = 0; i < documents; i++) {
		csvFile << probabilities[i] << "," << counts[i] << "\n";
	}

	csvFile.close();

	std::cout << "CSV文件已创建并写入数据。" << std::endl;
	//if (solve.isConnected()) 
	//	cout << "第一行和最后一行的0是连通的" << endl;
	//else 
	//	cout << "第一行和最后一行的0不是连通的" << endl;

	// 记录结束时间  
	auto end = std::chrono::high_resolution_clock::now();

	// 计算时间差  
	std::chrono::duration<double, std::milli> elapsed = end - start;

	// 输出运行时间（以毫秒为单位）  
	std::cout << "程序运行时间: " << elapsed.count() << " 毫秒" << std::endl;
	return 0;
}

csv文件部分截图

最后使用python读取这个文件，并绘图

确保你的python安装了pandas和matplotlib库

pip install pandas -i https://pypi.tuna.tsinghua.edu.cn/simple
pip install matplotlib -i https://pypi.tuna.tsinghua.edu.cn/simple

import pandas as pd
import matplotlib.pyplot as plt

# 读取CSV文件
path = r"your_path"
data = pd.read_csv(path)

# 绘图
# 长10，宽5
plt.figure(figsize=(10, 6))
plt.plot(data['Probability'], data['Connected Count'], marker='o')
plt.title('Connectivity Counts vs Probability')
plt.xlabel('Probability')
plt.ylabel('Connected Count')
plt.grid(True)
plt.show()

结果截图(200 * 200 , 0.01, 1000)