[Travelr-commits] r9 - in pre-alpha: . travelr travelr/R travelr/data travelr/man travelr/src travelr/tests travelr/tests/Network-Source
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Apr 26 04:15:29 CEST 2010
Author: jrawbits
Date: 2010-04-26 04:15:14 +0200 (Mon, 26 Apr 2010)
New Revision: 9
Added:
pre-alpha/travelr/
pre-alpha/travelr/DESCRIPTION
pre-alpha/travelr/NAMESPACE
pre-alpha/travelr/R/
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/R/tntp.R
pre-alpha/travelr/R/zzz.R
pre-alpha/travelr/data/
pre-alpha/travelr/data/SiouxFalls.rda
pre-alpha/travelr/data/datalist
pre-alpha/travelr/man/
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/intercept.set.Rd
pre-alpha/travelr/man/ipf.Rd
pre-alpha/travelr/man/node.map.Rd
pre-alpha/travelr/man/package.Rd
pre-alpha/travelr/man/sample.networks.Rd
pre-alpha/travelr/man/shortest.paths.Rd
pre-alpha/travelr/man/tntp.Rd
pre-alpha/travelr/man/utilities.Rd
pre-alpha/travelr/plan.txt
pre-alpha/travelr/src/
pre-alpha/travelr/src/build_path_R.cc
pre-alpha/travelr/src/build_path_R.h
pre-alpha/travelr/src/build_r_shlib.bat
pre-alpha/travelr/src/elements_R.h
pre-alpha/travelr/src/makevars.win
pre-alpha/travelr/src/path_templates_R.h
pre-alpha/travelr/src/short_path_R.cc
pre-alpha/travelr/tests/
pre-alpha/travelr/tests/Network-Source/
pre-alpha/travelr/tests/Network-Source/Richmond.VA_Meta.txt
pre-alpha/travelr/tests/Network-Source/Richmond.VA_Network.csv
pre-alpha/travelr/tests/Network-Source/Richmond.VA_Nodes.csv
pre-alpha/travelr/tests/Network-Source/Richmond.VA_OD.csv
pre-alpha/travelr/tests/Network-Source/Richmond.VA_Penalties.txt
pre-alpha/travelr/tests/Network-Source/SiouxFalls-Assignment.s
pre-alpha/travelr/tests/Network-Source/SiouxFalls_Cube.net
pre-alpha/travelr/tests/Network-Source/SiouxFalls_CubeLoaded.net
pre-alpha/travelr/tests/Network-Source/SiouxFalls_LoadTarget.dbf
pre-alpha/travelr/tests/Network-Source/SiouxFalls_Meta.txt
pre-alpha/travelr/tests/Network-Source/SiouxFalls_Network.csv
pre-alpha/travelr/tests/Network-Source/SiouxFalls_Nodes.csv
pre-alpha/travelr/tests/Network-Source/SiouxFalls_OD.csv
pre-alpha/travelr/tests/Network-Source/SiouxFalls_OD.mat
pre-alpha/travelr/tests/Network-Source/SiouxFalls_OD_Cube.csv
pre-alpha/travelr/tests/Network-Source/SiouxFalls_PathTrace.csv
pre-alpha/travelr/tests/RTC_Links.csv
pre-alpha/travelr/tests/RTC_Network.RData
pre-alpha/travelr/tests/RTC_Nodes.csv
pre-alpha/travelr/tests/RTC_OD.csv
pre-alpha/travelr/tests/RTC_Turn_Penalties.txt
pre-alpha/travelr/tests/Test_01_Highway_Network.R
pre-alpha/travelr/tests/Test_02_Assignment_Classes.R
pre-alpha/travelr/tests/Test_03_Highway_Assignment.R
pre-alpha/travelr/tests/Test_04_ipf.R
pre-alpha/travelr/tests/data/
pre-alpha/travelr/tests/install.travelr.bat
pre-alpha/travelr/tests/run.tests.bat
pre-alpha/travelr_0.1-1.zip
Log:
The initial pre-alpha release of TravelR. There's a bunch that probably doesn't work (especially, for example, the skim.paths function), the documentation sucks, and I haven't tested the latest round of multi-class assignment code with more than class, nor have I tested the MSA algorithm, nor have I copied the Frank-Wolfe algorithm into this tree and set it up for multi-class and intercepts. Whew! That said, it's amazingly close to done, at least until we look at the performance on a big network and decide that more needs to be moved to C++!
Added: pre-alpha/travelr/DESCRIPTION
===================================================================
--- pre-alpha/travelr/DESCRIPTION (rev 0)
+++ pre-alpha/travelr/DESCRIPTION 2010-04-26 02:15:14 UTC (rev 9)
@@ -0,0 +1,16 @@
+Package: travelr
+Version: 0.1-1
+Date: 2009-12-09
+Title: Functions for Travel Demand Modeling
+Author: Jeremy Raw
+Maintainer: Jeremy Raw <jeremy.raw at earthlink.net>
+Depends: R (>= 2.10.0)
+Suggests: igraph, sp
+Description: Provides a suite of functions for building travel demand models, with C/C++ support.
+ Includes fast routines for common operations that are specific to travel demand modeling including
+ gravity model, iterative proportional fitting (Fratar/Furness), network path building, and highway
+ assignment (several common and contemporary equilibrium assignment algorithms are implemented).
+Classification/JEL: R.41
+License: GPL (>= 2)
+LazyData: yes
+URL: http://www.opentravelmodels.org
Added: pre-alpha/travelr/NAMESPACE
===================================================================
--- pre-alpha/travelr/NAMESPACE (rev 0)
+++ pre-alpha/travelr/NAMESPACE 2010-04-26 02:15:14 UTC (rev 9)
@@ -0,0 +1,67 @@
+# travelr namespace
+
+# DLL load instruction
+useDynLib(travelr,
+ shortest_paths,
+ load_paths,
+ build_and_load_paths,
+ skim_paths,
+ walk_paths,
+ walk_pairs,
+ intercept_paths)
+
+# Highway Network functions
+export(map.highway.nodes)
+S3method(summary,highway.node.map)
+S3method(print,highway.node.map)
+export(as.highway.net)
+S3method(summary,highway.net)
+S3method(print,highway.net)
+S3method(print,sum.highway.net)
+
+# 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)
+export(new.assignment.set)
+export(make.assignment.class)
+export(add.assignment.class)
+export(build.BPR.function)
+export(new.intercept.set)
+
+# High-level algorithms
+export(highway.assign)
+
+# Algorithm Configuration Utilities
+export(parse.control)
+export(build.equil.stats.function)
+export(build.convergence.test)
+export(build.intercept.function)
+
+# Assignment Algorithm Helpers
+export(weighted.update.diff)
+export(weighted.update)
+export(weighted.update.intercept)
+export(build.lambda.function)
+export(build.lambda.search)
+
+# Utility functions
+export(rmse)
+export(ipf)
+S3method(print,iterative.fit)
+
+# Low-level path functions
+export(.build.network.set)
+export(.shortest.paths)
+export(.load.paths)
+export(.skim.paths)
+export(.intercept.paths)
+export(.walk.paths)
Added: pre-alpha/travelr/R/Assignment_Set.R
===================================================================
--- pre-alpha/travelr/R/Assignment_Set.R (rev 0)
+++ pre-alpha/travelr/R/Assignment_Set.R 2010-04-26 02:15:14 UTC (rev 9)
@@ -0,0 +1,267 @@
+####################################################################################################
+# travelr\r\Assignment_Set.R by Jeremy Raw Copyright (C) 2010
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# A copy of the GNU General Public License is available at
+# http://www.r-project.org/Licenses/
+# and included in the R distribution (in directory ‘share/licenses’).
+####################################################################################################
+
+# This file defines the Assignment Set mechanism for handling single and multi-class assignment
+# and presenting an abstract interface for assignment algorithms.
+
+# The assignment.set structures define the network and demand that will be fitted to each other
+# during the assignment process. The basic operations used in the assignment algorithms are
+# the following (implemented as generic functions):
+# 1. cost.function (uses assignment.set cost.function; generates cost.set from volume.set)
+# 2. build.paths (uses assignment.set network.set; generates path.set from cost.set)
+# 3. load.demand (uses assignment.set demands; generates volume.set from path.set)
+
+
+# free.flow will produce a volume set of zero flow for use in computing free-flow costs
+
+free.flow<-function(aset)
+ as.data.frame(
+ lapply(aset$classes,
+ FUN=function(aclass,numLinks) vector("numeric",numLinks),
+ numLinks=aset$network$numLinks)
+ )
+
+# 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,
+ USE.NAMES=TRUE, # Keep assignment class names on resulting list of path structures
+ SIMPLIFY=FALSE) # and leave results as a list
+
+# load.paths will map demand matrices onto path sets
+
+load.paths<-function(aset,paths)
+ as.data.frame(mapply(function(ac,paths) .load.paths(paths,ac$demand),
+ aset$classes, paths,
+ USE.NAMES=TRUE, # Keep assignment class names on resulting volume set list
+ SIMPLIFY=FALSE) # and leave the results as a list
+ )
+
+# build.and.load.paths provides a slight efficiency boost by leaving more of the
+# 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,
+ 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
+ return( list(paths=paths, volumes=volumes, result=result) )
+}
+
+# intercept.paths returns a list of OD matrices ("od") whose paths intercept selected links
+# we get one intercept structure per assignment class
+intercept.paths<-function( paths, links )
+ lapply( paths,FUN=.intercept.paths,links=links )
+
+# 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" ) {
+ if ( FUN=="sum" )
+ return( .skim.paths(paths,costs,empty.val) )
+ else
+ stop("Not implemented: Skims using functions other than 'sum'")
+}
+
+# Assignment Set
+# This factory function builds a suitable object of one of the assignment set types
+# The assignment sets themselves will need usable functions for computing diagnostic functions
+
+# Structure of an assignment set:
+
+# list with these elements:
+# network
+# cost.function (optional; default is to call a cost function for each class)
+# objective.function (optional; default is to do vector inner product of cost set and volume set)
+# ff.vol (free flow volume - zero flow in all classes)
+# ff.cost (cost function evaluated at zero flow in all classes)
+# classes (named list of assignment classes, each with these elements):
+# network.set -- routable subset of the network for assignment
+# 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, ... ) {
+ # 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.
+ # Suggestion: if any class provides a cost function, all of them receive a data.frame
+ # in that case, for classes with no cost function, the assignment set cost function should
+ # have a default process that takes a volume vector and returns a cost vector the unequipped class
+ # If only the assignment set provides a cost function, then the function receives a vector
+
+ # Check network
+ if ( is.null(network) || class(network) != "highway.net" )
+ stop("Assignment set requires a highway network")
+ aset<-list(network=network)
+
+ # 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]
+ if ( ! (cost.volume.type%in%c("vector","matrix") ) )
+ stop("Assignment set: Need vector or matrix cost function type.")
+ if (cost.volume.type=="vector")
+ cost.volume<-function(volume) rowSums(volume) # volume should be a matrix or data.frame
+ else
+ cost.volume<-function(volume) volume
+
+ # Install cost function
+ if ( is.null(cost.function) ) {
+ # if there is no function for the assignment set, there must be one for each assignment class
+ # cost functions for the assignment class get only a single parameter (volume)
+ # cost functions for the assignment set are passed two parameters (volume,aset)
+ if ( any(sapply( classes, function(x) { if (is.null(x$cost.function)) TRUE; FALSE } ) ) )
+ stop("Assignment set is missing required cost.functions")
+ # the master cost function will then apply the individual class functions
+ aset$cost.function<-function(volume,aset) {
+ as.data.frame(lapply(aset$classes,FUN=function(x,VOL) x$cost.function(cost.volume(VOL)),VOL=volume))
+ }
+ } else {
+ # There is great flexibility in how the cost function works
+ # It can call individual cost functions for only certain classes, build a partial matrix of
+ # volumes with only certain classes, etc. (Set the cost.volume.type to "matrix" in that case)
+ aset$cost.function<-function(volume,aset) {
+ as.data.frame(cost.function(cost.volume(volume),aset))
+ }
+ }
+
+ # 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
+
+ # 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,...)
+
+ # Mark this as a known structure
+ class(aset)<-"highway.assignment.set"
+
+ return(aset)
+}
+
+# The following convenience functions make it easy to build a set of assignment classes, once
+# you know the parameters 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)
+# class.list<-add.assignment.class(class.list,aclass) # preserves name, unlike c(...)
+# aclass <- make.assignment.class(network,"Class.2",demand.2,links.2,penalty.subset.2,cost.function.2)
+# class.list<-add.assignment.class(class.list,aclass)
+# ## etc...
+
+make.assignment.class <- function( network, name, demand, link.subset=TRUE, penalty.subset=NULL, cost.function=NULL ) {
+ 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
+ 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 ) {
+ if ( all(match(aclass$name,names(classes),nomatch=0)==0) ) {
+ classes[[aclass$name]]<-aclass
+ } else {
+ warning("Cannot overwrite existing class",aclass$name)
+ }
+ return(classes)
+}
+
+# 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) {
+ # 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 ) )} ) )
+}
+
+# Intercept Set, describes select link processing
+
+# The intercept set is built at the level of the highway assignment
+ # intercept is an intercept set that will assemble the intercept link volumes
+ # typically, this will be applied to all groups and we'll end up with a demand
+ # array with Z x Z x numSets dimensions, and a links matrix with L x numSets
+ # rows and columns attached to the intercept set
+
+new.intercept.set <- function( links, filter.od=NULL ) {
+ # links is a logical or numeric vector indexing the network links -- any path between
+ # selected OD pairs that crosses one of these links is marked for inclusion by the
+ # intercept.paths function.
+ # filter.od takes a logical OD matrix and keeps only those OD pairs that are interesting
+ # (default is to keep all origins and destinations)
+ if ( mode(links)!="logical" ) stop("links parameter for intercept.set must be a logical vector")
+ if (is.null(filter.od)) filter.od<-function(od) return(od)
+ iset <- list(links=links,filter.od=filter.od)
+ attr(iset,"class")<-"intercept.set"
+ return(iset)
+}
+
Added: pre-alpha/travelr/R/Highway_Assignment.R
===================================================================
--- pre-alpha/travelr/R/Highway_Assignment.R (rev 0)
+++ pre-alpha/travelr/R/Highway_Assignment.R 2010-04-26 02:15:14 UTC (rev 9)
@@ -0,0 +1,280 @@
+####################################################################################################
+# travelr\r\Highway_Assignment.R by Jeremy Raw Copyright (C) 2010
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# A copy of the GNU General Public License is available at
+# http://www.r-project.org/Licenses/
+# and included in the R distribution (in directory ‘share/licenses’).
+####################################################################################################
+
+# This file defines the basic highway assignment mechanism and a series of algorithms, using the
+# single/multi-class assignment sets to provide a set of abstract capabilities for evaluating
+# the algorithms.
+
+highway.assign <- function( aset, method=c("AON","MSA","Frank.Wolfe","ParTan"), control=vector("list",0) ) {
+ # simple error-checking driver for highway assignment methods
+ # aset: data for the assignment (network, demand, cost function)
+ # method: Certain standard methods are defined in this module, and if the value of method (a character
+ # string) does not find a function with a name of "highway.assign.<<method>>", it will
+ # attempt to locate a built-in to perform the actual assignment (the built-ins are named
+ # ".highway.assign.<<method>>" and are hidden in the package namespace.
+ # control: a list of control parameters:
+ # intercept=NULL # an intercept.set object; intercept results will be returned, if supported by 'method'
+ # min.relative.gap=1e-4 # value below which algorithm has "converged"
+ # max.iter=100 # maximum number of assignment iterations
+ # max.time=600 # maximum algorithm runtime in seconds (for the impatient)
+ # opt.tol=.Machine$double.eps^0.5 # tolerance for internal line search
+ # verbose=0:5 # how much detail to report interactively during the run (bigger number ==> more wordy)
+ # log=FALSE # if true, accumulates and returns a data.frame of iteration statistics
+ # the implementation method should return a HighwayAssignment object (q.v.)
+
+ # First check the aset parameter
+ if (class(aset)!="highway.assignment.set")
+ stop("Must provide an assignment.set object as first parameter (aset)")
+ # Then find the appropriate method
+ if (length(method)>1) {
+ method<-method[1] # which means the default method is the first in the built-in list
+ cat(sprintf(gettext("Using default assignment method '%s'\n"),as.character(method)))
+ }
+ method<-paste("highway.assignment.",as.character(method),sep="") # user can override by defining a suitable function
+ a.func<-NULL
+ if ( exists(method,mode="function") ) {
+ a.func <- get(method,mode="function")
+ } else {
+ method<-paste(".",method,sep="")
+ if ( exists(method,mode="function") ) a.func<-get(method,mode="function")
+ }
+ if ( is.null(a.func) )
+ stop(paste(gettext("Could not find assignment method implementation"),sprintf("'%s'",method)))
+ return(c(a.func(aset,control),list(method=method)))
+}
+
+# Parse control variable, including global defaults
+# 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
+ # 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
+ )
+
+parse.control<-function( control, control.element, default.value=NULL ) {
+ element<-control[[control.element]]
+ if ( !is.null(element) ) {
+ result<-element
+ default<-FALSE
+ } else {
+ default<-TRUE
+ if ( !is.null(default.value) ) result<-default.value
+ else result<-control.defaults[[control.element]]
+ }
+ if ( !is.null(result) ) attr(result,"default")<-default # can't put attribute on NULL
+ return(result)
+}
+
+# Build a function to compute equilibrium statistics from user-supplied objective function
+# 'flow' is a scalar representing the sum of all demands
+# Computations based on Boyce, Ralevic-Dekic, and Bar-Gera, presentation at
+# 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 )
+ start.time=proc.time()["elapsed"]
+ 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), ... ) {
+ # 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)
+ 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)
+ 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,...) )
+ }
+}
+
+# To go with the equilibrium stats function, we need a function to test for convergence
+# This function is extensible, and will use any control parameters that have prefixes
+# 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))]
+ 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])
+ function(results) {
+ all( results[min.names] < test.min ) &&
+ all( results[max.names] > test.max )
+ }
+}
+
+# Built-In Assignment Methods
+
+# 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
+.highway.assignment.AON <- function(aset,control) {
+ verbose <- parse.control(control,"verbose")
+ if ( verbose>0 )
+ cat("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
+.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
+
+ verbose <- parse.control(control,"verbose")
+ if ( verbose>0 )
+ cat("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))
+ } else if ( verbose>0 ) {
+ cat(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))
+ 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
+ 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
+
+ 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)
+
+ # 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
+ 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))
+}
+
+# Pre-build the intercept function (to avoid if-test inside inner loop)
+build.intercept.function <- function( iset, aset ) {
+ if ( is.null(iset) ) f <- function(paths) NULL
+ else {
+ classes <- aset$classes # list of matrices, one for each assignment class
+ demand <- lapply(classes,function(x) x$demand)
+ cat("Demand while building intercept function:\n")
+ links <- iset$links # vector of links of interest
+ f <- function(paths) {
+ od <- intercept.paths(paths,links) # each OD here is just 0/1
+ od <- mapply( function(od,d)d*od, # Expand to list of OD matrices
+ od,demand,
+ USE.NAMES=TRUE, SIMPLIFY=FALSE )
+ volumes <- mapply(.load.paths,
+ paths,demand,
+ USE.NAMES=TRUE, SIMPLIFY=FALSE )
+ return( list(od=od,volumes=volumes) ) # OD is list of matrices, per assignment class
+ # volumes is list of volumes, per assignment class
+ }
+ }
+ return( f )
+}
+
+# Helper functions for updating equilibrium iterations
+
+# Weighted update functions combine iteration values
+
+weighted.update.diff <- function( weight, base, new.diff ) base + weight * new.diff
+
+weighted.update <- function( weight, base, new ) {
+ if ( is.null(base) ) result <- new # ignoring weight in this case
+ else result <- base*(1-weight) + weight * new
+ return( result )
+}
+
+inner.weighted.intercept.function <- function(b,n,weight) {
+ list( od = weighted.update(weight,b$od,n$od),
+ volumes = weighted.update(weight,b$volumes,n$volumes) )
+}
+
+weighted.update.intercept<-function( weight, base, new ) {
+ if ( is.null(base) || is.null(new) ) return(new)
+ return( mapply(base,new,inner.weighted.intercept.function,vol=,SIMPLIFY=FALSE,weight=weight) )
+}
+
+# line search helper functions
+# Do like this before starting any assignment iterations:
+# lambda.func<-make.lambda.func(cost.function,objective.function)
+# lambda.search<-make.lambda.search(lambda.func,opt.tol=1e-8)
+# Then inside the iterations, when you need to find a lambda value, do this:
+# search.result<-lambda.search(volume,difference)
+# updated.volume <- volume + search.result$lambda * diff
+# updated.objective <- search.result$objective
+# Look at the Frank-Wolfe or ParTan algorithms to see applications
+
+build.lambda.function<-function(cost.function,objective.function)
+ function( lambda, vol, diff ) objective.function(cost.function(vol+lambda*diff),diff)
+
+build.lambda.search<-function(lambda.func,opt.tol=.Machine$double.eps^0.25) {
+ function(vol,diff) {
+ search.result<-uniroot(lambda.func,interval=c(0,1),tol=opt.tol,vol=vol,diff=diff)
+ return(list(lambda=search.result$root,objective=search.result$f.root))
+ }
+}
Added: pre-alpha/travelr/R/Highway_Network.R
===================================================================
--- pre-alpha/travelr/R/Highway_Network.R (rev 0)
+++ pre-alpha/travelr/R/Highway_Network.R 2010-04-26 02:15:14 UTC (rev 9)
@@ -0,0 +1,145 @@
+####################################################################################################
+# travelr\r\Highway_Network.R by Jeremy Raw Copyright (C) 2010
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# A copy of the GNU General Public License is available at
+# http://www.r-project.org/Licenses/
+# and included in the R distribution (in directory ‘share/licenses’).
+####################################################################################################
+
+# Highway network functions
+
+# map.highway.nodes constructs a node map object for the network
+# map<-map.highway.net(links)
+# links$A.mapped<-map$map(links$A)
+# all(links$A==map$unmap(links$A.mapped))
+
+# links should be a dataset or matrix with one row per network link
+# nodes is a list of node identifiers (integer) used to
+# identify ends of the links
+
+# Link.fields identify the "From" and "To" node identifiers. These should
+# be scalar values suitable for use as a list index.
+
+# Why NOT just use factors for the node identifiers, and make the "mapped" node fields into
+# the integer factor code?
+
+map.highway.nodes <- function( nodes, numZones ) {
+
+ # Extract nodes in use from the network
+ if ( any(nodes[1:numZones]!=c(1:numZones)) ) stop("Must have all zones (1 to",numZones,") present in node list")
+ # Path building will be unhappy and issue a warning if there are destination zones
+ # that don't exist and thus can't be reached; plus you'll waste a lot memory and processing time
+ # building path trees that have nothing in them.
+ node.map<-sort(unique(c(1:numZones,nodes[nodes>numZones])))
+ base.index<-0
+
+ # Construct functions with suitable closure for node mapping and unmapping
+ map.env<-new.env(parent=parent.frame(2)) # skip this function's environment
+ assign("node.map",node.map,envir=map.env) # but preserve necessary values
+ assign("base.index",base.index,envir=map.env)
+ map<-function(n) match(n,node.map,nomatch=base.index)
+ unmap<-function(n) node.map[n]
+ environment(map)<-map.env
+ environment(unmap)<-map.env
+
+ nmap<-list(map=map,unmap=unmap)
+ attr(nmap,"numNodes")<-length(node.map)
+ attr(nmap,"maxRawNode")<-max(node.map)
+ attr(nmap,"class")<-"highway.node.map"
+ return(nmap)
+}
+summary.highway.node.map<-function(object,...) return(unlist(attributes(object)[c("numNodes","maxRawNode")]))
+print.highway.node.map<-function(x,...) print(summary(x))
+length.highway.node.map<-function(x) return(attr(x,"maxRawNode"))
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/travelr -r 9
More information about the Travelr-commits
mailing list