[Travelr-commits] r10 - in pre-alpha: . travelr travelr/R travelr/man travelr/src travelr/tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Apr 30 03:10:30 CEST 2010
Author: jrawbits
Date: 2010-04-30 03:10:30 +0200 (Fri, 30 Apr 2010)
New Revision: 10
Added:
pre-alpha/travelr/man/assignment.class.Rd
pre-alpha/travelr/man/bpr.function.Rd
pre-alpha/travelr/man/free.flow.Rd
pre-alpha/travelr/man/skim.paths.Rd
pre-alpha/travelr/tests/Test_Cost_Integrator.R
pre-alpha/travelr_0.1-2.zip
Removed:
pre-alpha/travelr/plan.txt
pre-alpha/travelr_0.1-1.zip
Modified:
pre-alpha/ReadMe.txt
pre-alpha/travelr/DESCRIPTION
pre-alpha/travelr/NAMESPACE
pre-alpha/travelr/R/Assignment_Set.R
pre-alpha/travelr/R/Highway_Assignment.R
pre-alpha/travelr/R/Highway_Network.R
pre-alpha/travelr/R/Iterative_Fitting.R
pre-alpha/travelr/R/Path_Functions.R
pre-alpha/travelr/man/assignment.functions.Rd
pre-alpha/travelr/man/assignment.set.Rd
pre-alpha/travelr/man/build.intercepts.Rd
pre-alpha/travelr/man/build.network.set.Rd
pre-alpha/travelr/man/build.paths.Rd
pre-alpha/travelr/man/cost.functions.Rd
pre-alpha/travelr/man/highway.assign.Rd
pre-alpha/travelr/man/highway.net.Rd
pre-alpha/travelr/man/node.map.Rd
pre-alpha/travelr/man/shortest.paths.Rd
pre-alpha/travelr/man/utilities.Rd
pre-alpha/travelr/src/build_path_R.h
pre-alpha/travelr/src/path_templates_R.h
pre-alpha/travelr/tests/Test_02_Assignment_Classes.R
pre-alpha/travelr/tests/Test_03_Highway_Assignment.R
pre-alpha/travelr/tests/run.tests.bat
Log:
Release 0.1-2
Frank-Wolfe works on the small test network.
ParTan is there, but may not work
The cost.integrator does much better, but could use some additional testing.
Multi-class still needs testing with more than one class.
Intercepts need testing with iterative algorithms and multiple classes
Skims need testing (see sketch in skim.paths documentation)
Need to document and clean up cost functions -- drop the per-class cost function (too hard to manage in practice).
Modified: pre-alpha/ReadMe.txt
===================================================================
--- pre-alpha/ReadMe.txt 2010-04-26 02:15:14 UTC (rev 9)
+++ pre-alpha/ReadMe.txt 2010-04-30 01:10:30 UTC (rev 10)
@@ -7,7 +7,7 @@
for detailed information on what's in that package, where to get it, and how to
install it).
-Development (as of 4/25/2010) has been done with R 2.10.1 (even though 2.11.0 has
+Development (as of 4/29/2010) has been done with R 2.10.1 (even though 2.11.0 has
just been released) and a consistent Rtools package, and the package should
build with those tools.
Modified: pre-alpha/travelr/DESCRIPTION
===================================================================
--- pre-alpha/travelr/DESCRIPTION 2010-04-26 02:15:14 UTC (rev 9)
+++ pre-alpha/travelr/DESCRIPTION 2010-04-30 01:10:30 UTC (rev 10)
@@ -1,6 +1,6 @@
Package: travelr
-Version: 0.1-1
-Date: 2009-12-09
+Version: 0.1-2
+Date: 2009-04-26
Title: Functions for Travel Demand Modeling
Author: Jeremy Raw
Maintainer: Jeremy Raw <jeremy.raw at earthlink.net>
Modified: pre-alpha/travelr/NAMESPACE
===================================================================
--- pre-alpha/travelr/NAMESPACE 2010-04-26 02:15:14 UTC (rev 9)
+++ pre-alpha/travelr/NAMESPACE 2010-04-30 01:10:30 UTC (rev 10)
@@ -18,25 +18,37 @@
S3method(summary,highway.net)
S3method(print,highway.net)
S3method(print,sum.highway.net)
+export(tntp.network)
+export(tntp.od)
# Assignment Set operation functions
-export(free.flow)
export(build.paths)
export(load.paths)
export(build.and.load.paths)
-export(intercept.paths)
-export(skim.paths)
-# Assignment Set constrution functions
-export(cost.integrator)
-export(build.serious.objective.function)
-export(build.easy.objective.function)
+# Assignment Set construction functions
export(new.assignment.set)
export(make.assignment.class)
export(add.assignment.class)
-export(build.BPR.function)
+S3method(add.assignment.class,default)
+S3method(add.assignment.class,list)
+S3method(add.assignment.class,highway.assignment.set)
+export(hwy.update.class)
+export(hwy.update.penalties)
+export(hwy.update.demand)
+
+# Skims and intercepts
+export(skim.paths)
export(new.intercept.set)
+export(intercept.paths)
+# Cost function management
+export(free.flow)
+export(build.BPR.cost.function)
+export(build.BPR.objective.function)
+export(cost.integrator)
+export(build.general.objective.function)
+
# High-level algorithms
export(highway.assign)
@@ -54,6 +66,7 @@
export(build.lambda.search)
# Utility functions
+export(hwy.gamma.function)
export(rmse)
export(ipf)
S3method(print,iterative.fit)
@@ -62,6 +75,7 @@
export(.build.network.set)
export(.shortest.paths)
export(.load.paths)
+export(.build.and.load.paths)
export(.skim.paths)
export(.intercept.paths)
export(.walk.paths)
Modified: pre-alpha/travelr/R/Assignment_Set.R
===================================================================
--- pre-alpha/travelr/R/Assignment_Set.R 2010-04-26 02:15:14 UTC (rev 9)
+++ pre-alpha/travelr/R/Assignment_Set.R 2010-04-30 01:10:30 UTC (rev 10)
@@ -39,8 +39,8 @@
# build.paths will produce a list of path sets from route costs
build.paths<-function(aset,costs)
- mapply(function(ac,cost) .shortest.paths(ac$network,cost),
- aset$classes,costs,
+ mapply(function(ac,cost,penalties) .shortest.paths(ac$network,cost,penalties),
+ aset$classes,costs,MoreArgs=(penalties=aset$penalties),
USE.NAMES=TRUE, # Keep assignment class names on resulting list of path structures
SIMPLIFY=FALSE) # and leave results as a list
@@ -57,13 +57,13 @@
# logic at C-level. It's useful in link-based assignment, but not in path-based assignment
# where one is equilibrating path sets, rather than volume sets
build.and.load.paths<-function(aset,costs) {
- result=mapply( function(ac,cost).build.and.load.paths(ac$network,cost,ac$demand),
- aset$classes,costs,
+ result=mapply( function(ac,cost,penalties).build.and.load.paths(ac$network,cost,ac$demand,penalties),
+ aset$classes,costs,MoreArgs=list(penalties=aset$penalties),
USE.NAMES=TRUE, # Keep assignment class names on resulting volume set list
SIMPLIFY=FALSE) # and leave results as a list (each element of which is
# a list with "paths" and "volumes" elements)
paths = lapply(result,function(p)p$paths) # Paths for all classes in one list
- volumes = lapply(result,function(p)p$volumes) # Volumes for all classes in one list
+ volumes = as.data.frame(lapply(result,function(p)p$volumes)) # Volumes for all classes in one list
return( list(paths=paths, volumes=volumes, result=result) )
}
@@ -74,11 +74,11 @@
# skim paths returns a demand matrix by performing a function on a set of numeric values
# corresponding to the links on each path between origin and destination
-skim.paths<-function( paths, costs, empty.val=0.0, FUN="sum" ) {
+skim.paths<-function( paths, costs, empty.val=0.0, FUN="sum", ... ) {
if ( FUN=="sum" )
return( .skim.paths(paths,costs,empty.val) )
else
- stop("Not implemented: Skims using functions other than 'sum'")
+ stop("Skimming with functions other than 'sum' is not yet implemented")
}
# Assignment Set
@@ -98,38 +98,9 @@
# demand (matrix)
# cost.function (optional if assignment set cost function is defined)
-cost.integrator <- function( cost.function, ff.vol, ff.cost, cong.vol, cong.cost, depth=0, tol=1e-4, max.depth=7, ... ) {
- # Integrate the cost function using adaptive quadrature with inside-interval tolerance testing
- # The tolerance indicates how close the mid-segment interpolated cost must be to its true value in order to avoid
- # recursing into each of the sub-intervals
- # May do as many as 2^depth+1 calls to cost.function
- mid.cost <- (cong.cost-ff.cost)/2
- vol.ivl <- (cong.vol-ff.vol)/2
- mid.vol <- ff.vol+vol.ivl
- true.cost <- cost.function(mid.vol)
- test.diff <- abs(mid.cost-true.cost)/true.cost
- if ( test.diff < tol || depth>=max.depth ) # limit on recursion
- result <- sum(vol.ivl*((ff.cost+true.cost)+(true.cost+cong.cost)))/2 # area for this part
- else {
- depth <- depth+1
- result<- Recall( cost.function, ff.cost, ff.vol, mid.vol, true.cost, tol=tol, depth=depth, ... ) +
- Recall( cost.function, true.cost, mid.vol, cong.vol, cong.cost, tol=tol, depth=depth, ... )
- }
- attr(result,"max.depth")<-depth
- return(result)
-}
-
-# "correct" objective function
-build.serious.objective.function<-function(cost.function, ff.vol, ff.cost, tol=1e-4, ...)
- function(cost, volume) cost.integrator(cost.function, ff.vol, ff.cost, volume, cost, tol, ...)
-
-# expedient objective function
-build.easy.objective.function<-function(cost.function, ff.vol, ff.cost, ...)
- function(cost,volume) sum( cost*volume )
-
# Construct an assignment.set with all the right pieces in all the right places
new.assignment.set <- function( network,classes,cost.volume.type=c("vector","matrix"),cost.function=NULL,
- make.objective.function=NULL, obj.tol=1e-4, ... ) {
+ objective.function=NULL, obj.tol=1e-4, ... ) {
# We need cost functions to handle either a list/data.frame of volumes or a numeric vector
# But we would like to figure out how to determine whether to boil down the volumes to a vector
# or pass the entire data.frame.
@@ -142,7 +113,15 @@
if ( is.null(network) || class(network) != "highway.net" )
stop("Assignment set requires a highway network")
aset<-list(network=network)
+ if ( is.null(network$Penalty.fields) || is.null(network$Penalty[[network$Penalty.fields[["Penalty.value"]]]]) )
+ aset$penalties<-NULL
+ else
+ aset$penalties<-network$Penalty[[network$Penalty.fields[["Penalty.value"]]]] # may be NULL
+ # Check classes
+ if ( any(sapply(classes,function(c) class(c) )!="highway.assignment.class") )
+ stop("Must use make.assignment.class and add.assignment.classes to build the 'classes' list")
+
# Check cost.volume.type and construct a function to transform the volumes appropriately
if ( length(cost.volume.type)>1 ) # Use default: "vector"
cost.volume.type<-cost.volume.type[1]
@@ -173,30 +152,19 @@
}
}
- # Check classes
- if ( !is.list(classes) )
- stop("Assignment set cannot be built from non-list 'classes'")
-
# Validate presence of required class elements
- bad.names <-names(classes)!=sapply(classes,function(x){ ifelse(is.null(x$name),"", x$name)})
- if ( any(bad.names) )
- stop("Assignment set classes do not have consistent names:\n",paste(names(classes)[bad.names],collapse=", "))
- missing.cost.function<-sapply( classes, function(x) { if (is.null(x$cost.function)) TRUE; FALSE } )
- if ( !all(missing.cost.function) && any(missing.cost.function) )
- message("Assignment set classes have incomplete class functions:",paste(names(classes)[missing.cost.function],collapse=", "))
- missing.demand<-sapply( classes, function(x) { if (is.null(x$demand)) TRUE; FALSE } )
- if ( any(missing.demand) )
- stop("Assignment set classes are missing required demand matrix:",paste(names(classes)[missing.demand],collapse=", "))
- missing.network.set<-sapply( classes, function(x) { if (is.null(x$network.set)||class(x$network.set)!="highway.network.set") TRUE; FALSE } )
- if ( any(missing.network.set) )
- stop("Assignment set classes are missing suitable network set:",paste(names(classes)[missing.demand],collapse=", "))
- aset$classes<-classes
+ if ( validate.assignment.classes(aset,classes) )
+ aset$classes<-classes
+ else
+ stop("Invalid assignment class.")
# Construct free-flow volume and cost, and build objective function for assignment
aset$ff.vol<-free.flow(aset)
- aset$ff.cost<-aset$cost.function(aset$ff.vol)
- if ( is.null(make.objective.function) ) make.objective.function<-build.easy.objective.function
- aset$objective.function<-make.objective.function(cost.function,aset$ff.vol, aset$ff.cost, obj.tol,...)
+ aset$ff.cost<-aset$cost.function(aset$ff.vol,aset)
+ if ( is.null(objective.function) )
+ aset$objective.function<-build.general.objective.function( aset, aset$ff.vol, aset$ff.cost, tol=obj.tol)
+ else
+ aset$objective.function<-objective.function
# Mark this as a known structure
class(aset)<-"highway.assignment.set"
@@ -204,8 +172,37 @@
return(aset)
}
+# Convenience function for validating assignment classes
+validate.assignment.classes <- function( aset, classes ) {
+ bad.names <-names(classes)!=sapply(classes,function(x){ ifelse(is.null(x$name),"", x$name)})
+ if ( any(bad.names) ) {
+ warning("Assignment set classes do not have consistent names:\n",paste(names(classes)[bad.names],collapse=", "))
+ return(FALSE)
+ }
+ # note that the following line works, even if aset doesn't have any classes yet (is.null(aset$classes))
+ missing.cost.function<-sapply( c(aset$classes,classes), function(x) { if (is.null(x$cost.function)) TRUE; FALSE } )
+ if ( !all(missing.cost.function) && any(missing.cost.function) ) {
+ message("Assignment set classes have incomplete class functions:",
+ paste(names(classes)[missing.cost.function],collapse=", "))
+ return(FALSE)
+ }
+ missing.demand<-sapply( classes, function(x) { if (is.null(x$demand)) TRUE; FALSE } )
+ if ( any(missing.demand) ) {
+ warning("Assignment set classes are missing required demand matrix:",
+ paste(names(classes)[missing.demand],collapse=", "))
+ return(FALSE)
+ }
+ missing.network.set<-sapply( classes, function(x) { if (is.null(x$network.set)||class(x$network.set)!="highway.network.set") TRUE; FALSE } )
+ if ( any(missing.network.set) ) {
+ warning("Assignment set classes are missing suitable network set:",
+ paste(names(classes)[missing.demand],collapse=", "))
+ return(FALSE)
+ }
+ return(TRUE)
+}
+
# The following convenience functions make it easy to build a set of assignment classes, once
-# you know the parameters for each class.
+# you know the parameters (network, demand, link and penalty subsets, and cost.function) for each class.
# Example:
# class.list <- vector("list")
# aclass <- make.assignment.class(network,"Class.1",demand.1,links.1,penalty.subset.1,cost.function.1)
@@ -218,32 +215,194 @@
aclass<-list( name=name,demand=demand )
aclass$network.set<-.build.network.set(network,link.subset,penalty.subset)
if (!is.null(cost.function)) aclass$cost.function<-cost.function
+ class(aclass)<-"highway.assignment.class"
return(aclass)
}
# Create a named list entry with the same name as the class
# Do not permit a new class to overwrite an existing class with the same name
# There are easy ways to remove the existing name first (e.g. class.list[[replace.name]]<-NULL)
-add.assignment.class <- function( classes, aclass ) {
+add.assignment.class <- function( classes, aclass ) UseMethod("add.assignment.class",classes)
+
+add.assignment.class.default <- function( classes, aclass )
+ stop("Can only add assignment classes to a list of classes or to an assignment.set")
+
+add.assignment.class.highway.assignment.set <- function( classes, aclass ) {
+ aset <- classes
+ if ( class(aclass)!="highway.assignment.class" )
+ stop("Must use make.assignment.class and add.assignment.classes to build an assignment class")
+ if ( validate.assignment.classes(aset,added.class ) ) {
+ if ( all(match(aclass$name,names(aset$classes),nomatch=0)==0) ) {
+ added.class <- list(aclass)
+ names(added.class) <- aclass$name
+ aset$classes<-c(aset$classes,added.class)
+ } else {
+ warning("Cannot overwrite existing class ",aclass$name)
+ }
+ }
+ return(aset)
+}
+
+add.assignment.class.list <- function( classes, aclass ) {
if ( all(match(aclass$name,names(classes),nomatch=0)==0) ) {
classes[[aclass$name]]<-aclass
} else {
- warning("Cannot overwrite existing class",aclass$name)
+ warning("Cannot overwrite existing class ",aclass$name)
}
return(classes)
}
+# Function to update an entire assignment class... Can't imagine why, but here it is!
+hwy.update.class <- function( aset, name, aclass ) {
+ if ( is.null(name) || is.null(class) )
+ stop("Cannot update assignment class unless name and new data are provided")
+ if ( class(aclass) != "highway.assignment.class" )
+ stop("Use make.assignment.class to construct the new class")
+ if ( is.null(aset[[name]]) )
+ message("Adding assignment class named ",name)
+ aset$classes[[name]] <- aclass
+}
+
+# Function to update the demand for an assignment class once it's in the assignment set
+hwy.update.demand <- function( aset, name, demand ) {
+ if ( is.null(aset[[name]]) )
+ warning("No class named ",name)
+ else if (!is.null(demand))
+ aset[[name]]$demand<-demand
+ else
+ warning("No new demand: assignment set was unchanged")
+ return(aset)
+}
+
+# Function to update the penalty values for an assignment class
+hwy.update.penalties <- function( aset, penalties ) {
+ if ( is.null(aset$penalties) )
+ stop("Number of penalties must be defined when network is built; only the values may change.")
+ else if (!is.null(penalties) && length(penalties)==length(aset$penalties))
+ aset$penalties<-penalties
+ else
+ stop(sprintf(gettext("Must provide %d penalties"),length(aset$penalties)))
+ return(aset)
+}
+
# Cost functions
+
# Here is a helper to make a BPR-type function, using some static data
# This particular version wants a vector of volumes, which is suitable
# if we need a vector-based cost function
-build.BPR.function <- function(cost.data) {
+build.BPR.cost.function <- function(cost.data) {
# cost.data is either an environment, a list, or a data.frame
# All-cap elements must be provided in cost.data
# (either as scalars or as vectors that supply a value for each link in the network)
- return( with(cost.data, function(volume){TIME * ( 1 + (ALPHA/(CAPACITY^BETA)) * ( volume^BETA ) )} ) )
+ FACTOR <- with(cost.data,ALPHA/(CAPACITY^BETA))
+ return( with(cost.data, function(volume,...){TIME * ( 1 + FACTOR * ( volume^BETA ) )} ) )
}
+# Objective Functions
+build.BPR.objective.function <- function(cost.data) {
+ # cost.data is either an environment, a list, or a data.frame
+ # All-cap elements must be provided in cost.data
+ # (either as scalars or as vectors that supply a value for each link in the network)
+ BETA.OBJ<-with(cost.data,BETA+1)
+ FACTOR <- with(cost.data,ALPHA/(CAPACITY^BETA.OBJ))
+ with(cost.data, function(volume) { sum(volume * (TIME * ( 1 + FACTOR * ( volume^BETA ) ) ) ) } )
+}
+
+# "correct" objective function when there is no easy analytic form
+# Analytic form is to be preferred since it may be much faster
+build.general.objective.function<-function(aset, ff.vol, ff.cost, tol=1e-4, max.depth=14 ) {
+ cf <- function(volume) aset$cost.function(volume,aset)
+ function(volume) cost.integrator(cf, ff.vol, ff.cost, volume, cf(volume), tol, max.depth)
+}
+
+# From Wikipedia:
+#
+# Here is an implementation of the adaptive Simpson's method in C99 that avoids
+# redundant evaluations of f and quadrature computations, but also suffers from
+# memory runaway. Suggested improvement would be to use a for loop in
+# adaptiveSimpsonsAux:
+#
+# #include <math.h> // include file for fabs and sin
+# #include <stdio.h> // include file for printf
+#
+# //
+# // Recursive auxiliary function for adaptiveSimpsons() function below
+# //
+# double adaptiveSimpsonsAux(double (*f)(double), double a, double b, double epsilon,
+# double S, double fa, double fb, double fc, int bottom) {
+# double c = (a + b)/2, h = b - a;
+# double d = (a + c)/2, e = (c + b)/2;
+# double fd = f(d), fe = f(e);
+# double Sleft = (h/12)*(fa + 4*fd + fc);
+# double Sright = (h/12)*(fc + 4*fe + fb);
+# double S2 = Sleft + Sright;
+# if (bottom <= 0 || fabs(S2 - S) <= 15*epsilon)
+# return S2 + (S2 - S)/15;
+# return adaptiveSimpsonsAux(f, a, c, epsilon/2, Sleft, fa, fc, fd, bottom-1) +
+# adaptiveSimpsonsAux(f, c, b, epsilon/2, Sright, fc, fb, fe, bottom-1);
+# }
+#
+# //
+# // Adaptive Simpson's Rule
+# //
+# double adaptiveSimpsons(double (*f)(double), // ptr to function
+# double a, double b, // interval [a,b]
+# double epsilon, // error tolerance
+# int maxRecursionDepth) { // recursion cap
+# double c = (a + b)/2, h = b - a;
+# double fa = f(a), fb = f(b), fc = f(c);
+# double S = (h/6)*(fa + 4*fc + fb);
+# return adaptiveSimpsonsAux(f, a, b, epsilon, S, fa, fb, fc, maxRecursionDepth);
+# }
+#
+#
+# int main(){
+# double I = adaptiveSimpsons(sin, 0, 1, 0.000000001, 10); // compute integral of sin(x)
+# // from 0 to 1 and store it in
+# // the new variable I
+# printf("I = %lf\n",I); // print the result
+# return 0;
+# }
+#
+# The "runaway memory" problem could potentially be an issue in the R code, but the
+# Wikipedia-suggested solution of unrolling the recursion into a loop is overkill if the cost
+# function is well-behaved. Time and experience will tell if we need to revisit this
+
+# Note: cost integrator can handle the special requirement that the volume and cost parameters be
+# vectors (and that the integral of the whole be the sum of the integrals of each element from "free
+# flow" volume (0) to the corresponding element in "congested"). This function should also handle
+# everything happily even if ff.vol and cong.vol are data.frames or matrices.
+
+# The standard R 'integrate' function can't do any of that...
+
+cost.integrator <- function( cost.function, ff.vol, ff.cost, cong.vol, cong.cost, tol=1e-8, max.depth=14 ) {
+ # Integrate the cost function using adaptive quadrature with inside-interval tolerance testing
+ ivl <- (cong.vol-ff.vol)/2
+ mid.vol <- ff.vol+ivl
+ mid.cost <- cost.function(mid.vol)
+ result <- (ivl/3)*(ff.cost+4*mid.cost+cong.cost)
+ inner.integrator( cost.function, ff.vol, ff.cost, mid.vol, mid.cost, cong.vol, cong.cost, result, max.depth, tol )
+}
+
+inner.integrator <- function( cost.function, ff.vol,ff.cost,mid.vol,mid.cost,cong.vol,cong.cost,result,depth, tol) {
+ ivl <- (cong.vol-ff.vol)/4
+ left.vol <- ff.vol + ivl
+ left.cost <- cost.function(left.vol)
+ right.vol <- cong.vol - ivl
+ right.cost <- cost.function(right.vol)
+ h.6 <- ivl/3
+ result.left <- sum(h.6*(ff.cost+4*left.cost+mid.cost))
+ result.right<- sum(h.6*(mid.cost+4*right.cost+cong.cost))
+ result.2 <- result.left + result.right
+ if ( depth<=0 || (abs(result.2-result)<=(tol*15)) ) { # if splitting doesn't improve result, return best guess
+ result <- result.2 + (result.2-result)/15
+ } else { # otherwise, process each half recursively seeking a better fit
+ result <- Recall(cost.function,ff.vol,ff.cost,left.vol,left.cost,mid.vol,mid.cost,result.left,depth-1,tol/2)+
+ Recall(cost.function,mid.vol,mid.cost,right.vol,right.cost,cong.vol,cong.cost,result.right,depth-1,tol/2)
+ }
+ result
+ }
+
# Intercept Set, describes select link processing
# The intercept set is built at the level of the highway assignment
@@ -264,4 +423,3 @@
attr(iset,"class")<-"intercept.set"
return(iset)
}
-
Modified: pre-alpha/travelr/R/Highway_Assignment.R
===================================================================
--- pre-alpha/travelr/R/Highway_Assignment.R 2010-04-26 02:15:14 UTC (rev 9)
+++ pre-alpha/travelr/R/Highway_Assignment.R 2010-04-30 01:10:30 UTC (rev 10)
@@ -62,19 +62,20 @@
# You'll get the value explicitly set, or the local default in the algorithm, or the global default
# from the control.defaults list
# Note that some algorithms like MSA may provide alternate "global" defaults
-1
+
control.defaults <- list(
intercept=NULL,
- min.rgap = 1e-4, # stop if relative.gap goes below this value
- max.iter = 100, # for MSA, it's overridden to 4 and other stopping conditions are ignored
+ min.relative.gap = 1e-4, # stop if relative.gap goes below this value
+ max.iter = 100, # for MSA, it's overridden to 4 and other stopping conditions are ignored
# set max.iter=0 for open-ended running
- max.elapsed = 3600, # number of seconds - set to zero to run forever
- opt.tol = .Machine$double.eps^0.5, # tolerance for internal line searches
- verbose = 1, # depending on the algorithm, bigger numbers will produce more tracking output
- log = FALSE # if true, return 'log' attribute (a data.frame) with assignment result details
+ max.elapsed = 3600, # number of seconds - set to zero to run forever
+ opt.tol = .Machine$double.eps^0.5, # tolerance for internal line searches
+ verbose = 1, # depending on the algorithm, bigger numbers will produce more tracking output
+ log = FALSE # if true, return 'log' attribute (a data.frame) with assignment result details
)
parse.control<-function( control, control.element, default.value=NULL ) {
+ cat("Parsing",control.element,"\n")
element<-control[[control.element]]
if ( !is.null(element) ) {
result<-element
@@ -94,30 +95,33 @@
# at the 16th Annual International EMME/2 Users’ Group Conference, Albuquerque, NM,
# March 18-20, 2002.
build.equil.stats.function <- function( objective.function, flow=0 ) {
- avg.excess.cost.func <- ifelse( all.equal(flow,0), function(gap)NULL, function(gap)gap/flow )
+
+ if ( isTRUE(all.equal(flow,0)) ) { avg.excess.cost.func <- function(gap)NULL }
+ else { avg.excess.cost.func <- function(gap)gap/flow }
+
start.time=proc.time()["elapsed"]
- run.time <- function() as.numeric(proc.time["elapsed"]-start.time)
+ run.time <- function() as.numeric(proc.time()["elapsed"]-start.time)
# add additional values to the result list when calling this function
# by providing named numeric scalar values in the ... parameter
- function( iter, cost, vol, diff, best.lower.bound=as.numeric(NA), ... ) {
+ function( iter, cost, vol, diff, best.lower.bound=as.numeric(NA), extras=list() ) {
# vol is the vector of Equilibrium Path Volumes
# diff is the vector of Shortest Path Volumes - Equilibrium Path Volumes
# cost is the vector of link costs given current equilibrium path volumes
iter <- iter
elapsed <- run.time()
- objective.value <- objective.function(cost,vol)
+ objective.value <- objective.function(vol)
gap <- -sum(cost*diff)
lower.bound <- objective.value-gap
best.lower.bound <- max( best.lower.bound, lower.bound, na.rm=TRUE )
- relative.gap <- -gap/abs(best.lower.bound)
+ relative.gap <- (objective.value-best.lower.bound)/abs(best.lower.bound)
avg.excess.cost <- avg.excess.cost.func(gap)
# The return form here produces a named vector of numeric values
- return ( c(iter=iter,elapsed=elapsed,
- objective=objective.value,gap=gap,relative.gap=relative.gap,
- lower.bound=lower.bound,best.lower.bound=best.lower.bound,
- avg.excess.cost=avg.excess.cost,...) )
+ return ( c(list(iter=iter,elapsed=elapsed,
+ objective=objective.value,gap=gap,relative.gap=relative.gap,
+ lower.bound=lower.bound,best.lower.bound=best.lower.bound,
+ avg.excess.cost=avg.excess.cost),extras) )
}
}
@@ -126,15 +130,20 @@
# of "min." or "max." with the rest of their name matching the equilibrium statistic names
build.convergence.test<-function( control, test=c("min.relative.gap","max.iter","max.elapsed") ) {
# pre-compute a setup for convergence testing
- test.min.names <- test[which(grep("^min\\.",test))]
- test.max.names <- test[which(grep("^max\\.",test))]
+ test.min.names <- test[grep("^min\\.",test)]
+ test.max.names <- test[grep("^max\\.",test)]
min.names<-gsub("^min\\.","",test.min.names)
max.names<-gsub("^max\\.","",test.max.names)
- test.min<-unlist(control[test.min.names])
- test.max<-unlist(control[test.max.names])
+ test.min <- numeric(0)
+ for ( i in test.min.names ) test.min <- c(test.min,parse.control(control,i))
+ test.max <- numeric(0)
+ for ( i in test.max.names ) test.max <- c(test.max,parse.control(control,i))
+ cat("Convergence Test::\n")
+ cat("Min:",paste(test.min.names,sep=","),"Max:",paste(test.max.names,sep=","),"\n")
+ cat( "Min:",paste(test.min,collapse=","),"Max:",paste(test.max,collapse=","), "\n")
function(results) {
- all( results[min.names] < test.min ) &&
- all( results[max.names] > test.max )
+ any( unlist(results[min.names]) <= test.min ) ||
+ any( unlist(results[max.names]) >= test.max )
}
}
@@ -143,84 +152,195 @@
# Not much error checking on built-in functions, since they are hidden in the namespace
# and should only be called through the highway.assignment driver function, which
# should do more checking
+
+build.log.function <- function(control) {
+ log <- parse.control(control,"log")
+ if ( log ) {
+ log.function<-function(new.results,log,verbose) {
+ if ( is.data.frame(log) ) {
+ nr <- as.data.frame(new.results)
+ cat( paste(format(nr[1,]),collapse=", "), "\n" )
+ return( rbind(log,nr) )
+ } else {
+ cat( paste(names(new.results),collapse=", "),"\n")
+ return( as.data.frame(new.results) )
+ }
+ }
+ } else
+ log.function<-function(new.results,log,verbose) FALSE
+ return(log.function)
+}
+
+# All-Or-Nothing
.highway.assignment.AON <- function(aset,control) {
verbose <- parse.control(control,"verbose")
if ( verbose>0 )
- cat("All-or-Nothing (AON) highway assignment\n")
+ message("All-or-Nothing (AON) highway assignment\n")
assign.results<-build.and.load.paths(aset,aset$ff.cost)
- cat("Assignment results in AON:\n")
intercept<-(build.intercept.function(iset<-parse.control(control,"intercept"),aset))(assign.results$paths)
return(list(aset=aset,costs=aset$ff.cost,
paths=assign.results$paths,volumes=assign.results$volumes,
iset=iset,intercept=intercept))
}
-# Built-In Assignment Methods
+# Multiple Successive Averages
.highway.assignment.MSA <- function(aset,control) {
-
- log <- parse.control(control,"log")
- if ( log ) {
- log.function<-function(new.results,log,...) {
- if ( is.data.frame(log) )
- return( rbind(log,as.data.frame(new.results)) )
- else
- return( new.results )
- }
- } else
- log.function<-function(new.results,log,...) FALSE
-
+ log.function<-build.log.function(control)
verbose <- parse.control(control,"verbose")
if ( verbose>0 )
- cat("Multiple Successive Averages (MSA) highway assignment\n")
+ message("Multiple Successive Averages (MSA) highway assignment\n")
max.iter <- parse.control(control,"max.iter",4)
if (attr(max.iter,"default")) {
- cat(sprintf(gettext("Maximum iterations for MSA assignment set to default (%d)\n"),max.iter))
+ message(sprintf(gettext("Maximum iterations for MSA assignment set to default (%d)\n"),max.iter))
} else if ( verbose>0 ) {
- cat(sprintf("Maximum iterations for MSA assignment set to %d\n"),max.iter)
+ message(sprintf("Maximum iterations for MSA assignment set to %d\n"),max.iter)
}
control$max.iter<-max.iter
# Construct helper values and functions
- flow<-sum(sapply(aset$aclasses,function(ac)sum(ac$demand),SIMPLIFY=TRUE,USE.NAMES=FALSE))
+ flow<-sum(sapply(aset$classes,function(ac)sum(ac$demand),USE.NAMES=FALSE))
build.intercepts<-build.intercept.function(iset<-parse.control(control,"intercept"),aset)
converged<-build.convergence.test(control,c("max.iter"))
equil.stats<-build.equil.stats.function(aset$objective.function, flow)
iter<-1
- vol <- aset$ff.vol
cost <- aset$ff.cost
+ load <- build.and.load.paths(aset,cost)
+ vol <- load$volumes
+ cost <- aset$cost.function(vol,aset)
intercept<-NULL
best.lower.bound<-as.numeric(NA)
while (TRUE) {
load<-build.and.load.paths(aset,cost)
vol.new<-load$volumes
vol.diff <- vol.new-vol
- new.results <- equil.stats(iter,cost,vol,vol.diff,best.lower.bound)
- best.lower.bound<-new.results$best.lower.bound
+# cat(iter,"Sum cost:",sum(aset$cost.function(vol,aset)),"\n")
+# cat(iter,"Sum shortest-path cost:",sum(aset$cost.function(vol.new,aset)),"\n")
+ new.results <- equil.stats(iter,cost,vol,vol.diff,best.lower.bound,extras=list(cost=cost))
+ best.lower.bound<-new.results[["best.lower.bound"]]
iter<-iter+1
iter.weight <- 1/(iter+1)
vol <- weighted.update.diff( iter.weight, vol, vol.diff )
intercept<- weighted.update.intercept( iter.weight, intercept, build.intercepts(load$paths) )
- cost<-aset$cost.function(vol)
+ cost<-aset$cost.function(vol,aset)
- # Still not a happy logging process...
- log <- log.function( new.results, log )
- if ( nrow(log)>0 && verbose>1 )
- print(log,nrow(log)) # TODO: a prettier job on this output
+ log <- log.function( new.results, log, verbose )
if ( converged(new.results) )
break
}
return(list(aset=aset,costs=cost,paths=load$paths,volumes=load$volumes,iset=iset,intercept=intercept,results=new.results,log=log))
}
+# Frank.Wolfe Algorithm
+.highway.assignment.Frank.Wolfe <- function(aset,control) {
+ log.function<-build.log.function(control)
+ verbose <- parse.control(control,"verbose")
+ if ( verbose>0 )
+ message("Frank-Wolfe highway assignment\n")
+
+ # Construct helper values and functions
+ flow<-sum(sapply(aset$classes,function(ac)sum(ac$demand),USE.NAMES=FALSE))
+ build.intercepts<-build.intercept.function(iset<-parse.control(control,"intercept"),aset)
+ converged<-build.convergence.test(control,c("min.relative.gap","max.iter","max.elapsed"))
+ equil.stats<-build.equil.stats.function(aset$objective.function, flow)
+
+ iter<-1
+ cost <- aset$ff.cost
+ load <- build.and.load.paths(aset,cost)
+ vol <- load$volumes
+ intercept<-NULL
+ best.lower.bound<-as.numeric(NA)
+
+ lambda.func <- build.lambda.function(aset)
+ lambda.search <- build.lambda.search( lambda.func )
+
+ while (TRUE) {
+ cost <- aset$cost.function(vol,aset)
+ load<-build.and.load.paths(aset,cost)
+ vol.new<-load$volumes
+ vol.diff <- vol.new-vol
+ lambda <- lambda.search(vol,vol.diff)$lambda
+
+ new.results <- equil.stats(iter,cost,vol,vol.diff,best.lower.bound,extras=list(lambda=lambda))
+ best.lower.bound<-new.results[["best.lower.bound"]]
+
+ vol <- weighted.update.diff( lambda, vol, vol.diff )
+ intercept<- weighted.update.intercept( lambda, intercept, build.intercepts(load$paths) )
+
+ log <- log.function( new.results, log, verbose )
+ if ( converged(new.results) )
+ break
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/travelr -r 10
More information about the Travelr-commits
mailing list