网络导通概率的研究

　　最近老师给了一个题目，说是研究一个正常矩阵任意概率置点概率下，双向导通（x,y）的概率（要求有自然边界条件，也就是可以从0->length-1），用代码敲了一下demo，结果发现有个好有趣的结果：不同大小的矩阵，导通概率从某一个概率上升，点越多，上升程度越快（斜率越大），但是都会在（0.592...，35%个点导通）左右相交，符合ziff论文的结果。

　　其实这个结果是统计学的一个结果，不过我们还没学概率论，这个神奇的结论我还不能证明（不过这个证明也之能证明二维的，对于三维以上只能通过模拟得出）。

　　这个结果对于不同的图都有不同的结果，我模拟的这个是最简单的图（自然边界条件矩形图）。

　　下面是demo，没有什么特别的算法，也就是DFS稍微优化一下。

#include <iostream>
#include <functional>
#include <algorithm>
#include <time.h>
#define MARTIX_SIZE 256
#define NP 6

using namespace std;

typedef struct _tempset
{
    bool flag;//是否被访问过
    bool dir_x;//规定x一个方向
    bool dir_y;//规定x一个方向
    bool is_dot;//是否是一个节点
}Set;
typedef struct _group
{
    int particle_sum;
    bool Link_px;
    bool Link_py;
}Group;
typedef struct _temp_recur
{
    int x, y, dir;
}T_Group;
typedef int Position;

void Inivilize(Set ***const, T_Group **const,Position **const,Position **const);
Set **ReFresh_Martix(Set **const, int *const, Position *const, Position *const, int);
void Draw_The_Gragh(Set **const);
void Destory(Set **const, T_Group *, Position *, Position *);
void If_Cls(void);
void Show(const int, const int, const int, const double,FILE *);
void DFS_Martix(Set **, const int, const int, T_Group *const, Group *, Position *const, Position *const);

static double poss[NP] = { 0.590, 0.5925, 0.5926, 0.5927, 0.5928, 0.5930};

static pair<int,int>direction[4];

void Inivilize(Set ***const martix_tmp,
    T_Group **const G_Stack, Position **const Mark_x, Position **const Mark_y)
{
    /*函数名:       Inivilize
    *调用者函数:    void main(void)
    *函数形参:      martix_tmp（矩阵指针）
                    G_Stack（伪递归栈堆指针）
                    Mark_x，Mark_y（GPA优化的断层标记数组指针）
    *函数功能:      初始化
    *返回值：       无
    */
    *martix_tmp = (Set **)malloc(sizeof(Set *)*MARTIX_SIZE);
    Set *Zone = (Set *)malloc(sizeof(Set)*MARTIX_SIZE*MARTIX_SIZE);
    for (int i = 0; i < MARTIX_SIZE; i++)
        (*martix_tmp)[i] = &Zone[i*MARTIX_SIZE];

    *G_Stack = (T_Group *)malloc(sizeof(T_Group)*MARTIX_SIZE*MARTIX_SIZE);//整形最大值
    *Mark_x = (Position *)malloc(sizeof(Position)*MARTIX_SIZE);
    *Mark_y = (Position *)malloc(sizeof(Position)*MARTIX_SIZE);

    //...............................................
    srand((unsigned)time(NULL));
    direction[0].first = -1; direction[0].second = 0;
    direction[1].first = 1; direction[1].second = 0;
    direction[2].first = 0; direction[2].second = -1;
    direction[3].first = 0; direction[3].second = 1;
}

Set **ReFresh_Martix(Set **const Martix, int *const sum_particle,
    Position *const Mark_x, Position *const Mark_y, int k)
{
    /*函数名:       ReFresh_Martix
    *调用者函数:    void main(void)
    *函数形参:      Martix（矩阵指针）
                    sum_particle（新产生的粒子数）
                    Mark_x，Mark_y（GPA优化的断层标记数组指针）
                    k（以第k个概率产生点）
    *函数功能:      给矩阵产生点，以及赋值
    *返回值：       Martix（矩阵指针）
    */
    int tmp;
    memset(Mark_x, 0, sizeof(Position)*MARTIX_SIZE);
    memset(Mark_y, 0, sizeof(Position)*MARTIX_SIZE);
    for (int i = 0; i < MARTIX_SIZE; i++)
    {
        for (int j = 0; j < MARTIX_SIZE; j++)
        {
            tmp = rand();
            Martix[i][j].flag = tmp < RAND_MAX * poss[k] ? 1 : 0;
            Martix[i][j].dir_y = 0; Martix[i][j].dir_x = 0;//强行规定一个方向
            Martix[i][j].is_dot = 0;
            if (Martix[i][j].flag)
            {
                (*sum_particle)++;
                Mark_x[j]++; Mark_y[i]++;
                Martix[i][j].is_dot = 1; 
            }
        }
    }
    return Martix;
}

void Destory(Set **const Martix, T_Group *G_Stack, Position *Mark_x, Position *Mark_y)
{
    /*函数名:       Destory
    *调用者函数:    void main(void)
    *函数形参:      Martix（矩阵指针）
                    G_Stack（伪递归栈堆指针）
                    Mark_x，Mark_y（GPA优化的断层标记数组指针）
    *函数功能:      销毁
    *返回值：       无
    */
    free(Martix[0]);//去掉集体分配的首地址就可以了
    free(Martix);
    free(G_Stack);//销毁递归栈堆
    free(Mark_x); free(Mark_y);//销毁点集
}

void If_Cls(void)
{
    /*函数名:       If_Cls
    *调用者函数:    void main(void)
    *函数形参:      无
    *函数功能:      询问是否清屏
    *返回值：       无
    */
    char choice;
    printf("Do you want to clean the screen?(Y or N)");
    while (1)
    {
        choice = getchar();
        if (choice == 'Y')
        {
            system("cls");
            return;
        }
        else if (choice == 'N')
            return;
        else
        {
            puts("Error Input!Please enter again!\n");
            fflush(stdin);
        }
    }
}

void Show(const int k, const int link_group_sum, const int loop, const double max_poss, FILE *fp)
{
    /*函数名:       Show
    *调用者函数:    void main(void)
    *函数形参:      k（以第k个概率产生点）
                    link_group_sum（联通集团个数）
                    loop（循环次数）
                    max（联通最大点个数）
                    p_sum（所有联通集团的粒子总数）
                    fp（指向data.txt的文件指针）
    *函数功能:      显示当前概率集团信息，并且把数据写入data.txt中
    *返回值：       无
    */
    puts("===============================================================================\n");
    printf("\a%cEcho %d:\nPossibility of %.4f\nThe linked graphs` sum is %d\nAccount for %.2f%% in %d graghs\n",
        16, k + 1, poss[k], link_group_sum, 100 * (double)link_group_sum / (double)loop, loop);
    printf("The possibility of(max linked points/sum points)accounts for %.2f%%\n", 100 * max_poss);
    printf("\n");

    if (fp != NULL)
    {
        fprintf(fp, "===============================================================================\n");
        fprintf(fp, "%cEcho %d:\nPossibility of %.4f\nThe linked graphs` sum is %d\nAccount for %.2f%% in %d graghs\n",
            7, k + 1, poss[k], link_group_sum, 100 * (double)link_group_sum / (double)loop, loop);
        fprintf(fp, "The possibility of(max linked points/sum points)accounts for %.2f%%\n", 100 * max_poss);
        fprintf(fp, "\n");
    }
}

#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)\n(The program will stop until you input zero)\n", 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\n");if ((fp = fopen("data.txt", "r")) == NULL)puts("""data.txt"" doesn`t exsit,the program will create a new one\n");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\n\n", 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;
}