# R source code accompanying Einbeck & Dwyer, Transportmetrica

#
setwd("/home/stats/dma0je/public_html/lpc/traffic/")   # ignore!

setwd("../traffic")  # complete here the path name where you have extracted traffic.zip

# LPC source code
source("lpc032a.r")
source("lpctools.proj3.r")

# Additional R packages, from CRAN
library(princurve)
library(pspline)

# 1 mile = 1.609344 kilometers
v <-  1.609344

# Data

# VDS 716450
data06262007    <-  read.table("06262007.dat", header=TRUE)
data07032007    <-  read.table("07032007.dat", header=TRUE)
data07092007    <-  read.table("07092007.dat", header=TRUE)
data07102007    <-  read.table("07102007.dat", header=TRUE)
data07112007    <-  read.table("07112007.dat", header=TRUE)
data07122007    <-  read.table("07122007.dat", header=TRUE)
data07132007    <-  read.table("07132007.dat", header=TRUE)
data07142007    <-  read.table("07142007.dat", header=TRUE)
data07172007    <-  read.table("07172007.dat", header=TRUE)
data07242007    <-  read.table("07242007.dat", header=TRUE)
data07312007    <-  read.table("07312007.dat", header=TRUE)
data07102007Occ <-  read.table("07102007occ.dat", header=TRUE)

data07102007    <- cbind(data07102007, data07102007Occ[,"Lane1Occ"])
names(data07102007)[7]<- "Lane1Occ"
data07102007[,"kmSpeed"]<- v*data07102007Occ$Lane1Speed

## Figure 1  (Fundamental diagram)
par(mfrow=c(1,2))
plot(data07102007$Lane1Flow,data07102007$kmSpeed, col=1, xlab="Flow (veh/5 min)", ylab="Average Speed (km/h)")  # q(v)
plot(data07102007$Lane1Occ,data07102007$kmSpeed, col=1, xlab= "Occupancy", ylab="Average Speed (km/h)")# k(v)

## Figure 2 (LPC and projections)
h <-12
lpc07102007    <- lpc(cbind(data07102007$Lane1Flow, data07102007$kmSpeed), h=h, x0=c(150,50), plotlpc=0) # fits LPC
spline07102007 <- lpc.splinefun(lpc07102007)  # creates a spline function through local centers of mass
fit07102007    <- lpc.fit.spline(spline07102007, add.spline=FALSE) # evaluates spline function 
proj07102007   <- lpc.project.spline(spline07102007, cbind(data07102007$Lane1Flow, data07102007$kmSpeed), num.knots=100, optimize=TRUE, add.proj=FALSE) # computes projections

par(mfrow=c(2,2))
plot(data07102007$Lane1Flow,data07102007$kmSpeed, xlab="Flow (veh/5 min)", ylab="Average Speed (km/h)", col="dark grey")
points(lpc07102007$LPC[,1],lpc07102007$LPC[,2], pch="+")
lines(lpc07102007$LPC[,1],lpc07102007$LPC[,2])
plot(data07102007$Lane1Flow, data07102007$kmSpeed, xlab="Flow (veh/5 min)", ylab="Average Speed (km/h)", col="dark grey")
lines(lpc07102007$LPC[,1], lpc07102007$LPC[,2], lwd=1.5, lty=1)
for (i in 1:nrow(data07102007)){
        lines(rbind( c(data07102007[i,"Lane1Flow"], data07102007[i,"kmSpeed"]), c(proj07102007$closest.coords[i,1],proj07102007$closest.coords[i,2])),lty=2,col=4)
      }



## Figure 3 (Calibration curves)
par(mfrow=c(1,2))
ratio07102007 <-  12*lpc07102007$LPC[,1]/lpc07102007$LPC[,2]
plot(lpc07102007$Par[,1],ratio07102007,type="l",xlab="t", ylab="q(t)/v(t)") # flow and density re-scaled
plot(proj07102007$closest.pi, data07102007$Lane1Occ,type="p",xlab="t", ylab="occ(t)", col="dark grey") # occupancy
sm07102007 <- sm.spline(proj07102007$closest.pi,data07102007$Lane1Occ)
lines(sm07102007)

## Prediction example in Section 2.2
cal <- cbind(lpc07102007$Par[,1],ratio07102007)
cal[c(98,99),]
             ratio07102007
[1,] 35.6269      41.92800
[2,] 44.8267      38.86881

# for instance, density = 40 veh/km   
(41.92800 -40)/(41.92800- 38.86881 )   # linear interpolation
[1] 0.6302322                         
 35.6269+  0.6302322 *( 44.8267 -35.6269 )
[1] 41.42491                         # t
lpc.spline.eval(spline07102007,  41.42491 )
[1] 137.9369 41.49645                # q, v

# Test: k=q/v=
 (12*137.9369)/41.49645
 39.88878 # quite close to original value of 40!





## Figure 5  (6 consecutive weekdays)

par(mfrow=c(3,2))
h <- 12
lpc07092007   <- lpc(cbind(data07092007$Lane1Flow, v*data07092007$Lane1Speed), h=h, x0=c(150,50), plotlpc=2)
lpc07102007   <- lpc(cbind(data07102007$Lane1Flow, v*data07102007$Lane1Speed), h=h, x0=c(150,50), plotlpc=2)
lpc07112007   <- lpc(cbind(data07112007$Lane1Flow, v*data07112007$Lane1Speed), h=h, x0=c(150,50), plotlpc=2)
lpc07122007   <- lpc(cbind(data07122007$Lane1Flow, v*data07122007$Lane1Speed), h=h, x0=c(150,50), plotlpc=2) 
lpc07132007   <- lpc(cbind(data07132007$Lane1Flow, v*data07132007$Lane1Speed), h=h, x0=c(150,50), plotlpc=2) 
lpc07142007   <- lpc(cbind(data07142007$Lane1Flow, v*data07142007$Lane1Speed), h=h, x0=c(150,130), plotlpc=2) 

hastie07092007 <- principal.curve(cbind(data07092007$Lane1Flow/1.5, v*data07092007$Lane1Speed))  
hastie07102007 <- principal.curve(cbind(data07102007$Lane1Flow/1.5, v*data07102007$Lane1Speed))
hastie07112007 <- principal.curve(cbind(data07112007$Lane1Flow/1.5, v*data07112007$Lane1Speed))
hastie07122007 <- principal.curve(cbind(data07122007$Lane1Flow/1.5, v*data07122007$Lane1Speed))
hastie07132007 <- principal.curve(cbind(data07132007$Lane1Flow/1.5, v*data07132007$Lane1Speed))
hastie07142007 <- principal.curve(cbind(data07142007$Lane1Flow/1.5, v*data07142007$Lane1Speed))  

lmg07092007 <- lm(Lane1Flow~ -1+I(v*Lane1Speed)+ I(v*Lane1Speed^2), data= data07092007)
lmg07102007 <- lm(Lane1Flow~ -1+I(v*Lane1Speed)+ I(v*Lane1Speed^2), data= data07102007)
lmg07112007 <- lm(Lane1Flow~ -1+I(v*Lane1Speed)+ I(v*Lane1Speed^2), data= data07112007)
lmg07122007 <- lm(Lane1Flow~ -1+I(v*Lane1Speed)+ I(v*Lane1Speed^2), data= data07122007)
lmg07132007 <- lm(Lane1Flow~ -1+I(v*Lane1Speed)+ I(v*Lane1Speed^2), data= data07132007)
lmg07142007 <- lm(Lane1Flow~ -1+I(v*Lane1Speed)+ I(v*Lane1Speed^2), data= data07142007)

par(mfrow=c(3,2))
plot(data07092007$Lane1Flow, v*data07092007$Lane1Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (km/h)", col="grey", main="09/07/2007")
lines(lpc07092007$LPC[,1], lpc07092007$LPC[,2] , lwd=1.5, lty=1)
lines(1.5*hastie07092007$s[,1][order(hastie07092007$lambda)],hastie07092007$s[,2][order(hastie07092007$lambda)], lwd=1.5,  col=1, lty=2)
lines(fitted(lmg07092007)[order(data07092007$Lane1Speed)], v*data07092007$Lane1Speed[order(data07092007$Lane1Speed)], lwd=1.5, lty=4)
legend(10,80, c("LPC", "HS", "GM"), lwd=c(1.5,1.5,1.5), lty=c(1,2,4))

plot(data07102007$Lane1Flow, v*data07102007$Lane1Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (km/h)", col="grey", main="10/07/2007")
lines(lpc07102007$LPC[,1],lpc07102007$LPC[,2] , lwd=1.5, lty=1)
lines(1.5*hastie07102007$s[,1][order(hastie07102007$lambda)],hastie07102007$s[,2][order(hastie07102007$lambda)], lwd=1.5,  col=1, lty=2)
lines(fitted(lmg07102007)[order(data07102007$Lane1Speed)], v*data07102007$Lane1Speed[order(data07102007$Lane1Speed)], lwd=1.5, lty=4)

plot(data07112007$Lane1Flow, v*data07112007$Lane1Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (km/h)", col="grey", main="11/07/2007")
lines(lpc07112007$LPC[,1],lpc07112007$LPC[,2] , lwd=1.5, lty=1)
lines(1.5*hastie07112007$s[,1][order(hastie07112007$lambda)],hastie07112007$s[,2][order(hastie07112007$lambda)], lwd=1.5,  col=1, lty=2)
lines(fitted(lmg07112007)[order(data07112007$Lane1Speed)], v*data07112007$Lane1Speed[order(data07112007$Lane1Speed)], lwd=1.5, lty=4)

plot(data07122007$Lane1Flow, v*data07122007$Lane1Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (km/h)", col="grey", main="12/07/2007")
lines(lpc07122007$LPC[,1],lpc07122007$LPC[,2] , lwd=1.5, lty=1)
lines(1.5*hastie07122007$s[,1][order(hastie07122007$lambda)],hastie07122007$s[,2][order(hastie07122007$lambda)], lwd=1.5,  col=1, lty=2)
lines(fitted(lmg07122007)[order(data07122007$Lane1Speed)], v*data07122007$Lane1Speed[order(data07122007$Lane1Speed)], lwd=1.5, lty=4)

plot(data07132007$Lane1Flow, v*data07132007$Lane1Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (km/h)", col="grey", main="13/07/2007")
lines(lpc07132007$LPC[,1],lpc07132007$LPC[,2] , lwd=1.5, lty=1)
lines(1.5*hastie07132007$s[,1][order(hastie07132007$lambda)],hastie07132007$s[,2][order(hastie07132007$lambda)], lwd=1.5,  col=1, lty=2)
lines(fitted(lmg07132007)[order(data07132007$Lane1Speed)], v*data07132007$Lane1Speed[order(data07132007$Lane1Speed)], lwd=1.5, lty=4)

plot(data07142007$Lane1Flow, v*data07142007$Lane1Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (km/h)", col="grey", main="14/07/2007", ylim=c(20,130))
lines(lpc07142007$LPC[,1],lpc07142007$LPC[,2] , lwd=1.5, lty=1)
lines(1.5*hastie07142007$s[,1][order(hastie07142007$lambda)],hastie07142007$s[,2][order(hastie07142007$lambda)], lwd=1.5,  col=1, lty=2)
lines(fitted(lmg07142007)[order(data07142007$Lane1Speed)], v*data07142007$Lane1Speed[order(data07142007$Lane1Speed)], lwd=1.5, lty=4)


## Fig. 6 (six consecutive Tuesdays)

par(mfrow=c(3,2))
h<-12
lpc06262007   <- lpc(cbind(data06262007$Lane1Flow, v*data06262007$Lane1Speed), h=h, x0=c(150,50), plotlpc=2)
lpc07032007   <- lpc(cbind(data07032007$Lane1Flow, v*data07032007$Lane1Speed), h=h, x0=c(150,50), plotlpc=2)
lpc07102007   <- lpc(cbind(data07102007$Lane1Flow, v*data07102007$Lane1Speed), h=h, x0=c(150,50), plotlpc=2)
lpc07172007   <- lpc(cbind(data07172007$Lane1Flow, v*data07172007$Lane1Speed), h=h, x0=c(150,50), plotlpc=2)
lpc07242007   <- lpc(cbind(data07242007$Lane1Flow, v*data07242007$Lane1Speed), h=h, x0=c(150,50), plotlpc=2) 
lpc07312007   <- lpc(cbind(data07312007$Lane1Flow, v*data07312007$Lane1Speed), h=h, x0=c(150,50), plotlpc=2) 

hastie06262007 <- principal.curve(cbind(data06262007$Lane1Flow/1.6, v*data07092007$Lane1Speed))  
hastie07032007 <- principal.curve(cbind(data07032007$Lane1Flow/1.5, v*data07032007$Lane1Speed))
hastie07102007 <- principal.curve(cbind(data07102007$Lane1Flow/1.5, v*data07102007$Lane1Speed))
hastie07172007 <- principal.curve(cbind(data07172007$Lane1Flow/1.5, v*data07172007$Lane1Speed))
hastie07242007 <- principal.curve(cbind(data07242007$Lane1Flow/1.5, v*data07242007$Lane1Speed))
hastie07312007 <- principal.curve(cbind(data07312007$Lane1Flow/1.5, v*data07312007$Lane1Speed))  

lmg06262007 <- lm(Lane1Flow~ -1+I(v*Lane1Speed)+ I(v*Lane1Speed^2), data= data06262007)
lmg07032007 <- lm(Lane1Flow~ -1+I(v*Lane1Speed)+ I(v*Lane1Speed^2), data= data07032007)
lmg07102007 <- lm(Lane1Flow~ -1+I(v*Lane1Speed)+ I(v*Lane1Speed^2), data= data07102007)
lmg07172007 <- lm(Lane1Flow~ -1+I(v*Lane1Speed)+ I(v*Lane1Speed^2), data= data07172007)
lmg07242007 <- lm(Lane1Flow~ -1+I(v*Lane1Speed)+ I(v*Lane1Speed^2), data= data07242007)
lmg07312007 <- lm(Lane1Flow~ -1+I(v*Lane1Speed)+ I(v*Lane1Speed^2), data= data07312007)

par(mfrow=c(3,2))
plot(data06262007$Lane1Flow, v*data06262007$Lane1Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (mph)", col="grey", main ="26/06/2007")
lines(lpc06262007$LPC[,1],lpc06262007$LPC[,2] , lwd=1.5, lty=1)
lines(1.5*hastie06262007$s[,1][order(hastie06262007$lambda)],hastie06262007$s[,2][order(hastie06262007$lambda)], lwd=1.5,  col=1, lty=2)
lines(fitted(lmg06262007)[order(data06262007$Lane1Speed)], v*data06262007$Lane1Speed[order(data06262007$Lane1Speed)], lwd=1.5, lty=4)
legend(10,80, c("LPC", "HS", "GM"), lwd=c(1.5,1.5,1.5), lty=c(1,2,4))

plot(data07032007$Lane1Flow, v*data07032007$Lane1Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (mph)", col="grey", main="03/07/2007")
lines(lpc07032007$LPC[,1],lpc07032007$LPC[,2] , lwd=1.5, lty=1)
lines(1.5*hastie07032007$s[,1][order(hastie07032007$lambda)],hastie07032007$s[,2][order(hastie07032007$lambda)], lwd=1.5,  col=1, lty=2)
lines(fitted(lmg07032007)[order(data07032007$Lane1Speed)], v*data07032007$Lane1Speed[order(data07032007$Lane1Speed)], lwd=1.5, lty=4)

plot(data07102007$Lane1Flow, v*data07102007$Lane1Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (mph)", col="grey", main="10/07/2007")
lines(lpc07102007$LPC[,1],lpc07102007$LPC[,2] , lwd=1.5, lty=1)
lines(1.5*hastie07102007$s[,1][order(hastie07102007$lambda)],hastie07102007$s[,2][order(hastie07102007$lambda)], lwd=1.5,  col=1, lty=2)
lines(fitted(lmg07102007)[order(data07102007$Lane1Speed)], v*data07102007$Lane1Speed[order(data07102007$Lane1Speed)], lwd=1.5, lty=4)
lines(lpc07102007right$LPC[300:398,1],lpc07102007right$LPC[300:398,2] , lwd=0.75, lty=1)

plot(data07172007$Lane1Flow, v*data07172007$Lane1Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (mph)", col="grey", main="17/07/2007")
lines(lpc07172007$LPC[,1],lpc07172007$LPC[,2] , lwd=1.5, lty=1)
lines(1.5*hastie07172007$s[,1][order(hastie07172007$lambda)],hastie07172007$s[,2][order(hastie07172007$lambda)], lwd=1.5,  col=1, lty=2)
lines(fitted(lmg07172007)[order(data07172007$Lane1Speed)], v*data07172007$Lane1Speed[order(data07172007$Lane1Speed)], lwd=1.5, lty=4)

plot(data07242007$Lane1Flow, v*data07242007$Lane1Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (mph)", col="grey", main="24/07/2007")
lines(lpc07242007$LPC[,1],lpc07242007$LPC[,2] , lwd=1.5, lty=1)
lines(1.5*hastie07242007$s[,1][order(hastie07242007$lambda)],hastie07242007$s[,2][order(hastie07242007$lambda)], lwd=1.5,  col=1, lty=2)
lines(fitted(lmg07242007)[order(data07242007$Lane1Speed)], v*data07242007$Lane1Speed[order(data07242007$Lane1Speed)], lwd=1.5, lty=4)

plot(data07312007$Lane1Flow, v*data07312007$Lane1Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (mph)", col="grey", main="31/07/2007")
lines(lpc07312007$LPC[,1],lpc07312007$LPC[,2] , lwd=1.5, lty=1)
lines(1.5*hastie07312007$s[,1][order(hastie07312007$lambda)],hastie07312007$s[,2][order(hastie07312007$lambda)], lwd=1.5,  col=1, lty=2)
lines(fitted(lmg07312007)[order(data07312007$Lane1Speed)], v*data07312007$Lane1Speed[order(data07312007$Lane1Speed)], lwd=1.5, lty=4)



# Additional branch in Fig. 6 (1/07/2007)
lpc07102007right  <- lpc(cbind(data07102007$Lane1Flow, v*data07102007$Lane1Speed), h=12, x0=c(150,50, 150,120),  mult=2, plotlpc=2, rho=0.4)
plot(data07102007$Lane1Flow, data07102007$kmSpeed, xlab="Flow (veh/5 min)", ylab="Average Speed (km/h)", col="grey", main="10/07/2007")
lines(lpc07102007right$LPC[1:199,1],lpc07102007right$LPC[1:199,2] , lwd=1.5, lty=1)
lines(lpc07102007right$LPC[300:398,1],lpc07102007right$LPC[300:398,2] , lwd=1.5, lty=2)


## Figure 7 (consecutive weekdays and Tuesdays in each one plot)
par(mfrow=c(1,2))
plot(lpc07092007$LPC, type="l", xlim=c(0,200), ylim=c(0,130), xlab="Flow", ylab="Speed")
lines(lpc07102007$LPC, col=1)
lines(lpc07112007$LPC, col=1)
lines(lpc07122007$LPC, col=1)
lines(lpc07132007$LPC, col=1)
lines(lpc07142007$LPC, col=1)

plot(lpc07312007$LPC, type="l", xlim=c(0,200), ylim=c(0,130), xlab="Flow", ylab= "Speed")
lines(lpc06262007$LPC, col=1)
lines(lpc07032007$LPC, col=1)
lines(lpc07102007$LPC, col=1)
lines(lpc07172007$LPC, col=1)
lines(lpc07242007$LPC, col=1)

##  Table 1 (Capacities)
capacity1 <- matrix(0,6,3)
capacity1[1,1] <- max(lpc06262007$LPC[,1])
capacity1[2,1] <- max(lpc07032007$LPC[,1])
capacity1[3,1] <- max(lpc07102007$LPC[,1])
capacity1[4,1] <- max(lpc07172007$LPC[,1])
capacity1[5,1] <- max(lpc07242007$LPC[,1])
capacity1[6,1] <- max(lpc07312007$LPC[,1])

capacity1[1,2] <- max(1.5*hastie06262007$s[,1])
capacity1[2,2] <- max(1.5*hastie07032007$s[,1])
capacity1[3,2] <- max(1.5*hastie07102007$s[,1])
capacity1[4,2] <- max(1.5*hastie07172007$s[,1])
capacity1[5,2] <- max(1.5*hastie07242007$s[,1])
capacity1[6,2] <- max(1.5*hastie07312007$s[,1])

capacity1[1,3] <- predict(lmg06262007, newdata=data.frame(Lane1Speed= -lmg06262007$coef[1]/(2*lmg06262007$coef[2])))
capacity1[2,3] <- predict(lmg07032007, newdata=data.frame(Lane1Speed= -lmg07032007$coef[1]/(2*lmg07032007$coef[2])))
capacity1[3,3] <- predict(lmg07102007, newdata=data.frame(Lane1Speed= -lmg07102007$coef[1]/(2*lmg07102007$coef[2])))
capacity1[4,3] <- predict(lmg07172007, newdata=data.frame(Lane1Speed= -lmg07172007$coef[1]/(2*lmg07242007$coef[2])))
capacity1[5,3] <- predict(lmg07242007, newdata=data.frame(Lane1Speed= -lmg07242007$coef[1]/(2*lmg07242007$coef[2])))
capacity1[6,3] <- predict(lmg07312007, newdata=data.frame(Lane1Speed= -lmg07312007$coef[1]/(2*lmg07312007$coef[2])))

capacity2 <- capacity1*12  # rescale to veh/hour
apply(capacity2, 2,mean)
[1] 1993.839 2048.888 2122.837
apply(capacity2, 2,sd)
[1] 46.18938 67.40080 52.86674




## Fig. 8 (Selected representative patterns)
par(mfrow=c(3,2))

# Two regimes
h <-12
data07112007VDS1202053<- read.table("07112007VDS1202053.dat", header=TRUE)
lpc07112007VDS1202053 <- lpc(cbind(data07112007VDS1202053$Lane1Flow, v*data07112007VDS1202053$Lane1Speed), 12,plotlpc=0,xlab="Flow (veh/5 min)", ylab="Average Speed (mph)", x0=c(30,120,100,20), mult=2)
lpc07112007VDS1202053a <- lpc(cbind(data07112007VDS1202053$Lane1Flow, v*data07112007VDS1202053$Lane1Speed), 6,plotlpc=0,xlab="Flow (veh/5 min)", ylab="Average Speed (mph)", x0=c(30,120,100,20), mult=2)

plot(data07112007VDS1202053$Lane1Flow, v*data07112007VDS1202053$Lane1Speed,  xlab="Flow (veh/5 min)", ylab="Average Speed (km/h)", col="grey", main="(a) 11/07/2007, VDS 1202053")
lines(lpc07112007VDS1202053a$LPC[1:199,1], lpc07112007VDS1202053a$LPC[1:199,2], lwd=1.5, lty=1)
lines(lpc07112007VDS1202053a$LPC[200:398,1],lpc07112007VDS1202053a$LPC[200:398,2], lwd=1.5, lty=1)
lines(lpc07112007VDS1202053$LPC[300:398,1], lpc07112007VDS1202053$LPC[300:398,2], lwd=1, lty=2)

# Night time drivers
data07102007VDS1202263 <- read.table("07102007VDS1202263.dat", header=TRUE)
lpc19<-lpc(cbind(data07102007VDS1202263$Lane5Flow, v*data07102007VDS1202263$Lane5Speed),12, plotlpc=0,xlab="Flow (veh/5 min)", ylab="Average Speed (mph)", x0=c(120,50))

plot(data07102007VDS1202263$Lane5Flow, v*data07102007VDS1202263$Lane5Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (km/h)", col="grey", main="(b) 10/07/2007, VDS 1202263")
lines(lpc19$LPC[,1],lpc19$LPC[,2], lwd=1.5, lty=1)

# Undertaking
data29 <-   read.table("data29", header=TRUE)
lpc29<-lpc(cbind(data29$Lane5Flow, v*data29$Lane5Speed),6, x0=c(30,20,35,50),mult=2, plotlpc=0,xlab="Flow (veh/5 min)", ylab="Average Speed (mph)")

plot(data29$Lane5Flow, v*data29$Lane5Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (km/h)", col="grey", main="(c) 16/07/2007, VDS 401513")
lines(lpc29$LPC[1:199,1],lpc29$LPC[1:199,2], lwd=1.5, lty=1)
lines(lpc29$LPC[300:398,1],lpc29$LPC[300:398,2], lwd=1, lty=1)

# Adverse weather
data12142006  <-  read.table("12142006VDS717672.dat", header=TRUE)
lpc12142006   <- lpc(cbind(data12142006$Lane1Flow, v*data12142006$Lane1Speed), h=11, x0=c(60,80,30,130,160,50), mult=3, plotlpc=0)

plot(data12142006$Lane1Flow, v*data12142006$Lane1Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (km/h)", col="grey", main="(d) 14/12/2006, VDS 717672")
lines(lpc12142006$LPC[1:199,1], lpc12142006$LPC[1:199,2], lwd=1.5, lty=1)
lines(lpc12142006$LPC[200:398,1],lpc12142006$LPC[200:398,2], lwd=1.5, lty=1)
lines(lpc12142006$LPC[399:597,1],lpc12142006$LPC[399:597,2], lwd=1.5, lty=1)

# Slow lanes
data25 <-   read.table("data25", header=TRUE)
lpc25 <-lpc(cbind(data25$Lane5Flow, v*data25$Lane5Speed),12,  plotlpc=0,xlab="Flow (veh/5 min)", ylab="Average Speed (mph)", x0=c(30,20))
plot(data25$Lane5Flow, v*data25$Lane5Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (km/h)", col="grey", main="(e) 09/07/2007, VDS 1213624")
lines(lpc25$LPC[,1],lpc25$LPC[,2], lwd=1.5, lty=1)

# Faulty
data24 <-   read.table("data24", header=TRUE)
lpc24<-lpc(cbind(data24$Lane1Flow, v*data24$Lane1Speed),12, plotlpc=0, x0=c(0, 140, 120,80), mult=2,  xlab="Flow (veh/5 min)", ylab="Average Speed (mph)")

plot(data24$Lane1Flow, v*data24$Lane1Speed, xlab="Flow (veh/5 min)", ylab="Average Speed (km/h)", col="grey", main="(f) 09/07/2007, VDS 1202322")
lines(lpc24$LPC[1:199,1], lpc24$LPC[1:199,2], lwd=1.5, lty=1)
lines(lpc24$LPC[200:398,1],lpc24$LPC[200:398,2], lwd=1.5, lty=1)