~?!i&32768
~o=311
~$>end_of_copyright
~/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
~/Copyright (C) Heinz Spiess, CH-2558 Aegerten, 1991-2001. 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) The comments in this macros must not be removed and any additions or
~/      modification must be appropriately identified as such and give at least
~/      date, name and the reason of the modification.
~/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
~:end_of_copyright
~/************************************************************************
~/****   DEMADJT 1.8 - TRANSIT DEMAND ADJUSTMENT USING GRADIENT METHOD ***
~/****   Copyright (C) Heinz Spiess, CH-2558 Aegerten, 1991-2001       ***
~/
~/!!!! Not part of EMME/2 software licence - no support nor warranty.  !!!
~/ 
~/This macro implements the gradient method for the O-D matrix adjustment
~/problem.  Using observed volumes on a subset of segments 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:
~/    - transit assignment to get the segment volumes      5.11/5.31
~/    - computation segment derivatives and obj. funct.    2.41
~/    - generate scattergram observed vs predicted volumes 2.43
~/    - addl options assignment to compute gradient matrix 5.11/5.31
~/    - transit assignment to assign gradient matrix       5.11/5.31
~/    - 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:
~/@stmpx: Extra segment attribute reserved for temporary storage.  It must not
~/        yet exist, and it will be removed upon successful termination of the
~/        macro.
~/@ltmpx: Extra link attribute reserved for temporary storage.  It must not
~/        yet exist, and it will be removed upon successful termination of the
~/        macro.  Only used if link counts are specified.
~/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 <gstep> stopping criterion at each
~/        iteration.
~/ms91:   Scalar to hold the value of the objective function Z.
~/ms92:   Scalar to hold the  maximum gradient.
~/ms93:   Scalar to hold optimal step length.
~/ms94:   Scalar to hold temporary data.
~/
~/Other requirements:
~/1)This macro must be started at the main menu.
~/2)This macro needs a release 7.0 or later of the EMME/2 system.
~/3)A file "summary" is produced containing for each iteration
~/  a one line summary with the value of the objective function Z,
~/  the maximum gradient M and the step length L.  For other systems, the same
~/  information is written as comments into the log book.
~/4)The assignment parameters and weight factors are coded directly into
~/  macro.  They would need to be modified by editing the macro file itself.
~/5)For each gradient step <N>, a separate plot file scatt<N>.plt is created
~/  which contains the "observed-vs-assigned" scattergram.
~/
~/Usage: ~<demadjt <counts> <gpq> <dpq> <gstep>
~/
~/Macro Parameters:
~/<counts>:  Observed volumes (0=no count available), e.g. "us1"
~/           Can be either a segment or a link attribute.  If a link
~/           attribute is used, the segment gradient is computed by
~/           using the link difference for all lines using a link.
~/<gpq>:     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.
~/<dpq>:     Temporary matrix of type mf which is used to hold the gradient
~/           direction matrix.  Must exist and not be write protected.
~/<gstep>:   Number of gradient steps to be performed (adjustment
~/           iterations).
~/
~/*****************************************************************************
~x=%0%
~?x<4
~$STOP
~t2=%1%
~%
~t5=link
2.43 /##### test if counts are given as link or segment attribute
~+|2|%t2%|~?e|~t5=segment||q
~/%t5% counts stored in %t2% 
~/current demand stored in %1%
~/demand direction stored in %2%
~/number of gradient steps: %3%  (continue: %4%)
~>>summary
~"
~"----------------------------------------------------------------- 
~"Output generated by macro demadjt (1.8)
~"command line: ~<demadjt %t2% %1% %2% %3% %4%
~"%t5% counts stored in attribute %t2% 
~"current demand stored in matrix %1%
~"demand direction stored in matrix %2%
~"number of gradient steps: %3%
~"Iter: ObjFunction: MaxGradient: StepLength:
~>
c=demadjt %t2% %1% %2% %3% %4%
~t1=.001  
~/scale factor %t1% - change if needed to avoid output overflow (*******)
~/2.42 - create extra segment attribute @stmpx used for temporary storage
2.42 /##### create extra attribute @stmpx
2          /create attribute
4          /at segment level
~?e
~+|~/No space for extra segment attribute!|q|~$STOP
@stmpx     
~?e
~+|~/Extra attribute @stmpx exists already!||q|~$STOP
temporary attr. used by macro DEMADJT
0
~?t5=segment
~$>extra_attributes_created
2          /create attribute
2          /at link level
~?e
~+|~/No space for extra link attribute!|q|~$STOP
@ltmpx     
~?e
~+|~/Extra attribute @ltmpx exists already!||q|~$STOP
temporary attr. used by macro DEMADJT
0
~:extra_attributes_created
q
~/3.12 - check if matrices exist and create them if necessary 
3.12 /##### check if matrices exist and create them if necessary ##############
1          /create scalar ms90, if it does not already exist
ms90
~?e
~$CHECK MS91
step
matrix adjustment step counter
0
1          /create scalar ms91, if it does not already exist
~:CHECK MS91
ms91
~?e
~$CHECK MS92
Z-0
value of objective function Z
0
1          /create scalar ms92, if it does not already exist
~:CHECK MS92
ms92
~?e
~$CHECK MS93
M-0
max abs gradient
0
1          /create scalar ms93, if it does not already exist
~:CHECK MS93
ms93
~?e
~$CHECK MS94
L-0
maximum step length
0
1          /create scalar ms94, if it does not already exist
~:CHECK MS94
ms94
~?e
~$CHECK DPQ
temp
temporary scalar
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
~:NEXT STEP
~/------------------------- DEMADJT STEP %ms90% ------------------------------
~/5.11 - prepare transit assignment to compute transit volumes 
5.11 /##### prepare transit assignment to compute transit volumes ##########
2         / fixed demand transit assignment
~?q=2
2         / new assignment
%1%       / current demand matrix
          / no impedances
          / no in-vehicle time matrix
          / no auxiliary time matrix
          / no waiting time matrix
          / no 1st waiting time matrix
          / no boarding time matrix
          / no no. of boardings matrix
'*'       / active modes
1         / actual headways
1         / same boarding times
2.5       / minutes boarding time
1         / same wait time factor
0.5       / wait time factor
2.0       / wait time weight
2.0       / aux time weight
2.0       / boarding time weight
n         / no additional options
~/5.31 - perform transit assignment to compute transit volumes
5.31 /##### perform transit assignment to compute transit volumes ############
~?b=1
2
~?t5=link
~$>use_link_counts
~/2.41 - compute gradient and objective function
2.41 /##### compute gradient and objective function ###########################
1         / network calculation - compute segment gradient directly
y
@stmpx
y
Demadjt: observed-predicted iter%ms90%
%t1%*(%t2%>0)*(%t2%-voltr)

all
all
5         / store scalars (if faster, no pasue)
1         / back
1         / network calculation - compute objective function
n
(%t2%>0)*(%t2%-voltr)*(%t2%-voltr)*%t1%

all
all
5         / save summary results in scalars
4         / save sum of result values
91        / in scalar ms91
Z-%ms90%  / scalar name
'Dem adj. obj. function - step %ms90%'
q
~/       objective function is %ms91%
on=1
plots=scatt%ms90%.plt
~/2.43 - plot scattergram observed vs predicted volumes (scatt%ms90%.plt)
2.43 /##### plot scattergram observed vs predicted volumes ####################
4         / segment scattergram
%t2%      / observed volumes vs
voltr     / predicted volumes at step %ms90%
n
~?q=1
n         / no symbol index
~?q=1
n         / no color index
all
i=%ms90%
not 0
  
.001,9999999 /threshold - only segments with %t2%>0

y         / compute regression
0         / default x range minimum
0         / default y range minimum
~?q=2
2

q
off=1
plots=^
~$>compute_gradient_matrix
~:use_link_counts
~/2.41 - compute gradient and objective function
2.41 /##### compute gradient and objective function ###########################
1         / network calculation - compute total transit volumes on links
y
@ltmpx
y
link totals of transit segment volumes
voltr

4         / sum
all
all
5         / store scalars (if faster, no pasue)
1         / back
1         / network calculation - compute projected segment gradient
y
@stmpx
y
Demadjt: observed-predicted iter%ms90%
%t1%*(%t2%>0)*(%t2%-@ltmpx)

all
all
5         / store scalars (if faster, no pasue)
1         / back
1         / network calculation - compute objective function
n
(%t2%>0)*(%t2%-@ltmpx)*(%t2%-@ltmpx)*%t1%

all
5         / save summary results in scalars
4         / save sum of result values
91        / in scalar ms91
Z-%ms90%  / scalar name
'Dem adj. obj. function - step %ms90%'
q
~/       objective function is %ms91%
on=1
plots=scatt%ms90%.plt
~/2.43 - plot scattergram observed vs predicted volumes (scatt%ms90%.plt)
2.43 /##### plot scattergram observed vs predicted volumes ####################
2         / link scattergram
%t2%      / observed link volumes
@ltmpx    / predicted total transit volumes on links at step %ms90%
n
~?q=1
n         / no symbol index
~?q=1
n         / no color index
i=%ms90%
not 0
not %t2%=0

y         / compute regression
0         / default x range minimum
0         / default y range minimum
~?q=2
2

q
off=1
plots=^
~:compute_gradient_matrix
~x=%3%
~?x<%ms90%
~$DONE
~/5.11 - prepare add option assignment to compute gradient matrix 
5.11 /##### prepare add option assignment to compute gradient matrix ##########
2         / fixed demand transit assignment
~?q=2
2         / new assignment
%1%       / current demand matrix
          / no impedance matrix
          / no in-vehicle time matrix
          / no auxiliary time matrix
          / no waiting time matrix
          / no 1st waiting time matrix
          / no boarding time matrix
          / no no. of boardings matrix
'*'       / active modes
1         / actual headways
1         / same boarding times
2.5       / minutes boarding time
1         / same wait time factor
0.5       / wait time factor
2.0       / wait time weight
2.0       / aux time weight
2.0       / boarding time weight
y         / additional options
          / no additional boarding attribute
@stmpx    / additional segment attribute contains segment counts
          / no additional alighting attribute
          / no additional auxiliary transit attribute
1         / complete strategies
+         / path operator
          / default sub-strategy operator (average)
          / no strategy threshold
%2%       / addl attribute matrix  (patch! will be changed in 1.11)
y
D-%ms90%  / matrix name
demand gradient step %ms90%
~?q=1
y         / initialize matrix
0         / to zero
          / do not store active additional demand
~/5.31 - perform add. opt. assignment to obtain gradient matrix 
5.31 /##### perform add. opt. assignment to obtain gradient matrix ############
~?q=2
2         / Step %ms90%-B: %ms91% %ms92% %ms93%
~/3.21 - compute weighted gradient matrix
3.21 /##### 
1
y
%2%
n
%2%*(abs(%2%)<99999999)*%1%


n
~?q=2
2
~?m=321
q
~/5.11 - prepare assignment to get volume direction
5.11 /##### prepare assignment to get volume direction ###########
~?v<704
~+|q|1.11|3|1,-%s%,30,30|0|q|5.11 / reset MCADT to avoid ``transposition bug'' 
2         / fixed demand transit assignment
~?q=2
2         / new assignment
%2%       / current gradient matrix
          / no impedances
          / no in-vehicle time matrix
          / no auxiliary time matrix
          / no waiting time matrix
          / no 1st waiting time matrix
          / no boarding time matrix
          / no no. of boardings matrix
'*'       / active modes
1         / actual headways
1         / same boarding times
2.5       / minutes boarding time
1         / same wait time factor
0.5       / wait time factor
2.0       / wait time weight
2.0       / aux time weight
2.0       / boarding time weight
n         / no additional options
~/5.31 - perform assignment to get volume direction
5.31 /##### perform addl options assignment to get volume direction ###########
~?q=2
2         / Step %ms90%-C: %ms91% %ms92% %ms93%
~/3.21 - compute maximum gradient
3.21 /##### compute maximum gradient ##########################################
1         / matrix calculation
y         / save result
ms92      / in scalar 92
y
M-%ms90%  / scalar name
'maximum abs gradient step %ms90%' / and description
~?q=1
y
0
abs(%2%/%1%)

%1%
0,0,ex
n
.max.
.max.
~?q=2
2
~/       max gradient is %ms92%
~?m=321
q
~?t5=link
~$>step_link_counts
~/2.41 - compute optimal step length for segment counts
2.41 /##### compute optimal step length #######################################
1         / compute denominator of optimal step
n
(%t2%>0)*voltr*voltr/ms(92)

all
all
5
4
94        / save temporarily in scalar 94
T-%ms90%
SUM ( voltr*voltr ) .... step %ms90%
r
1         / compute optimal step length now
n
(%t2%>0)*(@stmpx/%t1%)*voltr/ms(94)

all
all
5
4
93        / save step length in scalar 93
L-%ms90%
Optimal step length .... step %ms90%
q
~$>step_computed
~:step_link_counts
~/2.41 - compute optimal step length for link counts
2.41 /##### compute optimal step length #######################################
1
y
tmpl1
voltr

4
all
all
5
r
1         / compute denominator of optimal step
n
(%t2%>0)*tmpl1*tmpl1/ms(92)

all
5
4
94        / save temporarily in scalar 94
T-%ms90%
SUM ( lvoltr*lvoltr ) .... step %ms90%
r
1         / compute optimal step length now
n
(%t2%>0)*(%t2%-@ltmpx)*tmpl1/ms(94)

all
5
4
93        / save step length in scalar 93
L-%ms90%
Optimal step length .... step %ms90%
q
~:step_computed
~/       step length is %ms93%
~/3.21 - update demand matrix 
3.21 /##### update demand matrix ##############################################
1         / matrix calculations
y         / save result 
%1%       / in <gpq>
y         / modify header
gpq%ms90%
Demand at step %ms90%
~?q=1
n
%1%+(ms93.min.1)*%2%/ms92
~/       %1% := %1%+(%ms93%.min.1)*%2%/%ms92%


n
~?q=2
2
~?m=321
q
~p=2004
~?p=1 (UNIX)
~>>summary
~"%ms90_5%   %ms91_>10%   %ms92_>10%  %ms93.6_>10%
~>
~/3.21 - increment step counter ms90 
3.21 /##### increment step counter ms90 #######################################
1         / matrix calculations
y         / save result
ms90      / in scalar 90
n         / don't change header
ms90+1


~?q=2
2
~?m=321
q
~?x>0
~$NEXT STEP
~:DONE
~>>summary
~"----------------------------------------------------------------- 
~"
~>
~/2.42 - delete extra attributes @stmpx/@ltmpx used for temporary storage
2.42 /##### delete extra attribute @stmpx and @ltmpx if applicable
~+|3|@stmpx|yes
~?t5=link
~+|3|@ltmpx|yes
q
off=1
~?b=2
q
~/************************************************************************
~/Macro 'DEMADJT %t2% %1% %2% %3% %4% ' terminated normally.
~/************************************************************************
~:STOP
~o=6