传播模型（SIR）

参考SIR模型原理：

①https://blog.csdn.net/Feng512275/article/details/82859526/

②https://zhuanlan.zhihu.com/p/104072104?app=zhihulite

代码参考：

https://github.com/pholme/sir

参考章节：

一.互联网络模型

构造了两种互连子网。一个是通过随机或优先连接两个相同的子网络形成的，包括scale-free-scale网络和e-mail-e-mail网络。这种互联网络可以用来表示现实世界中连接不同社区网络所形成的网络。互连密度是用参数γ来测量的，定义为γ=L/N。L表示互连的个数，N表示一个子网的大小。我们构造的另一种互连网络是将一个网络随机分成两个大小相同的互连子网。这种互联网络的结构表明，包含不同类型节点的网络被划分为不同的互联子网。每个子网包含相同类型的节点。我们使用的单一网络包括友谊网络和自动系统网络。

二.传染病传播模型

传染病传播模型本文采用的传染病传播模型是最基本、研究最充分的传染病传播模型。网络的元素可以分为三个部分包括易感者（感染者）、感染者（感染者）和康复者（康复者）。在每个时间步，如果易受感染的节点直接连接到一个受感染的节点，则易受感染的节点被感染的概率为λ。参数λ称为扩展率。同时，受感染的节点会出现被移除的节点，其概率δ称为恢复率。但是对于两个互连网络（A-B网络），我们需要分别指定这些过程，表示λ^A（λ^B）网络中节点之间的扩展速率，λ^AB（λ^BA）网络中节点A（B）到节点B（A）的扩展速率.δ^A（δ^B）表示网络中节点之间的恢复速率。在我们的研究中，在没有失去一般性的情况下，我们让δ^A=δ^B=1，但是我们需要使用相对小的数值作为传播率（λ^A，λ^B，λ^AB，λ^BA）和连接密度γ。因此感染率仍然很小，如果传播可以达到人口的很大一部分，那么单个节点的作用就不再重要，传播将覆盖几乎所有的网络，而与网络的起源无关

三.编码思路：

总共有S（0）、 I（1）、 R（2）三类节点。初始化所有节点为I即nodeState为1，然后SIR函数中处理I（1）节点，分别为感染其他节点和恢复

int SIR(ALGraph* G,int a[],int arrayLength,double inf,double rec)
//传入的分别为网络，感染节点数组,数组长度,感染率，恢复率 
{
    int infTatal=arrayLength;//感染节点总数 
    int Inf[G->vexnum];
    int newInf[G->vexnum];
    int i=0;
    //给感染数组赋初值 
    for(i=0;i<infTatal;i++)
    {
        Inf[i]=a[i];
    
        G->adjlist[i].nodeState=0;//传入的数组为感染态 
    }
    for(i=0;i<infTatal;i++)
    {
        newInf[i]=0;
    }
    double infection=inf;//感染概率 
    int count=infTatal;//当前网络中的感染个数
    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)
            {
                double test=rand()/(double)RAND_MAX;//rand()产生随机数为[1,32767],RAND_MAX设置为32767，那么test范围[0.1] 
                int nodej=p->adjvex;//记录当前连接的节点
                
                
            
                //如果随机数比感染概率小（能感染），且节点状态为易感染，就感染该节点 
                if(test<infection&&G->adjlist[nodej].nodeState==1)
                {
                newInf[newInfLength]=nodej;
                G->adjlist[nodej].nodeState=0;//被感染 
                newInfLength++;    
                }
                p=p->next; 
            } 
        }
        //感染节点恢复(不包括上一步新感染的)
         for(i=0;i<count;i++)
         {
             double recovRate=rec;
             double test_1=rand()/(double)RAND_MAX;
             //恢复分两种情况：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;//记录当前新感染的个数，作为继续循环的依据 
    }
    int recnum=0;//统计传播结束后，处于恢复状态的节点个数
    for(i=0;i<G->vexnum;i++)
    {
        if(G->adjlist[i].nodeState==2)
        {
        
        recnum++;    
        }
     } 
    return recnum;
}

加上之前的全部代码：（这里使用是一个感染节点数组作为传播源）

#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;
//返回文件行数（网络边数）,有换行符"
"就为一行
int lines(char* str)
{
int c;
FILE* fp;
int lines=0;
fp=fopen(str,"r");
if(fp)
    {
    while((c=fgetc(fp))!=EOF)
        if(c=='
')
            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[256];
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("
");
        }
}
//所属层插入
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
",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.传播完成的标志为不能再感染新的节点，并返回处于恢复节点的个数 
*/

int SIR(ALGraph* G,int a[],int arrayLength,double inf,double rec)
//传入的分别为网络，感染节点数组,数组长度,感染率，恢复率 
{
    int infTatal=arrayLength;//感染节点总数 
    int Inf[G->vexnum];
    int newInf[G->vexnum];
    int i=0;
    //给感染数组赋初值 
    for(i=0;i<infTatal;i++)
    {
        Inf[i]=a[i];
    
        G->adjlist[i].nodeState=0;//传入的数组为感染态 
    }
    for(i=0;i<infTatal;i++)
    {
        newInf[i]=0;
    }
    double infection=inf;//感染概率 
    int count=infTatal;//当前网络中的感染个数
    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)
            {
                double test=rand()/(double)RAND_MAX;//rand()产生随机数为[1,32767],RAND_MAX设置为32767，那么test范围[0.1] 
                int nodej=p->adjvex;//记录当前连接的节点
                
                
            
                //如果随机数比感染概率小（能感染），且节点状态为易感染，就感染该节点 
                if(test<infection&&G->adjlist[nodej].nodeState==1)
                {
                newInf[newInfLength]=nodej;
                G->adjlist[nodej].nodeState=0;//被感染 
                newInfLength++;    
                }
                p=p->next; 
            } 
        }
        //感染节点恢复(不包括上一步新感染的)
         for(i=0;i<count;i++)
         {
             double recovRate=rec;
             double test_1=rand()/(double)RAND_MAX;
             //恢复分两种情况：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;//记录当前新感染的个数，作为继续循环的依据 
    }
    int recnum=0;//统计传播结束后，处于恢复状态的节点个数
    for(i=0;i<G->vexnum;i++)
    {
        if(G->adjlist[i].nodeState==2)
        {
        
        recnum++;    
        }
     } 
    return recnum;
}


int main()
{
char* str1="E:\data_set\netsci1.txt";
char* str2="E:\data_set\netsci2.txt";
char* str3="E:\data_set\netsci3.txt";
//G1、G2为两个网络,G3为他们的连接
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(G2,2);
//计算特征向量中心性 
eigenvector_centrality(G1);
int a[1]={3};
int arrayLenth=sizeof(a)/sizeof(int);
//SIR
printf("recnum=%d",SIR(G1,a,arrayLenth,0.6,0.2));
    return 0;
}