# April 22, 1999 # May 6, 1999 (in Dallas) setdate <- function(month,day,year,hour,minute,ampm) { dys _ c(31,28,31,30,31,30,31,31,30,31,30,31) mon _ c(0, cumsum(dys[-12])) # days in year before xx/01 names(mon) _ c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec") apm <- c(0,12); names(apm) <- c("am","pm") day _ 365*(as.integer(year)-1998) + mon[month] + ( -1 + as.integer(day) ) hrs _ ( as.integer(hour) %% 12 ) + as.integer(minute)/60 + apm[ampm] ans <- day + hrs/24; names(ans) <- NULL ans } # Estimate demand parameters over time for NY1-->OK run. # fgrep NY1 data.txt | fgrep OK > hivol.dat (then edit out OK-->NY1 rows) # replace ":" with space "PM" with " pm" "AM" with " am" truck <- matrix(scan("hivol.txt",what=character()),ncol=21,byrow=T) # 747 x 21 date.pkup _ setdate(truck[,4],truck[,5],truck[,6],truck[,7],truck[,8],truck[,9]) date.wind _ setdate(truck[,10],truck[,11],truck[,12],truck[,13],truck[,14],truck[,15]) date.book _ setdate(truck[,16],truck[,17],truck[,18],truck[,19],truck[,20],truck[,21]) bt1 <- function() { brks_seq(115,450,1); col_0 # brks_seq(115,450,7); col_1 # brks_seq(115,450,14); col_1 # brks_seq(115,470,28); col_1 hist(date.book,breaks=brks,col=col); abline(v=seq(115,450,7),col=6) hist(date.wind,breaks=brks,col=col); abline(v=seq(115,450,7),col=6) hist(date.pkup,breaks=brks,col=col); abline(v=seq(115,450,7),col=6) hist(date.wind-date.book,40); abline(v=1,col=6) hist(date.pkup-date.book,40); abline(v=1,col=6) hist(date.pkup-date.wind,40); abline(v=1,col=6) hist(date.wind-date.book,40,xlim=c(-5,30),cex=1.0); abline(v=1,col=6) hist(date.pkup-date.book,40,xlim=c(-5,30),cex=1.0); abline(v=1,col=6) hist(date.pkup-date.wind,40,xlim=c(-5,30),cex=1.0); abline(v=1,col=6) diff21 _ date.wind-date.book diff31 _ date.pkup-date.book plot(diff21,diff31,pch="o") } bt2 <- function() { par(mfrow=c(3,1)) brks_seq(115,450,1); col_0 # brks_seq(115,450,7); col_1 # brks_seq(115,450,14); col_1 # brks_seq(115,470,28); col_1 diff21 _ date.wind-date.book; print(range(diff21)) diff31 _ date.pkup-date.book; print(range(diff31)) diff32 _ date.pkup-date.wind; print(range(diff32)) hist(date.wind-date.book,40); abline(v=1,col=6) hist(date.pkup-date.book,40); abline(v=1,col=6) hist(date.pkup-date.wind,40); abline(v=1,col=6) hist(date.wind-date.book,40,xlim=c(-5,30),cex=1.0); abline(v=1,col=6) hist(date.pkup-date.book,40,xlim=c(-5,30),cex=1.0); abline(v=1,col=6) hist(date.pkup-date.wind,40,xlim=c(-5,30),cex=1.0); abline(v=1,col=6) nd <- c("Thu","Fri","Sat","Sun","Mon","Tue","Wed") hist( (date.book)%%7 ,breaks=0:7,main="booking date"); axis(1,0:7) text(seq(.5,6.5),rep(mean(par("usr")[3:4]),7),nd) hist( (date.wind)%%7 ,breaks=0:7,main="pkup window date"); axis(1,0:7) text(seq(.5,6.5),rep(mean(par("usr")[3:4]),7),nd) hist( (date.pkup)%%7 ,breaks=0:7,main="pickup date"); axis(1,0:7) text(seq(.5,6.5),rep(mean(par("usr")[3:4]),7),nd) bk2 <- seq(0,7,.5) hist( (date.book)%%7 ,breaks=bk2,main="booking date"); axis(1,0:7) text(seq(.5,6.5),rep(mean(par("usr")[3:4]),7),nd) hist( (date.wind)%%7 ,breaks=bk2,main="pkup window date"); axis(1,0:7) text(seq(.5,6.5),rep(mean(par("usr")[3:4]),7),nd) hist( (date.pkup)%%7 ,breaks=bk2,main="pickup date"); axis(1,0:7) text(seq(.5,6.5),rep(mean(par("usr")[3:4]),7),nd) bk3 <- seq(0,7,.25) hist( (date.book)%%7 ,breaks=bk3,main="booking date"); axis(1,0:7) text(seq(.5,6.5),rep(mean(par("usr")[3:4]),7),nd) hist( (date.wind)%%7 ,breaks=bk3,main="pkup window date"); axis(1,0:7) text(seq(.5,6.5),rep(mean(par("usr")[3:4]),7),nd) hist( (date.pkup)%%7 ,breaks=bk3,main="pickup date"); axis(1,0:7) text(seq(.5,6.5),rep(mean(par("usr")[3:4]),7),nd) } bt3 <- function(dow=F) { # dow=day of week breakdown mybt <- function(x,hd1,hd2=" ") { # x=date hd1,hd2 = header characters tms <- as.character(c(12,1:12,1:12)); tcks <- seq(0,by=1/24,l=25) brks <- seq(-.0001,1,by=1/24); hdr <- paste(hd1,hd2,sep=" ") # print("in mybt"); browser() hist( x-trunc(x),breaks=brks,main=hdr) axis(1,tcks,tms,cex=.8*par("cex"),line=-2,ticks=F) if(hd2 == " ") { axis(3,tcks,tms,cex=.8*par("cex")) } } n <- length(date.book); nd <- c("Thu","Fri","Sat","Sun","Mon","Tue","Wed") if(dow) {par(mfrow=c(4,2))} else {par(mfrow=c(1,1))} mybt(date.book,"booking hour") if(dow) { for(i in 0:6) { j <- seq(n)[(trunc(date.book)%%7) == i] mybt(date.book[j],"booking hour",nd[i+1]) } } mybt(date.wind,"pickup window hour") if(dow) { for(i in 0:6) { j <- seq(n)[(trunc(date.wind)%%7) == i] mybt(date.wind[j],"pickup window hour",nd[i+1]) } } mybt(date.pkup,"pickup hour") if(dow) { for(i in 0:6) { j <- seq(n)[(trunc(date.pkup)%%7) == i] mybt(date.pkup[j],"pickup hour",nd[i+1]) } } } my.est <- function() { # Sunday Jan 4, 1998 at midnight is t=3.0 (or 12th time period) ## distribution of demand at t=0 from previous (uncalled) periods (2 weeks) p.b <- 1 + trunc(4*date.book) # integer book period ( 1 = jan 1 98 12-6 a.m. ) p.p <- 1 + trunc(4*date.pkup) # integer pkup period ( 1 = jan 1 98 12-6 a.m. ) pw.b <- 1 + trunc((4*(date.book-3))%%28) # integer book period of week (1--28) # 1 = Sunday 12-6 a.m. 28 = Sat 6-12 p.m. mij <- matrix(0,28,56) # two week window (56 periods) nskip <- 0 for(k in seq(date.pkup)) { i <- pw.b[k]; j <- p.p[k]-p.b[k]+1 # careful with time difference [0,55]+1 if(j>56) {nskip <- nskip+1} else {mij[i,j] <- mij[i,j] + 1} } mij <- mij / diff(range(p.p)) # average number of trucks per 6-hr time period print(paste(" number of runs skipped over two weeks = ",nskip)) ### probability of pickup at t=0 from previously scheduled ("window") p.w <- 1 + trunc(4*date.wind) # integer wind period ( 1 = jan 1 98 12-6 a.m. ) p.p <- 1 + trunc(4*date.pkup) # integer pkup period ( 1 = jan 1 98 12-6 a.m. ) pw.w <- 1 + trunc((4*(date.wind-3))%%28) # integer wind period of week (1--28) # 1 = Sunday 12-6 a.m. 28 = Sat 6-12 p.m. pij <- matrix(0,28,24) # six day window (24 periods) nskip <- 0 for(k in seq(date.pkup)) { i <- pw.w[k]; j <- max(1,p.p[k]-p.w[k]+1) # careful with time difference if(j>24) {nskip <- nskip+1} else {pij[i,j] <- pij[i,j] + 1} } for(i in 1:28) { tot <- sum(pij[i,]); pij[i,] <- pij[i,]/max(tot,1) } # probability truck picked up in 6-hr time period relative to window print(paste(" number of runs skipped over two weeks = ",nskip)) assign("mij",mij,where=1) assign("pij",pij,where=1) } # predict mean/variance of demand in future period [l1,l2] # pt0 is the period number (1--28) of the bin after t=0 # l1 ( >= 0) is the beginning of the forecast period # l2 ( >= l1) is the end of the forecast period requested # (default is "next" 12 hours [0,1]) # nt contains counts of "open" window requests # nto contains index in nt pointing to (beginning of) current 6-hr period # mij is 28x56 matrix i = 1...28 day period of week (bin #) of booking # j = 0...55 number of periods ahead that pickup occurs # assumes no seasonal trends (pickup occurs within 2 weeks) # pij is 28x24 matrix i = 1...28 day period of week (bin #) of pickup "window" # j = 0...23 number of periods ahead that pickup occurs # (pickup occurs within 6 days) my.pred <- function(pt0,l1=0,l2=1,nt,nto,mij,pij,bro=F) { period <- function(p,pt) { # compute time period number (1--28) p intervals after # t=0 time period which has number pt 1 + (p+(pt-1))%%28 # shift pt to [0,27] then add, mod, and increment } if( pt0<1 | pt0>28 | length(nt)<1 ) stop(" pt0 nt wrong --- check input") if( l1<0 | l2length(nt) ) stop(" l1 l2 nt nto wrong --- check input") m <- 0; v <- 0 for (k in l2:0) { # estimate demand after t=0 to be picked up in [l1,l2] kp <- period(k,pt0) # integer from 1-28 of bin "k" k1 <- max(l1,k); k2 <- l2 # portion of [l1,l2] beyond bin k j1 <- 1+k1-k; j2 <- 1+k2-k # indices to mij array # booking originates in bin k picked up in l1:l2 j <- j1:j2; kj <- j[j>=1 & j<=56] if(length(kj)>0) {mv <- sum(mij[kp,kj]) m <- m + mv # i.e. l-k = # time intervals ahead v <- v + mv } if(bro) {browser()} } for (k in 0:(23+l2-l1)) { # estimate pickups from "open" "windowed" schedule # go back at least 23 periods over [l1,l2] periods # starting at t = l2, back to l1-23 inlk <- l2-k+nto # index in vector nt of "pickup window" of period l2-k if(inlk<1 | inlk>length(nt)) {nlk <- 0} # assume none open else {nlk <- nt[inlk]} # number with "pickup window" in period l2-k if( nlk==0 ) next # nothing in this window period kp <- period(l2-k,pt0) # day period of window l2-k pden <- 1; if(l2-k<0) {pden <- 1-sum(pij[kp,1:(k-l2)])} # a window "open" before # t=0 but not picked up yet -- conditional prob pk <- 0 # conditional probability will be picked up if(pden>1.e-10) { j1 <- max(l1,l2-k); kj <- (j1:l2)+(k+1-l2) pk <- sum(pij[kp,kj]) / pden } # in periods l1--l2 (k periods ahead) # if pden == 0 this is really an error? m <- m + nlk*pk v <- v + nlk*pk*(1-pk) if(bro) {browser()} } list(m=m,v=v,l12=c(l1,l2)) } my.test <- function(book,wind,pkup,st=150,en=425,by=1.0,len=1.0) { # make 1-day predictions # st = day (Jan 1, 1998 is "zero") # en = last day # by = increment in prediction step (24-hour periods) # len = length of prediction interval ( 24 hour periods) k_sum(pkup<=wind);if(k>0){print("test data have some pickups before window N=");print(k)} no.p <- rep(0,(en-st+1)/by) # predictions (mean) no.v <- no.p # prediction variance no.t <- no.p # true counts l1 <- 0; l2 <- floor(4*len)-1 # prediction interval six-hour [0,3] is default for(i in seq(no.p)) { if(i%%10==0) cat(".") ti <- st + (i-1)*by # beginning time of prediction interval no.t[i] <- length(pkup[ pkup>=ti & pkup=wind ]) j <- seq(wind)[wind < ti+len & pkup>ti] # booked but not picked up if( length(j)==0 ) { nt <- 0; nto <- 1 } # nothing booked else { tj <- floor( 4*(wind[j]-ti) ) # 6-hour periods j1 <- min(0,min(tj)); j2 <- max(0,max(tj)) nt <- rep(0, j2-j1+1); nto <- 1-j1 # current time bin starts at nt[nto] for(k in seq(tj)) { kk <- tj[k]+nto; nt[kk] <- nt[kk] + 1 } } # nt vector now completed per.ti<- 1 + (floor(4*(4+ti))%%28) # integer from 1-28 which six-hour period # specific to 1/1/98 as time origin !!! ans <- my.pred(per.ti,l1,l2,nt,nto,mij,pij) no.p[i] <- ans$m; no.v[i] <- ans$v } browser() list(pred=no.p, var=no.v, true=no.t) } # 0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;