ClassSimulation.R

library(Rcpp)
library(igraph)
library(compiler)

cppFunction('
LogicalMatrix thresh_diffusion(int times, LogicalVector infected, NumericVector thresholds, IntegerMatrix el, double friendshipeffect){
    int n = infected.size();
    int ninfected = 0;
    int nedges = el.nrow();
    int friendsinfected;
    int friendsuninfected;
    double infectedeffect;
    double uninfectedeffect;
    double totaleffect;
    LogicalVector newinfected(n);
    LogicalMatrix infectedout(n, (times + 1));

    for (int i = 0; i < n; i++){
        if (infected[i]){
        ninfected++;
        newinfected[i] = true;
        }
    else{
        newinfected[i] = false;
        }
        infectedout(i, 0) = infected[i];
        }
    for (int i = 1; i < (times + 1); i++){
        for (int j = 0; j < n; j++){
        infectedout(j, i) = infectedout(j, (i - 1));
            if (infectedout(j, i) == false){
                friendsinfected = 0;
            friendsuninfected = 0;
            for (int k = 0; k < nedges; k++){
                if (el(k, 0) == j){
                    if (newinfected[el(k, 1)]){
                    friendsinfected++;
                    }
                else {
                    friendsuninfected++;
                    }
                    }
                if (el(k, 1) == j){
                    if (newinfected[el(k, 0)]){
                friendsinfected++;
                    friendsinfected++;
                    }
                else {
                    friendsuninfected++;
                    }
                    }
                }
            infectedeffect = ninfected + ((friendshipeffect - 1) * friendsinfected);
            uninfectedeffect = n - ninfected + ((friendshipeffect - 1) * friendsuninfected);
            totaleffect = infectedeffect / (infectedeffect + uninfectedeffect);
            infectedout(j, i) = totaleffect > thresholds[j];
            if (infectedout(j, i)){
            newinfected[j] = true;
                ninfected++;
                }
            }
            }
        }
    return infectedout;
    }
')

# I have tried to write the function below in C++ instead, though it's sometimes much slower.
vazquez.game.el<-function(N, u){
  #_N_ is the number of nodes
  #_u_ is the triadic closure parameter. Needs to be between 0 and 1.
  #Alexei Vazquez.  2003.  "Growing network with local rules: Preferential attachment, clustering hierarchy, and degree correlations."  Physical Review E 67, 056104
  add.node.edge<-function()
    return(add.node=c(rbinom(1,1,1-u), close.tri=rbinom(1,1,u)))
  vneigh<-function(v, el){
    receivers<-el[which(el[,1]==v), 2]
    senders<-el[which(el[,2]==v), 1]
    return(unique(c(receivers, senders)))
    }
  g<-matrix(NA, 0, 2)
  potentials<-matrix(NA, 0, 2)
  ncount<-1
  while(ncount<N){
    a<-add.node.edge()
    if (a[1]==1){
      source.v<-ncount+1
      target.v<-ifelse(ncount>1, sample.int(ncount, 1), 1)
      target.neighbors<-vneigh(target.v, g)
      g<-rbind(g, c(source.v, target.v))
      n.neighbors<-length(target.neighbors)
      if (n.neighbors>0){
        new.potentials<-matrix(NA, n.neighbors, 2)
        new.potentials[,1]<-source.v
        new.potentials[,2]<-target.neighbors
        potentials<-rbind(potentials, new.potentials)
        rm(new.potentials)
        }
      rm(source.v, target.v, target.neighbors, n.neighbors)
      ncount<-ncount+1
      }
    if ((a[2]==1)&(ncount>2)&(length(potentials)>=2)){
      if (length(potentials)==2){
        g<-rbind(g, c(potentials[1], potentials[2]))
        potentials<-matrix(NA, 0, 2)
        }
      else{
        potentials.ind<-sample.int(nrow(potentials), 1)
        adders<-potentials[potentials.ind, ]
        g<-rbind(g, c(adders[1], adders[2]))
        potentials<-potentials[-potentials.ind,]
        }
      }
    }
  return(g)
  }; vazquez.game.el<-cmpfun(vazquez.game.el)

vazquez.game.noise <- function(N, u, p){
    # _N_ is the number of nodes
    # _u_ is the triadic closure parameter. Needs to be between 0 and 1.
    # _p_ is the proportion of edges to rewire
    g <- graph.edgelist(vazquez.game.el(N, u), directed = FALSE)
    g <- rewire(g, mode = "simple", niter = round(p * ecount(g)))
    }; vazquez.game.noise <- cmpfun(vazquez.game.noise)

thresh_diffusion_ig = function(g, infected, thresholds, iters, friendshipeffect = 2){
    # _g_ is a simple, undirected igraph graph object
    # _infected_ is the initial conditions as to who is infected (TRUE) and who is not (FALSE)
    # _thresholds_ are the individual thresholds. If the proportion of the population exceeds an individual's threshold, that individual will participate. It should be between 0 and 1.
    # _iters_ is how many waves to go through the threshold experiment.
    # _friendshipeffect_ is how much each friend should count in an individual's decision. Granovetter says "2." That's the default, yet it's arbitrary and can be altered.
    el <- get.edgelist(g, names = FALSE) - 1
    return(thresh_diffusion(iters, infected, thresholds, el, friendshipeffect))
    }; thresh_diffusion_ig = cmpfun(thresh_diffusion_ig)

# A pretty nice example of the simulation
g <- vazquez.game.noise(100, .35, .20)
g.lo <- layout.fruchterman.reingold(g, params=list(niter = 2500, area = 100^3))
infvec <- sample(c(TRUE, FALSE), 100, c(.05, .95), replace = TRUE)
thresh <- runif(100, 0, 1)
diffuse.out <- thresh_diffusion_ig(g, infvec, thresh, 5, friendshipeffect = 2)
par(mfrow = c(2, 3), mar = c(0, 0, 0, 0))
for (i in 1:6){
    plot(g, vertex.size = 5 + 3 * log(1/sqrt(thresh)), vertex.color = c("light blue", "pink")[diffuse.out[, i] + 1], vertex.label = "", layout = g.lo, margin = c(0, 0, 0, 0))
    }
plot(x = c(0:5), y = 100 * apply(diffuse.out, 2, mean), xlab = "Time", ylab = "Engaged (%)", main = "Threshold Model of Collective Action", type = "l")

dev.off()
class_sim <- function(N, u, p.rewire, friendship = 2, trials = 5, p.initial = .05, threshdist = NULL){
    if (is.null(threshdist)){
        thresholds <- runif(N, min = 0, max = 1)
    }
    infvec <- sample(c(TRUE, FALSE), N, c(p.initial, 1 - p.initial), replace = TRUE)
    g <- vazquez.game.noise(N, u, p.rewire)
    infections <- thresh_diffusion_ig(g, infvec, thresholds, trials, friendshipeffect = friendship)
    infectedcases <- round(100 * apply(infections, 2, mean), digits = 3)
    infectionslope <- mean(diff(infectedcases))
    return(list(graph = g, initial = infvec, thresholds = thresholds, infections = infections, cumulative.cases = infectedcases, slope = infectionslope, parameters = list(N = N, u = u, p.rewire = p.rewire, friendship = friendship, trials = trials, p.initial = p.initial, threshdist = threshdist)))
    }; class_sim <- cmpfun(class_sim)

# Let's test increasing the density/transitivity of the network
f1 <- function(){
    temp = class_sim(100, runif(1, min = 0, max = 1), .20, 5)
    return(c(u = temp$parameters$u, slope = temp$slope))
    }; f1 <- cmpfun(f1)
utest <- t(replicate(100, f1()))
summary(lm(utest[,2] ~ utest[,1]))
# Yes, positive effect between u and slope. As density/transitivity increases, so too does the rate of diffusion.

# Let's test increasing the network's randomization
f2 <- function(){
    temp = class_sim(100, .35, runif(1, 0, 1), 5)
    return(c(rw = temp$parameters$p.rewire, slope = temp$slope))
    }; f2 <- cmpfun(f2)
rwtest <- t(replicate(100, f2()))
summary(lm(rwtest[,2] ~ rwtest[,1]))
# Seemingly no effect. Interesting--it would suggest the findings from early reflect density, not transitivity

# Let's test adjusting the friendship constant
f3 <- function(){
    temp = class_sim(100, .35, .20, 5, friendship = runif(1, 1, 5))
    return(c(friendship = temp$parameters$friendship, slope = temp$slope))
    }; f3 <- cmpfun(f3)
friendshiptest <- t(replicate(100, f3()))
summary(lm(friendshiptest[,2] ~ friendshiptest[,1]))
# Nothing. So long as there's a friendship constant, it doesn't seem to matter.

# Let's test both the density and the rewiring effect.
f4 <- function(){
    temp = class_sim(100, runif(1, min = 0, max = 1), runif(1, 0, 1), 5)
    return(c(u = temp$parameters$u, rw = temp$parameters$p.rewire, slope = temp$slope))
    }; f4 <- cmpfun(f4)
urwtest <- t(replicate(100, f4()))
summary(lm(urwtest[,3] ~ urwtest[,1] + urwtest[,2]))
# Again, density effects through u, but nothing from the rewiring test.
# I have increased the number of replications and the results are the same.