.Rhistory

# Run Simulations
for (sim in 1:reps) {
####################
# Seed structures and intial matrices
####################
# Set initial probability matrix (P_g)
P_g <- matrix(data = rep(0, n * m), ncol = m)
# Seed task (external) stimuli
stimMat <- seed_stimuls(intitial_stim = InitialStim,
Tsteps = Tsteps)
# Seed internal thresholds
threshMat <- seed_thresholds(n = n,
m = m,
threshold_means = ThreshM,
threshold_sds = ThreshSD)
# Start task performance
X_g <- matrix(data = rep(0, length(P_g)), ncol = ncol(P_g))
# Create cumulative task performance matrix
X_tot <- X_g
# Create cumulative adjacency matrix
g_tot <-  matrix(data = rep(0, n * n), ncol = n)
colnames(g_tot) <- paste0("v-", 1:n)
rownames(g_tot) <- paste0("v-", 1:n)
####################
# Simulate individual run
####################
# Run simulation
for (t in 1:Tsteps) {
# Current timestep is actually t+1 in this formulation, because first row is timestep 0
# Update stimuli
stimMat <- update_stim(stim_matrix = stimMat,
deltas = deltas,
alpha = alpha,
state_matrix = X_g,
time_step = t)
# Calculate task demand based on global stimuli
P_g <- calc_determ_thresh(time_step        = t + 1, # first row is generation 0
threshold_matrix = threshMat,
stimulus_matrix  = stimMat)
# Update task performance
X_g <- update_task_performance(task_probs   = P_g,
state_matrix = X_g,
quit_prob    = quitP)
# Update social network (previously this was before probability/task update)
g_adj <- temporalNetwork(X_sub_g = X_g,
prob_interact = p,
bias = beta)
g_tot <- g_tot + g_adj
# Adjust thresholds
threshMat <- adjust_thresholds_social_capped(social_network = g_adj,
threshold_matrix = threshMat,
state_matrix = X_g,
epsilon = epsilon,
threshold_max = 100)
# Update total task performance profile
X_tot <- X_tot + X_g
}
####################
# Post run calculations
####################
# Calculate Entropy
entropy <- mutualEntropy(TotalStateMat = X_tot)
# Calculate total task distribution
totalTaskDist <- X_tot / Tsteps
# Create tasktally table
stimMat <- cbind(stimMat, 0:(nrow(stimMat) - 1))
colnames(stimMat)[ncol(stimMat)] <- "t"
# Add total task distributions, entropy values, and graphs to lists
ens_taskDist[[sim]]    <- totalTaskDist
ens_entropy[[sim]]     <- entropy
ens_stim[[sim]]        <- stimMat
ens_thresh[[sim]]      <- threshMat
ens_graphs[[sim]]      <- g_tot / Tsteps
}
# Add to list of lists
groups_taskDist[[i]]    <- ens_taskDist
groups_stim[[i]]        <- ens_stim
groups_thresh[[i]]      <- ens_thresh
groups_entropy[[i]]     <- ens_entropy
groups_graphs[[i]]      <- ens_graphs
}
# Set threshold max/min
thresh_limit <- 100
# Set group size and replicate
size <- 35
size <- size/5
replicate <- 1
example_graph <- groups_graphs[[1]][[1]]
example_thresh <- as.data.frame(groups_thresh[[1]][[1]])
example_thresh$n <- 35
example_thresh$sim <- 0
example_thresh$chunk <- 0
example_thresh$Id <- row.names(example_thresh)
example_thresh$ThreshBias <- example_thresh$Thresh1 - example_thresh$Thresh2
example_thresh$ThreshBiasBounded <- example_thresh$ThreshBias
example_thresh$ThreshBiasBounded[example_thresh$ThreshBiasBounded < -thresh_limit] <- -thresh_limit
example_thresh$ThreshBiasBounded[example_thresh$ThreshBiasBounded > thresh_limit] <- thresh_limit
# If no node reaches upper or lower limits, add for coloring purposes in gephi
if (sum(example_thresh$ThreshBias == thresh_limit) == 0) {
max_row <- data.frame(Thresh1 = 50, Thresh2 = 50,
n = size * 5, sim = 0, chunk = 0,
Id = "Max", ThreshBias = thresh_limit, ThreshBiasBounded = thresh_limit)
example_thresh <- rbind(example_thresh, max_row)
}
if (sum(example_thresh$ThreshBias == -thresh_limit) == 0) {
min_row <- data.frame(Thresh1 = 50, Thresh2 = 50,
n = size * 5, sim = 0, chunk = 0,
Id = "Min", ThreshBias = -thresh_limit, ThreshBiasBounded = -thresh_limit)
example_thresh <- rbind(example_thresh, min_row)
}
not_chosen <- 1 - (( 1 / (nrow(example_graph) - 1)) * p)
expected_random <-  1 - not_chosen^2
example_graph[example_graph <= expected_random] <- 0
# Or take in those in top X percentile
# percentiles <- quantile(example_graph, na.rm = TRUE)
# fiftypercent <- percentiles[3]
# seventyfivepercent <- percentiles[4]
# example_graph[example_graph <= fiftypercent] <- 0
# diag(example_graph) <- 0
# Turn into graph object to get edgelist
g <- graph_from_adjacency_matrix(example_graph, mode = "undirected", weighted = TRUE)
edgelist <- get.edgelist(g)
edgelist <- as.data.frame(edgelist)
names(edgelist) <- c("Source", "Target")
edgelist$Weight <- E(g)$weight
# Write
write.csv(edgelist, file = paste0("output/Networks/ExampleNetworks/GroupSize_LowBiasLowSocInf", 5*size, "_AboveRandom_edgelist.csv"), row.names = FALSE)
write.csv(example_thresh, file = paste0("output/Networks/ExampleNetworks/GroupSize_LowIntBiasLowSocInf", 5*size, "nodelist.csv"), row.names = FALSE)
rm(list = ls())
####################
# Source necessary scripts/libraries
####################
source("scripts/util/__Util__MASTER.R")
####################
# Set global variables
####################
# Initial paramters: Free to change
# Base parameters
Ns             <- c(35) #vector of number of individuals to simulate
m              <- 2 #number of tasks
Tsteps         <- 50000 #number of time steps to run simulation
reps           <- 1 #number of replications per simulation (for ensemble)
# Threshold Parameters
ThreshM        <- rep(50, m) #population threshold means
ThreshSD       <- ThreshM * 0 #population threshold standard deviations
InitialStim    <- rep(0, m) #intital vector of stimuli
deltas         <- rep(0.8, m) #vector of stimuli increase rates
alpha          <- m #efficiency of task performance
quitP          <- 0.2 #probability of quitting task once active
# Social Network Parameters
p              <- 1 #baseline probablity of initiating an interaction per time step
epsilon        <- 0.1 #relative weighting of social interactions for adjusting thresholds
beta           <- 1.07 #probability of interacting with individual in same state relative to others
####################
# Run ensemble simulation
####################
# Prep meta-lists for collection of group size simulations
groups_taskDist    <- list()
groups_stim        <- list()
groups_thresh      <- list()
groups_entropy     <- list()
groups_graphs      <- list()
# Loop through group sizes
for (i in 1:length(Ns)) {
# Set group size
n <- Ns[i]
# Prep lists for collection of simulation outputs from this group size
ens_taskDist    <- list()
ens_entropy     <- list()
ens_stim        <- list()
ens_thresh      <- list()
ens_graphs      <- list()
# Run Simulations
for (sim in 1:reps) {
####################
# Seed structures and intial matrices
####################
# Set initial probability matrix (P_g)
P_g <- matrix(data = rep(0, n * m), ncol = m)
# Seed task (external) stimuli
stimMat <- seed_stimuls(intitial_stim = InitialStim,
Tsteps = Tsteps)
# Seed internal thresholds
threshMat <- seed_thresholds(n = n,
m = m,
threshold_means = ThreshM,
threshold_sds = ThreshSD)
# Start task performance
X_g <- matrix(data = rep(0, length(P_g)), ncol = ncol(P_g))
# Create cumulative task performance matrix
X_tot <- X_g
# Create cumulative adjacency matrix
g_tot <-  matrix(data = rep(0, n * n), ncol = n)
colnames(g_tot) <- paste0("v-", 1:n)
rownames(g_tot) <- paste0("v-", 1:n)
####################
# Simulate individual run
####################
# Run simulation
for (t in 1:Tsteps) {
# Current timestep is actually t+1 in this formulation, because first row is timestep 0
# Update stimuli
stimMat <- update_stim(stim_matrix = stimMat,
deltas = deltas,
alpha = alpha,
state_matrix = X_g,
time_step = t)
# Calculate task demand based on global stimuli
P_g <- calc_determ_thresh(time_step        = t + 1, # first row is generation 0
threshold_matrix = threshMat,
stimulus_matrix  = stimMat)
# Update task performance
X_g <- update_task_performance(task_probs   = P_g,
state_matrix = X_g,
quit_prob    = quitP)
# Update social network (previously this was before probability/task update)
g_adj <- temporalNetwork(X_sub_g = X_g,
prob_interact = p,
bias = beta)
g_tot <- g_tot + g_adj
# Adjust thresholds
threshMat <- adjust_thresholds_social_capped(social_network = g_adj,
threshold_matrix = threshMat,
state_matrix = X_g,
epsilon = epsilon,
threshold_max = 100)
# Update total task performance profile
X_tot <- X_tot + X_g
}
####################
# Post run calculations
####################
# Calculate Entropy
entropy <- mutualEntropy(TotalStateMat = X_tot)
# Calculate total task distribution
totalTaskDist <- X_tot / Tsteps
# Create tasktally table
stimMat <- cbind(stimMat, 0:(nrow(stimMat) - 1))
colnames(stimMat)[ncol(stimMat)] <- "t"
# Add total task distributions, entropy values, and graphs to lists
ens_taskDist[[sim]]    <- totalTaskDist
ens_entropy[[sim]]     <- entropy
ens_stim[[sim]]        <- stimMat
ens_thresh[[sim]]      <- threshMat
ens_graphs[[sim]]      <- g_tot / Tsteps
}
# Add to list of lists
groups_taskDist[[i]]    <- ens_taskDist
groups_stim[[i]]        <- ens_stim
groups_thresh[[i]]      <- ens_thresh
groups_entropy[[i]]     <- ens_entropy
groups_graphs[[i]]      <- ens_graphs
}
example_graph <- groups_graphs[[1]][[1]]
example_thresh <- as.data.frame(groups_thresh[[1]][[1]])
View(example_thresh)
example_thresh$n <- 35
example_thresh$sim <- 0
example_thresh$chunk <- 0
thresh_limit <- 100
# Set group size and replicate
size <- 35
size <- size/5
replicate <- 1
example_thresh$Id <- row.names(example_thresh)
example_thresh$ThreshBias <- example_thresh$Thresh1 - example_thresh$Thresh2
example_thresh$ThreshBiasBounded <- example_thresh$ThreshBias
example_thresh$ThreshBiasBounded[example_thresh$ThreshBiasBounded < -thresh_limit] <- -thresh_limit
example_thresh$ThreshBiasBounded[example_thresh$ThreshBiasBounded > thresh_limit] <- thresh_limit
# If no node reaches upper or lower limits, add for coloring purposes in gephi
if (sum(example_thresh$ThreshBias == thresh_limit) == 0) {
max_row <- data.frame(Thresh1 = 50, Thresh2 = 50,
n = size * 5, sim = 0, chunk = 0,
Id = "Max", ThreshBias = thresh_limit, ThreshBiasBounded = thresh_limit)
example_thresh <- rbind(example_thresh, max_row)
}
if (sum(example_thresh$ThreshBias == -thresh_limit) == 0) {
min_row <- data.frame(Thresh1 = 50, Thresh2 = 50,
n = size * 5, sim = 0, chunk = 0,
Id = "Min", ThreshBias = -thresh_limit, ThreshBiasBounded = -thresh_limit)
example_thresh <- rbind(example_thresh, min_row)
}
# example_thresh$ThreshRatio <- log(example_thresh$Thresh1 / example_thresh$Thresh2)
# example_thresh$ThreshRatioBounded <- example_thresh$ThreshRatio
# example_thresh$ThreshRatioBounded[example_thresh$ThreshRatioBounded < -thresh_limit] <- -thresh_limit
# example_thresh$ThreshRatioBounded[example_thresh$ThreshRatioBounded > thresh_limit] <- thresh_limit
# # If no node reaches upper or lower limits, add for coloring purposes in gephi
# if (sum(example_thresh$ThreshRatioBounded == thresh_limit) == 0) {
#   max_row <- data.frame(Thresh1 = NA, Thresh2 = NA,
#                         n = NA, sim = NA, chunk = NA,
#                         Id = "Max", ThreshRatio = thresh_limit, ThreshRatioBounded = thresh_limit)
#   example_thresh <- rbind(example_thresh, max_row)
# }
# if (sum(example_thresh$ThreshRatioBounded == -thresh_limit) == 0) {
#   min_row <- data.frame(Thresh1 = NA, Thresh2 = NA,
#                         n = NA, sim = NA, chunk = NA,
#                         Id = "Min", ThreshRatio = -thresh_limit, ThreshRatioBounded = -thresh_limit)
#   example_thresh <- rbind(example_thresh, min_row)
# }
# # Calculate values expected
# Zero out interactions equal to or less than random
not_chosen <- 1 - (( 1 / (nrow(example_graph) - 1)) * p)
expected_random <-  1 - not_chosen^2
example_graph[example_graph <= expected_random] <- 0
# Or take in those in top X percentile
# percentiles <- quantile(example_graph, na.rm = TRUE)
# fiftypercent <- percentiles[3]
# seventyfivepercent <- percentiles[4]
# example_graph[example_graph <= fiftypercent] <- 0
# diag(example_graph) <- 0
# Turn into graph object to get edgelist
g <- graph_from_adjacency_matrix(example_graph, mode = "undirected", weighted = TRUE)
edgelist <- get.edgelist(g)
edgelist <- as.data.frame(edgelist)
names(edgelist) <- c("Source", "Target")
edgelist$Weight <- E(g)$weight
# Write
write.csv(edgelist, file = paste0("output/Networks/ExampleNetworks/GroupSize_LowBiasLowSocInf", 5*size, "_AboveRandom_edgelist.csv"), row.names = FALSE)
write.csv(example_thresh, file = paste0("output/Networks/ExampleNetworks/GroupSize_LowIntBiasLowSocInf", 5*size, "nodelist.csv"), row.names = FALSE)
View(example_thresh)
####################
# Source necessary scripts/libraries
####################
source("scripts/util/__Util__MASTER.R")
####################
# Set global variables
####################
# Initial paramters: Free to change
# Base parameters
Ns             <- c(35) #vector of number of individuals to simulate
m              <- 2 #number of tasks
Tsteps         <- 50000 #number of time steps to run simulation
reps           <- 1 #number of replications per simulation (for ensemble)
# Threshold Parameters
ThreshM        <- rep(50, m) #population threshold means
ThreshSD       <- ThreshM * 0 #population threshold standard deviations
InitialStim    <- rep(0, m) #intital vector of stimuli
deltas         <- rep(0.8, m) #vector of stimuli increase rates
alpha          <- m #efficiency of task performance
quitP          <- 0.2 #probability of quitting task once active
# Social Network Parameters
p              <- 1 #baseline probablity of initiating an interaction per time step
epsilon        <- 0.07 #relative weighting of social interactions for adjusting thresholds
beta           <- 1.07 #probability of interacting with individual in same state relative to others
####################
# Run ensemble simulation
####################
# Prep meta-lists for collection of group size simulations
groups_taskDist    <- list()
groups_stim        <- list()
groups_thresh      <- list()
groups_entropy     <- list()
groups_graphs      <- list()
# Loop through group sizes
for (i in 1:length(Ns)) {
# Set group size
n <- Ns[i]
# Prep lists for collection of simulation outputs from this group size
ens_taskDist    <- list()
ens_entropy     <- list()
ens_stim        <- list()
ens_thresh      <- list()
ens_graphs      <- list()
# Run Simulations
for (sim in 1:reps) {
####################
# Seed structures and intial matrices
####################
# Set initial probability matrix (P_g)
P_g <- matrix(data = rep(0, n * m), ncol = m)
# Seed task (external) stimuli
stimMat <- seed_stimuls(intitial_stim = InitialStim,
Tsteps = Tsteps)
# Seed internal thresholds
threshMat <- seed_thresholds(n = n,
m = m,
threshold_means = ThreshM,
threshold_sds = ThreshSD)
# Start task performance
X_g <- matrix(data = rep(0, length(P_g)), ncol = ncol(P_g))
# Create cumulative task performance matrix
X_tot <- X_g
# Create cumulative adjacency matrix
g_tot <-  matrix(data = rep(0, n * n), ncol = n)
colnames(g_tot) <- paste0("v-", 1:n)
rownames(g_tot) <- paste0("v-", 1:n)
####################
# Simulate individual run
####################
# Run simulation
for (t in 1:Tsteps) {
# Current timestep is actually t+1 in this formulation, because first row is timestep 0
# Update stimuli
stimMat <- update_stim(stim_matrix = stimMat,
deltas = deltas,
alpha = alpha,
state_matrix = X_g,
time_step = t)
# Calculate task demand based on global stimuli
P_g <- calc_determ_thresh(time_step        = t + 1, # first row is generation 0
threshold_matrix = threshMat,
stimulus_matrix  = stimMat)
# Update task performance
X_g <- update_task_performance(task_probs   = P_g,
state_matrix = X_g,
quit_prob    = quitP)
# Update social network (previously this was before probability/task update)
g_adj <- temporalNetwork(X_sub_g = X_g,
prob_interact = p,
bias = beta)
g_tot <- g_tot + g_adj
# Adjust thresholds
threshMat <- adjust_thresholds_social_capped(social_network = g_adj,
threshold_matrix = threshMat,
state_matrix = X_g,
epsilon = epsilon,
threshold_max = 100)
# Update total task performance profile
X_tot <- X_tot + X_g
}
####################
# Post run calculations
####################
# Calculate Entropy
entropy <- mutualEntropy(TotalStateMat = X_tot)
# Calculate total task distribution
totalTaskDist <- X_tot / Tsteps
# Create tasktally table
stimMat <- cbind(stimMat, 0:(nrow(stimMat) - 1))
colnames(stimMat)[ncol(stimMat)] <- "t"
# Add total task distributions, entropy values, and graphs to lists
ens_taskDist[[sim]]    <- totalTaskDist
ens_entropy[[sim]]     <- entropy
ens_stim[[sim]]        <- stimMat
ens_thresh[[sim]]      <- threshMat
ens_graphs[[sim]]      <- g_tot / Tsteps
}
# Add to list of lists
groups_taskDist[[i]]    <- ens_taskDist
groups_stim[[i]]        <- ens_stim
groups_thresh[[i]]      <- ens_thresh
groups_entropy[[i]]     <- ens_entropy
groups_graphs[[i]]      <- ens_graphs
}
example_graph <- groups_graphs[[1]][[1]]
example_thresh <- as.data.frame(groups_thresh[[1]][[1]])
example_thresh$n <- 35
example_thresh$sim <- 0
example_thresh$chunk <- 0
thresh_limit <- 100
# Set group size and replicate
size <- 35
size <- size/5
replicate <- 1
# Get graph
example_graph <- soc_networks[[size]][[replicate]]
example_thresh <- as.data.frame(thresh_data[[size]][[replicate]])
example_thresh$Id <- row.names(example_thresh)
example_thresh$ThreshBias <- example_thresh$Thresh1 - example_thresh$Thresh2
example_thresh$ThreshBiasBounded <- example_thresh$ThreshBias
example_thresh$ThreshBiasBounded[example_thresh$ThreshBiasBounded < -thresh_limit] <- -thresh_limit
example_thresh$ThreshBiasBounded[example_thresh$ThreshBiasBounded > thresh_limit] <- thresh_limit
# If no node reaches upper or lower limits, add for coloring purposes in gephi
if (sum(example_thresh$ThreshBias == thresh_limit) == 0) {
max_row <- data.frame(Thresh1 = 50, Thresh2 = 50,
n = size * 5, sim = 0, chunk = 0,
Id = "Max", ThreshBias = thresh_limit, ThreshBiasBounded = thresh_limit)
example_thresh <- rbind(example_thresh, max_row)
}
if (sum(example_thresh$ThreshBias == -thresh_limit) == 0) {
min_row <- data.frame(Thresh1 = 50, Thresh2 = 50,
n = size * 5, sim = 0, chunk = 0,
Id = "Min", ThreshBias = -thresh_limit, ThreshBiasBounded = -thresh_limit)
example_thresh <- rbind(example_thresh, min_row)
}
# example_thresh$ThreshRatio <- log(example_thresh$Thresh1 / example_thresh$Thresh2)
# example_thresh$ThreshRatioBounded <- example_thresh$ThreshRatio
# example_thresh$ThreshRatioBounded[example_thresh$ThreshRatioBounded < -thresh_limit] <- -thresh_limit
# example_thresh$ThreshRatioBounded[example_thresh$ThreshRatioBounded > thresh_limit] <- thresh_limit
# # If no node reaches upper or lower limits, add for coloring purposes in gephi
# if (sum(example_thresh$ThreshRatioBounded == thresh_limit) == 0) {
#   max_row <- data.frame(Thresh1 = NA, Thresh2 = NA,
#                         n = NA, sim = NA, chunk = NA,
#                         Id = "Max", ThreshRatio = thresh_limit, ThreshRatioBounded = thresh_limit)
#   example_thresh <- rbind(example_thresh, max_row)
# }
# if (sum(example_thresh$ThreshRatioBounded == -thresh_limit) == 0) {
#   min_row <- data.frame(Thresh1 = NA, Thresh2 = NA,
#                         n = NA, sim = NA, chunk = NA,
#                         Id = "Min", ThreshRatio = -thresh_limit, ThreshRatioBounded = -thresh_limit)
#   example_thresh <- rbind(example_thresh, min_row)
# }
# # Calculate values expected
# Zero out interactions equal to or less than random
not_chosen <- 1 - (( 1 / (nrow(example_graph) - 1)) * p)
expected_random <-  1 - not_chosen^2
example_graph[example_graph <= expected_random] <- 0
# Or take in those in top X percentile
# percentiles <- quantile(example_graph, na.rm = TRUE)
# fiftypercent <- percentiles[3]
# seventyfivepercent <- percentiles[4]
# example_graph[example_graph <= fiftypercent] <- 0
# diag(example_graph) <- 0
# Turn into graph object to get edgelist
g <- graph_from_adjacency_matrix(example_graph, mode = "undirected", weighted = TRUE)
edgelist <- get.edgelist(g)
edgelist <- as.data.frame(edgelist)
names(edgelist) <- c("Source", "Target")
edgelist$Weight <- E(g)$weight
# Write
write.csv(edgelist, file = paste0("output/Networks/ExampleNetworks/GroupSize_LowBiasLowSocInf", 5*size, "_AboveRandom_edgelist.csv"), row.names = FALSE)
write.csv(example_thresh, file = paste0("output/Networks/ExampleNetworks/GroupSize_LowIntBiasLowSocInf", 5*size, "nodelist.csv"), row.names = FALSE)