~o|16
~?!i&32768
~o=39
~$>end_of_copyright
~/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
~/Copyright (C) Heinz Spiess, CH-2558 Aegerten, 1989-99. All rights reserved.
~/
~/The right to use this macro is granted to all EMME/2 users, provided the
~/following conditions are met:
~/ 1) The macro cannot be sold for a fee (but it can be used and distributed
~/ without charge within consulting projects).
~/ 2) The user is aware that this macro is not a part of the EMME/2 software
~/ licence and there is no explicit or implied warranty or support
~/ provided with this macro.
~/ 3) This macro must not be modified without also changing its name.
~/ 4) The comments in this macros must not be removed and any additions or
~/ modification must be appropriately identified as such and give at least
~/ the date, the name and the reason of the modification.
~/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
~:end_of_copyright
~#
~# Revision history:
~# 99-07-14 Version 2.2:
~# - Correction of a small bug which slipped into version 2.1 and made
~# the macro crash when being used without turn counts. (This bug
~# was found and reported by Zvi Leve - Thank you, Zvi!)
~# - Correction of two minor bugs regarding the case when a turn attribute
~# is specified which does not contain any counts.
~#
~# 99-07-13 Version 2.1:
~# - Support for counts on turns added. In order to keep the calling
~# sequence compatible with earlier versions, the turn attribute
~# containing the counts on turns is also specified in the
~# argument, separated from the link count attribute by a "/".
~# - R^2 is now reported on screen and in the summary file. Note that the
~# formula for R^2 is the one without correction for the mean (no constant
~# term!), thus it is *not* the same as the R^2 given in module 2.43.
~# - The summary file is now named "demadj.sum". It no longer uses command
~# espace commands ~!..., but the new write file feature.
~# - The number of link and turn counts is displayed and a tests were added
~# to detect negative values in the count data.
~# - Customization information can now be passed to DEMADJ in the registers
~# t3-t7. This avoids the need to modify the macro itself - instead it
~# is simply called via a wrapper macro which set the registers accordingly
~# before calling macro DEMADJ.
~# - Scalars ms91 - ms94 are no longer needed, floating point registers
~# r91 - r94 are now used instead.
~# - The macro will no longer work with releases older than Rel.9.
~#
~/*****************************************************************************
~/********* DEMADJ 2.2 - DEMAND ADJUSTMENT USING GRADIENT METHOD ***********
~/********* Copyright (C) Heinz Spiess, CH-2558 Aegerten, 1989-99 ***********
~/
~/ This macro implements the gradient method for the O-D matrix adjustment
~/ problem. Using observed volumes on a subset of links and a starting O-D
~/ matrix, the model will adjust the matrix to better fit with the observed
~/ volumes. The "steepest descent" property of the gradient approach ensures
~/ that the O-D matrix is not changed more than necessary. In a nut shell,
~/ here is a summary of the computations performed by this macro at each
~/ gradient step:
~/ - auto assignment to get the link volumes 5.11/5.21
~/ - computation link derivatives and objective funct. 2.41
~/ - generate scattergram observed vs predicted volumes 2.43
~/ - addl options assignment to compute gradient matrix 5.11/5.21
~/ - multiply gradient matrix with weights (optional) 3.21
~/ - addl options assignment to assign gradient matrix 5.11/5.21
~/ - compute maximal gradient 3.21
~/ - compute optimal step length 2.41
~/ - update demand matrix and step counter 3.21
~/
~/ This method has shown to work very well and, using this macro, is very
~/ simple to apply. However, as all matrix adjustment methods, it is extremely
~/ important to be sure that the differences between observed and assigned
~/ volumes are indeed due to the demand matrix and not a result of a coding
~/ errors in the network, badly calibrated volume delay functions or wrong
~/ counts. Very careful analysis of the data before and after the adjustment
~/ is of utmost importance!!!
~/
~/ Reference: Spiess, H., "A GRADIENT APPROACH FOR THE O-D MATRIX ADJUSTMENT
~/ PROBLEM", Publication 693, CRT, University of Montreal, 1990.
~/
~/ Data requirements:
~/ @ltmp, @ptmp: The extra attributes @ltmp and @ptmp are used for
~/ temporary storage of the link, resp. turn volume differences.
~/ If these attributes already exist, their contents will be lost.
~/ ms90: This scalar contains the current adjustment step number. If it
~/ exists, it must be initialized to 0 before a new adjustment
~/ is started. This scalar will be incremented after each adjustment
~/ step. It is compared with the stopping criterion at each
~/ gradient step.
~/
~/ Other requirements:
~/ 1)This macro must be started at the main menu.
~/ 2)This macro needs a release 9.0 or later of the EMME/2 system.
~/ 3)This macro produces for each gradient step a plot file named "scatNN.plt"
~/ containing the observed vs. assigned scattergram.
~/ 4)For each gradient step a one-line summary is written to the file
~/ demadj.sum containing the value of the objective function Z, the maximum
~/ gradient M, the optimal step length L and the R^2 measure of fit.
~/ 5)The demand adjustments can be restricted to a subset of O-D pairs by
~/ specifying an adjustability matrix in register t4 which contains a 1
~/ if the demand of an OD-pair can be adjusted and 0 if it should not be
~/ changed. (Value of t4 can be defined in an upper level wrapper macro.)
~/ 6)If needed, optional weight factors for the counts can be specified
~/ as an arbitrary link expression in text registers t6/t7.
~/ (Values of t6/t7 can be defined in an upper level wrapper macro.)
~/ 7)The R^2 measure reported in the progress report and the summary file is
~/ based on a linear regression *without* a constant term.
~/
~/ Usage: ~ []
~/
~/ Macro Parameters:
~/ : Link attribute containing link counts (0=no count), e.g. "ul1".
~/ An optional turn attribute containing counts on turns can be
~/ specified optionally after "/" (no blanks!), e.g. "ul2/up2".
~/ If counts are on turns only, use "0" for the link attribute.
~/ : Demand matrix to be adjusted (e.g. "mf6"). Since the initial
~/ contents of this matrix will be overwritten with the resulting
~/ adjusted matrix, it is recommended to first do a copy.
~/ This matrix must not be write protected.
~/ : Temporary matrix of type mf which is used to hold the gradient
~/ direction matrix. Must not be write protected. If this matrix
~/ already exists, its original contents will be lost.
~/ : Number of iteration stopping criteria used in assignment.
~/ : Relative gap stopping criteria used in assignment.
~/ : Normalized gap stopping criteria used in assignment.
~/ : Number of gradient steps to be performed (adjustment
~/ iterations).
~/ : Optional argument. If set to 1, the first "half step" (reg.
~/ assignment to compute volumes) is skipped, assuming that the
~/ assignment was already performed. Use only when continuing
~/ with more step on the same adjustment (i.e. higher
~/ value).
~/
~x=%0%
~?x<7
~$STOP
~p=50
~?p<9
~+|~/EMME/2 Release %p% - this version of DEMADJ needs Release 9 or later!|~$STOP
~?!m=000
~+|~/Macro DEMADJ must be started from the main menu of EMME/2!|~$STOP
~?!q=0
~+|~/Macro DEMADJ must be started from the main menu of EMME/2!|~$STOP
~t1=%1%
~z=0
~:SPLIT-NEXT
~+|~t2=%%%t1.-%z%%%%|~t2=%%%t2.1%%%|~?t2=|~$SPLIT-DONE|~z+1|~?!t2=/|~$SPLIT-NEXT
~+|~t2=%%%t1.-%z%%%%|~z-1
~+|~t1=%%%t1.%z%%%%|~?z=0|~t1=|
~:SPLIT-DONE
~?t1=
~+|~/No link counts specified!|~$STOP
~z=1
~+|~t9=%%%t1.-%z%
~# Registers t3 to t7 contain various configuration parameters. Usually,
~# they do not need to be changed. Should a change ever become necessary,
~# please do not change the default in the macro itself, rather set the
~# corresponding registers in a wrapper macro which then calls DEMADJ,
~# e.g.:
~# ~t4=mf4.ne.0
~# ~t6=@lwght
~# ~t7=@pwght
~# ~||"
~#
~# Set assignment type to default if it is not defined in calling macro:
~?t5=
~t5=1
~# Registers t6/t7 contain optional link/turn weight factors for count posts.
~# By default these registers are empty, which implies that all count posts
~# have an implicit weight of 1. If needed t6 (t7) can be set to an arbitrary
~# link (turn) expression which provides positive weight factors for all count
~# post links (e.g. "t6=@cwght" if the weights are stored in the extra link
~# attribute @cwght).
~#
~# Set link weights to default value if they are not defined in calling macro:
~?t6=
~t6=
~# Set turn weights to default value if they are not defined in calling macro:
~?t7=
~t7=
~%
~/*****************************************************************************
~+|2.41|1|n|put(get(1)+((%t1%)>0))+put(get(2)+((%t1%)<0))|
~?e
~+|~/Invalid link counts (%t1%)!||q|~$STOP
~+|all|5|~r10=%%%g1%%%|~r11=%%%g2%%%|q
~/Link counts: %t1% (contains counts on %r10% links)
~?r11>0
~+|~/%r11% negative link counts detected!|~$STOP
~?t2=
~$>NO TURN COUNTS
~+|2.41|1|n|0*tpf+put(get(1)+((%t2%)>0))+put(get(2)+((%t1%)<0))|
~?e
~+|~/Invalid turn counts (%t2%)!||q|~$STOP
~+|all|all|5|~r12=%%%g1%%%|~r13=%%%g2%%%|q
~/Turn counts: %t2% (contains counts on %r12% turns)
~?r13>0
~+|~/%r13% negative turn counts detected!|~$STOP
~:NO TURN COUNTS
~x=%r10%
~x+%r12%
~/Current demand: %1%
~/Demand direction: %2%
~?!t4=
~/Adjustability factor: %t4%
~?!t6=
~/Countpost weight factors: %t6% (on links)
~?!t7=
~/ %t7% (on turn)
~/Assignment: %3% iterations, %4% % rel.gap, %5% norm.gap
~/Gradient steps: %ms90% - %6% (continue: %7%)
c=demadj %t1%/%t2% %1% %2% %3% %4% %5% %6% %7%
~x=%7% 0
~?x>0
~$CONTINUE STEP
3.12 /##### check if matrices exist and create them if necessary ##############
1 /create scalar ms90, if it does not already exist
ms90
~?e
~$>CHECK DPQ
step
matrix adjustment step counter
0
1 /create full matrix dpq, if it does not already exist
~:CHECK DPQ
%2%
~?e
~$>MATRIX CHECKS DONE
dpq-0
demand gradient
0
~:MATRIX CHECKS DONE
~?q=0
q
2.42 /##### check if extra attribute @ltmp exists
3 / remove @ltmp if it exists
@ltmp
~?e
~?q=1
yes
2
2
~?e
~+|q|~/No space for extra link attribute @ltmp!|~$STOP
@ltmp
temporary link attribute (DEMADJ.MAC)
0
~?t2=
~$>NO PTMP
3 / remove @ptmp if it exists
@ptmp
~?e
~?q=1
yes
2
5
~?e
~+|q|~/No space for extra turn attribute @ptmp!|~$STOP
@ptmp
temporary turn attribute (DEMADJ.MAC)
0
~:NO PTMP
q
~:NEXT STEP
~/
~/Gradient step %ms90%:
5.11 /##### prepare auto assignment to compute auto volumes ###################
~/.... Assign demand to get link and turn volumes
1 / fixed demand auto assignment
~?q=2
2 / new assignment
%t5% / type of assignment (time based or generalized cost)
1 / no additional volumes
%1% / demand
/ no auto occupancy matrix
/ no addl demand
/ no auto times saved
~?q=1
y / initialize auto volumes
%3% / max iterations
%4% / % gap
%5% / time diff.
5.21 /##### perform auto assignment to compute auto volumes ###################
~?b=1
2
2.41 /##### compute R^2 observed vs predicted volumes #########################
1 / network calculation - compute R^2
n
put(get(1)+%t3%*(%t1%>0)*(%t1%)*volau)+
put(get(2)+%t3%*(%t1%>0)*(%t1%)*(%t1%))+
put(get(3)+%t3%*(%t1%>0)*volau*volau)
all
4 / none
~?t2=
~$>NO TURN R2
1 / network calculation - compute turn gradient
n
put(get(1)+%t3%*(%t2%>0)*(%t2%)*pvolau)+
put(get(2)+%t3%*(%t2%>0)*(%t2%)*(%t2%))+
put(get(3)+%t3%*(%t2%>0)*pvolau*pvolau)
all
all
4 / none
~:NO TURN R2
~# g1=%g1% g2=%g2% g3=%g3%
~r2=%g1%
~r2*%g1%
~r2/%g2%
~r2/%g3%
~/.... R^2 observed vs predicted volumes is %r2.5%
q
2.41 /##### compute gradient and objective function ###########################
1 / network calculation - compute link gradient
y
@ltmp
y
link gradient - macro demadj (it %ms90%)
%t3%*(%t1%>0)*(%t1%-volau)
~?!t6=
*(%t6%)
all
4 / none
~?t2=
~$>NO TURN GRADIENT
1 / network calculation - compute turn gradient
y
@ptmp
y
turn gradient - macro demadj (it %ms90%)
%t3%*(%t2%>0)*(%t2%-pvolau)
~?!t7=
*(%t7%)
all
all
4 / none
~:NO TURN GRADIENT
1 / network calculation - compute objective function
n
put(get(puti(91))+(%t1%>0)*(%t1%-volau)*(%t1%-volau)*%t3%
~?!t6=
*(%t6%)
)
all
4 / no report
~r3=%g91%
~?t2=
~$>NO TURN OBJECTIVE
1 / network calculation - compute objective function
n
put(get(puti(91))+
(%t2%>0)*(%t2%-pvolau)*(%t2%-pvolau)*%t3%
~?!t7=
*(%t7%)
)
all
all
4 / no report
~:NO TURN OBJECTIVE
~r91=%g91%
q
~r4=%r91%
~r4-%r3%
~/.... Objective function at step %ms90% is %r91%
~?!t2=
~/ (%r3% on links, %r4% on turns)
on=1
plots=scat%ms90%.plt
2.43 /##### plot scattergram observed vs predicted volumes ####################
2 / link scattergram
%t1% / observed volumes vs
~?e
~+||~$>NO LINK SCATTERGRAM
volau / predicted volumes at step %ms90%
n
~?q=1
n / no symbol index
~?q=1
n / no color index
i=%ms90%
and 0
not %t1%=0
y / compute regression
~?q>0
~$>NO LINK SCATTERGRAM
0 / default x range maximum
0 / default y range maximum
~?q=2
2
~:NO LINK SCATTERGRAM
~?t2=
~$>NO TURN SCATTERGRAM
5 / turn scattergram
%t2% / observed volumes vs
~?e
~+||~$>NO TURN SCATTERGRAM
pvolau / predicted volumes at step %ms90%
n
~?q=1
n / no symbol index
~?q=1
n / no color index
i=%ms90%
or not 0
i=%ms90%
or not 0
.00001,99999999
y / compute regression
~?q>0
~$>NO TURN SCATTERGRAM
0 / default x range maximum
0 / default y range maximum
~?q=2
2
~:NO TURN SCATTERGRAM
q
off=1
plots=^
~:CONTINUE STEP
~x=%6%
~?x<%ms90%
~$STOP
5.11 /##### prepare add option assignment to compute gradient matrix ##########
~/.... Compute gradient matrix
1 / fixed demand auto assignment
~?q=2
2 / new assignment
%t5% / type of assignment (time based or generalized cost)
5 / additional options
%1% / demand matrix
/ no auto occupancy
/ no add. demand = auto demand
/ no auto times
~?q=1
y / initialize volumes
~?t2=
~+|6|@ltmp / addl path attribute just on links
~?!t2=
~+|7|@ltmp|@ptmp / addl path attribute on links and on turns
+ / addl path operator
-99999,99999 / addl path threshold
%2% / addl attribute matrix (4.3 and later)
y
D-%ms90% / matrix name
demand gradient step %ms90%
~?q=1
y / initialize matrix
0 / to zero
3 / save weighted addl path attributes in matrix
~?q=1
y / initialize volumes
%3% / number of iterations
%4% / relative gap
%5% / normalized gap
5.21 /##### perform add. opt. assignment to obtain gradient matrix ############
~?q=2
2
~?t4=
~$>No adjustability factor specified in register t4
3.21 /##### multiply gradient matrix with adjustability factor - if applicable
1 / matrix calculations
y / save result
%2% / gradient matrix
n / don't change header
%2%*(%t4%)
/ no constraint matrix
n / no submatrix
~?q=2
2 / report to print file
q
~:No adjustability factor specified in register t4
5.11 /##### prepare addl options assignment to get volume direction ###########
~/.... Assign gradient matrix to obtain volume direction
1 / fixed demand auto assignment
2 / new assignment
%t5% / type of assignment (time based or generalized cost)
5 / addl options
%1% / demand matrix
/ no auto occupancy
%2% / demand gradient matrix
/ no auto times
~?q=1
y / initialize auto volumes
5 / no addl link attribute
%3% / number of iterations
%4% / relative gap
%5% / normalized gap
5.21 /##### perform addl options assignment to get volume direction ###########
~?q=2
2
3.21 /##### compute maximum gradient ##########################################
1 / matrix calculation
n / save result
put(get(puti(92)).max.abs(%2%/%1%))
%1%
0,0,ex
n
.max.
.max.
~?q=2
2
~r92=%g92%
q
~/.... Maximum gradient is %r92%
2.41 /##### compute denominator of optimal step length ##################
1 / compute denominator of optimal step (link part)
n
put(get(puti(94))+(%t1%>0)*volad*volad/%r92%)
all
4
~?t2=
~$>NO TURNS 1
1 / compute denominator of optimal step (turn part)
n
put(get(puti(94))+(%t2%>0)*pvolad*pvolad/%r92%)
all
all
4
~:NO TURNS 1
1 / compute optimal step length now
n
put(get(puti(93))+(%t1%>0)*(%t1%-volau)*volad/get(94)
~?!t6=
*(%t6%)
)
all
4
~?t2=
~$>NO TURNS 2
1 / compute optimal step length now
n
put(get(puti(93))+
(%t2%>0)*(%t2%-pvolau)*pvolad/get(94)
~?!t7=
*(%t7%)
)
all
all
4
~:NO TURNS 2
~r93=%g93%
q
~/.... Optimal step length is %r93%
3.21 /##### update demand matrix ##############################################
1 / matrix calculations
y / save result
%1% / in
y / modify header
gpq%ms90%
Demand at step %ms90%
~?q=1
n
%1%+(%r93%.min.1)*%2%/%r92%
~/.... Update demand: %1% := %1%+(%r93%.min.1)*%2%/%r92%
n
~?q=2
2
~?m=321
q
3.21 /##### increment step counter ms90 #######################################
1 / matrix calculations
~>>demadj.sum
~"I%ms90%: Z%r91% M%r92% L%r93% (R^2=%r2%)
~>
y / save result
ms90 / in scalar 90
n / don't change header
ms90+1
~?q=2
2
~?m=321
q
c= I%ms90%: Z%r91% M%r92% L%r93%
~$NEXT STEP
~:STOP
~/*****************************************************************************
~o=6
~?b=2
q