std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<> dis(0, N-1);

No. of rec = 99
No. of rec = 101
No. of rec = 101
No. of rec = 103
No. of rec = 102

No. of rec =  98
No. of rec =  309
No. of rec =  1433
No. of rec =  3730
No. of rec =  6590

import networkx as nx
import numpy as np
import matplotlib.pyplot as plt
import random

N = 10000
m = 2
epsilon = 0.001
eta = 1.
gamma = -2.1


def power_law(x0, x1, gamma):
    pl = []
    for i in range (N):
        pl.append(((x1**(gamma+1) - x0**(gamma+1))*random.uniform(0,1)  + x0**(gamma+1.0))**(1.0/(gamma + 1.0)))
    return pl


def graph_activity(act,m,N,G):

    active_nodes=[]

    c=0
    for i in act:
        if random.uniform(0,1)<i:
            active_nodes.append(c)
        G.add_node(c)
        c+=1

    for i in active_nodes:

        count=0
        while count<m:
            target=random.randint(0,N-1)
            if target!=i and target not in G.neighbors(i):
                    G.add_edge(i,target)
                    count+=1

    return 0



def SIR(status,l,mu,G):

    list_inf=set()
    list_rec=set()


    for i in range(status.size):
        if status[i]==1:
            for j in G.neighbors(i):
                if status[j]==0:
                    if random.uniform(0,1)<l:
                        list_inf.add(j)


            if random.uniform(0,1)<mu:
                list_rec.add(i)

    for i in list_inf:
        status[i]=1
    for i in list_rec:
        status[i]=2


def seeds(N,status,p):
    number=int(N*p)

    list_nodes=[]
    for i in range(N):
        list_nodes.append(i)

    random.shuffle(list_nodes)


    for i in range(number):
        status[list_nodes[i]]=1



list_l=list(np.arange(0,1,0.2))

recovered=[]
for l in list_l:
    act=power_law(0.001,1,gamma)
    status=np.zeros(N,int)
    seeds(N,status,0.01)

    for i in range(500):
        G=nx.Graph()
        graph_activity(act,m,N,G)
        SIR(status, l, 0.01, G)
        #print (2*G.number_of_edges()/G.number_of_nodes(), 2*m*average_a)


    count=0
    for j in range(N):
            if status[j]==2:
                count+=1

    print ("No. of rec = ", count)

    recovered.append(count)

#include <iostream>
#include <igraph/igraph.h>
#include <vector>
#include <random>
#include <ctime>
using namespace std;

int N = 10000;
int m = 2;
double epsilon = 0.001;
double eta = 1.0;
double g = -2.1;

double random_number(){
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<> dis(0, 1);//uniform distribution between 0 and 1
    return dis(gen); // random number between 0 and 1
}

vector<double> power_law(double x0, double x1, double _gamma){
    vector<double> pl;
    for (int i = 0; i < N; i++){
        double r = random_number();
        pl.push_back( pow((((pow(x1, g+1)) - (pow(x0,g+1)))*r + (pow(x0, g+1))),(1/(g+1))));
    }
    return pl;
}

vector<int> neighbours(igraph_vector_t *v) {
    vector<int> neigh;
    long int i;
    for (i=0; i<igraph_vector_size(v); i++) {
        neigh.push_back(VECTOR(*v)[i]);
    }
    return neigh;
}

igraph_t activate_graph(vector<double> activity, igraph_t g){
    vector<int> active_nodes;
    igraph_vector_t v;
    igraph_vector_init(&v, 8);

    for (int i = 0; i < activity.size(); i++){
        if (random_number() < activity[i]){
            active_nodes.push_back(i);
        }
    }

    for (int i = 0; i < active_nodes.size(); i++){
        int count = 0;
        while (count < m){
            int target = -1;
            while (target < 0) {
                std::random_device rd;
                std::mt19937 gen(rd());
                std::uniform_real_distribution<> dis(0, N-1);//uniform distribution between 0 and 1
                target = dis(gen);
            }
                igraph_neighbors(&g, &v, active_nodes[i], IGRAPH_ALL);
                vector<int> n = neighbours(&v);
                if (target != active_nodes[i] && !(std::find(n.begin(), n.end(), active_nodes[i]) != n.end())) {
                    igraph_add_edge(&g, active_nodes[i], target);
                    count++;
                }
        }
    }
    return g;
}

vector<int> seeds(int N, double p, vector<int> status){
    int number = (int)(N*p);
    vector<int> list_nodes;
    for (int i = 1; i < N-1; i++){
        list_nodes.push_back(i);
    }
    std::random_device rd; // Warning: implementation of this is allowed to be stupid and
    // return the same value every time. Does not work in mingw 4.8
    // Can't speak for 4.9 or greater
    std::mt19937 randomizer(rd());

    shuffle ( list_nodes.begin(), list_nodes.end(), randomizer); // shuffle the vector
    list_nodes.resize(number);

    for (int i = 0; i < list_nodes.size(); i++){
        status[list_nodes[i]] = 1;
    }

    return status;
}

vector<int> sir(vector<int> status_list, igraph_t g, double Lambda){
    set<int> list_inf;
    set<int> list_rec;

    igraph_vector_t v;
    igraph_vector_init(&v, 8);

    double mu = 0.01;
    for (int i = 0; i < status_list.size(); i++) {
            if (status_list[i] == 1) {
                igraph_neighbors(&g, &v, i, IGRAPH_ALL);
                vector<int> neigh = neighbours(&v);
                for (int j = 0; j < neigh.size(); j++) {
                    if (status_list[neigh[j]] == 0) {
                        if (random_number() < Lambda) {
                            list_inf.insert(j);
                        }
                    }
                }
                if (random_number() < mu) {
                    list_rec.insert(i);
                }
            }

        }
        srand((NULL));
        vector<int> vec_inf;
        vector<int> vec_rec;

        std::copy(list_inf.begin(), list_inf.end(), std::back_inserter(vec_inf));
        std::copy(list_rec.begin(), list_rec.end(), std::back_inserter(vec_rec));

        for (int i = 0; i < vec_inf.size(); i++) {
            status_list[vec_inf[i]] = 1;
        }
        for (int i = 0; i < vec_rec.size(); i++) {
            status_list[vec_rec[i]] = 2;
        }
    return status_list;
}

int main() {
    static vector<int> status_list;
    double Lambda = 0;
    vector<double> recovered;
    igraph_t graph;
    while (Lambda < 1.0) {

        vector<double> act = power_law(epsilon, 1, g);
        status_list = vector<int>(N, 0);
        status_list = seeds(N, 0.01, status_list);

        for (int i = 0; i < 500; i++) {
            igraph_empty(&graph, N, false);
            graph = activate_graph(act, graph);
            status_list = sir(status_list, graph, Lambda);
        }
        int count = 0;
        for (int j = 0; j < N; j++) {
            if (status_list[j] == 2) {
                count++;
            }
        }

        recovered.push_back(count);
        cout << "No. of rec = " << count << endl;
        igraph_destroy(&graph);
        Lambda += 0.2;
    }
    return 0;
}