.Rhistory

#' to only keep distance skims for those O-D pairs that are in the combined trips.
#' The diagonal value of the distance matrix is set to 0 to avoid making that
#' cell receive any values during balancing
dist <- read.csv("C:/Projects/Province-wide/Empties/od_road_distances.csv",
stringsAsFactors = FALSE) %>%
transform(., Flag = paste0(Orig, Dest))
#' only keep those o-D pairs that match the all_obs truck trip table
dist1 <- subset(dist, dist$Flag %in% all_obs$Flag)
#' wide format and set NAs to zero
dist1_cast <- dcast(dist1, Orig ~ Dest, value.var = "Dist")
rownames(dist1_cast) <- dist1_cast[, 1]
dist1_cast <- dist1_cast[, 2:ncol(dist1_cast)]
dist1_cast[is.na(dist1_cast)] <- 0
#' set diagonal to global value of 10
#diag(dist1_cast) <- 0
####################### COST MATRIX ###########################################
#' generate a cost matrix using distance between CSDs and set the beta to -0.01
#' Rolf uses a very small beta of -0.001. This will have a tendency to require
#' longer to calibrate and also it does not create enough variation in the
#' value of the various interchanges, especially those that are far apart.
cost <- dist1_cast
initbeta <- -0.019
################ COST FUNCTION MATRIX #########################################
#' Cost Function Matrix. We are using an exponential function here. Other formulations
#' can be tested as well, if desired. All those cells that had a value of 0
#' would automatically get set to 1 after taking the exponent. This would
#' incorrectly make these interchanges a candidate to receive trips thereby
#' making it difficult to calibrate and converge the betas. These cells are
#' forced to zero.
cfunc <- exp(initbeta*cost)
cfunc[cfunc == 1] <- 0
#' Estimate observed average trip lengths
tlfd_obs <- sum(dist1_cast * all_obs_cast)/totals
print(tlfd_obs)
############################ FURNESSING ##############################
#' This is the furnessing function, which is the backbone of the calibration
#' procedure.
furness <- function(cmat){
"
This function conducts furnessing or two dimensional balancing
Inputs: observed truck matrix
Arguments: cost function
Results: balanced matrix
"
#' row sum of cost function matrix
skim_trow <- as.data.frame(t(colSums(cmat)))
#' row sum of observed matrix
obs_tatt <- as.data.frame(t(colSums(all_obs_cast)))
#' row balancing factors
rowbal <- obs_tatt/skim_trow
rowbal1 <- do.call("rbind", replicate(nrow(all_obs_cast), rowbal,
simplify = FALSE))
#first iteration
cfunc1 <- cmat*rowbal1
#' col sum of cost function matrix
skim_tcol <- as.data.frame(rowSums(cfunc1))
#' col sum of observed matrix
obs_tprod <- as.data.frame(rowSums(all_obs_cast))
#' column balancing factors
colbal <- obs_tprod/skim_tcol
colbal1 <- do.call("cbind", replicate(ncol(all_obs_cast), colbal,
simplify = FALSE))
# next iteration
cfunc1 <- cfunc1*colbal1
return(cfunc1)
}
############################### CALIBRATION PROCESS ###########################
#' the outer loop is to calculate revised beta values. The user can set the
#' number of outer loops one wants to evaluate.
outer_iter <- 20
#' Start outer loop
for (j in 1:outer_iter) {
################## BETA COEFFICIENTS ########################################
#' set initial beta coefficient and cost matrix
if(j == 1){
beta <- initbeta
print(beta)
#' compute cfunc again
cfunc <- exp(initbeta*cost)
cfunc[cfunc == 1] <- 0
} else {
beta <- beta1
print(beta1)
#' compute cfunc again
cfunc <- exp(beta*cost)
cfunc[cfunc == 1] <- 0
}
#' run the loop with user defined iterations
maxiter <- 15
i <- 1
#' run the first iteration of furnessing
for (i in 1:1){
t <- furness(cfunc)
}
#' now loop furnessing through till the maxiterations
for (i in 2:maxiter){
t1 <- furness(t)
t <- t1
#print(i)
}
###############################################################################
#' Estimate simulated average trip lengths
tlfd_sim <- sum(dist1_cast * t1)/totals
#' calculate revised beta value to feed back to outer loop
if (j == 1){
beta1 <- initbeta * (tlfd_sim/tlfd_obs)
} else {
beta1 <- beta1 * (tlfd_sim/tlfd_obs)
}
}
write.csv(t, "/Projects/Province-wide/Empties/t.csv")
knitr::opts_chunk$set(echo = TRUE)
library(dplyr); library(knitr); library(grid); library(ggplot2); library(gridExtra)
clean <- function(truck, distdf){
"
Function to bring the observed truck matrix and O-D distance flows
together and generate a distance matrix that can be used to display
a trip length frequncy distribution and the cost function
Arguments: observed truck matrix, and distance by O-D
Inputs: same as argument
Retur: distance matrix
"
all_obs <- transform(truck, Flag = paste0(Orig,Dest)) %>%
subset(., Orig != Dest)
dist <- transform(distdf, Flag = paste0(Orig, Dest))
#' only keep those o-D pairs that match the all_obs truck trip table
dist1 <- subset(dist, dist$Flag %in% all_obs$Flag)
#' add in bands
dist1 <- merge(dist1, all_obs, by = "Flag") %>%
arrange(., Dist) %>%
transform(., DFlag = ifelse(Dist <= 50, 1,
ifelse(Dist > 50 & Dist <= 100, 2,
ifelse(Dist > 200 & Dist <= 500, 4,
ifelse(Dist > 500 & Dist <= 1000, 5,
ifelse(Dist > 1000 & Dist <= 2000, 6,
ifelse(Dist > 100 & Dist <= 200, 3, 7)))))))
return(dist1)
}
#' read observed truck trip matrix and combine
obs_hvy_cb <- read.csv("C:/Projects/Province-wide/Empties/empty_heavy_truck_matrix_int.csv",
stringsAsFactors = FALSE)
obs_med_cb <- read.csv("C:/Projects/Province-wide/Empties/empty_straight_truck_matrix_int.csv",
stringsAsFactors = FALSE)
obs_hvy <- read.csv("C:/Projects/Province-wide/Empties/empty_heavy_truck_matrix.csv",
stringsAsFactors = FALSE)
obs_med <- read.csv("C:/Projects/Province-wide/Empties/empty_straight_truck_matrix.csv",
stringsAsFactors = FALSE)
dist_cb <- read.csv("C:/Projects/Province-wide/Empties/od_road_distances_all.csv",
stringsAsFactors = FALSE)
dist <- read.csv("C:/Projects/Province-wide/Empties/od_road_distances.csv",
stringsAsFactors = FALSE)
dist1_cb_hvy <- clean(obs_hvy_cb, dist_cb)
dist1_cb_med <- clean(obs_med_cb, dist_cb)
dist1_hvy <- clean(obs_hvy, dist)
dist1_med <- clean(obs_med, dist)
#' Domestic Flows
oh <- ggplot(dist1_hvy, aes(DFlag, Trips)) + geom_bar(stat = "identity", aes(fill = factor(DFlag))) +
ggtitle("TLFD - Domestic Heavies") + scale_fill_discrete(name="Distance Bands",
labels=c("<50km", "50-100", "100-200", "200-500", "500-1000", "1000-2000", ">2000")) +
theme(plot.title = element_text(size=10)) + theme(legend.title = element_text(size = 8)) +
theme(legend.text = element_text(size = 6))
print("The average distance travelled by a heavy truck for domestic flows is 137.7 kms")
om <- ggplot(dist1_med, aes(DFlag, Trips)) + geom_bar(stat = "identity", aes(fill = factor(DFlag))) +
ggtitle("TLFD - Domestic Mediums") + scale_fill_discrete(name="Distance Bands",
labels=c("<50km", "50-100", "100-200", "200-500", "500-1000", "1000-2000", ">2000")) +
theme(plot.title = element_text(size=10))  + theme(legend.title = element_text(size = 8)) +
theme(legend.text = element_text(size = 6))
print("The average distance travelled by a medium truck for domestic flows is 77.3 kms")
#' International Flows
ch <- ggplot(dist1_cb_hvy, aes(DFlag, Trips)) + geom_bar(stat = "identity", aes(fill = factor(DFlag))) +
ggtitle("TLFD - International Heavies") + scale_fill_discrete(name="Distance Bands",
labels=c("<50km", "50-100", "100-200", "200-500", "500-1000", "1000-2000", ">2000")) +
theme(plot.title = element_text(size=10))  + theme(legend.title = element_text(size = 8)) +
theme(legend.text = element_text(size = 6))
print("The average distance travelled by a heavy truck for international flows is 403.9 kms")
cm <- ggplot(dist1_cb_med, aes(DFlag, Trips)) + geom_bar(stat = "identity", aes(fill = factor(DFlag))) +
ggtitle("TLFD - International Mediums") + scale_fill_discrete(name="Distance Bands",
labels=c("<50km", "50-100", "100-200", "200-500", "500-1000", "1000-2000", ">2000")) +
theme(plot.title = element_text(size=10))  + theme(legend.title = element_text(size = 8)) +
theme(legend.text = element_text(size = 6))
print("The average distance travelled by a medium truck for international flows is 195.9 kms")
grid.arrange(oh, om, ch, cm)
#' Domestic Cost FUnctions
d_h <- -0.00332
d_m <- -0.009442
dist1_hvy <- transform(dist1_hvy, Cf = exp(d_h*dist1_hvy$Dist))
dist1_med <- transform(dist1_med, Cf = exp(d_m*dist1_med$Dist))
#' International Const FUnctions
i_h <- -0.006311
i_m <- -0.015388
dist1_cb_hvy <- transform(dist1_cb_hvy, Cf = exp(i_h*dist1_cb_hvy$Dist))
dist1_cb_med <- transform(dist1_cb_med, Cf = exp(i_m*dist1_cb_med$Dist))
p <- ggplot() +
geom_line(data = dist1_hvy, aes(x = Dist, y = Cf), color="red", linetype = 2)+geom_point(color="red") +
geom_line(data = dist1_med, aes(x = Dist, y = Cf), color="red", linetype = 1)+geom_point(color="red") +
geom_line(data = dist1_cb_hvy, aes(x = Dist, y = Cf), color="blue", linetype = 2)+geom_point(color="blue") +
geom_line(data = dist1_cb_med, aes(x = Dist, y = Cf), color="blue", linetype = 1)+geom_point(color="blue") +
ggtitle("Cost Functions for Domestic and International Flows categorized by Medium and Heavy Trucks")
p
library(dplyr)
library(reshape2)
library(tidyr)
library(foreign)
###############################################################################
#' Title: "Calibrating a Doubly Constrained Gravity Model - Ontario Only"
#' Programmer: "Mausam Duggal, Systems Analysis Group, WSP|PB"
#' Date: "January 11, 2017"
###############################################################################
#' This code was written to calibrate the beta value in a doubly constrained
#' for the empty truck model. The data used was the 2012 CVS data. Bryce
#' prepared the empty truck matrices for both heavies and medium trucks.
#' The trucks have been currently combined because the commodity flow model is
#' not currently distinguishing between heavy and medium truck tonnage. Further,
#' while Rolf's Empty Truck model can be developed for different truck types it
#' was not deemed necessary in this version of the model.
###############################################################################
########################### INPUTS ############################################
#' read observed truck trip matrix and combine
obs_hvy <- read.csv("C:/Projects/Province-wide/Empties/empty_heavy_truck_matrix.csv",
stringsAsFactors = FALSE)
#' Combine all empties and also cast it in a format desired and getting rid of
#' intraCSD trips.
all_obs <- transform(obs_hvy, Flag = paste0(Orig,Dest)) %>%
subset(., Orig != Dest)
#' wide format and set NAs to zero
all_obs_cast <- dcast(all_obs, Orig ~ Dest, value.var = "Trips")
rownames(all_obs_cast) <- all_obs_cast[, 1]
all_obs_cast <- all_obs_cast[, 2:ncol(all_obs_cast)]
all_obs_cast[is.na(all_obs_cast)] <- 0
#' get total trips
totals <- sum(apply(all_obs_cast, 1, function(x) sum(x)))
library(dplyr)
library(reshape2)
library(tidyr)
library(foreign)
###############################################################################
#' Title: "Calibrating a Doubly Constrained Gravity Model - Ontario Only"
#' Programmer: "Mausam Duggal, Systems Analysis Group, WSP|PB"
#' Date: "January 11, 2017"
###############################################################################
#' This code was written to calibrate the beta value in a doubly constrained
#' for the empty truck model. The data used was the 2012 CVS data. Bryce
#' prepared the empty truck matrices for both heavies and medium trucks.
#' The trucks have been currently combined because the commodity flow model is
#' not currently distinguishing between heavy and medium truck tonnage. Further,
#' while Rolf's Empty Truck model can be developed for different truck types it
#' was not deemed necessary in this version of the model.
###############################################################################
########################### INPUTS ############################################
#' read observed truck trip matrix and combine
obs_hvy <- read.csv("C:/Projects/Province-wide/Empties/empty_heavy_truck_matrix.csv",
stringsAsFactors = FALSE)
#' Combine all empties and also cast it in a format desired and getting rid of
#' intraCSD trips.
all_obs <- transform(obs_hvy, Flag = paste0(Orig,Dest)) %>%
subset(., Orig != Dest)
#' wide format and set NAs to zero
all_obs_cast <- dcast(all_obs, Orig ~ Dest, value.var = "Trips")
rownames(all_obs_cast) <- all_obs_cast[, 1]
all_obs_cast <- all_obs_cast[, 2:ncol(all_obs_cast)]
all_obs_cast[is.na(all_obs_cast)] <- 0
#' get total trips
totals <- sum(apply(all_obs_cast, 1, function(x) sum(x)))
View(all_obs_cast)
write.csv(all_obs_cast, "C:/Projects/Province-wide/Empties/dome_hv.csv", row.names = FALSE)
dist <- read.csv("C:/Projects/Province-wide/Empties/od_road_distances.csv",
stringsAsFactors = FALSE) %>%
transform(., Flag = paste0(Orig, Dest))
#' only keep those o-D pairs that match the all_obs truck trip table
dist1 <- subset(dist, dist$Flag %in% all_obs$Flag)
#' wide format and set NAs to zero
dist1_cast <- dcast(dist1, Orig ~ Dest, value.var = "Dist")
rownames(dist1_cast) <- dist1_cast[, 1]
dist1_cast <- dist1_cast[, 2:ncol(dist1_cast)]
dist1_cast[is.na(dist1_cast)] <- 0
library(dplyr)
library(reshape2)
library(tidyr)
library(foreign)
###############################################################################
#' Title: "Calibrating a Doubly Constrained Gravity Model - International"
#' Programmer: "Mausam Duggal, Systems Analysis Group, WSP|PB"
#' Date: "January 11, 2017"
###############################################################################
#' This code was written to calibrate the beta value in a doubly constrained
#' for the empty truck model. The data used was the 2012 CVS data. Bryce
#' prepared the empty truck matrices for both heavies and medium trucks.
#' The trucks have been currently combined because the commodity flow model is
#' not currently distinguishing between heavy and medium truck tonnage. Further,
#' while Rolf's Empty Truck model can be developed for different truck types it
#' was not deemed necessary in this version of the model.
###############################################################################
########################### INPUTS ############################################
#' read observed truck trip matrix and combine
obs_hvy <- read.csv("C:/Projects/Province-wide/Empties/empty_heavy_truck_matrix_int.csv",
stringsAsFactors = FALSE)
#' Combine all empties and also cast it in a format desired and getting rid of
#' intraCSD trips.
all_obs <- transform(obs_hvy, Flag = paste0(Orig,Dest)) %>%
subset(., Orig != Dest)
#' wide format and set NAs to zero
all_obs_cast <- dcast(all_obs, Orig ~ Dest, value.var = "Trips")
rownames(all_obs_cast) <- all_obs_cast[, 1]
all_obs_cast <- all_obs_cast[, 2:ncol(all_obs_cast)]
all_obs_cast[is.na(all_obs_cast)] <- 0
#' get total trips
totals <- sum(apply(all_obs_cast, 1, function(x) sum(x)))
#' batch in distance skims and cast it in a format desired. There is also a need
#' to only keep distance skims for those O-D pairs that are in the combined trips.
#' The diagonal value of the distance matrix is set to 0 to avoid making that
#' cell receive any values during balancing
dist <- read.csv("C:/Projects/Province-wide/Empties/od_road_distances_all.csv",
stringsAsFactors = FALSE) %>%
transform(., Flag = paste0(Orig, Dest))
#' only keep those o-D pairs that match the all_obs truck trip table
dist1 <- subset(dist, dist$Flag %in% all_obs$Flag)
#' wide format and set NAs to zero
dist1_cast <- dcast(dist1, Orig ~ Dest, value.var = "Dist")
rownames(dist1_cast) <- dist1_cast[, 1]
dist1_cast <- dist1_cast[, 2:ncol(dist1_cast)]
dist1_cast[is.na(dist1_cast)] <- 0
#' set diagonal to global value of 10
#diag(dist1_cast) <- 0
####################### COST MATRIX ###########################################
#' generate a cost matrix using distance between CSDs and set the beta to -0.01
#' Rolf uses a very small beta of -0.001. This will have a tendency to require
#' longer to calibrate and also it does not create enough variation in the
#' value of the various interchanges, especially those that are far apart.
cost <- dist1_cast
initbeta <- -0.005
################ COST FUNCTION MATRIX #########################################
#' Cost Function Matrix. We are using an exponential function here. Other formulations
#' can be tested as well, if desired. All those cells that had a value of 0
#' would automatically get set to 1 after taking the exponent. This would
#' incorrectly make these interchanges a candidate to receive trips thereby
#' making it difficult to calibrate and converge the betas. These cells are
#' forced to zero.
cfunc <- exp(initbeta*cost)
cfunc[cfunc == 1] <- 0
#' Estimate observed average trip lengths
tlfd_obs <- sum(dist1_cast * all_obs_cast)/totals
print(tlfd_obs)
############################ FURNESSING ##############################
#' This is the furnessing function, which is the backbone of the calibration
#' procedure.
furness <- function(cmat){
"
This function conducts furnessing or two dimensional balancing
Inputs: observed truck matrix
Arguments: cost function
Results: balanced matrix
"
#' row sum of cost function matrix
skim_trow <- as.data.frame(t(colSums(cmat)))
#' row sum of observed matrix
obs_tatt <- as.data.frame(t(colSums(all_obs_cast)))
#' row balancing factors
rowbal <- obs_tatt/skim_trow
rowbal1 <- do.call("rbind", replicate(nrow(all_obs_cast), rowbal,
simplify = FALSE))
#first iteration
cfunc1 <- cmat*rowbal1
#' col sum of cost function matrix
skim_tcol <- as.data.frame(rowSums(cfunc1))
#' col sum of observed matrix
obs_tprod <- as.data.frame(rowSums(all_obs_cast))
#' column balancing factors
colbal <- obs_tprod/skim_tcol
colbal1 <- do.call("cbind", replicate(ncol(all_obs_cast), colbal,
simplify = FALSE))
# next iteration
cfunc1 <- cfunc1*colbal1
return(cfunc1)
}
############################### CALIBRATION PROCESS ###########################
#' the outer loop is to calculate revised beta values. The user can set the
#' number of outer loops one wants to evaluate.
outer_iter <- 30
#' Start outer loop
for (j in 1:outer_iter) {
################## BETA COEFFICIENTS ########################################
#' set initial beta coefficient and cost matrix
if(j == 1){
beta <- initbeta
print(beta)
#' compute cfunc again
cfunc <- exp(initbeta*cost)
cfunc[cfunc == 1] <- 0
} else {
beta <- beta1
print(beta1)
#' compute cfunc again
cfunc <- exp(beta*cost)
cfunc[cfunc == 1] <- 0
}
#' run the loop with user defined iterations
maxiter <- 15
i <- 1
#' run the first iteration of furnessing
for (i in 1:1){
t <- furness(cfunc)
}
#' now loop furnessing through till the maxiterations
for (i in 2:maxiter){
t1 <- furness(t)
t <- t1
#print(i)
}
###############################################################################
#' Estimate simulated average trip lengths
tlfd_sim <- sum(dist1_cast * t1)/totals
#' calculate revised beta value to feed back to outer loop
if (j == 1){
beta1 <- initbeta * (tlfd_sim/tlfd_obs)
} else {
beta1 <- beta1 * (tlfd_sim/tlfd_obs)
}
}
write.csv(t, "/Projects/Province-wide/Empties/t.csv")
library(dplyr); library(knitr); library(grid); library(ggplot2); library(gridExtra)
clean <- function(truck, distdf){
"
Function to bring the observed truck matrix and O-D distance flows
together and generate a distance matrix that can be used to display
a trip length frequncy distribution and the cost function
Arguments: observed truck matrix, and distance by O-D
Inputs: same as argument
Retur: distance matrix
"
all_obs <- transform(truck, Flag = paste0(Orig,Dest)) %>%
subset(., Orig != Dest)
dist <- transform(distdf, Flag = paste0(Orig, Dest))
#' only keep those o-D pairs that match the all_obs truck trip table
dist1 <- subset(dist, dist$Flag %in% all_obs$Flag)
#' add in bands
dist1 <- merge(dist1, all_obs, by = "Flag") %>%
arrange(., Dist) %>%
transform(., DFlag = ifelse(Dist <= 50, 1,
ifelse(Dist > 50 & Dist <= 100, 2,
ifelse(Dist > 200 & Dist <= 500, 4,
ifelse(Dist > 500 & Dist <= 1000, 5,
ifelse(Dist > 1000 & Dist <= 2000, 6,
ifelse(Dist > 100 & Dist <= 200, 3, 7)))))))
return(dist1)
}
#' read observed truck trip matrix and combine
obs_hvy_cb <- read.csv("C:/Projects/Province-wide/Empties/empty_heavy_truck_matrix_int.csv",
stringsAsFactors = FALSE)
obs_med_cb <- read.csv("C:/Projects/Province-wide/Empties/empty_straight_truck_matrix_int.csv",
stringsAsFactors = FALSE)
obs_hvy <- read.csv("C:/Projects/Province-wide/Empties/empty_heavy_truck_matrix.csv",
stringsAsFactors = FALSE)
obs_med <- read.csv("C:/Projects/Province-wide/Empties/empty_straight_truck_matrix.csv",
stringsAsFactors = FALSE)
dist_cb <- read.csv("C:/Projects/Province-wide/Empties/od_road_distances_all.csv",
stringsAsFactors = FALSE)
dist <- read.csv("C:/Projects/Province-wide/Empties/od_road_distances.csv",
stringsAsFactors = FALSE)
dist1_cb_hvy <- clean(obs_hvy_cb, dist_cb)
dist1_cb_med <- clean(obs_med_cb, dist_cb)
dist1_hvy <- clean(obs_hvy, dist)
dist1_med <- clean(obs_med, dist)
#' Domestic Flows
oh <- ggplot(dist1_hvy, aes(DFlag, Trips)) + geom_bar(stat = "identity", aes(fill = factor(DFlag))) +
ggtitle("TLFD - Domestic Heavies") + scale_fill_discrete(name="Distance Bands",
labels=c("<50km", "50-100", "100-200", "200-500", "500-1000", "1000-2000", ">2000")) +
theme(plot.title = element_text(size=10)) + theme(legend.title = element_text(size = 8)) +
theme(legend.text = element_text(size = 6))
print("The average distance travelled by a heavy truck for domestic flows is 137.7 kms")
om <- ggplot(dist1_med, aes(DFlag, Trips)) + geom_bar(stat = "identity", aes(fill = factor(DFlag))) +
ggtitle("TLFD - Domestic Mediums") + scale_fill_discrete(name="Distance Bands",
labels=c("<50km", "50-100", "100-200", "200-500", "500-1000", "1000-2000", ">2000")) +
theme(plot.title = element_text(size=10))  + theme(legend.title = element_text(size = 8)) +
theme(legend.text = element_text(size = 6))
print("The average distance travelled by a medium truck for domestic flows is 77.3 kms")
#' International Flows
ch <- ggplot(dist1_cb_hvy, aes(DFlag, Trips)) + geom_bar(stat = "identity", aes(fill = factor(DFlag))) +
ggtitle("TLFD - International Heavies") + scale_fill_discrete(name="Distance Bands",
labels=c("<50km", "50-100", "100-200", "200-500", "500-1000", "1000-2000", ">2000")) +
theme(plot.title = element_text(size=10))  + theme(legend.title = element_text(size = 8)) +
theme(legend.text = element_text(size = 6))
print("The average distance travelled by a heavy truck for international flows is 403.9 kms")
cm <- ggplot(dist1_cb_med, aes(DFlag, Trips)) + geom_bar(stat = "identity", aes(fill = factor(DFlag))) +
ggtitle("TLFD - International Mediums") + scale_fill_discrete(name="Distance Bands",
labels=c("<50km", "50-100", "100-200", "200-500", "500-1000", "1000-2000", ">2000")) +
theme(plot.title = element_text(size=10))  + theme(legend.title = element_text(size = 8)) +
theme(legend.text = element_text(size = 6))
print("The average distance travelled by a medium truck for international flows is 195.9 kms")
grid.arrange(oh, om, ch, cm)
#' Domestic Cost FUnctions
d_h <- -0.00332
d_m <- -0.009442
dist1_hvy <- transform(dist1_hvy, Cf = exp(d_h*dist1_hvy$Dist))
dist1_med <- transform(dist1_med, Cf = exp(d_m*dist1_med$Dist))
#' International Const FUnctions
i_h <- -0.006311
i_m <- -0.015388
dist1_cb_hvy <- transform(dist1_cb_hvy, Cf = exp(i_h*dist1_cb_hvy$Dist))
dist1_cb_med <- transform(dist1_cb_med, Cf = exp(i_m*dist1_cb_med$Dist))
p <- ggplot() +
geom_line(data = dist1_hvy, aes(x = Dist, y = Cf), color="red", linetype = 2)+geom_point(color="red") +
geom_line(data = dist1_med, aes(x = Dist, y = Cf), color="red", linetype = 1)+geom_point(color="red") +
geom_line(data = dist1_cb_hvy, aes(x = Dist, y = Cf), color="blue", linetype = 2)+geom_point(color="blue") +
geom_line(data = dist1_cb_med, aes(x = Dist, y = Cf), color="blue", linetype = 1)+geom_point(color="blue") +
ggtitle("Cost Functions for Domestic and International Flows categorized by Medium and Heavy Trucks")
p
readme.md
readme.md
setwd("C:/Projects/R Scripts/R-Codes")
readme.md
README.md
Sys.which("git")