library(igraph)
##
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
## The following object is masked from 'package:base':
##
## union
g <- sample_gnm(12, 20)
plot(g) # A simple plot of the network - we'll talk more about plots later
dists = distances(g)
#Setting the defender node and states of all other nodes
d <- 2
n <- vcount(g) # sample size
state_vec <- sample(c(-1,0,1), replace=TRUE, size=n)
state_vec[d] <- 1 #Hardcoding the state of the defender to 1
#This is the function that takes as input the parent node and outputs
#the nodes in a 1-hop neighborhood i.e., the child nodes
theta1nodes <- function(x, y, z){
#X is the parent node
#Y is the distance matrix
#Z is the defender node
#W is the nodes in the current hop that must not be counted
idx <- y[x,]
idx_theta1 <- which(idx == 1)
idx_theta1 <- idx_theta1[!idx_theta1 %in% z]
#idx_theta1 <- idx_theta1[!idx_theta1 %in% w]
return(idx_theta1)
}
#Function for calculating expected states
#Expected states for child nodes X
#Using probability conversion matrix information
#For parent node Z
#Under scenario s
exputil <- function(x, z, s, p, q){
#Expected states for child nodes X
#Using probability conversion matrix information
#For parent node Z
#Under scenario s
s1 <- vector(mode = "numeric", length = length(x))
#Scenario 1 - Parent induces child to have parent's opinion
if(s==1){
muc <- p
mu2c <- q
for (i in 1:length(x)) {
#First term is tied to probability of change and second term is tied to probability of no change
if(z ==1){
if(x[i] == 1){
s1[i] <- 1
}
if(x[i] == 0){
s1[i] <- muc*(z) + (1-muc)*(x[i])
}
if(x[i] == -1){
s1[i] <- mu2c*(z) + (1-mu2c)*(x[i])
}
}
if(z ==0){
if(x[i] == 1){
s1[i] <- (muc)*(z) + (1-muc)*x[i]
}
if(x[i] == 0){
s1[i] <- 0
}
if(x[i] == -1){
s1[i] <- mu2c*(0) + (1-mu2c)*(x[i])
}
}
if(z ==-1){
if(x[i] == 1){
s1[i] <- mu2c*(z) + (1-mu2c)*x[i]
}
if(x[i] == 0){
s1[i] <- muc*(z) + (1-muc)*(x[i])
}
if(x[i] == -1){
s1[i] <- -1
}
}
}
}
#Scenario 2 - Parent induces child to have a positive opinion
if(s==2){
muc <- p
mu2c <- q
for (i in 1:length(x)) {
#First term is tied to probability of change and second term is tied to probability of no change
if(z ==1){
if(x[i] == 1){
s1[i] <- 1
}
if(x[i] == 0){
s1[i] <- muc*(x[i]+1) + (1-muc)*(x[i])
}
if(x[i] == -1){
s1[i] <- mu2c*(x[i]+1) + (1-mu2c)*(x[i])
}
}
if(z==0){
if(x[i] == 1){
s1[i] <- 1
}
if(x[i] == 0){
s1[i] <- 0
}
if(x[i] == -1){
s1[i] <- mu2c*(x[i]+1) + (1-mu2c)*(x[i])
}
}
if(z ==-1){
if(x[i] == 1){
s1[i] <- 1
}
if(x[i] == 0){
s1[i] <- muc*(x[i]+1) + (1-muc)*(x[i])
}
if(x[i] == -1){
s1[i] <- -1
}
}
}
}
#Scenario 3 - Parent induces child to have a negative opinion
if(s==3){
muc <- p
mu2c <- q
for (i in 1:length(x)) {
#First term is tied to probability of change and second term is tied to probability of no change
if(z ==1){
if(x[i] == 1){
s1[i] <- 1
}
if(x[i] == 0){
s1[i] <- muc*(x[i]-1) + (1-muc)*(x[i])
}
if(x[i] == -1){
s1[i] <- -1
}
}
if(z ==0){
if(x[i] == 1){
s1[i] <- (muc)*(x[i]-1) + (1-muc)*x[i]
}
if(x[i] == 0){
s1[i] <- 0
}
if(x[i] == -1){
s1[i] <- -1
}
}
if(z ==-1){
if(x[i] == 1){
s1[i] <- mu2c*(x[i]-1) + (1-mu2c)*x[i]
}
if(x[i] == 0){
s1[i] <- muc*(x[i]-1) + (1-muc)*(x[i])
}
if(x[i] == -1){
s1[i] <- -1
}
}
}
}
return(s1)
}
#Function for enumerating the neighborhoods (k=1,2,3...) and their nodes
#This function outputs a list with: (1) max number of hops given graph and defender node, (2) List of nodes per hop, (3) Total
#number of nodes per hop
khop <- function(d,dmat){
#d is the defender node
#dmat is the distances matrix
#First establishing the max k - limit of neighborhoods
idx <- dmat[d,]
midx <- max(idx) #The maximum distance from the defender - the k in the k-hop
s1 <- vector(mode = "numeric", length = midx) #Vector for holding number of nodes per hop
s2 <- list() #List for holding the exact nodes per hop
for (i in 1:midx) {
if(length(which(idx==i))>0){
s2[[i]] <- which(idx==i)
s1[i] <- length(which(idx==i))
}
}
#Creating easy return list
rs <- list()
rs[["hops"]] <- midx #Max number of possible hops
rs[["nodesperhop"]] <- s2 #List of nodes per hop - List
rs[["numnodesperhop"]] <- s1 #Total number of nodes per hop - Vector
return(rs)
}
#Main Script
#Getting the distance (hop) matrix and the adjacency matrix
hopmat <- distances(g)
adjmat <- data.matrix(get.adjacency(g))
#Starting the main loop
#First, finding out the number of hops and max hops
iter_data <- khop(d,hopmat)
nhops_d <- iter_data$hops
#Second, creating matrix to store the expected state data
expstate <- matrix(data = 0, nrow = nrow(hopmat), ncol = ncol(hopmat))
#Third, starting off with defender node
st_node <- d
istart <- st_node
if(length(istart) == 1 && istart == d){
#Extracting child nodes for defender node
t1 <- theta1nodes(istart, hopmat, d)
#Calculate and store the expected state
expstate[d, t1] <- expstate[d, t1] + exputil(state_vec[t1], state_vec[istart], 1, 0.75, 0.25)
#Going on to next hop
st_node <- t1
}
for (i in 1:(nhops_d-1)) {
istart <- st_node
if(length(istart)>1){
for (j in 1:length(istart)) {
#Extracting child nodes for each node in parent nodes set
t1 <- theta1nodes(istart[j], hopmat, d)
if(length(t1)>0){
expstate[istart[j], t1] <- expstate[istart[j], t1] + exputil(state_vec[t1], state_vec[istart[j]], 1, 0.75, 0.25)
}
}
st_node <- t1 #Initializing the next hop
}
}
#Initializing Cost Matrix for immediate theta 1 neighborhood
cij <- vector(mode = "numeric", length = vcount(g))
cij[iter_data$nodesperhop[[1]]] <- runif(length(iter_data$nodesperhop[[1]]), 0, 1)
#All nodes except defender
d_nodes <- iter_data$nodesperhop[[1]]
#Enumerating the overall path of impact when selecting one or the other node
impact_list <- list()
path_list <- list()
#Computing sum within and across hops
sumutil <- function(x, y, z){
s1 <- 0
if(length(x) == 1){
ix <- which(z[[1]] == x)
s1 <- s1 + sum(y[z[[1]][ix], z[[2]]])
}
else{
s1 <- s1 + sum(y[x, z[[2]]])
}
if(length(z)>2){
for (j in 2:(length(z)-1)) {
s1 <- s1 + sum(y[z[[j]]], z[[j+1]])
}
}
return(s1)
}
#Power Set function
pset <- function(set) {
n <- length(set)
masks <- 2^(1:n-1)
lapply( 1:2^n-1, function(u) set[ bitwAnd(u, masks) != 0 ] )
}
#Generating the power set
all_possible_choices <- pset(d_nodes)
all_possible_choices <- all_possible_choices[-1]
#Matrix for benefits
decision_mat_utility <- matrix(data = 0, nrow = length(all_possible_choices), ncol = 1)
for (i in 1:length(all_possible_choices)) {
decision_mat_utility[i,1] <- sumutil(all_possible_choices[[i]], expstate, iter_data$nodesperhop)
}
#Decision not considering cost
dk1 <- which(decision_mat_utility == max(decision_mat_utility))
for (i in 1:length(dk1)) {
print("Choice")
print(all_possible_choices[[dk1[i]]])
print("Cost:")
print(sum(cij[all_possible_choices[[dk1[i]]]]))
print("________________")
}
## [1] "Choice"
## [1] 1 10 11
## [1] "Cost:"
## [1] 0.8303753
## [1] "________________"
## [1] "Choice"
## [1] 1 3 10 11
## [1] "Cost:"
## [1] 1.601867
## [1] "________________"