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

代码由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 + 1UnionFind 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)