MST_Script.R

library(igraph)

## 
## Attaching package: 'igraph'

## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum

## The following object is masked from 'package:base':
## 
##     union

library(leaflet)
library(geosphere)
library(optrees)

# *****************************************************************************
#                       R Script For Trasportation Network
# *****************************************************************************

# clears all objects in "global environment"
rm(list=ls())

set.seed(123)

# ************************ User defined functions ******************************
# Tests whether row vectors are within a matrix; identifies rows in matrix accordingly 
vectorsInVectorSeq <- function(testVectorMatrix, seqVectorMatrix) {
  r <- rep(FALSE,nrow(seqVectorMatrix))
  for (k in 1:nrow(testVectorMatrix)) {
    testElement = testVectorMatrix[k,]
    r <- r | apply(seqVectorMatrix, 1, function(x, test) isTRUE(all.equal(x, test)), testElement)
  }
  return(r)
}

# create distance matrix
mat2list <- function(D) {
  n = dim(D)[1]
  k <- 1
  e <- matrix(ncol = 3,nrow = n*(n-1)/2)
  for (i in 1:(n-1)) {
    for (j in (i+1):n) {
      e[k,] = c(i,j,D[i,j])
      k<-k+1
    }
  }
  return(e)
}

# ************************ End of user defined functions ******************************

# *************************************************************************************
#                               Start of main function
# *************************************************************************************

# 1.) identify lontitude and latitude of distributed centers
# latitude and lontitude is referenced from https://www.latlong.net/
lonlat <- read.table(textConnection( "lon, lat 
                                     9.189982,45.464203
                                     9.993682,53.551086 
                                     -0.127758,51.507351
                                     5.369780,43.296482
                                     2.173404,41.385063
                                     19.455982,51.759247
                                     4.477733,51.924419
                                     -9.139337,38.722252
                                     4.402464,51.219448
                                     11.974560,57.708870 
                                     "),header=TRUE,strip.white = TRUE, sep=",")

nname <- c("Milan","Hamburg","London","Marseille","Barcelona", 
          "Lodz","Rotterdam","Lisbon", 
          "Antwerp","Gothenburg")

v <- data.frame( ids = 1:10, 
                name = nname, 
                x = lonlat$lon, 
                y = lonlat$lat) 

leaflet(data=v[1:10,]) %>% addTiles() %>%
  addMarkers(lng=v$x, lat=v$y, popup=v$name)

leaflet(data=v[1:10,]) %>% addTiles() %>%
  addMarkers(lng=v$x, lat=v$y, popup=v$name) %>%
  addPolylines(lng=v$x, lat=v$y, popup=v$name)

# 2.) create distance in map 

D <- distm(lonlat, lonlat, fun=distVincentyEllipsoid) # in meters
eD <- mat2list(D/1000) #km
eD <- eD[order(eD[,3]),] #sort edges

# 3.) create Minimum Spanning Tree (Kruskal's algroithm)
mst <- getMinimumSpanningTree(v$name, round(eD), algorithm = "Kruskal",
                              show.data = TRUE, show.graph = TRUE, check.graph = FALSE)

## 
##  Minimum Cost Spanning Tree 
##  Algorithm: Kruskal 
##  Stages:  22 | Time:  0.001 
##  ------------------------------
##       ept1     ept2    weight 
##          7        9        79
##          3        9       317
##          4        5       338
##          1        4       388
##          2        7       414
##          2       10       479
##          2        6       670
##          1        9       731
##          5        8      1009
##  ------------------------------
##                    Total = 4425

# 4.) create a graph of MST
net <- graph.data.frame(eD[,1:2],directed = FALSE, vertices = v)
E(net)$weight <- eD[,3]
mstgraph <- minimum.spanning.tree(net)
par(mfrow=c(1,2), mar=c(0,1,0.75,0)) 
plot(net,vertex.label=nname)
plot(mstgraph,vertex.shape="none",edge.label=round(E(mstgraph)$weight))

# 5.) create k cluster from full mst graph
k <- 3
e.mst  = get.edges(mstgraph,1:ecount(mstgraph)) # matrix of edges of MST
e.clust = get.edges(mstgraph,1:(ecount(mstgraph)-k+1)) # matrix of edges of k cluster
clust_idx = vectorsInVectorSeq(e.clust, e.mst) # row indicis cluster in MST
ecol <- rep("grey80", ecount(mstgraph)) # default edge colour
ecol[clust_idx] <- "green" # colour for clust
ew <- rep(2,ecount(mstgraph)) # default edge width
ew[clust_idx] <- 5 # width for clust edges
par(mfrow=c(1,1), mar=c(0,1,0.75,0)) 
plot(mstgraph, vertex.shape="none",vertex.label.cex=0.85,
     edge.color=ecol, edge.width=ew,edge.label.cex=0.85,
     edge.label=round(E(mstgraph)$weight))

print("end")

## [1] "end"