C++中退火网络的SIS模型

退火网络介绍

退火网络与静态网络不同，由于退火网络的网络演化速率远高于动力学的演化速率，因此在每一个时刻，网络中各节点都在重新与其它节点进行连边，即每个感染或恢复的事件发生前，网络中各节点的邻居节点都在发生变化。这意味着退火网络并没有实质的网络连接关系，通常需要在动力学演化中过程中不断地进行生成。

因此，退火网络在连续时间SIS动力学中进行演化时，在每一个事件发生前，节点都应当重新选择其他节点作为新的邻居节点。

C++代码实现

举个例子，度分布满足无标度无关联网络的退火网络的SIS（易感-感染-易感）动力学过程，C++代码实现如下。

#include <iostream>
#include <ctime>
#include <stdlib.h>
#include <fstream>
#include <cmath>
#include "mt19937ar.c"
#include<algorithm>
#include <vector>

using namespace std;

int main()
{
    int Network_Size = 5000;  //网络尺度
    double r = 0.4;  //感染速率
    double g = 0.5;  //恢复速率
    int I0 = 100;    //初始感染种子
    double tend = 100;   //动力学演化最长时间
    double counter = 30;   //动力学模拟次数
    double wt = 0;
    double dt = 0;
    double t = 0;
    int I = 0;
    double nt = 0;
    int n_intervals = tend*10;   
    double Sum_degree = 0;
    double interval = tend/(1.0*n_intervals);  //每多少间隔记录一次数据
    vector<int> degree(Network_Size);
    vector<vector<int> > adjlist(Network_Size);
    vector<int> X_vec_ind(n_intervals + 1,0);
    vector<vector<int> > X_vec;
    unsigned long idum;
    idum=(unsigned long)time(NULL);
    init_genrand(idum);
    ofstream outfile_C_ext;
    outfile_C_ext.open("1.txt");  //输出文本
    ifstream infile;
    infile.open("UCM_degree.txt");   //导入网络中各节点的度，本代码采用的无标度网络度分布
    for(int i = 0; i < Network_Size; i++)
    {
        infile >> degree[i];
        Sum_degree = Sum_degree + degree[i];
    }

    for(int i = 0; i < counter; i++)   //动力学演化开始
    {
        double count1 = 1;
        t = 0;
        wt = 0;
        I = I0;
        nt = 0;
        X_vec_ind[0] = I0;
        vector<int> I_list;
        vector<int> I_node(Network_Size,0);
        vector<double> ni(Network_Size,0);
        vector<double> probability;
        for(;;)
        {
            if(I == 0)
            {
                break;
            }
            int random_1 = genrand_int32()%Network_Size;
            if(I_node[random_1] == 0)
            {
                I_node[random_1] = 1;
                I_list.push_back(random_1);
                I--;
            }
            else
                continue;
        }
        I = I0;
        for(int i = 0; i < I_list.size(); i++)
        {
            ni[I_list[i]] = g + r*degree[I_list[i]];;
            wt = wt + ni[I_list[i]];
        }
        dt = 1.0/wt;

        for(;;)
        {
            vector<double> selection_probability(I_list.size(),0);
            t = t + dt;
            while((t > count1*interval)&&(count1 <= n_intervals))
            {
                X_vec_ind[count1] = I;
                count1 ++;
            }

            if(t >= tend)
            {
                break;
            }
            double temp4 = 0;
            for(int i = 0; i < I_list.size(); i++)
            {
                temp4 = temp4 + ni[I_list[i]]/wt;
                selection_probability[i] = temp4;
            }
            double u = genrand_real2();
            for(int i = 0; i < I_list.size(); i++)
            {
                if( u <= selection_probability[i] )
                {
                    double v = genrand_real2();
                    if(v < g/ni[I_list[i]])
                    {
                        I_node[I_list[i]] = 0;
                        ni[I_list[i]] = 0;
                        double tempx;
                        tempx = I_list[i];
                        I_list[i] = I_list.back();
                        I_list.back() = tempx;
                        I_list.pop_back();
                        I--;
                        break;
                    }
                    else
                    {
                        double f1 = genrand_real2();
                        double temp5 = 0;
                        double temp6 = Sum_degree - degree[I_list[i]];
                        for(int j = 0; j < Network_Size; j++)
                        {
                            if(j == I_list[i])
                            {
                                continue;
                            }
                            temp5 = temp5 + degree[j];
                            probability.push_back(temp5/temp6);
                        }
                        int x = 0;
                        for(int k = 0; k < Network_Size; k++)
                        {
                            if(k == I_list[i])
                            {
                                continue;
                            }
                            if( f1 <= probability[x])
                            {
                                if(I_node[k] == 0)
                                {
                                    I_node[k] = 1;
                                    I_list.push_back(k);
                                    ni[k] = g + r*degree[k];
                                    I++;
                                    break;
                                }
                                else
                                    break;
                            }
                            x++;
                        }
                        probability.clear();
                        break;
                    }
                }
            }
            if(I == 0)
            {
                while (count1 <= n_intervals)
                {
                    X_vec_ind[count1] = 0;
                    count1 ++;
                }
                break;
            }
            wt = 0;
            for(int i = 0; i < I_list.size(); i++)
            {
                wt = wt + ni[I_list[i]];
            }
            dt = 1.0/wt;
        }
        X_vec.push_back(X_vec_ind);
    }
    double EXt;   //对所有模拟结果取平均
    for(int i=0; i <= n_intervals; i++)
    {
        EXt = 0;
        for(int j=0; j < counter; j++)
        {
            EXt += X_vec[j][i];
        }
        EXt /= 1.0*counter*Network_Size;
        outfile_C_ext << EXt << "\n";
    }
    outfile_C_ext << endl;




    return 0;
}