Nothing
R/functions_Metropolis.R
In ERPM: Exponential Random Partition Models
Defines functions
sample_new_partition_p3_restricted
compute_size_neighborhood_p3_restricted
sample_new_partition_p2_restricted
compute_size_neighborhood_p2_restricted
sample_new_partition_p1_restricted
compute_size_neighborhood_p1_restricted
reachable_p3
sample_new_partition_p3
compute_size_neighborhood_p3
reachable_p2
sample_new_partition_p2
compute_size_neighborhood_p2
reachable_p1
sample_new_partition_p1
compute_size_neighborhood_p1
reachable
sample_new_partition
compute_size_neighborhood
step_recalculate
recalculate_statistics
draw_step_multiple
draw_step_single
draw_Metropolis_multiple
draw_Metropolis_single
Documented in
draw_Metropolis_multiple
draw_Metropolis_single
######################################################################
## Simulation and estimation of Exponential Random Partition Models ##
## Function implementing the Metropolis HAstings algorithm to ##
## through partitions given a certain model specification ##
## Author: Marion Hoffman ##
######################################################################
## --- FUNCTIONS TO DRAW CHAINS ----
#' Draw Metropolis single
#'
#' Function to sample the model with a Markov chain (single partition procedure).
#'
#'
#' @param theta model parameters
#' @param first.partition, starting partition for the Markov chain
#' @param nodes nodeset (data frame)
#' @param effects effects/sufficient statistics (list with a vector "names", and a vector "objects")
#' @param objects objects used for statistics calculation (list with a vector "name", and a vector "object")
#' @param burnin integer for the number of burn-in steps before sampling
#' @param thining integer for the number of thining steps between sampling
#' @param num.steps number of samples
#' @param neighborhood = c(0.7,0.3,0), way of choosing partitions: probability vector (2 actors swap, merge/division, single actor move, single pair move, 2 pairs swap, 2 groups reshuffle)
#' @param numgroups.allowed = NULL, # vector containing the number of groups allowed in the partition (now, it only works with vectors like num_min:num_max)
#' @param numgroups.simulated = NULL, # vector containing the number of groups simulated
#' @param sizes.allowed = NULL, vector of group sizes allowed in sampling (now, it only works for vectors like size_min:size_max)
#' @param sizes.simulated = NULL, vector of group sizes allowed in the Markov chain but not necessraily sampled (now, it only works for vectors like size_min:size_max)
#' @param return.all.partitions = FALSE option to return the sampled partitions on top of their statistics (for GOF)
#' @return A list
#' @importFrom stats runif
#' @export
#' @examples
#' # define an arbitrary set of n = 6 nodes with attributes, and an arbitrary covariate matrix
#' n <- 6
#' nodes <- data.frame(label = c("A","B","C","D","E","F"),
#' gender = c(1,1,2,1,2,2),
#' age = c(20,22,25,30,30,31))
#' friendship <- matrix(c(0, 1, 1, 1, 0, 0,
#' 1, 0, 0, 0, 1, 0,
#' 1, 0, 0, 0, 1, 0,
#' 1, 0, 0, 0, 0, 0,
#' 0, 1, 1, 0, 0, 1,
#' 0, 0, 0, 0, 1, 0), 6, 6, TRUE)
#'
#' # choose the effects to be included (see manual for all effect names)
#' effects <- list(names = c("num_groups","same","diff","tie"),
#' objects = c("partition","gender","age","friendship"))
#' objects <- list()
#' objects[[1]] <- list(name = "friendship", object = friendship)
#'
#' # set parameter values for each of these effects
#' parameters <- c(-0.2, 0.2, -0.1, 0.5)
#'
#' \donttest{
#' # generate simulated sample, by setting the desired additional parameters for the
#' # Metropolis sampler and choosing a starting point for the chain (first.partition)
#' nsteps <- 100
#' sample <- draw_Metropolis_single(theta = parameters,
#' first.partition = c(1,1,2,2,3,3),
#' nodes = nodes,
#' effects = effects,
#' objects = objects,
#' burnin = 100,
#' thining = 10,
#' num.steps = nsteps,
#' neighborhood = c(0,1,0),
#' numgroups.allowed = 1:n,
#' numgroups.simulated = 1:n,
#' sizes.allowed = 1:n,
#' sizes.simulated = 1:n,
#' return.all.partitions = TRUE)
#'
#'
#' # or: simulate an estimated model
#' partition <- c(1,1,2,2,2,3) # the partition already defined for the (previous) estimation
#' nsimulations <- 1000
#' simulations <- draw_Metropolis_single(theta = estimation$results$est,
#' first.partition = partition,
#' nodes = nodes,
#' effects = effects,
#' objects = objects,
#' burnin = 100,
#' thining = 20,
#' num.steps = nsimulations,
#' neighborhood = c(0,1,0),
#' sizes.allowed = 1:n,
#' sizes.simulated = 1:n,
#' return.all.partitions = TRUE)
#' }
#'
draw_Metropolis_single <- function(theta,
first.partition,
nodes,
effects,
objects,
burnin,
thining,
num.steps,
neighborhood = c(0.7,0.3,0),
numgroups.allowed = NULL,
numgroups.simulated = NULL,
sizes.allowed = NULL,
sizes.simulated = NULL,
return.all.partitions = FALSE) {
num.nodes <- nrow(nodes)
num.effects <- length(effects$names)
# turn the neighborhood weights into probabilities if needed
if(sum(neighborhood) != 1){
neighborhood <- neighborhood / sum(neighborhood)
}
# check whether there are constraints
constraints <- FALSE
if(!is.null(numgroups.allowed) || !is.null(sizes.allowed)) {
constraints <- TRUE
if(is.null(numgroups.allowed)) numgroups.allowed <- 1:nrow(nodes)
if(is.null(numgroups.simulated)) numgroups.simulated <- numgroups.allowed
if(is.null(sizes.allowed)) sizes.allowed <- 1:nrow(nodes)
if(is.null(sizes.simulated)) sizes.simulated <- sizes.allowed
}
# instantiate with the starting network
current.partition <- first.partition
current.z <- computeStatistics(current.partition, nodes, effects, objects)
current.logit <- theta * current.z
# store the statistics collected for all networks simulated after the burn in
all.z <-c()
# store the partitions if needed (for GOF)
if(return.all.partitions){
all.partitions <- c()
}
end.walk <- FALSE
cpt_burnin <- 0
cpt_thining <- 0
cpt_steps <- 0
while(!end.walk){
new.step <- draw_step_single(theta,
current.partition,
current.logit,
current.z,
nodes,
effects,
objects,
neighborhood,
numgroups.simulated,
sizes.simulated)
proba.change <- min(1,new.step$hastings.ratio)
change.made <- (runif(1,0,1) <= proba.change)
# update partition if needed
old.partition <- current.partition
if(change.made) {
current.partition <- new.step$new.partition
current.z <- new.step$new.z
current.logit <- new.step$new.logit
}
cpt_burnin <- cpt_burnin + 1
if(cpt_burnin > burnin) cpt_thining <- cpt_thining + 1
# store the results if we are out of burnin
if(cpt_burnin >= burnin && cpt_thining == thining) {
store <- TRUE
if(constraints) store <- check_sizes(current.partition,sizes.allowed,numgroups.allowed)
if(store) {
all.z <- rbind(all.z,current.z)
cpt_thining <- 0
cpt_steps <- cpt_steps + 1
if(return.all.partitions){ # store all partitions if needed
all.partitions <- rbind(all.partitions,current.partition)
}
} else {
cpt_thining <- thining - 1
}
}
# stop the walk if number of steps reached
end.walk <- (cpt_steps >= num.steps)
}
# TODO: change this hack
row.names(all.z) <- rep("", dim(all.z)[1])
# compute the average statistics and the final network generated
if(!return.all.partitions) {
return( list("draws" = all.z,
"last.partition" = current.partition,
"all.partitions" = NULL))
} else {
return( list("draws" = all.z,
"last.partition" = current.partition,
"all.partitions" = all.partitions))
}
}
# MULTIPLE PARTITIONS PROCEDURE
#' Draw Metropolis multiple
#'
#' Function to sample the model with a Markov chain (single partition procedure).
#'
#' @param theta model parameters
#' @param first.partitions starting partition for the Markov chain
#' @param presence.tables matrix indicating which actors were present for each observations (mandatory)
#' @param nodes node set (data frame)
#' @param effects effects/sufficient statistics (list with a vector "names", and a vector "objects")
#' @param objects objects used for statistics calculation (list with a vector "name", and a vector "object")
#' @param burnin integer for the number of burn-in steps before sampling
#' @param thining integer for the number of thining steps between sampling
#' @param num.steps number of samples
#' @param neighborhood = c(0.7,0.3,0), way of choosing partitions: probability vector (2 actors swap, merge/division, single actor move, single pair move, 2 pairs swap, 2 groups reshuffle)
#' @param numgroups.allowed = NULL, # vector containing the number of groups allowed in the partition (now, it only works with vectors like num_min:num_max)
#' @param numgroups.simulated = NULL, # vector containing the number of groups simulated
#' @param sizes.allowed = NULL, vector of group sizes allowed in sampling (now, it only works for vectors like size_min:size_max)
#' @param sizes.simulated = NULL, vector of group sizes allowed in the Markov chain but not necessraily sampled (now, it only works for vectors like size_min:size_max)
#' @param return.all.partitions = FALSE, option to return the sampled partitions on top of their statistics (for GOF)
#' @param verbose logical: should intermediate results during the estimation be printed or not? Defaults to FALSE.
#' @return A list
#' @importFrom stats runif
#' @export
#' @examples
#' # define an arbitrary set of n = 6 nodes with attributes, and an arbitrary covariate matrix
#' n <- 6
#' nodes <- data.frame(label = c("A","B","C","D","E","F"),
#' gender = c(1,1,2,1,2,2),
#' age = c(20,22,25,30,30,31))
#' friendship <- matrix(c(0, 1, 1, 1, 0, 0,
#' 1, 0, 0, 0, 1, 0,
#' 1, 0, 0, 0, 1, 0,
#' 1, 0, 0, 0, 0, 0,
#' 0, 1, 1, 0, 0, 1,
#' 0, 0, 0, 0, 1, 0), 6, 6, TRUE)
#'
#' # specify whether nodes are present at different points of time
#' presence.tables <- matrix(c(1, 1, 1, 1, 1, 1,
#' 0, 1, 1, 1, 1, 1,
#' 1, 0, 1, 1, 1, 1), 6, 3)
#'
#' # choose effects to be included in the estimated model
#' effects_multiple <- list(names = c("num_groups","same","diff","tie","inertia_1"),
#' objects = c("partitions","gender","age","friendship","partitions"),
#' objects2 = c("","","","",""))
#' objects_multiple <- list()
#' objects_multiple[[1]] <- list(name = "friendship", object = friendship)
#'
#' # set parameter values for each of these effects
#' parameters <- c(-0.2,0.2,-0.1,0.5,1)
#'
#' # set a starting point for the simulation
#' first.partitions <- matrix(c(1, 1, 2, 2, 2, 3,
#' NA, 1, 1, 2, 2, 2,
#' 1, NA, 2, 3, 3, 1), 6, 3)
#'
#' \donttest{
#' # generate the simulated sample
#' nsteps <- 50
#' sample <- draw_Metropolis_multiple(theta = parameters,
#' first.partitions = first.partitions,
#' nodes = nodes,
#' presence.tables = presence.tables,
#' effects = effects_multiple,
#' objects = objects_multiple,
#' burnin = 100,
#' thining = 100,
#' num.steps = nsteps,
#' neighborhood = c(0,1,0),
#' numgroups.allowed = 1:n,
#' numgroups.simulated = 1:n,
#' sizes.allowed = 1:n,
#' sizes.simulated = 1:n,
#' return.all.partitions = TRUE)
#' }
#'
draw_Metropolis_multiple <- function(theta,
first.partitions,
presence.tables,
nodes,
effects,
objects,
burnin,
thining,
num.steps,
neighborhood = c(0.7,0.3,0),
numgroups.allowed,
numgroups.simulated,
sizes.allowed,
sizes.simulated,
return.all.partitions = FALSE,
verbose = FALSE)
{
num.nodes <- nrow(nodes)
num.effects <- length(effects$names)
num.obs <- ncol(presence.tables)
# turn the neighborhood weights into probabilities if needed
if(sum(neighborhood) != 1){
neighborhood <- neighborhood / sum(neighborhood)
}
# check whether there are constraints
constraints <- FALSE
if(!is.null(numgroups.allowed) || !is.null(sizes.allowed)) {
constraints <- TRUE
if(is.null(numgroups.allowed)) numgroups.allowed <- 1:nrow(nodes)
if(is.null(numgroups.simulated)) numgroups.simulated <- numgroups.allowed
if(is.null(sizes.allowed)) sizes.allowed <- 1:nrow(nodes)
if(is.null(sizes.simulated)) sizes.simulated <- sizes.allowed
}
# instantiate with the starting network
current.partitions <- first.partitions
current.z.contributions <- computeStatistics_multiple(current.partitions, presence.tables, nodes, effects, objects)
current.z <- rowSums(current.z.contributions)
current.logit <- theta * current.z
# store the statistics collected for all networks simulated after the burn in
all.z <-c()
# store the partitions if needed (for GOF)
if(return.all.partitions){
all.partitions <- list()
cpt_all <- 0
}
end.walk <- FALSE
cpt_burnin <- 0
cpt_thining <- 0
cpt_steps <- 0
while(!end.walk){
new.step <- draw_step_multiple(theta,
current.partitions,
current.logit,
current.z.contributions,
current.z,
presence.tables,
nodes,
effects,
objects,
neighborhood,
numgroups.simulated,
sizes.simulated)
proba.change <- min(1,new.step$hastings.ratio)
change.made <- (runif(1,0,1) <= proba.change)
move <- new.step$move
# update partition if needed
old.partitions <- current.partitions
if(change.made) {
current.partitions <- new.step$new.partitions
current.z.contributions <- new.step$new.z.contributions
current.z <- new.step$new.z
current.logit <- new.step$new.logit
}
cpt_burnin <- cpt_burnin + 1
if(cpt_burnin > burnin) cpt_thining <- cpt_thining + 1
if(cpt_burnin %% 10000 == 0 && verbose) cat(cpt_burnin, "\n")
# store the results if we are out of burnin
if(cpt_burnin >= burnin && cpt_thining == thining) {
store <- TRUE
if(constraints) { # check all partitions one by one
for(o in 1:num.obs) store <- store && check_sizes(current.partitions[as.logical(presence.tables[,o]),o],sizes.allowed,numgroups.allowed)
}
if(store) {
all.z <- rbind(all.z,current.z)
cpt_thining <- 0
cpt_steps <- cpt_steps+1
# store all partitions if needed
if(return.all.partitions){
all.partitions[[cpt_steps]] <- current.partitions
}
} else {
cpt_thining <- thining-1
}
}
# stop the walk if number of steps reached
end.walk <- (cpt_steps >= num.steps)
}
# TODO: change this hack
row.names(all.z) <- rep("", dim(all.z)[1])
# compute the average statistics and the final network generated
if(!return.all.partitions) {
return( list("draws" = all.z,
"last.partitions" = current.partitions,
"all.partitions" = NULL))
} else {
return( list("draws" = all.z,
"last.partitions" = current.partitions,
"all.partitions" = all.partitions))
}
}
# MULTIPLE PARTITIONS PROCEDURE
## --- FUNCTIONS TO DRAW ONE STEP IN THE CHAIN ----
# function to draw next partition and calculate HAstings ratio (one step in the Metropolis algorithm)
draw_step_single <- function(theta,
current.partition,
current.logit,
current.z,
nodes,
effects,
objects,
neighborhood,
numgroups.simulated,
sizes.simulated){
# pick a neighborhood
# 1 = swap two actors,
# 2 = merge 2 groups or split a group in two,
# 3 = move one actor,
## 4 = move one pair of nodes (FOR NOW REMOVED)
## 5 = swap two pairs of nodes (FOR NOW REMOVED)
## 6 = reshuffle the members of two groups (FOR NOW REMOVED)
n_neighborhood <- length(neighborhood)
move <- sample(1:n_neighborhood,1,prob=neighborhood)
current.sizes <- rep(0,n_neighborhood)
new.sizes <- rep(0,n_neighborhood)
# calculate current size of neighborhood and pick new partition for the chosen move
current.size <- compute_size_neighborhood(move, current.partition, numgroups.simulated, sizes.simulated)
if(current.size$total > 0) {
new.partition <- sample_new_partition(move, current.partition, current.size, numgroups.simulated, sizes.simulated)
new.size <- compute_size_neighborhood(move, new.partition, numgroups.simulated, sizes.simulated)
} else {
new.partition <- current.partition
new.size <- current.size
}
current.sizes[move] <- current.size$total
new.sizes[move] <- new.size$total
# if change, calculate sizes for all potential moves (depending on available neighborhoods)
if(current.size$total > 0) {
for(i in 1:n_neighborhood){
if(neighborhood[i] > 0 && i != move && reachable(i,current.partition,new.partition)) {
cs <- compute_size_neighborhood(i, current.partition, numgroups.simulated, sizes.simulated)
ns <- compute_size_neighborhood(i, new.partition, numgroups.simulated, sizes.simulated)
current.sizes[i] <- cs$total
new.sizes[i] <- ns$total
}
}
}
# compute new statistics only if it changed
if(current.size$total > 0) {
new.z <- computeStatistics(new.partition, nodes, effects, objects)
} else { new.z <- current.z }
new.logit <- theta * new.z
# calculate acceptance ratio if needed
if(current.size$total > 0) {
indexes <- current.sizes > 0
neighborhoods.ratio <- sum(neighborhood[indexes] / new.sizes[indexes]) /
sum(neighborhood[indexes] / current.sizes[indexes])
} else { neighborhoods.ratio <- 1 }
hastings.ratio <- (exp(sum(new.logit) - sum(current.logit))) * neighborhoods.ratio
return(list("hastings.ratio" = hastings.ratio,
"new.partition" = new.partition,
"new.z" = new.z,
"new.logit" = new.logit))
}
# function to draw next partition and calculate HAstings ratio (one step in the Metropolis algorithm)
draw_step_multiple <- function(theta,
current.partitions,
current.logit,
current.z.contributions,
current.z,
presence.tables,
nodes,
effects,
objects,
neighborhood,
numgroups.simulated,
sizes.simulated){
num.nodes <- nrow(nodes)
num.obs <- ncol(presence.tables)
new.partitions <- current.partitions
rand.o <- sample(1:num.obs,1)
nodes.rand.o <- as.logical(presence.tables[,rand.o])
current.partition <- current.partitions[nodes.rand.o,rand.o]
# calculate sizes for all potential moves (depending on available neighborhoods)
current.sizes <- rep(0,length(neighborhood))
# pick a neighborhood
# 1 = swap two actors,
# 2 = merge 2 groups or split a group in two,
# 3 = move one actor,
# 4 = move one pair of nodes (FOR NOW REMOVED)
# 5 = swap two pairs of nodes (FOR NOW REMOVED)
# 6 = reshuffle the members of two groups (FOR NOW REMOVED)
# 7 = reproduce a pair from the previous partition (empty for first time point) (FOR NOW REMOVED)
n_neighborhood <- length(neighborhood)
move <- sample(1:n_neighborhood,1,prob=neighborhood)
# calculate current size of neighborhood and pick new partition for neighborhoods 1 to 6
for(i in 1:n_neighborhood){
if(neighborhood[i] > 0){
current.size <- compute_size_neighborhood(i, current.partition, numgroups.simulated, sizes.simulated)
current.sizes[i] <- current.size$total
if(move == i) {
if(current.size$total > 0) new.partition <- sample_new_partition(i, current.partition, current.size, numgroups.simulated, sizes.simulated)
if(current.size$total == 0) new.partition <- current.partition
new.partitions[,rand.o] <- rep(NA,num.nodes)
new.partitions[nodes.rand.o,rand.o] <- new.partition
}
}
}
# # this should be for neighborhood 7
# if(neighborhood[7] > 0){
# if(rand.o == 1) current.partition2 <- 1:num.nodes
# else current.partition2 <- current.partitions[,rand.o-1]
# current.size <- compute_size_neighborhood(7, current.partitions[,rand.o], partition2 = current.partition2, numgroups.simulated, sizes.simulated)
# current.sizes[7] <- current.size$total
# if(move == 7){
# if(current.size$total > 0) new.partition <- sample_new_partition(7, current.partitions[,rand.o], current.size, partition2= current.partition2, numgroups.simulated, sizes.simulated)
# if(current.size$total == 0) new.partition <- current.partition
# new.partitions[,rand.o] <- new.partition
# }
# }
# compute new statistics only if it changed
if(all(current.partitions == new.partitions, na.rm = TRUE)) {
new.z.contributions <- current.z.contributions
new.z <- current.z
} else {
recalculated.stats <- recalculate_statistics(new.partitions, rand.o, nodes.rand.o, nodes, effects, objects, current.z.contributions)
new.z.contributions <- recalculated.stats$new.z.contributions
new.z <- recalculated.stats$new.z
}
new.logit <- theta * new.z
# chose whether to change or not
new.sizes <- rep(0,n_neighborhood)
for(i in 1:n_neighborhood) {
if(neighborhood[i] > 0){
new.size <- compute_size_neighborhood(i, new.partition, numgroups.simulated, sizes.simulated)
new.sizes[i] <- new.size$total
}
}
# # this should be for neighborhood 7
# if(neighborhood[7] > 0){
# new.size <- compute_size_neighborhood(7, new.partition, current.partition2, numgroups.simulated, sizes.simulated)
# new.sizes[7] <- new.size$total
# }
neighborhoods.ratio <- sum(current.sizes * neighborhood) /
sum(new.sizes * neighborhood)
hastings.ratio <- (exp(sum(new.logit) - sum(current.logit))) * neighborhoods.ratio
return(list("hastings.ratio" = hastings.ratio,
"new.partitions" = new.partitions,
"new.z.contributions" = new.z.contributions,
"new.z" = new.z,
"new.logit" = new.logit,
"move" = move))
}
# Recalculate statistics for a changed partition in multiple estimation
recalculate_statistics <- function(new.partitions, rand.o, nodes.rand.o, nodes, effects, objects, current.z.contributions) {
# store new statistics
new.z.contributions <- current.z.contributions
# calculate separately each effect
for(e in 1:length(effects$names)) {
new.z.contributions[e,rand.o] <- step_recalculate(new.partitions, rand.o, nodes.rand.o, nodes, effects, objects, e)
}
new.z <- rowSums(new.z.contributions)
return(list(new.z.contributions = new.z.contributions,
new.z = new.z))
}
step_recalculate <- function(new.partitions, rand.o, nodes.rand.o, nodes, effects, objects, e){
num.nodes <- nrow(nodes)
num.obs <- ncol(new.partitions)
object.name <- effects$objects[e]
# TIE EFFECT: keep only present nodes
if(effects$names[e] == "tie"){
effects.temp <- list(names="tie",objects="net.temp")
for(ob in 1:length(objects)){
if(objects[[ob]][[1]] == object.name){
net <- objects[[ob]][[2]]
objects.temp <- list(list(name="net.temp",object=net[nodes.rand.o,nodes.rand.o]))
}
}
nodes.temp <- nodes[nodes.rand.o,]
new.z.contribution <- as.numeric(computeStatistics(new.partitions[nodes.rand.o,rand.o], nodes.temp, effects.temp, objects.temp))
# SAME_VAR EFFECT: adapt the varying covariate
} else if(effects$names[e] == "same_var"){
effects.temp <- list(names="same",objects="var.temp")
for(ob in 1:length(objects)){
if(objects[[ob]][[1]] == object.name){
atts <- objects[[ob]][[2]]
}
}
objects.temp <- list()
nodes.temp <- nodes[nodes.rand.o,]
nodes.temp$var.temp <- atts[nodes.rand.o,rand.o]
new.z.contribution <- as.numeric(computeStatistics(new.partitions[nodes.rand.o,rand.o], nodes.temp, effects.temp, objects.temp))
# TIE_VAR EFFECT: adapt the varying network
} else if(effects$names[e] == "tie_var"){
effects.temp <- list(names="tie",objects="net.temp")
for(ob in 1:length(objects)){
if(objects[[ob]][[1]] == object.name){
nets <- objects[[ob]][[2]]
net <- nets[[rand.o]]
objects.temp <- list(list(name="net.temp",object=net[nodes.rand.o,nodes.rand.o]))
}
}
nodes.temp <- nodes[nodes.rand.o,]
new.z.contribution <- as.numeric(computeStatistics(new.partitions[nodes.rand.o,rand.o], nodes.temp, effects.temp, objects.temp))
# TIE_VAR_X_DIFF EFFECT: adapt the varying network
} else if(effects$names[e] == "tie_var_X_diff"){
effects.temp <- list(names="tie_X_diff",objects="net.temp",objects2=effects$objects2[e])
for(ob in 1:length(objects)){
if(objects[[ob]][[1]] == object.name){
nets <- objects[[ob]][[2]]
net <- nets[[rand.o]]
objects.temp <- list(list(name="net.temp",object=net[nodes.rand.o,nodes.rand.o]))
}
}
nodes.temp <- nodes[nodes.rand.o,]
new.z.contribution <- as.numeric(computeStatistics(new.partitions[nodes.rand.o,rand.o], nodes.temp, effects.temp, objects.temp))
# INERTIA_1 EFFECT: need to recalculate the partition before and after
} else if(effects$names[e] == "inertia_1"){
aff.temp <- as.matrix(table(data.frame(actor = 1:num.nodes, group= new.partitions[,rand.o])))
adj.temp <- aff.temp %*% t(aff.temp)
diag(adj.temp) <- 0
if(rand.o > 1) { # removed && length(groups[[rand.o]]) > 0
aff.temp2 <- as.matrix(table(data.frame(actor = 1:num.nodes, group= new.partitions[,rand.o-1])))
adj.temp2 <- aff.temp2 %*% t(aff.temp2)
diag(adj.temp2) <- 0
adj12 <- adj.temp * adj.temp2
new.z.contribution <- sum(1/2 * adj12)
} else {
new.z.contribution <- 0
}
if(rand.o < num.obs) { # removed && length(groups[[rand.o]]) > 0
aff.temp2 <- as.matrix(table(data.frame(actor = 1:num.nodes, group= new.partitions[,rand.o+1])))
adj.temp2 <- aff.temp2 %*% t(aff.temp2)
diag(adj.temp2) <- 0
adj12 <- adj.temp * adj.temp2
new.z.contribution <- sum(1/2 * adj12)
}
# ANY OTHER EFFECT: the effect also exists in the single partition version
} else {
effects.temp <- list(names=effects$names[e],objects=effects$objects[e])
objects.temp <- objects
nodes.temp <- nodes[nodes.rand.o,]
new.z.contribution <- as.numeric(computeStatistics(new.partitions[nodes.rand.o,rand.o], nodes.temp, effects.temp, objects.temp))
}
return(new.z.contribution)
}
## --- FUNCTIONS FOR THE PROPOSALS MADE AT EACH STEP ----
## GENERIC FUNCTION TO COMPUTE NEIGHBORHOOD SIZE
compute_size_neighborhood <- function(i,
partition,
numgroups.simulated = NULL,
sizes.simulated = NULL){#,partition2 = NULL, # for neighborhood 7, now deprecated){
# check whether there are constraints
constraints <- FALSE
if(!is.null(numgroups.simulated) || !is.null(sizes.simulated)) {
constraints <- TRUE
}
if(i == 1) {
if(!constraints) return(compute_size_neighborhood_p1(partition))
else return(compute_size_neighborhood_p1_restricted(partition, numgroups.simulated, sizes.simulated))
}
if(i == 2) {
if(!constraints) return(compute_size_neighborhood_p2(partition))
else return(compute_size_neighborhood_p2_restricted(partition, numgroups.simulated, sizes.simulated))
}
if(i == 3) {
if(!constraints) return(compute_size_neighborhood_p3(partition))
else return(compute_size_neighborhood_p3_restricted(partition, numgroups.simulated, sizes.simulated))
}
# For now, other neighborhoods removed
# if(i == 4) {
# if(is.null(sizes.simulated)) return(compute_size_neighborhood_p4(partition))
# else return(compute_size_neighborhood_p4_restricted(partition, numgroups.simulated, sizes.simulated))
# }
# if(i == 5) {
# if(is.null(sizes.simulated)) return(compute_size_neighborhood_p5(partition))
# else return(compute_size_neighborhood_p5_restricted(partition, numgroups.simulated, sizes.simulated))
# }
# if(i == 6) {
# if(is.null(sizes.simulated)) return(compute_size_neighborhood_p6(partition))
# else return(compute_size_neighborhood_p6_restricted(partition, numgroups.simulated, sizes.simulated))
# }
# if(i == 7) {
# if(is.null(sizes.simulated)) return(compute_size_neighborhood_p7(partition))
# else return(compute_size_neighborhood_p7_restricted(partition, numgroups.simulated, sizes.simulated))
# }
}
## GENERIC FUNCTION TO COMPUTE NEIGHBORHOOD SIZE
sample_new_partition <- function(i,
current.partition,
size_neighborhood,
numgroups.simulated = NULL,
sizes.simulated = NULL){#,current.partition2 = NULL, # for neighborhood 7, now deprecated){
# check whether there are constraints
constraints <- FALSE
if(!is.null(numgroups.simulated) || !is.null(sizes.simulated)) {
constraints <- TRUE
}
if(i == 1) {
if(!constraints) return(sample_new_partition_p1(current.partition, size_neighborhood))
else return(sample_new_partition_p1_restricted(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated))
}
if(i == 2) {
if(!constraints) return(sample_new_partition_p2(current.partition, size_neighborhood))
else return(sample_new_partition_p2_restricted(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated))
}
if(i == 3) {
if(!constraints) return(sample_new_partition_p3(current.partition, size_neighborhood))
else return(sample_new_partition_p3_restricted(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated))
}
# for now, other neighborhoods removed
# if(i == 4) {
# if(is.null(sizes.simulated)) return(sample_new_partition_p4(current.partition, size_neighborhood))
# else return(sample_new_partition_p4_restricted(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated))
# }
# if(i == 5) {
# if(is.null(sizes.simulated)) return(sample_new_partition_p5(current.partition, size_neighborhood))
# else return(sample_new_partition_p5_restricted(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated))
# }
# if(i == 6) {
# if(is.null(sizes.simulated)) return(sample_new_partition_p6(current.partition, size_neighborhood))
# else return(sample_new_partition_p6_restricted(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated))
# }
# if(i == 7) {
# if(is.null(sizes.simulated)) return(sample_new_partition_p7(current.partition, size_neighborhood, current.partition2))
# else return(sample_new_partition_p7_restricted(current.partition, size_neighborhood, current.partition2, numgroups.simulated, sizes.simulated))
# }
}
## GENERIC FUNCTION TO TEST WHETHER A PARTITION IS REACHABLE WITH A GIVEN NEIGHBORHOOD
reachable <- function(i,partition1,partition2){
if(i == 1) {
return(reachable_p1(partition1,partition2))
} else if(i == 2) {
return(reachable_p2(partition1,partition2))
} else if(i == 3) {
return(reachable_p3(partition1,partition2))
}
# for now, other nieghborhoods, removed
# else if(i == 4) {
# return(reachable_p4(partition1,partition2))
# } else if(i == 5) {
# return(reachable_p5(partition1,partition2))
# } else if(i == 6) {
# return(reachable_p6(partition1,partition2))
# }
}
## --- SPECIFIC PROPOSALS ----
## NEIGHBORHOOD PI 1: only swaps of two nodes (careful, cannot be used alone)
compute_size_neighborhood_p1 <- function(partition){
# find isolates, pairs and groups>2
num.groups <- max(partition)
num.nodes <- length(partition)
sizes <- table(partition)
isolates <- as.vector(which(sizes == 1))
pairs <- as.vector(which(sizes == 2))
others <- as.vector(which(sizes > 2))
nums.swaps <- matrix(0,num.nodes,num.nodes)
for(i in 1:(num.nodes-1)){
for(j in (i+1):num.nodes){
g_i <- partition[i]
g_j <- partition[j]
allowed <- TRUE
if(g_i == g_j) allowed <- FALSE # two members of the same group cannot swap, it wouldn't change anything
if(g_i %in% isolates && g_j %in% isolates){
allowed <- FALSE # two isolates cannot swap, it wouldn't change anything
}
if(allowed) nums.swaps[i,j] <- 1
}
}
return(list(isolates = isolates,
pairs = pairs,
others = others,
nums.swaps = nums.swaps,
num.swaps = sum(nums.swaps),
total = sum(nums.swaps)))
}
# sample partition neighborhood 1
sample_new_partition_p1 <- function(current.partition, size_neighborhood){
# calculate the number of neighbor partitions
num.nodes <- length(current.partition)
num.groups <- max(current.partition)
size <- size_neighborhood
isolates <- size$isolates
pairs <- size$pairs
others <- size$others
nums.swaps <- size$nums.swaps
num.swaps <- size$num.swaps
total <- size$total
if(total == 0) return(current.partition)
pick.1 <- sample(total,1)
# decide which actors to swap
new.partition <- current.partition
swap <- which(nums.swaps == 1)[pick.1]
j <- ((swap-1) %/% num.nodes) + 1
i <- (swap-1) %% num.nodes + 1
new.partition[i] <- current.partition[j]
new.partition[j] <- current.partition[i]
new.partition <- order_groupids(new.partition)
return(new.partition)
}
# reachable neighborhood 1
reachable_p1 <- function(partition1,partition2){
# they are not neighbors if they have a different number of groups
if(max(partition1) != max(partition2)) return(FALSE)
num.nodes <- length(partition1)
check <- FALSE
i <- 1
j <- i+1
# try remove all pairs of nodes
while(!check && i <= num.nodes){
newp1 <- order_groupids( partition1[-c(i,j)] )
newp2 <- order_groupids( partition2[-c(i,j)] )
# check if they swapped
if(all(newp1 == newp2)) {
othersi1 <- which(partition1 == partition1[i])
othersi1 <- othersi1[othersi1 != i]
othersj1 <- which(partition1 == partition1[j])
othersj1 <- othersj1[othersj1 != j]
othersi2 <- which(partition2 == partition2[i])
othersi2 <- othersi2[othersi2 != i]
othersj2 <- which(partition2 == partition2[j])
othersj2 <- othersj2[othersj2 != j]
check <- (setequal(othersi1,othersj2) && setequal(othersi2,othersj1))
if(check) break
}
j <- j+1
if(j > num.nodes){
i <- i+1
j <- i+1
}
}
return(check)
}
## NEIGHBORHOOD PI 2: only merges and divisions of 2 groups
compute_size_neighborhood_p2 <- function(partition){
# calculate the number of neighbor partitions
num.nodes <- length(partition)
num.groups <- max(partition)
num.merges <- num.groups*(num.groups-1)/2
nums.divisions <- rep(0,num.groups)
num.divisions <- 0
for(k in 1:num.groups){
sk <- length(which(partition == k))
nums.divisions[k] <- 2^(sk-1) - 1
num.divisions <- num.divisions + nums.divisions[k]
}
return(list(num.merges = num.merges,
num.divisions = num.divisions,
nums.divisions = nums.divisions,
total = num.merges + num.divisions))
}
# sample partition neighborhood 2
sample_new_partition_p2 <- function(current.partition, size_neighborhood){
# calculate the number of neighbor partitions
num.nodes <- length(current.partition)
num.groups <- max(current.partition)
size <- size_neighborhood
num.merges <- size$num.merges
num.divisions <- size$num.divisions
nums.divisions <- size$nums.divisions
total <- size$total
if(total == 0) return(current.partition)
# decide between merge or division (or self loop)
pick.1 <- sample(total,1)
new.partition <- current.partition
# if merge
if(pick.1 <= num.merges) {
# pick 2 groups to merge
old_gs <- sample(num.groups,2,replace = FALSE)
# reassign one of the groups and remove useless ids
new.partition[which(new.partition == old_gs[2])] <- old_gs[1]
new.partition <- order_groupids(new.partition)
}else if(pick.1 <= (num.merges+num.divisions)){
# pick group to divide
pick.2 <- sample(num.divisions,1)
sums <- unlist(lapply(1:num.groups,FUN = function(x){sum(nums.divisions[1:x])}))
starts <- c(0,sums[1:(num.groups-1)])
ends <- sums[1:num.groups]
old_g <- which(starts < pick.2 & ends >= pick.2)
# reassign one of the groups and remove useless ids
new.groups <- sample(c(old_g,num.groups+1),length(which(new.partition==old_g)),replace = TRUE)
found <- (length(unique(new.groups)) > 1)
while(!found){
new.groups <- sample(c(old_g,num.groups+1),length(which(new.partition==old_g)),replace = TRUE)
found <- (length(unique(new.groups)) > 1)
}
new.partition[which(new.partition == old_g)] <- new.groups
new.partition <- order_groupids(new.partition)
}
return(new.partition)
}
# reachable neighborhood 2
reachable_p2 <- function(partition1,partition2){
# they are not neighbors if one does not have one more group
if(abs(max(partition1) - max(partition2))!=1) return(FALSE)
num.nodes <- length(partition1)
num.groups1 <- max(partition1)
num.groups2 <- max(partition2)
check <- FALSE
g1 <- 1
g2 <- g1+1
# try remove groups (to check divisions) and pairs of groups (to check merges)
while(!check && g1 <= num.groups1){
# first check whether g1 is not divided in partition2
membersg1 <- which(partition1 == g1)
if(length(membersg1) == num.nodes) {
newp1 <- c()
newp2 <- c()
} else {
newp1 <- order_groupids( partition1[-membersg1] )
newp2 <- order_groupids( partition2[-membersg1] )
}
if(all(newp1 == newp2)) {
check <- length(unique(partition2[membersg1])) == 2
if(check) break
}
# second check whether g1 and g2 are merged in partition2
membersg2 <- which(partition1 == g2)
if(length(c(membersg1,membersg2)) == num.nodes) {
newp1 <- c()
newp2 <- c()
} else {
newp1 <- order_groupids( partition1[-c(membersg1,membersg2)] )
newp2 <- order_groupids( partition2[-c(membersg1,membersg2)] )
}
if(all(newp1 == newp2)) {
check <- length(unique(partition2[c(membersg1,membersg2)])) == 1
if(check) break
}
g2 <- g2+1
if(g2 > num.nodes){
g1 <- g1+1
g2 <- g1+1
}
}
return(check)
}
## NEIGHBORHOOD PI 3: only swaps of one node
compute_size_neighborhood_p3 <- function(partition){
# find isolates, pairs and groups>2
num.groups <- max(partition)
num.nodes <- length(partition)
sizes <- table(partition)
isolates <- as.vector(which(sizes == 1))
pairs <- as.vector(which(sizes == 2))
others <- as.vector(which(sizes > 2))
nums.swaps <- rep(0,num.nodes)
done.pairs <- rep(0,length(pairs))
done.isolates <- rep(0,length(isolates))
for(i in 1:num.nodes){
g <- partition[i]
if(g %in% isolates){
done.isolates[isolates == g] <- 1
nums.swaps[i] <- num.groups - sum(done.isolates) # can join any other group, except other isolates before (because already done) and itself
}
if(g %in% pairs){
nums.swaps[i] <- num.groups - done.pairs[pairs == g] # can join any other group or create its isolate, except if the possibilty of splitting the pair in 2 isolates is already counted
done.pairs[pairs == g] <- 1
}
if(g %in% others){
nums.swaps[i] <- num.groups # can join any other group or create its own isolate
}
}
return(list(isolates = isolates,
pairs = pairs,
others = others,
nums.swaps = nums.swaps,
num.swaps = sum(nums.swaps),
total = sum(nums.swaps)))
}
# sample partition neighborhood 3
sample_new_partition_p3 <- function(current.partition, size_neighborhood){
# calculate the number of neighbor partitions
num.nodes <- length(current.partition)
num.groups <- max(current.partition)
size <- size_neighborhood
isolates <- size$isolates
pairs <- size$pairs
others <- size$others
nums.swaps <- size$nums.swaps
num.swaps <- size$num.swaps
total <- size$total
if(total == 0) return(current.partition)
pick.1 <- sample(total,1)
# decide which actor to swap
all.groups <- 1:num.groups
new.partition <- current.partition
done.pairs <- rep(0,length(pairs))
done.isolates <- rep(0,length(isolates))
for(i in 1:num.nodes){
if(i == 1) {start <- 0} else {start <- sum(nums.swaps[1:i-1])}
if(i == num.nodes) {end <- num.swaps} else {end <- sum(nums.swaps[1:i])}
g <- current.partition[i]
if(g %in% isolates) done.isolates[isolates == g] <- 1
if(pick.1 > start && pick.1 <= end) {
g <- current.partition[i]
# an isolate is randomly joining another group except isolates before i
if(g %in% isolates){
previousisolates <- isolates[done.isolates == 1]
tosample <- all.groups[!all.groups %in% previousisolates]
if(length(tosample) == 1) newg <- tosample
if(length(tosample) >= 2) newg <- sample(tosample,1)
new.partition[i] <- newg
new.partition <- order_groupids(new.partition)
}
# a member of a pair is randomly joining another group
# or creating its isolate if it's the first of the pair to appear
if(g %in% pairs){
if(done.pairs[pairs == g] == 0){
newg <- sample(all.groups,1)
if(newg == g){
newg <- num.groups + 1
}
} else {
tosample <- all.groups[all.groups != g]
if(length(tosample) == 1) newg <- tosample
if(length(tosample) >= 2) newg <- sample(tosample,1)
}
new.partition[i] <- newg
new.partition <- order_groupids(new.partition)
}
# a member of a bigger group is randomly joining another group
# or creating its isolate
if(g %in% others){
newg <- sample(all.groups,1)
if(newg == g){
newg <- num.groups + 1
}
new.partition[i] <- newg
new.partition <- order_groupids(new.partition)
}
}
if(g %in% pairs) done.pairs[pairs == g] <- 1
}
return(new.partition)
}
# reachable neighborhood 3
reachable_p3 <- function(partition1,partition2){
# they are not neighbors if one doesn't have the same number of groups or one more
if(abs(max(partition1) - max(partition2))>1) return(FALSE)
num.nodes <- length(partition1)
check <- FALSE
i <- 1
# try remove all nodes
while(!check && i <= num.nodes){
newp1 <- order_groupids( partition1[-i] )
newp2 <- order_groupids( partition2[-i] )
# check if it moved
if(all(newp1 == newp2)) {
othersi1 <- which(partition1[-i] == partition1[i])
othersi2 <- which(partition2[-i] == partition2[i])
check <- !setequal(othersi1,othersi2)
if(check) break
}
i <- i+1
}
return(check)
}
## --- SPECIFIC PROPOSALS RESTRICTED VERSION ----
## NEIGHBORHOOD PI 1 RESTRICTED: only swaps of two nodes (careful, cannot be used alone)
# WARNING!!!!: for now it only works if sizes allowed are an interval ([size_min,size_max])
compute_size_neighborhood_p1_restricted <- function(partition, numgroups.simulated, sizes.simulated){
# check if current partition is allowed
if(!check_sizes(partition, sizes.simulated,numgroups.simulated)) stop("The partition we are in is not allowed.")
return(compute_size_neighborhood_p1(partition))
}
# sample partition neighborhood 1 restricted
sample_new_partition_p1_restricted <- function(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated){
return(sample_new_partition_p1(current.partition, size_neighborhood))
}
## NEIGHBORHOOD PI 2 RESTRICTED: only merges and divisions of 2 groups
compute_size_neighborhood_p2_restricted <- function(partition, numgroups.simulated, sizes.simulated){
# check if current partition is allowed
if(!check_sizes(partition, sizes.simulated,numgroups.simulated)) stop("The partition we are in is not allowed.")
# find minimum and maximum size allowed
smax <- max(sizes.simulated)
smin <- min(sizes.simulated)
# calculate the number of neighbor partitions
num.nodes <- length(partition)
num.groups <- max(partition)
sizes <- table(partition)
num.merges <- 0
merges <- matrix(0,num.groups,num.groups)
nums.divisions <- rep(0,num.groups)
num.divisions <- 0
# merges
if(num.groups>1 && (num.groups-1) %in% numgroups.simulated){ # we count merges if we have more than one group and having one less group is allowed
for(g1 in 1:(num.groups-1)){
for(g2 in (g1+1):num.groups){
if(length(which(partition == g1)) + length(which(partition == g2)) <= smax) { # merges are only possible if the new group size is allowed
num.merges <- num.merges + 1
merges[g1,g2] <- 1
}
}
}
}
# divisions
if((num.groups+1) %in% numgroups.simulated){ # we count divisions if having one more group is allowed
for(k in 1:num.groups){
sg <- sizes[k]
if(smin>1){
extras <- sum(unlist(lapply(1:(smin-1),FUN=function(x){choose(sg,x)})))
if(sg%%2 == 0 && sg/2 < smin) extras <- extras - choose(sg,sg/2)/2
} else {
extras <- 0
}
nums.divisions[k] <- 2^(sg-1) - 1 - extras
num.divisions <- num.divisions + nums.divisions[k]
}
}
return(list(num.merges = num.merges,
merges = merges,
num.divisions = num.divisions,
nums.divisions = nums.divisions,
total = num.merges + num.divisions))
}
# sample partition neighborhood 2 restricted
sample_new_partition_p2_restricted <- function(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated){
# calculate the number of neighbor partitions
num.nodes <- length(current.partition)
num.groups <- max(current.partition)
size <- size_neighborhood
num.merges <- size$num.merges
merges <- size$merges
num.divisions <- size$num.divisions
nums.divisions <- size$nums.divisions
total <- size$total
if(total == 0) return(current.partition)
# sizes allowed
smax <- max(sizes.simulated)
smin <- min(sizes.simulated)
# decide between merge or division (or self loop)
pick.1 <- sample(total,1)
all.groups <- 1:num.groups
new.partition <- current.partition
# if merge
if(pick.1 <= num.merges) {
# pick 2 groups to merge
allindexes <- which(merges == 1)
p <- allindexes[pick.1]
indexes <- arrayInd(p, dim(merges))
old_g1 <- indexes[,1]
old_g2 <- indexes[,2]
# reassign one of the groups and remove useless ids
new.partition[which(new.partition == old_g2)] <- old_g1
new.partition <- order_groupids(new.partition)
}
# if division
else if(pick.1 <= (num.merges+num.divisions)){
# pick group to divide
pick.2 <- sample(num.divisions,1)
old_g1 <- 1
cpt2 <- 1
g <- 1
found <- FALSE
while(g<=num.groups && !found){
if(pick.2 >= cpt2 && pick.2 <= sum(nums.divisions[1:g])){
old_g1 <- g
found <- TRUE
}
cpt2 <- cpt2 + nums.divisions[g]
g <- g+1
}
old_g2 <- 0
new_g1 <- old_g1
new_g2 <- num.groups + 1
# reassign one of the groups and remove useless ids
new.groups <- sample(c(old_g1,num.groups+1),length(which(new.partition==old_g1)),replace = TRUE)
found <- (length(unique(new.groups)) > 1) && min(table(new.groups)) >= smin
while(!found){
new.groups <- sample(c(old_g1,num.groups+1),length(which(new.partition==old_g1)),replace = TRUE)
found <- (length(unique(new.groups)) > 1) && min(table(new.groups)) >= smin
}
new.partition[which(new.partition == old_g1)] <- new.groups
new.partition <- order_groupids(new.partition)
}
return(new.partition)
}
## NEIGHBORHOOD PI 3 RESTRICTED: only swaps one node
# WARNING!!!!: for now it only works if sizes allowed are an interval ([size_min,size_max])
compute_size_neighborhood_p3_restricted <- function(partition, numgroups.simulated, sizes.simulated){
# check if current partition is allowed
if(!check_sizes(partition, sizes.simulated, numgroups.simulated)) stop("The partition we are in is not allowed.")
# find maximum size allowed
smax <- max(sizes.simulated)
smin <- min(sizes.simulated)
sizes <- table(partition)
# find isolates, pairs and groups>2
num.groups <- max(partition)
num.nodes <- length(partition)
isolates <- as.vector(which(sizes == 1))
pairs <- as.vector(which(sizes == 2))
others <- as.vector(which(sizes > 2))
abovemin_groups <- as.vector(which(sizes > smin))
belowmax_groups <- as.vector(which(sizes < smax))
nums.swaps <- rep(0,num.nodes)
done.pairs <- rep(0,length(pairs))
done.isolates <- rep(0,length(isolates))
for(i in 1:num.nodes){
g <- partition[i]
if(g %in% isolates && (num.groups-1) %in% numgroups.simulated){
# i can join any other group (<smax), except other isolates done before (because already done) and this one (that would not change anything)
# since it's an isolate, it means that isolates are allowed in that context so we don't worry about smin (=1)
done.isolates[isolates == g] <- 1
nums.swaps[i] <- length(belowmax_groups) - sum(done.isolates)
} # to move an isolate, a partition with one group less must be allowed
if(g %in% pairs && g %in% abovemin_groups){
# if isolates are allowed, then i can join any other group (smax), and splitting if it's not already counted
# here we don't worry about smin (=1) because 2 > smin
gnotbelowmax <- !(g %in% belowmax_groups)
if((num.groups+1) %in% numgroups.simulated && done.pairs[pairs == g] == 0) { # we check whether it's possible to create an isolate
nums.swaps[i] <- length(belowmax_groups) + gnotbelowmax
} else {
nums.swaps[i] <- length(belowmax_groups) + gnotbelowmax - 1
}
done.pairs[pairs == g] <- 1
} # to break a pair, size 2 should not be the minimum allowed
if(g %in% others && g %in% abovemin_groups){
# can join any other group (<smax) or create its own isolate (if it's allowed)
# we don't care here about checking that other groups are above min size, if they exist it means they are already allowed
gnotbelowmax <- !(g %in% belowmax_groups)
if(smin == 1 && (num.groups+1) %in% numgroups.simulated) { # we check whether it's possible to create an isolate
nums.swaps[i] <- length(belowmax_groups) + gnotbelowmax
} else {
nums.swaps[i] <- length(belowmax_groups) + gnotbelowmax - 1
}
} # to break a group, it should not be of minimal size
}
return(list(isolates = isolates,
pairs = pairs,
others = others,
abovemin_groups = abovemin_groups,
belowmax_groups = belowmax_groups,
nums.swaps = nums.swaps,
num.swaps = sum(nums.swaps),
total = sum(nums.swaps)))
}
# sample partition neighborhood 3 restricted
sample_new_partition_p3_restricted <- function(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated){
# calculate the number of neighbor partitions
num.nodes <- length(current.partition)
num.groups <- max(current.partition)
size <- size_neighborhood
isolates <- size$isolates
pairs <- size$pairs
others <- size$others
abovemin_groups <- size$abovemin_groups
belowmax_groups <- size$belowmax_groups
nums.swaps <- size$nums.swaps
num.swaps <- size$num.swaps
total <- size$total
if(total == 0) return(current.partition)
# allowed sizes
smax <- max(sizes.simulated)
smin <- min(sizes.simulated)
pick.1 <- sample(total,1)
# decide which actor to swap
all.groups <- 1:num.groups
new.partition <- current.partition
done.pairs <- rep(0,length(pairs))
done.isolates <- rep(0,length(isolates))
for(i in 1:num.nodes){
if(i == 1) {start <- 0} else {start <- sum(nums.swaps[1:i-1])}
if(i == num.nodes) {end <- num.swaps} else {end <- sum(nums.swaps[1:i])}
g <- current.partition[i]
if(g %in% isolates) done.isolates[isolates == g] <- 1
if(pick.1 > start && pick.1 <= end) {
# an isolate can join any other group (<smax), except other isolates done before (because already done) and this one (that would not change anything)
# if we are here, this implies that a partition with one less group is allowed
if(g %in% isolates){
previousisolates <- isolates[done.isolates == 1]
tosample <- all.groups[(all.groups %in% belowmax_groups) & !(all.groups %in% previousisolates)]
if(length(tosample) == 1) newg <- tosample
if(length(tosample) >= 2) newg <- sample(tosample,1)
new.partition[i] <- newg
new.partition <- order_groupids(new.partition)
}
# a member of a pair can join any other group (<smax), and splitting if it's not already counted, if isolates are allowed
if(g %in% pairs && g %in% abovemin_groups){
if((num.groups+1) %in% numgroups.simulated && done.pairs[pairs == g] == 0){ # if a split of the pair is possible
tosample <- all.groups[(all.groups %in% belowmax_groups) | (all.groups == g)]
if(length(tosample) == 1) newg <- tosample
if(length(tosample) >= 2) newg <- sample(tosample,1)
if(newg == g){
newg <- num.groups + 1
}
} else {
tosample <- all.groups[(all.groups %in% belowmax_groups) & (all.groups != g)]
if(length(tosample) == 1) newg <- tosample
if(length(tosample) >= 2) newg <- sample(tosample,1)
}
new.partition[i] <- newg
new.partition <- order_groupids(new.partition)
}
# a member of a bigger group (>smin) can join any other group (<smax) or create its own isolate (if it's allowed)
if(g %in% others && g %in% abovemin_groups){
if(smin == 1 && (num.groups+1) %in% numgroups.simulated) { # if we can create a isolate
tosample <- all.groups[(all.groups %in% belowmax_groups) | (all.groups == g)]
if(length(tosample) == 1) newg <- tosample
if(length(tosample) >= 2) newg <- sample(tosample,1)
if(newg == g){
newg <- num.groups + 1
}
} else {
tosample <- all.groups[(all.groups %in% belowmax_groups) & (all.groups != g)]
if(length(tosample) == 1) newg <- tosample
if(length(tosample) >= 2) newg <- sample(tosample,1)
}
new.partition[i] <- newg
new.partition <- order_groupids(new.partition)
}
}
if(g %in% pairs) done.pairs[pairs == g] <- 1
}
return(new.partition)
}
### Deprecated: unused neighborhoods
### Can be reintroduced, but they don't work with constraints on number of groups (only on sizes)
# ## NEIGHBORHOOD PI 4: swaps a pair of nodes
#
# compute_size_neighborhood_p4 <- function(partition){
#
# # find isolates, pairs and groups>2
# num.groups <- max(partition)
# num.nodes <- length(partition)
# sizes <- table(partition)
# isolates <- as.vector(which(sizes == 1))
# pairs <- as.vector(which(sizes == 2))
# others <- as.vector(which(sizes > 2))
#
# nums.swaps <- rep(0,max(partition))
#
# for(g in 1:num.groups){
#
# #if(g %in% isolates){}# nothing happens, there is no pair to swap here
#
# if(g %in% pairs){
# nums.swaps[g] <- length(pairs[pairs>g]) + length(isolates) + length(others)
# }
# if(g %in% others){
# sizeg <- sum(partition == g)
# nums.swaps[g] <- (length(pairs) + length(others) + length(isolates) ) * sizeg * (sizeg-1) / 2
# }
# }
#
# return(list(isolates = isolates,
# pairs = pairs,
# others = others,
# nums.swaps = nums.swaps,
# num.swaps = sum(nums.swaps),
# total = sum(nums.swaps)))
#
# }
#
#
# sample_new_partition_p4 <- function(current.partition, size_neighborhood){
#
# # calculate the number of neighbor partitions
# num.nodes <- length(current.partition)
# num.groups <- max(current.partition)
# size <- size_neighborhood
# isolates <- size$isolates
# pairs <- size$pairs
# others <- size$others
# nums.swaps <- size$nums.swaps
# num.swaps <- size$num.swaps
# total <- size$total
#
# if(total == 0) return(current.partition)
#
# pick.1 <- sample(total,1)
#
# # decide which group to pick a pair from
# all.groups <- 1:num.groups
# new.partition <- current.partition
# #done.pairs <- rep(0,length(pairs))
#
# sums <- unlist(lapply(1:num.groups,FUN = function(x){sum(nums.swaps[1:x])}))
# starts <- c(0,sums[1:(num.groups-1)])
# ends <- sums[1:num.groups]
# old_g <- which(starts < pick.1 & ends >= pick.1)
#
#
# # a pair can join other available groups, unless it's a pair that has already been "joined"
# if(old_g %in% pairs){
# members <- which(current.partition == old_g)
# tosample <- all.groups[!(all.groups %in% pairs & all.groups <= old_g)]
# if(length(tosample) == 1) newg <- tosample else newg <- sample(tosample,1)
# new.partition[members] <- newg
# new.partition <- order_groupids(new.partition)
# }
#
# # a pair in a bigger group can join any other group or create its own isolated pair
# if(old_g %in% others){
# members <- which(current.partition == old_g)
# sizeg <- length(members)
# pairg <- sample(members,2)
# tosample <- all.groups
# if(length(tosample) == 1) newg <- tosample else newg <- sample(tosample,1)
# if(newg == old_g) newg <- num.groups + 1
# new.partition[pairg] <- newg
# new.partition <- order_groupids(new.partition)
# }
#
# return(new.partition)
#
# }
#
#
# reachable_p4 <- function(partition1,partition2){
#
# # they are not neighbors if one doesn't have the same number of groups or one more
# if(abs(max(partition1) - max(partition2))>1) return(FALSE)
#
# num.nodes <- length(partition1)
# check <- FALSE
# i <- 1
# j <- i+1
#
# # try remove all pairs of nodes
# while(!check && i <= num.nodes){
#
# if(partition1[i] == partition1[j] && partition2[i] == partition2[j]) {
# newp1 <- order_groupids( partition1[-c(i,j)] )
# newp2 <- order_groupids( partition2[-c(i,j)] )
#
# # check if the pair moved together
# if(all(newp1 == newp2)) {
# others1 <- which(partition1[-c(i,j)] == partition1[i])
# others2 <- which(partition2[-c(i,j)] == partition2[i])
# check <- !setequal(others1,others2)
# if(check) break
# }
# }
#
# j <- j+1
# if(j > num.nodes){
# i <- i+1
# j <- i+1
# }
# }
#
# return(check)
# }
#
#
# ## NEIGHBORHOOD PI 5: swap two pairs of nodes
#
# compute_size_neighborhood_p5 <- function(partition){
#
# # find isolates, pairs and groups>2
# num.groups <- max(partition)
# num.nodes <- length(partition)
# sizes <- table(partition)
# isolates <- as.vector(which(sizes == 1))
# pairs <- as.vector(which(sizes == 2))
# others <- as.vector(which(sizes > 2))
#
# aff <- as.matrix(table(data.frame(actor = 1:num.nodes, group = partition)))
# net <- aff %*% t(aff)
# net[upper.tri(net)] <- 0
# diag(net) <- 0
# paired.actors <- which(net == 1)
# num.paired.actors <- sum(net)
#
# groups.paired.actors <- rep(0,num.paired.actors)
# pairs.paired.actors <- rep(0,num.paired.actors)
#
# impossible.swaps <- 0
# for(g in 1:num.groups){
# members <- which(partition == g)
#
# if(g %in% pairs){
# i <- members[1]
# j <- members[2]
# p <- (i-1)*num.nodes+j
# pairs.paired.actors[which(paired.actors == p)] <- 1
# groups.paired.actors[which(paired.actors == p)] <- g
# }
#
# if(g %in% others){
# for(i in 1:(length(members)-1)){
# for(j in (i+1):length(members)){
# p <- (members[i]-1)*num.nodes+members[j]
# groups.paired.actors[which(paired.actors == p)] <- g
# }
# }
# num.pairs.g <- length(members)*(length(members)-1)/2
# impossible.swaps <- impossible.swaps + num.pairs.g*(num.pairs.g-1)/2
# }
#
# }
#
# num.swaps <- num.paired.actors * (num.paired.actors-1) / 2 - impossible.swaps - sum(pairs.paired.actors) * (sum(pairs.paired.actors)-1) / 2
#
# return(list(isolates = isolates,
# pairs = pairs,
# others = others,
# paired.actors = paired.actors,
# groups.paired.actors = groups.paired.actors,
# pairs.paired.actors = pairs.paired.actors,
# num.swaps = num.swaps,
# total = num.swaps))
#
# }
#
#
# sample_new_partition_p5 <- function(current.partition, size_neighborhood){
#
# # calculate the number of neighbor partitions
# num.nodes <- length(current.partition)
# num.groups <- max(current.partition)
# size <- size_neighborhood
# isolates <- size$isolates
# pairs <- size$pairs
# others <- size$others
# paired.actors <- size$paired.actors
# groups.paired.actors <- size$groups.paired.actors
# pairs.paired.actors <- size$pairs.paired.actors
# num.swaps <- size$num.swaps
# total <- size$total
#
# if(total == 0) return(current.partition)
#
# new.partition <- current.partition
#
# # decide which actors to swap
# if(num.swaps > 0){
# found <- FALSE
# while(!found){
# pick.1 <- sample(paired.actors,1)
# pick.2 <- sample(paired.actors,1)
# if(groups.paired.actors[which(paired.actors == pick.1)] != groups.paired.actors[which(paired.actors == pick.2)] &&
# !(pairs.paired.actors[which(paired.actors == pick.1)] && pairs.paired.actors[which(paired.actors == pick.2)])) {
# found <- TRUE
# }
# }
#
# j1 <- ((pick.1-1) %/% num.nodes) + 1
# i1 <- (pick.1-1) %% num.nodes + 1
# j2 <- ((pick.2-1) %/% num.nodes) + 1
# i2 <- (pick.2-1) %% num.nodes + 1
# new.partition[i1] <- current.partition[i2]
# new.partition[j1] <- current.partition[j2]
# new.partition[i2] <- current.partition[i1]
# new.partition[j2] <- current.partition[j1]
# new.partition <- order_groupids(new.partition)
# }
#
#
# return(new.partition)
#
# }
#
#
# reachable_p5 <- function(partition1,partition2){
#
# # they are not neighbors if they have a different number of groups
# if(max(partition1) != max(partition2)) return(FALSE)
#
# num.nodes <- length(partition1)
# check <- FALSE
# i <- 1
# j <- i+1
# k <- 1
# l <- k+1
#
# # try remove all pairs of pairs
# while(!check && i <= num.nodes){
#
# if(partition1[i] == partition1[j] && partition2[i] == partition2[j] &&
# partition1[k] == partition1[l] && partition2[k] == partition2[l]) {
# newp1 <- order_groupids( partition1[-c(i,j,k,l)] )
# newp2 <- order_groupids( partition2[-c(i,j,k,l)] )
#
# # check if the pairs were swapped
# if(all(newp1 == newp2)) {
# othersij1 <- which(partition1[-c(i,j)] == partition1[i])
# othersij2 <- which(partition2[-c(i,j)] == partition2[i])
# otherskl1 <- which(partition1[-c(k,l)] == partition1[k])
# otherskl2 <- which(partition2[-c(k,l)] == partition2[k])
# check <- (setequal(othersij1,otherskl2) && setequal(othersij2,otherskl1))
# if(check) break
# }
# }
#
# l <- l+1
# if(l > num.nodes){
# k <- k+1
# l <- k+1
# if(k > num.nodes){
# j <- j+1
# if(j > num.nodes){
# i <- i+1
# j <- i+1
# }
# }
# }
# }
#
# return(check)
# }
#
#
# ## NEIGHBORHOOD PI 6: pick two groups, re attribute nodes within them (careful, cannot be used alone)
#
# compute_size_neighborhood_p6 <- function(partition){
#
# # find isolates, pairs and groups>2
# num.groups <- max(partition)
# num.nodes <- length(partition)
# sizes <- table(partition)
# isolates <- as.vector(which(sizes == 1))
# pairs <- as.vector(which(sizes == 2))
# others <- as.vector(which(sizes > 2))
#
# group.pairs <- matrix(0,num.groups,num.groups)
# nums.swaps <- matrix(0,num.groups,num.groups)
#
# if(num.groups > 1) {
#
# for(g1 in 1:(num.groups-1)){
# for(g2 in (g1+1):num.groups){
# size.1 <- sum(partition == g1)
# size.2 <- sum(partition == g2)
#
# if(size.1 == 1 && size.2 == 1){
# nums.swaps[g1,g2] <- 0
# } else {
# k <- size.1
# n <- size.1 + size.2
# nums.swaps[g1,g2] <- factorial(n) / (factorial(k)*factorial(n-k))
# if(size.1 == size.2) nums.swaps[g1,g2] <- nums.swaps[g1,g2] /2
# nums.swaps[g1,g2] <- nums.swaps[g1,g2] -1
# }
# }
# }
# }
#
# return(list(isolates = isolates,
# pairs = pairs,
# others = others,
# nums.swaps = nums.swaps,
# num.swaps = sum(nums.swaps),
# total = sum(nums.swaps)))
#
# }
#
#
# sample_new_partition_p6 <- function(current.partition, size_neighborhood){
#
# # calculate the number of neighbor partitions
# num.nodes <- length(current.partition)
# num.groups <- max(current.partition)
# size <- size_neighborhood
# isolates <- size$isolates
# pairs <- size$pairs
# others <- size$others
# nums.swaps <- size$nums.swaps
# num.swaps <- size$num.swaps
# total <- size$total
#
# if(total == 0) return(current.partition)
#
# # decide which actors to swap
# new.partition <- current.partition
# pick <- sample(1:num.swaps,1)
# start <- 0
# end <- 0
#
# if(num.groups > 1){
# for(g1 in 1:(num.groups-1)){
# for(g2 in (g1+1):num.groups){
#
# end <- end + nums.swaps[g1,g2]
#
# if(start < pick && pick <= end){
# found <- FALSE
# while(!found){
# tosample <- c(which(current.partition==g1),which(current.partition==g2))
# new.g1 <- sample(tosample,sum(current.partition==g1))
# new.g2 <- tosample[-which(tosample %in% new.g1)]
# if( !(setequal(new.g1,current.partition[g1]) && setequal(new.g2,current.partition[g2])) &&
# !(setequal(new.g1,current.partition[g2]) && setequal(new.g2,current.partition[g1])) ) found <- TRUE
# }
# new.partition[new.g1] <- g1
# new.partition[new.g2] <- g2
# new.partition <- order_groupids(new.partition)
# }
#
# start <- start + nums.swaps[g1,g2]
# }
# }
# }
#
# return(new.partition)
#
# }
#
#
# reachable_p6 <- function(partition1,partition2){
#
# # they are not neighbors if they have a different number of groups
# if(max(partition1) != max(partition2)) return(FALSE)
#
# num.nodes <- length(partition1)
# num.groups1 <- max(partition1)
# num.groups2 <- max(partition2)
# check <- FALSE
# g1 <- 1
# g2 <- g1+1
#
# # try remove groups (to check divisions) and pairs of groups (to check merges)
# while(!check && g <= num.groups1){
#
# # first check whether g1 is not divided in partition2
# membersg1 <- which(partition1 == g1)
# newp1 <- order_groupids( partition1[-membersg1] )
# newp2 <- order_groupids( partition2[-membersg1] )
#
# if(all(newp1 == newp2)) {
# check <- length(unique(partition2[membersg1])) == 2
# if(check) break
# }
#
# # second check whether g1 and g2 are merged in partition2
# membersg2 <- which(partition1 == g2)
# newp1 <- order_groupids( partition1[-c(membersg1,membersg2)] )
# newp2 <- order_groupids( partition2[-c(membersg1,membersg2)] )
#
# if(all(newp1 == newp2)) {
# check <- length(unique(partition2[c(membersg1,membersg2)])) == 1
# if(check) break
# }
#
# g2 <- g2+1
# if(g2 > num.nodes){
# g1 <- g1+1
# g2 <- g2+1
# }
# }
#
# return(check)
# }
#
# ## NEIGHBORHOOD PI 7: take two nodes together in partition 2, bring them together (single move)
# ## otherwise separate then
# ## careful need isolates
#
# compute_size_neighborhood_p7 <- function(partition, partition2){
#
# num.nodes <- length(partition)
# present <- !is.na(partition)
# present2 <- !is.na(partition2)
#
# num.groups <- max(partition, na.rm = TRUE)
# sizes <- table(partition)
# isolates <- as.vector(which(sizes == 1))
# pairs <- as.vector(which(sizes == 2))
#
# affiliation <- as.matrix(table(data.frame(actor = 1:num.nodes, group= partition)))
# affiliation2 <- as.matrix(table(data.frame(actor = 1:num.nodes, group= partition2)))
# adjacency <- affiliation %*% t(affiliation)
# adjacency2 <- affiliation2 %*% t(affiliation2)
# diag(adjacency) <- 0
# diag(adjacency2) <- 0
#
# nums.swaps <- matrix(0,num.nodes,num.nodes)
# for(i in 1:num.nodes) {
# for(j in 1:num.nodes){
# ad <- adjacency[i,j]
# ad2 <- adjacency2[i,j]
# if(!is.na(ad) && !is.na(ad2) && ad2 == 1) {
# if(ad == 0) { # if separated join if possible (two isolates will only join once)
# if(partition[i] %in% isolates && partition[j] %in% isolates && j < i) {
# nums.swaps[i,j] <- 0
# } else {
# nums.swaps[i,j] <- 1
# }
# }
# if(ad == 1) { # if together separate if possible (a pair will only separate once)
# if(partition[i] %in% pairs && j < i) {
# nums.swaps[i,j] <- num.groups - 1
# } else {
# nums.swaps[i,j] <- num.groups
# }
#
# }
# }
# }
# }
#
# return(list(adjacency = adjacency,
# pairs = pairs,
# nums.swaps = nums.swaps,
# num.swaps = sum(nums.swaps),
# total = sum(nums.swaps)))
# }
#
#
# sample_new_partition_p7 <- function(current.partition, current.partition2, size_neighborhood){
#
# num.nodes <- length(current.partition)
# num.groups <- max(current.partition,na.rm = TRUE)
# adjacency <- size_neighborhood$adjacency
# pairs <- size_neighborhood$pairs
# nums.swaps <- size_neighborhood$nums.swaps
# num.swaps <- size_neighborhood$num.swaps
# total <- size$total
#
# if(total == 0) return(current.partition1)
#
# present <- !(is.na(current.partition))
# new.partition <- current.partition
#
# probas_pairs <- as.vector(nums.swaps) / num.swaps
# pick <- sample(1:(num.nodes*num.nodes),1,prob=probas_pairs)
# i <- (pick-1) %% num.nodes + 1
# j <- ((pick-1) %/% num.nodes) + 1
#
# if(adjacency[i,j] == 0) {
# new.partition[i] <- current.partition[j]
# } else {
# if(current.partition[i] %in% pairs && j < i) {
# tosample <- (1:num.groups)[(1:num.groups) != current.partition[i]]
# if(length(tosample) == 1) pickgroup <- tosample
# else pickgroup <- sample(tosample,1)
# if(pickgroup == current.partition[i]) new.partition[i] <- num.groups + 1
# else new.partition[i] <- pickgroup
# }else {
# tosample <- 1:num.groups
# if(length(tosample) == 1) pickgroup <- tosample
# else pickgroup <- sample(tosample,1)
# if(pickgroup == current.partition[i]) new.partition[i] <- num.groups + 1
# else new.partition[i] <- pickgroup
# }
# }
# new.partition[present] <- order_groupids(new.partition[present])
#
# return(new.partition)
#
# }
#
#
#
# ## NEIGHBORHOOD PI 4 RESTRICTED: swaps a pair of nodes
# # WARNING!!!!: for now it only works if sizes allowed are an interval ([size_min,size_max])
# compute_size_neighborhood_p4_restricted <- function(partition, sizes.simulated){
#
# # check if current partition is allowed
# if(!check_sizes(partition, sizes.simulated)) {
# print(partition)
# #stop("The partition we are in is not allowed.")
# }
#
# # find maximum size allowed
# smax <- max(sizes.simulated)
# smin <- min(sizes.simulated)
#
# # find isolates, pairs and groups>2
# num.groups <- max(partition)
# num.nodes <- length(partition)
# sizes <- table(partition)
# isolates <- as.vector(which(sizes == 1))
# pairs <- as.vector(which(sizes == 2))
# others <- as.vector(which(sizes > 2))
# abovemin_groups <- as.vector(which(sizes > (smin+1)))
# belowmax_groups <- as.vector(which(sizes < (smax-1)))
#
# nums.swaps <- rep(0,max(partition))
# done.pairs <- rep(0,length(pairs))
# for(g in 1:num.groups){
#
# if(g %in% isolates){
# # nothing happens, there is no pair to swap here
# }
# if(g %in% pairs){
# # then the pair can move to a group <smax-1, unless it's another pair and the combination of the two pairs is already counted
# # is 2 <smax-1 we have to remove the pair from the count of possible other groups
# if(2 < smax-1) nums.swaps[g] <- length(belowmax_groups) - 1 - sum(done.pairs) + done.pairs[pairs==g]
# if(2 >= smax-1) nums.swaps[g] <- length(belowmax_groups) # then it's only isolates
# done.pairs[pairs == g] <- 1
# }
# if(g %in% others && g %in% abovemin_groups){
# # then any pair inside the group can be separated and added to an available group that is <smax-1, unless the group is not >smin+1
# sizeg <- sum(partition == g)
# gbelowmax <- (g %in% belowmax_groups)
# # if it can create its own isolated pair
# if(gbelowmax && 2 %in% sizes.simulated) nums.swaps[g] <- length(belowmax_groups) * sizeg * (sizeg-1) / 2
# if(!gbelowmax && 2 %in% sizes.simulated) nums.swaps[g] <- (length(belowmax_groups) + 1) * sizeg * (sizeg-1) / 2
# # if not
# if(gbelowmax && !(2 %in% sizes.simulated)) nums.swaps[g] <- (length(belowmax_groups)-1) * sizeg * (sizeg-1) / 2
# if(!gbelowmax && !(2 %in% sizes.simulated)) nums.swaps[g] <- length(belowmax_groups) * sizeg * (sizeg-1) / 2
# } # we can never break a group that is already of minimal size
# }
#
# return(list(isolates = isolates,
# pairs = pairs,
# others = others,
# abovemin_groups = abovemin_groups,
# belowmax_groups = belowmax_groups,
# nums.swaps = nums.swaps,
# num.swaps = sum(nums.swaps),
# total = sum(nums.swaps)))
#
# }
#
# sample_new_partition_p4_restricted <- function(current.partition, size_neighborhood, sizes.simulated){
#
# # calculate the number of neighbor partitions
# num.nodes <- length(current.partition)
# num.groups <- max(current.partition)
# size <- size_neighborhood
# isolates <- size$isolates
# pairs <- size$pairs
# others <- size$others
# abovemin_groups <- size$abovemin_groups
# belowmax_groups <- size$belowmax_groups
# nums.swaps <- size$nums.swaps
# num.swaps <- size$num.swaps
# total <- size$total
#
# # allowed sizes
# smax <- max(sizes.simulated)
# smin <- min(sizes.simulated)
#
# pick.1 <- sample(total,1)
#
# # decide which group to pick a pair from
# all.groups <- 1:num.groups
# new.partition <- current.partition
# done.pairs <- rep(0,length(pairs))
#
# for(g in 1:num.groups){
#
# if(g == 1) {start <- 0} else {start <- sum(nums.swaps[1:g-1])}
# if(g == num.groups) {end <- num.swaps} else {end <- sum(nums.swaps[1:g])}
#
# if(pick.1 > start && pick.1 <= end) {
#
# # a pair can join other available groups (<smax-1), unless it's a pair that has already been "joined"
# if(g %in% pairs){
# members <- which(current.partition == g)
# previous.pairs <- pairs[done.pairs == 1]
# if(2 < smax-1){
# tosample <- all.groups[(all.groups %in% belowmax_groups) & (all.groups != g) & !(all.groups %in% previous.pairs)]
# if(length(tosample) == 1) newg <- tosample
# if(length(tosample) >= 2) newg <- sample(tosample,1)
# }else{
# tosample <- all.groups[(all.groups %in% belowmax_groups)]
# if(length(tosample) == 1) newg <- tosample
# if(length(tosample) >= 2) newg <- sample(tosample,1)
# }
# new.partition[members] <- newg
# new.partition <- order_groupids(new.partition)
# }
#
# # a pair in a bigger group (>smin+1) can join any other group (<smax-1) or create its isolate
# if(g %in% others && g %in% abovemin_groups){
# members <- which(current.partition == g)
# sizeg <- length(members)
# pairg <- sample(members,2)
# if(2 %in% sizes.simulated) tosample <- all.groups[(all.groups %in% belowmax_groups) || (all.groups == g)]
# else tosample <- all.groups[(all.groups %in% belowmax_groups) & (all.groups != g)]
# if(length(tosample) == 1) newg <- tosample
# if(length(tosample) >= 2) newg <- sample(tosample,1)
# if(newg == g) newg <- num.groups + 1
# new.partition[pairg] <- newg
# new.partition <- order_groupids(new.partition)
# }
#
# }
#
# if(g %in% pairs) done.pairs[pairs == g] <- 1
# }
#
# return(new.partition)
#
# }
#
#
# ## NEIGHBORHOOD PI 5 RESTRICTED: only swaps of two pairs of nodes (careful, cannot be used alone)
# compute_size_neighborhood_p5_restricted <- function(partition, sizes.simulated){
# # check if current partition is allowed
# if(!check_sizes(partition, sizes.simulated)) stop("The partition we are in is not allowed.")
# return(compute_size_neighborhood_p5(partition))
# }
#
# sample_new_partition_p5_restricted <- function(current.partition, size_neighborhood, sizes.simulated){
# return(sample_new_partition_p5(current.partition, size_neighborhood))
# }
#
#
# ## NEIGHBORHOOD PI 6 RESTRICTED: exchanges of two groups (careful, cannot be used alone)
# compute_size_neighborhood_p6_restricted <- function(partition, sizes.simulated){
# # check if current partition is allowed
# if(!check_sizes(partition, sizes.simulated)) stop("The partition we are in is not allowed.")
# return(compute_size_neighborhood_p6(partition))
# }
#
# sample_new_partition_p6_restricted <- function(current.partition, size_neighborhood, sizes.simulated){
# return(sample_new_partition_p6(current.partition, size_neighborhood))
# }
#
#
# ## NEIGHBORHOOD PI 7 RESTRICTED: take two nodes together in partition 2, bring them together (single move)
# ## otherwise separate then
# ## careful need isolates
# compute_size_neighborhood_p7_restricted <- function(partition, partition2, sizes.simulated){
#
# # find maximum size allowed
# smax <- max(sizes.simulated)
# smin <- min(sizes.simulated)
#
# num.nodes <- length(partition)
# present <- !is.na(partition)
# present2 <- !is.na(partition2)
#
# num.groups <- max(partition,na.rm = TRUE)
# sizes <- table(partition)
# isolates <- as.vector(which(sizes == 1))
# pairs <- as.vector(which(sizes == 2))
# abovemin_groups <- as.vector(which(sizes > smin))
# belowmax_groups <- as.vector(which(sizes < smax))
# max_groups <- as.vector(which(sizes == smax))
#
# affiliation <- as.matrix(table(data.frame(actor = 1:num.nodes, group= partition)))
# affiliation2 <- as.matrix(table(data.frame(actor = 1:num.nodes, group= partition2)))
# adjacency <- affiliation %*% t(affiliation)
# adjacency2 <- affiliation2 %*% t(affiliation2)
# diag(adjacency) <- 0
# diag(adjacency2) <- 0
#
# nums.swaps <- matrix(0,num.nodes,num.nodes)
# for(i in 1:num.nodes) {
# for(j in 1:num.nodes){
# ad <- adjacency[i,j]
# ad2 <- adjacency2[i,j]
# if(!is.na(ad) && !is.na(ad2) && ad2 == 1) {
# if(ad == 0) { # if separated join if possible (two isolates will only join once)
# if(partition[i] %in% isolates && partition[j] %in% isolates && j < i) {
# nums.swaps[i,j] <- 0
# } else {
# nums.swaps[i,j] <- (partition[i] %in% c(abovemin_groups,isolates) ) & (partition[j] %in% belowmax_groups)
# }
# }
# if(ad == 1) { # if together separate if possible (a pair will only separate once)
# if(partition[i] %in% abovemin_groups) {
# if(partition[i] %in% pairs && j < i){
# nums.swaps[i,j] <- num.groups - length(max_groups) + (partition[i] %in% max_groups) - !(1 %in% sizes.simulated) - 1
# } else {
# nums.swaps[i,j] <- num.groups - length(max_groups) + (partition[i] %in% max_groups) - !(1 %in% sizes.simulated)
# }
# }
# }
# }
# }
# }
#
# return(list(abovemin_groups = abovemin_groups,
# belowmax_groups = belowmax_groups,
# adjacency = adjacency,
# pairs = pairs,
# nums.swaps = nums.swaps,
# num.swaps = sum(nums.swaps),
# total = sum(nums.swaps)))
# }
#
# sample_new_partition_p7_restricted <- function(current.partition, size_neighborhood, current.partition2, sizes.simulated){
#
# # find maximum size allowed
# smax <- max(sizes.simulated)
# smin <- min(sizes.simulated)
#
# num.nodes <- length(current.partition)
# num.groups <- max(current.partition,na.rm = TRUE)
# abovemin_groups <- size_neighborhood$abovemin_groups
# belowmax_groups <- size_neighborhood$belowmax_groups
# adjacency <- size_neighborhood$adjacency
# pairs <- size_neighborhood$pairs
# nums.swaps <- size_neighborhood$nums.swaps
# num.swaps <- size_neighborhood$num.swaps
# total <- size$total
#
# if(total == 0) return(current.partition)
#
# present <- !(is.na(current.partition))
# new.partition <- current.partition
#
# probas_pairs <- as.vector(nums.swaps) / num.swaps
# pick <- sample(1:(num.nodes*num.nodes),1,prob=probas_pairs)
# i <- (pick-1) %% num.nodes + 1
# j <- ((pick-1) %/% num.nodes) + 1
#
# if(adjacency[i,j] == 0) {
# new.partition[i] <- current.partition[j]
# }
# if(adjacency[i,j] == 1) {
# if(current.partition[i] %in% pairs && j < i){
# tosample <- belowmax_groups[belowmax_groups != current.partition[i]]
# } else {
# if(1 %in% sizes.simulated) tosample <- unique(c(belowmax_groups,current.partition[i]))
# else tosample <- belowmax_groups[belowmax_groups != current.partition[i]]
# }
# if(length(tosample) == 1) pickgroup <- tosample
# else pickgroup <- sample(tosample,1)
# if(pickgroup == current.partition[i]) new.partition[i] <- num.groups + 1
# else new.partition[i] <- pickgroup
# }
# new.partition[present] <- order_groupids(new.partition[present])
#
# return(new.partition)
#
# }
#
#
##############################################################################
### Deprecated 2: neighborhoods for the multiple partition estimation
# ## NEIGHBORHOOD PI 7: pick two partitions and exchange then (only for multiple partition estimation)
# # since we force one change, and it's reversible (bijection), we don't have to care about the size of the neighborhood
# # CAREFUL, only works when singletons are allowed
# compute_size_neighborhood_p7 <- function(partition1, partition2){
#
# num.nodes <- length(partition1)
# present1 <- !is.na(partition1)
# present2 <- !is.na(partition2)
#
# # calculate the number of ways to reconstruct p1 after having done the exchange with p2
# p1withoutabsentp2 <- partition1[present1 & present2]
# max1 <- length(unique(p1withoutabsentp2))
# nodesunknown1 <- sum(present1 & !present2)
# if(nodesunknown1 == 0) {
# possibilities1 <- 1
# } else {
# possibilities1 <- 0
# for(k in 1:nodesunknown1) {
# possibilities1 <- possibilities1 + Stirling(nodesunknown1,k)*(max1+1)^k
# }
# }
#
# # same the other way around
# p2withoutabsentp1 <- partition2[present2 & present1]
# max2 <- length(unique(p2withoutabsentp1))
# nodesunknown2 <- sum(present2 & !present1)
# if(nodesunknown2 == 0) {
# possibilities2 <- 1
# } else {
# possibilities2 <- 0
# for(k in 1:nodesunknown2) {
# possibilities2 <- possibilities2 + Stirling(nodesunknown2,k)*(max2+1)^k
# }
# }
#
# return(possibilities1 * possibilities2)
# }
#
# sample_new_partition_p7 <- function(current.partition1, current.partition2){
#
# num.nodes <- length(current.partition1)
# present1 <- which(!is.na(current.partition1))
# present2 <- which(!is.na(current.partition2))
#
# new.partition1 <- rep(0,num.nodes)
# new.partition2 <- rep(0,num.nodes)
#
# # fill in the first with the second, and create isolates for people not present in the second
# new.partition1[present1] <- current.partition2[present1]
# max1 <- max(new.partition1, na.rm = TRUE)
# singletons1 <- which(is.na(new.partition1))
# if(length(singletons1) > 0) { # assign randomly the ones who are not present in the other one
# for(s in singletons1){
# newgroup <- sample(1:(max1+1),1)
# new.partition1[s] <- newgroup
# if(newgroup == (max1+1)) max1 <- max1 + 1
# }
# }
# new.partition1[present1] <- order_groupids(new.partition1[present1])
# new.partition1[new.partition1 == 0] <- NA
#
# # same the other way around
# new.partition2[present2] <- current.partition1[present2]
# max2 <- max(new.partition2, na.rm = TRUE)
# singletons2 <- which(is.na(new.partition2))
# if(length(singletons2) > 0) { # assign randomly the ones who are not present in the other one
# for(s in singletons2){
# newgroup <- sample(1:(max2+1),1)
# new.partition2[s] <- newgroup
# if(newgroup == (max2+1)) max2 <- max2 + 1
# }
# }
# new.partition2[present2] <- order_groupids(new.partition2[present2])
# new.partition2[new.partition2 == 0] <- NA
#
# return(list(new.partition1 = new.partition1,
# new.partition2 = new.partition2))
#
# }
#
#
# ## NEIGHBORHOOD PI 8: pick another partition, a group in it, and reproduce it in the current partition, or reorganizes it
# # CAREFUL, only works when singletons are allowed
# compute_size_neighborhood_p8 <- function(partition1, partition2){
#
# num.nodes <- length(partition1)
# present1 <- !is.na(partition1)
# present2 <- !is.na(partition2)
#
# # count number of groups from the second partition
# groups.toreorganize <- unique(partition2[present1 & present2])
# nums.reorganizations <- rep(0,length(groups.toreorganize))
#
# # count for each of them number of moves
# for(g in 1:length(groups.toreorganize)){
# group <- groups.toreorganize[g]
# members <- which(partition2 == group)
# max1 <- length(unique(partition1[present1 & partition1 != group]))
# for(k in 1:length(members)) {
# nums.reorganizations[g] <- nums.reorganizations[g] + Stirling(length(members),k)*(max1+1)^k
# }
# nums.reorganizations[g] <- nums.reorganizations[g] + 1 # when we just copy the group
# }
#
# return(list(groups.toreorganize = groups.toreorganize,
# nums.reorganizations = nums.reorganizations,
# num.reorganizations = sum(nums.reorganizations)))
# }
#
# sample_new_partition_p8 <- function(current.partition1, current.partition2, size_neighborhood){
#
# groups.toreorganize <- size_neighborhood$groups.toreorganize
# nums.reorganizations <- size_neighborhood$nums.reorganizations
# num.reorganizations <- size_neighborhood$num.reorganizations
#
# num.nodes <- length(current.partition1)
# present1 <- !is.na(current.partition1)
# present2 <- !is.na(current.partition2)
#
# new.partition <- current.partition1
#
# # pick
# probas_groups <- nums.reorganizations / num.reorganizations
# pick <- sample(c(0,groups.toreorganize),1,prob=c(1/ num.reorganizations, probas_groups))
#
# #randomize group
# if(pick > 0) {
# group <- pick
# members <- which(current.partition2 == group && present1)
# max1 <- length(unique(current.partition1[present1 & current.partition1 != group]))
# for(a in members){
# newgroup <- sample(1:(max1+1),1)
# new.partition[a] <- newgroup
# if(newgroup == (max1+1)) max1 <- max1 + 1
# }
# new.partition[present1] <- order_groupids(new.partition[present1])
# }
#
# # copy group
# if(pick ==0) {
# # reproduce the group
# new.partition[members] <- max(new.partition,na.rm = TRUE) + 1
# #reorder
# new.partition[present1] <- order_groupids(new.partition[present1])
# }
#
# return(new.partition)
#
# }
#
# #sample_new_partition_p8 <- function(current.partition1, current.partition2){
# #
# # num.nodes <- length(current.partition1)
# # present1 <- !is.na(current.partition1)
# # present2 <- !is.na(current.partition2)
# #
# # new.partition <- current.partition1
# #
# # # pick a group to reproduce
# # topick <- unique(current.partition2[present1 & present2])
# # g <- sample(topick,1)
# # members <- which(current.partition2 == g && present1)
# #
# # # reproduce the group
# # new.partition[members] <- max(new.partition,na.rm = TRUE) + 1
# #
# # #reorder
# # new.partition[present1] <- order_groupids(new.partition[present1])
# #
# # return(new.partition)
# #
# #}
#
# ## NEIGHBORHOOD PI 8: pick another partition, a group in it, and reproduce it in the current partition, or reorganizes it
# # CAREFUL, only works when singletons are allowed
# compute_size_neighborhood_p8_restricted <- function(partition1, partition2, sizes.simulated){
#
# num.nodes <- length(partition1)
# present1 <- !is.na(partition1)
# present2 <- !is.na(partition2)
# smax <- max(sizes.simulated)
#
# # count number of groups from the second partition
# groups.toreorganize <- unique(partition2[present1 & present2])
# combinations.pergroup <- list()
# nums.partitions.pergroup <- list()
# nums.reorganizations.pergroup <- list()
# nums.reorganizations <- rep(0,length(groups.toreorganize))
#
# # count for each of them number of moves
# for(g in 1:length(groups.toreorganize)){
# group <- groups.toreorganize[g]
# members <- which(partition2 == group)
# members <- which(1:num.nodes %in% members & present1)
#
# # if the group is not present, we will create it
# if(length(unique(partition1[members])) == 1 && sum(partition1[-members] == (partition1[members])[1], na.rm = TRUE) == 0) {
#
# nums.reorganizations.pergroup[[g]] <- NULL
# nums.reorganizations[g] <- 1
#
# # if the group is present, we will reshuffle
# } else {
# nmembers <- length(members)
# sizes <- table(partition1[-members])
#
# # find ways to divide the group
# combs <- parts(nmembers)
# combinations.pergroup[[g]] <- combs
#
# nums.p <- rep(1,ncol(combs))
# nums.r <- rep(1,ncol(combs))
# for(c in 1:ncol(combs)){
#
# # find number partitions of each
# ncpt <- nmembers
# for(s in 1:nrow(combs)){
# ncs <- sum(combs[,c] == s)
# nums.p[c] <- nums.p[c] * (choose(ncpt,s))^ncs / factorial(ncs)
# ncpt <- ncpt - ncs*s
# }
#
# # find ways to reorganize
# for(s in 1:nrow(combs)){
# ncs <- sum(combs[,c] == s)
# ncpt <- sum(sizes <= (smax-s))
# for(k in 1:ncs){
# nums.r[c] <- nums.r[c] * ncpt
# ncpt <- ncpt - 1
# }
# }
# }
#
# nums.partitions.pergroup[[g]] <- nums.p
# nums.reorganizations.pergroup[[g]] <- nums.r
# nums.reorganizations[g] <- nums.p * nums.r
# }
#
# }
#
# return(list(groups.toreorganize = groups.toreorganize,
# combinations.pergroup = combinations.pergroup,
# nums.partitions.pergroup = nums.partitions.pergroup,
# nums.reorganizations.pergroup = nums.reorganizations.pergroup,
# nums.reorganizations = nums.reorganizations,
# num.reorganizations = sum(nums.reorganizations) ))
# }
#
# sample_new_partition_p8_restricted <- function(current.partition1, current.partition2, size_neighborhood, sizes.simulated){
#
# groups.toreorganize <- size_neighborhood$groups.toreorganize
# combinations.pergroup <- size_neighborhood$combinations.pergroup
# nums.partitions.pergroup <- size_neighborhood$nums.partitions.pergroup
# nums.reorganizations.pergroup <- size_neighborhood$nums.reorganizations.pergroup
# nums.reorganizations <- size_neighborhood$nums.reorganizations
# num.reorganizations <- size_neighborhood$num.reorganizations
#
# num.nodes <- length(current.partition1)
# present1 <- !is.na(current.partition1)
# present2 <- !is.na(current.partition2)
#
# new.partition <- current.partition1
#
# # pick
# probas_groups <- nums.reorganizations / num.reorganizations
# group <- sample(groups.toreorganize,1,prob=probas_groups)
# members <- which(current.partition2 == group)
# members <- which(1:num.nodes %in% members & present1)
#
# # if the group is not present, we will create it
# if(length(unique(current.partition1[members])) == 1 && sum(current.partition1[-members] == (current.partition1[members])[1], na.rm = TRUE) == 0) {
#
# new.partition[members] <- max(new.partition,na.rm = TRUE) + 1
# new.partition[present1] <- order_groupids(new.partition[present1])
#
# # if the group is present, we will reshuffle
# } else {
#
# nmembers <- length(members)
# gleft <- unique(current.partition1[-members])
# gleft <- gleft[!is.na(gleft)]
# sizes <- apply(gleft, FUN = function(x) {return(sum(current.partition1[-members] == x))} )
#
# g <- which(groups.toreorganize == group)
# combs <- combinations.pergroup[[g]]
#
# # sample a combination
# probas_combs <- nums.partitions.pergroup[[g]] * nums.reorganizations.pergroup[[g]] / nums.reorganizations[g]
# c <- sample(1:ncol(combs),1,prob=probas_combs)
#
# # create a combination
# ntaken <- rep(F,nmembers)
# gtaken <- rep(F,length(gleft))
# for(s in nrow(combs)) {
# ncs <- sum(combs[,c] == s)
# for(k in 1:ncs) {
# ns <- choose(members[!ntaken], s)
# topick <- gleft[sizes <= (smax-s) & !gtaken]
# newg <- sample(c(0,topick),1)
# if(newg == 0){
# new.partition[ns] <- max(new.partition,na.rm = TRUE) + 1
# } else {
# new.partition[ns] <- newg
# gtaken[gleft == newg] <- TRUE
# }
# ntaken[members == ns] <- TRUE
# }
# }
# newnodes <- sample(1:nmembers,nmembers)
#
# new.partition[present1] <- order_groupids(new.partition[present1])
#
# }
#
# return(new.partition)
#
# }
Try the ERPM package in your browser
Any scripts or data that you put into this service are public.
ERPM documentation built on May 29, 2024, 10:05 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions sample_new_partition_p3_restricted compute_size_neighborhood_p3_restricted sample_new_partition_p2_restricted compute_size_neighborhood_p2_restricted sample_new_partition_p1_restricted compute_size_neighborhood_p1_restricted reachable_p3 sample_new_partition_p3 compute_size_neighborhood_p3 reachable_p2 sample_new_partition_p2 compute_size_neighborhood_p2 reachable_p1 sample_new_partition_p1 compute_size_neighborhood_p1 reachable sample_new_partition compute_size_neighborhood step_recalculate recalculate_statistics draw_step_multiple draw_step_single draw_Metropolis_multiple draw_Metropolis_single
Documented in draw_Metropolis_multiple draw_Metropolis_single
######################################################################
## Simulation and estimation of Exponential Random Partition Models ##
## Function implementing the Metropolis HAstings algorithm to ##
## through partitions given a certain model specification ##
## Author: Marion Hoffman ##
######################################################################
## --- FUNCTIONS TO DRAW CHAINS ----
#' Draw Metropolis single
#'
#' Function to sample the model with a Markov chain (single partition procedure).
#'
#'
#' @param theta model parameters
#' @param first.partition, starting partition for the Markov chain
#' @param nodes nodeset (data frame)
#' @param effects effects/sufficient statistics (list with a vector "names", and a vector "objects")
#' @param objects objects used for statistics calculation (list with a vector "name", and a vector "object")
#' @param burnin integer for the number of burn-in steps before sampling
#' @param thining integer for the number of thining steps between sampling
#' @param num.steps number of samples
#' @param neighborhood = c(0.7,0.3,0), way of choosing partitions: probability vector (2 actors swap, merge/division, single actor move, single pair move, 2 pairs swap, 2 groups reshuffle)
#' @param numgroups.allowed = NULL, # vector containing the number of groups allowed in the partition (now, it only works with vectors like num_min:num_max)
#' @param numgroups.simulated = NULL, # vector containing the number of groups simulated
#' @param sizes.allowed = NULL, vector of group sizes allowed in sampling (now, it only works for vectors like size_min:size_max)
#' @param sizes.simulated = NULL, vector of group sizes allowed in the Markov chain but not necessraily sampled (now, it only works for vectors like size_min:size_max)
#' @param return.all.partitions = FALSE option to return the sampled partitions on top of their statistics (for GOF)
#' @return A list
#' @importFrom stats runif
#' @export
#' @examples
#' # define an arbitrary set of n = 6 nodes with attributes, and an arbitrary covariate matrix
#' n <- 6
#' nodes <- data.frame(label = c("A","B","C","D","E","F"),
#' gender = c(1,1,2,1,2,2),
#' age = c(20,22,25,30,30,31))
#' friendship <- matrix(c(0, 1, 1, 1, 0, 0,
#' 1, 0, 0, 0, 1, 0,
#' 1, 0, 0, 0, 1, 0,
#' 1, 0, 0, 0, 0, 0,
#' 0, 1, 1, 0, 0, 1,
#' 0, 0, 0, 0, 1, 0), 6, 6, TRUE)
#'
#' # choose the effects to be included (see manual for all effect names)
#' effects <- list(names = c("num_groups","same","diff","tie"),
#' objects = c("partition","gender","age","friendship"))
#' objects <- list()
#' objects[[1]] <- list(name = "friendship", object = friendship)
#'
#' # set parameter values for each of these effects
#' parameters <- c(-0.2, 0.2, -0.1, 0.5)
#'
#' \donttest{
#' # generate simulated sample, by setting the desired additional parameters for the
#' # Metropolis sampler and choosing a starting point for the chain (first.partition)
#' nsteps <- 100
#' sample <- draw_Metropolis_single(theta = parameters,
#' first.partition = c(1,1,2,2,3,3),
#' nodes = nodes,
#' effects = effects,
#' objects = objects,
#' burnin = 100,
#' thining = 10,
#' num.steps = nsteps,
#' neighborhood = c(0,1,0),
#' numgroups.allowed = 1:n,
#' numgroups.simulated = 1:n,
#' sizes.allowed = 1:n,
#' sizes.simulated = 1:n,
#' return.all.partitions = TRUE)
#'
#'
#' # or: simulate an estimated model
#' partition <- c(1,1,2,2,2,3) # the partition already defined for the (previous) estimation
#' nsimulations <- 1000
#' simulations <- draw_Metropolis_single(theta = estimation$results$est,
#' first.partition = partition,
#' nodes = nodes,
#' effects = effects,
#' objects = objects,
#' burnin = 100,
#' thining = 20,
#' num.steps = nsimulations,
#' neighborhood = c(0,1,0),
#' sizes.allowed = 1:n,
#' sizes.simulated = 1:n,
#' return.all.partitions = TRUE)
#' }
#'
draw_Metropolis_single <- function(theta,
first.partition,
nodes,
effects,
objects,
burnin,
thining,
num.steps,
neighborhood = c(0.7,0.3,0),
numgroups.allowed = NULL,
numgroups.simulated = NULL,
sizes.allowed = NULL,
sizes.simulated = NULL,
return.all.partitions = FALSE) {
num.nodes <- nrow(nodes)
num.effects <- length(effects$names)
# turn the neighborhood weights into probabilities if needed
if(sum(neighborhood) != 1){
neighborhood <- neighborhood / sum(neighborhood)
}
# check whether there are constraints
constraints <- FALSE
if(!is.null(numgroups.allowed) || !is.null(sizes.allowed)) {
constraints <- TRUE
if(is.null(numgroups.allowed)) numgroups.allowed <- 1:nrow(nodes)
if(is.null(numgroups.simulated)) numgroups.simulated <- numgroups.allowed
if(is.null(sizes.allowed)) sizes.allowed <- 1:nrow(nodes)
if(is.null(sizes.simulated)) sizes.simulated <- sizes.allowed
}
# instantiate with the starting network
current.partition <- first.partition
current.z <- computeStatistics(current.partition, nodes, effects, objects)
current.logit <- theta * current.z
# store the statistics collected for all networks simulated after the burn in
all.z <-c()
# store the partitions if needed (for GOF)
if(return.all.partitions){
all.partitions <- c()
}
end.walk <- FALSE
cpt_burnin <- 0
cpt_thining <- 0
cpt_steps <- 0
while(!end.walk){
new.step <- draw_step_single(theta,
current.partition,
current.logit,
current.z,
nodes,
effects,
objects,
neighborhood,
numgroups.simulated,
sizes.simulated)
proba.change <- min(1,new.step$hastings.ratio)
change.made <- (runif(1,0,1) <= proba.change)
# update partition if needed
old.partition <- current.partition
if(change.made) {
current.partition <- new.step$new.partition
current.z <- new.step$new.z
current.logit <- new.step$new.logit
}
cpt_burnin <- cpt_burnin + 1
if(cpt_burnin > burnin) cpt_thining <- cpt_thining + 1
# store the results if we are out of burnin
if(cpt_burnin >= burnin && cpt_thining == thining) {
store <- TRUE
if(constraints) store <- check_sizes(current.partition,sizes.allowed,numgroups.allowed)
if(store) {
all.z <- rbind(all.z,current.z)
cpt_thining <- 0
cpt_steps <- cpt_steps + 1
if(return.all.partitions){ # store all partitions if needed
all.partitions <- rbind(all.partitions,current.partition)
}
} else {
cpt_thining <- thining - 1
}
}
# stop the walk if number of steps reached
end.walk <- (cpt_steps >= num.steps)
}
# TODO: change this hack
row.names(all.z) <- rep("", dim(all.z)[1])
# compute the average statistics and the final network generated
if(!return.all.partitions) {
return( list("draws" = all.z,
"last.partition" = current.partition,
"all.partitions" = NULL))
} else {
return( list("draws" = all.z,
"last.partition" = current.partition,
"all.partitions" = all.partitions))
}
}
# MULTIPLE PARTITIONS PROCEDURE
#' Draw Metropolis multiple
#'
#' Function to sample the model with a Markov chain (single partition procedure).
#'
#' @param theta model parameters
#' @param first.partitions starting partition for the Markov chain
#' @param presence.tables matrix indicating which actors were present for each observations (mandatory)
#' @param nodes node set (data frame)
#' @param effects effects/sufficient statistics (list with a vector "names", and a vector "objects")
#' @param objects objects used for statistics calculation (list with a vector "name", and a vector "object")
#' @param burnin integer for the number of burn-in steps before sampling
#' @param thining integer for the number of thining steps between sampling
#' @param num.steps number of samples
#' @param neighborhood = c(0.7,0.3,0), way of choosing partitions: probability vector (2 actors swap, merge/division, single actor move, single pair move, 2 pairs swap, 2 groups reshuffle)
#' @param numgroups.allowed = NULL, # vector containing the number of groups allowed in the partition (now, it only works with vectors like num_min:num_max)
#' @param numgroups.simulated = NULL, # vector containing the number of groups simulated
#' @param sizes.allowed = NULL, vector of group sizes allowed in sampling (now, it only works for vectors like size_min:size_max)
#' @param sizes.simulated = NULL, vector of group sizes allowed in the Markov chain but not necessraily sampled (now, it only works for vectors like size_min:size_max)
#' @param return.all.partitions = FALSE, option to return the sampled partitions on top of their statistics (for GOF)
#' @param verbose logical: should intermediate results during the estimation be printed or not? Defaults to FALSE.
#' @return A list
#' @importFrom stats runif
#' @export
#' @examples
#' # define an arbitrary set of n = 6 nodes with attributes, and an arbitrary covariate matrix
#' n <- 6
#' nodes <- data.frame(label = c("A","B","C","D","E","F"),
#' gender = c(1,1,2,1,2,2),
#' age = c(20,22,25,30,30,31))
#' friendship <- matrix(c(0, 1, 1, 1, 0, 0,
#' 1, 0, 0, 0, 1, 0,
#' 1, 0, 0, 0, 1, 0,
#' 1, 0, 0, 0, 0, 0,
#' 0, 1, 1, 0, 0, 1,
#' 0, 0, 0, 0, 1, 0), 6, 6, TRUE)
#'
#' # specify whether nodes are present at different points of time
#' presence.tables <- matrix(c(1, 1, 1, 1, 1, 1,
#' 0, 1, 1, 1, 1, 1,
#' 1, 0, 1, 1, 1, 1), 6, 3)
#'
#' # choose effects to be included in the estimated model
#' effects_multiple <- list(names = c("num_groups","same","diff","tie","inertia_1"),
#' objects = c("partitions","gender","age","friendship","partitions"),
#' objects2 = c("","","","",""))
#' objects_multiple <- list()
#' objects_multiple[[1]] <- list(name = "friendship", object = friendship)
#'
#' # set parameter values for each of these effects
#' parameters <- c(-0.2,0.2,-0.1,0.5,1)
#'
#' # set a starting point for the simulation
#' first.partitions <- matrix(c(1, 1, 2, 2, 2, 3,
#' NA, 1, 1, 2, 2, 2,
#' 1, NA, 2, 3, 3, 1), 6, 3)
#'
#' \donttest{
#' # generate the simulated sample
#' nsteps <- 50
#' sample <- draw_Metropolis_multiple(theta = parameters,
#' first.partitions = first.partitions,
#' nodes = nodes,
#' presence.tables = presence.tables,
#' effects = effects_multiple,
#' objects = objects_multiple,
#' burnin = 100,
#' thining = 100,
#' num.steps = nsteps,
#' neighborhood = c(0,1,0),
#' numgroups.allowed = 1:n,
#' numgroups.simulated = 1:n,
#' sizes.allowed = 1:n,
#' sizes.simulated = 1:n,
#' return.all.partitions = TRUE)
#' }
#'
draw_Metropolis_multiple <- function(theta,
first.partitions,
presence.tables,
nodes,
effects,
objects,
burnin,
thining,
num.steps,
neighborhood = c(0.7,0.3,0),
numgroups.allowed,
numgroups.simulated,
sizes.allowed,
sizes.simulated,
return.all.partitions = FALSE,
verbose = FALSE)
{
num.nodes <- nrow(nodes)
num.effects <- length(effects$names)
num.obs <- ncol(presence.tables)
# turn the neighborhood weights into probabilities if needed
if(sum(neighborhood) != 1){
neighborhood <- neighborhood / sum(neighborhood)
}
# check whether there are constraints
constraints <- FALSE
if(!is.null(numgroups.allowed) || !is.null(sizes.allowed)) {
constraints <- TRUE
if(is.null(numgroups.allowed)) numgroups.allowed <- 1:nrow(nodes)
if(is.null(numgroups.simulated)) numgroups.simulated <- numgroups.allowed
if(is.null(sizes.allowed)) sizes.allowed <- 1:nrow(nodes)
if(is.null(sizes.simulated)) sizes.simulated <- sizes.allowed
}
# instantiate with the starting network
current.partitions <- first.partitions
current.z.contributions <- computeStatistics_multiple(current.partitions, presence.tables, nodes, effects, objects)
current.z <- rowSums(current.z.contributions)
current.logit <- theta * current.z
# store the statistics collected for all networks simulated after the burn in
all.z <-c()
# store the partitions if needed (for GOF)
if(return.all.partitions){
all.partitions <- list()
cpt_all <- 0
}
end.walk <- FALSE
cpt_burnin <- 0
cpt_thining <- 0
cpt_steps <- 0
while(!end.walk){
new.step <- draw_step_multiple(theta,
current.partitions,
current.logit,
current.z.contributions,
current.z,
presence.tables,
nodes,
effects,
objects,
neighborhood,
numgroups.simulated,
sizes.simulated)
proba.change <- min(1,new.step$hastings.ratio)
change.made <- (runif(1,0,1) <= proba.change)
move <- new.step$move
# update partition if needed
old.partitions <- current.partitions
if(change.made) {
current.partitions <- new.step$new.partitions
current.z.contributions <- new.step$new.z.contributions
current.z <- new.step$new.z
current.logit <- new.step$new.logit
}
cpt_burnin <- cpt_burnin + 1
if(cpt_burnin > burnin) cpt_thining <- cpt_thining + 1
if(cpt_burnin %% 10000 == 0 && verbose) cat(cpt_burnin, "\n")
# store the results if we are out of burnin
if(cpt_burnin >= burnin && cpt_thining == thining) {
store <- TRUE
if(constraints) { # check all partitions one by one
for(o in 1:num.obs) store <- store && check_sizes(current.partitions[as.logical(presence.tables[,o]),o],sizes.allowed,numgroups.allowed)
}
if(store) {
all.z <- rbind(all.z,current.z)
cpt_thining <- 0
cpt_steps <- cpt_steps+1
# store all partitions if needed
if(return.all.partitions){
all.partitions[[cpt_steps]] <- current.partitions
}
} else {
cpt_thining <- thining-1
}
}
# stop the walk if number of steps reached
end.walk <- (cpt_steps >= num.steps)
}
# TODO: change this hack
row.names(all.z) <- rep("", dim(all.z)[1])
# compute the average statistics and the final network generated
if(!return.all.partitions) {
return( list("draws" = all.z,
"last.partitions" = current.partitions,
"all.partitions" = NULL))
} else {
return( list("draws" = all.z,
"last.partitions" = current.partitions,
"all.partitions" = all.partitions))
}
}
# MULTIPLE PARTITIONS PROCEDURE
## --- FUNCTIONS TO DRAW ONE STEP IN THE CHAIN ----
# function to draw next partition and calculate HAstings ratio (one step in the Metropolis algorithm)
draw_step_single <- function(theta,
current.partition,
current.logit,
current.z,
nodes,
effects,
objects,
neighborhood,
numgroups.simulated,
sizes.simulated){
# pick a neighborhood
# 1 = swap two actors,
# 2 = merge 2 groups or split a group in two,
# 3 = move one actor,
## 4 = move one pair of nodes (FOR NOW REMOVED)
## 5 = swap two pairs of nodes (FOR NOW REMOVED)
## 6 = reshuffle the members of two groups (FOR NOW REMOVED)
n_neighborhood <- length(neighborhood)
move <- sample(1:n_neighborhood,1,prob=neighborhood)
current.sizes <- rep(0,n_neighborhood)
new.sizes <- rep(0,n_neighborhood)
# calculate current size of neighborhood and pick new partition for the chosen move
current.size <- compute_size_neighborhood(move, current.partition, numgroups.simulated, sizes.simulated)
if(current.size$total > 0) {
new.partition <- sample_new_partition(move, current.partition, current.size, numgroups.simulated, sizes.simulated)
new.size <- compute_size_neighborhood(move, new.partition, numgroups.simulated, sizes.simulated)
} else {
new.partition <- current.partition
new.size <- current.size
}
current.sizes[move] <- current.size$total
new.sizes[move] <- new.size$total
# if change, calculate sizes for all potential moves (depending on available neighborhoods)
if(current.size$total > 0) {
for(i in 1:n_neighborhood){
if(neighborhood[i] > 0 && i != move && reachable(i,current.partition,new.partition)) {
cs <- compute_size_neighborhood(i, current.partition, numgroups.simulated, sizes.simulated)
ns <- compute_size_neighborhood(i, new.partition, numgroups.simulated, sizes.simulated)
current.sizes[i] <- cs$total
new.sizes[i] <- ns$total
}
}
}
# compute new statistics only if it changed
if(current.size$total > 0) {
new.z <- computeStatistics(new.partition, nodes, effects, objects)
} else { new.z <- current.z }
new.logit <- theta * new.z
# calculate acceptance ratio if needed
if(current.size$total > 0) {
indexes <- current.sizes > 0
neighborhoods.ratio <- sum(neighborhood[indexes] / new.sizes[indexes]) /
sum(neighborhood[indexes] / current.sizes[indexes])
} else { neighborhoods.ratio <- 1 }
hastings.ratio <- (exp(sum(new.logit) - sum(current.logit))) * neighborhoods.ratio
return(list("hastings.ratio" = hastings.ratio,
"new.partition" = new.partition,
"new.z" = new.z,
"new.logit" = new.logit))
}
# function to draw next partition and calculate HAstings ratio (one step in the Metropolis algorithm)
draw_step_multiple <- function(theta,
current.partitions,
current.logit,
current.z.contributions,
current.z,
presence.tables,
nodes,
effects,
objects,
neighborhood,
numgroups.simulated,
sizes.simulated){
num.nodes <- nrow(nodes)
num.obs <- ncol(presence.tables)
new.partitions <- current.partitions
rand.o <- sample(1:num.obs,1)
nodes.rand.o <- as.logical(presence.tables[,rand.o])
current.partition <- current.partitions[nodes.rand.o,rand.o]
# calculate sizes for all potential moves (depending on available neighborhoods)
current.sizes <- rep(0,length(neighborhood))
# pick a neighborhood
# 1 = swap two actors,
# 2 = merge 2 groups or split a group in two,
# 3 = move one actor,
# 4 = move one pair of nodes (FOR NOW REMOVED)
# 5 = swap two pairs of nodes (FOR NOW REMOVED)
# 6 = reshuffle the members of two groups (FOR NOW REMOVED)
# 7 = reproduce a pair from the previous partition (empty for first time point) (FOR NOW REMOVED)
n_neighborhood <- length(neighborhood)
move <- sample(1:n_neighborhood,1,prob=neighborhood)
# calculate current size of neighborhood and pick new partition for neighborhoods 1 to 6
for(i in 1:n_neighborhood){
if(neighborhood[i] > 0){
current.size <- compute_size_neighborhood(i, current.partition, numgroups.simulated, sizes.simulated)
current.sizes[i] <- current.size$total
if(move == i) {
if(current.size$total > 0) new.partition <- sample_new_partition(i, current.partition, current.size, numgroups.simulated, sizes.simulated)
if(current.size$total == 0) new.partition <- current.partition
new.partitions[,rand.o] <- rep(NA,num.nodes)
new.partitions[nodes.rand.o,rand.o] <- new.partition
}
}
}
# # this should be for neighborhood 7
# if(neighborhood[7] > 0){
# if(rand.o == 1) current.partition2 <- 1:num.nodes
# else current.partition2 <- current.partitions[,rand.o-1]
# current.size <- compute_size_neighborhood(7, current.partitions[,rand.o], partition2 = current.partition2, numgroups.simulated, sizes.simulated)
# current.sizes[7] <- current.size$total
# if(move == 7){
# if(current.size$total > 0) new.partition <- sample_new_partition(7, current.partitions[,rand.o], current.size, partition2= current.partition2, numgroups.simulated, sizes.simulated)
# if(current.size$total == 0) new.partition <- current.partition
# new.partitions[,rand.o] <- new.partition
# }
# }
# compute new statistics only if it changed
if(all(current.partitions == new.partitions, na.rm = TRUE)) {
new.z.contributions <- current.z.contributions
new.z <- current.z
} else {
recalculated.stats <- recalculate_statistics(new.partitions, rand.o, nodes.rand.o, nodes, effects, objects, current.z.contributions)
new.z.contributions <- recalculated.stats$new.z.contributions
new.z <- recalculated.stats$new.z
}
new.logit <- theta * new.z
# chose whether to change or not
new.sizes <- rep(0,n_neighborhood)
for(i in 1:n_neighborhood) {
if(neighborhood[i] > 0){
new.size <- compute_size_neighborhood(i, new.partition, numgroups.simulated, sizes.simulated)
new.sizes[i] <- new.size$total
}
}
# # this should be for neighborhood 7
# if(neighborhood[7] > 0){
# new.size <- compute_size_neighborhood(7, new.partition, current.partition2, numgroups.simulated, sizes.simulated)
# new.sizes[7] <- new.size$total
# }
neighborhoods.ratio <- sum(current.sizes * neighborhood) /
sum(new.sizes * neighborhood)
hastings.ratio <- (exp(sum(new.logit) - sum(current.logit))) * neighborhoods.ratio
return(list("hastings.ratio" = hastings.ratio,
"new.partitions" = new.partitions,
"new.z.contributions" = new.z.contributions,
"new.z" = new.z,
"new.logit" = new.logit,
"move" = move))
}
# Recalculate statistics for a changed partition in multiple estimation
recalculate_statistics <- function(new.partitions, rand.o, nodes.rand.o, nodes, effects, objects, current.z.contributions) {
# store new statistics
new.z.contributions <- current.z.contributions
# calculate separately each effect
for(e in 1:length(effects$names)) {
new.z.contributions[e,rand.o] <- step_recalculate(new.partitions, rand.o, nodes.rand.o, nodes, effects, objects, e)
}
new.z <- rowSums(new.z.contributions)
return(list(new.z.contributions = new.z.contributions,
new.z = new.z))
}
step_recalculate <- function(new.partitions, rand.o, nodes.rand.o, nodes, effects, objects, e){
num.nodes <- nrow(nodes)
num.obs <- ncol(new.partitions)
object.name <- effects$objects[e]
# TIE EFFECT: keep only present nodes
if(effects$names[e] == "tie"){
effects.temp <- list(names="tie",objects="net.temp")
for(ob in 1:length(objects)){
if(objects[[ob]][[1]] == object.name){
net <- objects[[ob]][[2]]
objects.temp <- list(list(name="net.temp",object=net[nodes.rand.o,nodes.rand.o]))
}
}
nodes.temp <- nodes[nodes.rand.o,]
new.z.contribution <- as.numeric(computeStatistics(new.partitions[nodes.rand.o,rand.o], nodes.temp, effects.temp, objects.temp))
# SAME_VAR EFFECT: adapt the varying covariate
} else if(effects$names[e] == "same_var"){
effects.temp <- list(names="same",objects="var.temp")
for(ob in 1:length(objects)){
if(objects[[ob]][[1]] == object.name){
atts <- objects[[ob]][[2]]
}
}
objects.temp <- list()
nodes.temp <- nodes[nodes.rand.o,]
nodes.temp$var.temp <- atts[nodes.rand.o,rand.o]
new.z.contribution <- as.numeric(computeStatistics(new.partitions[nodes.rand.o,rand.o], nodes.temp, effects.temp, objects.temp))
# TIE_VAR EFFECT: adapt the varying network
} else if(effects$names[e] == "tie_var"){
effects.temp <- list(names="tie",objects="net.temp")
for(ob in 1:length(objects)){
if(objects[[ob]][[1]] == object.name){
nets <- objects[[ob]][[2]]
net <- nets[[rand.o]]
objects.temp <- list(list(name="net.temp",object=net[nodes.rand.o,nodes.rand.o]))
}
}
nodes.temp <- nodes[nodes.rand.o,]
new.z.contribution <- as.numeric(computeStatistics(new.partitions[nodes.rand.o,rand.o], nodes.temp, effects.temp, objects.temp))
# TIE_VAR_X_DIFF EFFECT: adapt the varying network
} else if(effects$names[e] == "tie_var_X_diff"){
effects.temp <- list(names="tie_X_diff",objects="net.temp",objects2=effects$objects2[e])
for(ob in 1:length(objects)){
if(objects[[ob]][[1]] == object.name){
nets <- objects[[ob]][[2]]
net <- nets[[rand.o]]
objects.temp <- list(list(name="net.temp",object=net[nodes.rand.o,nodes.rand.o]))
}
}
nodes.temp <- nodes[nodes.rand.o,]
new.z.contribution <- as.numeric(computeStatistics(new.partitions[nodes.rand.o,rand.o], nodes.temp, effects.temp, objects.temp))
# INERTIA_1 EFFECT: need to recalculate the partition before and after
} else if(effects$names[e] == "inertia_1"){
aff.temp <- as.matrix(table(data.frame(actor = 1:num.nodes, group= new.partitions[,rand.o])))
adj.temp <- aff.temp %*% t(aff.temp)
diag(adj.temp) <- 0
if(rand.o > 1) { # removed && length(groups[[rand.o]]) > 0
aff.temp2 <- as.matrix(table(data.frame(actor = 1:num.nodes, group= new.partitions[,rand.o-1])))
adj.temp2 <- aff.temp2 %*% t(aff.temp2)
diag(adj.temp2) <- 0
adj12 <- adj.temp * adj.temp2
new.z.contribution <- sum(1/2 * adj12)
} else {
new.z.contribution <- 0
}
if(rand.o < num.obs) { # removed && length(groups[[rand.o]]) > 0
aff.temp2 <- as.matrix(table(data.frame(actor = 1:num.nodes, group= new.partitions[,rand.o+1])))
adj.temp2 <- aff.temp2 %*% t(aff.temp2)
diag(adj.temp2) <- 0
adj12 <- adj.temp * adj.temp2
new.z.contribution <- sum(1/2 * adj12)
}
# ANY OTHER EFFECT: the effect also exists in the single partition version
} else {
effects.temp <- list(names=effects$names[e],objects=effects$objects[e])
objects.temp <- objects
nodes.temp <- nodes[nodes.rand.o,]
new.z.contribution <- as.numeric(computeStatistics(new.partitions[nodes.rand.o,rand.o], nodes.temp, effects.temp, objects.temp))
}
return(new.z.contribution)
}
## --- FUNCTIONS FOR THE PROPOSALS MADE AT EACH STEP ----
## GENERIC FUNCTION TO COMPUTE NEIGHBORHOOD SIZE
compute_size_neighborhood <- function(i,
partition,
numgroups.simulated = NULL,
sizes.simulated = NULL){#,partition2 = NULL, # for neighborhood 7, now deprecated){
# check whether there are constraints
constraints <- FALSE
if(!is.null(numgroups.simulated) || !is.null(sizes.simulated)) {
constraints <- TRUE
}
if(i == 1) {
if(!constraints) return(compute_size_neighborhood_p1(partition))
else return(compute_size_neighborhood_p1_restricted(partition, numgroups.simulated, sizes.simulated))
}
if(i == 2) {
if(!constraints) return(compute_size_neighborhood_p2(partition))
else return(compute_size_neighborhood_p2_restricted(partition, numgroups.simulated, sizes.simulated))
}
if(i == 3) {
if(!constraints) return(compute_size_neighborhood_p3(partition))
else return(compute_size_neighborhood_p3_restricted(partition, numgroups.simulated, sizes.simulated))
}
# For now, other neighborhoods removed
# if(i == 4) {
# if(is.null(sizes.simulated)) return(compute_size_neighborhood_p4(partition))
# else return(compute_size_neighborhood_p4_restricted(partition, numgroups.simulated, sizes.simulated))
# }
# if(i == 5) {
# if(is.null(sizes.simulated)) return(compute_size_neighborhood_p5(partition))
# else return(compute_size_neighborhood_p5_restricted(partition, numgroups.simulated, sizes.simulated))
# }
# if(i == 6) {
# if(is.null(sizes.simulated)) return(compute_size_neighborhood_p6(partition))
# else return(compute_size_neighborhood_p6_restricted(partition, numgroups.simulated, sizes.simulated))
# }
# if(i == 7) {
# if(is.null(sizes.simulated)) return(compute_size_neighborhood_p7(partition))
# else return(compute_size_neighborhood_p7_restricted(partition, numgroups.simulated, sizes.simulated))
# }
}
## GENERIC FUNCTION TO COMPUTE NEIGHBORHOOD SIZE
sample_new_partition <- function(i,
current.partition,
size_neighborhood,
numgroups.simulated = NULL,
sizes.simulated = NULL){#,current.partition2 = NULL, # for neighborhood 7, now deprecated){
# check whether there are constraints
constraints <- FALSE
if(!is.null(numgroups.simulated) || !is.null(sizes.simulated)) {
constraints <- TRUE
}
if(i == 1) {
if(!constraints) return(sample_new_partition_p1(current.partition, size_neighborhood))
else return(sample_new_partition_p1_restricted(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated))
}
if(i == 2) {
if(!constraints) return(sample_new_partition_p2(current.partition, size_neighborhood))
else return(sample_new_partition_p2_restricted(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated))
}
if(i == 3) {
if(!constraints) return(sample_new_partition_p3(current.partition, size_neighborhood))
else return(sample_new_partition_p3_restricted(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated))
}
# for now, other neighborhoods removed
# if(i == 4) {
# if(is.null(sizes.simulated)) return(sample_new_partition_p4(current.partition, size_neighborhood))
# else return(sample_new_partition_p4_restricted(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated))
# }
# if(i == 5) {
# if(is.null(sizes.simulated)) return(sample_new_partition_p5(current.partition, size_neighborhood))
# else return(sample_new_partition_p5_restricted(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated))
# }
# if(i == 6) {
# if(is.null(sizes.simulated)) return(sample_new_partition_p6(current.partition, size_neighborhood))
# else return(sample_new_partition_p6_restricted(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated))
# }
# if(i == 7) {
# if(is.null(sizes.simulated)) return(sample_new_partition_p7(current.partition, size_neighborhood, current.partition2))
# else return(sample_new_partition_p7_restricted(current.partition, size_neighborhood, current.partition2, numgroups.simulated, sizes.simulated))
# }
}
## GENERIC FUNCTION TO TEST WHETHER A PARTITION IS REACHABLE WITH A GIVEN NEIGHBORHOOD
reachable <- function(i,partition1,partition2){
if(i == 1) {
return(reachable_p1(partition1,partition2))
} else if(i == 2) {
return(reachable_p2(partition1,partition2))
} else if(i == 3) {
return(reachable_p3(partition1,partition2))
}
# for now, other nieghborhoods, removed
# else if(i == 4) {
# return(reachable_p4(partition1,partition2))
# } else if(i == 5) {
# return(reachable_p5(partition1,partition2))
# } else if(i == 6) {
# return(reachable_p6(partition1,partition2))
# }
}
## --- SPECIFIC PROPOSALS ----
## NEIGHBORHOOD PI 1: only swaps of two nodes (careful, cannot be used alone)
compute_size_neighborhood_p1 <- function(partition){
# find isolates, pairs and groups>2
num.groups <- max(partition)
num.nodes <- length(partition)
sizes <- table(partition)
isolates <- as.vector(which(sizes == 1))
pairs <- as.vector(which(sizes == 2))
others <- as.vector(which(sizes > 2))
nums.swaps <- matrix(0,num.nodes,num.nodes)
for(i in 1:(num.nodes-1)){
for(j in (i+1):num.nodes){
g_i <- partition[i]
g_j <- partition[j]
allowed <- TRUE
if(g_i == g_j) allowed <- FALSE # two members of the same group cannot swap, it wouldn't change anything
if(g_i %in% isolates && g_j %in% isolates){
allowed <- FALSE # two isolates cannot swap, it wouldn't change anything
}
if(allowed) nums.swaps[i,j] <- 1
}
}
return(list(isolates = isolates,
pairs = pairs,
others = others,
nums.swaps = nums.swaps,
num.swaps = sum(nums.swaps),
total = sum(nums.swaps)))
}
# sample partition neighborhood 1
sample_new_partition_p1 <- function(current.partition, size_neighborhood){
# calculate the number of neighbor partitions
num.nodes <- length(current.partition)
num.groups <- max(current.partition)
size <- size_neighborhood
isolates <- size$isolates
pairs <- size$pairs
others <- size$others
nums.swaps <- size$nums.swaps
num.swaps <- size$num.swaps
total <- size$total
if(total == 0) return(current.partition)
pick.1 <- sample(total,1)
# decide which actors to swap
new.partition <- current.partition
swap <- which(nums.swaps == 1)[pick.1]
j <- ((swap-1) %/% num.nodes) + 1
i <- (swap-1) %% num.nodes + 1
new.partition[i] <- current.partition[j]
new.partition[j] <- current.partition[i]
new.partition <- order_groupids(new.partition)
return(new.partition)
}
# reachable neighborhood 1
reachable_p1 <- function(partition1,partition2){
# they are not neighbors if they have a different number of groups
if(max(partition1) != max(partition2)) return(FALSE)
num.nodes <- length(partition1)
check <- FALSE
i <- 1
j <- i+1
# try remove all pairs of nodes
while(!check && i <= num.nodes){
newp1 <- order_groupids( partition1[-c(i,j)] )
newp2 <- order_groupids( partition2[-c(i,j)] )
# check if they swapped
if(all(newp1 == newp2)) {
othersi1 <- which(partition1 == partition1[i])
othersi1 <- othersi1[othersi1 != i]
othersj1 <- which(partition1 == partition1[j])
othersj1 <- othersj1[othersj1 != j]
othersi2 <- which(partition2 == partition2[i])
othersi2 <- othersi2[othersi2 != i]
othersj2 <- which(partition2 == partition2[j])
othersj2 <- othersj2[othersj2 != j]
check <- (setequal(othersi1,othersj2) && setequal(othersi2,othersj1))
if(check) break
}
j <- j+1
if(j > num.nodes){
i <- i+1
j <- i+1
}
}
return(check)
}
## NEIGHBORHOOD PI 2: only merges and divisions of 2 groups
compute_size_neighborhood_p2 <- function(partition){
# calculate the number of neighbor partitions
num.nodes <- length(partition)
num.groups <- max(partition)
num.merges <- num.groups*(num.groups-1)/2
nums.divisions <- rep(0,num.groups)
num.divisions <- 0
for(k in 1:num.groups){
sk <- length(which(partition == k))
nums.divisions[k] <- 2^(sk-1) - 1
num.divisions <- num.divisions + nums.divisions[k]
}
return(list(num.merges = num.merges,
num.divisions = num.divisions,
nums.divisions = nums.divisions,
total = num.merges + num.divisions))
}
# sample partition neighborhood 2
sample_new_partition_p2 <- function(current.partition, size_neighborhood){
# calculate the number of neighbor partitions
num.nodes <- length(current.partition)
num.groups <- max(current.partition)
size <- size_neighborhood
num.merges <- size$num.merges
num.divisions <- size$num.divisions
nums.divisions <- size$nums.divisions
total <- size$total
if(total == 0) return(current.partition)
# decide between merge or division (or self loop)
pick.1 <- sample(total,1)
new.partition <- current.partition
# if merge
if(pick.1 <= num.merges) {
# pick 2 groups to merge
old_gs <- sample(num.groups,2,replace = FALSE)
# reassign one of the groups and remove useless ids
new.partition[which(new.partition == old_gs[2])] <- old_gs[1]
new.partition <- order_groupids(new.partition)
}else if(pick.1 <= (num.merges+num.divisions)){
# pick group to divide
pick.2 <- sample(num.divisions,1)
sums <- unlist(lapply(1:num.groups,FUN = function(x){sum(nums.divisions[1:x])}))
starts <- c(0,sums[1:(num.groups-1)])
ends <- sums[1:num.groups]
old_g <- which(starts < pick.2 & ends >= pick.2)
# reassign one of the groups and remove useless ids
new.groups <- sample(c(old_g,num.groups+1),length(which(new.partition==old_g)),replace = TRUE)
found <- (length(unique(new.groups)) > 1)
while(!found){
new.groups <- sample(c(old_g,num.groups+1),length(which(new.partition==old_g)),replace = TRUE)
found <- (length(unique(new.groups)) > 1)
}
new.partition[which(new.partition == old_g)] <- new.groups
new.partition <- order_groupids(new.partition)
}
return(new.partition)
}
# reachable neighborhood 2
reachable_p2 <- function(partition1,partition2){
# they are not neighbors if one does not have one more group
if(abs(max(partition1) - max(partition2))!=1) return(FALSE)
num.nodes <- length(partition1)
num.groups1 <- max(partition1)
num.groups2 <- max(partition2)
check <- FALSE
g1 <- 1
g2 <- g1+1
# try remove groups (to check divisions) and pairs of groups (to check merges)
while(!check && g1 <= num.groups1){
# first check whether g1 is not divided in partition2
membersg1 <- which(partition1 == g1)
if(length(membersg1) == num.nodes) {
newp1 <- c()
newp2 <- c()
} else {
newp1 <- order_groupids( partition1[-membersg1] )
newp2 <- order_groupids( partition2[-membersg1] )
}
if(all(newp1 == newp2)) {
check <- length(unique(partition2[membersg1])) == 2
if(check) break
}
# second check whether g1 and g2 are merged in partition2
membersg2 <- which(partition1 == g2)
if(length(c(membersg1,membersg2)) == num.nodes) {
newp1 <- c()
newp2 <- c()
} else {
newp1 <- order_groupids( partition1[-c(membersg1,membersg2)] )
newp2 <- order_groupids( partition2[-c(membersg1,membersg2)] )
}
if(all(newp1 == newp2)) {
check <- length(unique(partition2[c(membersg1,membersg2)])) == 1
if(check) break
}
g2 <- g2+1
if(g2 > num.nodes){
g1 <- g1+1
g2 <- g1+1
}
}
return(check)
}
## NEIGHBORHOOD PI 3: only swaps of one node
compute_size_neighborhood_p3 <- function(partition){
# find isolates, pairs and groups>2
num.groups <- max(partition)
num.nodes <- length(partition)
sizes <- table(partition)
isolates <- as.vector(which(sizes == 1))
pairs <- as.vector(which(sizes == 2))
others <- as.vector(which(sizes > 2))
nums.swaps <- rep(0,num.nodes)
done.pairs <- rep(0,length(pairs))
done.isolates <- rep(0,length(isolates))
for(i in 1:num.nodes){
g <- partition[i]
if(g %in% isolates){
done.isolates[isolates == g] <- 1
nums.swaps[i] <- num.groups - sum(done.isolates) # can join any other group, except other isolates before (because already done) and itself
}
if(g %in% pairs){
nums.swaps[i] <- num.groups - done.pairs[pairs == g] # can join any other group or create its isolate, except if the possibilty of splitting the pair in 2 isolates is already counted
done.pairs[pairs == g] <- 1
}
if(g %in% others){
nums.swaps[i] <- num.groups # can join any other group or create its own isolate
}
}
return(list(isolates = isolates,
pairs = pairs,
others = others,
nums.swaps = nums.swaps,
num.swaps = sum(nums.swaps),
total = sum(nums.swaps)))
}
# sample partition neighborhood 3
sample_new_partition_p3 <- function(current.partition, size_neighborhood){
# calculate the number of neighbor partitions
num.nodes <- length(current.partition)
num.groups <- max(current.partition)
size <- size_neighborhood
isolates <- size$isolates
pairs <- size$pairs
others <- size$others
nums.swaps <- size$nums.swaps
num.swaps <- size$num.swaps
total <- size$total
if(total == 0) return(current.partition)
pick.1 <- sample(total,1)
# decide which actor to swap
all.groups <- 1:num.groups
new.partition <- current.partition
done.pairs <- rep(0,length(pairs))
done.isolates <- rep(0,length(isolates))
for(i in 1:num.nodes){
if(i == 1) {start <- 0} else {start <- sum(nums.swaps[1:i-1])}
if(i == num.nodes) {end <- num.swaps} else {end <- sum(nums.swaps[1:i])}
g <- current.partition[i]
if(g %in% isolates) done.isolates[isolates == g] <- 1
if(pick.1 > start && pick.1 <= end) {
g <- current.partition[i]
# an isolate is randomly joining another group except isolates before i
if(g %in% isolates){
previousisolates <- isolates[done.isolates == 1]
tosample <- all.groups[!all.groups %in% previousisolates]
if(length(tosample) == 1) newg <- tosample
if(length(tosample) >= 2) newg <- sample(tosample,1)
new.partition[i] <- newg
new.partition <- order_groupids(new.partition)
}
# a member of a pair is randomly joining another group
# or creating its isolate if it's the first of the pair to appear
if(g %in% pairs){
if(done.pairs[pairs == g] == 0){
newg <- sample(all.groups,1)
if(newg == g){
newg <- num.groups + 1
}
} else {
tosample <- all.groups[all.groups != g]
if(length(tosample) == 1) newg <- tosample
if(length(tosample) >= 2) newg <- sample(tosample,1)
}
new.partition[i] <- newg
new.partition <- order_groupids(new.partition)
}
# a member of a bigger group is randomly joining another group
# or creating its isolate
if(g %in% others){
newg <- sample(all.groups,1)
if(newg == g){
newg <- num.groups + 1
}
new.partition[i] <- newg
new.partition <- order_groupids(new.partition)
}
}
if(g %in% pairs) done.pairs[pairs == g] <- 1
}
return(new.partition)
}
# reachable neighborhood 3
reachable_p3 <- function(partition1,partition2){
# they are not neighbors if one doesn't have the same number of groups or one more
if(abs(max(partition1) - max(partition2))>1) return(FALSE)
num.nodes <- length(partition1)
check <- FALSE
i <- 1
# try remove all nodes
while(!check && i <= num.nodes){
newp1 <- order_groupids( partition1[-i] )
newp2 <- order_groupids( partition2[-i] )
# check if it moved
if(all(newp1 == newp2)) {
othersi1 <- which(partition1[-i] == partition1[i])
othersi2 <- which(partition2[-i] == partition2[i])
check <- !setequal(othersi1,othersi2)
if(check) break
}
i <- i+1
}
return(check)
}
## --- SPECIFIC PROPOSALS RESTRICTED VERSION ----
## NEIGHBORHOOD PI 1 RESTRICTED: only swaps of two nodes (careful, cannot be used alone)
# WARNING!!!!: for now it only works if sizes allowed are an interval ([size_min,size_max])
compute_size_neighborhood_p1_restricted <- function(partition, numgroups.simulated, sizes.simulated){
# check if current partition is allowed
if(!check_sizes(partition, sizes.simulated,numgroups.simulated)) stop("The partition we are in is not allowed.")
return(compute_size_neighborhood_p1(partition))
}
# sample partition neighborhood 1 restricted
sample_new_partition_p1_restricted <- function(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated){
return(sample_new_partition_p1(current.partition, size_neighborhood))
}
## NEIGHBORHOOD PI 2 RESTRICTED: only merges and divisions of 2 groups
compute_size_neighborhood_p2_restricted <- function(partition, numgroups.simulated, sizes.simulated){
# check if current partition is allowed
if(!check_sizes(partition, sizes.simulated,numgroups.simulated)) stop("The partition we are in is not allowed.")
# find minimum and maximum size allowed
smax <- max(sizes.simulated)
smin <- min(sizes.simulated)
# calculate the number of neighbor partitions
num.nodes <- length(partition)
num.groups <- max(partition)
sizes <- table(partition)
num.merges <- 0
merges <- matrix(0,num.groups,num.groups)
nums.divisions <- rep(0,num.groups)
num.divisions <- 0
# merges
if(num.groups>1 && (num.groups-1) %in% numgroups.simulated){ # we count merges if we have more than one group and having one less group is allowed
for(g1 in 1:(num.groups-1)){
for(g2 in (g1+1):num.groups){
if(length(which(partition == g1)) + length(which(partition == g2)) <= smax) { # merges are only possible if the new group size is allowed
num.merges <- num.merges + 1
merges[g1,g2] <- 1
}
}
}
}
# divisions
if((num.groups+1) %in% numgroups.simulated){ # we count divisions if having one more group is allowed
for(k in 1:num.groups){
sg <- sizes[k]
if(smin>1){
extras <- sum(unlist(lapply(1:(smin-1),FUN=function(x){choose(sg,x)})))
if(sg%%2 == 0 && sg/2 < smin) extras <- extras - choose(sg,sg/2)/2
} else {
extras <- 0
}
nums.divisions[k] <- 2^(sg-1) - 1 - extras
num.divisions <- num.divisions + nums.divisions[k]
}
}
return(list(num.merges = num.merges,
merges = merges,
num.divisions = num.divisions,
nums.divisions = nums.divisions,
total = num.merges + num.divisions))
}
# sample partition neighborhood 2 restricted
sample_new_partition_p2_restricted <- function(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated){
# calculate the number of neighbor partitions
num.nodes <- length(current.partition)
num.groups <- max(current.partition)
size <- size_neighborhood
num.merges <- size$num.merges
merges <- size$merges
num.divisions <- size$num.divisions
nums.divisions <- size$nums.divisions
total <- size$total
if(total == 0) return(current.partition)
# sizes allowed
smax <- max(sizes.simulated)
smin <- min(sizes.simulated)
# decide between merge or division (or self loop)
pick.1 <- sample(total,1)
all.groups <- 1:num.groups
new.partition <- current.partition
# if merge
if(pick.1 <= num.merges) {
# pick 2 groups to merge
allindexes <- which(merges == 1)
p <- allindexes[pick.1]
indexes <- arrayInd(p, dim(merges))
old_g1 <- indexes[,1]
old_g2 <- indexes[,2]
# reassign one of the groups and remove useless ids
new.partition[which(new.partition == old_g2)] <- old_g1
new.partition <- order_groupids(new.partition)
}
# if division
else if(pick.1 <= (num.merges+num.divisions)){
# pick group to divide
pick.2 <- sample(num.divisions,1)
old_g1 <- 1
cpt2 <- 1
g <- 1
found <- FALSE
while(g<=num.groups && !found){
if(pick.2 >= cpt2 && pick.2 <= sum(nums.divisions[1:g])){
old_g1 <- g
found <- TRUE
}
cpt2 <- cpt2 + nums.divisions[g]
g <- g+1
}
old_g2 <- 0
new_g1 <- old_g1
new_g2 <- num.groups + 1
# reassign one of the groups and remove useless ids
new.groups <- sample(c(old_g1,num.groups+1),length(which(new.partition==old_g1)),replace = TRUE)
found <- (length(unique(new.groups)) > 1) && min(table(new.groups)) >= smin
while(!found){
new.groups <- sample(c(old_g1,num.groups+1),length(which(new.partition==old_g1)),replace = TRUE)
found <- (length(unique(new.groups)) > 1) && min(table(new.groups)) >= smin
}
new.partition[which(new.partition == old_g1)] <- new.groups
new.partition <- order_groupids(new.partition)
}
return(new.partition)
}
## NEIGHBORHOOD PI 3 RESTRICTED: only swaps one node
# WARNING!!!!: for now it only works if sizes allowed are an interval ([size_min,size_max])
compute_size_neighborhood_p3_restricted <- function(partition, numgroups.simulated, sizes.simulated){
# check if current partition is allowed
if(!check_sizes(partition, sizes.simulated, numgroups.simulated)) stop("The partition we are in is not allowed.")
# find maximum size allowed
smax <- max(sizes.simulated)
smin <- min(sizes.simulated)
sizes <- table(partition)
# find isolates, pairs and groups>2
num.groups <- max(partition)
num.nodes <- length(partition)
isolates <- as.vector(which(sizes == 1))
pairs <- as.vector(which(sizes == 2))
others <- as.vector(which(sizes > 2))
abovemin_groups <- as.vector(which(sizes > smin))
belowmax_groups <- as.vector(which(sizes < smax))
nums.swaps <- rep(0,num.nodes)
done.pairs <- rep(0,length(pairs))
done.isolates <- rep(0,length(isolates))
for(i in 1:num.nodes){
g <- partition[i]
if(g %in% isolates && (num.groups-1) %in% numgroups.simulated){
# i can join any other group (<smax), except other isolates done before (because already done) and this one (that would not change anything)
# since it's an isolate, it means that isolates are allowed in that context so we don't worry about smin (=1)
done.isolates[isolates == g] <- 1
nums.swaps[i] <- length(belowmax_groups) - sum(done.isolates)
} # to move an isolate, a partition with one group less must be allowed
if(g %in% pairs && g %in% abovemin_groups){
# if isolates are allowed, then i can join any other group (smax), and splitting if it's not already counted
# here we don't worry about smin (=1) because 2 > smin
gnotbelowmax <- !(g %in% belowmax_groups)
if((num.groups+1) %in% numgroups.simulated && done.pairs[pairs == g] == 0) { # we check whether it's possible to create an isolate
nums.swaps[i] <- length(belowmax_groups) + gnotbelowmax
} else {
nums.swaps[i] <- length(belowmax_groups) + gnotbelowmax - 1
}
done.pairs[pairs == g] <- 1
} # to break a pair, size 2 should not be the minimum allowed
if(g %in% others && g %in% abovemin_groups){
# can join any other group (<smax) or create its own isolate (if it's allowed)
# we don't care here about checking that other groups are above min size, if they exist it means they are already allowed
gnotbelowmax <- !(g %in% belowmax_groups)
if(smin == 1 && (num.groups+1) %in% numgroups.simulated) { # we check whether it's possible to create an isolate
nums.swaps[i] <- length(belowmax_groups) + gnotbelowmax
} else {
nums.swaps[i] <- length(belowmax_groups) + gnotbelowmax - 1
}
} # to break a group, it should not be of minimal size
}
return(list(isolates = isolates,
pairs = pairs,
others = others,
abovemin_groups = abovemin_groups,
belowmax_groups = belowmax_groups,
nums.swaps = nums.swaps,
num.swaps = sum(nums.swaps),
total = sum(nums.swaps)))
}
# sample partition neighborhood 3 restricted
sample_new_partition_p3_restricted <- function(current.partition, size_neighborhood, numgroups.simulated, sizes.simulated){
# calculate the number of neighbor partitions
num.nodes <- length(current.partition)
num.groups <- max(current.partition)
size <- size_neighborhood
isolates <- size$isolates
pairs <- size$pairs
others <- size$others
abovemin_groups <- size$abovemin_groups
belowmax_groups <- size$belowmax_groups
nums.swaps <- size$nums.swaps
num.swaps <- size$num.swaps
total <- size$total
if(total == 0) return(current.partition)
# allowed sizes
smax <- max(sizes.simulated)
smin <- min(sizes.simulated)
pick.1 <- sample(total,1)
# decide which actor to swap
all.groups <- 1:num.groups
new.partition <- current.partition
done.pairs <- rep(0,length(pairs))
done.isolates <- rep(0,length(isolates))
for(i in 1:num.nodes){
if(i == 1) {start <- 0} else {start <- sum(nums.swaps[1:i-1])}
if(i == num.nodes) {end <- num.swaps} else {end <- sum(nums.swaps[1:i])}
g <- current.partition[i]
if(g %in% isolates) done.isolates[isolates == g] <- 1
if(pick.1 > start && pick.1 <= end) {
# an isolate can join any other group (<smax), except other isolates done before (because already done) and this one (that would not change anything)
# if we are here, this implies that a partition with one less group is allowed
if(g %in% isolates){
previousisolates <- isolates[done.isolates == 1]
tosample <- all.groups[(all.groups %in% belowmax_groups) & !(all.groups %in% previousisolates)]
if(length(tosample) == 1) newg <- tosample
if(length(tosample) >= 2) newg <- sample(tosample,1)
new.partition[i] <- newg
new.partition <- order_groupids(new.partition)
}
# a member of a pair can join any other group (<smax), and splitting if it's not already counted, if isolates are allowed
if(g %in% pairs && g %in% abovemin_groups){
if((num.groups+1) %in% numgroups.simulated && done.pairs[pairs == g] == 0){ # if a split of the pair is possible
tosample <- all.groups[(all.groups %in% belowmax_groups) | (all.groups == g)]
if(length(tosample) == 1) newg <- tosample
if(length(tosample) >= 2) newg <- sample(tosample,1)
if(newg == g){
newg <- num.groups + 1
}
} else {
tosample <- all.groups[(all.groups %in% belowmax_groups) & (all.groups != g)]
if(length(tosample) == 1) newg <- tosample
if(length(tosample) >= 2) newg <- sample(tosample,1)
}
new.partition[i] <- newg
new.partition <- order_groupids(new.partition)
}
# a member of a bigger group (>smin) can join any other group (<smax) or create its own isolate (if it's allowed)
if(g %in% others && g %in% abovemin_groups){
if(smin == 1 && (num.groups+1) %in% numgroups.simulated) { # if we can create a isolate
tosample <- all.groups[(all.groups %in% belowmax_groups) | (all.groups == g)]
if(length(tosample) == 1) newg <- tosample
if(length(tosample) >= 2) newg <- sample(tosample,1)
if(newg == g){
newg <- num.groups + 1
}
} else {
tosample <- all.groups[(all.groups %in% belowmax_groups) & (all.groups != g)]
if(length(tosample) == 1) newg <- tosample
if(length(tosample) >= 2) newg <- sample(tosample,1)
}
new.partition[i] <- newg
new.partition <- order_groupids(new.partition)
}
}
if(g %in% pairs) done.pairs[pairs == g] <- 1
}
return(new.partition)
}
### Deprecated: unused neighborhoods
### Can be reintroduced, but they don't work with constraints on number of groups (only on sizes)
# ## NEIGHBORHOOD PI 4: swaps a pair of nodes
#
# compute_size_neighborhood_p4 <- function(partition){
#
# # find isolates, pairs and groups>2
# num.groups <- max(partition)
# num.nodes <- length(partition)
# sizes <- table(partition)
# isolates <- as.vector(which(sizes == 1))
# pairs <- as.vector(which(sizes == 2))
# others <- as.vector(which(sizes > 2))
#
# nums.swaps <- rep(0,max(partition))
#
# for(g in 1:num.groups){
#
# #if(g %in% isolates){}# nothing happens, there is no pair to swap here
#
# if(g %in% pairs){
# nums.swaps[g] <- length(pairs[pairs>g]) + length(isolates) + length(others)
# }
# if(g %in% others){
# sizeg <- sum(partition == g)
# nums.swaps[g] <- (length(pairs) + length(others) + length(isolates) ) * sizeg * (sizeg-1) / 2
# }
# }
#
# return(list(isolates = isolates,
# pairs = pairs,
# others = others,
# nums.swaps = nums.swaps,
# num.swaps = sum(nums.swaps),
# total = sum(nums.swaps)))
#
# }
#
#
# sample_new_partition_p4 <- function(current.partition, size_neighborhood){
#
# # calculate the number of neighbor partitions
# num.nodes <- length(current.partition)
# num.groups <- max(current.partition)
# size <- size_neighborhood
# isolates <- size$isolates
# pairs <- size$pairs
# others <- size$others
# nums.swaps <- size$nums.swaps
# num.swaps <- size$num.swaps
# total <- size$total
#
# if(total == 0) return(current.partition)
#
# pick.1 <- sample(total,1)
#
# # decide which group to pick a pair from
# all.groups <- 1:num.groups
# new.partition <- current.partition
# #done.pairs <- rep(0,length(pairs))
#
# sums <- unlist(lapply(1:num.groups,FUN = function(x){sum(nums.swaps[1:x])}))
# starts <- c(0,sums[1:(num.groups-1)])
# ends <- sums[1:num.groups]
# old_g <- which(starts < pick.1 & ends >= pick.1)
#
#
# # a pair can join other available groups, unless it's a pair that has already been "joined"
# if(old_g %in% pairs){
# members <- which(current.partition == old_g)
# tosample <- all.groups[!(all.groups %in% pairs & all.groups <= old_g)]
# if(length(tosample) == 1) newg <- tosample else newg <- sample(tosample,1)
# new.partition[members] <- newg
# new.partition <- order_groupids(new.partition)
# }
#
# # a pair in a bigger group can join any other group or create its own isolated pair
# if(old_g %in% others){
# members <- which(current.partition == old_g)
# sizeg <- length(members)
# pairg <- sample(members,2)
# tosample <- all.groups
# if(length(tosample) == 1) newg <- tosample else newg <- sample(tosample,1)
# if(newg == old_g) newg <- num.groups + 1
# new.partition[pairg] <- newg
# new.partition <- order_groupids(new.partition)
# }
#
# return(new.partition)
#
# }
#
#
# reachable_p4 <- function(partition1,partition2){
#
# # they are not neighbors if one doesn't have the same number of groups or one more
# if(abs(max(partition1) - max(partition2))>1) return(FALSE)
#
# num.nodes <- length(partition1)
# check <- FALSE
# i <- 1
# j <- i+1
#
# # try remove all pairs of nodes
# while(!check && i <= num.nodes){
#
# if(partition1[i] == partition1[j] && partition2[i] == partition2[j]) {
# newp1 <- order_groupids( partition1[-c(i,j)] )
# newp2 <- order_groupids( partition2[-c(i,j)] )
#
# # check if the pair moved together
# if(all(newp1 == newp2)) {
# others1 <- which(partition1[-c(i,j)] == partition1[i])
# others2 <- which(partition2[-c(i,j)] == partition2[i])
# check <- !setequal(others1,others2)
# if(check) break
# }
# }
#
# j <- j+1
# if(j > num.nodes){
# i <- i+1
# j <- i+1
# }
# }
#
# return(check)
# }
#
#
# ## NEIGHBORHOOD PI 5: swap two pairs of nodes
#
# compute_size_neighborhood_p5 <- function(partition){
#
# # find isolates, pairs and groups>2
# num.groups <- max(partition)
# num.nodes <- length(partition)
# sizes <- table(partition)
# isolates <- as.vector(which(sizes == 1))
# pairs <- as.vector(which(sizes == 2))
# others <- as.vector(which(sizes > 2))
#
# aff <- as.matrix(table(data.frame(actor = 1:num.nodes, group = partition)))
# net <- aff %*% t(aff)
# net[upper.tri(net)] <- 0
# diag(net) <- 0
# paired.actors <- which(net == 1)
# num.paired.actors <- sum(net)
#
# groups.paired.actors <- rep(0,num.paired.actors)
# pairs.paired.actors <- rep(0,num.paired.actors)
#
# impossible.swaps <- 0
# for(g in 1:num.groups){
# members <- which(partition == g)
#
# if(g %in% pairs){
# i <- members[1]
# j <- members[2]
# p <- (i-1)*num.nodes+j
# pairs.paired.actors[which(paired.actors == p)] <- 1
# groups.paired.actors[which(paired.actors == p)] <- g
# }
#
# if(g %in% others){
# for(i in 1:(length(members)-1)){
# for(j in (i+1):length(members)){
# p <- (members[i]-1)*num.nodes+members[j]
# groups.paired.actors[which(paired.actors == p)] <- g
# }
# }
# num.pairs.g <- length(members)*(length(members)-1)/2
# impossible.swaps <- impossible.swaps + num.pairs.g*(num.pairs.g-1)/2
# }
#
# }
#
# num.swaps <- num.paired.actors * (num.paired.actors-1) / 2 - impossible.swaps - sum(pairs.paired.actors) * (sum(pairs.paired.actors)-1) / 2
#
# return(list(isolates = isolates,
# pairs = pairs,
# others = others,
# paired.actors = paired.actors,
# groups.paired.actors = groups.paired.actors,
# pairs.paired.actors = pairs.paired.actors,
# num.swaps = num.swaps,
# total = num.swaps))
#
# }
#
#
# sample_new_partition_p5 <- function(current.partition, size_neighborhood){
#
# # calculate the number of neighbor partitions
# num.nodes <- length(current.partition)
# num.groups <- max(current.partition)
# size <- size_neighborhood
# isolates <- size$isolates
# pairs <- size$pairs
# others <- size$others
# paired.actors <- size$paired.actors
# groups.paired.actors <- size$groups.paired.actors
# pairs.paired.actors <- size$pairs.paired.actors
# num.swaps <- size$num.swaps
# total <- size$total
#
# if(total == 0) return(current.partition)
#
# new.partition <- current.partition
#
# # decide which actors to swap
# if(num.swaps > 0){
# found <- FALSE
# while(!found){
# pick.1 <- sample(paired.actors,1)
# pick.2 <- sample(paired.actors,1)
# if(groups.paired.actors[which(paired.actors == pick.1)] != groups.paired.actors[which(paired.actors == pick.2)] &&
# !(pairs.paired.actors[which(paired.actors == pick.1)] && pairs.paired.actors[which(paired.actors == pick.2)])) {
# found <- TRUE
# }
# }
#
# j1 <- ((pick.1-1) %/% num.nodes) + 1
# i1 <- (pick.1-1) %% num.nodes + 1
# j2 <- ((pick.2-1) %/% num.nodes) + 1
# i2 <- (pick.2-1) %% num.nodes + 1
# new.partition[i1] <- current.partition[i2]
# new.partition[j1] <- current.partition[j2]
# new.partition[i2] <- current.partition[i1]
# new.partition[j2] <- current.partition[j1]
# new.partition <- order_groupids(new.partition)
# }
#
#
# return(new.partition)
#
# }
#
#
# reachable_p5 <- function(partition1,partition2){
#
# # they are not neighbors if they have a different number of groups
# if(max(partition1) != max(partition2)) return(FALSE)
#
# num.nodes <- length(partition1)
# check <- FALSE
# i <- 1
# j <- i+1
# k <- 1
# l <- k+1
#
# # try remove all pairs of pairs
# while(!check && i <= num.nodes){
#
# if(partition1[i] == partition1[j] && partition2[i] == partition2[j] &&
# partition1[k] == partition1[l] && partition2[k] == partition2[l]) {
# newp1 <- order_groupids( partition1[-c(i,j,k,l)] )
# newp2 <- order_groupids( partition2[-c(i,j,k,l)] )
#
# # check if the pairs were swapped
# if(all(newp1 == newp2)) {
# othersij1 <- which(partition1[-c(i,j)] == partition1[i])
# othersij2 <- which(partition2[-c(i,j)] == partition2[i])
# otherskl1 <- which(partition1[-c(k,l)] == partition1[k])
# otherskl2 <- which(partition2[-c(k,l)] == partition2[k])
# check <- (setequal(othersij1,otherskl2) && setequal(othersij2,otherskl1))
# if(check) break
# }
# }
#
# l <- l+1
# if(l > num.nodes){
# k <- k+1
# l <- k+1
# if(k > num.nodes){
# j <- j+1
# if(j > num.nodes){
# i <- i+1
# j <- i+1
# }
# }
# }
# }
#
# return(check)
# }
#
#
# ## NEIGHBORHOOD PI 6: pick two groups, re attribute nodes within them (careful, cannot be used alone)
#
# compute_size_neighborhood_p6 <- function(partition){
#
# # find isolates, pairs and groups>2
# num.groups <- max(partition)
# num.nodes <- length(partition)
# sizes <- table(partition)
# isolates <- as.vector(which(sizes == 1))
# pairs <- as.vector(which(sizes == 2))
# others <- as.vector(which(sizes > 2))
#
# group.pairs <- matrix(0,num.groups,num.groups)
# nums.swaps <- matrix(0,num.groups,num.groups)
#
# if(num.groups > 1) {
#
# for(g1 in 1:(num.groups-1)){
# for(g2 in (g1+1):num.groups){
# size.1 <- sum(partition == g1)
# size.2 <- sum(partition == g2)
#
# if(size.1 == 1 && size.2 == 1){
# nums.swaps[g1,g2] <- 0
# } else {
# k <- size.1
# n <- size.1 + size.2
# nums.swaps[g1,g2] <- factorial(n) / (factorial(k)*factorial(n-k))
# if(size.1 == size.2) nums.swaps[g1,g2] <- nums.swaps[g1,g2] /2
# nums.swaps[g1,g2] <- nums.swaps[g1,g2] -1
# }
# }
# }
# }
#
# return(list(isolates = isolates,
# pairs = pairs,
# others = others,
# nums.swaps = nums.swaps,
# num.swaps = sum(nums.swaps),
# total = sum(nums.swaps)))
#
# }
#
#
# sample_new_partition_p6 <- function(current.partition, size_neighborhood){
#
# # calculate the number of neighbor partitions
# num.nodes <- length(current.partition)
# num.groups <- max(current.partition)
# size <- size_neighborhood
# isolates <- size$isolates
# pairs <- size$pairs
# others <- size$others
# nums.swaps <- size$nums.swaps
# num.swaps <- size$num.swaps
# total <- size$total
#
# if(total == 0) return(current.partition)
#
# # decide which actors to swap
# new.partition <- current.partition
# pick <- sample(1:num.swaps,1)
# start <- 0
# end <- 0
#
# if(num.groups > 1){
# for(g1 in 1:(num.groups-1)){
# for(g2 in (g1+1):num.groups){
#
# end <- end + nums.swaps[g1,g2]
#
# if(start < pick && pick <= end){
# found <- FALSE
# while(!found){
# tosample <- c(which(current.partition==g1),which(current.partition==g2))
# new.g1 <- sample(tosample,sum(current.partition==g1))
# new.g2 <- tosample[-which(tosample %in% new.g1)]
# if( !(setequal(new.g1,current.partition[g1]) && setequal(new.g2,current.partition[g2])) &&
# !(setequal(new.g1,current.partition[g2]) && setequal(new.g2,current.partition[g1])) ) found <- TRUE
# }
# new.partition[new.g1] <- g1
# new.partition[new.g2] <- g2
# new.partition <- order_groupids(new.partition)
# }
#
# start <- start + nums.swaps[g1,g2]
# }
# }
# }
#
# return(new.partition)
#
# }
#
#
# reachable_p6 <- function(partition1,partition2){
#
# # they are not neighbors if they have a different number of groups
# if(max(partition1) != max(partition2)) return(FALSE)
#
# num.nodes <- length(partition1)
# num.groups1 <- max(partition1)
# num.groups2 <- max(partition2)
# check <- FALSE
# g1 <- 1
# g2 <- g1+1
#
# # try remove groups (to check divisions) and pairs of groups (to check merges)
# while(!check && g <= num.groups1){
#
# # first check whether g1 is not divided in partition2
# membersg1 <- which(partition1 == g1)
# newp1 <- order_groupids( partition1[-membersg1] )
# newp2 <- order_groupids( partition2[-membersg1] )
#
# if(all(newp1 == newp2)) {
# check <- length(unique(partition2[membersg1])) == 2
# if(check) break
# }
#
# # second check whether g1 and g2 are merged in partition2
# membersg2 <- which(partition1 == g2)
# newp1 <- order_groupids( partition1[-c(membersg1,membersg2)] )
# newp2 <- order_groupids( partition2[-c(membersg1,membersg2)] )
#
# if(all(newp1 == newp2)) {
# check <- length(unique(partition2[c(membersg1,membersg2)])) == 1
# if(check) break
# }
#
# g2 <- g2+1
# if(g2 > num.nodes){
# g1 <- g1+1
# g2 <- g2+1
# }
# }
#
# return(check)
# }
#
# ## NEIGHBORHOOD PI 7: take two nodes together in partition 2, bring them together (single move)
# ## otherwise separate then
# ## careful need isolates
#
# compute_size_neighborhood_p7 <- function(partition, partition2){
#
# num.nodes <- length(partition)
# present <- !is.na(partition)
# present2 <- !is.na(partition2)
#
# num.groups <- max(partition, na.rm = TRUE)
# sizes <- table(partition)
# isolates <- as.vector(which(sizes == 1))
# pairs <- as.vector(which(sizes == 2))
#
# affiliation <- as.matrix(table(data.frame(actor = 1:num.nodes, group= partition)))
# affiliation2 <- as.matrix(table(data.frame(actor = 1:num.nodes, group= partition2)))
# adjacency <- affiliation %*% t(affiliation)
# adjacency2 <- affiliation2 %*% t(affiliation2)
# diag(adjacency) <- 0
# diag(adjacency2) <- 0
#
# nums.swaps <- matrix(0,num.nodes,num.nodes)
# for(i in 1:num.nodes) {
# for(j in 1:num.nodes){
# ad <- adjacency[i,j]
# ad2 <- adjacency2[i,j]
# if(!is.na(ad) && !is.na(ad2) && ad2 == 1) {
# if(ad == 0) { # if separated join if possible (two isolates will only join once)
# if(partition[i] %in% isolates && partition[j] %in% isolates && j < i) {
# nums.swaps[i,j] <- 0
# } else {
# nums.swaps[i,j] <- 1
# }
# }
# if(ad == 1) { # if together separate if possible (a pair will only separate once)
# if(partition[i] %in% pairs && j < i) {
# nums.swaps[i,j] <- num.groups - 1
# } else {
# nums.swaps[i,j] <- num.groups
# }
#
# }
# }
# }
# }
#
# return(list(adjacency = adjacency,
# pairs = pairs,
# nums.swaps = nums.swaps,
# num.swaps = sum(nums.swaps),
# total = sum(nums.swaps)))
# }
#
#
# sample_new_partition_p7 <- function(current.partition, current.partition2, size_neighborhood){
#
# num.nodes <- length(current.partition)
# num.groups <- max(current.partition,na.rm = TRUE)
# adjacency <- size_neighborhood$adjacency
# pairs <- size_neighborhood$pairs
# nums.swaps <- size_neighborhood$nums.swaps
# num.swaps <- size_neighborhood$num.swaps
# total <- size$total
#
# if(total == 0) return(current.partition1)
#
# present <- !(is.na(current.partition))
# new.partition <- current.partition
#
# probas_pairs <- as.vector(nums.swaps) / num.swaps
# pick <- sample(1:(num.nodes*num.nodes),1,prob=probas_pairs)
# i <- (pick-1) %% num.nodes + 1
# j <- ((pick-1) %/% num.nodes) + 1
#
# if(adjacency[i,j] == 0) {
# new.partition[i] <- current.partition[j]
# } else {
# if(current.partition[i] %in% pairs && j < i) {
# tosample <- (1:num.groups)[(1:num.groups) != current.partition[i]]
# if(length(tosample) == 1) pickgroup <- tosample
# else pickgroup <- sample(tosample,1)
# if(pickgroup == current.partition[i]) new.partition[i] <- num.groups + 1
# else new.partition[i] <- pickgroup
# }else {
# tosample <- 1:num.groups
# if(length(tosample) == 1) pickgroup <- tosample
# else pickgroup <- sample(tosample,1)
# if(pickgroup == current.partition[i]) new.partition[i] <- num.groups + 1
# else new.partition[i] <- pickgroup
# }
# }
# new.partition[present] <- order_groupids(new.partition[present])
#
# return(new.partition)
#
# }
#
#
#
# ## NEIGHBORHOOD PI 4 RESTRICTED: swaps a pair of nodes
# # WARNING!!!!: for now it only works if sizes allowed are an interval ([size_min,size_max])
# compute_size_neighborhood_p4_restricted <- function(partition, sizes.simulated){
#
# # check if current partition is allowed
# if(!check_sizes(partition, sizes.simulated)) {
# print(partition)
# #stop("The partition we are in is not allowed.")
# }
#
# # find maximum size allowed
# smax <- max(sizes.simulated)
# smin <- min(sizes.simulated)
#
# # find isolates, pairs and groups>2
# num.groups <- max(partition)
# num.nodes <- length(partition)
# sizes <- table(partition)
# isolates <- as.vector(which(sizes == 1))
# pairs <- as.vector(which(sizes == 2))
# others <- as.vector(which(sizes > 2))
# abovemin_groups <- as.vector(which(sizes > (smin+1)))
# belowmax_groups <- as.vector(which(sizes < (smax-1)))
#
# nums.swaps <- rep(0,max(partition))
# done.pairs <- rep(0,length(pairs))
# for(g in 1:num.groups){
#
# if(g %in% isolates){
# # nothing happens, there is no pair to swap here
# }
# if(g %in% pairs){
# # then the pair can move to a group <smax-1, unless it's another pair and the combination of the two pairs is already counted
# # is 2 <smax-1 we have to remove the pair from the count of possible other groups
# if(2 < smax-1) nums.swaps[g] <- length(belowmax_groups) - 1 - sum(done.pairs) + done.pairs[pairs==g]
# if(2 >= smax-1) nums.swaps[g] <- length(belowmax_groups) # then it's only isolates
# done.pairs[pairs == g] <- 1
# }
# if(g %in% others && g %in% abovemin_groups){
# # then any pair inside the group can be separated and added to an available group that is <smax-1, unless the group is not >smin+1
# sizeg <- sum(partition == g)
# gbelowmax <- (g %in% belowmax_groups)
# # if it can create its own isolated pair
# if(gbelowmax && 2 %in% sizes.simulated) nums.swaps[g] <- length(belowmax_groups) * sizeg * (sizeg-1) / 2
# if(!gbelowmax && 2 %in% sizes.simulated) nums.swaps[g] <- (length(belowmax_groups) + 1) * sizeg * (sizeg-1) / 2
# # if not
# if(gbelowmax && !(2 %in% sizes.simulated)) nums.swaps[g] <- (length(belowmax_groups)-1) * sizeg * (sizeg-1) / 2
# if(!gbelowmax && !(2 %in% sizes.simulated)) nums.swaps[g] <- length(belowmax_groups) * sizeg * (sizeg-1) / 2
# } # we can never break a group that is already of minimal size
# }
#
# return(list(isolates = isolates,
# pairs = pairs,
# others = others,
# abovemin_groups = abovemin_groups,
# belowmax_groups = belowmax_groups,
# nums.swaps = nums.swaps,
# num.swaps = sum(nums.swaps),
# total = sum(nums.swaps)))
#
# }
#
# sample_new_partition_p4_restricted <- function(current.partition, size_neighborhood, sizes.simulated){
#
# # calculate the number of neighbor partitions
# num.nodes <- length(current.partition)
# num.groups <- max(current.partition)
# size <- size_neighborhood
# isolates <- size$isolates
# pairs <- size$pairs
# others <- size$others
# abovemin_groups <- size$abovemin_groups
# belowmax_groups <- size$belowmax_groups
# nums.swaps <- size$nums.swaps
# num.swaps <- size$num.swaps
# total <- size$total
#
# # allowed sizes
# smax <- max(sizes.simulated)
# smin <- min(sizes.simulated)
#
# pick.1 <- sample(total,1)
#
# # decide which group to pick a pair from
# all.groups <- 1:num.groups
# new.partition <- current.partition
# done.pairs <- rep(0,length(pairs))
#
# for(g in 1:num.groups){
#
# if(g == 1) {start <- 0} else {start <- sum(nums.swaps[1:g-1])}
# if(g == num.groups) {end <- num.swaps} else {end <- sum(nums.swaps[1:g])}
#
# if(pick.1 > start && pick.1 <= end) {
#
# # a pair can join other available groups (<smax-1), unless it's a pair that has already been "joined"
# if(g %in% pairs){
# members <- which(current.partition == g)
# previous.pairs <- pairs[done.pairs == 1]
# if(2 < smax-1){
# tosample <- all.groups[(all.groups %in% belowmax_groups) & (all.groups != g) & !(all.groups %in% previous.pairs)]
# if(length(tosample) == 1) newg <- tosample
# if(length(tosample) >= 2) newg <- sample(tosample,1)
# }else{
# tosample <- all.groups[(all.groups %in% belowmax_groups)]
# if(length(tosample) == 1) newg <- tosample
# if(length(tosample) >= 2) newg <- sample(tosample,1)
# }
# new.partition[members] <- newg
# new.partition <- order_groupids(new.partition)
# }
#
# # a pair in a bigger group (>smin+1) can join any other group (<smax-1) or create its isolate
# if(g %in% others && g %in% abovemin_groups){
# members <- which(current.partition == g)
# sizeg <- length(members)
# pairg <- sample(members,2)
# if(2 %in% sizes.simulated) tosample <- all.groups[(all.groups %in% belowmax_groups) || (all.groups == g)]
# else tosample <- all.groups[(all.groups %in% belowmax_groups) & (all.groups != g)]
# if(length(tosample) == 1) newg <- tosample
# if(length(tosample) >= 2) newg <- sample(tosample,1)
# if(newg == g) newg <- num.groups + 1
# new.partition[pairg] <- newg
# new.partition <- order_groupids(new.partition)
# }
#
# }
#
# if(g %in% pairs) done.pairs[pairs == g] <- 1
# }
#
# return(new.partition)
#
# }
#
#
# ## NEIGHBORHOOD PI 5 RESTRICTED: only swaps of two pairs of nodes (careful, cannot be used alone)
# compute_size_neighborhood_p5_restricted <- function(partition, sizes.simulated){
# # check if current partition is allowed
# if(!check_sizes(partition, sizes.simulated)) stop("The partition we are in is not allowed.")
# return(compute_size_neighborhood_p5(partition))
# }
#
# sample_new_partition_p5_restricted <- function(current.partition, size_neighborhood, sizes.simulated){
# return(sample_new_partition_p5(current.partition, size_neighborhood))
# }
#
#
# ## NEIGHBORHOOD PI 6 RESTRICTED: exchanges of two groups (careful, cannot be used alone)
# compute_size_neighborhood_p6_restricted <- function(partition, sizes.simulated){
# # check if current partition is allowed
# if(!check_sizes(partition, sizes.simulated)) stop("The partition we are in is not allowed.")
# return(compute_size_neighborhood_p6(partition))
# }
#
# sample_new_partition_p6_restricted <- function(current.partition, size_neighborhood, sizes.simulated){
# return(sample_new_partition_p6(current.partition, size_neighborhood))
# }
#
#
# ## NEIGHBORHOOD PI 7 RESTRICTED: take two nodes together in partition 2, bring them together (single move)
# ## otherwise separate then
# ## careful need isolates
# compute_size_neighborhood_p7_restricted <- function(partition, partition2, sizes.simulated){
#
# # find maximum size allowed
# smax <- max(sizes.simulated)
# smin <- min(sizes.simulated)
#
# num.nodes <- length(partition)
# present <- !is.na(partition)
# present2 <- !is.na(partition2)
#
# num.groups <- max(partition,na.rm = TRUE)
# sizes <- table(partition)
# isolates <- as.vector(which(sizes == 1))
# pairs <- as.vector(which(sizes == 2))
# abovemin_groups <- as.vector(which(sizes > smin))
# belowmax_groups <- as.vector(which(sizes < smax))
# max_groups <- as.vector(which(sizes == smax))
#
# affiliation <- as.matrix(table(data.frame(actor = 1:num.nodes, group= partition)))
# affiliation2 <- as.matrix(table(data.frame(actor = 1:num.nodes, group= partition2)))
# adjacency <- affiliation %*% t(affiliation)
# adjacency2 <- affiliation2 %*% t(affiliation2)
# diag(adjacency) <- 0
# diag(adjacency2) <- 0
#
# nums.swaps <- matrix(0,num.nodes,num.nodes)
# for(i in 1:num.nodes) {
# for(j in 1:num.nodes){
# ad <- adjacency[i,j]
# ad2 <- adjacency2[i,j]
# if(!is.na(ad) && !is.na(ad2) && ad2 == 1) {
# if(ad == 0) { # if separated join if possible (two isolates will only join once)
# if(partition[i] %in% isolates && partition[j] %in% isolates && j < i) {
# nums.swaps[i,j] <- 0
# } else {
# nums.swaps[i,j] <- (partition[i] %in% c(abovemin_groups,isolates) ) & (partition[j] %in% belowmax_groups)
# }
# }
# if(ad == 1) { # if together separate if possible (a pair will only separate once)
# if(partition[i] %in% abovemin_groups) {
# if(partition[i] %in% pairs && j < i){
# nums.swaps[i,j] <- num.groups - length(max_groups) + (partition[i] %in% max_groups) - !(1 %in% sizes.simulated) - 1
# } else {
# nums.swaps[i,j] <- num.groups - length(max_groups) + (partition[i] %in% max_groups) - !(1 %in% sizes.simulated)
# }
# }
# }
# }
# }
# }
#
# return(list(abovemin_groups = abovemin_groups,
# belowmax_groups = belowmax_groups,
# adjacency = adjacency,
# pairs = pairs,
# nums.swaps = nums.swaps,
# num.swaps = sum(nums.swaps),
# total = sum(nums.swaps)))
# }
#
# sample_new_partition_p7_restricted <- function(current.partition, size_neighborhood, current.partition2, sizes.simulated){
#
# # find maximum size allowed
# smax <- max(sizes.simulated)
# smin <- min(sizes.simulated)
#
# num.nodes <- length(current.partition)
# num.groups <- max(current.partition,na.rm = TRUE)
# abovemin_groups <- size_neighborhood$abovemin_groups
# belowmax_groups <- size_neighborhood$belowmax_groups
# adjacency <- size_neighborhood$adjacency
# pairs <- size_neighborhood$pairs
# nums.swaps <- size_neighborhood$nums.swaps
# num.swaps <- size_neighborhood$num.swaps
# total <- size$total
#
# if(total == 0) return(current.partition)
#
# present <- !(is.na(current.partition))
# new.partition <- current.partition
#
# probas_pairs <- as.vector(nums.swaps) / num.swaps
# pick <- sample(1:(num.nodes*num.nodes),1,prob=probas_pairs)
# i <- (pick-1) %% num.nodes + 1
# j <- ((pick-1) %/% num.nodes) + 1
#
# if(adjacency[i,j] == 0) {
# new.partition[i] <- current.partition[j]
# }
# if(adjacency[i,j] == 1) {
# if(current.partition[i] %in% pairs && j < i){
# tosample <- belowmax_groups[belowmax_groups != current.partition[i]]
# } else {
# if(1 %in% sizes.simulated) tosample <- unique(c(belowmax_groups,current.partition[i]))
# else tosample <- belowmax_groups[belowmax_groups != current.partition[i]]
# }
# if(length(tosample) == 1) pickgroup <- tosample
# else pickgroup <- sample(tosample,1)
# if(pickgroup == current.partition[i]) new.partition[i] <- num.groups + 1
# else new.partition[i] <- pickgroup
# }
# new.partition[present] <- order_groupids(new.partition[present])
#
# return(new.partition)
#
# }
#
#
##############################################################################
### Deprecated 2: neighborhoods for the multiple partition estimation
# ## NEIGHBORHOOD PI 7: pick two partitions and exchange then (only for multiple partition estimation)
# # since we force one change, and it's reversible (bijection), we don't have to care about the size of the neighborhood
# # CAREFUL, only works when singletons are allowed
# compute_size_neighborhood_p7 <- function(partition1, partition2){
#
# num.nodes <- length(partition1)
# present1 <- !is.na(partition1)
# present2 <- !is.na(partition2)
#
# # calculate the number of ways to reconstruct p1 after having done the exchange with p2
# p1withoutabsentp2 <- partition1[present1 & present2]
# max1 <- length(unique(p1withoutabsentp2))
# nodesunknown1 <- sum(present1 & !present2)
# if(nodesunknown1 == 0) {
# possibilities1 <- 1
# } else {
# possibilities1 <- 0
# for(k in 1:nodesunknown1) {
# possibilities1 <- possibilities1 + Stirling(nodesunknown1,k)*(max1+1)^k
# }
# }
#
# # same the other way around
# p2withoutabsentp1 <- partition2[present2 & present1]
# max2 <- length(unique(p2withoutabsentp1))
# nodesunknown2 <- sum(present2 & !present1)
# if(nodesunknown2 == 0) {
# possibilities2 <- 1
# } else {
# possibilities2 <- 0
# for(k in 1:nodesunknown2) {
# possibilities2 <- possibilities2 + Stirling(nodesunknown2,k)*(max2+1)^k
# }
# }
#
# return(possibilities1 * possibilities2)
# }
#
# sample_new_partition_p7 <- function(current.partition1, current.partition2){
#
# num.nodes <- length(current.partition1)
# present1 <- which(!is.na(current.partition1))
# present2 <- which(!is.na(current.partition2))
#
# new.partition1 <- rep(0,num.nodes)
# new.partition2 <- rep(0,num.nodes)
#
# # fill in the first with the second, and create isolates for people not present in the second
# new.partition1[present1] <- current.partition2[present1]
# max1 <- max(new.partition1, na.rm = TRUE)
# singletons1 <- which(is.na(new.partition1))
# if(length(singletons1) > 0) { # assign randomly the ones who are not present in the other one
# for(s in singletons1){
# newgroup <- sample(1:(max1+1),1)
# new.partition1[s] <- newgroup
# if(newgroup == (max1+1)) max1 <- max1 + 1
# }
# }
# new.partition1[present1] <- order_groupids(new.partition1[present1])
# new.partition1[new.partition1 == 0] <- NA
#
# # same the other way around
# new.partition2[present2] <- current.partition1[present2]
# max2 <- max(new.partition2, na.rm = TRUE)
# singletons2 <- which(is.na(new.partition2))
# if(length(singletons2) > 0) { # assign randomly the ones who are not present in the other one
# for(s in singletons2){
# newgroup <- sample(1:(max2+1),1)
# new.partition2[s] <- newgroup
# if(newgroup == (max2+1)) max2 <- max2 + 1
# }
# }
# new.partition2[present2] <- order_groupids(new.partition2[present2])
# new.partition2[new.partition2 == 0] <- NA
#
# return(list(new.partition1 = new.partition1,
# new.partition2 = new.partition2))
#
# }
#
#
# ## NEIGHBORHOOD PI 8: pick another partition, a group in it, and reproduce it in the current partition, or reorganizes it
# # CAREFUL, only works when singletons are allowed
# compute_size_neighborhood_p8 <- function(partition1, partition2){
#
# num.nodes <- length(partition1)
# present1 <- !is.na(partition1)
# present2 <- !is.na(partition2)
#
# # count number of groups from the second partition
# groups.toreorganize <- unique(partition2[present1 & present2])
# nums.reorganizations <- rep(0,length(groups.toreorganize))
#
# # count for each of them number of moves
# for(g in 1:length(groups.toreorganize)){
# group <- groups.toreorganize[g]
# members <- which(partition2 == group)
# max1 <- length(unique(partition1[present1 & partition1 != group]))
# for(k in 1:length(members)) {
# nums.reorganizations[g] <- nums.reorganizations[g] + Stirling(length(members),k)*(max1+1)^k
# }
# nums.reorganizations[g] <- nums.reorganizations[g] + 1 # when we just copy the group
# }
#
# return(list(groups.toreorganize = groups.toreorganize,
# nums.reorganizations = nums.reorganizations,
# num.reorganizations = sum(nums.reorganizations)))
# }
#
# sample_new_partition_p8 <- function(current.partition1, current.partition2, size_neighborhood){
#
# groups.toreorganize <- size_neighborhood$groups.toreorganize
# nums.reorganizations <- size_neighborhood$nums.reorganizations
# num.reorganizations <- size_neighborhood$num.reorganizations
#
# num.nodes <- length(current.partition1)
# present1 <- !is.na(current.partition1)
# present2 <- !is.na(current.partition2)
#
# new.partition <- current.partition1
#
# # pick
# probas_groups <- nums.reorganizations / num.reorganizations
# pick <- sample(c(0,groups.toreorganize),1,prob=c(1/ num.reorganizations, probas_groups))
#
# #randomize group
# if(pick > 0) {
# group <- pick
# members <- which(current.partition2 == group && present1)
# max1 <- length(unique(current.partition1[present1 & current.partition1 != group]))
# for(a in members){
# newgroup <- sample(1:(max1+1),1)
# new.partition[a] <- newgroup
# if(newgroup == (max1+1)) max1 <- max1 + 1
# }
# new.partition[present1] <- order_groupids(new.partition[present1])
# }
#
# # copy group
# if(pick ==0) {
# # reproduce the group
# new.partition[members] <- max(new.partition,na.rm = TRUE) + 1
# #reorder
# new.partition[present1] <- order_groupids(new.partition[present1])
# }
#
# return(new.partition)
#
# }
#
# #sample_new_partition_p8 <- function(current.partition1, current.partition2){
# #
# # num.nodes <- length(current.partition1)
# # present1 <- !is.na(current.partition1)
# # present2 <- !is.na(current.partition2)
# #
# # new.partition <- current.partition1
# #
# # # pick a group to reproduce
# # topick <- unique(current.partition2[present1 & present2])
# # g <- sample(topick,1)
# # members <- which(current.partition2 == g && present1)
# #
# # # reproduce the group
# # new.partition[members] <- max(new.partition,na.rm = TRUE) + 1
# #
# # #reorder
# # new.partition[present1] <- order_groupids(new.partition[present1])
# #
# # return(new.partition)
# #
# #}
#
# ## NEIGHBORHOOD PI 8: pick another partition, a group in it, and reproduce it in the current partition, or reorganizes it
# # CAREFUL, only works when singletons are allowed
# compute_size_neighborhood_p8_restricted <- function(partition1, partition2, sizes.simulated){
#
# num.nodes <- length(partition1)
# present1 <- !is.na(partition1)
# present2 <- !is.na(partition2)
# smax <- max(sizes.simulated)
#
# # count number of groups from the second partition
# groups.toreorganize <- unique(partition2[present1 & present2])
# combinations.pergroup <- list()
# nums.partitions.pergroup <- list()
# nums.reorganizations.pergroup <- list()
# nums.reorganizations <- rep(0,length(groups.toreorganize))
#
# # count for each of them number of moves
# for(g in 1:length(groups.toreorganize)){
# group <- groups.toreorganize[g]
# members <- which(partition2 == group)
# members <- which(1:num.nodes %in% members & present1)
#
# # if the group is not present, we will create it
# if(length(unique(partition1[members])) == 1 && sum(partition1[-members] == (partition1[members])[1], na.rm = TRUE) == 0) {
#
# nums.reorganizations.pergroup[[g]] <- NULL
# nums.reorganizations[g] <- 1
#
# # if the group is present, we will reshuffle
# } else {
# nmembers <- length(members)
# sizes <- table(partition1[-members])
#
# # find ways to divide the group
# combs <- parts(nmembers)
# combinations.pergroup[[g]] <- combs
#
# nums.p <- rep(1,ncol(combs))
# nums.r <- rep(1,ncol(combs))
# for(c in 1:ncol(combs)){
#
# # find number partitions of each
# ncpt <- nmembers
# for(s in 1:nrow(combs)){
# ncs <- sum(combs[,c] == s)
# nums.p[c] <- nums.p[c] * (choose(ncpt,s))^ncs / factorial(ncs)
# ncpt <- ncpt - ncs*s
# }
#
# # find ways to reorganize
# for(s in 1:nrow(combs)){
# ncs <- sum(combs[,c] == s)
# ncpt <- sum(sizes <= (smax-s))
# for(k in 1:ncs){
# nums.r[c] <- nums.r[c] * ncpt
# ncpt <- ncpt - 1
# }
# }
# }
#
# nums.partitions.pergroup[[g]] <- nums.p
# nums.reorganizations.pergroup[[g]] <- nums.r
# nums.reorganizations[g] <- nums.p * nums.r
# }
#
# }
#
# return(list(groups.toreorganize = groups.toreorganize,
# combinations.pergroup = combinations.pergroup,
# nums.partitions.pergroup = nums.partitions.pergroup,
# nums.reorganizations.pergroup = nums.reorganizations.pergroup,
# nums.reorganizations = nums.reorganizations,
# num.reorganizations = sum(nums.reorganizations) ))
# }
#
# sample_new_partition_p8_restricted <- function(current.partition1, current.partition2, size_neighborhood, sizes.simulated){
#
# groups.toreorganize <- size_neighborhood$groups.toreorganize
# combinations.pergroup <- size_neighborhood$combinations.pergroup
# nums.partitions.pergroup <- size_neighborhood$nums.partitions.pergroup
# nums.reorganizations.pergroup <- size_neighborhood$nums.reorganizations.pergroup
# nums.reorganizations <- size_neighborhood$nums.reorganizations
# num.reorganizations <- size_neighborhood$num.reorganizations
#
# num.nodes <- length(current.partition1)
# present1 <- !is.na(current.partition1)
# present2 <- !is.na(current.partition2)
#
# new.partition <- current.partition1
#
# # pick
# probas_groups <- nums.reorganizations / num.reorganizations
# group <- sample(groups.toreorganize,1,prob=probas_groups)
# members <- which(current.partition2 == group)
# members <- which(1:num.nodes %in% members & present1)
#
# # if the group is not present, we will create it
# if(length(unique(current.partition1[members])) == 1 && sum(current.partition1[-members] == (current.partition1[members])[1], na.rm = TRUE) == 0) {
#
# new.partition[members] <- max(new.partition,na.rm = TRUE) + 1
# new.partition[present1] <- order_groupids(new.partition[present1])
#
# # if the group is present, we will reshuffle
# } else {
#
# nmembers <- length(members)
# gleft <- unique(current.partition1[-members])
# gleft <- gleft[!is.na(gleft)]
# sizes <- apply(gleft, FUN = function(x) {return(sum(current.partition1[-members] == x))} )
#
# g <- which(groups.toreorganize == group)
# combs <- combinations.pergroup[[g]]
#
# # sample a combination
# probas_combs <- nums.partitions.pergroup[[g]] * nums.reorganizations.pergroup[[g]] / nums.reorganizations[g]
# c <- sample(1:ncol(combs),1,prob=probas_combs)
#
# # create a combination
# ntaken <- rep(F,nmembers)
# gtaken <- rep(F,length(gleft))
# for(s in nrow(combs)) {
# ncs <- sum(combs[,c] == s)
# for(k in 1:ncs) {
# ns <- choose(members[!ntaken], s)
# topick <- gleft[sizes <= (smax-s) & !gtaken]
# newg <- sample(c(0,topick),1)
# if(newg == 0){
# new.partition[ns] <- max(new.partition,na.rm = TRUE) + 1
# } else {
# new.partition[ns] <- newg
# gtaken[gleft == newg] <- TRUE
# }
# ntaken[members == ns] <- TRUE
# }
# }
# newnodes <- sample(1:nmembers,nmembers)
#
# new.partition[present1] <- order_groupids(new.partition[present1])
#
# }
#
# return(new.partition)
#
# }
Try the ERPM package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.