For our last course on flows and networks, I have uploaded additional material. Slides are now online

A markdown-based webpage is now online, with all the codes.

For our last course on flows and networks, I have uploaded additional material. Slides are now online

A markdown-based webpage is now online, with all the codes.

This morning, we had our first practicals on network flows, using an example mentioned in some papers published by Noraini Abdullah and Ting Kien Hua, max flow min cut theorem to minimize traffic congestion in Kota Kinabalu and application of the Shortest Path and Maximum Flow with Bottleneck in Traffic Flow of Kota Kinabalu. From the roads mentioned in the articles, I did try my best to locate the nodes on a map,

`m=matrix(c(0,5.995910, 116.105520,`

1,5.992737, 116.093718,

2,5.992066, 116.109883,

3,5.976947, 116.095760,

4,5.985766, 116.091580,

5,5.988940, 116.080112,

6,5.968318, 116.080764,

7,5.977454, 116.075460,

8,5.974226, 116.073604,

9,5.969651, 116.073753,

10,5.972341, 116.069270,

11,5.978818, 116.072880),3,12)

we can be visualized below

`library(OpenStreetMap)`

map = openmap(c(lat= 6.000, lon= 116.06),

c(lat= 5.960, lon= 116.12))

map=openproj(map)

plot(map)

points(t(m[3:2,]),col="black", pch=19, cex=3 )

text(t(m[3:2,]),c("s",1:10,"t"),col="white")

If the source is realistic (up north), I do not feel very confortable with the location of the sink (on the west). But let’s pretend it’s find (to do the maths, at least).

To extract information about edge capacity, on that network use the following code that will extract the three tables from the paper

`library(devtools)`

install_github("ropensci/tabulizer")

library(tabulizer)

location <- 'http://www.jistm.com/PDF/JISTM-2017-04-06-02.pdf'

out <- extract_tables(location)

with Windows, it seems to be necessary to download another package first

`library(devtools)`

install_github("ropensci/tabulizerjars")

install_github("ropensci/tabulizer")

library(tabulizer)

location <- 'http://www.jistm.com/PDF/JISTM-2017-04-06-02.pdf'

out <- extract_tables(location)

Now we can get out data frame with capacities

`B1=as.data.frame(out[[2]])`

B2=as.data.frame(out[[3]])

E=data.frame(from=B1[3:20,"V3"],

to=B1[3:20,"V4"])

E=E[-c(6,8),]

capacity=as.character(B2$V3[-1])

capacity[6]="843"

capacity[4]="2913"

E$capacity=as.numeric(capacity)

We can add those edges on our map (without the arrows to indicate the direction, it would be to heavy to read)

`plot(map)`

points(t(m[3:2,]),col="black", pch=19, cex=3 )

B=data.frame(i=as.character(c("s",paste("V",1:10,sep=""),"t")),

x=m[3,],y=m[2,])

for(i in 1:nrow(E)){

i1=which(B$i==as.character(E$from[i]))

i2=which(B$i==as.character(E$to[i]))

segments(B[i1,"x"],B[i1,"y"],B[i2,"x"],B[i2,"y"],lwd=3)

}

text(t(m[3:2,]),c("s",1:10,"t"),col="white")

To get the graph with capacities, an alternative is to use

`library(igraph)`

g=graph_from_data_frame(E)

E(g)$label=E$capacity

plot(g)

but it does not respect geographical locations of nodes. It can actually be done using

`plot(g, layout=as.matrix(B[,c("x","y")]))`

To get a better understanding of the capacities of the road, use

`plot(g, layout=as.matrix(B[,c("x","y")]),`

edge.width=E$capacity/200)

From that network with capacities, the goal is to determine maximum flow on that network, from the source to the sink. This can be done with R using

`> (m=max_flow(graph=g, source="s", target="t"))`

$value

[1] 2571

`$flow`

[1] 1191 1380 1422 1380 231 0 231 0 1149 1422 1149 0 0 1149 1422

[16] 1149

Our maximum flow is here 2571, which is different from was is actually claimed both in the two papers max flow min cut theorem to… and application of the Shortest Path… (“*the maximum flow for the capacitated network with 12 nodes and 16 edges of the selected scope in this study was 2598 vehicles per hour*“) where there are clearly typos since values in the table and on the graph are different. Here I did use the ones from the tables.

`E$flux1=m$flow`

E(g)$label=E$flux1

plot(g, layout=as.matrix(B[,c("x","y")]),

edge.width=E$flux1/200)

That is nice, but rather odd. Actually, a much simpler flow can be considered, but the same global value

`E$flux2=c(1422,1149,1422,1149,0,0,0,0,`

1149,1422,1149,0,0,1149,1422,1149)

E(g)$label=E$flux2

plot(g, layout=as.matrix(B[,c("x","y")]),

edge.width=E$flux2/200)

Nice, isn’t it. It is actually possible to do exactly the same on another paper they have, on the same city, traffic congestion problem of road networks in Kota Kinabalu.

`location <- 'http://www.worldresearchlibrary.org/up_proc/pdf/999-150486366625-30.pdf'`

out <- extract_tables(location)

dim(out[[3]])

B1=as.data.frame(out[[3]])

E=data.frame(from=B1[2:61,"V2"],

to=B1[2:61,"V3"],

capacity=B1[2:61,"V4"])

E$capacity=as.numeric(

as.character(E$capacity))

library(igraph)

g=graph_from_data_frame(E)

m=max_flow(graph=g,

source="S",

target="T")

E$flux1=m$flow

E(g)$label=E$flux1

plot(g,

edge.width=E$flux1/200,

edge.arrow.size=0.15)

Here the value of the maximal flow is 4017, just as they found in the original paper

For practicals on networks and flows, we will use the R package flows dedicated to flows on networks

`library(flows)`

data(nav)

myflows <- prepflows(mat = nav, i = "i", j = "j", fij = "fij")

diag(myflows) <- 0

Select flows that represent at least 20% of the sum of outgoing flows for each urban area.

`flowSel1 <- firstflows(mat = myflows/rowSums(myflows)*100, method = "xfirst",k = 20)`

Then select the dominant flows (incoming flows criterion)

`flowSel2 <- domflows(mat = myflows, w = colSums(myflows), k = 1)`

flowSel <- myflows * flowSel1 * flowSel2

inflows <- data.frame(id = colnames(myflows), w = colSums(myflows))

and finally plot dominant flows map

`opar <- par(mar = c(0,0,2,0))`

sp::plot(GE, col = "#cceae7", border = NA)

plotMapDomFlows(mat = flowSel, spdf = UA, spdfid = "ID", w = inflows, wid = "id",wvar = "w", wcex = 0.05, add = TRUE,legend.flows.pos = "topright",legend.flows.title = "Nb. of commuters")

title("Dominant Flows of Commuters")

The code to get the background map is based on the GE object, defined in that package.

To go further on dominant flows read Nystuen & Dacey (1961)

We will discuss in the last course, next week two extensions that were not mentioned in the course. The first one is about congestion models. The second one is a nice application of flow to discuss sports issues in NBA (or NHL).

For the second practicals of our course on networks and flows, we will study traffic flow of Kota Kinabalu (Malaysia), following several papers published by Noraini Abdullah and Ting Kien Hua, such as max flow min cut theorem to minimize traffic congestion in Kota Kinabalu, traffic congestion problem of road networks in Kota Kinabalu and application of the Shortest Path and Maximum Flow with Bottleneck in Traffic Flow of Kota Kinabalu.

After the Traveling Salesman part, we will see tools used to study flows in networks this week. Slides are now online (from slide 42).

Slides are probably full of typos, not to say mistakes. All apologies.

I will also include some slides on the simplex method (since I am not sure that all students have seen it)

Demain matin, nous aurons un TP pour le cours sur les réseaux, les flux et les transports. En particulier, en nous inspirant des travaux de Sybil Derrible, nous allons commencer par étudier la centralité dans les différents systèmes de métro, mais aussi la robustesse. Les matrices d’adjacence d’une trentaine de métros dans le monde sont en ligne dans un fichier xls. Histoire de gagner un peu de temps, le code pour créer une matrice d’adjacence peut être le suivant

`loc="/data/Metro_Networks_Adjacency.xls"`

library(xlsx)

E=read.xlsx(loc,"StPetersburg")

n=nrow(E)

nom=as.character(E[3:(n-2),1])

Adj=E[3:(n-2),(4:ncol(E)-1)]

Adj[is.na(Adj)]=0

Adj=as.matrix(Adj)

colnames(Adj)=rownames(Adj)=nom

On est ensuite prêt à manipuler le réseau,

`library(igraph)`

iflo=graph_from_adjacency_matrix(Adj,mode = "undirected")

plot(iflo)

On va utiliser les notions vues en cours, sur la centralité, mais surtout, on travaillera sur Quantifying the robustness of metro networks, inspiré de The complexity and robustness of metro networks. Plusieurs fonctions utiles sont déjà programmées dans R, comme l’assortativité.

In the second part of the course on graphs and networks, we will focus on economic applications, and flows. The first series of slides are on the traveling salesman problem. Slides are available online.

In order to practice with network data with R, we have been playing with the Padgett (1994) Florentine’s wedding dataset (discussed in the lecture). The dataset is available from

`> library(network)`

> data(flo)

> nflo=network(flo,directed=FALSE)

> plot(nflo, displaylabels = TRUE,

+ boxed.labels =

+ FALSE)

The next step was to move from the network package to igraph. Since we have the adjacency matrix, we can use it

`> library(igraph)`

> iflo=graph_from_adjacency_matrix(flo,

+ mode = "undirected")

> plot(iflo)

The good thing is that a lot of functions are available, for instance we can get shortest paths, between two specific nodes. And we can give appropriate colors to the nodes that we’ll cross

`> AP=all_shortest_paths(iflo,`

+ from="Peruzzi",

+ to="Ginori")

> L=AP$res[[1]]

> V(iflo)$color="yellow"

> V(iflo)$color[L[2:4]]="light blue"

> V(iflo)$color[L[c(1,5)]]="blue"

> plot(iflo)

We can also visualize edges, but I found it slightly more complicated (to extract edges from the output)

`> liens=c(paste(as.character(L)[1:4],`

+ "--",

+ as.character(L)[2:5],sep=""),

+ paste(as.character(L)[2:5],

+ "--",

+ as.character(L)[1:4],sep=""))

> df=as.data.frame(ends(iflo,E(iflo)))

> names(df)=c("src","target")

> lstn=sort(unique(c(as.character(df[,1]),as.character(df[,2]),"Pucci")))

> Eliens=paste(as.numeric(factor(df[,1],levels=lstn)),"--",

+ as.numeric(factor(df[,2],levels=lstn)),sep="")

> EU=unlist(lapply(Eliens,function(x) x%in%liens))

> E(iflo)$color=c("grey","black")[1+EU]

> plot(iflo)

But it works. It is also possible to use some D3js visualization

`> library( networkD3 )`

> simpleNetwork (df)

Then the next question was to add a vertice to the network. The most simple way to do it is probability through the adjacency matrix

`> flo2=flo`

> flo2["Pucci","Bischeri"]=1

> flo2["Bischeri","Pucci"]=1

> nflo2=network(flo2,directed=FALSE)

> plot(nflo2, displaylabels = TRUE,

+ boxed.labels =

+ FALSE)

Then, we’ve been playing with centrality measures.

`> plot(iflo,vertex.size=betweenness(iflo))`

The goal was to see how related they were. Here, for all of them, “Medici” is the central node. But what about the others?

`> B=betweenness(iflo)`

> C=closeness(iflo)

> D=degree(iflo)

> E=eigen_centrality(iflo)$vector

> base=data.frame(betw=B,close=C,deg=D,eig=E)

> cor(base)

betw close deg eig

betw 1.0000000 0.5763487 0.8333763 0.6737162

close 0.5763487 1.0000000 0.7572778 0.7989789

deg 0.8333763 0.7572778 1.0000000 0.9404647

eig 0.6737162 0.7989789 0.9404647 1.0000000

Those measures are quite correlated. It is also possible to use a hierarchical graph to visualize how close those centrality measures can be

`> H=hclust(dist(t(base)),`

+ method="ward")

> plot(H)

Instead of looking at values of centrality measures, it is possible to looks are ranks

`> rbase=base`

> for(i in 1:4) rbase[,i]=rank(base[,i])

> H=hclust(dist(t(rbase)),

+ method="ward")

> plot(H)

Here the eigenvector measure is very close to the degree of vertices.

Finally, it is possible to seek clusters (in the context of coalition here, in case a war should start between those families)

`> kc <- fastgreedy.community ( iflo )`

Here we have 3 classes (+1 for the node that is disconnected from the other families)

`> V(iflo)$color=c("yellow","orange",`

+ "light blue")[membership ( kc )]

> plot(iflo)

`> plot(kc,iflo)`

This week, we will start the first graduate course on graphs and networks. Slides are available online.

Later on, we will see additional applications, in the context of flows, matching, transportation, etc.