最近老师给了一个题目，说是研究一个正常矩阵任意概率置点概率下，双向导通（x,y）的概率（要求有自然边界条件，也就是可以从0->length-1），用代码敲了一下demo，结果发现有个好有趣的结果：不同大小的矩阵，导通概率从某一个概率上升，点越多，上升程度越快（斜率越大），但是都会在（0.592...，35%个点导通）左右相交，符合ziff论文的结果。
其实这个结果是统计学的一个结果，不过我们还没学概率论，这个神奇的结论我还不能证明（不过这个证明也之能证明二维的，对于三维以上只能通过模拟得出）。
这个结果对于不同的图都有不同的结果，我模拟的这个是最简单的图（自然边界条件矩形图）。
下面是demo，没有什么特别的算法，也就是DFS稍微优化一下。
1 #include <iostream> 2 #include <functional> 3 #include <algorithm> 4 #include <time.h> 5 #define MARTIX_SIZE 256 6 #define NP 6 7 8 using namespace std; 9 10 typedef struct _tempset 11 { 12 bool flag;//是否被访问过 13 bool dir_x;//规定x一个方向 14 bool dir_y;//规定x一个方向 15 bool is_dot;//是否是一个节点 16 }Set; 17 typedef struct _group 18 { 19 int particle_sum; 20 bool Link_px; 21 bool Link_py; 22 }Group; 23 typedef struct _temp_recur 24 { 25 int x, y, dir; 26 }T_Group; 27 typedef int Position; 28 29 void Inivilize(Set ***const, T_Group **const,Position **const,Position **const); 30 Set **ReFresh_Martix(Set **const, int *const, Position *const, Position *const, int); 31 void Draw_The_Gragh(Set **const); 32 void Destory(Set **const, T_Group *, Position *, Position *); 33 void If_Cls(void); 34 void Show(const int, const int, const int, const double,FILE *); 35 void DFS_Martix(Set **, const int, const int, T_Group *const, Group *, Position *const, Position *const); 36 37 static double poss[NP] = { 0.590, 0.5925, 0.5926, 0.5927, 0.5928, 0.5930}; 38 39 static pair<int,int>direction[4]; 40 41 void Inivilize(Set ***const martix_tmp, 42 T_Group **const G_Stack, Position **const Mark_x, Position **const Mark_y) 43 { 44 /*函数名: Inivilize 45 *调用者函数: void main(void) 46 *函数形参: martix_tmp（矩阵指针） 47 G_Stack（伪递归栈堆指针） 48 Mark_x，Mark_y（GPA优化的断层标记数组指针） 49 *函数功能: 初始化 50 *返回值： 无 51 */ 52 *martix_tmp = (Set **)malloc(sizeof(Set *)*MARTIX_SIZE); 53 Set *Zone = (Set *)malloc(sizeof(Set)*MARTIX_SIZE*MARTIX_SIZE); 54 for (int i = 0; i < MARTIX_SIZE; i++) 55 (*martix_tmp)[i] = &Zone[i*MARTIX_SIZE]; 56 57 *G_Stack = (T_Group *)malloc(sizeof(T_Group)*MARTIX_SIZE*MARTIX_SIZE);//整形最大值 58 *Mark_x = (Position *)malloc(sizeof(Position)*MARTIX_SIZE); 59 *Mark_y = (Position *)malloc(sizeof(Position)*MARTIX_SIZE); 60 61 //............................................... 62 srand((unsigned)time(NULL)); 63 direction[0].first = -1; direction[0].second = 0; 64 direction[1].first = 1; direction[1].second = 0; 65 direction[2].first = 0; direction[2].second = -1; 66 direction[3].first = 0; direction[3].second = 1; 67 } 68 69 Set **ReFresh_Martix(Set **const Martix, int *const sum_particle, 70 Position *const Mark_x, Position *const Mark_y, int k) 71 { 72 /*函数名: ReFresh_Martix 73 *调用者函数: void main(void) 74 *函数形参: Martix（矩阵指针） 75 sum_particle（新产生的粒子数） 76 Mark_x，Mark_y（GPA优化的断层标记数组指针） 77 k（以第k个概率产生点） 78 *函数功能: 给矩阵产生点，以及赋值 79 *返回值： Martix（矩阵指针） 80 */ 81 int tmp; 82 memset(Mark_x, 0, sizeof(Position)*MARTIX_SIZE); 83 memset(Mark_y, 0, sizeof(Position)*MARTIX_SIZE); 84 for (int i = 0; i < MARTIX_SIZE; i++) 85 { 86 for (int j = 0; j < MARTIX_SIZE; j++) 87 { 88 tmp = rand(); 89 Martix[i][j].flag = tmp < RAND_MAX * poss[k] ? 1 : 0; 90 Martix[i][j].dir_y = 0; Martix[i][j].dir_x = 0;//强行规定一个方向 91 Martix[i][j].is_dot = 0; 92 if (Martix[i][j].flag) 93 { 94 (*sum_particle)++; 95 Mark_x[j]++; Mark_y[i]++; 96 Martix[i][j].is_dot = 1; 97 } 98 } 99 } 100 return Martix; 101 } 102 103 void Destory(Set **const Martix, T_Group *G_Stack, Position *Mark_x, Position *Mark_y) 104 { 105 /*函数名: Destory 106 *调用者函数: void main(void) 107 *函数形参: Martix（矩阵指针） 108 G_Stack（伪递归栈堆指针） 109 Mark_x，Mark_y（GPA优化的断层标记数组指针） 110 *函数功能: 销毁 111 *返回值： 无 112 */ 113 free(Martix[0]);//去掉集体分配的首地址就可以了 114 free(Martix); 115 free(G_Stack);//销毁递归栈堆 116 free(Mark_x); free(Mark_y);//销毁点集 117 } 118 119 void If_Cls(void) 120 { 121 /*函数名: If_Cls 122 *调用者函数: void main(void) 123 *函数形参: 无 124 *函数功能: 询问是否清屏 125 *返回值： 无 126 */ 127 char choice; 128 printf("Do you want to clean the screen?(Y or N)"); 129 while (1) 130 { 131 choice = getchar(); 132 if (choice == 'Y') 133 { 134 system("cls"); 135 return; 136 } 137 else if (choice == 'N') 138 return; 139 else 140 { 141 puts("Error Input!Please enter again! "); 142 fflush(stdin); 143 } 144 } 145 } 146 147 void Show(const int k, const int link_group_sum, const int loop, const double max_poss, FILE *fp) 148 { 149 /*函数名: Show 150 *调用者函数: void main(void) 151 *函数形参: k（以第k个概率产生点） 152 link_group_sum（联通集团个数） 153 loop（循环次数） 154 max（联通最大点个数） 155 p_sum（所有联通集团的粒子总数） 156 fp（指向data.txt的文件指针） 157 *函数功能: 显示当前概率集团信息，并且把数据写入data.txt中 158 *返回值： 无 159 */ 160 puts("=============================================================================== "); 161 printf("a%cEcho %d: Possibility of %.4f The linked graphs` sum is %d Account for %.2f%% in %d graghs ", 162 16, k + 1, poss[k], link_group_sum, 100 * (double)link_group_sum / (double)loop, loop); 163 printf("The possibility of(max linked points/sum points)accounts for %.2f%% ", 100 * max_poss); 164 printf(" "); 165 166 if (fp != NULL) 167 { 168 fprintf(fp, "=============================================================================== "); 169 fprintf(fp, "%cEcho %d: Possibility of %.4f The linked graphs` sum is %d Account for %.2f%% in %d graghs ", 170 7, k + 1, poss[k], link_group_sum, 100 * (double)link_group_sum / (double)loop, loop); 171 fprintf(fp, "The possibility of(max linked points/sum points)accounts for %.2f%% ", 100 * max_poss); 172 fprintf(fp, " "); 173 } 174 }
#include "plug.h" static bool x_full, y_full; int main(void)//森林火灾项目的计算连通性的概率 { int loop, sum_particle, link_group_sum, max_particle_sum; double max_poss, poss_tmp; Set **martix = NULL; //矩阵 Group tmp; T_Group *G_Stack; //范围太大，不用递归，手动直接开大栈提高效率 Position *Mark_x = NULL, *Mark_y = NULL; //Gpa优化标记数组 FILE *fp = NULL; //文件指针 //...................................................................................... Inivilize(&martix, &G_Stack, &Mark_x, &Mark_y); //初始化 while (1) { printf("Please enter the times of the loops(The martix`s size is %d*%d) (The program will stop until you input zero) ", MARTIX_SIZE, MARTIX_SIZE); while (1) { scanf("%d", &loop); if (loop >= 0)break; puts("DO NOT try to input a negative loop"); system("cls"); puts("Please enter the time of the loop(break until you input 0)"); } if (loop == 0) break; fflush(stdin); puts("The Date will be written in ""data.txt"" in the folder of this program "); if ((fp = fopen("data.txt", "r")) == NULL) puts("""data.txt"" doesn`t exsit,the program will create a new one "); fp = fopen("data.txt", "a+"); time_t c1 = clock(); //这里先得到一个循环开始的时间 for (int k = 0; k < NP; k++) { link_group_sum = 0; max_poss = 0; /*DFS: 递归每一个概率，和每一张图，然后每一个点深度优先搜索看是否联通 *GPA优化：如果Mark数组出现断层，则此图一定不连通(x_full和y_full) */ for (int i = 0; i < loop; i++) { sum_particle = 0; x_full = 1; y_full = 1; martix = ReFresh_Martix(martix, &sum_particle, Mark_x, Mark_y, k); max_particle_sum = 0; for (int y = 0; y < MARTIX_SIZE && x_full && y_full; y++) { for (int x = 0; x < MARTIX_SIZE && x_full && y_full; x++) { if (martix[y][x].flag) { DFS_Martix(martix, x, y, G_Stack, &tmp, Mark_x, Mark_y); if (tmp.Link_px && tmp.Link_py) { link_group_sum++; if (tmp.particle_sum > max_particle_sum) max_particle_sum = tmp.particle_sum; } } } } poss_tmp = (double)max_particle_sum / (double)sum_particle; max_poss = poss_tmp > max_poss ? poss_tmp : max_poss; } Show(k, link_group_sum, loop, max_poss, fp); //联通的概率（loop个图），联通的时候最大联通集团粒子数 } time_t c2 = clock(); //得到循环结束的时间（ms） printf("The program has run %lldms ", c2 - c1); fclose(fp); //记得关闭文件流 system("pause"); If_Cls(); } Destory(martix, G_Stack, Mark_x, Mark_y); //销毁栈组 fflush(stdin); system("cls"); system("pause"); return 0; } void DFS_Martix(Set **martix, const int st_x, const int st_y, T_Group *const G_Stack, Group *Part_Group, Position *const Mark_x, Position *const Mark_y) { /*函数名: DFS_Martix *调用者函数: void main(void) *函数形参: martix（矩阵） st_x，st_y（开始坐标横纵坐标） G_Stack（伪递归栈堆） Part_Group（集团参数指针） Mark_x，Mark_y（GPA优化的断层标记数组） *函数功能: DFS搜索联通集团 *返回值： 无 */ bool If_Push, //出入栈标志 if_link_y = 0, if_link_x = 0; //图是否联通标志 bool round_dir_x = 0, round_dir_y = 0, //环图指向标记：规定沿着某个轴沿一个方向突然发生反转时，某个方向就会被标记1 tmp_dir_x, tmp_dir_y; int dir_tmp, particle_sum = 0, //粒子数 top = 0; //栈顶位 Position dir_x, dir_y; T_Group tmp; //...................................................................................... //初始化栈.......... martix[st_y][st_x].flag = 0; martix[st_y][st_x].dir_x = 0; martix[st_y][st_x].dir_y = 0; Mark_x[st_x]--; Mark_y[st_y]--; tmp.x = st_x, tmp.y = st_y; tmp.dir = dir_tmp = 0; G_Stack[top++] = tmp; particle_sum++; //.................. Part_Group->Link_px = 0; Part_Group->Link_py = 0; while (top != 0) { If_Push = 0; tmp = G_Stack[top - 1]; //读取栈顶元素,以及栈顶元素环图两个指向标记 round_dir_x = martix[tmp.y][tmp.x].dir_x; round_dir_y = martix[tmp.y][tmp.x].dir_y; //每弹出一个点，就把该节点的x和y做标记,直接从上一个搜索方向之后开始搜索 for (int i = dir_tmp; i < 4; i++) { dir_x = tmp.x + direction[i].first; tmp_dir_x = (dir_x < 0 || dir_x == MARTIX_SIZE) ? !round_dir_x : round_dir_x; dir_x = (dir_x + MARTIX_SIZE) % MARTIX_SIZE; //循环坐标 dir_y = tmp.y + direction[i].second; tmp_dir_y = (dir_y < 0 || dir_y == MARTIX_SIZE) ? !round_dir_y : round_dir_y; dir_y = (dir_y + MARTIX_SIZE) % MARTIX_SIZE; //循环坐标 if (martix[dir_y][dir_x].is_dot && !martix[dir_y][dir_x].flag) { //如果环图指向标记出现相反，则图已经联通 if (!if_link_x) if_link_x = martix[dir_y][dir_x].dir_x == tmp_dir_x ? 0 : 1; if (!if_link_y) if_link_y = martix[dir_y][dir_x].dir_y == tmp_dir_y ? 0 : 1; } if (martix[dir_y][dir_x].flag) { If_Push = 1; Mark_x[dir_x]--; Mark_y[dir_y]--; round_dir_x = tmp_dir_x; round_dir_y = tmp_dir_y; martix[dir_y][dir_x].flag = 0; martix[dir_y][dir_x].dir_x = round_dir_x; martix[dir_y][dir_x].dir_y = round_dir_y; if (!Mark_x[dir_x]) x_full = 0; if (!Mark_y[dir_y]) y_full = 0; tmp.x = dir_x; tmp.y = dir_y; tmp.dir = i; G_Stack[top++] = tmp; particle_sum++; break; } } dir_tmp = 0; //这里一定要初始化为0，否则接下来的点就无法搜索其他方向 if (!If_Push) //如果没有可以入栈的，那说明开始递归了，pop栈顶 { top--; dir_tmp = tmp.dir; //定义为当前点上一次出去的点的方向 } } if (particle_sum >= MARTIX_SIZE) { Part_Group->Link_px = if_link_x; Part_Group->Link_py = if_link_y; } Part_Group->particle_sum = particle_sum; }