多层网络的SIR实现

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<malloc.h>
#include<math.h>
#include<time.h>
#define MaxVertexNum 90000
#define RAND_MAX 0x7fff 
//边表节点
typedef struct node
{
int adjvex;
struct node *next;
int visit;
}EdgeNode;
//顶点表节点
typedef struct vnode
{
int vertex;
int KS;//k-core
float ec;//特征向量中心性
int is_infected;
int nodeState;
int is_visit;
int layer;//层
int degree;//该节点的度
EdgeNode* firstedge;//指向边表的指针
}VertexNode;
typedef VertexNode AdjList[MaxVertexNum];
//图的结构
typedef struct
{
AdjList adjlist;//顶点表（类似于数组的）
int vexnum;//顶点个数
int arcnum;//边个数
}ALGraph;
//返回文件行数（网络边数）,有换行符"\n"就为一行
int lines(char* str)
{
int c;
FILE* fp;
int lines=0;
fp=fopen(str,"r");
if(fp)
    {
    while((c=fgetc(fp))!=EOF)
        if(c=='\n')
            lines++;
        fclose(fp);
    }
return lines;
}
//返回文件最大数（网络节点数）
int max(char* str)
{
FILE* fp;
char* p;
int line=lines(str);
int i=0;
int a=0;
int b=0;
fp=fopen(str,"r");
char buf[1024];
if((fp=fopen(str,"r"))==NULL)
    {
    perror("fail to read");
        exit(1);
    }
//把文件的内容给buf

while(fgets(buf,line,fp)!=NULL)
    {


//p就没有    
    p=buf;
sscanf(p,"%d %d",&a,&b);//输入源为p

    //i始终为最大的
    if(a>i)
        i=a;
    if(b>i)
        i=b;
    }

return i;
}
//创建图
void createAlgraph(ALGraph*  G,char* str)
{
FILE* fp;
int line=lines(str);
int node=max(str);//其中最大数
G->vexnum=node+1;//因为是从0开始，所以+1，多了一个0 
G->arcnum=line;
fp=fopen(str,"r");
char buf[1024];
int len;
int m;
int n;
EdgeNode* s;
char* p;
int a=0;
int b=0;
int i=0;
char StrLine[1024];
//每个节点的顶点表（vn1(i),vn2(i)）
int vn1[line];//这里本来要用vn[line],如果是vc++就不能通过编译,有多少行就有多少（i,j）
int vn2[line];
//顶点录入
for(int j=0;j<G->vexnum;j++)
    {
    G->adjlist[j].vertex=j;
    G->adjlist[j].firstedge=NULL;
    G->adjlist[j].nodeState=1;
    }
if((fp=fopen(str,"r"))==NULL)
    {
    perror("faile to read");
    exit(1);
    }

while(!feof(fp))//因为行数等于边数，则读取行数个就可以把其他的节点的连接读完
    {

fgets(StrLine,1024,fp);
sscanf(StrLine,"%d%d",&a,&b);

    vn1[i]=a;
    vn2[i]=b;
    i++;
    }

//边节点放入链表
//一行就是一个坐标，有多少行就有多少坐标 
for(int k=0;k<line;k++)//有多少行就有多少节点，每个节点对应一个边链
    {
    m=vn1[k],n=vn2[k];
    int ii=0;
    EdgeNode* p;
    p=G->adjlist[m].firstedge;
    while(p!=NULL)
        {
        if(p->adjvex==n)
            {
            ii=1;
            break;
            }
        p=p->next;
        }
    if(ii!=1)
        {
        s=(EdgeNode*)malloc(sizeof(EdgeNode));
        s->adjvex=n;//相连接的顶点
        s->next=NULL;
        s->next=G->adjlist[m].firstedge;//类似于自己写的链表 
        G->adjlist[m].firstedge=s;
        //无向图 有来回
        s=(EdgeNode*)malloc(sizeof(EdgeNode));
        s->adjvex=m;
        s->next=NULL;
        s->next=G->adjlist[n].firstedge;
        G->adjlist[n].firstedge=s;
        }
    }
//深度为每个节点后面连接的链长度
EdgeNode* q;
for( i=0;i<G->vexnum;i++)
    {
    int k=0;
    q=G->adjlist[i].firstedge;
    while(q!=NULL)
        {
        k++;
        q=q->next;
        }
    G->adjlist[i].degree=k;    
    }
    
}
//打印邻接表
void printGraph(ALGraph* G)
{
        EdgeNode* s;
    for(int i=0;i<G->vexnum;i++)
        {
        
        s=G->adjlist[i].firstedge;//s为一个带adjvex，next指针的边表节点 
    
        while(s)
            {
            printf("(%d,%d)",G->adjlist[i].vertex,s->adjvex);    
            s=s->next;
            }
        printf("\n");
        }
}
//所属层插入
void insertLayer(ALGraph* G,int layer)
{
for(int i=0;i<G->vexnum;i++)
    {
    G->adjlist[i].layer=layer;
    }
}
//打印度中心性
void printDegreeCentrality(ALGraph* G)
{
    
for(int i=0;i<G->vexnum;i++)
    {
    
    printf("node %d dgree centrality is:%d\n",i,G->adjlist[i].degree);
    
    }
}


//计算特征向量中心性
void eigenvector_centrality(ALGraph *G)
{
    float e[G->vexnum];//记录上一次的指标（最终的特征向量中心性指标 ，因为会把最终的计算赋值给e）；下面都用指标代表特征向量指标
    float e1[G->vexnum];//记录这一次的指标 
    float max = 0;//这一次的最大指标 
    float max1 = 0;//记录上一次最大指标 
    int  flag=0;//当flag=1时，代表找到各个指标 
    for(int i=0; i<G->vexnum; i++)
    {
        e[i]=1;//将每个点初始化为1 
        e1[i]=0;
    }
    EdgeNode *p;
    p=(EdgeNode*)malloc(sizeof(EdgeNode));
    //循环开始 
    while(flag==0)
    {
        max1=max;//max1为上一次的最大值 
        max=0;
        for (int i=0; i<G->vexnum; i++)
        {
            p=G->adjlist[i].firstedge;
            while(p!=NULL)
            {
                e1[i]+=e[p->adjvex];//第一次的计算结果为他们各自的度 
                p=p->next;
            }
            if(e1[i]>max)
                max=e1[i];//记录本次的最大指标 
        }
        for(int i=0; i<G->vexnum; i++)
        {
            if(e[i]!=e1[i])
                break;
            if(i==G->vexnum-1)
                flag=1;//两次计算结果相同结束循环 
        }
        if((1.0/max1-1.0/max)<0.01&&(1.0/max1-1.0/max)>-0.01)
            flag=1;//当差值较小时也可结束循环 
            //保留这次的结果到e中，并且将ei重置为0，方便下次计算 
        for(int i=0; i<G->vexnum; i++)
        {
            e[i]=e1[i]; 
            e1[i]=0;
        }
    }
    for(int i=0; i<G->vexnum; i++)
    {
        e[i]=e[i]/max;
        G->adjlist[i].ec=e[i];
    }
}
/*
1.SIR传播模型，改变node.state(0为感染，1为易感染，2为恢复。所有节点初始化为1) 
2.SIR函数主要操作感染点，感染点会做两件事：①感染易感节点；②感染节点恢复 
3.传播完成的标志为不能再感染新的节点，并返回处于恢复节点的个数 
*/
void SIR1(ALGraph* G,int a,double inf,double rec,char* str)
//传入的分别为网络,感染节点,感染率，恢复率,写入的文件 
{
    double rate;//传入节点作为感染节点的感染规模 
    
    int i=0;
    
    int sum=0;//统计100次传播后的结果 
    FILE* fp;
    
    fp=fopen(str,"at"); 
    //所有节点最开始为易感态 
    for(int k=0;k<100;k++)
    {    
        int Inf[G->vexnum];
        int newInf[G->vexnum];
        int recnum=0;//统计传播结束后，处于恢复状态的节点个数
        for(i=0;i<G->vexnum;i++)
        {
            G->adjlist[i].nodeState=0;
        } 
//给感染节点赋值 
    G->adjlist[a].nodeState=1;//传入的j节点a为感染态 
    //把感染节点交给数组 
    Inf[0]=a;
    double infection=inf;//感染概率 
    int count=1;//当前网络中的感染个数
    srand((unsigned)time(NULL)); //设置种子，用于随机数产生 
    while(count>0)//还能继续感染 
    {
        int newInfLength=0;//表示新感染节点的个数 
        //节点感染 
        for(i=0;i<count;i++)
        {
            int vert=Inf[i];//当前的感染点 
            
            EdgeNode* p;
            
            p=G->adjlist[vert].firstedge;
            //用当前节点去感染其他节点
            while(p!=NULL)
            {
                int n=0; 
                double infect_rate;//感染的概率为1-（1-λ）^n;其中λ为感染率，n为周围节点是感染者的个数 
                double test=rand()/(double)RAND_MAX;//rand()产生随机数为[1,32767],RAND_MAX设置为32767，那么test范围[0.1] ;
                 //计算n
                int nodej=p->adjvex;//记录当前连接的节点
                //用s查看当前连接节点的周围有多少感染者 
                EdgeNode* s=G->adjlist[nodej].firstedge;
                while(s!=NULL)
                {
                    if(G->adjlist[s->adjvex].nodeState==1)
                    {
                        n++;
                    }
                    s=s->next;
                }
                 //计算感染率infect_rate
                 infect_rate=1.0-pow(1.0-infection,n);
                
                //如果随机数比感染概率小（能感染），且节点状态为易感染，就感染该节点 
                if(test<=infect_rate&&G->adjlist[nodej].nodeState==0)
                {
                newInf[newInfLength]=nodej;
                G->adjlist[nodej].nodeState=1;//被感染 
                newInfLength++;    
                }
                p=p->next; 
            } 
        }
        //感染节点恢复(不包括上一步新感染的)
         for(i=0;i<count;i++)
         {
             double recovRate=rec;
             double test_1=rand()/(double)RAND_MAX;//rand()产生随机数为[1,32767],RAND_MAX设置为32767，那么test范围[0.1] 
             //此处当恢复率设置为1时所有感染节点都能恢复 
             //恢复分两种情况：1.能恢复,改变nodeState为2；2.不能恢复，放入新感染数组 
             if(test_1<=recovRate)
             {
                 G->adjlist[Inf[i]].nodeState=2;
              } 
              else
              {
                  newInf[newInfLength]=Inf[i];
                  newInfLength++;
              }
         }
         //newInf数组中元素两个来源：1.易感染节点被感染；2.感染节点未恢复 
         //再把新感染的数组newInf交给Inf进行下一次循环
         for(i=0;i<newInfLength;i++)
         {
             Inf[i]=newInf[i];
             
          } 
         
          count=newInfLength;//记录当前新感染的个数，作为继续循环的依据 
    }
    for(i=0;i<G->vexnum;i++)
    {
        if(G->adjlist[i].nodeState==2)
        {
        
        recnum++;    
        }
     }
     sum+=recnum;
//     //重置所有会节点状态 ，作为下一轮感染准备 
//     recnum=0;
////     Inf[G->vexnum]={0};
////    newInf[G->vexnum]={0};
//    for(i=0;i<G->vexnum;i++)
//    {
//        G->adjlist[i].nodeState=;
//     } 
         } 
      rate=(sum*1.0)/(G->vexnum*100); 
    fprintf(fp,"%d    %lf     %lf\n",G->adjlist[a].degree,G->adjlist[a].ec,rate);
    //重置所有节点的node.state
    
    fclose(fp);
    return;
    
}

void SIR(ALGraph *G,char* str,double beta)
{
    //0:易感染 1:感染 2：恢复 
    double gama=1.0;
    srand(time( NULL ));
    int sum_of_recovered;
    int sum_of_infected;
    FILE *fp;
    fp=fopen(str,"w");
    for(int node=0;node<G->vexnum;node++)
    {
        printf("%d  ",node);
        int sum=0;
        for(int k=0; k<100; k++)
        {
            sum_of_recovered =0;
            //先把所有节点设置为易感染 
            for(int i=0; i<G->vexnum; i++)
                G->adjlist[i].is_infected=0;
                //最初那个节点设置为感染 
            G->adjlist[node].is_infected=1;
            sum_of_infected=1;
            //下面的while中一轮一轮的把所有节点进行感染 
            while(1)
            {
                //记录所有节点的感染状态 
                int temp[G->vexnum];
                for(int i=0; i<G->vexnum; i++)
                    temp[i]=G->adjlist[i].is_infected;
                    //下面把所有节点进行判断是否感染或易感 
                for(int i=0; i<G->vexnum; i++)
                {
                    if(G->adjlist[i].is_infected==0)
                    {
                        EdgeNode *p;
                        p=G->adjlist[i].firstedge;
                        int count=0;//统计易感旁边感染的个数，用于计算感染的概率 
                        while(p!=NULL)
                        {
                            if(temp[p->adjvex]==1)
                                count++;
                            p=p->next;
                        }
                        double infect_rate=1-pow(1.0-beta,count);
                        if(rand()/(double)RAND_MAX<=infect_rate)
                        {
                            G->adjlist[i].is_infected=1;
                            sum_of_infected++;
                        }
                    }
                    else if(G->adjlist[i].is_infected==1)
                    {
                        if(rand()/(double)RAND_MAX<=gama)
                        {
                            G->adjlist[i].is_infected=2;
                            sum_of_infected--;
                            sum_of_recovered++;
                        }
                    }
                }
                if(sum_of_infected==0)
                    break;
            }
            sum+=sum_of_recovered;
        }
        double rate=(sum*1.0)/(G->vexnum*100.0);
        fprintf(fp,"%d    %lf     %lf\n",G->adjlist[node].degree,G->adjlist[node].ec,rate);
    }
    fclose(fp);
}

void DoubleNetworkSIR(ALGraph *G1,ALGraph* G2,ALGraph* G3,char* str,double beta1,double beta2) 
{    
    //传入参数依次为（网络G1.连接边G2，网络G3，写入的文件，G1中的传播率，G2中的传播率） 
    //0:易感染 1:感染 2：恢复 
    double gama=1.0;
    srand(time( NULL ));
    int sum_of_recovered;
    int sum_of_infected;
    FILE *fp;
    fp=fopen(str,"w");
    for(int node=0;node<(G1->vexnum+G3->vexnum);node++)
    {
        printf("%d  ",node);
        int sum=0;
        for(int k=0; k<100; k++)
        {
            sum_of_recovered =0;
            //先把所有节点设置为易感染 
            for(int i=0; i<G1->vexnum; i++)
                G1->adjlist[i].is_infected=0;
            for(int i=0; i<G3->vexnum; i++)
                G3->adjlist[i].is_infected=0;
                //最初那个节点设置为感染(需要区分两个不同的网络)
            if(node<G1->vexnum) 
            G1->adjlist[node].is_infected=1;
            else
            G3->adjlist[node-(G1->vexnum)].is_infected=1;
            sum_of_infected=1;
            //下面的while中一轮一轮的把所有节点进行感染 
            while(1)
            {
                
                
                //记录所有节点的感染状态 
                int temp1[G1->vexnum];//G1中的感染点 
                int temp2[G3->vexnum];//G3中的感染点 
                for(int i=0; i<G1->vexnum; i++)
                    temp1[i]=G1->adjlist[i].is_infected;
                for(int i=0; i<G3->vexnum; i++)
                    temp2[i]=G3->adjlist[i].is_infected;
                    //下面把所有节点进行判断是否感染或易感 
                    //G1中的点 
                for(int i=0; i<G1->vexnum; i++)
                {
                    if(G1->adjlist[i].is_infected==0)
                    {
                        EdgeNode *p;
                        p=G1->adjlist[i].firstedge;
                        int count1=0;//统计易感旁边感染的个数，用于计算感染的概率 （先统计自己网络的，再统计相连接网络的） ,这里为G1中的邻居 
                        int count2=0;//统计G3中的邻居 
                        while(p!=NULL)
                        {
                            if(temp1[p->adjvex]==1)
                                count1++;
                            p=p->next;
                        }
                        //需要判断这个点是否与G3中的点相连（根据这个点是否在G2判断） 
                         EdgeNode *s;
                         s=G2->adjlist[i].firstedge;
                         //如果相连，看相连G3中的点是否感染 
                         while(s!=NULL)
                         {
                             if(temp2[s->adjvex]==1)
                             count2++;
                             s=s->next;
                         }
                        
                        double infect_rate=1-pow(1.0-beta1,count1)*pow(1.0-beta2,count2);
                        if(rand()/(double)RAND_MAX<=infect_rate)
                        {
                            G1->adjlist[i].is_infected=1;
                            sum_of_infected++;
                        }
                    }
                    else if(G1->adjlist[i].is_infected==1)
                    {
                        if(rand()/(double)RAND_MAX<=gama)
                        {
                            G1->adjlist[i].is_infected=2;
                            sum_of_infected--;
                            sum_of_recovered++;
                        }
                    }
                }
                //G3中的点 
                for(int i=0; i<G3->vexnum; i++)
                {
                    if(G3->adjlist[i].is_infected==0)
                    {
                        EdgeNode *p;
                        p=G3->adjlist[i].firstedge;
                         int count1=0;//统计易感旁边感染的个数，用于计算感染的概率 
                         int count2=0;//统计G3的邻居 
                        while(p!=NULL)
                        {
                            if(temp2[p->adjvex]==1)
                                count2++;
                            p=p->next;
                        }
                        //需要判断这个点是否与G3中的点相连（根据这个点是否在G2判断） 
                         EdgeNode *s;
                         s=G2->adjlist[i].firstedge;
                         //如果相连，看相连G3中的点是否感染 
                         while(s!=NULL)
                         {
                             if(temp1[s->adjvex]==1)
                             count1++;
                             s=s->next;
                         }
                        
                        double infect_rate=1-pow(1.0-beta1,count1)*pow(1.0-beta2,count2);
                        if(rand()/(double)RAND_MAX<=infect_rate)
                        {
                            G3->adjlist[i].is_infected=1;
                            sum_of_infected++;
                        }
                    }
                    else if(G3->adjlist[i].is_infected==1)
                    {
                        if(rand()/(double)RAND_MAX<=gama)
                        {
                            G3->adjlist[i].is_infected=2;
                            sum_of_infected--;
                            sum_of_recovered++;
                        }
                    }
                }
                
                if(sum_of_infected==0)
                    break;
            }
            sum+=sum_of_recovered;
        }
        double rate=(sum*1.0)/((G1->vexnum+G3->vexnum)*100.0);
        if(node<G1->vexnum)
        fprintf(fp,"%d    %lf     %lf\n",G1->adjlist[node].degree,G1->adjlist[node].ec,rate);
        else
        fprintf(fp,"%d    %lf     %lf\n",G3->adjlist[node-(G1->vexnum)].degree,G3->adjlist[node-(G1->vexnum)].ec,rate);
    }
    fclose(fp);    
}

//创建连接的边，并写入文件
void createLink(char* str)
{
    FILE* fp;
    fp=fopen(str,"w");
    srand((unsigned)time(NULL));
    for (int i=0;i<2000;i++)
    {
        fprintf(fp,"%d    %d\n",rand()%2000,rand()%2000); 
    //fpintf(fp,"%d    %d",rand()%2000,rand()%2000);
    }
    fclose(fp);
 } 

int main()
{

char* str1="E:\\data_set\\netsci1.txt";//G1，网络 
char* str2="E:\\data_set\\netsci2.txt";//G2，连边 
char* str3="E:\\data_set\\netsci3.txt";//G3,网络 
 createLink(str2);//创建连边 

//G1、G3为两个网络,G2为他们的连接
ALGraph* G1;
ALGraph* G2;
ALGraph* G3;
G1=(ALGraph*)malloc(sizeof(ALGraph));
G2=(ALGraph*)malloc(sizeof(ALGraph));
G3=(ALGraph*)malloc(sizeof(ALGraph));
//创建三个表的信息 
createAlgraph(G1,str1);//分别插入图的地址，连接顶点信息
createAlgraph(G2,str2);
createAlgraph(G3,str3);
//插入层数
insertLayer(G1,1);
insertLayer(G3,2);
//计算特征向量中心性 
eigenvector_centrality(G1);
eigenvector_centrality(G3);
//SIR，此处为单个节点推动 
char* str4="E:\\data_set\\result.txt";
//{
//printf("%d ",i);    
//SIR1(G1,i,0.08,1,str4);    
//}
//SIR(G1,str4,0.08); 
DoubleNetworkSIR(G1,G2,G3,str4,0.15,0.1);


    return 0;
}