devtools::load_all (".", export_all = FALSE)## Loading osmprob
#devtools::load_all (".")
#devtools::document (".")
#devtools::check (".")First set up a simple network and call the Rcpp function, which in this version dumps the initial and final transition probability matrices, and returns the shortest path (calculated by simply Dijkstra).
netdf <- data.frame (
xfr = c (rep (0, 3), rep (1, 3), rep (2, 4),
rep (3, 3), rep (4, 2), rep (5, 3)),
xto = c (1, 2, 5, 0, 2, 3, 0, 1, 3, 5,
1, 2, 4, 3, 5, 0, 2, 4),
d = c (7., 9., 14., 7., 10., 15., 9., 10., 11., 2.,
15., 11., 6., 6., 9., 14., 2., 9.))Then a routine to visualise the network
make_polygon <- function (netdf)
{
# arrange polygon points around a circle
n <- length (unique (netdf$xfr)) - 1
xy <- data.frame ('x' = cos (2 * pi * seq (n) / n),
'y' = sin (2 * pi * seq (n) / n))
# then insert third point in the middle
rbind (xy [1:2,], c (0, 0), xy [3:n,])
}
plot_net <- function (netdf, lines=TRUE, mt='Network weights', lwdmax=10)
{
xy <- make_polygon (netdf)
#plot.new ()
par (mar=c(0, 0, 1, 0))
plot (xy$x, xy$y, pch=1, xlim = c (-1.1, 1.1), ylim = c (-1.1, 1.1),
xaxt="n", yaxt="n", xlab="", ylab="", frame=FALSE)
title (main=mt)
text (xy$x, xy$y, labels=0:5, pos=4, cex=2)
xy1 <- xy [netdf$xfr + 1, ]
xy2 <- xy [netdf$xto + 1, ]
lwd <- netdf$d * lwdmax / max (netdf$d)
if (lines)
for (i in seq (nrow (netdf)))
lines (c (xy1 [i, 1], xy2 [i, 1]),
c (xy1 [i, 2], xy2 [i, 2]), lwd = lwd [i])
idx <- c (1, 5) # start and end nodes
points (xy$x [idx], xy$y [idx], pch=19, col=c ("green", "red"), cex=2)
}plot_net (netdf)
This is Algorithm#1 of Saerens et al (Neural Computation 2009).
dmat <- array (NA, dim=c(6, 6))
dmat [netdf$xf * 6 + netdf$xto + 1] <- netdf$d
# add start node as first column
dmat <- cbind (rep (NA, 6), dmat)
dmat <- rbind (c (NA, 1, rep (NA, 5)), dmat)
qmat <- dmat
qmat [!is.na (qmat)] <- 1
qmat [is.na (qmat)] <- 0
qmat <- qmat / rowSums (qmat)
# Add the link to the absorbing destination state:
dest <- 6 # node#4 with 0-index is 5, plus one row = 6
n <- length (which (qmat [dest, ] > 0))
qmat [dest, ] <- qmat [dest, ] * n / (n + 1)
nmat <- solve (diag (nrow (qmat)) - qmat)The final nmat calculated above holds the expected number of trips along each edge, which this function then plots (after removing the first row and first column of nmat).
plot_ntrips <- function (netdf, nmat, mt='Number of trips', lwdmax=10)
{
n <- length (unique (netdf$xfr))
nmat_sm <- nmat [2:(n + 1), 2:(n + 1)]
nmat_sm <- nmat_sm * lwdmax / max (nmat_sm)
plot_net (netdf, lines=FALSE, mt=mt)
xy <- make_polygon (netdf)
xy1 <- xy [netdf$xfr + 1, ]
xy2 <- xy [netdf$xto + 1, ]
for (i in seq (nrow (netdf)))
lines (c (xy1 [i, 1], xy2 [i, 1]),
c (xy1 [i, 2], xy2 [i, 2]),
lwd = nmat_sm [netdf$xfr [i] + 1, netdf$xto [i] + 1])
idx <- c (1, 5) # start and end nodes
points (xy$x [idx], xy$y [idx], pch=19, col=c ("green", "red"), cex=2)
}
par (mfrow=c(1,2))
plot_net (netdf)
plot_ntrips (netdf, nmat)
Then back on to the subsequent steps of the routing algorithm.
lq <- t (-log (qmat))
lq [!is.finite (lq)] <- 0
hvec <- diag (qmat %*% lq)
xvec <- nmat %*% hvecCalculation of the v-vector then requires multiplying Q by D, but dmat has NA values for non-finite costs. These first have to be re-set to 0 prior to the multiplication. Also note that the final component of Eq.(7) of Saerens et al. (2009) is
Why is is then included at all in the Algorithm?
dtemp <- dmat
dtemp [is.na (dtemp)] <- 0
vvec <- nmat %*% diag (qmat %*% t (dtemp))That finally suffices to iterate the transition probability matrix, which can be done on the rows of qmat
qmat0 <- qmat
eta <- 8 # The entropy parameter
for (i in seq (nrow (qmat)))
{
indx <- which (qmat [i, ] > 0)
rowi <- exp (-(1/eta) * (qmat [i, indx] + vvec [indx]) + xvec [indx])
qmat [i, indx] <- rowi / sum (rowi, na.rm=TRUE)
}Then compare the two probability matrices
qmat0## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [2,] 0 0.0000000 0.3333333 0.3333333 0.0000000 0.0000000 0.3333333
## [3,] 0 0.3333333 0.0000000 0.3333333 0.3333333 0.0000000 0.0000000
## [4,] 0 0.2500000 0.2500000 0.0000000 0.2500000 0.0000000 0.2500000
## [5,] 0 0.0000000 0.3333333 0.3333333 0.0000000 0.3333333 0.0000000
## [6,] 0 0.0000000 0.0000000 0.0000000 0.3333333 0.0000000 0.3333333
## [7,] 0 0.3333333 0.0000000 0.3333333 0.0000000 0.3333333 0.0000000
qmat## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [2,] 0 0.0000000 0.2414607 0.4079540 0.0000000 0.0000000 0.3505853
## [3,] 0 0.2986676 0.0000000 0.4353070 0.2660255 0.0000000 0.0000000
## [4,] 0 0.2496312 0.2153484 0.0000000 0.2223484 0.0000000 0.3126720
## [5,] 0 0.0000000 0.2495663 0.4216485 0.0000000 0.3287852 0.0000000
## [6,] 0 0.0000000 0.0000000 0.0000000 0.4155887 0.0000000 0.5844113
## [7,] 0 0.2782418 0.0000000 0.4055365 0.0000000 0.3162217 0.0000000
The above steps can then be pulled together as the steps necessary to iterate until convergence of qmat
delta <- sum (abs (qmat - qmat0))
tol <- 1e-3
n <- 2 # already done 1 step above
while (delta > tol)
{
qmat_old <- qmat
lq <- t (-log (qmat))
lq [!is.finite (lq)] <- 0
hvec <- diag (qmat %*% lq)
xvec <- nmat %*% hvec
vvec <- nmat %*% diag (qmat %*% t (dtemp))
for (i in seq (nrow (qmat)))
{
indx <- which (qmat [i, ] > 0)
rowi <- exp (-(1/eta) * (qmat [i, indx] + vvec [indx]) + xvec [indx])
qmat [i, indx] <- rowi / sum (rowi, na.rm=TRUE)
}
delta <- sum (abs (qmat - qmat_old))
n <- n + 1
}
message ('converged in ', n, ' steps')## converged in 11 steps
#png (filename="probnet.png", width=800, height=400)
par (mfrow=c(1,2))
plot_net (netdf)
plot_ntrips (netdf, qmat, mt='Probability')
#graphics.off ()Finally, compare the above output with the result of the Rcpp function, which in this version dumps the initial and final transition probability matrices, and returns the shortest path (calculated by simply Dijkstra).
qmat## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [2,] 0 0.0000000 0.2330359 0.4110607 0.0000000 0.0000000 0.3559033
## [3,] 0 0.2915579 0.0000000 0.4541876 0.2542545 0.0000000 0.0000000
## [4,] 0 0.2434057 0.2153019 0.0000000 0.2121036 0.0000000 0.3291888
## [5,] 0 0.0000000 0.2569678 0.4522957 0.0000000 0.2907365 0.0000000
## [6,] 0 0.0000000 0.0000000 0.0000000 0.3946412 0.0000000 0.6053588
## [7,] 0 0.2808723 0.0000000 0.4378515 0.0000000 0.2812762 0.0000000
osm_router (netdf, 0, 4, eta=8)## ---converged in 19 loops
## ------ Q1_MAT ------
## S 0 1 2 3 4 5 E
## 0 1.0000 0 0 0 0 0
## 0 0 0.2330 0.4111 0 0 0.3559
## 0 0.2916 0 0.4542 0.2542 0 0
## 0 0.2434 0.2153 0 0.2121 0 0.3292
## 0 0 0.2570 0.4523 0 0.2908 0
## 0 0 0 0 0.3946 0 0.6054
## 0 0.2809 0 0.4378 0 0.2813 0
## [,1] [,2]
## [1,] 0 7
## [2,] 2 9
## [3,] 5 20
## [4,] 4 20
