######################################################################
# STATNET SCRIPT - March 21st 2023                                       #
######################################################################


# Getting Started with network and Igraph packages (1st part)
######################################################################
## 1) Finding statnet built-in data sets

library(statnet)
data(package="network")
data(package="sna")

data(flo) # Load a built-in data set; see ?flo for more
flo # Examine the flo adjacency matrix
class(flo)
?flo #more info
?emon

data("faux.magnolia.high")
fmh <- faux.magnolia.high
?faux.magnolia.high

######################################################################
## 2) Importing Relational Data
#Be sure to be in the directory where you stored the data
#You can use setwd() to change it, or platform specific shortcuts
getwd() # Check which working dir are you pointing
list.files() # Check what is the working dir content

#Read an adjacency matrix (R stores it as a data frame by default)
relations <- read.csv("relationalData.csv",header=FALSE,stringsAsFactors=FALSE, sep=",")
relations
is.data.frame(relations)

# the matrix class
relations <-as.matrix(relations) # Convert in matrix

# Import attributes (this can be left as dataframe)
nodeInfo <- read.csv("vertexAttributes.csv",header=TRUE,stringsAsFactors=FALSE, sep=",")
nodeInfo

# relational data has no row/column names by default, it is better to set them
rownames(relations)<-nodeInfo$name
colnames(relations)<-nodeInfo$name
relations

#relations is weighted. Let dichotomize it (for now)

relations[relations>1] <- 1

######################################################################
## 3) objects of class "network"

nempty <- network.initialize(5) # Create an empty graph with 5 vertices
nempty

randomG <- rgraph(n=10, m=1, tprob=0.5, mode="digraph") # Create a random graph with n=10 nodes and tie prob = 0.5

# Let's transform the imported data into a network object
nrelations<-network(relations,directed=FALSE) 
nrelations # Get a quick description of the data
is.network(nrelations)
class(nrelations)

######################################################################
## Initializing networks
##
net <- network.initialize(5)
net

nmat <- matrix(rbinom(25, 1, 0.5), nr = 5, nc = 5)
net <- network(nmat, loops = TRUE)
net
net[,]
summary(net)

net <- as.network(nmat, loops = FALSE)
all(nmat == net[,])

net <- as.network(nmat, loops = TRUE)
all(nmat == net[,])

# recover relational data from matrices within a network object
net <- network.initialize(5, loops = TRUE)
net[,]
net[nmat>0] <- 1
all(nmat == net[,])

net[,] <- 0
net[,] <- nmat
all(nmat == net[,])

net[,] <- 0
network::add.edges(net, row(nmat)[nmat>0], col(nmat)[nmat>0])
?add.edges
all(nmat == net[,])

net[,] <- 0
net[,] <- nmat[,]
all(nmat == net[,])

######################################################################
##### Basic counting

summary(nrelations) # Get an overall summary
is.directed(nrelations)
network.dyadcount(nrelations) # How many dyads in nrelations? 
network.edgecount(nrelations) # How many edges are present?
network.size(nrelations) # How many nodes?
class(nrelations)
sna::dyad.census(nrelations)
sna::triad.census(nrelations)

##############################################################################
# Some elementary operations on networks
# pointing, deleting, subsetting

#let's fill the empty network. First actor is connected with everyone

nempty[1,] <- 1
#which network is this?

#gplot(nempty)

nrelations
nrelations[,]
nrelations[8,9] # UNIFI to UNIPI?
nrelations[8,] # Entire UNIFI row
nrelations[1:4,5:8] # Subsets
nrelations[-9,-9] # deleting UNIPI row


# Add vertex attributes to a network object

nrelations %v% "id" <- nodeInfo$id # Add in our vertex attributes
nrelations %v% "numProj" <- nodeInfo$NumProject
nrelations %v% "type" <- nodeInfo$Type
nrelations %v% "prov" <- nodeInfo$Prov
nrelations %v% "firstProj" <- nodeInfo$firstProj

#Listing attributes

network::list.vertex.attributes(nrelations) # List all vertex attributes

#Retrieving attributes

nrelations %v% "vertex.names" # Retrieve vertex name
nrelations %v% "prov" # Retrieve vertex prov
nrelations %v% "firstProj" # Retrieve vertex first project
nrelations %v% "type" # Retrieve vertex type

nrelations

# Get induced subgraphs

nrelations %s% which(nrelations %v% "type"=="PRO")


sub_nrelations <- nrelations %s% which(nrelations %v% "type"=="PRO")

plot(sub_nrelations, displaylabels=TRUE)

plot(nrelations, displaylabels=TRUE)

sub_nrelations <- nrelations %s% which(nrelations %v% "type"=="PRO" | nrelations %v% "type"=="I")

View(data.frame(nrelations %v% "vertex.names", nrelations %v% "type"))

######################################################################
##### Edgelist class

#The sna package also supports a special kind of matrix called an 'sna edgelist'. 
#These are three-column matrices, each row of which represents an edge (via its sender, recipient, and value, respectively). 
#These 'sna edgelists' have special attributes that indicate their size, vertex names (if any), and bipartite status (if applicable).

erelations <-as.edgelist.sna(nrelations) # Coerce nrelations to an sna edgelist
erelations

get.dyads.eids(nrelations, tails=7, heads=9)
get.dyads.eids(nrelations, tails=9, heads=7)

attr(erelations,"n") # How many vertices are there?
attr(erelations,"vnames") # Are there vertex names?
b<-as.sociomatrix.sna(erelations) # Can transform back with as.sociomatrix.sna
class(b)

# 'sna edgelists' can be handy with large data sets (as a simple alternative to network objects). 

# (induced) subgraph with more than 1 attribute (works in a ego-network sense)

sub_nrelations2 <- nrelations %s% list(which(nrelations %v% "type"=="PRO"), which(nrelations %v% "type"=="I"))

plot(sub_nrelations2, displaylabels=TRUE)

sub_nrelations2 <- get.inducedSubgraph(nrelations, v=c(3,8,9))

plot(sub_nrelations2, displaylabels=TRUE)

sub_nrelations2 <- get.inducedSubgraph(nrelations, eid = 1:3)

plot(sub_nrelations2, displaylabels=TRUE)
plot(nrelations, displaylabels=TRUE)

######################################################################
##### Import 2-mode networks or incidence matrices

davis<-read.table("Davis.txt",header = T, row.names = 1)
davis <-as.matrix(davis)
gplot(davis, gmode="twomode",displaylabels = T, label.cex = 0.7)
?gplot

Arow <- davis %*% t(davis) #transform into a one-mode network
Arow <- tcrossprod(davis) #faster function
Acol <- tcrossprod(t(davis))

gplot(Arow, gmode="graph",displaylabels = T, label.cex = 0.7, edge.lwd=Arow)

#dichotomize
Arow[Arow>1] <- 1

## load GOT data

load("netgot.RData")

ls()
netgot
netgot[,]
netgot %v% "vertex.names"
network::list.vertex.attributes(netgot) # List all vertex attributes

Ego <- ego.extract(netgot, ego =1)
gplot(as.network(Ego$`Waymar Royce`), displaylabels = T)