New nested choice function added

README.md

@@ -11,7 +11,7 @@ Python modules for typical travel demand modeling calculations
 <p class=MsoNormal>License URI: https://opensource.org/licenses/MIT</p>
 
 <p class=MsoNormal>Description: These are some python scripts that can be used
-as modules for performing typicaly calculations in a travel demand model. The
+as modules for performing typical calculations in a travel demand model. The
 expected usage is to put the py file in the site packages or other directory
 and then to use the import command in python to load the functions in the
 modules. The documentation for usage of each of the functions is shown in the
@@ -36,6 +36,8 @@ trip matrices with multiple target production vectors and one attraction vector<
 set of friction matrices with multiple production vectors and one target
 attraction vector</p>
 
+<p class=MsoNormal>f) CalcGravityShadow : Implements attraction balancing by scaling attractions instead of furnessing flows, this method is more 'correct'</p>
+
 <p class=MsoNormal>Choice functions:</p>
 
 <p class=MsoNormal>a) CalcMultinomialChoice : Calculates a multinomial choice
@@ -48,9 +50,15 @@ probability given base utilities, current utilities and base proabilities</p>
 probabilities given dictionary with tree definition, matrix references and
 number of zones</p>
 
-<p class=MsoNormal>d) TODO: CalcNestedChoiceFlat : Calculate nested choice on
-flat array so it can be used for stuff like microsim ABM etc... c above can in
-general be easily modified for this</p>
+<p class=MsoNormal>d) CalcNestedChoiceFlat : Calculate nested choice on
+flat array so it can be used for stuff like microsim ABM etc. usage is same as c) 
+but inputs are flat arrays instead of square matrices and length of vector/s instead of 
+number of zones</p>
+
+<p class=MsoNormal>Matrix estimation (ODME):</p>
+
+<p class=MsoNormal>a) MatEstimateGradient : Performs synthetic matrix estimation using a least squares formulation. The solution algorithm is gradient descent (see Spiess, H., "A GRADIENT APPROACH FOR THE O-D MATRIX ADJUSTMENT PROBLEM", Publication 693, CRT, University of Montreal, 1990.) </p>
 
-<p class=MsoNormal>Have fun!</p>
+<p class=MsoNormal>If you use some of the components or code in this repo, please consider citing as shown below. Have fun!</p>
 
+<p class=MsoNormal>Joshi. C, python-tdm, (2015), GitHub repository, https://github.com/joshchea/python-tdm#python-tdm</p>

scripts/CalcDistribution.py

@@ -8,13 +8,14 @@
 #                  friction factor matrix (P and A should be balanced before usage, if not then A is scaled to P)
 #               d) CalcMultiFratar : Applies fratar model to given set of trip matrices with multiple target production vectors and one attraction vector
 #               e) CalcMultiDistribute : Applies gravity model to a given set of frication matrices with multiple production vectors and one target attraction vector
+#               f) CalcGravityShadow : Implements attraction balancing by scaling attractions instead of furnessing flows, this method is more 'correct'
+#
+#              **All input vectors are expected to be numpy arrays
 #
-#              **All input vectors are expected to be numpy arrays  
-#               
 # Author:      Chetan Joshi, Portland OR
 # Dependencies:numpy [www.numpy.org]
 # Created:     5/14/2015
-#              
+#
 # Copyright:   (c) Chetan Joshi 2015
 # Licence:     Permission is hereby granted, free of charge, to any person obtaining a copy
 #              of this software and associated documentation files (the "Software"), to deal
@@ -34,7 +35,7 @@
 #              OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 #              SOFTWARE.
 #--------------------------------------------------------------------------------------------------------------------------------------------------------------------#
-import numpy
+import numpy as np
 
 def CalcFratar(ProdA, AttrA, Trips1, maxIter = 10):
     '''Calculates fratar trip distribution
@@ -45,23 +46,23 @@ def CalcFratar(ProdA, AttrA, Trips1, maxIter = 10):
        Returns fratared trip table
     '''
     print 'Checking production, attraction balancing:'
-    sumP = sum(ProdA)
-    sumA = sum(AttrA)
+    sumP = ProdA.sum()
+    sumA = AttrA.sum()
     print 'Production: ', sumP
     print 'Attraction: ', sumA
-    if sumP <> sumA:
+    if sumP != sumA:
         print 'Productions and attractions do not balance, attractions will be scaled to productions!'
         AttrA = AttrA*(sumP/sumA)
     else:
         print 'Production, attraction balancing OK.'
     #Run 2D balancing --->
-    for balIter in xrange(0, maxIter):
+    for balIter in range(0, maxIter):
         ComputedProductions = Trips1.sum(1)
         ComputedProductions[ComputedProductions==0]=1
         OrigFac = (ProdA/ComputedProductions)
-        Trips1 = Trips1*OrigFac[:, numpy.newaxis]
+        Trips1 = Trips1*OrigFac[:, np.newaxis]
 
-        ComputedAttractions = Trips1.sum(0) 
+        ComputedAttractions = Trips1.sum(0)
         ComputedAttractions[ComputedAttractions==0]=1
         DestFac = (AttrA/ComputedAttractions)
         Trips1 = Trips1*DestFac
@@ -86,10 +87,10 @@ def CalcDoublyConstrained(ProdA, AttrA, F, maxIter = 10):
     maxIter (optional) = maximum iterations, default is 10
     Returns trip table
     '''
-    Trips1 = numpy.zeros((len(ProdA),len(ProdA)))
+    Trips1 = np.zeros((len(ProdA),len(ProdA)))
     print 'Checking production, attraction balancing:'
-    sumP = sum(ProdA)
-    sumA = sum(AttrA)
+    sumP = ProdA.sum()
+    sumA = AttrA.sum()
     print 'Production: ', sumP
     print 'Attraction: ', sumA
     if sumP <> sumA:
@@ -102,12 +103,12 @@ def CalcDoublyConstrained(ProdA, AttrA, F, maxIter = 10):
         AttrT = AttrA.copy()
         ProdT = ProdA.copy()
 
-    for balIter in xrange(0, maxIter):
-        for i in range(0,len(ProdA)):
-            Trips1[i,:] = ProdA[i]*AttrA*F[i,:]/max(0.000001, sum(AttrA * F[i,:]))
-        
+    for balIter in range(0, maxIter):
+        for i in range(0, len(ProdA)):
+            Trips1[i,:] = ProdA[i]*AttrA*F[i,:]/max(0.000001, np.sum(AttrA * F[i,:]))
+
         #Run 2D balancing --->
-        ComputedAttractions = Trips1.sum(0) 
+        ComputedAttractions = Trips1.sum(0)
         ComputedAttractions[ComputedAttractions==0]=1
         AttrA = AttrA*(AttrT/ComputedAttractions)
 
@@ -117,10 +118,53 @@ def CalcDoublyConstrained(ProdA, AttrA, F, maxIter = 10):
 
     for i in range(0,len(ProdA)):
             Trips1[i,:] = ProdA[i]*AttrA*F[i,:]/max(0.000001, sum(AttrA * F[i,:]))
-            
+
     return Trips1
 
-def CalcMultiFratar(Prods, Attr, TripMatrices, maxIter =10):
+def CalcGravityShadow(ProdA, AttrA, F, maxIter = 10):
+    '''Calculates doubly constrained trip distribution for a given friction factor matrix,
+    uses shadow pricing at attraction end
+    ProdA = Production array
+    AttrA = Attraction array (Target attractions)
+    F = Friction factor matrix
+    maxIter (optional) = maximum iterations, default is 10
+    Returns trip table
+    '''
+    T = np.zeros(F.shape)
+    AttrA = AttrA*ProdA.sum() / AttrA.sum()  #in case P and A totals don't match - balance A to P
+    AttrA[AttrA<0.000001] = 0.0001 #avoid divide by zero
+    Attr = AttrA.copy()
+    F[F<0.000001] = 0.0001
+
+    for k in range(maxIter):
+        if k > 0:
+            A_calc = T.sum(0)
+            print('A_calc:' + str(A_calc))
+            Attr = Attr * AttrA / A_calc
+
+        for i in range(ProdA.shape[0]):
+            T[i,:] = ProdA[i] * Attr * F[i, :] / (Attr * F[i, :]).sum()
+
+    return T
+
+def CalcGravity(P, A, F, maxIter=10):
+    '''A more vectorized verion of doubly constrinaed gravity model
+    P = Array of zone Productions
+    A = Array of zone Attractions | also Target attractions
+    F = Friction factor or transformed utility matrix
+    '''
+    T = A*F*P[:, np.newaxis]/np.maximum(np.sum(A*F, axis=1), 0.00001)[:, np.newaxis]
+    for i in range(maxIter):
+        cA = T.sum(0) #sum of calculated attractions
+        factor = np.where(A > 0, A/cA, 0)
+        F = factor*F
+        T = A*F*P[:, np.newaxis]/np.maximum(np.sum(A*F, axis=1), 0.00001)[:, np.newaxis]
+    cA = T.sum(0)
+    diff = np.absolute(A - cA).max()  #maximum absolute difference between target and calculated
+    print ('final max abs diff: {}'.format(diff))
+    return T
+
+def CalcMultiFratar(Prods, Attr, TripMatrices, maxIter=10):
     '''Applies fratar model to given set of trip matrices with target productions and one attraction vector
     Prods = Array of Productions (n production segments)
     AttrAtt = Array of Attraction ( 1 attraction segment)
@@ -130,25 +174,25 @@ def CalcMultiFratar(Prods, Attr, TripMatrices, maxIter =10):
     '''
     numZones = len(Attr)
     numTripMats = len(TripMatrices)
-    TripMatrices = numpy.zeros((numTripMats,numZones,numZones))
-   
+    TripMatrices = np.zeros((numTripMats,numZones,numZones))
+
     ProdOp = Prods.copy()
     AttrOp = Attr.copy()
 
     #Run 2D balancing --->
-    for Iter in xrange(0, maxIter):        
+    for Iter in range(0, maxIter):
         #ComputedAttractions = numpy.ones(numZones)
         ComputedAttractions = TripMatrices.sum(1).sum(0)
         ComputedAttractions[ComputedAttractions==0]=1
         DestFac = Attr/ComputedAttractions
-        
-        for k in xrange(0, len(numTripMats)):
+
+        for k in range(0, len(numTripMats)):
             TripMatrices[k]=TripMatrices[k]*DestFac
             ComputedProductions = TripMatrices[k].sum(1)
             ComputedProductions[ComputedProductions==0]=1
             OrigFac = Prods[:,k]/ComputedProductions #P[i, k1, k2, k3]...
-            TripMatrices[k]=TripMatrices[k]*OrigFac[:, numpy.newaxis]
-            
+            TripMatrices[k]=TripMatrices[k]*OrigFac[:, np.newaxis]
+
     return TripMatrices
 
 def CalcMultiDistribute(Prods, Attr, FricMatrices, maxIter = 10):
@@ -159,32 +203,26 @@ def CalcMultiDistribute(Prods, Attr, FricMatrices, maxIter = 10):
        Returns N-Dim array of trip matrices corresponding to each production segment
     '''
     numZones = len(Attr)
-    TripMatrices = numpy.zeros(FricMatrices.shape)
+    TripMatrices = np.zeros(FricMatrices.shape)
     numFricMats = len(FricMatrices)
 
     ProdOp = Prods.copy()
     AttrOp = Attr.copy()
+    #Initial trip distribution --->
+    for k in range(0, numFricMats):
+        for i in range(0, numZones):
+            if ProdOp[i, k] > 0:
+                TripMatrices[k, i, :] = ProdOp[i, k] * AttrOp * FricMatrices[k, i, :] / max(0.000001, np.sum(AttrOp * FricMatrices[k, i, :]))
 
-    for Iter in xrange(0, maxIter):
-        #Distribution --->
-        for k in xrange(0, numFricMats):
-            for i in xrange(0, numZones):
-                if ProdOp[i, k] > 0:
-                    TripMatrices[k, i, :] = ProdOp[i, k] * AttrOp * FricMatrices[k, i, :] / max(0.000001, sum(AttrOp * FricMatrices[k, i, :]))
-        #Balancing --->        
+    for Iter in range(0, maxIter):
+        #Balancing --->
         ComputedAttractions = TripMatrices.sum(1).sum(0)
         ComputedAttractions[ComputedAttractions==0]=1
         AttrOp = AttrOp*(Attr/ComputedAttractions)
-        for k in xrange(0, len(fricmatnos)):
-            ComputedProductions = TripMatrices[k].sum(1)
-            ComputedProductions[ComputedProductions==0]=1
-            OrigFac = Prods[:,k]/ComputedProductions
-            ProdOp[:,k] = OrigFac*ProdOp[:,k]
-    #Final Distribution --->
-    for k in xrange(0, numFricMats):
-        for i in xrange(0, numZones):
-            if ProdOp[i, k] > 0:
-                TripMatrices[k, i, :] = ProdOp[i, k] * AttrOp * FricMatrices[k, i, :] / max(0.000001, sum(AttrOp * FricMatrices[k, i, :]))
-
-    return TripMatrices
+        #Distribution --->
+        for k in range(0, numFricMats):
+            for i in range(0, numZones):
+                if ProdOp[i, k] > 0:
+                    TripMatrices[k, i, :] = ProdOp[i, k] * AttrOp * FricMatrices[k, i, :] / max(0.000001, np.sum(AttrOp * FricMatrices[k, i, :]))
 
+    return TripMatrices

scripts/CalcLogitChoice.py

@@ -3,8 +3,8 @@
 # Purpose:     Utilities for various calculations of different types of choice models.
 #               a) CalcMultinomialChoice : Calculates a multinomial choice model probability given a dictionary of mode utilities 
 #               b) CalcPivotPoint : Calculates pivot point choice probability given base utilities, current utilities and base proabilities
-#               e) CalcNestedChoice : Calculates n-level nested mode choice probabilities given dictionary with tree definition, matrix references and number of zones
-#               f) TODO: CalcNestedChoiceFlat : Calculate nested choice on flat array so it can be used for stuff like microsim ABM etc... e) can in general be easily modified for this 
+#               c) CalcNestedChoice : Calculates n-level nested mode choice probabilities given dictionary with tree definition, matrix references and number of zones
+#               d) CalcNestedChoiceFlat : Calculate nested choice on flat array so it can be used for stuff like microsim ABM etc... e) can in general be easily modified for this 
 #              **All input vectors are expected to be numpy arrays  
 #               
 # Author:      Chetan Joshi, Portland OR
@@ -61,12 +61,9 @@ def CalcMultinomialChoice(Utils, getLogSumAccess = 0):
     else:
         return Probs, lnSumAccess
 
-def CalcPivotPoint(BaseUtils, Utils, Po):
+def CalcPivotPoint(Utils, Po):
     '''
-       BaseUtils = Base utility matrices in a dictionary
-       ex. BaseUtils = {'auto':mat1, 'transit':mat2, 'bike':mat3, 'walk':mat4}
-
-       Utils = Updated utility matrices in a dictionary
+       Utils = Updated delta utility matrices in a dictionary i.e delta of Uk (k = mode)
        ex. Utils = {'auto':mat1, 'transit':mat2, 'bike':mat3, 'walk':mat4}
 
        Po = Base probabilities in a dictionary
@@ -162,11 +159,92 @@ def CalcNestedChoice(TreeDefn, MatRefs, numZn, getLogSumAccess = 0):
         return ProbMats
     else:
         return ProbMats, MatRefs['ROOT']
+
+
+def CalcNestedChoiceFlat(TreeDefn, MatRefs, vecLen, getLogSumAccess = 0):
+    '''
+    #TreeDefn = {(0,'ROOT'):[1.0,['AU', 'TR', 'AC']],
+    #            (1,'AU'):[0.992,['CD', 'CP']],
+    #            (1,'TR'):[0.992,['TB', 'TP']],
+    #            (1,'AC'):[0.992,['BK', 'WK']]}
+    #
+    #Key-> (Level ID, Level Code): Values-> (LogSum Parameter enters as: 1/lambda, SubLevel IDs)
+    #       ROOT should always be ID = 0 and Code = 'ROOT'                            
+    #                             ROOT
+    #                            / |  \
+    #                           /  |   \
+    #	                       /   |    \
+    #	                     AU    TR    AC(logsum parameter)
+    #	                    /\     /\     /\ 
+    #	                  CD CP  TB TP  BK WK         
+    #             		              
+    #MatRefs = {'ROOT': 1.0, 'AU':0, 'TR':0, 'AC':0,
+    #           'CD':Ucd), 'CP':Ucp),
+    #           'TB':Utb), 'TP':Utp),
+    #           'BK':Ubk), 'WK':Uwk)} Stores utilities in dict of vectors, base level utilities are pre-specified!!
+    #
+    #vecLen = number of od pairs being evaluated 
+    #
+    #getLogSumAccess (optional, accessibility log sum) 0=no, <>0=yes
+    ''' 
+    #ProbMats = {'ROOT': 1.0, 'AU':0, 'TR':0, 'AC':0, 'CD':0, 'CP':0, 'TB':0, 'TP':0, 'BK':0, 'WK':0}   #Stores probabilities at each level
+    #TripMat = GetMatrixRaw(Visum, tripmatno) #--> Input trip distribution matrix
+    #numZn = Visum.Net.Zones.Count
+    ProbMats = dict(zip(MatRefs.keys(), numpy.zeros(len(MatRefs.keys()))))
+    ProbMats['ROOT'] = 1.0
+    #Utility calculator going up...
+    #print 'Getting logsums and utilities...'
+    for key in sorted(TreeDefn.keys(), reverse= True):
+        #print key, TreeDefn[key]
+        sumExp = numpy.zeros(vecLen)   
+        sublevelmat_codes = TreeDefn[key][1]  #produces --> ex. ['WB', 'WX', 'DX']
+
+        for code in sublevelmat_codes:
+            #print ([code, TreeDefn[key][0]])
+            MatRefs[code] = MatRefs[code]/TreeDefn[key][0]  #---> scale the utility
+            sumExp+=numpy.exp(MatRefs[code])
+        lnSum = sumExp.copy() #Maybe there is a better way of doing the next 4 steps in 1 shot
+        lnSum[sumExp == 0] = 0.000000001 
+        lnSum = numpy.log(lnSum)
+        lnSum[sumExp == 0] = -999
+        MatRefs[key[1]] = TreeDefn[key][0]*lnSum #---> Get ln sum of sublevel
+
+    #Probability going down...
+    #print 'Getting probabilities...'
+    for key in sorted(TreeDefn.keys()):
+        #print key, TreeDefn[key]
+        eU_total = numpy.zeros(vecLen)
+        sublevelmat_codes = TreeDefn[key][1] #1st set--> ROOT : AU, TR
+        for code in sublevelmat_codes:
+            #print ([code, TreeDefn[key][0]])
+            eU_total+=numpy.exp(MatRefs[code])
+
+        eU_total[eU_total == 0] = 0.0001  #Avoid divide by 0 error
+##        for code in sublevelmat_codes:
+##            ProbMats[code] = ProbMats[key[1]]*numpy.exp(MatRefs[code])/eU_total
+        nSublevels = len(sublevelmat_codes)
+        cumProb = 0 
+        for i in xrange(nSublevels - 1):
+            code = sublevelmat_codes[i]
+            temp = numpy.exp(MatRefs[code])/eU_total
+            ProbMats[code] = ProbMats[key[1]]*temp
+            cumProb+=temp
+        code = sublevelmat_codes[i+1]
+        ProbMats[code] = ProbMats[key[1]]*(1.0-cumProb) 
+
+    if getLogSumAccess == 0:
+        return ProbMats
+    else:
+        return ProbMats, MatRefs['ROOT']
+
 #some generic utilities for reading and writing numpy arrays to disk..
 
 def GetMatrix(fn, numZn):
     return numpy.fromfile(fn).reshape((numZn, numZn))
 
+def GetMatrixFlat(fn):
+    return numpy.fromfile(fn)
+
 def PushMatrix(fn, mat):
     mat.tofile(fn)
 
@@ -205,4 +283,4 @@ def PushMatrix(fn, mat):
 ##print 'Calculation completed.'
 ##print 'Time taken(secs): ', time.time()-start   
 
-
+